function [h] = fscatter3(X,Y,Z,C,cmap);
% [h] = fscatter3(x,y,z,C,cmap);
% Plots point cloud data in cmap color classes and 3 Dimensions,
% much faster and very little memory usage compared to scatter3 !
% x,y,z,C are vectors of the same length
% with C being used as index into colormap (can be any values though)
% cmap is optional colourmap to be used
% h are handles to the line objects

% Felix Morsdorf, Jan 2003, Remote Sensing Laboratory Zuerich

%  Copyright (c) 2010, Felix Morsdorf
%  All rights reserved.
%
%  Redistribution and use in source and binary forms, with or without
%  modification, are permitted provided that the following conditions are
%  met:
%
%      * Redistributions of source code must retain the above copyright
%        notice, this list of conditions and the following disclaimer.
%      * Redistributions in binary form must reproduce the above copyright
%        notice, this list of conditions and the following disclaimer in
%        the documentation and/or other materials provided with the distribution
%
%  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
%  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%  POSSIBILITY OF SUCH DAMAGE.

if nargin == 3
  C = Z;
end
if nargin <= 4
  cmap = [];
end
filled = 0;
if strfind(cmap,'filled')
  cmap=[];
  filled = 1;
elseif strfind(cmap,'bubble')
  cmap=[];
  filled = 2;
elseif ischar(cmap)
  cmap='';
end
if isempty(cmap)
  cmap = jet(256);
end
numclass = max(size(cmap));
if numclass == 1
  cmap = jet(256);
  numclass = 256;
end

% avoid too many calculations
if ~isreal(C), C=abs(C); end
index=find(~isfinite(C(:)));
C(index)=0;
mins = min(C(:));
maxs = max(C(:));
if mins == maxs, h=[]; return; end
minz = min(Z);
maxz = max(Z);
minx = min(X);
maxx = max(X);
miny = min(Y);
maxy = max(Y);
C(index)=nan;

% construct colormap :

col = cmap;

% determine index into colormap
ii = round(interp1([mins maxs],[1 numclass],C));
ish = ishold;
hold on
colormap(cmap);

% plot each color class in a loop

if ~filled
  marker = max(2, 10-log10(length(X(:))));
end

k = 0;
for j = 1:numclass
  jj = find(ii == j);
  if ~isempty(jj)
    k = k + 1;
    if filled
      marker = ceil(20*sqrt((C(jj(1))-mins)/(maxs-mins))+2);
    end
    if filled == 2
      h(k) = plot3(X(jj),Y(jj),Z(jj),'o','color',col(j,:), ...
		 'markersize',marker/2);
    else
      h(k) = plot3(X(jj),Y(jj),Z(jj),'.','color',col(j,:), ...
		 'markersize',marker);
	  end
  end  
end

if ~ish, hold off; end
