Changeset 18859
- Timestamp:
- Aug 1, 2008, 1:58:14 PM (18 years ago)
- Location:
- trunk/psModules/src/detrend
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/detrend/pmDark.c
r18163 r18859 98 98 } 99 99 100 101 bool pmDarkCombine(pmCell *output, const psArray *inputs, psArray *ordinates, const char *normConcept, 102 int iter, float rej, psMaskType maskVal) 103 { 104 PS_ASSERT_PTR_NON_NULL(output, false); 105 PS_ASSERT_ARRAY_NON_NULL(inputs, false); 106 PS_ASSERT_ARRAY_NON_NULL(ordinates, false); 107 PS_ASSERT_INT_NONNEGATIVE(iter, false); 108 PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false); 109 110 bool inRange = true; 111 112 // Extract fitting orders 113 int numOrdinates = ordinates->n; // Number of ordinates 114 int numInputs = inputs->n; // Number of inputs 100 // this creates and saves: values, roMask, norm, poly, counts, sigma, and saves the on output->analysis 101 bool pmDarkCombinePrepare(pmCell *output, const psArray *inputs, psArray *ordinates, const char *normConcept) 102 { 103 psArray *values = psArrayAlloc(inputs->n); 104 psVector *roMask = psVectorAlloc(inputs->n, PS_TYPE_U8); // Mask for bad readouts 105 psVector *norm = normConcept ? psVectorAlloc(inputs->n, PS_TYPE_F32) : NULL; // Normalizations for each 106 psVector *orders = psVectorAlloc(ordinates->n, PS_TYPE_U8); // Orders for each concept 107 108 psVectorInit(roMask, 0); 109 110 bool inRange = false; 115 111 int numBadInputs = 0; // Number of bad inputs 116 psArray *values = psArrayAlloc(numInputs); 117 psVector *roMask = psVectorAlloc(numInputs, PS_TYPE_U8); // Mask for bad readouts 118 psVectorInit(roMask, 0); 119 psVector *norm = normConcept ? psVectorAlloc(numInputs, PS_TYPE_F32) : NULL; // Normalisations for each 120 for (int i = 0; i < numInputs; i++) { 121 values->data[i] = psVectorAlloc(numOrdinates, PS_TYPE_F32); 122 if (norm) { 123 pmReadout *ro = inputs->data[i]; // Readout of interest 124 float normValue; // Normalisation value 125 if (!ordinateLookup(&normValue, &inRange, normConcept, false, NAN, NAN, ro)) { 126 psWarning("Unable to find value of %s for readout %d", normConcept, i); 127 roMask->data.U8[i] = 0xff; 128 norm->data.F32[i] = NAN; 129 numBadInputs++; 130 continue; 131 } 132 if (normValue == 0.0) { 133 psWarning("Normalisation value (%s) for readout %d is zero", normConcept, i); 134 roMask->data.U8[i] = 0xff; 135 norm->data.F32[i] = NAN; 136 numBadInputs++; 137 continue; 138 } 139 norm->data.F32[i] = 1.0 / normValue; 140 } 141 } 142 psVector *orders = psVectorAlloc(numOrdinates, PS_TYPE_U8); // Orders for each concept 143 for (int i = 0; i < numOrdinates; i++) { 112 113 // build the 'norm' vector and the 'values' vectors, count the number of bad inputs 114 for (int i = 0; i < inputs->n; i++) { 115 values->data[i] = psVectorAlloc(ordinates->n, PS_TYPE_F32); 116 if (!norm) continue; 117 118 pmReadout *readout = inputs->data[i]; // Readout of interest 119 float normValue; // Normalisation value 120 if (!ordinateLookup(&normValue, &inRange, normConcept, false, NAN, NAN, readout)) { 121 psWarning("Unable to find value of %s for readout %d", normConcept, i); 122 roMask->data.U8[i] = 0xff; 123 norm->data.F32[i] = NAN; 124 numBadInputs++; 125 continue; 126 } 127 if (normValue == 0.0) { 128 psWarning("Normalisation value (%s) for readout %d is zero", normConcept, i); 129 roMask->data.U8[i] = 0xff; 130 norm->data.F32[i] = NAN; 131 numBadInputs++; 132 continue; 133 } 134 norm->data.F32[i] = 1.0 / normValue; 135 } 136 137 // build the 'orders' vector and set the array of 'values' 138 for (int i = 0; i < ordinates->n; i++) { 144 139 pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate information 145 140 if (ord->order <= 0) { … … 148 143 psFree(roMask); 149 144 psFree(orders); 145 psFree(norm); 150 146 return false; 151 147 } 152 148 orders->data.U8[i] = ord->order; 153 149 154 // Mask the readout and move on 155 #define MASK_READOUT_VALUE { \ 156 roMask->data.U8[j] = 0xff; \ 157 val->data.F32[i] = NAN; \ 158 numBadInputs++; \ 159 continue; \ 160 } 161 162 for (int j = 0; j < numInputs; j++) { 150 for (int j = 0; j < inputs->n; j++) { 163 151 psVector *val = values->data[j]; // Value vector for readout 164 152 if (roMask->data.U8[j]) { … … 167 155 } 168 156 169 pmReadout *r o= inputs->data[j]; // Readout of interest157 pmReadout *readout = inputs->data[j]; // Readout of interest 170 158 float value = NAN; // Value of ordinate 171 if (!ordinateLookup(&value, &inRange, ord->name, ord->scale, ord->min, ord->max, r o)) {159 if (!ordinateLookup(&value, &inRange, ord->name, ord->scale, ord->min, ord->max, readout)) { 172 160 roMask->data.U8[j] = 0xff; 173 161 val->data.F32[i] = NAN; … … 186 174 187 175 if (psTraceGetLevel("psModules.detrend") > 9) { 188 for (int i = 0; i < numInputs; i++) {176 for (int i = 0; i < inputs->n; i++) { 189 177 psVector *val = values->data[i]; 190 178 (void) val; // avoid unused variable message when tracing is compiled out 191 for (int j = 0; j < numOrdinates; j++) {179 for (int j = 0; j < ordinates->n; j++) { 192 180 psTrace("psModules.detrend", 9, "Image %d, ordinate %d: %f\n", i, j, val->data.F32[j]); 193 181 } 194 182 } 195 183 } 196 197 184 psPolynomialMD *poly = psPolynomialMDAlloc(orders); // Polynomial for fitting 198 185 psFree(orders); 186 199 187 int numTerms = poly->coeff->n; // Number of terms in polynomial 200 if (numTerms > numInputs- numBadInputs) {201 psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Insufficient inputs (% d) to fit polynomial terms (%d).",202 numInputs- numBadInputs, numTerms);188 if (numTerms > inputs->n - numBadInputs) { 189 psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Insufficient inputs (%ld) to fit polynomial terms (%d).", 190 inputs->n - numBadInputs, numTerms); 203 191 psFree(values); 204 192 psFree(roMask); … … 207 195 } 208 196 209 // Set up output 197 // determine the output image size based on the input images 198 int row0, col0, numCols, numRows; 199 if (!pmReadoutStackSetOutputSize(&col0, &row0, &numCols, &numRows, inputs)) { 200 psError(PS_ERR_UNKNOWN, false, "problem setting output readout size."); 201 return false; 202 } 203 204 // the output is potentially a cube (depending on the dimensionality of the fit) 205 if (output->readouts->n != numTerms) { 206 output->readouts = psArrayRealloc(output->readouts, numTerms); 207 } 208 209 // generate the required output images based on the specified sizes 210 for (int i = 0; i < numTerms; i++) { 211 pmReadout *readout = output->readouts->data[i]; // Readout to update 212 if (!readout) { 213 readout = output->readouts->data[i] = pmReadoutAlloc(output); 214 } 215 216 pmReadoutStackDefineOutput(readout, col0, row0, numCols, numRows, false, false, 0); 217 psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", col0, row0); 218 } 219 220 // these calls allocate and save the requested images on the output analysis metadata 221 psImage *counts = pmReadoutSetAnalysisImage(output->readouts->data[0], PM_READOUT_STACK_ANALYSIS_COUNT, numCols, numRows, PS_TYPE_U16, 0); 222 if (!counts) { 223 return false; 224 } 225 psImage *sigma = pmReadoutSetAnalysisImage(output->readouts->data[0], PM_READOUT_STACK_ANALYSIS_SIGMA, numCols, numRows, PS_TYPE_F32, NAN); 226 if (!sigma) { 227 return false; 228 } 229 230 psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_ORDINATES, PS_DATA_ARRAY | PS_META_REPLACE, "Dark ordinates", ordinates); 231 psMetadataAddStr(output->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_NORM, PS_META_REPLACE, "Dark normalisation", normConcept); 232 233 psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "DARK.VALUES", PS_DATA_ARRAY | PS_META_REPLACE, "Dark values", values); 234 psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "DARK.RO.MASK", PS_DATA_VECTOR | PS_META_REPLACE, "Dark Readout Mask", roMask); 235 psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "DARK.NORM", PS_DATA_VECTOR | PS_META_REPLACE, "Dark norm", norm); 236 psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "DARK.POLY", PS_DATA_UNKNOWN | PS_META_REPLACE, "Dark poly", poly); 237 238 for (int i = 0; i < numTerms; i++) { 239 pmReadout *readout = output->readouts->data[i]; // Readout to update 240 readout->data_exists = true; 241 } 242 output->data_exists = true; 243 output->parent->data_exists = true; 244 245 psFree(norm); 246 psFree(roMask); 247 psFree(poly); 248 psFree(values); 249 250 return true; 251 } 252 253 // do the combine work for this portion of the output (range is set by input data) 254 bool pmDarkCombine(pmCell *output, const psArray *inputs, int iter, float rej, psMaskType maskVal) 255 { 256 PS_ASSERT_PTR_NON_NULL(output, false); 257 PS_ASSERT_PTR_NON_NULL(output->readouts, false); 258 PS_ASSERT_INT_NONNEGATIVE(output->readouts->n, false); 259 PS_ASSERT_ARRAY_NON_NULL(inputs, false); 260 PS_ASSERT_INT_NONNEGATIVE(iter, false); 261 PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false); 262 263 pthread_t id = pthread_self(); 264 char name[64]; 265 sprintf (name, "%x", (unsigned int) id); 266 psTimerStart (name); 267 268 bool mdok = false; 269 270 // retrieve the required parameter vectors 271 psArray *values = psMetadataLookupPtr(&mdok, output->analysis, "DARK.VALUES"); psAssert (values, "values not supplied"); 272 psVector *roMask = psMetadataLookupPtr(&mdok, output->analysis, "DARK.RO.MASK"); psAssert (roMask, "roMask not supplied"); 273 psPolynomialMD *poly = psMetadataLookupPtr(&mdok, output->analysis, "DARK.POLY"); psAssert (poly, "orders not supplied"); 274 275 // retrieve the norm vector, if supplied 276 psVector *norm = psMetadataLookupPtr(&mdok, output->analysis, "DARK.NORM"); 277 278 // retrieve the 'counts' and 'sigma' images 279 psImage *counts = pmReadoutGetAnalysisImage(output->readouts->data[0], PM_READOUT_STACK_ANALYSIS_COUNT); 280 if (!counts) { 281 return false; 282 } 283 psImage *sigma = pmReadoutGetAnalysisImage(output->readouts->data[0], PM_READOUT_STACK_ANALYSIS_SIGMA); 284 if (!sigma) { 285 return false; 286 } 287 210 288 int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine 211 289 int xSize, ySize; // Size of the output image … … 213 291 inputs)) { 214 292 psError(PS_ERR_UNKNOWN, false, "No valid input readouts."); 215 psFree(values); 216 psFree(roMask); 217 psFree(norm); 218 return false; 219 } 220 if (output->readouts->n != numTerms) { 221 output->readouts = psArrayRealloc(output->readouts, numTerms); 222 } 223 int outRow0 = 0, outCol0 = 0; // Output row0 and col0 224 for (int i = 0; i < numTerms; i++) { 225 pmReadout *ro = output->readouts->data[i]; // Readout to update 226 if (!ro) { 227 ro = output->readouts->data[i] = pmReadoutAlloc(output); 228 } 229 230 pmReadoutUpdateSize(ro, minInputCols, minInputRows, xSize, ySize, false, false, 0); 231 if (i == 0) { 232 outRow0 = ro->row0; 233 outCol0 = ro->col0; 234 } else { 235 assert(ro->row0 == outRow0); 236 assert(ro->col0 == outCol0); 237 } 238 } 239 240 psImage *counts = pmReadoutAnalysisImage(output->readouts->data[0], PM_READOUT_STACK_ANALYSIS_COUNT, 241 xSize, ySize, PS_TYPE_U16, 0); 242 if (!counts) { 243 return false; 244 } 245 psImage *sigma = pmReadoutAnalysisImage(output->readouts->data[0], PM_READOUT_STACK_ANALYSIS_SIGMA, 246 xSize, ySize, PS_TYPE_F32, NAN); 247 if (!sigma) { 248 psFree(counts); 249 return false; 250 } 293 return false; 294 } 295 296 pmReadout *outReadout = output->readouts->data[0]; 251 297 252 298 // Iterate over pixels, fitting polynomial 253 psVector *pixels = psVectorAlloc( numInputs, PS_TYPE_F32); // Stack of pixels254 psVector *mask = psVectorAlloc( numInputs, PS_TYPE_MASK); // Mask for stack299 psVector *pixels = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of pixels 300 psVector *mask = psVectorAlloc(inputs->n, PS_TYPE_MASK); // Mask for stack 255 301 for (int i = minInputRows; i < maxInputRows; i++) { 256 int yOut = i - outRow0; // y position on output readout 257 #ifdef SHOW_BUSY 302 int yOut = i - outReadout->row0; // y position on output readout 303 304 # ifdef SHOW_BUSY 258 305 if (psTraceGetLevel("psModules.detrend") > 9) { 259 306 printf("Processing row %d\r", i); 260 307 fflush(stdout); 261 308 } 262 # endif309 # endif 263 310 264 311 for (int j = minInputCols; j < maxInputCols; j++) { 265 int xOut = j - out Col0; // x position on output readout312 int xOut = j - outReadout->col0; // x position on output readout 266 313 267 314 psVectorInit(mask, 0); … … 303 350 psVectorInit(poly->coeff, NAN); 304 351 } 305 for (int k = 0; k < numTerms; k++) {352 for (int k = 0; k < poly->coeff->n; k++) { 306 353 pmReadout *ro = output->readouts->data[k]; // Readout of interest 307 354 ro->image->data.F32[yOut][xOut] = poly->coeff->data.F64[k]; … … 312 359 } 313 360 314 psFree(norm);315 psFree(roMask);316 psFree(poly);317 361 psFree(pixels); 318 362 psFree(mask); 319 psFree(values); 320 321 psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_ORDINATES, 322 PS_DATA_ARRAY | PS_META_REPLACE, "Dark ordinates", ordinates); 323 psMetadataAddStr(output->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_NORM, 324 PS_META_REPLACE, "Dark normalisation", normConcept); 325 326 for (int i = 0; i < numTerms; i++) { 327 pmReadout *ro = output->readouts->data[i]; // Readout to update 328 ro->data_exists = true; 329 } 330 output->data_exists = true; 331 output->parent->data_exists = true; 332 363 364 fprintf (stderr, "done with combine %x : %f sec\n", (unsigned int) id, psTimerMark (name)); 333 365 return true; 334 366 } 335 336 337 367 338 368 bool pmDarkApply(pmReadout *readout, const pmCell *dark, psMaskType bad) -
trunk/psModules/src/detrend/pmDark.h
r18163 r18859 25 25 26 26 27 // Combine darks 27 // Combine darks -- preparation step 28 bool pmDarkCombinePrepare(pmCell *output, // Output cell; readouts will be attached 29 const psArray *inputs, // Input readouts for combination 30 psArray *ordinates, // Ordinates for fitting 31 const char *normConcept // Concept name to use to divide input pixel values 32 ); 33 34 // Combine darks -- do the actual work 28 35 bool pmDarkCombine(pmCell *output, // Output cell; readouts will be attached 29 36 const psArray *inputs, // Input readouts for combination 30 psArray *ordinates, // Ordinates for fitting31 const char *normConcept, // Concept name to use to divide input pixel values32 37 int iter, // Number of rejection iterations 33 38 float rej, // Rejection threshold (standard deviations)
Note:
See TracChangeset
for help on using the changeset viewer.
