基于曲波变换的医学图像噪声抑制

对传统去噪算法、小波去噪算法和曲波去噪算法的理论进行了介绍，采用Matlab编程对两种算法进行了测试。对原图加不同类型和不同强度的噪声，采用小波去噪和曲波去噪算法对这些图片进行去噪仿真实验。计算不同噪声图片处理后的信噪比以及人眼视觉直观感受两个方面来评价去噪的效果。通过观察得到的去噪图片以及信噪比数据，小波去噪和曲波去噪算法对高斯噪声都能起到良好的去噪效果。而对于椒盐噪声，则会导致图片的失真，但曲波去噪算法处理椒盐噪声后得到的图片能够更加符合人眼的视觉感受。

function C = fdct_wrapping（x is_real finest nbscales nbangles_coarse）

% fdct_wrapping.m - Fast Discrete Curvelet Transform via wedge wrapping - Version 1.0
%
% Inputs
%   x           M-by-N matrix
%
% Optional Inputs
%   is_real     Type of the transform
%                   0: complex-valued curvelets
%                   1: real-valued curvelets
%               [default set to 0]
%   finest      Chooses one of two possibilities for the coefficients at the
%               finest level:
%                   1: curvelets
%                   2: wavelets
%               [default set to 2]
%   nbscales    number of scales including the coarsest wavelet level
%               [default set to ceil（log2（min（MN）） - 3）]
%   nbangles_coarse
%               number of angles at the 2nd coarsest level minimum 8
%               must be a multiple of 4. [default set to 16]
%
%
%
% See also ifdct_wrapping.m fdct_wrapping_param.m
%
% By Laurent Demanet 2004

X = fftshift（fft2（ifftshift（x）））/sqrt（prod（size（x）））;
[N1N2] = size（X）;
if nargin < 2 is_real = 0; end;
if nargin < 3 finest = 2; end;
if nargin < 4 nbscales = ceil（log2（min（N1N2）） - 3）; end;
if nargin < 5 nbangles_coarse = 16; end;

% Initialization: data structure
nbangles = [1 nbangles_coarse .* 2.^（ceil（（nbscales-（nbscales:-1:2））/2））];
if finest == 2 nbangles（nbscales） = 1; end;
C = cell（1nbscales）;
for j = 1:nbscales
    C{j} = cell（1nbangles（j））;
end;

% Loop: pyramidal scale decomposition
M1 = N1/3;
M2 = N2/3;
if finest == 1

    % Initialization: smooth periodic extension of high frequencies
    bigN1 = 2*floor（2*M1）+1;
    bigN2 = 2*floor（2*M2）+1;
    equiv_index_1 = 1+mod（floor（N1/2）-floor（2*M1）+（1:bigN1）-1N1）;
    equiv_index_2 = 1+mod（floor（N2/2）-floor（2*M2）+（1:bigN2）-1N2）;
    X = X（equiv_index_1equiv_index_2）;
        % Invariant: equiv_index_1（floor（2*M1）+1） == （N1 + 2 - mod（N12））/2
        % is the center in frequency. Same for M2 N2.
    window_length_1 = floor（2*M1） - floor（M1） - 1 - （mod（N13）==0）;
    window_length_2 = floor（2*M2） - floor（M2） - 1 - （mod（N23）==0）;
        % Invariant: floor（M1） + floor（2*M1） == N1 - （mod（M13）~=0）
        % Same for M2 N2.
    coord_1 = 0:（1/window_length_1）:1;
    coord_2 = 0:（1/window_length_2）:1;
    [wl_1wr_1] = fdct_wrapping_window（coord_1）;
    [wl_2wr_2] = fdct_wrapping_window（coord_2）;
    lowpass_1 = [wl_1 ones（12*floor（M1）+1） wr_1];
    if mod（N13）==0 lowpass_1 = [0 lowpass_1 0]; end;
    lowpass_2 = [wl_2 ones（12*floor（M2）+1） wr_2];
    if mod（N23）==0 lowpass_2 = [0 lowpass_2 0]; end;
    lowpass = lowpass_1‘*lowpass_2;
    Xlow = X .* lowpass;

    scales = nbscales:-1:2;

else
    
    M1 = M1/2;
    M2 = M2/2;
    window_length_1 = floor（2*M1） - floor（M1） - 1;
    window_length_2 = floor（2*M2） - floor（M2） - 1;
    coord_1 = 0:（1/window_length_1）:1;
    coord_2 = 0:（1/window_length_2）:1;
    [wl_1wr_1] = fdct_wrapping_windo

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     目录           0  2018-10-10 15:39  源码\
     文件      325120  2018-10-13 14:59  源码\小波去噪和曲波去噪.doc
     目录           0  2017-03-29 16:05  源码\曲波程序及结果\
     文件       15003  2018-10-13 14:57  源码\曲波程序及结果\fdct_wrapping.m
     文件         751  2014-04-27 11:56  源码\曲波程序及结果\fdct_wrapping_window.m
     文件       16277  2012-03-13 03:21  源码\曲波程序及结果\ifdct_wrapping.m
     文件       32087  2012-03-13 03:21  源码\曲波程序及结果\Lena.jpg
     文件       39936  2017-02-22 12:36  源码\曲波程序及结果\Thumbs.db
     文件        5998  2018-10-13 14:57  源码\曲波程序及结果\wavedenoise.m
     目录           0  2018-10-10 15:39  演示\
     文件    42717693  2017-03-24 18:08  演示\10545.wmv
     文件          64  2018-10-10 16:34  说明.txt
     目录           0  2018-10-10 15:38  需求\
     文件          84  2017-03-06 16:28  需求\新建文本文档 （2）.txt