2019-04-11 17:22:57 Tifa_Best 阅读数 189

要求

  1. 均值滤波
    具体内容:利用 OpenCV 对灰度图像像素进行操作,分别利用算术均值滤波器. 几何均值滤波器. 谐波和逆谐波均值滤波器进行图像去噪。模板大小为5*5。(注:请分别为图像添加高斯噪声. 胡椒噪声. 盐噪声和椒盐噪声,并观察滤波效果)
  2. 中值滤波
    具体内容:利用 OpenCV 对灰度图像像素进行操作,分别利用 5*5 和 9*9尺寸的模板对图像进行中值滤波。(注:请分别为图像添加胡椒噪声. 盐噪声和椒盐噪声,并观察滤波效果)
  3. 自适应均值滤波。
    具体内容:利用 OpenCV 对灰度图像像素进行操作,设计自适应局部降低噪声滤波器去噪算法。模板大小 7*7 (对比该算法的效果和均值滤波器的效果)
  4. 自适应中值滤波
    具体内容:利用 OpenCV 对灰度图像像素进行操作,设计自适应中值滤波算法对椒盐图像进行去噪。模板大小 7*7(对比中值滤波器的效果)
  5. 彩色图像均值滤波
    具体内容:利用 OpenCV 对彩色图像 RGB 三个通道的像素进行操作,利用算术均值滤波器和几何均值滤波器进行彩色图像去噪。模板大小为 5*5。

过程

calcHist

使用的是


void cv::calcHist	(	const Mat * 	images,
int 	nimages,
const int * 	channels,
InputArray 	mask,
OutputArray 	hist,
int 	dims,
const int * 	histSize,
const float ** 	ranges,
bool 	uniform = true,
bool 	accumulate = false 
)

channels

当参数channels为1时,效果和np.histogram差不多。当channels为2时,效果和np.histogram2d差不多。

注意:当channels>1时,函数不是分别计算各个channel的一维直方图

ranges

range是个二维float数组。其中的各项(一维数组)指定了直方图各个维度的bin边界(取值范围)。当直方图中bin的边界是均匀分割得到的时候(unifor = true),对于每个维度i,只需要指定该维度的第一个bin(索引0)的下边界(包含在内)L0L_0和该维度的最后一个bin(索引为histSize[i]-1)的上边界UhistSize[i]1U_{\texttt{histSize}[i]-1}(排除在外)。在这种情况下,ranges的各项(一维数组)只需要包含2个元素(float)就行了。在另外的情况下(uniform=false),ranges[i]包含histSize[i]+1个元素L0,U0=L1,U1=L2,...,UhistSize[i]2=LhistSize[i]1,UhistSize[i]1L_0, U_0=L_1, U_1=L_2, ..., U_{\texttt{histSize[i]}-2}=L_{\texttt{histSize[i]}-1}, U_{\texttt{histSize[i]}-1}。图片中没有被包含在L0L_0UhistSize[i]1U_{\texttt{histSize[i]}-1}之间的,不会被统计进直方图中。

STL和Mat

vector是深复制。

高效的使用STL 当对象很大时,建立指针的容器而不是对象的容器

带你深入理解STL之迭代器和Traits技法

cppMat中对应的是

template<typename _Tp>
class cv::DataType< _Tp >
Template "trait" class for OpenCV primitive data types.

static local var

使用静态局部变量实现带状态的函数

nth_element

STL中用来寻找在全排序中位列第n个元素,实际上并不需要全排序,局部排序即可。应该是借鉴快排的思路了。

ROI

locateROI的函数用来寻找ROI(如果ROI是从Mat中截取出来的)在原Mat中的位置。

adjustROI的函数用来调整ROI的边界。

结构设计

我不是太清楚这些设计的正确命名

调用处理框架

特定函数,首先实现其特定的小函数,然后将特定小函数作为回调函数传递给并调用公共的处理函数(类似于处理框架)。主体是特定函数

函数表驱动

将特定函数部分实现后放进函数表里。公共处理部分(框架)主动调用函数表里的函数。主体为公共处理部分

如果将公共部分抽象出来成函数,将特定的小函数和公共处理函数封装成函数,就成了实验一中的特定函数。估计在公共处理部分只是整个程序的一小部分时才会这么做吧,比如说库。opencv的特定滤波器好像就这么实现的

代码

#include <opencv2/core.hpp>
#include <opencv2/imgproc.hpp>
#include <opencv2/highgui.hpp>

#include <numeric>
#include <iostream>
#include <vector>
#include <map>
#include <functional>
#include <string>

using namespace cv;
using namespace std;
using namespace std::placeholders;

float gauss_mean = 0, gauss_stdev = 15, pulse_prcnt = 0.01;

// line's point's x-aixs gap is depend on point number and image's width.
void plot_poly(Mat &imgIO, Mat &pts, const Scalar color, bool reset_scale=true)
{
    static float scale = 1.;
    Mat points;
    pts.convertTo(points, CV_32F);
    int width = imgIO.cols, height = imgIO.rows, pts_num = points.total();

    if (reset_scale) {
        double max_val;
        minMaxLoc(points, NULL, &max_val);
        scale = (double)height / max_val;
    }
    Mat tmpM(imgIO.size(), imgIO.type(), Scalar::all(0));

    for (int i = 1; i < pts_num; ++i)
        //Round operator should be done at last.
        line(tmpM,
            Point(cvRound((float)(i-1) / pts_num * width), height - cvRound(points.at<float>(i-1) * scale)),
            Point(cvRound((float)(i) / pts_num * width), height - cvRound(points.at<float>(i) * scale)),
            color, 2, LINE_AA, 0);

    imgIO += tmpM;
}

void add_gaussian_noise(const Mat & imgI, Mat & imgO, float mean, float stdev)
{
    cout << "Add gussian noise.\tmean:" << mean << "\tstdev:" << stdev << endl;
    Mat noise(imgI.size(), imgI.type());
    randn(noise, Scalar::all(mean), Scalar::all(stdev));
    add(imgI, noise, imgO);
}


void add_pulse_noise(const Mat &imgI, Mat &imgO, float percent, uchar value)
{
    CV_Assert(imgI.type() == CV_8U);

    cout << "Add pulse noise.\tpulse value:" << (int)value << endl;
    RNG rng(getTickCount());

    imgO = imgI.clone();
    int rows = imgI.rows, cols = imgI.cols;
    int counts = rows * cols * percent;
    for (int i = 0; i < counts; ++i)
        imgO.at<uchar>(rng.uniform(0, rows), rng.uniform(0, cols)) = value;
}

