Changeset 7059 for trunk/psModules/src/imcombine/pmReadoutCombine.c
- Timestamp:
- May 3, 2006, 4:38:20 PM (20 years ago)
- File:
-
- 1 edited
-
trunk/psModules/src/imcombine/pmReadoutCombine.c (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/pmReadoutCombine.c
r6873 r7059 5 5 * @author GLG, MHPCC 6 6 * 7 * @version $Revision: 1. 7$ $Name: not supported by cvs2svn $8 * @date $Date: 2006-0 4-17 18:10:08$7 * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2006-05-04 02:38:20 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 12 12 */ 13 13 14 #include<stdio.h> 15 #include<math.h> 16 #include "pslib.h" 14 #include <stdio.h> 15 #include <assert.h> 16 #include <pslib.h> 17 18 #include "pmFPA.h" 19 #include "pmMaskBadPixels.h" 20 17 21 #include "pmReadoutCombine.h" 18 22 23 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 24 // File-static (private) functions 25 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 26 27 // Return the statistic of interest 28 static double getStat(const psStats *stats, // Statistics structure 29 psStatsOptions option // Statistics option 30 ) 31 { 32 switch (option) { 33 case PS_STAT_SAMPLE_MEAN: 34 return stats->sampleMean; 35 case PS_STAT_SAMPLE_MEDIAN: 36 return stats->sampleMedian; 37 case PS_STAT_SAMPLE_STDEV: 38 return stats->sampleStdev; 39 case PS_STAT_ROBUST_MEDIAN: 40 return stats->robustMedian; 41 case PS_STAT_ROBUST_STDEV: 42 return stats->robustStdev; 43 case PS_STAT_FITTED_MEAN: 44 return stats->fittedMean; 45 case PS_STAT_FITTED_STDEV: 46 return stats->fittedStdev; 47 case PS_STAT_CLIPPED_MEAN: 48 return stats->clippedMean; 49 case PS_STAT_CLIPPED_STDEV: 50 return stats->clippedStdev; 51 case PS_STAT_MAX: 52 return stats->max; 53 case PS_STAT_MIN: 54 return stats->min; 55 default: 56 psAbort(__func__, "Bad option: %x\n", option); 57 } 58 return NAN; 59 } 60 61 // Mask for combination --- used only locally 62 typedef enum { 63 PM_READOUT_COMBINE_CLEAR = 0x00, // No problem 64 PM_READOUT_COMBINE_NO_IMAGE = 0x01, // No image available 65 PM_READOUT_COMBINE_MASKED = 0x02, // Pixel is masked 66 PM_READOUT_COMBINE_BAD = 0x03, // Pixel is bad (NO_IMAGE | MASKED) 67 PM_READOUT_COMBINE_CLIPPED = 0x04 // Pixel has been clipped 68 } pmReadoutCombineMask; 69 70 71 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 72 // Public functions 73 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 74 75 // Allocator for pmCombineParams 76 pmCombineParams *pmCombineParamsAlloc(psStatsOptions combine) 77 { 78 pmCombineParams *params = psAlloc(sizeof(pmCombineParams)); 79 80 params->combine = combine; 81 params->maskVal = 0; 82 params->nKeep = 0; 83 params->fracHigh = 0.0; 84 params->fracHigh = 0.0; 85 params->iter = 1; 86 params->rej = 3.0; 87 88 return params; 89 } 90 91 19 92 /****************************************************************************** 20 DetermineNumBits(data): This routine takes an enum psStatsOptions as an 21 argument and returns the number of non-zero bits. 93 XXX: Maybe add support for S16 and S32 types. Currently, only F32 supported. 22 94 *****************************************************************************/ 23 static psS32 DetermineNumBits(psStatsOptions data) 95 bool pmReadoutCombine(pmReadout *output, 96 const psArray *inputs, 97 const psVector *zero, 98 const psVector *scale, 99 pmCombineParams *params 100 ) 24 101 { 25 psS32 i; 26 psU64 tmpData = data; 27 psS32 numBits = 0; 28 29 for (i=0;i<(8 * sizeof(psStatsOptions));i++) { 30 if (0x0001 & tmpData) { 31 numBits++; 32 } 33 tmpData = tmpData >> 1; 34 } 35 return(numBits); 102 // Check inputs 103 PS_ASSERT_PTR_NON_NULL(output, false); 104 PS_ASSERT_PTR_NON_NULL(inputs, false); 105 PS_ASSERT_PTR_NON_NULL(params, false); 106 if (zero) { 107 PS_ASSERT_VECTOR_TYPE(zero, PS_TYPE_F32, false); 108 if (zero->n != inputs->n) { 109 psError(PS_ERR_UNKNOWN, true, "Zero vector has incorrect size (%d vs %d).\n", 110 zero->n, inputs->n); 111 return false; 112 } 113 } 114 if (scale) { 115 PS_ASSERT_VECTOR_TYPE(scale, PS_TYPE_F32, false); 116 if (scale->n != inputs->n) { 117 psError(PS_ERR_UNKNOWN, true, "Scale vector has incorrect size (%d vs %d).\n", 118 scale->n, inputs->n); 119 return false; 120 } 121 } 122 assert(params->fracLow >= 0.0 && params->fracLow < 1.0); 123 assert(params->fracHigh >= 0.0 && params->fracHigh < 1.0); 124 assert(params->combine == PS_STAT_SAMPLE_MEAN || 125 params->combine == PS_STAT_SAMPLE_MEDIAN || 126 params->combine == PS_STAT_ROBUST_MEDIAN || 127 params->combine == PS_STAT_FITTED_MEAN || 128 params->combine == PS_STAT_CLIPPED_MEAN); 129 130 psStats *stats = psStatsAlloc(params->combine); // The statistics to use in the combination 131 if (params->combine == PS_STAT_CLIPPED_MEAN) { 132 stats->clipSigma = params->rej; 133 stats->clipIter = params->iter; 134 } 135 136 // Step through each readout in the input image list to determine how big of an output image is needed to 137 // combine these input images. 138 long maxInputCols = 0; // The largest input column value 139 long maxInputRows = 0; // The largest input row value 140 long minInputCols = LONG_MAX; // The smallest input column value 141 long minInputRows = LONG_MAX; // The smallest input row value 142 psVector *rowLower = psVectorAlloc(inputs->n, PS_TYPE_U32); // The lower y bound for each image 143 psVector *rowUpper = psVectorAlloc(inputs->n, PS_TYPE_U32); // The upper y bound for each image 144 psVector *colLower = psVectorAlloc(inputs->n, PS_TYPE_U32); // The lower x bound for each image 145 psVector *colUpper = psVectorAlloc(inputs->n, PS_TYPE_U32); // The upper x bound for each image 146 psVector *mask = psVectorAlloc(inputs->n, PS_TYPE_U8); // Mask for stack 147 psVectorInit(mask, 0); 148 bool haveWeights = false; // Do we have weight images? 149 bool valid = false; // Do we have a single valid input? 150 for (int i = 0; i < inputs->n; i++) { 151 pmReadout *readout = inputs->data[i]; // Readout of interest 152 if (!readout || !readout->image) { 153 mask->data.U8[i] = PM_READOUT_COMBINE_NO_IMAGE; 154 continue; 155 } 156 157 if (readout->weight) { 158 if (valid && !haveWeights) { 159 psLogMsg(__func__, PS_LOG_WARN, "Readout %d has a weight map, but others don't --- " 160 "weights ignored.\n", i); 161 } else { 162 haveWeights = true; 163 } 164 } else if (haveWeights) { 165 psLogMsg(__func__, PS_LOG_WARN, "Readout %d doesn't have a weight map, but others do --- " 166 "weights ignored.\n", i); 167 haveWeights = false; 168 } 169 valid = true; 170 171 // Size of output image 172 minInputRows = PS_MIN(minInputRows, readout->row0); 173 maxInputRows = PS_MAX(maxInputRows, readout->row0 + readout->image->numRows); 174 minInputCols = PS_MIN(minInputCols, readout->col0); 175 maxInputCols = PS_MAX(maxInputCols, readout->col0 + readout->image->numCols); 176 // Bounds of input image 177 rowLower->data.U32[i] = readout->row0; 178 colLower->data.U32[i] = readout->col0; 179 rowUpper->data.U32[i] = readout->row0 + readout->image->numRows; 180 colUpper->data.U32[i] = readout->col0 + readout->image->numCols; 181 } 182 183 // Reset output readout components 184 *(psS32 *) &(output->col0) = minInputCols; 185 *(psS32 *) &(output->row0) = minInputRows; 186 output->image = psImageRecycle(output->image, maxInputCols - minInputCols, maxInputRows - minInputRows, 187 PS_TYPE_F32); 188 output->mask = psImageRecycle(output->mask, maxInputCols - minInputCols, maxInputRows - minInputRows, 189 PS_TYPE_U8); 190 psImageInit(output->mask, 0); 191 psStatsOptions combineStdev = 0; // Statistics option for weights 192 if (haveWeights) { 193 output->weight = psImageRecycle(output->weight, maxInputCols - minInputCols, 194 maxInputRows - minInputRows, PS_TYPE_F32); 195 switch (params->combine) { 196 case PS_STAT_SAMPLE_MEAN: 197 case PS_STAT_SAMPLE_MEDIAN: 198 combineStdev = PS_STAT_SAMPLE_STDEV; 199 break; 200 case PS_STAT_ROBUST_MEDIAN: 201 combineStdev = PS_STAT_ROBUST_STDEV; 202 break; 203 case PS_STAT_FITTED_MEAN: 204 combineStdev = PS_STAT_FITTED_STDEV; 205 break; 206 case PS_STAT_CLIPPED_MEAN: 207 combineStdev = PS_STAT_CLIPPED_STDEV; 208 break; 209 default: 210 psAbort(__func__, "Should never get here --- checked params->combine before.\n"); 211 } 212 stats->options |= combineStdev; 213 } 214 215 // We loop through each pixel in the output image. We loop through each input readout. We determine if 216 // that output pixel is contained in the image from that readout. If so, we save it in psVector 217 // tmpPixels. If not, we set a mask for that element in tmpPixels. Then, we mask off pixels not between 218 // fracLow and fracHigh. Then we call the vector stats routine on those pixels/mask. Then we set the 219 // output pixel value to the result of the stats call. 220 221 psVector *pixels = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of pixels 222 pixels->n = inputs->n; 223 psVector *weights = NULL; // Stack of weights 224 if (haveWeights) { 225 weights = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of weights 226 weights->n = inputs->n; 227 } 228 psVector *index = NULL; // The indices to sort the pixels 229 230 float keepFrac = 1.0 - params->fracLow - params->fracHigh; // Fraction of pixels to keep 231 psMaskType maskVal = params->maskVal; // The mask value 232 233 for (int i = output->row0; i < output->row0 + output->image->numRows; i++) { 234 for (int j = output->col0; j < output->col0 + output->image->numCols; j++) { 235 236 int numValid = 0; // Number of valid pixels in the stack 237 for (int r = 0; r < inputs->n; r++) { 238 // Check existence and bounds 239 if (mask->data.U8[r] & PM_READOUT_COMBINE_NO_IMAGE || 240 i < colLower->data.U32[r] || 241 i >= colUpper->data.U32[r] || 242 j < rowLower->data.U32[r] || 243 j >= rowUpper->data.U32[r]) { 244 continue; 245 } 246 247 pmReadout *readout = inputs->data[r]; // Input readout 248 int y = i - readout->row0; // y position on input readout 249 int x = j - readout->col0; // x position on input readout 250 if (readout->mask && readout->mask->data.U8[y][x] & maskVal) { 251 mask->data.U8[r] &= PM_READOUT_COMBINE_MASKED; 252 continue; 253 } 254 255 pixels->data.F32[r] = readout->image->data.F32[y][x]; 256 257 // Apply zero and scale 258 if (zero) { 259 pixels->data.F32[r] -= zero->data.F32[r]; 260 } 261 if (scale) { 262 pixels->data.F32[r] /= scale->data.F32[r]; 263 } 264 265 if (haveWeights) { 266 weights->data.F32[r] = readout->weight->data.F32[y][x]; 267 if (scale) { 268 weights->data.F32[r] /= scale->data.F32[r] * scale->data.F32[r]; 269 } 270 } 271 numValid++; 272 } 273 274 if (numValid == 0) { 275 output->mask->data.U8[i][j] = PM_MASK_FLAT; 276 continue; 277 } 278 279 // Apply fracLow,fracHigh if there are enough pixels 280 if (numValid * keepFrac >= params->nKeep) { 281 index = psVectorSortIndex(index, pixels); 282 int numLow = numValid * params->fracLow; // Number of low pixels to clip 283 int numHigh = numValid * params->fracHigh; // Number of high pixels to clip 284 // Low pixels 285 for (int k = 0, numMasked = 0; numMasked < numLow; k++) { 286 // Don't count the ones that are already masked 287 if (! mask->data.U8[index->data.S32[k]] & PM_READOUT_COMBINE_BAD) { 288 mask->data.U8[index->data.S32[k]] |= PM_READOUT_COMBINE_CLIPPED; 289 numMasked++; 290 } 291 } 292 // High pixels 293 for (int k = pixels->n, numMasked = 0; numMasked < numHigh; k++) { 294 // Don't count the ones that are already masked 295 if (! mask->data.U8[index->data.S32[k]] & PM_READOUT_COMBINE_BAD) { 296 mask->data.U8[index->data.S32[k]] |= PM_READOUT_COMBINE_CLIPPED; 297 numMasked++; 298 } 299 } 300 psFree(index); 301 } 302 303 // Combination 304 psVectorStats(stats, pixels, weights, mask, PM_READOUT_COMBINE_BAD | PM_READOUT_COMBINE_CLIPPED); 305 output->image->data.F32[i][j] = getStat(stats, params->combine); 306 if (haveWeights) { 307 output->weight->data.F32[i][j] = getStat(stats, combineStdev); 308 output->weight->data.F32[i][j] *= output->weight->data.F32[i][j]; // Squared for variance 309 } 310 311 // Clear the clipping 312 psBinaryOp(mask, mask, "&", psScalarAlloc(~PM_READOUT_COMBINE_CLIPPED, PS_TYPE_U8)); 313 314 } 315 } 316 317 psFree(rowLower); 318 psFree(rowUpper); 319 psFree(colLower); 320 psFree(colUpper); 321 psFree(pixels); 322 psFree(mask); 323 psFree(weights); 324 325 psFree(stats); 326 327 return true; 36 328 } 37 329 38 static void pmCombineParamsFree (pmCombineParams *params)39 {40 41 if (params == NULL)42 return;43 44 psFree (params->stats);45 return;46 }47 48 pmCombineParams *pmCombineParamsAlloc (psStatsOptions statsOptions)49 {50 51 pmCombineParams *params = psAlloc (sizeof(pmCombineParams));52 psMemSetDeallocator(params, (psFreeFunc) pmCombineParamsFree);53 54 params->stats = psStatsAlloc (statsOptions);55 params->maskVal = 0;56 params->fracHigh = 0.25;57 params->fracHigh = 0.25;58 params->nKeep = 3;59 60 return (params);61 }62 63 /******************************************************************************64 XXX: Must add support for S16 and S32 types. F32 currently supported.65 *****************************************************************************/66 psImage *pmReadoutCombine(psImage *output,67 const psArray *inputs,68 const psVector *zero,69 const psVector *scale,70 pmCombineParams *params,71 bool applyZeroScale,72 psF32 gain,73 psF32 readnoise)74 {75 PS_ASSERT_PTR_NON_NULL(inputs, NULL);76 PS_ASSERT_PTR_NON_NULL(params, NULL);77 PS_ASSERT_PTR_NON_NULL(params->stats, NULL);78 if (zero != NULL) {79 PS_ASSERT_VECTOR_TYPE(zero, PS_TYPE_F32, NULL);80 // PS_ASSERT_VECTOR_TYPE_S16_S32_F32(zero, NULL);81 }82 if (scale != NULL) {83 PS_ASSERT_VECTOR_TYPE(scale, PS_TYPE_F32, NULL);84 // PS_ASSERT_VECTOR_TYPE_S16_S32_F32(scale, NULL);85 }86 if ((zero != NULL) && (scale != NULL)) {87 PS_ASSERT_VECTOR_TYPE_EQUAL(zero, scale, NULL);88 // PS_ASSERT_VECTOR_TYPE_S16_S32_F32(scale, NULL);89 }90 91 psStats *stats = params->stats;92 psS32 maxInputCols = 0;93 psS32 maxInputRows = 0;94 psS32 minInputCols = PS_MAX_S32;95 psS32 minInputRows = PS_MAX_S32;96 pmReadout *tmpReadout = NULL;97 psS32 tmpI;98 psElemType outputType = PS_TYPE_F32;99 100 if (DetermineNumBits(stats->options) != 1) {101 psError(PS_ERR_UNKNOWN, true,102 "Multiple statistical options have been requested. Returning NULL.\n");103 return(NULL);104 }105 106 // We step through each readout in the input image list. If any readout is107 // NULL, empty, or has the wrong type, we generate an error and return108 // NULL. We determine how big of an output image is needed to combine109 // these input images. We do this by taking the110 // max(readout->col0 + readout->numCols + image->col0 + image->numCols)111 // max(readout->row0 + readout->numRows + image->row0 + image->numRows)112 //113 for (int i = 0; i < inputs->n; i++) {114 tmpReadout = inputs->data[i];115 PS_ASSERT_READOUT_NON_NULL(tmpReadout, output);116 PS_ASSERT_READOUT_NON_EMPTY(tmpReadout, output);117 PS_ASSERT_READOUT_TYPE(tmpReadout, PS_TYPE_F32, output);118 119 minInputRows = PS_MIN(minInputRows, (tmpReadout->row0 + tmpReadout->image->row0));120 tmpI = tmpReadout->row0 +121 tmpReadout->image->row0 +122 tmpReadout->image->numRows;123 maxInputRows = PS_MAX(maxInputRows, tmpI);124 125 minInputCols = PS_MIN(minInputCols, (tmpReadout->col0 + tmpReadout->image->col0));126 tmpI = tmpReadout->col0 +127 tmpReadout->image->col0 +128 tmpReadout->image->numCols;129 maxInputCols = PS_MAX(maxInputCols, tmpI);130 }131 132 // We ensure that the zero vector is of the proper size.133 if (zero != NULL) {134 PS_ASSERT_VECTOR_TYPE(zero, PS_TYPE_F32, NULL);135 if (zero->n < inputs->n) {136 psError(PS_ERR_UNKNOWN, true, "zero vector has incorrect size (%d). Returning NULL.\n", zero->n);137 return(NULL);138 } else if (zero->n > inputs->n) {139 // XXX EAM : abort on this condition? is probably an error140 psLogMsg(__func__, PS_LOG_WARN,141 "WARNING: the zero vector too many elements (%d)\n", zero->n);142 }143 }144 145 // We ensure that the scale vector is of the proper size.146 if (scale != NULL) {147 PS_ASSERT_VECTOR_TYPE(scale, PS_TYPE_F32, NULL);148 if (scale->n < inputs->n) {149 psError(PS_ERR_UNKNOWN, true, "scale vector has incorrect size (%d). Returning NULL.\n", scale->n);150 return(NULL);151 } else if (scale->n > inputs->n) {152 // XXX EAM : abort on this condition? is probably an error153 psLogMsg(__func__, PS_LOG_WARN,154 "WARNING: the scale vector has too many elements (%d)\n", scale->n);155 }156 }157 158 // At this point, the following variables have been computed:159 // maxInputRows: the largest input row value, in output image space.160 // maxInputCols: the largest input column value, in output image space.161 // minInputRows: the smallest input row value, in output image space.162 // minInputCols: the smallest input column value, in output image space.163 //164 if (output == NULL) {165 output = psImageAlloc(maxInputCols-minInputCols, maxInputRows-minInputRows, outputType);166 *(psS32 *) &(output->col0) = minInputCols;167 *(psS32 *) &(output->row0) = minInputRows;168 } else {169 170 // XXX EAM : recycle the existing output image? why would we care about the existing pixels?171 PS_ASSERT_IMAGE_TYPE(output, PS_TYPE_F32, NULL);172 if (((output->col0 + output->numCols) < maxInputCols) ||173 ((output->row0 + output->numRows) < maxInputRows)) {174 psError(PS_ERR_UNKNOWN, true,175 "Output image (%d, %d) is too small to hold combined images. Returning NULL.\n",176 output->row0 + output->numRows,177 output->col0 + output->numCols);178 return(NULL);179 }180 181 // reset output origin using logic of above182 *(psS32 *) &(output->col0) = minInputCols;183 *(psS32 *) &(output->row0) = minInputRows;184 }185 186 psVector *tmpPixels = psVectorAlloc(inputs->n, PS_TYPE_F32);187 psVector *tmpPixelsKeep = psVectorAlloc(inputs->n, PS_TYPE_F32);188 psVector *outRowLower = psVectorAlloc(inputs->n, PS_TYPE_U32);189 psVector *outRowUpper = psVectorAlloc(inputs->n, PS_TYPE_U32);190 psVector *outColLower = psVectorAlloc(inputs->n, PS_TYPE_U32);191 psVector *outColUpper = psVectorAlloc(inputs->n, PS_TYPE_U32);192 193 // For each input readout, we store the min/max pixel indices for that readout, in detector coordinates,194 // in the psVectors (outRowLower, outColLower, outRowUpper, outColUpper).195 for (int i = 0; i < inputs->n; i++) {196 tmpReadout = (pmReadout *) inputs->data[i];197 outRowLower->data.U32[i] = tmpReadout->row0 + tmpReadout->image->row0;198 outColLower->data.U32[i] = tmpReadout->col0 + tmpReadout->image->col0;199 outRowUpper->data.U32[i] = tmpReadout->row0 +200 tmpReadout->image->row0 +201 tmpReadout->image->numRows;202 outColUpper->data.U32[i] = tmpReadout->col0 +203 tmpReadout->image->col0 +204 tmpReadout->image->numCols;205 }206 207 // We loop through each pixel in the output image. We loop through each208 // input readout. We determine if that output pixel is contained in the209 // image from that readout. If so, we save it in psVector tmpPixels.210 // If not, we set a mask for that element in tmpPixels. Then, we mask off211 // pixels not between fracLow and fracHigh. Then we call the vector212 // stats routine on those pixels/mask. Then we set the output pixel value213 // to the result of the stats call.214 215 int nx, ny;216 int nKeep, nMin;217 float keepFrac = 1.0 - params->fracLow - params->fracHigh;218 float value = 0;219 psF32 *saveVector = tmpPixelsKeep->data.F32;220 221 for (int j = output->row0; j < (output->row0 + output->numRows) ; j++) {222 if (j % 10 == 0)223 fprintf (stderr, ".");224 for (int i = output->col0; i < (output->col0 + output->numCols) ; i++) {225 int nPix = 0;226 for (int r = 0; r < inputs->n; r++) {227 tmpReadout = (pmReadout *) inputs->data[r];228 229 // psTrace (__func__, 6, "[%d][%d]: [%d][%d] to [%d][%d]\n", i, j, outColLower->data.U32[r], outRowLower->data.U32[r], outColUpper->data.U32[r], outRowUpper->data.U32[r]);230 if (i < outColLower->data.U32[r])231 continue;232 if (i >= outColUpper->data.U32[r])233 continue;234 if (j < outRowLower->data.U32[r])235 continue;236 if (j >= outRowUpper->data.U32[r])237 continue;238 239 nx = i - (tmpReadout->col0 + tmpReadout->image->col0);240 ny = j - (tmpReadout->row0 + tmpReadout->image->row0);241 242 if (tmpReadout->mask != NULL) {243 if (tmpReadout->mask->data.U8[ny][nx] && params->maskVal)244 continue;245 }246 247 tmpPixels->data.F32[nPix] = tmpReadout->image->data.F32[ny][nx];248 // psTrace (__func__, 6, "readout[%d], image [%d][%d] is %f\n", r, i, j, tmpPixels->data.F32[nPix]);249 nPix ++;250 }251 tmpPixels->n = nPix;252 253 // are there enough valid pixels to apply fracLow,fracHigh?254 nKeep = nPix * keepFrac;255 if (nKeep >= params->nKeep) {256 psVectorSort (tmpPixels, tmpPixels);257 nMin = nPix * params->fracLow;258 tmpPixelsKeep->data.F32 = &tmpPixels->data.F32[nMin];259 tmpPixelsKeep->n = nKeep;260 } else {261 tmpPixelsKeep->data.F32 = tmpPixels->data.F32;262 tmpPixelsKeep->n = nPix;263 }264 265 // tmpPixelsKeep is already sorted. sample mean and median are very easy266 if (stats->options & PS_STAT_SAMPLE_MEAN) {267 value = 0;268 for (int r = 0; r < tmpPixelsKeep->n; r++) {269 value += tmpPixelsKeep->data.F32[r];270 }271 if (tmpPixelsKeep->n == 0) {272 value = 0;273 } else {274 value = value / tmpPixelsKeep->n;275 }276 }277 if (stats->options & PS_STAT_SAMPLE_MEDIAN) {278 int r = tmpPixelsKeep->n / 2;279 if (tmpPixelsKeep->n == 0) {280 value = 0;281 goto got_value;282 }283 if (tmpPixelsKeep->n % 2 == 1) {284 int r = 0.5*tmpPixelsKeep->n;285 value = tmpPixelsKeep->data.F32[r];286 goto got_value;287 }288 if (tmpPixelsKeep->n % 2 == 0) {289 value = 0.5*(tmpPixelsKeep->data.F32[r] +290 tmpPixelsKeep->data.F32[r-1]);291 goto got_value;292 }293 }294 got_value:295 output->data.F32[j-output->row0][i-output->col0] = value;296 }297 }298 tmpPixelsKeep->data.F32 = saveVector;299 300 psFree(tmpPixels);301 psFree(tmpPixelsKeep);302 psFree(outRowLower);303 psFree(outRowUpper);304 psFree(outColLower);305 psFree(outColUpper);306 307 return(output);308 }309 310 /******************************************************************************311 XXX: Must add support for S16 and S32 types. F32 currently supported.312 *****************************************************************************/313 psImage *pmReadoutCombine_OLD(psImage *output,314 const psList *inputs,315 pmCombineParams *params,316 const psVector *zero,317 const psVector *scale,318 bool applyZeroScale,319 psF32 gain,320 psF32 readnoise)321 {322 PS_ASSERT_PTR_NON_NULL(inputs, NULL);323 PS_ASSERT_PTR_NON_NULL(params, NULL);324 PS_ASSERT_PTR_NON_NULL(params->stats, NULL);325 if (zero != NULL) {326 PS_ASSERT_VECTOR_TYPE(zero, PS_TYPE_F32, NULL);327 // PS_ASSERT_VECTOR_TYPE_S16_S32_F32(zero, NULL);328 }329 if (scale != NULL) {330 PS_ASSERT_VECTOR_TYPE(scale, PS_TYPE_F32, NULL);331 // PS_ASSERT_VECTOR_TYPE_S16_S32_F32(scale, NULL);332 }333 if ((zero != NULL) && (scale != NULL)) {334 PS_ASSERT_VECTOR_TYPE_EQUAL(zero, scale, NULL);335 // PS_ASSERT_VECTOR_TYPE_S16_S32_F32(scale, NULL);336 }337 338 psStats *stats = params->stats;339 psS32 i;340 psS32 j;341 psS32 maxInputCols = 0;342 psS32 maxInputRows = 0;343 psS32 minInputCols = PS_MAX_S32;344 psS32 minInputRows = PS_MAX_S32;345 psListElem *tmpInput = NULL;346 pmReadout *tmpReadout = NULL;347 psS32 numInputs = 0;348 psS32 tmpI;349 psElemType outputType = PS_TYPE_F32;350 351 if (1 < DetermineNumBits(params->stats->options)) {352 psError(PS_ERR_UNKNOWN, true,353 "Multiple statistical options have been requested. Returning NULL.\n");354 return(NULL);355 }356 357 //358 // We step through each readout in the input image list. If any readout is359 // NULL, empty, or has the wrong type, we generate an error and return360 // NULL. We determine how big of an output image is needed to combine361 // these input images. We do this by taking the362 // max(readout->col0 + readout->numCols + image->col0 + image->numCols)363 // max(readout->row0 + readout->numRows + image->row0 + image->numRows)364 // We then compare that to365 // output->col0 + output->numCols366 // output->row0 + output->numRows367 // to determine if the output image actually stores that pixel. A similar368 // thing is done for the minimum row and column.369 //370 tmpInput = (psListElem *) inputs->head;371 while (NULL != tmpInput) {372 tmpReadout = (pmReadout *) tmpInput->data;373 PS_ASSERT_READOUT_NON_NULL(tmpReadout, output);374 PS_ASSERT_READOUT_NON_EMPTY(tmpReadout, output);375 PS_ASSERT_READOUT_TYPE(tmpReadout, PS_TYPE_F32, output);376 377 outputType = tmpReadout->image->type.type;378 379 minInputRows = PS_MIN(minInputRows,380 (tmpReadout->row0 + tmpReadout->image->row0));381 tmpI = tmpReadout->row0 +382 tmpReadout->image->row0 +383 tmpReadout->image->numRows;384 maxInputRows = PS_MAX(maxInputRows, tmpI);385 386 minInputCols = PS_MIN(minInputCols,387 (tmpReadout->col0 + tmpReadout->image->col0));388 tmpI = tmpReadout->col0 +389 tmpReadout->image->col0 +390 tmpReadout->image->numCols;391 maxInputCols = PS_MAX(maxInputCols, tmpI);392 tmpInput = tmpInput->next;393 numInputs++;394 }395 396 // We ensure that the zero vector is of the proper size.397 if (zero != NULL) {398 PS_ASSERT_VECTOR_TYPE(zero, PS_TYPE_F32, NULL);399 if (numInputs > zero->n) {400 psError(PS_ERR_UNKNOWN, true, "zero vector has incorrect size (%d). Returning NULL.\n", zero->n);401 return(NULL);402 } else if (numInputs < zero->n) {403 psLogMsg(__func__, PS_LOG_WARN,404 "WARNING: the zero vector too many elements (%d)\n", zero->n);405 }406 }407 408 // We ensure that the scale vector is of the proper size.409 if (scale != NULL) {410 PS_ASSERT_VECTOR_TYPE(scale, PS_TYPE_F32, NULL);411 if (numInputs > scale->n) {412 psError(PS_ERR_UNKNOWN, true, "scale vector has incorrect size (%d). Returning NULL.\n", scale->n);413 return(NULL);414 } else if (numInputs < scale->n) {415 psLogMsg(__func__, PS_LOG_WARN,416 "WARNING: the scale vector has too many elements (%d)\n", scale->n);417 }418 }419 420 // At this point, the following variables have been computed:421 // maxInputRows: the largest input row value, in output image space.422 // maxInputCols: the largest input column value, in output image space.423 // minInputRows: the smallest input row value, in output image space.424 // minInputCols: the smallest input column value, in output image space.425 //426 if (output == NULL) {427 output = psImageAlloc(maxInputCols-minInputCols,428 maxInputRows-minInputRows, outputType);429 *(psS32 *) &(output->col0) = minInputCols;430 *(psS32 *) &(output->row0) = minInputRows;431 } else {432 PS_ASSERT_IMAGE_TYPE(output, PS_TYPE_F32, NULL);433 if (((output->col0 + output->numCols) < maxInputCols) ||434 ((output->row0 + output->numRows) < maxInputRows)) {435 psError(PS_ERR_UNKNOWN, true,436 "Output image (%d, %d) is too small to hold combined images. Returning NULL.\n",437 output->row0 + output->numRows,438 output->col0 + output->numCols);439 return(NULL);440 }441 442 if ((output->col0 > minInputCols) || (output->row0 > minInputRows)) {443 psError(PS_ERR_UNKNOWN, true,444 "Output image offset is larger then input image offset. Returning NULL.\n");445 return(NULL);446 }447 }448 449 psVector *tmpPixels = psVectorAlloc(numInputs, PS_TYPE_F32);450 tmpPixels->n = tmpPixels->nalloc;451 psVector *tmpPixelErrors = psVectorAlloc(numInputs, PS_TYPE_F32);452 tmpPixelErrors->n = tmpPixelErrors->nalloc;453 psVector *tmpPixelMask = psVectorAlloc(numInputs, PS_TYPE_U8);454 tmpPixelMask->n = tmpPixelMask->nalloc;455 psVector *tmpPixelMaskNKeep = psVectorAlloc(numInputs, PS_TYPE_U8);456 tmpPixelMaskNKeep->n = tmpPixelMaskNKeep->nalloc;457 psVector *outRowLower = psVectorAlloc(numInputs, PS_TYPE_U32);458 outRowLower->n = outRowLower->nalloc;459 psVector *outRowUpper = psVectorAlloc(numInputs, PS_TYPE_U32);460 outRowUpper->n = outRowUpper->nalloc;461 psVector *outColLower = psVectorAlloc(numInputs, PS_TYPE_U32);462 outColLower->n = outColLower->nalloc;463 psVector *outColUpper = psVectorAlloc(numInputs, PS_TYPE_U32);464 outColUpper->n = outColUpper->nalloc;465 pmReadout **tmpReadouts = (pmReadout **) psAlloc(numInputs * sizeof(pmReadout *));466 467 // For each input readout, we create a pointer to that readout in468 // "tmpReadouts[]", and we store the min/max pixel indices for that469 // readout, in output image coordinates, in the psVectors470 // (outRowLower, outColLower, outRowUpper, outColUpper).471 i = 0;472 tmpInput = (psListElem *) inputs->head;473 while (NULL != tmpInput) {474 tmpReadouts[i] = (pmReadout *) tmpInput->data;475 outRowLower->data.U32[i] = tmpReadouts[i]->row0 + tmpReadouts[i]->image->row0;476 outColLower->data.U32[i] = tmpReadouts[i]->col0 + tmpReadouts[i]->image->col0;477 outRowUpper->data.U32[i] = tmpReadouts[i]->row0 +478 tmpReadouts[i]->image->row0 +479 tmpReadouts[i]->image->numRows;480 outColUpper->data.U32[i] = tmpReadouts[i]->col0 +481 tmpReadouts[i]->image->col0 +482 tmpReadouts[i]->image->numCols;483 tmpInput = tmpInput->next;484 i++;485 }486 487 // We loop through each pixel in the output image. We loop through each488 // input readout. We determine if that output pixel is contained in the489 // image from that readout. If so, we save it in psVector tmpPixels.490 // If not, we set a mask for that element in tmpPixels. Then, we mask off491 // pixels not between fracLow and fracHigh. Then we call the vector492 // stats routine on those pixels/mask. Then we set the output pixel value493 // to the result of the stats call.494 495 for (i = output->row0; i < (output->row0 + output->numRows) ; i++) {496 for (j = output->col0; j < (output->col0 + output->numCols) ; j++) {497 for (psS32 r = 0; r < numInputs ; r++) {498 // printf("[%d][%d]: [%d][%d] to [%d][%d]\n", i, j, outRowLower->data.U32[r], outColLower->data.U32[r], outRowUpper->data.U32[r], outColUpper->data.U32[r]);499 if ((outRowLower->data.U32[r] <= i) &&500 (outColLower->data.U32[r] <= j) &&501 (outRowUpper->data.U32[r] > i) &&502 (outColUpper->data.U32[r] > j)) {503 504 psS32 imageRow = i - (tmpReadouts[r]->row0 +505 tmpReadouts[r]->image->row0);506 psS32 imageCol = j - (tmpReadouts[r]->col0 +507 tmpReadouts[r]->image->col0);508 509 if ((NULL == tmpReadouts[r]->mask) ||510 !(params->maskVal && tmpReadouts[r]->mask->data.U8[imageRow][imageCol])) {511 tmpPixels->data.F32[r] = tmpReadouts[r]->image->data.F32[imageRow][imageCol];512 tmpPixelMask->data.U8[r] = 0;513 } else {514 tmpPixels->data.F32[r] = 0.0;515 tmpPixelMask->data.U8[r] = 1;516 }517 } else {518 tmpPixels->data.F32[r] = 0.0;519 tmpPixelMask->data.U8[r] = 1;520 }521 // printf("readout[%d], image [%d][%d] is %f\n", r, i, j, tmpPixels->data.F32[r]);522 }523 // At this point, we have scanned all input readouts for this524 // one output pixel.525 // for (psS32 r = 0; r < numInputs ; r++) printf("(0)tmpPixels->data.F32[%d] is %f\n", r, tmpPixels->data.F32[r]);526 527 // Determine how many pixels lie between fracLow and fracHigh.528 psS32 pixelCount = 0;529 for (psS32 r = 0; r < numInputs ; r++) {530 if (tmpPixelMask->data.U8[r] == 0) {531 if ((params->fracLow <= tmpPixels->data.F32[r]) &&532 (params->fracHigh >= tmpPixels->data.F32[r])) {533 pixelCount++;534 }535 }536 }537 538 // If more than params->nKeep pixels lie between the valid range,539 // then loop through the pixels, and mask away any pixels outside540 // that range.541 if (pixelCount >= params->nKeep) {542 for (psS32 r = 0; r < numInputs ; r++) {543 if (tmpPixelMask->data.U8[r] == 0) {544 if ((params->fracLow <= tmpPixels->data.F32[r]) &&545 (params->fracHigh >= tmpPixels->data.F32[r])) {546 tmpPixelMaskNKeep->data.U8[r] = 0;547 } else {548 tmpPixelMaskNKeep->data.U8[r] = 1;549 }550 }551 }552 }553 554 if ((gain > 0.0) && (readnoise >= 0.0)) {555 psF32 x;556 psF32 sigma;557 if (applyZeroScale == true) {558 for (psS32 r = 0; r < numInputs ; r++) {559 if (zero != NULL) {560 x = zero->data.F32[r];561 } else {562 x = 0.0;563 }564 if (scale != NULL) {565 x+= tmpPixels->data.F32[r] * scale->data.F32[r];566 } else {567 x+= tmpPixels->data.F32[r];568 }569 sigma = sqrtf((readnoise*readnoise) + gain * x) / gain;570 571 tmpPixelErrors->data.F32[r] = sigma;572 tmpPixels->data.F32[r]= x;573 }574 } else {575 for (psS32 r = 0; r < numInputs ; r++) {576 x= tmpPixels->data.F32[r];577 578 if (zero != NULL) {579 sigma = zero->data.F32[r];580 } else {581 sigma = 0.0;582 }583 if (scale != NULL) {584 sigma+= tmpPixels->data.F32[r] * scale->data.F32[r];585 } else {586 sigma+= tmpPixels->data.F32[r];587 }588 sigma = sqrtf((readnoise*readnoise) + (gain * sigma)) / gain;589 590 tmpPixelErrors->data.F32[r] = sigma;591 tmpPixels->data.F32[r]= x;592 }593 }594 // Calculate the specified statistic on the stack of pixels.595 // for (psS32 r = 0; r < numInputs ; r++) printf("(1)tmpPixels->data.F32[%d] is %f\n", r, tmpPixels->data.F32[r]);596 psStats *rc = psVectorStats(stats,597 tmpPixels,598 tmpPixelErrors,599 tmpPixelMaskNKeep,600 1);601 if (rc == NULL) {602 psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation. Returning NULL.\n");603 return(NULL);604 }605 } else {606 if (scale != NULL) {607 for (psS32 r = 0; r < numInputs ; r++) {608 tmpPixels->data.F32[r]*= scale->data.F32[r];609 }610 }611 612 // We add the zero vector, if non-NULL.613 if (zero != NULL) {614 for (psS32 r = 0; r < numInputs ; r++) {615 tmpPixels->data.F32[r]+= zero->data.F32[r];616 }617 }618 619 // Calculate the specified statistic on the stack of pixels.620 // for (psS32 r = 0; r < numInputs ; r++) printf("(2)tmpPixels->data.F32[%d] is %f\n", r, tmpPixels->data.F32[r]);621 psStats *rc = psVectorStats(stats,622 tmpPixels,623 NULL,624 tmpPixelMaskNKeep,625 1);626 if (rc == NULL) {627 psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation. Returning NULL.\n");628 return(NULL);629 }630 }631 632 633 // Set the pixel value in the output image to the stat value.634 double statValue;635 if (!p_psGetStatValue(stats, &statValue)) {636 psError(PS_ERR_UNKNOWN, true, "Could not determine stats value. Returning NULL.\n");637 return(NULL);638 } else {639 output->data.F32[i-output->row0][j-output->col0] = (psF32) statValue;640 }641 }642 }643 644 psFree(tmpPixels);645 psFree(tmpPixelErrors);646 psFree(tmpPixelMask);647 psFree(tmpPixelMaskNKeep);648 psFree(outRowLower);649 psFree(outRowUpper);650 psFree(outColLower);651 psFree(outColUpper);652 psFree(tmpReadouts);653 654 return(output);655 }656 657 658 /* This function measures the robust median at each of the minimum and maximum659 * coordinates and determines the difference and mean of the two values. The size660 * of the box used to make the measurement at each point is specified by the661 * configuration variable FRINGE_SQUARE_RADIUS. From the collection of662 * differences, the robust median is calculated, and returned as part of the663 * fringe statistics. For each fringe point, the values of delta and midValue are664 * also assigned and available to the user on return.665 */666 667 psStats *pmFringeStats(668 psArray *fringePoints,669 psImage *image,670 psMetadata *config)671 {672 PS_ASSERT_PTR_NON_NULL(fringePoints, NULL);673 // for (psS32 i = 0 ; i < fringePoints->n ; i++) {674 // if (!psMemCheckFringePoint((pmFringePoint *) fringePoints->data[i])) {675 // psError(PS_ERR_UNKNOWN, true, "fringePoints->data[%d] is not of type pmFringePoint.\n");676 // return(NULL);677 // }678 // }679 PS_ASSERT_IMAGE_NON_NULL(image, NULL);680 PS_ASSERT_IMAGE_NON_EMPTY(image, NULL);681 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);682 PS_ASSERT_PTR_NON_NULL(config, NULL);683 684 psBool rc;685 psF32 frSquareRadius = psMetadataLookupF32(&rc, config, "FRINGE_SQUARE_RADIUS");686 if (!rc) {687 psError(PS_ERR_UNKNOWN, true, "Could not determing the fringe radius from the metadata.\n");688 }689 690 psRegion minRegion;691 psRegion maxRegion;692 psStats *minStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);693 psStats *maxStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);694 psStats *diffStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);695 psVector *diffs = psVectorAlloc(fringePoints->n, PS_TYPE_F32);696 diffs->n = diffs->nalloc;697 698 //699 // Loop through each fringe point. Determine the robust mean around700 // the min and max for that fringe point. Save the difference between701 // those two numbers in psVector diffs.702 //703 // XXX: Ensure you have the radius correct. (add 1 to x1 and y1?)704 //705 for (psS32 i = 0 ; i < fringePoints->n ; i++) {706 pmFringePoint *fp = (pmFringePoint *) fringePoints->data[i];707 minRegion.x0 = fp->xMin - frSquareRadius;708 minRegion.x1 = fp->xMin + frSquareRadius;709 minRegion.y0 = fp->yMin - frSquareRadius;710 minRegion.y1 = fp->yMin + frSquareRadius;711 psImage *minSubImage = psImageSubset(image, minRegion);712 minStats = psImageStats(minStats, minSubImage, NULL, 0);713 714 maxRegion.x0 = fp->xMax - frSquareRadius;715 maxRegion.x1 = fp->xMax + frSquareRadius;716 maxRegion.y0 = fp->yMax - frSquareRadius;717 maxRegion.y1 = fp->yMax + frSquareRadius;718 psImage *maxSubImage = psImageSubset(image, maxRegion);719 maxStats = psImageStats(maxStats, maxSubImage, NULL, 0);720 721 if ((minStats == NULL) || (maxStats == NULL)) {722 psError(PS_ERR_UNKNOWN, true, "Could not determine robust mean on subimage.\n");723 psFree(minStats);724 psFree(maxStats);725 return(NULL);726 }727 728 fp->midValue = 0.5 * (maxStats->robustMedian + minStats->robustMedian);729 fp->delta = maxStats->robustMedian - minStats->robustMedian;730 diffs->data.F32[i] = fp->delta;731 }732 psFree(minStats);733 psFree(maxStats);734 735 diffStats = psVectorStats(diffStats, diffs, NULL, NULL, 0);736 psFree(diffs);737 if (diffStats == NULL) {738 psError(PS_ERR_UNKNOWN, true, "Could not determine robust median of the differences.\n");739 return(NULL);740 }741 return(diffStats);742 }743 744 /**745 *746 * The input array fluxLevels consists of Ni vectors, one per mosaic image.747 * Each vector consists of Nj elements, each a measurement of the input748 * flat-field image flux levels. All of these vectors must be constructed with749 * the same number of elements, or the function will return an error. If a chip750 * is missing from a particular image, that element should be set to NaN. The751 * vector chipGains supplies initial guesses for the chip gains. If the vector752 * contains the values 0.0 or NaN for any of the elements, the gain is set to the753 * mean of the valid values. If the vector length does not match the number of754 * chips, an warning is raised, all chip gain guesses will be set to 1.0, and the755 * vector length modified to match the number of chips defined by the supplied756 * fluxLevels. The sourceFlux input vector must be allocated (not NULL), but the757 * routine will set the vector length to the number of source images regardless758 * of the initial state of the vector. All vectors used by this function must be759 * of type PS_DATA_F64.760 *761 762 fluxLevels(i, j): for each flat field image i, this psArray contains a vector763 with an elemenmt for each chip j. So, fluxLevels(i, j) corresponds to the764 measured flux M_(i, j) for flat image i, chip j.765 766 chipGains[]: has j elements, one for each chip.767 768 769 They have the observed flux levels for each chip of each image. They want to770 solve for the actual flux levels and the gain of each chip.771 772 Okay, they want to solve for source fluxes and chip gains.773 774 *775 */776 bool pmFlatNormalization(777 psVector *sourceFlux,778 psVector *chipGains,779 psArray *fluxLevels)780 {781 PS_ASSERT_PTR_NON_NULL(fluxLevels, false);782 psS32 numImages = fluxLevels->n;783 psS32 numChips = ((psVector *) fluxLevels->data[0])->n;784 for (psS32 i = 0 ; i < numImages ; i++) {785 psVector *tmpVec = (psVector *) fluxLevels->data[i];786 PS_ASSERT_VECTOR_NON_NULL(tmpVec, false);787 PS_ASSERT_VECTOR_TYPE(tmpVec, PS_TYPE_F64, false);788 PS_ASSERT_VECTOR_SIZE(tmpVec, numChips, false);789 }790 791 //792 // Ensure that *localChipGains points to a vector of the same length as numImages.793 //794 PS_ASSERT_PTR_NON_NULL(chipGains, false);795 PS_ASSERT_VECTOR_TYPE(chipGains, PS_TYPE_F64, false);796 psVector *localChipGains = chipGains;797 if (numChips != chipGains->n) {798 psLogMsg(__func__, PS_LOG_WARN,799 "WARNING: the chipGains vector length does not match the number of chips.\n");800 localChipGains = psVectorAlloc(numChips, PS_TYPE_F64);801 localChipGains->n = localChipGains->nalloc;802 psBool rc = psVectorInit(localChipGains, 1.0);803 if (rc == false) {804 printf("XXX: gen error\n");805 }806 }807 808 //809 // If the chipGains vector contains the values 0.0 or NaN for any of the elements,810 // the gain is set to the mean of the valid values.811 //812 psBool meanFlag = false;813 psVector *chipGainsMask = psVectorAlloc(chipGains->n, PS_TYPE_U8);814 chipGainsMask->n = chipGainsMask->nalloc;815 for (psS32 i = 0 ; i < chipGains->n ; i++) {816 if ((fabs(chipGains->data.F64[i]) < FLT_EPSILON) ||817 (isnan(chipGains->data.F64[i]))) {818 chipGainsMask->data.U8[i] = 1;819 meanFlag = true;820 }821 }822 // Must calculate the mean.823 if (meanFlag == true) {824 psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);825 stats = psVectorStats(stats, chipGains, NULL, chipGainsMask, 1);826 if (stats == NULL) {827 printf("XXX: gen error\n");828 }829 psF64 mean;830 psBool rc = p_psGetStatValue(stats, &mean);831 if (rc == false) {832 printf("XXX: gen error\n");833 }834 // Set the gain to this mean for chips with a gain of 0.0 or NAN835 836 for (psS32 i = 0 ; i < chipGains->n ; i++) {837 if ((fabs(chipGains->data.F64[i]) < FLT_EPSILON) ||838 (isnan(chipGains->data.F64[i]))) {839 chipGains->data.F64[i] = mean;840 }841 }842 }843 844 //845 // Assert that sourceFlux is non-NULL, correct type, correct size.846 //847 PS_ASSERT_PTR_NON_NULL(sourceFlux, false);848 PS_ASSERT_VECTOR_TYPE(sourceFlux, PS_TYPE_F64, false);849 psVectorRealloc(sourceFlux, numImages);850 851 // psFree(psVector);852 if (numImages != chipGains->n) {853 psFree(localChipGains);854 }855 856 return(true);857 }
Note:
See TracChangeset
for help on using the changeset viewer.
