%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. 
%
%   [I,V] = FINDCOORD(B) returns also a 1-D image containin the vaules of 
%   the non-zero pixels in B.
%   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. :{
% 19 November 2008:  Added return of pixel values (BR)

function [C,V] = findcoord(in)
if prod(builtin('size',in)) ~= 1
   error('Parameter "in" is an array of images.')
else
   if in.dims < 1 | isempty(in.data)
      error('Input image is empty or 0-D.');
   end
   try
      I = builtin('find',in.data);
      if nargout>1
         V = dip_image(in.data(I),in.dip_type);
      end
      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
end
