function h = cdata_direct(h, clev, varargin)
% function h = cdata_direct(h, clev, <i0>, <'data', z> )
%
% Converts patches (such as made by contourf) or
% surfaces (such as made by pcolor) to
% indexed colors ('CDataMapping', 'direct').
%
% Inputs:
%  (required)
%   h     Vector of patch and/or surface handles.  For example,
%            this is the second output argument of contourf.
%   clev  Vector of color boundaries, such as the
%             contour levels (boundaries) *used* in
%             the call to contourf.  It should start and
%             end with extremes, so that all data fall
%             within its limits.
%  (optional)
%   i0    Index in the colormap where the first range in
%             clev will fall.  Default is 1.
%             Other values are used to pick a section of
%             a concatenated colormap.
%         Or, if i0 is a vector, it will be interpreted
%             directly as the set of indices in the colormap
%             that correspond to the ranges in clev.  In
%             this case, you must have
%             length(i0) = length(clev) - 1
%
%   'data', [zmin, zmax] or <full data array>
%             zmin, zmax are min and max of contoured data;
%             or pass in the original data array and zmin
%             and zmax will be calculated from it.
%             If this keyword-value pair is present, patch
%             objects for clev values outside this range will
%             be deleted. Patches with CData = NaN are retained,
%             so that internal data holes are blanked out.
%
% Output:
%   h     Vector of patch handles, as updated by deletion
%             of patch objects if the 'datarange' option
%             was used.
%
% See also: cbar_direct, a companion routine to make a colorbar
%             using the same indexing scheme.
%
% EF 2001/10/17, generalized from cf_color, which it replaces.
%
% EF 2001/11/15, added 'data', [zmin, zmax] or <data array> argument pair.

drange = [];
i0 = 1;

if nargin < 2,
   help cdata_direct
   error('Need at least 2 arguments.');
end
if nargin == 3 | nargin == 5,
   i0 = varargin{1};
   varargin(1) = [];          % Chop out i0 leaving nothing, or 2 args.
end
if length(varargin) == 2
   if ischar(varargin{1}) & strcmp(varargin{1}, 'data'),
      z = varargin{2};
      drange = [min(z(:)), max(z(:))];
   else
      help cdata_direct
      error('Expected datarange keyword-value pair.')
   end
end

nlev = length(clev);
if length(i0) == 1,
   indices = (i0 - 1) + (1:(nlev-1));
elseif length(i0) == (nlev-1),
   indices = i0;
else
   help cdata_direct
   error('length of i0 must be either 1 or length(clev)-1')
end

cd = get(h, 'CData');
nh = length(h);
if nh == 1,
   set(h, 'CData', to_index(cd, clev, indices));
else
   for ii = 1:nh,
      set(h(ii), 'CData', to_index(cd{ii}, clev, indices));
         % This must be in the loop; for a vector
         % of handles, a property can only be set to a
         % scalar, even though one can "get" a cell
         % array of values.
   end
end

set(h, 'CDataMapping', 'direct');

if length(drange) == 2 & iscell(cd),
   ibad = logical(ones(nh, 1));
   for ii = 1:nh,
      if isnan(cd{ii}) | (cd{ii} >= drange(1) & cd{ii} <= drange(2)),
         ibad(ii) = 0;
      end
   end
   delete(h(ibad));
   h(ibad) = [];
end

%----------------------------------------
function ind = to_index(cd, clev, indices)
nlev = length(clev);
ind = NaN * cd;
for ilev = 1:(nlev-1),
   ii = find(cd >= clev(ilev) & cd < clev(ilev+1));
   ind(ii) = indices(ilev);
end


