IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 21, 2024, 8:50:38 PM (2 years ago)
Author:
hgao
Message:

implement 2D overscan subtraction per cell as configured in OVERSCAN.2D.SUBSET

Location:
branches/2dbias/ppImage/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/2dbias/ppImage/src/ppImage.h

    r42382 r42679  
    188188
    189189bool ppImageDetrendFree(pmConfig *config, pmFPAview *view);
     190bool ppImageDoPatternForView(bool *doit, const pmConfig *config, const pmChip *chip, const pmFPAview *view, const char *recipeName, const char *recipeValue);
    190191bool ppImageFringeFree(pmConfig *config, pmFPAview *view);
    191192

  • branches/2dbias/ppImage/src/ppImageDetrendPattern.c

    r42382 r42679  
    55#include "ppImage.h"
    66
    7 #define ESCAPE(STATUS,...) {                            \
    8         psError(PS_ERR_UNKNOWN, STATUS, __VA_ARGS__);   \
    9         psFree(view);                                   \
    10         return false;                                   \
    11     }
    12 
    13 static bool doPatternForView (bool *doit, const pmConfig *config, const pmChip *chip, const pmFPAview *view, const char *recipename, const char *recipevalue);
     7#define ESCAPE(STATUS, ...)                           \
     8    {                                                 \
     9        psError(PS_ERR_UNKNOWN, STATUS, __VA_ARGS__); \
     10        psFree(view);                                 \
     11        return false;                                 \
     12    }
     13
     14bool ppImageDoPatternForView(bool *doit, const pmConfig *config, const pmChip *chip, const pmFPAview *view, const char *recipename, const char *recipevalue);
    1415
    1516bool ppImageDetrendPatternRowApply(pmConfig *config, pmChip *chip, const pmFPAview *inputView, ppImageOptions *options);
    … …  
    2627    assert(inputView->readout == -1);
    2728
    28     if (!ppImageDetrendPatternRowApply (config, chip, inputView, options)) {
    29         return false;
    30     }
    31     if (!ppImageDetrendPatternContinuityApply (config, chip, inputView, options)) {
    32         return false;
    33     }
    34     if (!ppImageDetrendPatternCellApply (config, chip, inputView, options)) {
    35         return false;
    36     }
    37     if (!ppImageDetrendPatternDeadCellsApply (config, chip, inputView, options)) {
    38         return false;
    39     }
    40     return(true);
     29    if (!ppImageDetrendPatternRowApply(config, chip, inputView, options))
     30    {
     31        return false;
     32    }
     33    if (!ppImageDetrendPatternContinuityApply(config, chip, inputView, options))
     34    {
     35        return false;
     36    }
     37    if (!ppImageDetrendPatternCellApply(config, chip, inputView, options))
     38    {
     39        return false;
     40    }
     41    if (!ppImageDetrendPatternDeadCellsApply(config, chip, inputView, options))
     42    {
     43        return false;
     44    }
     45    return (true);
    4146}
    4247
    … …  
    4550    bool status;
    4651
    47     if (!options->doPatternRow) return true;
    48    
     52    if (!options->doPatternRow)
     53        return true;
     54
    4955    assert(inputView->chip != -1);
    5056    assert(inputView->cell == -1);
    … …  
    7379
    7480    pmHDU *hdu = pmHDUFromChip(chip);
    75     if (psMetadataLookupBool(&status, hdu->header, "PTRN_ROW")) {
    76         psLogMsg("ppImage", PS_LOG_INFO, "Not performing row pattern correction as it has already been done.");
    77         return true;
     81    if (psMetadataLookupBool(&status, hdu->header, "PTRN_ROW"))
     82    {
     83        psLogMsg("ppImage", PS_LOG_INFO, "Not performing row pattern correction as it has already been done.");
     84        return true;
    7885    }
    7986
    … …  
    8592
    8693    // grab the PATTERN.ROW.AMP file
    87     pmFPAfile *PRAfile = psMetadataLookupPtr (&status, config->files, "PPIMAGE.PATTERN.ROW.AMP");
    88     if (!PRAfile) {
    89         psLogMsg("ppImage", PS_LOG_INFO, "Pattern Row Amplitude file not found, applying to all (with limits from ROW.SUBSET) ");
    90     }
    91        
     94    pmFPAfile *PRAfile = psMetadataLookupPtr(&status, config->files, "PPIMAGE.PATTERN.ROW.AMP");
     95    if (!PRAfile)
     96    {
     97        psLogMsg("ppImage", PS_LOG_INFO, "Pattern Row Amplitude file not found, applying to all (with limits from ROW.SUBSET) ");
     98    }
     99
    92100    pmFPAview *view = pmFPAviewAlloc(0); // View for local processing
    93101    *view = *inputView;
    94102
    95103    pmCell *cell = NULL;
    96     while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) {
    97         // const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME");
    98         // const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME");
    99         // psWarning ("Looping through %s, %s\n", chipName, cellName);
    100 
    101         if (!cell->process || !cell->file_exists) {
    102             continue;
    103         }
    104 
    105         // this forces pmFPAfileRead of the PATTERN.ROW.AMP file (XXX but is this needed?)
    106         if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    107             ESCAPE(false, "load failure for Cell");
    108         }
    109            
    110         if (!cell->data_exists) {
    111             continue;
    112         }
    113 
    114         if (cell->readouts->n > 1) {
    115             psWarning ("Skipping Video Cell for ppImageDetrendPatternApply");
    116             continue;
    117         }
    118 
    119         // grab the corresponding cell
    120         if (PRAfile) {
    121             pmCell *PRAcell = pmFPAviewThisCell (view, PRAfile->fpa);
    122             psAssert (PRAcell, "found Pattern Row Amplitude file, but not cell?");
    123 
    124             // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file)
    125             float amplitude = psMetadataLookupF32 (&status, PRAcell->analysis, "PTN.ROW.AMP");
    126             if (!status) amplitude = NAN;
    127            
    128             // put the value on the science cell
    129             psMetadataAddF32 (cell->analysis, PS_LIST_TAIL, "PTN.ROW.AMP", PS_META_REPLACE, "", amplitude);
    130         }
    131 
    132         bool doPattern = false;
    133         if (!doPatternForView(&doPattern, config, chip, view, RECIPE_NAME, "PATTERN.ROW.SUBSET")) {
    134             ESCAPE(false, "Unable to determine whether row pattern matching should be applied.");
    135         }
    136         if (!doPattern) continue;
    137 
    138         // switch to test threaded version
    139         if (true) {
    140             // I need to allocate a view here to be freed by the
    141             // called function below. 
    142             pmFPAview *myView = pmFPAviewAlloc(0); // View for local processing
    143             *myView = *view;
    144 
    145             // allocate a job, construct the arguments for this job
    146             psThreadJob *job = psThreadJobAlloc("PPIMAGE_PATTERN_ROW_CELL");
    147             psArrayAdd(job->args, 1, config);
    148             psArrayAdd(job->args, 1, input->fpa);
    149             psArrayAdd(job->args, 1, chip);
    150             psArrayAdd(job->args, 1, cell);
    151             psArrayAdd(job->args, 1, myView);
    152             psArrayAdd(job->args, 1, options);
    153             if (!psThreadJobAddPending(job)) {
    154                 return false;
    155             }
    156         } else {
    157             // bump the counter since it must be freed by the function below. 
    158             psMemIncrRefCounter (view); // View for local processing
    159             if (!ppImageDetrendPatternApplyCell (config, input->fpa, chip, cell, view, options)) {
    160                 return false;
    161             }
    162         }
     104    while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL)
     105    {
     106        // const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME");
     107        // const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME");
     108        // psWarning ("Looping through %s, %s\n", chipName, cellName);
     109
     110        if (!cell->process || !cell->file_exists)
     111        {
     112            continue;
     113        }
     114
     115        // this forces pmFPAfileRead of the PATTERN.ROW.AMP file (XXX but is this needed?)
     116        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE))
     117        {
     118            ESCAPE(false, "load failure for Cell");
     119        }
     120
     121        if (!cell->data_exists)
     122        {
     123            continue;
     124        }
     125
     126        if (cell->readouts->n > 1)
     127        {
     128            psWarning("Skipping Video Cell for ppImageDetrendPatternApply");
     129            continue;
     130        }
     131
     132        // grab the corresponding cell
     133        if (PRAfile)
     134        {
     135            pmCell *PRAcell = pmFPAviewThisCell(view, PRAfile->fpa);
     136            psAssert(PRAcell, "found Pattern Row Amplitude file, but not cell?");
     137
     138            // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file)
     139            float amplitude = psMetadataLookupF32(&status, PRAcell->analysis, "PTN.ROW.AMP");
     140            if (!status)
     141                amplitude = NAN;
     142
     143            // put the value on the science cell
     144            psMetadataAddF32(cell->analysis, PS_LIST_TAIL, "PTN.ROW.AMP", PS_META_REPLACE, "", amplitude);
     145        }
     146
     147        bool doPattern = false;
     148        if (!ppImageDoPatternForView(&doPattern, config, chip, view, RECIPE_NAME, "PATTERN.ROW.SUBSET"))
     149        {
     150            ESCAPE(false, "Unable to determine whether row pattern matching should be applied.");
     151        }
     152        if (!doPattern)
     153            continue;
     154
     155        // switch to test threaded version
     156        if (true)
     157        {
     158            // I need to allocate a view here to be freed by the
     159            // called function below.
     160            pmFPAview *myView = pmFPAviewAlloc(0); // View for local processing
     161            *myView = *view;
     162
     163            // allocate a job, construct the arguments for this job
     164            psThreadJob *job = psThreadJobAlloc("PPIMAGE_PATTERN_ROW_CELL");
     165            psArrayAdd(job->args, 1, config);
     166            psArrayAdd(job->args, 1, input->fpa);
     167            psArrayAdd(job->args, 1, chip);
     168            psArrayAdd(job->args, 1, cell);
     169            psArrayAdd(job->args, 1, myView);
     170            psArrayAdd(job->args, 1, options);
     171            if (!psThreadJobAddPending(job))
     172            {
     173                return false;
     174            }
     175        }
     176        else
     177        {
     178            // bump the counter since it must be freed by the function below.
     179            psMemIncrRefCounter(view); // View for local processing
     180            if (!ppImageDetrendPatternApplyCell(config, input->fpa, chip, cell, view, options))
     181            {
     182                return false;
     183            }
     184        }
    163185    }
    164186
    165187    // wait here for the threaded jobs to finish
    166     // if no threads are allocated, this
    167     if (!psThreadPoolWait(true, true)) {
    168         psError(PS_ERR_UNKNOWN, false, "Unable to apply bias correction.");
    169         return false;
    170     }
    171 
    172     psMetadataAddBool(hdu->header, PS_LIST_TAIL, "PTRN_ROW",PS_META_REPLACE,"PATTERN.ROW correction applied",true);
     188    // if no threads are allocated, this
     189    if (!psThreadPoolWait(true, true))
     190    {
     191        psError(PS_ERR_UNKNOWN, false, "Unable to apply bias correction.");
     192        return false;
     193    }
     194
     195    psMetadataAddBool(hdu->header, PS_LIST_TAIL, "PTRN_ROW", PS_META_REPLACE, "PATTERN.ROW correction applied", true);
    173196    psFree(view);
    174    
     197
    175198    return true;
    176199}
    … …  
    180203    bool status;
    181204
    182     if (!options->doPatternContinuity) return true;
     205    if (!options->doPatternContinuity)
     206        return true;
    183207
    184208    assert(inputView->chip != -1);
    … …  
    188212    // see the comment for PATTERN.ROW; the same rules apply for PATTERN.CELL
    189213
    190     int numCells = chip->cells->n;       // Number of cells
     214    int numCells = chip->cells->n;                         // Number of cells
    191215    psVector *tweak = psVectorAlloc(numCells, PS_TYPE_U8); // Tweak cell?
    192     pmFPAview *view = pmFPAviewAlloc(0); // View for local processing
     216    pmFPAview *view = pmFPAviewAlloc(0);                   // View for local processing
    193217    *view = *inputView;
    194218
    195219    pmHDU *hdu = pmHDUFromChip(chip);
    196     if (psMetadataLookupBool(&status, hdu->header, "PTRN_CON")) {
    197         psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell continuity correction as it has already been done.");
    198         return true;
     220    if (psMetadataLookupBool(&status, hdu->header, "PTRN_CON"))
     221    {
     222        psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell continuity correction as it has already been done.");
     223        return true;
    199224    }
    200225
    … …  
    202227    psLogMsg("ppImage", PS_LOG_INFO, "Performing pattern continuity correction for %s\n", chipName);
    203228
    204     for (int i = 0; i < chip->cells->n; i++) {
    205         view->cell = i;
    206 
    207         pmCell *cell = chip->cells->data[i]; // Cell of interest
    208 
    209         if (cell->readouts->n > 1) {
    210             psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell continuity correction on video cell.");
    211             continue;
    212         }
    213 
    214         bool doPattern = false;
    215         if (!doPatternForView(&doPattern, config, chip, view, RECIPE_NAME, "PATTERN.CONTINUITY.SUBSET")) {
    216             ESCAPE(false, "Unable to determine whether row pattern matching should be applied.");
    217         }
    218 
    219         if (doPattern) {
    220             tweak->data.U8[i] = 0xFF;
    221         }
     229    for (int i = 0; i < chip->cells->n; i++)
     230    {
     231        view->cell = i;
     232
     233        pmCell *cell = chip->cells->data[i]; // Cell of interest
     234
     235        if (cell->readouts->n > 1)
     236        {
     237            psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell continuity correction on video cell.");
     238            continue;
     239        }
     240
     241        bool doPattern = false;
     242        if (!ppImageDoPatternForView(&doPattern, config, chip, view, RECIPE_NAME, "PATTERN.CONTINUITY.SUBSET"))
     243        {
     244            ESCAPE(false, "Unable to determine whether row pattern matching should be applied.");
     245        }
     246
     247        if (doPattern)
     248        {
     249            tweak->data.U8[i] = 0xFF;
     250        }
    222251    }
    223252
    224253    // Tweak the cells
    225254    if (!pmPatternContinuity(chip, tweak, options->patternCellBG, options->patternCellMean,
    226                              options->maskValue, options->darkMask,options->patternContinuityEdgeWidth)) {
    227         psFree(tweak);
    228         psFree(view);
    229         return false;
     255                             options->maskValue, options->darkMask, options->patternContinuityEdgeWidth))
     256    {
     257        psFree(tweak);
     258        psFree(view);
     259        return false;
    230260    }
    231261    psFree(tweak);
    232262    psFree(view);
    233263
    234     psMetadataAddBool(hdu->header, PS_LIST_TAIL, "PTRN_CON",PS_META_REPLACE,"PATTERN.CONTINUITY correction applied",true);
     264    psMetadataAddBool(hdu->header, PS_LIST_TAIL, "PTRN_CON", PS_META_REPLACE, "PATTERN.CONTINUITY correction applied", true);
    235265    return true;
    236266}
    … …  
    240270    bool status;
    241271
    242     if (!options->doPatternCell) return true;
     272    if (!options->doPatternCell)
     273        return true;
    243274
    244275    assert(inputView->chip != -1);
    … …  
    246277    assert(inputView->readout == -1);
    247278
    248     int numCells = chip->cells->n;       // Number of cells
     279    int numCells = chip->cells->n;                         // Number of cells
    249280    psVector *tweak = psVectorAlloc(numCells, PS_TYPE_U8); // Tweak cell?
    250     pmFPAview *view = pmFPAviewAlloc(0); // View for local processing
     281    pmFPAview *view = pmFPAviewAlloc(0);                   // View for local processing
    251282    *view = *inputView;
    252283
    253284    pmHDU *hdu = pmHDUFromChip(chip);
    254     if (psMetadataLookupBool(&status, hdu->header, "PTRN_CEL")) {
    255         psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell pattern correction as it has already been done.");
    256         return true;
     285    if (psMetadataLookupBool(&status, hdu->header, "PTRN_CEL"))
     286    {
     287        psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell pattern correction as it has already been done.");
     288        return true;
    257289    }
    258290
    … …  
    260292    psLogMsg("ppImage", PS_LOG_INFO, "Performing cell pattern correction for %s\n", chipName);
    261293
    262     for (int i = 0; i < chip->cells->n; i++) {
    263         view->cell = i;
    264 
    265         pmCell *cell = chip->cells->data[i]; // Cell of interest
    266 
    267         if (cell->readouts->n > 1) {
    268             psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell pattern correction on video cell.");
    269             continue;
    270         }
    271 
    272         bool doPattern = false;
    273         if (!doPatternForView(&doPattern, config, chip, view, RECIPE_NAME, "PATTERN.CELL.SUBSET")) {
    274             ESCAPE(false, "Unable to determine whether row pattern matching should be applied.");
    275         }
    276         if (doPattern) {
    277             tweak->data.U8[i] = 0xFF;
    278         }
     294    for (int i = 0; i < chip->cells->n; i++)
     295    {
     296        view->cell = i;
     297
     298        pmCell *cell = chip->cells->data[i]; // Cell of interest
     299
     300        if (cell->readouts->n > 1)
     301        {
     302            psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell pattern correction on video cell.");
     303            continue;
     304        }
     305
     306        bool doPattern = false;
     307        if (!ppImageDoPatternForView(&doPattern, config, chip, view, RECIPE_NAME, "PATTERN.CELL.SUBSET"))
     308        {
     309            ESCAPE(false, "Unable to determine whether row pattern matching should be applied.");
     310        }
     311        if (doPattern)
     312        {
     313            tweak->data.U8[i] = 0xFF;
     314        }
    279315    }
    280316
    281317    // Tweak the cells
    282318    if (!pmPatternCell(chip, tweak, options->patternCellBG, options->patternCellMean,
    283                        options->maskValue, options->darkMask)) {
    284         psFree(tweak);
    285         psFree(view);
    286         return false;
     319                       options->maskValue, options->darkMask))
     320    {
     321        psFree(tweak);
     322        psFree(view);
     323        return false;
    287324    }
    288325    psFree(tweak);
    289326    psFree(view);
    290327
    291     psMetadataAddBool(hdu->header, PS_LIST_TAIL, "PTRN_CEL",PS_META_REPLACE,"PATTERN.CELL correction applied",true);
     328    psMetadataAddBool(hdu->header, PS_LIST_TAIL, "PTRN_CEL", PS_META_REPLACE, "PATTERN.CELL correction applied", true);
    292329    return true;
    293330}
    294331
    295 bool ppImageDetrendPatternDeadCellsMask (pmChip *chip, psImageMaskType maskVal) {
    296 
    297     int numCells = chip->cells->n;       // Number of cells
     332bool ppImageDetrendPatternDeadCellsMask(pmChip *chip, psImageMaskType maskVal)
     333{
     334
     335    int numCells = chip->cells->n; // Number of cells
    298336
    299337    // now mask bad cells
    300     for (int i = 0; i < numCells; i++) {
    301         pmCell *cell = chip->cells->data[i]; // Cell of interest
     338    for (int i = 0; i < numCells; i++)
     339    {
     340        pmCell *cell = chip->cells->data[i];     // Cell of interest
    302341        pmReadout *ro = cell->readouts->data[0]; // Readout of interest
    303342
    304         psImage *mask = ro->mask; // mask of interest
    305         int numCols = mask->numCols, numRows = mask->numRows; // Size of image
    306        
    307         for (int y = 0; y < numRows; y++) {
    308             for (int x = 0; x < numCols; x++) {
    309                 mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskVal;
    310             }
    311         }
     343        psImage *mask = ro->mask;                             // mask of interest
     344        int numCols = mask->numCols, numRows = mask->numRows; // Size of image
     345
     346        for (int y = 0; y < numRows; y++)
     347        {
     348            for (int x = 0; x < numCols; x++)
     349            {
     350                mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskVal;
     351            }
     352        }
    312353    }
    313354    return true;
    314355}
    315356
    316 bool ppImageDetrendPatternDeadCellsXY42 (pmChip *chip, psImageMaskType maskVal) {
     357bool ppImageDetrendPatternDeadCellsXY42(pmChip *chip, psImageMaskType maskVal)
     358{
    317359
    318360    bool status;
    319361
    320362    // extract the MEDIAN_CELL_BACKGROUND values to check for problems
    321     int numCells = chip->cells->n;       // Number of cells
    322     psVector *cellBackground = psVectorAllocEmpty (numCells, PS_DATA_F32);
    323 
    324     for (int i = 0; i < numCells; i++) {
    325         pmCell *cell = chip->cells->data[i];;
    326 
    327         if (!cell->process || !cell->file_exists || !cell->data_exists) continue;
    328        
    329         // select the HDU for this cell
    330         pmHDU *cellHDU = pmHDUFromCell(cell);  // HDU of interest
    331        
    332         psF32 BackValue = psMetadataLookupF32(&status, cellHDU->header, "BACK_VAL");
    333         if (!status) continue;
    334 
    335         psF32 BackError = psMetadataLookupF32(&status, cellHDU->header, "BACK_ERR");
    336         if (!status) continue;
    337 
    338         if (BackValue == 0.0) continue;
    339         if (BackError == 0.0) continue;
    340 
    341         psVectorAppend (cellBackground, BackValue / BackError);
    342     }
    343 
    344     if (cellBackground->n < 0.375*numCells) {
    345         // Chip is bad, mask the whole thing
    346         ppImageDetrendPatternDeadCellsMask (chip, maskVal);
    347         psFree (cellBackground);
    348         return true;
     363    int numCells = chip->cells->n; // Number of cells
     364    psVector *cellBackground = psVectorAllocEmpty(numCells, PS_DATA_F32);
     365
     366    for (int i = 0; i < numCells; i++)
     367    {
     368        pmCell *cell = chip->cells->data[i];
     369        ;
     370
     371        if (!cell->process || !cell->file_exists || !cell->data_exists)
     372            continue;
     373
     374        // select the HDU for this cell
     375        pmHDU *cellHDU = pmHDUFromCell(cell); // HDU of interest
     376
     377        psF32 BackValue = psMetadataLookupF32(&status, cellHDU->header, "BACK_VAL");
     378        if (!status)
     379            continue;
     380
     381        psF32 BackError = psMetadataLookupF32(&status, cellHDU->header, "BACK_ERR");
     382        if (!status)
     383            continue;
     384
     385        if (BackValue == 0.0)
     386            continue;
     387        if (BackError == 0.0)
     388            continue;
     389
     390        psVectorAppend(cellBackground, BackValue / BackError);
     391    }
     392
     393    if (cellBackground->n < 0.375 * numCells)
     394    {
     395        // Chip is bad, mask the whole thing
     396        ppImageDetrendPatternDeadCellsMask(chip, maskVal);
     397        psFree(cellBackground);
     398        return true;
    349399    }
    350400
    351401    // Second, calculate the median
    352     psVectorSortInPlace (cellBackground);
     402    psVectorSortInPlace(cellBackground);
    353403    int midPt = cellBackground->n / 2.0;
    354     float median = cellBackground->n % 2 ? cellBackground->data.F32[midPt] : 0.5*(cellBackground->data.F32[midPt] + cellBackground->data.F32[midPt-1]);
    355     psFree (cellBackground);
     404    float median = cellBackground->n % 2 ? cellBackground->data.F32[midPt] : 0.5 * (cellBackground->data.F32[midPt] + cellBackground->data.F32[midPt - 1]);
     405    psFree(cellBackground);
    356406
    357407    // XXX hardwired number
    358     if (median < 3.0) {
    359         // Chip is bad, mask the whole thing
    360         ppImageDetrendPatternDeadCellsMask (chip, maskVal);
    361     }
    362     psFree (cellBackground);
     408    if (median < 3.0)
     409    {
     410        // Chip is bad, mask the whole thing
     411        ppImageDetrendPatternDeadCellsMask(chip, maskVal);
     412    }
     413    psFree(cellBackground);
    363414    return true;
    364415}
    365416
    366417// to apply the dead cell pattern, we need to transfer the pattern for this chip from the
    367 // pmFPAfile for the pattern to the one for the image being processed. 
     418// pmFPAfile for the pattern to the one for the image being processed.
    368419bool ppImageDetrendPatternDeadCellsApply(pmConfig *config, pmChip *chip, const pmFPAview *inputView, ppImageOptions *options)
    369420{
    370421    bool status;
    371422
    372     if (!options->doPatternDeadCells) return true;
     423    if (!options->doPatternDeadCells)
     424        return true;
    373425
    374426    assert(inputView->chip != -1);
    … …  
    380432
    381433    pmHDU *hdu = pmHDUFromChip(chip);
    382     if (psMetadataLookupBool(&status, hdu->header, "PTRN_DED")) {
    383         psLogMsg("ppImage", PS_LOG_INFO, "Not performing dead cell pattern correction as it has already been done.");
    384         return true;
    385     }
    386 
    387     if (!strcmp (chipName, "XY40") || !strcmp (chipName, "XY42")) {
    388         // special case : check for BACK_VAL / BACK_ERR > 3 or < 3
    389         psLogMsg("ppImage", PS_LOG_INFO, "Using special case for XY40 and XY42");
    390         ppImageDetrendPatternDeadCellsXY42 (chip, options->blankMask);
    391         return true;
     434    if (psMetadataLookupBool(&status, hdu->header, "PTRN_DED"))
     435    {
     436        psLogMsg("ppImage", PS_LOG_INFO, "Not performing dead cell pattern correction as it has already been done.");
     437        return true;
     438    }
     439
     440    if (!strcmp(chipName, "XY40") || !strcmp(chipName, "XY42"))
     441    {
     442        // special case : check for BACK_VAL / BACK_ERR > 3 or < 3
     443        psLogMsg("ppImage", PS_LOG_INFO, "Using special case for XY40 and XY42");
     444        ppImageDetrendPatternDeadCellsXY42(chip, options->blankMask);
     445        return true;
    392446    }
    393447
    394448    pmFPAfile *pattern = psMetadataLookupPtr(&status, config->files, "PPIMAGE.PATTERN.DEAD.CELLS");
    395     if (!pattern) {
    396         psLogMsg("ppImage", PS_LOG_INFO, "Pattern Dead Cells file not found, skipping");
    397         return true;
     449    if (!pattern)
     450    {
     451        psLogMsg("ppImage", PS_LOG_INFO, "Pattern Dead Cells file not found, skipping");
     452        return true;
    398453    }
    399454
    … …  
    401456    *view = *inputView;
    402457
    403     pmChip *patternChip = pmFPAviewThisChip (view, pattern->fpa);
     458    pmChip *patternChip = pmFPAviewThisChip(view, pattern->fpa);
    404459
    405460    // grab the pattern from the input pattern file chip
    406     psImage *deadCellPattern = (psImage *) psMetadataLookupPtr (&status, patternChip->analysis, "PTN.DEAD.CELL");
    407     if (!deadCellPattern) {
     461    psImage *deadCellPattern = (psImage *)psMetadataLookupPtr(&status, patternChip->analysis, "PTN.DEAD.CELL");
     462    if (!deadCellPattern)
     463    {
    408464        psLogMsg("psModules.detrend", PS_LOG_DETAIL, "No DEAD CELL pattern for chip, skipping\n");
    409         psFree (view);
    410         return true;
     465        psFree(view);
     466        return true;
    411467    }
    412468
    413469    // copy the pattern pointer to the chip of the image being processed
    414     psMetadataAddImage (chip->analysis, PS_LIST_TAIL, "PTN.DEAD.CELL", PS_META_REPLACE, "", deadCellPattern);
     470    psMetadataAddImage(chip->analysis, PS_LIST_TAIL, "PTN.DEAD.CELL", PS_META_REPLACE, "", deadCellPattern);
    415471
    416472    // extract the MEDIAN_CELL_BACKGROUND values to check for problems
    417     int numCells = chip->cells->n;       // Number of cells
    418     psVector *cellBackground = psVectorAllocEmpty (numCells, PS_DATA_F32);
     473    int numCells = chip->cells->n; // Number of cells
     474    psVector *cellBackground = psVectorAllocEmpty(numCells, PS_DATA_F32);
    419475
    420476    pmCell *cell = NULL;
    421     while ((cell = pmFPAviewNextCell(view, chip->parent, 1)) != NULL) {
    422         if (!cell->process || !cell->file_exists || !cell->data_exists) {
    423             psVectorAppend (cellBackground, NAN);
    424             continue;
    425         }
    426        
    427         // select the HDU for this cell
    428         pmHDU *cellHDU = pmHDUFromCell(cell);  // HDU of interest
    429        
    430         psF32 value = psMetadataLookupF32(&status, cellHDU->header, "BACK_VAL");
    431         if (!status) {
    432             psVectorAppend (cellBackground, NAN);
    433             continue;
    434         }
    435         psVectorAppend (cellBackground, value);
     477    while ((cell = pmFPAviewNextCell(view, chip->parent, 1)) != NULL)
     478    {
     479        if (!cell->process || !cell->file_exists || !cell->data_exists)
     480        {
     481            psVectorAppend(cellBackground, NAN);
     482            continue;
     483        }
     484
     485        // select the HDU for this cell
     486        pmHDU *cellHDU = pmHDUFromCell(cell); // HDU of interest
     487
     488        psF32 value = psMetadataLookupF32(&status, cellHDU->header, "BACK_VAL");
     489        if (!status)
     490        {
     491            psVectorAppend(cellBackground, NAN);
     492            continue;
     493        }
     494        psVectorAppend(cellBackground, value);
    436495    }
    437496
    438497    // match cellBackground pattern to registered patterns and mask as needed
    439     if (!pmPatternDeadCells(chip, cellBackground, options->blankMask)) {
    440         psFree(cellBackground);
    441         psFree(view);
    442         return false;
     498    if (!pmPatternDeadCells(chip, cellBackground, options->blankMask))
     499    {
     500        psFree(cellBackground);
     501        psFree(view);
     502        return false;
    443503    }
    444504    psFree(cellBackground);
    … …  
    449509}
    450510
    451 static bool doPatternForView (bool *doit, const pmConfig *config, const pmChip *chip, const pmFPAview *view, const char *recipeName, const char *recipeValue) {
     511bool ppImageDoPatternForView(bool *doit, const pmConfig *config, const pmChip *chip, const pmFPAview *view, const char *recipeName, const char *recipeValue)
     512{
    452513
    453514    *doit = false;
    … …  
    456517
    457518    doPattern = pmConfigRecipeValueByView(config, recipeName, recipeValue, chip->parent, view);
    458    
    459     if (!doPattern) {
     519
     520    if (!doPattern)
     521    {
    460522        psError(PS_ERR_UNKNOWN, false, "Unable to determine whether row pattern matching should be applied.");
    461523        return false;
    462524    }
    463     if (doPattern->type == PS_DATA_BOOL) {
     525    if (doPattern->type == PS_DATA_BOOL)
     526    {
    464527        *doit = doPattern->data.B;
    465528        return true;
    466529    }
    467     if (doPattern->type == PS_DATA_STRING) {
     530    if (doPattern->type == PS_DATA_STRING)
     531    {
    468532        // expect a string of the form "000110001001" with at least view->cell entries
    469533        char *string = doPattern->data.str;
    470         if (strlen(string) < view->cell) {
    471             psError(PS_ERR_UNKNOWN, true, "error in PATTERN.ROW.SUBSET chip string (too few elements %d)", (int) strlen(string));
     534        if (strlen(string) < view->cell)
     535        {
     536            psError(PS_ERR_UNKNOWN, true, "error in PATTERN.ROW.SUBSET chip string (too few elements %d)", (int)strlen(string));
    472537            return false;
    473538        }
    474         switch (string[view->cell]) {
    475           case '0':
    476           case 'f':
    477           case 'F':
    478           case 'n':
    479           case 'N':
     539        switch (string[view->cell])
     540        {
     541        case '0':
     542        case 'f':
     543        case 'F':
     544        case 'n':
     545        case 'N':
    480546            *doit = false;
    481547            return true;
    482           case '1':
    483           case 't':
    484           case 'T':
    485           case 'y':
    486           case 'Y':
     548        case '1':
     549        case 't':
     550        case 'T':
     551        case 'y':
     552        case 'Y':
    487553            *doit = true;
    488554            return true;
    489           default:
     555        default:
    490556            psError(PS_ERR_UNKNOWN, true, "error in PATTERN.ROW.SUBSET chip string %s (unknown value %c))", string, string[view->cell]);
    491557            return false;
    … …  
    497563}
    498564
    499 // thread safety :
    500 bool ppImageDetrendPatternApplyCell (pmConfig *config, pmFPA *fpa, pmChip *chip, pmCell *cell, pmFPAview *view, ppImageOptions *options) {
    501  
     565// thread safety :
     566bool ppImageDetrendPatternApplyCell(pmConfig *config, pmFPA *fpa, pmChip *chip, pmCell *cell, pmFPAview *view, ppImageOptions *options)
     567{
     568
    502569    // process each of the readouts
    503     pmReadout *readout;         // Readout from cell
    504     while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
    505         if (!readout->data_exists) {
    506             continue;
    507         }
    508 
    509         // Perform pattern correction
    510         if (!pmPatternRow(readout, options->patternRowOrder, options->patternRowIter,
    511                           options->patternRowRej, options->patternRowThresh, options->patternRowMean,
    512                           options->patternRowStdev, options->maskValue, options->darkMask)) {
    513             psFree (view);
    514             return false;
    515         }
    516     }
    517     psFree (view); // supplied view must be freeable
     570    pmReadout *readout; // Readout from cell
     571    while ((readout = pmFPAviewNextReadout(view, fpa, 1)) != NULL)
     572    {
     573        if (!readout->data_exists)
     574        {
     575            continue;
     576        }
     577
     578        // Perform pattern correction
     579        if (!pmPatternRow(readout, options->patternRowOrder, options->patternRowIter,
     580                          options->patternRowRej, options->patternRowThresh, options->patternRowMean,
     581                          options->patternRowStdev, options->maskValue, options->darkMask))
     582        {
     583            psFree(view);
     584            return false;
     585        }
     586    }
     587    psFree(view); // supplied view must be freeable
    518588    return true;
    519589}
    520590
    521 psVector *ppImageDetrendPatternCellFailures(pmFPAfile *input, const pmFPAview *inputView) {
    522    
     591psVector *ppImageDetrendPatternCellFailures(pmFPAfile *input, const pmFPAview *inputView)
     592{
     593
    523594    bool status = false;
    524595    pmCell *cell = NULL;
    525    
     596
    526597    pmFPAview *view = pmFPAviewAlloc(0); // View for local processing
    527598    *view = *inputView;
    528599
    529600    // extract the MEDIAN_CELL_BACKGROUND values to check for problems
    530     psVector *cellBackground = psVectorAllocEmpty (64, PS_DATA_F32);
    531     while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) {
    532         if (!cell->process || !cell->file_exists || !cell->data_exists) {
    533             psVectorAppend (cellBackground, NAN);
    534             continue;
    535         }
    536        
    537         // select the HDU for this cell
    538         pmHDU *hdu = pmHDUFromCell(cell);  // HDU of interest
    539        
    540         psF32 value = psMetadataLookupF32(&status, hdu->header, "BACK_VAL");
    541         psVectorAppend (cellBackground, value);
     601    psVector *cellBackground = psVectorAllocEmpty(64, PS_DATA_F32);
     602    while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL)
     603    {
     604        if (!cell->process || !cell->file_exists || !cell->data_exists)
     605        {
     606            psVectorAppend(cellBackground, NAN);
     607            continue;
     608        }
     609
     610        // select the HDU for this cell
     611        pmHDU *hdu = pmHDUFromCell(cell); // HDU of interest
     612
     613        psF32 value = psMetadataLookupF32(&status, hdu->header, "BACK_VAL");
     614        psVectorAppend(cellBackground, value);
    542615    }
    543616    return cellBackground;

  • branches/2dbias/ppImage/src/ppImageDetrendReadout.c

    r42382 r42679  
    55#include "ppImage.h"
    66
     7#define ESCAPE(STATUS, ...)                           \
     8    {                                                 \
     9        psError(PS_ERR_UNKNOWN, STATUS, __VA_ARGS__); \
     10        psFree(view);                                 \
     11        return false;                                 \
     12    }
     13
    714bool ppImageDetrendReadout(pmConfig *config, ppImageOptions *options, pmFPAview *view)
    815{
    9   // psTimerStart("detrend.readout");
     16    // psTimerStart("detrend.readout");
    1017
    1118    // construct a view for the detrend images (which have only one readout)
    … …  
    1825
    1926    // Check that the gain is set (this is used by both pmReadoutGenerateMask and pmReadoutGenerateVariance)
    20     {
    21       float gain = psMetadataLookupF32(NULL, input->parent->concepts, "CELL.GAIN"); // Gain for cell
    22       if (!isfinite(gain)) {
    23         psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, RECIPE_NAME); // Recipe
    24         psAssert(recipe, "Should be there!");
    25         bool override = psMetadataLookupBool(NULL, recipe, "GAIN.OVERRIDE"); // Override the bad gain?
    26         if (override) {
    27           psWarning("CELL.GAIN is not set for readout (%d,%d,%d) --- setting to unity.", view->chip, view->cell, view->readout);
    28           psMetadataItem *item = psMetadataLookup(input->parent->concepts, "CELL.GAIN"); // Gain item
    29           psAssert(item, "Should be there!");
    30           item->data.F32 = 1.0;
    31 
    32           // for unity gain, there is no modification for the readnoise, note that it has (effectively) been updated
    33           psMetadataRemoveKey(input->parent->concepts, "CELL.READNOISE.UPDATE");
    34         } else {
    35           psWarning("CELL.GAIN is NAN for readout (%d,%d,%d), image will be masked.", view->chip, view->cell, view->readout);
    36         }
    37       }
     27    {
     28        float gain = psMetadataLookupF32(NULL, input->parent->concepts, "CELL.GAIN"); // Gain for cell
     29        if (!isfinite(gain))
     30        {
     31            psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, RECIPE_NAME); // Recipe
     32            psAssert(recipe, "Should be there!");
     33            bool override = psMetadataLookupBool(NULL, recipe, "GAIN.OVERRIDE"); // Override the bad gain?
     34            if (override)
     35            {
     36                psWarning("CELL.GAIN is not set for readout (%d,%d,%d) --- setting to unity.", view->chip, view->cell, view->readout);
     37                psMetadataItem *item = psMetadataLookup(input->parent->concepts, "CELL.GAIN"); // Gain item
     38                psAssert(item, "Should be there!");
     39                item->data.F32 = 1.0;
     40
     41                // for unity gain, there is no modification for the readnoise, note that it has (effectively) been updated
     42                psMetadataRemoveKey(input->parent->concepts, "CELL.READNOISE.UPDATE");
     43            }
     44            else
     45            {
     46                psWarning("CELL.GAIN is NAN for readout (%d,%d,%d), image will be masked.", view->chip, view->cell, view->readout);
     47            }
     48        }
    3849    }
    3950    // psLogMsg ("ppImage", 6, "check gain: %f sec\n", psTimerMark ("detrend.readout"));
    … …  
    4354    bool hasVideo = false;
    4455    {
    45       // XXX test for GPC1? or CAN_HAVE_VIDEO in camera config?
    46       bool status = false;
    47       if (!input) goto done_video_check;
    48       if (!input->parent) goto done_video_check;
    49       if (!input->parent->parent) goto done_video_check;
    50       if (!input->parent->parent->hdu) goto done_video_check;
    51       if (!input->parent->parent->hdu->header) goto done_video_check;
    52       char *ptr = psMetadataLookupStr(&status,input->parent->parent->hdu->header,"CELLMODE");
    53       if (status) {
    54         psLogMsg ("ppImage.detrend", PS_LOG_DETAIL, "VIDEO: %d %d %d\n",(int) options->hasVideo,(int) options->useVideoDark, (int) options->useVideoMask);
    55         char *Vptr = strchr(ptr,'V');
    56         if (Vptr) {
    57           hasVideo = options->hasVideo = true;
    58           psLogMsg ("ppImage.detrend", PS_LOG_INFO, "VIDEO: %d %d %d\n",(int) options->hasVideo,(int) options->useVideoDark, (int) options->useVideoMask);
    59         }
    60       }
    61     }
    62     done_video_check:
     56        // XXX test for GPC1? or CAN_HAVE_VIDEO in camera config?
     57        bool status = false;
     58        if (!input)
     59            goto done_video_check;
     60        if (!input->parent)
     61            goto done_video_check;
     62        if (!input->parent->parent)
     63            goto done_video_check;
     64        if (!input->parent->parent->hdu)
     65            goto done_video_check;
     66        if (!input->parent->parent->hdu->header)
     67            goto done_video_check;
     68        char *ptr = psMetadataLookupStr(&status, input->parent->parent->hdu->header, "CELLMODE");
     69        if (status)
     70        {
     71            psLogMsg("ppImage.detrend", PS_LOG_DETAIL, "VIDEO: %d %d %d\n", (int)options->hasVideo, (int)options->useVideoDark, (int)options->useVideoMask);
     72            char *Vptr = strchr(ptr, 'V');
     73            if (Vptr)
     74            {
     75                hasVideo = options->hasVideo = true;
     76                psLogMsg("ppImage.detrend", PS_LOG_INFO, "VIDEO: %d %d %d\n", (int)options->hasVideo, (int)options->useVideoDark, (int)options->useVideoMask);
     77            }
     78        }
     79    }
     80done_video_check:
    6381    // psLogMsg ("ppImage", 6, "check video: %f sec\n", psTimerMark ("detrend.readout"));
    6482
    … …  
    7492
    7593    // Masking on the basis of pixel value needs to be done before anything else, so the values are pristine.
    76     if (options->doMaskBuild) {
     94    if (options->doMaskBuild)
     95    {
    7796        psImageMaskType satMask = options->doMaskSat ? options->satMask : 0;
    7897        psImageMaskType lowMask = options->doMaskLow ? options->lowMask : 0;
    79         if (!pmReadoutGenerateMask(input, satMask, lowMask)) {
    80           psError(PS_ERR_UNKNOWN, false, "Unable to generate a mask.");
    81           psFree(detview);
    82           return false;
    83         }
    84         // psLogMsg ("ppImage", 6, "generate mask: %f sec\n", psTimerMark ("detrend.readout"));
     98        if (!pmReadoutGenerateMask(input, satMask, lowMask))
     99        {
     100            psError(PS_ERR_UNKNOWN, false, "Unable to generate a mask.");
     101            psFree(detview);
     102            return false;
     103        }
     104        // psLogMsg ("ppImage", 6, "generate mask: %f sec\n", psTimerMark ("detrend.readout"));
    85105    }
    86106    // apply the externally supplied mask to the input->mask pixels
    87     if (options->doMask) {
    88       pmReadout *mask;
    89       if ((options->useVideoMask)&&(hasVideo)) {
    90         mask = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.VIDEOMASK");
    91       }
    92       else {
    93         mask = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.MASK");
    94       }
    95       pmMaskBadPixels(input, mask, options->maskValue);
    96       // psLogMsg ("ppImage", 6, "apply mask: %f sec\n", psTimerMark ("detrend.readout"));
     107    if (options->doMask)
     108    {
     109        pmReadout *mask;
     110        if ((options->useVideoMask) && (hasVideo))
     111        {
     112            mask = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.VIDEOMASK");
     113        }
     114        else
     115        {
     116            mask = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.MASK");
     117        }
     118        pmMaskBadPixels(input, mask, options->maskValue);
     119        // psLogMsg ("ppImage", 6, "apply mask: %f sec\n", psTimerMark ("detrend.readout"));
    97120    }
    98121
    … …  
    101124        // extern bool ppImageBurntoolApply(pmConfig *, ppImageOptions *, pmFPAview *, pmReadout *);
    102125        ppImageBurntoolApply(config, options, view, input);
    103         // psLogMsg ("ppImage", 6, "apply burntool: %f sec\n", psTimerMark ("detrend.readout"));
    104     }
    105    
    106     if (options->doMaskBurntool) {
    107         if (options->doApplyBurntool) {
     126        // psLogMsg ("ppImage", 6, "apply burntool: %f sec\n", psTimerMark ("detrend.readout"));
     127    }
     128
     129    if (options->doMaskBurntool)
     130    {
     131        if (options->doApplyBurntool)
     132        {
    108133            // build burntool mask from data input burntool table
    109             ppImageBurntoolMaskFromTable(config,options,view,input);
    110         } else {
     134            ppImageBurntoolMaskFromTable(config, options, view, input);
     135        }
     136        else
     137        {
    111138            // build burntool mask from data in the input image's fits extension
    112             ppImageBurntoolMask(config,options,view,input);
    113         }
    114         // psLogMsg ("ppImage", 6, "apply burntool mask: %f sec\n", psTimerMark ("detrend.readout"));
     139            ppImageBurntoolMask(config, options, view, input);
     140        }
     141        // psLogMsg ("ppImage", 6, "apply burntool mask: %f sec\n", psTimerMark ("detrend.readout"));
    115142    }
    116143
    117144    // Subtract the overscan
    118     if (options->doOverscan) {
    119       if (!pmOverscanSubtract (input, options->overscan)) {
    120         psError(PS_ERR_UNKNOWN, false, "Unable to subtract overscan.");
    121         psFree(detview);
    122         return false;
    123       }
    124       // psLogMsg ("ppImage", 6, "subtract overscan: %f sec\n", psTimerMark ("detrend.readout"));
     145    if (options->doOverscan)
     146    {
     147        bool doTwoDOverscan = false;
     148        if (options->overscan->TwoD)
     149        {
     150            pmChip *chip = pmFPAfileThisChip(config->files, view, "PPIMAGE.INPUT");
     151            if (!ppImageDoPatternForView(&doTwoDOverscan, config, chip, view, RECIPE_NAME, "OVERSCAN.2D.SUBSET"))
     152            {
     153                ESCAPE(false, "Unable to determine whether 2D Overscan subtraction should be applied.");
     154            }
     155        }
     156        if (!pmOverscanSubtract(input, options->overscan, doTwoDOverscan))
     157        {
     158            psError(PS_ERR_UNKNOWN, false, "Unable to subtract overscan.");
     159            psFree(detview);
     160            return false;
     161        }
     162        // psLogMsg ("ppImage", 6, "subtract overscan: %f sec\n", psTimerMark ("detrend.readout"));
    125163    }
    126164
    127165    // measure the overscan-subtracted readoutBackground here (or subtract the overscan value?)
    128166    // XXX this is the measurements and should be independent of the pattern masking
    129     if (options->doPatternDeadCells) {
    130         psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    131         psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
    132         psImageBackground (stats, NULL, input->image, NULL, 0xffff, rng);
    133 
    134         // save this value somewhere
    135         pmHDU *hdu = pmHDUFromReadout(input);  // HDU of interest
    136         psMetadataAddF32(hdu->header, PS_LIST_TAIL, "BACK_VAL", PS_META_REPLACE, "Median cell background", stats->robustMedian);
    137         psMetadataAddF32(hdu->header, PS_LIST_TAIL, "BACK_ERR", PS_META_REPLACE, "Stdev of cell background", stats->robustStdev);
    138         psFree (stats);
    139         psFree (rng);
     167    if (options->doPatternDeadCells)
     168    {
     169        psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     170        psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
     171        psImageBackground(stats, NULL, input->image, NULL, 0xffff, rng);
     172
     173        // save this value somewhere
     174        pmHDU *hdu = pmHDUFromReadout(input); // HDU of interest
     175        psMetadataAddF32(hdu->header, PS_LIST_TAIL, "BACK_VAL", PS_META_REPLACE, "Median cell background", stats->robustMedian);
     176        psMetadataAddF32(hdu->header, PS_LIST_TAIL, "BACK_ERR", PS_META_REPLACE, "Stdev of cell background", stats->robustStdev);
     177        psFree(stats);
     178        psFree(rng);
    140179    }
    141180
    142181    // Non-linearity correction
    143     if (options->doNonLin) {
    144       if (!ppImageDetrendNonLinear(input,detview,config)) {
    145         psError(PS_ERR_UNKNOWN, false, "Unable to correct NonLinearity");
    146         psFree(detview);
    147         return(false);
    148       }
    149       // psLogMsg ("ppImage", 6, "nonlinear correction: %f sec\n", psTimerMark ("detrend.readout"));
     182    if (options->doNonLin)
     183    {
     184        if (!ppImageDetrendNonLinear(input, detview, config))
     185        {
     186            psError(PS_ERR_UNKNOWN, false, "Unable to correct NonLinearity");
     187            psFree(detview);
     188            return (false);
     189        }
     190        // psLogMsg ("ppImage", 6, "nonlinear correction: %f sec\n", psTimerMark ("detrend.readout"));
    150191    }
    151192    // New Non-linearity correction (exclusive of the above)
    152     if (options->doNewNonLin) {
    153       if (!ppImageDetrendNewNonLinear(input, detview, config)) {
    154         psError(PS_ERR_UNKNOWN, false, "Unable to correct Non-Linearity with new version (2023)");
    155         psFree(detview);
    156         return(false);
    157       }
    158       // psLogMsg ("ppImage", 6, "nonlinear correction: %f sec\n", psTimerMark ("detrend.readout"));
     193    if (options->doNewNonLin)
     194    {
     195        if (!ppImageDetrendNewNonLinear(input, detview, config))
     196        {
     197            psError(PS_ERR_UNKNOWN, false, "Unable to correct Non-Linearity with new version (2023)");
     198            psFree(detview);
     199            return (false);
     200        }
     201        // psLogMsg ("ppImage", 6, "nonlinear correction: %f sec\n", psTimerMark ("detrend.readout"));
    159202    }
    160203
    161204    // set up the dark and bias
    162     pmCell *dark = NULL;                // Multi-dark
    163     pmReadout *oldDark = NULL;          // Old-fashioned dark
    164     pmReadout  *bias = NULL;
    165     if (options->doBias) {
     205    pmCell *dark = NULL;       // Multi-dark
     206    pmReadout *oldDark = NULL; // Old-fashioned dark
     207    pmReadout *bias = NULL;
     208    if (options->doBias)
     209    {
    166210        bias = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.BIAS");
    167211    }
    168     if (options->doDark) {
    169         bool mdok;                      // Status of MD lookup
    170         psMetadata *recipe = psMetadataLookupPtr (&mdok, config->recipes, RECIPE_NAME);
     212    if (options->doDark)
     213    {
     214        bool mdok; // Status of MD lookup
     215        psMetadata *recipe = psMetadataLookupPtr(&mdok, config->recipes, RECIPE_NAME);
    171216        assert(mdok && recipe);
    172217
    173         if ((options->useVideoDark)&&(hasVideo)) {
    174           if (psMetadataLookupBool(&mdok, recipe, "OLDDARK")) {
    175             oldDark = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.VIDEODARK");
    176           } else {
    177             dark = pmFPAfileThisCell(config->files, detview, "PPIMAGE.VIDEODARK");
    178           }
    179         }
    180         else {
    181           if (psMetadataLookupBool(&mdok, recipe, "OLDDARK")) {
    182             oldDark = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.DARK");
    183           } else {
    184             dark = pmFPAfileThisCell(config->files, detview, "PPIMAGE.DARK");
    185           }
    186         }
     218        if ((options->useVideoDark) && (hasVideo))
     219        {
     220            if (psMetadataLookupBool(&mdok, recipe, "OLDDARK"))
     221            {
     222                oldDark = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.VIDEODARK");
     223            }
     224            else
     225            {
     226                dark = pmFPAfileThisCell(config->files, detview, "PPIMAGE.VIDEODARK");
     227            }
     228        }
     229        else
     230        {
     231            if (psMetadataLookupBool(&mdok, recipe, "OLDDARK"))
     232            {
     233                oldDark = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.DARK");
     234            }
     235            else
     236            {
     237                dark = pmFPAfileThisCell(config->files, detview, "PPIMAGE.DARK");
     238            }
     239        }
    187240    }
    188241
    189242    // Bias and temperature-independent-dark subtraction are merged.
    190     if (options->doBias) {
    191         if (!pmBiasSubtract(input, bias, oldDark, view)) {
     243    if (options->doBias)
     244    {
     245        if (!pmBiasSubtract(input, bias, oldDark, view))
     246        {
    192247            psError(PS_ERR_UNKNOWN, false, "Unable to subtract bias.");
    193248            psFree(detview);
    194249            return false;
    195250        }
    196         // psLogMsg ("ppImage", 6, "apply bias: %f sec\n", psTimerMark ("detrend.readout"));
    197     }
    198    
     251        // psLogMsg ("ppImage", 6, "apply bias: %f sec\n", psTimerMark ("detrend.readout"));
     252    }
     253
    199254    // Weight on the basis of pixel value needs to be done after the overscan has been subtracted
    200     if (options->doVarianceBuild) {
     255    if (options->doVarianceBuild)
     256    {
    201257        // create the target mask and variance images
    202258        psImage *noiseImage = NULL;
    203         if (options->doNoiseMap) {
     259        if (options->doNoiseMap)
     260        {
    204261            // XXX convert the noiseMap image to a binned image
    205262            pmReadout *noiseMap = NULL;
    206263            noiseMap = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.NOISEMAP");
    207             noiseImage = psImageCopy (NULL, input->image, PS_TYPE_F32);
    208             psImageInit (noiseImage, 0.0);
     264            noiseImage = psImageCopy(NULL, input->image, PS_TYPE_F32);
     265            psImageInit(noiseImage, 0.0);
    209266
    210267            // XXX this works, but is not really quite right: the model shoud include the
    211268            // offset information, we are not really getting exactly the right mapping from the
    212269            // original file.
    213             // CZW 2012-03-21: I do not believe this is true anymore.  In any case, this seems to be what
    214             //                 ppImageReplaceBackground does to do sky subtraction.
     270            // CZW 2012-03-21: I do not believe this is true anymore.  In any case, this seems to be what
     271            //                 ppImageReplaceBackground does to do sky subtraction.
    215272            psImageBinning *binning = psImageBinningAlloc();
    216273            binning->nXruff = noiseMap->image->numCols;
    … …  
    220277            psImageBinningSetScale(binning, PS_IMAGE_BINNING_LEFT);
    221278
    222             psImageUnbin (noiseImage, noiseMap->image, binning);
    223             psFree (binning);
    224         }
    225 
    226         if (!pmReadoutGenerateVariance(input, noiseImage, true)) {
    227             psError(PS_ERR_UNKNOWN, false, "Unable to generate a variance image.");
    228             psFree(detview);
    229             return false;
    230         }
    231         psFree (noiseImage);
    232         // psLogMsg ("ppImage", 6, "generate variance: %f sec\n", psTimerMark ("detrend.readout"));
    233     }
    234 
    235     if (options->doDark && dark) {
    236         if (!pmDarkApply(input, dark, options->darkMask)) {
     279            psImageUnbin(noiseImage, noiseMap->image, binning);
     280            psFree(binning);
     281        }
     282
     283        if (!pmReadoutGenerateVariance(input, noiseImage, true))
     284        {
     285            psError(PS_ERR_UNKNOWN, false, "Unable to generate a variance image.");
     286            psFree(detview);
     287            return false;
     288        }
     289        psFree(noiseImage);
     290        // psLogMsg ("ppImage", 6, "generate variance: %f sec\n", psTimerMark ("detrend.readout"));
     291    }
     292
     293    if (options->doDark && dark)
     294    {
     295        if (!pmDarkApply(input, dark, options->darkMask))
     296        {
    237297            psError(PS_ERR_UNKNOWN, false, "Unable to subtract dark.");
    238298            psFree(detview);
    239299            return false;
    240300        }
    241         // psLogMsg ("ppImage", 6, "apply dark: %f sec\n", psTimerMark ("detrend.readout"));
    242     }
    243 
    244     if (options->doRemnance) {
     301        // psLogMsg ("ppImage", 6, "apply dark: %f sec\n", psTimerMark ("detrend.readout"));
     302    }
     303
     304    if (options->doRemnance)
     305    {
    245306        if (!pmRemnance(input, options->maskValue, options->lowMask,
    246                         options->remnanceSize, options->remnanceThresh)) {
     307                        options->remnanceSize, options->remnanceThresh))
     308        {
    247309            psError(PS_ERR_UNKNOWN, false, "Unable to mask remnance.");
    248310            psFree(detview);
    249311            return false;
    250312        }
    251         // psLogMsg ("ppImage", 6, "mask remnance: %f sec\n", psTimerMark ("detrend.readout"));
     313        // psLogMsg ("ppImage", 6, "mask remnance: %f sec\n", psTimerMark ("detrend.readout"));
    252314    }
    253315
    254316    // Shutter correction
    255     if (options->doShutter) {
     317    if (options->doShutter)
     318    {
    256319        pmReadout *shutter = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.SHUTTER");
    257         if (!pmShutterCorrectionApply(input, shutter, pmConfigMaskGet("FLAT", config))) {
    258             psFree(detview);
    259             return false;
    260         }
    261         // psLogMsg ("ppImage", 6, "shutter correction: %f sec\n", psTimerMark ("detrend.readout"));
     320        if (!pmShutterCorrectionApply(input, shutter, pmConfigMaskGet("FLAT", config)))
     321        {
     322            psFree(detview);
     323            return false;
     324        }
     325        // psLogMsg ("ppImage", 6, "shutter correction: %f sec\n", psTimerMark ("detrend.readout"));
    262326    }
    263327
    264328    // Flat-field correction (no options used?)
    265     if (options->doFlat) {
     329    if (options->doFlat)
     330    {
    266331        pmReadout *flat = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.FLAT");
    267         if (!pmFlatField(input, flat, options->flatMask)) {
    268             psFree(detview);
    269             return false;
    270         }
    271         // psLogMsg ("ppImage", 6, "apply flat: %f sec\n", psTimerMark ("detrend.readout"));
    272     }
    273 
    274 /*     // Pattern noise correction */
    275 /*     if (options->doPattern) { */
    276 /*         if (!pmPatternRow(input, options->patternOrder, options->patternIter, options->patternRej, */
    277 /*                           options->patternThresh, options->patternMean, options->patternStdev, */
    278 /*                           options->maskValue, options->darkMask)) { */
    279 /*             psFree(detview); */
    280 /*             return false; */
    281 /*         } */
    282 /*     } */
     332        if (!pmFlatField(input, flat, options->flatMask))
     333        {
     334            psFree(detview);
     335            return false;
     336        }
     337        // psLogMsg ("ppImage", 6, "apply flat: %f sec\n", psTimerMark ("detrend.readout"));
     338    }
     339
     340    /*     // Pattern noise correction */
     341    /*     if (options->doPattern) { */
     342    /*         if (!pmPatternRow(input, options->patternOrder, options->patternIter, options->patternRej, */
     343    /*                           options->patternThresh, options->patternMean, options->patternStdev, */
     344    /*                           options->maskValue, options->darkMask)) { */
     345    /*             psFree(detview); */
     346    /*             return false; */
     347    /*         } */
     348    /*     } */
    283349
    284350    // Normalization by a single (known) constant
    285     bool mdok;                          // Status of MD lookup
     351    bool mdok; // Status of MD lookup
    286352    float norm = psMetadataLookupF32(&mdok, config->arguments, "NORMALIZATION");
    287     if (mdok && isfinite(norm) && norm != 1.0) {
     353    if (mdok && isfinite(norm) && norm != 1.0)
     354    {
    288355        pmHDU *hdu = pmHDUFromReadout(input); // HDU of interest
    289         psString comment = NULL;        // Comment to add
     356        psString comment = NULL;              // Comment to add
    290357        psStringAppend(&comment, "Normalization: %f", norm);
    291358        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
    … …  
    293360
    294361        psBinaryOp(input->image, input->image, "*", psScalarAlloc(norm, PS_TYPE_F32));
    295         // psLogMsg ("ppImage", 6, "renormalize: %f sec\n", psTimerMark ("detrend.readout"));
    296     }
    297 
    298 # if (1)
     362        // psLogMsg ("ppImage", 6, "renormalize: %f sec\n", psTimerMark ("detrend.readout"));
     363    }
     364
     365#if (1)
    299366    // Normalization by per-class values
    300367    psMetadata *normlist = psMetadataLookupMetadata(&mdok, config->arguments, "NORMALIZATION.TABLE");
    301     if (normlist) {
     368    if (normlist)
     369    {
    302370        pmFPAfile *inputFile = psMetadataLookupPtr(&mdok, config->files, "PPIMAGE.INPUT");
    303371
    304372        // get the menu of class IDs
    305373        psMetadata *menu = psMetadataLookupMetadata(&mdok, inputFile->camera, "CLASSID");
    306         if (!menu) {
     374        if (!menu)
     375        {
    307376            psError(PS_ERR_IO, false, "Unable to find CLASSID metadata in camera configuration");
    308377            psFree(detview);
    … …  
    311380        // get the rule for class_id for the desired class
    312381        const char *rule = psMetadataLookupStr(&mdok, menu, options->normClass);
    313         if (!rule) {
     382        if (!rule)
     383        {
    314384            psError(PS_ERR_IO, false, "Unable to find NORM.CLASS value %s in CLASSID in camera configuration", options->normClass);
    315385            psFree(detview);
    … …  
    318388        // get the class_id from the rule
    319389        char *classID = pmFPAfileNameFromRule(rule, inputFile, view);
    320         if (!classID) {
     390        if (!classID)
     391        {
    321392            psError(PS_ERR_IO, false, "error converting CLASSID rule %s to name\n", rule);
    322393            psFree(detview);
    … …  
    325396
    326397        // get normalization from the class_id
    327         float norm = psMetadataLookupF32 (&mdok, normlist, classID);
    328         if (!mdok) {
     398        float norm = psMetadataLookupF32(&mdok, normlist, classID);
     399        if (!mdok)
     400        {
    329401            psError(PS_ERR_IO, false, "failed to find class ID %s in normalization table\n", classID);
    330402            psFree(detview);
    … …  
    333405
    334406        pmHDU *hdu = pmHDUFromReadout(input); // HDU of interest
    335         psString comment = NULL;        // Comment to add
     407        psString comment = NULL;              // Comment to add
    336408        psStringAppend(&comment, "Normalization: %f", norm);
    337409        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
    … …  
    341413        psBinaryOp(input->image, input->image, "*", psScalarAlloc(norm, PS_TYPE_F32));
    342414
    343         psFree (classID);
    344 
    345         // psLogMsg ("ppImage", 6, "renormalize by class: %f sec\n", psTimerMark ("detrend.readout"));
    346     }
    347 # endif
    348 
    349     if (options->doFringe) {
     415        psFree(classID);
     416
     417        // psLogMsg ("ppImage", 6, "renormalize by class: %f sec\n", psTimerMark ("detrend.readout"));
     418    }
     419#endif
     420
     421    if (options->doFringe)
     422    {
    350423        pmCell *fringe = pmFPAfileThisCell(config->files, detview, "PPIMAGE.FRINGE");
    351         if (!ppImageDetrendFringeMeasure(input, fringe, false, options)) {
    352             psFree(detview);
    353             return false;
    354         }
    355         // psLogMsg ("ppImage", 6, "measure fringe: %f sec\n", psTimerMark ("detrend.readout"));
     424        if (!ppImageDetrendFringeMeasure(input, fringe, false, options))
     425        {
     426            psFree(detview);
     427            return false;
     428        }
     429        // psLogMsg ("ppImage", 6, "measure fringe: %f sec\n", psTimerMark ("detrend.readout"));
    356430    }
    357431
Note: See TracChangeset for help on using the changeset viewer.