IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 14, 2006, 7:13:32 PM (20 years ago)
Author:
Paul Price
Message:

Allowing use of an FPAview to restrict the analysis

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppStats/src/ppStatsLoop.c

    r8338 r8346  
    11#include <stdio.h>
     2#include <assert.h>
    23#include <string.h>
    34#include <pslib.h>
    … …  
    1011static void getMetadata(psMetadata *target, // Target for metadata
    1112                        psMetadata *source, // Source for metadata
    12                         psListIterator *iterator // Iterator for keywords
    13     )
    14 {
    15     psListIteratorSet(iterator, PS_LIST_HEAD);
     13                        psList *list    // List containing keywords
     14    )
     15{
     16    psListIterator *iterator = psListIteratorAlloc(list, PS_LIST_HEAD, false); // Iterator
    1617    psString name;                      // Name from iteration
    1718    while ((name = psListGetAndIncrement(iterator))) {
    … …  
    2122        }
    2223    }
     24    psFree(iterator);
    2325    return;
    2426}
    … …  
    6264
    6365
     66static void cellStats(psMetadata *chipResults, // Metadata holding the chip results
     67                      pmCell *cell,     // Cell for which to get statistics
     68                      psFits *fits,     // FITS file handle
     69                      ppStatsData *data,// The data
     70                      const pmConfig *config // Configuration
     71    )
     72{
     73    assert(chipResults);
     74    assert(cell);
     75    assert(data);
     76    assert(config);
     77
     78    const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell
     79
     80    // Check to see if this is a cell of interest
     81    if (!doThis(data->cells, cellName)) {
     82        return;
     83    }
     84
     85    // Cell-level results
     86    bool mdok;                          // Status of MD lookup
     87    psMetadata *cellResults = psMemIncrRefCounter(psMetadataLookupMD(&mdok, chipResults, cellName));
     88    if (!mdok || !cellResults) {
     89        cellResults = psMetadataAlloc();
     90    }
     91
     92    if (psListLength(data->headers) > 0 && cell->hdu) {
     93        if (fits && !pmCellReadHeader(cell, fits)) {
     94            goto cellDone;
     95        }
     96        pmHDU *hdu = cell->hdu;     // HDU for headers
     97        getMetadata(cellResults, hdu->header, data->headers);
     98    }
     99    if (psListLength(data->concepts) > 0) {
     100        if (fits && !pmCellReadHeader(cell, fits)) {
     101            goto cellDone;
     102        }
     103        pmConceptsReadCell(cell, PM_CONCEPT_SOURCE_ALL, false, config->database);
     104        getMetadata(cellResults, cell->concepts, data->concepts);
     105    }
     106
     107    if (!data->doStats) {
     108        // Nothing further to do --- don't want to waste our time reading the data
     109        if (psListLength(cellResults->list) > 0) {
     110            psMetadataAdd(chipResults, PS_LIST_TAIL, cellName, PS_DATA_METADATA,
     111                          "Results for cell", cellResults);
     112        }
     113        goto cellDone;
     114    }
     115
     116    pmHDU *hdu = pmHDUFromCell(cell); // HDU for cell
     117    if (!hdu || hdu->blankPHU) {
     118        goto cellDone;
     119    }
     120
     121    if (fits && !pmCellRead(cell, fits, config->database)) {
     122        psLogMsg(__func__, PS_LOG_WARN, "Unable to read cell %s\n", cellName);
     123        goto cellDone;
     124    }
     125
     126    psArray *readouts = cell->readouts; // Array of component readouts
     127    if (readouts->n == 0) {
     128        goto cellDone;
     129    }
     130    if (readouts->n > 1) {
     131        psLogMsg(__func__, PS_LOG_WARN, "Multiple readouts (%ld) present in cell %s --- "
     132                 "using only the first.\n", readouts->n, cellName);
     133    }
     134    pmReadout *readout = readouts->data[0]; // The readout of interest
     135    if (!readout->image) {
     136        psLogMsg(__func__, PS_LOG_WARN, "No image associated with readout in cell %s --- "
     137                 "ignored.\n", cellName);
     138        goto cellDone;
     139    }
     140
     141    // Do the statistics
     142    if (data->sample <= 0.0) {
     143        if (!psImageStats(data->stats, readout->image, readout->mask, data->maskVal)) {
     144            psLogMsg(__func__, PS_LOG_WARN, "Unable to perform statistics on cell %s --- "
     145                     "ignored.\n", cellName);
     146            goto cellDone;
     147        }
     148    } else {
     149        // Apply sampling
     150        psImage *image = readout->image; // The image of interest
     151        psImage *mask = readout->mask; // The mask image
     152        int numSamples = data->sample * image->numCols * image->numRows; // Number of samples
     153        int sampleSpace = 1.0 / data->sample; // Space between samples
     154        psVector *sampleValues = psVectorAlloc(numSamples, PS_TYPE_F32); // Vector of samples
     155        sampleValues->n = numSamples;
     156        psVector *sampleMask = NULL;  // Corresponding mask
     157        if (mask) {
     158            sampleMask = psVectorAlloc(numSamples, PS_TYPE_U8);
     159            sampleMask->n = numSamples;
     160        }
     161        for (int i = 0; i < numSamples; i++) {
     162            int j = i * sampleSpace;
     163            int y = j / image->numRows;
     164            int x = j % image->numRows;
     165            sampleValues->data.F32[i] = image->data.F32[y][x];
     166            if (mask) {
     167                sampleMask->data.U8[i] = mask->data.U8[y][x];
     168            }
     169        }
     170        if (!psVectorStats(data->stats, sampleValues, NULL, sampleMask, data->maskVal)) {
     171            psLogMsg(__func__, PS_LOG_WARN, "Unable to perform statistics on cell %s --- "
     172                     "ignored.\n", cellName);
     173            psFree(sampleValues);
     174            psFree(sampleMask);
     175            goto cellDone;
     176        }
     177        psFree(sampleValues);
     178        psFree(sampleMask);
     179    }
     180
     181#define WRITE_STAT(SYMBOL, NAME, SOURCE) \
     182    if (data->stats->options & SYMBOL) { \
     183        psMetadataAddF32(cellResults, PS_LIST_TAIL, NAME, 0, NULL, data->stats->SOURCE); \
     184    }
     185
     186    WRITE_STAT(PS_STAT_SAMPLE_MEAN,     "SAMPLE_MEAN",   sampleMean);
     187    WRITE_STAT(PS_STAT_SAMPLE_MEDIAN,   "SAMPLE_MEDIAN", sampleMedian);
     188    WRITE_STAT(PS_STAT_SAMPLE_STDEV,    "SAMPLE_STDEV",  sampleStdev);
     189    WRITE_STAT(PS_STAT_SAMPLE_QUARTILE, "SAMPLE_LQ",     sampleLQ);
     190    WRITE_STAT(PS_STAT_SAMPLE_QUARTILE, "SAMPLE_UQ",     sampleUQ);
     191    WRITE_STAT(PS_STAT_ROBUST_MEDIAN,   "ROBUST_MEDIAN", robustMedian);
     192    WRITE_STAT(PS_STAT_ROBUST_STDEV,    "ROBUST_STDEV",  robustStdev);
     193    WRITE_STAT(PS_STAT_ROBUST_QUARTILE, "ROBUST_LQ",     robustLQ);
     194    WRITE_STAT(PS_STAT_ROBUST_QUARTILE, "ROBUST_UQ",     robustUQ);
     195    WRITE_STAT(PS_STAT_ROBUST_QUARTILE, "ROBUST_N50",    robustN50);
     196    WRITE_STAT(PS_STAT_FITTED_MEAN,     "FITTED_MEAN",   fittedMean);
     197    WRITE_STAT(PS_STAT_FITTED_STDEV,    "FITTED_STDEV",  fittedStdev);
     198    WRITE_STAT(PS_STAT_CLIPPED_MEAN,    "CLIPPED_MEAN",  clippedMean);
     199    WRITE_STAT(PS_STAT_CLIPPED_STDEV,   "CLIPPED_STDEV", clippedStdev);
     200
     201cellDone:
     202    // Add the cell results to the chip
     203    addToHierarchy(cellResults, chipResults, cellName, "Results for cell");
     204    if (fits) {
     205        pmCellFreeData(cell);
     206    }
     207    return;
     208}
     209
     210static bool chipStats(psMetadata *fpaResults, // Metadata holding the fpa results
     211                      pmChip *chip,     // Chip for which to get statistics
     212                      psFits *fits,     // FITS file handle
     213                      pmFPAview *view,  // View for analysis
     214                      ppStatsData *data,// The data
     215                      const pmConfig *config // Configuration
     216    )
     217{
     218    assert(fpaResults);
     219    assert(chip);
     220    assert(view);
     221    assert(data);
     222    assert(config);
     223
     224    const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip
     225
     226    // Check to see if this is a chip of interest
     227    if (!doThis(data->chips, chipName)) {
     228        return true;
     229    }
     230
     231    psArray *cells = chip->cells;       // Array of cells
     232    if (view->cell >= cells->n) {
     233        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Desired cell view (%d) doesn't match "
     234                "number of cells (%d)\n", view->cell, cells->n);
     235        return false;
     236    }
     237
     238    // Chip-level results
     239    bool mdok;                          // Status of MD lookup
     240    psMetadata *chipResults = psMemIncrRefCounter(psMetadataLookupMD(&mdok, fpaResults, chipName));
     241    if (!mdok || !chipResults) {
     242        chipResults = psMetadataAlloc();
     243    }
     244
     245    if (psListLength(data->headers) > 0 && chip->hdu) {
     246        if (fits && !pmChipReadHeader(chip, fits)) {
     247            goto chipDone;
     248        }
     249        pmHDU *hdu = chip->hdu;     // HDU for headers
     250        getMetadata(chipResults, hdu->header, data->headers);
     251    }
     252    if (psListLength(data->concepts) > 0) {
     253        if (fits && !pmChipReadHeader(chip, fits)) {
     254            goto chipDone;
     255        }
     256        pmConceptsReadChip(chip, PM_CONCEPT_SOURCE_ALL, false, false, config->database);
     257        getMetadata(chipResults, chip->concepts, data->concepts);
     258    }
     259
     260    if (view->cell >= 0) {
     261        pmCell *cell = cells->data[view->cell]; // Cell of interest
     262        cellStats(chipResults, cell, fits, data, config);
     263        goto chipDone;
     264    }
     265
     266    // Iterate over cells
     267    for (int i = 0; i < cells->n; i++) {
     268        pmCell *cell = cells->data[i];  // Cell of interest
     269        cellStats(chipResults, cell, fits, data, config);
     270    }
     271
     272chipDone:
     273    addToHierarchy(chipResults, fpaResults, chipName, "Results for chip");
     274    if (fits) {
     275        pmChipFreeData(chip);
     276    }
     277
     278    return true;
     279}
     280
     281
    64282psMetadata *ppStatsLoop(psMetadata *fpaResults, // Metadata to hold the FPA results
    65283                        ppStatsData *data, // The data
    … …  
    72290    PS_ASSERT_PTR_NON_NULL(fpa, NULL);
    73291
     292    pmFPAview *view = psMemIncrRefCounter(data->view); // View for analysis
     293    if (!view) {
     294        view = pmFPAviewAlloc(0);
     295    }
     296
    74297    if (!fpaResults) {
    75298        fpaResults = psMetadataAlloc();
    76299    }
    77 
    78     bool mdok;                          // Status of MD lookup
    79 
    80     // Iterators for the headers and concepts
    81     psListIterator *headersIter = psListIteratorAlloc(data->headers, PS_LIST_HEAD, false); // Headers
    82     psListIterator *conceptsIter = psListIteratorAlloc(data->concepts, PS_LIST_HEAD, false); // Concepts
    83300
    84301    // Iterate through the FPA
    … …  
    88305        }
    89306        pmHDU *hdu = fpa->hdu;          // HDU for headers
    90         getMetadata(fpaResults, hdu->header, headersIter);
     307        getMetadata(fpaResults, hdu->header, data->headers);
    91308    }
    92309    if (psListLength(data->concepts) > 0) {
    … …  
    95312        }
    96313        pmConceptsReadFPA(fpa, PM_CONCEPT_SOURCE_ALL, false, config->database);
    97         getMetadata(fpaResults, fpa->concepts, conceptsIter);
    98     }
    99     psArray *chips = fpa->chips;        // Array of component chips
    100     for (long i = 0; i < chips->n; i++) {
    101         pmChip *chip = chips->data[i];  // The chip of interest
    102         if (!chip) {
    103             continue;
    104         }
    105         const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip
    106 
    107         // Check to see if this is a chip of interest
    108         if (!doThis(data->chips, chipName)) {
    109             continue;
    110         }
    111 
    112         // Chip-level results
    113         psMetadata *chipResults = psMemIncrRefCounter(psMetadataLookupMD(&mdok, fpaResults, chipName));
    114         if (!mdok || !chipResults) {
    115             chipResults = psMetadataAlloc();
    116         }
    117 
    118         if (psListLength(data->headers) > 0 && chip->hdu) {
    119             if (fits && !pmChipReadHeader(chip, fits)) {
    120                 continue;
    121             }
    122             pmHDU *hdu = chip->hdu;     // HDU for headers
    123             getMetadata(chipResults, hdu->header, headersIter);
    124         }
    125         if (psListLength(data->concepts) > 0) {
    126             if (fits && !pmChipReadHeader(chip, fits)) {
    127                 continue;
    128             }
    129             pmConceptsReadChip(chip, PM_CONCEPT_SOURCE_ALL, false, false, config->database);
    130             getMetadata(chipResults, chip->concepts, conceptsIter);
    131         }
    132 
    133         psArray *cells = chip->cells;   // Array of component cells
    134         for (long j = 0; j < cells->n; j++) {
    135             pmCell *cell = cells->data[j]; // The cell of interest
    136             if (!cell) {
    137                 continue;
    138             }
    139             const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell
    140 
    141             // Check to see if this is a cell of interest
    142             if (!doThis(data->cells, cellName)) {
    143                 continue;
    144             }
    145 
    146             // Cell-level results
    147             psMetadata *cellResults = psMemIncrRefCounter(psMetadataLookupMD(&mdok, chipResults, cellName));
    148             if (!mdok || !cellResults) {
    149                 cellResults = psMetadataAlloc();
    150             }
    151 
    152             if (psListLength(data->headers) > 0 && cell->hdu) {
    153                 if (fits && !pmCellReadHeader(cell, fits)) {
    154                     addToHierarchy(cellResults, chipResults, cellName, "Results for cell");
    155                     continue;
    156                 }
    157                 pmHDU *hdu = cell->hdu;     // HDU for headers
    158                 getMetadata(cellResults, hdu->header, headersIter);
    159             }
    160             if (psListLength(data->concepts) > 0) {
    161                 if (fits && !pmCellReadHeader(cell, fits)) {
    162                     addToHierarchy(cellResults, chipResults, cellName, "Results for cell");
    163                     continue;
    164                 }
    165                 pmConceptsReadCell(cell, PM_CONCEPT_SOURCE_ALL, false, config->database);
    166                 getMetadata(cellResults, cell->concepts, conceptsIter);
    167             }
    168 
    169             if (!data->doStats) {
    170                 // Nothing further to do --- don't want to waste our time reading the data
    171                 if (psListLength(cellResults->list) > 0) {
    172                     psMetadataAdd(chipResults, PS_LIST_TAIL, cellName, PS_DATA_METADATA,
    173                                   "Results for cell", cellResults);
    174                 }
    175                 addToHierarchy(cellResults, chipResults, cellName, "Results for cell");
    176                 continue;
    177             }
    178 
    179             pmHDU *hdu = pmHDUFromCell(cell); // HDU for cell
    180             if (!hdu || hdu->blankPHU) {
    181                 addToHierarchy(cellResults, chipResults, cellName, "Results for cell");
    182                 continue;
    183             }
    184 
    185             if (fits && !pmCellRead(cell, fits, config->database)) {
    186                 psLogMsg(__func__, PS_LOG_WARN, "Unable to read chip %s cell %s\n", chipName, cellName);
    187                 pmCellFreeData(cell);
    188                 addToHierarchy(cellResults, chipResults, cellName, "Results for cell");
    189                 continue;
    190             }
    191 
    192             psArray *readouts = cell->readouts; // Array of component readouts
    193             if (readouts->n == 0) {
    194                 if (fits) {
    195                     pmCellFreeData(cell);
    196                 }
    197                 addToHierarchy(cellResults, chipResults, cellName, "Results for cell");
    198                 continue;
    199             }
    200             if (readouts->n > 1) {
    201                 psLogMsg(__func__, PS_LOG_WARN, "Multiple readouts (%ld) present in chip %s cell %s --- "
    202                          "using only the first.\n", readouts->n, chipName, cellName);
    203             }
    204             pmReadout *readout = readouts->data[0]; // The readout of interest
    205             if (!readout->image) {
    206                 psLogMsg(__func__, PS_LOG_WARN, "No image associated with readout in chip %s cell %s --- "
    207                          "ignored.\n", chipName, cellName);
    208                 if (fits) {
    209                     pmCellFreeData(cell);
    210                 }
    211                 addToHierarchy(cellResults, chipResults, cellName, "Results for cell");
    212                 continue;
    213             }
    214 
    215             // Do the statistics
    216             if (data->sample <= 0.0) {
    217                 if (!psImageStats(data->stats, readout->image, readout->mask, data->maskVal)) {
    218                     psLogMsg(__func__, PS_LOG_WARN, "Unable to perform statistics on chip %s cell %s --- "
    219                              "ignored.\n", chipName, cellName);
    220                     if (fits) {
    221                         pmCellFreeData(cell);
    222                     }
    223                     addToHierarchy(cellResults, chipResults, cellName, "Results for cell");
    224                     continue;
    225                 }
    226             } else {
    227                 // Apply sampling
    228                 psImage *image = readout->image; // The image of interest
    229                 psImage *mask = readout->mask; // The mask image
    230                 int numSamples = data->sample * image->numCols * image->numRows; // Number of samples
    231                 int sampleSpace = 1.0 / data->sample; // Space between samples
    232                 psVector *sampleValues = psVectorAlloc(numSamples, PS_TYPE_F32); // Vector of samples
    233                 sampleValues->n = numSamples;
    234                 psVector *sampleMask = NULL;  // Corresponding mask
    235                 if (mask) {
    236                     sampleMask = psVectorAlloc(numSamples, PS_TYPE_U8);
    237                     sampleMask->n = numSamples;
    238                 }
    239                 for (int i = 0; i < numSamples; i++) {
    240                     int j = i * sampleSpace;
    241                     int y = j / image->numRows;
    242                     int x = j % image->numRows;
    243                     sampleValues->data.F32[i] = image->data.F32[y][x];
    244                     if (mask) {
    245                         sampleMask->data.U8[i] = mask->data.U8[y][x];
    246                     }
    247                 }
    248                 if (!psVectorStats(data->stats, sampleValues, NULL, sampleMask, data->maskVal)) {
    249                     psLogMsg(__func__, PS_LOG_WARN, "Unable to perform statistics on chip %s cell %s --- "
    250                              "ignored.\n", chipName, cellName);
    251                     psFree(sampleValues);
    252                     psFree(sampleMask);
    253                     if (fits) {
    254                         pmCellFreeData(cell);
    255                     }
    256                     addToHierarchy(cellResults, chipResults, cellName, "Results for cell");
    257                     continue;
    258                 }
    259                 psFree(sampleValues);
    260                 psFree(sampleMask);
    261             }
    262 
    263             #define WRITE_STAT(SYMBOL, NAME, SOURCE) \
    264             if (data->stats->options & SYMBOL) { \
    265                 psMetadataAddF32(cellResults, PS_LIST_TAIL, NAME, 0, NULL, data->stats->SOURCE); \
    266             }
    267 
    268             WRITE_STAT(PS_STAT_SAMPLE_MEAN,     "SAMPLE_MEAN",   sampleMean);
    269             WRITE_STAT(PS_STAT_SAMPLE_MEDIAN,   "SAMPLE_MEDIAN", sampleMedian);
    270             WRITE_STAT(PS_STAT_SAMPLE_STDEV,    "SAMPLE_STDEV",  sampleStdev);
    271             WRITE_STAT(PS_STAT_SAMPLE_QUARTILE, "SAMPLE_LQ",     sampleLQ);
    272             WRITE_STAT(PS_STAT_SAMPLE_QUARTILE, "SAMPLE_UQ",     sampleUQ);
    273             WRITE_STAT(PS_STAT_ROBUST_MEDIAN,   "ROBUST_MEDIAN", robustMedian);
    274             WRITE_STAT(PS_STAT_ROBUST_STDEV,    "ROBUST_STDEV",  robustStdev);
    275             WRITE_STAT(PS_STAT_ROBUST_QUARTILE, "ROBUST_LQ",     robustLQ);
    276             WRITE_STAT(PS_STAT_ROBUST_QUARTILE, "ROBUST_UQ",     robustUQ);
    277             WRITE_STAT(PS_STAT_ROBUST_QUARTILE, "ROBUST_N50",    robustN50);
    278             WRITE_STAT(PS_STAT_FITTED_MEAN,     "FITTED_MEAN",   fittedMean);
    279             WRITE_STAT(PS_STAT_FITTED_STDEV,    "FITTED_STDEV",  fittedStdev);
    280             WRITE_STAT(PS_STAT_CLIPPED_MEAN,    "CLIPPED_MEAN",  clippedMean);
    281             WRITE_STAT(PS_STAT_CLIPPED_STDEV,   "CLIPPED_STDEV", clippedStdev);
    282 
    283             // Add the cell results to the chip
    284             addToHierarchy(cellResults, chipResults, cellName, "Results for cell");
    285             if (fits) {
    286                 pmCellFreeData(cell);
    287             }
    288         }
    289 
    290         addToHierarchy(chipResults, fpaResults, chipName, "Results for chip");
    291         if (fits) {
    292             pmChipFreeData(chip);
    293         }
    294 
    295     }
     314        getMetadata(fpaResults, fpa->concepts, data->concepts);
     315    }
     316
     317    psArray *chips = fpa->chips;        // Array of chips
     318    if (view->chip >= 0) {
     319        pmChip *chip = chips->data[view->chip]; // Chip of interest
     320        chipStats(fpaResults, chip, fits, view, data, config);
     321        goto fpaDone;
     322    }
     323
     324    // Iterate over cells
     325    for (int i = 0; i < chips->n; i++) {
     326        pmChip *chip = chips->data[i];  // Chip of interest
     327        chipStats(fpaResults, chip, fits, view, data, config);
     328    }
     329
     330fpaDone:
    296331    if (fits) {
    297332        pmFPAFreeData(fpa);
    298333    }
    299     psFree(headersIter);
    300     psFree(conceptsIter);
     334    psFree(view);
    301335
    302336    return fpaResults;
Note: See TracChangeset for help on using the changeset viewer.