The following example is based around a simple algorithm which detects prominent pixels (e.g. data spikes) in a 2-dimensional image and replaces them with bad pixels. It is typical of applications which take a single NDF as input and produce a new NDF with the same size as output. It illustrates the use of propagation (NDF_PROP) in producing the new output NDF using the input as a template. Note that this application modifies the data array but does not handle the variance array, which will therefore become invalid and is not propagated.
This example also illustrates how bad pixels might be handled in a reasonably realistic image-processing algorithm. No attempt is made here to distinguish cases where bad pixels are present from those where they are not, and we do not really know afterwards if there are any bad pixels in the output image (although a check for this could easily be added). The output bad-pixel flag is therefore left with its default value of .TRUE..
SUBROUTINE ZAPPIX( STATUS )
*+
* Name:
* ZAPPIX
* Purpose:
* Zap prominent pixels.
* Description:
* This routine "zaps" prominent pixels in a 2-dimensional image
* (stored in the data array of an NDF). It searches for pixels
* which deviate by more than a specified amount from the mean of
* their nearest neighbours, and replaces them with bad pixels.
* ADAM Parameters:
* IN = NDF (Read)
* Input NDF data structure.
* OUT = NDF (Write)
* The output NDF data structure.
* THRESH = _REAL (Read)
* Threshold for zapping pixels; pixels will be set bad if they
* deviate by more than this amount from the mean of their
* nearest neighbours (the absolute value of THRESH is used).
* Implementation Status:
* This routine correctly handles bad pixels but does not handle NDF
* variance arrays. Real arithmetic is used.
*-
* Type Definitions:
IMPLICIT NONE ! No implicit typing
* Global Constants:
INCLUDE 'SAE_PAR' ! Standard SAE constants
* Status:
INTEGER STATUS ! Global status
* Local Variables:
INTEGER DIM( 2 ) ! Image dimension sizes
INTEGER EL ! Number of mapped values
INTEGER INDF1 ! Input NDF identifier
INTEGER INDF2 ! Output NDF identifier
INTEGER NDIM ! Number of NDF dimensions (junk)
INTEGER PNTR1( 1 ) ! Pointer to mapped input values
INTEGER PNTR2( 1 ) ! Pointer to mapped output values
REAL THRESH ! Threshold for zapping pixels
*.
* Check inherited global status.
IF ( STATUS .NE. SAI__OK ) RETURN
* Begin an NDF context.
CALL NDF_BEGIN
* Obtain the input NDF and obtain its first two dimension sizes.
CALL NDF_ASSOC( 'IN', 'READ', INDF1, STATUS )
CALL NDF_DIM( INDF1, 2, DIM, NDIM, STATUS )
* Obtain a threshold value.
CALL PAR_GET0R( 'THRESH', THRESH, STATUS )
* Create an output NDF based on the input one. Propagate the AXIS,
* QUALITY and UNITS components.
CALL NDF_PROP( INDF1, 'Axis,Quality,Units', 'OUT', INDF2, STATUS )
* Map the input and output data arrays for reading and writing
* respectively.
CALL NDF_MAP( INDF1, 'Data', '_REAL', 'READ', PNTR1, EL, STATUS )
CALL NDF_MAP( INDF2, 'Data', '_REAL', 'WRITE', PNTR2, EL, STATUS )
* Process the input array, writing the new values to the output array.
CALL ZAPIT( ABS( THRESH ), DIM( 1 ), DIM( 2 ), %VAL( PNTR1( 1 ) ),
: %VAL( PNTR2( 1 ) ), STATUS )
* End the NDF context (this cleans everything up).
CALL NDF_END( STATUS )
* If an error occurred, then report a contextual message.
IF ( STATUS .NE. SAI__OK ) THEN
CALL ERR_REP( 'ZAPPIX_ERR',
: 'ZAPPIX: Error zapping prominent pixels in an image.',
: STATUS )
END IF
END
SUBROUTINE ZAPIT( THRESH, NX, NY, A, B, STATUS )
*+
* Name:
* ZAPIT
* Purpose:
* Zap prominent pixels in an image.
* Invocation:
* CALL ZAPIT( THRESH, NX, NY, A, B, STATUS )
* Description:
* The routine finds all pixels in a 2-dimensional image which
* deviate by more than a specified amount from the mean of their
* nearest neighbours and replaces them with the bad pixel value.
* Bad pixels in the input image are correctly handled.
* Arguments:
* THRESH = REAL (Given)
* The threshold for zapping pixels.
* NX = INTEGER (Given)
* X dimension of image.
* NY = INTEGER (Given)
* Y dimension of image.
* A( NX, NY ) = REAL (Given)
* The input image.
* B( NX, NY ) = REAL (Returned)
* The output image.
* STATUS = INTEGER (Given and Returned)
* The global status.
*-
* Type Definitions:
IMPLICIT NONE ! No implicit typing
* Global Constants:
INCLUDE 'SAE_PAR' ! Standard SAE constants
INCLUDE 'PRM_PAR' ! Define VAL__BADR constant
* Arguments Given:
REAL THRESH
INTEGER NX
INTEGER NY
REAL A( NX, NY )
* Arguments Returned:
REAL B( NX, NY )
* Status:
INTEGER STATUS ! Global status
* Local Variables:
INTEGER IIX ! Loop counter for neighbours
INTEGER IIY ! Loop counter for neighbours
INTEGER IX ! Loop counter for image pixels
INTEGER IY ! Loop counter for image pixels
INTEGER N ! Number of good neighbours
REAL DIFF ! Deviation from mean of neighbours
REAL S ! Sum of good neighbours
*.
* Check inherited global status.
IF ( STATUS .NE. SAI__OK ) RETURN
* Loop through all the pixels in the image.
DO 4 IY = 1, NY
DO 3 IX = 1, NX
* If the input pixel is bad, then so is the output pixel.
IF ( A( IX, IY ) .EQ. VAL__BADR ) THEN
B( IX, IY ) = VAL__BADR
* Otherwise, loop to find the average of the nearest neighbours.
ELSE
S = 0.0
N = 0
DO 2 IIY = MAX( 1, IY - 1 ), MIN( NY, IY + 1 )
DO 1 IIX = MAX( 1, IX - 1 ), MIN( NX, IX + 1 )
* Only count neighbours which are not bad themselves.
IF ( A( IIX, IIY ) .NE. VAL__BADR ) THEN
S = S + A( IIX, IIY )
N = N + 1
END IF
1 CONTINUE
2 CONTINUE
* If all the neighbours were bad, then just copy the central pixel.
IF ( N .EQ. 0 ) THEN
B( IX, IY ) = A( IX, IY )
* Otherwise, see if the central pixel deviates by more than THRESH from
* the average. If not, copy it. If so, set it bad.
ELSE
DIFF = A( IX, IY ) - ( S / REAL( N ) )
IF ( ABS( DIFF ) .LE. THRESH ) THEN
B( IX, IY ) = A( IX, IY )
ELSE
B( IX, IY ) = VAL__BADR
END IF
END IF
END IF
3 CONTINUE
4 CONTINUE
END
The following is an example ADAM interface file (zappix.ifl) for the application above.
interface ZAPPIX
parameter IN # Input NDF
position 1
prompt 'Input NDF'
endparameter
parameter OUT # Output NDF
position 3
prompt 'Output NDF'
endparameter
parameter THRESH # Zapping threshold
position 2
type _REAL
prompt 'Threshold'
endparameter
endinterface