基于偏微分方程的图像处理

《基于偏微分方程的图像处理》
1．图像的基本操作
(1)把一幅彩色图像分解为R、G、B三副单色图像；
clear;
image_I=imread('rgbtest2.bmp');
subplot(2,2,1);
imshow(image_I);
matrix_R(:,:,1)=image_I(:,:,1);
matrix_R(:,:,2)=0;
matrix_R(:,:,3)=0;
subplot(2,2,2);
imshow(matrix_R);
title('R分量');
matrix_G(:,:,2)=image_I(:,:,2);
matrix_G(:,:,1)=0;
matrix_G(:,:,3)=0;
subplot(2,2,3);
imshow(matrix_G);
title('G分量');
matrix_B(:,:,3)=image_I(:,:,3);
matrix_B(:,:,1)=0;
matrix_B(:,:,2)=0;
subplot(2,2,4);
imshow(matrix_B);
title('B分量');
R分量
G分量B分量
(2)把一幅灰度图像分别沿x轴和y轴做反射，扩展为四倍大小；
clear;
Image=imread('graytest2.bmp');
imshow(Image);
[m,n]=size(Image);
image11=Image;
for i=1:m
for j=1:n
image12(i,j)=image11(i,n-j+1);
end
end
for j=1:n
for i=1:m
image21(i,j)=image11(m-i+1,j);
image22(i,j)=image12(m-i+1,j);
end
end
image1=[image11,image12];
image2=[image21,image22];
image=[image1;image2];
figure
imshow(image)
原图像：
扩展图像：
2．把一幅灰度图像的像素值都变换为原来的1/2使图像质量变
差，然后利用Matlab图像处理工具箱中的直方图均衡化函数histeq 对图像进行增强，输出原图像、质量变差图像和增强后的图像以及它们的直方图。

clear;
Image=imread('graytest2.bmp');
Image1=0.5*double(Image);
subplot(3,2,1);
imshow(Image);
subplot(3,2,2);
imhist(Image);
subplot(3,2,3);
imshow(uint8(Image1));
image1=uint8(Image1);
imwrite(image1,'graytest2_1.bmp')
subplot(3,2,4);
imhist(image1);
Image_hist=histeq(image1);
subplot(3,2,5);
imshow(Image_hist);
imwrite(Image_hist,'graytest2_2.bmp')
subplot(3,2,6);
imhist(Image_hist);
0100
200
5000
100
200
0500
10000
100
200
1000
3． 图像的水平集、水平线和等高线的显示和处理
(1) 显示一幅灰度图像在某一水平的水平集、水平线和等高线；
function level(image_name,N) %求一图像的水平集、线
%输入参数：image_name －－－图像文件名
% ：N －－－－－－－－阀值（0～255） if nargin<1
disp('请输入文件名!!!'); return; elseif nargin<2 N=128; end
image_input=imread(image_name); subplot(2,2,1);
imshow(image_input); title('原图像');
[size_m,size_n]=size(image_input); matrix_temp=zeros(size_m,size_n); % 求水平集 for row=1:size_m for col=1:size_n
if image_input(row,col) > N
matrix_temp(row,col)=1;
end
end
end
subplot(2,2,2);
imshow(matrix_temp,[]);
title('图像的水平集');
imwrite(matrix_temp,'level_setzhan.bmp');
% 图像矩阵扩展赋值便于处理边界
matrix_ex=zeros(size_m+2,size_n+2);
for row=1:size_m
for col=1:size_n
matrix_ex(row+1,col+1)=matrix_temp(row,col);
end
end
% 四邻域反填充得水平线
matrix_new=matrix_temp;
for row=2:size_m+1
for col=2:size_n+1
if matrix_ex(row+1,col)==0 & matrix_ex(row-1,col)==0 & matrix_ex(row,col+1)==0 & matrix_ex(row,col-1)==0
matrix_new(row-1,col-1)=1;
end
end
end
subplot(2,2,3);
imshow(matrix_new,[]);
title('图像的水平线');
%imwrite(matrix_new,'level_line.bmp');
% 求图像的等高线
contour=zeros(size_m,size_n);
for row=1:size_m
for col=1:size_n
if image_input(row,col)==N
contour(row,col)=1;
end
contour(row,col)=1-contour(row,col);
end
end
subplot(2,2,4);
imshow(contour,[]);
title('图像的等高线');
%imwrite(contour,'contour.bmp');
function level_line_set(image_name)
% 求一图像的阀值N取各值所对应的水平线的叠加
% 输入参数：image_name－－－图像文件名
image_input=imread(image_name); %不变量
[size_m,size_n]=size(image_input); %不变量
matrix_out=zeros(size_m,size_n); %循环中不需重新初始化的变量
% 输出矩阵初始化（白）
for row=1:size_m
for col=1:size_n
matrix_out(row,col)=1;
end
end
matrix_ex=zeros(size_m+2,size_n+2); % 扩展矩阵
subplot(1,2,1);
imshow(image_input);
title('原图像');
for N=1:32:255 % N为相应阀值
matrix_temp=zeros(size_m,size_n); % 重新初始化（黑）
for row=1:size_m
for col=1:size_n
if image_input(row,col) >=N
matrix_temp(row,col)=1; % 求水平集－－－－－
end
end
end
% 图像矩阵扩展赋值便于处理边界
for row=1:size_m
for col=1:size_n
matrix_ex(row+1,col+1)=matrix_temp(row,col);
end
end
% 四邻域（黑区域填白）反填充得水平线
matrix_new=matrix_temp;
for row=2:size_m+1
for col=2:size_n+1
if matrix_ex(row+1,col)==0 & matrix_ex(row-1,col)==0 & matrix_ex(row,col+1)==0 & matrix_ex(row,col-1)==0
matrix_new(row-1,col-1)=1; % 填白
end
end
end
matrix_out=min(matrix_out,matrix_new); %各级水平线叠加end
matrix_out(1,1)=0;
matrix_out(2,1)=1;
subplot(1,2,2);
imshow(matrix_out,[]);
title('图像部分水平集的叠加');
level('graytest2.bmp',100);
原图像图像的水平集
图像的水平线图像的等高线
(2)显示图像并显示它的一族水平线（要求阀值 为32, 64, 96, 128, 160, 192, 224）。

