function [pt_array2,pt2i,pt_array1,pt1i,dist] = closest(array1,array2);
%CLOSEST - Find point in one array closest to point in another array
% function [pt_array2,pt2i,pt_array1,pt1i,dist] = closest(array1,array2);
% find point the row in array2 closest in a Euclidean sense to locations in 
% array1.  return as pt_array2.
% Optionally, pt_array1 is point in array1 closest to array2.
% Each row one multidimensional location.
% Optionally, return index pt2i and pt1i giving corresponding row number
% of closest elements.
% dist is the Euclidean distance between the two points
%
% (Note reversal of output arguments compared with input.  This reversal
% makes this function compatible with earlier versions of "closest".
% Can be called as closest([x y z],array), returning member in array.)
%
% See also CLUSTER, NEIGHBOR

%<autobegin> ---------------------- 26-May-2004 11:29:53 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Utility - Numeric
%
% At Check-in: $Author: Mosher $  $Revision: 14 $  $Date: 5/26/04 9:59a $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 24-May-2004
% 
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
%    University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
%    Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
%    CNRS, Hopital de la Salpetriere, Paris, France
% 
% See BrainStorm website at http://neuroimage.usc.edu for further information.
% 
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed  under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%<autoend> ------------------------ 26-May-2004 11:29:53 -----------------------

% ----------------------------- Script History ---------------------------------
% JCM        1991  Creation
% JCM 20-Jul-1993  Handle multiple points
% modified JCM 10/29/99 to use repmat function instead of kron for speed
%  use sum(x,2) instead of sum(x',2);
%  probably should consider shift to column vectors, but legacy code here.
% JCM 19-May-2004  Comments cleaning
% ----------------------------- Script History ---------------------------------


[m2,n2] = size(array2);
[m1,n1] = size(array1);
if(n2 ~= n1),
  error(sprintf(...
  'CLOSEST: Arrays are not of same column size: %.0f, %.0f',n1,n2));
end

if(0) % old code
   array2 = kron(ones(m1,1),array2);
   array1 = kron(array1,ones(m2,1));
else % new code 10/29/99
   if(m1 > 1), % multiple rows in the first array
      array2 = repmat(array2,m1,1);
      array1 = repmat(array1,[1 1 m2]);
      array1 = shiftdim(array1,2);
      array1 = reshape(array1(:),[m1*m2 n1]); % same as above kron without the multiply
   else % only a single row in array1
      array1 = array1(ones(m2,1),:); % simple rep out to match array2
   end
end

darray = array2 - array1;  % all combinations
if(n1 > 1),
  darray = sum(darray.^2,2); % squared Euclidean distance, no need for sqrt
  [dist,i] = min(darray);  % closest point
  dist = sqrt(dist);  % Euclidean distance
else
  darray = abs(darray); % one dimensional vector
  [dist,i] = min(darray);  % closest point
end

pt_array2 = array2(i,:);

if(nargout > 1),
   pt2i = rem((i-1),m2)+1;  % index in orginal arrays
else
   return % only wanted the row vector, not the index
end

if(nargout > 2), % also want the location in the first array
   pt_array1 = array1(i,:);
   pt1i = floor((i-1)/m2)+1;
end
   
return
