IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 28, 2008, 5:08:31 PM (18 years ago)
Author:
Paul Price
Message:

Merging in pap_branch_080320 --- modernised version of ppMerge.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppMerge/src/ppMergeScaleZero.c

    r16842 r17227  
    1010
    1111#include "ppMerge.h"
    12 #include "ppMergeScaleZero.h"
    1312
    1413// Get the scale and zero for each chip of each input
    15 bool ppMergeScaleZero(psImage **scales, // The scales for each integration/cell
    16                       psImage **zeros, // The zeroes for each integration/cell
    17                       psArray **shutters, // The shutter correction data for each cell
    18                       ppMergeData *data,// The data
    19                       const ppMergeOptions *options, // The options
    20                       const pmConfig *config // The configuration
    21     )
     14bool ppMergeScaleZero(pmConfig *config)
    2215{
    23     assert(data);
    24     assert(options);
    2516    assert(config);
    2617
    27     if (!options->scale && !options->zero && !options->shutter) {
    28         return true;                    // We did everything we were asked for
     18    ppMergeType type = psMetadataLookupS32(NULL, config->arguments, "TYPE"); // Type of frame
     19    int numInputs = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
     20    int numCells = psMetadataLookupS32(NULL, config->arguments, "INPUTS.CELLS"); // Number of cells
     21    psStatsOptions meanStat = psMetadataLookupS32(NULL, config->arguments, "MEAN"); // Statistic for mean
     22    psStatsOptions stdevStat = psMetadataLookupS32(NULL, config->arguments, "STDEV"); // Statistic for stdev
     23    psMaskType maskVal = psMetadataLookupU8(NULL, config->arguments, "MASKVAL"); // Value to mask
     24    int shutterSize = psMetadataLookupS32(NULL, config->arguments, "SHUTTER.SIZE"); // Size of shutter region
     25
     26    psVector *gains = NULL;             // Gains for each cell
     27    psArray *shutters = NULL;           // Shutter data for each cell
     28    psStats *stats = NULL;              // Statistics for background
     29    psImage *background = NULL;         // Background measurements per cell per file
     30
     31    switch (type) {
     32      case PPMERGE_TYPE_BIAS:
     33      case PPMERGE_TYPE_DARK:
     34        // Nothing to measure
     35        return true;
     36      case PPMERGE_TYPE_FLAT:
     37      case PPMERGE_TYPE_FRINGE:
     38        gains = psVectorAlloc(numCells, PS_TYPE_F32);
     39        background = psImageAlloc(numCells, numInputs, PS_TYPE_F32);
     40        psImageInit(background, NAN);
     41        stats = psStatsAlloc(meanStat);
     42        break;
     43      case PPMERGE_TYPE_SHUTTER:
     44        shutters = psArrayAlloc(numCells);
     45        break;
     46      case PPMERGE_TYPE_MASK:
     47      default:
     48        break;
    2949    }
    30 
    31     assert(config->camera);             // Need the camera configuration
    32     assert(config->arguments);          // Need the list of files
    33 
    34     psArray *filenames = psMetadataLookupPtr(NULL, config->arguments, "INPUT"); // The input file names
    35     assert(filenames);                  // It should be here --- it's put here in ppMergeConfig
    36 
    37     // Sanity checks
    38     assert(!options->scale || scales);
    39     assert(!scales || !*scales || ((*scales)->type.type == PS_TYPE_F32 &&
    40                                    (*scales)->numCols == data->numCells &&
    41                                    (*scales)->numRows == filenames->n));
    42     assert(!options->zero || zeros);
    43     assert(!zeros || !*zeros || ((*zeros)->type.type == PS_TYPE_F32 &&
    44                                  (*zeros)->numCols == data->numCells &&
    45                                  (*zeros)->numRows == filenames->n));
    46     assert(!options->shutter || shutters);
    47     assert(!shutters || !*shutters || (*shutters)->n == data->numCells);
    48 
    49     // Allocate the outputs
    50     if (options->scale) {
    51         if (*scales) {
    52             psFree(*scales);
    53         }
    54         *scales = psImageAlloc(data->numCells, filenames->n, PS_TYPE_F32);
     50    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
     51    pmFPAview *view = NULL;             // View into FPA
     52
     53    for (int i = 0; i < numInputs; i++) {
     54        pmFPAfileActivate(config->files, false, NULL);
     55        psArray *files = ppMergeFileActivateSingle(config, PPMERGE_FILES_INPUT, true, i); // Activated files
     56        pmFPAfile *input = files->data[0]; // Representative file; should be the image (not mask or weight)
     57        pmFPA *fpa = input->fpa;        // FPA of interest
     58        view = pmFPAviewAlloc(0);       // View to component of interest
     59        int cellNum = 0;                // Index for cell
     60        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
     61            goto ERROR;
     62        }
     63        pmChip *chip;                   // Chip of interest
     64        while ((chip = pmFPAviewNextChip(view, fpa, 1))) {
     65            if (!chip->file_exists) {
     66                continue;
     67            }
     68            if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
     69                goto ERROR;
     70            }
     71
     72            pmCell *cell;               // Cell of interest
     73            while ((cell = pmFPAviewNextCell(view, fpa, 1))) {
     74                if (!cell->file_exists) {
     75                    continue;
     76                }
     77                if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
     78                    goto ERROR;
     79                }
     80
     81                if (!cell->data_exists) {
     82                    continue;
     83                }
     84
     85                if (cell->readouts->n > 1) {
     86                    psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     87                            "File %d chip %d cell %d contains more than one readout (%ld)",
     88                            i, view->chip, view->cell, cell->readouts->n);
     89                    goto ERROR;
     90                }
     91                pmReadout *readout = cell->readouts->data[0]; // Readout of interest
     92
     93                switch (type) {
     94                  case PPMERGE_TYPE_FLAT:
     95                  case PPMERGE_TYPE_FRINGE: {
     96                      // Extract the gain
     97                      float gain = psMetadataLookupF32(NULL, cell->concepts, "CELL.GAIN"); // Cell gain
     98                      if (!isfinite(gain)) {
     99                          psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     100                                  "CELL.GAIN for file %d chip %d cell %d is not set.",
     101                                  i, view->chip, view->cell);
     102                          goto ERROR;
     103                      }
     104                      gains->data.F32[cellNum] = gain;
     105
     106                      // Measure the background
     107                      if (!psImageBackground(stats, NULL, readout->image, readout->mask, maskVal, rng)) {
     108                          psError(PS_ERR_UNKNOWN, false,
     109                                  "Unable to get statistics for file %d chip %d cell %d",
     110                                  i, view->chip, view->cell);
     111                          goto ERROR;
     112                      }
     113                      background->data.F32[i][cellNum] = psStatsGetValue(stats, meanStat);
     114                      break;
     115                  }
     116                  case PPMERGE_TYPE_SHUTTER: {
     117                      pmShutterCorrectionData *shutter = shutters->data[cellNum]; // Shutter correction data
     118                      if (!shutter) {
     119                          shutter = pmShutterCorrectionDataAlloc(readout->image->numCols,
     120                                                                 readout->image->numRows,
     121                                                                 shutterSize);
     122                          shutters->data[cellNum] = shutter;
     123                      }
     124                      if (!pmShutterCorrectionAddReadout(shutter, readout, meanStat, stdevStat,
     125                                                         maskVal, rng)) {
     126                          psError(PS_ERR_UNKNOWN, false,
     127                                  "Can't add file %d chip %d cell %d to shutter correction.",
     128                                  i, view->chip, view->cell);
     129                          goto ERROR;
     130                      }
     131                      break;
     132                  }
     133                  case PPMERGE_TYPE_MASK:
     134                  default:
     135                    psAbort("Should never get here.");
     136                }
     137
     138                cellNum++;
     139
     140                if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
     141                    goto ERROR;
     142                }
     143            }
     144
     145            if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
     146                goto ERROR;
     147            }
     148        }
     149
     150        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
     151            goto ERROR;
     152        }
     153
     154        psFree(view);
     155
     156#if 0
     157        // Reset files for reading again
     158        for (int i = 0; i < files->n; i++) {
     159            pmFPAfile *file = files->data[i]; // File of interest
     160        }
     161#endif
     162        psFree(files);
    55163    }
    56     if (options->zero) {
    57         if (*zeros) {
    58             psFree(*zeros);
    59         }
    60         *zeros = psImageAlloc(data->numCells, filenames->n, PS_TYPE_F32);
     164
     165    psFree(rng); rng = NULL;
     166    psFree(stats); stats = NULL;
     167
     168    // Store results
     169    switch (type) {
     170      case PPMERGE_TYPE_FRINGE:
     171        psMetadataAddImage(config->arguments, PS_LIST_TAIL, "ZEROS", 0,
     172                           "Zero to subtract from each input cell", background);
     173        // Flow through
     174      case PPMERGE_TYPE_FLAT: {
     175          // Need to normalize over the focal plane
     176          if (psTraceGetLevel("ppMerge") > 9) {
     177              for (int i = 0; i < gains->n; i++) {
     178                  psTrace("ppMerge", 10, "Gain for cell %d is %f\n", i, gains->data.F32[i]);
     179              }
     180          }
     181          psVector *fluxes = NULL;        // Solution to fluxes
     182          if (!pmFlatNormalize(&fluxes, &gains, background)) {
     183              psError(PS_ERR_UNKNOWN, false, "Normalisation failed to converge --- continuing anyway.");
     184              psFree(fluxes);
     185              goto ERROR;
     186          }
     187
     188          psMetadataAddVector(config->arguments, PS_LIST_TAIL, "SCALES", 0,
     189                              "Scale to divide into each input file", fluxes);
     190          psFree(fluxes);               // Drop reference
     191          break;
     192      }
     193      case PPMERGE_TYPE_SHUTTER:
     194        psMetadataAddArray(config->arguments, PS_LIST_TAIL, "SHUTTER", 0,
     195                           "Shutter data", shutters);
     196        break;
     197      case PPMERGE_TYPE_MASK:
     198      default:
     199        psAbort("Should never get here.");
    61200    }
    62     psRandom *rng = NULL;               // Random number generator
    63     if (options->shutter) {
    64         if (*shutters) {
    65             psFree(*shutters);
    66         }
    67         *shutters = psArrayAlloc(data->numCells);
    68         rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
    69     }
    70 
    71     bool fromConcepts = false;          // Do we get the scale and zero points from the concepts?
    72     bool first = true;                  // Are we on the first cell (that sets the standard for the rest)?
    73     bool done = false;                  // Are we done going through the list?
    74     bool mdok = true;                   // Status of MD lookup
    75     for (long i = 0; i < data->in->n && !done; i++) {
    76         pmFPA *fpa = data->in->data[i]; // The FPA
    77         if (!fpa) {
    78             continue;
    79         }
    80         long cellNum = -1;              // Number of the cell
    81         psArray *chips = fpa->chips;    // The array of chips
    82         for (long j = 0; j < chips->n && !done; j++) {
    83             pmChip *chip = chips->data[j]; // The chip
    84             if (!chip) {
    85                 continue;
    86             }
    87             psArray *cells = chip->cells; // The array of cells
    88             for (long k = 0; k < cells->n && !done; k++) {
    89                 pmCell *cell = cells->data[k]; // The cell
    90                 if (!cell) {
    91                     continue;
    92                 }
    93                 cellNum++;
    94 
    95                 if (options->scale) {
    96                     float scale = psMetadataLookupF32(&mdok, cell->concepts, "PPMERGE.SCALE"); // The scale
    97                     if (mdok && !isnan(scale)) {
    98                         if (!first && !fromConcepts) {
    99                             psWarning("PPMERGE.SCALE and PPMERGE.ZERO have been set "
    100                                      "for some, but not all cells --- we will re-measure it for all cells.");
    101                             done = true;
    102                             continue;
    103                         }
    104                         fromConcepts = true;
    105                         (*scales)->data.F32[i][cellNum] = scale;
    106                         psTrace("ppMerge", 9, "Scale for input %ld, chip %ld, cell %ld: %f\n",
    107                                 i, j, k, scale);
    108                     } else if (!first && fromConcepts) {
    109                         psWarning("PPMERGE.SCALE and PPMERGE.ZERO have been set "
    110                                  "for some, but not all cells --- we will re-measure it for all cells.");
    111                         fromConcepts = false;
    112                         done = true;
    113                         continue;
    114                     }
    115                 }
    116 
    117                 if (options->zero) {
    118                     float zero = psMetadataLookupF32(&mdok, cell->concepts, "PPMERGE.ZERO"); // The zero
    119                     if (mdok && !isnan(zero)) {
    120                         if (!first && !fromConcepts) {
    121                             psWarning("PPMERGE.SCALE and PPMERGE.ZERO have been set "
    122                                      "for some, but not all cells --- we will re-measure it for all cells.");
    123                             done = true;
    124                             continue;
    125                         }
    126                         fromConcepts = true;
    127                         (*zeros)->data.F32[i][cellNum] = zero;
    128                         psTrace("ppMerge", 9, "Zero for input %ld, chip %ld, cell %ld: %f\n", i, j, k, zero);
    129                     } else if (!first && fromConcepts) {
    130                         psWarning("PPMERGE.SCALE and PPMERGE.ZERO have been set "
    131                                  "for some, but not all cells --- we will re-measure it for all cells.");
    132                         fromConcepts = false;
    133                         done = true;
    134                         continue;
    135                     }
    136                 }
    137 
    138                 first = false;
    139             }
    140         }
    141     }
    142 
    143     if (fromConcepts) {
    144         // We've already done everything we need to
    145         psFree(rng);
    146         return true;
    147     }
    148 
    149     psImage *background = psImageAlloc(data->numCells, filenames->n, PS_TYPE_F32); // Background measurements
    150     psImageInit(background, NAN);
    151     psStats *bgStats = psStatsAlloc(options->mean); // Statistic to measure the background
    152     psVector *gains = NULL;             // The gains for each cell
    153     if (options->scale) {
    154         gains = psVectorAlloc(data->numCells, PS_TYPE_F32);
    155     }
    156 
    157     pmFPAview *view = pmFPAviewAlloc(0);
    158 
    159     bool status = true;                 // Status of getting the scale and zero --- did everything go right?
    160     for (int i = 0; i < filenames->n; i++) {
    161         psString name = filenames->data[i]; // The name of the file
    162         if (!name || strlen(name) == 0) {
    163             continue;
    164         }
    165         psTrace("ppMerge", 9, "Opening %s to get background...\n", name);
    166         psFits *inFile = data->files->data[i]; // The FITS file to read
    167         pmFPA *fpa = data->in->data[i]; // The FPA for this input
    168         int cellNum = -1;               // Number of the cell
    169         psArray *chips = fpa->chips;    // Array of chips
    170         for (int j = 0; j < chips->n; j++) {
    171             pmChip *chip = chips->data[j]; // The chip of interest
    172             if (!chip) {
    173                 continue;
    174             }
    175             psArray *cells = chip->cells; // Array of cells
    176             for (int k = 0; k < cells->n; k++) {
    177                 pmCell *cell = cells->data[k]; // The cell of interest
    178                 if (!cell) {
    179                     continue;
    180                 }
    181                 cellNum++;
    182 
    183                 if (!pmCellReadHeader(cell, inFile)) {
    184                     continue;
    185                 }
    186 
    187                 // Scaling by the background
    188                 if (options->scale || options->zero) {
    189                     if (!pmCellRead(cell, inFile, config->database)) {
    190                         // Nothing here
    191                         pmCellFreeData(cell);
    192                         continue;
    193                     }
    194 
    195                     if (cell->readouts->n > 1) {
    196                         psWarning("File %s chip %d cell %d contains more than one "
    197                                  "readout --- ignoring all but the first.\n", name, j, k);
    198                         status = false;
    199                     }
    200 
    201                     pmReadout *readout = cell->readouts->data[0]; // The readout of interest
    202                     psImage *image = readout->image; // The pixels of interest
    203                     if (!image) {
    204                         pmCellFreeData(cell);
    205                         continue;
    206                     }
    207 
    208                     // Get the gain
    209                     if (options->scale) {
    210                         bool mdok = true;   // Status of MD lookup
    211                         gains->data.F32[cellNum] = psMetadataLookupF32(&mdok, cell->concepts, "CELL.GAIN");
    212                         if (!mdok || isnan(gains->data.F32[cellNum])) {
    213                             psWarning("CELL.GAIN for file %s chip %d cell %d is not "
    214                                      "set.\n", name, j, k);
    215                             gains->data.F32[cellNum] = NAN;
    216                             status = false;
    217                         }
    218                     }
    219 
    220                     // Get the background
    221                     int sampleSize = (image->numCols * image->numRows) / options->sample; // Size of sample
    222                     psVector *sample = psVectorAlloc(sampleSize, PS_TYPE_F32); // Sample of the image
    223                     psVector *sampleMask = NULL; // Mask for sample
    224                     if (readout->mask) {
    225                         sampleMask = psVectorAlloc(sampleSize, PS_TYPE_U8);
    226                     }
    227                     psImage *mask = readout->mask; // The mask image
    228                     for (long i = 0; i < sampleSize; i++) {
    229                         int j = i * options->sample; // Index into image
    230                         int x = j % image->numCols; // x index
    231                         int y = j / image->numCols; // y index
    232                         sample->data.F32[i] = image->data.F32[y][x];
    233                         if (readout->mask) {
    234                             sampleMask->data.U8[i] = mask->data.U8[y][x];
    235                         }
    236                     }
    237                     status = psVectorStats(bgStats, sample, sampleMask, NULL, options->combine->maskVal);
    238                     if (!status) {
    239                       psTrace("ppMerge", 3, "failed to get stats for for %s, cell %d is %f\n", name, cellNum,
    240                               background->data.F32[i][cellNum]);
    241                       psErrorClear();
    242                     }
    243                     psFree(sample);
    244                     psFree(sampleMask);
    245                     background->data.F32[i][cellNum] = psStatsGetValue(bgStats, options->mean);
    246                     psTrace("ppMerge", 3, "Background for %s, cell %d is %f\n", name, cellNum,
    247                             background->data.F32[i][cellNum]);
    248                 }
    249 
    250                 // Shutter correction
    251                 if (options->shutter) {
    252                     if (!pmCellRead(cell, inFile, config->database)) {
    253                         // Nothing here
    254                         pmCellFreeData(cell);
    255                         continue;
    256                     }
    257 
    258                     if (cell->readouts->n > 1) {
    259                         psWarning("File %s chip %d cell %d contains more than one "
    260                                   "readout --- ignoring all but the first.\n", name, j, k);
    261                         status = false;
    262                     }
    263                     pmReadout *readout = cell->readouts->data[0]; // The readout of interest
    264 
    265                     pmShutterCorrectionData *shutter = (*shutters)->data[cellNum]; // Shutter correction data
    266                     if (!shutter) {
    267                         shutter = pmShutterCorrectionDataAlloc(readout->image->numCols,
    268                                                                readout->image->numRows,
    269                                                                options->shutterSize);
    270                         (*shutters)->data[cellNum] = shutter;
    271                     }
    272                     if (!pmShutterCorrectionAddReadout(shutter, readout, options->mean, options->stdev,
    273                                                        options->combine->maskVal, rng)) {
    274                         psWarning("Can't add file %s chip %d cell %d to shutter correction --- ignored.",
    275                                   name, j, k);
    276                         status = false;
    277                     }
    278                 }
    279 
    280 
    281                 pmCellFreeData(cell);
    282             }
    283             pmChipFreeData(chip);
    284         }
    285         pmFPAFreeData(fpa);
    286     }
     201
     202    psFree(background);
     203    psFree(shutters);
     204    return true;
     205
     206ERROR:
     207    // Common path for errors
     208    psFree(gains);
     209    psFree(background);
     210    psFree(shutters);
    287211    psFree(rng);
    288     psFree(bgStats);
     212    psFree(stats);
    289213    psFree(view);
    290 
    291     if (options->scale) {
    292         // Need to normalize over the focal plane
    293         if (psTraceGetLevel("ppMerge") > 9) {
    294             for (int i = 0; i < gains->n; i++) {
    295                 psTrace("ppMerge", 10, "Gain for cell %d is %f\n", i, gains->data.F32[i]);
    296             }
    297         }
    298         psVector *fluxes = NULL;        // Solution to fluxes
    299         if (!pmFlatNormalize(&fluxes, &gains, background)) {
    300             psWarning("Normalisation failed to converge --- continuing anyway.\n");
    301             status = false;
    302         }
    303         psFree(gains);
    304 
    305         psImage *scalesDeref = *scales; // Dereference the pointer
    306 
    307         for (int i = 0; i < scalesDeref->numRows; i++) {
    308             psF32 bg = fluxes->data.F32[i];
    309             for (int j = 0; j < scalesDeref->numCols; j++) {
    310                 scalesDeref->data.F32[i][j] = bg;
    311             }
    312         }
    313     }
    314 
    315     if (options->zero) {
    316         if (!*zeros) {
    317             // This is much faster than copying!
    318             *zeros = psMemIncrRefCounter(background);
    319         } else {
    320             *zeros = psImageCopy(*zeros, background, PS_TYPE_F32);
    321         }
    322     }
    323 
    324     // Diagnostic stuff
    325     if (scales && *scales && psTraceGetLevel("ppMerge") > 9) {
    326         psImage *scalesDeref = *scales; // Dereference the pointer
    327         for (int i = 0; i < scalesDeref->numRows; i++) {
    328             for (int j = 0; j < scalesDeref->numCols; j++) {
    329                 psTrace("ppMerge", 9, "Scale for exposure %d, cell %d is: %f\n", i, j,
    330                         scalesDeref->data.F32[i][j]);
    331             }
    332         }
    333     }
    334     if (zeros && *zeros && psTraceGetLevel("ppMerge") > 9) {
    335         psImage *zerosDeref = *zeros; // Dereference the pointer
    336         for (int i = 0; i < zerosDeref->numRows; i++) {
    337             for (int j = 0; j < zerosDeref->numCols; j++) {
    338                 psTrace("ppMerge", 9, "Zero for exposure %d, cell %d is: %f\n", i, j,
    339                         zerosDeref->data.F32[i][j]);
    340             }
    341         }
    342     }
    343 
    344     psFree(background);
    345 
    346     return status;
     214    return false;
    347215}
    348216
Note: See TracChangeset for help on using the changeset viewer.