function level_line_set(image_name)
% 求一图像的阀值N取各值所对应的水平线的叠加
% 输入参数：image_name－－－图像文件名
image_input=imread(image_name); %不变量
[size_m,size_n]=size(image_input); %不变量
matrix_out=zeros(size_m,size_n); %循环中不需重新初始化的变量
% 输出矩阵初始化（白）
for row=1:size_m
for col=1:size_n
matrix_out(row,col)=1;
end
end
matrix_ex=zeros(size_m+2,size_n+2); % 扩展矩阵
subplot(1,2,1);
imshow(image_input);
title('原图像');
for N=1:32:255 % N为相应阀值（要求阀值 为32, 64, 96, 128, 160, 192, 224）。

matrix_temp=zeros(size_m,size_n); % 重新初始化（黑）
for row=1:size_m
for col=1:size_n
if image_input(row,col) >=N
matrix_temp(row,col)=1; % 求水平集－－－－－
end
end
end
% 图像矩阵扩展赋值便于处理边界
for row=1:size_m
for col=1:size_n
matrix_ex(row+1,col+1)=matrix_temp(row,col);
end
end
% 四邻域（黑区域填白）反填充得水平线
matrix_new=matrix_temp;
for row=2:size_m+1
for col=2:size_n+1
if matrix_ex(row+1,col)==0 & matrix_ex(row-1,col)==0 & matrix_ex(row,col+1)==0 & matrix_ex(row,col-1)==0
matrix_new(row-1,col-1)=1; % 填白
end
end
end
matrix_out=min(matrix_out,matrix_new); %各级水平线叠加end
matrix_out(1,1)=0;
matrix_out(2,1)=1;
subplot(1,2,2);
imshow(matrix_out,[]);
title('图像部分水平集的叠加');
原图
像图像部分水平集的叠加
原图
像图像部分水平集的叠加
原图像
图像部分水平集的叠加
原图像图像部分水平集的叠加
原图像图像部分水平集的叠加
原图像图像部分水平集的叠加
4．显示一幅灰度图像及其加高斯白噪声（01.0
μ）和椒盐
=σ
,02=
噪声（噪声密度02
d）后的图像及直方图。

