function [cout,map] = imapprox(a,cm,n,flag)
%IMAPPROX Approximate indexed image by an image with less colors.
%   [Y,NEWMAP] = IMAPPROX(X,MAP,N) approximates the indexed image X
%   with colormap MAP by another indexed image Y with colormap
%   NEWMAP.  NEWMAP will have N colors which will be chosen
%   to best approximate the colors in MAP.
%
%   [Y,NEWMAP] = IMAPPROX(X,MAP,TOL) approximates the colors in the
%   indexed image by choosing colors via uniform quantization 
%   in RGB space. NEWMAP will contain at most (FLOOR(1/TOL)+1)^3 
%   colors.  TOL must be between 0 and 1.0. 
%
%   B = IMAPPROX(X,MAP,NEWMAP) approximates the colors in MAP by
%   the closest match in NEWMAP.
%
%   [Y,NEWMAP] = IMAPPROX(X,MAP) uniformly removes colors
%   from the colormap if necessary so that NEWMAP has fewer
%   than 256 colors.  The removed colors are approximated
%   by other colors in the colormap.
%
%   [...] = IMAPPROX(...,'dither_option') allows you to enable
%   or disable dithering.  Possiblities are 'dither' or 'nodither'.
%   'dither' is the default.
%
%   Output image is of class uint8 if the number of colors is not
%   greater than 256; otherwise it is double.
%
%   See also RGB2IND, CMUNIQUE.

%   Clay M. Thompson 2-25-93
%   updated by Chris Griffin 8-28-96
%   Copyright (c) 1993 by The MathWorks, Inc.
%   $Revision: 5.5 $  $Date: 1996/10/23 19:28:53 $

if nargin<2, error('Not enough input arguments.'); end

if nargin==2, % imapprox(a,cm) (remove colors)
  if size(cm,1)<=256,
    c = a; map = cm;
  else
    n = size(cm,1)-256; 
    ndx = floor((1:n)/(n+1)*size(cm,1)+1);
    map = cm; map(ndx,:) = [];
    rgb = ind2rgb(a,cm);
    c = dither(rgb,map);
  end

elseif isstr(n), % imapprox(a,cm,'dither_option')
  flag = [lower(n) ' ']; % Protect against short string
  if size(cm,1)<=256,
    c = a; map = cm;
  elseif flag(1)=='n', % 'nodither'
    n = size(cm,1)-256; 
    ndx = floor((1:n)/(n+1)*size(cm,1)+1);
    map = cm; map(ndx,:) = [];

    % Manhattan Distance to colors in reduced colormap
    M = size(map,1); N = length(ndx);
    [v,loc] = min(abs(repmat(map(:,1),1,N)-repmat(cm(ndx,1)',M,1)) + ...
                  abs(repmat(map(:,2),1,N)-repmat(cm(ndx,2)',M,1)) + ...
                  abs(repmat(map(:,3),1,N)-repmat(cm(ndx,3)',M,1)));
    dpos = ones(size(cm,1),1); 
    dpos(ndx) = zeros(size(ndx));
    pos = cumsum(dpos); pos(ndx) = loc;
    c = zeros(size(a));
    if isa(a, 'uint8')   % promotion to double is necessary until we have
        a = double(a)+1; % the utility function for uint8 indexing.
    end
    c(:) = pos(a);
  else
    n = size(cm,1)-256; 
    ndx = floor((1:n)/(n+1)*size(cm,1)+1);
    map = cm; map(ndx,:) = [];
    rgb = ind2rgb(a,cm);
    c = dither(rgb,map);
  end

elseif prod(size(n))==1, % imapprox(a,cm,tol) or imapprox(a,cm,n)
  if nargin<4, flag = 'dither'; end
  if n>size(cm,1),
    c = a; map = cm;
  else
    % Convert colormap to RGB so that rgb2ind can do all the work
    rgb(:,:,1) = cm(:,1); rgb(:,:,2) = cm(:,2); rgb(:,:,3) = cm(:,3);

    [ndx,map] = rgb2ind(rgb,n,'n');
    if flag(1)=='n', % 'nodither'
      c = zeros(size(a));
      if isa(a, 'uint8')   % promotion to double is necessary until we have
          a = double(a)+1; % the utility function for uint8 indexing.
      end
      c(:) = ndx(a);
    else
      rgb = ind2rgb(a,cm);
      c = dither(rgb,map);
    end
  end

elseif size(n,2)==3, % imapprox(a,cm,map)
  if nargin<4, flag = 'dither'; end
  flag = [lower(flag) ' ']; % Protect against short string
  map = n;
  if flag(1)=='n', % 'nodither'
    % Manhattan Distance to colors in new colormap
    M = size(map,1); N = size(cm,1);
    
    [v,loc] = min(abs(repmat(map(:,1),1,N)-repmat(cm(:,1)',M,1)) + ...
                  abs(repmat(map(:,2),1,N)-repmat(cm(:,2)',M,1)) + ...
                  abs(repmat(map(:,3),1,N)-repmat(cm(:,3)',M,1)));
    c = zeros(size(a));
    if isa(a, 'uint8')   % promotion to double is necessary until we have
        a = double(a)+1; % the utility function for uint8 indexing.
    end
    c(:) = loc(a);
  else
    rgb = ind2rgb(a,cm);
    c = dither(rgb,map);
  end    
end

if nargout==0,
  imshow(c,map)
  return
end

%if isa(c, 'double') & size(map,1) <= 256
%    c = uint8(c-1);
%end

cout = c;






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 1997-1999. The Regents of the University
% of California (Regents). All Rights Reserved.
% 
% Permission to use, copy, modify, and distribute this software and its
% documentation for educational, research, and not-for-profit purposes,
% without fee and without a signed licensing agreement, is hereby
% granted, provided that the above copyright notice, this paragraph and
% the following two paragraphs appear in all copies, modifications, and
% distributions. Contact The Office of Technology Licensing, UC Berkeley,
% 2150 Shattuck Avenue, Suite 510, Berkeley, CA 94720-1620, (510)
% 643-7201, for commercial licensing opportunities. Created by
% Chad Carson, Serge Belongie, and Jitendra Malik, Department of Electrical
% Engineering and Computer Sciences, University of California, Berkeley.
% For technical information, contact carson@eecs.berkeley.edu.
% 
% IN NO EVENT SHALL REGENTS BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT,
% SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS,
% ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF
% REGENTS HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
% 
% REGENTS SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
% PARTICULAR PURPOSE. THE SOFTWARE AND ACCOMPANYING DOCUMENTATION, IF
% ANY, PROVIDED HEREUNDER IS PROVIDED "AS IS". REGENTS HAS NO OBLIGATION
% TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR
% MODIFICATIONS.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

