Subject: Astronomical image compression

From: jwright at cfht.hawaii.edu (Jim Wright)
Date: Wed, 22 Jan 1992 03:31:53 GMT

This article is to announce the first release of the program Compfits.
Compfits was developed at Canada-France-Hawaii Telescope Corp. for
lossless (or lossy) compression of 16 bit astronomical images
(FITS or non-FITS format).
Compfits is a public domain program and is available via anonymous
ftp from uwila.cfht.hawaii.edu (128.171.80.50).
To get it, log in as anonymous, retrieve the compressed tar file
compfits.tar.Z in the pub/compfits directory and untar it with
the command  "zcat compfits.tar.Z | tar xf -".
What follows is the manpage of Compfits. Just cut after the "------"
and use the command nroff -man to read it.
Feel free to test this program and to mail any comments, criticism or
suggestions to compfits at cfht.hawaii.edu .

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

----------------------------------------------------------------------

.         at (#)Compfits.1	1.0 01/14/92
.TH COMPFITS 1 "January 1992" "Public Domain"
.SH NAME
compfits \- compress/uncompress 16 bit images (lossless or lossy)
.SH SYNOPSIS
.SH Compression:
.PP
.HP 5
.B compfits
[\ \fB\-cinqvwL\ \fR] 
\ \fB\-F\ \fIorigfile\ \fR
[\ \fB\-f\ \fIcompfile\ \fR]
[\ \fB\-P\ \fIpathext\ \fR]
[\ \fB\-p\ \fIargext\ \fR]
[\ \fB\-N\ \fInuser\ \fR]
[\ \fB\-e\ \fIepsilon\ \fR]
[\ \fB\-m\ \fImsbbuf\ \fR]
[\ \fB\-l\ \fIlsbbuf\ \fR]
[\ \fB\-t\ \fItesbuf\ \fR]
.PP
.SH Decompression:
.PP
.HP 5
.B compfits
\ \fB\-d\ \fR
[\ \fB\-inqv\ \fR]
\ \fB\-f\ \fIcompfile\ \fR
[\ \fB\-F\ \fIdecfile\ \fR]
[\ \fB\-P\ \fIpathext\ \fR]
[\ \fB\-p\ \fIargext\ \fR]
[\ \fB\-m\ \fImsbbuf\ \fR]
[\ \fB\-l\ \fIlsbbuf\ \fR]
[\ \fB\-t\ \fItesbuf\ \fR]
.SH DESCRIPTION
\fICompfits\fP 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.
.PP
\fICompfits\fP doesn't compress anything by itself and therefore needs an
external compression program to work.
\fICompfits\fP 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
\fBLSB image\fP (containing the n least significant bits of the original
image) and the \fBMSB image\fP (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.
\fICompfits\fP 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.
\fICompfits\fP 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.
\fICompfits\fP makes a guess for the optimal n trying to estimate
the entropy of the different bit planes, unless the user specifies
a partition.
.PP
The compressed file has the following structure:
.nf
.I HEADER
.I LSB image
.I Compressed MSB image
.sp
.fi
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.
.PP
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 \fBCOMPART\fP.
.PP
If the original file is a FITS file the value of the \fBCOMPART\fP card
is a long integer:
\fIn+17*size_diff\fP where \fIsize_diff\fP is the difference in bytes between
the actual size of the file and the expected size according to the
FITS header (\fIsize_diff\fP should normally be zero).
.PP
If the original file is not a FITS file the value of the \fBCOMPART\fP card
is a long integer:
\fIn+17*size_file\fP where \fIsize_file\fP is the size in bytes of
the original file.
.PP
Another card named \fBCOMPRESS\fP is added to indicate that the file is
compressed and to give the name of the compression algorithm.
When \fICompfits\fP is used, the value of the \fBCOMPRESS\fP card is set to
\fBCOMPFITS\fP for a lossless compression and to \fBLOSSFITS\fP for
a lossy compression.
When it is not a FITS image, the added header contains both the
\fBCOMPRESS\fP and the \fBCOMPART\fP card following a \fBSIMPLE\fP card
with the logical value \fB'F'\fP (false), to denote that although
the header may be recognized as FITS, the data does not conform
to the FITS standard.
.PP
Because it needs to first analyse the data before starting the
compression process, \fICompfits\fP 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, \fICompfits\fP needs
an output disk file.
In other words, \fICompfits\fP cannot be used as a filter.
The recommended suffix for files compressed with \fICompfits\fP is \fB.CF\fP.
.SH OPTIONS
All numbers indicating sizes may be suffixed by \fBb\fP (512 bytes blocks),
\fBk\fP (kilobytes), or \fBm\fP (megabytes).
.TP 12
.B \-c
Compress data.
.TP
.B \-d
Decompress data.
.TP
.BI \-e \ epsilon
Use this value of \fIepsilon\fP instead of the defined one.
.TP
.BI \-f \ filename
Name of the compressed file.
For compression, the default is \fIoriginal\fP\fB.CF\fP.
For decompression, there is no default; this option must be given.
.TP
.BI \-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 \fIcompressed\fP\fB.REC\fP.
.TP
.B \-i
Ignore garbage at the end of file if any (this may happen if the
file is blocked for storage).
.TP
.B \-L
Lossy compression: the LSB image is not written which allows the
program to achieve particularly good compression ratios.
Reconstruction is not exact.
.TP
.BI \-l \ lsbbuf
Size of the buffer for the LSB image.
Data is written to or read from the disk
in blocks of this size.
.TP
.BI \-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.
.TP
.BI \-N \ nuser
Force the partion value \fIn\fP to be \fInuser\fP and not the one determined
by the program.
.TP
.B \-n
Indicates that the file to be compressed or uncompressed is not a FITS file.
.TP
.BI \-P \ pathext
Path to the external compression or decompression program.
Default is \fI/usr/local/bin/compact_1.0\fP for compression.
Default is \fI/usr/local/bin/uncompact_1.0\fP for decompression.
.TP
.BI \-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 "\fB-q\fP" (quiet mode) for compact.
.TP
.B \-q
Quiet mode: no compression or decompression statistics will be
printed on stderr.
.TP
.BI \-t \ tesbuf
Size of the central test zone for determining the best partition.
Default is 512*512*16 bits.
.TP
.B \-v
Verbose: prints some information on the image to stderr.
.TP
.B \-w
Print the entropy of all the bit planes.
.SH THEORY OF OPERATION
For \fICompfits\fP, the critical operation is to determine the
optimal partition of LSB/MSB (i.e. determine the optimal \fIn\fP).
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 \fIi\fP, 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 \fIi-entropy(i)>epsilon\fP.
Thus the optimal \fIn\fP is chosen such that \fIn-entropy(n)<epsilon\fP and
\fI(n+1)-entropy(n+1)>epsilon\fP.
This comes simply from the fact that the low order bit planes are
supposed to contain almost only noise and therefore \fIentropy(i)\fP is
very close to \fIi\fP, whereas high order bit planes mainly contain
information so \fIentropy(i)\fP is much lower than \fIi\fP.
With \fICompfits\fP plus \fICompact 1.0\fP, the experimental best value of
\fIepsilon\fP is 0.01 (used as the default \fIepsilon\fP).
Usually the function \fIf(i)=i-entropy(i)\fP is non-decreasing with a
sharp transition for \fIi=n\fP, 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.
.SH PERFORMANCE
The compression ratio achieved depends of course on the external
compression program.
As \fICompfits\fP 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 \fICompact 1.0\fP,
using the Lempel-Ziv algorithm.
The default value of \fIepsilon\fP is set for use with this program.
.PP
\fICompfits\fP plus \fICompact 1.0\fP
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
\fIUnix Compress\fP or \fICompact 1.0\fP 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 \fICompfits\fP
plus \fICompact 1.0\fP.
.PP
On these astronomical images, and with the default value of \fIepsilon\fP,
the partition found by \fICompfits\fP is almost always optimal; any
other partition leads to worse results.
Note that a partition with \fIn=0\fP (0 LSB) is the same as using the
external program alone.
The more noise there is in the image, the higher \fIn\fP is and the greater
the difference between the use of \fICompfits\fP and the simple use of the
external program is.
.PP
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.
.PP
.B WARNING:
This progam will work well on any noisy 16 bits/word file, but
may yield poor results on any other kind of file.
.SH DIAGNOSTIC
\fICompfits\fP prints error messages and warning messages on stderr.
It differentiates between "error due to operating system", typically
producing an \fIerrno\fP code and "error due to program", (all the others).
.PP
.nf
\fICompfits\fP exits with a 0 status on normal completion.
.sp
The exit status is negative on error:
-1 for a system error;
-2 for a program error.
.sp
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).
.SH EXAMPLES
.B compfits -F image
will compress the file image, using \fI/usr/local/bin/compact_1.0\fP as
the external compression program and with all the default parameters.
The compressed file will be \fIimage.CF\fP.
.PP
.B compfits -d -f image.CF
will do the reverse, uncompressing \fIimage.CF\fP using
\fI/usr/local/bin/uncompact_1.0\fP as the external program.
The reconstructed file will be \fIimage.CF.REC\fP and will be
identical to \fIimage\fP.
.PP
.B compfits -F image -f compimage
.B -P /usr/bin/compress -p -b16
will compress \fIimage\fP to \fIcompimage\fP, using \fI/usr/bin/compress\fP
with options \fI-b16\fP as the external compression program.
All other parameters will be default ones.
.PP
.B compfits -d -f compimage -F recimage
.B -P /usr/bin/uncompress -p -b16
will do the reverse.
\fIrecimage\fP and \fIimage\fP will be identical.
.PP
.B compfits -F image -N 4 -w
(reverse:
.B compfits -d -f image.CF)
will compress \fIimage\fP to \fIimage.CF\fP 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 \fIepsilon\fP is optimal.
.PP
.B compfits -F image -L
(reverse:
.B compfits -d -f image.CF)
will compress \fIimage\fP to \fIimage.CF\fP.
The compression will be lossy.
Therefore after the decompression, \fIimage\fP and \fIimage.CF.REC\fP
won't be identical.
.SH AUTHOR
.nf
Jean-Pierre Veran
Canada-France-Hawaii Telescope Corporation /
   Ecole Nationale Superieure des Telecommunications de Paris
.sp
Copyright (C) 1992    Canada-France-Hawaii Telescope Corp.
.fi
.PP
All the programs in the compfits package are distributed WITHOUT any
warranty, and are under the terms of the GNU General Public License.
.PP
.nf
Please, send any questions or comments to
"compfits at cfht.hawaii.edu".
.fi
--
Jim Wright
jwright at cfht.hawaii.edu
Canada-France-Hawaii Telescope Corp.