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

《基于偏微分方程的图像处理》 《基于偏微分方程的图像处理》 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分量'); (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); 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．​ 显示一幅灰度图像及其加高斯白噪声（ ）和椒盐噪声（噪声密度 ）后的图像及直方图。 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); 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') 算术平均滤波： 加权的线性光滑滤波: 热传导方程滤波: