主动轮廓法的算法实现

主动轮廓法于1987年被提出，现在轮廓提取、对象跟踪等许多方面得到了广泛的应用。本资源实现了主动轮廓法的算法。

function DeformSnake
%轮廓变形
I = imread（‘imc.bmp‘）;

lstCtrlPt = InitSnake; %控制点

imshow（I）;
hold on
plot（lstCtrlPt（1:）lstCtrlPt（2:）‘go‘）;

[noUse NCtrlPt] = size（lstCtrlPt）;  %NCtrlPt:轮廓控制点数

[imHeight imWidth noUse] = size（I）;

param = [0.30.10.5]‘; %[alphabetagamma]‘

THRESHOLD1 = 100;    %梯度平方最小差距
THRESHOLD2 = 0.05;   %此次移动点数点总控制点数最小比例
rPtMoved = 1;       %此次移动点数点总控制点数比例

I = I（::1）;   %转换成8位灰度图
I= double（I）;
I = gaussMask（I）;   %高斯模糊
Igrad = edge_detect（I）;%求图像梯度
nCount = 0;
while rPtMoved >= THRESHOLD2
%%****************初始化轮廓****************%%
    lstCtrlPt = snakeInterp（lstCtrlPt）;
    [noUse NCtrlPt] = size（lstCtrlPt）;  %NCtrlPt:轮廓控制点数  
%%************************轮廓变形*************************%%
    nPtMoved = 0;
    i = 1;
    while i <= NCtrlPt
        
        gradIadj = Igrad（lstCtrlPt（2i）-2:lstCtrlPt（2i）+2lstCtrlPt（1i）-2:lstCtrlPt（1i）+2）;
        
        Max = max（gradIadj.^2）;
        Min = min（gradIadj.^2）;
        %保证 Max-Min >= THRESHOLD1
        if Max（1）-Min（1） < THRESHOLD1
            Min（1） = Max（1）-THRESHOLD1;
        end
        %dist:相邻两点的距离和davg:平均距离
        dist = sqrt（sum（（[lstCtrlPt（:2:NCtrlPt） lstCtrlPt（:1）]-lstCtrlPt）.^2））;
        davg = sum（dist）/NCtrlPt;
        
        %计算邻域内每个点的能量
        for yOffset = -2:2
            for xOffset = -2:2
                %当前邻域中的点:pt
                pt = [lstCtrlPt（1i）+xOffsetlstCtrlPt（2i）+yOffset]‘;
                subList = [lstCtrlPt（:pre（iNCtrlPt）） pt lstCtrlPt（:next（iNCtrlPt））];
                dist = sqrt（sum（（subList（:2:3）-subList（:1:2））.^2））;
                %计算弹性能量：Eelastic
                Eelastic = （davg-dist（1））^2;

                %计算弯曲能量：Ebending
                Ebending = sum（（subList（:1）-2*pt+subList（:3））.^2）;

                %计算图像能量：Eimage
                Eimage = （Min（1）-gradIadj（3+yOffset3+xOffset）^2）/（Max（1）-Min（1））;

                %将三个量保存起来:Energy
                Energy（（3+yOffset-1）*5+3+xOffset:） = [EelasticEbendingEimage];
            end
        end

        MaxEelastic = max（Energy（:1））;
        MaxEbending = max（Energy（:2））;

        %单位化能量：
        Energy（:1） = Energy（:1）/MaxEelastic;
        Energy（:2） = Energy（:2）/MaxEbending;

        Etotal = Energy*param;
        
        %找能量最小点的坐标
        index = find（Etotal == min（Etotal））;
        yOffset = ceil（index（1）/5）-3;
        xOffset = index（1） -（2+yOffset）*5-3;
        ptNewCtrl = [lstCtrlPt（1i）+xOffsetlstCtrlPt（2i）+yOffset]‘;
        
        %更新控制点
        if sum（ptNewCtrl ~= lstCtrlPt（:i））
            nPtMoved = nPtMoved+1;
            lstCtrlPt（:i）=ptNewCtrl;
        end
        i = i+1;
    end
    nCount = nCount+1;
    rPtMoved = nPtMoved/NCtrlPt;
end
plot（lstCtrlPt（1:）lstCtrlPt（2:）‘b*‘）;
%显示参数
legend（‘origional CtrlPt‘‘final CtrlPt‘）;
paramStr{11} = [‘\alpha = ‘ num2str（param（11））];
paramStr{21} = [‘\beta = ‘ num2str（param（21））];
paramStr{31} =

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件        646  2008-08-13 21:30  运行方式.txt

     文件        714  2008-08-15 20:04  matlab\InitSnake.m

     文件     463816  2008-05-12 20:46  matlab\ima.bmp

     文件        284  2008-05-19 10:12  matlab\edge_detect.m

     文件        231  2008-08-15 20:03  matlab\isInRegion.m

     文件          0  2008-05-16 14:39  matlab\ToGrad.m

     文件        499  2008-05-17 09:23  matlab\insertPt.m

     文件     265198  2008-05-17 16:10  matlab\1_1_1p5.bmp

     文件        120  2008-05-18 10:03  matlab\pre.m

     文件        126  2008-05-18 10:01  matlab\next.m

     文件     156120  2008-05-18 10:41  matlab\imb.bmp

     文件       3338  2008-05-19 16:37  matlab\DeformSnake.m

     文件        299  2008-05-19 08:59  matlab\gaussMask.m

     文件     265198  2008-05-19 10:09  matlab\ima_1p5_1_3.bmp

     文件     265198  2008-05-19 10:10  matlab\imb_1p5_1_3.bmp

     文件     156120  2008-05-19 16:08  matlab\imc.bmp

     文件     265198  2008-05-19 16:13  matlab\Snake不能收敛到凹陷示例.bmp

     文件        810  2008-05-20 17:20  matlab\snakeInterp.m

    ..A.SH.        98  2008-06-02 22:59  matlab\desktop.ini

    ...D..R         0  2008-08-15 19:58  matlab

----------- ---------  ---------- -----  ----

              1844013                    20