IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 23, 2021, 1:48:02 PM (5 years ago)
Author:
eugene
Message:

adding PATTERN.ROW.AMP I/O functions, clean up pmPatternRow code

Location:
branches/eam_branches/ipp-dev-20210817/psModules/src
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-dev-20210817/psModules/src/camera/pmFPAfile.c

    r41535 r41856  
    551551    if (!strcasecmp(type, "KH.CORRECT"))     {
    552552        return PM_FPA_FILE_KH_CORRECT;
     553    }
     554    if (!strcasecmp(type, "PATTERN.ROW.AMP"))     {
     555        return PM_FPA_FILE_PATTERN_ROW_AMP;
    553556    }
    554557    if (!strcasecmp(type, "SUBKERNEL"))     {
    … …  
    606609      case PM_FPA_FILE_KH_CORRECT:
    607610        return ("KH.CORRECT");
     611      case PM_FPA_FILE_PATTERN_ROW_AMP:
     612        return ("PATTERN.ROW.AMP");
    608613      case PM_FPA_FILE_SUBKERNEL:
    609614        return ("SUBKERNEL");

  • branches/eam_branches/ipp-dev-20210817/psModules/src/camera/pmFPAfile.h

    r36834 r41856  
    5252    PM_FPA_FILE_SRCTEXT,
    5353    PM_FPA_FILE_PATTERN,
     54    PM_FPA_FILE_PATTERN_ROW_AMP,
    5455    PM_FPA_FILE_LINEARITY,
    5556    PM_FPA_FILE_EXPNUM,

  • branches/eam_branches/ipp-dev-20210817/psModules/src/camera/pmFPAfileFitsIO.c

    r36834 r41856  
    152152      case PM_FPA_FILE_ASTROM_REFSTARS:
    153153      case PM_FPA_FILE_KH_CORRECT:
     154      case PM_FPA_FILE_PATTERN_ROW_AMP:
    154155        {
    155156          pmHDU *hdu = pmFPAviewThisHDU(view, fpa);

  • branches/eam_branches/ipp-dev-20210817/psModules/src/camera/pmFPAfileIO.c

    r41843 r41856  
    253253        status = pmPatternRead(view, file, config);
    254254        break;
     255      case PM_FPA_FILE_PATTERN_ROW_AMP:
     256        status = pmPatternRowAmpRead(view, file, config);
     257        break;
    255258      case PM_FPA_FILE_SX:
    256259      case PM_FPA_FILE_RAW:
    … …  
    363366      case PM_FPA_FILE_ASTROM_REFSTARS:
    364367      case PM_FPA_FILE_KH_CORRECT:
     368      case PM_FPA_FILE_PATTERN_ROW_AMP:
    365369      case PM_FPA_FILE_JPEG:
    366370      case PM_FPA_FILE_KAPA:
    … …  
    451455      return true;
    452456    }
     457    if (file->type == PM_FPA_FILE_PATTERN_ROW_AMP) {
     458      psTrace("psModules.camera", 6, "skip write for %s, no write function defined", file->name);
     459      return true;
     460    }
    453461
    454462    // open the file if not yet opened
    … …  
    547555      case PM_FPA_FILE_KH_CORRECT:
    548556        psError(PS_ERR_IO, true, "cannot write type KH.CORRECT (%s)", file->name);
     557        break;
     558
     559      case PM_FPA_FILE_PATTERN_ROW_AMP:
     560        psError(PS_ERR_IO, true, "cannot write type PATTERN.ROW.AMP (%s)", file->name);
    549561        break;
    550562
    … …  
    624636      case PM_FPA_FILE_ASTROM_REFSTARS:
    625637      case PM_FPA_FILE_KH_CORRECT:
     638      case PM_FPA_FILE_PATTERN_ROW_AMP:
    626639      case PM_FPA_FILE_LINEARITY:
    627640      case PM_FPA_FILE_EXPNUM:
    … …  
    702715      case PM_FPA_FILE_ASTROM_REFSTARS:
    703716      case PM_FPA_FILE_KH_CORRECT:
     717      case PM_FPA_FILE_PATTERN_ROW_AMP:
    704718      case PM_FPA_FILE_EXPNUM:
    705719        psTrace ("psModules.camera", 6, "NOT freeing %s (%s) : save for further analysis\n", file->filename, file->name);
    … …  
    865879      case PM_FPA_FILE_ASTROM_REFSTARS:
    866880      case PM_FPA_FILE_KH_CORRECT:
     881      case PM_FPA_FILE_PATTERN_ROW_AMP:
    867882      case PM_FPA_FILE_LINEARITY:
    868883      case PM_FPA_FILE_EXPNUM:
    … …  
    10691084      case PM_FPA_FILE_ASTROM_MODEL:
    10701085      case PM_FPA_FILE_KH_CORRECT:
     1086      case PM_FPA_FILE_PATTERN_ROW_AMP:
    10711087      case PM_FPA_FILE_SX:
    10721088      case PM_FPA_FILE_RAW:

  • branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmDetrendDB.c

    r38038 r41856  
    112112        DETREND_STRING_CASE(AUXMASK);
    113113        DETREND_STRING_CASE(KH_CORRECT);
     114        DETREND_STRING_CASE(PATTERN_ROW_AMP);
    114115    default:
    115116        return NULL;

  • branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmDetrendDB.h

    r36834 r41856  
    4141    PM_DETREND_TYPE_AUXMASK,
    4242    PM_DETREND_TYPE_KH_CORRECT,
     43    PM_DETREND_TYPE_PATTERN_ROW_AMP,
    4344} pmDetrendType;
    4445

  • branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmPattern.c

    r41851 r41856  
    160160
    161161
     162// XXX allow user choice of binned vs unbinned analysis?
    162163bool pmPatternRow(pmReadout *ro, int order, int iter, float rej, float thresh,
    163164                  psStatsOptions clipMean, psStatsOptions clipStdev,
    164165                  psImageMaskType maskVal, psImageMaskType maskBad) {
    165166
    166   bool status = pmPatternRowBinned(ro, order, iter, rej, thresh, clipMean, clipStdev, maskVal, maskBad);
    167   // bool status = pmPatternRowUnbinned(ro, order, iter, rej, thresh, clipMean, clipStdev, maskVal, maskBad);
    168     return status;
     167  bool status = false;
     168  if (true) {
     169    status = pmPatternRowBinned(ro, order, iter, rej, thresh, clipMean, clipStdev, maskVal, maskBad);
     170  } else {
     171    status = pmPatternRowUnbinned(ro, order, iter, rej, thresh, clipMean, clipStdev, maskVal, maskBad);
     172  }
     173  return status;
    169174}
    170175
    171 // USE_BACKGROUND_STDEV: if true, the analysis will use the measured robust stdev to clip the out-of-range pixles
    172 // if false, the stdev will be estimated based on the Poisson statistics of the median (sky level).
     176// USE_BACKGROUND_STDEV: if TRUE, the analysis will use the measured robust stdev to clip the out-of-range pixles
     177// if FALSE, the stdev will be estimated based on the Poisson statistics of the median (sky level).
    173178# define USE_BACKGROUND_STDEV 0
    174179
    … …  
    184189    PS_ASSERT_FLOAT_LARGER_THAN(thresh, 0.0, false);
    185190
     191    bool mdok;                          // Status of MD lookup
     192
     193    pmCell *cell = ro->parent;
     194    pmChip *chip = cell->parent;
     195    const char *chipName = psMetadataLookupStr(&mdok, chip->concepts, "CHIP.NAME"); // Name of chip
     196    const char *cellName = psMetadataLookupStr(&mdok, cell->concepts, "CELL.NAME"); // Name of cell
     197
    186198    psImage *image = ro->image;         // Image to correct
    187199    psImage *mask = ro->mask;           // Mask for image
    … …  
    191203    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
    192204    if (!psImageBackground(stats, NULL, ro->image, ro->mask, maskVal, rng)) {
    193         psWarning("Unable to calculate statistics on readout; masking entire readout.");
     205        psWarning("Unable to calculate statistics on readout; skipping pattern correction for %s, %s.", chipName, cellName);
    194206        psErrorClear();
    195207        psFree(stats);
    196208        psFree(rng);
    197         psImageInit(image, NAN);
    198         if (mask) {
    199             psBinaryOp(mask, mask, "|", psScalarAlloc(maskBad, PS_TYPE_IMAGE_MASK));
    200         }
    201         if (ro->variance) {
    202             psImageInit(image, NAN);
    203         }
    204209        return true;
    205210    }
    … …  
    226231    // so smaller by a factor of 1.4 than what we predict, which is not very large
    227232
    228     # define READNOISE 10
     233    // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file)
     234    float nominalAmplitude = psMetadataLookupF32 (&mdok, cell->analysis, "PTN.ROW.AMP");
     235    if (!mdok) nominalAmplitude = 20; // XXX EAM : somewhat arbitrary number
     236    // If we cannot determine the nominal amplitude, we fall-back on the worst case
     237
     238    // XXX retrieve noise and gain from 'concepts' if possible
     239# define READNOISE 10 /* arbitrary number */
    229240    float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE));
    230     float lower = stats->robustMedian - thresh * sigma; // Lower bound for data
    231     float upper = stats->robustMedian + thresh * sigma; // Upper bound for data
     241    float delta = PS_MIN (thresh * sigma, 2*nominalAmplitude);
     242    float lower = stats->robustMedian - delta; // Lower bound for data
     243    float upper = stats->robustMedian + delta; // Upper bound for data
    232244    float background = stats->robustMedian;
     245
     246    // if the noise from the background is too large
     247    float significance = nominalAmplitude / sigma;
     248   
     249    // XXX EAM : arbitrary number
     250    if (isfinite(nominalAmplitude) && (significance < 1.0/6.0)) {
     251      psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
     252      return true;
     253    }
     254    fprintf (stderr, "correcting pattern row background %s: %f - %f - %f : %f %f %f\n", cellName, lower, background, upper, sigma, nominalAmplitude, significance);
     255    psLogMsg("ppImage", PS_LOG_INFO, "Performing row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
    233256# endif   
    234257   
    235     pmCell *cell = ro->parent;
    236     const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell
    237     fprintf (stderr, "pattern row background %s: %f - %f - %f\n", cellName,lower, background, upper);
    238 
    239 
    240     // XXX add code to skip fit if background stdev is large compared to expected signal
    241     // XXX add a constraint to the fit so the |amplitude| < x (x ~ 20 - 25)
    242 
    243258    psFree(stats);
    244259    psFree(rng);
    … …  
    458473}
    459474
    460 # define NPIX 31   
     475# define NPIX 15
    461476
    462477// bin by NPIX in the x-direction to reduce the number of calculations needed to measure
    … …  
    473488    PS_ASSERT_FLOAT_LARGER_THAN(thresh, 0.0, false);
    474489
     490    bool mdok;                          // Status of MD lookup
     491
     492    pmCell *cell = ro->parent;
     493    pmChip *chip = cell->parent;
     494    const char *chipName = psMetadataLookupStr(&mdok, chip->concepts, "CHIP.NAME"); // Name of chip
     495    const char *cellName = psMetadataLookupStr(&mdok, cell->concepts, "CELL.NAME"); // Name of cell
     496
    475497    psImage *image = ro->image;         // Image to correct
    476498    psImage *mask = ro->mask;           // Mask for image
    477499    int numCols = image->numCols, numRows = image->numRows; // Size of image
    478500
    479     // use the background and typical amplitude value to choose to correct or not
    480     // pmCell *cell = ro->parent;
    481 
    482501    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    483502    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
    484503    if (!psImageBackground(stats, NULL, ro->image, ro->mask, maskVal, rng)) {
    485         psWarning("Unable to calculate statistics on readout; skipping pattern correction.");
     504        psWarning("Unable to calculate statistics on readout; skipping pattern correction for %s, %s.", chipName, cellName);
    486505        psErrorClear();
    487506        psFree(stats);
    … …  
    525544    // so smaller by a factor of 1.4 than what we predict, which is not very large
    526545
     546    // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file)
     547    float nominalAmplitude = psMetadataLookupF32 (&mdok, cell->analysis, "PTN.ROW.AMP");
     548    if (!mdok) nominalAmplitude = 20; // XXX EAM : somewhat arbitrary number
     549    if (!isfinite(nominalAmplitude)) nominalAmplitude = 20; // XXX EAM : somewhat arbitrary number
     550    // If we cannot determine the nominal amplitude, we fall-back on the worst case
     551
    527552    // retrieve noise and gain from 'concepts' if possible
    528553# define READNOISE 10 /* arbitrary number */
    529554    float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE));
    530     float delta = PS_MIN (thresh * sigma, 40); // XXX EAM : arbitrary number
     555    float delta = PS_MIN (thresh * sigma, 5*nominalAmplitude);
    531556    float lower = stats->robustMedian - delta; // Lower bound for data
    532557    float upper = stats->robustMedian + delta; // Upper bound for data
    533558    float background = stats->robustMedian;
    534559
    535     // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file)
    536     float nominalAmplitude = 20; // XXX EAM : arbitrary number
    537 
    538560    // if the noise from the background is too large
    539561    float significance = nominalAmplitude / sigma;
    540     if (significance < 0.5) return true; // XXX EAM : arbitrary number
    541    
    542     pmCell *cell = ro->parent;
    543     const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell
    544     fprintf (stderr, "pattern row background %s: %f - %f - %f\n", cellName, lower, background, upper);
     562   
     563    // XXX EAM : arbitrary number
     564    if (isfinite(nominalAmplitude) && (significance < 1.0/6.0)) {
     565      psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
     566      psFree(stats);
     567      psFree(rng);
     568      return true;
     569    }
     570    psLogMsg("ppImage", PS_LOG_INFO, "Performing row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
    545571# endif   
    546572
    … …  
    643669        }
    644670
    645         // XXX how much time is spent in the fitting? to test: make this always true
     671        // If not enough points are valid, skip
    646672        if (num < validNmin) {
    647             // Not enough points to fit
    648             // EAM 20211011 : we used to mask rows which could not be fit, but this is excessive.  just skip.
    649             // patternMaskRow(ro, y, maskBad);
    650673            // Ignore this row in our subsequent fits, because the fit failed.
    651674            yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF;
    … …  
    656679            psWarning("Unable to fit polynomial to row %d", y);
    657680            psErrorClear();
    658             // EAM 20211011 : we used to mask rows which could not be fit, but this is excessive.  just skip.
    659             // patternMaskRow(ro, y, maskBad);
    660681            // Ignore this row in our subsequent fits, because the fit failed.
    661682            yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF;
    … …  
    672693            psWarning("Unable to evaluate polynomial for row %d", y);
    673694            psErrorClear();
    674             // EAM 20211011 : we used to mask rows which could not be fit, but this is excessive.  just skip.
    675             // patternMaskRow(ro, y, maskBad);
    676695            yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF;
    677696            continue;

  • branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmPatternIO.c

    r41851 r41856  
    88#include "pmFPAview.h"
    99#include "pmFPAfile.h"
     10#include "pmFPAUtils.h"
    1011#include "pmFPAfileFitsIO.h"
    1112#include "pmFPAHeader.h"
    … …  
    455456/**************** PatternRowAmp(litude) I/O *************************/
    456457
    457 # if (0)
    458 /********************* Read Data functions *****************************/
    459 
    460 bool pmPatternRowAmpReadForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
     458bool pmPatternRowAmpRead (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
    461459    {
    462460        // read the full model in one pass: require the level to be FPA
    … …  
    499497    }
    500498
     499    int nGood = 0;
     500    int nTotal = 0;
    501501    while (haveData) {
    502502
    … …  
    511511        // determine the chip:
    512512        char *extname = psMetadataLookupStr (&status, header, "EXTNAME");
    513         psLogMsg ("psModules.detrend", 4, "read %ld rows from Pattern Row Amplitude file, extname %s\n", model->n, extname);
     513        psLogMsg ("psModules.detrend", 8, "read %ld rows from Pattern Row Amplitude file, extname %s\n", model->n, extname);
     514        nTotal += model->n;
    514515
    515516        // I expect to find a name of the form: chipName.ptn (eg, XY01.ptn)
    516517        // where chipName like 'XY01'
    517         psAssert (strlen(extname) == 10, "invalid extension %s", extname);
     518        psAssert (strlen(extname) == 8, "invalid extension %s", extname);
    518519        psAssert (extname[5] == 'p', "invalid extension %s", extname);
    519520        psAssert (extname[6] == 't', "invalid extension %s", extname);
    … …  
    530531        for (int i = 0; i < model->n; i++) {
    531532            psMetadata *row = model->data[i];
    532             char *cellName = psMetadataLookupStr(&status, row, "cellname");
    533             float *amplitude = psMetadataLookupF32(&status, row, "medianValue");
     533
     534            char *cellName = psMetadataLookupStr(&status, row, "CELL_NAME");
     535            if (!cellName) continue;
     536
     537            float amplitude = psMetadataLookupF32(&status, row, "VALUE_MEDIAN");
    534538
    535539            int cellNumber = pmChipFindCell (chip, cellName);
    536             pmCell *cell = chip->data[cellNumber];
     540            pmCell *cell = chip->cells->data[cellNumber];
     541            if (!cell) continue;
    537542
    538543            psMetadataAddF32 (cell->analysis, PS_LIST_TAIL, "PTN.ROW.AMP", PS_META_REPLACE, "", amplitude);
     544            nGood ++;
    539545        }
     546
    540547        psFree (model);
    541548        psFree (header);
    … …  
    544551        haveData = psFitsMoveExtNum (file->fits, 1, true);
    545552    }
    546 
    547     return true;
    548 }
    549 # endif
     553    psLogMsg ("psModules.detrend", 4, "read %d of %d rows from Pattern Row Amplitude file\n", nGood, nTotal);
     554
     555    return true;
     556}

  • branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmPatternIO.h

    r26893 r41856  
    3535    );
    3636
     37
     38bool pmPatternRowAmpRead (const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
     39bool pmPatternRowAmpReadFPA (pmFPAfile *file);
     40bool pmPatternRowAmpReadChips (pmFPAfile *file);
     41
    3742#endif
Note: See TracChangeset for help on using the changeset viewer.