function points = cps_warp_1d(points, warp_centre, r, d)
% CPS_WARP_1D: Apply a Clamped plate spline warp to points
%
% Code written by Katherine Smith, 2003
%
%    GENERAL
%
%      points = cps_warp_1d(points, warp_centre, r, d)
%
%    INPUT/S
%
%      -points:
%           The coordinates of the points before warping.
%
%      -warp centre:
%           co-ordinate of the centre of the warp, 
%
%      -r:
%           Radius of the region to be warped
%
%      -d:
%           The magnitude of the warp (in the same units as
%           the points)
%           
%    OUTPUT/S
%
%      -points:
%           The coordinates of the points after warping.
%
%    PENDING WORK
%
%      -Rework these comments. Many of them must be inaccurate.
%
%    KNOWN BUG/S
%
%      -
%
%    COMMENT/S
%
%      -
%
%    RELATED FUNCTION/S
%
%      
%
%    ABOUT
%
%      -Created:     November 23rd, 2003
%      -Last update: Novermber 25th, 2003
%      -Revision:    0.0.2
%      -Author:      R. S. Schestowitz, University of Manchester
% ==============================================================

% maximum displacement which is diffeomorphism
max_d = 0.67;
%d = d/r;
if abs(d) > max_d
  error('Non-diffeomorphic shift');
end
%if (warp_centre-r <0) | (warp_centre+r > max(points))
if (warp_centre-r <-1) | (warp_centre+r > 1)
  error(['Warp out of bounds warp_centre: ',num2str(warp_centre),' r: ',num2str(r)]);
end

% Smith: (RSS: is this correct?) normalise to -1 to +1
points = (points - warp_centre)/r;

% for each point, find and apply G (which scales the shift)
points_to_warp = points((points>= -1) & (points<= 1));
G = zeros(size(points));
G((points>= (-1)) & (points<= 1)) = ...
  ((ones(size(points_to_warp))*1/3 ...
  - points_to_warp.^2 ...
  + 2/3*(abs(points_to_warp).^3)));

points = points + G*3*d;

% and rescale
points = points*r + warp_centre;