=
.0
function noise(filename)
Image=imread(filename);
subplot(2,2,1);
imshow(Image,[]); % 原图像
title('原图像');
noise_g=imnoise(Image,'gaussian'); % 增加高斯白噪声
subplot(2,2,2);
imshow(noise_g,[]);
%imwrite(noise_g,'noise_g.bmp')
title('高斯白噪声');
noise_s=imnoise(Image,'salt & pepper');
subplot(2,2,3);
imshow(noise_s,[]);
title('黑白象素点噪声(椒盐)');
%imwrite(noise_s,'noise_s.bmp')
noise_p=imnoise(Image,'speckle');
subplot(2,2,4);
imshow(noise_p,[]);
%imwrite(noise_p,'noise_p.bmp')
title('乘积性噪声');
Image=imread('graytest2.bmp'); subplot(3,2,1);
imshow(Image,[]); % 原图像 title('原图像'); subplot(3,2,2); imhist(Image);
noise_g=imnoise(Image,'gaussian',0,0.01); % 增加高斯白噪声 subplot(3,2,3);
imshow(noise_g,[]); title('高斯白噪声');
imwrite(noise_g,'graytest2_g.bmp'); subplot(3,2,4); imhist(noise_g);
noise_s=imnoise(Image,'salt & pepper',0.02); subplot(3,2,5);
imshow(noise_s,[]);
title('黑白象素点噪声(椒盐)'); imwrite(noise_s,'graytest2_s.bmp') subplot(3,2,6); imhist(noise_s);
原图像
100
200
500高斯白噪声
100
200
500
黑白象素点噪声(椒盐)
100
200
500
5． 分别利用算术平均滤波、加权的线性光滑和热传导方程滤波
方法对4题中得到的两幅噪声图像进行去噪处理。

