%FINDCOORD   Find coordinates of nonzero elements.
%   I = FINDCOORD(B) returns the coordinates of the image B that
%   are non-zero. I(N,:) is a vector with coordinates for non-zero
%   pixel number N.
%
%   FINDCOORD is similar to FIND, in that it returns the same
%   list of pixels, but in a different form.
%
%   Note that the output of FINDCOORD cannot be directly used to
%   index an image. Each of the coordinates should be used separately,
%   i.e., B(I(N,1),I(N,2)), or by computing the pixel indices:
%      sz = size(b);
%      strides = cumprod([1,sz([2,1,3:end])]);
%      strides = strides([2,1,3:end-1])';
%      b(I*strides)
%   Here, I*strides is the same as what is returned by the function
%   FIND. Note that the re-ordering of the size and stride arrays is
%   required because the data is stored column-wise in MATLAB.
%
%   See also FIND, COORD2IMAGE.

% (C) Copyright 1999-2007               Pattern Recognition Group
%     All rights reserved               Faculty of Applied Physics
%                                       Delft University of Technology
%                                       Lorentzweg 1
%                                       2628 CJ Delft
%                                       The Netherlands
%
% Cris Luengo, September 2000.
% 22 June 2001:      Fixed bug caused by FIND on 1D data.
%                    Also, 1D images produced two coordinates.
% 27 August 2001:    Fixed important bug: coordinates returned were bogus
%                    if x and y sizes were not the same.
%                    Fixed bug in documentation.
% 11 September 2001: Fixed previous bug correctly. :{

function C = findcoord(in)
if ~isscalar(in), error('Parameter "in" is an array of images.'); end
if in.dims < 1 | isempty(in.data)
   error('Input image is empty or 0-D.');
end
try
   I = find(in.data);
   if isempty(I)
      C = [];
   else
      I = I-1;
      I = reshape(I,length(I),1);  % For 1D data, FIND returns a row vector!
      % translate indices in I to coordinates C.
      if in.dims==1
         C = I;
      else
         C = zeros(length(I),in.dims);
         sz = size(in.data);
         sz = [sz,ones(1,in.dims-length(sz))];
         k = [1,cumprod(sz(1:end-1))];
         for ii=in.dims:-1:2
            C(:,ii) = floor(I/k(ii));
            I = I-(C(:,ii)*k(ii));
         end
         % Switch the X and Y coordinates!
         C(:,1) = C(:,2);
         C(:,2) = I;
      end
   end
catch
   error(di_firsterr)
end
