Subject: FITS File Compression / CCD Image Compression

From: jim at galileo.ifa.hawaii.edu (Bubba Mersault)
Date: Tue, 21 Apr 1992 08:56:31 GMT

I am searching for suggestions on FITS file compression.  I
have experimented with the standard UNIX compress utility
and found it not to be very effective. I have also used W.
Press's FITSPRESS program, which is effective on some
images but slow (19 minutes to compress a 2048x2048 image
on a SPARC Server) for my purposes. It is, however,
ineffective on on some images I have with only a couple of
bright stars in the field.

My needs dictate that a lossless compression system be
used.  Also, I am dealing with large format CCDs -
2048x2048 now, and larger ones in the near future.

If anyone could suggest other software, or even a
hardware solution, that I could evaluate, I would
appreciate it.

Kevin T.C. Jim            Voice: (808) 956-9867
Institute for Astronomy   FAX:   (808) 988-2790        
2680 Woodlawn Drive       Email: jim at galileo.ifa.hawaii.edu

---

From: dwells at azalea.cx.nrao.edu (Don Wells)
Date: Tue, 21 Apr 1992 14:07:26 GMT

In article <1992Apr21.085631.20342 at news.Hawaii.Edu>
jim at galileo.ifa.hawaii.edu (Bubba Mersault) writes:

 BM> Summary: Request for info on image compression software
 BM> Keywords: FITS,compression,CCD

 BM> "I am searching for suggestions on FITS file compression.. [have]
 BM> used W. Press's FITSPRESS program, which is effective on some
 BM> images but slow.. [need] a lossless compression system.. CCDs..
 BM> 2048x2048..."  

This subject was discussed in the alt.sci.astro.fits newsgroup (also
"fitsbits" exploder) in September-91 and January-92. The traffic is in
the files:

fits.cx.nrao.edu[192.33.115.8]:/FITS/fitsbits/fitsbits.{9109,9110,9201}

I think that the discussion in the 9109 and 9110 files is quite
complete; many, many people gave their opinions and ideas.  Take a
look at the announcement of program COMPFITS in the 9201 file.

			-*-

You will find my opinions in several of the messages in the archives,
but I will take this opportunity to state them once more:

\begin{soapbox}

Truly random noise is {\em incompressible}. Read that sentence
again---it is fundamental! If any portion of your bitstream is random
noise, then that portion is incompressible.

Suppose you have sky-limited images with 16-bit pixels and sky level
around 1000 (out of the 65536 total range) with negligible amplifier
noise and dark value. Almost all of your pixels have approximately 10
non-zero bits (1000 ~= 2^10). The Poisson noise is sqrt(value), so the
bottom 5 bits of each pixel are random noise. The top six bits are
zero, so they will compress with high efficiency. The middle five bits
are approximately constant (for a flattened image or a chip of good
cosmetic quality) and so they too will compress well. The bottom five
bits will not compress. In summary, this trivial back-of-the-envelope
calculation tells us that an optimum non-destructive algorithm can
reduce such an image from 16 bits/pixel to ~5 bits/pixel. This is a
size reduction by a factor of ~3.2.

Now you should ask yourself how many of those 5 bits you really need.
I contend that for even the most demanding applications (e.g.
detecting low surface brightness in the outermost parts of galaxies)
you only need ~3 of them. Indeed, ~2 of them might be
enough.\footnote{Use Monte Carlo trials of the demanding
application(s) to decide the question.} THIS MEANS THAT ANOTHER FACTOR
OF {\em TWO} IS AVAILABLE IF YOU ARE WILLING TO MASK EXTRA NOISE BITS!

But wait---there's more. Probably the direct imagery is sbomewhat
oversampled. By that I mean that probably the camera samples more than
2 pixels per seeing FWHM (this is the recommended practice in
imaging).  I.e., the power spectrum of the image probably extends
beyond the power spectrum of the signal that came through the
atmosphere and imaging system. (The cutoff of the signal must be
evaluated in those portions of the image which have interesting signal
composed of high spatial frequencies which the atmospheric turbulence
has limited.) If a Wiener noise-limiting filter is applied to such an
image the redundant noise at high spatial frequencies will be
eliminated but the valid signal will remain. The effects on the
signals produced by the linear filter operator are precisely
calculable. Suppose that the effective cutoff of the signal is about
2/3 of the Nyquist frequency.  THIS MEANS THAT ANOTHER FACTOR OF 1.5
IS AVAILABLE IF YOU ARE WILLING TO APPLY A LINEAR NOISE-LIMITING
FILTER!\footnote{I believe that the reduction is actually 1.5^2
because the filter is applied in the y-direction too, but I will
ignore this additional 1.5x for this discussion.}

