function points = find_target_warp(reference_image, input_image, image_width)
% FIND_TARGET_WARP: Finds the warp curve that transforms the input image
% into the reference space
%
%    GENERAL
%
%      points = find_target_warp(reference_image, input_image)
%
%    INPUT/S
%
%      -reference_image:
%           The refrence image
%
%      -input_image:
%           The raw image to be warped to some reference
%
%      -input_width:
%           The width of the images
%           
%    OUTPUT/S
%
%      -points:
%           The warp curve points after correct warp has been found
%
%    PENDING WORK
%
%      -March 2004: Works for just image width 50 and needs to be extended.
%      --FIXED March 2004
%
%    KNOWN BUG/S
%
%      -None.
%
%    COMMENT/S
%
%      -DOES NOT YET WORK WITH OFFSETS. IT ASSUME BUMPS FLAT AT 0.
%
%    RELATED FUNCTION/S
%
%      LINEAR_WARP
%
%    ABOUT
%
%      -Created:     March 7th, 2004
%      -Last update: March 22nd, 2004
%      -Revision:    0.0.4
%      -Author:      R. S. Schestowitz, University of Manchester
% ==============================================================

% firstly find the edges. Will fail is smoothing is applied because it is a
% very data-specific term.

for i=1:image_width-1,
  if (input_image(i) == 0 & input_image(i + 1) ~= 0)
    lo = i;
  elseif (input_image(i) ~= 0 & input_image(i + 1) == 0)
    hi = i;
  end
  if (reference_image(i) == 0 & reference_image(i + 1) ~= 0)
    ref_lo = i;
  elseif (reference_image(i) ~= 0 & reference_image(i + 1) == 0)
    ref_hi = i;
  end
end
% lo is the lower point where the bump begins (from left to right) and hi
% is the point where the bump goes flat

if ((hi - lo) <= 0) 
    % RSS: Check for exceptions here
  disp('Invalid values of high and low. Low must be lower than high');
  return; % break out of function
end

% Smith: line up the edge points, do linear stretch in between!

% Smith: first section - between 1 and lo
% scale_factor1 = (ref_lo-1)/(lo-1);

points(1:image_width) = 1:image_width;

% set up an initial curve

scale_factor1 = (lo - 1) / (ref_lo - 1);

% factors are the a in f = ax + b

points(1:ref_lo) = (points(1:ref_lo) - 1) * scale_factor1 + 1;

% second section of the curve

scale_factor2 = (hi - lo)/(ref_hi - ref_lo);
points(ref_lo + 1:ref_hi - 1) = (points(ref_lo + 1:ref_hi - 1) - ref_lo) * scale_factor2 + lo;

% scale_factor2 = (ref_hi - ref_lo)/(hi-lo);
% points(:,lo+1:hi-1) = (points(:,lo+1:hi-1) - lo)*scale_factor2 + ref_lo;

% third section of the curve
scale_factor3 = (image_width - hi) / (image_width - ref_hi);
points(ref_hi:image_width) = ...
  (points(ref_hi:image_width) - ref_hi) * scale_factor3 + hi;

%scale_factor3 = (image_width-ref_hi)/(image_width-hi);
%points(:,hi:image_width) = ...
%  (points(:,hi:image_width) - hi)*scale_factor3 + ref_hi;