首页 文章

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

提问于
浏览
2

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

Input image

Result image

加载图片:

input_image = imread('test.gif');

步骤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);

步骤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

步骤3:逐个像素地扫描标记图 . 如果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);

STEP5:

通过传统的分水岭算法根据标记对距离图进行分割,得到二值图像的分割 .

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

步骤6:通过将分水线的末端与直线连接并沿着直线重新分类像素来使分水岭线变直 .

I don't know how to implement this step

1 回答

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

    enter image description here

相关问题