IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 13, 2024, 11:43:14 AM (2 years ago)
Author:
eugene
Message:

updated version of pmOverscan which can support 2D bias subtraction

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/psModules.20240412/src/detrend/pmOverscan.c

    r42653 r42654  
    1919#define SMOOTH_NSIGMA 4.0               // Number of Gaussian sigma the smoothing kernel extends
    2020
    21 static void pmOverscanOptionsFree(pmOverscanOptions *options)
    22 {
    23     psFree(options->stat);
    24     psFree(options->poly);
    25     psFree(options->spline);
    26 }
    27 
    28 pmOverscanOptions *pmOverscanOptionsAlloc(bool single, pmFit fitType, unsigned int order, psStats *stat,
    29                                           int boxcar, float gauss)
    30 {
    31     pmOverscanOptions *opts = psAlloc(sizeof(pmOverscanOptions));
    32     psMemSetDeallocator(opts, (psFreeFunc)pmOverscanOptionsFree);
    33 
    34     // Inputs
    35     opts->single = single;
    36     opts->constant = false;
    37     opts->fitType = fitType;
    38     opts->order = order;
    39     opts->stat = psMemIncrRefCounter(stat);
    40 
    41     opts->minValid = 0.0; // default value if not defined
    42     opts->maxValid = (float) 0x10000; // default value if not defined
    43     opts->maskVal = 0x0001; // default value if not defined
    44 
    45     // Smoothing
    46     opts->boxcar = boxcar;
    47     opts->gauss = gauss;
    48 
    49     // Outputs
    50     opts->poly = NULL;
    51     opts->spline = NULL;
    52 
    53     return opts;
    54 }
    55 
    56 // Produce an overscan vector from an array of pixels
    57 psVector *pmOverscanVector(float *chi2, // chi^2 from fit
    58                            pmOverscanOptions *overscanOpts, // Overscan options
    59                            const psArray *pixels, // Array of vectors containing the pixel values
    60                            psStats *myStats // Statistic to use in reducing the overscan
    61     )
    62 {
    63     assert(overscanOpts);
    64     assert(pixels);
    65     assert(myStats);
    66 
    67     psStatsOptions statistic = psStatsSingleOption(myStats->options); // Statistic to use
    68     assert(statistic != 0);
    69 
    70     // Reduce the overscans
    71     psVector *reduced = psVectorAlloc(pixels->n, PS_TYPE_F32); // Overscan for each row
    72     psVector *ordinate = psVectorAlloc(pixels->n, PS_TYPE_F32); // Ordinate
    73     psVector *mask = psVectorAlloc(pixels->n, PS_TYPE_VECTOR_MASK); // Mask for fitting
    74 
    75     for (int i = 0; i < pixels->n; i++) {
    76         psVector *values = pixels->data[i]; // Vector with overscan values
    77         if (values->n > 0) {
    78             mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
    79             ordinate->data.F32[i] = 2.0*(float)i/(float)pixels->n - 1.0; // Scale to [-1,1]
    80             if (!psVectorStats(myStats, values, NULL, NULL, 0)) {
    81                 psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    82                 return false;
    83             }
    84             reduced->data.F32[i] = psStatsGetValue(myStats, statistic);
    85         } else if (overscanOpts->fitType == PM_FIT_NONE) {
    86             psError(PS_ERR_UNKNOWN, true, "The overscan is not supplied for all points on the "
    87                     "image, and no fit is requested.\n");
    88             return NULL;
    89         } else {
    90             // We'll fit this one out
    91             mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 1;
    92         }
    93     }
    94 
    95     // Smooth the reduced vector
    96     if (overscanOpts->boxcar > 0) {
    97         psVector *smoothed = psVectorBoxcar(NULL, reduced, overscanOpts->boxcar); // Smoothed vector
    98         psFree(reduced);
    99         reduced = smoothed;
    100     }
    101     if (isfinite(overscanOpts->gauss) && overscanOpts->gauss > 0) {
    102         if (overscanOpts->boxcar > 0) {
    103             psWarning("Gaussian smoothing the boxcar smoothed overscan --- you asked for it.");
    104         }
    105         psVector *smoothed = psVectorSmooth(NULL, reduced, overscanOpts->gauss, SMOOTH_NSIGMA);
    106         psFree(reduced);
    107         reduced = smoothed;
    108     }
    109 
    110     // Fit the overscan, if required
    111     psVector *fitted = NULL;                   // Fitted overscan values
    112     switch (overscanOpts->fitType) {
    113       case PM_FIT_NONE:
    114         // No fitting --- that's easy.
    115         fitted = psMemIncrRefCounter(reduced);
    116         break;
    117       case PM_FIT_POLY_ORD:
    118         if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order ||
    119                                    overscanOpts->poly->type != PS_POLYNOMIAL_ORD)) {
    120             psFree(overscanOpts->poly);
    121             overscanOpts->poly = NULL;
    122         }
    123         if (! overscanOpts->poly) {
    124             overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, overscanOpts->order);
    125         }
    126         psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate);
    127         fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
    128         break;
    129       case PM_FIT_POLY_CHEBY:
    130         if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order ||
    131                                    overscanOpts->poly->type != PS_POLYNOMIAL_CHEB)) {
    132             psFree(overscanOpts->poly);
    133             overscanOpts->poly = NULL;
    134         }
    135         if (! overscanOpts->poly) {
    136             overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_CHEB, overscanOpts->order);
    137         }
    138         psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate);
    139         fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
    140         break;
    141       case PM_FIT_SPLINE:
    142 
    143         // XXX I don't think psSpline1D is up to scratch yet --- it has no mask, and it assumes
    144         // a knot for every input point.  it needs an argument like 'number of knots' for the
    145         // output spline.  EAM: still true 2023.01.22
    146 
    147         // overscanOpts->spline = psVectorFitSpline1D(reduced, ordinate);
    148         // fitted = psSpline1DEvalVector(overscanOpts->spline, ordinate);
    149         psError(PS_ERR_UNKNOWN, true, "Spline overscan fitting is broken\n");
    150         break;
    151       default:
    152         psError(PS_ERR_UNKNOWN, true, "Unknown value for the fitting type: %d\n", overscanOpts->fitType);
    153         return NULL;
    154         break;
    155     }
    156 
    157     if (chi2) {
    158         *chi2 = 0.0;                    // chi^2 (sort of)
    159         for (int i = 0; i < reduced->n; i++) {
    160             *chi2 += PS_SQR(fitted->data.F32[i] - reduced->data.F32[i]);
    161         }
    162     }
    163 
    164     psFree(reduced);
    165     psFree(ordinate);
    166     psFree(mask);
    167 
    168     return fitted;
    169 }
    170 
    171 bool pmOverscanUpdateHeader (pmHDU *hdu, pmOverscanOptions *overscanOpts, float chi2) {
    172 
    173     psString comment = NULL;    // Comment to add
    174 
    175     switch (overscanOpts->fitType) {
    176       case PM_FIT_POLY_ORD:
    177       case PM_FIT_POLY_CHEBY: {
    178           psStringAppend(&comment, "Overscan fit (chi2: %.2f): ", chi2);
    179           psPolynomial1D *poly = overscanOpts->poly; // The polynomial
    180           for (int i = 0; i < poly->nX; i++) {
    181               psStringAppend(&comment, "%.1f ", poly->coeff[i]);
    182           }
    183           psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
    184           psFree(comment);
    185           comment = NULL;
    186 
    187           // write metadata header value
    188           psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE,
    189                            "Overscan value", poly->coeff[0]);
    190           psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE,
    191                            "Overscan stdev", poly->coeffErr[0]);
    192           break;
    193       }
    194       case PM_FIT_SPLINE: {
    195         /*
    196           psSpline1D *spline = overscanOpts->spline; // The spline
    197           for (int i = 0; i < spline->n; i++) {
    198               psStringAppend(&comment, "Overscan fit (chi2: %.2f) %d:", chi2, i);
    199               psPolynomial1D *poly = spline->spline[i]; // i-th polynomial
    200               for (int j = 0; j < poly->nX; j++) {
    201                   psStringAppend(&comment, "%.1f ", poly->coeff[i]);
    202               }
    203               psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,
    204                                comment, "");
    205               psFree(comment);
    206               comment = NULL;
    207           }
    208         */
    209           // write metadata header value
    210           psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE,
    211                            "Overscan value", NAN);
    212           psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE,
    213                            "Overscan stdev", NAN);
    214           break;
    215       }
    216       case PM_FIT_NONE:
    217         break;
    218       default:
    219         psAbort("Should never get here!!!\n");
    220     }
    221     return true;
    222 }
    223 
    224 static int saveVectors = 1;
     21static void pmOverscanOptionsFree(pmOverscanOptions *options);
     22static void pmOverscanStatOptionsFree(pmOverscanStatOptions *options);
     23bool pmOverscanUpdateHeaderVector (pmReadout *input, pmHDU *hdu, pmOverscanOptions *overscanOpts, psVector *reduced);
    22524
    22625bool pmOverscanSubtract (pmReadout *input, pmOverscanOptions *overscanOpts) {
    … …  
    23635    if (overscanOpts->constant) {
    23736        // write metadata header value
    238         psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value",
    239                          overscanOpts->value);
     37        psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", overscanOpts->value);
    24038        psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", NAN);
    24139
    … …  
    24947
    25048    // Check for an unallowable pmFit.
    251     if (overscanOpts->fitType != PM_FIT_NONE && overscanOpts->fitType != PM_FIT_POLY_ORD &&
    252         overscanOpts->fitType != PM_FIT_POLY_CHEBY && overscanOpts->fitType != PM_FIT_SPLINE) {
    253         psError(PS_ERR_UNKNOWN, true, "Invalid fit type (%d).  Returning original image.\n",
    254                 overscanOpts->fitType);
    255         return false;
     49    if (overscanOpts->primary) {
     50        if (overscanOpts->primary->fitType != PM_FIT_NONE &&
     51            overscanOpts->primary->fitType != PM_FIT_POLY_ORD &&
     52            overscanOpts->primary->fitType != PM_FIT_POLY_CHEBY &&
     53            overscanOpts->primary->fitType != PM_FIT_SPLINE) {
     54            psError(PS_ERR_UNKNOWN, true, "Invalid fit type (%d).  Returning original image.\n",
     55                    overscanOpts->primary->fitType);
     56            return false;
     57        }
     58    }
     59    if (overscanOpts->secondary) {
     60        if (overscanOpts->secondary->fitType != PM_FIT_NONE &&
     61            overscanOpts->secondary->fitType != PM_FIT_POLY_ORD &&
     62            overscanOpts->secondary->fitType != PM_FIT_POLY_CHEBY &&
     63            overscanOpts->secondary->fitType != PM_FIT_SPLINE) {
     64            psError(PS_ERR_UNKNOWN, true, "Invalid fit type (2D) (%d).  Returning original image.\n",
     65                    overscanOpts->secondary->fitType);
     66            return false;
     67        }
    25668    }
    25769
    25870    psList *overscans = input->bias; // List of the overscan images
    259 
    260     psStatsOptions statistic = psStatsSingleOption(overscanOpts->stat->options); // Statistic to use
    261     if (statistic == 0) {
    262         psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics options set: %p\n",
    263                 overscanOpts->stat);
    264         return false;
    265     }
    266     psStats *stats = psStatsAlloc(statistic); // A new psStats, to avoid clobbering original
    267 
    268     psString comment = NULL;    // Comment to add
    269     psStringAppend(&comment, "Subtracting overscan (stat %x; type %x; order %d)",
    270                    statistic, overscanOpts->fitType, overscanOpts->order);
    271     psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,
    272                      comment, "");
    273     psTrace ("psModules.detrend", 4, "%s\n", comment);
    274     psFree(comment);
    27571
    27672    // Reduce all overscan pixels to a single value
    27773    if (overscanOpts->single) {
     74
     75        // extract overscan pixels to a single vector
    27876        psVector *pixels = psVectorAlloc(0, PS_TYPE_F32);
    27977        psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
    … …  
    29189        psFree(iter);
    29290
     91        // statistic to be calculated
     92        psStatsOptions statistic = psStatsSingleOption(overscanOpts->primary->stat->options); // Statistic to use
     93        if (!statistic) {
     94            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics options set: %p\n",
     95                    overscanOpts->primary->stat);
     96            return false;
     97        }
     98        psStats *stats = psStatsAlloc(statistic); // A new psStats, to avoid clobbering original
     99
    293100        if (!psVectorStats(stats, pixels, NULL, NULL, 0)) {
    294101            psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    … …  
    304111
    305112        // write metadata header value
    306         psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value",
    307                         reduced);
     113        // XXX EAM : this could / should write the stdev of the overscan region
     114        psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", reduced);
    308115        psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", NAN);
    309116
    … …  
    315122        overscan = NULL;   // Overscan image from iterator
    316123        while ((overscan = psListGetAndIncrement(iter))) {
    317           psBinaryOp(overscan, overscan, "-", psMemIncrRefCounter(reducedScalar)); // NOTE: psBinaryOp frees arg2 if it a scalar, so we need to bump to re-use
     124            psBinaryOp(overscan, overscan, "-", psMemIncrRefCounter(reducedScalar)); // NOTE: psBinaryOp frees arg2 if it a scalar, so we need to bump to re-use
    318125        }
    319126        psFree(iter);
    … …  
    341148    // We are performing a row-by-row overscan subtraction
    342149    int cellreaddir = psMetadataLookupS32(&mdok, input->parent->concepts, "CELL.READDIR"); // Read direction
    343     if ((cellreaddir != -1) && (cellreaddir != 1) && (cellreaddir != 2)) {
    344         psError(PS_ERR_UNKNOWN, true, "CELL.READDIR must be 1 (rows) or 2 (cols) or 0 (2D)\n");
     150    if ((cellreaddir != 1) && (cellreaddir != 2)) {
     151        psError(PS_ERR_UNKNOWN, true, "CELL.READDIR must be 1 (rows) or 2 (cols)\n");
    345152        return false;
    346153    }
    … …  
    349156
    350157    // adjust operation depending on the read direction : need to re-org pixels for columns
    351     if (cellreaddir == 1) {
     158    if (!overscanOpts->TwoD && (cellreaddir == 1)) {
    352159        // The read direction is rows
    353160        psArray *pixels = psArrayAlloc(image->numRows); // Array of vectors containing pixels
    … …  
    381188
    382189        // Reduce the overscans
    383         psVector *reduced = pmOverscanVector(&chi2, overscanOpts, pixels, stats);
     190        psVector *reduced = pmOverscanVector(&chi2, overscanOpts->primary, pixels);
    384191        psFree(pixels);
    385192        if (! reduced) {
    386193            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate overscan vector.\n");
    387             psFree(stats);
    388             return false;
    389         }
    390 
    391         // generate stats of overscan vector for header
    392         {
    393             psString comment = NULL;    // Comment to add
    394             psStats *vectorStats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
    395             if (!psVectorStats (vectorStats, reduced, NULL, NULL, 0)) {
    396                 psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    397                 return false;
    398             }
    399             psStringAppend(&comment, "Mean Overscan value: %f", vectorStats->sampleMean);
    400             psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
    401             psFree(comment);
    402 
    403             // write metadata header value
    404             psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan mean", vectorStats->sampleMean);
    405             psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", vectorStats->sampleStdev);
    406 
    407             // EAM 2022.03.29 : if the calculated overscan value is below the threshold,
    408             // declare the readout dead and mask
    409          
    410             if ((vectorStats->sampleMean < overscanOpts->minValid) || (vectorStats->sampleMean > overscanOpts->maxValid)) {
    411                 fprintf (stderr, "bad overscan (2) %f, masking readout\n", vectorStats->sampleMean);
    412                 psImage *mask = input->mask;
    413                 for (int y = 0; y < mask->numRows; y++) {
    414                     for (int x = 0; x < mask->numCols; x++) {
    415                         mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= overscanOpts->maskVal;
    416                     }
    417                 }
    418             }
    419 
    420             psFree (vectorStats);
     194            return false;
     195        }
     196
     197        if (!pmOverscanUpdateHeaderVector (input, hdu, overscanOpts, reduced)) {
     198            psError(PS_ERR_UNKNOWN, false, "failure to update header");
     199            return false;
    421200        }
    422201
    … …  
    430209        // subtract from the overscan regions
    431210        {
    432           psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
    433           psImage *overscan = NULL; // Overscan image from iterator
    434           while ((overscan = psListGetAndIncrement(iter))) {
    435             // the overscan and image might not be aligned.
    436             int diff = overscan->row0 - image->row0; // Offset between the two regions
    437             for (int i = PS_MAX(0,diff); i < PS_MIN(image->numRows, overscan->numRows + diff); i++) {
    438               int j = i - diff;
    439               // i is row on image
    440               // j is row on overscan
    441               for (int k = 0; k < overscan->numCols; k++) {
    442                 overscan->data.F32[j][k] -= reduced->data.F32[j];
    443               }
    444             }
    445           }
    446           psFree(iter);
     211            psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
     212            psImage *overscan = NULL; // Overscan image from iterator
     213            while ((overscan = psListGetAndIncrement(iter))) {
     214                // the overscan and image might not be aligned.
     215                int diff = overscan->row0 - image->row0; // Offset between the two regions
     216                for (int i = PS_MAX(0,diff); i < PS_MIN(image->numRows, overscan->numRows + diff); i++) {
     217                    int j = i - diff;
     218                    // i is row on image
     219                    // j is row on overscan
     220                    for (int k = 0; k < overscan->numCols; k++) {
     221                        overscan->data.F32[j][k] -= reduced->data.F32[j];
     222                    }
     223                }
     224            }
     225            psFree(iter);
    447226        }
    448227        psFree(reduced);
    449228    }
    450     if (cellreaddir == 2) {
     229
     230    if (!overscanOpts->TwoD && (cellreaddir == 2)) {
    451231        // The read direction is columns
    452232        psArray *pixels = psArrayAlloc(image->numCols); // Array of vectors containing pixels
    … …  
    480260
    481261        // Reduce the overscans
    482         psVector *reduced = pmOverscanVector(&chi2, overscanOpts, pixels, stats);
     262        psVector *reduced = pmOverscanVector(&chi2, overscanOpts->primary, pixels);
    483263        psFree(pixels);
    484264        if (! reduced) {
    485265            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate overscan vector.\n");
    486             psFree(stats);
    487             return false;
    488         }
    489 
    490         // generate stats of overscan vector for header
    491         {
    492           psString comment = NULL;    // Comment to add
    493           psStats *vectorStats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
    494           if (!psVectorStats (vectorStats, reduced, NULL, NULL, 0)) {
    495               psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    496               return false;
    497           }
    498           psStringAppend(&comment, "Mean Overscan value: %f", vectorStats->sampleMean);
    499           psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
    500           psFree(comment);
    501 
    502           // write metadata header value
    503           psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan mean", vectorStats->sampleMean);
    504           psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", vectorStats->sampleStdev);
    505 
    506           // EAM 2022.03.29 : if the calculated overscan value is below the threshold,
    507           // declare the readout dead and mask
    508          
    509           if ((vectorStats->sampleMean < overscanOpts->minValid) || (vectorStats->sampleMean > overscanOpts->maxValid)) {
    510               fprintf (stderr, "bad overscan (3) %f, masking readout\n", vectorStats->sampleMean);
    511               psImage *mask = input->mask;
    512               for (int y = 0; y < mask->numRows; y++) {
    513                   for (int x = 0; x < mask->numCols; x++) {
    514                       mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= overscanOpts->maskVal;
    515                   }
    516               }
    517           }
    518           psFree (vectorStats);
     266            return false;
     267        }
     268
     269        if (!pmOverscanUpdateHeaderVector (input, hdu, overscanOpts, reduced)) {
     270            psError(PS_ERR_UNKNOWN, false, "failure to update header");
     271            return false;
    519272        }
    520273
    521274        // Subtract column by column
    522275        for (int j = 0; j < image->numRows; j++) {
    523           for (int i = 0; i < image->numCols; i++) {
     276            for (int i = 0; i < image->numCols; i++) {
    524277                image->data.F32[j][i] -= reduced->data.F32[i];
    525278            }
    … …  
    528281        // subtract from the overscan regions
    529282        {
    530           psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
    531           psImage *overscan = NULL; // Overscan image from iterator
    532           while ((overscan = psListGetAndIncrement(iter))) {
    533             // the overscan and image might not be aligned.
    534             int diff = overscan->col0 - image->col0; // Offset between the two regions
    535             for (int i = PS_MAX(0,diff); i < PS_MIN(image->numCols, overscan->numCols + diff); i++) {
    536               int j = i - diff;
    537               // i is col on image
    538               // j is col on overscan
    539               for (int k = 0; i < overscan->numRows; k++) {
    540                 overscan->data.F32[k][j] -= reduced->data.F32[j];
    541               }
    542             }
    543           }
    544           psFree(iter);
     283            psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
     284            psImage *overscan = NULL; // Overscan image from iterator
     285            while ((overscan = psListGetAndIncrement(iter))) {
     286                // the overscan and image might not be aligned.
     287                int diff = overscan->col0 - image->col0; // Offset between the two regions
     288                for (int i = PS_MAX(0,diff); i < PS_MIN(image->numCols, overscan->numCols + diff); i++) {
     289                    int j = i - diff;
     290                    // i is col on image
     291                    // j is col on overscan
     292                    for (int k = 0; i < overscan->numRows; k++) {
     293                        overscan->data.F32[k][j] -= reduced->data.F32[j];
     294                    }
     295                }
     296            }
     297            psFree(iter);
    545298        }
    546299
    … …  
    552305    // readout direction from the top overscan region and a second
    553306    // 1D pattern in the cross direction from the overscan
    554     if (cellreaddir == -1) {
     307    if (overscanOpts->TwoD && (cellreaddir == 1)) {
     308        psAssert (overscanOpts->secondary, "2D overscan subtraction requires OVERSCAN.2D parameters");
    555309        // we require 2 overscan regions, and they must match these directions:
    556310        if (overscans->n != 2) {
    … …  
    558312            psLogMsg (__func__, PS_LOG_ERROR, "e.g.: CELL.BIASSEC STR [591:624,1:608],[1:590,598:608]");
    559313            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to generate overscan vector.\n");
    560             psFree(stats);
    561314            return false;
    562315        }
    … …  
    620373        // Reduce the overscans
    621374        // XXX need to save 2 different chi-square values
    622         psVector *xReduced = pmOverscanVector(&chi2, overscanOpts, xscanPixels, stats);
     375        psVector *yReduced = pmOverscanVector(&chi2, overscanOpts->primary, yscanPixels);
     376        psFree(yscanPixels);
     377        if (!yReduced) {
     378            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate y-dir overscan vector.\n");
     379            return false;
     380        }
     381        if (!pmOverscanUpdateHeaderVector (input, hdu, overscanOpts, yReduced)) {
     382            psError(PS_ERR_UNKNOWN, false, "failure to update header");
     383            return false;
     384        }
     385
     386        // Reduce the overscans
     387        // XXX need to save 2 different chi-square values
     388        psVector *xReduced = pmOverscanVector(&chi2, overscanOpts->secondary, xscanPixels);
    623389        psFree(xscanPixels);
    624390        if (!xReduced) {
    625391            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate x-dir overscan vector.\n");
    626             psFree(stats);
    627             return false;
    628         }
    629 
    630         // Reduce the overscans
    631         // XXX need to save 2 different chi-square values
    632         psVector *yReduced = pmOverscanVector(&chi2, overscanOpts, yscanPixels, stats);
    633         psFree(yscanPixels);
    634         if (!yReduced) {
    635             psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate y-dir overscan vector.\n");
    636             psFree(stats);
    637             return false;
    638         }
    639 
    640 
    641         if (saveVectors) {
    642             FILE *foutX = fopen ("xdir.dat", "w");
    643             int fdX = fileno (foutX);
    644             p_psVectorPrint (fdX, xReduced, "xscan");
    645             fclose (foutX);
    646             FILE *foutY = fopen ("ydir.dat", "w");
    647             int fdY = fileno (foutY);
    648             p_psVectorPrint (fdY, yReduced, "yscan");
    649             fclose (foutY);
    650             saveVectors = 0;
    651         }
    652 
    653         // write info about stats of y-dir overscan vector in header
    654         {
    655             psString comment = NULL;    // Comment to add
    656             psStats *vectorStats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
    657             if (!psVectorStats (vectorStats, yReduced, NULL, NULL, 0)) {
    658                 psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    659                 return false;
    660             }
    661             psStringAppend(&comment, "Mean Overscan value: %f", vectorStats->sampleMean);
    662             psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
    663             psFree(comment);
    664 
    665             // write metadata header value
    666             psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan mean", vectorStats->sampleMean);
    667             psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", vectorStats->sampleStdev);
    668 
    669             // EAM 2022.03.29 : if the calculated overscan value is below the threshold,
    670             // declare the readout dead and mask
    671          
    672             if ((vectorStats->sampleMean < overscanOpts->minValid) || (vectorStats->sampleMean > overscanOpts->maxValid)) {
    673                 fprintf (stderr, "bad overscan (2) %f, masking readout\n", vectorStats->sampleMean);
    674                 psImage *mask = input->mask;
    675                 for (int y = 0; y < mask->numRows; y++) {
    676                     for (int x = 0; x < mask->numCols; x++) {
    677                         mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= overscanOpts->maskVal;
    678                     }
    679                 }
    680             }
    681             psFree (vectorStats);
     392            return false;
    682393        }
    683394
    … …  
    719430        psFree(yReduced);
    720431    }
    721 
    722     pmOverscanUpdateHeader (hdu, overscanOpts, chi2);
     432    // pmOverscanUpdateHeader (hdu, overscanOpts, chi2);
     433    return true;
     434
     435} // End of overscan subtraction
     436
     437static void pmOverscanOptionsFree(pmOverscanOptions *options)
     438{
     439    psFree(options->primary);
     440    psFree(options->secondary);
     441}
     442
     443static void pmOverscanStatOptionsFree(pmOverscanStatOptions *options)
     444{
     445    psFree(options->stat);
     446    psFree(options->poly);
     447    psFree(options->spline);
     448}
     449
     450// Globally defined
     451pmOverscanStatOptions *pmOverscanStatOptionsAlloc(void)
     452{
     453    pmOverscanStatOptions *opts = psAlloc(sizeof(pmOverscanStatOptions));
     454    psMemSetDeallocator(opts, (psFreeFunc)pmOverscanStatOptionsFree);
     455
     456    // Inputs
     457    opts->fitType = PM_FIT_NONE;
     458    opts->order   = 0;
     459    opts->stat    = NULL;
     460
     461    // Smoothing
     462    opts->boxcar = 0;
     463    opts->gauss  = 0.0;
     464
     465    // Outputs
     466    opts->poly = NULL;
     467    opts->spline = NULL;
     468
     469    return opts;
     470}
     471
     472// Globally defined
     473pmOverscanOptions *pmOverscanOptionsAlloc(void)
     474{
     475    pmOverscanOptions *opts = psAlloc(sizeof(pmOverscanOptions));
     476    psMemSetDeallocator(opts, (psFreeFunc)pmOverscanOptionsFree);
     477
     478    // Inputs
     479    opts->single   = false;
     480    opts->constant = false;
     481    opts->TwoD     = false;
     482
     483    opts->value    = 0.0;
     484
     485    opts->minValid = 0.0; // default value if not defined
     486    opts->maxValid = (float) 0x10000; // default value if not defined
     487    opts->maskVal  = 0x0001; // default value if not defined
     488
     489    // stat options
     490    opts->primary = NULL;
     491    opts->secondary = NULL;
     492
     493    return opts;
     494}
     495
     496// Produce an overscan vector from an array of pixels
     497psVector *pmOverscanVector(float *chi2, // chi^2 from fit
     498                           pmOverscanStatOptions *overscanOpts, // Overscan statistic options
     499                           const psArray *pixels // Array of vectors containing the pixel values
     500    )
     501{
     502    assert(overscanOpts);
     503    assert(pixels);
     504   
     505    // statisctic to be calculated
     506    psStatsOptions statistic = psStatsSingleOption(overscanOpts->stat->options); // Statistic to use
     507    if (!statistic) {
     508        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics options set: %p\n", overscanOpts->stat);
     509        return false;
     510    }
     511    psStats *stats = psStatsAlloc(statistic); // A new psStats, to avoid clobbering original
     512
     513    // Reduce the overscans
     514    psVector *reduced = psVectorAlloc(pixels->n, PS_TYPE_F32); // Overscan for each row
     515    psVector *ordinate = psVectorAlloc(pixels->n, PS_TYPE_F32); // Ordinate
     516    psVector *mask = psVectorAlloc(pixels->n, PS_TYPE_VECTOR_MASK); // Mask for fitting
     517
     518    for (int i = 0; i < pixels->n; i++) {
     519        psVector *values = pixels->data[i]; // Vector with overscan values
     520        if (values->n > 0) {
     521            mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
     522            ordinate->data.F32[i] = 2.0*(float)i/(float)pixels->n - 1.0; // Scale to [-1,1]
     523            if (!psVectorStats(stats, values, NULL, NULL, 0)) {
     524                psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     525                goto escape;
     526            }
     527            reduced->data.F32[i] = psStatsGetValue(stats, statistic);
     528        } else {
     529            if (overscanOpts->fitType == PM_FIT_NONE) {
     530                psError(PS_ERR_UNKNOWN, true, "The overscan is not supplied for all points on the image, and no fit is requested.\n");
     531                goto escape;
     532            } else {
     533                // We'll fit this one out
     534                mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 1;
     535            }
     536        }
     537    }
     538    // Smooth the reduced vector
     539    if (overscanOpts->boxcar > 0) {
     540        psVector *smoothed = psVectorBoxcar(NULL, reduced, overscanOpts->boxcar); // Smoothed vector
     541        psFree(reduced);
     542        reduced = smoothed;
     543    }
     544    if (isfinite(overscanOpts->gauss) && overscanOpts->gauss > 0) {
     545        if (overscanOpts->boxcar > 0) {
     546            psWarning("Gaussian smoothing the boxcar smoothed overscan --- you asked for it.");
     547        }
     548        psVector *smoothed = psVectorSmooth(NULL, reduced, overscanOpts->gauss, SMOOTH_NSIGMA);
     549        psFree(reduced);
     550        reduced = smoothed;
     551    }
     552
     553    // Fit the overscan, if required
     554    psVector *fitted = NULL;                   // Fitted overscan values
     555    switch (overscanOpts->fitType) {
     556      case PM_FIT_NONE:
     557        // No fitting --- that's easy.
     558        fitted = psMemIncrRefCounter(reduced);
     559        break;
     560      case PM_FIT_POLY_ORD:
     561        if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order ||
     562                                   overscanOpts->poly->type != PS_POLYNOMIAL_ORD)) {
     563            psFree(overscanOpts->poly);
     564            overscanOpts->poly = NULL;
     565        }
     566        if (! overscanOpts->poly) {
     567            overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, overscanOpts->order);
     568        }
     569        psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate);
     570        fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
     571        break;
     572      case PM_FIT_POLY_CHEBY:
     573        if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order ||
     574                                   overscanOpts->poly->type != PS_POLYNOMIAL_CHEB)) {
     575            psFree(overscanOpts->poly);
     576            overscanOpts->poly = NULL;
     577        }
     578        if (! overscanOpts->poly) {
     579            overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_CHEB, overscanOpts->order);
     580        }
     581        psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate);
     582        fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
     583        break;
     584      case PM_FIT_SPLINE:
     585
     586        // XXX I don't think psSpline1D is up to scratch yet --- it has no mask, and it assumes
     587        // a knot for every input point.  it needs an argument like 'number of knots' for the
     588        // output spline.  EAM: still true 2023.01.22
     589
     590        // overscanOpts->spline = psVectorFitSpline1D(reduced, ordinate);
     591        // fitted = psSpline1DEvalVector(overscanOpts->spline, ordinate);
     592        psError(PS_ERR_UNKNOWN, true, "Spline overscan fitting is broken\n");
     593        break;
     594      default:
     595        psError(PS_ERR_UNKNOWN, true, "Unknown value for the fitting type: %d\n", overscanOpts->fitType);
     596        goto escape;
     597    }
     598
     599    if (chi2) {
     600        *chi2 = 0.0;                    // chi^2 (sort of)
     601        for (int i = 0; i < reduced->n; i++) {
     602            *chi2 += PS_SQR(fitted->data.F32[i] - reduced->data.F32[i]);
     603        }
     604    }
     605
     606    psFree(reduced);
     607    psFree(ordinate);
     608    psFree(mask);
    723609    psFree(stats);
    724 
     610    return fitted;
     611
     612escape:
     613    psFree(reduced);
     614    psFree(ordinate);
     615    psFree(mask);
     616    psFree(stats);
     617    return NULL;
     618}
     619
     620// XXX EAM 2024.04.13 : this function seems poorly conceived.  it is replacing the stats from
     621// the reduced overscan vector with the 0-order element from the polynomial fit.  Not clear is
     622// adding any useful information
     623bool pmOverscanUpdateHeader (pmHDU *hdu, pmOverscanStatOptions *overscanOpts, float chi2) {
     624
     625    psString comment = NULL;    // Comment to add
     626
     627    switch (overscanOpts->fitType) {
     628      case PM_FIT_POLY_ORD:
     629      case PM_FIT_POLY_CHEBY: {
     630          psStringAppend(&comment, "Overscan fit (chi2: %.2f): ", chi2);
     631          psPolynomial1D *poly = overscanOpts->poly; // The polynomial
     632          for (int i = 0; i < poly->nX; i++) {
     633              psStringAppend(&comment, "%.1f ", poly->coeff[i]);
     634          }
     635          psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
     636          psFree(comment);
     637          comment = NULL;
     638
     639          // write metadata header value
     640          psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", poly->coeff[0]);
     641          psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", poly->coeffErr[0]);
     642          break;
     643      }
     644      case PM_FIT_SPLINE: {
     645          /*
     646            psSpline1D *spline = overscanOpts->spline; // The spline
     647            for (int i = 0; i < spline->n; i++) {
     648            psStringAppend(&comment, "Overscan fit (chi2: %.2f) %d:", chi2, i);
     649            psPolynomial1D *poly = spline->spline[i]; // i-th polynomial
     650            for (int j = 0; j < poly->nX; j++) {
     651            psStringAppend(&comment, "%.1f ", poly->coeff[i]);
     652            }
     653            psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,
     654            comment, "");
     655            psFree(comment);
     656            comment = NULL;
     657            }
     658          */
     659          // write metadata header value
     660          psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE,
     661                           "Overscan value", NAN);
     662          psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE,
     663                           "Overscan stdev", NAN);
     664          break;
     665      }
     666      case PM_FIT_NONE:
     667        break;
     668      default:
     669        psAbort("Should never get here!!!\n");
     670    }
    725671    return true;
    726 
    727 } // End of overscan subtraction
    728    
     672}
     673
     674// generate stats of overscan vector for header
     675// reduced: 1D vector with overscan stats
     676bool pmOverscanUpdateHeaderVector (pmReadout *input, pmHDU *hdu, pmOverscanOptions *overscanOpts, psVector *reduced) {
     677
     678    psString comment = NULL;    // Comment to add
     679    psStats *vectorStats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
     680    if (!psVectorStats (vectorStats, reduced, NULL, NULL, 0)) {
     681        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     682        return false;
     683    }
     684    psStringAppend(&comment, "Mean Overscan value: %f", vectorStats->sampleMean);
     685    psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
     686    psFree(comment);
     687
     688    // write metadata header value
     689    psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan mean", vectorStats->sampleMean);
     690    psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", vectorStats->sampleStdev);
     691
     692    // EAM 2022.03.29 : if the calculated overscan value is below the threshold,
     693    // declare the readout dead and mask
     694         
     695    if ((vectorStats->sampleMean < overscanOpts->minValid) || (vectorStats->sampleMean > overscanOpts->maxValid)) {
     696        fprintf (stderr, "bad overscan (2) %f, masking readout\n", vectorStats->sampleMean);
     697        psImage *mask = input->mask;
     698        for (int y = 0; y < mask->numRows; y++) {
     699            for (int x = 0; x < mask->numCols; x++) {
     700                mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= overscanOpts->maskVal;
     701            }
     702        }
     703    }
     704
     705    psFree (vectorStats);
     706    return true;
     707}
     708
Note: See TracChangeset for help on using the changeset viewer.