心脏图像去噪

这篇博文围绕如下问题展开

问题：

分析心脏B超图的噪声类型并对其进行去噪处理。

心脏B超图用如下软件打开

对于这个问题，首先找了一篇论文，其中提出可以通过离散小波变换识别噪声的类型。论文标题如下：

论文提出，不同类型噪声在图像小波子带HH系数的分布大不相同。其系数直方图有以下特点:

1)对称于零点并在该点处 达到最大值;

2)高斯噪声的不显著小波系数较多, 而椒盐噪声的不显著小波系数较少。

根据以上特点, 可以用系数能量对噪声类型进行识别。

HH系数矩阵为{D(l,k) (l,k=1,2,…,N/2)}，按系数幅值（绝对值）的大小可将不显著的系数提出：

    其中T为确定系数是否显著的阈值，进一步定义不显著系数的能量与系数总能量与系数总能量的比为：

   (2)

    式中num为总的HH系数个数，而为满足(1)式的系数个数。

    在此基础上就可以判断图像的噪声类型。先计算和并记录显著个数和总系数个数num。

接下来就可以根据公式(2)计算ER的值，论文提出，当ER>0.5时，为高斯噪声；当ER<0.5时，为椒盐噪声。

题目提供的心脏图片的ER值为0.7973，因此噪声类型为高斯噪声，可以用均值模板去除心脏图片的噪声。去除噪声后的心脏图片如下图所示：

如左下角所示的图片，虽然去除了一部分噪声，但是效果并不好。除此之外，还分析心脏图片的直方图分布，如下图所示：

    从结果来看，类似于卡方分布或者瑞利分布，这种分布可以对其做一个小波分解，之后对小波系数做一个软阈值去噪。但是去噪图片会特别模糊，类似于频率域低通滤波。

    利用频率域去除噪声，虽然可以去除部分噪声，但是图像会因为失去边缘、细节、纹理等高频信息而变得模糊，因此使用gamma矫正，结果如下：

    利用伽马矫正，得到的去噪图像不仅保留了边缘、细节、纹理等高频信息，而且还去除了噪声，效果还可以。

matlab代码

%%%%%%%二、均值滤波、锐化滤波与伽马矫正
nii = load_nii( 'D:\课程\数字图像处理\课程设计\image\echocardiogram.nii' );  % 装载.nii数据
img = nii.img;  %文件包含img和head，img是图像数据
image = img(:,:,1);
figure;
imshow(image,[]);title('原图');

% 获取图像数据
img_data = double(image);
% 绘制直方图
figure;
histogram(img_data, 'BinWidth', 1, 'FaceColor', 'b');
title('心脏图片直方图');
xlabel('像素值');
ylabel('像素数量');

% 将图像数据转换为灰度图像
% 读取.nii文件
nii_data = niftiread('D:\课程\数字图像处理\课程设计\image\echocardiogram.nii');
gray_image = mean(nii_data, 4);

a=zeros(size(image));
[row,col]=size(image);
for i=1:row
    for j=1:col
        if(image(i,j)>5)
            a(i,j)=1;
        else
            a(i,j)=0;
        end
    end
end

dd2=xiaobo(image);
ERR = calculateER(image, 20);
disp(['Energy Ratio (ERR): ' num2str(ERR)]);

%用均值滤波去除高斯噪声
avg_filter_image=imfilter(image,fspecial('average',3),'replicate');
avg_filter_image1=imfilter(avg_filter_image,fspecial('average',3),'replicate');
avg_filter_image2=imfilter(avg_filter_image1,fspecial('average',3),'replicate');
avg_filter_image2=uint8(avg_filter_image2);

% %定义拉普拉斯锐化滤波模板
% H1=[0,1,0;1,-4,1;0,1,0];
% H2=[1,1,1;1,-8,1;1,1,1];
% %用H1对原图进行拉普拉斯锐化滤波
% sharpImage1=imfilter(double(avg_filter_image1),H1);