The two techniques I have just advocated are capable of making big
gains in compression ratio without destroying the real signals in the
data. What does the term "lossless" mean in light of this argument?
If you are only a programmer and know nothing about the {\em meaning}
of the bits, you will say that the specification "lossless" means that
the de-compressed bitstream must be {\em identical} to the original
bitstream. But if you are aware that the atmosphere, imaging system
and random noise processes have already corrupted the signal in the
original bitstream, you can decide that the specification "lossless"
will have a different meaning in your compression application.

For the example numbers (nbits=16, sky=1000, cutoff=0.67), the basic
gain will eliminate about 2/3 of the bits, about 1/2 of the remaining
1/3 of the bits can be removed because they over-sample the noise, and
finally about 1/3 of the remaining 1/6 of the bits can be removed
because they over-sample the signal. IF ALL THREE TECHNIQUES ARE USED
THE IMAGE CAN BE REPRESENTED BY ABOUT 1.7 BITS/PIXEL, OR A COMPRESSION
RATIO OF ABOUT 9.5.\footnote{If the extra 1.5x Wiener filter factor is
available then the numbers are ~1.1 bits/pixels, compression ~14x.}

If compression really matters to you and your users it will be
productive to think long and hard about using these techniques. In
particular, masking the redundant (Poisson) noise bits is cheap and
easy to implement and will make a {\em big} gain in many cases.

\end{soapbox}

--

Donald C. Wells             Associate Scientist        dwells at nrao.edu
National Radio Astronomy Observatory                   +1-804-296-0277
520 Edgemont Road                                 Fax= +1-804-296-0278
Charlottesville, Virginia 22903-2475 USA            78:31.1W, 38:02.2N

---

From: steve at cfht.hawaii.edu (Steven S Smith)
Date: Tue, 21 Apr 1992 20:00:55 GMT

>I am searching for suggestions on FITS file compression.  I
...

FYI, I have enclosed the man page for some code we have here.
I am not the author, but in the spirit of aloha...