void filter2Dfun(const Mat &imgI, Mat &imgO, function<uchar(const Mat&)> pf, int m=3, int n=0)
{
    if (!n) n = m;
    cout << "\n\tsize(mxn):" << m << "x" << n;

    CV_Assert(m % 2 == 1);
    CV_Assert(n % 2 == 1);

    imgO = imgI.clone();
    // imgI.convertTo(imgO, CV_64F);

    int cols = imgI.cols, rows = imgI.rows;
    // top == (n - 1) / 2 == n / 2;  left == (m - 1) / 2 == m / 2;
    int top = n / 2, bottom = rows - top, left = m / 2, right = cols - left;
    for (int y = top; y < bottom; ++y) {
        for (int x = left; x < right; ++x) {
            Mat roi(imgI, Rect(x - left, y - top, m, n));
            imgO.at<uchar>(Point(x, y)) = pf(roi);
        }
    }
}

void apply_arithmetic_mean_filter(const Mat &imgI, Mat &imgO, int m=3, int n=0)
{
    if (!n) n = m;
    cout << "apply arithmetic mean filter.";
    auto fun = [=](const Mat &roi) -> uchar { return mean(roi)[0]; };
    filter2Dfun(imgI, imgO, fun, m, n);
    cout << endl;
}

void apply_geometric_mean_filter(const Mat &imgI, Mat &imgO, int m=3, int n=0)
{
    if (!n) n = m;
    cout << "apply geometric mean filter.";
    auto fun = [=](const Mat &roi) -> uchar {
        Mat tmpM;
        roi.convertTo(tmpM, CV_64F);
        // cv::pow(tmpM, (double)1/(m*n), tmpM);

// TODO: split 2-D into 1-D's may be better
        auto product = std::accumulate(tmpM.begin<double>(), tmpM.end<double>(), (double)1., [](double lhs, double rhs) -> double {
            if (!lhs) lhs = 1;
            if (!rhs) rhs = 1;
            return lhs * rhs;});
        product = pow(product, 1./m/n);
        // auto product = std::accumulate(tmpM.begin<double>(), tmpM.end<double>(), (double)1., multiplies<double>());
        return saturate_cast<uchar>(product);
    };
    filter2Dfun(imgI, imgO, fun, m, n);
    cout << endl;
}

void apply_harmonic_mean_filter(const Mat &imgI, Mat &imgO, int m=3, int n=0)
{
    if (!n) n = m;
    cout << "apply harmonic mean filter.";
    auto fun = [=](const Mat &roi) -> uchar {
        Mat tmpM;
        roi.convertTo(tmpM, CV_64F);
        divide(1, tmpM, tmpM);
        return saturate_cast<uchar>((double)m * n / sum(tmpM)[0]);
    };
    filter2Dfun(imgI, imgO, fun, m, n);
    cout << endl;
}

void apply_inverse_harmonic_mean_filter(const Mat &imgI, Mat &imgO, int Q, int m=3, int n=0)
{
    if (!n) n = m;
    cout << "apply inverse harmonic mean filter.\tQ:" << Q;
    auto fun = [=](const Mat &roi) -> uchar {
        Mat tmpM, tmpM_Qp1, tmpM_Q;
        roi.convertTo(tmpM, CV_64F);
        pow(tmpM, Q, tmpM_Q);
        pow(tmpM, Q+1, tmpM_Qp1);
        return saturate_cast<uchar>(sum(tmpM_Qp1)[0] / sum(tmpM_Q)[0]);
    };
    filter2Dfun(imgI, imgO, fun, m, n);
    cout << endl;
}

void apply_median_filter(const Mat &imgI, Mat &imgO, int m=3, int n=0)
{
    if (!n) n = m;
    cout << "apply median filter.";
    auto fun = [=](const Mat &roi) -> uchar {
        Mat tmpM = roi.clone();
        nth_element(tmpM.begin<uchar>(), tmpM.begin<uchar>() + tmpM.total()/2, tmpM.end<uchar>());
        return tmpM.at<uchar>(tmpM.total()/2);

    };
    filter2Dfun(imgI, imgO, fun, m, n);
    cout << endl;
}


void apply_adaptive_mean_filtering(const Mat &imgI, Mat &imgO, int m=3, int n=0)
{
    if (!n) n = m;
    cout << "apply adaptive mean filter.";
    float var_g = gauss_stdev;
    auto fun = [=](const Mat &roi) -> uchar {
        Mat mean_l, stdev_l;
        meanStdDev(roi, mean_l, stdev_l);
        double var_l = stdev_l.at<double>(0) * stdev_l.at<double>(0);
        if (var_g > var_l) {
            return saturate_cast<uchar>(mean_l.at<double>(0));
        } else {
            uchar val = roi.at<uchar>(roi.total()/2);
            return saturate_cast<uchar>(val - var_g / var_l * (val - mean_l.at<double>(0)));
        }
    };
    filter2Dfun(imgI, imgO, fun, m, n);
    cout << endl;
}

void apply_adaptive_median_filtering(const Mat &imgI, Mat &imgO, int S_max=9)
{
    int m = 3, r_max = (sqrt(S_max)-1)/2;
    cout << "apply adaptive median filter.";
    auto fun = [=](const Mat &roi0) -> uchar {
        uchar z_med, z_min, z_max, z_xy = roi0.at<uchar>(roi0.total()/2);
        Size sz;
        Point ofs;
        roi0.locateROI(sz, ofs);
        int w = sz.width, h = sz.height, cx = ofs.x + roi0.cols/2, cy = ofs.y + roi0.rows/2;
        Mat roi = roi0;

        // cout << endl << "cx:" << cx << "\tcy:" << cy << endl;

        for (int r = 1; r<=cx && r<=cy && r<w-cx && r<h-cy && r<=r_max; ++r) {
            auto z_minmax = minmax_element(roi.begin<uchar>(),roi.end<uchar>());
            z_min = *z_minmax.first, z_max = *z_minmax.second;

            Mat tmpM = roi.clone();
            nth_element(tmpM.begin<uchar>(), tmpM.begin<uchar>() + tmpM.total()/2, tmpM.end<uchar>());
            z_med =  tmpM.at<uchar>(tmpM.total()/2);

            if (z_min < z_med && z_med < z_max) {

                // cout << "return form program B." << endl;

                return (z_min < z_xy  && z_xy < z_max) ? z_xy : z_med;
            } else {
                roi.adjustROI(1, 1, 1, 1);

                // cout << "roi.size():" << roi.size() << endl;

            }
        }

        // cout << "size can't be bigger.return from program A." << endl;

        return z_med;
    };
    filter2Dfun(imgI, imgO, fun, m);
    cout << "~" << r_max *2 + 1 << "x" << r_max * 2 + 1 << endl;
}

