IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5516


Ignore:
Timestamp:
Nov 15, 2005, 10:09:03 AM (20 years ago)
Author:
gusciora
Message:

SubtractBias was recoded. Significant mods to removeBadPixels.
Additional mods to other files.

Location:
trunk/psModules
Files:
15 edited

Legend:

Unmodified
Added
Removed
  • trunk/psModules/src/config/pmConfig.c

    r5435 r5516  
    33 *  @author PAP, IfA
    44 *
    5  *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    6  *  @date $Date: 2005-10-20 23:06:24 $
     5 *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
     6 *  @date $Date: 2005-11-15 20:09:03 $
    77 *
    88 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    1111#include <stdio.h>
    1212#include <string.h>
    13 #include <assert.h>
    1413#include "pslib.h"
    15 //#include "psAdditionals.h"
    1614#include "pmConfig.h"
    1715
    1816#define PS_SITE "PS_SITE"               // Name of the environment variable containing the site config file
    19 #define PS_DEFAULT_SITE "ipprc.config"     // Default site config file
    20 
    21 
     17#define PS_DEFAULT_SITE "ipprc.config"  // Default site config file
    2218
    2319/** readConfig
     
    3228    const char *description)            // Description of file
    3329{
    34     int numBadLines = 0;  // Number of bad lines in config file
     30    int numBadLines = 0;
     31
    3532    psLogMsg(__func__, PS_LOG_INFO, "Loading %s configuration from file %s\n",
    3633             description, name);
     
    5653must take precedence override the values set here.  This must be, somehow,
    5754coded.
     55 
     56XXX: Must load camera and recipe configuration if specified in the command line.
    5857 *****************************************************************************/
    5958bool pmConfigRead(
     
    119118    // file and store in psMetadata struct site.
    120119    //
     120
    121121    if (!readConfig(site, siteName, "site")) {
    122122        return false;
     
    249249    }
    250250
     251    //
     252    // Allow command line options to override defaults for logging.
     253    // XXX: Is it appropriate to use the ArgVerbosity function for this?
     254    //   A: it removes the options from the command line.
     255    //   B: will the pmConfigRead function always be called on initialization.
     256    //
     257    psS32 saveLogLevel = psLogGetLevel();
     258    psArgumentVerbosity(argc, argv);
     259    // XXX: substitute the string for the default log level for "2".
     260    if (2 == psLogGetLevel()) {
     261        psLogSetLevel(saveLogLevel);
     262    }
    251263    return(true);
    252264}
     
    266278    // Apply the rules
    267279    psMetadataIterator *ruleIter = psMetadataIteratorAlloc(rule, PS_LIST_HEAD, NULL); // Rule iterator
    268     psMetadataItem *ruleItem = NULL; // Item from the metadata
    269     bool match = true;   // Does it match?
     280    psMetadataItem *ruleItem = NULL;    // Item from the metadata
     281    bool match = true;                  // Does it match?
    270282    while ((ruleItem = psMetadataGetAndIncrement(ruleIter)) && match) {
    271283        // Check for the existence of the rule
     
    389401    const char *recipeName)
    390402{
    391     assert(camera);
    392     assert(recipeName);
    393 
    394     psMetadata *recipe = NULL; // Recipe to read
    395     bool mdok = true;   // Status of MD lookup
     403    PS_ASSERT_PTR_NON_NULL(camera, false);
     404    PS_ASSERT_PTR_NON_NULL(recipeName, false);
     405
     406    psMetadata *recipe = NULL;          // Recipe to read
     407    bool mdok = true;                   // Status of MD lookup
    396408    psMetadata *recipes = psMetadataLookupMD(&mdok, camera, "RECIPES"); // The list of recipes
    397409    if (! mdok || ! recipes) {
  • trunk/psModules/src/detrend/pmFlatField.c

    r5435 r5516  
    1818 *  @author Ross Harman, MHPCC
    1919 *
    20  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    21  *  @date $Date: 2005-10-20 23:06:24 $
     20 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     21 *  @date $Date: 2005-11-15 20:09:03 $
    2222 *
    2323 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    3636#include "pmMaskBadPixels.h"
    3737#include "pmFlatFieldErrors.h"
     38#include "pmSubtractBias.h"
    3839
    3940
     
    4243    const pmReadout *flat)
    4344{
     45    // XXX: Use the proper image and readout asserts.
     46    PS_ASSERT_PTR_NON_NULL(in, true);
     47    PS_ASSERT_PTR_NON_NULL(in->image, false);
     48    PS_ASSERT_PTR_NON_NULL(in->mask, false);
     49    PS_ASSERT_PTR_NON_NULL(flat, false);
     50    PS_ASSERT_PTR_NON_NULL(flat->image, false);
     51    if (in == NULL)
     52        printf("XXX: NULL\n");
     53
    4454    // XXX: Not sure if this is correct.  Must consult with IfA.
    4555    PS_ASSERT_PTR_NON_NULL(in->mask, false);
     56
     57    PS_WARN_PTR_NON_NULL(in->parent);
     58    if (in->parent != NULL) {
     59        PS_WARN_PTR_NON_NULL(in->parent->concepts);
     60    }
    4661    int i = 0;
    4762    int j = 0;
     
    5166    psElemType flatType;
    5267    psElemType maskType;
    53     psImage *inImage = NULL;
    5468    psImage *inMask = NULL;
    5569    psImage *flatImage = NULL;
    5670
    57 
    58     // Check for nulls
    59     if (in == NULL) {
    60         return true;       // Readout might have data in it
    61     } else if(flat==NULL) {
    62         psError( PS_ERR_BAD_PARAMETER_NULL, true,
    63                  PS_ERRORTEXT_pmFlatField_NULL_FLAT_READOUT);
    64         return false;
    65     }
    66 
    67     inImage = in->image;
     71    //
     72    // Determine trimmed image from metadata.
     73    //
     74    psImage *trimmedImg = p_psDetermineTrimmedImage(in);
    6875    flatImage = flat->image;
    69     if (inImage == NULL) {
    70         psError( PS_ERR_BAD_PARAMETER_NULL, true,
    71                  PS_ERRORTEXT_pmFlatField_NULL_INPUT_IMAGE);
    72         return false;
    73     } else if(flatImage == NULL) {
    74         psError( PS_ERR_BAD_PARAMETER_NULL, true,
    75                  PS_ERRORTEXT_pmFlatField_NULL_FLAT_IMAGE);
    76         return false;
    77     }
    7876    inMask = in->mask;
    7977
    8078    // Check input image and its mask are not larger than flat image
    8179
    82     if (inImage->numRows>flatImage->numRows || inImage->numCols>flatImage->numCols) {
     80    if (trimmedImg == NULL)
     81        printf("XXX: 00\n");
     82    if (flatImage == NULL)
     83        printf("XXX 01\n");
     84
     85    if (trimmedImg->numRows>flatImage->numRows || trimmedImg->numCols>flatImage->numCols) {
    8386        psError( PS_ERR_BAD_PARAMETER_SIZE, true,
    8487                 PS_ERRORTEXT_pmFlatField_SIZE_INPUT_IMAGE,
    85                  inImage->numRows, inImage->numCols, flatImage->numRows, flatImage->numCols);
     88                 trimmedImg->numRows, trimmedImg->numCols, flatImage->numRows, flatImage->numCols);
    8689        return false;
    8790    }
     
    9497
    9598    // Determine total offset based on image offset with chip offset
    96     totOffCol = inImage->col0 + in->col0;
    97     totOffRow = inImage->row0 + in->row0;
     99    totOffCol = trimmedImg->col0 + in->col0;
     100    totOffRow = trimmedImg->row0 + in->row0;
    98101
    99102    // Check that offsets are within image limits
     
    103106                 totOffRow, totOffCol, flatImage->numRows, flatImage->numCols);
    104107        return false;
    105     } else if(totOffRow>=inImage->numRows || totOffCol>=inImage->numCols) {
     108    } else if(totOffRow>=trimmedImg->numRows || totOffCol>=trimmedImg->numCols) {
    106109        psError( PS_ERR_BAD_PARAMETER_SIZE, true,
    107110                 PS_ERRORTEXT_pmFlatField_OFFSET_INPUT_IMAGE,
    108                  totOffRow, totOffCol, inImage->numRows, inImage->numCols);
     111                 totOffRow, totOffCol, trimmedImg->numRows, trimmedImg->numCols);
    109112        return false;
    110113    } else if(totOffRow>=inMask->numRows || totOffCol>=inMask->numCols) {
     
    116119
    117120    // Check for incorrect types
    118     inType = inImage->type.type;
     121    inType = trimmedImg->type.type;
    119122    flatType = flatImage->type.type;
    120123    maskType = inMask->type.type;
     
    145148case PS_TYPE_##TYPE:                                                                                         \
    146149    /* Per Eugene's request, use two sets of loops: first to fill mask, second to avoid div with bad pix */  \
    147     for(j = totOffRow; j < inImage->numRows; j++) {                                                          \
    148         for(i = totOffCol; i < inImage->numCols; i++) {                                                      \
     150    for(j = totOffRow; j < trimmedImg->numRows; j++) {                                                          \
     151        for(i = totOffCol; i < trimmedImg->numCols; i++) {                                                      \
    149152            if(flatImage->data.TYPE[j][i] <= 0.0) {                                                          \
    150153                /* Negative or zero flat pixels shall be masked in input image as  PM_MASK_FLAT */           \
     
    154157        }                                                                                                    \
    155158    }                                                                                                        \
    156     for(j = totOffRow; j < inImage->numRows; j++) {                                                          \
    157         for(i = totOffCol; i < inImage->numCols; i++) {                                                      \
     159    for(j = totOffRow; j < trimmedImg->numRows; j++) {                                                          \
     160        for(i = totOffCol; i < trimmedImg->numCols; i++) {                                                      \
    158161            if(!inMask->data.PS_TYPE_MASK_DATA[j][i]) {                                                      \
    159162                /* Module shall divide the input image by the flat-fielded image */                          \
    160                 inImage->data.TYPE[j][i] /= flatImage->data.TYPE[j][i];                                      \
     163                trimmedImg->data.TYPE[j][i] /= flatImage->data.TYPE[j][i];                                      \
    161164            }                                                                                                \
    162165        }                                                                                                    \
  • trunk/psModules/src/detrend/pmMaskBadPixels.c

    r5170 r5516  
    1919 *  @author Ross Harman, MHPCC
    2020 *
    21  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    22  *  @date $Date: 2005-09-28 20:43:52 $
     21 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
     22 *  @date $Date: 2005-11-15 20:09:03 $
    2323 *
    2424 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    3535#include "pmMaskBadPixels.h"
    3636#include "pmMaskBadPixelsErrors.h"
     37#include "pmSubtractBias.h"
     38
     39/******************************************************************************
     40GrowPixel(inMask, row, col, radius, maskVal): This private routine takes an
     41input image mask and a pixel location, then sets (logical or) all pixels with
     42parameter radius if that pixel to maskVal parameter.
     43 *****************************************************************************/
     44psBool GrowPixel(
     45    psImage *inMask,
     46    psS32 col,
     47    psS32 row,
     48    psS32 radius,
     49    psU32 maskVal)
     50{
     51    psS32 rowMin = PS_MAX(row-radius, 0);
     52    psS32 rowMax = PS_MIN(row+radius+1, inMask->numRows);
     53    psS32 colMin = PS_MAX(col-radius, 0);
     54    psS32 colMax = PS_MIN(col+radius+1, inMask->numCols);
     55    psF32 squareRadius = PS_SQR(radius);
     56
     57
     58    for(psS32 i=rowMin; i<rowMax; i++) {
     59        for(psS32 j=colMin; j<colMax; j++) {
     60            // Old code:
     61            //            psF32 rRound = 0.5 + sqrtf((psF32) (((col-j)*(col-j)) + ((row-i)*(row-i))));
     62            //            if(rRound <= radius) {
     63            psF32 squareDis = (psF32) (((col-j)*(col-j)) + ((row-i)*(row-i)));
     64            if (squareDis <= squareRadius) {
     65                inMask->data.U8[i][j] |= maskVal;
     66            }
     67        }
     68    }
     69    return(true);
     70}
     71
     72/******************************************************************************
     73GetRadius(inImg, col, row, sat, growVal, grow): This private routine takes an
     74input image and pixel location and determines what radius that pixel must grow
     75by.
     76 
     77//XXX: Inline this or macro it for speed.
     78 *****************************************************************************/
     79static psS32 GetRadius(
     80    psImage *inImg,
     81    psS32 col,
     82    psS32 row,
     83    psF32 sat,
     84    psU32 growVal,
     85    psS32 grow)
     86{
     87    psS32 growRadius = 0;
     88    if (inImg->type.type == PS_TYPE_F32) {
     89        if(inImg->data.F32[row][col] == growVal) {
     90            growRadius = grow;
     91        }
     92        if (inImg->data.F32[row][col] > sat) {
     93            growRadius = PS_MAX(growRadius+1, 2);
     94        }
     95    } else if (inImg->type.type == PS_TYPE_S32) {
     96        if(inImg->data.S32[row][col] == growVal) {
     97            growRadius = grow;
     98        }
     99        if (inImg->data.S32[row][col] > sat) {
     100            growRadius = PS_MAX(growRadius+1, 2);
     101        }
     102    } else if (inImg->type.type == PS_TYPE_U16) {
     103        if(inImg->data.U16[row][col] == growVal) {
     104            growRadius = grow;
     105        }
     106        if (inImg->data.U16[row][col] > sat) {
     107            growRadius = PS_MAX(growRadius+1, 2);
     108        }
     109    } else {
     110        psError( PS_ERR_BAD_PARAMETER_TYPE, true,
     111                 PS_ERRORTEXT_pmMaskBadPixels_TYPE_MASK_IMAGE,
     112                 inImg->type.type);
     113    }
     114    return(growRadius);
     115}
     116
    37117
    38118bool pmMaskBadPixels(
    39119    pmReadout *in,
    40     const psImage *mask,
     120    const pmReadout *mask,
    41121    unsigned int maskVal,
    42122    float sat,
     
    44124    int grow)
    45125{
     126    // XXX: Review this code then put all asserts at the top.
     127    PS_ASSERT_PTR_NON_NULL(in, true);
     128    PS_ASSERT_PTR_NON_NULL(in->image, false);
     129    PS_WARN_PTR_NON_NULL(in->parent);
     130    if (in->parent != NULL) {
     131        PS_WARN_PTR_NON_NULL(in->parent->concepts);
     132    }
     133    PS_ASSERT_PTR_NON_NULL(mask, false);
    46134    int i = 0;
    47135    int j = 0;
    48     int jj = 0;
    49     int ii = 0;
    50     int rRound = 0;
    51     int rowMin = 0;
    52     int rowMax = 0;
    53     int colMin = 0;
    54     int colMax = 0;
    55136    int totOffCol = 0;
    56137    int totOffRow = 0;
    57     float r = 0.0f;
    58138    psElemType inType;
    59139    psElemType maskType;
    60     psImage *inImage = NULL;
    61140    psImage *inMask = NULL;
    62141
    63 
    64     // Check for nulls
    65     if (in == NULL) {
    66         return true;   // Readout may not have data in it
    67     } else if(mask==NULL) {
    68         psError( PS_ERR_BAD_PARAMETER_NULL, true,
    69                  PS_ERRORTEXT_pmMaskBadPixels_NULL_MASK_IMAGE);
    70         return false;
    71     }
    72 
    73     inImage = in->image;
    74     if (inImage == NULL) {
    75         psError( PS_ERR_BAD_PARAMETER_NULL, true,
    76                  PS_ERRORTEXT_pmMaskBadPixels_NULL_INPUT_IMAGE);
    77         return false;
    78     } else if(in->mask == NULL) {
    79         in->mask = psImageAlloc(inImage->numCols, inImage->numRows, PS_TYPE_MASK);
    80         memset(in->mask->data.V[0], 0, inImage->numCols*inImage->numRows*PSELEMTYPE_SIZEOF(PS_TYPE_MASK));
     142    //
     143    // Determine trimmed image from metadata.
     144    //
     145    psImage *trimmedImg = p_psDetermineTrimmedImage(in);
     146    if(in->mask == NULL) {
     147        in->mask = psImageAlloc(trimmedImg->numCols, trimmedImg->numRows, PS_TYPE_MASK);
     148        PS_IMAGE_SET_U8(in->mask, 0);
    81149    }
    82150    inMask = in->mask;
    83151
    84152    // Check input image and its mask are not larger than mask
    85     if(inImage->numRows > mask->numRows || inImage->numCols > mask->numCols) {
     153    if(trimmedImg->numRows > mask->image->numRows || trimmedImg->numCols > mask->image->numCols) {
    86154        psError( PS_ERR_BAD_PARAMETER_SIZE, true,
    87155                 PS_ERRORTEXT_pmMaskBadPixels_SIZE_INPUT_IMAGE,
    88                  inImage->numRows, inImage->numCols, mask->numRows, mask->numCols);
    89         return false;
    90     } else if(inMask->numRows>mask->numRows || inMask->numCols>mask->numCols) {
     156                 trimmedImg->numRows, trimmedImg->numCols, mask->image->numRows, mask->image->numCols);
     157        return false;
     158    } else if(inMask->numRows > mask->image->numRows || inMask->numCols > mask->image->numCols) {
    91159        psError( PS_ERR_BAD_PARAMETER_SIZE, true,
    92160                 PS_ERRORTEXT_pmMaskBadPixels_SIZE_MASK_IMAGE,
    93                  inMask->numRows, inMask->numCols, mask->numRows, mask->numCols);
     161                 inMask->numRows, inMask->numCols, mask->image->numRows, mask->image->numCols);
    94162        return false;
    95163    }
    96164
    97165    // Determine total offset based on image offset with chip offset
    98     totOffCol = inImage->col0 + in->col0;
    99     totOffRow = inImage->row0 + in->row0;
     166    totOffCol = trimmedImg->col0 + in->col0;
     167    totOffRow = trimmedImg->row0 + in->row0;
    100168
    101169    // Check that offsets are within image limits
    102     if(totOffRow>=mask->numRows || totOffCol>=mask->numCols) {
     170    if(totOffRow>=mask->image->numRows || totOffCol>=mask->image->numCols) {
    103171        psError( PS_ERR_BAD_PARAMETER_SIZE, true,
    104172                 PS_ERRORTEXT_pmMaskBadPixels_OFFSET_MASK_IMAGE,
    105                  totOffRow, totOffCol, mask->numRows, mask->numCols);
    106         return false;
    107     } else if(totOffRow>=inImage->numRows || totOffCol>=inImage->numCols) {
     173                 totOffRow, totOffCol, mask->image->numRows, mask->image->numCols);
     174        return false;
     175    } else if(totOffRow>=trimmedImg->numRows || totOffCol>=trimmedImg->numCols) {
    108176        psError( PS_ERR_BAD_PARAMETER_SIZE, true,
    109177                 PS_ERRORTEXT_pmMaskBadPixels_OFFSET_INPUT_IMAGE,
    110                  totOffRow, totOffCol, inImage->numRows, inImage->numCols);
     178                 totOffRow, totOffCol, trimmedImg->numRows, trimmedImg->numCols);
    111179        return false;
    112180    } else if(totOffRow>=inMask->numRows || totOffCol>=inMask->numCols) {
     
    118186
    119187    // Check for incorrect types
    120     inType = inImage->type.type;
    121     maskType = mask->type.type;
     188    inType = trimmedImg->type.type;
     189    maskType = mask->image->type.type;
    122190    if(PS_IS_PSELEMTYPE_COMPLEX(inType)) {
    123191        psError( PS_ERR_BAD_PARAMETER_TYPE, true,
     
    132200    }
    133201
    134     // Macro for all PS types
    135     #define PM_BAD_PIXELS(TYPE)                                                                              \
    136 case PS_TYPE_##TYPE:                                                                                         \
    137     for(j=totOffRow; j<inImage->numRows; j++) {                                                              \
    138         for(i=totOffCol; i<inImage->numCols; i++) {                                                          \
    139             \
    140             /* Pixels with flux greater than sat shall be masked */                                          \
    141             if(inImage->data.TYPE[j][i] > sat) {                                                             \
    142                 inMask->data.PS_TYPE_MASK_DATA[j][i] |= PM_MASK_SAT;                                         \
    143             }                                                                                                \
    144             \
    145             /* Pixels which satisfy maskVal shall be masked */                                               \
    146             inMask->data.PS_TYPE_MASK_DATA[j][i] |= (mask->data.PS_TYPE_MASK_DATA[j][i]&maskVal);            \
    147             \
    148             /* Pixels which satisfy growVal and within the grow radius shall be masked */                    \
    149             if(mask->data.PS_TYPE_MASK_DATA[j][i] & growVal) {                                               \
    150                 rowMin = MAX(j-grow, 0);                                                                     \
    151                 rowMax = MIN(j+grow+1, inImage->numRows);                                                    \
    152                 colMin = MAX(i-grow, 0);                                                                     \
    153                 colMax = MIN(i+grow+1, inImage->numCols);                                                    \
    154                 for(jj=rowMin; jj<rowMax; jj++) {                                                            \
    155                     for(ii=colMin; ii<colMax; ii++) {                                                        \
    156                         r = sqrtf((ii-i)*(ii-i)+(jj-j)*(jj-j));                                              \
    157                         rRound = r + 0.5;                                                                    \
    158                         if(rRound <= grow) {                                                                 \
    159                             inMask->data.PS_TYPE_MASK_DATA[jj][ii] |=                                        \
    160                                     (mask->data.PS_TYPE_MASK_DATA[j][i]&growVal);                            \
    161                         }                                                                                    \
    162                     }                                                                                        \
    163                 }                                                                                            \
    164             }                                                                                                \
    165         }                                                                                                    \
    166     }                                                                                                        \
    167     break;
    168 
    169     // Switch to call bad pixel masking macro defined above
    170     switch(inType) {
    171         PM_BAD_PIXELS(U8);
    172         PM_BAD_PIXELS(U16);
    173         PM_BAD_PIXELS(U32);
    174         PM_BAD_PIXELS(U64);
    175         PM_BAD_PIXELS(S8);
    176         PM_BAD_PIXELS(S16);
    177         PM_BAD_PIXELS(S32);
    178         PM_BAD_PIXELS(S64);
    179         PM_BAD_PIXELS(F32);
    180         PM_BAD_PIXELS(F64);
    181     default:
    182         psError( PS_ERR_BAD_PARAMETER_TYPE, true,
    183                  PS_ERRORTEXT_pmMaskBadPixels_TYPE_UNSUPPORTED,
    184                  inType);
    185     }
    186 
    187     return false;
     202    //
     203    // This is the main loop which examines each pixel and masks pixels if
     204    //    A: The mask matches the maskValue
     205    //    B: The pixel is larger than the saturation value
     206    //    C: The pixel equals the grow value (in which case a circle is masked)
     207    //
     208    for(i=totOffRow; i<trimmedImg->numRows; i++) {
     209        for(j=totOffCol; j<trimmedImg->numCols; j++) {
     210            //
     211            // A: Pixels which satisfy maskVal shall be masked
     212            //
     213            if (mask->image->data.U8[i][j] == maskVal) {
     214                in->mask->data.U8[i][j] |= maskVal;
     215            }
     216
     217            //
     218            // We first determine how much we need to grow by and store this in
     219            // growRadius.
     220            //
     221            psS32 growRadius = GetRadius(trimmedImg, j, i, sat, growVal, grow);
     222
     223            //
     224            // Grow the pixel mask if necessary.
     225            //
     226            if (growRadius != 0) {
     227                GrowPixel(in->mask, j, i, growRadius, maskVal);
     228            }
     229        }
     230    }
     231
     232    return true;
    188233}
  • trunk/psModules/src/detrend/pmMaskBadPixels.h

    r5170 r5516  
    1919 *  @author Ross Harman, MHPCC
    2020 *
    21  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    22  *  @date $Date: 2005-09-28 20:43:52 $
     21 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
     22 *  @date $Date: 2005-11-15 20:09:03 $
    2323 *
    2424 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    3636} pmMaskValue;
    3737
     38// XXX: Use PS_MIN, PS_MAX
    3839/** Macro to find maximum of two numbers */
    3940#define MAX(A,B)((A)>=(B)?(A):(B))
     
    5354bool pmMaskBadPixels(
    5455    pmReadout *in,          ///< Readout containing input image data.
    55     const psImage *mask,    ///< Mask data to be added to readout mask data.
     56    const pmReadout *mask,   ///< Mask data to be added to readout mask data.
    5657    unsigned int maskVal,   ///< Mask value to determine what to add to input mask.
    5758    float sat,              ///< Saturation limit to mask bad pixels.
  • trunk/psModules/src/detrend/pmNonLinear.c

    r5435 r5516  
    55 *  @author GLG, MHPCC
    66 *
    7  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-10-20 23:06:24 $
     7 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2005-11-15 20:09:03 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    2222
    2323#include "pmNonLinear.h"
     24#include "pmSubtractBias.h"
    2425
    2526/******************************************************************************
     
    3738    PS_ASSERT_IMAGE_TYPE(inputReadout->image, PS_TYPE_F32, NULL);
    3839    PS_ASSERT_PTR_NON_NULL(input1DPoly, NULL);
     40    PS_WARN_PTR_NON_NULL(inputReadout->parent);
     41    if (inputReadout->parent != NULL) {
     42        PS_WARN_PTR_NON_NULL(inputReadout->parent->concepts);
     43    }
     44    //
     45    // Determine trimmed image from metadata.
     46    //
     47    psImage *trimmedImg = p_psDetermineTrimmedImage(inputReadout);
    3948
    40     psS32 i;
    41     psS32 j;
    42 
    43     for (i=0;i<inputReadout->image->numRows;i++) {
    44         for (j=0;j<inputReadout->image->numCols;j++) {
    45             inputReadout->image->data.F32[i][j] = psPolynomial1DEval(input1DPoly, inputReadout->image->data.F32[i][j]);
     49    for (psS32 i=0;i<trimmedImg->numRows;i++) {
     50        for (psS32 j=0;j<trimmedImg->numCols;j++) {
     51            trimmedImg->data.F32[i][j] = psPolynomial1DEval(input1DPoly,
     52                                         trimmedImg->data.F32[i][j]);
    4653        }
    4754    }
     
    6875    PS_ASSERT_PTR_NON_NULL(inputReadout->image,NULL);
    6976    PS_ASSERT_IMAGE_TYPE(inputReadout->image, PS_TYPE_F32, NULL);
     77    PS_WARN_PTR_NON_NULL(inputReadout->parent);
     78    if (inputReadout->parent != NULL) {
     79        PS_WARN_PTR_NON_NULL(inputReadout->parent->concepts);
     80    }
     81    //
     82    // Determine trimmed image from metadata.
     83    //
     84    psImage *trimmedImg = p_psDetermineTrimmedImage(inputReadout);
     85
    7086    psLookupTable *tmpLT = psLookupTableAlloc(filename, "%f %f", 0);
    7187    psS32 numLines = psLookupTableRead(tmpLT);
     
    7591                 "WARNING: Lookup Table is too small.  Returning original pmReadout.\n");
    7692    } else {
    77         for (psS32 i=0;i<inputReadout->image->numRows;i++) {
    78             for (psS32 j=0;j<inputReadout->image->numCols;j++) {
    79                 psF64 tmpD = psLookupTableInterpolate(tmpLT, inputReadout->image->data.F32[i][j], 1);
     93        for (psS32 i=0;i<trimmedImg->numRows;i++) {
     94            for (psS32 j=0;j<trimmedImg->numCols;j++) {
     95                psF64 tmpD = psLookupTableInterpolate(tmpLT, trimmedImg->data.F32[i][j], 1);
    8096                if (!isnan(tmpD)) {
    81                     inputReadout->image->data.F32[i][j] = tmpD;
     97                    trimmedImg->data.F32[i][j] = tmpD;
    8298                } else {
    8399                    numPixels++;
  • trunk/psModules/src/imcombine/pmReadoutCombine.c

    r5170 r5516  
    55 *  @author GLG, MHPCC
    66 *
    7  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-09-28 20:43:52 $
     7 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2005-11-15 20:09:03 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    374374    return(output);
    375375}
     376
     377
     378/* This function measures the robust median at each of the minimum and maximum
     379 * coordinates and determines the difference and mean of the two values. The size
     380 * of the box used to make the measurement at each point is specified by the
     381 * configuration variable FRINGE_SQUARE_RADIUS. From the collection of
     382 * differences, the robust median is calculated, and returned as part of the
     383 * fringe statistics. For each fringe point, the values of delta and midValue are
     384 * also assigned and available to the user on return.
     385 */
     386
     387psStats *pmFringeStats(
     388    psArray *fringePoints,
     389    psImage *image,
     390    psMetadata *config)
     391{
     392    PS_ASSERT_PTR_NON_NULL(fringePoints, NULL);
     393    //    for (psS32 i = 0 ; i < fringePoints->n ; i++) {
     394    //        if (!psMemCheckFringePoint((pmFringePoint *) fringePoints->data[i])) {
     395    //            psError(PS_ERR_UNKNOWN, true, "fringePoints->data[%d] is not of type pmFringePoint.\n");
     396    //            return(NULL);
     397    //        }
     398    //    }
     399    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
     400    PS_ASSERT_IMAGE_NON_EMPTY(image, NULL);
     401    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);
     402    PS_ASSERT_PTR_NON_NULL(config, NULL);
     403
     404    psBool rc;
     405    psF32 frSquareRadius = psMetadataLookupF32(&rc, config, "FRINGE_SQUARE_RADIUS");
     406    if (!rc) {
     407        psError(PS_ERR_UNKNOWN, true, "Could not determing the fringe radius from the metadata.\n");
     408    }
     409
     410    psRegion minRegion;
     411    psRegion maxRegion;
     412    psStats *minStats = psStatsAlloc(PS_STAT_ROBUST_MEAN);
     413    psStats *maxStats = psStatsAlloc(PS_STAT_ROBUST_MEAN);
     414    psStats *diffStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
     415    psVector *diffs = psVectorAlloc(fringePoints->n, PS_TYPE_F32);
     416
     417    //
     418    // Loop through each fringe point.  Determine the robust mean around
     419    // the min and max for that fringe point.  Save the difference between
     420    // those two numbers in psVector diffs.
     421    //
     422    // XXX: Ensure you have the radius correct.  (add 1 to x1 and y1?)
     423    //
     424    for (psS32 i = 0 ; i < fringePoints->n ; i++) {
     425        pmFringePoint *fp = (pmFringePoint *) fringePoints->data[i];
     426        minRegion.x0 = fp->xMin - frSquareRadius;
     427        minRegion.x1 = fp->xMin + frSquareRadius;
     428        minRegion.y0 = fp->yMin - frSquareRadius;
     429        minRegion.y1 = fp->yMin + frSquareRadius;
     430        psImage *minSubImage = psImageSubset(image, minRegion);
     431        minStats = psImageStats(minStats, minSubImage, NULL, 0);
     432
     433        maxRegion.x0 = fp->xMax - frSquareRadius;
     434        maxRegion.x1 = fp->xMax + frSquareRadius;
     435        maxRegion.y0 = fp->yMax - frSquareRadius;
     436        maxRegion.y1 = fp->yMax + frSquareRadius;
     437        psImage *maxSubImage = psImageSubset(image, maxRegion);
     438        maxStats = psImageStats(maxStats, maxSubImage, NULL, 0);
     439
     440        if ((minStats == NULL) || (maxStats == NULL)) {
     441            psError(PS_ERR_UNKNOWN, true, "Could not determine robust mean on subimage.\n");
     442            psFree(minStats);
     443            psFree(maxStats);
     444            return(NULL);
     445        }
     446
     447        fp->midValue = 0.5 * (maxStats->robustMean + minStats->robustMean);
     448        fp->delta = maxStats->robustMean - minStats->robustMean;
     449        diffs->data.F32[i] = fp->delta;
     450    }
     451    psFree(minStats);
     452    psFree(maxStats);
     453
     454    diffStats = psVectorStats(diffStats, diffs, NULL, NULL, 0);
     455    psFree(diffs);
     456    if (diffStats == NULL) {
     457        psError(PS_ERR_UNKNOWN, true, "Could not determine robust median of the differences.\n");
     458        return(NULL);
     459    }
     460    return(diffStats);
     461}
     462
  • trunk/psModules/src/imcombine/pmReadoutCombine.h

    r5170 r5516  
    55 *  @author GLG, MHPCC
    66 *
    7  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-09-28 20:43:52 $
     7 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2005-11-15 20:09:03 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    4444                          float readnoise);
    4545
     46/**
     47 *
     48 * This function measures the robust median at each of the minimum and maximum
     49 * coordinates and determines the difference and mean of the two values. The size
     50 * of the box used to make the measurement at each point is specified by the
     51 * configuration variable FRINGE_SQUARE_RADIUS. From the collection of
     52 * differences, the robust median is calculated, and returned as part of the
     53 * fringe statistics. For each fringe point, the values of delta and midValue are
     54 * also assigned and available to the user on return.
     55 *
     56 */
     57psStats *pmFringeStats(
     58    psArray *fringePoints,
     59    psImage *image,
     60    psMetadata *config
     61);
     62
     63typedef struct
     64{
     65    psF64 xMin;
     66    psF64 yMin;
     67    psF64 xMax;
     68    psF64 yMax;
     69    psF64 delta;
     70    psF64 midValue;
     71}
     72pmFringePoint;
     73
    4674#endif
  • trunk/psModules/src/imsubtract/pmSubtractBias.c

    r5435 r5516  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-10-20 23:06:24 $
     8 *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2005-11-15 20:09:03 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1212 *
    1313 */
    14 
     14/*****************************************************************************/
     15/* INCLUDE FILES                                                             */
     16/*****************************************************************************/
     17#include <stdio.h>
     18#include <math.h>
     19#include <string.h>
     20#include "pslib.h"
    1521#if HAVE_CONFIG_H
    1622#include <config.h>
    1723#endif
    18 
    1924#include "pmSubtractBias.h"
    2025
    21 #define PM_SUBTRACT_BIAS_POLYNOMIAL_ORDER 2
    22 #define PM_SUBTRACT_BIAS_SPLINE_ORDER 3
    23 
     26/*****************************************************************************/
     27/* DEFINE STATEMENTS                                                         */
     28/*****************************************************************************/
    2429// XXX: put these in psConstants.h
    25 void PS_POLY1D_PRINT(psPolynomial1D *poly)
     30void PS_POLY1D_PRINT(
     31    psPolynomial1D *poly)
    2632{
    2733    printf("-------------- PS_POLY1D_PRINT() --------------\n");
     
    5157}\
    5258
    53 /******************************************************************************
    54 psSubtractFrame(): this routine will take as input a readout for the input
    55 image and a readout for the bias image.  The bias image is subtracted in
    56 place from the input image.
     59/*****************************************************************************/
     60/* TYPE DEFINITIONS                                                          */
     61/*****************************************************************************/
     62
     63/*****************************************************************************/
     64/* GLOBAL VARIABLES                                                          */
     65/*****************************************************************************/
     66psS32 currentId = 0;                // XXX: remove
     67psS32 memLeaks = 0;                 // XXX: remove
     68//PRINT_MEMLEAKS(8); XXX
     69/*****************************************************************************/
     70/* FILE STATIC VARIABLES                                                     */
     71/*****************************************************************************/
     72
     73/*****************************************************************************/
     74/* FUNCTION IMPLEMENTATION - LOCAL                                           */
     75/*****************************************************************************/
     76
     77/******************************************************************************
     78psSubtractFrame(): this routine will take as input the pmReadout for the input
     79image and a pmReadout for the bias image.  The bias image is subtracted in
     80place from the input image.  We assume that sizes and types are checked
     81elsewhere.
     82 
     83XXX: Verify that the image and readout offsets are being used the right way.
     84 
     85XXX: Ensure that it does the correct thing with image size.
    5786*****************************************************************************/
    58 static pmReadout *SubtractFrame(pmReadout *in,
    59                                 const pmReadout *bias)
    60 {
    61     psS32 i;
    62     psS32 j;
    63 
    64     if (bias == NULL) {
    65         psLogMsg(__func__, PS_LOG_WARN,
    66                  "WARNING: pmSubtractBias.c: SubtractFrame(): bias frame is NULL.  Returning original image.\n");
    67         return(in);
    68     }
    69 
    70 
    71     if ((in->image->numRows + in->row0 - bias->row0) > bias->image->numRows) {
    72         psError(PS_ERR_UNKNOWN,true, "bias image does not have enough rows.  Returning in image\n");
    73         return(in);
    74     }
    75     if ((in->image->numCols + in->col0 - bias->col0) > bias->image->numCols) {
    76         psError(PS_ERR_UNKNOWN,true, "bias image does not have enough columns.  Returning in image\n");
    77         return(in);
    78     }
    79 
    80     for (i=0;i<in->image->numRows;i++) {
    81         for (j=0;j<in->image->numCols;j++) {
     87static pmReadout *SubtractFrame(
     88    pmReadout *in,
     89    const pmReadout *bias)
     90{
     91    // XXX: When did the ->row0 and ->col0 offsets get coded?
     92    for (psS32 i=0;i<in->image->numRows;i++) {
     93        for (psS32 j=0;j<in->image->numCols;j++) {
    8294            in->image->data.F32[i][j]-=
    8395                bias->image->data.F32[i+in->row0-bias->row0][j+in->col0-bias->col0];
     96
    8497            if ((in->mask != NULL) && (bias->mask != NULL)) {
    8598                (in->mask->data.U8[i][j])|=
     
    92105}
    93106
     107
     108/******************************************************************************
     109psSubtractDarkFrame(): this routine will take as input the pmReadout for the
     110input image and a pmReadout for the dark image.  The dark image is scaled and
     111subtracted in place from the input image.
     112 
     113XXX: Verify that the image and readout offsets are being used the right way.
     114 
     115XXX: Ensure that it does the correct thing with image size.
     116*****************************************************************************/
     117static pmReadout *SubtractDarkFrame(
     118    pmReadout *in,
     119    const pmReadout *dark,
     120    psF32 scale)
     121{
     122    // XXX: When did the ->row0 and ->col0 offsets get coded?
     123    if (fabs(scale) > FLT_EPSILON) {
     124        for (psS32 i=0;i<in->image->numRows;i++) {
     125            for (psS32 j=0;j<in->image->numCols;j++) {
     126                in->image->data.F32[i][j]-=
     127                    (scale * dark->image->data.F32[i+in->row0-dark->row0][j+in->col0-dark->col0]);
     128
     129                if ((in->mask != NULL) && (dark->mask != NULL)) {
     130                    (in->mask->data.U8[i][j])|=
     131                        dark->mask->data.U8[i+in->row0-dark->row0][j+in->col0-dark->col0];
     132                }
     133            }
     134        }
     135    } else {
     136        for (psS32 i=0;i<in->image->numRows;i++) {
     137            for (psS32 j=0;j<in->image->numCols;j++) {
     138                in->image->data.F32[i][j]-=
     139                    dark->image->data.F32[i+in->row0-dark->row0][j+in->col0-dark->col0];
     140
     141                if ((in->mask != NULL) && (dark->mask != NULL)) {
     142                    (in->mask->data.U8[i][j])|=
     143                        dark->mask->data.U8[i+in->row0-dark->row0][j+in->col0-dark->col0];
     144                }
     145            }
     146        }
     147    }
     148
     149    return(in);
     150}
     151
    94152/******************************************************************************
    95153ImageSubtractScalar(): subtract a scalar from the input image.
    96154 
    97 XXX: Use a psLib function for this.
    98  
    99 XXX: This should
    100  *****************************************************************************/
    101 static psImage *ImageSubtractScalar(psImage *image,
    102                                     psF32 scalar)
     155XXX: Is there a psLib function for this?
     156 *****************************************************************************/
     157static psImage *ImageSubtractScalar(
     158    psImage *image,
     159    psF32 scalar)
    103160{
    104161    for (psS32 i=0;i<image->numRows;i++) {
     
    164221
    165222    if (numOptions == 0) {
    166         psError(PS_ERR_UNKNOWN,true, "No statistics options have been specified.\n");
     223        psError(PS_ERR_UNKNOWN,true, "No allowable statistics options have been specified.\n");
    167224    }
    168225    if (numOptions != 1) {
     
    173230}
    174231
    175 
     232/******************************************************************************
     233Polynomial1DCopy(): This private function copies the members of the existing
     234psPolynomial1D "in" into the existing psPolynomial1D "out".  The previous
     235members of the existing psPolynomial1D "out" are psFree'ed.
     236 *****************************************************************************/
     237static psBool Polynomial1DCopy(
     238    psPolynomial1D *out,
     239    psPolynomial1D *in)
     240{
     241    psFree(out->coeff);
     242    psFree(out->coeffErr);
     243    psFree(out->mask);
     244
     245    out->type = in->type;
     246    out->nX = in->nX;
     247
     248    out->coeff = (psF64 *) psAlloc((in->nX + 1) * sizeof(psF64));
     249    // XXX: use memcpy
     250    for (psS32 i = 0 ; i < (in->nX + 1) ; i++) {
     251        out->coeff[i] = in->coeff[i];
     252    }
     253
     254    out->coeffErr = (psF64 *) psAlloc((in->nX + 1) * sizeof(psF64));
     255    // XXX: use memcpy
     256    for (psS32 i = 0 ; i < (in->nX + 1) ; i++) {
     257        out->coeffErr[i] = in->coeffErr[i];
     258    }
     259
     260    out->mask = (psMaskType *) psAlloc((in->nX + 1) * sizeof(psMaskType));
     261    // XXX: use memcpy
     262    for (psS32 i = 0 ; i < (in->nX + 1) ; i++) {
     263        out->mask[i] = in->mask[i];
     264    }
     265
     266    return(true);
     267}
     268
     269/******************************************************************************
     270Polynomial1DDup(): This private function duplicates and then returns the input
     271psPolynomial1D "in".
     272 *****************************************************************************/
     273static psPolynomial1D *Polynomial1DDup(
     274    psPolynomial1D *in)
     275{
     276    psPolynomial1D *out = psPolynomial1DAlloc(in->nX, in->type);
     277    Polynomial1DCopy(out, in);
     278    return(out);
     279}
     280
     281
     282/******************************************************************************
     283SplineCopy(): This private function copies the members of the existing
     284psSpline in into the existing psSpline out.
     285 *****************************************************************************/
     286static psBool SplineCopy(
     287    psSpline1D *out,
     288    psSpline1D *in)
     289{
     290    PS_ASSERT_PTR_NON_NULL(out, false);
     291    PS_ASSERT_PTR_NON_NULL(in, false);
     292
     293    for (psS32 i = 0 ; i < out->n ; i++) {
     294        psFree(out->spline[i]);
     295    }
     296    psFree(out->spline);
     297    psFree(out->knots);
     298    psFree(out->p_psDeriv2);
     299
     300    out->n = in->n;
     301    out->spline = (psPolynomial1D **) psAlloc(in->n * sizeof(psPolynomial1D *));
     302    for (psS32 i = 0 ; i < in->n ; i++) {
     303        out->spline[i] = Polynomial1DDup(in->spline[i]);
     304    }
     305
     306    out->knots = psVectorCopy(out->knots, in->knots, in->knots->type.type);
     307
     308    out->p_psDeriv2 = (psF32 *) psAlloc((in->n + 1) * sizeof(psF32));
     309    // XXX: use memcpy
     310    for (psS32 i = 0 ; i < (in->n + 1) ; i++) {
     311        out->p_psDeriv2[i] = in->p_psDeriv2[i];
     312    }
     313
     314    return(true);
     315}
    176316
    177317/******************************************************************************
     
    181321    PM_FIT_POLYNOMIAL: fit a polynomial to the entire input vector data.
    182322    PM_FIT_SPLINE: fit splines to the input vector data.
    183 XXX: Doesn't it make more sense to do polynomial interpolation on a few
    184 elements of the input vector, rather than fit a polynomial to the entire
    185 vector?
    186  *****************************************************************************/
    187 static psVector *ScaleOverscanVector(psVector *overscanVector,
    188                                      psS32 n,
    189                                      void *fitSpec,
    190                                      pmFit fit)
     323The resulting spline or polynomial is set in the fitSpec argument.
     324 *****************************************************************************/
     325static psVector *ScaleOverscanVector(
     326    psVector *overscanVector,
     327    psS32 n,
     328    void *fitSpec,
     329    pmFit fit)
    191330{
    192331    psTrace(".psModule.pmSubtracBias.ScaleOverscanVector", 4,
    193332            "---- ScaleOverscanVector() begin (%d -> %d) ----\n", overscanVector->n, n);
    194     //    PS_VECTOR_PRINT_F32(overscanVector);
    195333
    196334    if (NULL == overscanVector) {
     
    205343    //
    206344    if (n == overscanVector->n) {
    207         for (psS32 i = 0 ; i < n ; i++) {
    208             newVec->data.F32[i] = overscanVector->data.F32[i];
    209         }
    210         return(newVec);
    211     }
    212     psPolynomial1D *myPoly;
    213     psSpline1D *mySpline;
     345        return(psVectorCopy(newVec, overscanVector, PS_TYPE_F32));
     346    }
    214347    psF32 x;
    215     psS32 i;
    216348
    217349    if (fit == PM_FIT_POLYNOMIAL) {
    218350        // Fit a polynomial to the old overscan vector.
    219         myPoly = (psPolynomial1D *) fitSpec;
     351        psPolynomial1D *myPoly = (psPolynomial1D *) fitSpec;
    220352        PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
     353        PS_ASSERT_POLY1D(myPoly, NULL);
    221354        myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL);
    222355        if (myPoly == NULL) {
    223             psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(1): Could not fit a polynomial to the psVector.\n");
     356            psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the psVector.\n");
    224357            return(NULL);
    225358        }
     
    228361        // of the old vector, use the fitted polynomial to determine the
    229362        // interpolated value at that point, and set the new vector.
    230         for (i=0;i<n;i++) {
     363        for (psS32 i=0;i<n;i++) {
    231364            x = ((psF32) i) * ((psF32) overscanVector->n) / ((psF32) n);
    232365            newVec->data.F32[i] = psPolynomial1DEval(myPoly, x);
    233366        }
    234367    } else if (fit == PM_FIT_SPLINE) {
    235         psS32 mustFreeSpline = 0;
    236         // Fit a spline to the old overscan vector.
    237         mySpline = (psSpline1D *) fitSpec;
    238         // XXX: Does it make any sense to have a psSpline argument?
    239         if (mySpline == NULL) {
    240             mustFreeSpline = 1;
    241         }
    242 
    243368        //
    244369        // NOTE: Since the X arg in the psVectorFitSpline1D() function is NULL,
     
    246371        // properly when doing the spline eval.
    247372        //
    248         //        mySpline = psVectorFitSpline1D(mySpline, NULL, overscanVector, NULL);
    249         mySpline = psVectorFitSpline1D(NULL, overscanVector);
     373        psSpline1D *mySpline = psVectorFitSpline1D(NULL, overscanVector);
    250374        if (mySpline == NULL) {
    251             psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(2): Could not fit a spline to the psVector.\n");
     375            psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to the psVector.\n");
    252376            return(NULL);
    253377        }
    254         //        PS_PRINT_SPLINE(mySpline);
    255378
    256379        // For each element of the new vector, convert the x-ordinate to that
    257         // of the old vector, use the fitted polynomial to determine the
     380        // of the old vector, use the fitted spline to determine the
    258381        // interpolated value at that point, and set the new vector.
    259         for (i=0;i<n;i++) {
     382        for (psS32 i=0;i<n;i++) {
    260383            // Scale to [0 : overscanVector->n - 1]
    261384            x = ((psF32) i) * ((psF32) (overscanVector->n-1)) / ((psF32) n);
    262385            newVec->data.F32[i] = psSpline1DEval(mySpline, x);
    263386        }
    264         if (mustFreeSpline ==1) {
    265             psFree(mySpline);
    266         }
    267         //        PS_VECTOR_PRINT_F32(newVec);
    268 
    269 
     387
     388        psSpline1D *ptrSpline = (psSpline1D *) fitSpec;
     389        if (ptrSpline != NULL) {
     390            // Copy the resulting spline fit into ptrSpline.
     391            PS_ASSERT_SPLINE(ptrSpline, NULL);
     392            SplineCopy(ptrSpline, mySpline);
     393        }
     394        psFree(mySpline);
    270395    } else {
    271396        psError(PS_ERR_UNKNOWN, true, "unknown fit type.  Returning NULL.\n");
     
    280405
    281406/******************************************************************************
    282 XXX: The SDRS does not specify type support.  F32 is implemented here.
     407 *****************************************************************************/
     408static psS32 GetOverscanSize(
     409    psImage *inImg,
     410    pmOverscanAxis overScanAxis)
     411{
     412    if (overScanAxis == PM_OVERSCAN_ROWS) {
     413        return(inImg->numCols);
     414    } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
     415        return(inImg->numRows);
     416    } else if (overScanAxis == PM_OVERSCAN_ALL) {
     417        return(1);
     418    }
     419    return(0);
     420}
     421
     422/******************************************************************************
     423GetOverscanAxis(in) this private routine determines the appropiate overscan
     424axis from the parent cell metadata.
     425 
     426XXX: Verify the READDIR corresponds with my overscan axis.
     427 *****************************************************************************/
     428static pmOverscanAxis GetOverscanAxis(pmReadout *in)
     429{
     430    psBool rc;
     431    if ((in->parent != NULL) && (in->parent->concepts)) {
     432        psS32 dir = psMetadataLookupS32(&rc, in->parent->concepts, "CELL.READDIR");
     433        if (rc == true) {
     434            if (dir == 1) {
     435                return(PM_OVERSCAN_ROWS);
     436            } else if (dir == 2) {
     437                return(PM_OVERSCAN_COLUMNS);
     438            } else if (dir == 3) {
     439                return(PM_OVERSCAN_ALL);
     440            }
     441        }
     442    }
     443
     444    psLogMsg(__func__, PS_LOG_WARN,
     445             "WARNING: pmSubtractBias.(): could not determine CELL.READDIR from in->parent metadata.  Setting overscan axis to PM_OVERSCAN_NONE.\n");
     446    return(PM_OVERSCAN_NONE);
     447}
     448
     449/******************************************************************************
     450psListLength(list): determine the length of a psList.
     451 
     452XXX: Put this elsewhere.
     453 *****************************************************************************/
     454static psS32 psListLength(
     455    psList *list)
     456{
     457    psS32 length = 0;
     458    psListElem *tmpElem = (psListElem *) list->head;
     459    while (NULL != tmpElem) {
     460        tmpElem = tmpElem->next;
     461        length++;
     462    }
     463    return(length);
     464}
     465
     466/******************************************************************************
     467Note: this isn't needed anymore as of psModule SDRS 12-09.
     468 *****************************************************************************/
     469static psBool OverscanReducePixel(
     470    psImage *in,
     471    psList *bias,
     472    psStats *myStats)
     473{
     474    PS_ASSERT_PTR_NON_NULL(in, NULL);
     475    PS_ASSERT_PTR_NON_NULL(bias, NULL);
     476    PS_ASSERT_PTR_NON_NULL(bias->head, NULL);
     477    PS_ASSERT_PTR_NON_NULL(myStats, NULL);
     478
     479    // Allocate a psVector with one element per overscan image.
     480    psS32 numOverscanImages = psListLength(bias);
     481    psVector *statsAll = psVectorAlloc(numOverscanImages, PS_TYPE_F32);
     482    psListElem *tmpOverscan = (psListElem *) bias->head;
     483    psS32 i = 0;
     484    psF64 statValue;
     485    //
     486    // We loop through each overscan image, calculating the specified
     487    // statistic on that image.
     488    //
     489    while (NULL != tmpOverscan) {
     490        psImage *myOverscanImage = (psImage *) tmpOverscan->data;
     491
     492        PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL);
     493        myStats = psImageStats(myStats, myOverscanImage, NULL, (psMaskType)0xffffffff);
     494        if (myStats == NULL) {
     495            psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation.  Returning in image.\n");
     496            psFree(statsAll);
     497            return(false);
     498        }
     499        if (false == p_psGetStatValue(myStats, &statValue)) {
     500            psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
     501            psFree(statsAll);
     502            return(false);
     503        }
     504        statsAll->data.F32[i] = statValue;
     505        i++;
     506        tmpOverscan = tmpOverscan->next;
     507    }
     508
     509    //
     510    // We reduce the individual stats for each overscan image to
     511    // a single psF32.
     512    //
     513    myStats = psVectorStats(myStats, statsAll, NULL, NULL, 0);
     514    if (myStats == NULL) {
     515        psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation.  Returning in image.\n");
     516        psFree(statsAll);
     517        return(false);
     518    }
     519    if (false == p_psGetStatValue(myStats, &statValue)) {
     520        psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
     521        psFree(statsAll);
     522        return(false);
     523    }
     524
     525    //
     526    // Subtract the result and return.
     527    //
     528    ImageSubtractScalar(in, statValue);
     529    psFree(statsAll);
     530    return(in);
     531}
     532
     533/******************************************************************************
     534ReduceOverscanImageToCol(overscanImage, myStats): This private routine reduces
     535a single psImage to a column by combining all pixels from each row into a
     536single pixel via requested statistic in myStats.
     537 *****************************************************************************/
     538static psVector *ReduceOverscanImageToCol(
     539    psImage *overscanImage,
     540    psStats *myStats)
     541{
     542    psF64 statValue;
     543    psVector *tmpRow = psVectorAlloc(overscanImage->numCols, PS_TYPE_F32);
     544    psVector *tmpCol = psVectorAlloc(overscanImage->numRows, PS_TYPE_F32);
     545
     546    //
     547    // For each row, we store all pixels in that row into a temporary psVector,
     548    // then we run psVectorStats() on that vector.
     549    //
     550    for (psS32 i=0;i<overscanImage->numRows;i++) {
     551        for (psS32 j=0;j<overscanImage->numCols;j++) {
     552            tmpRow->data.F32[j] = overscanImage->data.F32[i][j];
     553        }
     554
     555        psStats *rc = psVectorStats(myStats, tmpRow, NULL, NULL, 0);
     556        if (rc == NULL) {
     557            psError(PS_ERR_UNKNOWN, true, "psVectorStats() could not perform requested statistical operation.  Returning in image.\n");
     558            return(NULL);
     559        }
     560
     561        if (false ==  p_psGetStatValue(rc, &statValue)) {
     562            psError(PS_ERR_UNKNOWN, true, "p_psGetStatValue() could not determine result from requested statistical operation.  Returning in image.\n");
     563            return(NULL);
     564        }
     565
     566        tmpCol->data.F32[i] = (psF32) statValue;
     567    }
     568    psFree(tmpRow);
     569
     570    return(tmpCol);
     571}
     572
     573/******************************************************************************
     574ReduceOverscanImageToCol(overscanImage, myStats): This private routine reduces
     575a single psImage to a row by combining all pixels from each column into a
     576single pixel via requested statistic in myStats.
     577 *****************************************************************************/
     578static psVector *ReduceOverscanImageToRow(
     579    psImage *overscanImage,
     580    psStats *myStats)
     581{
     582    psF64 statValue;
     583    psVector *tmpRow = psVectorAlloc(overscanImage->numCols, PS_TYPE_F32);
     584    psVector *tmpCol = psVectorAlloc(overscanImage->numRows, PS_TYPE_F32);
     585
     586    //
     587    // For each column, we store all pixels in that column into a temporary psVector,
     588    // then we run psVectorStats() on that vector.
     589    //
     590    for (psS32 i=0;i<overscanImage->numCols;i++) {
     591        for (psS32 j=0;j<overscanImage->numRows;j++) {
     592            tmpCol->data.F32[j] = overscanImage->data.F32[j][i];
     593        }
     594
     595        psStats *rc = psVectorStats(myStats, tmpCol, NULL, NULL, 0);
     596        if (rc == NULL) {
     597            psError(PS_ERR_UNKNOWN, true, "psVectorStats() could not perform requested statistical operation.  Returning in image.\n");
     598            return(NULL);
     599        }
     600
     601        if (false ==  p_psGetStatValue(rc, &statValue)) {
     602            psError(PS_ERR_UNKNOWN, true, "p_psGetStatValue() could not determine result from requested statistical operation.  Returning in image.\n");
     603            return(NULL);
     604        }
     605
     606        tmpRow->data.F32[i] = (psF32) statValue;
     607    }
     608    psFree(tmpCol);
     609
     610    return(tmpRow);
     611}
     612
     613/******************************************************************************
     614OverscanReduce(vecSize, bias, myStats): This private routine takes a psList of
     615overscan images (in bias) and reduces them to a single psVector via the
     616specified psStats struct.  The vector is then scaled to the length or the
     617row/column in inImg.
     618 *****************************************************************************/
     619static psVector* OverscanReduce(
     620    psImage *inImg,
     621    pmOverscanAxis overScanAxis,
     622    psList *bias,
     623    void *fitSpec,
     624    pmFit fit,
     625    psStats *myStats)
     626{
     627    if ((overScanAxis != PM_OVERSCAN_ROWS) && (overScanAxis != PM_OVERSCAN_COLUMNS)) {
     628        psError(PS_ERR_UNKNOWN, true, "overScanAxis must be PM_OVERSCAN_ROWS or PM_OVERSCAN_COLUMNS\n");
     629        return(NULL);
     630    }
     631    PS_ASSERT_PTR_NON_NULL(inImg, NULL);
     632    PS_ASSERT_PTR_NON_NULL(bias, NULL);
     633    PS_ASSERT_PTR_NON_NULL(bias->head, NULL);
     634    PS_ASSERT_PTR_NON_NULL(myStats, NULL);
     635    //
     636    // Allocate a psVector for the output of this routine.
     637    //
     638    psS32 vecSize = GetOverscanSize(inImg, overScanAxis);
     639    psVector *overscanVector = psVectorAlloc(vecSize, PS_TYPE_F32);
     640
     641    //
     642    // Allocate an array of psVectors with one psVector per element of the
     643    // final oversan column vector.  These psVectors will be used with
     644    // psStats to reduce the multiple elements from each overscan column
     645    // vector to a single final column vector.
     646    //
     647    psS32 numOverscanImages = psListLength(bias);
     648    psVector **overscanVectors = (psVector **) psAlloc(numOverscanImages * sizeof(psVector *));
     649    for (psS32 i = 0 ; i < numOverscanImages ; i++) {
     650        overscanVectors[i] = NULL;
     651    }
     652
     653    //
     654    // We iterate through the list of overscan images.  For each image,
     655    // we reduce it to a single column or row.  Save the overscan vector
     656    // in overscanVectors[].
     657    //
     658    psListElem *tmpOverscan = (psListElem *) bias->head;
     659    psS32 overscanID = 0;
     660    while (tmpOverscan != NULL) {
     661        psImage *tmpOverscanImage = (psImage *) tmpOverscan->data;
     662        if (overScanAxis == PM_OVERSCAN_ROWS) {
     663            overscanVectors[overscanID] = ReduceOverscanImageToRow(tmpOverscanImage, myStats);
     664        } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
     665            overscanVectors[overscanID] = ReduceOverscanImageToCol(tmpOverscanImage, myStats);
     666        }
     667
     668        tmpOverscan = tmpOverscan->next;
     669        overscanID++;
     670    }
     671
     672    //
     673    // For each overscan vector, if necessary, we scale that column or
     674    // row to vecSize.  Note: we should have already ensured that the
     675    // fit is poly or spline.
     676    //
     677    for (psS32 i = 0 ; i < numOverscanImages ; i++) {
     678        psVector *tmpOverscanVector = overscanVectors[i];
     679
     680        if (tmpOverscanVector->n != vecSize) {
     681            overscanVectors[i] = ScaleOverscanVector(tmpOverscanVector, vecSize, fitSpec, fit);
     682            if (overscanVectors[i] == NULL) {
     683                psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector.\n");
     684                for (psS32 i = 0 ; i < numOverscanImages ; i++) {
     685                    psFree(overscanVectors[i]);
     686                }
     687                psFree(overscanVectors);
     688                psFree(tmpOverscanVector);
     689                return(NULL);
     690            }
     691            psFree(tmpOverscanVector);
     692        }
     693    }
     694
     695    //
     696    // We collect all elements in the overscan vectors for the various
     697    // overscan images into a single psVector (tmpVec).  Then we call
     698    // psStats on that vector to determine the final values for the
     699    // overscan vector.
     700    //
     701    psVector *tmpVec = psVectorAlloc(numOverscanImages, PS_TYPE_F32);
     702    psF64 statValue;
     703    for (psS32 i = 0 ; i < vecSize ; i++) {
     704        // Collect the i-th elements from each overscan vector into a single vector.
     705        for (psS32 j = 0 ; j < numOverscanImages ; j++) {
     706            tmpVec->data.F32[j] = overscanVectors[j]->data.F32[i];
     707        }
     708
     709        if (NULL == psVectorStats(myStats, tmpVec, NULL, NULL, 0)) {
     710            psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
     711            for (psS32 i = 0 ; i < numOverscanImages ; i++) {
     712                psFree(overscanVectors[i]);
     713            }
     714            psFree(overscanVectors);
     715            psFree(tmpVec);
     716            return(NULL);
     717        }
     718        if (false == p_psGetStatValue(myStats, &statValue)) {
     719            psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
     720            for (psS32 i = 0 ; i < numOverscanImages ; i++) {
     721                psFree(overscanVectors[i]);
     722            }
     723            psFree(overscanVectors);
     724            psFree(tmpVec);
     725            return(NULL);
     726        }
     727
     728        overscanVector->data.F32[i] = (psF32) statValue;
     729    }
     730
     731    //
     732    // We're done.  Free the intermediate overscan vectors.
     733    //
     734    psFree(tmpVec);
     735    for (psS32 i = 0 ; i < numOverscanImages ; i++) {
     736        psFree(overscanVectors[i]);
     737    }
     738    psFree(overscanVectors);
     739
     740    //
     741    // Return the computed overscanVector
     742    //
     743    return(overscanVector);
     744}
     745
     746/******************************************************************************
     747RebinOverscanVector(overscanVector, nBinOrig, myStats): this private routine
     748takes groups of nBinOrig elements in the input vector, combines them into a
     749single pixel via myStats and psVectorStats(), and then outputs a vector of
     750those pixels.
     751 *****************************************************************************/
     752static psS32 RebinOverscanVector(
     753    psVector *overscanVector,
     754    psS32 nBinOrig,
     755    psStats *myStats)
     756{
     757    psF64 statValue;
     758    psS32 nBin;
     759    if ((nBinOrig > 1) && (nBinOrig < overscanVector->n)) {
     760        psS32 numBins = 1+((overscanVector->n)/nBinOrig);
     761        psVector *myBin = psVectorAlloc(numBins, PS_TYPE_F32);
     762        psVector *binVec = psVectorAlloc(nBinOrig, PS_TYPE_F32);
     763
     764        for (psS32 i=0;i<numBins;i++) {
     765            for(psS32 j=0;j<nBinOrig;j++) {
     766                if (overscanVector->n > ((i*nBinOrig)+j)) {
     767                    binVec->data.F32[j] = overscanVector->data.F32[(i*nBinOrig)+j];
     768                } else {
     769                    // XXX: we get here if nBinOrig does not evenly divide
     770                    // the overscanVector vector.  This is the last bin.  Should
     771                    // we change the binVec->n to acknowledge that?
     772                    binVec->n = j;
     773                }
     774            }
     775            psStats *rc = psVectorStats(myStats, binVec, NULL, NULL, 0);
     776            if (rc == NULL) {
     777                psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
     778                return(-1);
     779            }
     780            if (false ==  p_psGetStatValue(rc, &statValue)) {
     781                psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
     782                return(-1);
     783            }
     784            myBin->data.F32[i] = statValue;
     785        }
     786
     787        // Change the effective size of overscanVector.
     788        overscanVector->n = numBins;
     789        for (psS32 i=0;i<numBins;i++) {
     790            overscanVector->data.F32[i] = myBin->data.F32[i];
     791        }
     792        psFree(binVec);
     793        psFree(myBin);
     794        nBin = nBinOrig;
     795    } else {
     796        nBin = 1;
     797    }
     798
     799    return(nBin);
     800}
     801
     802/******************************************************************************
     803FitOverscanVectorAndUnbin(inImg, overscanVector, overScanAxis, fitSpec, fit,
     804nBin):  this private routine fits a psPolynomial or psSpline to the overscan
     805vector.  It then creates a new vector, with a size determined by the input
     806image, evaluates the psPolynomial or psSpline at each element in that vector,
     807then returns that vector.
     808 *****************************************************************************/
     809static psVector *FitOverscanVectorAndUnbin(
     810    psImage *inImg,
     811    psVector *overscanVector,
     812    pmOverscanAxis overScanAxis,
     813    void *fitSpec,
     814    pmFit fit,
     815    psS32 nBin)
     816{
     817    psPolynomial1D* myPoly = NULL;
     818    psSpline1D *mySpline = NULL;
     819    //
     820    // Fit a polynomial or spline to the overscan vector.
     821    //
     822    if (fit == PM_FIT_POLYNOMIAL) {
     823        myPoly = (psPolynomial1D *) fitSpec;
     824        PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
     825        PS_ASSERT_POLY1D(myPoly, NULL);
     826        myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL);
     827        if (myPoly == NULL) {
     828            psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to overscan vector.  Returning NULL.\n");
     829            return(NULL);
     830        }
     831    } else if (fit == PM_FIT_SPLINE) {
     832        mySpline = psVectorFitSpline1D(NULL, overscanVector);
     833        if (mySpline == NULL) {
     834            psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to overscan vector.  Returning NULL.\n");
     835            return(NULL);
     836        }
     837        if (fitSpec != NULL) {
     838            // Copy the resulting spline fit into fitSpec.
     839            psSpline1D *ptrSpline = (psSpline1D *) fitSpec;
     840            PS_ASSERT_SPLINE(ptrSpline, NULL);
     841            SplineCopy(ptrSpline, mySpline);
     842        }
     843    }
     844
     845    //
     846    // Evaluate the poly/spline at each pixel in the overscan row/column.
     847    //
     848    psS32 vecSize = GetOverscanSize(inImg, overScanAxis);
     849    psVector *newVec = psVectorAlloc(vecSize, PS_TYPE_F32);
     850    if ((nBin > 1) && (nBin < overscanVector->n)) {
     851        for (psS32 i = 0 ; i < vecSize ; i++) {
     852            if (fit == PM_FIT_POLYNOMIAL) {
     853                newVec->data.F32[i] = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin));
     854            } else if (fit == PM_FIT_SPLINE) {
     855                newVec->data.F32[i] = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin));
     856            }
     857        }
     858    } else {
     859        for (psS32 i = 0 ; i < vecSize ; i++) {
     860            if (fit == PM_FIT_POLYNOMIAL) {
     861                newVec->data.F32[i] = psPolynomial1DEval(myPoly, (psF32) i);
     862            } else if (fit == PM_FIT_SPLINE) {
     863                newVec->data.F32[i] = psSpline1DEval(mySpline, (psF32) i);
     864            }
     865        }
     866    }
     867
     868    psFree(mySpline);
     869    psFree(overscanVector);
     870    return(newVec);
     871}
     872
     873
     874
     875/******************************************************************************
     876UnbinOverscanVector(inImg, overscanVector, overScanAxis, nBin):  this private
     877routine takes a psVector overscanVector that was previously binned by a factor
     878of nBin, and then expands it to its original size, duplicated elements nBin
     879times for each element in the input vector overscanVector.
     880 *****************************************************************************/
     881static psVector *UnbinOverscanVector(
     882    psImage *inImg,
     883    psVector *overscanVector,
     884    pmOverscanAxis overScanAxis,
     885    psS32 nBin)
     886{
     887    psS32 vecSize;
     888
     889    if (overScanAxis == PM_OVERSCAN_ROWS) {
     890        vecSize = inImg->numCols;
     891    } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
     892        vecSize = inImg->numRows;
     893    }
     894
     895    psVector *newVec = psVectorAlloc(vecSize, PS_TYPE_F32);
     896    for (psS32 i = 0 ; i < vecSize ; i++) {
     897        newVec->data.F32[i] = overscanVector->data.F32[i/nBin];
     898    }
     899
     900    psFree(overscanVector);
     901    return(newVec);
     902}
     903
     904
     905/******************************************************************************
     906SubtractVectorFromImage(inImg, overscanVector, overScanAxis):  this private
     907routine subtracts the overscanVector column-wise or row-wise from inImg.
     908 *****************************************************************************/
     909static psImage *SubtractVectorFromImage(
     910    psImage *inImg,
     911    psVector *overscanVector,
     912    pmOverscanAxis overScanAxis)
     913{
     914    //
     915    // Subtract overscan vector row-wise from the image.
     916    //
     917    if (overScanAxis == PM_OVERSCAN_ROWS) {
     918        for (psS32 i=0;i<inImg->numCols;i++) {
     919            for (psS32 j=0;j<inImg->numRows;j++) {
     920                inImg->data.F32[j][i]-= overscanVector->data.F32[i];
     921            }
     922        }
     923    }
     924
     925    //
     926    // Subtract overscan vector column-wise from the image.
     927    //
     928    if (overScanAxis == PM_OVERSCAN_COLUMNS) {
     929        for (psS32 i=0;i<inImg->numRows;i++) {
     930            for (psS32 j=0;j<inImg->numCols;j++) {
     931                inImg->data.F32[i][j]-= overscanVector->data.F32[i];
     932            }
     933        }
     934    }
     935
     936    return(inImg);
     937}
     938
     939
     940
     941typedef enum {
     942    PM_ERROR_NO_SUBTRACTION,
     943    PM_WARNING_NO_SUBTRACTION,
     944    PM_ERROR_NO_BIAS_SUBTRACT,
     945    PM_WARNING_NO_BIAS_SUBTRACT,
     946    PM_ERROR_NO_DARK_SUBTRACT,
     947    PM_WARNING_NO_DARK_SUBTRACT,
     948    PM_OKAY
     949} pmSubtractBiasAssertStatus;
     950/******************************************************************************
     951AssertCodeOverscan(....) this private routine verifies that the various input
     952parameters to pmSubtractBias() are correct for overscan subtraction.
     953 *****************************************************************************/
     954pmSubtractBiasAssertStatus AssertCodeOverscan(
     955    pmReadout *in,
     956    void *fitSpec,
     957    pmFit fit,
     958    bool overscan,
     959    psStats *stat,
     960    int nBinOrig,
     961    const pmReadout *bias,
     962    const pmReadout *dark)
     963{
     964
     965    PS_ASSERT_READOUT_NON_NULL(in, PM_ERROR_NO_SUBTRACTION);
     966    PS_ASSERT_READOUT_NON_EMPTY(in, PM_ERROR_NO_SUBTRACTION);
     967    PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, PM_ERROR_NO_SUBTRACTION);
     968    PS_WARN_PTR_NON_NULL(in->parent);
     969    if (in->parent != NULL) {
     970        PS_WARN_PTR_NON_NULL(in->parent->concepts);
     971    }
     972
     973    if (overscan == true) {
     974        pmOverscanAxis overScanAxis = GetOverscanAxis(in);
     975        PS_ASSERT_PTR_NON_NULL(stat, PM_ERROR_NO_SUBTRACTION);
     976        PS_ASSERT_PTR_NON_NULL(in->bias, PM_ERROR_NO_SUBTRACTION);
     977        PS_ASSERT_PTR_NON_NULL(in->bias->head, PM_ERROR_NO_SUBTRACTION);
     978        //
     979        // Check the type, size of each bias image.
     980        //
     981        psListElem *tmpOverscan = (psListElem *) in->bias->head;
     982        psS32 numOverscans = 0;
     983        while (NULL != tmpOverscan) {
     984            numOverscans++;
     985            psImage *myOverscanImage = (psImage *) tmpOverscan->data;
     986            PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, PM_ERROR_NO_SUBTRACTION);
     987            // XXX: Get this right with the rows and columns.
     988            if (overScanAxis == PM_OVERSCAN_ROWS) {
     989                if (myOverscanImage->numRows != in->image->numRows) {
     990                    psLogMsg(__func__, PS_LOG_WARN,
     991                             "WARNING: pmSubtractBias.(): overscan image (# %d) has %d rows, input image has %d rows\n",
     992                             numOverscans, myOverscanImage->numCols, in->image->numRows);
     993                    if (fit == PM_FIT_NONE) {
     994                        psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vectors.  Set fit to PM_FIT_POLYNOMIAL or PM_FIT_SPLINE.\n");
     995                        return(PM_ERROR_NO_SUBTRACTION);
     996                    }
     997                }
     998            } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
     999                if (myOverscanImage->numCols != in->image->numCols) {
     1000                    psLogMsg(__func__, PS_LOG_WARN,
     1001                             "WARNING: pmSubtractBias.(): overscan image (# %d) has %d columns, input image has %d columns\n",
     1002                             numOverscans, myOverscanImage->numCols, in->image->numCols);
     1003                    if (fit == PM_FIT_NONE) {
     1004                        psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vectors.  Set fit to PM_FIT_POLYNOMIAL or PM_FIT_SPLINE.\n");
     1005                        return(PM_ERROR_NO_SUBTRACTION);
     1006                    }
     1007                }
     1008            } else if (overScanAxis != PM_OVERSCAN_ALL) {
     1009                psError(PS_ERR_UNKNOWN, true, "Must specify and overscan axis.\n");
     1010                return(PM_ERROR_NO_SUBTRACTION);
     1011            }
     1012            tmpOverscan = tmpOverscan->next;
     1013        }
     1014    } else {
     1015        if (fit != PM_FIT_NONE) {
     1016            psLogMsg(__func__, PS_LOG_WARN,
     1017                     "WARNING: pmSubtractBias.(): overscan is FALSE and fit is not PM_FIT_NONE.\n");
     1018            return(PM_WARNING_NO_SUBTRACTION);
     1019        }
     1020    }
     1021
     1022    // XXX: I do not like the following spec since it's useless to specify
     1023    // a psSpline as the fitSpec.
     1024    if (0) {
     1025        if ((fitSpec == NULL) &&
     1026                ((fit != PM_FIT_NONE) || (overscan == true))) {
     1027            psError(PS_ERR_UNKNOWN, true, "fitSpec is NULL and fit is not PM_FIT_NONE or overscan is TRUE.\n");
     1028            return(PM_ERROR_NO_SUBTRACTION);
     1029        }
     1030    }
     1031
     1032    return(PM_OKAY);
     1033}
     1034
     1035/******************************************************************************
     1036AssertCodeBias(....) this private routine verifies that the various input
     1037parameters to pmSubtractBias() are correct for bias subtraction.
     1038 *****************************************************************************/
     1039static pmSubtractBiasAssertStatus AssertCodeBias(
     1040    pmReadout *in,
     1041    void *fitSpec,
     1042    pmFit fit,
     1043    bool overscan,
     1044    psStats *stat,
     1045    int nBinOrig,
     1046    const pmReadout *bias,
     1047    const pmReadout *dark)
     1048{
     1049    if ((in->image->numRows + in->row0 - bias->row0) > bias->image->numRows) {
     1050        psError(PS_ERR_UNKNOWN,true, "bias image does not have enough rows.  Returning in image\n");
     1051        return(PM_ERROR_NO_BIAS_SUBTRACT);
     1052    }
     1053    if ((in->image->numCols + in->col0 - bias->col0) > bias->image->numCols) {
     1054        psError(PS_ERR_UNKNOWN,true, "bias image does not have enough columns.  Returning in image\n");
     1055        return(PM_ERROR_NO_BIAS_SUBTRACT);
     1056    }
     1057
     1058    if (bias != NULL) {
     1059        PS_ASSERT_READOUT_NON_EMPTY(bias, PM_ERROR_NO_BIAS_SUBTRACT);
     1060        PS_ASSERT_READOUT_TYPE(bias, PS_TYPE_F32, PM_ERROR_NO_DARK_SUBTRACT);
     1061    }
     1062    return(PM_OKAY);
     1063}
     1064
     1065/******************************************************************************
     1066AssertCodeDark(....) this private routine verifies that the various input
     1067parameters to pmSubtractBias() are correct for dark subtraction.
     1068 *****************************************************************************/
     1069pmSubtractBiasAssertStatus AssertCodeDark(
     1070    pmReadout *in,
     1071    void *fitSpec,
     1072    pmFit fit,
     1073    bool overscan,
     1074    psStats *stat,
     1075    int nBinOrig,
     1076    const pmReadout *bias,
     1077    const pmReadout *dark)
     1078{
     1079    if ((in->image->numRows + in->row0 - dark->row0) > dark->image->numRows) {
     1080        psError(PS_ERR_UNKNOWN, true, "dark image does not have enough rows.  Returning in image\n");
     1081        return(PM_ERROR_NO_DARK_SUBTRACT);
     1082    }
     1083    if ((in->image->numCols + in->col0 - dark->col0) > dark->image->numCols) {
     1084        psError(PS_ERR_UNKNOWN, true, "dark image does not have enough columns.  Returning in image\n");
     1085        return(PM_ERROR_NO_DARK_SUBTRACT);
     1086    }
     1087
     1088    if (dark != NULL) {
     1089        PS_ASSERT_READOUT_NON_EMPTY(dark, PM_ERROR_NO_DARK_SUBTRACT);
     1090        PS_ASSERT_READOUT_TYPE(dark, PS_TYPE_F32, PM_ERROR_NO_DARK_SUBTRACT);
     1091    }
     1092    return(PM_OKAY);
     1093}
     1094
     1095/******************************************************************************
     1096p_psDetermineTrimmedImage(): global routine: determines the region of the
     1097input pmReadout which will be operated on by the various detrend modules.  It
     1098does a metadata fetch on "CELL.TRIMSEC" for the parent cell of the pmReadout.
     1099 
     1100Use it this way:
     1101    PS_WARN_PTR_NON_NULL(in->parent);
     1102    if (in->parent != NULL) {
     1103        PS_WARN_PTR_NON_NULL(in->parent->concepts);
     1104    }
     1105    //
     1106    // Determine trimmed image from metadata.
     1107    //
     1108    psImage *trimmedImg = p_psDetermineTrimmedImage(in);
     1109 
     1110XXX: Create a pmUtils.c file and put this routine there.
     1111 *****************************************************************************/
     1112psImage *p_psDetermineTrimmedImage(pmReadout *in)
     1113{
     1114    if ((in->parent == NULL) || (in->parent->concepts == NULL)) {
     1115        psLogMsg(__func__, PS_LOG_WARN,
     1116                 "WARNING: could not determine CELL.TRIMSEC from parent cell Metadata (NULL).\n");
     1117        return(in->image);
     1118    }
     1119
     1120    psBool rc = false;
     1121    psImage *trimmedImg = NULL;
     1122    psRegion *trimRegion = psMetadataLookupPtr(&rc, in->parent->concepts,
     1123                           "CELL.TRIMSEC");
     1124    if (rc == false) {
     1125        psLogMsg(__func__, PS_LOG_WARN,
     1126                 "WARNING: could not determine CELL.TRIMSEC from parent cell Metadata.\n");
     1127        trimmedImg = in->image;
     1128    } else {
     1129        trimmedImg = psImageSubset(in->image, *trimRegion);
     1130    }
     1131
     1132    return(trimmedImg);
     1133}
     1134
     1135
     1136
     1137
     1138
     1139
     1140
     1141
     1142
     1143
     1144
     1145
     1146
     1147
     1148
     1149
     1150
     1151
     1152/******************************************************************************
     1153pmSubtractBias(....): see SDRS for complete specification.
     1154 
     1155XXX: Code and assert type support: U16, S32, F32.
     1156XXX: Add trace messages.
    2831157 *****************************************************************************/
    2841158pmReadout *pmSubtractBias(
     
    2901164    int nBin,
    2911165    const pmReadout *bias,
    292     const pmReadout *dark
    293 )
     1166    const pmReadout *dark)
    2941167{
    2951168    psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4,
    2961169            "---- pmSubtractBias() begin ----\n");
    297     PS_ASSERT_READOUT_NON_NULL(in, NULL);
    298     PS_ASSERT_READOUT_NON_EMPTY(in, NULL);
    299     PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, NULL);
    300 
    301     //
    302     // If the overscans != NULL, then check the type of each image.
    303     //
    304     if (overscans != NULL) {
    305         psListElem *tmpOverscan = (psListElem *) overscans->head;
    306         while (NULL != tmpOverscan) {
    307             psImage *myOverscanImage = (psImage *) tmpOverscan->data;
    308             PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL);
    309             tmpOverscan = tmpOverscan->next;
    310         }
    311     }
    312 
    313     if ((overscans == NULL) && (overScanAxis != PM_OVERSCAN_NONE)) {
    314         psError(PS_ERR_UNKNOWN,true, "(overscans == NULL) && (overScanAxis != PM_OVERSCAN_NONE).  Returning in image\n");
     1170    //
     1171    // Check input parameters, generate warnings and errors.
     1172    //
     1173    if (PM_OKAY != AssertCodeOverscan(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) {
    3151174        return(in);
    3161175    }
    317 
    318     // Check for an unallowable pmFit.
    319     if ((fit != PM_OVERSCAN_NONE) &&
    320             (fit != PM_OVERSCAN_ROWS) &&
    321             (fit != PM_OVERSCAN_COLUMNS) &&
    322             (fit != PM_OVERSCAN_ALL)) {
    323         psError(PS_ERR_UNKNOWN, true, "fit is unallowable (%d).  Returning in image.\n", fit);
    324         return(in);
    325     }
    326     // Check for an unallowable pmOverscanAxis.
    327     if ((overScanAxis != PM_OVERSCAN_NONE) &&
    328             (overScanAxis != PM_OVERSCAN_ROWS) &&
    329             (overScanAxis != PM_OVERSCAN_COLUMNS) &&
    330             (overScanAxis != PM_OVERSCAN_ALL)) {
    331         psError(PS_ERR_UNKNOWN, true, "overScanAxis is unallowable (%d).  Returning in image.\n", overScanAxis);
    332         return(in);
    333     }
    334     psS32 i;
    335     psS32 j;
    336     psS32 numBins = 0;
    337     static psVector *overscanVector = NULL;
    338     psVector *tmpRow = NULL;
    339     psVector *tmpCol = NULL;
    340     psVector *myBin = NULL;
    341     psVector *binVec = NULL;
    342     psListElem *tmpOverscan = NULL;
    343     double statValue;
    344     psImage *myOverscanImage = NULL;
    345     psPolynomial1D *myPoly = NULL;
    346     psSpline1D *mySpline = NULL;
    347     psS32 nBin;
    348 
    349     //
    350     //  Create a static stats data structure and determine the highest
    351     //  priority stats option.
    352     //
    353     static psStats *myStats = NULL;
    354     if (myStats == NULL) {
    355         myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
    356         p_psMemSetPersistent(myStats, true);
    357     }
    358     if (stat != NULL) {
    359         myStats->options = GenNewStatOptions(stat);
    360     }
    361 
    362 
    363     if (overScanAxis == PM_OVERSCAN_NONE) {
    364         if (fit != PM_FIT_NONE) {
    365             psLogMsg(__func__, PS_LOG_WARN,
    366                      "WARNING: pmSubtractBias.(): overScanAxis equals NONE, and fit does not equal NONE.  Proceeding to full fram subtraction.\n");
    367         }
    368 
    369         if (overscans != NULL) {
    370             psLogMsg(__func__, PS_LOG_WARN,
    371                      "WARNING: pmSubtractBias.(): overScanAxis equals NONE and overscans does not equal NULL.  Proceeding to full fram subtraction.\n");
    372         }
    373         return(SubtractFrame(in, bias));
    374     }
    375 
    376     if ((overScanAxis == PM_OVERSCAN_ALL) && (fit != PM_FIT_NONE)) {
    377         psLogMsg(__func__, PS_LOG_WARN,
    378                  "WARNING: pmSubtractBias.(): overScanAxis equals ALL, and fit does not equal NONE.  Proceeding with the rest of the module.\n");
    379     }
    380 
    381 
    382     //
    383     // We subtract each overscan region from the image data.
    384     // If we get here we know that overscans != NULL.
    385     //
    386 
    387     if (overScanAxis == PM_OVERSCAN_ALL) {
    388         tmpOverscan = (psListElem *) overscans->head;
    389         while (NULL != tmpOverscan) {
    390             myOverscanImage = (psImage *) tmpOverscan->data;
    391 
    392             PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL);
    393             psStats *rc = psImageStats(myStats, myOverscanImage, NULL, (psMaskType)0xffffffff);
    394             if (rc == NULL) {
    395                 psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation.  Returning in image.\n");
     1176    //
     1177    // Determine trimmed image from metadata.
     1178    //
     1179    psImage *trimmedImg = p_psDetermineTrimmedImage(in);
     1180
     1181    //
     1182    // Subtract overscan frames if necessary.
     1183    //
     1184    if (overscan == true) {
     1185        pmOverscanAxis overScanAxis = GetOverscanAxis(in);
     1186
     1187        //
     1188        //  Create a psStats data structure and determine the highest
     1189        //  priority stats option.
     1190        //
     1191        psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
     1192        if (stat != NULL) {
     1193            myStats->options = GenNewStatOptions(stat);
     1194        }
     1195
     1196        //
     1197        // Reduce overscan images to a single pixel, then subtract.
     1198        // This code is no longer required as of SDRS 12-09.
     1199        //
     1200        if (overScanAxis == PM_OVERSCAN_ALL) {
     1201            if (false == OverscanReducePixel(trimmedImg, in->bias, myStats)) {
    3961202                return(in);
    3971203            }
    398             if (false == p_psGetStatValue(myStats, &statValue)) {
    399                 psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
     1204            psFree(myStats);
     1205        } else {
     1206            //
     1207            // Reduce the overscan images to a single overscan vector.
     1208            //
     1209            psVector *overscanVector = OverscanReduce(in->image, overScanAxis,
     1210                                       in->bias, fitSpec,
     1211                                       fit, myStats);
     1212            if (overscanVector == NULL) {
     1213                psError(PS_ERR_UNKNOWN, false, "Could not reduce overscan images to a single overscan vector.  Returning in image\n");
     1214                psFree(myStats);
    4001215                return(in);
    4011216            }
    402             ImageSubtractScalar(in->image, statValue);
    403 
    404             tmpOverscan = tmpOverscan->next;
    405         }
    406         return(in);
    407     }
    408 
    409     // This check is redundant with above code.
    410     if (!((overScanAxis == PM_OVERSCAN_ROWS) || (overScanAxis == PM_OVERSCAN_COLUMNS))) {
    411         psError(PS_ERR_UNKNOWN, true, "overScanAxis is unallowable (%d).\nReturning in image.\n", overScanAxis);
    412         return(in);
    413     }
    414 
    415     tmpOverscan = (psListElem *) overscans->head;
    416     while (NULL != tmpOverscan) {
    417         //        PS_IMAGE_PRINT_F32_HIDEF(in->image);
    418         myOverscanImage = (psImage *) tmpOverscan->data;
    419 
    420         if (overScanAxis == PM_OVERSCAN_ROWS) {
    421             if (myOverscanImage->numCols != (in->image)->numCols) {
    422                 psLogMsg(__func__, PS_LOG_WARN,
    423                          "WARNING: pmSubtractBias.(): overscan image has %d columns, input image has %d columns\n",
    424                          myOverscanImage->numCols, in->image->numCols);
    425             }
    426 
    427             // We create a row vector and subtract this vector from image.
    428             // XXX: Is there a better way to extract a psVector from a psImage without
    429             // having to copy every element in that vector?
    430             overscanVector = psVectorAlloc(myOverscanImage->numCols, PS_TYPE_F32);
    431             for (i=0;i<overscanVector->n;i++) {
    432                 overscanVector->data.F32[i] = 0.0;
    433             }
    434             tmpRow = psVectorAlloc(myOverscanImage->numRows, PS_TYPE_F32);
    435 
    436             // For each column of the input image, loop through every row,
    437             // collect the pixel in that row, then performed the specified
    438             // statistical op on those pixels.  Store this in overscanVector.
    439             for (i=0;i<myOverscanImage->numCols;i++) {
    440                 for (j=0;j<myOverscanImage->numRows;j++) {
    441                     tmpRow->data.F32[j] = myOverscanImage->data.F32[j][i];
    442                 }
    443                 psStats *rc = psVectorStats(myStats, tmpRow, NULL, NULL, 0);
    444                 if (rc == NULL) {
    445                     psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
     1217
     1218            //
     1219            // Rebin the overscan vector if necessary.
     1220            //
     1221            psS32 newBin = RebinOverscanVector(overscanVector, nBin, myStats);
     1222            if (newBin < 0) {
     1223                psError(PS_ERR_UNKNOWN, false, "Could rebin the overscan vector.  Returning in image\n");
     1224                psFree(myStats);
     1225                return(in);
     1226            }
     1227
     1228            //
     1229            // If necessary, fit a psPolynomial or psSpline to the overscan vector.
     1230            // Then, unbin the overscan vector to appropriate length for the in image.
     1231            //
     1232            if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) {
     1233                overscanVector = FitOverscanVectorAndUnbin(trimmedImg, overscanVector, overScanAxis, fitSpec, fit, newBin);
     1234                if (overscanVector == NULL) {
     1235                    psError(PS_ERR_UNKNOWN, false, "Could not fit the polynomial or spline to the overscan vector.  Returning in image\n");
     1236                    psFree(myStats);
    4461237                    return(in);
    4471238                }
    448                 if (false ==  p_psGetStatValue(rc, &statValue)) {
    449                     psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
    450                     return(in);
     1239            } else {
     1240                overscanVector = UnbinOverscanVector(trimmedImg, overscanVector, overScanAxis, newBin);
     1241            }
     1242
     1243            //
     1244            // Subtract the overscan vector from the input image.
     1245            //
     1246            SubtractVectorFromImage(trimmedImg, overscanVector, overScanAxis);
     1247            psFree(myStats);
     1248            psFree(overscanVector);
     1249        }
     1250    }
     1251
     1252    //
     1253    // Perform bias subtraction if necessary.
     1254    //
     1255    if (bias != NULL) {
     1256        if (PM_OKAY == AssertCodeBias(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) {
     1257            SubtractFrame(in, bias);
     1258        }
     1259    }
     1260
     1261    //
     1262    // Perform dark subtraction if necessary.
     1263    //
     1264    if (dark != NULL) {
     1265        if (PM_OKAY == AssertCodeDark(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) {
     1266            psBool rc;
     1267            psF32 scale = 0.0;
     1268            if (in->parent != NULL) {
     1269                scale = psMetadataLookupS32(&rc, in->parent->concepts, "CELL.DARKTIME");
     1270                if (rc == false) {
     1271                    psLogMsg(__func__, PS_LOG_WARN,
     1272                             "WARNING: pmSubtractBias.(): could not determine CELL.FARKTIME from in->parent metadata.\n");
    4511273                }
    452                 overscanVector->data.F32[i] = statValue;
    453             }
    454             psFree(tmpRow);
    455 
    456             // Scale the overscan vector to the size of the input image.
    457             if (overscanVector->n != in->image->numCols) {
    458                 if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) {
    459                     psVector *newVec = ScaleOverscanVector(overscanVector,
    460                                                            in->image->numCols,
    461                                                            fitSpec, fit);
    462                     if (newVec == NULL) {
    463                         psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector.  Returning in image.\n");
    464                         return(in);
    465                     }
    466                     psFree(overscanVector);
    467                     overscanVector = newVec;
    468                 } else {
    469                     psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vector.  Set fit to PM_FIT_SPLINE or PM_FIT_POLYNOMIAL.  Returning in image.\n");
    470                     psFree(overscanVector);
    471                     return(in);
    472                 }
    473             }
    474         }
    475 
    476         if (overScanAxis == PM_OVERSCAN_COLUMNS) {
    477             if (myOverscanImage->numRows != (in->image)->numRows) {
    478                 psLogMsg(__func__, PS_LOG_WARN,
    479                          "WARNING: pmSubtractBias.(): overscan image has %d rows, input image has %d rows\n",
    480                          myOverscanImage->numRows, in->image->numRows);
    481             }
    482 
    483             // We create a column vector and subtract this vector from image.
    484             overscanVector = psVectorAlloc(myOverscanImage->numRows, PS_TYPE_F32);
    485             for (i=0;i<overscanVector->n;i++) {
    486                 overscanVector->data.F32[i] = 0.0;
    487             }
    488             tmpCol = psVectorAlloc(myOverscanImage->numCols, PS_TYPE_F32);
    489 
    490             // For each row of the input image, loop through every column,
    491             // collect the pixel in that row, then performed the specified
    492             // statistical op on those pixels.  Store this in overscanVector.
    493             for (i=0;i<myOverscanImage->numRows;i++) {
    494                 for (j=0;j<myOverscanImage->numCols;j++) {
    495                     tmpCol->data.F32[j] = myOverscanImage->data.F32[i][j];
    496                 }
    497                 psStats *rc = psVectorStats(myStats, tmpCol, NULL, NULL, 0);
    498                 if (rc == NULL) {
    499                     psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
    500                     return(in);
    501                 }
    502                 if (false ==  p_psGetStatValue(rc, &statValue)) {
    503                     psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
    504                     return(in);
    505                 }
    506                 overscanVector->data.F32[i] = statValue;
    507             }
    508             psFree(tmpCol);
    509 
    510             // Scale the overscan vector to the size of the input image.
    511             if (overscanVector->n != in->image->numRows) {
    512                 if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) {
    513                     psVector *newVec = ScaleOverscanVector(overscanVector,
    514                                                            in->image->numRows,
    515                                                            fitSpec, fit);
    516                     if (newVec == NULL) {
    517                         psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector.  Returning in image.\n");
    518                         return(in);
    519                     }
    520                     psFree(overscanVector);
    521                     overscanVector = newVec;
    522                 } else {
    523                     psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vector.  Set fit to PM_FIT_SPLINE or PM_FIT_POLYNOMIAL.  Returning in image.\n");
    524                     psFree(overscanVector);
    525                     return(in);
    526                 }
    527             }
    528         }
    529 
    530         //
    531         // Re-bin the overscan vector (change its length).
    532         //
    533         // Only if nBinOrig > 1.
    534         if ((nBinOrig > 1) && (nBinOrig < overscanVector->n)) {
    535             numBins = 1+((overscanVector->n)/nBinOrig);
    536             myBin = psVectorAlloc(numBins, PS_TYPE_F32);
    537             binVec = psVectorAlloc(nBinOrig, PS_TYPE_F32);
    538 
    539             for (i=0;i<numBins;i++) {
    540                 for(j=0;j<nBinOrig;j++) {
    541                     if (overscanVector->n > ((i*nBinOrig)+j)) {
    542                         binVec->data.F32[j] = overscanVector->data.F32[(i*nBinOrig)+j];
    543                     } else {
    544                         // XXX: we get here if nBinOrig does not evenly divide
    545                         // the overscanVector vector.  This is the last bin.  Should
    546                         // we change the binVec->n to acknowledge that?
    547                         binVec->n = j;
    548                     }
    549                 }
    550                 psStats *rc = psVectorStats(myStats, binVec, NULL, NULL, 0);
    551                 if (rc == NULL) {
    552                     psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
    553                     return(in);
    554                 }
    555                 if (false ==  p_psGetStatValue(rc, &statValue)) {
    556                     psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
    557                     return(in);
    558                 }
    559                 myBin->data.F32[i] = statValue;
    560             }
    561 
    562             // Change the effective size of overscanVector.
    563             overscanVector->n = numBins;
    564             for (i=0;i<numBins;i++) {
    565                 overscanVector->data.F32[i] = myBin->data.F32[i];
    566             }
    567             psFree(binVec);
    568             psFree(myBin);
    569             nBin = nBinOrig;
    570         } else {
    571             nBin = 1;
    572         }
    573 
    574         // At this point the number of data points in overscanVector should be
    575         // equal to the number of rows/columns (whatever is appropriate) in the
    576         // image divided by numBins.
    577         //
    578 
    579 
    580         //
    581         // This doesn't seem right.  The only way to do a spline fit is if,
    582         // by SDRS requirements, fitSpec is not-NULL>  But in order for it
    583         // to be non-NULL, someone must have called psSpline1DAlloc() with
    584         // the min, max, and number of splines.
    585         //
    586         if (!((fitSpec == NULL) || (fit == PM_FIT_NONE))) {
    587             //
    588             // Fit a polynomial or spline to the overscan vector.
    589             //
    590             if (fit == PM_FIT_POLYNOMIAL) {
    591                 myPoly = (psPolynomial1D *) fitSpec;
    592                 myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL);
    593                 if (myPoly == NULL) {
    594                     psError(PS_ERR_UNKNOWN, false, "(3) Could not fit a polynomial to overscan vector.  Returning in image.\n");
    595                     psFree(overscanVector);
    596                     return(in);
    597                 }
    598             } else if (fit == PM_FIT_SPLINE) {
    599                 // XXX: This makes no sense
    600                 // XXX: must free mySpline?
    601                 mySpline = (psSpline1D *) fitSpec;
    602                 mySpline = psVectorFitSpline1D(NULL, overscanVector);
    603                 if (mySpline == NULL) {
    604                     psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to overscan vector.  Returning in image.\n");
    605                     psFree(overscanVector);
    606                     return(in);
    607                 }
    608             }
    609 
    610             //
    611             // Subtract fitted overscan vector row-wise from the image.
    612             //
    613             if (overScanAxis == PM_OVERSCAN_ROWS) {
    614                 for (i=0;i<(in->image)->numCols;i++) {
    615                     psF32 tmpF32 = 0.0;
    616                     if (fit == PM_FIT_POLYNOMIAL) {
    617                         tmpF32 = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin));
    618                     } else if (fit == PM_FIT_SPLINE) {
    619                         tmpF32 = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin));
    620                     }
    621                     for (j=0;j<(in->image)->numRows;j++) {
    622                         (in->image)->data.F32[j][i]-= tmpF32;
    623                     }
    624                 }
    625             }
    626 
    627             //
    628             // Subtract fitted overscan vector column-wise from the image.
    629             //
    630             if (overScanAxis == PM_OVERSCAN_COLUMNS) {
    631                 for (i=0;i<(in->image)->numRows;i++) {
    632                     psF32 tmpF32 = 0.0;
    633                     if (fit == PM_FIT_POLYNOMIAL) {
    634                         tmpF32 = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin));
    635                     } else if (fit == PM_FIT_SPLINE) {
    636                         tmpF32 = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin));
    637                     }
    638 
    639                     for (j=0;j<(in->image)->numCols;j++) {
    640                         (in->image)->data.F32[i][j]-= tmpF32;
    641                     }
    642                 }
    643             }
    644         } else {
    645             //
    646             // If we get here, then no polynomials were fit to the overscan
    647             // vector.  We simply subtract it, taking into account binning,
    648             // from the image.
    649             //
    650 
    651             //
    652             // Subtract overscan vector row-wise from the image.
    653             //
    654             if (overScanAxis == PM_OVERSCAN_ROWS) {
    655                 for (i=0;i<(in->image)->numCols;i++) {
    656                     for (j=0;j<(in->image)->numRows;j++) {
    657                         (in->image)->data.F32[j][i]-= overscanVector->data.F32[i/nBin];
    658                     }
    659                 }
    660             }
    661 
    662             //
    663             // Subtract overscan vector column-wise from the image.
    664             //
    665             if (overScanAxis == PM_OVERSCAN_COLUMNS) {
    666                 for (i=0;i<(in->image)->numRows;i++) {
    667                     for (j=0;j<(in->image)->numCols;j++) {
    668                         (in->image)->data.F32[i][j]-= overscanVector->data.F32[i/nBin];
    669                     }
    670                 }
    671             }
    672         }
    673 
    674         psFree(overscanVector);
    675 
    676         tmpOverscan = tmpOverscan->next;
    677     }
    678 
     1274            }
     1275            SubtractDarkFrame(in, dark, scale);
     1276        }
     1277    }
     1278
     1279    //
     1280    // All done.
     1281    //
    6791282    psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4,
    6801283            "---- pmSubtractBias() exit ----\n");
    681 
    682     if (bias != NULL) {
    683         return(SubtractFrame(in, bias));
    684     }
    6851284    return(in);
    6861285}
  • trunk/psModules/src/imsubtract/pmSubtractBias.h

    r5435 r5516  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-10-20 23:06:24 $
     8 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2005-11-15 20:09:03 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    4747    const pmReadout *bias,              ///< A possibly NULL bias pmReadout which is to be subtracted
    4848    const pmReadout *dark               ///< A possibly NULL bias pmReadout which is to be subtracted
    49 )
     49);
    5050// OLD: remove this
    5151//    const psList *overscans,      ///< A psList of overscan images
    5252//    pmOverscanAxis overScanAxis,  ///< Defines how overscans are applied
    5353
     54/******************************************************************************
     55DetermineTrimmedImageg() This private routine determines the region of the
     56input pmReadout which will be operated on by the various detrend modules.  It
     57does a metadata fetch on "CELL.TRIMSEC" for the parent cell of the pmReadout.
     58 
     59XXX: Consider making a pmUtils.c file and put this routine there.
     60 *****************************************************************************/
     61psImage *p_psDetermineTrimmedImage(
     62    pmReadout *in
     63);
     64
    5465#endif
  • trunk/psModules/src/imsubtract/pmSubtractSky.c

    r5294 r5516  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-10-12 21:02:04 $
     8 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2005-11-15 20:09:03 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    1818#include<math.h>
    1919#include "pslib.h"
    20 #include "psConstants.h"
    2120#include "pmSubtractSky.h"
     21
     22// XXX: Get rid of the.  Create pmUtils.h
     23psImage *p_psDetermineTrimmedImage(
     24    pmReadout *in
     25);
    2226
    2327/******************************************************************************
     
    152156 
    153157XXX: Use your brain and figure out the analytical expression.
     158 
     159XXX: Why isn't it simply (xOrder+1) * (yOrder+1)?
    154160 *****************************************************************************/
    155161static psS32 CalculatePolyTerms(psS32 xOrder, psS32 yOrder)
     
    170176        }
    171177    }
    172 
    173178    psTrace("SubtractSky.CalculatePolyTerms", 4,
    174179            "Exiting CalculatePolyTerms(%d, %d) -> %d\n", xOrder, yOrder, localPolyTerms);
    175180    return(localPolyTerms);
     181
     182    //    return((xOrder+1) * (yOrder+1));
    176183}
    177184
     
    283290XXX: Different trace message facilities in use here.
    284291 *****************************************************************************/
    285 static psPolynomial2D *ImageFitPolynomial(psPolynomial2D *myPoly,
    286         psImage *dataImage,
    287         psImage *maskImage)
     292static psPolynomial2D *ImageFitPolynomial(
     293    psPolynomial2D *myPoly,
     294    psImage *dataImage,
     295    psImage *maskImage)
    288296{
    289297    psTrace("SubtractSky.ImageFitPolynomial", 4,
     
    471479    PS_ASSERT_READOUT_NON_EMPTY(in, NULL);
    472480    PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, NULL);
     481    PS_WARN_PTR_NON_NULL(in->parent);
     482    if (in->parent != NULL) {
     483        PS_WARN_PTR_NON_NULL(in->parent->concepts);
     484    }
    473485    psTrace(".psModule.pmSubtractSky", 4,
    474486            "---- pmSubtractSky() begin ----\n");
     
    492504        return(in);
    493505    }
    494     psImage *origImage = in->image;
     506
     507    //
     508    // Determine trimmed image from metadata.
     509    //
     510
     511    psImage *trimmedImg = p_psDetermineTrimmedImage(in);
    495512    psImage *binnedImage = NULL;
    496513    psPolynomial2D *myPoly = NULL;
     
    539556        // No binning is required here.  Simply create a copy of the image
    540557        // and a mask.
    541         binnedImage = psImageCopy(binnedImage, origImage, PS_TYPE_F32);
     558        binnedImage = psImageCopy(binnedImage, trimmedImg, PS_TYPE_F32);
    542559        if (binnedImage == NULL) {
    543560            psError(PS_ERR_UNKNOWN, false, "psImageCopy() returned NULL.  Returning in image.\n");
     
    559576        }
    560577    } else {
    561         binnedImage = psImageRebin(NULL, origImage, in->mask, 0, binFactor, stats);
     578        binnedImage = psImageRebin(NULL, trimmedImg, in->mask, 0, binFactor, stats);
    562579        if (binnedImage == NULL) {
    563580            psError(PS_ERR_UNKNOWN, false, "psImageRebin() returned NULL.  Returning in image.\n");
     
    677694    if (binFactor <= 1) {
    678695        // The binned image is the same size as the original image.
    679         for (psS32 row = 0; row < origImage->numRows ; row++) {
    680             for (psS32 col = 0; col < origImage->numCols ; col++) {
    681                 origImage->data.F32[row][col]-= binnedImage->data.F32[row][col];
     696        for (psS32 row = 0; row < trimmedImg->numRows ; row++) {
     697            for (psS32 col = 0; col < trimmedImg->numCols ; col++) {
     698                trimmedImg->data.F32[row][col]-= binnedImage->data.F32[row][col];
    682699            }
    683700        }
    684701    } else {
    685         for (psS32 row = 0; row < origImage->numRows ; row++) {
    686             for (psS32 col = 0; col < origImage->numCols ; col++) {
     702        for (psS32 row = 0; row < trimmedImg->numRows ; row++) {
     703            for (psS32 col = 0; col < trimmedImg->numCols ; col++) {
    687704                // We calculate the F32 value of the pixel coordinates in the
    688705                // binned image and then use a pixel interpolation routine to
     
    700717                                     binnedImage, binColF64, binRowF64,
    701718                                     NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
    702                 origImage->data.F32[row][col]-= binPixel;
     719                trimmedImg->data.F32[row][col]-= binPixel;
    703720
    704721                psTrace(".psModule.pmSubtractSky", 8,
  • trunk/psModules/test/detrend/tst_pmFlatField.c

    r5435 r5516  
    2121 *  @author Ross Harman, MHPCC
    2222 *
    23  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    24  *  @date $Date: 2005-10-20 23:06:24 $
     23 *  XXX: I added the CELL.TRIMSEC region code but there are not tests for it.
     24 *
     25 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     26 *  @date $Date: 2005-11-15 20:09:03 $
    2527 *
    2628 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  • trunk/psModules/test/detrend/tst_pmMaskBadPixels.c

    r5169 r5516  
    1717 *  @author Ross Harman, MHPCC
    1818 *
    19  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    20  *  @date $Date: 2005-09-28 20:42:52 $
     19 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
     20 *  @date $Date: 2005-11-15 20:09:03 $
    2121 *
    2222 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    3232printf(STRING);                                                                                              \
    3333printf("\n");                                                                                                \
    34 for(int i=IMAGE->numRows-1; i>-1; i--) {                                                                     \
    35     for(int j=0; j<IMAGE->numCols; j++) {                                                                    \
    36         if(PS_IS_PSELEMTYPE_COMPLEX(IMAGE->type.type)) {                                                     \
    37             printf("%f+%fi ", creal(IMAGE->data.TYPE[i][j]), cimag(IMAGE->data.TYPE[i][j]));                 \
    38         } else if(PS_IS_PSELEMTYPE_INT(IMAGE->type.type)) {                                                  \
    39             printf("%d", (int)IMAGE->data.TYPE[i][j]);                                                       \
     34for(int i=(IMAGE)->numRows-1; i>-1; i--) {                                                                     \
     35    for(int j=0; j<(IMAGE)->numCols; j++) {                                                                    \
     36        if(PS_IS_PSELEMTYPE_COMPLEX((IMAGE)->type.type)) {                                                     \
     37            printf("%f+%fi ", creal((IMAGE)->data.TYPE[i][j]), cimag((IMAGE)->data.TYPE[i][j]));                 \
     38        } else if(PS_IS_PSELEMTYPE_INT((IMAGE)->type.type)) {                                                  \
     39            printf("%d", (int)(IMAGE)->data.TYPE[i][j]);                                                       \
    4040        } else {                                                                                             \
    41             printf("%f", (double)IMAGE->data.TYPE[i][j]);                                                    \
     41            printf("%f", (double)(IMAGE)->data.TYPE[i][j]);                                                    \
    4242        }                                                                                                    \
    4343    }                                                                                                        \
     
    4848
    4949#define CREATE_AND_SET_IMAGE(NAME,TYPE,VALUE,NROWS,NCOLS)                                                    \
    50 psImage *NAME = (psImage*)psImageAlloc(NCOLS,NROWS,PS_TYPE_##TYPE);                                          \
    51 for(int i=0; i<NAME->numRows; i++) {                                                                         \
    52     for(int j=0; j<NAME->numCols; j++) {                                                                     \
    53         NAME->data.TYPE[i][j] = VALUE;                                                                       \
     50(NAME) = (psImage*)psImageAlloc(NCOLS,NROWS,PS_TYPE_##TYPE);                                             \
     51for(int i=0; i<(NAME)->numRows; i++) {                                                                         \
     52    for(int j=0; j<(NAME)->numCols; j++) {                                                                     \
     53        (NAME)->data.TYPE[i][j] = VALUE;                                                                       \
    5454    }                                                                                                        \
    5555}
    56 
    5756
    5857static int testMaskBadPixels1(void);
     
    8079                              {testMaskBadPixels9, 885, "pmMaskBadPixels - Attempt to use offset greater than input image", 0, false},
    8180                              {testMaskBadPixels10, 885, "pmMaskBadPixels - Attempt to use complex input image", 0, false},
    82                               {testMaskBadPixels11, 885, "pmMaskBadPixels - Attempt to use non-mask type mask image", 0, false},
     81                              {testMaskBadPixels11, 885, "pmMaskBadPixels - Attempt to use non-mask type mask image", 0, false        },
    8382                              {NULL}
    8483                          };
    8584
     85/*
     86    #define PS_TYPE_MASK PS_TYPE_U8
     87    #define PS_TYPE_MASK_DATA U8
     88    #define PS_TYPE_MASK_NAME "psU8"
     89*/
    8690
    8791int main(int argc, char* argv[])
     
    9296
    9397
     98#define NUM_ROWS 50
     99#define NUM_COLS 50
     100#define DEFAULT_IMAGE_VAL 0.0
     101#define DEFAULT_MASK_VAL 0
     102#define MASK_VAL 1
     103#define SAT_VAL  100.0
     104#define GROW_VAL 1
     105#define GROW_RAD 10
    94106int testMaskBadPixels1( void )
    95107{
     108    //
    96109    // Test A - Create mask based on maskVal argument
    97     CREATE_AND_SET_IMAGE(inImage1,F64,0,50,50);
    98     CREATE_AND_SET_IMAGE(mask1,U8,0,50,50)
    99     //    pmReadout *inReadout = pmReadoutAlloc(0, 0, inImage1);
    100     pmReadout *inReadout = pmReadoutAlloc(NULL);
    101     inReadout->row0 = 0;
    102     inReadout->col0 = 0;
    103     inReadout->image = inImage1;
    104     mask1->data.PS_TYPE_MASK_DATA[24][24]=1;
    105     PRINT_MATRIX(mask1, U8, "Data mask:");
    106     pmMaskBadPixels(inReadout, mask1, 1, 100.0, 0, 0);
     110    //
     111
     112    pmReadout *inReadout = pmReadoutAlloc(NULL);
     113    CREATE_AND_SET_IMAGE(inReadout->image, F32, DEFAULT_IMAGE_VAL, NUM_ROWS, NUM_COLS);
     114    inReadout->image->row0 = 0;
     115    inReadout->image->col0 = 0;
     116    inReadout->row0 = 0;
     117    inReadout->col0 = 0;
     118
     119    pmReadout *maskReadout = pmReadoutAlloc(NULL);
     120    CREATE_AND_SET_IMAGE(maskReadout->image, U8, DEFAULT_MASK_VAL, NUM_ROWS, NUM_COLS)
     121    maskReadout->image->data.PS_TYPE_MASK_DATA[NUM_ROWS/2][NUM_COLS/2]=1;
     122    maskReadout->image->row0 = 0;
     123    maskReadout->image->col0 = 0;
     124
     125    PRINT_MATRIX(maskReadout->image, U8, "Data mask:");
     126
     127    pmMaskBadPixels(inReadout, maskReadout, MASK_VAL, SAT_VAL, GROW_VAL, GROW_RAD);
     128    PRINT_MATRIX(inReadout->mask, PS_TYPE_MASK_DATA, "Resulting mask:");
     129
     130    psFree(inReadout);
     131    psFree(maskReadout);
     132
     133    return 0;
     134}
     135
     136int testMaskBadPixels2( void )
     137{
     138    //
     139    // Test B - Create mask based on saturation argument
     140    //
     141
     142    pmReadout *inReadout = pmReadoutAlloc(NULL);
     143    CREATE_AND_SET_IMAGE(inReadout->image, F32, DEFAULT_IMAGE_VAL, NUM_ROWS, NUM_COLS);
     144    inReadout->image->row0 = 0;
     145    inReadout->image->col0 = 0;
     146    inReadout->row0 = 0;
     147    inReadout->col0 = 0;
     148    inReadout->image->data.F32[NUM_ROWS/2][NUM_COLS/2] = 150.0;
     149
     150    pmReadout *maskReadout = pmReadoutAlloc(NULL);
     151    CREATE_AND_SET_IMAGE(maskReadout->image, U8, DEFAULT_MASK_VAL, NUM_ROWS, NUM_COLS)
     152
     153    //PS_IMAGE_PRINT_F32(inReadout->image);
     154    PRINT_MATRIX(maskReadout->image, U8, "Data mask:");
     155    pmMaskBadPixels(inReadout, maskReadout, MASK_VAL, SAT_VAL, GROW_VAL, GROW_RAD);
    107156    PRINT_MATRIX(inReadout->mask, U8, "Resulting mask:");
    108     psFree(mask1);
    109     psFree(inReadout);
    110 
    111     return 0;
    112 }
    113 
    114 int testMaskBadPixels2( void )
    115 {
    116     // Test B - Create mask based on saturation argument
    117     CREATE_AND_SET_IMAGE(inImage2,F64,150.0,50,50);
    118     CREATE_AND_SET_IMAGE(mask2,U8,0,50,50)
    119     //    pmReadout *inReadout2 = pmReadoutAlloc(0, 0, inImage2);
    120     pmReadout *inReadout2 = pmReadoutAlloc(NULL);
    121     inReadout2->row0 = 0;
    122     inReadout2->col0 = 0;
    123     inReadout2->image = inImage2;
    124     PRINT_MATRIX(mask2, U8, "Data mask:");
    125     pmMaskBadPixels(inReadout2, mask2, 0, 100.0, 0, 0);
    126     PRINT_MATRIX(inReadout2->mask, U8, "Resulting mask:");
    127     psFree(mask2);
    128     psFree(inReadout2);
     157
     158    psFree(inReadout);
     159    psFree(maskReadout);
    129160
    130161    return 0;
     
    133164int testMaskBadPixels3( void )
    134165{
     166    //
    135167    // Test C - Create mask based on growVal and grow arguments
    136     CREATE_AND_SET_IMAGE(inImage3,F64,50.0,50,50);
    137     CREATE_AND_SET_IMAGE(mask3,U8,0,50,50)
    138     //    pmReadout *inReadout3 = pmReadoutAlloc(0, 0, inImage3);
    139     pmReadout *inReadout3 = pmReadoutAlloc(NULL);
    140     inReadout3->row0 = 0;
    141     inReadout3->col0 = 0;
    142     inReadout3->image = inImage3;
    143     mask3->data.PS_TYPE_MASK_DATA[24][24]=1;
    144     mask3->data.PS_TYPE_MASK_DATA[4][3]=1;
    145     mask3->data.PS_TYPE_MASK_DATA[4][46]=1;
    146     PRINT_MATRIX(mask3, U8, "Data mask:");
    147     pmMaskBadPixels(inReadout3, mask3, 0, 100.0, 1, 10);
    148     PRINT_MATRIX(inReadout3->mask, U8, "Resulting mask:");
    149     psFree(mask3);
    150     psFree(inReadout3);
     168    //
     169
     170    pmReadout *inReadout = pmReadoutAlloc(NULL);
     171    CREATE_AND_SET_IMAGE(inReadout->image, F32, DEFAULT_IMAGE_VAL, NUM_ROWS, NUM_COLS);
     172    inReadout->image->row0 = 0;
     173    inReadout->image->col0 = 0;
     174    inReadout->row0 = 0;
     175    inReadout->col0 = 0;
     176    inReadout->image->data.F32[NUM_ROWS/2][NUM_COLS/2]=GROW_VAL;
     177    inReadout->image->data.F32[NUM_ROWS/4][NUM_COLS/4]=GROW_VAL;
     178    inReadout->image->data.F32[NUM_ROWS/4][NUM_COLS-(NUM_COLS/4)]=GROW_VAL;
     179
     180    pmReadout *maskReadout = pmReadoutAlloc(NULL);
     181    CREATE_AND_SET_IMAGE(maskReadout->image, U8, DEFAULT_MASK_VAL, NUM_ROWS, NUM_COLS);
     182
     183    PRINT_MATRIX(maskReadout->image, U8, "Data mask:");
     184    //PS_IMAGE_PRINT_F32(inReadout->image);
     185    pmMaskBadPixels(inReadout, maskReadout, MASK_VAL, SAT_VAL, GROW_VAL, GROW_RAD);
     186    PRINT_MATRIX(inReadout->mask, U8, "Resulting mask:");
     187
     188    psFree(inReadout);
     189    psFree(maskReadout);
    151190
    152191    return 0;
     
    155194int testMaskBadPixels4( void )
    156195{
     196    //
    157197    // Test D - Auto Create mask based on maskVal argument
    158     CREATE_AND_SET_IMAGE(inImage4,F64,50.0,50,50);
    159     CREATE_AND_SET_IMAGE(mask4,U8,0,50,50)
    160     CREATE_AND_SET_IMAGE(mask4i,U8,0,50,50)
    161     //    pmReadout *inReadout4 = pmReadoutAlloc(0, 0, inImage4);
    162     pmReadout *inReadout4 = pmReadoutAlloc(NULL);
    163     inReadout4->row0 = 0;
    164     inReadout4->col0 = 0;
    165     inReadout4->image = inImage4;
    166     inReadout4->mask = mask4i;
    167     mask4->data.PS_TYPE_MASK_DATA[24][24]=1;
    168     PRINT_MATRIX(mask4, U8, "Data mask:");
    169     pmMaskBadPixels(inReadout4, mask4, 0, 100.0, 1, 10);
    170     PRINT_MATRIX(inReadout4->mask, U8, "Resulting mask:");
    171     psFree(mask4);
    172     psFree(inReadout4);
     198    //
     199
     200    pmReadout *inReadout = pmReadoutAlloc(NULL);
     201    CREATE_AND_SET_IMAGE(inReadout->image, F32, DEFAULT_IMAGE_VAL, NUM_ROWS, NUM_COLS);
     202    inReadout->image->row0 = 0;
     203    inReadout->image->col0 = 0;
     204    inReadout->row0 = 0;
     205    inReadout->col0 = 0;
     206
     207    pmReadout *maskReadout = pmReadoutAlloc(NULL);
     208    CREATE_AND_SET_IMAGE(maskReadout->image, U8, DEFAULT_MASK_VAL, NUM_ROWS, NUM_COLS)
     209    maskReadout->image->data.PS_TYPE_MASK_DATA[NUM_ROWS/2][NUM_COLS/2]=1;
     210
     211    PRINT_MATRIX(maskReadout->image, U8, "Data mask:");
     212    pmMaskBadPixels(inReadout, maskReadout, MASK_VAL, SAT_VAL, GROW_VAL, GROW_RAD);
     213    PRINT_MATRIX(inReadout->mask, U8, "Resulting mask:");
     214
     215    psFree(inReadout);
     216    psFree(maskReadout);
    173217
    174218    return 0;
     
    177221int testMaskBadPixels5( void )
    178222{
     223    //
    179224    // Test E - Attempt to use null mask
    180     CREATE_AND_SET_IMAGE(inImage5,F64,50.0,50,50);
    181     //    pmReadout *inReadout5 = pmReadoutAlloc(0, 0, inImage5);
    182     pmReadout *inReadout5 = pmReadoutAlloc(NULL);
    183     inReadout5->row0 = 0;
    184     inReadout5->col0 = 0;
    185     inReadout5->image = inImage5;
    186     pmMaskBadPixels(inReadout5, NULL, 0, 100.0, 1, 10);
    187     psFree(inReadout5);
     225    //
     226
     227    pmReadout *inReadout = pmReadoutAlloc(NULL);
     228    CREATE_AND_SET_IMAGE(inReadout->image, F32, DEFAULT_IMAGE_VAL, NUM_ROWS, NUM_COLS);
     229    inReadout->image->row0 = 0;
     230    inReadout->image->col0 = 0;
     231    inReadout->row0 = 0;
     232    inReadout->col0 = 0;
     233
     234    pmMaskBadPixels(inReadout, NULL, MASK_VAL, SAT_VAL, GROW_VAL, GROW_RAD);
     235    psFree(inReadout);
    188236
    189237    return 0;
     
    192240int testMaskBadPixels6( void )
    193241{
     242    //
    194243    // Test F - Attempt tp use null input image
    195     CREATE_AND_SET_IMAGE(mask6,U8,0,50,50)
    196     //    pmReadout *inReadout6 = pmReadoutAlloc(0, 0, NULL);
    197     pmReadout *inReadout6 = pmReadoutAlloc(NULL);
    198     inReadout6->row0 = 0;
    199     inReadout6->col0 = 0;
    200     inReadout6->mask = mask6;
    201     pmMaskBadPixels(inReadout6, mask6, 0, 100.0, 1, 10);
    202     psFree(inReadout6);
     244    //
     245
     246    pmReadout *inReadout = pmReadoutAlloc(NULL);
     247    inReadout->row0 = 0;
     248    inReadout->col0 = 0;
     249
     250    pmReadout *maskReadout = pmReadoutAlloc(NULL);
     251    CREATE_AND_SET_IMAGE(maskReadout->image, U8, DEFAULT_MASK_VAL, NUM_ROWS, NUM_COLS)
     252
     253    pmMaskBadPixels(inReadout, maskReadout, MASK_VAL, SAT_VAL, GROW_VAL, GROW_RAD);
     254
     255    psFree(inReadout);
     256    psFree(maskReadout);
    203257
    204258    return 0;
     
    207261int testMaskBadPixels7( void )
    208262{
     263    //
    209264    // Test G - Attempt to use input image bigger than mask
    210     CREATE_AND_SET_IMAGE(inImage7,F64,0.0,60,60);
    211     CREATE_AND_SET_IMAGE(mask7,U8,0,50,50)
    212     CREATE_AND_SET_IMAGE(mask7i,U8,0,50,50)
    213     //    pmReadout *inReadout7 = pmReadoutAlloc(0, 0, inImage7);
    214     pmReadout *inReadout7 = pmReadoutAlloc(NULL);
    215     inReadout7->row0 = 0;
    216     inReadout7->col0 = 0;
    217     inReadout7->image = inImage7;
    218     inReadout7->mask = mask7i;
    219     pmMaskBadPixels(inReadout7, mask7, 0, 100.0, 1, 10);
    220     psFree(mask7);
    221     psFree(inReadout7);
     265    //
     266
     267    pmReadout *inReadout = pmReadoutAlloc(NULL);
     268    CREATE_AND_SET_IMAGE(inReadout->image, F32, DEFAULT_IMAGE_VAL, NUM_ROWS+10, NUM_COLS+10);
     269    inReadout->image->row0 = 0;
     270    inReadout->image->col0 = 0;
     271    inReadout->row0 = 0;
     272    inReadout->col0 = 0;
     273
     274    pmReadout *maskReadout = pmReadoutAlloc(NULL);
     275    CREATE_AND_SET_IMAGE(maskReadout->image, U8, DEFAULT_MASK_VAL, NUM_ROWS, NUM_COLS)
     276
     277    pmMaskBadPixels(inReadout, maskReadout, MASK_VAL, SAT_VAL, GROW_VAL, GROW_RAD);
     278
     279    psFree(inReadout);
     280    psFree(maskReadout);
    222281
    223282    return 0;
     
    226285int testMaskBadPixels8( void )
    227286{
    228     // Test H - Attempt to use input image mask bigger than mask
    229     CREATE_AND_SET_IMAGE(inImage8,F64,0.0,50,50);
    230     CREATE_AND_SET_IMAGE(mask8,U8,0,50,50)
    231     CREATE_AND_SET_IMAGE(mask8i,U8,0,60,60)
    232     //    pmReadout *inReadout8 = pmReadoutAlloc(0, 0, inImage8);
    233     pmReadout *inReadout8 = pmReadoutAlloc(NULL);
    234     inReadout8->row0 = 0;
    235     inReadout8->col0 = 0;
    236     inReadout8->image = inImage8;
    237     inReadout8->mask = mask8i;
    238     pmMaskBadPixels(inReadout8, mask8, 0, 100.0, 1, 10);
    239     psFree(mask8);
    240     psFree(inReadout8);
     287    //
     288    // Test H - Attempt to use mask bigger than image
     289    //
     290
     291    pmReadout *inReadout = pmReadoutAlloc(NULL);
     292    CREATE_AND_SET_IMAGE(inReadout->image, F32, DEFAULT_IMAGE_VAL, NUM_ROWS, NUM_COLS);
     293    inReadout->image->row0 = 0;
     294    inReadout->image->col0 = 0;
     295    inReadout->row0 = 0;
     296    inReadout->col0 = 0;
     297
     298    pmReadout *maskReadout = pmReadoutAlloc(NULL);
     299    CREATE_AND_SET_IMAGE(maskReadout->image, U8, DEFAULT_MASK_VAL, NUM_ROWS+10, NUM_COLS+10)
     300
     301    pmMaskBadPixels(inReadout, maskReadout, MASK_VAL, SAT_VAL, GROW_VAL, GROW_RAD);
     302
     303    psFree(inReadout);
     304    psFree(maskReadout);
    241305
    242306    return 0;
     
    245309int testMaskBadPixels9( void )
    246310{
     311    //
    247312    // Test I - Attempt to use offset greater than input image
    248     CREATE_AND_SET_IMAGE(inImage9,F64,0.0,50,50);
    249     CREATE_AND_SET_IMAGE(mask9,U8,0,50,50)
    250     CREATE_AND_SET_IMAGE(mask9i,U8,0,50,50)
    251     //    pmReadout *inReadout9 = pmReadoutAlloc(0, 0, inImage9);
    252     pmReadout *inReadout9 = pmReadoutAlloc(NULL);
    253     inReadout9->image = inImage9;
    254     inReadout9->mask = mask9i;
    255     inReadout9->row0 = 0;
    256     inReadout9->col0 = 0;
    257     *(int*)&inReadout9->col0 = 150;
    258     *(int*)&inReadout9->row0 = 150;
    259     pmMaskBadPixels(inReadout9, mask9, 0, 100.0, 1, 10);
    260     psFree(mask9);
    261     psFree(inReadout9);
     313    //
     314
     315    pmReadout *inReadout = pmReadoutAlloc(NULL);
     316    CREATE_AND_SET_IMAGE(inReadout->image, F32, DEFAULT_IMAGE_VAL, NUM_ROWS, NUM_COLS);
     317    inReadout->image->row0 = 0;
     318    inReadout->image->col0 = 0;
     319    inReadout->row0 = 0;
     320    inReadout->col0 = 0;
     321    *(int*)&inReadout->image->col0 = 150;
     322    *(int*)&inReadout->image->row0 = 150;
     323
     324    pmReadout *maskReadout = pmReadoutAlloc(NULL);
     325    CREATE_AND_SET_IMAGE(maskReadout->image, U8, DEFAULT_MASK_VAL, NUM_ROWS, NUM_COLS)
     326    pmMaskBadPixels(inReadout, maskReadout, MASK_VAL, SAT_VAL, GROW_VAL, GROW_RAD);
     327
     328    psFree(inReadout);
     329    psFree(maskReadout);
    262330
    263331    return 0;
     
    266334int testMaskBadPixels10( void )
    267335{
     336    //
    268337    // Test J - Attempt to use complex input image
    269     CREATE_AND_SET_IMAGE(inImage10,C64,50.0,50,50);
    270     CREATE_AND_SET_IMAGE(mask10,U8,0,50,50)
    271     CREATE_AND_SET_IMAGE(mask10i,U8,0,50,50)
    272     //    pmReadout *inReadout10 = pmReadoutAlloc(0, 0, inImage10);
    273     pmReadout *inReadout10 = pmReadoutAlloc(NULL);
    274     inReadout10->row0 = 0;
    275     inReadout10->col0 = 0;
    276     inReadout10->image = inImage10;
    277     inReadout10->mask = mask10i;
    278     pmMaskBadPixels(inReadout10, mask10, 0, 100.0, 1, 10);
    279     psFree(mask10);
    280     psFree(inReadout10);
     338    //
     339
     340    pmReadout *inReadout = pmReadoutAlloc(NULL);
     341    CREATE_AND_SET_IMAGE(inReadout->image, C64, DEFAULT_IMAGE_VAL, NUM_ROWS, NUM_COLS);
     342    inReadout->image->row0 = 0;
     343    inReadout->image->col0 = 0;
     344    inReadout->row0 = 0;
     345    inReadout->col0 = 0;
     346
     347    pmReadout *maskReadout = pmReadoutAlloc(NULL);
     348    CREATE_AND_SET_IMAGE(maskReadout->image, U8, DEFAULT_MASK_VAL, NUM_ROWS, NUM_COLS)
     349
     350    pmMaskBadPixels(inReadout, maskReadout, MASK_VAL, SAT_VAL, GROW_VAL, GROW_RAD);
     351
     352    psFree(inReadout);
     353    psFree(maskReadout);
    281354
    282355    return 0;
     
    285358int testMaskBadPixels11( void )
    286359{
    287     // Test K - Attempt to use non-mask type mask image
    288     CREATE_AND_SET_IMAGE(inImage11,F64,50.0,50,50);
    289     CREATE_AND_SET_IMAGE(mask11,F64,0,50,50)
    290     CREATE_AND_SET_IMAGE(mask11i,U8,0,50,50)
    291     //    pmReadout *inReadout11 = pmReadoutAlloc(0, 0, inImage11);
    292     pmReadout *inReadout11 = pmReadoutAlloc(NULL);
    293     inReadout11->row0 = 0;
    294     inReadout11->col0 = 0;
    295     inReadout11->image = inImage11;
    296     inReadout11->mask = mask11i;
    297     pmMaskBadPixels(inReadout11, mask11, 0, 100.0, 1, 10);
    298     psFree(mask11);
    299     psFree(inReadout11);
    300 
    301     return 0;
    302 }
    303 
     360    //
     361    // Test K - Attempt to use mask image with wrong data type.
     362    //
     363
     364    pmReadout *inReadout = pmReadoutAlloc(NULL);
     365    CREATE_AND_SET_IMAGE(inReadout->image, F32, DEFAULT_IMAGE_VAL, NUM_ROWS, NUM_COLS);
     366    inReadout->image->row0 = 0;
     367    inReadout->image->col0 = 0;
     368    inReadout->row0 = 0;
     369    inReadout->col0 = 0;
     370
     371    pmReadout *maskReadout = pmReadoutAlloc(NULL);
     372    CREATE_AND_SET_IMAGE(maskReadout->image, F64, DEFAULT_MASK_VAL, NUM_ROWS, NUM_COLS)
     373
     374    pmMaskBadPixels(inReadout, maskReadout, MASK_VAL, SAT_VAL, GROW_VAL, GROW_RAD);
     375
     376    psFree(inReadout);
     377    psFree(maskReadout);
     378
     379    return 0;
     380}
     381
  • trunk/psModules/test/imsubtract/tst_pmSubtractBias.c

    r5188 r5516  
    33 *  @brief Contains the tests for pmSubtractBias.c:
    44 *
    5  * test00: This code will subtract full bias frames from the input image.
    6  * test01: Multiple overscan regions, calculate a scalar statistic for
    7  *  each, then subtract from the input image.
    8  * test02: Calculate a column overscan vector and subtract it from each
    9  *  column in the input image.
    10  * test03: Calculate a row overscan vector and subtract it from each
     5 * test00a: This code will subtract full bias frames from the input image.
     6 * XXX: Must test:
     7 *  Various image offsets.
     8 *  Various image size combinations.
     9 *  Various data types for the bias and input images.
     10 *  Ensure code works when CELL.TRIMSEC is not set.
     11 * test00b: This code will subtract full dark frames from the input image.
     12 * XXX: Must test:
     13 *  Various image offsets.
     14 *  Various image size combinations.
     15 *  Various data types for the bias and input images.
     16 *  Code properly determines CELL.DARKTIME from cell metadata.
     17 *  Ensure code works when CELL.DARKTIME is not set.
     18 *  Ensure code works when CELL.TRIMSEC is not set.
     19 *  test03: Calculate a row overscan vector and subtract it from each
    1120 *  row in the input image.
    12  * test04:
     21 * test05:
    1322 *
    1423 *  @author GLG, MHPCC
     
    1625 *  XXX: Memory leaks are not being detected.
    1726 *
    18  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    19  *  @date $Date: 2005-09-29 21:57:31 $
     27 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     28 *  @date $Date: 2005-11-15 20:09:03 $
    2029 *
    2130 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    2534#include "pslib.h"
    2635#include "pmSubtractBias.h"
    27 static int test00(void);
    28 static int test01(void);
    29 static int test02(void);
    30 static int test03(void);
    31 static int test04(void);
    32 static int testX(void);
     36static int test00a(void);
     37static int test00b(void);
     38//static int test01(void);
     39//static int test02(void);
     40//static int test03(void);
     41//static int test04(void);
     42static int test05(void);
    3343testDescription tests[] = {
    34                               {test00, 000, "pmSubtractBias", 0, true},
    35                               {test01, 000, "pmSubtractBias", 0, true},
    36                               {test02, 000, "pmSubtractBias", 0, true},
    37                               {test03, 000, "pmSubtractBias", 0, false},
    38                               {test04, 000, "pmSubtractBias", 0, true},
    39                               {testX,  000, "pmSubtractBias", 0, true},
     44                              {test00a, 000, "doSubtractBiasFullFrame", 0, false},
     45                              {test00b, 000, "doSubtractDarkFullFrame", 0, false},
     46                              //                              {test01, 000, "pmSubtractBias", 0, true},
     47                              //                              {test02, 000, "pmSubtractBias", 0, true},
     48                              //                              {test03, 000, "pmSubtractBias", 0, false},
     49                              //                              {test04, 000, "pmSubtractBias", 0, true},
     50                              {test05, 000, "pmSubtractBias", 0, false},
    4051                              {NULL}
    4152                          };
     
    4556{
    4657    psLogSetFormat("HLNM");
     58
     59    psTraceSetLevel(".", 0);
     60    psTraceSetLevel("spline1DFree", 0);
     61    psTraceSetLevel("calculateSecondDerivs", 0);
     62    psTraceSetLevel("vectorBinDisectF32", 0);
     63    psTraceSetLevel("vectorBinDisectF64", 0);
     64    psTraceSetLevel("p_psVectorBinDisect", 0);
     65    psTraceSetLevel("psSpline1DAlloc", 0);
     66    psTraceSetLevel("psVectorFitSpline1D", 0);
     67    psTraceSetLevel("psSpline1DEval", 0);
     68    psTraceSetLevel("psSpline1DEvalVector", 0);
     69
    4770    return !runTestSuite(stderr, "Test Point Driver", tests, argc, argv);
    4871}
     
    5073#define NUM_ROWS 8
    5174#define NUM_COLS 8
     75#define MAX_HEADER_MSG_LENGTH 1000
     76#define POLYNOMIAL_FIT_ORDER 2
     77#define NUM_OVERSCANS 2
    5278/******************************************************************************
    5379doSubtractBiasFullFrame(): a sample pmReadout as well as a bias image are
     
    6389    psImage *tmpImage1 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    6490    psImage *tmpImage2 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    65     //    pmReadout *myReadout = pmReadoutAlloc(numCols, numRows, tmpImage1);
    66     //    pmReadout *myBias = pmReadoutAlloc(numCols, numRows, tmpImage2);
    6791    pmReadout *myReadout = pmReadoutAlloc(NULL);
    6892    pmReadout *myBias = pmReadoutAlloc(NULL);
     
    7094    myBias->image = tmpImage2;
    7195
    72 
    73     printPositiveTestHeader(stdout, "pmSubtractBias", "doSubtractBiasFullFrame");
     96    char *HeaderMessageStr = (char *) psAlloc(MAX_HEADER_MSG_LENGTH);
     97    sprintf(HeaderMessageStr, "doSubtractBiasFullFrame(%d, %d)", numRows, numCols);
     98    printPositiveTestHeader(stdout, "pmSubtractBias", HeaderMessageStr);
    7499    for (i=0;i<numRows;i++) {
    75100        for (j=0;j<numCols;j++) {
     
    79104    }
    80105
    81     myReadout = pmSubtractBias(myReadout, NULL, NULL, PM_OVERSCAN_NONE, NULL,
    82                                0, PM_FIT_NONE, myBias);
     106    myReadout = pmSubtractBias(myReadout, NULL, PM_FIT_NONE, false,
     107                               NULL, 0, myBias, NULL);
    83108
    84109    for (i=0;i<numRows;i++) {
     
    96121    psFree(myReadout);
    97122    psFree(myBias);
    98     printFooter(stdout, "pmSubtractBias", "doSubtractBiasFullFrame", true);
     123    printFooter(stdout, "pmSubtractBias", HeaderMessageStr, true);
     124    psFree(HeaderMessageStr);
    99125    return(testStatus);
    100126}
    101127
    102128
    103 int test00( void )
     129int test00a( void )
    104130{
    105131    int testStatus = 0;
     
    112138}
    113139
     140
    114141/******************************************************************************
    115 doSubtractFullOverscans(): a sample pmReadout as well as several overscan
    116 images of the same size are created.  The overscan images are then subtracted
    117 from the pmReadout.
     142doSubtractDarkFullFrame(): a sample pmReadout as well as a dark image are
     143created and the dark image is subtracted from the pmReadout.
    118144 *****************************************************************************/
    119 int doSubtractFullOverscans(int numCols, int numRows)
     145int doSubtractDarkFullFrame(int numCols, int numRows)
    120146{
    121147    int i;
     
    126152    psImage *tmpImage1 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    127153    psImage *tmpImage2 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    128     psImage *tmpImage3 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    129     psImage *tmpImage4 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    130     //    pmReadout *myReadout = pmReadoutAlloc(numCols, numRows, tmpImage1);
    131154    pmReadout *myReadout = pmReadoutAlloc(NULL);
     155    pmReadout *myDark = pmReadoutAlloc(NULL);
    132156    myReadout->image = tmpImage1;
    133 
    134 
    135     printPositiveTestHeader(stdout, "pmSubtractBias", "doSubtractFullOverscans");
    136     psList *list;
    137     psStats *stat = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
    138 
     157    myDark->image = tmpImage2;
     158
     159    char *HeaderMessageStr = (char *) psAlloc(MAX_HEADER_MSG_LENGTH);
     160    sprintf(HeaderMessageStr, "doSubtractDarkFullFrame(%d, %d)", numRows, numCols);
     161    printPositiveTestHeader(stdout, "pmSubtractBias", HeaderMessageStr);
    139162    for (i=0;i<numRows;i++) {
    140163        for (j=0;j<numCols;j++) {
    141164            myReadout->image->data.F32[i][j] = (float) (i + j);
    142             tmpImage2->data.F32[i][j] = 3.0;
    143             tmpImage3->data.F32[i][j] = 4.0;
    144             tmpImage4->data.F32[i][j] = 5.0;
    145         }
    146     }
    147     list = psListAlloc(tmpImage2);
    148     psListAdd(list, PS_LIST_HEAD, tmpImage3);
    149     psListAdd(list, PS_LIST_HEAD, tmpImage4);
    150 
    151     myReadout = pmSubtractBias(myReadout, NULL, list, PM_OVERSCAN_ALL, stat,
    152                                0, PM_FIT_NONE, NULL);
     165            myDark->image->data.F32[i][j] = 1.0;
     166        }
     167    }
     168
     169    myReadout = pmSubtractBias(myReadout, NULL, PM_FIT_NONE, false,
     170                               NULL, 0, NULL, myDark);
    153171
    154172    for (i=0;i<numRows;i++) {
    155173        for (j=0;j<numCols;j++) {
    156             expect = ((float) (i + j)) - 12.0;
     174            expect = ((float) (i + j)) - 1.0;
    157175            actual = myReadout->image->data.F32[i][j];
    158176            if (FLT_EPSILON < fabs(expect - actual)) {
     
    165183
    166184    psFree(myReadout);
    167     psFree(tmpImage2);
    168     psFree(tmpImage3);
    169     psFree(tmpImage4);
    170     psFree(stat);
    171     psFree(list);
    172 
    173     printFooter(stdout, "pmSubtractBias", "doSubtractFullOverscans", true);
     185    psFree(myDark);
     186    printFooter(stdout, "pmSubtractBias", HeaderMessageStr, true);
     187    psFree(HeaderMessageStr);
    174188    return(testStatus);
    175189}
    176190
    177 int test01( void )
     191
     192int test00b( void )
    178193{
    179194    int testStatus = 0;
    180195
    181     testStatus |= doSubtractFullOverscans(1, 1);
    182     testStatus |= doSubtractFullOverscans(1, NUM_ROWS);
    183     testStatus |= doSubtractFullOverscans(NUM_COLS, 1);
    184     testStatus |= doSubtractFullOverscans(NUM_COLS, NUM_ROWS);
    185 
     196    testStatus |= doSubtractDarkFullFrame(1, 1);
     197    testStatus |= doSubtractDarkFullFrame(NUM_COLS, 1);
     198    testStatus |= doSubtractDarkFullFrame(1, NUM_ROWS);
     199    testStatus |= doSubtractDarkFullFrame(NUM_COLS, NUM_ROWS);
    186200    return(testStatus);
    187201}
    188202
    189 /******************************************************************************
    190 doSubtractFullOverscans(): a sample pmReadout as well as several overscan
    191 images of the same size are created.  The overscan images are collected
    192 pixel-by-pixel then subtracted column-wise from the pmReadout.
    193  *****************************************************************************/
    194 int doSubtractFullOverscanColumns(int numCols, int numRows)
    195 {
    196     int i;
    197     int j;
    198     float actual;
    199     float expect;
    200     int testStatus = 0;
    201     psImage *tmpImage1 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    202     psImage *tmpImage2 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    203     psImage *tmpImage3 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    204     psImage *tmpImage4 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    205     pmReadout *myReadout = pmReadoutAlloc(NULL);
    206     myReadout->image = tmpImage1;
    207     psList *list;
    208     psStats *stat = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
    209 
    210     for (i=0;i<numRows;i++) {
    211         for (j=0;j<numCols;j++) {
    212             myReadout->image->data.F32[i][j] = (float) (i + j);
    213             tmpImage2->data.F32[i][j] = 3.0;
    214             tmpImage3->data.F32[i][j] = 4.0;
    215             tmpImage4->data.F32[i][j] = 5.0;
    216         }
    217     }
    218     list = psListAlloc(tmpImage2);
    219     psListAdd(list, PS_LIST_HEAD, tmpImage3);
    220     psListAdd(list, PS_LIST_HEAD, tmpImage4);
    221 
    222     printPositiveTestHeader(stdout, "pmSubtractBias", "Column Overscans");
    223 
    224     myReadout = pmSubtractBias(myReadout, NULL, list, PM_OVERSCAN_COLUMNS, stat,
    225                                0, PM_FIT_NONE, NULL);
    226 
    227     for (i=0;i<numRows;i++) {
    228         for (j=0;j<numCols;j++) {
    229             expect = ((float) (i + j)) - 12.0;
    230             actual = myReadout->image->data.F32[i][j];
    231             if (FLT_EPSILON < fabs(expect - actual)) {
    232                 printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
    233                 testStatus = 1;
    234             }
    235         }
    236     }
    237 
    238 
    239     psFree(myReadout);
    240     psFree(tmpImage2);
    241     psFree(tmpImage3);
    242     psFree(tmpImage4);
    243     psFree(stat);
    244     psFree(list);
    245 
    246     printFooter(stdout, "pmSubtractBias", "Column Overscans", true);
    247     return(testStatus);
    248 }
    249 /******************************************************************************
    250 doSubtractFullOverscans(): a sample pmReadout as well as several overscan
    251 images of the same size are created.  The overscan images are collected
    252 pixel-by-pixel then subtracted column-wise from the pmReadout.
    253  *****************************************************************************/
    254 int doSubtractFullOverscanColumnsPoly(int numCols, int numRows)
    255 {
    256     int i;
    257     int j;
    258     float actual;
    259     float expect;
    260     int testStatus = 0;
    261     psImage *tmpImage1 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    262     psImage *tmpImage2 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    263     psImage *tmpImage3 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    264     psImage *tmpImage4 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    265     pmReadout *myReadout = pmReadoutAlloc(NULL);
    266     myReadout->image = tmpImage1;
    267     psList *list;
    268     psStats *stat = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
    269     psPolynomial1D *myPoly = psPolynomial1DAlloc(1, PS_POLYNOMIAL_ORD);
    270 
    271     for (i=0;i<numRows;i++) {
    272         for (j=0;j<numCols;j++) {
    273             myReadout->image->data.F32[i][j] = (float) (i + j);
    274             tmpImage2->data.F32[i][j] = 3.0;
    275             tmpImage3->data.F32[i][j] = 4.0;
    276             tmpImage4->data.F32[i][j] = 5.0;
    277         }
    278     }
    279     list = psListAlloc(tmpImage2);
    280     psListAdd(list, PS_LIST_HEAD, tmpImage3);
    281     psListAdd(list, PS_LIST_HEAD, tmpImage4);
    282 
    283     printPositiveTestHeader(stdout, "pmSubtractBias", "Column Overscans");
    284 
    285     myReadout = pmSubtractBias(myReadout, myPoly, list, PM_OVERSCAN_COLUMNS, stat,
    286                                0, PM_FIT_POLYNOMIAL, NULL);
    287 
    288     for (i=0;i<numRows;i++) {
    289         for (j=0;j<numCols;j++) {
    290             expect = ((float) (i + j)) - 12.0;
    291             actual = myReadout->image->data.F32[i][j];
    292             if (FLT_EPSILON < fabs(expect - actual)) {
    293                 printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
    294                 testStatus = 1;
    295             }
    296         }
    297     }
    298 
    299 
    300     psFree(myReadout);
    301     psFree(tmpImage2);
    302     psFree(tmpImage3);
    303     psFree(tmpImage4);
    304     psFree(stat);
    305     psFree(list);
    306     psFree(myPoly);
    307 
    308     printFooter(stdout, "pmSubtractBias", "Column Overscans", true);
    309     return(testStatus);
    310 }
    311 /******************************************************************************
    312 doSubtractFullOverscansSmall(): a sample pmReadout as well as several overscan
    313 images of smaller size are created.  The overscan images are collected
    314 pixel-by-pixel then subtracted column-wise from the pmReadout.
    315  *****************************************************************************/
    316 int doSubtractFullOverscanColumnsSmall(int numCols, int numRows)
    317 {
    318     int i;
    319     int j;
    320     float actual;
    321     float expect;
    322     int testStatus = 0;
    323     psImage *tmpImage1 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    324     pmReadout *myReadout = pmReadoutAlloc(NULL);
    325     myReadout->image = tmpImage1;
    326     psImage *tmpImage2 = psImageAlloc(numCols/2, numRows/2, PS_TYPE_F32);
    327     psImage *tmpImage3 = psImageAlloc(numCols/2, numRows/2, PS_TYPE_F32);
    328     psImage *tmpImage4 = psImageAlloc(numCols/2, numRows/2, PS_TYPE_F32);
    329     psList *list;
    330     psStats *stat = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
    331 
    332     for (i=0;i<numRows;i++) {
    333         for (j=0;j<numCols;j++) {
    334             myReadout->image->data.F32[i][j] = (float) (i + j);
    335         }
    336     }
    337     for (i=0;i<numRows/2;i++) {
    338         for (j=0;j<numCols/2;j++) {
    339             tmpImage2->data.F32[i][j] = 3.0;
    340             tmpImage3->data.F32[i][j] = 4.0;
    341             tmpImage4->data.F32[i][j] = 5.0;
    342         }
    343     }
    344     list = psListAlloc(tmpImage2);
    345     psListAdd(list, PS_LIST_HEAD, tmpImage3);
    346     psListAdd(list, PS_LIST_HEAD, tmpImage4);
    347 
    348     printPositiveTestHeader(stdout, "pmSubtractBias", "Column Overscans");
    349 
    350     myReadout = pmSubtractBias(myReadout, NULL, list, PM_OVERSCAN_COLUMNS, stat,
    351                                0, PM_FIT_POLYNOMIAL, NULL);
    352 
    353     for (i=0;i<numRows;i++) {
    354         for (j=0;j<numCols;j++) {
    355             expect = ((float) (i + j)) - 12.0;
    356             actual = myReadout->image->data.F32[i][j];
    357             if (FLT_EPSILON < fabs(expect - actual)) {
    358                 printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
    359                 testStatus = 1;
    360             }
    361         }
    362     }
    363 
    364 
    365     psFree(myReadout);
    366     psFree(tmpImage2);
    367     psFree(tmpImage3);
    368     psFree(tmpImage4);
    369     psFree(stat);
    370     psFree(list);
    371 
    372     printFooter(stdout, "pmSubtractBias", "Column Overscans", true);
    373     return(testStatus);
    374 }
    375 
    376 int test02( void )
    377 {
    378     int testStatus = 0;
    379 
    380     testStatus |= doSubtractFullOverscanColumns(1, 1);
    381     testStatus |= doSubtractFullOverscanColumns(1, NUM_ROWS);
    382     testStatus |= doSubtractFullOverscanColumns(NUM_COLS, 1);
    383     testStatus |= doSubtractFullOverscanColumns(NUM_COLS, NUM_ROWS);
    384     /* These tests do not make sense until the SDRS is clarified.
    385         testStatus |= doSubtractFullOverscanColumnsSmall(1, 1);
    386         testStatus |= doSubtractFullOverscanColumnsSmall(1, NUM_ROWS);
    387         testStatus |= doSubtractFullOverscanColumnsSmall(NUM_COLS, 1);
    388         testStatus |= doSubtractFullOverscanColumnsSmall(NUM_COLS, NUM_ROWS);
    389     */
    390 
    391     return(testStatus);
    392 }
    393 
    394 /******************************************************************************
    395 doSubtractFullOverscans(): a sample pmReadout as well as several overscan
    396 images of the same size are created.  The overscan images are collected
    397 pixel-by-pixel then subtracted row-wise from the pmReadout.
    398  *****************************************************************************/
    399 int doSubtractFullOverscanRows(int numCols, int numRows)
    400 {
    401     int i;
    402     int j;
    403     float actual;
    404     float expect;
    405     int testStatus = 0;
    406     psImage *tmpImage1 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    407     psImage *tmpImage2 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    408     psImage *tmpImage3 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    409     psImage *tmpImage4 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    410     pmReadout *myReadout = pmReadoutAlloc(NULL);
    411     myReadout->image = tmpImage1;
    412     psList *list;
    413     psStats *stat = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
    414 
    415     for (i=0;i<numRows;i++) {
    416         for (j=0;j<numCols;j++) {
    417             myReadout->image->data.F32[i][j] = (float) (i + j);
    418             tmpImage2->data.F32[i][j] = 3.0;
    419             tmpImage3->data.F32[i][j] = 4.0;
    420             tmpImage4->data.F32[i][j] = 5.0;
    421         }
    422     }
    423     list = psListAlloc(tmpImage2);
    424     psListAdd(list, PS_LIST_HEAD, tmpImage3);
    425     psListAdd(list, PS_LIST_HEAD, tmpImage4);
    426 
    427     printPositiveTestHeader(stdout, "pmSubtractBias", "Row Overscans");
    428 
    429     myReadout = pmSubtractBias(myReadout, NULL, list, PM_OVERSCAN_ROWS, stat,
    430                                0, PM_FIT_NONE, NULL);
    431 
    432     for (i=0;i<numRows;i++) {
    433         for (j=0;j<numCols;j++) {
    434             expect = ((float) (i + j)) - 12.0;
    435             actual = myReadout->image->data.F32[i][j];
    436             if (FLT_EPSILON < fabs(expect - actual)) {
    437                 printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
    438                 testStatus = 1;
    439             }
    440         }
    441     }
    442 
    443 
    444     psFree(myReadout);
    445     psFree(tmpImage2);
    446     psFree(tmpImage3);
    447     psFree(tmpImage4);
    448     psFree(stat);
    449     psFree(list);
    450 
    451     printFooter(stdout, "pmSubtractBias", "Row Overscans", true);
    452     return(testStatus);
    453 }
    454 
    455 /******************************************************************************
    456 doSubtractFullOverscansSmall(): a sample pmReadout as well as several overscan
    457 images of smaller size are created.  The overscan images are collected
    458 pixel-by-pixel then subtracted row-wise from the pmReadout.
    459  *****************************************************************************/
    460 int doSubtractFullOverscanRowsSmall(int numCols, int numRows)
    461 {
    462     int i;
    463     int j;
    464     float actual;
    465     float expect;
    466     int testStatus = 0;
    467     psS32 OSnumRows = numRows-1;
    468     if (OSnumRows == 0) {
    469         OSnumRows = 1;
    470     }
    471     psS32 OSnumCols = numCols-1;
    472     if (OSnumCols == 0) {
    473         OSnumCols = 1;
    474     }
    475     psImage *tmpImage1 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    476     psImage *tmpImage2 = psImageAlloc(OSnumCols, OSnumRows, PS_TYPE_F32);
    477     psImage *tmpImage3 = psImageAlloc(OSnumCols, OSnumRows, PS_TYPE_F32);
    478     psImage *tmpImage4 = psImageAlloc(OSnumCols, OSnumRows, PS_TYPE_F32);
    479     pmReadout *myReadout = pmReadoutAlloc(NULL);
    480     myReadout->image = tmpImage1;
    481     psList *list;
    482     psStats *stat = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
    483 
    484     for (i=0;i<numRows;i++) {
    485         for (j=0;j<numCols;j++) {
    486             myReadout->image->data.F32[i][j] = (float) (i + j);
    487         }
    488     }
    489     for (i=0;i<OSnumRows;i++) {
    490         for (j=0;j<OSnumCols;j++) {
    491             tmpImage2->data.F32[i][j] = 3.0;
    492             tmpImage3->data.F32[i][j] = 4.0;
    493             tmpImage4->data.F32[i][j] = 5.0;
    494         }
    495     }
    496     list = psListAlloc(tmpImage2);
    497     psListAdd(list, PS_LIST_HEAD, tmpImage3);
    498     psListAdd(list, PS_LIST_HEAD, tmpImage4);
    499 
    500     printPositiveTestHeader(stdout, "pmSubtractBias", "Row Overscans");
    501 
    502     myReadout = pmSubtractBias(myReadout, NULL, list, PM_OVERSCAN_ROWS, stat,
    503                                0, PM_FIT_SPLINE, NULL);
    504 
    505     for (i=0;i<numRows;i++) {
    506         for (j=0;j<numCols;j++) {
    507             expect = ((float) (i + j)) - 12.0;
    508             actual = myReadout->image->data.F32[i][j];
    509             if (FLT_EPSILON < fabs(expect - actual)) {
    510                 printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
    511                 testStatus = 1;
    512             }
    513         }
    514     }
    515 
    516 
    517     psFree(myReadout);
    518     psFree(tmpImage2);
    519     psFree(tmpImage3);
    520     psFree(tmpImage4);
    521     psFree(stat);
    522     psFree(list);
    523 
    524     printFooter(stdout, "pmSubtractBias", "Row Overscans", true);
    525     return(testStatus);
    526 }
    527 
    528 
    529 // XXX: HEY
    530 int test03( void )
    531 {
    532     int testStatus = 0;
    533 
    534     //    testStatus |= doSubtractFullOverscanRows(1, 1);
    535     //    testStatus |= doSubtractFullOverscanRows(1, NUM_ROWS);
    536     //    testStatus |= doSubtractFullOverscanRows(NUM_COLS, 1);
    537     //    testStatus |= doSubtractFullOverscanRows(NUM_COLS, NUM_ROWS);
    538     //    testStatus |= doSubtractFullOverscanRowsSmall(1, 1);
    539     //    testStatus |= doSubtractFullOverscanRowsSmall(1, NUM_ROWS);
    540     //    testStatus |= doSubtractFullOverscanRowsSmall(NUM_COLS, 1);
    541     testStatus |= doSubtractFullOverscanRowsSmall(NUM_COLS, NUM_ROWS);
    542 
    543     return(testStatus);
    544 }
    545 
    546 
    547 
     203
     204/*
    548205int doSubtractOverscansTestInputCases(int numCols, int numRows)
    549206{
     
    568225    myBias->image = tmpImage5;
    569226    printPositiveTestHeader(stdout, "pmSubtractBias", "Testing input parameter error conditions");
    570 
     227 
    571228    psImage *tmpImage5ShortRows = psImageAlloc(numCols, numRows-1, PS_TYPE_F32);
    572229    pmReadout *myBiasShortRows = pmReadoutAlloc(NULL);
     
    575232    pmReadout *myBiasShortCols = pmReadoutAlloc(NULL);
    576233    myBiasShortCols->image = tmpImage5ShortCols;
    577 
     234 
    578235    for (i=0;i<numRows;i++) {
    579236        for (j=0;j<numCols;j++) {
     
    588245    psListAdd(list, PS_LIST_HEAD, tmpImage3);
    589246    psListAdd(list, PS_LIST_HEAD, tmpImage4);
    590 
     247 
    591248    for (i=0;i<numRows-1;i++) {
    592249        for (j=0;j<numCols-1;j++) {
     
    609266        }
    610267    }
    611 
    612 
     268 
     269 
    613270    printf("------------------------------------------------------------------\n");
    614271    printf("Calling pmSubtractBias() with NULL overscan list and PM_OVERSCAN_ALL.  Should generate error.\n");
     
    618275        testStatus = false;
    619276    }
    620 
     277 
    621278    printf("------------------------------------------------------------------\n");
    622279    printf("Calling pmSubtractBias() with NULL overscan list and PM_OVERSCAN_ROWS.  Should generate error.\n");
     
    627284        psFree(rc);
    628285    }
    629 
     286 
    630287    printf("------------------------------------------------------------------\n");
    631288    printf("Calling pmSubtractBias() with NULL overscan list and PM_OVERSCAN_COLUMNS.  Should generate error.\n");
     
    636293        psFree(rc);
    637294    }
    638 
     295 
    639296    printf("------------------------------------------------------------------\n");
    640297    printf("Calling pmSubtractBias() with non-NULL overscan list and PM_OVERSCAN_NONE.  Should generate warning.\n");
    641298    rc = pmSubtractBias(myReadout, NULL, list, PM_OVERSCAN_NONE, stat, 0, PM_FIT_NONE, myBias);
    642 
     299 
    643300    for (i=0;i<numRows;i++) {
    644301        for (j=0;j<numCols;j++) {
     
    649306                testStatus = 1;
    650307            }
    651 
     308 
    652309            // Restore myReadout for next test.
    653310            myReadout->image->data.F32[i][j] = (float) (i + j);
    654311        }
    655312    }
    656 
    657     /* XXX: This does not seem to be a requirement.
     313 
     314    // XXX: This does not seem to be a requirement.
     315    if (0) {
    658316        printf("------------------------------------------------------------------\n");
    659317        printf("Calling pmSubtractBias() with NULL overscan list and PM_OVERSCAN_NONE.  Should generate warning.\n");
     
    674332            }
    675333        }
    676     */
    677 
     334    }
     335 
    678336    printf("------------------------------------------------------------------\n");
    679337    printf("Calling pmSubtractBias() with PM_OVERSCAN_NONE and PM_FIT_POLYNOMIAL.  Should generate Warning.\n");
     
    687345                testStatus = 1;
    688346            }
    689 
     347 
    690348            // Restore myReadout for next test.
    691349            myReadout->image->data.F32[i][j] = (float) (i + j);
    692350        }
    693351    }
    694 
    695 
     352 
     353 
    696354    printf("------------------------------------------------------------------\n");
    697355    printf("Calling pmSubtractBias() with PM_OVERSCAN_ALL and PM_FIT_SPLINE.  Should generate Warning.\n");
     
    702360        psFree(rc);
    703361    }
    704 
     362 
    705363    for (i=0;i<numRows;i++) {
    706364        for (j=0;j<numCols;j++) {
     
    711369                testStatus = 1;
    712370            }
    713 
     371 
    714372            // Restore myReadout for next test.
    715373            myReadout->image->data.F32[i][j] = (float) (i + j);
    716374        }
    717375    }
    718 
    719 
     376 
     377 
    720378    printf("------------------------------------------------------------------\n");
    721379    printf("Calling pmSubtractBias() with multiple stats->options.  Should generate Warning.\n");
     
    734392    }
    735393    stat->options = PS_STAT_SAMPLE_MEAN;
    736 
     394 
    737395    printf("------------------------------------------------------------------\n");
    738396    printf("Calling pmSubtractBias() undersize overscans (PM_OVERSCAN_ROWS).  Should generate Warning.\n");
    739397    rc = pmSubtractBias(myReadout, NULL, listShort, PM_OVERSCAN_ROWS, stat,
    740398                        0, PM_FIT_NONE, NULL);
    741     /*
     399    if (0) {
    742400        for (i=0;i<numRows;i++) {
    743401            for (j=0;j<numCols;j++) {
     
    750408            }
    751409        }
    752     */
    753 
     410    }
     411 
    754412    printf("------------------------------------------------------------------\n");
    755413    printf("Calling pmSubtractBias() undersize overscans (PM_OVERSCAN_COLUMNS).  Should generate Warning.\n");
    756414    rc = pmSubtractBias(myReadout, NULL, listShort, PM_OVERSCAN_COLUMNS, stat,
    757415                        0, PM_FIT_NONE, NULL);
    758     /*
     416    if (0) {
    759417        for (i=0;i<numRows;i++) {
    760418            for (j=0;j<numCols;j++) {
     
    767425            }
    768426        }
    769     */
    770 
     427    }
     428 
    771429    printf("------------------------------------------------------------------\n");
    772430    printf("Calling pmSubtractBias() undersize bias image (short rows).  Should generate Error.\n");
     
    778436        psFree(rc);
    779437    }
    780 
    781 
     438 
     439 
    782440    printf("------------------------------------------------------------------\n");
    783441    printf("Calling pmSubtractBias() undersize bias image (short columns).  Should generate Error.\n");
     
    789447        psFree(rc);
    790448    }
    791 
     449 
    792450    printf("------------------------------------------------------------------\n");
    793451    printf("Calling pmSubtractBias() with bogus PM_FIT.  Should generate Error.\n");
     
    799457        psFree(rc);
    800458    }
    801 
     459 
    802460    printf("------------------------------------------------------------------\n");
    803461    printf("Calling pmSubtractBias() with bogus overScanAxis.  Should generate Error.\n");
     
    809467        psFree(rc);
    810468    }
    811 
    812     /*
     469 
     470    if (0) {
    813471        for (i=0;i<numRows;i++) {
    814472            for (j=0;j<numCols;j++) {
     
    821479            }
    822480        }
    823     */
    824 
     481    }
     482 
    825483    printf("------------------------------------------------------------------\n");
    826484    psFree(myReadout);
     
    837495    psFree(list);
    838496    psFree(listShort);
    839 
     497 
    840498    printFooter(stdout, "pmSubtractBias", "Testing input parameter error conditions", true);
    841499    return(testStatus);
    842500}
    843 
     501 
    844502int test04( void )
    845503{
    846504    int testStatus = 0;
    847 
     505 
    848506    testStatus |= doSubtractOverscansTestInputCases(NUM_COLS, NUM_ROWS);
    849507    return(testStatus);
    850508}
     509*/
    851510
    852511/******************************************************************************
    853 doSubtractFullOverscansSmall(): a sample pmReadout as well as several overscan
    854 images of smaller size are created.  The overscan images are collected
    855 pixel-by-pixel then subtracted column-wise from the pmReadout.
     512doSubtractFullOverscanColumnsGeneric(): This is a general version of the
     513bias subtraction tests which allows the various parameters to be specified
     514as arguments.
    856515 *****************************************************************************/
    857 int doSubtractFullOverscanColumnsGeneric(int imageNumCols,
    858         int imageNumRows,
    859         int overscanNumCols,
    860         int overscanNumRows,
    861         pmOverscanAxis overscanaxis,
    862         pmFit fit,
    863         psS32 nBin)
     516int doSubtractFullOverscanColumnsGeneric(
     517    int imageNumCols,
     518    int imageNumRows,
     519    int overscanNumCols,
     520    int overscanNumRows,
     521    int numOverscans,
     522    pmOverscanAxis overscanaxis,
     523    pmFit fit,
     524    psS32 nBin)
    864525{
    865526    int i;
     
    869530    int testStatus = 0;
    870531
    871     psImage *tmpImage1 = psImageAlloc(imageNumCols, imageNumRows, PS_TYPE_F32);
    872     //    pmReadout *myReadout = pmReadoutAlloc(imageNumCols, imageNumRows, tmpImage1);
    873     pmReadout *myReadout = pmReadoutAlloc(NULL);
    874     myReadout->image = tmpImage1;
     532    printPositiveTestHeader(stdout, "pmSubtractBias", "PUT COMMENT HERE");
     533    printf("---- doSubtractFullOverscanColumnsGeneric() ----\n");
     534    printf("    Image size: %d by %d\n", imageNumRows, imageNumCols);
     535    printf("    Overscan size: %d by %d\n", overscanNumRows, overscanNumCols);
     536    printf("    Total Overscans: %d\n", numOverscans);
     537    printf("    Binning factor: %d\n", nBin);
     538    if (overscanaxis == PM_OVERSCAN_ROWS)
     539        printf("    Overscan axis: PM_OVERSCAN_ROWS\n");
     540    if (overscanaxis == PM_OVERSCAN_COLUMNS)
     541        printf("    Overscan axis: PM_OVERSCAN_COLUMNS\n");
     542    if (overscanaxis == PM_OVERSCAN_ALL)
     543        printf("    Overscan axis: PM_OVERSCAN_ALL\n");
     544    if (overscanaxis == PM_OVERSCAN_NONE)
     545        printf("    Overscan axis: PM_OVERSCAN_NONE\n");
     546    if (fit == PM_FIT_NONE)
     547        printf("    Fit type: PM_FIT_NONE\n");
     548    if (fit == PM_FIT_POLYNOMIAL)
     549        printf("    Fit type: PM_FIT_POLYNOMIAL\n");
     550    if (fit == PM_FIT_SPLINE)
     551        printf("    Fit type: PM_FIT_SPLINE\n");
     552
     553    //
     554    // Create and initialize input image, FPA hierarchy.
     555    //
     556    const psMetadata *camera = psMetadataAlloc();
     557    pmFPA* fpa = pmFPAAlloc(camera);
     558
     559    if (fpa == NULL) {
     560        psLogMsg(__func__,PS_LOG_ERROR, "TEST ERROR: pmFPAAlloc returned a NULL.\n");
     561        return 1;
     562    }
     563
     564    pmChip *chip = pmChipAlloc(fpa, "ChipName");
     565    if (chip == NULL) {
     566        psLogMsg(__func__,PS_LOG_ERROR, "TEST ERROR: pmChipAlloc returned a NULL.\n");
     567        return 2;
     568    }
     569
     570    pmCell *cell = pmCellAlloc(chip, (psMetadata *) camera, "CellName");
     571    if (cell == NULL) {
     572        psLogMsg(__func__,PS_LOG_ERROR, "TEST ERROR: pmCellAlloc returned a NULL.\n");
     573        return 3;
     574    }
     575
     576    pmReadout *myReadout = pmReadoutAlloc(cell);
     577    myReadout->image = psImageAlloc(imageNumCols, imageNumRows, PS_TYPE_F32);
     578    for (i=0;i<myReadout->image->numRows;i++) {
     579        for (j=0;j<myReadout->image->numCols;j++) {
     580            myReadout->image->data.F32[i][j] = (float) (i + j);
     581        }
     582    }
     583
     584    //
     585    // Set overscan axis in the metadata.
     586    //
     587    psBool rc;
     588    if (overscanaxis == PM_OVERSCAN_ROWS) {
     589        rc = psMetadataAddS32(myReadout->parent->concepts, PS_LIST_HEAD, "CELL.READDIR", 0, NULL, 1);
     590    } else if (overscanaxis == PM_OVERSCAN_COLUMNS) {
     591        rc = psMetadataAddS32(myReadout->parent->concepts, PS_LIST_HEAD, "CELL.READDIR", 0, NULL, 2);
     592    } else if (overscanaxis == PM_OVERSCAN_ALL) {
     593        rc = psMetadataAddS32(myReadout->parent->concepts, PS_LIST_HEAD, "CELL.READDIR", 0, NULL, 3);
     594    } else if (overscanaxis == PM_OVERSCAN_NONE) {
     595        rc = psMetadataAddS32(myReadout->parent->concepts, PS_LIST_HEAD, "CELL.READDIR", 0, NULL, 0);
     596    }
     597    if (rc == false) {
     598        printf("TEST ERROR: Could not set CELL.READDIR metadata.\n");
     599        testStatus = 1;
     600    }
     601
     602    psStats *stat = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
     603    psPolynomial1D *myPoly = psPolynomial1DAlloc(POLYNOMIAL_FIT_ORDER, PS_POLYNOMIAL_ORD);
     604    psSpline1D *mySpline = NULL;
     605
     606
     607    if (0) {
     608        if (overscanNumRows <= 0) {
     609            overscanNumRows = 1;
     610        }
     611        if (overscanNumCols <= 0) {
     612            overscanNumCols = 1;
     613        }
     614    }
     615    psF32 oAverage = 0.0;
     616    myReadout->bias = NULL;
     617    for (psS32 i = 0 ; i < numOverscans ; i++) {
     618        psImage *tmpImage = psImageAlloc(overscanNumCols, overscanNumRows, PS_TYPE_F32);
     619        psF32 oValue = (float) (i + 3);
     620        PS_IMAGE_SET_F32(tmpImage, oValue);
     621        oAverage += oValue;
     622        if (myReadout->bias == NULL) {
     623            myReadout->bias = psListAlloc(tmpImage);
     624        } else {
     625            psListAdd(myReadout->bias, PS_LIST_HEAD, tmpImage);
     626        }
     627    }
     628    oAverage/= (psF32) numOverscans;
     629
     630
     631    if (fit == PM_FIT_NONE) {
     632        myReadout = pmSubtractBias(myReadout, NULL, PM_FIT_NONE, overscanaxis,
     633                                   stat, nBin, NULL, NULL);
     634    } else if (fit == PM_FIT_POLYNOMIAL) {
     635        myReadout = pmSubtractBias(myReadout, myPoly, PM_FIT_POLYNOMIAL, overscanaxis,
     636                                   stat, nBin, NULL, NULL);
     637    } else if (fit == PM_FIT_SPLINE) {
     638        myReadout = pmSubtractBias(myReadout, mySpline, PM_FIT_SPLINE, overscanaxis,
     639                                   stat, nBin, NULL, NULL);
     640    }
     641
    875642    for (i=0;i<imageNumRows;i++) {
    876643        for (j=0;j<imageNumCols;j++) {
    877             myReadout->image->data.F32[i][j] = (float) (i + j);
    878         }
    879     }
    880 
    881     psStats *stat = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
    882     psPolynomial1D *myPoly = psPolynomial1DAlloc(1, PS_POLYNOMIAL_ORD);
    883     psSpline1D *mySpline = NULL;
    884 
    885 
    886     /*
    887         if (overscanNumRows <= 0) {
    888             overscanNumRows = 1;
    889         }
    890         if (overscanNumCols <= 0) {
    891             overscanNumCols = 1;
    892         }
    893     */
    894     psImage *tmpImage2 = psImageAlloc(overscanNumCols, overscanNumRows, PS_TYPE_F32);
    895     psImage *tmpImage3 = psImageAlloc(overscanNumCols, overscanNumRows, PS_TYPE_F32);
    896     psImage *tmpImage4 = psImageAlloc(overscanNumCols, overscanNumRows, PS_TYPE_F32);
    897     psList *list;
    898     for (i=0;i<overscanNumRows;i++) {
    899         for (j=0;j<overscanNumCols;j++) {
    900             tmpImage2->data.F32[i][j] = 3.0;
    901             tmpImage3->data.F32[i][j] = 4.0;
    902             tmpImage4->data.F32[i][j] = 5.0;
    903         }
    904     }
    905     list = psListAlloc(tmpImage2);
    906     psListAdd(list, PS_LIST_HEAD, tmpImage3);
    907     psListAdd(list, PS_LIST_HEAD, tmpImage4);
    908 
    909     printPositiveTestHeader(stdout, "pmSubtractBias", "Column Overscans");
    910 
    911     if (fit == PM_FIT_NONE) {
    912         myReadout = pmSubtractBias(myReadout, NULL, list, overscanaxis, stat,
    913                                    nBin, fit, NULL);
    914     } else if (fit == PM_FIT_POLYNOMIAL) {
    915         myReadout = pmSubtractBias(myReadout, myPoly, list, overscanaxis, stat,
    916                                    nBin, fit, NULL);
    917     } else if (fit == PM_FIT_SPLINE) {
    918         myReadout = pmSubtractBias(myReadout, mySpline, list, overscanaxis, stat,
    919                                    nBin, fit, NULL);
    920     }
    921 
    922     for (i=0;i<imageNumRows;i++) {
    923         for (j=0;j<imageNumCols;j++) {
    924             expect = ((float) (i + j)) - 12.0;
     644            expect = ((float) (i + j)) - oAverage;
    925645            actual = myReadout->image->data.F32[i][j];
    926646            if (FLT_EPSILON < fabs(expect - actual)) {
    927647                printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
    928648                testStatus = 1;
     649            } else {
     650                //printf("COOL: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
    929651            }
    930652        }
    931653    }
    932 
    933 
     654    psFree(fpa);
     655    psFree(chip);
     656    psFree(cell);
     657    psListElem *tmpElem = (psListElem *) myReadout->bias->head;
     658    while (NULL != tmpElem) {
     659        psFree((psImage *) tmpElem->data);
     660        tmpElem = tmpElem->next;
     661    }
    934662    psFree(myReadout);
    935     psFree(tmpImage2);
    936     psFree(tmpImage3);
    937     psFree(tmpImage4);
     663    psFree(camera);
    938664    psFree(stat);
    939     psFree(list);
    940665    psFree(myPoly);
    941666    psFree(mySpline);
     
    945670}
    946671
    947 int testX( void )
     672
     673
     674
     675
     676
     677
     678
     679
     680
     681
     682
     683
     684
     685
     686
     687
     688
     689
     690
     691
     692
     693
     694
     695
     696
     697
     698
     699
     700
     701
     702
     703
     704
     705/******************************************************************************
     706test05a() The following combinations are tested here:
     707 Overscan images are same size, no fit, bin factor is 1.
     708 Overscan images are same size, no fit, bin factor is 2.
     709 *****************************************************************************/
     710int test05a(
     711    psS32 imageNumCols,
     712    psS32 imageNumRows,
     713    psS32 overscanNumCols,
     714    psS32 overscanNumRows)
    948715{
    949716    int testStatus = 0;
     
    955722    // Overscan images are same size, no fit, bin factor is 1.
    956723    //
    957     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    958                   NUM_COLS, NUM_ROWS,
     724    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     725                  overscanNumCols, overscanNumRows, NUM_OVERSCANS,
    959726                  PM_OVERSCAN_ALL,
    960727                  PM_FIT_NONE, 1);
    961728
    962     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    963                   NUM_COLS, NUM_ROWS,
     729    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     730                  overscanNumCols, overscanNumRows, NUM_OVERSCANS,
    964731                  PM_OVERSCAN_COLUMNS,
    965732                  PM_FIT_NONE, 1);
    966733
    967     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    968                   NUM_COLS, NUM_ROWS,
     734    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     735                  overscanNumCols, overscanNumRows, NUM_OVERSCANS,
    969736                  PM_OVERSCAN_ROWS,
    970737                  PM_FIT_NONE, 1);
     
    973740    // Overscan images are same size, no fit, bin factor is 2.
    974741    //
    975     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    976                   NUM_COLS, NUM_ROWS,
     742    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     743                  overscanNumCols, overscanNumRows, NUM_OVERSCANS,
    977744                  PM_OVERSCAN_ALL,
    978745                  PM_FIT_NONE, 2);
    979746
    980     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    981                   NUM_COLS, NUM_ROWS,
     747    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     748                  overscanNumCols, overscanNumRows, NUM_OVERSCANS,
    982749                  PM_OVERSCAN_COLUMNS,
    983750                  PM_FIT_NONE, 2);
    984751
    985     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    986                   NUM_COLS, NUM_ROWS,
     752    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     753                  overscanNumCols, overscanNumRows, NUM_OVERSCANS,
    987754                  PM_OVERSCAN_ROWS,
    988755                  PM_FIT_NONE, 2);
    989     //
    990     // Overscan images are too small, spline fit, bin factor is 1.
    991     //
    992     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    993                   NUM_COLS-1, NUM_ROWS-1,
    994                   PM_OVERSCAN_ALL,
    995                   PM_FIT_SPLINE, 1);
    996 
    997     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    998                   NUM_COLS-1, NUM_ROWS-1,
    999                   PM_OVERSCAN_COLUMNS,
    1000                   PM_FIT_SPLINE, 1);
    1001 
    1002     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    1003                   NUM_COLS-1, NUM_ROWS-1,
    1004                   PM_OVERSCAN_ROWS,
    1005                   PM_FIT_SPLINE, 1);
    1006 
    1007     //
    1008     // Overscan images are too small, spline fit, bin factor is 2.
    1009     //
    1010     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    1011                   NUM_COLS-1, NUM_ROWS-1,
    1012                   PM_OVERSCAN_ALL,
    1013                   PM_FIT_SPLINE, 2);
    1014 
    1015     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    1016                   NUM_COLS-1, NUM_ROWS-1,
    1017                   PM_OVERSCAN_COLUMNS,
    1018                   PM_FIT_SPLINE, 2);
    1019 
    1020     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    1021                   NUM_COLS-1, NUM_ROWS-1,
    1022                   PM_OVERSCAN_ROWS,
    1023                   PM_FIT_SPLINE, 2);
    1024 
    1025     //
    1026     // Overscan images are same size, spline fit, bin factor is 1.
    1027     //
    1028     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    1029                   NUM_COLS, NUM_ROWS,
    1030                   PM_OVERSCAN_ALL,
    1031                   PM_FIT_SPLINE, 1);
    1032 
    1033     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    1034                   NUM_COLS, NUM_ROWS,
    1035                   PM_OVERSCAN_COLUMNS,
    1036                   PM_FIT_SPLINE, 1);
    1037 
    1038     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    1039                   NUM_COLS, NUM_ROWS,
    1040                   PM_OVERSCAN_ROWS,
    1041                   PM_FIT_SPLINE, 1);
    1042     //
    1043     // Overscan images are same size, spline fit, bin factor is 2.
    1044     //
    1045     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    1046                   NUM_COLS, NUM_ROWS,
    1047                   PM_OVERSCAN_ALL,
    1048                   PM_FIT_SPLINE, 2);
    1049 
    1050     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    1051                   NUM_COLS, NUM_ROWS,
    1052                   PM_OVERSCAN_COLUMNS,
    1053                   PM_FIT_SPLINE, 2);
    1054 
    1055     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    1056                   NUM_COLS, NUM_ROWS,
    1057                   PM_OVERSCAN_ROWS,
    1058                   PM_FIT_SPLINE, 2);
    1059 
    1060     //
    1061     // Overscan images are same size, polynomial fit, bin factor is 1.
    1062     //
    1063     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    1064                   NUM_COLS, NUM_ROWS,
    1065                   PM_OVERSCAN_ALL,
    1066                   PM_FIT_POLYNOMIAL, 1);
    1067 
    1068     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    1069                   NUM_COLS, NUM_ROWS,
    1070                   PM_OVERSCAN_COLUMNS,
    1071                   PM_FIT_POLYNOMIAL, 1);
    1072 
    1073     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    1074                   NUM_COLS, NUM_ROWS,
    1075                   PM_OVERSCAN_ROWS,
    1076                   PM_FIT_POLYNOMIAL, 1);
    1077     //
    1078     // Overscan images are same size, polynomial fit, bin factor is 2.
    1079     //
    1080     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    1081                   NUM_COLS, NUM_ROWS,
    1082                   PM_OVERSCAN_ALL,
    1083                   PM_FIT_POLYNOMIAL, 2);
    1084 
    1085     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    1086                   NUM_COLS, NUM_ROWS,
    1087                   PM_OVERSCAN_COLUMNS,
    1088                   PM_FIT_POLYNOMIAL, 2);
    1089 
    1090     testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
    1091                   NUM_COLS, NUM_ROWS,
    1092                   PM_OVERSCAN_ROWS,
    1093                   PM_FIT_POLYNOMIAL, 2);
    1094756
    1095757    return(testStatus);
    1096758}
    1097759
     760
     761/******************************************************************************
     762test05b() The following combinations are tested here:
     763 Overscan images are too small, spline fit, bin factor is 1.
     764 Overscan images are too small, spline fit, bin factor is 2.
     765 Overscan images are same size, spline fit, bin factor is 1.
     766 Overscan images are same size, spline fit, bin factor is 2.
     767 Overscan images are too big,   spline fit, bin factor is 1.
     768 Overscan images are too big,   spline fit, bin factor is 2.
     769 A single overscan image of the same size, spline fit, bin factor is 1.
     770 
     771 Overscan images are too small, polynomial fit, bin factor is 1.
     772 Overscan images are too small, polynomial fit, bin factor is 2.
     773 Overscan images are same size, polynomial fit, bin factor is 1.
     774 Overscan images are same size, polynomial fit, bin factor is 2.
     775 Overscan images are too big,   polynomial fit, bin factor is 1.
     776 Overscan images are too big,   polynomial fit, bin factor is 2.
     777 A single overscan image of the same size, polynomial fit, bin factor is 1.
     778 
     779XXX: Must add M-by-N image size tests.
     780 *****************************************************************************/
     781int test05b(
     782    psS32 imageNumCols,
     783    psS32 imageNumRows,
     784    psS32 overscanNumCols,
     785    psS32 overscanNumRows)
     786{
     787    int testStatus = 0;
     788
     789    //
     790    // Overscan images are too small, spline fit, bin factor is 1.
     791    //
     792    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     793                  overscanNumCols-2, overscanNumRows-2, NUM_OVERSCANS,
     794                  PM_OVERSCAN_ALL,
     795                  PM_FIT_SPLINE, 1);
     796
     797    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     798                  overscanNumCols-2, overscanNumRows-2, NUM_OVERSCANS,
     799                  PM_OVERSCAN_COLUMNS,
     800                  PM_FIT_SPLINE, 1);
     801
     802    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     803                  overscanNumCols-2, overscanNumRows-2, NUM_OVERSCANS,
     804                  PM_OVERSCAN_ROWS,
     805                  PM_FIT_SPLINE, 1);
     806
     807    //
     808    // Overscan images are too small, spline fit, bin factor is 2.
     809    //
     810    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     811                  overscanNumCols-2, overscanNumRows-2, NUM_OVERSCANS,
     812                  PM_OVERSCAN_ALL,
     813                  PM_FIT_SPLINE, 2);
     814
     815    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     816                  overscanNumCols-2, overscanNumRows-2, NUM_OVERSCANS,
     817                  PM_OVERSCAN_COLUMNS,
     818                  PM_FIT_SPLINE, 2);
     819
     820    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     821                  overscanNumCols-2, overscanNumRows-2, NUM_OVERSCANS,
     822                  PM_OVERSCAN_ROWS,
     823                  PM_FIT_SPLINE, 2);
     824
     825
     826    //
     827    // Overscan images are same size, spline fit, bin factor is 1.
     828    //
     829    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     830                  overscanNumCols, overscanNumRows, NUM_OVERSCANS,
     831                  PM_OVERSCAN_ALL,
     832                  PM_FIT_SPLINE, 1);
     833
     834    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     835                  overscanNumCols, overscanNumRows, NUM_OVERSCANS,
     836                  PM_OVERSCAN_COLUMNS,
     837                  PM_FIT_SPLINE, 1);
     838
     839    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     840                  overscanNumCols, overscanNumRows, NUM_OVERSCANS,
     841                  PM_OVERSCAN_ROWS,
     842                  PM_FIT_SPLINE, 1);
     843
     844    //
     845    // Overscan images are same size, spline fit, bin factor is 2.
     846    //
     847    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     848                  overscanNumCols, overscanNumRows, NUM_OVERSCANS,
     849                  PM_OVERSCAN_ALL,
     850                  PM_FIT_SPLINE, 2);
     851
     852    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     853                  overscanNumCols, overscanNumRows, NUM_OVERSCANS,
     854                  PM_OVERSCAN_COLUMNS,
     855                  PM_FIT_SPLINE, 2);
     856
     857    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     858                  overscanNumCols, overscanNumRows, NUM_OVERSCANS,
     859                  PM_OVERSCAN_ROWS,
     860                  PM_FIT_SPLINE, 2);
     861
     862    //
     863    // Overscan images are too big, spline fit, bin factor is 1.
     864    //
     865    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     866                  overscanNumCols+2, overscanNumRows+2, NUM_OVERSCANS,
     867                  PM_OVERSCAN_ALL,
     868                  PM_FIT_SPLINE, 1);
     869
     870    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     871                  overscanNumCols+2, overscanNumRows+2, NUM_OVERSCANS,
     872                  PM_OVERSCAN_COLUMNS,
     873                  PM_FIT_SPLINE, 1);
     874
     875    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     876                  overscanNumCols+2, overscanNumRows+2, NUM_OVERSCANS,
     877                  PM_OVERSCAN_ROWS,
     878                  PM_FIT_SPLINE, 1);
     879
     880    //
     881    // Overscan images are too big, spline fit, bin factor is 2.
     882    //
     883    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     884                  overscanNumCols+2, overscanNumRows+2, NUM_OVERSCANS,
     885                  PM_OVERSCAN_ALL,
     886                  PM_FIT_SPLINE, 2);
     887
     888    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     889                  overscanNumCols+2, overscanNumRows+2, NUM_OVERSCANS,
     890                  PM_OVERSCAN_COLUMNS,
     891                  PM_FIT_SPLINE, 2);
     892
     893    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     894                  overscanNumCols+2, overscanNumRows+2 , NUM_OVERSCANS,
     895                  PM_OVERSCAN_ROWS,
     896                  PM_FIT_SPLINE, 2);
     897
     898
     899    //
     900    // A single overscan image of the same size, spline fit, bin factor is 1.
     901    //
     902    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     903                  overscanNumCols, overscanNumRows, 1,
     904                  PM_OVERSCAN_ALL,
     905                  PM_FIT_SPLINE, 1);
     906
     907    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     908                  overscanNumCols, overscanNumRows, 1,
     909                  PM_OVERSCAN_COLUMNS,
     910                  PM_FIT_SPLINE, 1);
     911
     912    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     913                  overscanNumCols, overscanNumRows, 1,
     914                  PM_OVERSCAN_ROWS,
     915                  PM_FIT_SPLINE, 1);
     916
     917
     918    //
     919    // Overscan images are too small, polynomial fit, bin factor is 1.
     920    //
     921    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     922                  overscanNumCols-2, overscanNumRows-2, NUM_OVERSCANS,
     923                  PM_OVERSCAN_ALL,
     924                  PM_FIT_POLYNOMIAL, 1);
     925
     926    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     927                  overscanNumCols-2, overscanNumRows-2, NUM_OVERSCANS,
     928                  PM_OVERSCAN_COLUMNS,
     929                  PM_FIT_POLYNOMIAL, 1);
     930
     931    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     932                  overscanNumCols-2, overscanNumRows-2, NUM_OVERSCANS,
     933                  PM_OVERSCAN_ROWS,
     934                  PM_FIT_POLYNOMIAL, 1);
     935
     936    //
     937    // Overscan images are too small, polynomial fit, bin factor is 2.
     938    //
     939    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     940                  overscanNumCols-2, overscanNumRows-2, NUM_OVERSCANS,
     941                  PM_OVERSCAN_ALL,
     942                  PM_FIT_POLYNOMIAL, 2);
     943
     944    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     945                  overscanNumCols-2, overscanNumRows-2, NUM_OVERSCANS,
     946                  PM_OVERSCAN_COLUMNS,
     947                  PM_FIT_POLYNOMIAL, 2);
     948
     949    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     950                  overscanNumCols-2, overscanNumRows-2, NUM_OVERSCANS,
     951                  PM_OVERSCAN_ROWS,
     952                  PM_FIT_POLYNOMIAL, 2);
     953
     954
     955    //
     956    // Overscan images are same size, polynomial fit, bin factor is 1.
     957    //
     958    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     959                  overscanNumCols, overscanNumRows, NUM_OVERSCANS,
     960                  PM_OVERSCAN_ALL,
     961                  PM_FIT_POLYNOMIAL, 1);
     962
     963    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     964                  overscanNumCols, overscanNumRows, NUM_OVERSCANS,
     965                  PM_OVERSCAN_COLUMNS,
     966                  PM_FIT_POLYNOMIAL, 1);
     967
     968    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     969                  overscanNumCols, overscanNumRows, NUM_OVERSCANS,
     970                  PM_OVERSCAN_ROWS,
     971                  PM_FIT_POLYNOMIAL, 1);
     972
     973    //
     974    // Overscan images are same size, polynomial fit, bin factor is 2.
     975    //
     976    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     977                  overscanNumCols, overscanNumRows, NUM_OVERSCANS,
     978                  PM_OVERSCAN_ALL,
     979                  PM_FIT_POLYNOMIAL, 2);
     980
     981    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     982                  overscanNumCols, overscanNumRows, NUM_OVERSCANS,
     983                  PM_OVERSCAN_COLUMNS,
     984                  PM_FIT_POLYNOMIAL, 2);
     985
     986    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     987                  overscanNumCols, overscanNumRows, NUM_OVERSCANS,
     988                  PM_OVERSCAN_ROWS,
     989                  PM_FIT_POLYNOMIAL, 2);
     990
     991    //
     992    // Overscan images are too big, polynomial fit, bin factor is 1.
     993    //
     994    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     995                  overscanNumCols+2, overscanNumRows+2, NUM_OVERSCANS,
     996                  PM_OVERSCAN_ALL,
     997                  PM_FIT_POLYNOMIAL, 1);
     998
     999    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     1000                  overscanNumCols+2, overscanNumRows+2, NUM_OVERSCANS,
     1001                  PM_OVERSCAN_COLUMNS,
     1002                  PM_FIT_POLYNOMIAL, 1);
     1003
     1004    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     1005                  overscanNumCols+2, overscanNumRows+2, NUM_OVERSCANS,
     1006                  PM_OVERSCAN_ROWS,
     1007                  PM_FIT_POLYNOMIAL, 1);
     1008
     1009    //
     1010    // Overscan images are too big, polynomial fit, bin factor is 2.
     1011    //
     1012    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     1013                  overscanNumCols+2, overscanNumRows+2, NUM_OVERSCANS,
     1014                  PM_OVERSCAN_ALL,
     1015                  PM_FIT_POLYNOMIAL, 2);
     1016
     1017    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     1018                  overscanNumCols+2, overscanNumRows+2, NUM_OVERSCANS,
     1019                  PM_OVERSCAN_COLUMNS,
     1020                  PM_FIT_POLYNOMIAL, 2);
     1021
     1022    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     1023                  overscanNumCols+2, overscanNumRows+2, NUM_OVERSCANS,
     1024                  PM_OVERSCAN_ROWS,
     1025                  PM_FIT_POLYNOMIAL, 2);
     1026
     1027    //
     1028    // A single overscan image of the same size, polynomial fit, bin factor is 1.
     1029    //
     1030    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     1031                  overscanNumCols, overscanNumRows, 1,
     1032                  PM_OVERSCAN_ALL,
     1033                  PM_FIT_POLYNOMIAL, 1);
     1034
     1035    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     1036                  overscanNumCols, overscanNumRows, 1,
     1037                  PM_OVERSCAN_COLUMNS,
     1038                  PM_FIT_POLYNOMIAL, 1);
     1039
     1040    testStatus |= doSubtractFullOverscanColumnsGeneric(imageNumCols, imageNumRows,
     1041                  overscanNumCols, overscanNumRows, 1,
     1042                  PM_OVERSCAN_ROWS,
     1043                  PM_FIT_POLYNOMIAL, 1);
     1044
     1045    return(testStatus);
     1046}
     1047
     1048
     1049
     1050#define LOW_COLS 3
     1051#define LOW_ROWS 3
     1052
     1053/******************************************************************************
     1054test05(): See test05a() and test05b().
     1055 
     1056We run the tests in test05b() starting with all possible combinations of
     1057sizes.
     1058 
     1059XXX: Must add M-by-N image size tests.
     1060 *****************************************************************************/
     1061int test05()
     1062{
     1063    int testStatus = 0;
     1064    testStatus = test05a(NUM_COLS, NUM_ROWS, NUM_COLS, NUM_ROWS);
     1065
     1066    //    testStatus|= test05b(LOW_COLS, LOW_ROWS, LOW_COLS, LOW_ROWS);
     1067    //    testStatus|= test05b(LOW_COLS, LOW_ROWS, LOW_COLS, NUM_ROWS);
     1068    //    testStatus|= test05b(LOW_COLS, LOW_ROWS, NUM_COLS, LOW_ROWS);
     1069    //    testStatus|= test05b(LOW_COLS, LOW_ROWS, NUM_COLS, NUM_ROWS);
     1070
     1071    //    testStatus|= test05b(LOW_COLS, NUM_ROWS, LOW_COLS, LOW_ROWS);
     1072    //    testStatus|= test05b(LOW_COLS, NUM_ROWS, LOW_COLS, NUM_ROWS);
     1073    //    testStatus|= test05b(LOW_COLS, NUM_ROWS, NUM_COLS, LOW_ROWS);
     1074    //    testStatus|= test05b(LOW_COLS, NUM_ROWS, NUM_COLS, NUM_ROWS);
     1075
     1076    //    testStatus|= test05b(NUM_COLS, LOW_ROWS, LOW_COLS, LOW_ROWS);
     1077    //    testStatus|= test05b(NUM_COLS, LOW_ROWS, LOW_COLS, NUM_ROWS);
     1078    //    testStatus|= test05b(NUM_COLS, LOW_ROWS, NUM_COLS, LOW_ROWS);
     1079    //    testStatus|= test05b(NUM_COLS, LOW_ROWS, NUM_COLS, NUM_ROWS);
     1080
     1081    //    testStatus|= test05b(NUM_COLS, NUM_ROWS, LOW_COLS, LOW_ROWS);
     1082    //    testStatus|= test05b(NUM_COLS, NUM_ROWS, LOW_COLS, NUM_ROWS);
     1083    //    testStatus|= test05b(NUM_COLS, NUM_ROWS, NUM_COLS, LOW_ROWS);
     1084    testStatus|= test05b(NUM_COLS, NUM_ROWS, NUM_COLS, NUM_ROWS);
     1085
     1086
     1087    return(testStatus);
     1088}
     1089
     1090
  • trunk/psModules/test/imsubtract/tst_pmSubtractSky.c

    r5169 r5516  
    77 *  @author GLG, MHPCC
    88 *
    9  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2005-09-28 20:42:52 $
     9 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2005-11-15 20:09:03 $
     11 *
     12 *  XXX: I added the CELL.TRIMSEC region code but there are not tests for it.
    1113 *
    1214 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1315 */
    14 
    15 
    16 
    1716#include "psTest.h"
    1817#include "pslib.h"
     
    6766    myReadout = pmSubtractSky(myReadout, (void *) myPoly, PM_FIT_POLYNOMIAL,
    6867                              binFactor, myStats, 10.0);
    69     for (i=0;i<numRows;i++) {
    70         for (j=0;j<numCols;j++) {
    71             if (ERROR_TOLERANCE < fabs(myReadout->image->data.F32[i][j])) {
    72                 printf("TEST ERROR: image[%d][%d] is %f, should be 0.0\n", i, j, myReadout->image->data.F32[i][j]);
    73                 testStatus = 1;
     68    if (myReadout == NULL) {
     69        printf("TEST ERROR: pmSubtractSky() returned NULL.\n");
     70        testStatus = 1;
     71    } else {
     72        for (i=0;i<numRows;i++) {
     73            for (j=0;j<numCols;j++) {
     74                if (ERROR_TOLERANCE < fabs(myReadout->image->data.F32[i][j])) {
     75                    printf("TEST ERROR: image[%d][%d] is %f, should be 0.0\n", i, j, myReadout->image->data.F32[i][j]);
     76                    testStatus = 1;
     77                }
    7478            }
    7579        }
     
    117121    myReadout = pmSubtractSky(myReadout, (void *) myPoly, PM_FIT_POLYNOMIAL,
    118122                              binFactor, myStats, 2.0);
    119     for (i=0;i<numRows;i++) {
    120         for (j=0;j<numCols;j++) {
    121             if (errorTolerance < fabs(myReadout->image->data.F32[i][j] - trueImage->data.F32[i][j])) {
    122                 printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j,
    123                        myReadout->image->data.F32[i][j], trueImage->data.F32[i][j]);
    124                 testStatus = 1;
     123    if (myReadout == NULL) {
     124        printf("TEST ERROR: pmSubtractSky() returned NULL.\n");
     125        testStatus = 1;
     126    } else {
     127        for (i=0;i<numRows;i++) {
     128            for (j=0;j<numCols;j++) {
     129                if (errorTolerance < fabs(myReadout->image->data.F32[i][j] - trueImage->data.F32[i][j])) {
     130                    printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j,
     131                           myReadout->image->data.F32[i][j], trueImage->data.F32[i][j]);
     132                    testStatus = 1;
     133                }
    125134            }
    126135        }
  • trunk/psModules/test/objects/tst_pmObjects01.c

    r5435 r5516  
    2727    most of psObjects.c is not tested
    2828 *
    29  *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    30  *  @date $Date: 2005-10-20 23:06:24 $
     29 *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
     30 *  @date $Date: 2005-11-15 20:09:03 $
    3131 *
    3232 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    12871287
    12881288
    1289 // this code will
    1290 
    1291 
    1292 
Note: See TracChangeset for help on using the changeset viewer.