algorithm - 在 Matlab 中实现自适应分水岭分割

Question

我想在 Matlab 中实现“自适应分水岭分割”。这个算法有六个步骤。输入是图(a)，结果是图(d)。请你帮我检查一下我的代码有没有错误，我不知道如何实现第六步。非常感谢!

加载图片:

input_image = imread('test.gif');

第一步:计算每个(x,y)处的D(x,y)，得到二值图像的欧几里得距离图，并将M(x,y)的每个值赋为0。

DT = bwdist(input_image,'euclidean'); % Trandform distance:Euclidian distance 
[h,w]=size(DT);
M = zeros(h,w);

第二步:使用高斯滤波器对距离图进行平滑处理，合并相邻的极大值，如果D(x,y)为局部极大值，则将M(x,y)设为1，得到距离图的标记图.

H = fspecial('gaussian');
gfDT = imfilter(DT,H); 
M = imregionalmax(gfDT); % maker map, M = local maximum of gfDT

第三步:逐像素扫描标记图。如果 M(x0,y0) 为 1，则在其半径为 D(x ,y ) 的邻域中寻找伪最大值。当 M(x,y) 等于 1 且 sqr((x − x0)^2 + (y − y0)^2 ) ≤ D(x0, y0) , 如果 D(x,y) < D(x0,y0) 则设置 M(x,y) 为 0。

for x0 = 1:h
    for y0 = 1:w
        if M(x0,y0) == 1
            r = ceil(gfDT(x0,y0));

            % range begin:(x0-r,y0-r) end:(x0+r,y0+r)
            xb = x0-r;
            if xb <= 0
                xb =1;
            end

            yb = y0-r;
            if yb <= 0
                yb =1;
            end

            xe = x0+r;
            if xe > w
            xe = w;
            end

            ye = y0+r;
            if ye > h
                ye = h;
            end

            for x = yb:ye
                for y = xb:xe
                    if M(x,y)==1 
                        Pos = [x0,y0 ;x,y];
                        Dis = pdist(Pos,'euclidean');
                        IFA = Dis<= (gfDT(x0,y0));
                        IFB = gfDT(x,y)<gfDT(x0,y0);
                        if ( IFA && IFB)
                            M(x,y) = 0;
                        end
                    end
                end
            end
        end
    end
end

第四步:

计算距离图的倒数，局部最大值就是局部最小值。

igfDT = -(gfDT);

第五步:

利用常规分水岭算法根据标记对距离图进行分割，得到二值图像的分割结果。

I2 = imimposemin(igfDT,M);
L = watershed(I2);
igfDT (L==0)=0;

第 6 步:通过用直线连接分水岭线的末端并沿直线重新分类像素来拉直分水岭线。

我不知道如何实现这一步

Answer 1

先尝试距离变换，然后再进行分水岭变换。

im=imread('n6BRI.gif');

imb=bwdist(im);

sigma=3;
kernel = fspecial('gaussian',4*sigma+1,sigma);
im2=imfilter(imb,kernel,'symmetric');

L = watershed(max(im2(:))-im2);
[x,y]=find(L==0);

lblImg = bwlabel(L&~im);

figure,imshow(label2rgb(lblImg,'jet','k','shuffle'));