算术平均滤波：
function aver_filter( Image_name,N,out_filename )
%非线性平滑滤波器G1 平均
%输入参数文件名N为模板大小(正奇数)
%读取图像－－－－－－－－－－－－－
if nargin<3
disp('输入参数错误！！！');
return;
end
image=imread(Image_name);
if isrgb(image)
image=rgb2gray(image);
end
% 图像扩展－－－－－－－－－－－－－－
[size_m,size_n]=size(image);
image_ex=zeros(size_m+N-1,size_n+N-1);
% 赋值
for row=1:size_m
for col=1:size_n
image_ex(row+(N-1)/2,col+(N-1)/2)=image(row,col);
end
end
% 对每个象素－－－－－－－－－－－－－－－－－－
pattern_plate=zeros(N,N); %模板
image_new=zeros(size_m,size_n); %滤后图像矩阵
for row=1:size_m
for col=1:size_n
% 取模板的元素
for row_p=1:N
for col_p=1:N
pattern_plate(row_p,col_p)=image_ex(row+row_p-1,col+col_p-1);
end
end
% 求模板的平均值得滤后图像矩阵
image_new(row,col)=mean2(pattern_plate);
end
end
% 输出－－－－－－－－－－－－－－－
imshow(image_new,[]);
title('滤后图像');
imwrite( uint8(image_new),out_filename )
加权线性光滑：
function weight_filter( Image_name,flag,out_filename )
% 用各模板进行去噪和锐化
%输入参数Image_name 文件名
% flag 模板类型标志1~4
if nargin<3
disp('请按参数格式输入参数！！！');
return;
end
if flag~=1 & flag~=2 & flag~=3 & flag~=4
disp('参数flag输入错误！！！');
return;
end
N=3;
pattern_plate=zeros(N,N); %模板
% 各类模板
if flag==1
pattern_plate=[1/9 1/9 1/9; 1/9 1/9 1/9; 1/9 1/9 1/9];
elseif flag==2
pattern_plate=[1/16 1/8 1/16; 1/8 1/4 1/8; 1/16 1/8 1/16]; % 平滑滤波器
elseif flag==3
pattern_plate=[0 -1 0; -1 5 -1; 0 -1 0]; % 反锐化掩蔽滤波器elseif flag==4
A=1.7;
pattern_plate=[0 -1 0; -1 A+5 -1; 0 -1 0]; % 高提升滤波器end
% 读取图像－－－－－－－－－－
image_I=imread(Image_name);
if isrgb(image_I)
image_I=rgb2gray(image_I);
end
subplot(1,2,1);
imshow(image_I,[]);
title('原图像');
%图像扩展－－－－－－－－－－－－－－
[size_m,size_n]=size(image_I);
image_ex=zeros(size_m+N-1,size_n+N-1);
%赋值
for row=1:size_m
for col=1:size_n
image_ex(row+(N-1)/2,col+(N-1)/2)=image_I(row,col);
end
end
%对每个象素－－－－－－－－－－－－－－－－－－
temp_plate=zeros(N,N); %模板
image_new=zeros(size_m,size_n); %滤后图像矩阵
for row=1:size_m
for col=1:size_n
%取模板的元素
for row_p=1:N
for col_p=1:N
temp_plate(row_p,col_p)=image_ex(row+row_p-1,col+col_p-1);
end
end
image_new(row,col)=sum(sum(temp_plate.*pattern_plate));
end
end
%输出－－－－－－－－－－－－－－－
subplot(1,2,2);
imshow(uint8(image_new));
title('滤波后图像');
imwrite(uint8(image_new),out_filename);
热传导方程滤波：
function Heat_Equation( image_name,t,out_filename )
% 热传导方程[ Ut=U0+t*Delta(U0) ]
% image_name : 原图像文件名
% t : 迭代次数
% out_filename: 输出文件名
if nargin<3
disp('请按参数格式输入参数！！！');
return;
end
image_I=imread(image_name);
subplot(1,2,1);
imshow(image_I);
title('原图像');
mkdir('Heat_Eq'); % 新建文件夹以存放过程文件
cd('Heat_Eq');
[size_r,size_c]=size(image_I);
image_I=double(image_I);
if(t>0)
N=t/0.1;
else
N=(-t)/0.1% 每次叠代的步长为0.1
end
for index=1:N
image_ex=matrix_expand(image_I); % 边界拓展
diff_two_order=zeros(size_r,size_c);
for row=2:size_r+1
for col=2:size_c+1
diff_two_order(row-1,col-1)=-4*image_ex(row,col)+image_ex(row-1,col)+image_ex(row+1,col) +image_ex(row,col-1)+image_ex(row-1,col+1);
end
end
if(t>0)
image_I=image_I+0.1*diff_two_order;
else
image_I=image_I-0.1*diff_two_order;
end
imwrite( uint8(image_I),strcat('Heat_Eq_',num2str(index),'.bmp'));
end
cd('..')
subplot(1,2,2);
imshow(uint8(image_I));
if t>0
title('热传导方程光滑后图像');
else
title('Gabor锐化后的图像');
end
imwrite( uint8(image_I),out_filename );
% －－－－－－－－－－－－－－－－－－－－－－－－－end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%
function matrix_ex=matrix_expand(matrix)
% 对矩阵进行延拓
[size_r,size_c]=size(matrix);
matrix_ex=zeros(size_r+2,size_c+2);
% 中
for row=1:size_r
for col=1:size_c
matrix_ex(row+1,col+1)=matrix(row,col);
end
end
% 左，右
for row=1:size_r
matrix_ex(row+1,1)=matrix(row,1);
matrix_ex(row+1,size_c+2)=matrix(row,size_c);
end
% 上，下
for col=1:size_c
matrix_ex(1,col+1)=matrix(1,col);
matrix_ex(size_r+2,col+1)=matrix(size_r,col);
end
matrix_ex(1,1)=matrix(1,1);
matrix_ex(1,size_c+2)=matrix(1,size_c);
matrix_ex(size_r+2,1)=matrix(size_r,1);
matrix_ex(size_r+2,size_c+2)=matrix(size_r,size_c);
% －－－－－－－－－－－－－－－－－－－－－－－－－－运行命令：
aver_filter('graytest2_g.bmp',5,'averfilter_g.bmp')
figure
aver_filter('graytest2_s.bmp',5,'averfilter_s.bmp')
figure
weight_filter('graytest2_g.bmp',2,'weightfilter_g.bmp')
figure
weight_filter('graytest2_s.bmp',2,'weightfilter_s.bmp')
figure
Heat_Equation('graytest2_g.bmp',3,'Heat_Equation_g.bmp')
figure
Heat_Equation('graytest2_s.bmp',3,'Heat_Equation_s.bmp')
算术平均滤波：
滤后图像
滤后图像
加权的线性光滑滤波:
原图像滤波后图像
原图像滤波后图像
热传导方程滤波:
原图像热传导方程光滑后图像
如有帮助，欢迎下载。

21 原图
像
热传导方程光滑后图像。