void update_display(vector<string> wndwnms, vector<Mat*> images, string hist_wn="hist_img")
{
    namedWindow(hist_wn);
    Mat hist_img(images[0]->rows, images[0]->cols, CV_8UC3, Scalar::all(0));

    // origin image pdf, noise image pdf, processed image pdf
    vector<Mat> pdfs(3);
    float range[] = {0, 256};
    const float * hist_range = {range};
    int hist_size = 256;

    // cout << images[1] << endl;
    for (int i = 0; i < 3; ++i) {
        calcHist(images[i], 1, 0, Mat(), pdfs[i], 1, &hist_size, &hist_range);
        // plot_poly(hist_img,pdfs[i],Scalar(255*(0==i),255*(1==i),255*(2==i)));
        plot_poly(hist_img, pdfs[i], Scalar(255*(0==i),255*(1==i),255*(2==i)), i==0);
       
        imshow(wndwnms[i], *images[i]);
    }

    imshow(hist_wn, hist_img);
    waitKey();
    destroyWindow(hist_wn);
}

int main(int argc, char const *argv[])
{
    const char *image_path = (argc > 1) ? argv[1] : "img_test.jpg";

    Mat img, image, image_color, image_gray, image_noise, image_dst;
    img = imread(image_path, IMREAD_COLOR);
    if (img.empty()) {
        cout << "Cannot read image: " << image_path << std::endl;
        return -1;
    }

    image_color = img;
    cvtColor(image_color, image_gray, CV_BGR2GRAY);
    image = image_gray;

    vector<Mat*> imgs {&image, &image_noise, &image_dst};
    vector<string> wndnms {"image_src", "image_noise", "image_processed"};
    for (auto it : wndnms)
        namedWindow(it, CV_WINDOW_AUTOSIZE);

    vector<function<void(const Mat&, Mat&)>> add_noise_funs {
        bind(add_gaussian_noise, _1, _2, gauss_mean, gauss_stdev), 
        bind(add_pulse_noise, _1, _2, pulse_prcnt, 0), 
        bind(add_pulse_noise, _1, _2, pulse_prcnt, 255), 
        // TODO: [add_pulse_noise]
        [&](const Mat& imgI, Mat& imgO) -> void {
            add_pulse_noise(imgI, imgO, pulse_prcnt, 0); 
            add_pulse_noise(imgO, imgO, pulse_prcnt, 255);}
    };
    
    vector<function<void(const Mat&, Mat&)>> apply_mean_filters {
        bind(apply_arithmetic_mean_filter, _1, _2, 5, 5),
        bind(apply_geometric_mean_filter, _1, _2, 5, 5),
        bind(apply_harmonic_mean_filter, _1, _2, 5, 5),
        bind(apply_inverse_harmonic_mean_filter, _1, _2, 1, 5, 5)
    };


    for (auto apply_mean_filter : apply_mean_filters) {
        for (auto add_noise : add_noise_funs) {
            add_noise(image, image_noise);
            apply_mean_filter(image_noise, image_dst);
            update_display(wndnms, imgs);
            cout << endl;
        }
    }


    for (size_t i = 1; i < add_noise_funs.size(); ++i) {
        add_noise_funs[i](image, image_noise);

        apply_median_filter(image_noise, image_dst, 5);
        update_display(wndnms, imgs);
        cout << endl;

        apply_median_filter(image_noise, image_dst, 9);
        update_display(wndnms, imgs);
        cout << endl;
    }


    add_gaussian_noise(image, image_noise, gauss_mean, gauss_stdev);
    apply_arithmetic_mean_filter(image_noise, image_dst);
    update_display(wndnms, imgs);
    apply_adaptive_mean_filtering(image_noise, image_dst, 7);
    update_display(wndnms, imgs);
    cout << endl;


    add_pulse_noise(image, image_noise, pulse_prcnt, 255);
    add_pulse_noise(image_noise, image_noise, pulse_prcnt, 0);

    apply_median_filter(image_noise, image_dst, 7);
    update_display(wndnms, imgs);

    apply_adaptive_median_filtering(image_noise, image_dst, 7*7);
    update_display(wndnms, imgs);
    cout << endl;


    image = image_color;
    vector<Mat> bgr_planes;
    split(image, bgr_planes);
    for (auto &plane : bgr_planes) {
        add_gaussian_noise(plane, plane, gauss_mean, gauss_stdev);
        add_pulse_noise(plane, plane, pulse_prcnt, 0);
        add_pulse_noise(plane, plane, pulse_prcnt, 255);
    }
    merge(bgr_planes, image_noise);

    split(image_noise, bgr_planes);
    for(auto &plane : bgr_planes) {
        apply_arithmetic_mean_filter(plane, plane, 5);
    }
    merge(bgr_planes, image_dst);
    for (size_t i = 0; i < imgs.size(); ++i)
        imshow(wndnms[i], *imgs[i]);
    waitKey();

    split(image_noise, bgr_planes);
    for(auto &plane : bgr_planes) {
        apply_geometric_mean_filter(plane, plane, 5);
    }
    merge(bgr_planes, image_dst);
    for (size_t i = 0; i < imgs.size(); ++i)
        imshow(wndnms[i], *imgs[i]);
    waitKey();


    return 0;
}
2016-05-29 16:12:26 EbowTang 阅读数 4059

叠加平均去噪:反复对同一张干净的图像添加随机噪声(椒盐,乘性,高斯等噪声,此处噪声不叠加,具体见代码),进行200张噪声图像的平均。

小波阈值去噪:利用改进了阈值的小波阈值去噪对图像进行去噪处理,其中半阈值是改进了的阈值函数。

图像信息:经典的lena图,大小256*256,256位灰度图

噪声浓度均为0.05


1,添加'speckle'噪声时,两种算法的去砸效果对比


利用客观参数psnr评估可知:

pnsr_averge =40.9679
pnsr_noise =18.9210
pnsr_hard =23.9319
pnsr_soft =24.4425
pnsr_half =25.0121

结论1,清晰灰度图像下的乘性噪声污染,不管是主观还是客观参数,都说明叠加平均的去噪能力远优于小波去噪,叠加平均效果非常好。


2,添加'salt & pepper'噪声时,两种算法的去砸效果对比


pnsr_averge =38.0104
pnsr_noise =18.3476
pnsr_hard =20.2537
pnsr_soft =24.0672
pnsr_half =22.5696

结论2,清晰灰度图像下的椒盐噪声污染,不管是主观还是客观参数,都说明叠加平均的去噪能力远优于小波去噪,叠加平均效果非常好。


3,添加'gaussian'噪声时,两种算法的去砸效果对比


pnsr_averge =25.9784
pnsr_noise =19.0491
pnsr_hard =22.6575
pnsr_soft =22.1231
pnsr_half =23.0135


结论3:结论还是一致的,叠加平均效果非常好。


