本文系《数字图像处理原理与实践（MATLAB版）》一书之代码系列的Part4（由于之前发布顺序调整，请读者注意页码标注而不要仅仅依据系列文章的标题编号），辑录该书第186至第280页之代码，供有需要读者下载研究使用。至此全书代码发布已经过半。代码执行结果请参见原书配图。
-------------------------------------------
P186
A = rgb2gray(imread('circle.png'));B = edge(A, 'canny');
[centers, radii, metric] = imfindcircles(B,[22 65]);
imshow(A);
viscircles(centers, radii,'EdgeColor','b');
P195
BW = imread('contour.bmp');imshow(BW,[]);hold ons=size(BW);for row = 2:55:s(1) for col=1:s(2) if BW(row,col), break; end end contour = bwtraceboundary(BW, [row, col], 'W', 8); if(~isempty(contour)) plot(contour(:,2),contour(:,1),'g','LineWidth',2); endend
P197
I = im2bw(imread('penguins.bmp'), 0.38);BW = 1-I;B = bwboundaries(BW,8,'noholes');imshow(I)hold onfor k = 1:length(B) boundary = B{k}; plot(boundary(:,2), boundary(:,1), 'g', 'LineWidth', 2)end
-------------------------------------------
补充一点小技巧。先前在写Demo的时候没想过这个问题，后来是因为要为图书做插图，所以就需要把处理结果的白边去掉，下面这段代码实现了这种结果。与图像处理无关，这种方法也没有出现在书里，只是关于MATLAB保存图像时的一点小技巧而已。
I = im2bw(imread('penguins.bmp'), 0.38);BW = 1-I;B = bwboundaries(BW,8,'noholes');imshow(I,'border','tight');hold onfor k = 1:length(B) boundary = B{k}; plot(boundary(:,2), boundary(:,1), 'g', 'LineWidth', 2)endsaveas(gcf,'pengs3.bmp');
-------------------------------------------
P203
I = imread('nums.bmp');locs =[57 64;47 103;81 224;94 274;11 365;85 461;86 540];BW = imfill(I, locs, 4);imshow(BW);
P204
I = imread('nums.bmp');BW2 = imfill(I,'holes');imshow(BW2);
P205
I = imread('tire.tif');I2 = imfill(I);figure, imshow(I), figure, imshow(I2)
P206
I = imread('eight.tif');c = [222 272 300 270 221 194];r = [21 21 75 121 121 75];J = roifill(I,c,r);imshow(I)figure, imshow(J)
P207
function J = regiongrowing(I,x,y,threshold)if(exist('threshold','var')==0), threshold=0.2; endJ = zeros(size(I)); % 用来标记输出结果的二值矩阵[m n] = size(I); % 输入图像的尺寸reg_mean = I(x,y); % 被分割区域的灰度均值reg_size = 1; % 区域中像素的数目% 用以存储被分割出来的区域的邻域点的堆栈neg_free = 10000; neg_pos=0;neg_list = zeros(neg_free,3);delta=0; % 最新被引入的像素与区域灰度均值的差值% 区域生长直至满足终止条件while(delta<threshold && reg_size<numel(I)) % 检测邻域像素，并判读是否将其划入区域 for i = -1:1 for j = -1:1 xn = x + i; yn = y + j; % 计算邻域点的坐标 % 检查邻域像素是否越界 indicator = (xn >= 1)&&(yn >= 1)&&(xn <= m)&&(yn <= n); % 如果邻域像素还不属于被分割区域则加入堆栈 if(indicator && (J(xn,yn)==0)) neg_pos = neg_pos+1; neg_list(neg_pos,:) = [xn yn I(xn,yn)]; J(xn,yn)=1; end end end if(neg_pos+10>neg_free), % 如果堆栈空间不足，则对其进行扩容 neg_free=neg_free+10000; neg_list((neg_pos+1):neg_free,:)=0; end % 将那些灰度值最接近区域均值的像素加入到区域中去 dist = abs(neg_list(1:neg_pos,3)-reg_mean); [delta, index] = min(dist); J(x,y)=2; reg_size=reg_size+1; % 计算新区域的均值 reg_mean = (reg_mean*reg_size + neg_list(index,3))/(reg_size+1); % 保存像素坐标，然后将像素从堆栈中移除 x = neg_list(index,1); y = neg_list(index,2); neg_list(index,:)=neg_list(neg_pos,:); neg_pos=neg_pos-1;end% 将由区域生长得到的分割区域以二值矩阵的形式返回J=J>1;
P208
I = im2double(rgb2gray(imread('penguins.bmp')));x = 244; y = 679;J = regiongrowing(I,x,y,0.2);figure, imshow(I+J);
P213
I = imread('liftingbody.png');S = qtdecomp(I,.27);blocks = repmat(uint8(0),size(S));for dim = [512 256 128 64 32 16 8 4 2 1]; numblocks = length(find(S==dim)); if (numblocks > 0) values = repmat(uint8(1),[dim dim numblocks]); values(2:dim,2:dim,:) = 0; blocks = qtsetblk(blocks,S,dim,values); endendblocks(end,1:end) = 1;blocks(1:end,end) = 1;imshow(I), figure, imshow(blocks,[])
P219
rgb = imread('potatos.jpg');I = rgb2gray(rgb);hy = fspecial('sobel');hx = hy';Iy = imfilter(double(I), hy, 'replicate');Ix = imfilter(double(I), hx, 'replicate');gradmag = sqrt(Ix.^2 + Iy.^2);L = watershed(gradmag);Lrgb = label2rgb(L);figuresubplot(1, 2, 1); imshow(gradmag,[]), title('梯度幅值图像')subplot(1, 2, 2); imshow(Lrgb); title('对梯度图直接做分水岭分割')
P221-P224
rgb = imread('potatos.jpg');I = rgb2gray(rgb);hy = fspecial('sobel');hx = hy';Iy = imfilter(double(I), hy, 'replicate');Ix = imfilter(double(I), hx, 'replicate');gradmag = sqrt(Ix.^2 + Iy.^2);se = strel('disk', 12);Ie = imerode(I, se);Iobr = imreconstruct(Ie, I);Iobrd = imdilate(Iobr, se);Iobrcbr = imreconstruct(imcomplement(Iobrd), imcomplement(Iobr));Iobrcbr = imcomplement(Iobrcbr);fgm = imregionalmax(Iobrcbr);It1 = rgb(:, :, 1);It2 = rgb(:, :, 2);It3 = rgb(:, :, 3);It1(fgm) = 255; It2(fgm) = 0; It3(fgm) = 0;I2 = cat(3, It1, It2, It3);figuresubplot(1, 2, 1); imshow(fgm, []); title('局部极大值图像');subplot(1, 2, 2); imshow(I2); title('局部极大值叠加图像');se2 = strel(ones(15,15));fgm2 = imclose(fgm, se2);fgm3 = imerode(fgm2, se2);fgm4 = bwareaopen(fgm3, 400);bw = im2bw(Iobrcbr, graythresh(Iobrcbr));D = bwdist(bw);DL = watershed(D);bgm = DL == 0;gradmag2 = imimposemin(gradmag, bgm | fgm4);L = watershed(gradmag2);%第一种显示方法Lrgb = label2rgb(L, 'jet', 'w', 'shuffle');figuresubplot(1,2,1), imshow(Lrgb), title('分水岭分割结果显示1');%第二种显示方法subplot(1, 2, 2); imshow(rgb, []), title('分水岭分割结果显示2');hold on;himage = imshow(Lrgb);set(himage, 'AlphaData', 0.3);
P245
I = imread('lena.png');fcoef=fft2(double(I)); %FFT变换tmp1 =log(1+abs(fcoef));spectrum = fftshift(fcoef); %调整中心tmp2 = log(1+abs(spectrum));ifcoef = ifft2(fcoef); %逆变换figure %显示处理结果subplot(2,2,1), imshow(I), title('source image');subplot(2,2,2), imshow(tmp1,[]), title('FFT image');subplot(2,2,3), imshow(tmp2,[]), title('shift FFT image');subplot(2,2,4), imshow(ifcoef,[]), title('IFFT image');
P251
J= double(imread('lena.bmp'));K = dct2(J);figure, imshow(K,[0 255])
P252-1
J= double(imread('lena.bmp'));K = dct2(J);figure, imshow(K,[0 255]);K_i = idct2(K);figure, imshow(K_i,[0 255])
P252-2
J= double(imread('lena.bmp'));A = J(1:8,1:8);D = dctmtx(8);dct_1 = D*A;dct_2 = D'*dct_1;
P252-3
J= double(imread('lena.bmp'));A = J(1:8,1:8);D = dctmtx(8);dct_1 = D*A*D';dct_2 = dct2(A);
P253
I = imread('cameraman.tif');I = im2double(I);T = dctmtx(8);dct = @(block_struct) T * block_struct.data * T';B = blockproc(I,[8 8],dct);mask = [1 1 1 1 0 0 0 0 1 1 1 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];B2 = blockproc(B,[8 8],@(block_struct) mask .* block_struct.data);invdct = @(block_struct) T' * block_struct.data * T;I2 = blockproc(B2,[8 8],invdct);imshow(I), figure, imshow(I2)
P262
a = [0 0 1 1 0 0 1 1];b = fwht(a);
P263
I = imread('baboon.bmp');I1 = double(I);T = hadamard(8);myFun1 = @(block_struct)T*block_struct.data*T/64;H = blockproc(I1, [8 8], myFun1);H(abs(H)<3.5)=0;myFun2 = @(block_struct)T*block_struct.data*T;I2 = blockproc(H, [8 8], myFun2);subplot(121), imshow(I1,[]), title('original image');subplot(122), imshow(I2,[]), title('zipped image');
P264
I = imread('baboon.bmp');I1 = double(I);[m n] =size(I);sizi = 8;num = 16;%分块进行离散沃尔什变换T = hadamard(sizi);myFun1 = @(block_struct)T*block_struct.data*T/(sizi.^2);hdcoe = blockproc(I1, [sizi, sizi], myFun1);%重新排列系数coe = im2col(hdcoe, [sizi, sizi], 'distinct');coe_t = abs(coe);[Y, ind] = sort(coe_t);%舍去绝对值较小的系数[m_c, n_c] = size(coe);for i = 1:n_ccoe(ind(1:num, i), i)=0;end%重建图像re_hdcoe = col2im(coe, [sizi, sizi], [m, n], 'distinct');myFun2 = @(block_struct)T*block_struct.data*T;re_s = blockproc(re_hdcoe, [sizi, sizi], myFun2);subplot(121), imshow(I1,[]), title('original image');subplot(122), imshow(re_s,[]), title('compressed image');
P268
dim1 = [1 1 1 2 2 2 3 3 3];dim2 = [1 2 3 1 2 3 1 2 3];dim3 = [63 75 78 50 56 65 70 71 80];sum( (dim1-mean(dim1)) .* (dim2-mean(dim2)) ) / ( 9-1 ) % 0sum( (dim1-mean(dim1)) .* (dim3-mean(dim3)) ) / ( 9-1 ) % 0.625sum( (dim2-mean(dim2)) .* (dim3-mean(dim3)) ) / (9-1 ) % 5std(dim1)^2 % 0.75std(dim2)^2 % 0.75std(dim3)^2 % 100.7778
P274
X = [2 2; 2 3; 3 4; 4 3; 5 4; 5 5];[COEFF,SCORE,latent,tsquare] = princomp(X);
P275-1
X0=X-repmat(mean(X),6,1);SCORE_1 = X0*COEFF;
P275-2
X = [2 2; 2 3; 3 4; 4 3; 5 4; 5 5];V = cov(X);[COEFF,latent] = pcacov(V)
P277
I = imread('baboon.bmp');x = double(I)/255;[m,n]=size(x);y =[];%拆解图像for i = 1:m/8; for j = 1:n/8; ii = (i-1)*8+1; jj = (j-1)*8+1; y_app = reshape(x(ii:ii+7,jj:jj+7),1,64); y=[y;y_app]; endend%KL变换[COEFF,SCORE,latent] = princomp(y);kl = y * COEFF;kl1 = kl;kl2 = kl;kl3 = kl;%置零压缩过程kl1(:, 33:64)=0;kl2(:, 17:64)=0;kl3(:, 9:64)=0;%KL逆变换kl_i = kl*COEFF';kl1_i = kl1*COEFF';kl2_i = kl2*COEFF';kl3_i = kl3*COEFF';image = ones(256,256);image1 = ones(256,256);image2 = ones(256,256);image3 = ones(256,256);k=1;%重组图像for i = 1:m/8; for j = 1:n/8; y = reshape(kl_i(k, 1:64),8,8); y1 = reshape(kl1_i(k, 1:64),8,8); y2 = reshape(kl2_i(k, 1:64),8,8); y3 = reshape(kl3_i(k, 1:64),8,8); ii = (i-1)*8+1; jj = (j-1)*8+1; image(ii:ii+7,jj:jj+7) = y; image1(ii:ii+7,jj:jj+7) = y1; image2(ii:ii+7,jj:jj+7) = y2; image3(ii:ii+7,jj:jj+7) = y3; k=k+1; endend
（代码发布未完，请待后续...）
