function imwrite_compressed(source_image, filename)
% imwrite_compressed(A, FILENAME)
% writes the image A to FILENAME,
% in experimental format (using Huffman compression).
% much like imwrite.
% Assumes A is a integer-valued matrix.
% (A can have 2, 4, 8, 12, or more bits per pixel).
% (A can have any number of dimensions --
% a 2D grayscale image, a 3D grayscale movie
% a 4D color movie, etc.)
% Usage:
%	cameraman = double(imread('cameraman.tif'));
%	imwrite_compressed(cameraman, 'temp.dat');
%	X = imread_compressed('temp.dat');
%	if( isequal(cameraman,X) )
%		disp('perfect lossless reconstruction')
%	else,
%		error(' I messed up somewhere ... ')
%	end;
%
% You can get better compression by prefiltering;
%   t = cameraman - 128;
%   diff_cameraman = ...
%     [ double(t(1,:)); diff(double(t)) ];
%   imwrite_compressed(diff_cameraman+256, 'tempdiff.dat');
%   diffX = imread_compressed('tempdiff.dat') - 256;
%   X = 128 + cumsum(diffX);
%   if( isequal(cameraman,X) )
%       disp('perfect lossless reconstruction')
%   else,
%		error(' I messed up somewhere ... ')
%   end;
%
% I wish I could use this with
%    imwrite(A, FILENAME, 'dave');
%    X = imread(FILENAME, 'dave');
% See also IMREAD_COMPRESSED, HUFFMAN_COMPRESS, IMWRITE.

% Change Log
% 1999-07-19:DAV: Threw out needless complexity.
% 1999-06-26:DAV: David Cary <d.cary at ieee.org> started.

% FIXME: this is really crude. I'd rather
% base it on the "PNG" format. See
% http://www.cdrom.com/pub/png/spec/PNG-CRCAppendix.html

% FIXME: Do funky wavelet transform here.
plaintext = source_image;

% header = \211   P   N   G  \r  \n \032 \n
header = uint8(double( [ 137  80  78  71  13  10  26  10 ] ));

shape = size(plaintext); % *not* packed into uint8 !
[compressed_data, bitlengths, encoded_data_length] = huffman_compress(plaintext(:));

max_length = max(bitlengths(:))
if( max_length < 23 )
   bitlengths = uint8(bitlengths);
else,
   max_length
   error('Sorry, code lengths greater than 22 bits not yet implemented.')
end;

save( filename,  'header', 'shape', 'bitlengths', 'compressed_data' );

if(1)
   
   length(compressed_data)
   disp(' bytes of compressed data ')
   length(bitlengths)
   disp(' distinct values ')
   
   D = dir(filename);
   if( 1 == length(D) ),
      
      filesize = D(1).bytes
      % When saving just one variable of type uint8,
      % Matlab apparently prepends a (exactly) 200 byte header.
      if( length(compressed_data) + 1000 < filesize )
         warning(' Unusually long header -- over 1000 bytes. ')
      end;
      
      disp('Compressed to:')
      pixels = prod(size(source_image));
      compressed_data_bpp = encoded_data_length/pixels;
      total_bpp = (filesize*8)/pixels;
      header_bpp = total_bpp - compressed_data_bpp;
      
      header_bpp
      compressed_data_bpp
      total_bpp
      
      disp('bpp.')
      
      
      
      
      if(1)
         restored_data = imread_compressed(filename);
         if( isequal(source_image,restored_data) )
            disp('perfect lossless reconstruction')
         else,
            error(' I messed up somewhere ... ')
         end;
         
      end;
      
   else,
      ls
      D
      filename
      warning(' I can''t seem to find the file ')
      warning( filename )
      warning(' that I just saved. Perhaps Matlab appended a ".mat" extension ?')
   end;
   
end;
% end imwrite_compressed.m