4,添加浓度(0.1)更高更复杂的混合(三种噪声的混合)噪声,两种算法的去砸效果对比


pnsr_averge =20.7570
pnsr_noise =12.5467
pnsr_hard =18.3766
pnsr_soft =18.7248
pnsr_half =18.4966



最终结论:

本文简单地以小波阈值去噪为辅助参考效果,并采用客观评价指标得出叠加平均去噪法对于随机噪声具有非常强的去噪能力。该方法适合对同一景物的多次采集图像相加后取平均值来消除噪声。在实际运算过程中不需要图像的统计特性,可以克服线性滤波器所带来的图像细节模糊,而且对随机的滤波脉冲干扰和散斑图像噪声最有效。设想,如果同一景物的多次采集图像由于不可抗拒因素存在微小的位移,那么此时就需要用配准来校准了。



参考代码:

clc;
clear;
close all;

src_img=imread('lena.bmp');
%I=rgb2gray(I);
figure(1);
subplot(231)
imshow(src_img);title('原图像');  
[a,b]=size(src_img);
dst_img=zeros(a,b);

var=0.05;
%noise_type='salt & pepper';
noise_type='speckle';
%noise_type='gaussian';
img_noise=imnoise(src_img,noise_type,var);

%%%%%%%%进行叠加平均%%%%%%%%%%%%%%%%%%%%
subplot(232)
imshow(img_noise);title('噪声图像');
for n=1:200
    img_noise=imnoise(src_img,noise_type,var);
    dst_img=dst_img+double(img_noise);%求和
end;
dst_img=dst_img/n;%取平均
%求取叠加平均除燥后的psnr结果
pnsr_averge = myPsnr(dst_img,src_img)
%求取加燥后的psnr结果
pnsr_noise = myPsnr(img_noise,src_img)

dst_img=uint8(dst_img);
subplot(233)
imshow(dst_img);title('叠加平均');

%%%%%%%%进行小波去噪%%%%%%%%%%%%%%%%%%%%
% 获取输入参数    
w    = 'db6';%小波类型    
n    = 3;%分解层数     
  
sorh1 = 'hard';%硬阈值    
sorh2 = 'soft';%软阈值    
sorh3 = 'half';%半阈值     
% 对图像进行小波分解    
[c,l] = wavedec2(img_noise,n,w);    
    
  
%求取阈值    
N = numel(src_img);    
[chd1,cvd1,cdd1] = detcoef2('all',c,l,1);     
cdd1=cdd1(:)';    
sigma = median(abs(cdd1))/0.6745;%提取细节系数求中值并除以0.6745    
thr = sigma*sqrt(2*log(N))/1.414; %对阈值做了改进  
  
  
% 对小波系数全局阈值处理    
cxchard = c;% 保留近似系数    
cxcsoft = c;% 保留近似系数    
cxchalf = c;% 保留近似系数    
justdet = prod(l(1,:))+1:length(c);%截取细节系数(不处理近似系数)    
    
% 阈值处理细节系数    
cxchard(justdet) = myWthresh(cxchard(justdet),sorh1,thr);%硬阈值去噪    
cxcsoft(justdet) = myWthresh(cxcsoft(justdet),sorh2,thr);%软阈值去噪    
cxchalf(justdet) = myWthresh(cxchalf(justdet),sorh3,thr);%软阈值去噪    
  
%小波重建    
xchard = waverec2(cxchard,l,w);    
xcsoft = waverec2(cxcsoft,l,w);    
xchalf = waverec2(cxchalf,l,w);     
  
subplot(234)   
imshow(uint8(xchard));title('硬阈值去噪图像');      
pnsr_hard = myPsnr(xchard,src_img)

subplot(235)   
imshow(uint8(xcsoft));title('软阈值去噪图像');     
pnsr_soft = myPsnr(xcsoft,src_img)

subplot(236)   
imshow(uint8(xchalf));title('半阈值去噪图像');   
pnsr_half = myPsnr(xchalf,src_img)


pnsr程序

function pnsr_result = myPsnr(img_ref,img_den)
[a,b]=size(img_ref);
XX2=double(img_ref) - double(img_den);
mse_value2 = sum(sum( XX2.^2 ))/(a*b);
pnsr_result = 10*log10( 255*255 / mse_value2 );


2014-10-26 21:53:38 EbowTang 阅读数 14952

本文代码的实现严重依赖前面的两篇文章:


1,一维信号的小波阈值去噪


2,小波变换一维Mallat算法的C++实现


注本文的大部分文字提取于参考论文
重要说明:本文代码是学习小波变换时写的,文中的代码有较严重的性能问题(但是运算结果是正确的),如你需要本代码,需要自行优化或者更改(一维阈值去噪那篇文章中的性能挺快的)

一,小波阈值去噪基本理论

      图像在获取或传输过程中会因各种噪声的干扰使质量下降,这将对后续图像的处理产生不利影响.所以必须对图像进行去噪处理,而去噪所要达到的目的就是在较好去除噪声的基础上,良好的保持图像的边缘等重要细节.在图像去噪领域得到广泛的应用.本博文根据小波的分解与重构原理,实现了基于硬阈值和软阈值函数的小波阈值去噪的C++版本,最终结果与matlab库函数运算结果完全一致。

1,小波阈值处理

小波阈值收缩法是Donoho和Johnstone提出的,一下便是养活不少学者的三篇基础论文:

【1】 Donoho D L. De-noising by soft-thresholding. IEEE Trans- actions on Information Theory, 1995, 41(3): 613−627
【2】 Donoho D L, Johnstone I M. Adapting to unknown smooth- ness via wavelet shrinkage. Journal of the American Statistic
Association, 1995, 90(432): 1200−1224
【3】 Donoho D L, Johnstone I M, Kerkyacharian G, Picard D. Wavelet shrinkage: asymptopia? Journal of Royal Statisti-
cal Society Series B, 1995, 57(2): 301−369

小波阈值去噪其主要理论依据是,小波变换具有很强的去数据相关性,它能够使信号的能量在小波域集中在一些大的小波系数中;而噪声的能量却分布于整个小波域内.因此,经小波分解后,信号的小波系数幅值要大于噪声的系数幅值.可以认为,幅值比较大的小波系数一般以信号为主,而幅值比较小的系数在很大程度上是噪声.于是,采用阈值的办法可以把信号系数保留,而使大部分噪声系数减小至零.小波阈值收缩法去噪的具体处理过程为:将含噪信号在各尺度上进行小波分解,设定一个阈值,幅值低于该阈值的小波系数置为0,高于该阈值的小波系数或者完全保留,或者做相应的“收缩(shrinkage)”处理.最后将处理后获得的小波系数用逆小波变换进行重构,得到去噪后的图像.


