function a = mrto8bit(a)
% a = mrto8bit(a)
%
% Intelligently convert 16-bit integer image A into 256 levels (uint8).
%
% Dynamic range compression is similar to a histogram equalization,
% except that no integer level in the input image is ever spread out
% to take up more than one output level.  Since most MRI pixels have
% intensity less than 100, the 256-level output image looks visually
% similar to the input image.
%
% Disadvantage: the function is rather slow.
%
% Mark Hasegawa-Johnson, 4/6/2002

if prod(size(a))>0,
   % First, convert a into integer values if it is not already
   a = round(a);
   
   % Get useful statistics about image: histogram, npix
   HIST = hist(a(:),min(0,min(a(:))):max(a(:)));
   NPIX = length(a(:));

   % Image dynamic range compression.
   % NLU of the input levels map to single output levels.
   % These levels account for NLP pixels.
   % Therefore the set of input levels that map to single output levels
   % consists of the levels with more than (NPIX-NLP)/(256-NLU) pixels.
   NLU = 0;
   SINGLEOUTPUTS = find(HIST > NPIX/256);
   while length(SINGLEOUTPUTS)>NLU,
      NLU=length(SINGLEOUTPUTS);
      SINGLEOUTPUTS = find(HIST > (NPIX-sum(HIST(SINGLEOUTPUTS)))/(256-NLU));
   end
   
  % Number of output levels per input level is min(1,inputpixperlevel/outputpixperlevel)
  OUTPUTPIXPERLEVEL = (NPIX-sum(HIST(SINGLEOUTPUTS)))/(256-NLU);
  OUTPUTSPERINPUT = min(1,HIST ./ OUTPUTPIXPERLEVEL);
  
  % Map input image levels to output image levels using a CDF
  IOMAP = max(0, min(255, ceil(cumsum(OUTPUTSPERINPUT)) - 1));
  a = IOMAP(a+1+min(0,min(a(:))));
end

