function Q = closestquadraticpoint(C,P)
% CLOSESTQUADRATICPOINT Find the minimum-distance point on a quadratic 
% 
%  Q = CLOSESTQUADRATICPOINT(C,P) where C = [c b a] for a quadratic
%  curve y=a*x^2 + b*x + c and P=[x y] is a point in the plane. Q
%  is the closest point on C from P.

% Set up the cubic equation derived from the solution to the
% Lagrange multipler problem
% min D(Q,P) = norm(P,Q)^2
% st  dot(C,phi(P(1)) - P(2)

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

% $Id: closestquadraticpoint.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/>.


L = [ 4*C(3)^2 ...                 % x^3
      4*C(3) * C(2) ...            % x^2
      4*C(3) * (C(3)-P(2)) + 2 ... % x
      -2 * P(1) ];

R = roots(L)';

% Eliminate complex roots!
R(imag(R)~=0) = [];

% Candidate points
Q = [R ; C*[ones(size(R)) ; R ; R.^2]]';

% Find closest
if length(R)>1
  Q = Q(argmin( sum( (P(ones(size(Q,1),1),:)-Q).^2, 2) ), :);
end;
  
