第四章考试

4.1使用均匀随机数生成具有规定分布的随机数题目：实验代码：imnoise2.m函数：function R = imnoise2(type, M, N, a, b)if nargin == 1a = 0; b = 1;M = 1; N = 1;elseif nargin == 3a = 0; b = 1;endswitch lower(type)case ‘uniform’R = a + (b - a)rand(M, N);case ‘gaussian’R = a + brandn(M, N);case ‘salt & pepper’if nargin <= 3a = 0.05; b = 0.05;endif (a + b) > 1error(‘The sum Pa + Pb must not exceed 1.’)endR(1:M, 1:N) = 0.5;X = rand(M, N);c = find(X <= a);R[表情] = 0;u = a + b;c = find(X > a & X <= u);R[表情] = 1;case ‘lognormal’if nargin <= 3a = 1; b = 0.25;endR = aexp(brandn(M, N));case ‘rayleigh’R = a + (-blog(1 - rand(M, N))).^0.5;case ‘exponential’if nargin <= 3a = 1;endif a <= 0error(‘Parameter a must be positive for exponential type.’)endk = -1/a;R = klog(1 - rand(M, N));case ‘erlang’if nargin <= 3a = 2; b = 5;endif (b ~= round(b) | b <= 0)error(‘Param b must be a positive integer for Erlang.’)endk = -1/a;R = zeros(M, N);for j = 1:bR = R + klog(1 - rand(M, N));endotherwiseerror(‘Unknown distribution type.’)end例4.2 使用函数imnois2所生成数据的直防图题目：代码及结果：代码：r = imnoise2(‘gaussian’,100000,1,0,1);bins = 50;hist(r,bins)title(‘gaussian’)r = imnoise2(‘uniform’,100000,1,0,1);bins = 50;figure,hist(r,bins)title(‘uniform’)r = imnoise2(‘salt & pepper’,1000,1,0.1,0.27);bins = 50;figure,hist(r,bins)title(‘salt & pepper’)r = imnoise2(‘lognormal’,100000,1);bins = 50;figure,hist(r,bins)title(‘lognormal’)r = imnoise2(‘rayleigh’,100000,1,0,1);bins = 50;figure,hist(r,bins)title(‘rayleigh’)r = imnoise2(‘exponential’,100000,1);bins = 50;figure,hist(r,bins)title(‘exponential’)r = imnoise2(‘erlang’,100000,1);bins = 50;figure,hist(r,bins)title(‘erlang’)结果：IMnoise3.mfunction [r, R, S] = imnoise3(M, N, C, A, B)[K, n] = size[表情];if nargin == 3A(1:K) = 1.0;B(1:K, 1:2) = 0;elseif nargin == 4B(1:K, 1:2) = 0;endR = zeros(M, N);for j = 1:Ku1 = M/2 + 1 + C(j, 1); v1 = N/2 + 1 + C(j, 2);R(u1, v1) = i * (A(j)/2) * exp(i2piC(j, 1) * B(j, 1)/M);u2 = M/2 + 1 - C(j, 1); v2 = N/2 + 1 - C(j, 2);R(u2, v2) = -i * (A(j)/2) * exp(i2piC(j, 2) * B(j, 2)/N);endS = abs[表情];r = real(ifft2(ifftshift[表情]));例4.3 使用函数imnoise3代码及实验结果：代码：clcclearC = [0 64; 0 128; 32 32; 64 0; 128 0; -32 32];[r,R,S] = imnoise3(512, 512, C);C1 = C/2;[r1,R1,S1] = imnoise3(512, 512, C1);C2 = [6 32; -2 2];[r2,R2,S2] = imnoise3(512, 512, C2);A = [1 5];[r3,R3,S3] = imnoise3(512, 512, C2, A);imshow(S3,[]) %有两个不清楚的点，因为其振幅较小subplot(3, 2, 1), imshow(S,[]),title(‘[6个]指定冲击的正弦噪声周期频谱[1]’)subplot(3, 2, 2), imshow(r,[]),title(‘[6个]相应的正弦噪声周期模式[1]’)subplot(3, 2, 3), imshow(S1,[]),title(‘[6个]指定冲击的正弦噪声周期频谱[2]’)subplot(3, 2, 4), imshow(r1,[]),title(‘[6个]相应的正弦噪声周期模式[2]’)subplot(3, 2, 5), imshow(r2,[]),title(‘[2个][使用非默认的不同振幅]指定冲击的正弦噪声周期频谱[4]’)subplot(3, 2, 6), imshow(r3,[]),title(‘[2个][使用非默认的不同振幅]相应的正弦噪声周期模式[4]’);结果：创建histroi.m函数function [p, npix] = histroi(f, c, r)B = roipoly(f, c, r);p = imhist(f(B));if nargout > 1npix = sum(B([表情]);end例4.4 估计噪声参数题目代码及运行结果代码：f = imread(‘D:\BaiduNetdiskDownload\Matlab 2020b\bin\数字图像处理\第四章考试\44.tif’);imshow(f)title(‘原始含噪声图像’)[B,c,r] = roipoly(f);figure,imshow(B)[h,npix] = histroi(f,c,r);figure,bar(h,1)title(‘交互式选取区域产生的直方图’)axis tight[v,unv] = statmoments(h,2) ;X=imnoise2(‘gaussian’,npix,1, unv(1), sqrt(unv(2)) );figure,hist(X,130)title(‘使用函数[imnoise2]产生的高斯数据的直方图’)%axis([0 300 0 140])axis tight运行结果：创建spfilt.m函数function f = spfilt(g, type, m, n, parameter)if nargin == 2m = 3; n = 3; Q = 1.5; d = 2;elseif nargin == 5Q = parameter; d = parameter;elseif nargin == 4Q = 1.5; d = 2;elseerror(‘Wrong number of inputs.’);endswitch typecase ‘amean’w = fspecial(‘average’, [m n]);f = imfilter(g, w, ‘replicate’);case ‘gmean’f = gmean(g, m, n);case ‘hmean’f = harmean(g, m, n);case ‘chmean’f = charmean(g, m, n, Q);case ‘median’f = medfilt2(g, [m n], ‘symmetric’);case ‘max’f = ordfilt2(g, mn, ones(m, n), ‘symmetric’);case ‘min’f = ordfilt2(g, 1, ones(m, n), ‘symmetric’);case ‘midpoint’f1 = ordfilt2(g, 1, ones(m, n), ‘symmetric’);f2 = ordfilt2(g, mn, ones(m, n), ‘symmetric’);f = imlincomb(0.5, f1, 0.5, f2);case ‘atrimmed’if (d <= 0) | (d/2 ~= round(d/2))error(‘d must be a positive, even integer.’)endf = alphatrim(g, m, n, d);otherwiseerror(‘Unknown filter type.’)end例4.5使用spfilt函数题目代码及运行结果代码：clcclearf = imread(‘D:\BaiduNetdiskDownload\Matlab 2020b\bin\数字图像处理\第四章考试\45.tif’);imshow(f),title(‘原始图像’)[M,N] = size(f);R = imnoise2(‘salt & pepper’,M,N,0.1,0);c = find(R == 0);gp = f;gp[表情] = 0;imshow(gp),title(‘被概率为0.1的胡椒噪声污染的图像’)R = imnoise2(‘salt & pepper’,M,N,0,0.1);c = find(R == 1);gs = f;gs[表情] = 255;imshow(gs),title(‘被概率为0.1的盐粒噪声污染的图像’)fp = spfilt(gp,‘chmean’,3,3,1.5);imshow(fp),title(‘用阶为Q=1.5的33反调和滤波器对[被概率为0.1的胡椒噪声污染的图像]滤波的结果’)fs = spfilt(gs,‘chmean’,3,3,-1.5);imshow(fs),title(‘用阶为Q=-1.5的33反调和滤波器对[被概率为0.1的盐粒噪声污染的图像]滤波的结果’)fpmax = spfilt(gp,‘max’,3,3);imshow(fpmax),title(‘用33最大滤波器对[被概率为0.1的胡椒噪声污染的图像]滤波的结果’)fsmin = spfilt(gs,‘min’,3,3);imshow(fsmin),title(‘用33最小滤波器对[被概率为0.1的盐粒噪声污染的图像]滤波的结果’)subplot(3, 3, 1), imshow(f),title(‘原始图像’)subplot(3, 3, 2), imshow(gp),title(‘被概率为0.1的胡椒噪声污染的图像’)subplot(3, 3, 3), imshow(gs),title(‘被概率为0.1的盐粒噪声污染的图像’)subplot(3, 3, 4), imshow(fp),title(‘用阶为Q=1.5的33反调和滤波器对[被概率为0.1的胡椒噪声污染的图像]滤波的结果’)subplot(3, 3, 5),imshow(fs),title('用阶为Q=-1.5的33反调和滤波器对[被概率为0.1的盐粒噪声污染的图像]滤波的结果’)subplot(3, 3, 6), imshow(fpmax),title(‘用33最大滤波器对[被概率为0.1的胡椒噪声污染的图像]滤波的结果’)subplot(3, 3, 7), imshow(fsmin),title('用33最小滤波器对[被概率为0.1的盐粒噪声污染的图像]滤波的结果’)运行结果：例4.6自适应中值滤波题目代码及运行结果代码：f = imread(‘D:\BaiduNetdiskDownload\Matlab 2020b\bin\数字图像处理\第四章考试\45.tif’);imshow(f)title(‘原始图像’)g = imnoise(f,‘salt & pepper’,0.25);% 噪声点有黑有白imshow(g)title(‘被概率为0.25椒盐噪声污染的图像’)f1 = medfilt2(g,[7 7],‘symmetric’);imshow(f1)title(‘用77中值滤波器对[被概率为0.25椒盐噪声污染的图像]滤波的结果’)f2 = adpmedian(g,7);imshow(f2)title(‘用Smax=7的自适应中值滤波器对[被概率为0.25椒盐噪声污染的图像]滤波的结果’)subplot(2, 2, 1), imshow(f),title(‘原始图像’)subplot(2, 2, 2), imshow(g),title(‘被概率为0.25椒盐噪声污染的图像’)subplot(2, 2, 3), imshow(f1),title(‘用77中值滤波器对滤波的结果’)subplot(2, 2, 4), imshow(f2),title(‘用Smax=7的自适应中值滤波器滤波的结果’)运行结果：例4.7模糊的带噪图像的建模题目代码及运行结果代码：f = checkerboard(8);imshow(f),title(‘原始图像’)PSF = fspecial(‘motion’,7,45); % sum(PSF([表情]) = 1gb = imfilter(f,PSF,‘circular’);imshow(gb),title(‘使用 PSF = fspecial(motion,7,45) 模糊后的图像’)noise = imnoise(zeros(size(f)),‘gaussian’,0,0.001);imshow(noise,[]),title(‘高斯纯噪声图像’)g = gb + noise;imshow(g,[]),title(‘模糊加噪声的图像’)% imshowMy(pixeldup(f,8),[])subplot(2, 2, 1), imshow(f),title(‘原始图像’)subplot(2, 2, 2), imshow(gb),title(‘使用 PSF = fspecial(motion,7,45) 模糊后的图像’)subplot(2, 2, 3), imshow(noise,[]),title(‘高斯纯噪声图像’)subplot(2, 2, 4), imshow(g,[]),title(‘模糊加噪声的图像’)运行结果：例4.8 使用函数deconvwnr复原模糊带噪图像题目代码及运行结果代码：f = checkerboard(8);PSF = fspecial(‘motion’,7,45);gb = imfilter(f,PSF,‘circular’);noise = imnoise(zeros(size(f)),‘gaussian’,0,0.001);g = gb + noise;imshow(g,[]),title(‘模糊加噪声的图像’)fr1 = deconvwnr(g,PSF);imshow(fr1,[]),title(‘简单的维纳滤波（逆滤波）后的结果’)Sn = abs(fft2(noise)).^2;nA = sum(Sn([表情])/prod(size(noise));Sf = abs(fft2(f)).^2;fA = sum(Sf([表情])/prod(size(f));R = nA/fA;fr2 = deconvwnr(g,PSF,R);imshow(fr2,[]),title(‘使用常数比率的维纳滤波后的结果’)NCORR = fftshift(real(ifft(Sn)));ICORR = fftshift(real(ifft(Sf)));fr3 = deconvwnr(g,PSF,NCORR,ICORR);imshow(fr3,[]),title(‘使用自相关函数的维纳滤波后的结果’)% imshowMy(pixeldup(fr3,8))subplot(2, 2, 1), imshow(g,[]),title(‘模糊加噪声的图像’)subplot(2, 2, 2), imshow(fr1,[]),title(‘简单的维纳滤波（逆滤波）后的结果’)subplot(2, 2, 3), imshow(fr2,[]),title(‘使用常数比率的维纳滤波后的结果’)subplot(2, 2, 4), imshow(fr3,[]),title(‘使用自相关函数的维纳滤波后的结果’)运行结果：例4.9使用radon。题目代码及运行结果代码：g1=zeros(600,600);g1(100:500,250:350)=1;g2=phantom(‘Modified Shepp-Logan’,600);imshow(g1);figure,imshow(g2);theta=0:0.5:179.5;[R1,xp1]=radon(g1,theta);[R2,xp2]=radon(g2,theta);R1=flipud(R1’);R2=flipud(R2’);figure,imshow(R1,[],‘XData’,xp1([1 end]),‘YData’,[179.5 0]);axis xy;axis on;xlabel(‘\rho’),ylabel(‘\theta’);figure,imshow(R2,[],‘XData’,xp2([1 end]),‘YData’,[179.5 0]);axis xy;axis on;xlabel(‘\rho’),ylabel(‘\theta’);subplot(2, 2, 1), imshow(g1)subplot(2, 2, 2), imshow(g2)subplot(2, 2, 3), imshow(R1,[],‘XData’,xp1([1 end]),‘YData’,[179.5 0]);subplot(2, 2, 4), imshow(R2,[],‘XData’,xp2([1 end]),‘YData’,[179.5 0]);运行结果：4.10函数iradon的使用题目;R2=radon(g2, theta);f1=iradon(R1, theta,‘none’);f2=iradon(R2, theta,‘none’);f1_ram=iradon(R1, theta);f2_ram=iradon(R2, theta);f1_hamm=iradon(R1, theta,‘hamming’);f2_hamm=iradon(R2, theta,‘hamming’);f1_near=iradon(R1,theta,‘nearest’);f1_lin=iradon(R1,theta,‘linear’);f1_cub=iradon(R1,theta,‘cubic’);subplot(3, 4, 1), imshow(g1)subplot(3, 4, 2), imshow(g2)subplot(3, 4, 3), imshow(f1,[])subplot(3, 4, 4), imshow(f2,[])subplot(3, 4, 5), imshow (f1_ram,[])subplot(3, 4, 6), imshow (f2_ram,[])subplot(3, 4, 7), imshow (f1_hamm,[])subplot(3, 4, 8), imshow (f2_hamm,[])subplot(3, 4, 9), imshow (f1_near,[])subplot(3, 4, 10), imshow (f1_lin,[])subplot(3, 4, 11), imshow (f1_cub,[])运行结果：例4.11使用fanbeam。题目代码及运行结果代码：g1 =zeros(600,600);g1( 100:500,250:350) = 1;g2=phantom(‘Modified Shepp-Logan’,600);D = 1.5hypot(size(g1,1), size(g1,2))/2;B1_line= fanbeam( g1,D,‘FanSensorGeometry’,‘line’,…‘FanSensorSpacing’ ,1 ,‘FanRotationIncrement’ ,0.5);B1_line= flipud(B1_line’);B2_line= fanbeam( g2,D,‘FanSensorGeometry’,‘line’,…‘FanSensorSpacing’ ,1 ,‘FanRotationIncrement’ ,0.5);B2_line= flipud(B2_line’);imshow(B1_line,[]);figure, imshow(B2_line,[]);B1_arc= fanbeam( g1,D,‘FanSensorGeometry’,‘arc’,…‘FanSensorSpacing’ ,.08,‘FanRotationIncrement’ ,0.5);B2_arc=fanbeam( g2,D,‘FanSensorGeometry’,‘arc’,…‘FanSensorSpacing’ ,.08,‘FanRotationIncrement’ ,0.5);figure,imshow( flipud(B1_arc’),[]);figure,imshow( flipud(B2_arc’),[]);subplot(2, 2, 1), imshow(B1_line,[]);subplot(2, 2, 2), imshow(B2_line,[]);subplot(2, 2, 3), imshow( flipud(B1_arc’),[]);subplot(2, 2, 4), imshow( flipud(B2_arc’),[]);运行结果：例4.12使用函数ifanbeam。题目代码及运行结果：代码：g=phantom(‘Modified Shepp-Logan’,600);D = 1.5hypot(size(g,1), size(g,2))/2;B1 = fanbeam(g,D);f1 =ifanbeam(B1,D);figure, imshow(f1,[])B2 = fanbeam(g,D,‘FanRotationIncrement’,0.5,…‘FanSensorSpacing’,0.5);f2 = ifanbeam(B2,D,‘FanRotationIncrement’,0.5,…‘FanSensorSpacing’,0.5,‘Filter’,‘Hamming’);figure, imshow(f2,[])B3 = fanbeam(g,D,‘FanRotationIncrement’,0.5,…‘FanSensorSpacing’,0.05) ;f3 = ifanbeam(B3,D,‘FanRotationIncrement’,0.5,…?‘FanSensorSpacing’,0.05,‘Filter’,‘Hamming’);figure, imshow(f3,[])subplot(1, 3, 1), imshow(f1,[])subplot(1, 3, 2), imshow(f2,[])subplot(1, 3, 3), imshow(f3,[])运行结果：创建函数adpmedianfunction f = adpmedian(g, Smax)%ADPMEDIAN Perform adaptive median filtering.% F = ADPMEDIAN(G, SMAX) performs adaptive median filtering of% image G. The median filter starts at size 3-by-3 and iterates up% to size SMAX-by-SMAX. SMAX must be an odd integer greater than 1.% Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins% Digital Image Processing Using MATLAB, Prentice-Hall, 2004% $Revision: 1.5 $ $Date: 2003/11/21 14:19:05 $% SMAX must be an odd, positive integer greater than 1.if (Smax <= 1) | (Smax/2 == round(Smax/2)) | (Smax ~= round(Smax))error(‘SMAX must be an odd integer > 1.’)end[M, N] = size(g);% Initial setup.f = g;f([表情] = 0;alreadyProcessed = false(size(g));% Begin filtering.for k = 3:2:Smaxzmin = ordfilt2(g, 1, ones(k, k), ‘symmetric’);zmax = ordfilt2(g, k * k, ones(k, k), ‘symmetric’);zmed = medfilt2(g, [k k], ‘symmetric’);processUsingLevelB = (zmed > zmin) & (zmax > zmed) & …~alreadyProcessed;zB = (g > zmin) & (zmax > g);outputZxy = processUsingLevelB & zB;outputZmed = processUsingLevelB & ~zB;f(outputZxy) = g(outputZxy);f(outputZmed) = zmed(outputZmed);alreadyProcessed = alreadyProcessed | processUsingLevelB;if all(alreadyProcessed([表情])break;endend% Output zmed for any remaining unprocessed pixels. Note that this% zmed was computed using a window of size Smax-by-Smax, which is% the final value of k in the loop.f(~alreadyProcessed) = zmed(~alreadyProcessed);创建函数changeclassfunction image = changeclass(class, varargin)switch classcase ‘uint8’image = im2uint8(varargin{:});case ‘uint16’image = im2uint16(varargin{:});case ‘double’image = im2double(varargin{:});otherwiseerror(‘Unsupported IPT data class.’);end创建函数statmomentsfunction [v, unv] = statmoments(p, n)Lp = length§;if (Lp ~= 256) & (Lp ~= 65536)error(‘P must be a 256- or 65536-element vector.’);endG = Lp - 1;p = p/sum§; p = p([表情];z = 0:G;z = z./G;m = zp;z = z - m;v = zeros(1, n);v(1) = m;for j = 2:nv(j) = (z.^j)*p;endif nargout > 1unv = zeros(1, n);unv(1)=m.G;for j = 2:nunv(j) = ((zG).^j)p;endend例4.13使用fan2pare题目代码及运行结果代码：g1 =zeros(600,600);g1( 100:500,250:350) = 1 ;g2=phantom(‘Modified Shepp-Logan’,600);D =1.5hypot(size(g1,1), size(g1,2))/2;B1_line= fanbeam( g1,D,‘FanSensorGeometry’,…‘line’,‘FanSensorSpacing’,1,…‘FanRotationIncrement’ ,0.5);B2_arc= fanbeam( g2,D,‘FanSensorGeometry’,‘arc’,…‘FanSensorSpacing’,.08,‘FanRotationIncrement’,0.5);P1_line=fan2para(B1_line,D,‘FanRotationIncrement’,0.5,…‘FanSensorGeometry’,‘line’ ,…‘FanSensorSpacing’ ,1,…‘ParallelCoverage’,‘halfcycle’,…‘ParallelRotationIncrement’,0.5,…‘ParallelSensorSpacing’,1);P2_arc=fan2para(B2_arc,D,‘FanRotationIncrement’,0.5,…‘FanSensorGeometry’,‘arc’,…‘FanSensorSpacing’ ,0.08,…‘ParallelCoverage’,‘halfcycle’,…‘ParallelRotationIncrement’,0.5,…‘ParallelSensorSpacing’,1);P1_line= flipud(P1_line’);P2_arc= flipud(P2_arc’);figure, imshow(P1_line,[])figure, imshow(P2_arc,[])subplot(2, 1, 1), imshow(P1_line,[])subplot(2, 1, 2), imshow(P2_arc,[])