IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 17, 2006, 7:13:42 AM (20 years ago)
Author:
magnier
Message:

bulk merge of eam_rel9_p0 onto this branch

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/rel10_ifa/psModules/src/imsubtract/pmSubtractBias.c

    r6325 r6448  
     1//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     2// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the
     3// one that was being worked on.
     4//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     5
    16/** @file  pmSubtractBias.c
    27 *
    … …  
    611 *  @author GLG, MHPCC
    712 *
    8  *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2006-02-06 21:03:25 $
     13 *  @version $Revision: 1.9.4.1 $ $Name: not supported by cvs2svn $
     14 *  @date $Date: 2006-02-17 17:13:42 $
    1015 *
    1116 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1217 *
    1318 */
    14 /*****************************************************************************/
    15 /* INCLUDE FILES                                                             */
    16 /*****************************************************************************/
    17 #include <stdio.h>
    18 #include <math.h>
    19 #include <string.h>
    20 #include "pslib.h"
     19
    2120#if HAVE_CONFIG_H
    2221#include <config.h>
    2322#endif
     23
     24#include <assert.h>
    2425#include "pmSubtractBias.h"
    2526
    26 /*****************************************************************************/
    27 /* DEFINE STATEMENTS                                                         */
    28 /*****************************************************************************/
     27#define PM_SUBTRACT_BIAS_POLYNOMIAL_ORDER 2
     28#define PM_SUBTRACT_BIAS_SPLINE_ORDER 3
     29
     30
     31#define MAX(a,b) ((a) > (b) ? (a) : (b))
     32#define MIN(a,b) ((a) < (b) ? (a) : (b))
     33
     34
    2935// XXX: put these in psConstants.h
    30 void PS_POLY1D_PRINT(
    31     psPolynomial1D *poly)
     36void PS_POLY1D_PRINT(psPolynomial1D *poly)
    3237{
    3338    printf("-------------- PS_POLY1D_PRINT() --------------\n");
    … …  
    5762}\
    5863
    59 /*****************************************************************************/
    60 /* TYPE DEFINITIONS                                                          */
    61 /*****************************************************************************/
    62 
    63 /*****************************************************************************/
    64 /* GLOBAL VARIABLES                                                          */
    65 /*****************************************************************************/
    66 psS32 currentId = 0;                // XXX: remove
    67 psS32 memLeaks = 0;                 // XXX: remove
    68 //PRINT_MEMLEAKS(8); XXX
    69 /*****************************************************************************/
    70 /* FILE STATIC VARIABLES                                                     */
    71 /*****************************************************************************/
    72 
    73 /*****************************************************************************/
    74 /* FUNCTION IMPLEMENTATION - LOCAL                                           */
    75 /*****************************************************************************/
     64
     65void overscanOptionsFree(pmOverscanOptions *options)
     66{
     67    psFree(options->stat);
     68    psFree(options->poly);
     69    psFree(options->spline);
     70}
     71
     72pmOverscanOptions *pmOverscanOptionsAlloc(bool single, pmFit fitType, unsigned int order, psStats *stat)
     73{
     74    pmOverscanOptions *opts = psAlloc(sizeof(pmOverscanOptions));
     75    psMemSetDeallocator(opts, (psFreeFunc)overscanOptionsFree);
     76
     77    // Inputs
     78    opts->single = single;
     79    opts->fitType = fitType;
     80    opts->order = order;
     81    opts->stat = psMemIncrRefCounter(stat);
     82
     83    // Outputs
     84    opts->poly = NULL;
     85    opts->spline = NULL;
     86
     87    return opts;
     88}
     89
    7690
    7791/******************************************************************************
    78 psSubtractFrame(): this routine will take as input the pmReadout for the input
    79 image and a pmReadout for the bias image.  The bias image is subtracted in
    80 place from the input image.  We assume that sizes and types are checked
    81 elsewhere.
    82  
    83 XXX: Verify that the image and readout offsets are being used the right way.
    84  
    85 XXX: Ensure that it does the correct thing with image size.
     92psSubtractFrame(): this routine will take as input a readout for the input
     93image and a readout for the bias image.  The bias image is subtracted in
     94place from the input image.
    8695*****************************************************************************/
    87 static 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++) {
    94             in->image->data.F32[i][j]-=
    95                 bias->image->data.F32[i+in->row0-bias->row0][j+in->col0-bias->col0];
    96 
    97             if ((in->mask != NULL) && (bias->mask != NULL)) {
    98                 (in->mask->data.U8[i][j])|=
    99                     bias->mask->data.U8[i+in->row0-bias->row0][j+in->col0-bias->col0];
     96static bool SubtractFrame(pmReadout *in,// Input readout
     97                          const pmReadout *sub, // Readout to be subtracted from input
     98                          float scale   // Scale to apply before subtracting
     99                         )
     100{
     101    assert(in);
     102    assert(sub);
     103
     104    psImage *inImage  = in->image;      // The input image
     105    psImage *inMask   = in->mask;       // The input mask
     106    psImage *subImage = sub->image;     // The image to be subtracted
     107    psImage *subMask  = sub->mask;      // The mask for the subtraction image
     108
     109    // Offsets of the cells
     110    int x0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.X0");
     111    int y0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.Y0");
     112    int x0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.X0");
     113    int y0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.Y0");
     114
     115    if ((inImage->numCols + x0in - x0sub) > subImage->numCols) {
     116        psError(PS_ERR_UNKNOWN, true, "Image does not have enough columns for subtraction.\n");
     117        return false;
     118    }
     119    if ((inImage->numRows + y0in - y0sub) > subImage->numRows) {
     120        psError(PS_ERR_UNKNOWN, true, "Image does not have enough rows for subtraction.\n");
     121        return false;
     122    }
     123
     124    if (scale == 1.0) {
     125        for (int i = 0; i < inImage->numRows; i++) {
     126            for (int j = 0; j < inImage->numCols; j++) {
     127                inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub];
     128                if (inMask && subMask) {
     129                    inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub];
     130                }
    100131            }
    101132        }
    102     }
    103 
    104     return(in);
    105 }
    106 
    107 
    108 /******************************************************************************
    109 psSubtractDarkFrame(): this routine will take as input the pmReadout for the
    110 input image and a pmReadout for the dark image.  The dark image is scaled and
    111 subtracted in place from the input image.
    112  
    113 XXX: Verify that the image and readout offsets are being used the right way.
    114  
    115 XXX: Ensure that it does the correct thing with image size.
    116 *****************************************************************************/
    117 static 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];
     133    } else {
     134        for (int i = 0; i < inImage->numRows; i++) {
     135            for (int j = 0; j < inImage->numCols; j++) {
     136                inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub] * scale;
     137                if (inMask && subMask) {
     138                    inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub];
    132139                }
    133140            }
    134141        }
    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 
     142    }
     143
     144    return true;
     145}
     146
     147
     148#if 0
    152149/******************************************************************************
    153150ImageSubtractScalar(): subtract a scalar from the input image.
    154151 
    155 XXX: Is there a psLib function for this?
     152XXX: Use a psLib function for this.
     153 
     154XXX: This should
    156155 *****************************************************************************/
    157 static psImage *ImageSubtractScalar(
    158     psImage *image,
    159     psF32 scalar)
     156static psImage *ImageSubtractScalar(psImage *image,
     157                                    psF32 scalar)
    160158{
    161159    for (psS32 i=0;i<image->numRows;i++) {
    … …  
    166164    return(image);
    167165}
     166#endif
    168167
    169168/******************************************************************************
    … …  
    179178    psStatsOptions opt = 0;
    180179
    181     /*
    182         if (stat->options & PS_STAT_ROBUST_MODE) {
    183             if (numOptions == 0) {
    184                 opt = PS_STAT_ROBUST_MODE;
    185             }
    186             numOptions++;
    187         }
    188     */
    189180    if (stat->options & PS_STAT_ROBUST_MEDIAN) {
    190181        if (numOptions == 0) {
    … …  
    194185    }
    195186
    196     if (stat->options & PS_STAT_FITTED_MEAN) {
    197         if (numOptions == 0) {
    198             opt = PS_STAT_FITTED_MEAN;
    199         }
    200         numOptions++;
    201     }
    202 
    203187    if (stat->options & PS_STAT_CLIPPED_MEAN) {
    204188        if (numOptions == 0) {
    … …  
    222206
    223207    if (numOptions == 0) {
    224         psError(PS_ERR_UNKNOWN,true, "No allowable statistics options have been specified.\n");
     208        psError(PS_ERR_UNKNOWN,true, "No statistics options have been specified.\n");
    225209    }
    226210    if (numOptions != 1) {
    … …  
    231215}
    232216
    233 /******************************************************************************
    234 Polynomial1DCopy(): This private function copies the members of the existing
    235 psPolynomial1D "in" into the existing psPolynomial1D "out".  The previous
    236 members of the existing psPolynomial1D "out" are psFree'ed.
    237  *****************************************************************************/
    238 static psBool Polynomial1DCopy(
    239     psPolynomial1D *out,
    240     psPolynomial1D *in)
    241 {
    242     psFree(out->coeff);
    243     psFree(out->coeffErr);
    244     psFree(out->mask);
    245 
    246     out->type = in->type;
    247     out->nX = in->nX;
    248 
    249     out->coeff = (psF64 *) psAlloc((in->nX + 1) * sizeof(psF64));
    250     // XXX: use memcpy
    251     for (psS32 i = 0 ; i < (in->nX + 1) ; i++) {
    252         out->coeff[i] = in->coeff[i];
    253     }
    254 
    255     out->coeffErr = (psF64 *) psAlloc((in->nX + 1) * sizeof(psF64));
    256     // XXX: use memcpy
    257     for (psS32 i = 0 ; i < (in->nX + 1) ; i++) {
    258         out->coeffErr[i] = in->coeffErr[i];
    259     }
    260 
    261     out->mask = (psMaskType *) psAlloc((in->nX + 1) * sizeof(psMaskType));
    262     // XXX: use memcpy
    263     for (psS32 i = 0 ; i < (in->nX + 1) ; i++) {
    264         out->mask[i] = in->mask[i];
    265     }
    266 
    267     return(true);
    268 }
    269 
    270 /******************************************************************************
    271 Polynomial1DDup(): This private function duplicates and then returns the input
    272 psPolynomial1D "in".
    273  *****************************************************************************/
    274 static psPolynomial1D *Polynomial1DDup(
    275     psPolynomial1D *in)
    276 {
    277     psPolynomial1D *out = psPolynomial1DAlloc(in->type, in->nX);
    278     Polynomial1DCopy(out, in);
    279     return(out);
    280 }
    281 
    282 
    283 /******************************************************************************
    284 SplineCopy(): This private function copies the members of the existing
    285 psSpline in into the existing psSpline out.
    286  *****************************************************************************/
    287 static psBool SplineCopy(
    288     psSpline1D *out,
    289     psSpline1D *in)
    290 {
    291     PS_ASSERT_PTR_NON_NULL(out, false);
    292     PS_ASSERT_PTR_NON_NULL(in, false);
    293 
    294     for (psS32 i = 0 ; i < out->n ; i++) {
    295         psFree(out->spline[i]);
    296     }
    297     psFree(out->spline);
    298     psFree(out->knots);
    299     psFree(out->p_psDeriv2);
    300 
    301     out->n = in->n;
    302     out->spline = (psPolynomial1D **) psAlloc(in->n * sizeof(psPolynomial1D *));
    303     for (psS32 i = 0 ; i < in->n ; i++) {
    304         out->spline[i] = Polynomial1DDup(in->spline[i]);
    305     }
    306 
    307     // XXX: use psVectorCopy if they get it working.
    308     out->knots = psVectorAlloc(in->knots->n, in->knots->type.type);
    309     for (psS32 i = 0 ; i < in->knots->n ; i++) {
    310         out->knots->data.F32[i] = in->knots->data.F32[i];
    311     }
    312     /*
    313         out->knots = psVectorCopy(out->knots, in->knots, in->knots->type.type);
    314     */
    315 
    316     out->p_psDeriv2 = (psF32 *) psAlloc((in->n + 1) * sizeof(psF32));
    317     // XXX: use memcpy
    318     for (psS32 i = 0 ; i < (in->n + 1) ; i++) {
    319         out->p_psDeriv2[i] = in->p_psDeriv2[i];
    320     }
    321 
    322     return(true);
    323 }
    324 
     217
     218
     219#if 0
    325220/******************************************************************************
    326221ScaleOverscanVector(): this routine takes as input an arbitrary vector,
    … …  
    329224    PM_FIT_POLYNOMIAL: fit a polynomial to the entire input vector data.
    330225    PM_FIT_SPLINE: fit splines to the input vector data.
    331 The resulting spline or polynomial is set in the fitSpec argument.
     226XXX: Doesn't it make more sense to do polynomial interpolation on a few
     227elements of the input vector, rather than fit a polynomial to the entire
     228vector?
    332229 *****************************************************************************/
    333 static psVector *ScaleOverscanVector(
    334     psVector *overscanVector,
    335     psS32 n,
    336     void *fitSpec,
    337     pmFit fit)
     230static psVector *ScaleOverscanVector(psVector *overscanVector,
     231                                     psS32 n,
     232                                     void *fitSpec,
     233                                     pmFit fit)
    338234{
    339235    psTrace(".psModule.pmSubtracBias.ScaleOverscanVector", 4,
    340236            "---- ScaleOverscanVector() begin (%d -> %d) ----\n", overscanVector->n, n);
     237    //    PS_VECTOR_PRINT_F32(overscanVector);
    341238
    342239    if (NULL == overscanVector) {
    … …  
    351248    //
    352249    if (n == overscanVector->n) {
    353         return(psVectorCopy(newVec, overscanVector, PS_TYPE_F32));
    354     }
     250        for (psS32 i = 0 ; i < n ; i++) {
     251            newVec->data.F32[i] = overscanVector->data.F32[i];
     252        }
     253        return(newVec);
     254    }
     255    psPolynomial1D *myPoly;
     256    psSpline1D *mySpline;
    355257    psF32 x;
    356 
     258    psS32 i;
    357259    if (fit == PM_FIT_POLYNOMIAL) {
    358260        // Fit a polynomial to the old overscan vector.
    359         psPolynomial1D *myPoly = (psPolynomial1D *) fitSpec;
     261        myPoly = (psPolynomial1D *) fitSpec;
    360262        PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
    361         PS_ASSERT_POLY1D(myPoly, NULL);
    362263        myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL);
    363264        if (myPoly == NULL) {
    364             psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the psVector.\n");
     265            psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(1): Could not fit a polynomial to the psVector.\n");
    365266            return(NULL);
    366267        }
    … …  
    369270        // of the old vector, use the fitted polynomial to determine the
    370271        // interpolated value at that point, and set the new vector.
    371         for (psS32 i=0;i<n;i++) {
     272        for (i=0;i<n;i++) {
    372273            x = ((psF32) i) * ((psF32) overscanVector->n) / ((psF32) n);
    373274            newVec->data.F32[i] = psPolynomial1DEval(myPoly, x);
    374275        }
    375276    } else if (fit == PM_FIT_SPLINE) {
     277        psS32 mustFreeSpline = 0;
     278        // Fit a spline to the old overscan vector.
     279        mySpline = (psSpline1D *) fitSpec;
     280        // XXX: Does it make any sense to have a psSpline argument?
     281        if (mySpline == NULL) {
     282            mustFreeSpline = 1;
     283        }
     284
    376285        //
    377286        // NOTE: Since the X arg in the psVectorFitSpline1D() function is NULL,
    … …  
    379288        // properly when doing the spline eval.
    380289        //
    381         psSpline1D *mySpline = psVectorFitSpline1D(NULL, overscanVector);
     290        //        mySpline = psVectorFitSpline1D(mySpline, NULL, overscanVector, NULL);
     291        mySpline = psVectorFitSpline1D(NULL, overscanVector);
    382292        if (mySpline == NULL) {
    383             psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to the psVector.\n");
     293            psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(2): Could not fit a spline to the psVector.\n");
    384294            return(NULL);
    385295        }
     296        //        PS_PRINT_SPLINE(mySpline);
    386297
    387298        // For each element of the new vector, convert the x-ordinate to that
    388         // of the old vector, use the fitted spline to determine the
     299        // of the old vector, use the fitted polynomial to determine the
    389300        // interpolated value at that point, and set the new vector.
    390         for (psS32 i=0;i<n;i++) {
     301        for (i=0;i<n;i++) {
    391302            // Scale to [0 : overscanVector->n - 1]
    392303            x = ((psF32) i) * ((psF32) (overscanVector->n-1)) / ((psF32) n);
    393304            newVec->data.F32[i] = psSpline1DEval(mySpline, x);
    394305        }
    395 
    396         psSpline1D *ptrSpline = (psSpline1D *) fitSpec;
    397         if (ptrSpline != NULL) {
    398             // Copy the resulting spline fit into ptrSpline.
    399             PS_ASSERT_SPLINE(ptrSpline, NULL);
    400             SplineCopy(ptrSpline, mySpline);
    401         }
    402         psFree(mySpline);
     306        if (mustFreeSpline ==1) {
     307            psFree(mySpline);
     308        }
     309        //        PS_VECTOR_PRINT_F32(newVec);
     310
     311
    403312    } else {
    404313        psError(PS_ERR_UNKNOWN, true, "unknown fit type.  Returning NULL.\n");
    … …  
    412321}
    413322
     323#endif
     324
     325// Produce an overscan vector from an array of pixels
     326static psVector *overscanVector(pmOverscanOptions *overscanOpts, // Overscan options
     327                                const psArray *pixels, // Array of vectors containing the pixel values
     328                                psStats *myStats // Statistic to use in reducing the overscan
     329                               )
     330{
     331    // Reduce the overscans
     332    psVector *reduced = psVectorAlloc(pixels->n, PS_TYPE_F32); // Overscan for each row
     333    psVector *ordinate = psVectorAlloc(pixels->n, PS_TYPE_F32); // Ordinate
     334    psVector *mask = psVectorAlloc(pixels->n, PS_TYPE_U8); // Mask for fitting
     335    for (int i = 0; i < pixels->n; i++) {
     336        psVector *values = pixels->data[i]; // Vector with overscan values
     337        if (values->n > 0) {
     338            mask->data.U8[i] = 0;
     339            ordinate->data.F32[i] = 2.0*(float)i/(float)pixels->n - 1.0; // Scale to [-1,1]
     340            psVectorStats(myStats, values, NULL, NULL, 0);
     341            double reducedVal = NAN; // Result of statistics
     342            if (! p_psGetStatValue(myStats, &reducedVal)) {
     343                psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result "
     344                        "of statistics on row %d.\n", i);
     345                return NULL;
     346            }
     347            reduced->data.F32[i] = reducedVal;
     348        } else if (overscanOpts->fitType == PM_FIT_NONE) {
     349            psError(PS_ERR_UNKNOWN, true, "The overscan is not supplied for all points on the "
     350                    "image, and no fit is requested.\n");
     351            return NULL;
     352        } else {
     353            // We'll fit this one out
     354            mask->data.U8[i] = 1;
     355        }
     356    }
     357
     358    // Fit the overscan, if required
     359    switch (overscanOpts->fitType) {
     360    case PM_FIT_NONE:
     361        // No fitting --- that's easy.
     362        break;
     363    case PM_FIT_POLY_ORD:
     364        if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order ||
     365                                   overscanOpts->poly->type != PS_POLYNOMIAL_ORD)) {
     366            psFree(overscanOpts->poly);
     367            overscanOpts->poly = NULL;
     368        }
     369        if (! overscanOpts->poly) {
     370            overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, overscanOpts->order);
     371        }
     372        psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate);
     373        psFree(reduced);
     374        reduced = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
     375        break;
     376    case PM_FIT_POLY_CHEBY:
     377        if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order ||
     378                                   overscanOpts->poly->type != PS_POLYNOMIAL_CHEB)) {
     379            psFree(overscanOpts->poly);
     380            overscanOpts->poly = NULL;
     381        }
     382        if (! overscanOpts->poly) {
     383            overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_CHEB, overscanOpts->order);
     384        }
     385        psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate);
     386        psFree(reduced);
     387        reduced = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
     388        break;
     389    case PM_FIT_SPLINE:
     390        // XXX I don't think psSpline1D is up to scratch yet --- it has no mask, and requires an
     391        // input spline
     392        overscanOpts->spline = psVectorFitSpline1D(reduced, ordinate);
     393        psFree(reduced);
     394        reduced = psSpline1DEvalVector(overscanOpts->spline, ordinate);
     395        break;
     396    default:
     397        psError(PS_ERR_UNKNOWN, true, "Unknown value for the fitting type: %d\n", overscanOpts->fitType);
     398        return NULL;
     399        break;
     400    }
     401
     402    psFree(ordinate);
     403    psFree(mask);
     404
     405    return reduced;
     406}
     407
     408
     409
    414410/******************************************************************************
     411XXX: The SDRS does not specify type support.  F32 is implemented here.
    415412 *****************************************************************************/
    416 static psS32 GetOverscanSize(
    417     psImage *inImg,
    418     pmOverscanAxis overScanAxis)
    419 {
    420     if (overScanAxis == PM_OVERSCAN_ROWS) {
    421         return(inImg->numCols);
    422     } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
    423         return(inImg->numRows);
    424     } else if (overScanAxis == PM_OVERSCAN_ALL) {
    425         return(1);
    426     }
    427     return(0);
    428 }
    429 
    430 /******************************************************************************
    431 GetOverscanAxis(in) this private routine determines the appropiate overscan
    432 axis from the parent cell metadata.
    433  
    434 XXX: Verify the READDIR corresponds with my overscan axis.
    435  *****************************************************************************/
    436 static pmOverscanAxis GetOverscanAxis(pmReadout *in)
    437 {
    438     psBool rc;
    439     if ((in->parent != NULL) && (in->parent->concepts)) {
    440         psS32 dir = psMetadataLookupS32(&rc, in->parent->concepts, "CELL.READDIR");
    441         if (rc == true) {
    442             if (dir == 1) {
    443                 return(PM_OVERSCAN_ROWS);
    444             } else if (dir == 2) {
    445                 return(PM_OVERSCAN_COLUMNS);
    446             } else if (dir == 3) {
    447                 return(PM_OVERSCAN_ALL);
    448             }
    449         }
    450     }
    451 
    452     psLogMsg(__func__, PS_LOG_WARN,
    453              "WARNING: pmSubtractBias.(): could not determine CELL.READDIR from in->parent metadata.  Setting overscan axis to PM_OVERSCAN_NONE.\n");
    454     return(PM_OVERSCAN_NONE);
    455 }
    456 
    457 /******************************************************************************
    458 psListLength(list): determine the length of a psList.
    459  
    460 XXX: Put this elsewhere.
    461  *****************************************************************************/
    462 static psS32 psListLength(
    463     psList *list)
    464 {
    465     psS32 length = 0;
    466     psListElem *tmpElem = (psListElem *) list->head;
    467     while (NULL != tmpElem) {
    468         tmpElem = tmpElem->next;
    469         length++;
    470     }
    471     return(length);
    472 }
    473 
    474 /******************************************************************************
    475 Note: this isn't needed anymore as of psModule SDRS 12-09.
    476  *****************************************************************************/
    477 static psBool OverscanReducePixel(
    478     psImage *in,
    479     psList *bias,
    480     psStats *myStats)
    481 {
    482     PS_ASSERT_PTR_NON_NULL(in, NULL);
    483     PS_ASSERT_PTR_NON_NULL(bias, NULL);
    484     PS_ASSERT_PTR_NON_NULL(bias->head, NULL);
    485     PS_ASSERT_PTR_NON_NULL(myStats, NULL);
    486 
    487     // Allocate a psVector with one element per overscan image.
    488     psS32 numOverscanImages = psListLength(bias);
    489     psVector *statsAll = psVectorAlloc(numOverscanImages, PS_TYPE_F32);
    490     psListElem *tmpOverscan = (psListElem *) bias->head;
    491     psS32 i = 0;
    492     psF64 statValue;
    493     //
    494     // We loop through each overscan image, calculating the specified
    495     // statistic on that image.
    496     //
    497     while (NULL != tmpOverscan) {
    498         psImage *myOverscanImage = (psImage *) tmpOverscan->data;
    499 
    500         PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL);
    501         myStats = psImageStats(myStats, myOverscanImage, NULL, (psMaskType)0xffffffff);
    502         if (myStats == NULL) {
    503             psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation.  Returning in image.\n");
    504             psFree(statsAll);
    505             return(false);
    506         }
    507         if (false == p_psGetStatValue(myStats, &statValue)) {
    508             psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
    509             psFree(statsAll);
    510             return(false);
    511         }
    512         statsAll->data.F32[i] = statValue;
    513         i++;
    514         tmpOverscan = tmpOverscan->next;
    515     }
    516 
    517     //
    518     // We reduce the individual stats for each overscan image to
    519     // a single psF32.
    520     //
    521     myStats = psVectorStats(myStats, statsAll, NULL, NULL, 0);
    522     if (myStats == NULL) {
    523         psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation.  Returning in image.\n");
    524         psFree(statsAll);
    525         return(false);
    526     }
    527     if (false == p_psGetStatValue(myStats, &statValue)) {
    528         psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
    529         psFree(statsAll);
    530         return(false);
    531     }
    532 
    533     //
    534     // Subtract the result and return.
    535     //
    536     ImageSubtractScalar(in, statValue);
    537     psFree(statsAll);
    538     return(in);
    539 }
    540 
    541 /******************************************************************************
    542 ReduceOverscanImageToCol(overscanImage, myStats): This private routine reduces
    543 a single psImage to a column by combining all pixels from each row into a
    544 single pixel via requested statistic in myStats.
    545  *****************************************************************************/
    546 static psVector *ReduceOverscanImageToCol(
    547     psImage *overscanImage,
    548     psStats *myStats)
    549 {
    550     psF64 statValue;
    551     psVector *tmpRow = psVectorAlloc(overscanImage->numCols, PS_TYPE_F32);
    552     psVector *tmpCol = psVectorAlloc(overscanImage->numRows, PS_TYPE_F32);
    553 
    554     //
    555     // For each row, we store all pixels in that row into a temporary psVector,
    556     // then we run psVectorStats() on that vector.
    557     //
    558     for (psS32 i=0;i<overscanImage->numRows;i++) {
    559         for (psS32 j=0;j<overscanImage->numCols;j++) {
    560             tmpRow->data.F32[j] = overscanImage->data.F32[i][j];
    561         }
    562 
    563         psStats *rc = psVectorStats(myStats, tmpRow, NULL, NULL, 0);
    564         if (rc == NULL) {
    565             psError(PS_ERR_UNKNOWN, true, "psVectorStats() could not perform requested statistical operation.  Returning in image.\n");
    566             return(NULL);
    567         }
    568 
    569         if (false ==  p_psGetStatValue(rc, &statValue)) {
    570             psError(PS_ERR_UNKNOWN, true, "p_psGetStatValue() could not determine result from requested statistical operation.  Returning in image.\n");
    571             return(NULL);
    572         }
    573 
    574         tmpCol->data.F32[i] = (psF32) statValue;
    575     }
    576     psFree(tmpRow);
    577 
    578     return(tmpCol);
    579 }
    580 
    581 /******************************************************************************
    582 ReduceOverscanImageToCol(overscanImage, myStats): This private routine reduces
    583 a single psImage to a row by combining all pixels from each column into a
    584 single pixel via requested statistic in myStats.
    585  *****************************************************************************/
    586 static psVector *ReduceOverscanImageToRow(
    587     psImage *overscanImage,
    588     psStats *myStats)
    589 {
    590     psF64 statValue;
    591     psVector *tmpRow = psVectorAlloc(overscanImage->numCols, PS_TYPE_F32);
    592     psVector *tmpCol = psVectorAlloc(overscanImage->numRows, PS_TYPE_F32);
    593 
    594     //
    595     // For each column, we store all pixels in that column into a temporary psVector,
    596     // then we run psVectorStats() on that vector.
    597     //
    598     for (psS32 i=0;i<overscanImage->numCols;i++) {
    599         for (psS32 j=0;j<overscanImage->numRows;j++) {
    600             tmpCol->data.F32[j] = overscanImage->data.F32[j][i];
    601         }
    602 
    603         psStats *rc = psVectorStats(myStats, tmpCol, NULL, NULL, 0);
    604         if (rc == NULL) {
    605             psError(PS_ERR_UNKNOWN, true, "psVectorStats() could not perform requested statistical operation.  Returning in image.\n");
    606             return(NULL);
    607         }
    608 
    609         if (false ==  p_psGetStatValue(rc, &statValue)) {
    610             psError(PS_ERR_UNKNOWN, true, "p_psGetStatValue() could not determine result from requested statistical operation.  Returning in image.\n");
    611             return(NULL);
    612         }
    613 
    614         tmpRow->data.F32[i] = (psF32) statValue;
    615     }
    616     psFree(tmpCol);
    617 
    618     return(tmpRow);
    619 }
    620 
    621 /******************************************************************************
    622 OverscanReduce(vecSize, bias, myStats): This private routine takes a psList of
    623 overscan images (in bias) and reduces them to a single psVector via the
    624 specified psStats struct.  The vector is then scaled to the length or the
    625 row/column in inImg.
    626  *****************************************************************************/
    627 static psVector* OverscanReduce(
    628     psImage *inImg,
    629     pmOverscanAxis overScanAxis,
    630     psList *bias,
    631     void *fitSpec,
    632     pmFit fit,
    633     psStats *myStats)
    634 {
    635     if ((overScanAxis != PM_OVERSCAN_ROWS) && (overScanAxis != PM_OVERSCAN_COLUMNS)) {
    636         psError(PS_ERR_UNKNOWN, true, "overScanAxis must be PM_OVERSCAN_ROWS or PM_OVERSCAN_COLUMNS\n");
    637         return(NULL);
    638     }
    639     PS_ASSERT_PTR_NON_NULL(inImg, NULL);
    640     PS_ASSERT_PTR_NON_NULL(bias, NULL);
    641     PS_ASSERT_PTR_NON_NULL(bias->head, NULL);
    642     PS_ASSERT_PTR_NON_NULL(myStats, NULL);
    643     //
    644     // Allocate a psVector for the output of this routine.
    645     //
    646     psS32 vecSize = GetOverscanSize(inImg, overScanAxis);
    647     psVector *overscanVector = psVectorAlloc(vecSize, PS_TYPE_F32);
    648 
    649     //
    650     // Allocate an array of psVectors with one psVector per element of the
    651     // final oversan column vector.  These psVectors will be used with
    652     // psStats to reduce the multiple elements from each overscan column
    653     // vector to a single final column vector.
    654     //
    655     psS32 numOverscanImages = psListLength(bias);
    656     psVector **overscanVectors = (psVector **) psAlloc(numOverscanImages * sizeof(psVector *));
    657     for (psS32 i = 0 ; i < numOverscanImages ; i++) {
    658         overscanVectors[i] = NULL;
    659     }
    660 
    661     //
    662     // We iterate through the list of overscan images.  For each image,
    663     // we reduce it to a single column or row.  Save the overscan vector
    664     // in overscanVectors[].
    665     //
    666     psListElem *tmpOverscan = (psListElem *) bias->head;
    667     psS32 overscanID = 0;
    668     while (tmpOverscan != NULL) {
    669         psImage *tmpOverscanImage = (psImage *) tmpOverscan->data;
    670         if (overScanAxis == PM_OVERSCAN_ROWS) {
    671             overscanVectors[overscanID] = ReduceOverscanImageToRow(tmpOverscanImage, myStats);
    672         } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
    673             overscanVectors[overscanID] = ReduceOverscanImageToCol(tmpOverscanImage, myStats);
    674         }
    675 
    676         tmpOverscan = tmpOverscan->next;
    677         overscanID++;
    678     }
    679 
    680     //
    681     // For each overscan vector, if necessary, we scale that column or
    682     // row to vecSize.  Note: we should have already ensured that the
    683     // fit is poly or spline.
    684     //
    685     for (psS32 i = 0 ; i < numOverscanImages ; i++) {
    686         psVector *tmpOverscanVector = overscanVectors[i];
    687 
    688         if (tmpOverscanVector->n != vecSize) {
    689             overscanVectors[i] = ScaleOverscanVector(tmpOverscanVector, vecSize, fitSpec, fit);
    690             if (overscanVectors[i] == NULL) {
    691                 psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector.\n");
    692                 for (psS32 i = 0 ; i < numOverscanImages ; i++) {
    693                     psFree(overscanVectors[i]);
    694                 }
    695                 psFree(overscanVectors);
    696                 psFree(tmpOverscanVector);
    697                 return(NULL);
    698             }
    699             psFree(tmpOverscanVector);
    700         }
    701     }
    702 
    703     //
    704     // We collect all elements in the overscan vectors for the various
    705     // overscan images into a single psVector (tmpVec).  Then we call
    706     // psStats on that vector to determine the final values for the
    707     // overscan vector.
    708     //
    709     psVector *tmpVec = psVectorAlloc(numOverscanImages, PS_TYPE_F32);
    710     psF64 statValue;
    711     for (psS32 i = 0 ; i < vecSize ; i++) {
    712         // Collect the i-th elements from each overscan vector into a single vector.
    713         for (psS32 j = 0 ; j < numOverscanImages ; j++) {
    714             tmpVec->data.F32[j] = overscanVectors[j]->data.F32[i];
    715         }
    716 
    717         if (NULL == psVectorStats(myStats, tmpVec, NULL, NULL, 0)) {
    718             psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
    719             for (psS32 i = 0 ; i < numOverscanImages ; i++) {
    720                 psFree(overscanVectors[i]);
    721             }
    722             psFree(overscanVectors);
    723             psFree(tmpVec);
    724             return(NULL);
    725         }
    726         if (false == p_psGetStatValue(myStats, &statValue)) {
    727             psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
    728             for (psS32 i = 0 ; i < numOverscanImages ; i++) {
    729                 psFree(overscanVectors[i]);
    730             }
    731             psFree(overscanVectors);
    732             psFree(tmpVec);
    733             return(NULL);
    734         }
    735 
    736         overscanVector->data.F32[i] = (psF32) statValue;
    737     }
    738 
    739     //
    740     // We're done.  Free the intermediate overscan vectors.
    741     //
    742     psFree(tmpVec);
    743     for (psS32 i = 0 ; i < numOverscanImages ; i++) {
    744         psFree(overscanVectors[i]);
    745     }
    746     psFree(overscanVectors);
    747 
    748     //
    749     // Return the computed overscanVector
    750     //
    751     return(overscanVector);
    752 }
    753 
    754 /******************************************************************************
    755 RebinOverscanVector(overscanVector, nBinOrig, myStats): this private routine
    756 takes groups of nBinOrig elements in the input vector, combines them into a
    757 single pixel via myStats and psVectorStats(), and then outputs a vector of
    758 those pixels.
    759  *****************************************************************************/
    760 static psS32 RebinOverscanVector(
    761     psVector *overscanVector,
    762     psS32 nBinOrig,
    763     psStats *myStats)
    764 {
    765     psF64 statValue;
    766     psS32 nBin;
    767     if ((nBinOrig > 1) && (nBinOrig < overscanVector->n)) {
    768         psS32 numBins = 1+((overscanVector->n)/nBinOrig);
    769         psVector *myBin = psVectorAlloc(numBins, PS_TYPE_F32);
    770         psVector *binVec = psVectorAlloc(nBinOrig, PS_TYPE_F32);
    771 
    772         for (psS32 i=0;i<numBins;i++) {
    773             for(psS32 j=0;j<nBinOrig;j++) {
    774                 if (overscanVector->n > ((i*nBinOrig)+j)) {
    775                     binVec->data.F32[j] = overscanVector->data.F32[(i*nBinOrig)+j];
    776                 } else {
    777                     // XXX: we get here if nBinOrig does not evenly divide
    778                     // the overscanVector vector.  This is the last bin.  Should
    779                     // we change the binVec->n to acknowledge that?
    780                     binVec->n = j;
    781                 }
    782             }
    783             psStats *rc = psVectorStats(myStats, binVec, NULL, NULL, 0);
    784             if (rc == NULL) {
    785                 psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
    786                 return(-1);
    787             }
    788             if (false ==  p_psGetStatValue(rc, &statValue)) {
    789                 psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
    790                 return(-1);
    791             }
    792             myBin->data.F32[i] = statValue;
    793         }
    794 
    795         // Change the effective size of overscanVector.
    796         overscanVector->n = numBins;
    797         for (psS32 i=0;i<numBins;i++) {
    798             overscanVector->data.F32[i] = myBin->data.F32[i];
    799         }
    800         psFree(binVec);
    801         psFree(myBin);
    802         nBin = nBinOrig;
    803     } else {
    804         nBin = 1;
    805     }
    806 
    807     return(nBin);
    808 }
    809 
    810 /******************************************************************************
    811 FitOverscanVectorAndUnbin(inImg, overscanVector, overScanAxis, fitSpec, fit,
    812 nBin):  this private routine fits a psPolynomial or psSpline to the overscan
    813 vector.  It then creates a new vector, with a size determined by the input
    814 image, evaluates the psPolynomial or psSpline at each element in that vector,
    815 then returns that vector.
    816  *****************************************************************************/
    817 static psVector *FitOverscanVectorAndUnbin(
    818     psImage *inImg,
    819     psVector *overscanVector,
    820     pmOverscanAxis overScanAxis,
    821     void *fitSpec,
    822     pmFit fit,
    823     psS32 nBin)
    824 {
    825     psPolynomial1D* myPoly = NULL;
    826     psSpline1D *mySpline = NULL;
    827     //
    828     // Fit a polynomial or spline to the overscan vector.
    829     //
    830     if (fit == PM_FIT_POLYNOMIAL) {
    831         myPoly = (psPolynomial1D *) fitSpec;
    832         PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
    833         PS_ASSERT_POLY1D(myPoly, NULL);
    834         myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL);
    835         if (myPoly == NULL) {
    836             psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to overscan vector.  Returning NULL.\n");
    837             return(NULL);
    838         }
    839     } else if (fit == PM_FIT_SPLINE) {
    840         mySpline = psVectorFitSpline1D(NULL, overscanVector);
    841         if (mySpline == NULL) {
    842             psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to overscan vector.  Returning NULL.\n");
    843             return(NULL);
    844         }
    845         if (fitSpec != NULL) {
    846             // Copy the resulting spline fit into fitSpec.
    847             psSpline1D *ptrSpline = (psSpline1D *) fitSpec;
    848             PS_ASSERT_SPLINE(ptrSpline, NULL);
    849             SplineCopy(ptrSpline, mySpline);
    850         }
    851     }
    852 
    853     //
    854     // Evaluate the poly/spline at each pixel in the overscan row/column.
    855     //
    856     psS32 vecSize = GetOverscanSize(inImg, overScanAxis);
    857     psVector *newVec = psVectorAlloc(vecSize, PS_TYPE_F32);
    858     if ((nBin > 1) && (nBin < overscanVector->n)) {
    859         for (psS32 i = 0 ; i < vecSize ; i++) {
    860             if (fit == PM_FIT_POLYNOMIAL) {
    861                 newVec->data.F32[i] = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin));
    862             } else if (fit == PM_FIT_SPLINE) {
    863                 newVec->data.F32[i] = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin));
    864             }
    865         }
    866     } else {
    867         for (psS32 i = 0 ; i < vecSize ; i++) {
    868             if (fit == PM_FIT_POLYNOMIAL) {
    869                 newVec->data.F32[i] = psPolynomial1DEval(myPoly, (psF32) i);
    870             } else if (fit == PM_FIT_SPLINE) {
    871                 newVec->data.F32[i] = psSpline1DEval(mySpline, (psF32) i);
    872             }
    873         }
    874     }
    875 
    876     psFree(mySpline);
    877     psFree(overscanVector);
    878     return(newVec);
    879 }
    880 
    881 
    882 
    883 /******************************************************************************
    884 UnbinOverscanVector(inImg, overscanVector, overScanAxis, nBin):  this private
    885 routine takes a psVector overscanVector that was previously binned by a factor
    886 of nBin, and then expands it to its original size, duplicated elements nBin
    887 times for each element in the input vector overscanVector.
    888  *****************************************************************************/
    889 static psVector *UnbinOverscanVector(
    890     psImage *inImg,
    891     psVector *overscanVector,
    892     pmOverscanAxis overScanAxis,
    893     psS32 nBin)
    894 {
    895     psS32 vecSize = 0;
    896 
    897     if (overScanAxis == PM_OVERSCAN_ROWS) {
    898         vecSize = inImg->numCols;
    899     } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
    900         vecSize = inImg->numRows;
    901     }
    902 
    903     psVector *newVec = psVectorAlloc(vecSize, PS_TYPE_F32);
    904     for (psS32 i = 0 ; i < vecSize ; i++) {
    905         newVec->data.F32[i] = overscanVector->data.F32[i/nBin];
    906     }
    907 
    908     psFree(overscanVector);
    909     return(newVec);
    910 }
    911 
    912 
    913 /******************************************************************************
    914 SubtractVectorFromImage(inImg, overscanVector, overScanAxis):  this private
    915 routine subtracts the overscanVector column-wise or row-wise from inImg.
    916  *****************************************************************************/
    917 static psImage *SubtractVectorFromImage(
    918     psImage *inImg,
    919     psVector *overscanVector,
    920     pmOverscanAxis overScanAxis)
    921 {
    922     //
    923     // Subtract overscan vector row-wise from the image.
    924     //
    925     if (overScanAxis == PM_OVERSCAN_ROWS) {
    926         for (psS32 i=0;i<inImg->numCols;i++) {
    927             for (psS32 j=0;j<inImg->numRows;j++) {
    928                 inImg->data.F32[j][i]-= overscanVector->data.F32[i];
    929             }
    930         }
    931     }
    932 
    933     //
    934     // Subtract overscan vector column-wise from the image.
    935     //
    936     if (overScanAxis == PM_OVERSCAN_COLUMNS) {
    937         for (psS32 i=0;i<inImg->numRows;i++) {
    938             for (psS32 j=0;j<inImg->numCols;j++) {
    939                 inImg->data.F32[i][j]-= overscanVector->data.F32[i];
    940             }
    941         }
    942     }
    943 
    944     return(inImg);
    945 }
    946 
    947 
    948 
    949 typedef enum {
    950     PM_ERROR_NO_SUBTRACTION,
    951     PM_WARNING_NO_SUBTRACTION,
    952     PM_ERROR_NO_BIAS_SUBTRACT,
    953     PM_WARNING_NO_BIAS_SUBTRACT,
    954     PM_ERROR_NO_DARK_SUBTRACT,
    955     PM_WARNING_NO_DARK_SUBTRACT,
    956     PM_OKAY
    957 } pmSubtractBiasAssertStatus;
    958 /******************************************************************************
    959 AssertCodeOverscan(....) this private routine verifies that the various input
    960 parameters to pmSubtractBias() are correct for overscan subtraction.
    961  *****************************************************************************/
    962 pmSubtractBiasAssertStatus AssertCodeOverscan(
    963     pmReadout *in,
    964     void *fitSpec,
    965     pmFit fit,
    966     bool overscan,
    967     psStats *stat,
    968     int nBinOrig,
    969     const pmReadout *bias,
    970     const pmReadout *dark)
    971 {
    972 
    973     PS_ASSERT_READOUT_NON_NULL(in, PM_ERROR_NO_SUBTRACTION);
    974     PS_ASSERT_READOUT_NON_EMPTY(in, PM_ERROR_NO_SUBTRACTION);
    975     PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, PM_ERROR_NO_SUBTRACTION);
    976     PS_WARN_PTR_NON_NULL(in->parent);
    977     if (in->parent != NULL) {
    978         PS_WARN_PTR_NON_NULL(in->parent->concepts);
    979     }
    980 
    981     if (overscan == true) {
    982         pmOverscanAxis overScanAxis = GetOverscanAxis(in);
    983         PS_ASSERT_PTR_NON_NULL(stat, PM_ERROR_NO_SUBTRACTION);
    984         PS_ASSERT_PTR_NON_NULL(in->bias, PM_ERROR_NO_SUBTRACTION);
    985         PS_ASSERT_PTR_NON_NULL(in->bias->head, PM_ERROR_NO_SUBTRACTION);
    986         //
    987         // Check the type, size of each bias image.
    988         //
    989         psListElem *tmpOverscan = (psListElem *) in->bias->head;
    990         psS32 numOverscans = 0;
    991         while (NULL != tmpOverscan) {
    992             numOverscans++;
    993             psImage *myOverscanImage = (psImage *) tmpOverscan->data;
    994             PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, PM_ERROR_NO_SUBTRACTION);
    995             // XXX: Get this right with the rows and columns.
    996             if (overScanAxis == PM_OVERSCAN_ROWS) {
    997                 if (myOverscanImage->numRows != in->image->numRows) {
    998                     psLogMsg(__func__, PS_LOG_WARN,
    999                              "WARNING: pmSubtractBias.(): overscan image (# %d) has %d rows, input image has %d rows\n",
    1000                              numOverscans, myOverscanImage->numCols, in->image->numRows);
    1001                     if (fit == PM_FIT_NONE) {
    1002                         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");
    1003                         return(PM_ERROR_NO_SUBTRACTION);
    1004                     }
    1005                 }
    1006             } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
    1007                 if (myOverscanImage->numCols != in->image->numCols) {
    1008                     psLogMsg(__func__, PS_LOG_WARN,
    1009                              "WARNING: pmSubtractBias.(): overscan image (# %d) has %d columns, input image has %d columns\n",
    1010                              numOverscans, myOverscanImage->numCols, in->image->numCols);
    1011                     if (fit == PM_FIT_NONE) {
    1012                         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");
    1013                         return(PM_ERROR_NO_SUBTRACTION);
    1014                     }
    1015                 }
    1016             } else if (overScanAxis != PM_OVERSCAN_ALL) {
    1017                 psError(PS_ERR_UNKNOWN, true, "Must specify and overscan axis.\n");
    1018                 return(PM_ERROR_NO_SUBTRACTION);
    1019             }
    1020             tmpOverscan = tmpOverscan->next;
    1021         }
    1022     } else {
    1023         if (fit != PM_FIT_NONE) {
    1024             psLogMsg(__func__, PS_LOG_WARN,
    1025                      "WARNING: pmSubtractBias.(): overscan is FALSE and fit is not PM_FIT_NONE.\n");
    1026             return(PM_WARNING_NO_SUBTRACTION);
    1027         }
    1028     }
    1029 
    1030     // XXX: I do not like the following spec since it's useless to specify
    1031     // a psSpline as the fitSpec.
    1032     if (0) {
    1033         if ((fitSpec == NULL) &&
    1034                 ((fit != PM_FIT_NONE) || (overscan == true))) {
    1035             psError(PS_ERR_UNKNOWN, true, "fitSpec is NULL and fit is not PM_FIT_NONE or overscan is TRUE.\n");
    1036             return(PM_ERROR_NO_SUBTRACTION);
    1037         }
    1038     }
    1039 
    1040     return(PM_OKAY);
    1041 }
    1042 
    1043 /******************************************************************************
    1044 AssertCodeBias(....) this private routine verifies that the various input
    1045 parameters to pmSubtractBias() are correct for bias subtraction.
    1046  *****************************************************************************/
    1047 static pmSubtractBiasAssertStatus AssertCodeBias(
    1048     pmReadout *in,
    1049     void *fitSpec,
    1050     pmFit fit,
    1051     bool overscan,
    1052     psStats *stat,
    1053     int nBinOrig,
    1054     const pmReadout *bias,
    1055     const pmReadout *dark)
    1056 {
    1057     if ((in->image->numRows + in->row0 - bias->row0) > bias->image->numRows) {
    1058         psError(PS_ERR_UNKNOWN,true, "bias image does not have enough rows.  Returning in image\n");
    1059         return(PM_ERROR_NO_BIAS_SUBTRACT);
    1060     }
    1061     if ((in->image->numCols + in->col0 - bias->col0) > bias->image->numCols) {
    1062         psError(PS_ERR_UNKNOWN,true, "bias image does not have enough columns.  Returning in image\n");
    1063         return(PM_ERROR_NO_BIAS_SUBTRACT);
    1064     }
    1065 
    1066     if (bias != NULL) {
    1067         PS_ASSERT_READOUT_NON_EMPTY(bias, PM_ERROR_NO_BIAS_SUBTRACT);
    1068         PS_ASSERT_READOUT_TYPE(bias, PS_TYPE_F32, PM_ERROR_NO_DARK_SUBTRACT);
    1069     }
    1070     return(PM_OKAY);
    1071 }
    1072 
    1073 /******************************************************************************
    1074 AssertCodeDark(....) this private routine verifies that the various input
    1075 parameters to pmSubtractBias() are correct for dark subtraction.
    1076  *****************************************************************************/
    1077 pmSubtractBiasAssertStatus AssertCodeDark(
    1078     pmReadout *in,
    1079     void *fitSpec,
    1080     pmFit fit,
    1081     bool overscan,
    1082     psStats *stat,
    1083     int nBinOrig,
    1084     const pmReadout *bias,
    1085     const pmReadout *dark)
    1086 {
    1087     if ((in->image->numRows + in->row0 - dark->row0) > dark->image->numRows) {
    1088         psError(PS_ERR_UNKNOWN, true, "dark image does not have enough rows.  Returning in image\n");
    1089         return(PM_ERROR_NO_DARK_SUBTRACT);
    1090     }
    1091     if ((in->image->numCols + in->col0 - dark->col0) > dark->image->numCols) {
    1092         psError(PS_ERR_UNKNOWN, true, "dark image does not have enough columns.  Returning in image\n");
    1093         return(PM_ERROR_NO_DARK_SUBTRACT);
    1094     }
    1095 
    1096     if (dark != NULL) {
    1097         PS_ASSERT_READOUT_NON_EMPTY(dark, PM_ERROR_NO_DARK_SUBTRACT);
    1098         PS_ASSERT_READOUT_TYPE(dark, PS_TYPE_F32, PM_ERROR_NO_DARK_SUBTRACT);
    1099     }
    1100     return(PM_OKAY);
    1101 }
    1102 
    1103 /******************************************************************************
    1104 p_psDetermineTrimmedImage(): global routine: determines the region of the
    1105 input pmReadout which will be operated on by the various detrend modules.  It
    1106 does a metadata fetch on "CELL.TRIMSEC" for the parent cell of the pmReadout.
    1107  
    1108 Use it this way:
    1109     PS_WARN_PTR_NON_NULL(in->parent);
    1110     if (in->parent != NULL) {
    1111         PS_WARN_PTR_NON_NULL(in->parent->concepts);
    1112     }
    1113     //
    1114     // Determine trimmed image from metadata.
    1115     //
    1116     psImage *trimmedImg = p_psDetermineTrimmedImage(in);
    1117  
    1118 XXX: Create a pmUtils.c file and put this routine there.
    1119  *****************************************************************************/
    1120 psImage *p_psDetermineTrimmedImage(pmReadout *in)
    1121 {
    1122     if ((in->parent == NULL) || (in->parent->concepts == NULL)) {
    1123         psLogMsg(__func__, PS_LOG_WARN,
    1124                  "WARNING: could not determine CELL.TRIMSEC from parent cell Metadata (NULL).\n");
    1125         return(in->image);
    1126     }
    1127 
    1128     psBool rc = false;
    1129     psImage *trimmedImg = NULL;
    1130     psRegion *trimRegion = psMetadataLookupPtr(&rc, in->parent->concepts,
    1131                            "CELL.TRIMSEC");
    1132     if (rc == false) {
    1133         psLogMsg(__func__, PS_LOG_WARN,
    1134                  "WARNING: could not determine CELL.TRIMSEC from parent cell Metadata.\n");
    1135         trimmedImg = in->image;
    1136     } else {
    1137         trimmedImg = psImageSubset(in->image, *trimRegion);
    1138     }
    1139 
    1140     return(trimmedImg);
    1141 }
    1142 
    1143 
    1144 /******************************************************************************
    1145 pmSubtractBias(....): see SDRS for complete specification.
    1146  
    1147 XXX: Code and assert type support: U16, S32, F32.
    1148 XXX: Add trace messages.
    1149  *****************************************************************************/
    1150 pmReadout *pmSubtractBias(
    1151     pmReadout *in,
    1152     void *fitSpec,
    1153     pmFit fit,
    1154     bool overscan,
    1155     psStats *stat,
    1156     int nBin,
    1157     const pmReadout *bias,
    1158     const pmReadout *dark)
     413pmReadout *pmSubtractBias(pmReadout *in, pmOverscanOptions *overscanOpts,
     414                          const pmReadout *bias, const pmReadout *dark)
    1159415{
    1160416    psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4,
    1161417            "---- pmSubtractBias() begin ----\n");
    1162     //
    1163     // Check input parameters, generate warnings and errors.
    1164     //
    1165     if (PM_OKAY != AssertCodeOverscan(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) {
    1166         return(in);
    1167     }
    1168     //
    1169     // Determine trimmed image from metadata.
    1170     //
    1171     psImage *trimmedImg = p_psDetermineTrimmedImage(in);
    1172 
    1173     //
    1174     // Subtract overscan frames if necessary.
    1175     //
    1176     if (overscan == true) {
    1177         pmOverscanAxis overScanAxis = GetOverscanAxis(in);
    1178         //
    1179         //  Create a psStats data structure and determine the highest
    1180         //  priority stats option.
    1181         //
    1182         psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
    1183         if (stat != NULL) {
    1184             myStats->options = GenNewStatOptions(stat);
    1185         }
    1186 
    1187         //
    1188         // Reduce overscan images to a single pixel, then subtract.
    1189         // This code is no longer required as of SDRS 12-09.
    1190         //
    1191         if (overScanAxis == PM_OVERSCAN_ALL) {
    1192             if (false == OverscanReducePixel(trimmedImg, in->bias, myStats)) {
     418    PS_ASSERT_READOUT_NON_NULL(in, NULL);
     419    PS_ASSERT_READOUT_NON_EMPTY(in, NULL);
     420    PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, NULL);
     421
     422    psImage *image = in->image;         // The input image
     423
     424    // Overscan processing
     425    if (overscanOpts) {
     426        // Check for an unallowable pmFit.
     427        if (overscanOpts->fitType != PM_FIT_NONE && overscanOpts->fitType != PM_FIT_POLY_ORD &&
     428                overscanOpts->fitType != PM_FIT_POLY_CHEBY && overscanOpts->fitType != PM_FIT_SPLINE) {
     429            psError(PS_ERR_UNKNOWN, true, "Invalid fit type (%d).  Returning original image.\n", overscanOpts->fitType);
     430            return(in);
     431        }
     432
     433        psList *overscans = in->bias; // List of the overscan images
     434
     435        psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // A new psStats, to avoid clobbering original
     436        myStats->options = GenNewStatOptions(overscanOpts->stat);
     437
     438        // Reduce all overscan pixels to a single value
     439        if (overscanOpts->single) {
     440            psVector *pixels = psVectorAlloc(0, PS_TYPE_F32);
     441            pixels->n = 0;
     442            psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
     443            psImage *overscan = NULL;   // Overscan image from iterator
     444            while ((overscan = psListGetAndIncrement(iter))) {
     445                int index = pixels->n;  // Index
     446                pixels = psVectorRealloc(pixels, pixels->n + overscan->numRows * overscan->numCols);
     447                // XXX Reimplement with memcpy
     448                for (int i = 0; i < overscan->numRows; i++) {
     449                    for (int j = 0; j < overscan->numCols; j++) {
     450                        pixels->data.F32[index++] = overscan->data.F32[i][j];
     451                    }
     452                }
     453
     454            }
     455            psFree(iter);
     456
     457            (void)psVectorStats(myStats, pixels, NULL, NULL, 0);
     458            double reduced = NAN;     // Result of statistics
     459            if (! p_psGetStatValue(myStats, &reduced)) {
     460                psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning input image.\n");
    1193461                return(in);
    1194462            }
    1195             psFree(myStats);
     463            (void)psBinaryOp(image, image, "-", psScalarAlloc((float)reduced, PS_TYPE_F32));
    1196464        } else {
    1197             //
    1198             // Reduce the overscan images to a single overscan vector.
    1199             //
    1200             psVector *overscanVector = OverscanReduce(in->image, overScanAxis,
    1201                                        in->bias, fitSpec,
    1202                                        fit, myStats);
    1203             if (overscanVector == NULL) {
    1204                 psError(PS_ERR_UNKNOWN, false, "Could not reduce overscan images to a single overscan vector.  Returning in image\n");
    1205                 psFree(myStats);
    1206                 return(in);
     465
     466            // We do the regular overscan subtraction
     467
     468            bool readRows = psMetadataLookupBool(NULL, in->parent->concepts, "CELL.READDIR");// Read direction
     469
     470            if (readRows) {
     471                // The read direction is rows
     472                psArray *pixels = psArrayAlloc(image->numRows); // Array of vectors containing pixels
     473                for (int i = 0; i < pixels->n; i++) {
     474                    psVector *values = psVectorAlloc(0, PS_TYPE_F32);
     475                    values->n = 0;
     476                    pixels->data[i] = values;
     477                }
     478
     479                // Pull the pixels out into the vectors
     480                psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
     481                psImage *overscan = NULL; // Overscan image from iterator
     482                while ((overscan = psListGetAndIncrement(iter))) {
     483                    int diff = image->row0 - overscan->row0; // Offset between the two regions
     484                    for (int i = MAX(0,diff); i < MIN(image->numRows, overscan->numRows + diff); i++) {
     485                        // i is row on overscan
     486                        // XXX Reimplement with memcpy
     487                        psVector *values = pixels->data[i];
     488                        int index = values->n; // Index in the vector
     489                        values = psVectorRealloc(values, values->n + overscan->numCols);
     490                        for (int j = 0; j < overscan->numCols; j++) {
     491                            values->data.F32[index++] = overscan->data.F32[i][j];
     492                        }
     493                        values->n += overscan->numCols;
     494                        pixels->data[i] = values; // Update the pointer in case it's moved
     495                    }
     496                }
     497                psFree(iter);
     498
     499                // Reduce the overscans
     500                psVector *reduced = overscanVector(overscanOpts, pixels, myStats);
     501                psFree(pixels);
     502                if (! reduced) {
     503                    return in;
     504                }
     505
     506                // Subtract row by row
     507                for (int i = 0; i < image->numRows; i++) {
     508                    for (int j = 0; j < image->numCols; j++) {
     509                        image->data.F32[i][j] -= reduced->data.F32[i];
     510                    }
     511                }
     512                psFree(reduced);
     513
     514            } else {
     515                // The read direction is columns
     516                psArray *pixels = psArrayAlloc(image->numCols); // Array of vectors containing pixels
     517                for (int i = 0; i < pixels->n; i++) {
     518                    psVector *values = psVectorAlloc(0, PS_TYPE_F32);
     519                    values->n = 0;
     520                    pixels->data[i] = values;
     521                }
     522
     523                // Pull the pixels out into the vectors
     524                psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
     525                psImage *overscan = NULL; // Overscan image from iterator
     526                while ((overscan = psListGetAndIncrement(iter))) {
     527                    int diff = image->col0 - overscan->col0; // Offset between the two regions
     528                    for (int i = MAX(0,diff); i < MIN(image->numCols, overscan->numCols + diff); i++) {
     529                        // i is column on overscan
     530                        // XXX Reimplement with memcpy
     531                        psVector *values = pixels->data[i];
     532                        int index = values->n; // Index in the vector
     533                        values = psVectorRealloc(values, values->n + overscan->numRows);
     534                        for (int j = 0; j < overscan->numRows; j++) {
     535                            values->data.F32[index++] = overscan->data.F32[i][j];
     536                        }
     537                        values->n += overscan->numRows;
     538                        pixels->data[i] = values; // Update the pointer in case it's moved
     539                    }
     540                }
     541                psFree(iter);
     542
     543                // Reduce the overscans
     544                psVector *reduced = overscanVector(overscanOpts, pixels, myStats);
     545                psFree(pixels);
     546                if (! reduced) {
     547                    return in;
     548                }
     549
     550                // Subtract column by column
     551                for (int i = 0; i < image->numCols; i++) {
     552                    for (int j = 0; j < image->numRows; j++) {
     553                        image->data.F32[j][i] -= reduced->data.F32[i];
     554                    }
     555                }
     556                psFree(reduced);
    1207557            }
    1208 
    1209             //
    1210             // Rebin the overscan vector if necessary.
    1211             //
    1212             psS32 newBin = RebinOverscanVector(overscanVector, nBin, myStats);
    1213             if (newBin < 0) {
    1214                 psError(PS_ERR_UNKNOWN, false, "Could rebin the overscan vector.  Returning in image\n");
    1215                 psFree(myStats);
    1216                 return(in);
    1217             }
    1218 
    1219             //
    1220             // If necessary, fit a psPolynomial or psSpline to the overscan vector.
    1221             // Then, unbin the overscan vector to appropriate length for the in image.
    1222             //
    1223             if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) {
    1224                 overscanVector = FitOverscanVectorAndUnbin(trimmedImg, overscanVector, overScanAxis, fitSpec, fit, newBin);
    1225                 if (overscanVector == NULL) {
    1226                     psError(PS_ERR_UNKNOWN, false, "Could not fit the polynomial or spline to the overscan vector.  Returning in image\n");
    1227                     psFree(myStats);
    1228                     return(in);
    1229                 }
    1230             } else {
    1231                 overscanVector = UnbinOverscanVector(trimmedImg, overscanVector, overScanAxis, newBin);
    1232             }
    1233 
    1234             //
    1235             // Subtract the overscan vector from the input image.
    1236             //
    1237             SubtractVectorFromImage(trimmedImg, overscanVector, overScanAxis);
    1238             psFree(myStats);
    1239             psFree(overscanVector);
    1240         }
    1241     }
    1242 
    1243     //
    1244     // Perform bias subtraction if necessary.
    1245     //
    1246     if (bias != NULL) {
    1247         if (PM_OKAY == AssertCodeBias(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) {
    1248             SubtractFrame(in, bias);
    1249         }
    1250     }
    1251 
    1252     //
    1253     // Perform dark subtraction if necessary.
    1254     //
    1255     if (dark != NULL) {
    1256         if (PM_OKAY == AssertCodeDark(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) {
    1257             psBool rc;
    1258             psF32 scale = 0.0;
    1259             if (in->parent != NULL) {
    1260                 scale = psMetadataLookupS32(&rc, in->parent->concepts, "CELL.DARKTIME");
    1261                 if (rc == false) {
    1262                     psLogMsg(__func__, PS_LOG_WARN,
    1263                              "WARNING: pmSubtractBias.(): could not determine CELL.FARKTIME from in->parent metadata.\n");
    1264                 }
    1265             }
    1266             SubtractDarkFrame(in, dark, scale);
    1267         }
    1268     }
    1269 
    1270     //
    1271     // All done.
    1272     //
    1273     psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4,
    1274             "---- pmSubtractBias() exit ----\n");
    1275     return(in);
    1276 }
    1277 
    1278 
     558        }
     559        psFree(myStats);
     560    } // End of overscan subtraction
     561
     562    // Bias frame subtraction
     563    if (bias) {
     564        SubtractFrame(in, bias, 1.0);
     565    }
     566
     567    if (dark) {
     568        // Get the scaling
     569        float inTime = psMetadataLookupF32(NULL, in->parent->concepts, "CELL.DARKTIME");
     570        float darkTime = psMetadataLookupF32(NULL, dark->parent->concepts, "CELL.DARKTIME");
     571        SubtractFrame(in, dark, inTime/darkTime);
     572    }
     573
     574    return in;
     575}
     576
     577
Note: See TracChangeset for help on using the changeset viewer.