8 Variances, bad pixels, and quality – The art of error propagation

 8.1 Processing a variance array
 8.2 Example 9 – A code which does not propagate the variance array
 8.3 Example 10 – A code which DOES propagate the variance array
 8.4 The Quality component
 8.5 Rules for propagation of variance arrays

NDFs contain a lot more than just the data values. They contain so called headers, along with values such as Quality and Variance associated with each pixel. For some NDFs there is also axes information.

The easiest way to take a look at what structures exist within an NDF is to use hdstrace. A sample trace might produce an output such as:

  F2  <NDF>
  
     DATA_ARRAY     <ARRAY>         {structure}
        DATA(1021,200)  <_REAL>        41.46541,69.59,51.00168,48.32741,
                                       ... 169.1622,151.1605,150.1822,166.3774
        ORIGIN(2)      <_INTEGER>      1,1
        BAD_PIXEL      <_LOGICAL>      FALSE
  
     TITLE          <_CHAR*8>       ’FEIGE 34’
     UNITS          <_CHAR*9>       ’ELECTRONS’
     VARIANCE       <ARRAY>         {structure}
        DATA(1021,200)  <_REAL>        521.5062,571.9399,534.9544,532.3478,
                                       ... 587.9011,556.7531,627.0024,590.9243
        ORIGIN(2)      <_INTEGER>      1,1
  
     MORE           <EXT>           {structure}
        FITS(174)      <_CHAR*80>      ’SIMPLE  =                    T’,’BI...’
                                       ... ’15TELE  =
  ...’,’PACKEND’,’END’
        CCDPACK        <CCDPACK_EXT>   {structure}
           DEBIAS         <_CHAR*24>      ’Tue Apr 29 21:11:30 1997’End of
  Trace.

We can see from this that there are arrays for both the data and the variances. There is also some FITS information, and the file has at some point been processed by CCDPACK. Let’s write a new code that processes the variances in a similar way to how clip processed the data array.

8.1 Processing a variance array

Let’s illustrate the way in which the NDF library can handle variance information by developing a simple application. In the following example, an image is searched for “spikes” which are potential cosmic ray events.

When spikes are detected, the pixels are set to what is known as a BAD value. BAD pixels arise in situations where the values of the pixels are either unknown, useless, or meaningless. In this simple routine, we assign BAD values to pixels whose information is corrupted by cosmic rays.

In this example, any variance array present in the NDF is left untouched. This means that the variance array is no longer valid. It is, therefore, not passed on to the output image. Guidelines for when components of an NDF should and should not be propagated are listed in SUN/33. We will deal with a case where the variance information is propagated later in this section.

8.2 Example 9 – A code which does not propagate the variance array

