本文系《数字图像处理原理与实践（MATLAB版）》一书之代码系列的Part7（由于之前发布顺序调整，请读者注意页码标注而不要仅仅依据系列文章的标题编号），辑录该书第186至第225页之代码，供有需要读者下载研究使用。至此全书代码发布已经过半。代码执行结果请参见原书配图，建议下载代码前阅读下文：
关于《数字图像处理原理与实践（MATLAB版）》一书代码发布的说明
http://blog.csdn.net/baimafujinji/article/details/40987807
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 on
s=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);
end
end
P197
I = im2bw(imread(‘penguins.bmp‘), 0.38);
BW = 1-I;
B = bwboundaries(BW,8,‘noholes‘);
imshow(I)
hold on
for 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 on
for k = 1:length(B)
boundary = B{k};
plot(boundary(:,2), boundary(:,1), ‘g‘, ‘LineWidth‘, 2)
end
saveas(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; end
J = 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);
end
end
blocks(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);
figure
subplot(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);
figure
subplot(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‘);
figure
subplot(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);
（代码发布未完，请待后续...）