IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 3, 2007, 12:40:38 PM (19 years ago)
Author:
eugene
Message:

splitting out ppStatsLoop.c into ppStatsChip, ppStatsCell, ppStatsReadout; adding complete readout-level iteration

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppStats/src/ppStatsLoop.c

    r13999 r14003  
    11# include "ppStatsInternal.h"
    2 
    3 static void getMetadata(psMetadata *target, // Target for metadata
    4                         psMetadata *source, // Source for metadata
    5                         psList *list    // List containing keywords
    6     )
    7 {
    8     psListIterator *iterator = psListIteratorAlloc(list, PS_LIST_HEAD, false); // Iterator
    9     psString name;                      // Name from iteration
    10     while ((name = psListGetAndIncrement(iterator))) {
    11         psMetadataItem *item = psMetadataLookup(source, name); // Item of interest, or NULL
    12         if (item) {
    13             psMetadataAddItem(target, item, PS_LIST_TAIL, 0);
    14         }
    15     }
    16     psFree(iterator);
    17     return;
    18 }
    19 
    20 static void getAnalysis(psMetadata *target, // Output Target for metadata
    21                         psList *headers,    // List containing desired keywords
    22                         psMetadata *source, // Input Source for metadata
    23                         psList *list        // List containing analysis blocks
    24     )
    25 {
    26     bool status;
    27 
    28     psListIterator *iterator = psListIteratorAlloc(list, PS_LIST_HEAD, false); // Iterator
    29     psString name;                      // Name from iteration
    30     while ((name = psListGetAndIncrement(iterator))) {
    31         psMetadata *folder = psMetadataLookupMetadata(&status, source, name); // Item of interest, or NULL
    32         if (folder) {
    33             getMetadata (target, folder, headers);
    34         }
    35     }
    36     psFree(iterator);
    37     return;
    38 }
    39 
    40 static bool doThis(psList *toDoList,    // List of things to do
    41                    const char *this     // The name of "this"
    42     )
    43 {
    44     if (psListLength(toDoList) == 0) {
    45         // No list --- do everything
    46         return true;
    47     }
    48 
    49     psListIterator *iterator = psListIteratorAlloc(toDoList, PS_LIST_HEAD, false); // Iterator
    50     psString test;                      // Test string, from iteration
    51     while ((test = psListGetAndIncrement(iterator))) {
    52         if (strcmp(this, test) == 0) {
    53             // It's in the list --- do it
    54             psFree(iterator);
    55             return true;
    56         }
    57     }
    58     psFree(iterator);
    59     // Couldn't find it --- don't do it
    60     return false;
    61 }
    62 
    63 static void addToHierarchy(psMetadata *source, // Source to add
    64                            psMetadata *target, // Target to which to add
    65                            const char *name, // Name of source
    66                            const char *comment // Comment for source
    67     )
    68 {
    69     if (psListLength(source->list) > 0 && !psMetadataLookup(target, name)) {
    70         psMetadataAdd(target, PS_LIST_TAIL, name, PS_DATA_METADATA,
    71                       comment, source);
    72     }
    73     return;
    74 }
    75 
    76 
    77 psExit ppStatsCell(psMetadata *chipResults, // Metadata holding the chip results
    78                    pmCell *cell,        // Cell for which to get statistics
    79                    psFits *fits,        // FITS file handle
    80                    ppStatsData *data,   // The data
    81                    const pmConfig *config // Configuration
    82     )
    83 {
    84     assert(chipResults);
    85     assert(cell);
    86     assert(data);
    87     assert(config);
    88 
    89     const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell
    90 
    91     // Check to see if this is a cell of interest
    92     if (!doThis(data->cells, cellName)) {
    93         return PS_EXIT_SUCCESS;
    94     }
    95 
    96     // select the header unit for this cell
    97     pmHDU *hdu = pmHDUFromCell(cell); // HDU for cell
    98     if (!hdu || hdu->blankPHU) {
    99         if (fits) {
    100             // No HDU means there's no data in this cell
    101             return PS_EXIT_SUCCESS;
    102         } else {
    103             psError(PS_ERR_UNKNOWN, false, "Can't find HDU for cell\n");
    104             return PS_EXIT_CONFIG_ERROR;
    105         }
    106     }
    107 
    108     /*** psphot and psastro put their results on the readout->analysis metadata (PSPHOT.HEADER,
    109          PSASTRO.HEADER).  we need to pull quantities of interest from those locations. to do
    110          this, we need to select the appropriate readout.  ***/
    111 
    112     // Select the readout
    113     psArray *readouts = cell->readouts; // Array of component readouts
    114     if (readouts->n == 0) {
    115         psLogMsg(__func__, PS_LOG_WARN, "No readouts present in cell %s --- skipping\n", cellName);
    116         goto cellDone;
    117     }
    118     if (readouts->n > 1) {
    119         psLogMsg(__func__, PS_LOG_WARN, "Multiple readouts (%ld) present in cell %s --- "
    120                  "using only the first.\n", readouts->n, cellName);
    121     }
    122     pmReadout *readout = readouts->data[0]; // The readout of interest
    123 
    124     // Extract Header and Concept values from the Cell and Readout->analysis level
    125     bool mdok;                          // Status of MD lookup
    126     psMetadata *cellResults = psMemIncrRefCounter(psMetadataLookupMetadata(&mdok, chipResults, cellName));
    127     if (!mdok || !cellResults) {
    128         cellResults = psMetadataAlloc();
    129     }
    130 
    131     // Extract Header values
    132     if (psListLength(data->headers) > 0 && cell->hdu) {
    133         if (fits && !pmCellReadHeader(cell, fits)) {
    134             psError (PS_ERR_IO, false, "trouble reading cell header\n");
    135             psFree(cellResults);
    136             return PS_EXIT_DATA_ERROR;
    137         }
    138         pmHDU *hdu = cell->hdu;     // HDU for headers
    139         getMetadata(cellResults, hdu->header, data->headers);
    140 
    141         // search in the readout->analysis metadata blocks listed in data->analysis
    142         if (psListLength(data->analysis) > 0) {
    143             getAnalysis (cellResults, data->headers, readout->analysis, data->analysis);
    144         }
    145     }
    146 
    147     // Extract Concept values
    148     if (psListLength(data->concepts) > 0) {
    149         if (fits && cell->hdu && !pmCellReadHeader(cell, fits)) {
    150             psError (PS_ERR_IO, false, "trouble reading cell header\n");
    151             psFree(cellResults);
    152             return PS_EXIT_DATA_ERROR;
    153         }
    154         pmConceptsReadCell(cell, PM_CONCEPT_SOURCE_ALL, false, config->database);
    155         getMetadata(cellResults, cell->concepts, data->concepts);
    156     }
    157 
    158     // Do we want to measure pixel statistics?
    159     if (!data->doStats && psListLength(data->summary)) {
    160         // Nothing further to do --- don't want to waste our time reading the data
    161         goto cellDone;
    162     }
    163 
    164     // Read the image pixel data
    165     if (fits && !pmCellRead(cell, fits, config->database)) {
    166         psError (PS_ERR_IO, false, "trouble reading cell data\n");
    167         psFree(cellResults);
    168         return PS_EXIT_DATA_ERROR;
    169     }
    170 
    171     if (!readout->image) {
    172         psLogMsg(__func__, PS_LOG_WARN, "No image associated with readout in cell %s --- ignored.\n", cellName);
    173         goto cellDone;
    174     }
    175 
    176     // Measure basic image statistics (means, stdevs, etc).
    177     if (data->sample <= 0.0) {
    178         if (!psImageStats(data->stats, readout->image, readout->mask, data->maskVal)) {
    179             psWarning("Unable to perform statistics on cell %s --- ignored.\n", cellName);
    180             goto statsDone;
    181         }
    182     } else {
    183         // Apply sampling
    184         psImage *image = readout->image; // The image of interest
    185         psImage *mask = readout->mask; // The mask image
    186         int numSamples = data->sample * image->numCols * image->numRows; // Number of samples
    187         int sampleSpace = 1.0 / data->sample; // Space between samples
    188         psVector *sampleValues = psVectorAlloc(numSamples, PS_TYPE_F32); // Vector of samples
    189         psVector *sampleMask = NULL;  // Corresponding mask
    190         if (mask) {
    191             sampleMask = psVectorAlloc(numSamples, PS_TYPE_U8);
    192         }
    193         for (int i = 0; i < numSamples; i++) {
    194             int j = i * sampleSpace;
    195             int y = j / image->numCols;
    196             int x = j % image->numCols;
    197             sampleValues->data.F32[i] = image->data.F32[y][x];
    198             if (mask) {
    199                 sampleMask->data.U8[i] = mask->data.U8[y][x];
    200             }
    201         }
    202         if (!psVectorStats(data->stats, sampleValues, NULL, sampleMask, data->maskVal)) {
    203             psLogMsg(__func__, PS_LOG_WARN, "Unable to perform statistics on cell %s --- "
    204                      "ignored.\n", cellName);
    205             psFree(sampleValues);
    206             psFree(sampleMask);
    207             goto statsDone;
    208         }
    209         psFree(sampleValues);
    210         psFree(sampleMask);
    211     }
    212 
    213 #define WRITE_STAT(SYMBOL, NAME, SOURCE) \
    214     if (data->stats->options & SYMBOL) { \
    215         psMetadataAddF32(cellResults, PS_LIST_TAIL, NAME, 0, NULL, data->stats->SOURCE); \
    216     }
    217 
    218     WRITE_STAT(PS_STAT_SAMPLE_MEAN,     "SAMPLE_MEAN",   sampleMean);
    219     WRITE_STAT(PS_STAT_SAMPLE_MEDIAN,   "SAMPLE_MEDIAN", sampleMedian);
    220     WRITE_STAT(PS_STAT_SAMPLE_STDEV,    "SAMPLE_STDEV",  sampleStdev);
    221     WRITE_STAT(PS_STAT_SAMPLE_QUARTILE, "SAMPLE_LQ",     sampleLQ);
    222     WRITE_STAT(PS_STAT_SAMPLE_QUARTILE, "SAMPLE_UQ",     sampleUQ);
    223     WRITE_STAT(PS_STAT_ROBUST_MEDIAN,   "ROBUST_MEDIAN", robustMedian);
    224     WRITE_STAT(PS_STAT_ROBUST_STDEV,    "ROBUST_STDEV",  robustStdev);
    225     WRITE_STAT(PS_STAT_ROBUST_QUARTILE, "ROBUST_LQ",     robustLQ);
    226     WRITE_STAT(PS_STAT_ROBUST_QUARTILE, "ROBUST_UQ",     robustUQ);
    227     WRITE_STAT(PS_STAT_ROBUST_QUARTILE, "ROBUST_N50",    robustN50);
    228     WRITE_STAT(PS_STAT_FITTED_MEAN,     "FITTED_MEAN",   fittedMean);
    229     WRITE_STAT(PS_STAT_FITTED_STDEV,    "FITTED_STDEV",  fittedStdev);
    230     WRITE_STAT(PS_STAT_CLIPPED_MEAN,    "CLIPPED_MEAN",  clippedMean);
    231     WRITE_STAT(PS_STAT_CLIPPED_STDEV,   "CLIPPED_STDEV", clippedStdev);
    232 
    233     // measure other types of statistics tests
    234 
    235 statsDone:
    236     // count saturated pixels
    237     if (psListLength(data->summary) > 0) {
    238         bool get_nSatPixels = false;
    239         bool get_fSatPixels = false;
    240 
    241         psListIterator *iterator = psListIteratorAlloc(data->summary, PS_LIST_HEAD, false);
    242         psString choice;
    243         while ((choice = psListGetAndIncrement(iterator))) {
    244             if (!strcasecmp (choice, "SAT_PIXEL_NUM"))  get_nSatPixels = true;
    245             if (!strcasecmp (choice, "SAT_PIXEL_FRAC")) get_fSatPixels = true;
    246         }
    247 
    248         if (!get_nSatPixels && !get_fSatPixels) {
    249             goto cellDone;
    250         }
    251 
    252         // Get the "concepts" of interest
    253         float saturation = psMetadataLookupF32(&mdok, cell->concepts, "CELL.SATURATION"); // Saturation level
    254         if (!mdok || isnan(saturation)) {
    255             psLogMsg(__func__, PS_LOG_WARN, "CELL.SATURATION is not set --- unable to measure N_SAT_PIXELS.\n");
    256             if (get_nSatPixels) psMetadataAddS32(cellResults, PS_LIST_TAIL, "SAT_PIXEL_NUM", 0, NULL, 0);
    257             if (get_fSatPixels) psMetadataAddF32(cellResults, PS_LIST_TAIL, "SAT_PIXEL_FRAC", 0, NULL, NAN);
    258             goto cellDone;
    259         }
    260 
    261         int nSatPixels = 0;
    262         for (int j = 0; j < readout->image->numRows; j++) {
    263             for (int i = 0; i < readout->image->numCols; i++) {
    264                 if (readout->image->data.F32[j][i] >= saturation) {
    265                     nSatPixels ++;
    266                 }
    267             }
    268         }
    269         if (get_nSatPixels) psMetadataAddS32(cellResults, PS_LIST_TAIL, "SAT_PIXEL_NUM", 0, NULL, nSatPixels);
    270         if (get_fSatPixels) psMetadataAddF32(cellResults, PS_LIST_TAIL, "SAT_PIXEL_FRAC", 0, NULL, nSatPixels / (double)(readout->image->numRows * readout->image->numCols));
    271     }
    272 
    273 cellDone:
    274     // Add the cell results to the chip
    275     addToHierarchy(cellResults, chipResults, cellName, "Results for cell");
    276     psFree (cellResults);
    277     if (fits) {
    278         pmCellFreeData(cell);
    279     }
    280     return PS_EXIT_SUCCESS;
    281 }
    282 
    283 psExit ppStatsChip(psMetadata *fpaResults, // Metadata holding the fpa results
    284                    pmChip *chip,     // Chip for which to get statistics
    285                    psFits *fits,     // FITS file handle
    286                    pmFPAview *view,  // View for analysis
    287                    ppStatsData *data,// The data
    288                    const pmConfig *config // Configuration
    289     )
    290 {
    291     assert(fpaResults);
    292     assert(chip);
    293     assert(view);
    294     assert(data);
    295     assert(config);
    296 
    297     psExit result = PS_EXIT_SUCCESS;
    298 
    299     const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip
    300 
    301     // Check to see if this is a chip of interest
    302     if (!doThis(data->chips, chipName)) {
    303         return PS_EXIT_SUCCESS;
    304     }
    305 
    306     psArray *cells = chip->cells;       // Array of cells
    307     if (view->cell >= cells->n) {
    308         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Desired cell view (%d) doesn't match "
    309                 "number of cells (%ld)\n", view->cell, cells->n);
    310         return PS_EXIT_CONFIG_ERROR;
    311     }
    312 
    313     // Chip-level results
    314     bool mdok;                          // Status of MD lookup
    315     psMetadata *chipResults = psMemIncrRefCounter(psMetadataLookupMetadata(&mdok, fpaResults, chipName));
    316     if (!mdok || !chipResults) {
    317         chipResults = psMetadataAlloc();
    318     }
    319 
    320     if (psListLength(data->headers) > 0 && chip->hdu) {
    321         if (fits && !pmChipReadHeader(chip, fits)) {
    322             psError (PS_ERR_IO, false, "trouble reading chip header\n");
    323             psFree(chipResults);
    324             return PS_EXIT_DATA_ERROR;
    325         }
    326         pmHDU *hdu = chip->hdu;     // HDU for headers
    327         getMetadata(chipResults, hdu->header, data->headers);
    328     }
    329     if (psListLength(data->concepts) > 0) {
    330         if (fits && chip->hdu && !pmChipReadHeader(chip, fits)) {
    331             psError (PS_ERR_IO, false, "trouble reading chip header\n");
    332             psFree(chipResults);
    333             return PS_EXIT_DATA_ERROR;
    334         }
    335         pmConceptsReadChip(chip, PM_CONCEPT_SOURCE_ALL, false, false, config->database);
    336         getMetadata(chipResults, chip->concepts, data->concepts);
    337     }
    338 
    339     if (view->cell >= 0) {
    340         pmCell *cell = cells->data[view->cell]; // Cell of interest
    341         result = ppStatsCell(chipResults, cell, fits, data, config);
    342         if (result != PS_EXIT_SUCCESS) {
    343             psError(PS_ERR_UNKNOWN, false, "trouble with cell stats for %d\n", view->cell);
    344         }
    345         addToHierarchy(chipResults, fpaResults, chipName, "Results for chip");
    346         psFree (chipResults);
    347         if (fits) {
    348             pmChipFreeData(chip);
    349         }
    350         return result;
    351     }
    352 
    353     // Iterate over cells
    354     for (int i = 0; i < cells->n; i++) {
    355         pmCell *cell = cells->data[i];  // Cell of interest
    356         result = ppStatsCell(chipResults, cell, fits, data, config);
    357         if (result != PS_EXIT_SUCCESS) {
    358             psError(PS_ERR_UNKNOWN, false, "trouble with cell stats for %d\n", i);
    359             if (fits) {
    360                 pmChipFreeData(chip);
    361             }
    362             psFree(chipResults);
    363             return result;
    364         }
    365     }
    366 
    367     addToHierarchy(chipResults, fpaResults, chipName, "Results for chip");
    368     psFree (chipResults);
    369     if (fits) {
    370         pmChipFreeData(chip);
    371     }
    372     return PS_EXIT_SUCCESS;
    373 }
    3742
    3753psMetadata *ppStatsLoop(psExit *result,
    … …  
    40129        }
    40230        pmHDU *hdu = fpa->hdu;          // HDU for headers
    403         getMetadata(newResults, hdu->header, data->headers);
     31        p_ppStatsGetMetadata(newResults, hdu->header, data->headers);
    40432    }
    40533    if (psListLength(data->concepts) > 0) {
    … …  
    41240        }
    41341        pmConceptsReadFPA(fpa, PM_CONCEPT_SOURCE_ALL, false, config->database);
    414         getMetadata(newResults, fpa->concepts, data->concepts);
     42        p_ppStatsGetMetadata(newResults, fpa->concepts, data->concepts);
    41543    }
    41644
     45    // check if we can legitimately iterate to the chip level
    41746    psArray *chips = fpa->chips;        // Array of chips
    418     if (view->chip >= 0) {
    419         pmChip *chip = chips->data[view->chip]; // Chip of interest
    420         *result = ppStatsChip(newResults, chip, fits, view, data, config);
    421         if (*result != PS_EXIT_SUCCESS) {
    422             psError(PS_ERR_UNKNOWN, false, "trouble with stats for cell %d\n", view->cell);
    423             psFree (view);
    424             psFree (newResults);
    425             return NULL;
    426         }
    427         if (fits) {
    428             pmFPAFreeData(fpa);
    429         }
    430         psFree(view);
    431         return newResults;
     47    if (view->chip >= chips->n) {
     48        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Desired chip view (%d) doesn't match "
     49                "number of chips (%ld)\n", view->chip, chips->n);
     50        return NULL;
    43251    }
    43352
    434     // Iterate over chips
     53    // Iterate over chips (if view->chip is set, skip all others)
    43554    for (int i = 0; i < chips->n; i++) {
     55        if ((view->chip >= 0) && (i != view->chip)) continue;
    43656        pmChip *chip = chips->data[i];  // Chip of interest
    43757        *result = ppStatsChip(newResults, chip, fits, view, data, config);
    … …  
    44767        pmFPAFreeData(fpa);
    44868    }
     69
    44970    psFree(view);
    450 
    45171    return newResults;
    45272}
Note: See TracChangeset for help on using the changeset viewer.