Code:

        SUBROUTINE ZAPPER(STATUS)
  *
  * This routine zaps all pixels whose value goes above a certain
  * threshold defined by the user. A new NDF is propagated from the
  * old one. In this example the variance array (if present) is
  * not modified.
  *
  * Parameters
  *
  * IN = NDF (Read)
  *  Input NDF
  *
  * OUT = NDF (Write)
  *  Output NDF
  *
  * THRESH = _REAL (Read)
  *  Threshold for zapping
  *
        IMPLICIT NONE     ! No implicit typing
        INCLUDE ’SAE_PAR’ ! SAE constants
        INTEGER STATUS    ! Global Status
  *
  * Local Variables
  *
        INTEGER DIM(2)    ! Dimensions of image
        INTEGER NPIX      ! Number of pixels
        INTEGER INDF      ! Input NDF identifier
        INTEGER ONDF      ! Output NDF identifier
        INTEGER IPNTR     ! Input NDF pointer
        INTEGER OPNTR     ! Output NDF pointer
        INTEGER NDIM      ! Number of dimensions
        REAL THRESH       ! Threshold value
  *
  * Check the "inherited status" before starting the code proper
  *
        IF (STATUS .NE. SAI__OK) RETURN
  *
  * Obtain the input NDF. Find out how big it is in both dimensions
  *
        CALL NDF_ASSOC(’IN’,’READ’,INDF,STATUS)
        CALL NDF_DIM(INDF,2,DIM,NDIM,STATUS)
  *
  * Read in the threshold value
  *
        CALL PAR_GETOR(’THRESH’,THRESH,STATUS)
  *
  * Create a new output NDF. Model it on the old one. Propagate and
  * DATA information. The variances will NOT be propagated.
  *
        CALL NDF_PROP(INDF,’Data,NoVariance’,’OUT’,ONDF,STATUS)
  *
  * Map the data arrays in both the input and the output files
  *
        CALL NDF_MAP(INDF,’Data’,’_REAL’,’READ’,IPNTR,NPIX,STATUS)
        CALL NDF_MAP(ONDF,’Data’,’_REAL’,’WRITE’,OPNTR,NPIX,STATUS)
  *
  * Call the main working subroutine. Write new values to output
  *
        CALL ZAP(THRESH,NPIX,%VAL(IPNTR),%VAL(OPNTR),STATUS)
  *
  * Close down
  *
        CALL NDF_END(STATUS)
        END
  
        SUBROUTINE ZAP(THRESH,NPIX,IMAGE,OUT,STATUS)
  *
  * Zap pixels with more counts than THRESH by assigning BAD values to them
  *
        IMPLICIT NONE       ! No implicit typing
        INCLUDE ’SAE_PAR’   ! Standard SAE constants
        INCLUDE ’PRM_PAR’   ! Define BAD constants
  *
  * (>) - Given,  (<) - Output
  *
        REAL THRESH         ! (>) Threshold value
        INTEGER NPIX        ! (>) Number of pixels in images
        REAL IMAGE(NPIX)    ! (>) Array of input pixel values
        REAL OUT(NPIX)      ! (<) Array of output pixel values
        INTEGER STATUS      ! (<>) Global Status
  *
  * Local variables
  *
        INTEGER N
  *
  * Go through image and find excess values
  *
        DO 1, N = 1, NPIX
          IF (IMAGE(N).NE.VAL__BADR .AND. IMAGE(N).GT.THRESH) THEN
            OUT(N) = VAL__BADR
          ELSE
            OUT(N) = IMAGE(N)
          ENDIF
   1    CONTINUE
        END

Interface file:

  interface ZAPPER
  
   parameter IN
    position 1
    prompt ’Input file’
   endparameter
  
   parameter OUT
    position 2
    prompt ’Output file’
   endparameter
  
   parameter THRESH
    position 3
    type _REAL
    prompt ’Threshold value’
   end parameter
  
  endinterface

Note how the interface file now includes a parameter which is a REAL number as well as filenames.

The hdstrace application reveals that if the input NDF contained a variance array, it does not exist in the output NDF. This is because the NDF_PROP call did not explicitly list it. In order to get the variance array propagated through to the output, let’s adapt the code and call it zapper2.

8.3 Example 10 – A code which DOES propagate the variance array

