next up previous
Next: ADD Add Two NDF Data Structures
Up: EXAMPLE APPLICATIONS
Previous: READIMG Read an image into an NDF

ZAPPIX -- ``Zap'' Prominent Pixels in an Image  

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



next up previous
Next: ADD Add Two NDF Data Structures
Up: EXAMPLE APPLICATIONS
Previous: READIMG Read an image into an NDF


Starlink User Note 33
R.F. Warren-Smith
11th January 2000
E-mail:rfws@star.rl.ac.uk

Copyright © 2000 Council for the Central Laboratory of the Research Councils