生物电子与影像技术实验指导书2013.doc

生物电子与影像技术实验指导生物电子与影像技术实验指导哈尔滨工业大学（威海）控制科学与工程系实验一：图像处理基本操作与灰度变换1. 基础知识1) Matlab支持的图像格式TIFF, JEPG, GIF, BMP, PNG, XWD (X Window Dump)，其中GIF不支持写。2) 与图像处理相关的最基本函数读：imread; 写：imwrite; 显示：imshow; 信息查看：imfinfo; 3) Matlab支持的数据类型double, unit8, int8, uint16, int16, uint32, int32, single, char (2 bytes per element), logical.4) Matlab支持的图像类型Intensity images, binary images, indexed images, RGB image5) 数据类及图像类型间的基本转换函数数据类转换：B = data_class_name(A);IPT图像数据类型转换函数：函数名输出图像数据类型输入图像数据类型im2uint8uint8logical, uint8, uint16, doubleim2uint16uint16logical, uint8, uint16, doublemat2graydouble(0, 1)doubleim2douledoublelogical, uint8, uint16, doubleim2bwlogical logical, uint8, uint16, double6) Matlab的灰度直方图函数显示 imhist, 均衡化 histeq7) Matlab的图像灰度变换图像滤波 imfilter (image, filter) 灰度变换 imadjust2. 图像信息显示与图像读写 imfinfo(parrots.bmp) %图像文件信息显示 Filename: parrots.bmp FileModDate: 27-Sep-2010 19:42:59 FileSize: 1179702 Format: bmp FormatVersion: Version 3 (Microsoft Windows 3.x) Width: 768 Height: 512 BitDepth: 24 ColorType: truecolor imfinfo(pollen.tif) %图像文件信息显示 Filename: pollen.tif FileModDate: 27-Sep-2010 19:42:59 FileSize: 250398 Format: tif FormatVersion: Width: 500 Height: 500 BitDepth: 8 ColorType: grayscale img=imread(parrots.bmp); %读入图像文件到内存变量 imshow(img); %图像显示 whos img %显示变量信息 Name Size Bytes Class img 512x768x3 1179648 uint8 arrayGrand total is 1179648 elements using 1179648 bytes imwrite(img, NewParrots.png); %以其它格式另外保存图像 img2=imread(pollen.tif); whos img2 Name Size Bytes Class img2 500x500 250000 uint8 arrayGrand total is 250000 elements using 250000 bytes figure, imshow(img2);3. 直方图显示与均衡化 im=imread(pollen.tif); %读入图像文件 size(im) %观察2维还是3维ans = % 如果是3维，使用 im=im(:,:,1); 转成2维 500 500 subplot(221), imshow(im) %显示原图像 subplot(222), imhist(im) %显示直方图 im2=histeq(im); %直方图均衡化 subplot(223), imshow(im2) %显示均衡化后的图像 subplot(224), imhist(im2) %显示均衡化后的图像的直方图4. 灰度变换函数形式： J = imadjust(I,low_in; high_in,low_out; high_out, Gamma)参数 Gamma用于指定变换函数形状，小于1时，输出更亮，大于1时，输出更暗。 srcimage=imread(chest_xray.tif); %读入图像文件 g1=imadjust(srcimage, 0 1, 0 1, 0.5); %灰度变换 g2=imadjust(srcimage, 0 1, 0 1, 2.5); subplot(1, 3, 1); imshow(srcimage); title (Source Image); subplot(1, 3, 2); imshow(g1); title (Gamma=0.5); subplot(1, 3, 3); imshow(g2); title (Gamma=2.5);实验二：频域图像处理1. 基础知识1) 二维离散Fourier变换对设图像f(x,y)的大小为(M,N), 其二维离散Fourier变换对2) 二维离散Fourier谱的计算与可视化为了便于观察，一般取Fourier谱的一个整周期，并且平移到坐标轴的正向，如下图所示。二维离散Fourier变换的FFT计算函数：F = fft2 (f, P, Q)F: 输出大小为P x Q, 通过对输入图像f添加需要数目的0实现二维离散Fourier变换的Fourier谱计算S = abs ( F )二维离散Fourier变换的原点移动Fc = fftshift (F)Log变换S2=log(1+abs(Fc)2. Fourier谱的计算与显示close allclear all im=imread(Xraychest.png); % 读入图像文件 imdouble=double(im); % fft2函数不处理uint8类型数据，转换成double F=fft2(imdouble); % 2D FFT计算 Fc=fftshift(F); % 平移 S2=log(1+abs(Fc); % Log变换 subplot(121);imshow(im); % 显示原图像 subplot(122);imshow(S2,); % 显示Fourier谱3. Fourier逆变换计算 F = ifftshift(Fc); % 对于Fc要使用平移复原， 对于F则不需要平移 f=real(ifft2(F); % 2D FFT逆变换 imifft=uint8(f); % 转换成uint8类型 figure, imshow(imifft); % 显示图像4. 频域滤波频域滤波基本步骤：输入图像-图像预处理-Fourier变换-滤波- Fourier反变换-图像后处理-输出图像F = fft2(f, size(H, 1), size(H, 2);g = real(ifft2(H.*F);g = g(1:size(f, 1), 1:size(f, 2);5. 空域和频域滤波转换% 显示原图像及其Fourier谱close allclear allf = imread(bld.tif);F = fft2(double(f);S = fftshift(log(1+abs(F);figure (1)subplot(221);imshow(f);title(Original image);subplot(222);imshow(S, );title( Fourier spectrum);% 生成空域和频域滤波器h = fspecial(sobel); % 使用Sobel 滤波器PQ = paddedsize(size(f); % 自编函数，计算填0后的图像大小H = freqz2(h, PQ(1), PQ(2); % 生成 Sobel 算子的频率滤波器H1 = ifftshift(H); % 频率滤波器原点移位到左上角figure(2);subplot(121);mesh(abs(H(1:20:end, 1:20:end); %绘制频率滤波器axis off;title(perspective plot of frequency filter);subplot(122);mesh(abs(H1(1:20:end, 1:20:end); %绘制移位后的频率滤波器axis off;title(frequency filter with shift);% 分别在空域和频域滤波gs = imfilter(double(f), h); %空域滤波F = fft2(double(f), size(H1, 1), size(H1, 2); %频域滤波gf = real(ifft2(H1.*F);gf = gf(1:size(f, 1), 1:size(f, 2);figure(1);subplot(223);imshow(gs, );title(Filtered image in spatial domain);subplot(224);imshow(gf, ); title(Filtered image in frequency domain);6. Gaussian低通和高通频域滤波Gaussian低通滤波器的传递函数为：其中是标准方差，是从点 (u, v) 到滤波器中心的距离。令，得到关于截止参数的表达式：Gaussian低通频域滤波close allclear allf = imread(FigMark.tif);PQ = paddedsize(size(f);U, V = dftuv(PQ(1), PQ(2); % 生成距离计算需要的网格D0 = 0.05*PQ(2); % D0取宽度的5%, 根据需要调整F = fft2(double(f), PQ(1), PQ(2); H = exp(-(U.2+V.2) / (2*(D02); % 滤波器计算g = real(ifft2(H.*F); % 频域滤波g = g(1:size(f, 1), 1:size(f, 2);figure % 结果显示subplot(221); imshow(f, );title(Original image);subplot(222); imshow(log(1 + abs(fftshift(F), );title( Fourier spectrum of original image);subplot(223); imshow(g, );title(Filtered image by Gaussian lowpass filter);F2= fft2(double(g), PQ(1), PQ(2);subplot(224); imshow(log(1 + abs(fftshift(F2), );title( Fourier spectrum of lowpass filtered image);记低通滤波器的传递函数为，则高通滤波器的传递函数为Gaussian高通频域滤波H2 =1-H; % 高通滤波器计算g2 = real(ifft2(H2.*F); % 频域滤波g2 = g2(1:size(f, 1), 1:size(f, 2);figure % 结果显示subplot(221); imshow(f, );title(Original image);subplot(222); imshow(log(1 + abs(fftshift(F), );title( Fourier spectrum of original image);subplot(223); imshow(g2, );title(Filtered image by Gaussian highpass filter);FH2= fft2(double(g2), PQ(1), PQ(2);subplot(224); imshow(log(1 + abs(fftshift(FH2), );title( Fourier spectrum of highpass filtered image);实验三：图像的恢复与分割1. 基础知识1) 图像恢复设真实图像为f(x,y)，其退化和恢复模型如下图所示。图像恢复就是利用某种先验知识对退化过程建模，并在某个质量准则下，应用其逆过程来得到原图像的最优估计。实际的恢复过程相当于设计一个滤波器，使其能从降质图像计算得到真实图像的估值，根据预先规定的误差准则，最大程度接近真实图像。2) 噪声模型对噪声的行为和效果进行仿真的能力是图像恢复处理中的重点。除了周期噪声，通常假定噪声与像素坐标无关。Matlab中噪声仿真函数为：g = imnoise(f, type, parameters)，常见的高斯、椒盐噪声的仿真实例如下：g = imnoise(f, gaussian, m, var);g = imnoise(f,salt & pepper, d); (其中 d 代表噪声密度)Matlab中两个常用的随机分布产生函数是：rand（均匀分布噪声）, randn（标准正态分布噪声）。常用的噪声模型包括：Gaussian噪声分布：Rayleigh噪声分布：脉冲噪声分布：噪声的估计主要是对参数均值和方差和的估计。可以通过选择图像的一块小的背景灰度近似为常数的区域估计均值和方差，再求出a和b。选择图像区域可以通过Matlab函数roipoly进行：BW = roipoly(I,c,r): 返回图像I中由c, r确定的区域，结果是一个二值图像，其大小和图像I相同，区域内为1，区域外为0。参数c和r保存组成该区域多边形的各个顶点的坐标对。BW, x, y=roipoly(I); 允许用户交互式地通过鼠标选择图像I的多边形区域，结果返回多边形坐标对x,y以及二值区域BW。均值和方差的估计可以通过自编函数statmoments完成。3) 仅有噪声的图像恢复对于仅有噪声影响的图像，可以采用空域滤波进行图像恢复。常用的空域滤波包括：算术均值滤波器：几何均值滤波器：谐波均值滤波器：逆谐波均值滤波器：，其中Q为滤波器的阶数，Q=0时成为算术均值滤波，Q0 时对于盐噪声处理效果较好，Q1 对于椒噪声处理效果较好，Q=-1时成为谐波均值滤波。自编函数spfilt可以完成常用的空域滤波工作：f=spfilt (img, type, m, n) img输入图像，m,n滤波窗口大小，如3x3 typeamean:算术均值滤波, gmean: 几何均值滤波器, hmean: 谐波均值滤波器, median: 中值滤波 对于逆谐波均值滤波器, type设置为chmean, 同时还应该加上参数Q, 用法如下f=spfilt (img, type, m, n, Q)2. 噪声的产生和分布的显示噪声可以直接用自编函数imnoise2来产生，分布可以使用直方图进行观察。g=imnoise2 (type, M, N, a, b )Type uniform, gaussian, salt & pepper, rayleigh, exponential, erlangM, N 数组大小a, b 对应type的参数g 返回的噪声数组例如：R=imnoise2(gaussian,500,500,0,1); %产生Gaussian噪声 500x500像素, 均值为0, 方差为1figure, hist(R(:), 50); % 显示直方图 50为统计的区间个数title(Gaussian Noise);R2=imnoise2(rayleigh, 100000, 1); % 产生Rayleigh噪声figure, hist(R2(:), 50);title(Rayleigh Noise);3. 仅有噪声影响图像的恢复噪声可以直接用自编函数imnoise2来产生，空域降噪可以采用函数spfilt完成。close allclear allim=imread(bld.tif);subplot(231);imshow(im); title (Source image);% 计算噪声图像imnoise=imnoise2(gaussian,size(im,1), size(im,2),0,0.1); % 噪声参数设置im2=double(im)+255.0*imnoise; % 加性噪声影响subplot(232);imshow(im2,); title (Noise image);% 进行空域滤波f=spfilt (im2,amean , 3, 3); % 算术均值滤波subplot(233);imshow(f,); title (Arithmetic mean);f=spfilt (im2,gmean , 3, 3); % 几何均值滤波subplot(234);imshow(f,); title (Geometric mean);f=spfilt (uint8(im2),hmean , 3, 3); % :谐波均值滤波subplot(235);imshow(f,); title (Harmonic mean);f=spfilt (im2,median , 3, 3); % 中值滤波subplot(236);imshow(f,); title (Median);4. 含椒盐噪声图像的去噪椒盐噪声中椒(黑色)和盐(白色)可以分别设置强度（即在图像中的占的比例，注意椒和盐的强度之和应该小于1），也可以仅出现椒噪声或者盐噪声。中值滤波和逆谐波滤波都经常用于椒盐噪声的消除。% 产生椒盐噪声图像的方法： %close allclear allf=imread(cktboard.tif);%产生椒盐噪声M, N=size(f);R=imnoise2(salt & pepper,M, N, 0.1, 0.1); % 噪声强度 20%-Peppers 20%-Saltc = find(R=0);g2 = f;g2(c)=0;c = find(R=1);g2(c)=255;%产生椒噪声R=imnoise2(salt & pepper,M, N, 0.2, 0); %噪声强度30%-Peppersc = find(R=0);gp = f;gp(c)=0;%产生盐噪声R=imnoise2(salt & pepper,M, N, 0,0.2); %噪声强度30%-Saltc = find(R=1);gs = f;gs(c)=255;figure(1)subplot(241); imshow(f); title(source);subplot(242); imshow(g2); title(salt & pepper);subplot(243); imshow(gp); title(pepper);subplot(244); imshow(gs); title(salt);% 中值滤波降噪 %f1=spfilt (g2,median , 3, 3); % 中值滤波,窗口大小3x3f2=spfilt (gp,median , 3, 3); f3=spfilt (gs,median , 3, 3); subplot(246); imshow(f1); title (Median Denoise Salt&Pepper);subplot(247); imshow(f2); title (Median Denoise Pepper);subplot(248); imshow(f3); title (Median Denoise Salt);% 不同Q值的逆谐波均值滤波降噪%figure(2)subplot(241); imshow(f); title(source);subplot(