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