Code: greedy_ice/greedyice

function [con, ptsmoved] = greedyice_movecontour(con, im, alpha, beta, gamma, delta);
% GREEDYICE_MOVECONTOUR Moves contour points to new local optima

%% constants
nr = 2; %% radius of square enclosing neighbouring points

%% init
ptsmoved = 0; % total points moved
n = length(con); % number of points in contour
imsz = size(im); % size of image matrix
mdist = meandist(con); % mean distance between points on contour
mpoint = sum(con) ./ length(con); % center of the contour

%% loop to move points to new locations
vlast = con(n,:); %% initial vi-1
for i = 1:n+1

    %% point 0 is the first and last
    if i > n, i = 1; end
    
    %% compute vectors needed in energy calc
    vi = con(i,:);
    if i == n, vnext = con(1,:);
    else vnext = con(i+1,:); end

    %% compute neighbours of vi
    neighbours = [];
    j = 1;
    for y = vi(1,2)-nr:vi(1,2)+nr
        for x = vi(1,1)-nr:vi(1,1)+nr
            %% check is inside image
            if x > 0 && y > 0 && x <= imsz(1) && y <= imsz(2)
                neighbours(j,1) = x;
                neighbours(j,2) = y;
                j = j + 1;
            end
        end
    end

    %% compute term values for energy functionals of neighbours
    econt = []; ecurv = []; eimage = []; ecent = [];
    econtmax = 0; ecurvmax = -inf;
    ecentmax = 0; ecentmin = inf;
    eimagemax = -inf; eimagemin = inf;
    for j = 1:length(neighbours)
        vj = neighbours(j,:);
        %% calc econt
        econt(j) = norm(mdist - norm(vj - vlast));
        if econt(j) > econtmax, econtmax = econt(j); end
        %% calc ecurv
        ecurv(j) = norm(vlast - (2 * vj) + vnext) ^ 2;
        if ecurv(j) > ecurvmax, ecurvmax = ecurv(j); end
        %% calc eimage
        eimage(j) = im(vj(1), vj(2));
        if eimage(j) > eimagemax, eimagemax = eimage(j); end
        if eimage(j) < eimagemin, eimagemin = eimage(j); end
        %% calc ecent
        ecent(j) = norm(mpoint - vj);
        if ecent(j) > ecentmax, ecentmax = ecent(j); end
        if ecent(j) < ecentmin, ecentmin = ecent(j); end
    end

    %% normalize energy terms
    econt = econt ./ econtmax;
    ecurv = ecurv ./ ecurvmax;
    o = ones(1, length(eimage));
    eimage = eimage - (o * eimagemin);
    if eimagemax - eimagemin < 0.02, eimagemin = eimagemax - 0.02; end
    eimage = eimage ./ (eimagemax - eimagemin);
    eimage = o - eimage;
    
    ecent = ecent - (o * ecentmin);
    ecent = ecent ./ (ecentmax - ecentmin);

    %% compute energy functional of each point
    emin = inf;
    jmin = 1;
    for j = 1:length(neighbours)
        %% compute energy function at j
        ej = (alpha(i) * econt(j)) + (beta(i) * ecurv(j)) + (gamma(i) * eimage(j)) + (delta(i) * ecent(j));
        %% check if energy at ej is less than at current min
        if ej < emin
            emin = ej;
            jmin = j;
        end
    end

    %% move point vi to location jmin
    vjmin = neighbours(jmin,:);
    con(i,:) = vjmin;
    if vi(1) ~= vjmin(1) || vi(2) ~= vjmin(2),
        ptsmoved = ptsmoved + 1;
    end;

    %% set previous point
    vlast = vi;
end
Image Contour Extraction

greedyice_movecontour.m
