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

Question

我想在 Matlab 中实现“自适应分水岭分割”。该算法有六个步骤。输入是图（a），结果是图（d）。请您帮我检查一下我的代码是否有错误，我不知道如何实现第六步。太感谢了！

加载图像：

input_image = imread('test.gif');

Step 1：计算每个(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);

Step 2：使用高斯滤波器对距离图进行平滑，合并相邻的最大值，如果D(x,y)是局部最大值，则将M(x,y)设为1，然后得到距离图的标记图。

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

Step3：逐像素扫描标记图。如果 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

第4步：

计算距离图的倒数，局部极大值变成局部极小值。

igfDT = -(gfDT);

第 5 步：

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

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