2,阈值函数的选取

  阈值去噪中,阈值函数体现了对超过和低于阈值的小波系数不同处理策略,是阈值去噪中关键的一步。设w表示小波系数,T为给定阈值,sign(*)为符号函数,常见的阈值函数有:


硬阈值函数:     (小波系数的绝对值低于阈值的置零,高于的保留不变)

     

软阈值函数:   (小波系数的绝对值低于阈值的置零,高于的系数shrinkage处理)

      


值得注意的是:

1) 硬阈值函数在阈值点是不连续的,在下图中已经用黑线标出。不连续会带来振铃,伪吉布斯效应等。

2) 软阈值函数,原系数和分解得到的小波系数总存在着恒定的偏差,这将影响重构的精度使得重构图像的边缘模糊等现象.

所以这里也添加一种简单的改进阈值函数,我们称之为半阈值

注意:图片中的公式有误,应该是sgn(w)*(|w|-P*T)(即将软函数中的阈值T缩小):

         

三种阈值处理策略见下图:



其实不少文章出现各种优秀的改进方案(于是有养活了不少学者,非本文重点):


3,阈值的确定

选取的阈值最好刚好大于噪声的最大水平,可以证明的是噪声的最大限度以非常高的概率低于,此阈值是Donoho和Johnstone提出的,其实我一直很想吐槽为什么阈值和信号的size相关呢?当然我的疑问也是大家的疑问,此问题有养活了一批学者,其中根号右边的这个参数(叫做sigma)就是估计出来的噪声标准偏差(第一级分解出的小波细节系数,即整个HH系数的绝对值的中值),本文将用此阈值去处理各尺度上的细节系数。


4,阈值策略

阈值策略有两种,一种是全局阈值策略,一种是分层阈值策略,从读论文了解到,全局阈值策略有他的缺陷:如果采用全局阈值进行处理,则会对所有尺度下的高频系数进行同一阈值处理,然而随着小波分解尺度的增加,有用信息分解后的小波系数将会越来越大,而噪声系数却会越来越小,所以为了防止高尺度下有用信息的分解系数被过度处理,分层阈值处理将会更合理。但是从我实际测试的结果来看------没什么区别!!!简便起见代码还是采用全局阈值。


为了不误导小伙伴们,请注意:这张图中阈值的估算不是来自CV1而是来自CD1,因为CD1是两次高通滤波的结果,里面含有最多的噪声,估计出来的噪声偏差也是最准确的。





5,小波阈值去噪过程


   有用信号经小波变换后, 其能量将集中在少数的小波系数上, 而噪声点的小波系数互不相关, 分布在各个尺度的所有时间轴上. 保留小波变换的各尺度下的模极大值点, 而将其他点置零或最大程度的减小, 然后将处理后的小波系数做小波逆变换, 即可达到抑制噪声的目的.阈值去噪是通过对变换域系数与阈值进行比较判断, 然后将处理后的系数进行逆变换重构去噪图像. 小波阈值去噪法的具体步骤如下:

步骤 1. 图像的小波分解: 确定小波函数和分解层次 N, 对图像进行 N 层的小波分解;
步骤 2. 阈值处理: 对分解得到的各层系数选择阈值, 并对细节系数进行阈值判断;
步骤 3. 图像重构: 对阈值处理后的系数通过小波逆变换重建图像.
   信号和噪声在小波域内具有不同的相关性. 信号在尺度间相应位置上的小波系数具有很强的相关性, 而噪声的小波系数则具有弱相关性或者不相关. 在阈值去噪中, 由于所选定的阈值通常固定, 不会随着小波系数的不同而变化, 这就不可避免地会对部分小波系数进行误判,于是又养活了一批学者。。。分析改进,分析改进,发文章!


1),小波的分解





2,小波的重构





二,Matlab库函数实现

1,核心库函数说明

1)wavedec2

图像的多级小波分解,将返回分解出来的小波系数以及小波系数的各级长度

2)waverec2

多级小波系数的重构,重构出原信号

3)wthresh函数

对系数进行指定类型(全局阈值或者分层阈值)的阈值去噪处理,软硬阈值处理函数。下图所示程序和运行结果可以比较清晰地看出该程序的执行过程。

clear;
y = linspace(-1,1,100);
thr = 0.4;
ythard = wthresh(y,'h',thr);
ytsoft = wthresh(y,'s',thr);
subplot(311);plot(y);
subplot(312);plot(ythard);
subplot(313);plot(ytsoft);



更具体的函数说明可以在,matlab里键入“doc 函数名”将得到很详细的说明,当然也可以百度

2,软,硬阈值处理效果:



哎哟,不错哦,半阈值去噪效果在视觉上的确有明显的改进效果




局部放大图像:

四幅图象均放大两倍,便于查看区别

这幅图像是源图像放大的效果:









3,完整的代码实现

说明:代码实现对图像一层分解后的细节系数进行全局阈值处理(即LHL,LH,HH均用同一阈值处理)。

% 获取输入参数  
w    = 'db3';%小波类型  
n    = 3;%分解层数   

sorh1 = 'hard';%硬阈值  
sorh2 = 'soft';%软阈值  
sorh3 = 'half';%半阈值   
% 对图像进行小波分解  
[c,l] = wavedec2(octimage,n,w);  
  

%求取阈值  
N = numel(octimage);  
[chd1,cvd1,cdd1] = detcoef2('all',c,l,1);   
cdd1=cdd1(:)';  
sigma = median(abs(cdd1))/0.6745;%提取细节系数求中值并除以0.6745  
thr = sigma*sqrt(2*log(N))/sqrt(1+sqrt(n)); %对阈值做了改进


% 对小波系数全局阈值处理  
cxchard = c;% 保留近似系数  
cxcsoft = c;% 保留近似系数  
cxchalf = c;% 保留近似系数  
justdet = prod(l(1,:))+1:length(c);%截取细节系数(不处理近似系数)  
  
% 阈值处理细节系数  
cxchard(justdet) = myWthresh(cxchard(justdet),sorh1,thr);%硬阈值去噪  
cxcsoft(justdet) = myWthresh(cxcsoft(justdet),sorh2,thr);%软阈值去噪  
cxchalf(justdet) = myWthresh(cxchalf(justdet),sorh3,thr);%软阈值去噪  

%小波重建  
xchard = waverec2(cxchard,l,w);  
xcsoft = waverec2(cxcsoft,l,w);  
xchalf = waverec2(cxchalf,l,w);   

figure(2);   
imshow(uint8(xchard(1:iCount, :)));title('硬阈值去噪图像');    
    
figure(3);   
imshow(uint8(xcsoft(1:iCount, :)));title('软阈值去噪图像');   

