From: awpaeth@watpix.uwaterloo.ca (Alan Wm Paeth)
Subject: Lossless (Re: Best) compression for 24 bit images
Organization: University of Waterloo
Date: Mon, 31 Aug 1992 22:10:20 GMT

>In article <1992Aug28.133552.16085@crd.ge.com>
>davidsen@ariel.crd.GE.COM (william E Davidsen) writes:
> 
What is the best way to do lossless compression of 24 bit images?

A lossless technique accommodating arbitrary precision without the use of
tables is described in Graphics Gems II. The method "whitens" 2D images so that
their sequential 1D raster representations are more amenable to compression by
tools such as Unix Compress or MSDOS ZIP.  It is intentionally nonadaptive, as
that work is left to the conventional (1D) coder which follows downstream. The
method is based on predictive-corrective methods and may be regarded as an
extension to the "use the pixel 'above' on the previous scanline" technique.

In a nutshell, given a raster presented in standard CRT-scan order then the
pixel in question (*) has three adjacent neighbors, a, b and c which have
already been seen:


     a  b
     
     c (*).
     
A trial prediction employs Robert's gradient finder to estimate the slope of
the NW/SE diagonal:

    slope = b + c - 2a

and adds this to the pixel at "a":

    trial = b + c - a 

which forms an exact estimate should all four pixels lie on a planar "wedge"
of intensity (linear slope in two dimensions).

The "twist" is the non-linear final step forming an estimate by selecting that
pixel from the trial set {a, b, c} having a value closest (ie, which minimizes
the absolute value) to the trial prediction. Since two pixels of dissimilar
value may be the closest to the trial estimate (ie, equidistant), a tie-break
order is required. It is a, b, c. The code for the last step resembles:

   if |trial - a| < |trial - b| and |trial - a| < |trial - c| return(a)
   if |trial - b| < |trial - c| return(b)
   return(c)

(The Gem pseudocode is more concise; the concluding C-code more efficient).

The motivation in terms of symmetry of design is treated at length in the text;
this non-linear step yields a number of benefits. Among them, it allows the
algorithm to hold up well when presented 1-bit data -- the final estimate is
never fractional. The estimates also match the "best" choices predicted by
empirical analysis of 1-bit data. The method also supports a fast encoding
given arbitrary precision, in much the same way that Median filtering of input
does. (BTW, Median filtering and many varients was also studied in the course
of deriving this filter; it could not be made to work as well). Also, forming a
final choice from the input seen guarentees the estimate lies within the right
domain. This means that any pixel-wise "range clipping" of the estimate (ie,
when negative) step may be removed from the inner loop of the algorithm. More
important, if the "colours" presented to the coder are "sparse" (e.g., only a
few levels occuring in any subregion, despite 24-bit precision), the coder
correctly locates estimates from within the "sparse" set, unlike most other
techniques. This is a further win should the filter encounter synthetically
rendered images.

Lest I steal my own thunder, the Gems article also describes requirements of
a suitable non-adaptive 2D coder and derives the algorithm from this. The
entry also includes empirical data, histograms and sample images, plus some
pseudocode, which models the simple equations listed above.

For those who enjoy hacking inner loops of image processing tools, the fastest
implementation I've been able to create is for use with Waterloo's IM Raster
Toolkit. (In fact, the method was invented with this code suite in mind: a
lossless method for use with composable Unix-based tools accommodating
arbitrary bit precision). I've minimized the number of variables employed to
allow use of the C-language "register" compiler directive for further
optimization given integer pixels of arbitrary precision. The variables and
then the code snipped follow the signature.

   /Alan Paeth VE3AWP/KD3XG
   Computer Graphics Laboratory,
   University of Waterloo


Prefacing Notes:

b1  -- a pointer to the previous scanline of (integer) input values
b2  -- a pointer to the current  "	"
ob  -- a pointer to the output scanline of identical precision
 t  -- temporary variable
col -- number of pixels in scanline arrays b1, b2 and ob
 g  -- estimated gradient
d1  -- trial distances to pixel a
d2  -- trial distances to pixel b
d3  -- trial distances to pixel c

 p  -- the current pixel (relates to the "correction" step)
mask-- a bit-use mask (relates to reuse of the tool for "decoding" of data)

-----
...

/*
 * We consider (implicitly) the equations:
 *
 * g  = (b+c-a) = Roberts' gradient
 * d1 = |g - a| = |b+c-2a|
 * d2 = |g - b| = |c-a|
 * d3 = |g - c| = |b-a|
 *
 * then pick the smallest of d1,d2,d3 in that order.
 */
	    while(--col>=0)
		{
		a = *b1++;
		c = *b2++;		/* for "b" use "*b1" */
		if ((m = *b1 + c - a - a) < 0) m = -m;	/* m = |b+c-2a| */
		p = a;
		if ((t = a - c) < 0) t = -t;		/* t = |a-c| */
		if (t < m)
		    {
		    m = t;
		    p = *b1;
		    }
		if ((t = *b1 - a) < 0) t = -t;		/* t = |b-a| */
		if (t < m) p = c;
		*ob++ = (p - *ib++) & mask;
		}
	    }
	writeline();
	}
    imclose();
    exit(0);
    }
Newsgroups: comp.compression