This code propagates the variance array. Note you will need to copy the zapper.ifl file to a zapper2.ifl file.

        SUBROUTINE ZAPPER2(STATUS)
  *
  * This routine zaps all pixels whose value goes above a certain
  * threshold defined by the user. A new NDF is propagated from the
  * old one.
  *
  * Parameters
  *
  * IN = NDF (Read)
  *  Input NDF
  *
  * OUT = NDF (Write)
  *  Output NDF
  *
  * THRESH = _REAL (Read)
  *  Threshold for zapping
  *
        IMPLICIT NONE     ! No implicit typing
        INCLUDE ’SAE_PAR’ ! SAE constants
        INTEGER STATUS    ! Global Status
  *
  * Local Variables
  *
        INTEGER DIM(2)    ! Dimensions of image
        INTEGER NPIX      ! Number of pixels
        INTEGER INDF      ! Input NDF identifier
        INTEGER ONDF      ! Output NDF identifier
        INTEGER IPNTR     ! Input NDF data pointer
        INTEGER OPNTR     ! Output NDF data pointer
        INTEGER VIPNTR    ! Input NDF variance pointer
        INTEGER VOPNTR    ! Output NDF variance pointer
        INTEGER NDIM      ! Number of dimensions
        REAL THRESH       ! Threshold value
  *
  * Check the "inherited status" before starting the code proper
  *
        IF (STATUS .NE. SAI__OK) RETURN
  *
  * Begin the NDF context
  *
        CALL NDF_BEGIN
  *
  * Obtain the input NDF. Find out how big it is in both dimensions
  *
        CALL NDF_ASSOC(’IN’,’READ’,INDF,STATUS)
        CALL NDF_DIM(INDF,2,DIM,NDIM,STATUS)
  *
  * Read in the threshold value
  *
        CALL PAR_GET0R(’THRESH’,THRESH,STATUS)
  *
  * Create a new output NDF. Model it on the old one. Propagate
  * DATA and VARIANCES.
  *
        CALL NDF_PROP(INDF,’Data,Variance’,’OUT’,ONDF,STATUS)
  *
  * Map the data arrays in both the input and the output files
  *
        CALL NDF_MAP(INDF,’Data’,’_REAL’,’READ’,IPNTR,NPIX,STATUS)
        CALL NDF_MAP(INDF,’Variance’,’_REAL’,’READ’,VIPNTR,NPIX,STATUS)
        CALL NDF_MAP(ONDF,’Data’,’_REAL’,’WRITE’,OPNTR,NPIX,STATUS)
        CALL NDF_MAP(ONDF,’Variance’,’_REAL’,’WRITE’,VOPNTR,NPIX,STATUS)
  *
  * Call the main working subroutine. Write new values to output
  *
        CALL ZAP(THRESH,NPIX,%VAL(IPNTR),%VAL(OPNTR),STATUS)
  *
  * Close down
  *
        CALL NDF_END(STATUS)
        END
  
        SUBROUTINE ZAP(THRESH,NPIX,IMAGE,OUT,STATUS)
  *
  * Zap pixels with more counts than THRESH by assigning BAD values to them
  *
        IMPLICIT NONE       ! No implicit typing
        INCLUDE ’SAE_PAR’   ! Standard SAE constants
        INCLUDE ’PRM_PAR’   ! Define BAD constants
  *
  * (>) - Given,  (<) - Output
  *
        REAL THRESH         ! (>) Threshold value
        INTEGER NPIX        ! (>) Number of pixels in images
        REAL IMAGE(NPIX)    ! (>) Array of input pixel values
        REAL OUT(NPIX)      ! (<) Array of output pixel values
        INTEGER STATUS      ! (<>) Global Status
  *
  * Local variables
  *
        INTEGER N
  *
  * Go through image and find excess values
  *
        DO 1, N = 1, NPIX
          IF (IMAGE(N).NE.VAL__BADR .AND. IMAGE(N).GT.THRESH) THEN
            OUT(N) = VAL__BADR
          ELSE
            OUT(N) = IMAGE(N)
          ENDIF
   1    CONTINUE
        END

Note that although the variance array was mapped, none of its values were changed in this simple application. The variance array was mapped and therefore allocated pointers (although this wasn’t necessary to force the variance propagation). The interested reader might wish to try to adapt the call to ZAP to fix the variance values using these pointers.

8.4 The Quality component

In addition to the Data and Variance components, some NDFs also contain a Quality component. The quality array is not usually accessed directly by the user. Instead, the processing of the Data and Variance arrays automatically takes into account the values in the Quality array unless the application explicitly accesses and maps the Quality array.

Why bother with a Quality array when there is a variance array? The Quality value of a pixel can be used to describe its suitability for performing a particular task. For example, imagine a situation where photometry is performed on an object. The quality array could be set to one value for the object and another for the sky pixels.

Manipulation of the Quality array is not commonly done in most applications and so there will be no further discussion of it here. SUN/33 contains a full discussion of how to access the Quality array.

8.5 Rules for propagation of variance arrays

It is not appropriate to propagate the values in the variance array in all circumstances. For example, adding two NDFs together will result in new variances for the result. If these new variances are not calculated, the output NDF must not be allowed to contain a variance array. On the other hand, if a constant is added to an NDF, the variances can be passed on to the output with no calculations necessary.

The default used by the routine NDF_PROP is not to propagate the data, quality, or variance arrays. This explains why is was necessary to explicitly list what we wanted to propagate in the last example.