figure(4);   
imshow(uint8(xchalf(1:iCount, :)));title('半阈值去噪图像'); 




三,C加加实现

说明:如同一维的阈值去噪一样,在执行自己编写的wavedec2函数时必须先初始化,初始化的目的是为了获取信号的长度,选择的是什么小波,以及分解的等级等信息,然后计算出未来的各种信息,比如每个等级的系数的size,其中共有变量m_msgCL2D记录了这些信息。二维小波分解的初始化函数如下:

//初始化二维图像的分解信息,保存未来需要的信息
bool  CWavelet::InitDecInfo2D(
	const int height,//预分解的图像的高度
	const int width,//预分解的图像的宽度
	const int Scale,//分解尺度
	const int dbn//db滤波器编号,有默认值
	)
{
	if (dbn != 3)
		SetFilter(dbn);

	if (height < m_dbFilter.filterLen - 1 || width < m_dbFilter.filterLen - 1)
	{
		cerr << "错误信息:滤波器长度大于信号的高度或者宽度!" << endl;
		return false;
	}

	int srcHeight = height;
	int srcWidth = width;
	m_msgCL2D.dbn = dbn;
	m_msgCL2D.Scale = Scale;
	m_msgCL2D.msgHeight.resize(Scale + 2);
	m_msgCL2D.msgWidth.resize(Scale + 2);
	//源图像的尺寸
	m_msgCL2D.msgHeight[0] = height;
	m_msgCL2D.msgWidth[0] = width;

	//每一尺度上的尺寸
	for (int i = 1; i <= Scale; i++)//注意:每个尺度的四个分量的长宽是一样的
	{
		int exHeight = (srcHeight + m_dbFilter.filterLen - 1) / 2;//对称拓延后系数的长度的一半
		srcHeight = exHeight;
		m_msgCL2D.msgHeight[i] = srcHeight;

		int exWidth = (srcWidth + m_dbFilter.filterLen - 1) / 2;//对称拓延后系数的长度一半
		srcWidth = exWidth;
		m_msgCL2D.msgWidth[i] = srcWidth;
	}
	m_msgCL2D.msgHeight[Scale + 1] = srcHeight;
	m_msgCL2D.msgWidth[Scale + 1] = srcWidth;

	//计算总的数据个数
	int tmpAllSize = 0;
	int curPartSize = 0;
	int prePartSize = 0;
	for (int i = 1; i <= Scale; i++)
	{
		curPartSize = m_msgCL2D.msgHeight[i] * m_msgCL2D.msgWidth[i];
		tmpAllSize += curPartSize * 4 - prePartSize;
		prePartSize = curPartSize;
	}
	m_msgCL2D.allSize = tmpAllSize;
	m_bInitFlag2D = true;
	return true;
}


1,核心函数的实现

1)二维信号的单次分解

说明:本函数建立在一维的小波分解函数基础上(DWT)

// 二维数据的小波分解
void  CWavelet::DWT2(
	double *pSrcImage,//源图像数据(存储成一维数据,行优先存储)
	int height,//图像的高
	int width,//图像的宽
	double *pDstCeof//分解出来的图像系数
	)
{

	if (!m_bInitFlag2D)
	{
		cerr << "错误信息:未初始化,无法对信号进行分解!" << endl;
		return;
	}

	if (pSrcImage == NULL || pDstCeof == NULL)
	{
		cerr << "错误信息:dwt2数据无内存" << endl;
		Sleep(3000);
		exit(1);
	}
	
	int exwidth = (width + m_dbFilter.filterLen - 1) / 2 * 2;//pImagCeof的宽度
	int exheight = (height + m_dbFilter.filterLen - 1) / 2 * 2;//pImagCeof的高度

	double *tempImage = new double[exwidth*height];

	
	// 对每一行进行行变换
	double *tempAhang = new double[width]; 
	double *tempExhang = new double[exwidth]; // 临时存放每一行的处理数据
	for (int i = 0; i < height; i++)
	{
		for (int j = 0; j < width; j++)
			tempAhang[j] = pSrcImage[i*width + j];//提取每一行的数据

		DWT(tempAhang, width, tempExhang);

		for (int j = 0; j < exwidth; j++)
			tempImage[i*exwidth + j] = tempExhang[j];
	}

	// 对每一列进行列变换
	double *tempAlie = new double[height]; // 临时存放每一列的转置数据
	double *tempexlie = new double[exheight]; // 临时存放每一列的处理数据
	for (int i = 0; i < exwidth; i++)
	{
		// 列转置
		for (int j = 0; j < height; j++)
			tempAlie[j] = tempImage[j*exwidth + i];//提取每一列数据

		//执行变换
		DWT(tempAlie, height, tempexlie);

		// 反转置
		for (int j = 0; j < exheight; j++)
			pDstCeof[j*exwidth + i] = tempexlie[j];
	}

	AdjustData(pDstCeof, exheight, exwidth);//调整数据

	delete[] tempAlie;
	tempAlie = NULL;
	delete[] tempexlie;
	tempexlie = NULL;

	delete[] tempAhang;
	tempAhang = NULL;
	delete[] tempExhang;
	tempExhang = NULL;

	delete[] tempImage;
	tempImage = NULL;
	
}



2)二维信号的单次重构

说明:

//二维小波反变换
void  CWavelet::IDWT2(
	double *pSrcCeof, //二维源图像系数数据
	int dstHeight,//重构出来后数据的高度
	int dstWidth,//重构出来后数据的宽度
	double *pDstImage//重构出来的图像
	)
{
	int srcHeight = (dstHeight + m_dbFilter.filterLen - 1) / 2 * 2;
	int srcWidth = (dstWidth + m_dbFilter.filterLen - 1) / 2 * 2;//pSrcCeof的高度
	IAdjustData(pSrcCeof, srcHeight, srcWidth);//调整成LL,HL,LH,HH

	double *tempAline = new double[srcHeight]; // 临时存放每一列的数据
	double *tempdstline = new double[dstHeight]; // 临时存放每一列的重构结果

	double *pTmpImage = new double[srcWidth*dstHeight];
	// 列重构
	for (int i = 0; i < srcWidth; i++)//每一列
	{
		// 列转置
		for (int j = 0; j<srcHeight; j++)
			tempAline[j] = pSrcCeof[j*srcWidth + i];//提取每一列

		IDWT(tempAline, dstHeight, tempdstline);

		// 反转置
		for (int j = 0; j < dstHeight; j++)
			pTmpImage[j*srcWidth + i] = tempdstline[j];
	}


	// 对每一行进行行变换
	double *tempAhang = new double[srcWidth];
	double *tempdsthang = new double[dstWidth]; // 临时存放每一行的处理数据
	for (int i = 0; i < dstHeight; i++)
	{
		for (int j = 0; j < srcWidth; j++)
			tempAhang[j] = pTmpImage[i*srcWidth + j];//提取每一行的数据

		IDWT(tempAhang, dstWidth, tempdsthang);

		for (int j = 0; j < dstWidth; j++)
			pDstImage[i*dstWidth + j] = tempdsthang[j];
	}


	delete[] tempAline;
	tempAline = NULL;
	delete[] tempdstline;
	tempdstline = NULL;

	delete[] tempAhang;
	tempAhang = NULL;
	delete[] tempdsthang;
	tempdsthang = NULL;

	delete[] pTmpImage;
	pTmpImage = NULL;
}


