function [bottoms_img,tops_img] = normalextremalocal(bin_img,coeff)
% NORMALEXTREMALOCAL Component local extrema oriented to normal of curve
%
%   [BOTTOMS_IMG,TOPS_IMG] = NORMALEXTREMALOCAL(BIN_IMG,COEFF)
%   For parameters BIN_IMG consisting of binarized text.  The
%   parameter COEFF contains a vector containing [c b a], the
%   coefficients of the quadratic model y=ax^2+bx+c. The curve
%   represented by COEFF is an estimate of the skew/curvature of the
%   text.  This function returns an image, BOTTOMS_IMG, consisting of
%   the bottom points of the text. BOTTOMS_IMG is the same size as
%   BIN_IMG. Likewise TOPS_IMG
%
%   For each column along the line, and for each connected component
%   of BIN_IMG, we find the farthest point in the component along the
%   normal below the curve COEFF and store that point in
%   BOTTOMS_IMG. Likewise points above COEFF for TOPS_IMG.
%   
%   Example Usage:
%
%   bin_img = gray_image>128;
%   coeff = initregress();
%   coeff = pruneregress(bin_img,coeff);
%   bottoms_img = normalextrema2(bin_img,coeff);
%
%   See also INITREGRESS, PRUNEREGRESS.

% Copyright (c) 2012
% Jerod Weinman
% jerod@acm.org

% $Id: normalextremalocal.m 595 2014-10-24 18:15:54Z 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/>.

% Total size of the image
sz = size(bin_img);

% Create a mask image for the tops and bottoms
bottoms_img = false(sz);
tops_img = false(sz);

% Run a connected components analysis
CC = bwconncomp(bin_img);

if(CC.NumObjects == 0)
    bottoms_img = bin_img;
    tops_img = bin_img;
    return;
end
 
B = coeff;

[rows,cols] = find(bin_img);

% Horizontal range of columns
cc = min(cols):max(cols);

% Vertical range of points
startRow = min(rows);
endRow = max(rows);

for col=cc(1:end) %% for each column along the curve %%

  row = dot(B,[1 col col^2]); % corresponding row

  % along normal, last column within region
  endCol = (row-endRow)*(2*B(3)*col+B(2)) + col;
  % along normal, first column within region
  startCol = (row-startRow)*(2*B(3)*col+B(2)) + col;

  %% Calculate normal line for tops
  xWidth = startCol-col; %run (should be non-positive)
  yWidth = startRow-row; %rise
  lnlen = sqrt(xWidth^2 + yWidth^2); % The length of our line
  slope = yWidth / xWidth; % Determine slope of normal line as rise over run.
 
  % Create enough x values to predict all of the y values.
  x = (xWidth):(-xWidth/(2*lnlen)):0;
  % Calculate y values from slope.
  y = x*slope;
  % Round to nearest pixel.
  x = round(x + col);
  y = round(y + row);
  x = max(min(x,sz(2)),1);
  y = max(min(y,sz(1)),1);
   
  indices_top = sub2ind(sz,y,x);
  

  %% Calculate normal line for bottoms
  xWidth = endCol-col; %run
  yWidth = endRow-row; %rise
  lnlen = sqrt(xWidth^2 + yWidth^2); % The length of our line
  slope = yWidth / xWidth; % Determine slope of normal line as rise over run.
 
  % Create enough x values to predict all of the y values.
  x = 0:(xWidth/(2*lnlen)):(xWidth);
  % Calculate y values from slope.
  y = x*slope;
  % Round to nearest pixel.
  x = round(x + col);
  y = round(y + row);
  x = max(min(x,sz(2)),1);
  y = max(min(y,sz(1)),1);
   
  indices_bot = sub2ind(sz,y,x);
  
  for comp=1:CC.NumObjects
    
    %% For each component, find farthest points of component along normal line

    comp_ndx = CC.PixelIdxList{comp};
    bot_intersect_ndx = myintersect(indices_bot, comp_ndx);
    top_intersect_ndx = myintersect(indices_top, comp_ndx);

    %Bottom
    if ~isempty(bot_intersect_ndx)
      % Convert indices to subscripts
      [pr,pc] = ind2sub(sz,bot_intersect_ndx);
      [cr,cc] = ind2sub(sz,comp_ndx);
      
      % Choose the curve/component intersection point with the greatest
      % (squared) distance from (row,col) - point on the curve
      pt = argmax( (pr-row).^2 + (pc-col).^2 );
      
      % Translate the coordinate system so that the farthest
      % intersecting point pt is the origin
      crt = cr-pr(pt);
      cct = cc-pc(pt);
      
      % Slope of tangent line
      slope = 2*col*B(3) + B(2);
      
      % Rotate the axes so that the normal is vertical
      ang = -atan(slope);
      cctr = round(cct*cos(ang) - crt*sin(ang));
      crtr = round(cct*sin(ang) + crt*cos(ang));
      
      % Trace the bottom contour of the rotated points
      cctru = unique(cctr); % Unique columns: x-coordinates
      crtrb = arrayfun( @(c) max( crtr(cctr==c) ), ... % Max row: y-coordinates
                        cctru); 

      % Because the farthest point is at the origin, we want both
      % neighbors to be non-positive (which NOT BELOW in ij coords)

      nbr = abs(cctru)<=5; % Neighborhood (5 pixels to either side)
      rnbr = crtrb(nbr);   % Row of pixels in neighborhood
      
       if ~any(rnbr > 0) % Bottom-most point in neighborhood
         % Set a pixel in the mask for that point
         bottoms_img(pr(pt),pc(pt)) = 1;
         bottom_pt = [pr(pt),pc(pt)];
      end;
    end;
    
    % Add Top point along curve
    if ~isempty(top_intersect_ndx)
      
      % Convert indices to subscripts
      [pr,pc] = ind2sub(sz,top_intersect_ndx);
      [cr,cc] = ind2sub(sz,comp_ndx);
      
      % Choose the curve/component intersection point with the greatest
      % (squared) distance from (row,col) - point on the curve
      pt = argmax( (pr-row).^2 + (pc-col).^2 );
      
      % Translate the coordinate system so that the farthest
      % intersecting point pt is the origin
      crt = cr-pr(pt);
      cct = cc-pc(pt);
      
      % Slope of tangent line
      slope = 2*col*B(3) + B(2);
      
      % Rotate the axes so that the normal is vertical
      ang = -atan(slope);
      cctr = round(cct*cos(ang) - crt*sin(ang));
      crtr = round(cct*sin(ang) + crt*cos(ang));
      
      % Trace the top contour of the rotated points
      cctru = unique(cctr); % Unique columns: x-coordinates
      crtrt = arrayfun( @(c) min( crtr(cctr==c) ), ... % Min row: y-coordinates
                        cctru); 

      % Because the farthest point is at the origin, we want both
      % neighbors to be non-negative (which is NOT ABOVE in ij coords)
      nbr = abs(cctru)<=5; % Neighborhood (5 pixels to either side)
      rnbr = crtrt(nbr);   % Row of pixels in neighborhood
      
      if ~any(rnbr < 0) % Top-most point in neighborhood
        % Set a pixel in the mask for that point
        tops_img(pr(pt),pc(pt)) = 1;
        top_pt = [pr(pt),pc(pt)];
      end
      
    end
    
  end % comp
end % col 
