function [B,pad] = tpsimtform(A,X,Y,interp,fill)
% TPSIMTFORM Perform a thin-plate-spline transformation using control points
%
% B = TPSIMTFORM(A,X,Y) where A is a grayscale double image, and X
% and Y are the 2xN control points in (x,y) format for the input
% (A) and output (B) image, respectively.
%
% TPSIMTFORM(...,interp) where interp is the interpolation option string
% given to IMTRANFORM, and may be one of bilinear, bicubic, or
% nearest. 
  

% Jerod Weinman
% jerod@acm.org
% Copyright (c) 2010
    
% $Id: tpsimtform.m 63 2013-03-16 15:41:50Z weinman $
 
% Jerod Weinman
% jerod@acm.org
% /home/data/svn/vision/tools/package/GPL_STATEMENT
 


    
    if nargin<4 || isempty(interp)
      interp = 'bicubic';
    end;
    
    if nargin<5 || isempty(fill)
      fill = 'none';
    end

    pad = [0;0];
    if ~strcmp(fill,'none') % Fill desired?
          
      % Check input value ranges, padding where input points are
      % out of the input image bounds
      if any(X(1,:) < 1 ) % Left
        pad(1) = 1+abs(floor(min(X(1,:))));
        A = padarray( A, [0 pad(1)],fill,'pre');
        X(1,:) = X(1,:) + pad(1);
      end

      if any(X(2,:) < 1 ) % Top
        pad(2) = 1+abs(floor(min(X(2,:))));
        A = padarray( A, [pad(2) 0],fill,'pre');
        X(2,:) = X(2,:) + pad(2);
      end
      
      if any(X(1,:) > size(A,2) ) % Right
        A = padarray( A, [0 ceil(max(X(1,:)))-size(A,2)],fill,'post');
      end

      if any(X(2,:) > size(A,1) ) % Bottom
        A = padarray( A, [ceil(max(X(2,:)))-size(A,1) 0],fill,'post');
      end
    end
      
      
    % Calculate the inverse thin-plane-spline function from control
    % points X (input) and Y (output)
    fn = tpaps(Y,X);
    
    % Determine integral size of output (nrows, ncols)
    sz = fliplr(ceil(max(Y')));

    % Construct custom image transform
    T = maketform('custom',2,2, [], @(Y,tt)fnval(fn,Y')',[]);

    % Apply transform
    B = imtransform(A,T,interp, ...
                    'XData',[1 sz(2)], ...
                    'YData',[1 sz(1)], ...
                    'XYScale', 1);
