Changeset 8346
- Timestamp:
- Aug 14, 2006, 7:13:32 PM (20 years ago)
- Location:
- trunk/ppStats/src
- Files:
-
- 5 edited
-
ppStats.c (modified) (2 diffs)
-
ppStats.h (modified) (1 diff)
-
ppStatsData.c (modified) (2 diffs)
-
ppStatsData.h (modified) (1 diff)
-
ppStatsLoop.c (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/ppStats/src/ppStats.c
r8337 r8346 7 7 psMetadata *ppStats(psMetadata *out, // Output metadata 8 8 pmFPA *fpa, // FPA for which to get statistics 9 pmFPAview *view, // View for analysis 9 10 pmConfig *config // Configuration 10 11 ) … … 18 19 data->fpa = psMemIncrRefCounter(fpa); 19 20 21 if (data->view) { 22 psFree(data->view); 23 } 24 data->view = psMemIncrRefCounter(view); 25 20 26 // Go through the FPA and do the hard work 21 27 out = ppStatsLoop(out, data, config); -
trunk/ppStats/src/ppStats.h
r8337 r8346 10 10 #include "ppStatsLoop.h" 11 11 12 // Perform the ppStats steps 13 psMetadata *ppStats(psMetadata *out, // Output metadata 14 pmFPA *fpa, // FPA for which to get statistics 15 pmFPAview *view, // View for analysis 16 pmConfig *config // Configuration 17 ); 18 12 19 #endif -
trunk/ppStats/src/ppStatsData.c
r8337 r8346 7 7 ) 8 8 { 9 // inName and region are not on the psLib memory system (they are from argv).10 psFree(data->fpa);11 9 if (data->fits) { 12 10 psFitsClose(data->fits); 13 11 data->fits = NULL; 14 12 } 13 psFree(data->fpa); 14 psFree(data->view); 15 15 psFree(data->headers); 16 16 psFree(data->concepts); … … 29 29 psMemSetDeallocator(data, (psFreeFunc)statsDataFree); 30 30 31 data->fits = NULL; 31 32 data->fpa = NULL; 32 data-> fits= NULL;33 data->view = NULL; 33 34 34 35 data->headers = psListAlloc(NULL); -
trunk/ppStats/src/ppStatsData.h
r8337 r8346 9 9 psFits *fits; // Input file handle 10 10 pmFPA *fpa; // FPA to analyse 11 pmFPAview *view; // View to analyse 11 12 // Stuff to output 12 13 psStats *stats; // Statistics to calculate -
trunk/ppStats/src/ppStatsLoop.c
r8338 r8346 1 1 #include <stdio.h> 2 #include <assert.h> 2 3 #include <string.h> 3 4 #include <pslib.h> … … 10 11 static void getMetadata(psMetadata *target, // Target for metadata 11 12 psMetadata *source, // Source for metadata 12 psList Iterator *iterator // Iterator forkeywords13 ) 14 { 15 psListIterator Set(iterator, PS_LIST_HEAD);13 psList *list // List containing keywords 14 ) 15 { 16 psListIterator *iterator = psListIteratorAlloc(list, PS_LIST_HEAD, false); // Iterator 16 17 psString name; // Name from iteration 17 18 while ((name = psListGetAndIncrement(iterator))) { … … 21 22 } 22 23 } 24 psFree(iterator); 23 25 return; 24 26 } … … 62 64 63 65 66 static 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 201 cellDone: 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 210 static 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 272 chipDone: 273 addToHierarchy(chipResults, fpaResults, chipName, "Results for chip"); 274 if (fits) { 275 pmChipFreeData(chip); 276 } 277 278 return true; 279 } 280 281 64 282 psMetadata *ppStatsLoop(psMetadata *fpaResults, // Metadata to hold the FPA results 65 283 ppStatsData *data, // The data … … 72 290 PS_ASSERT_PTR_NON_NULL(fpa, NULL); 73 291 292 pmFPAview *view = psMemIncrRefCounter(data->view); // View for analysis 293 if (!view) { 294 view = pmFPAviewAlloc(0); 295 } 296 74 297 if (!fpaResults) { 75 298 fpaResults = psMetadataAlloc(); 76 299 } 77 78 bool mdok; // Status of MD lookup79 80 // Iterators for the headers and concepts81 psListIterator *headersIter = psListIteratorAlloc(data->headers, PS_LIST_HEAD, false); // Headers82 psListIterator *conceptsIter = psListIteratorAlloc(data->concepts, PS_LIST_HEAD, false); // Concepts83 300 84 301 // Iterate through the FPA … … 88 305 } 89 306 pmHDU *hdu = fpa->hdu; // HDU for headers 90 getMetadata(fpaResults, hdu->header, headersIter);307 getMetadata(fpaResults, hdu->header, data->headers); 91 308 } 92 309 if (psListLength(data->concepts) > 0) { … … 95 312 } 96 313 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 330 fpaDone: 296 331 if (fits) { 297 332 pmFPAFreeData(fpa); 298 333 } 299 psFree(headersIter); 300 psFree(conceptsIter); 334 psFree(view); 301 335 302 336 return fpaResults;
Note:
See TracChangeset
for help on using the changeset viewer.
