数字图像处理实验报告.docx

数字图像处理实验报告班级：姓名:学号：日期：目录实验1 直方图均衡4问题描述4程序说明4运行结果5结果分析6小结6实验2 图像几何变换7问题描述7程序说明7运行结果9结果分析10小结10实验3 空间滤波11问题描述11程序说明11运行结果13结果分析13小结14实验4 傅里叶变换14问题描述14程序说明14运行结果18结果分析18小结18实验5 图像复原19问题描述19程序说明19运行结果21结果分析22小结23实验6 图像压缩23问题描述23程序说明23运行结果26结果分析26小结26实验7 边缘检测27问题描述27程序说明27运行结果30结果分析30小结31心得体会31致谢31参考文献31实验1 直方图均衡问题描述histogram equalization (test images: fig1.jpg, fig2.jpg)(a) write a computer program for computing the histogram of an image.(b) implement the histogram equalization technique.(c) your program must be general to allow any gray-level image as its input. as a minimum, your report should include the original image, a plot of its histogram, a plot of the histogram-equalization, transformation function, the enhanced image, and a plot of its histogram.程序说明在此问题中通过编写直方图均衡函数，读取输入图像的各像素点灰度值，计算直方图，再做直方图均衡操作并绘图，具体函数histogram如下：function orihist, trans, enhimg, enhhist = histogram( oriimg, graynum )%histogram 直方图均衡函数% oriimg 原图片% graynum 图片灰度级数% orihist 原图片直方图% trans 灰度变换数组% enhimg 均衡后图片% enhhist 均衡后直方图 height width = size(oriimg);orihist = zeros(1, graynum);pixelnum = height * width; for i = 1:graynum temp = oriimg = i - 1; orihist(i) = sum(temp(:) / pixelnum;endtrans = uint8(cumsum(orihist) * graynum);enhimg = trans(int16(oriimg) + 1);enhhist = full(sum( . sparse(1:graynum, double(trans) + 1, orihist, graynum, graynum); end算法说明：其中的for循环计算出原图片的直方图。下面一行通过累加（cumsum）计算出均衡化变换曲线。再下面一行通过变换曲线进行直方图均衡化。最后一行计算出均衡化后的直方图。最后，按题目要求绘出各图线。运行结果fig1处理结果对比：fig2处理结果对比：结果分析fig1中原图片的背景较暗，导致骨骼图片的细节无法显示，通过直方图均衡的方法，可以看到处理后的图片骨骼的细节很好的显示出来，同时可以看到处理后的直方图区域不再集中在灰度值较低的部分。fig2中原图片背景与需要观察的花粉粒区分度不够高，不能很好地观察到图片的细节，通过直方图均衡后得到高对比度的图片，可以看到花粉粒的细节被很好地突出出来。同时在直方图中也可以看到原图片的直方图分布主要集中在100附近，而均衡后的直方图较为均匀地分布在0255之间，故起到增强图像对比度的效果。小结直方图均衡化是图像处理领域中利用图像直方图对对比度进行调整的方法。这种方法通常用来增加许多图像的局部对比度，尤其是当图像的有用数据的对比度相当接近的时候。通过这种方法，亮度可以更好地在直方图上分布。这样就可以用于增强局部的对比度而不影响整体的对比度，直方图均衡化通过有效地扩展常用的亮度来实现这种功能。这种方法对于背景和前景都太亮或者太暗的图像非常有用，这种方法尤其是可以带来x光图像中更好的骨骼结构显示以及曝光过度或者曝光不足照片中更好的细节。在fig1中有很好的体现。另一方面，种方法的一个缺点是它对处理的数据不加选择，它可能会增加背景杂讯的对比度并且降低有用信号的对比度，同样在fig1处理后的图片中可以看到增加了很多黑斑，但并不影响骨骼的观测效果。实验2 图像几何变换问题描述geometric transform (test image: fig3.tif)develop a geometric transform program that will rotate, translate, and scale an image by specified amounts, using the nearest neighbor and bilinear interpolation methods, respectively.程序说明在此问题中的图像的翻转和拉伸变换均可以通过对图像像素矩阵的像素点变换加以实现，在这里可以乘上相应的翻转拉伸矩阵。并分别使用最近差值和双线性差值的方法得到几何变换后的图片。相应的变换矩阵生成函数如下：旋转变换：function a = rotatemat( angle )%rotatemat % angle a = cos(angle),-sin(angle),0;sin(angle),cos(angle),0;0,0,1; end拉伸变换：function a = scalemat( scalex, scaley )%scalemat % scalex % scaley a = scalex,0,0;0,scaley,0;0,0,1; end平移变换：function a = translatemat( offsetx, offsety )%translatemat % offsetx % offsety a = 1,0,offsetx;0,1,offsety;0,0,1; end最近差值和双向性差值算法：% 最近差值function value = nearestneighbor(oriimg, x, y) oriheight oriwidth = size(oriimg);if x oriwidth | y oriheight value = 0;else value = oriimg(round(y), round(x);end end % 双线性差值function value = bilinear(oriimg, x, y) x1 = floor(x);x2 = ceil(x);y1 = floor(y);y2 = ceil(y);up = (getvalue(oriimg, x2, y1) - getvalue(oriimg, x1, y1) * (x - x1) + getvalue(oriimg, x1, y1);down = (getvalue(oriimg, x2, y2) - getvalue(oriimg, x1, y2) * (x - x1) + getvalue(oriimg, x1, y2);value = (down - up) * double(y - y1) + up; end运行结果通过运行相应的变换得到如下结果：原图片：处理后图片：说明：上面两幅图片是用最近差的方法经过旋转平移和拉伸平移得到，下面两幅图片是用双线性差值的方法经过旋转平移和拉伸平移得到结果分析从上图可以看出，平移、拉伸、旋转变换都很正确。就图像平滑度方面来说，双线性差值的结果比最近邻差值的结果更平滑。由于在计算变换结果像素时，只在原图像中采样一个点，所以图像在放大的情况下，双线性差值的结果才比最近邻差值更平滑，而在图像缩小的情况下，两种算法结果类似。小结最邻近差值就是输入图像中与其最邻近的采样点的像素值，但是此方法会在图像中产生认为是加工的痕迹（锯齿形状）。对比左列两幅图片可以再边缘部分发现最近差值方法图片的边缘与双线性差值方法得到的图片相比也确实出现了锯齿状。而双向性差值会使图像显得平滑，然而也会使图像的细节退化。实验3 空间滤波问题描述spatial filtering (test image: fig4.jpg)(a).noise generators. 1. find (or develop) a program to add gaussian noise to an image. you must be able to specify the noise mean and variance; 2. find (or develop) a program to add salt-and-pepper (impulse) noise to an image. you must be able to specify the probabilities of each of the two noise components.(b) write program to perform spatial mean filtering, applying your filter to the noisy images obtained in (a). you can fix the size of the spatial mask at 3 x 3, but the coefficients need to be variables that can be input into your program.(c) modify the program that you developed in (b) to perform 3 x 3 median filtering. compare the differences between filtered images by mean and median filters.程序说明高斯噪声生成函数：function b = gaussiannoise( a, mean, var)% gaussiannoise 高斯噪声函数% a 输入信号% mean 噪声均值% var 噪声方差% b 输出信号 b = double(a) + randn(size(a) * sqrt(var) + mean; end椒盐噪声生成函数：function b = saltpeppernoise( a, s, p )% saltpeppernoise椒盐噪声生成函数% a 输入信号% s 盐噪声概率% p 椒噪声概率% b 输出信号 r = rand(size(a);b = double(a) + (r 1 - p) * 999999; end均值滤波函数：function b = meanfilter( a, m )% meanfilter 均值滤波函数% a 输入信号% m 3*3掩膜矩阵% b 输出信号 m n = size(a);temp = zeros(m+2,n+2);for i = -1:1 for j = -1:1 temp(2+i):(m+1+i),(2+j):(n+1+j) =. temp(2+i):(m+1+i),(2+j):(n+1+j) + double(a) * m(2-i,2-j); endendb = temp(2:m+1,2:n+1); end中值滤波函数：function b = medianfilter( a )%medianfilter 中值滤波% a 输入信号% b 输出信号 m n = size(a);temp = zeros(m+2,n+2,9);k = 1;for i = -1:1 for j = -1:1 temp(2+i):(m+1+i),(2+j):(n+1+j),k) = a; k = k + 1; endendtemp2 = median(temp, 3);b = temp2(2:m+1,2:n+1);在该程序中，通过运行高斯噪声函数和椒盐噪声函数添加方差500高斯噪声与概率(0.1, 0.1)的椒盐噪声，并通过均值滤波函数和中值滤波函数滤波。运行结果结果分析从上图中对比可以发现，对于高斯噪声，均值滤波与中值滤波的效果类似，都较好；而对于椒盐噪声，均值滤波效果较差，而中值滤波效果较好。小结对于高斯噪声，出现位置是一定的，但噪声幅值是随机的；而椒盐噪声出现位置是随机的，但噪声的幅值是基本相同的。而中值滤波是选择适当的点来代替污染点的值，所以在处理椒盐噪声时效果会比均值滤波好。实验4 傅里叶变换问题描述two-dimensional fast fourier transform (test image: fig5.jpg)consult a reference book for fast fourier transform (fft), and then develop a program that can compute and display the two-dimensional discrete fourier transform (amplitude and phase spectra) of a digital image. your implementation must have the capabilities to:(a) multiply the input image by (-1)x+y to center the transform.(b) compute and display the amplitude spectrum and phase spectrum.(c) compute the inverse fourier transform.(d) use the image fig5.jpg to illustrate your algorithm.程序说明这个算法的实现为fft2d函数。它先做中心化变换，然后分别对图像的行与列进行一维fft运算，随后就得到中心化的二维fft结果。先读入fig5.jpg图像，然后用fft2d函数对它进行二维fft运算，得到一个复数图像。通过傅里叶逆变换ifft2d函数对它进行二维fft逆运算，画出复原图像。最后分别画出这个复数矩阵的幅度矩阵与相位矩阵，即原图的幅度谱与相位谱。其中绘制幅度谱前将这个矩阵取了对数，这样避免了图像中明暗变化过大，导致只能看到图像中心的亮点，而看不出其余部分的明暗变化。这两幅图绘制时，都将图像亮度归一化到了0到1之间。二维快速傅里叶变换函数：function b = fft( a )%fft 快速傅里叶变换函数% a 输入矩阵% b 输出矩阵 n = length(a);m = log2(n); t = zeros(1, m);mark = zeros(1, n);for k = 0:1:n-1 h = k; for j = m-1:-1:0 temp = power(2, j); t(j+1) = floor(h / temp); h = rem(h , temp); end t2 = zeros(1, m); for j = 0:1:m-1 if t(j+1) = 1 t2(m-j) = 1; end end sum = 0; for j = m-1:-1:0 sum = sum*2 + t2(j+1); end if mark(sum+1) = 0 temp = a(k+1); a(k+1) = a(sum+1); a(sum+1) = temp; mark(k+1) = 1; endend for k = 1:m l = power(2, k); r = n / l; l2 = l / 2; for s = 0:l2-1 w = cos( 2*pi*s / l ) - 1i*sin( 2*pi*s / l); for j = 0:1:r-1 t = w * a(j*l + s + l2 + 1); a(j*l + s + l2 + 1) = a(j*l + s + 1) - t; a(j*l + s + 1) = a(j*l + s + 1) + t; end endend b = a;end快速傅里叶变换逆变换函数：function b = ifft2d( a )% ifft2d 二维快速傅里叶变换逆变换% a 输入矩阵% b 输出矩阵 m n = size(a);b = zeros(m, n); for j = 1:n b(:, j) = ifft(a(:, j);end for j = 1:m b(j, :) = ifft(b(j, :);end end function b = ifft( a ) n = length(a);m = log2(n); t = zeros(1, m);mark = zeros(1, n);for k = 0:1:n-1 h = k; for j = m-1:-1:0 temp = power(2, j); t(j+1) = floor(h / temp); h = rem(h , temp); end t2 = zeros(1, m); for j = 0:1:m-1 if t(j+1) = 1 t2(m-j) = 1; end end sum = 0; for j = m-1:-1:0 sum = sum*2 + t2(j+1); end if mark(sum+1) = 0 temp = a(k+1); a(k+1) = a(sum+1); a(sum+1) = temp; mark(k+1) = 1; endend for k = 1:m l = power(2, k); r = n / l; l2 = l / 2; for l = 0:1:l2-1 w = cos( 2*pi*l / l ) + 1i*sin( 2*pi*l / l); for j = 0:1:r-1 t = w*a(j*l + l + l2 + 1); a(j*l + l + l2 + 1) = a(j*l +l + 1) - t; a(j*l + l + 1) = a(j*l + l + 1) + t; end endend b = 1/n .* a; end运行结果结果分析从幅度谱中可以看出，有一道纵向的亮线，这应该对应原图像中的水平波纹噪声与水面横向纹理。其余线条可能与原图像中的三角架等明显直线边缘对应。小结快速傅里叶变换能克服时间域与频率域之间相互转换的计算障碍，在光谱、大气波谱分析、数字信号处理等方面有广泛应用。上图中正是据此得出图像的幅度谱与相位谱。实验5 图像复原问题描述image restoration(test image: fig6.jpg) suppose a blurring degradation function as (1)(a) implement a blurring filter using eq. (1).(b) blur the test image fig6.jpg using parameters a=b=0.1 and t = 1.(c) add gaussian noise of 0 mean and variance of 650 to the blurred image.(d) restore the blurred image and the blurred noisy image using the inverse filter, wiener deconvolution filter and the parametric wiener filter, respectively.(e) add gaussian noise of 0 and different variances to the blurred image and repeat (d), investigate the performance of the wiener deconvolution filter.程序说明h.m函数定义了上面那个传递函数。filter函数定义了频域中的正过滤与反过滤。wiener函数输入原图像的频域矩阵、噪声矩阵、原图像与算法参数，最后输出图像。运行程序可以演示上述算法。这个程序会读入fig6.jpg图像，然后进行过滤，加噪声，最后分别使用反过滤、维纳滤波、参数维纳滤波进行图片复原。这里选取了方差为0, 10, 650的高斯噪声分别进行了测试。过滤函数：function b = filter( a, inv )% a 输入信号% inv 0表示反过滤，其它表示正过滤 if exist(inv,var) inv = 0;end m n = size(a);b = zeros(m, n);if inv = 0 for i = 1:m for j = 1:n b(i, j) = a(i, j) * h(i, j, 0.1, 0.1, 1); end endelse for i = 1:m for j = 1:n b(i, j) = a(i, j) / h(i, j, 0.1, 0.1, 1); end endend end维纳滤波函数：function b = wiener( a, noise, ori, arg )%wiener 维纳滤波函数% a 输入信号% noise 噪声% ori 原图像% arg 参数% b 输出图像 m n = size(a);b = a;fnoise = fft2(noise);forig = fft2(double(ori);ratio = arg .* (conj(fnoise) .* fnoise) ./ (conj(forig) .* forig);for i = 1:m for j = 1:n temp = conj(h(i, j, 0.1, 0.1, 1) * h(i, j, 0.1, 0.1, 1); b(i, j) = b(i, j) .* conj(h(i, j, 0.1, 0.1, 1) ./ abs(temp + ratio(i, j); endend b = ifft2(b); end运行结果0方差噪声下，图片复原效果图：10方差噪声下，图片复原效果图：650方差噪声下，图片复原效果图：结果分析从结果看出，在0方差噪声时，三种方法均可以很好地复原图像。但是，反过滤对噪声的敏感度很大，稍有噪声则无法复原图像，而维纳滤波对噪声的敏感度较小，可以在一定噪声条件下复原图像。然而在噪声过大的情况下，复原后的图像会模糊不清。小结从本实验中可以看出维纳滤波较反过滤方法在有噪声的情况下有很好的效果，但维纳滤波要求输入过程是广义平稳的，并且.输入过程的统计特性是已知的。由于输入过程取决于外界的信号、干扰环境，这种环境的统计特性常常是未知的、变化的，因而难以满足上述两个要求。实验6 图像压缩问题描述image compression(a)develop a program to compute the two-dimensional discrete cosine transform of an image(b)compress the test image to different qualities by using the inverse discrete cosine transform with fewer transform coefficients.(c)compute and display the difference image.程序说明dct2函数实现了二维dct变换，idct2函数实现了二维dct逆变换。 compress函数实现了对图像的压缩。它将图像按88分成小块，对每个小块分别进行dct变换，然后将变换后的结果按照参数q进行量化。q越大，即量化单位越大，则压缩率越高，图片质量越差。二维dct变换函数：function b = dct2( a )m n = size(a);b = zeros(m, n); for i = 1:m b(i, :) = dct(a(i,:);end for i = 1:n b(:, i) = dct(b(:,i);end end function b = dct( a ) a = double(a);m n = size(a);if n = 1 a = a;else m = n;end; b = a(1:2:m-1),(m:-2:2);b = fft(b);temp = exp(-1i * (0:m-1) * pi / (2 * m);b = b .* temp;b = real(b);b(1) = b(1) * sqrt(1/m);b(2:m) = b(2:m) * sqrt(2/m); end图像压缩函数：function b = compress( a, q )%compress 图像压缩函数% a 被压缩图片% q 压缩质量% b 压缩后的信息矩阵 m n = size(a);b = double(a); for i = 1:8:m-7 for j = 1:8:n-7 s = b(i:i+7, j:j+7); b(i:i+7, j:j+7) = dct2(s); b(i:i+7, j:j+7) = round( b(i:i+7, j:j+7) ./ q); endend end解压缩函数：function b = extract( a, q )%extract 解压缩函数% a 压缩的图片信息% q 压缩质量% b 解压后的图片 m n = size(a);b = double(a); for i = 1:8:m-7 for j = 1:8:n-7 s = b(i:i+7, j:j+7); s = round( s .* q ); b(i:i+7, j:j+7) = idct2(s); endend end它首先读取出fig3.tif的图像，然后分别进行量化参数为10, 30, 80 的压缩与解压缩，最后显示出压缩结果。压缩结果中包括了：解压缩后的图片、原图片与解压图片的差异。这里差异图片用0.5灰色来表示没有差异，深色与浅色分别表示正负差异。运行结果结果分析从不同的参数的压缩结果来看，压缩率越高，图像的质量就越差，即与原图像的差距越大。差异主要在图像的细节部分。每个分块的主要明暗变化得以保留，细节变化被去除。小结图像压缩是指以较少的比特有损或无损地表示原来的像素矩阵的技术。在此实验中主要是有损压缩的方式，在压缩率较高时，则损失了较多信息，在解压缩时也出现较大的失真，但无损压缩有时有无法实现较高的压缩率。实验7 边缘检测问题描述edge detection (test image: fig3.tif)develop a program to implement the laplacian, roberts, sobel and the prewitt edge detection algorithms.程序说明由于这里的算法都是空间域中的掩膜滤波，在遍历图片像素前需要用extendimage函数对原图片进行扩展。extendimage函数会对输入图像外围增加一圈像素，并使新加的像素尽量接近周围的像素。laplacian算法使用的掩膜矩阵为：0-10-14-10-10其余三种算法使用的掩膜矩阵是成对的，计算时需要取绝对值的和。roberts算法使用的掩膜矩阵为：-1001, 0-110sobel算法使用的掩膜矩阵为：-1-2-1000121, -101-202-101prewitt算法使用的掩膜矩阵为：-1-1-1000111, -101-101-101laplacian算法函数：function b = laplacian( a )%laplacian laplacian函数 m n = size(a);c = extendimage(a);lapmask = 0 -1 0 -1 4 -1 0 -1 0 ;b = zeros(m, n);for i = 2:m+1 for j = 2:n+1 temp = lapmask .* c(i-1:i+1, j-1:j+1); b(i-1,j-1) = abs(sum(temp(:); endendendsobel算法函数：function b = sobel( a )%sobel sobel算法函数 m n = size(a);c = extendimage(a);sobelmask1 = -1 -2 -1 0 0 0 1 2 1 ;sobelmask2 = -1 0 1