snake.m

Tracking a moving object through several frames, provided changes from frame to frame are on the ord

function [snake_pnts,e] = snake(pnts,...
                                alpha, beta, gamma,...
                                max_delta_x, resol_x,...
                                max_delta_y, resol_y,...
                                feat_img,...
                                feat_img_gradient)
% Usage
%
% [snake_pnts,e] = snake(pnts, alpha, beta, ...
%                        max_delta_y, resol_y, max_delta_x, resol_x, feat_img)
%
% Inputs:
%   pnts          Starting contour. Each row is a [x,y] coordinate.
%   alpha         Energy contributed by the distance between control points.
%                 Set to zero if length of slope doesn't matter.
%   beta          Energy contributed by the curvature of the snake.  Larger
%                 values of beta cause bends in the snake to have a high cost
%                 and lead to smoother snakes.
%   max_delta_y   Max number of pixels to move each contour point vertically
%   resol_y       Contour points will be moved by multiples of resol_y
%   max_delta_x   Analog to max_delta_y
%   resol_x       Analog to resol_y
%   feat_img      2D-Array of the feature responses in the image.  For example
%                 it can contain the magnitude of the image gradients
%
% Outputs:
%   snake_pnts    New contour points.
%   e             Energy value of these new contour points
%

nPnts = size(pnts);

nSnakePnts = nPnts;
nSnakePnts(1) = nSnakePnts(1) + 2;
new_snake_pnts = zeros( nSnakePnts );

% taking care of the boundary cases
% needs an improvement
new_snake_pnts(2:nPnts(1)+1, :) = pnts;
new_snake_pnts(1, :) = pnts(2, :);
new_snake_pnts(nSnakePnts(1), :) = pnts(nPnts(1) - 1, :);

% set up indices for the square
% neighborhood

index_x = -max_delta_x:max_delta_x;
index_y = -max_delta_y:max_delta_y;

sz_index_x = size(index_x);
sz_index_y = size(index_y);

imSize = size( feat_img );

ab_index = 0:nPnts(1);

% loop through all points
for n = 2:nSnakePnts(1)-1
  % set E_min for this point to a very
  % big value
    E_min = 1000000;
  %  
    cur_aplha = alpha(ab_index(n));
    cur_beta  = beta(ab_index(n));
    cur_gamma = gamma(ab_index(n));
  % set the deltas corresponding to the min
  % energy to 0's
    delta_x_min = 0;
    delta_y_min = 0;
  % for each point in the neighborhood of the point new_snake_pnts(n)...
    for i = 1:sz_index_y(2)
        for j = 1:sz_index_x(2)
          % calculate deltas to move the point
          % in the neighorhood
            delta_y = index_y(i);
            delta_x = index_x(j);
          % move the current point according to deltas
            pnt = new_snake_pnts(n, :);
            pnt = pnt + [delta_x delta_y];

            savePnt = pnt;

            pnt(1) = clipValue(pnt(1), 1, imSize(2));
            pnt(2) = clipValue(pnt(2), 1, imSize(1));

            delta_x = (pnt(1) - savePnt(1)) + delta_x;
            delta_y = (pnt(2) - savePnt(2)) + delta_y;

          % calculate a new energy estimate
            E_j = cur_aplha*continuityEstimate(n,...
                                               new_snake_pnts,...
					       pnt,...
					       feat_img) +...
                  cur_beta*curvatureEstimate(n,...
					     new_snake_pnts,...
					     pnt,...
					     feat_img) +...
                  cur_gamma*imageForces(pnt,...
					10,...
					10,...
					feat_img,...
                                        feat_img_gradient);
          % greedy strategy
            if E_j <= E_min,
                E_min = E_j;
                delta_x_min = delta_x;
                delta_y_min = delta_y;
            end
        end
    end
  % change the current point according to the min
  % energy estimate
    new_snake_pnts(n, :) = new_snake_pnts(n, :) + [delta_x_min delta_y_min];
end

    snake_pnts = new_snake_pnts(2:nPnts(1)+1, :);
end