IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 27, 2007, 6:36:29 PM (18 years ago)
Author:
eugene
Message:

splitting up ppMergeMask, adding iteration and grow

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppMerge/src/ppMergeMask.c

    r15913 r15937  
    2222    )
    2323{
    24     psArray *filenames = psMetadataLookupPtr(NULL, config->arguments, "INPUT"); // The input file names
    25     assert(filenames);                  // It should be here --- it's put here in ppMergeConfig
    26 
    27     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
    28     pmFPA *fpaOut = data->out;          // Output FPA
    29 
    30     // Iterate over each file
    31     for (int i = 0; i < filenames->n; i++) {
    32         if (! filenames->data[i] || strlen(filenames->data[i]) == 0) {
    33             continue;
    34         }
    35         psFits *fits = data->files->data[i]; // FITS file handle
    36         if (!fits) {
    37             continue;
    38         }
    39         psTrace("ppMerge", 3, "File %d: %s\n", i, (const char*)filenames->data[i]);
    40 
    41         pmFPA *fpaIn = data->in->data[i]; // Input FPA
    42         pmFPAview *view = pmFPAviewAlloc(0); // View of FPA, for iteration
    43         pmChip *chipIn;                 // Input chip of interest
    44         while ((chipIn = pmFPAviewNextChip(view, fpaIn, 1))) {
    45 
    46             // handle chip vs cell statistics & avoid reading the data twice
    47 
    48             // load the data of all cells
    49             pmCell *cellIn;             // Input cell of interest
    50             while ((cellIn = pmFPAviewNextCell(view, fpaIn, 1))) {
    51                 if (!pmCellRead(cellIn, fits, config->database)) continue;
    52                 if (cellIn->readouts->n == 0) continue;
    53             }
    54 
    55             // calculate the readout statistics either for each readout, or across the entire chip
    56             if (options->statsByChip) {
    57                 ppMergeMaskChipStats (chipIn, options, rng);
    58             } else {
    59                 // calculate the stats for each cell independently
    60                 while ((cellIn = pmFPAviewNextCell(view, fpaIn, 1))) {
    61                     if (cellIn->readouts->n == 0) continue;
    62                     pmReadout *roIn;        // Input readout of interest
    63                     while ((roIn = pmFPAviewNextReadout(view, fpaIn, 1))) {
    64                         if (!roIn->image) continue;
    65                         psTrace("ppMerge", 4, "Measure statistics for chip %d, cell %d, ro %d\n",
    66                                 view->chip, view->cell, view->readout);
    67                         ppMergeMaskReadoutStats (roIn, options, rng);
    68                     }
    69                 }
    70             }
    71 
    72             // apply the measured statistics to determine the outliers to be masked
    73             while ((cellIn = pmFPAviewNextCell(view, fpaIn, 1))) {
    74                 if (cellIn->readouts->n == 0) continue;
    75 
    76                 pmCell *cellOut = pmFPAviewThisCell(view, fpaOut); // Output cell
    77                 // Suspect pixels image
    78                 psImage *suspect = psMetadataLookupPtr(NULL, cellOut->analysis, "MASK.SUSPECT");
    79                 bool first = suspect ? false : true;
    80 
    81                 pmReadout *roIn;        // Input readout of interest
    82                 while ((roIn = pmFPAviewNextReadout(view, fpaIn, 1))) {
    83                     if (!roIn->image) {
    84                         continue;
    85                     }
    86                     psTrace("ppMerge", 4, "Flagging suspect pixels in chip %d, cell %d, ro %d\n",
    87                             view->chip, view->cell, view->readout);
    88                     suspect = pmMaskFlagSuspectPixels(suspect, roIn, options->maskSuspect, options->combine->maskVal);
    89                 }
    90 
    91                 if (first) {
    92                     psMetadataAddImage(cellOut->analysis, PS_LIST_TAIL, "MASK.SUSPECT", 0,
    93                                        "Suspect pixels", suspect);
    94                     psFree(suspect);
    95                 }
    96                 pmCellFreeData(cellIn);
    97             }
    98             pmChipFreeData(chipIn);
    99         }
    100         pmFPAFreeData(fpaIn);
    101         psFree(view);
    102     }
    103     psFree(rng);
    104 
    105     pmFPAview *view = pmFPAviewAlloc(0);// View of FPA, for iteration
    106     if (fpaOut->hdu) {
    107         pmFPAUpdateNames(fpaOut, NULL, NULL);
    108     }
    109     pmFPAWriteMask(fpaOut, data->outFile, config->database, true, false); // Write header only
    110     pmChip *chipOut;                    // Output chip of interest
    111     while ((chipOut = pmFPAviewNextChip(view, fpaOut, 1))) {
    112         if (chipOut->hdu) {
    113             chipOut->data_exists = true;
    114             pmFPAUpdateNames(fpaOut, chipOut, NULL);
    115         }
    116         pmChipWriteMask(chipOut, data->outFile, config->database, true, false); // Write header only
    117         pmCell *cellOut;                   // Output cell of interest
    118         while ((cellOut = pmFPAviewNextCell(view, fpaOut, 1))) {
    119             if (cellOut->hdu) {
    120                 chipOut->data_exists = cellOut->data_exists = true;
    121                 pmFPAUpdateNames(fpaOut, chipOut, cellOut);
    122             }
    123             pmCellWriteMask(cellOut, data->outFile, config->database, true); // Write header only
    124 
    125             psImage *suspect = psMetadataLookupPtr(NULL, cellOut->analysis, "MASK.SUSPECT");
    126             if (! suspect) {
    127                 continue;
    128             }
    129 
    130             pmReadout *roOut = pmReadoutAlloc(cellOut); // Output readout
    131             roOut->mask = pmMaskIdentifyBadPixels(suspect, options->combine->maskVal, filenames->n, options->maskBad, options->maskMode);
    132             roOut->data_exists = cellOut->data_exists = chipOut->data_exists = true;
    133 
    134             // Get list of cells for concepts averaging
    135             {
    136                 psList *inCells = psListAlloc(NULL); // List of cells
    137                 for (int i = 0; i < filenames->n; i++) {
    138                     if (! filenames->data[i] || strlen(filenames->data[i]) == 0) {
    139                         continue;
    140                     }
    141                     pmCell *cellIn = pmFPAviewThisCell(view, data->in->data[i]); // Input cell
    142                     psListAdd(inCells, PS_LIST_TAIL, cellIn);
    143                 }
    144                 if (!pmConceptsAverageCells(cellOut, inCells, NULL, NULL, true)) {
    145                     psError(PS_ERR_UNKNOWN, false, "Unable to average cell concepts.");
    146                     psFree(inCells);
    147                     return false;
    148                 }
    149                 psFree(inCells);
    150             }
    151 
    152             // Statistics on the merged cell
    153             if (data->statsFile) {
    154                 if (!data->stats) {
    155                     data->stats = psMetadataAlloc();
    156                 }
    157 
    158                 // Build a fake image to do statistics
    159                 roOut->image = psImageAlloc(roOut->mask->numCols, roOut->mask->numRows, PS_TYPE_F32);
    160                 psImageInit(roOut->image, 1.0);
    161                 if (!ppStatsFPA(data->stats, data->out, view,
    162                                 options->combine->maskVal | pmConfigMask("BLANK", config),
    163                                 config)) {
    164                     psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to generate stats for image.\n");
    165                     return false;
    166                 }
    167                 psFree(roOut->image);
    168                 roOut->image = NULL;
    169             }
    170 
    171             psFree(roOut);              // Drop reference
    172 
    173             if (cellOut->hdu && !cellOut->hdu->blankPHU) {
    174                 psTrace("ppMerge", 5, "Writing out cell HDU.\n");
    175                 pmCellWriteMask(cellOut, data->outFile, config->database, false);
    176                 pmCellFreeData(cellOut);
    177             }
    178         }
    179 
    180         if (chipOut->hdu && !chipOut->hdu->blankPHU) {
    181             psTrace("ppMerge", 5, "Writing out chip HDU.\n");
    182             pmChipWriteMask(chipOut, data->outFile, config->database, false, false);
    183             pmChipFreeData(chipOut);
    184         }
    185     }
    186 
    187     // Get list of FPAs for concepts averaging
    188     {
    189         psList *inFPAs = psListAlloc(NULL); // List of FPAs
    190         for (int i = 0; i < filenames->n; i++) {
    191             if (! filenames->data[i] || strlen(filenames->data[i]) == 0) {
    192                 continue;
    193             }
    194             pmFPA *fpaIn = data->in->data[i]; // Input FPA
    195             psListAdd(inFPAs, PS_LIST_TAIL, fpaIn);
    196         }
    197 
    198         if (!pmConceptsAverageFPAs(fpaOut, inFPAs)) {
    199             psError(PS_ERR_UNKNOWN, false, "Unable to average FPA concepts.");
    200             psFree(inFPAs);
    201             return false;
    202         }
    203         psFree(inFPAs);
    204     }
    205 
    206     if (fpaOut->hdu && !fpaOut->hdu->blankPHU) {
    207         psTrace("ppMerge", 5, "Writing out FPA HDU.\n");
    208         pmFPAWriteMask(fpaOut, data->outFile, config->database, false, false);
    209     }
    210     pmFPAFreeData(fpaOut);
    211 
    212     psFree(view);
     24    for (int i = 0; i < 2; i++) {
     25        if (!ppMergeMaskSuspect (data, options, config)) {
     26            psError(PS_ERR_UNKNOWN, true, "failed on mask suspect.\n");
     27            return false;
     28        }
     29
     30        if (!ppMergeMaskBad (data, options, config)) {
     31            psError(PS_ERR_UNKNOWN, true, "failed on mask bad.\n");
     32            return false;
     33        }
     34    }
     35
     36    if (!ppMergeMaskAverageConcepts (data, options, config)) {
     37        psError(PS_ERR_UNKNOWN, true, "failed on average concepts.\n");
     38        return false;
     39    }
     40
     41    if (!ppMergeMaskGrow (data, options, config)) {
     42        psError(PS_ERR_UNKNOWN, true, "failed on mask grow.\n");
     43        return false;
     44    }
     45
     46    if (!ppMergeMaskWrite (data, options, config)) {
     47        psError(PS_ERR_UNKNOWN, true, "failed on mask write.\n");
     48        return false;
     49    }
    21350
    21451    return true;
    … …  
    21653
    21754bool ppMergeMaskReadoutStats(const pmReadout *readout,
     55                             const pmReadout *output,
    21856                             ppMergeOptions *options, // Options
    21957                             psRandom *rng)
    … …  
    22967    }
    23068
     69    psImage *mask = NULL;
    23170    psImage *image = readout->image;    // Image of interest
    232     psImage *mask = readout->mask;      // Corresponding mask
     71
     72    if (output) {
     73        mask = output->mask;      // Corresponding mask
     74    }
    23375
    23476    if (rng) {
    … …  
    263105
    264106bool ppMergeMaskChipStats (const pmChip *chip,
     107                           const pmChip *output,
    265108                           ppMergeOptions *options,
    266109                           psRandom *rng)
    … …  
    291134
    292135            psImage *image = readout->image;    // Image of interest
    293             psImage *mask = readout->mask;      // Corresponding mask
     136
     137            pmCell *cellOutput = output->cells->data[nCell];
     138            pmReadout *readoutOutput = NULL;
     139            psImage *mask = NULL;
     140            if (cellOutput->readouts->n > 0) {
     141                readoutOutput = cellOutput->readouts->data[0];
     142                mask = readoutOutput->mask;      // Corresponding mask
     143            }
    294144
    295145            // Size of image
    … …  
    318168    }
    319169
     170    // no valid data, skip the chip
     171    if (!values->n) {
     172        psFree(values);
     173        psFree(stats);
     174        psFree(rng);
     175        return true;
     176    }
     177
    320178    // calculate the statistics
    321179    if (!psVectorStats (stats, values, NULL, NULL, 0)) {
    … …  
    327185    }
    328186    if (!isfinite(stats->robustMedian) || !isfinite(stats->robustUQ) || !isfinite(stats->robustLQ)) {
    329         psError(PS_ERR_UNKNOWN, false, "invalide image statistics (nan).\n");
     187        psError(PS_ERR_UNKNOWN, false, "invalid image statistics (nan).\n");
    330188        psFree(values);
    331189        psFree(stats);
    … …  
    347205        }
    348206    }
     207
     208    psLogMsg ("ppMerge", PS_LOG_INFO, "statistics for chip: %f +/- %f\n", stats->robustMedian, stats->robustStdev);
    349209
    350210    psFree(values);
Note: See TracChangeset for help on using the changeset viewer.