Help for DESPIKE

PURPOSE:
DESPIKE identifies and removes single-pixel errors (e.g. telemetry bit errors)
from images.

EXECUTION:

   DESPIKE  INP  OUT
or DESPIKE  INP  OUT  SCALE=3  TOL=3

where

   INP and OUT are the input and output images 
   SCALE and TOL are optional scale and offset parameters for controlling the
   DN threshold.

The program converts the input to floating-point (real*4) and converts the 
results back to the format of the input.  Therefore, all VICAR formats are
accepted, but some loss of precision may occur with longword integer or
double precision formats.

REFERENCE:

  "A comparison of image despiking algorithms", Gary Yagi, IOM, June 1, 1999.

OPERATION:

DESPIKE identifies single-pixel errors by comparing the DN value of a pixel
with that of its eight surrounding neighbors.  If the pixel differs from its
neighbors by some specifiable threshold (see SCALE and TOL parameters) the
pixel is removed by interpolation.

DESPIKE will actually remove errors affecting multiple pixels as long as no
more than 3 noisy pixels are included in any 3x3 area.  Thus, small radiation
noise events and single line or column dropouts may be successfully removed.

See also programs REMNOISE and REMRAY.

METHOD:

The algorithm consists of the following steps:

  1) Use a median filter to identify all "suspicious" pixels in the 3x3 area.
  2) If the central pixel is one of the suspects, use only the reliable pixels
     in the area to determine if it is a spike.
  3) If the central pixel is determined to be a spike, replace it by fitting
     all the reliable pixels to a surface.

In the following diagram, let Do represent the DN value of the pixel being
tested, and D1, D2, D3, ..., D8 represent the DN values of the pixels in the
3x3 area surrounding it:


		D1  D2  D3			M1  M2  M3
			      median filter
		D8  Do  D4    ------------->	M8  Mo  M4

		D7  D6  D5			M7  M6  M5

A 3X3 median filter (see program MEDIAN) is applied (to the entire image),
resulting in median values Mo, M1, M2, ..., M8 for each of the nine pixels in
the area.  The median filter is used to provide an estimate of what the
"correct" value of each pixel would be if there were no noise.  A pixel Di
is considered suspicious if it differs from its median by an amount
significantly more than the mean difference in the area:

			 	 8
				---
			     S  \
		  |Di-Mi| > --- /  |Dj-Mj| + To			(1)
			     9  ---
                                j=0

where S and To are scale and offset constants specified via the SCALE and TOL
parameters.

If, after the above screening, the central pixel is one of the suspects, it is
subjected to a second test, based only on reliable pixels within the area:

			|Do - Eo| > S x A + To			(2)
where
    Eo is an estimate of the "correct" value of Do,
    A is a scene activity measure,
    S and To are scale and offset constants specified via the SCALE and TOL
      parameters.
    
The estimator Eo is computed by fitting the 3x3 DN surface to a plane.  The
scene activity measure A is the mean absolute difference between adjacent pixels
along the perimeter.  Only valid pixels are used in these calculation.  If the
central pixel is a spike, it is replaced by Eo.

USING THE SCALE AND TOL PARAMETERS:

As described above, the threshold for determining whether a pixel is a spike
is computed from some linear function of the local scene activity.  The scale
and offset terms (S and To) of this linear function are specified via the
SCALE and TOL parameters.

In general, the useful range for S is between 1.5 and 3.  The useful range for
To is between 2 and 6.  Setting To=0 or 1 is likely to cause large numbers of
valid pixels to be destroyed (this is similar to setting the lawn mower blades
too low and chewing up large chunks of root and dirt).  The following values
are recommended:

	SCALE=3  TOL=3  for bit-error rates of 10**-3 or less
	SCALE=2  TOL=2  for bit-error rates of 10**-2 or greater

Note that it is more fruitful to adjust SCALE rather than the TOL.  To get
some feel for what effect a particular bit error rate has on an image,
experiment with program ADDNOISE.

The effectiveness of the program in removing noise may be determined by
applying the following test on a noise-free image:

   addnoise clean noisy rate=1000	  !add 10**-3 BER of noise to image
   despike  noisy   d   scale=3  tol=3    !despike the noisy image
   f2 (clean,d) diff func="abs(in1-in2)"  !compare despiked and clean images
   hist  diff  'nohist			  !and print residual noise

In the test, we add noise to the clean image (Bit-Error-Rate=10**-3) and
attempt to remove the noise using despike.  By subtracting the despiked image
from the original image, we obtain a measure of the residual noise, i.e. the
noise that DESPIKE failed to remove.  By varying the SCALE and TOL parameters,
we can determine optimal values which minimize the residual noise.

The following data was obtained by applying the above test to an image of the
moon after adding 4.7 DN of noise:

	                residual noise (DN)
	      To   s=1.5    s=2    s=2.5    s=3
	       1  1.4160  1.0983  0.8762  0.8196
	       2  1.1476  0.9032  0.7785  0.7605
	       3  1.0062  0.8215  0.7463  0.7620
	       4  0.9190  0.7865 *0.7380  0.7602
	       5  0.8608  0.7674  0.7459  0.7779

In this example, the minimum occurs when SCALE=2.5 and TOL=4.  In general,
lowering the value of S or To causes more noise spikes to be detected and
removed, but causes more valid pixels to be destroyed as well.  The minumum
occurs approximately at that point where lowering either value any further
causes more valid pixels to be destroyed than noise spikes.

To answer the question of how much harm DESPIKE does to an image, we can
apply it to the clean image (assumed to be noise-free) at various combinations of
S and To.  The following data was obtained by applying DESPIKE to the Earth:

		    Valid pixels destroyed
		To  S=1.5   S=2  S=2.5  S=3
		1   7070   2161   598   232
		2   3075    769   162    48
		3   1498    329    66    18
		4    771    162    32    12
		5    394     75    24     7

Note that the above data is skewed by the fact that the original image
(assumed to be noise-free) contains radiation noise (Galileo was traveling
through the Van Allen radiation belt at the time).


PROGRAMMING NOTES:

The following comments are addressed to the programmer charged with maintaining
the program:

The median of each 3x3 area is computed by sorting the 9 pixels and selecting
the 5th lowest value.  The sort algorithm uses the fact that it is only
necessary to find the 5 lowest values to determine the median.  The resulting
algorithm is twice as fast as a conventional sort routine.

Because the filter window scans the image from left-to-right and from top-to-
bottom, at any given position of the window only the four pixels in the
lower right quadrant must be screened for suspicious pixels:

			D1  D2  D3

			D8  Do?	D4?

			D7  D6? D5?

The remaining pixels in the area have already received a more accurate
screening at previous positions of the filter window.

Note that the same SCALE and TOL values are used in two separate tests (see
equations 1 and 2 above).  Originally, independent SCALE and TOL values were
used for each test.  However, this provided only marginal improvement in
performance at the expense of greatly increasing the difficulty of determining
optimal values.

The current algorithm is not optimal.  However, further marginal improvements in
accuracy was only achieved at the expense of greatly increased computational
cost.  It is a matter of how hard you are willing to squeeze a lemon in order
to get out that last drop.


PROGRAM HISTORY:

Written by:  Gary Yagi, June 4, 1999
Cognizant programmer:  Gary Yagi
Revisions: 
  2003-08-14   NTT  Enabled despike for 3D images
  2005-10-11   lwk  converted operations to float, added POSONLY keyword