function [X,Y,B,H,L,R] = controlpts(bot_coefs,top_coefs,left, right, ...
                                    dim, numPoints)
% CONTROLPTS Calculate control points for text normalization
%
% [X,Y,B,H,L,R] = CONTROLPTS(bot_coefs,top_coefs,left, right)
% where bot_coefs and top_coefs are vectors of the form [c b a] for
% quadratic curves y=a*x^2+b*x+c, left and right are scalars
% representing the horizontal extent of the points to use in the
% original image.  X is a 2xN matrix of control
% points from the original image, and Y is the 2xN matrix of matching
% target points. B is the approximate vertical location for the
% baseline in the target, H is the approximate height of the text in
% the target, and L and R give the horizontal bounds of the text in
% the target.
%
% CONTROLPTS(...,dim) where dim is the size of the original image
% (used to limit the extent of control points).
%
% CONTROLPTS(...,numPoints) where the number of samples along the
% baseline is restricted to approximately numPoints, giving roughly
% N=4*numPoints total control points.
%
% See also TPSIMTFORM
    
% Jerod Weinman
% jerod@acm.org
% (c) 2011
    
% $Id: controlpts.m 594 2014-10-24 12:04:55Z weinman $
 
% Jerod Weinman, Zach Butler, Dugan Knoll, Alexander Turner, Yirong Jing
% jerod@acm.org
% Copyright 2011, 2012, 2013, 2014
 
%    This file is part of PIONEER WordNorm.
%
%    PIONEER WordNorm is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    PIONEER WordNorm is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with PIONEER WordNorm.  If not, see <http://www.gnu.org/licenses/>.

if nargin<5
  dim = [];
end

if nargin<6 || isempty(numPoints)
  numPoints = 100;
end

% Preserve average vertical distance as height 
H = vquaddist(bot_coefs,top_coefs,left, right);

    %%%%
    %%%%
    %%%% NOTE: All of the points below are (x,y) or (col,row) pairs
    %%%%
    %%%%

X = zeros(2,0);
Y = zeros(2,0);




if ~isempty(dim)
  first = max( round(left-2*H),  1);        % Extended starting column
  last  = min( round(right+2*H), dim(2));  % Extended ending column
else
  first = round(left-2*H);        % Extended starting column
  last = round(right+2*H);
end

hoff = 1;   % Accumulating horizontal offset
hasPrev = 0;

% Limit the number of control points for very large images
step = ceil( (last-first)/numPoints );

% Include the left and right bound points for both accuracy and
% calculating target locations
columns = sort([left first:step:last right]);

%% For each column in the source
for col=columns(:)'
  

  [pbot,ptop,plow,phi] = calcpts( col, bot_coefs, top_coefs, H);
  
  %% Update offset
  
  if hasPrev; % Not the first point
    
    % Add the distance between bottom points to the offset
    % (linearized segment length)
    hoff = hoff + norm(pbot - prevBot);
  else
    hasPrev = 1; % Only needs to be done once
  end
  
  %% Add point pairs

  % Source
  X = [X ...
       phi ...  % Ascender
       ptop ... % Top
       pbot ... % Base
       plow ];  % Descender
       
  % Target
  Y = [Y ...
       [ hoff ;     1 ] ... % Ascender
       [ hoff ;   H+1 ] ... % Top
       [ hoff ; 2*H+1 ] ... % Base
       [ hoff ; 3*H+1 ] ];  % Descender

    
  % Update previous bottom point
  prevBot = pbot;

  if col==left
    L = hoff;
  elseif col==right
    R = hoff;
  end
  
end; % for each column


B = 2*H+1; 

%% Filter out-of-bounds source points
if ~isempty(dim)
  ndx = X(1,:) >= 1      & ...
        X(1,:) <= dim(2) & ...
        X(2,:) >= 1      & ...
        X(2,:) <= dim(1);


  X = X(:,ndx);
  Y = Y(:,ndx);
end

%% Filter imaginary points (no intersection along normal)
ndx = imag(X(1,:))==0 & imag(X(2,:))==0;
X = X(:,ndx);
Y = Y(:,ndx);

end % function controlpts

function [pbot,ptop,plow,phi] = calcpts( col, bot_coefs, top_coefs, H)
% CALCPTS  Calculate control points for a column
%
% Inputs
%   col:       Column along bottom guide 
%   bot_coefs: Coefficients of bottom quadratic guide
%   top_coefs: Coefficients of bottom quadratic guide
%   H:         Amount of space to preserve between top and bottom
%              points and extend low and high control points through
%   
% Outputs
%   pbot: Coordinate of bottom guide at col
%   ptop: Corresponding coordinate of top guide
%   plow: Corresponding coordinate of extended point below (i.e., descender)
%   phi:  Corresponding coordinate of extended point above (i.e., ascender)

  %% Calculate the normal line to the bottom quadratic

  row = dot(bot_coefs, [1 col col^2]);
  
  
  % Normal-line slope, offset
  deriv = 2*bot_coefs(3) * col + bot_coefs(2);
  slope = -1/deriv;
  slopeAng = atan(slope);
  offset = col/deriv + row;

  if top_coefs(3) ~=0 % Quadratic?
    
    %% Find intersection(s) of normal with top quadratic
    
    % Line/quadratic intersection equation discriminant root
    %   y = a*x^2 + b*x + c
    %   y = slope*x + offset
    %   a*x^2 + (b-slope)*x + (c-offset) = 0
    %   x = ( -B +/- sqrt( B^2 - 4*A*C ) ) / 2*A
    %   where A = a        = top_coefs(3)
    %         B = b-slope  = top_coefs(2) - slope
    %         C = c-offset = top_coefs(1) - offset
    
    discr = (top_coefs(2) - slope)^2 - ...
            4 * top_coefs(3) * (top_coefs(1) - offset);
    
    root = sqrt(discr);
    
    % Intersection point one
    col0 = (slope - top_coefs(2)  + root) / ( 2 * top_coefs(3) );
    row0 = dot(top_coefs, [1 col0 col0^2]);
    
    % Intersection point two
    col1 = (slope - top_coefs(2)  - root) / ( 2 * top_coefs(3) );
    row1 = dot(top_coefs, [1 col1 col1^2]);
    
    
    
    %% Choose closest intersection
    
    % This assumes the quadratics are somewhat flat
    if norm([row0 col0]-[row col]) < norm([row1 col1]-[row col])
      rowt = row0;
      colt = col0;
    else
      rowt = row1;
      colt = col1;
    end
  else % Linear ==> only one intersection
    colt = (top_coefs(1)-offset)/(slope-top_coefs(2));
    rowt = dot(top_coefs, [1 colt colt^2]);
  end
  
  % Set extension distance/magnitude    
  del = H;
  
  if slope>0
    del = -del;
  end

  
  %% Construct descender extension points
  
  cold = col - del*cos(slopeAng);
  rowd = row - del*sin(slopeAng);

  %% Construct ascender extension points
  cole = colt + del*cos(slopeAng);
  rowe = rowt + del*sin(slopeAng);
  
  %% Final points
  pbot = [col  ; row];
  ptop = [colt ; rowt];
  plow = [cold ; rowd];
  phi  = [cole ; rowe];
end % function calcpts
  
