function [snake_pnts,e] = snake(pnts, alpha, beta, max_delta_x, resol_x, ...
                          max_delta_y, resol_y, feat_img)
% SNAKES
% By Chris Bregler and Malcolm Slaney, Interval Technical Report IRC 1995-017
% Copyright (c) 1995 Interval Research Corporation.
%
% Usage
%
% [snake_pnts,e] = snake(pnts, alpha, beta, ...
%                        max_delta_y, resol_y, max_delta_x, resol_x, feat_img)
%
% This function computes one iteration of the energy-minimization of
% active contour models described in the paper "Using Dynamic Programming
% for Solving Variational Problems in Vision" by Amir A. Amini, Terry E.
% Weymouth, and Ramesh C. Jain, IEEE Transactions on Pattern Analysis and
% Machine Intelligence, Vol. 12, No. 9, September 1990, pp 855-867.
%
% Snakes align a contour to some feature maxima in an image (for example)
% image boundaries) using dynamic programming.  The quality of the
% alignment is measured by an energy function consisting of an internal
% "contour-smoothness-term" and an external feature term (equation 50 in the
% paper).  Minimizing this energy term leads to a contour that trade offs
% these two criterias (smoothness + maximal feature responses) in some desired
% way.
%
% 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
%

% SNAKES - A MatLab MEX file to demonstrate snake contour-following.
% This Software was developed by Chris Bregler and Malcolm Slaney of
% Interval Research Corporation.
% Copyright (c) 1995 Interval Research Corporation.
%
% This is experimental software and is being provided to Licensee
% 'AS IS.'  Although the software has been tested on a PowerMac
% 8100 running version 4.2c of MatLab with MEX support and on an
% SGI running version 4.2c, Interval makes no warranties relating
% to the software's performance on these or any other platforms.
%
% Disclaimer
% THIS SOFTWARE IS BEING PROVIDED TO YOU 'AS IS.'  INTERVAL MAKES
% NO EXPRESS, IMPLIED OR STATUTORY WARRANTY OF ANY KIND FOR THE
% SOFTWARE INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY OF
% PERFORMANCE, MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
% IN NO EVENT WILL INTERVAL BE LIABLE TO LICENSEE OR ANY THIRD
% PARTY FOR ANY DAMAGES, INCLUDING LOST PROFITS OR OTHER INCIDENTAL
% OR CONSEQUENTIAL DAMAGES, EVEN IF INTERVAL HAS BEEN ADVISED OF
% THE POSSIBLITY THEREOF.
%
%   This software program is owned by Interval Research
% Corporation, but may be used, reproduced, modified and
% distributed by Licensee.  Licensee agrees that any copies of the
% software program will contain the same proprietary notices and
% warranty disclaimers which appear in this software program.

  if resol_y < 1; resol_y = 1; end;
  if resol_x < 1; resol_x = 1; end;

  n = size(pnts,1);
  [row,col] = size(feat_img);
  target = reshape(feat_img,row*col,1);
  scan_y = -max_delta_y:resol_y:max_delta_y;
  scan_x = -max_delta_x:resol_x:max_delta_x;
  num_scan_y = size(scan_y,2);
  num_scan_x = size(scan_x,2);
  num_states = num_scan_y * num_scan_x;

  fprintf('n = %d; num_states = %d; ',n,num_states);

  delta_x = ones(num_scan_y,1)*scan_x;
  delta_y = scan_y'*ones(1,num_scan_x);
  delta_x = reshape(delta_x,1,num_states);
  delta_y = reshape(delta_y,1,num_states);

  states_x = round(pnts(:,1))*ones(1,num_states) + ones(n,1)*delta_x;
  states_y = round(pnts(:,2))*ones(1,num_states) + ones(n,1)*delta_y;

  % take care of boundary cases
  states_x = min(max(states_x,1),col);
  states_y = min(max(states_y,1),row);

  states_i = (states_x-1)*row + states_y;

  Smat = zeros(n,num_states^2);
  Imat = zeros(n,num_states^2);

  % forward pass

  for v2 = 1:num_states,
    Smat(1,(v2-1)*num_states+1:v2*num_states) = ...
      -target(states_i(1,:))';
  end;

  for k = 2:n-1,
    fprintf('.');  % debug
    for v2 = 1:num_states, for v1 = 1:num_states,
      v0_domain = 1:num_states;
      [y,i] = min( Smat(k-1,(v1-1)*num_states+v0_domain) ...
              + alpha(k)*( (states_x(k,v1)-states_x(k-1,v0_domain)).^2 ...
                      + (states_y(k,v1)-states_y(k-1,v0_domain)).^2) ...
              + beta(k)*( (states_x(k+1,v2)-2*states_x(k,v1) ...
                       + states_x(k-1,v0_domain)).^2 ...
                     + (states_y(k+1,v2)-2*states_y(k,v1) ...
                       + states_y(k-1,v0_domain)).^2) );
      Smat(k,(v2-1)*num_states+v1) = ...
        y-target(states_i(k,v1));
      Imat(k,(v2-1)*num_states+v1) = i;
    end; end;
  end;

  for v1 = 1:num_states,
    v0_domain = 1:num_states;
    [y,i] = min( Smat(n-1,(v1-1)*num_states+v0_domain) ...
            + alpha(n)*( (states_x(n,v1)-states_x(n-1,v0_domain)).^2 ...
                    + (states_y(n,v1)-states_y(n-1,v0_domain)).^2));
    Smat(n,v1) = y-target(states_i(n,v1));
    Imat(n,v1) = i;
  end;

  [e,final_i] = min(Smat(n,1:num_states));


  % backward pass

  snake_pnts = zeros(n,2);

  snake_pnts(n,:) = [states_x(n,final_i),states_y(n,final_i)];
  v1 = final_i; v2 = 1;
  for k=n-1:-1:1,
    v = Imat(k+1,(v2-1)*num_states+v1);
    v2 = v1; v1 = v;
    snake_pnts(k,:) = [states_x(k,v1),states_y(k,v1)];
  end;

  fprintf('\n');