3)二维信号的多级分解

说明:对于每一级分解都将调用单次二维分解函数来实现,所以本函数是建立在函数IDW2基础上

// 二维小波多级分解,需要先初始化获取未来数据信息
bool CWavelet::WaveDec2(
	double *pSrcData,//源图像数据,存储为一维信号
	double *pDstCeof//分解后的系数,它的大小必须是m_msgCL2D.allSize
	)
{
	if (!m_bInitFlag2D)
	{
		cerr << "错误信息:未初始化,无法对图像进行分解!" << endl;
		return false;
	}
	if (pSrcData == NULL || pDstCeof == NULL)//错误:无内存
		return false;

	int height = m_msgCL2D.msgHeight[0];
	int width = m_msgCL2D.msgWidth[0];
	int scale = m_msgCL2D.Scale;

	// 临时变量,图像数据
	double *tempImage = new double[height*width];
	int maxCoefSize =4 * m_msgCL2D.msgHeight[1] * m_msgCL2D.msgWidth[1];
	double *tempDst = new double[maxCoefSize];

	for (int i = 0; i < height*width; i++)
		tempImage[i] = pSrcData[i];

	int gap = m_msgCL2D.allSize - maxCoefSize;
	for (int i = 1; i <= scale; i++)
	{
		DWT2(tempImage, height, width, tempDst);

		// 低频子图像的高和宽
		height = m_msgCL2D.msgHeight[i];
		width = m_msgCL2D.msgWidth[i];
		
		for (int j = 0; j < height*width; j++)
			tempImage[j] = tempDst[j];//提取低频系数(近似系数)
		//
		for (int j = 0, k = gap; j < 4 * height*width; j++, k++)
			pDstCeof[k] = tempDst[j];//所有系数
		gap -= 4 * m_msgCL2D.msgWidth[i + 1] * m_msgCL2D.msgHeight[i + 1] - height*width;
	}
	delete[] tempDst;
	tempDst = NULL;
	delete[] tempImage;
	tempImage = NULL;
	return true;
}


4)多级分解系数的重构

// 根据多级分解系数重构出二维信号,必须先初始化获取分解信息
bool CWavelet::WaveRec2(
	double *pSrcCoef,//多级分解出的源系数
	double *pDstData//重构出来的信号
	)
{
	if (!m_bInitFlag2D)
	{
		cerr << "错误信息:未初始化,无法对信号进行分解!" << endl;
		return false;
	}
	if (pSrcCoef == NULL || pDstData == NULL)//错误:无内存
		return false;

	int height = m_msgCL2D.msgHeight[0];
	int width = m_msgCL2D.msgWidth[0];
	int decLevel = m_msgCL2D.Scale;

	int maxCeofSize = 4 * m_msgCL2D.msgHeight[1] * m_msgCL2D.msgWidth[1];

	double *pTmpImage = new double[maxCeofSize];

	int minCeofSize = 4 * m_msgCL2D.msgHeight[decLevel] * m_msgCL2D.msgWidth[decLevel];
	for (int i = 0; i < minCeofSize; i++)
		pTmpImage[i] = pSrcCoef[i];

	int gap = minCeofSize;
	for (int i = decLevel; i >= 1; i--)
	{
		int nextheight = m_msgCL2D.msgHeight[i - 1];//重构出来的高度
		int nextwidth = m_msgCL2D.msgWidth[i - 1];//重构出来的宽度
		IDWT2(pTmpImage, nextheight, nextwidth, pDstData);

		if (i > 1)//i==1已经重构出来了,不再需要提取系数
		{
			for (int j = 0; j < nextheight*nextwidth; j++)
				pTmpImage[j] = pDstData[j];
			for (int j = 0; j < 3 * nextheight*nextwidth; j++)
				pTmpImage[nextheight*nextwidth + j] = pSrcCoef[gap + j];
			gap += 3 * nextheight*nextwidth;
		}
	}
	
	delete[] pTmpImage;
	pTmpImage = NULL;

	return true;
}


2,函数正确性验证

1),二维单次分解与重构测试



2)二维多级分解与重构测试

说明:对二维数据进行了5层分解,选取的是小波族db3



3,阈值去噪结果:

说明:本测试只是模拟测试,对图像的处理也是一样的(完全一致)

硬阈值去噪结果




软阈值去噪结果





实际的VC++图像处理结果为:

源噪声图像为:



注意以下是采用:db6,3层分解,软阈值去噪,阈值是在前文提及的阈值基础上缩小2.5倍得到的效果:



注意以下是采用:db6,3层分解,硬阈值去噪,阈值是在前文提及的阈值基础上缩小2.5倍得到的效果:




附带上述matlab验证程序

clc;  
clear all;  
close all;  
% %%%%%%%%%%%%%%%%%%%%%%%%通过matlab的函数来实现阈值去噪%%%%%%%%%%%%%%%%%%%%%%%%%% %  
X=[ 10, 12, 30, 4, 5, 61, 2, 3;
    41, 5, 6, 27, 3, 4, 15, 6;
    72, 8, 41, 5, 6, 7, 8, 9;
    5, 64, 7, 8, 9, 14, 6, 27;
    8, 9, 40, 31,10, 12, 30, 4;
    50, 61, 2, 3, 41, 5, 6, 27];

X=double(X);  

% 获取输入参数  
wname    = 'db3';%小波类型  
n    = 3;%分解层数  
sorh1 = 'h';%硬阈值  
sorh2 = 's';%软阈值  

% 对图像进行小波分解  
[c,l] = wavedec2(X,n,wname);

%求取阈值
N = numel(X);
[chd1,cvd1,cdd1] = detcoef2('all',c,l,1); 
cvd1=cvd1(:)';
sigma = median(abs(cvd1))/0.6745;%提取细节系数求中值并除以0.6745
thr = sigma*sqrt(2*log(N)); 

% 对小波系数全局阈值处理  
cxch = c;% 保留近似系数  
cxcs = c;% 保留近似系数  
justdet = prod(l(1,:))+1:length(c);%截取细节系数(不处理近似系数)  