enjoy,
steve

 COMPFITS(1)			Public Domain			 COMPFITS(1)
				January 1992

 NAME
      compfits - compress/uncompress 16 bit images (lossless or lossy)

 SYNOPSIS
 Compression:
      compfits [ -cinqvwL ]  -F origfile  [ -f compfile ] [ -P pathext ]
	   [ -p argext ] [ -N nuser ] [ -e epsilon ] [ -m msbbuf ]
	   [ -l lsbbuf ] [ -t tesbuf ]

 Decompression:
      compfits	-d  [ -inqv ]  -f compfile  [ -F decfile ] [ -P pathext ]
	   [ -p argext ] [ -m msbbuf ] [ -l lsbbuf ] [ -t tesbuf ]

 DESCRIPTION
      Compfits is a utility program developed at Canada-France-Hawaii
      Telescope Corporation for lossless compression of astronomical 16
      bits/pixel FITS (Flexible Image Transport System) format images.
      These images are usually very noisy and don't shrink well with the
      typical lossless compression algorithms.	Options are available to
      compress non-FITS files and to perform lossy compression.

      Compfits doesn't compress anything by itself and therefore needs an
      external compression program to work.  Compfits analyzes the bit
      planes of the original image to see which contain only noise (the low
      order planes) and which contain the actual information (the high order
      planes).	These two sets of bit planes, or partitions, are
      respectively called the LSB image (containing the n least significant
      bits of the original image) and the MSB image (containing the 16-n
      most significant bits of the original image).  The LSB image is
      supposed to contain only noise, therefore any attempt to compress it
      would actually increase its size.  Compfits writes it as it is in the
      compressed file, just packing the bits one after the other.  The MSB
      image contains the information and can usually compress fairly well.
      Compfits sends it to the chosen external compression program which
      must be a filter (i.e. read stdin and send compressed version to
      stdout) and writes the compressed MSB image in the compressed file.
      Of course, the best partition (i.e. the best value of n) depends a
      good deal on the image and slightly on the external compression
      program used.  Compfits makes a guess for the optimal n trying to
      estimate the entropy of the different bit planes, unless the user
      specifies a partition.

      The compressed file has the following structure:
      HEADER
      LSB image
      Compressed MSB image

      The header is the FITS header if the original is a FITS format image
      or a blank 2880 bytes FITS header if it is not.

				    - 1 -	  Formatted:  April 21, 1992

 COMPFITS(1)			Public Domain			 COMPFITS(1)
				January 1992

      For the decompression, it is necessary to know which partition was
      used during the compression.  This information is written in the
      header as a card named COMPART.

      If the original file is a FITS file the value of the COMPART card is a
      long integer: n+17*size_diff where size_diff is the difference in
      bytes between the actual size of the file and the expected size
      according to the FITS header (size_diff should normally be zero).

      If the original file is not a FITS file the value of the COMPART card
      is a long integer: n+17*size_file where size_file is the size in bytes
      of the original file.

      Another card named COMPRESS is added to indicate that the file is
      compressed and to give the name of the compression algorithm.  When
      Compfits is used, the value of the COMPRESS card is set to COMPFITS
      for a lossless compression and to LOSSFITS for a lossy compression.
      When it is not a FITS image, the added header contains both the
      COMPRESS and the COMPART card following a SIMPLE card with the logical
      value 'F' (false), to denote that although the header may be
      recognized as FITS, the data does not conform to the FITS standard.

      Because it needs to first analyse the data before starting the
      compression process, Compfits needs an input disk file.  Because it
      needs to write the compressed MSB and LSB images at the same time to
      different locations in the compressed file, Compfits needs an output
      disk file.  In other words, Compfits cannot be used as a filter.	The
      recommended suffix for files compressed with Compfits is .CF.

 OPTIONS
      All numbers indicating sizes may be suffixed by b (512 bytes blocks),
      k (kilobytes), or m (megabytes).

      -c	  Compress data.

      -d	  Decompress data.

      -e epsilon  Use this value of epsilon instead of the defined one.

      -f filename Name of the compressed file.	For compression, the default
		  is original.CF.  For decompression, there is no default;
		  this option must be given.

      -F filename Name of the original or reconstructed file.  For
		  compression, there is no default; this option must be
		  given.  For decompression, the default is compressed.REC.

      -i	  Ignore garbage at the end of file if any (this may happen
		  if the file is blocked for storage).

				    - 2 -	  Formatted:  April 21, 1992

 COMPFITS(1)			Public Domain			 COMPFITS(1)
				January 1992

      -L	  Lossy compression: the LSB image is not written which
		  allows the program to achieve particularly good
		  compression ratios.  Reconstruction is not exact.

      -l lsbbuf   Size of the buffer for the LSB image.  Data is written to
		  or read from the disk in blocks of this size.

      -m msbbuf   Size of the buffer for the MSB image.  Data is sent to or
		  received from the external compression program in blocks
		  of this size.

      -N nuser	  Force the partion value n to be nuser and not the one
		  determined by the program.

      -n	  Indicates that the file to be compressed or uncompressed
		  is not a FITS file.

      -P pathext  Path to the external compression or decompression program.
		  Default is /usr/local/bin/compact_1.0 for compression.
		  Default is /usr/local/bin/uncompact_1.0 for decompression.

      -p argext   Arguments of the external compression or decompression
		  program.  If multiple arguments with imbedded spaces are
		  present, this must be a quoted string.  Default is "-q"
		  (quiet mode) for compact.

      -q	  Quiet mode: no compression or decompression statistics
		  will be printed on stderr.

      -t tesbuf   Size of the central test zone for determining the best
		  partition.  Default is 512*512*16 bits.

      -v	  Verbose: prints some information on the image to stderr.

      -w	  Print the entropy of all the bit planes.

 THEORY OF OPERATION
      For Compfits, the critical operation is to determine the optimal
      partition of LSB/MSB (i.e. determine the optimal n).  To do so, it
      tries to compute the entropy of the different bit planes of a
      subimage.  For homogeneous 16 bits/pixel images, experiments show that
      there is no advantage in sampling a subimage bigger than 512*512.
      First the image is translated from the binary code into the Gray code
      to give more coherence to the bit planes.  Then for each bit plane i,
      starting with the zeroth order and increasing, the entropy is computed
      assuming that the bit plane is a binary discrete memoryless source.
      This continues until i-entropy(i)>epsilon.  Thus the optimal n is
      chosen such that n-entropy(n)<epsilon and (n+1)-entropy(n+1)>epsilon.
      This comes simply from the fact that the low order bit planes are
      supposed to contain almost only noise and therefore entropy(i) is very
      close to i, whereas high order bit planes mainly contain information

				    - 3 -	  Formatted:  April 21, 1992

 COMPFITS(1)			Public Domain			 COMPFITS(1)
				January 1992

      so entropy(i) is much lower than i.  With Compfits plus Compact 1.0,
      the experimental best value of epsilon is 0.01 (used as the default
      epsilon).  Usually the function f(i)=i-entropy(i) is non-decreasing
      with a sharp transition for i=n, which makes the determination of the
      optimal partition quite reliable.  Of course, considering each bit
      plane as a binary discrete memoryless source is a bit simplistic.
      However, as the computation is actually performed only on the low
      order bit planes (i.e. almost white noise most of the time), this
      model is quite accurate.

 PERFORMANCE
      The compression ratio achieved depends of course on the external
      compression program.  As Compfits assumes a 16 bits/pixel image, the
      best results are obtained when using an external compression program
      which assumes a 16 bit word input stream too.  An example of such a
      program is the public domain Compact 1.0, using the Lempel-Ziv
      algorithm.  The default value of epsilon is set for use with this
      program.

      Compfits plus Compact 1.0 were tested on a number of 8 megabyte FITS
      images (2048*2048 pixels, 16 bits/pixels) taken at the Canada-France-
      Hawaii Telescope.  Typically an object image (i.e. an actual sky image
      rather than a flat field or bias image) is compressed to 4 megabytes
      whereas usual lossless compression programs such as Unix Compress or
      Compact 1.0 rarely do better than 6 megabytes.  The difference is even
      more drammatic when the images are very noisy (e.g. flat field
      images).	Such images may not shrink at all with a usual lossless
      compression program whereas it is reduced to 6 megabytes with Compfits
      plus Compact 1.0.

      On these astronomical images, and with the default value of epsilon,
      the partition found by Compfits is almost always optimal; any other
      partition leads to worse results.  Note that a partition with n=0 (0
      LSB) is the same as using the external program alone.  The more noise
      there is in the image, the higher n is and the greater the difference
      between the use of Compfits and the simple use of the external program
      is.

      The visual result of lossy compression is usually quite good as the
      LSB image which is removed is supposed to contain only noise.
      Moreover, 256 color workstation displays usually can't show the whole
      range of intensity which hides the defects.  However a lossy
      compressed astronomical image can't possibly be considered as a
      scientific frame.

      WARNING: This progam will work well on any noisy 16 bits/word file,
      but may yield poor results on any other kind of file.

 DIAGNOSTIC
      Compfits prints error messages and warning messages on stderr.  It
      differentiates between "error due to operating system", typically

				    - 4 -	  Formatted:  April 21, 1992

 COMPFITS(1)			Public Domain			 COMPFITS(1)
				January 1992

      producing an errno code and "error due to program", (all the others).

      Compfits exits with a 0 status on normal completion.

      The exit status is negative on error:
      -1 for a system error;
      -2 for a program error.

      The exit status is positive on warning:
      1 for bad arguments;
      2 if the file is shorter than expected from examining the header;
      3 if the file larger than expected from examining the header;
      4 if the image is not 16 bits/pixel;
      10 if the image is not a FITS file (unexpected).

 EXAMPLES
      compfits -F image will compress the file image, using
      /usr/local/bin/compact_1.0 as the external compression program and
      with all the default parameters.	The compressed file will be
      image.CF.

      compfits -d -f image.CF will do the reverse, uncompressing image.CF
      using /usr/local/bin/uncompact_1.0 as the external program.  The
      reconstructed file will be image.CF.REC and will be identical to
      image.

      compfits -F image -f compimage -P /usr/bin/compress -p -b16 will
      compress image to compimage, using /usr/bin/compress with options -b16
      as the external compression program.  All other parameters will be
      default ones.

      compfits -d -f compimage -F recimage -P /usr/bin/uncompress -p -b16
      will do the reverse.  recimage and image will be identical.

      compfits -F image -N 4 -w (reverse: compfits -d -f image.CF) will
      compress image to image.CF forcing the partition to be 4 LSB / 12 MSB,
      and displaying the entropy of the different bit planes.  This is
      useful when testing if the default value of epsilon is optimal.

      compfits -F image -L (reverse: compfits -d -f image.CF) will compress
      image to image.CF.  The compression will be lossy.  Therefore after
      the decompression, image and image.CF.REC won't be identical.

 AUTHOR
      Jean-Pierre Veran
      Canada-France-Hawaii Telescope Corporation /
	 Ecole Nationale Superieure des Telecommunications de Paris

      Copyright (C) 1992    Canada-France-Hawaii Telescope Corp.

				    - 5 -	  Formatted:  April 21, 1992

 COMPFITS(1)			Public Domain			 COMPFITS(1)
				January 1992

      All the programs in the compfits package are distributed WITHOUT any
      warranty, and are under the terms of the GNU General Public License.

      Please, send any question or comment to
      "compfits at cfht.hawaii.edu".

				    - 6 -	  Formatted:  April 21, 1992

--
Steven S. Smith 		Canada-France-Hawaii Telescope Corp. 
steve at cfht.hawaii.edu 		(NeXT mail grudgingly accepted ;-)