北京航空航天大学数字图像处理实验报告.docx

成 绩 数字图像处理实验报告 院（系）名称电子信息工程学院 专业名称通信工程 学生学号 学生姓名2017年5月目录目录2实验一 图像变换实验3实验二 图像复原实验9实验三 图像分割处理实验14实验四 用Hough变换进行曲线的参数提取19实验五 手写数字识别26附录40实验一 图像变换实验1实验目的：学会对图像进行傅立叶等变换，在频谱上对图像进行分析，增进对图像频域上的感性认识，并用图像变换进行压缩。2实验内容：对Lena或cameraman图像进行傅立叶、离散余弦、哈达玛变换。在频域，对比他们的变换后系数矩阵的频谱情况，进一步，通过逆变换观察不同变换下的图像重建质量情况。3. 实验要求实验采用获取的图像，为灰度图像，该图像每象素由8比特表示。具体要求如下：（1）对图像进行傅立叶变换、获得变换后的系数矩阵；（2）将傅立叶变换后系数矩阵的频谱用图像输出，观察频谱；（3）通过设定门限，将系数矩阵中95%的（小值）系数置为0，对图像进行反变换，获得逆变换后图像；（4）观察逆变换后图像质量，并比较原始图像与逆变后的峰值信噪比（PSNR）。（5）对输入图像进行离散余弦、哈达玛变换，重复步骤1-5；（6）比较三种变换的频谱情况、以及逆变换后图像的质量（PSNR）。4实验结果图1.1傅里叶变换及逆变换图1.2离散余弦变换及逆变换图1.3哈达玛变换及逆变换观察结果可以发现，将系数矩阵中95%的（小值）系数置为0，对图像进行反变换之后获得逆变换的图像质量有所下降，信噪比提高。但是图像质量的下降在可以接受的范围内，并且过程中实现了对于图像的压缩。通过对于信噪比的比较可以发现：离散余弦变换的效果最好，傅里叶变换的效果最差。5实验代码：*傅里叶变换*function Fourier_Callback(hObject, eventdata, handles)% -显示原图像-figure=imread('lena512.bmp');subplot(1,3,1);imshow(figure);title('Lena 原图');% -傅里叶变换频谱图-fourier=fft2(figure);fourier=fftshift(fourier);R=real(fourier);I=imag(fourier);spectrum=abs(fourier);spectrum_norm =(spectrum-min(min(spectrum)/(max(max(spectrum)-min(min(spectrum)*255;subplot(1,3,2);imshow(spectrum_norm);title('傅里叶变换频谱图');% -傅里叶逆变换-height,width=size(spectrum);spectrum_adjust=sort(reshape(spectrum,1,height*width);threshold=spectrum_adjust(round(0.95*height*width);for i=1:height for j=1:width if spectrum(i,j)<=threshold R(i,j)=0; I(i,j)=0; end endendrefigure=ifft2(ifftshift(R+1i*I);subplot(1,3,3);imshow(uint8(refigure);title('逆变换图像');% -PSNR-PSNR = PSNR_cal(figure,refigure,8); set(handles.edit2,'string',PSNR);*离散余弦变换*function dct_Callback(hObject, eventdata, handles)% -显示原图像-figure=imread('lena512.bmp');subplot(1,3,1);imshow(figure);title('Lena 原图');% -离散余弦变换频谱图-DCT=dct2(figure);subplot(1,3,2);imshow(uint8(abs(DCT);title('离散余弦变换频谱图');% -离散余弦逆变换-Size=size(DCT);height=Size(1);width=Size(2);DCT_sort=sort(reshape(DCT,1,width*height);threhold=DCT_sort(round(width*height*0.95); DCT(abs(DCT)<threhold)=0;refigure=idct2(DCT);subplot(1,3,3);imshow(uint8(refigure);title('逆变换图像');% -PSNR-PSNR = PSNR_cal(figure,refigure,8); set(handles.edit2,'string',PSNR);*哈达玛变换*function Hadamard_Callback(hObject, eventdata, handles)% -显示原图像 -figure=imread('lena512.bmp');subplot(1,3,1);imshow(figure);title('Lena 原图');% -哈达玛变换频谱图-H=hadamard(512);figure=double(figure);DHT=(H*figure*H)./512;subplot(1,3,2);imshow(uint8(DHT);title('哈达玛变换频谱图');% -哈达玛逆变换-Size=size(figure);height=Size(1);width=Size(2);DHT_sort=sort(reshape(DHT,1,width*height);threhold=DHT_sort(round(width*height*0.95); DHT(abs(DHT)<threhold)=0;iH=H-1;refigure=(iH*DHT*iH)*512;subplot(1,3,3);imshow(uint8(refigure);title('逆变换图像');% -PSNR-PSNR = PSNR_cal(figure,refigure,8); set(handles.edit2,'string',PSNR);*PSNR计算函数*function PSNR = PSNR_cal( figure1,figure2,bit )img=double(figure2);imgn=double(figure1);Size=size(figure1);height=Size(1);width=Size(2);MAX=2bit-1; PMSE=sum(sum(img-imgn).2)/(height*width)/MAX2; PSNR=-10*log10(PMSE); end实验二 图像复原实验1实验目的：利用反向滤波和维纳滤波进行降质图像复原，比较不同参数选择对复原结果的影响。2实验内容：（1）利用反向滤波方法进行图像复原；（2）利用维纳滤波方法进行图像复原。3实验要求：（1）输入图像采用实验1所获取的图像，对输入图像采用运动降质模型，如下式所示与降值图像相关的参数是：；（2）对每一种方法通过计算复原出来的图像的峰值信噪比，进行最优参数的选择，包括反向滤波方法中进行复原的区域半径、维纳方法中的噪声对信号的频谱密度比值K；（3）将降质图像和利用最优参数恢复后的图像同时显示出来，以便比较。4实验结果：图2.1 反向滤波图2.2 维纳滤波反向滤波方法中进行复原的区域半径越大，复原效果越好，当r0=1时，能够完全复原降质图像；维纳滤波方法中的噪声对信号的频谱密度比值K值越小，复原效果越好，当K=0时，能够完全复原降质图像。5实验代码：*降质模型*lena=imread('lena512.bmp');subplot(2,3,1);imshow(lena);title('Lena 原图');x,y=size(lena);fftlena=fft2(lena);H1 = zeros(x,y);j=sqrt(-1); for uu=1:x; for vv=1:y; u=-x/2+uu-1; v=-y/2+vv-1; H1(uu,vv)=5/(pi*(u+v)*sin(pi*(u+v)*exp(-pi*j*(u+v); endendnnn=isnan(H1);for uu=1:x; for vv=1:y; if nnn(uu,vv)=1; u=-x/2+uu-1; v=-y/2+vv-1; H1(uu,vv)=T/(pi*(a*u+b*v+0.001)*sin(pi*(a*u+b*v+0.001)*exp(-pi*j*(a*u+b*v+0.001); end endendQ0=fftlena.*H1; q0=real(ifft2(Q0); q0=uint8(255.*mat2gray(q0); peaksnr=PSNR_cal(lena,q0,8);subplot(2,3,2);imshow(q0);title(sprintf('降质图像, PSNR=%.4f dB',peaksnr);*反向滤波*function inv_Callback(hObject, eventdata, handles)degradation_model;r=1,0.8,0.5,0.2; for i=1:4 inverse_filter(r(i),x,y,i+2,H1,Q0,lena);endfunction = inverse_filter( r , x, y, card, H1, Q0, lena)x_r=round(x*r);y_r=round(y*r);H_r=ones(x,y).*100000;H_r(1:x_r,1:y_r)=H1(1:x_r,1:y_r);Q0_inv=Q0./H_r;q0_inv=uint8(real(ifft2(Q0_inv);peaksnr=PSNR_cal(lena,q0_inv,8); % PSNR% show the imagesubplot(2,3,card);imshow(uint8(255.*mat2gray(q0_inv);title(sprintf('r_0=%.1f , PSNR=%.4f dB',r,peaksnr);end*维纳滤波*function wiener_Callback(hObject, eventdata, handles)degradation_model;HM=abs(H1.*H1);k=0 10(-35) 10(-31) 10(-28); % parameter for wiener filterW=cell(1,4); % wiener filter spectrumfor i=1:4 Wi=HM./(H1.*(HM+k(i);end% run wiener filterfor i=1:4 wiener_filter(i,k(i),Wi,Q0,lena);endfunction = wiener_filter( card, ki, W_k, Q0, lena )Q_k=Q0.*W_k;q_k=real(ifft2(Q_k);% q_k=uint8(255.*mat2gray(q_k);% q_k=uint8(q_k);peaksnr=PSNR_cal(lena,uint8(q_k),8); % PSNRsubplot(2,3,card+2);imshow(uint8(255.*mat2gray(q_k);title(sprintf('K=%f , PSNR=%.4f dB',ki,peaksnr);end实验三 图像分割处理实验1实验目的：（1）了解图像分割的基本原理，并利用图像分割算法进行图像分割处理；（2）掌握数学形态学的基本运算。2实验内容：（1）利用类间方差阈值算法实现图像的分割处理；（2）利用形态学处理进行处理结果修正。3实验要求：（1）实验用图如图3.1所示；图3.1 原始图像（2）对输入图像进行平滑处理，以减小噪声对分割处理的影响；（3）利用类间方差阈值算法对滤波处理后图像进行分割处理，获取分割图像；（4）利用数学形态学中的腐蚀和膨胀运算处理，剔除分割处理结果中的一些细小的残余误分割点，在进行腐蚀和膨胀运算时可采用半径为r的圆形结构元素，注意比较选取不同r值时的处理结果。4实验结果：3.2 类间方差阈值分割3.3 腐蚀运算3.4 膨胀运算可以看出通过类间方差阈值分割，有效的将目标区域和背景区域区分开。对于腐蚀运算，它将邻接背景的目标像素设置为背景像素来达到较小目标面积的目的，r越大，效果越明显。对于膨胀运算，它将邻接背景的背景像素设置为目标像素来达到增加目标面积的目的，r越大，效果越明显。5实验代码：*类间方差阈值分割*function DIVIDE_Callback(hObject, eventdata, handles)dock=rgb2gray(imread('dock.png');dock2=medfilt2(dock); t=Ostu(dock2);dock3=dock2;dock3(dock3<t)=0;dock3(dock3>=t)=255;subplot(1,3,1);imshow(dock);title('原图','color','k');subplot(1,3,2);imshow(dock2);title('中值滤波','color','k');subplot(1,3,3);imshow(dock3);title('类间方差阈值分割','color','k');*腐蚀运算 *function ERODE_Callback(hObject, eventdata, handles)dock=rgb2gray(imread('dock.png');dock2=medfilt2(dock); t=Ostu(dock2);dock3=dock2;dock3(dock3<t)=0;dock3(dock3>=t)=255;subplot(2,2,1);imshow(dock3);title('类间方差阈值分割','color','k');for card=1:3 erode_exp3(card,dock3);endfunction = erode_exp3( card,dock3 )ele=strel('disk',card,8);im=imerode(dock3,ele);subplot(2,2,card+1);imshow(im);title(sprintf('腐蚀: r=%d',card),'color','k');hold onend*膨胀运算 *function DILATE_Callback(hObject, eventdata, handles)dock=rgb2gray(imread('dock.png');dock2=medfilt2(dock); t=Ostu(dock2);dock3=dock2;dock3(dock3<t)=0;dock3(dock3>=t)=255;subplot(2,2,1);imshow(dock3);title('类间方差阈值分割 ','color','k');for card=1:3 dilate_exp3(card,dock3);endfunction = dilate_exp4( card,dock3 )ele=strel('disk',card,8);im=imdilate(dock3,ele);%subplot(2,2,card+1);imshow(im);title(sprintf('膨胀: r=%d',card),'color','k');end实验四 用Hough变换进行曲线的参数提取1实验目的：（1）了解边缘检测算子的原理，并利用边缘算子对图像进行检测；（2）掌握Hough变换的基本原理。2实验内容：（1）分别将原始图像及加高斯噪声、椒盐噪声后的图像中圆形边缘检测出来；（2）用Hough变换对边缘进行参数提取。3实验要求：（1）实验用图像文件：原始图像（houghorg.bmp）、加高斯噪声后图像（houghgau.bmp）和加椒盐噪声后图像（houghsalt.bmp）；图4.1 原始图像 （2）在含有噪声的背景下，先对图像中值滤波，再进行边缘检测；（3）将目标的边界提取出来。边缘检测算子可利用matlab自带函数实现，使用Robert、Sobel和Laplacian算子；（4）利用Hough变换提取的参数绘制曲线，并叠加在噪声图像上。4实验要点：（1）利用算子进行边缘检测：可先将加噪以后的图像进行平滑滤波，如采用9*9的掩膜模板进行中值滤波；为了对图像中图形边缘进行线性提取，可通过设置阈值将图像变为二值图像，再利用三种不同的算子（Robert、Sobel和Laplacian）来完成边缘的检测；（2）Hough变换进行曲线参数提取：在使用三种算子对加噪后图像进行边缘检测以后，使用Hough变换对检测后图像进行参数提取，并在提取成功以后，使用提取获得的参数进行图像的重建，最后将重建图像叠加到加噪图像中。注意在进行Hough变换时，对比观察获得图像与使用算子进行边缘检测获得图像之间的区别5实验结果：4.2 噪声产生4.3 高斯噪声图像的边缘检测及Hough重建4.4 椒盐噪声图像的边缘检测及Hough重建Hough变换检测圆：通过a=x-r*cos(angle),b=y-r*sin(angle)将原图像中的边缘点映射到参数空间（a,b,r）中，由于数字图像且采用极坐标，angle和r都取一定的范围和步长，这样通过两重循环将原图像空间中的点映射到参数空间中，再在参数空间中寻找圆心，然后求出半径坐标。观察实验结果，三种算子都能有效提取边界，但使用不同的算子提取出来的边缘质量有差别。 在噪声的影响下， roberts 算子和 sobel 算子提取边界的效果好于 laplacian 算子，提取结果更精确； laplacian 算子在Hough 变换后重建的边缘与原图像有一定偏离，这可能是由于laplacian算子对噪声点有一定的放大作用。由于此处用Hough 变换提取的是圆，也就将之前算子检测到的图像边缘的直线型边界滤除了，所以在图中仅显示红色部分的圆。6实验代码：*噪声产生*function NOISE_Callback(hObject, eventdata, handles)shape=rgb2gray(imread('houghorg.bmp');gauss=imnoise(shape,'gaussian',0,0.2);pepper=imnoise(shape,'salt & pepper',0.2);subplot(1,3,1);imshow(shape);title('原图','color','k');subplot(1,3,2);imshow(gauss);title('高斯噪声','color','k');subplot(1,3,3);imshow(pepper);title('椒盐噪声','color','k');*高斯噪声图像的边缘检测及Hough重建*function GAUSS_Callback(hObject, eventdata, handles)shape=rgb2gray(imread('houghorg.bmp');gauss=imnoise(shape,'gaussian',0,0.2);gs0=medfilt2(gauss,9,9); gs0(gs0>127)=255;gs0(gs0<128)=0;flag='高斯噪声' op_exp4(gs0,gauss,flag);*椒盐噪声图像的边缘检测及Hough重建*function PEPPERSALT_Callback(hObject, eventdata, handles)shape=rgb2gray(imread('houghorg.bmp');pepper=imnoise(shape,'salt & pepper',0.2);p0=medfilt2(pepper,9,9); % filtered% binarizep0(p0>127)=255;p0(p0<128)=0;flag='椒盐噪声'op_exp4(p0, pepper,flag);*边缘检测及Hough重建*function =op_exp4( p0 , original ,flag)% edge detectps=cell(1,3);ps1=edge(p0,'roberts'); subplot(2,3,1); imshow(ps1); title(flag,':','roberts算子','color','k');ps2=edge(p0,'sobel'); subplot(2,3,2); imshow(ps2); title(flag,':','sobel算子','color','k');ps3=edge(p0,'log'); subplot(2,3,3); imshow(ps3);title(flag,':','laplacian算子','color','k'); for i=1:3 restore_exp4(psi,i,original);endendfunction =restore_exp4( p, card, original )par1,par3=Hough(p);subplot(2,3,card+3);imshow(original);viscircles(par1,par3);if card=1 op='roberts算子'elseif card=2 op='sobel算子'else op='laplacian算子'endtitle('重建图像:',op,'color','k');endfunction par1, par3 = Hough( BW )r_max=100;r_min=40;step_r=1;step_angle=pi/20;p=0.5;m,n = size(BW);size_r = round(r_max-r_min)/step_r)+1;size_angle = round(2*pi/step_angle);hough_space = zeros(m,n,size_r);rows,cols = find(BW);ecount = size(rows);for i=1:ecount for r=1:size_r for k=1:size_angle a = round(rows(i)-(r_min+(r-1)*step_r)*cos(k*step_angle); b = round(cols(i)-(r_min+(r-1)*step_r)*sin(k*step_angle); if(a>0 && a<=m && b>0 && b<=n) hough_space(a,b,r) = hough_space(a,b,r)+1; end end endendmax_para = max(max(max(hough_space);index = find(hough_space>=max_para*p);length = size(index);hough_circle = false(m,n); for i=1:ecount for k=1:length par3 = floor(index(k)/(m*n)+1; par2 = floor(index(k)-(par3-1)*(m*n)/m)+1;%y par1 = index(k)-(par3-1)*(m*n)-(par2-1)*m; %x if(rows(i)-par1)2+(cols(i)-par2)2<(r_min+(par3-1)*step_r)2+5 && (rows(i)-par1)2+(cols(i)-par2)2>(r_min+(par3-1)*step_r)2-5) hough_circle(rows(i),cols(i) = true; end endendfor k=1:length par3 = floor(index(k)/(m*n)+1; par2 = floor(index(k)-(par3-1)*(m*n)/m)+1;%y par1 = index(k)-(par3-1)*(m*n)-(par2-1)*m; %x par3 = r_min+(par3-1)*step_r; %r para(:,k) = par1,par2,par3;endpar1=par2 par1;end实验五 手写数字识别1实验目的： 掌握分类、识别问题的实质，了解各种分类问题的机器学习方法，并至少掌握一种。2实验内容：对实验提供的手写数据库（包含数字“1”，“6”，“8”）进行训练和测试，最终能够较为准确的识别数据库中的手写体数字。3实验要求：编写一完整的MATLAB程序，选取一种合适的机器学习方法，对实验提供的手写数据库（包含数字“1”，“6”，“8”）进行训练和测试，最终能够较为准确的识别数据库中的手写体数字。要求程序输出：（1）该算法的精确度，并且不能低于75%；（2）测试集中所有的测试结果（可用fprintf显示出来）。提示：可以选取以下几种方法：（1）利用数字的集合几何形状的特点，计算每幅图的连通域，来进行分类识别；（2）逻辑回归算法；（3）支持向量机（SVM）；4实验结果：集体预测结果：数字观测结果：本实验通过神经网络算法设计。由实验结果可以看出，迭代训练的次数越多，网络构建越完善，正确识别率越高。5实验代码：*神经网络训练*function prepare_key_Callback(hObject, eventdata, handles)close(figure);NNtraindata;* NNtraindata*global T;global labda;input_layer_size = 784;hidden_layer_size = 25; num_labels = 10; load('traindata.mat');m = size(X, 1);X = double(X);sel = randperm(size(X, 1);sel = sel(1:100);displayData(X(sel, :);load('ex5weights.mat');nn_params = Theta1(:) ; Theta2(:);lambda = 0;J = nnCostFunction(nn_params, input_layer_size, hidden_layer_size, . num_labels, X, y, lambda);lambda = 1;J = nnCostFunction(nn_params, input_layer_size, hidden_layer_size, . num_labels, X, y, lambda);g = sigmoidGradient(1 -0.5 0 0.5 1);initial_Theta1 = randInitializeWeights(input_layer_size, hidden_layer_size);initial_Theta2 = randInitializeWeights(hidden_layer_size, num_labels);initial_nn_params = initial_Theta1(:) ; initial_Theta2(:);checkNNGradients;lambda = 3;checkNNGradients(lambda);% Also output the costFunction debugging valuesdebug_J = nnCostFunction(nn_params, input_layer_size, . hidden_layer_size, num_labels, X, y, lambda);% After you have completed the assignment, chan