% 阈值处理细节系数  
cxch(justdet) = wthresh(cxch(justdet),sorh1,thr);%硬阈值去噪  
cxcs(justdet) = wthresh(cxcs(justdet),sorh2,thr);%软阈值去噪  

%小波重建  
xch = waverec2(cxch,l,wname);  
xcs = waverec2(cxcs,l,wname);  


鉴于不少人要代码,出于不误导人,给一些我自己的建议:
1,对程序要有自己的思考,要有自己的验证和学习过程,写时当时我还是个C++新手!点击下载完整程序
2,二维小波变换不建议用,因为速度太慢了(自己徒手写的),一维可以用。
3,推荐你另外的小波库网站,http://wavelet2d.sourceforge.net/

注:本博文为EbowTang原创,后续可能继续更新本文。如果转载,请务必复制本条信息!

原文地址:http://blog.csdn.net/ebowtang/article/details/40481539

原作者博客:http://blog.csdn.net/ebowtang


参考资源:

【1】《数字图像处理》(冈萨雷斯matlab第二版)
【2】 Donoho D L. De-noising by soft-thresholding. IEEE Trans- actions on Information Theory, 1995, 41(3): 613−627
【3】 Donoho D L, Johnstone I M. Adapting to unknown smooth- ness via wavelet shrinkage. Journal of the American Statistic
Association, 1995, 90(432): 1200−1224
【4】 Donoho D L, Johnstone I M, Kerkyacharian G, Picard D. Wavelet shrinkage: asymptopia? Journal of Royal Statisti-
cal Society Series B, 1995, 57(2): 301−369
【5】杨恢先,王绪四,改进阈值与尺度间相关的小波红外图像去噪
【6】《小波分析及其应用》,孙延奎著
【7】杨建国.小波分析及其工程应用[M].北京:机械工业出版社.2005
【8】毛艳辉.小波去噪在语音识别预处理中的应用.上海交通大学硕士学位论文.2010
【9】matlab各种函数说明,及其内部函数实现
【10】http://www.bearcave.com/misl/misl_tech/wavelets/haar.html

2018-07-24 17:38:38 qq_29441995 阅读数 1178

在图像处理中,经常处理需要过滤处理一些噪声,opencv自带算法不一定能满足我们工程的要求,自己写了一个调暗处理的简单函数,记录下来:

//取ksize大小矩形区域内的最小值替代整个区域的像素值
void Ave2(Mat mat, Mat&ave,int ksize)
{
    Size size(ksize, ksize);
    ave = mat.clone();
    for (int i = 0; i < mat.rows; i++)
    {
        for (int j = 0; j < mat.cols; j++)
        {
            int min = 256;
            for (int m = 0; m < size.width; m++)
            {
                for (int n = 0; n < size.height; n++)
                {
                    int x = m - size.width / 2 + j;
                    int y = n - size.height / 2 + i;
                    if (x < 0 || x >= mat.cols || y < 0 || y >= mat.rows)
                        continue;
                    int data = mat.at<unsigned char>(y, x);
                    if (data < min)
                        min = data;
                }
            }
            ave.at<unsigned char>(i, j) = min;
        }
    }
}

//p点的size范围内如果非0值的数量小于threshold,则置为0
void myBlur(Mat mat,Mat& blur,Size size,int threshold)
{
    blur = mat.clone();
    for (int i = 0; i < mat.rows; i++)
    {
        for (int j = 0; j < mat.cols; j++)
        {
            int count = 0;
            for (int m = 0; m < size.width; m++)
            {
                for (int n = 0; n < size.height; n++)
                {
                    int x = m - size.width / 2 + j;
                    int y = n - size.height / 2 + i;
                    if (x < 0 || x >= mat.cols || y < 0 || y >= mat.rows)
                        continue;
                    int data = mat.at<unsigned char>(y, x);
                    if (data != 0)
                        count++;
                }
            }
            if (count <= threshold)
            {
                blur.at<unsigned char>(i, j)=0;
            }
        }
    }
}

2018-01-06 21:24:55 zhangquan2015 阅读数 42244

Github个人博客:https://joeyos.github.io

图像去噪常用方法

图像去噪处理方法可分为空间域法和变换域法两大类。

基于离散余弦变换的图像去噪

一般而言,我们认为图像的噪声在离散余弦变换结果中处在其高频部分,而高频部分的幅值一般很小,利用这一性质,就可以实现去噪。然而,同时会失去图像的部分细节。

%读取图像
X=imread('wangshi.jpg'); 
X=rgb2gray(X);
%读取图像尺寸
[m,n]=size(X); 
%给图像加噪
Xnoised=imnoise(X,'speckle',0.01); 
%输出加噪图像
subplot(121); 
imshow(Xnoised);
%DCT变换
Y=dct2(Xnoised); 
I=zeros(m,n);
%高频屏蔽
I(1:m/3,1:n/3)=1; 
Ydct=Y.*I;
%逆DCT变换
Y=uint8(idct2(Ydct)); 
%结果输出
subplot(122);
imshow(Y);

这里写图片描述

基于小波变换的图像去噪

小波去噪是小波变换较为成功的一类应用,其去噪的基本思路为:含噪图像-小波分解-分尺度去噪-小波逆变换-恢复图像。含噪信号经过预处理,然后利用小波变换把信号分解到各尺度中,在每一尺度下把属于噪声的小波系数去掉,保留并增强属于信号的小波系数,最后再经过小波逆变换恢复检测信号。比基于傅里叶变换的去噪方法好。

clear;                 
X=imread('life.jpg');            
X=rgb2gray(X);
subplot(221);          
imshow(X);             
title('原始图像');                  
% 生成含噪图像并图示
init=2055615866;       
randn('seed',init);      
X=double(X);
% 添加随机噪声
XX=X+8*randn(size(X));  
subplot(222);             
imshow(uint8(XX));              
title(' 含噪图像 ');       
%用小波函数coif2对图像XX进行2层
% 分解
[c,l]=wavedec2(XX,2,'coif2'); 
% 设置尺度向量
n=[1,2];                  
% 设置阈值向量 , 对高频小波系数进行阈值处理
p=[10.28,24.08]; 
nc=wthcoef2('h',c,l,n,p,'s');
% 图像的二维小波重构
X1=waverec2(nc,l,'coif2');   
subplot(223);              
imshow(uint8(X1));                
%colormap(map);            
title(' 第一次消噪后的图像 '); 
%再次对高频小波系数进行阈值处理
mc=wthcoef2('v',nc,l,n,p,'s');
% 图像的二维小波重构
X2=waverec2(mc,l,'coif2');  
subplot(224);             
imshow(uint8(X2));               
title(' 第二次消噪后的图像 ');   

这里写图片描述

图像去噪处理

阅读数 3152

没有更多推荐了,返回首页