IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 17, 2006, 8:10:08 AM (20 years ago)
Author:
magnier
Message:

fixed up conflicts

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/imsubtract/pmSubtractBias.c

    r6511 r6873  
     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.11 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2006-03-04 01:01:33 $
     13 *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
     14 *  @date $Date: 2006-04-17 18:10:08 $
    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     out->knots->n = out->knots->nalloc;
    310     for (psS32 i = 0 ; i < in->knots->n ; i++) {
    311         out->knots->data.F32[i] = in->knots->data.F32[i];
    312     }
    313     /*
    314         out->knots = psVectorCopy(out->knots, in->knots, in->knots->type.type);
    315     */
    316 
    317     out->p_psDeriv2 = (psF32 *) psAlloc((in->n + 1) * sizeof(psF32));
    318     // XXX: use memcpy
    319     for (psS32 i = 0 ; i < (in->n + 1) ; i++) {
    320         out->p_psDeriv2[i] = in->p_psDeriv2[i];
    321     }
    322 
    323     return(true);
    324 }
    325 
     217
     218
     219#if 0
    326220/******************************************************************************
    327221ScaleOverscanVector(): this routine takes as input an arbitrary vector,
    … …  
    330224    PM_FIT_POLYNOMIAL: fit a polynomial to the entire input vector data.
    331225    PM_FIT_SPLINE: fit splines to the input vector data.
    332 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?
    333229 *****************************************************************************/
    334 static psVector *ScaleOverscanVector(
    335     psVector *overscanVector,
    336     psS32 n,
    337     void *fitSpec,
    338     pmFit fit)
     230static psVector *ScaleOverscanVector(psVector *overscanVector,
     231                                     psS32 n,
     232                                     void *fitSpec,
     233                                     pmFit fit)
    339234{
    340235    psTrace(".psModule.pmSubtracBias.ScaleOverscanVector", 4,
    341236            "---- ScaleOverscanVector() begin (%d -> %d) ----\n", overscanVector->n, n);
     237    //    PS_VECTOR_PRINT_F32(overscanVector);
    342238
    343239    if (NULL == overscanVector) {
    … …  
    347243    // Allocate the new vector.
    348244    psVector *newVec = psVectorAlloc(n, PS_TYPE_F32);
    349     newVec->n = newVec->nalloc;
     245
    350246    //
    351247    // If the new vector is the same size as the old, simply copy the data.
    352248    //
    353249    if (n == overscanVector->n) {
    354         return(psVectorCopy(newVec, overscanVector, PS_TYPE_F32));
    355     }
     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;
    356257    psF32 x;
    357 
     258    psS32 i;
    358259    if (fit == PM_FIT_POLYNOMIAL) {
    359260        // Fit a polynomial to the old overscan vector.
    360         psPolynomial1D *myPoly = (psPolynomial1D *) fitSpec;
     261        myPoly = (psPolynomial1D *) fitSpec;
    361262        PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
    362         PS_ASSERT_POLY1D(myPoly, NULL);
    363263        myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL);
    364264        if (myPoly == NULL) {
    365             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");
    366266            return(NULL);
    367267        }
    … …  
    370270        // of the old vector, use the fitted polynomial to determine the
    371271        // interpolated value at that point, and set the new vector.
    372         for (psS32 i=0;i<n;i++) {
     272        for (i=0;i<n;i++) {
    373273            x = ((psF32) i) * ((psF32) overscanVector->n) / ((psF32) n);
    374274            newVec->data.F32[i] = psPolynomial1DEval(myPoly, x);
    375275        }
    376276    } 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
    377285        //
    378286        // NOTE: Since the X arg in the psVectorFitSpline1D() function is NULL,
    … …  
    380288        // properly when doing the spline eval.
    381289        //
    382         psSpline1D *mySpline = psVectorFitSpline1D(NULL, overscanVector);
     290        //        mySpline = psVectorFitSpline1D(mySpline, NULL, overscanVector, NULL);
     291        mySpline = psVectorFitSpline1D(NULL, overscanVector);
    383292        if (mySpline == NULL) {
    384             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");
    385294            return(NULL);
    386295        }
     296        //        PS_PRINT_SPLINE(mySpline);
    387297
    388298        // For each element of the new vector, convert the x-ordinate to that
    389         // of the old vector, use the fitted spline to determine the
     299        // of the old vector, use the fitted polynomial to determine the
    390300        // interpolated value at that point, and set the new vector.
    391         for (psS32 i=0;i<n;i++) {
     301        for (i=0;i<n;i++) {
    392302            // Scale to [0 : overscanVector->n - 1]
    393303            x = ((psF32) i) * ((psF32) (overscanVector->n-1)) / ((psF32) n);
    394304            newVec->data.F32[i] = psSpline1DEval(mySpline, x);
    395305        }
    396 
    397         psSpline1D *ptrSpline = (psSpline1D *) fitSpec;
    398         if (ptrSpline != NULL) {
    399             // Copy the resulting spline fit into ptrSpline.
    400             PS_ASSERT_SPLINE(ptrSpline, NULL);
    401             SplineCopy(ptrSpline, mySpline);
    402         }
    403         psFree(mySpline);
     306        if (mustFreeSpline ==1) {
     307            psFree(mySpline);
     308        }
     309        //        PS_VECTOR_PRINT_F32(newVec);
     310
     311
    404312    } else {
    405313        psError(PS_ERR_UNKNOWN, true, "unknown fit type.  Returning NULL.\n");
    … …  
    413321}
    414322
     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
    415410/******************************************************************************
     411XXX: The SDRS does not specify type support.  F32 is implemented here.
    416412 *****************************************************************************/
    417 static psS32 GetOverscanSize(
    418     psImage *inImg,
    419     pmOverscanAxis overScanAxis)
    420 {
    421     if (overScanAxis == PM_OVERSCAN_ROWS) {
    422         return(inImg->numCols);
    423     } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
    424         return(inImg->numRows);
    425     } else if (overScanAxis == PM_OVERSCAN_ALL) {
    426         return(1);
    427     }
    428     return(0);
    429 }
    430 
    431 /******************************************************************************
    432 GetOverscanAxis(in) this private routine determines the appropiate overscan
    433 axis from the parent cell metadata.
    434  
    435 XXX: Verify the READDIR corresponds with my overscan axis.
    436  *****************************************************************************/
    437 static pmOverscanAxis GetOverscanAxis(pmReadout *in)
    438 {
    439     psBool rc;
    440     if ((in->parent != NULL) && (in->parent->concepts)) {
    441         psS32 dir = psMetadataLookupS32(&rc, in->parent->concepts, "CELL.READDIR");
    442         if (rc == true) {
    443             if (dir == 1) {
    444                 return(PM_OVERSCAN_ROWS);
    445             } else if (dir == 2) {
    446                 return(PM_OVERSCAN_COLUMNS);
    447             } else if (dir == 3) {
    448                 return(PM_OVERSCAN_ALL);
    449             }
    450         }
    451     }
    452 
    453     psLogMsg(__func__, PS_LOG_WARN,
    454              "WARNING: pmSubtractBias.(): could not determine CELL.READDIR from in->parent metadata.  Setting overscan axis to PM_OVERSCAN_NONE.\n");
    455     return(PM_OVERSCAN_NONE);
    456 }
    457 
    458 /******************************************************************************
    459 my_psListLength(list): determine the length of a psList.
    460  
    461 XXX: Put this elsewhere.
    462  
    463 XXX: psList.h now has a version of this function.  Use that instead.
    464  *****************************************************************************/
    465 
    466 static psS32 my_psListLength(
    467     psList *list)
    468 {
    469     psS32 length = 0;
    470     psListElem *tmpElem = (psListElem *) list->head;
    471     while (NULL != tmpElem) {
    472         tmpElem = tmpElem->next;
    473         length++;
    474     }
    475     return(length);
    476 }
    477 
    478 /******************************************************************************
    479 Note: this isn't needed anymore as of psModule SDRS 12-09.
    480  *****************************************************************************/
    481 static psBool OverscanReducePixel(
    482     psImage *in,
    483     psList *bias,
    484     psStats *myStats)
    485 {
    486     PS_ASSERT_PTR_NON_NULL(in, NULL);
    487     PS_ASSERT_PTR_NON_NULL(bias, NULL);
    488     PS_ASSERT_PTR_NON_NULL(bias->head, NULL);
    489     PS_ASSERT_PTR_NON_NULL(myStats, NULL);
    490 
    491     // Allocate a psVector with one element per overscan image.
    492     psS32 numOverscanImages = my_psListLength(bias);
    493     psVector *statsAll = psVectorAlloc(numOverscanImages, PS_TYPE_F32);
    494     statsAll->n = statsAll->nalloc;
    495     psListElem *tmpOverscan = (psListElem *) bias->head;
    496     psS32 i = 0;
    497     psF64 statValue;
    498     //
    499     // We loop through each overscan image, calculating the specified
    500     // statistic on that image.
    501     //
    502     while (NULL != tmpOverscan) {
    503         psImage *myOverscanImage = (psImage *) tmpOverscan->data;
    504 
    505         PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL);
    506         myStats = psImageStats(myStats, myOverscanImage, NULL, (psMaskType)0xffffffff);
    507         if (myStats == NULL) {
    508             psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation.  Returning in image.\n");
    509             psFree(statsAll);
    510             return(false);
    511         }
    512         if (false == p_psGetStatValue(myStats, &statValue)) {
    513             psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
    514             psFree(statsAll);
    515             return(false);
    516         }
    517         statsAll->data.F32[i] = statValue;
    518         i++;
    519         tmpOverscan = tmpOverscan->next;
    520     }
    521 
    522     //
    523     // We reduce the individual stats for each overscan image to
    524     // a single psF32.
    525     //
    526     myStats = psVectorStats(myStats, statsAll, NULL, NULL, 0);
    527     if (myStats == NULL) {
    528         psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation.  Returning in image.\n");
    529         psFree(statsAll);
    530         return(false);
    531     }
    532     if (false == p_psGetStatValue(myStats, &statValue)) {
    533         psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
    534         psFree(statsAll);
    535         return(false);
    536     }
    537 
    538     //
    539     // Subtract the result and return.
    540     //
    541     ImageSubtractScalar(in, statValue);
    542     psFree(statsAll);
    543     return(in);
    544 }
    545 
    546 /******************************************************************************
    547 ReduceOverscanImageToCol(overscanImage, myStats): This private routine reduces
    548 a single psImage to a column by combining all pixels from each row into a
    549 single pixel via requested statistic in myStats.
    550  *****************************************************************************/
    551 static psVector *ReduceOverscanImageToCol(
    552     psImage *overscanImage,
    553     psStats *myStats)
    554 {
    555     psF64 statValue;
    556     psVector *tmpRow = psVectorAlloc(overscanImage->numCols, PS_TYPE_F32);
    557     psVector *tmpCol = psVectorAlloc(overscanImage->numRows, PS_TYPE_F32);
    558     tmpRow->n = tmpRow->nalloc;
    559     tmpCol->n = tmpCol->nalloc;
    560 
    561     //
    562     // For each row, we store all pixels in that row into a temporary psVector,
    563     // then we run psVectorStats() on that vector.
    564     //
    565     for (psS32 i=0;i<overscanImage->numRows;i++) {
    566         for (psS32 j=0;j<overscanImage->numCols;j++) {
    567             tmpRow->data.F32[j] = overscanImage->data.F32[i][j];
    568         }
    569 
    570         psStats *rc = psVectorStats(myStats, tmpRow, NULL, NULL, 0);
    571         if (rc == NULL) {
    572             psError(PS_ERR_UNKNOWN, true, "psVectorStats() could not perform requested statistical operation.  Returning in image.\n");
    573             return(NULL);
    574         }
    575 
    576         if (false ==  p_psGetStatValue(rc, &statValue)) {
    577             psError(PS_ERR_UNKNOWN, true, "p_psGetStatValue() could not determine result from requested statistical operation.  Returning in image.\n");
    578             return(NULL);
    579         }
    580 
    581         tmpCol->data.F32[i] = (psF32) statValue;
    582     }
    583     psFree(tmpRow);
    584 
    585     return(tmpCol);
    586 }
    587 
    588 /******************************************************************************
    589 ReduceOverscanImageToCol(overscanImage, myStats): This private routine reduces
    590 a single psImage to a row by combining all pixels from each column into a
    591 single pixel via requested statistic in myStats.
    592  *****************************************************************************/
    593 static psVector *ReduceOverscanImageToRow(
    594     psImage *overscanImage,
    595     psStats *myStats)
    596 {
    597     psF64 statValue;
    598     psVector *tmpRow = psVectorAlloc(overscanImage->numCols, PS_TYPE_F32);
    599     psVector *tmpCol = psVectorAlloc(overscanImage->numRows, PS_TYPE_F32);
    600     tmpRow->n = tmpRow->nalloc;
    601     tmpCol->n = tmpCol->nalloc;
    602 
    603     //
    604     // For each column, we store all pixels in that column into a temporary psVector,
    605     // then we run psVectorStats() on that vector.
    606     //
    607     for (psS32 i=0;i<overscanImage->numCols;i++) {
    608         for (psS32 j=0;j<overscanImage->numRows;j++) {
    609             tmpCol->data.F32[j] = overscanImage->data.F32[j][i];
    610         }
    611 
    612         psStats *rc = psVectorStats(myStats, tmpCol, NULL, NULL, 0);
    613         if (rc == NULL) {
    614             psError(PS_ERR_UNKNOWN, true, "psVectorStats() could not perform requested statistical operation.  Returning in image.\n");
    615             return(NULL);
    616         }
    617 
    618         if (false ==  p_psGetStatValue(rc, &statValue)) {
    619             psError(PS_ERR_UNKNOWN, true, "p_psGetStatValue() could not determine result from requested statistical operation.  Returning in image.\n");
    620             return(NULL);
    621         }
    622 
    623         tmpRow->data.F32[i] = (psF32) statValue;
    624     }
    625     psFree(tmpCol);
    626 
    627     return(tmpRow);
    628 }
    629 
    630 /******************************************************************************
    631 OverscanReduce(vecSize, bias, myStats): This private routine takes a psList of
    632 overscan images (in bias) and reduces them to a single psVector via the
    633 specified psStats struct.  The vector is then scaled to the length or the
    634 row/column in inImg.
    635  *****************************************************************************/
    636 static psVector* OverscanReduce(
    637     psImage *inImg,
    638     pmOverscanAxis overScanAxis,
    639     psList *bias,
    640     void *fitSpec,
    641     pmFit fit,
    642     psStats *myStats)
    643 {
    644     if ((overScanAxis != PM_OVERSCAN_ROWS) && (overScanAxis != PM_OVERSCAN_COLUMNS)) {
    645         psError(PS_ERR_UNKNOWN, true, "overScanAxis must be PM_OVERSCAN_ROWS or PM_OVERSCAN_COLUMNS\n");
    646         return(NULL);
    647     }
    648     PS_ASSERT_PTR_NON_NULL(inImg, NULL);
    649     PS_ASSERT_PTR_NON_NULL(bias, NULL);
    650     PS_ASSERT_PTR_NON_NULL(bias->head, NULL);
    651     PS_ASSERT_PTR_NON_NULL(myStats, NULL);
    652     //
    653     // Allocate a psVector for the output of this routine.
    654     //
    655     psS32 vecSize = GetOverscanSize(inImg, overScanAxis);
    656     psVector *overscanVector = psVectorAlloc(vecSize, PS_TYPE_F32);
    657     overscanVector->n = overscanVector->nalloc;
    658 
    659     //
    660     // Allocate an array of psVectors with one psVector per element of the
    661     // final oversan column vector.  These psVectors will be used with
    662     // psStats to reduce the multiple elements from each overscan column
    663     // vector to a single final column vector.
    664     //
    665     psS32 numOverscanImages = my_psListLength(bias);
    666     psVector **overscanVectors = (psVector **) psAlloc(numOverscanImages * sizeof(psVector *));
    667     //    (*overscanVectors)->n = (*overscanVectors)->nalloc;
    668     for (psS32 i = 0 ; i < numOverscanImages ; i++) {
    669         overscanVectors[i] = NULL;
    670     }
    671 
    672     //
    673     // We iterate through the list of overscan images.  For each image,
    674     // we reduce it to a single column or row.  Save the overscan vector
    675     // in overscanVectors[].
    676     //
    677     psListElem *tmpOverscan = (psListElem *) bias->head;
    678     psS32 overscanID = 0;
    679     while (tmpOverscan != NULL) {
    680         psImage *tmpOverscanImage = (psImage *) tmpOverscan->data;
    681         if (overScanAxis == PM_OVERSCAN_ROWS) {
    682             overscanVectors[overscanID] = ReduceOverscanImageToRow(tmpOverscanImage, myStats);
    683         } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
    684             overscanVectors[overscanID] = ReduceOverscanImageToCol(tmpOverscanImage, myStats);
    685         }
    686 
    687         tmpOverscan = tmpOverscan->next;
    688         overscanID++;
    689     }
    690 
    691     //
    692     // For each overscan vector, if necessary, we scale that column or
    693     // row to vecSize.  Note: we should have already ensured that the
    694     // fit is poly or spline.
    695     //
    696     for (psS32 i = 0 ; i < numOverscanImages ; i++) {
    697         psVector *tmpOverscanVector = overscanVectors[i];
    698 
    699         if (tmpOverscanVector->n != vecSize) {
    700             overscanVectors[i] = ScaleOverscanVector(tmpOverscanVector, vecSize, fitSpec, fit);
    701             if (overscanVectors[i] == NULL) {
    702                 psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector.\n");
    703                 for (psS32 i = 0 ; i < numOverscanImages ; i++) {
    704                     psFree(overscanVectors[i]);
    705                 }
    706                 psFree(overscanVectors);
    707                 psFree(tmpOverscanVector);
    708                 return(NULL);
    709             }
    710             psFree(tmpOverscanVector);
    711         }
    712     }
    713 
    714     //
    715     // We collect all elements in the overscan vectors for the various
    716     // overscan images into a single psVector (tmpVec).  Then we call
    717     // psStats on that vector to determine the final values for the
    718     // overscan vector.
    719     //
    720     psVector *tmpVec = psVectorAlloc(numOverscanImages, PS_TYPE_F32);
    721     tmpVec->n = tmpVec->nalloc;
    722     psF64 statValue;
    723     for (psS32 i = 0 ; i < vecSize ; i++) {
    724         // Collect the i-th elements from each overscan vector into a single vector.
    725         for (psS32 j = 0 ; j < numOverscanImages ; j++) {
    726             tmpVec->data.F32[j] = overscanVectors[j]->data.F32[i];
    727         }
    728 
    729         if (NULL == psVectorStats(myStats, tmpVec, NULL, NULL, 0)) {
    730             psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
    731             for (psS32 i = 0 ; i < numOverscanImages ; i++) {
    732                 psFree(overscanVectors[i]);
    733             }
    734             psFree(overscanVectors);
    735             psFree(tmpVec);
    736             return(NULL);
    737         }
    738         if (false == p_psGetStatValue(myStats, &statValue)) {
    739             psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
    740             for (psS32 i = 0 ; i < numOverscanImages ; i++) {
    741                 psFree(overscanVectors[i]);
    742             }
    743             psFree(overscanVectors);
    744             psFree(tmpVec);
    745             return(NULL);
    746         }
    747 
    748         overscanVector->data.F32[i] = (psF32) statValue;
    749     }
    750 
    751     //
    752     // We're done.  Free the intermediate overscan vectors.
    753     //
    754     psFree(tmpVec);
    755     for (psS32 i = 0 ; i < numOverscanImages ; i++) {
    756         psFree(overscanVectors[i]);
    757     }
    758     psFree(overscanVectors);
    759 
    760     //
    761     // Return the computed overscanVector
    762     //
    763     return(overscanVector);
    764 }
    765 
    766 /******************************************************************************
    767 RebinOverscanVector(overscanVector, nBinOrig, myStats): this private routine
    768 takes groups of nBinOrig elements in the input vector, combines them into a
    769 single pixel via myStats and psVectorStats(), and then outputs a vector of
    770 those pixels.
    771  *****************************************************************************/
    772 static psS32 RebinOverscanVector(
    773     psVector *overscanVector,
    774     psS32 nBinOrig,
    775     psStats *myStats)
    776 {
    777     psF64 statValue;
    778     psS32 nBin;
    779     if ((nBinOrig > 1) && (nBinOrig < overscanVector->n)) {
    780         psS32 numBins = 1+((overscanVector->n)/nBinOrig);
    781         psVector *myBin = psVectorAlloc(numBins, PS_TYPE_F32);
    782         psVector *binVec = psVectorAlloc(nBinOrig, PS_TYPE_F32);
    783         myBin->n = myBin->nalloc;
    784         binVec->n = binVec->nalloc;
    785 
    786         for (psS32 i=0;i<numBins;i++) {
    787             for(psS32 j=0;j<nBinOrig;j++) {
    788                 if (overscanVector->n > ((i*nBinOrig)+j)) {
    789                     binVec->data.F32[j] = overscanVector->data.F32[(i*nBinOrig)+j];
    790                 } else {
    791                     // XXX: we get here if nBinOrig does not evenly divide
    792                     // the overscanVector vector.  This is the last bin.  Should
    793                     // we change the binVec->n to acknowledge that?
    794                     binVec->n = j;
    795                 }
    796             }
    797             psStats *rc = psVectorStats(myStats, binVec, NULL, NULL, 0);
    798             if (rc == NULL) {
    799                 psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
    800                 return(-1);
    801             }
    802             if (false ==  p_psGetStatValue(rc, &statValue)) {
    803                 psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
    804                 return(-1);
    805             }
    806             myBin->data.F32[i] = statValue;
    807         }
    808 
    809         // Change the effective size of overscanVector.
    810         overscanVector->n = numBins;
    811         for (psS32 i=0;i<numBins;i++) {
    812             overscanVector->data.F32[i] = myBin->data.F32[i];
    813         }
    814         psFree(binVec);
    815         psFree(myBin);
    816         nBin = nBinOrig;
    817     } else {
    818         nBin = 1;
    819     }
    820 
    821     return(nBin);
    822 }
    823 
    824 /******************************************************************************
    825 FitOverscanVectorAndUnbin(inImg, overscanVector, overScanAxis, fitSpec, fit,
    826 nBin):  this private routine fits a psPolynomial or psSpline to the overscan
    827 vector.  It then creates a new vector, with a size determined by the input
    828 image, evaluates the psPolynomial or psSpline at each element in that vector,
    829 then returns that vector.
    830  *****************************************************************************/
    831 static psVector *FitOverscanVectorAndUnbin(
    832     psImage *inImg,
    833     psVector *overscanVector,
    834     pmOverscanAxis overScanAxis,
    835     void *fitSpec,
    836     pmFit fit,
    837     psS32 nBin)
    838 {
    839     psPolynomial1D* myPoly = NULL;
    840     psSpline1D *mySpline = NULL;
    841     //
    842     // Fit a polynomial or spline to the overscan vector.
    843     //
    844     if (fit == PM_FIT_POLYNOMIAL) {
    845         myPoly = (psPolynomial1D *) fitSpec;
    846         PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
    847         PS_ASSERT_POLY1D(myPoly, NULL);
    848         myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL);
    849         if (myPoly == NULL) {
    850             psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to overscan vector.  Returning NULL.\n");
    851             return(NULL);
    852         }
    853     } else if (fit == PM_FIT_SPLINE) {
    854         mySpline = psVectorFitSpline1D(NULL, overscanVector);
    855         if (mySpline == NULL) {
    856             psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to overscan vector.  Returning NULL.\n");
    857             return(NULL);
    858         }
    859         if (fitSpec != NULL) {
    860             // Copy the resulting spline fit into fitSpec.
    861             psSpline1D *ptrSpline = (psSpline1D *) fitSpec;
    862             PS_ASSERT_SPLINE(ptrSpline, NULL);
    863             SplineCopy(ptrSpline, mySpline);
    864         }
    865     }
    866 
    867     //
    868     // Evaluate the poly/spline at each pixel in the overscan row/column.
    869     //
    870     psS32 vecSize = GetOverscanSize(inImg, overScanAxis);
    871     psVector *newVec = psVectorAlloc(vecSize, PS_TYPE_F32);
    872     newVec->n = newVec->nalloc;
    873     if ((nBin > 1) && (nBin < overscanVector->n)) {
    874         for (psS32 i = 0 ; i < vecSize ; i++) {
    875             if (fit == PM_FIT_POLYNOMIAL) {
    876                 newVec->data.F32[i] = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin));
    877             } else if (fit == PM_FIT_SPLINE) {
    878                 newVec->data.F32[i] = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin));
    879             }
    880         }
    881     } else {
    882         for (psS32 i = 0 ; i < vecSize ; i++) {
    883             if (fit == PM_FIT_POLYNOMIAL) {
    884                 newVec->data.F32[i] = psPolynomial1DEval(myPoly, (psF32) i);
    885             } else if (fit == PM_FIT_SPLINE) {
    886                 newVec->data.F32[i] = psSpline1DEval(mySpline, (psF32) i);
    887             }
    888         }
    889     }
    890 
    891     psFree(mySpline);
    892     psFree(overscanVector);
    893     return(newVec);
    894 }
    895 
    896 
    897 
    898 /******************************************************************************
    899 UnbinOverscanVector(inImg, overscanVector, overScanAxis, nBin):  this private
    900 routine takes a psVector overscanVector that was previously binned by a factor
    901 of nBin, and then expands it to its original size, duplicated elements nBin
    902 times for each element in the input vector overscanVector.
    903  *****************************************************************************/
    904 static psVector *UnbinOverscanVector(
    905     psImage *inImg,
    906     psVector *overscanVector,
    907     pmOverscanAxis overScanAxis,
    908     psS32 nBin)
    909 {
    910     psS32 vecSize = 0;
    911 
    912     if (overScanAxis == PM_OVERSCAN_ROWS) {
    913         vecSize = inImg->numCols;
    914     } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
    915         vecSize = inImg->numRows;
    916     }
    917 
    918     psVector *newVec = psVectorAlloc(vecSize, PS_TYPE_F32);
    919     newVec->n = newVec->nalloc;
    920     for (psS32 i = 0 ; i < vecSize ; i++) {
    921         newVec->data.F32[i] = overscanVector->data.F32[i/nBin];
    922     }
    923 
    924     psFree(overscanVector);
    925     return(newVec);
    926 }
    927 
    928 
    929 /******************************************************************************
    930 SubtractVectorFromImage(inImg, overscanVector, overScanAxis):  this private
    931 routine subtracts the overscanVector column-wise or row-wise from inImg.
    932  *****************************************************************************/
    933 static psImage *SubtractVectorFromImage(
    934     psImage *inImg,
    935     psVector *overscanVector,
    936     pmOverscanAxis overScanAxis)
    937 {
    938     //
    939     // Subtract overscan vector row-wise from the image.
    940     //
    941     if (overScanAxis == PM_OVERSCAN_ROWS) {
    942         for (psS32 i=0;i<inImg->numCols;i++) {
    943             for (psS32 j=0;j<inImg->numRows;j++) {
    944                 inImg->data.F32[j][i]-= overscanVector->data.F32[i];
    945             }
    946         }
    947     }
    948 
    949     //
    950     // Subtract overscan vector column-wise from the image.
    951     //
    952     if (overScanAxis == PM_OVERSCAN_COLUMNS) {
    953         for (psS32 i=0;i<inImg->numRows;i++) {
    954             for (psS32 j=0;j<inImg->numCols;j++) {
    955                 inImg->data.F32[i][j]-= overscanVector->data.F32[i];
    956             }
    957         }
    958     }
    959 
    960     return(inImg);
    961 }
    962 
    963 
    964 
    965 typedef enum {
    966     PM_ERROR_NO_SUBTRACTION,
    967     PM_WARNING_NO_SUBTRACTION,
    968     PM_ERROR_NO_BIAS_SUBTRACT,
    969     PM_WARNING_NO_BIAS_SUBTRACT,
    970     PM_ERROR_NO_DARK_SUBTRACT,
    971     PM_WARNING_NO_DARK_SUBTRACT,
    972     PM_OKAY
    973 } pmSubtractBiasAssertStatus;
    974 /******************************************************************************
    975 AssertCodeOverscan(....) this private routine verifies that the various input
    976 parameters to pmSubtractBias() are correct for overscan subtraction.
    977  *****************************************************************************/
    978 pmSubtractBiasAssertStatus AssertCodeOverscan(
    979     pmReadout *in,
    980     void *fitSpec,
    981     pmFit fit,
    982     bool overscan,
    983     psStats *stat,
    984     int nBinOrig,
    985     const pmReadout *bias,
    986     const pmReadout *dark)
    987 {
    988 
    989     PS_ASSERT_READOUT_NON_NULL(in, PM_ERROR_NO_SUBTRACTION);
    990     PS_ASSERT_READOUT_NON_EMPTY(in, PM_ERROR_NO_SUBTRACTION);
    991     PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, PM_ERROR_NO_SUBTRACTION);
    992     PS_WARN_PTR_NON_NULL(in->parent);
    993     if (in->parent != NULL) {
    994         PS_WARN_PTR_NON_NULL(in->parent->concepts);
    995     }
    996 
    997     if (overscan == true) {
    998         pmOverscanAxis overScanAxis = GetOverscanAxis(in);
    999         PS_ASSERT_PTR_NON_NULL(stat, PM_ERROR_NO_SUBTRACTION);
    1000         PS_ASSERT_PTR_NON_NULL(in->bias, PM_ERROR_NO_SUBTRACTION);
    1001         PS_ASSERT_PTR_NON_NULL(in->bias->head, PM_ERROR_NO_SUBTRACTION);
    1002         //
    1003         // Check the type, size of each bias image.
    1004         //
    1005         psListElem *tmpOverscan = (psListElem *) in->bias->head;
    1006         psS32 numOverscans = 0;
    1007         while (NULL != tmpOverscan) {
    1008             numOverscans++;
    1009             psImage *myOverscanImage = (psImage *) tmpOverscan->data;
    1010             PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, PM_ERROR_NO_SUBTRACTION);
    1011             // XXX: Get this right with the rows and columns.
    1012             if (overScanAxis == PM_OVERSCAN_ROWS) {
    1013                 if (myOverscanImage->numRows != in->image->numRows) {
    1014                     psLogMsg(__func__, PS_LOG_WARN,
    1015                              "WARNING: pmSubtractBias.(): overscan image (# %d) has %d rows, input image has %d rows\n",
    1016                              numOverscans, myOverscanImage->numCols, in->image->numRows);
    1017                     if (fit == PM_FIT_NONE) {
    1018                         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");
    1019                         return(PM_ERROR_NO_SUBTRACTION);
    1020                     }
    1021                 }
    1022             } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
    1023                 if (myOverscanImage->numCols != in->image->numCols) {
    1024                     psLogMsg(__func__, PS_LOG_WARN,
    1025                              "WARNING: pmSubtractBias.(): overscan image (# %d) has %d columns, input image has %d columns\n",
    1026                              numOverscans, myOverscanImage->numCols, in->image->numCols);
    1027                     if (fit == PM_FIT_NONE) {
    1028                         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");
    1029                         return(PM_ERROR_NO_SUBTRACTION);
    1030                     }
    1031                 }
    1032             } else if (overScanAxis != PM_OVERSCAN_ALL) {
    1033                 psError(PS_ERR_UNKNOWN, true, "Must specify and overscan axis.\n");
    1034                 return(PM_ERROR_NO_SUBTRACTION);
    1035             }
    1036             tmpOverscan = tmpOverscan->next;
    1037         }
    1038     } else {
    1039         if (fit != PM_FIT_NONE) {
    1040             psLogMsg(__func__, PS_LOG_WARN,
    1041                      "WARNING: pmSubtractBias.(): overscan is FALSE and fit is not PM_FIT_NONE.\n");
    1042             return(PM_WARNING_NO_SUBTRACTION);
    1043         }
    1044     }
    1045 
    1046     // XXX: I do not like the following spec since it's useless to specify
    1047     // a psSpline as the fitSpec.
    1048     if (0) {
    1049         if ((fitSpec == NULL) &&
    1050                 ((fit != PM_FIT_NONE) || (overscan == true))) {
    1051             psError(PS_ERR_UNKNOWN, true, "fitSpec is NULL and fit is not PM_FIT_NONE or overscan is TRUE.\n");
    1052             return(PM_ERROR_NO_SUBTRACTION);
    1053         }
    1054     }
    1055 
    1056     return(PM_OKAY);
    1057 }
    1058 
    1059 /******************************************************************************
    1060 AssertCodeBias(....) this private routine verifies that the various input
    1061 parameters to pmSubtractBias() are correct for bias subtraction.
    1062  *****************************************************************************/
    1063 static pmSubtractBiasAssertStatus AssertCodeBias(
    1064     pmReadout *in,
    1065     void *fitSpec,
    1066     pmFit fit,
    1067     bool overscan,
    1068     psStats *stat,
    1069     int nBinOrig,
    1070     const pmReadout *bias,
    1071     const pmReadout *dark)
    1072 {
    1073     if ((in->image->numRows + in->row0 - bias->row0) > bias->image->numRows) {
    1074         psError(PS_ERR_UNKNOWN,true, "bias image does not have enough rows.  Returning in image\n");
    1075         return(PM_ERROR_NO_BIAS_SUBTRACT);
    1076     }
    1077     if ((in->image->numCols + in->col0 - bias->col0) > bias->image->numCols) {
    1078         psError(PS_ERR_UNKNOWN,true, "bias image does not have enough columns.  Returning in image\n");
    1079         return(PM_ERROR_NO_BIAS_SUBTRACT);
    1080     }
    1081 
    1082     if (bias != NULL) {
    1083         PS_ASSERT_READOUT_NON_EMPTY(bias, PM_ERROR_NO_BIAS_SUBTRACT);
    1084         PS_ASSERT_READOUT_TYPE(bias, PS_TYPE_F32, PM_ERROR_NO_DARK_SUBTRACT);
    1085     }
    1086     return(PM_OKAY);
    1087 }
    1088 
    1089 /******************************************************************************
    1090 AssertCodeDark(....) this private routine verifies that the various input
    1091 parameters to pmSubtractBias() are correct for dark subtraction.
    1092  *****************************************************************************/
    1093 pmSubtractBiasAssertStatus AssertCodeDark(
    1094     pmReadout *in,
    1095     void *fitSpec,
    1096     pmFit fit,
    1097     bool overscan,
    1098     psStats *stat,
    1099     int nBinOrig,
    1100     const pmReadout *bias,
    1101     const pmReadout *dark)
    1102 {
    1103     if ((in->image->numRows + in->row0 - dark->row0) > dark->image->numRows) {
    1104         psError(PS_ERR_UNKNOWN, true, "dark image does not have enough rows.  Returning in image\n");
    1105         return(PM_ERROR_NO_DARK_SUBTRACT);
    1106     }
    1107     if ((in->image->numCols + in->col0 - dark->col0) > dark->image->numCols) {
    1108         psError(PS_ERR_UNKNOWN, true, "dark image does not have enough columns.  Returning in image\n");
    1109         return(PM_ERROR_NO_DARK_SUBTRACT);
    1110     }
    1111 
    1112     if (dark != NULL) {
    1113         PS_ASSERT_READOUT_NON_EMPTY(dark, PM_ERROR_NO_DARK_SUBTRACT);
    1114         PS_ASSERT_READOUT_TYPE(dark, PS_TYPE_F32, PM_ERROR_NO_DARK_SUBTRACT);
    1115     }
    1116     return(PM_OKAY);
    1117 }
    1118 
    1119 /******************************************************************************
    1120 p_psDetermineTrimmedImage(): global routine: determines the region of the
    1121 input pmReadout which will be operated on by the various detrend modules.  It
    1122 does a metadata fetch on "CELL.TRIMSEC" for the parent cell of the pmReadout.
    1123  
    1124 Use it this way:
    1125     PS_WARN_PTR_NON_NULL(in->parent);
    1126     if (in->parent != NULL) {
    1127         PS_WARN_PTR_NON_NULL(in->parent->concepts);
    1128     }
    1129     //
    1130     // Determine trimmed image from metadata.
    1131     //
    1132     psImage *trimmedImg = p_psDetermineTrimmedImage(in);
    1133  
    1134 XXX: Create a pmUtils.c file and put this routine there.
    1135  *****************************************************************************/
    1136 psImage *p_psDetermineTrimmedImage(pmReadout *in)
    1137 {
    1138     if ((in->parent == NULL) || (in->parent->concepts == NULL)) {
    1139         psLogMsg(__func__, PS_LOG_WARN,
    1140                  "WARNING: could not determine CELL.TRIMSEC from parent cell Metadata (NULL).\n");
    1141         return(in->image);
    1142     }
    1143 
    1144     psBool rc = false;
    1145     psImage *trimmedImg = NULL;
    1146     psRegion *trimRegion = psMetadataLookupPtr(&rc, in->parent->concepts,
    1147                            "CELL.TRIMSEC");
    1148     if (rc == false) {
    1149         psLogMsg(__func__, PS_LOG_WARN,
    1150                  "WARNING: could not determine CELL.TRIMSEC from parent cell Metadata.\n");
    1151         trimmedImg = in->image;
    1152     } else {
    1153         trimmedImg = psImageSubset(in->image, *trimRegion);
    1154     }
    1155 
    1156     return(trimmedImg);
    1157 }
    1158 
    1159 
    1160 /******************************************************************************
    1161 pmSubtractBias(....): see SDRS for complete specification.
    1162  
    1163 XXX: Code and assert type support: U16, S32, F32.
    1164 XXX: Add trace messages.
    1165  *****************************************************************************/
    1166 pmReadout *pmSubtractBias(
    1167     pmReadout *in,
    1168     void *fitSpec,
    1169     pmFit fit,
    1170     bool overscan,
    1171     psStats *stat,
    1172     int nBin,
    1173     const pmReadout *bias,
    1174     const pmReadout *dark)
     413pmReadout *pmSubtractBias(pmReadout *in, pmOverscanOptions *overscanOpts,
     414                          const pmReadout *bias, const pmReadout *dark)
    1175415{
    1176416    psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4,
    1177417            "---- pmSubtractBias() begin ----\n");
    1178     //
    1179     // Check input parameters, generate warnings and errors.
    1180     //
    1181     if (PM_OKAY != AssertCodeOverscan(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) {
    1182         return(in);
    1183     }
    1184     //
    1185     // Determine trimmed image from metadata.
    1186     //
    1187     psImage *trimmedImg = p_psDetermineTrimmedImage(in);
    1188 
    1189     //
    1190     // Subtract overscan frames if necessary.
    1191     //
    1192     if (overscan == true) {
    1193         pmOverscanAxis overScanAxis = GetOverscanAxis(in);
    1194         //
    1195         //  Create a psStats data structure and determine the highest
    1196         //  priority stats option.
    1197         //
    1198         psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
    1199         if (stat != NULL) {
    1200             myStats->options = GenNewStatOptions(stat);
    1201         }
    1202 
    1203         //
    1204         // Reduce overscan images to a single pixel, then subtract.
    1205         // This code is no longer required as of SDRS 12-09.
    1206         //
    1207         if (overScanAxis == PM_OVERSCAN_ALL) {
    1208             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");
    1209461                return(in);
    1210462            }
    1211             psFree(myStats);
     463            (void)psBinaryOp(image, image, "-", psScalarAlloc((float)reduced, PS_TYPE_F32));
    1212464        } else {
    1213             //
    1214             // Reduce the overscan images to a single overscan vector.
    1215             //
    1216             psVector *overscanVector = OverscanReduce(in->image, overScanAxis,
    1217                                        in->bias, fitSpec,
    1218                                        fit, myStats);
    1219             if (overscanVector == NULL) {
    1220                 psError(PS_ERR_UNKNOWN, false, "Could not reduce overscan images to a single overscan vector.  Returning in image\n");
    1221                 psFree(myStats);
    1222                 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);
    1223557            }
    1224 
    1225             //
    1226             // Rebin the overscan vector if necessary.
    1227             //
    1228             psS32 newBin = RebinOverscanVector(overscanVector, nBin, myStats);
    1229             if (newBin < 0) {
    1230                 psError(PS_ERR_UNKNOWN, false, "Could rebin the overscan vector.  Returning in image\n");
    1231                 psFree(myStats);
    1232                 return(in);
    1233             }
    1234 
    1235             //
    1236             // If necessary, fit a psPolynomial or psSpline to the overscan vector.
    1237             // Then, unbin the overscan vector to appropriate length for the in image.
    1238             //
    1239             if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) {
    1240                 overscanVector = FitOverscanVectorAndUnbin(trimmedImg, overscanVector, overScanAxis, fitSpec, fit, newBin);
    1241                 if (overscanVector == NULL) {
    1242                     psError(PS_ERR_UNKNOWN, false, "Could not fit the polynomial or spline to the overscan vector.  Returning in image\n");
    1243                     psFree(myStats);
    1244                     return(in);
    1245                 }
    1246             } else {
    1247                 overscanVector = UnbinOverscanVector(trimmedImg, overscanVector, overScanAxis, newBin);
    1248             }
    1249 
    1250             //
    1251             // Subtract the overscan vector from the input image.
    1252             //
    1253             SubtractVectorFromImage(trimmedImg, overscanVector, overScanAxis);
    1254             psFree(myStats);
    1255             psFree(overscanVector);
    1256         }
    1257     }
    1258 
    1259     //
    1260     // Perform bias subtraction if necessary.
    1261     //
    1262     if (bias != NULL) {
    1263         if (PM_OKAY == AssertCodeBias(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) {
    1264             SubtractFrame(in, bias);
    1265         }
    1266     }
    1267 
    1268     //
    1269     // Perform dark subtraction if necessary.
    1270     //
    1271     if (dark != NULL) {
    1272         if (PM_OKAY == AssertCodeDark(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) {
    1273             psBool rc;
    1274             psF32 scale = 0.0;
    1275             if (in->parent != NULL) {
    1276                 scale = psMetadataLookupS32(&rc, in->parent->concepts, "CELL.DARKTIME");
    1277                 if (rc == false) {
    1278                     psLogMsg(__func__, PS_LOG_WARN,
    1279                              "WARNING: pmSubtractBias.(): could not determine CELL.FARKTIME from in->parent metadata.\n");
    1280                 }
    1281             }
    1282             SubtractDarkFrame(in, dark, scale);
    1283         }
    1284     }
    1285 
    1286     //
    1287     // All done.
    1288     //
    1289     psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4,
    1290             "---- pmSubtractBias() exit ----\n");
    1291     return(in);
    1292 }
    1293 
    1294 
     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.