图像处理（MATLAB版）第四章----4.1

4.1使用均匀随机数生成具有规定分布的随机数

题目：

实验代码：

imnoise2.m函数：
function R = imnoise2(type, M, N, a, b)
if nargin == 1
a = 0; b = 1;
M = 1; N = 1;
elseif nargin == 3
a = 0; b = 1;
end
switch lower(type)
case ‘uniform’
R = a + (b - a)rand(M, N);
case ‘gaussian’
R = a + brandn(M, N);
case ‘salt & pepper’
if nargin <= 3
a = 0.05; b = 0.05;
end

if (a + b) > 1
error(‘The sum Pa + Pb must not exceed 1.’)
end
R(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 <= 3
a = 1; b = 0.25;
end
R = aexp(brandn(M, N));
case ‘rayleigh’
R = a + (-blog(1 - rand(M, N))).^0.5;
case ‘exponential’
if nargin <= 3
a = 1;
end
if a <= 0
error(‘Parameter a must be positive for exponential type.’)
end
k = -1/a;
R = klog(1 - rand(M, N));
case ‘erlang’
if nargin <= 3
a = 2; b = 5;
end
if (b ~= round(b) | b <= 0)
error(‘Param b must be a positive integer for Erlang.’)
end
k = -1/a;
R = zeros(M, N);
for j = 1:b
R = R + k*log(1 - rand(M, N));
end
otherwise
error(‘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.m
function [r, R, S] = imnoise3(M, N, C, A, B)
[K, n] = size©;
if nargin == 3
A(1:K) = 1.0;
B(1:K, 1:2) = 0;
elseif nargin == 4
B(1:K, 1:2) = 0;
end

R = zeros(M, N);
for j = 1:K
u1 = 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);
end

S = abs®;
r = real(ifft2(ifftshift®));

例4.3 使用函数imnoise3

代码及实验结果：
代码：
clc
clear
C = [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 > 1
npix = 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 == 2
m = 3; n = 3; Q = 1.5; d = 2;
elseif nargin == 5
Q = parameter; d = parameter;
elseif nargin == 4
Q = 1.5; d = 2;
else
error(‘Wrong number of inputs.’);
end
switch type
case ‘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.’)
end
f = alphatrim(g, m, n, d);
otherwise
error(‘Unknown filter type.’)
end

例4.5使用spfilt函数

题目


代码及运行结果
代码：
clc
clear
f = 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的3*3反调和滤波器对[被概率为0.1的胡椒噪声污染的图像]滤波的结果’)

fs = spfilt(gs,‘chmean’,3,3,-1.5);
imshow(fs),title(‘用阶为Q=-1.5的3*3反调和滤波器对[被概率为0.1的盐粒噪声污染的图像]滤波的结果’)

fpmax = spfilt(gp,‘max’,3,3);
imshow(fpmax),title(‘用3*3最大滤波器对[被概率为0.1的胡椒噪声污染的图像]滤波的结果’)

fsmin = spfilt(gs,‘min’,3,3);
imshow(fsmin),title(‘用3*3最小滤波器对[被概率为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(‘用7*7中值滤波器对[被概率为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(‘用7*7中值滤波器对滤波的结果’)
subplot(2, 2, 4), imshow(f2),title(‘用Smax=7的自适应中值滤波器滤波的结果’)

运行结果：

例4.7模糊的带噪图像的建模

题目


代码及运行结果
代码：
f = checkerboard(8);
imshow(f),title(‘原始图像’)

PSF = fspecial(‘motion’,7,45); % sum(PSF(😃) = 1
gb = 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的使用

题目

![在这里插入图片描述](https://img-blog.csdnimg.cn/8c71889beb5c48a2a196d59fccdda4fb.png
代码及运行结果
代码：
theta =0:0.5:179.5;
R1=radon(g1, theta);
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.5*hypot(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.5*hypot(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,[])

运行结果：

创建函数adpmedian

function 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:Smax
zmin = 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;
end
end

% 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);

创建函数changeclass

function image = changeclass(class, varargin)
switch class
case ‘uint8’
image = im2uint8(varargin{:});
case ‘uint16’
image = im2uint16(varargin{:});
case ‘double’
image = im2double(varargin{:});
otherwise
error(‘Unsupported IPT data class.’);
end

创建函数statmoments

function [v, unv] = statmoments(p, n)
Lp = length§;
if (Lp ~= 256) & (Lp ~= 65536)
error(‘P must be a 256- or 65536-element vector.’);
end
G = Lp - 1;
p = p/sum§; p = p(😃;
z = 0:G;
z = z./G;
m = z*p;
z = z - m;
v = zeros(1, n);
v(1) = m;
for j = 2:n
v(j) = (z.^j)*p;
end

if nargout > 1
unv = zeros(1, n);
unv(1)=m.G;
for j = 2:n
unv(j) = ((zG).^j)*p;
end
end

例4.13使用fan2pare

题目

代码及运行结果
代码：
g1 =zeros(600,600);
g1( 100:500,250:350) = 1 ;
g2=phantom(‘Modified Shepp-Logan’,600);
D =1.5*hypot(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,[])

运行结果：