% H1=[-1,-2,-1;0,0,0;1,2,1];H2=[-1,0,1;-2,0,2;-1,0,1];
% H1=[-1,-1,-1;0,0,0;1,1,1];H2=[-1,0,1;-1,0,1;-1,0,1];
H1=[1,0;0,-1];H2=[0,1;-1,0];
R1=imfilter(double(avg_filter_image1),H1);
R2=imfilter(double(avg_filter_image1),H2);
edgeImage=abs(R1)+abs(R2);

minx=min(edgeImage(:));
maxx=max(edgeImage(:));
edgeImage=(edgeImage-minx)/(maxx-minx);
edgeImage=edgeImage.*255;
edgeImage=uint8(edgeImage);

sharpImage=uint8(avg_filter_image1)+edgeImage;

%显示均值滤波、锐化滤波结果
figure;
subplot(2,2,1);imshow(image,[]);title('原图');
subplot(2,2,2);imshow(edgeImage);title('边缘细节图像');
subplot(2,2,3);imshow(avg_filter_image2);title('处理后的图像');
subplot(2,2,4);imshow(sharpImage);title('边缘细节图像+处理后的图像');

% 显示矫正后的图像
% 伽马矫正参数
gamma_value = 1.5; % 可调整的伽马值，通常在1和2之间
% 应用伽马矫正
gamma_corrected_image = gray_image.^gamma_value;
%显示伽马矫正结果
figure;
imshow(gamma_corrected_image,[]);title('Gamma Corrected Image');

function ER = calculateER(image, threshold)
    % 读取图像
    img = image;
    
    % 小波分解
    [c, s] = wavedec2(img, 2, 'sym5');
    
    % HH 子带提取
    HH_coefficients = wrcoef2('d', c, s, 'sym5', 2);
    
    % 计算显著系数个数和总系数个数
    numT = sum(abs(HH_coefficients(:)) < threshold);
    num = numel(HH_coefficients);
    % 计算 ER
    numerator = sum(HH_coefficients(abs(HH_coefficients) < threshold).^2);
    denominator = sum(HH_coefficients.^2, 'all');

    disp(numerator);
    disp(denominator);
    if numT == 0 || denominator == 0
        ER = 0;  % 可以根据实际情况调整
    else
        ER = (1/numT) * numerator / ((1/num) * denominator);
    end
end

function [dd2]=xiaobo(image)
% X包含图像 % 图像分解尺度为2，采用sym5小波
[c,s] = wavedec2(image,2,'sym5'); 
%%高频分量
%HH
dd2 = wrcoef2('d',c,s,'sym5',2);
end

%自编均值滤波
function result=myMeanFilter(inputImage,template)
    [row,col]=size(inputImage);
    [t_row,t_col]=size(template);
    result=zeros(row,col);
    padding=floor(t_row/2);
    new_inputImage=zeros(row+padding*2,col+padding*2);
    new_inputImage(1+padding:row+padding*2-padding,1+padding:col+padding*2-padding)=inputImage(:,:);
    for i=1+padding:row+padding*2-padding
        for j=1+padding:col+padding*2-padding
            window=new_inputImage(i-padding:i+padding,j-padding:j+padding);
            result(i-padding,j-padding)=sum(sum(window.*template));
        end
    end
end

绘制心脏B超图直方图matlab代码

nii = load_nii( 'D:\课程\数字图像处理\课程设计\image\echocardiogram.nii.gz' );  % 装载.nii数据
img = nii.img;  %文件包含img和head，img是图像数据
Image = img(:,:,1);
figure;
imshow(Image,[]);title('原图');
% 获取图像数据
img_data = double(Image);
disp(size(img_data));

% 绘制直方图
figure;
histogram(img_data, 'BinWidth', 1, 'FaceColor', 'b');
title('NIfTI文件直方图');
xlabel('像素值');
ylabel('像素数量');