MATLAB图像平滑与滤波实验

一、实验题目

1.      拍摄的含有正弦波（ref1.bmp, obj1.bmp）和三角波条纹(T20r1.bmp, T20o1.bmp)投影图像

2.      对图像进行滤波去除噪声，噪声来源未知，滤波方法和滤波器可参考课程中的内容

3.      选择最优的滤波器和滤波器参数，要求滤波后的效果：能够基本恢复原来正弦波或三角波的形状

4.      要求提交实验报告作为平时成绩，占总成绩50%，包括源代码和滤波前后的截图，滤波器设计和使用的思路，评价，分析等。

例如：

图像中一个y截面滤波前三角波波形（左），滤波后三角波波形（右）

注：滤波三角波要防止三角波顶底部的特征点因为滤波后变得圆滑平缓

二、实验代码：

2.1 滤波器函数filt.m

%均值滤波
subplot(1,2,1);
img=imread('obj1.bmp');
imshow(img);
xlabel('(1) 原始图像');
A=fspecial('average',[7,7]);%3x3均值滤波
m1=imfilter(img,A);
subplot(1,2,2),imshow(m1,[]);
xlabel('(2) 1次7x7均值滤波');
imwrite(m1,'obj1_1_7x7avgfilt.bmp');

%中值滤波
subplot(1,2,1);
img=imread('obj1.bmp');
imshow(img);
xlabel('(1) 原始图像');
m1=medfilt2(img,[3,3]);
subplot(1,2,2),imshow(m1,[]);
xlabel('(2) 1次3x3中值滤波');
imwrite(m1,'obj1_1_3x3medfilt.bmp');

%维纳滤波
subplot(1,2,1);
img=imread('obj1.bmp');
imshow(img);
xlabel('(1) 原始图像');
m1=wiener2(img, [3 3]);%3x3自适应维纳滤波
subplot(1,2,2),imshow(m1,[]);
xlabel('(2) 1次3x3维纳滤波');
imwrite(m1,'obj1_1_3x3_wienerfilt.bmp');

%高斯滤波 
subplot(1,2,1);
img=imread('T20r1.bmp');
imshow(img);
xlabel('(1) 原始图像');
%k=input('请输入高斯滤波器的均值/n');  
%n=input('请输入高斯滤波器的方差/n');  
A1=fspecial('gaussian',3,5);      %生成高斯序列  
m1=filter2(A1,img)/255;             %用生成的高斯序列进行滤波
subplot(1,2,2),imshow(m1,[]);
xlabel('(2) 1次均值为3方差为0.1的高斯滤波');
imwrite(m1,'T20r1_3_50_gaussianfilt.bmp');

2.2 波形比较plotComp.m

%滤波前后对比
figure;
subplot(2,1,1);
h=imread('T20r1_3_50_gaussianfilt.bmp');
[m n]=size(h);
s=1;
arr=[];
    for i=1:m
        for j=1:n       
            if(j==50)
            arr(s)=h(i,j,1);
            s=s+1; 
           end
        end         
    end
plot(arr);

subplot(2,1,2);
h2=imread('obj1_1_3x3avgfilt.bmp');
[m n]=size(h2);
s=1;
arr=[];
    for i=1:m
        for j=1:n       
            if(j==50)
            arr(s)=h2(i,j,1);
            s=s+1; 
           end
        end         
    end
plot(arr);

2.3 滤波评价RMSE(均方根误差) PSNR（峰值信噪比）

% Demo to calculate PSNR of a gray scale image.
% http://en.wikipedia.org/wiki/PSNR
% Clean up.
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;

%------ GET DEMO IMAGES ----------------------------------------------------------
% Read in a standard MATLAB gray scale demo image.
refImage = imread('T20r1_3_50_gaussianfilt.bmp');

% Display the first image.
subplot(3, 2, 1);
imshow(refImage, []);
title('After filt ref Image', 'FontSize', fontSize);

% Cut ref Image 从上到下 从左到右
refImage_cut = refImage(20:1000,20:150);
[rows columns] = size(refImage_cut);
subplot(3, 2, 2);
imshow(refImage_cut, []);
title('Cut after filt ref Image', 'FontSize', fontSize);

% Get a second image by adding noise to the first image.
%noisyImage = imnoise(grayImage, 'gaussian', 0, 0.003);

% Read the second image
objImage = imread('T20o1_3_50_gaussianfilt.bmp');
% Display the second image.
subplot(3, 2, 3);
imshow(objImage, []);
title('After filt obj Image', 'FontSize', fontSize);

% Cut obj Image 从上到下 从左到右
objImage_cut = objImage(20:1000,20:150);
subplot(3, 2, 4);
imshow(objImage_cut, []);
title('Cut After filt obj Image', 'FontSize', fontSize);

%------ PSNR CALCULATION ----------------------------------------------------------
% Now we have our two images and we can calculate the PSNR.
% First, calculate the "square error" image.
% Make sure they're cast to floating point so that we can get negative differences.
% Otherwise two uint8's that should subtract to give a negative number
% would get clipped to zero and not be negative.
squaredErrorImage = (double(refImage_cut) - double(objImage_cut)) .^ 2;
% Display the squared error image.
subplot(3, 2, 5);
imshow(squaredErrorImage, []);
title('Squared Error Image', 'FontSize', fontSize);
% Sum the Squared Image and divide by the number of elements
% to get the Mean Squared Error.  It will be a scalar (a single number).
mse = sum(sum(squaredErrorImage)) / (rows * columns);
rmse = sqrt(mse);
% Calculate PSNR (Peak Signal to Noise Ratio) from the MSE according to the formula.
PSNR = 10 * log10( 256^2 / mse);
% Alert user of the answer.
message = sprintf('The mean square error is %.2f.\nThe PSNR = %.2f', rmse, PSNR);
msgbox(message);