Index: /trunk/ppStats/src/ppStats.c =================================================================== --- /trunk/ppStats/src/ppStats.c (revision 8345) +++ /trunk/ppStats/src/ppStats.c (revision 8346) @@ -7,4 +7,5 @@ psMetadata *ppStats(psMetadata *out, // Output metadata pmFPA *fpa, // FPA for which to get statistics + pmFPAview *view, // View for analysis pmConfig *config // Configuration ) @@ -18,4 +19,9 @@ data->fpa = psMemIncrRefCounter(fpa); + if (data->view) { + psFree(data->view); + } + data->view = psMemIncrRefCounter(view); + // Go through the FPA and do the hard work out = ppStatsLoop(out, data, config); Index: /trunk/ppStats/src/ppStats.h =================================================================== --- /trunk/ppStats/src/ppStats.h (revision 8345) +++ /trunk/ppStats/src/ppStats.h (revision 8346) @@ -10,3 +10,10 @@ #include "ppStatsLoop.h" +// Perform the ppStats steps +psMetadata *ppStats(psMetadata *out, // Output metadata + pmFPA *fpa, // FPA for which to get statistics + pmFPAview *view, // View for analysis + pmConfig *config // Configuration + ); + #endif Index: /trunk/ppStats/src/ppStatsData.c =================================================================== --- /trunk/ppStats/src/ppStatsData.c (revision 8345) +++ /trunk/ppStats/src/ppStatsData.c (revision 8346) @@ -7,10 +7,10 @@ ) { - // inName and region are not on the psLib memory system (they are from argv). - psFree(data->fpa); if (data->fits) { psFitsClose(data->fits); data->fits = NULL; } + psFree(data->fpa); + psFree(data->view); psFree(data->headers); psFree(data->concepts); @@ -29,6 +29,7 @@ psMemSetDeallocator(data, (psFreeFunc)statsDataFree); + data->fits = NULL; data->fpa = NULL; - data->fits = NULL; + data->view = NULL; data->headers = psListAlloc(NULL); Index: /trunk/ppStats/src/ppStatsData.h =================================================================== --- /trunk/ppStats/src/ppStatsData.h (revision 8345) +++ /trunk/ppStats/src/ppStatsData.h (revision 8346) @@ -9,4 +9,5 @@ psFits *fits; // Input file handle pmFPA *fpa; // FPA to analyse + pmFPAview *view; // View to analyse // Stuff to output psStats *stats; // Statistics to calculate Index: /trunk/ppStats/src/ppStatsLoop.c =================================================================== --- /trunk/ppStats/src/ppStatsLoop.c (revision 8345) +++ /trunk/ppStats/src/ppStatsLoop.c (revision 8346) @@ -1,3 +1,4 @@ #include +#include #include #include @@ -10,8 +11,8 @@ static void getMetadata(psMetadata *target, // Target for metadata psMetadata *source, // Source for metadata - psListIterator *iterator // Iterator for keywords - ) -{ - psListIteratorSet(iterator, PS_LIST_HEAD); + psList *list // List containing keywords + ) +{ + psListIterator *iterator = psListIteratorAlloc(list, PS_LIST_HEAD, false); // Iterator psString name; // Name from iteration while ((name = psListGetAndIncrement(iterator))) { @@ -21,4 +22,5 @@ } } + psFree(iterator); return; } @@ -62,4 +64,220 @@ +static void cellStats(psMetadata *chipResults, // Metadata holding the chip results + pmCell *cell, // Cell for which to get statistics + psFits *fits, // FITS file handle + ppStatsData *data,// The data + const pmConfig *config // Configuration + ) +{ + assert(chipResults); + assert(cell); + assert(data); + assert(config); + + const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell + + // Check to see if this is a cell of interest + if (!doThis(data->cells, cellName)) { + return; + } + + // Cell-level results + bool mdok; // Status of MD lookup + psMetadata *cellResults = psMemIncrRefCounter(psMetadataLookupMD(&mdok, chipResults, cellName)); + if (!mdok || !cellResults) { + cellResults = psMetadataAlloc(); + } + + if (psListLength(data->headers) > 0 && cell->hdu) { + if (fits && !pmCellReadHeader(cell, fits)) { + goto cellDone; + } + pmHDU *hdu = cell->hdu; // HDU for headers + getMetadata(cellResults, hdu->header, data->headers); + } + if (psListLength(data->concepts) > 0) { + if (fits && !pmCellReadHeader(cell, fits)) { + goto cellDone; + } + pmConceptsReadCell(cell, PM_CONCEPT_SOURCE_ALL, false, config->database); + getMetadata(cellResults, cell->concepts, data->concepts); + } + + if (!data->doStats) { + // Nothing further to do --- don't want to waste our time reading the data + if (psListLength(cellResults->list) > 0) { + psMetadataAdd(chipResults, PS_LIST_TAIL, cellName, PS_DATA_METADATA, + "Results for cell", cellResults); + } + goto cellDone; + } + + pmHDU *hdu = pmHDUFromCell(cell); // HDU for cell + if (!hdu || hdu->blankPHU) { + goto cellDone; + } + + if (fits && !pmCellRead(cell, fits, config->database)) { + psLogMsg(__func__, PS_LOG_WARN, "Unable to read cell %s\n", cellName); + goto cellDone; + } + + psArray *readouts = cell->readouts; // Array of component readouts + if (readouts->n == 0) { + goto cellDone; + } + if (readouts->n > 1) { + psLogMsg(__func__, PS_LOG_WARN, "Multiple readouts (%ld) present in cell %s --- " + "using only the first.\n", readouts->n, cellName); + } + pmReadout *readout = readouts->data[0]; // The readout of interest + if (!readout->image) { + psLogMsg(__func__, PS_LOG_WARN, "No image associated with readout in cell %s --- " + "ignored.\n", cellName); + goto cellDone; + } + + // Do the statistics + if (data->sample <= 0.0) { + if (!psImageStats(data->stats, readout->image, readout->mask, data->maskVal)) { + psLogMsg(__func__, PS_LOG_WARN, "Unable to perform statistics on cell %s --- " + "ignored.\n", cellName); + goto cellDone; + } + } else { + // Apply sampling + psImage *image = readout->image; // The image of interest + psImage *mask = readout->mask; // The mask image + int numSamples = data->sample * image->numCols * image->numRows; // Number of samples + int sampleSpace = 1.0 / data->sample; // Space between samples + psVector *sampleValues = psVectorAlloc(numSamples, PS_TYPE_F32); // Vector of samples + sampleValues->n = numSamples; + psVector *sampleMask = NULL; // Corresponding mask + if (mask) { + sampleMask = psVectorAlloc(numSamples, PS_TYPE_U8); + sampleMask->n = numSamples; + } + for (int i = 0; i < numSamples; i++) { + int j = i * sampleSpace; + int y = j / image->numRows; + int x = j % image->numRows; + sampleValues->data.F32[i] = image->data.F32[y][x]; + if (mask) { + sampleMask->data.U8[i] = mask->data.U8[y][x]; + } + } + if (!psVectorStats(data->stats, sampleValues, NULL, sampleMask, data->maskVal)) { + psLogMsg(__func__, PS_LOG_WARN, "Unable to perform statistics on cell %s --- " + "ignored.\n", cellName); + psFree(sampleValues); + psFree(sampleMask); + goto cellDone; + } + psFree(sampleValues); + psFree(sampleMask); + } + +#define WRITE_STAT(SYMBOL, NAME, SOURCE) \ + if (data->stats->options & SYMBOL) { \ + psMetadataAddF32(cellResults, PS_LIST_TAIL, NAME, 0, NULL, data->stats->SOURCE); \ + } + + WRITE_STAT(PS_STAT_SAMPLE_MEAN, "SAMPLE_MEAN", sampleMean); + WRITE_STAT(PS_STAT_SAMPLE_MEDIAN, "SAMPLE_MEDIAN", sampleMedian); + WRITE_STAT(PS_STAT_SAMPLE_STDEV, "SAMPLE_STDEV", sampleStdev); + WRITE_STAT(PS_STAT_SAMPLE_QUARTILE, "SAMPLE_LQ", sampleLQ); + WRITE_STAT(PS_STAT_SAMPLE_QUARTILE, "SAMPLE_UQ", sampleUQ); + WRITE_STAT(PS_STAT_ROBUST_MEDIAN, "ROBUST_MEDIAN", robustMedian); + WRITE_STAT(PS_STAT_ROBUST_STDEV, "ROBUST_STDEV", robustStdev); + WRITE_STAT(PS_STAT_ROBUST_QUARTILE, "ROBUST_LQ", robustLQ); + WRITE_STAT(PS_STAT_ROBUST_QUARTILE, "ROBUST_UQ", robustUQ); + WRITE_STAT(PS_STAT_ROBUST_QUARTILE, "ROBUST_N50", robustN50); + WRITE_STAT(PS_STAT_FITTED_MEAN, "FITTED_MEAN", fittedMean); + WRITE_STAT(PS_STAT_FITTED_STDEV, "FITTED_STDEV", fittedStdev); + WRITE_STAT(PS_STAT_CLIPPED_MEAN, "CLIPPED_MEAN", clippedMean); + WRITE_STAT(PS_STAT_CLIPPED_STDEV, "CLIPPED_STDEV", clippedStdev); + +cellDone: + // Add the cell results to the chip + addToHierarchy(cellResults, chipResults, cellName, "Results for cell"); + if (fits) { + pmCellFreeData(cell); + } + return; +} + +static bool chipStats(psMetadata *fpaResults, // Metadata holding the fpa results + pmChip *chip, // Chip for which to get statistics + psFits *fits, // FITS file handle + pmFPAview *view, // View for analysis + ppStatsData *data,// The data + const pmConfig *config // Configuration + ) +{ + assert(fpaResults); + assert(chip); + assert(view); + assert(data); + assert(config); + + const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip + + // Check to see if this is a chip of interest + if (!doThis(data->chips, chipName)) { + return true; + } + + psArray *cells = chip->cells; // Array of cells + if (view->cell >= cells->n) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Desired cell view (%d) doesn't match " + "number of cells (%d)\n", view->cell, cells->n); + return false; + } + + // Chip-level results + bool mdok; // Status of MD lookup + psMetadata *chipResults = psMemIncrRefCounter(psMetadataLookupMD(&mdok, fpaResults, chipName)); + if (!mdok || !chipResults) { + chipResults = psMetadataAlloc(); + } + + if (psListLength(data->headers) > 0 && chip->hdu) { + if (fits && !pmChipReadHeader(chip, fits)) { + goto chipDone; + } + pmHDU *hdu = chip->hdu; // HDU for headers + getMetadata(chipResults, hdu->header, data->headers); + } + if (psListLength(data->concepts) > 0) { + if (fits && !pmChipReadHeader(chip, fits)) { + goto chipDone; + } + pmConceptsReadChip(chip, PM_CONCEPT_SOURCE_ALL, false, false, config->database); + getMetadata(chipResults, chip->concepts, data->concepts); + } + + if (view->cell >= 0) { + pmCell *cell = cells->data[view->cell]; // Cell of interest + cellStats(chipResults, cell, fits, data, config); + goto chipDone; + } + + // Iterate over cells + for (int i = 0; i < cells->n; i++) { + pmCell *cell = cells->data[i]; // Cell of interest + cellStats(chipResults, cell, fits, data, config); + } + +chipDone: + addToHierarchy(chipResults, fpaResults, chipName, "Results for chip"); + if (fits) { + pmChipFreeData(chip); + } + + return true; +} + + psMetadata *ppStatsLoop(psMetadata *fpaResults, // Metadata to hold the FPA results ppStatsData *data, // The data @@ -72,13 +290,12 @@ PS_ASSERT_PTR_NON_NULL(fpa, NULL); + pmFPAview *view = psMemIncrRefCounter(data->view); // View for analysis + if (!view) { + view = pmFPAviewAlloc(0); + } + if (!fpaResults) { fpaResults = psMetadataAlloc(); } - - bool mdok; // Status of MD lookup - - // Iterators for the headers and concepts - psListIterator *headersIter = psListIteratorAlloc(data->headers, PS_LIST_HEAD, false); // Headers - psListIterator *conceptsIter = psListIteratorAlloc(data->concepts, PS_LIST_HEAD, false); // Concepts // Iterate through the FPA @@ -88,5 +305,5 @@ } pmHDU *hdu = fpa->hdu; // HDU for headers - getMetadata(fpaResults, hdu->header, headersIter); + getMetadata(fpaResults, hdu->header, data->headers); } if (psListLength(data->concepts) > 0) { @@ -95,208 +312,25 @@ } pmConceptsReadFPA(fpa, PM_CONCEPT_SOURCE_ALL, false, config->database); - getMetadata(fpaResults, fpa->concepts, conceptsIter); - } - psArray *chips = fpa->chips; // Array of component chips - for (long i = 0; i < chips->n; i++) { - pmChip *chip = chips->data[i]; // The chip of interest - if (!chip) { - continue; - } - const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip - - // Check to see if this is a chip of interest - if (!doThis(data->chips, chipName)) { - continue; - } - - // Chip-level results - psMetadata *chipResults = psMemIncrRefCounter(psMetadataLookupMD(&mdok, fpaResults, chipName)); - if (!mdok || !chipResults) { - chipResults = psMetadataAlloc(); - } - - if (psListLength(data->headers) > 0 && chip->hdu) { - if (fits && !pmChipReadHeader(chip, fits)) { - continue; - } - pmHDU *hdu = chip->hdu; // HDU for headers - getMetadata(chipResults, hdu->header, headersIter); - } - if (psListLength(data->concepts) > 0) { - if (fits && !pmChipReadHeader(chip, fits)) { - continue; - } - pmConceptsReadChip(chip, PM_CONCEPT_SOURCE_ALL, false, false, config->database); - getMetadata(chipResults, chip->concepts, conceptsIter); - } - - psArray *cells = chip->cells; // Array of component cells - for (long j = 0; j < cells->n; j++) { - pmCell *cell = cells->data[j]; // The cell of interest - if (!cell) { - continue; - } - const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell - - // Check to see if this is a cell of interest - if (!doThis(data->cells, cellName)) { - continue; - } - - // Cell-level results - psMetadata *cellResults = psMemIncrRefCounter(psMetadataLookupMD(&mdok, chipResults, cellName)); - if (!mdok || !cellResults) { - cellResults = psMetadataAlloc(); - } - - if (psListLength(data->headers) > 0 && cell->hdu) { - if (fits && !pmCellReadHeader(cell, fits)) { - addToHierarchy(cellResults, chipResults, cellName, "Results for cell"); - continue; - } - pmHDU *hdu = cell->hdu; // HDU for headers - getMetadata(cellResults, hdu->header, headersIter); - } - if (psListLength(data->concepts) > 0) { - if (fits && !pmCellReadHeader(cell, fits)) { - addToHierarchy(cellResults, chipResults, cellName, "Results for cell"); - continue; - } - pmConceptsReadCell(cell, PM_CONCEPT_SOURCE_ALL, false, config->database); - getMetadata(cellResults, cell->concepts, conceptsIter); - } - - if (!data->doStats) { - // Nothing further to do --- don't want to waste our time reading the data - if (psListLength(cellResults->list) > 0) { - psMetadataAdd(chipResults, PS_LIST_TAIL, cellName, PS_DATA_METADATA, - "Results for cell", cellResults); - } - addToHierarchy(cellResults, chipResults, cellName, "Results for cell"); - continue; - } - - pmHDU *hdu = pmHDUFromCell(cell); // HDU for cell - if (!hdu || hdu->blankPHU) { - addToHierarchy(cellResults, chipResults, cellName, "Results for cell"); - continue; - } - - if (fits && !pmCellRead(cell, fits, config->database)) { - psLogMsg(__func__, PS_LOG_WARN, "Unable to read chip %s cell %s\n", chipName, cellName); - pmCellFreeData(cell); - addToHierarchy(cellResults, chipResults, cellName, "Results for cell"); - continue; - } - - psArray *readouts = cell->readouts; // Array of component readouts - if (readouts->n == 0) { - if (fits) { - pmCellFreeData(cell); - } - addToHierarchy(cellResults, chipResults, cellName, "Results for cell"); - continue; - } - if (readouts->n > 1) { - psLogMsg(__func__, PS_LOG_WARN, "Multiple readouts (%ld) present in chip %s cell %s --- " - "using only the first.\n", readouts->n, chipName, cellName); - } - pmReadout *readout = readouts->data[0]; // The readout of interest - if (!readout->image) { - psLogMsg(__func__, PS_LOG_WARN, "No image associated with readout in chip %s cell %s --- " - "ignored.\n", chipName, cellName); - if (fits) { - pmCellFreeData(cell); - } - addToHierarchy(cellResults, chipResults, cellName, "Results for cell"); - continue; - } - - // Do the statistics - if (data->sample <= 0.0) { - if (!psImageStats(data->stats, readout->image, readout->mask, data->maskVal)) { - psLogMsg(__func__, PS_LOG_WARN, "Unable to perform statistics on chip %s cell %s --- " - "ignored.\n", chipName, cellName); - if (fits) { - pmCellFreeData(cell); - } - addToHierarchy(cellResults, chipResults, cellName, "Results for cell"); - continue; - } - } else { - // Apply sampling - psImage *image = readout->image; // The image of interest - psImage *mask = readout->mask; // The mask image - int numSamples = data->sample * image->numCols * image->numRows; // Number of samples - int sampleSpace = 1.0 / data->sample; // Space between samples - psVector *sampleValues = psVectorAlloc(numSamples, PS_TYPE_F32); // Vector of samples - sampleValues->n = numSamples; - psVector *sampleMask = NULL; // Corresponding mask - if (mask) { - sampleMask = psVectorAlloc(numSamples, PS_TYPE_U8); - sampleMask->n = numSamples; - } - for (int i = 0; i < numSamples; i++) { - int j = i * sampleSpace; - int y = j / image->numRows; - int x = j % image->numRows; - sampleValues->data.F32[i] = image->data.F32[y][x]; - if (mask) { - sampleMask->data.U8[i] = mask->data.U8[y][x]; - } - } - if (!psVectorStats(data->stats, sampleValues, NULL, sampleMask, data->maskVal)) { - psLogMsg(__func__, PS_LOG_WARN, "Unable to perform statistics on chip %s cell %s --- " - "ignored.\n", chipName, cellName); - psFree(sampleValues); - psFree(sampleMask); - if (fits) { - pmCellFreeData(cell); - } - addToHierarchy(cellResults, chipResults, cellName, "Results for cell"); - continue; - } - psFree(sampleValues); - psFree(sampleMask); - } - - #define WRITE_STAT(SYMBOL, NAME, SOURCE) \ - if (data->stats->options & SYMBOL) { \ - psMetadataAddF32(cellResults, PS_LIST_TAIL, NAME, 0, NULL, data->stats->SOURCE); \ - } - - WRITE_STAT(PS_STAT_SAMPLE_MEAN, "SAMPLE_MEAN", sampleMean); - WRITE_STAT(PS_STAT_SAMPLE_MEDIAN, "SAMPLE_MEDIAN", sampleMedian); - WRITE_STAT(PS_STAT_SAMPLE_STDEV, "SAMPLE_STDEV", sampleStdev); - WRITE_STAT(PS_STAT_SAMPLE_QUARTILE, "SAMPLE_LQ", sampleLQ); - WRITE_STAT(PS_STAT_SAMPLE_QUARTILE, "SAMPLE_UQ", sampleUQ); - WRITE_STAT(PS_STAT_ROBUST_MEDIAN, "ROBUST_MEDIAN", robustMedian); - WRITE_STAT(PS_STAT_ROBUST_STDEV, "ROBUST_STDEV", robustStdev); - WRITE_STAT(PS_STAT_ROBUST_QUARTILE, "ROBUST_LQ", robustLQ); - WRITE_STAT(PS_STAT_ROBUST_QUARTILE, "ROBUST_UQ", robustUQ); - WRITE_STAT(PS_STAT_ROBUST_QUARTILE, "ROBUST_N50", robustN50); - WRITE_STAT(PS_STAT_FITTED_MEAN, "FITTED_MEAN", fittedMean); - WRITE_STAT(PS_STAT_FITTED_STDEV, "FITTED_STDEV", fittedStdev); - WRITE_STAT(PS_STAT_CLIPPED_MEAN, "CLIPPED_MEAN", clippedMean); - WRITE_STAT(PS_STAT_CLIPPED_STDEV, "CLIPPED_STDEV", clippedStdev); - - // Add the cell results to the chip - addToHierarchy(cellResults, chipResults, cellName, "Results for cell"); - if (fits) { - pmCellFreeData(cell); - } - } - - addToHierarchy(chipResults, fpaResults, chipName, "Results for chip"); - if (fits) { - pmChipFreeData(chip); - } - - } + getMetadata(fpaResults, fpa->concepts, data->concepts); + } + + psArray *chips = fpa->chips; // Array of chips + if (view->chip >= 0) { + pmChip *chip = chips->data[view->chip]; // Chip of interest + chipStats(fpaResults, chip, fits, view, data, config); + goto fpaDone; + } + + // Iterate over cells + for (int i = 0; i < chips->n; i++) { + pmChip *chip = chips->data[i]; // Chip of interest + chipStats(fpaResults, chip, fits, view, data, config); + } + +fpaDone: if (fits) { pmFPAFreeData(fpa); } - psFree(headersIter); - psFree(conceptsIter); + psFree(view); return fpaResults;