Changeset 17228
- Timestamp:
- Mar 28, 2008, 5:10:17 PM (18 years ago)
- Location:
- trunk/psModules/src
- Files:
-
- 9 edited
-
camera/pmReadoutStack.c (modified) (1 diff)
-
camera/pmReadoutStack.h (modified) (2 diffs)
-
detrend/pmDark.c (modified) (2 diffs)
-
detrend/pmMaskBadPixels.c (modified) (7 diffs)
-
detrend/pmMaskBadPixels.h (modified) (3 diffs)
-
detrend/pmShutterCorrection.c (modified) (14 diffs)
-
detrend/pmShutterCorrection.h (modified) (9 diffs)
-
imcombine/pmReadoutCombine.c (modified) (10 diffs)
-
imcombine/pmReadoutCombine.h (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/camera/pmReadoutStack.c
r17009 r17228 3 3 #endif 4 4 5 #include <stdio.h> 6 #include <string.h> 5 7 #include <pslib.h> 6 8 7 9 #include "pmReadoutStack.h" 8 10 11 psImage *pmReadoutAnalysisImage(pmReadout *readout, // Readout containing image 12 const char *name, // Name of image in analysis metadata 13 int numCols, int numRows, // Expected size of image 14 psElemType type, // Expected type of image 15 double init // Initial value 16 ) 17 { 18 PS_ASSERT_PTR_NON_NULL(readout, false); 19 PS_ASSERT_STRING_NON_EMPTY(name, false); 20 21 bool mdok; // Status of MD lookup 22 psImage *image = psMetadataLookupPtr(&mdok, readout->analysis, name); 23 if (!image) { 24 image = psImageAlloc(numCols, numRows, type); 25 psMetadataAddImage(readout->analysis, PS_LIST_TAIL, name, 0, "Analysis image from " __FILE__, image); 26 psImageInit(image, init); 27 return image; 28 } 29 if (image->numCols != numCols || image->numRows != numRows) { 30 psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Analysis image %s has incorrect size (%dx%d vs %dx%d)", 31 name, image->numCols, image->numRows, numCols, numRows); 32 return NULL; 33 } 34 if (image->type.type != type) { 35 psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Analysis image %s has incorrect type (%x vs %x)", 36 name, image->type.type, type); 37 return NULL; 38 } 39 return psMemIncrRefCounter(image); 40 } 9 41 10 42 bool pmReadoutUpdateSize(pmReadout *readout, int minCols, int minRows, -
trunk/psModules/src/camera/pmReadoutStack.h
r16600 r17228 4 4 #include "pmHDU.h" 5 5 #include "pmFPA.h" 6 7 #define PM_READOUT_STACK_ANALYSIS_COUNT "STACK.COUNT" // Name for count image in analysis metadata 8 #define PM_READOUT_STACK_ANALYSIS_SIGMA "STACK.SIGMA" // Name for sigma image in analysis metadata 6 9 7 10 /// Update an output readout (for a stack) with the correct col0,row0 and the image size … … 21 24 ); 22 25 26 /// Return an image from analysis metadata, produced while stacking 27 psImage *pmReadoutAnalysisImage(pmReadout *readout, // Readout containing image 28 const char *name, // Name of image in analysis metadata 29 int numCols, int numRows, // Expected size of image 30 psElemType type, // Expected type of image 31 double init // Initial value 32 ); 23 33 24 34 #endif -
trunk/psModules/src/detrend/pmDark.c
r17058 r17228 227 227 } 228 228 229 psImage *counts = pmReadoutAnalysisImage(output->readouts->data[0], PM_READOUT_STACK_ANALYSIS_COUNT, 230 xSize, ySize, PS_TYPE_U16, 0); 231 if (!counts) { 232 return false; 233 } 234 psImage *sigma = pmReadoutAnalysisImage(output->readouts->data[0], PM_READOUT_STACK_ANALYSIS_SIGMA, 235 xSize, ySize, PS_TYPE_F32, NAN); 236 if (!sigma) { 237 psFree(counts); 238 return false; 239 } 240 229 241 // Iterate over pixels, fitting polynomial 230 242 psVector *pixels = psVectorAlloc(numInputs, PS_TYPE_F32); // Stack of pixels … … 269 281 ro->image->data.F32[yOut][xOut] = poly->coeff->data.F64[k]; 270 282 } 283 counts->data.U16[yOut][xOut] = poly->numFit; 284 sigma->data.F32[yOut][xOut] = poly->stdevFit; 271 285 } 272 286 } -
trunk/psModules/src/detrend/pmMaskBadPixels.c
r16413 r17228 81 81 82 82 83 psImage *pmMaskFlagSuspectPixels(psImage *out, const pmReadout *readout, float rej,84 83 bool pmMaskFlagSuspectPixels(pmReadout *output, const pmReadout *readout, float median, float stdev, 84 float rej, psMaskType maskVal) 85 85 { 86 PS_ASSERT_PTR_NON_NULL(readout, NULL);87 PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, NULL);88 PS_ASSERT_IMAGE_NON_NULL(readout->image, NULL);89 PS_ASSERT_IMAGE_NON_EMPTY(readout->image, NULL);90 PS_ASSERT_IMAGE_TYPE(readout->image, PS_TYPE_F32, NULL);86 PS_ASSERT_PTR_NON_NULL(readout, false); 87 PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false); 88 PS_ASSERT_IMAGE_NON_NULL(readout->image, false); 89 PS_ASSERT_IMAGE_NON_EMPTY(readout->image, false); 90 PS_ASSERT_IMAGE_TYPE(readout->image, PS_TYPE_F32, false); 91 91 if (readout->mask) { 92 PS_ASSERT_IMAGE_NON_EMPTY(readout->mask, NULL); 93 PS_ASSERT_IMAGES_SIZE_EQUAL(readout->image, readout->mask, NULL); 94 PS_ASSERT_IMAGE_TYPE(readout->mask, PS_TYPE_MASK, NULL); 95 } 96 if (out) { 97 PS_ASSERT_IMAGE_NON_EMPTY(out, NULL); 98 PS_ASSERT_IMAGE_TYPE(out, PS_TYPE_S32, NULL); 99 PS_ASSERT_IMAGES_SIZE_EQUAL(readout->image, out, NULL); 100 } 101 102 bool status; 92 PS_ASSERT_IMAGE_NON_EMPTY(readout->mask, false); 93 PS_ASSERT_IMAGES_SIZE_EQUAL(readout->image, readout->mask, false); 94 PS_ASSERT_IMAGE_TYPE(readout->mask, PS_TYPE_MASK, false); 95 } 96 PS_ASSERT_PTR_NON_NULL(output, false); 97 98 bool mdok; // Status of MD lookup 99 psImage *suspect = psMetadataLookupPtr(&mdok, output->analysis, PM_MASK_ANALYSIS_SUSPECT); // Suspect img 100 if (suspect) { 101 PS_ASSERT_IMAGE_NON_EMPTY(suspect, false); 102 PS_ASSERT_IMAGE_TYPE(suspect, PS_TYPE_S32, false); 103 PS_ASSERT_IMAGES_SIZE_EQUAL(readout->image, suspect, false); 104 psMemIncrRefCounter(suspect); 105 } else { 106 suspect = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_S32); 107 psImageInit(suspect, 0); 108 psMetadataAddImage(output->analysis, PS_LIST_TAIL, PM_MASK_ANALYSIS_SUSPECT, PS_META_REPLACE, 109 "Suspect pixels", suspect); 110 psMetadataAddS32(output->analysis, PS_LIST_TAIL, PM_MASK_ANALYSIS_NUM, PS_META_REPLACE, 111 "Number of input images", 0); 112 } 113 114 if (!isfinite(median) || !isfinite(stdev)) { 115 // If we get down here and the statistics are missing, then we should go and mask the entire image 116 psWarning("Missing statistics --- flagging entire image as suspect."); 117 return (psImage*)psBinaryOp(suspect, suspect, "+", psScalarAlloc(1.0, PS_TYPE_S32)); 118 } 119 103 120 psImage *image = readout->image; // Image of interest 104 121 psImage *mask = readout->mask; // Corresponding mask 105 122 106 if (!out) {107 out = psImageAlloc(image->numCols, image->numRows, PS_TYPE_S32);108 psImageInit(out, 0);109 }110 111 bool whole = false; // Mask whole image?112 float median = psMetadataLookupF32 (&status, readout->analysis, "READOUT.MEDIAN");113 if (!status || !isfinite(median)) {114 whole = true;115 }116 float stdev = psMetadataLookupF32 (&status, readout->analysis, "READOUT.STDEV");117 if (!status || !isfinite(stdev)) {118 whole = true;119 }120 121 122 123 psTrace ("psModules.detrend", 3, "suspect: %f +/- %f\n", median, stdev); 123 124 if (whole) {125 // If we get down here and the statistics are missing, then we should go and mask the entire image126 psWarning("Missing statistics --- flagging entire image as suspect.");127 return (psImage*)psBinaryOp(out, out, "+", psScalarAlloc(1.0, PS_TYPE_S32));128 }129 124 130 125 for (int y = 0; y < image->numRows; y++) { … … 132 127 if (fabs((image->data.F32[y][x] - median) / stdev) >= rej && 133 128 (!mask || !(mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal))) { 134 out->data.S32[y][x]++; 135 } 136 } 137 } 138 139 return out; 140 } 141 142 psImage *pmMaskIdentifyBadPixels(const psImage *suspects, psMaskType maskVal, int nTotal, float thresh, pmMaskIdentifyMode mode) 129 suspect->data.S32[y][x]++; 130 } 131 } 132 } 133 psFree(suspect); // Drop reference 134 135 psMetadataItem *numItem = psMetadataLookup(output->analysis, PM_MASK_ANALYSIS_NUM); // Item with number 136 assert(numItem); 137 numItem->data.S32++; 138 139 return true; 140 } 141 142 bool pmMaskIdentifyBadPixels(pmReadout *output, psMaskType maskVal, float thresh, pmMaskIdentifyMode mode) 143 143 { 144 PS_ASSERT_IMAGE_NON_NULL(suspects, NULL); 145 PS_ASSERT_IMAGE_NON_EMPTY(suspects, NULL); 146 PS_ASSERT_IMAGE_TYPE(suspects, PS_TYPE_S32, NULL); 144 PS_ASSERT_PTR_NON_NULL(output, false); 145 psImage *suspects = psMetadataLookupPtr(NULL, output->analysis, PM_MASK_ANALYSIS_SUSPECT); // Suspect img 146 if (!suspects) { 147 psError(PS_ERR_UNKNOWN, false, "Unable to find image with suspected bad pixels."); 148 return false; 149 } 150 PS_ASSERT_IMAGE_NON_EMPTY(suspects, false); 151 PS_ASSERT_IMAGE_TYPE(suspects, PS_TYPE_S32, false); 152 if (output->mask) { 153 PS_ASSERT_IMAGE_NON_EMPTY(output->mask, false); 154 PS_ASSERT_IMAGES_SIZE_EQUAL(output->mask, suspects, false); 155 PS_ASSERT_IMAGE_TYPE(output->mask, PS_TYPE_MASK, false); 156 } else { 157 output->mask = psImageAlloc(suspects->numCols, suspects->numRows, PS_TYPE_MASK); 158 } 159 int num = psMetadataLookupS32(NULL, output->analysis, PM_MASK_ANALYSIS_NUM); // Number of inputs 160 PS_ASSERT_INT_POSITIVE(num, false); 147 161 148 162 float limit = NAN; // Limit for masking 149 150 163 switch (mode) { 151 164 case PM_MASK_ID_VALUE: … … 154 167 155 168 case PM_MASK_ID_FRACTION: 156 limit = thresh * n Total;169 limit = thresh * num; 157 170 break; 158 171 … … 200 213 limit = max + 1.0 - thresh * sqrtf((float)max + 1.0); 201 214 202 psTrace ("psModules.detrend", 3, "bad: mode: %d, stdev: %f, limit: %f\n", max, sqrtf((float)max + 1.0), limit); 215 psTrace ("psModules.detrend", 3, "bad: mode: %d, stdev: %f, limit: %f\n", 216 max, sqrtf((float)max + 1.0), limit); 203 217 break; 204 218 } … … 207 221 return NULL; 208 222 } 209 210 psImage *badpix = psImageAlloc(suspects->numCols, suspects->numRows, PS_TYPE_MASK); // Bad pixel mask211 psImageInit(badpix, 0);212 223 213 224 if (psTraceGetLevel("psModules.detrend") > 9) { … … 227 238 psTrace ("psModules.detrend", 3, "bad pixel threshold: %f", limit); 228 239 240 psImage *badpix = output->mask; // Bad pixel mask 241 psImageInit(badpix, 0); 242 229 243 for (int y = 0; y < suspects->numRows; y++) { 230 244 for (int x = 0; x < suspects->numCols; x++) { … … 235 249 } 236 250 237 return badpix;251 return true; 238 252 } 239 253 240 254 pmMaskIdentifyMode pmMaskIdentifyModeFromString (const char *string) { 241 255 242 if (!strcasecmp 256 if (!strcasecmp(string, "VALUE")) { 243 257 return PM_MASK_ID_VALUE; 244 258 } 245 if (!strcasecmp 259 if (!strcasecmp(string, "FRACTION")) { 246 260 return PM_MASK_ID_FRACTION; 247 261 } 248 if (!strcasecmp 262 if (!strcasecmp(string, "SIGMA")) { 249 263 return PM_MASK_ID_SIGMA; 250 264 } 251 if (!strcasecmp 265 if (!strcasecmp(string, "POISSON")) { 252 266 return PM_MASK_ID_POISSON; 253 267 } -
trunk/psModules/src/detrend/pmMaskBadPixels.h
r15910 r17228 5 5 * @author Eugene Magnier, IfA 6 6 * 7 * @version $Revision: 1.1 5$ $Name: not supported by cvs2svn $8 * @date $Date: 200 7-12-24 21:10:28$7 * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2008-03-29 03:10:17 $ 9 9 * Copyright 2004 Institute for Astronomy, University of Hawaii 10 10 */ … … 15 15 /// @addtogroup detrend Detrend Creation and Application 16 16 /// @{ 17 18 #define PM_MASK_ANALYSIS_SUSPECT "MASK.SUSPECT" // Readout analysis metadata keyword for suspect image 19 #define PM_MASK_ANALYSIS_NUM "MASK.NUM" // Readout analysis metadata keyword for number of inputs 20 17 21 18 22 typedef enum { … … 47 51 /// image is of type S32. The relevant median and standard deviation must be supplied in the 48 52 /// readout->analysis metadata as READOUT.MEDIAN, READOUT.STDEVe 49 psImage *pmMaskFlagSuspectPixels(psImage *out, ///< Suspected bad pixels image, or NULL 50 const pmReadout *readout, ///< Readout to inspect 51 float rej, ///< Rejection threshold (standard deviations) 52 psMaskType maskVal ///< Mask value for statistics 53 ); 53 bool pmMaskFlagSuspectPixels(pmReadout *output, ///< Output readout, optionally with suspect pixels image 54 const pmReadout *readout, ///< Readout to inspect 55 float median, ///< Image median 56 float stdev, ///< Image standard deviation 57 float rej, ///< Rejection threshold (standard deviations) 58 psMaskType maskVal ///< Mask value for statistics 59 ); 54 60 55 61 /// Identify bad pixels from the suspect pixels image 56 62 /// 57 /// Bad pixels are identified from the suspect pixels image (accumulated over a large number of images). 58 /// Pixels marked as suspect in more than "thresh" standard deviations from the mean are identified as bad 59 /// pixels (output image). If "thresh" is negative, a Poisson is assumed. 60 psImage *pmMaskIdentifyBadPixels(const psImage *suspects, ///< Accumulated suspect pixels image 61 psMaskType maskVal, ///< Value to set for bad pixels 62 int nTotal, 63 float thresh, ///< Threshold for bad pixel (standard deviations) 64 pmMaskIdentifyMode mode 65 ); 63 /// Bad pixels are identified from the suspect pixels image (accumulated over a large number of images), 64 /// according to the chosen mode. 65 bool pmMaskIdentifyBadPixels(pmReadout *output, ///< Output readout, with suspect pixels imageOut 66 psMaskType maskVal, ///< Value to set for bad pixels 67 float thresh, ///< Threshold for bad pixel 68 pmMaskIdentifyMode mode ///< Mode for identifying bad pixels 69 ); 66 70 /// @} 67 71 #endif -
trunk/psModules/src/detrend/pmShutterCorrection.c
r16600 r17228 78 78 corr->offset = 0.0; 79 79 corr->offref = 0.0; 80 corr->num = 0; 81 corr->stdev = NAN; 80 82 81 83 return corr; … … 190 192 191 193 // linear fit to the counts and exptime, given a value for offref 192 pmShutterCorrection *pmShutterCorrectionLinFit(const psVector *exptime, 193 const psVector *counts, 194 const psVector *cntError, 195 const psVector *mask, 196 float offref, 197 int nIter, 198 float rej, 199 psMaskType maskVal 200 ) 194 pmShutterCorrection *pmShutterCorrectionLinFit(const psVector *exptime, const psVector *counts, 195 const psVector *cntError, const psVector *mask, float offref, 196 int nIter, float rej, psMaskType maskVal) 201 197 { 202 198 PS_ASSERT_VECTOR_NON_NULL(exptime, NULL); … … 249 245 return NULL; 250 246 } 251 psFree(stats);252 247 253 248 pmShutterCorrection *corr = pmShutterCorrectionAlloc(); … … 255 250 corr->scale = line->coeff[1][0]; 256 251 corr->offset = line->coeff[0][1] / line->coeff[1][0]; 257 252 corr->num = stats->clippedNvalues; 253 corr->stdev = stats->clippedStdev; 254 255 psFree(stats); 258 256 psFree(x); 259 257 psFree(y); … … 263 261 } 264 262 265 static psF32 pmShutterCorrectionModel(psVector *deriv, 266 const psVector *params, 267 const psVector *x) 263 static psF32 pmShutterCorrectionModel(psVector *deriv, const psVector *params, const psVector *x) 268 264 { 269 265 // This is in a tight loop, so we won't assert here. … … 283 279 284 280 // non-linear fit to the counts and exptime, given a guess for the three parameters 285 pmShutterCorrection *pmShutterCorrectionFullFit(const psVector *exptime, 286 const psVector *counts, 287 const psVector *cntError, 288 const pmShutterCorrection *guess 289 ) 281 pmShutterCorrection *pmShutterCorrectionFullFit(const psVector *exptime, const psVector *counts, 282 const psVector *cntError, const pmShutterCorrection *guess) 290 283 { 291 284 PS_ASSERT_VECTOR_NON_NULL(exptime, NULL); … … 638 631 639 632 640 bool pmShutterCorrectionApply(pmReadout *readout, const pmReadout *shutter )633 bool pmShutterCorrectionApply(pmReadout *readout, const pmReadout *shutter, psMaskType blank) 641 634 { 642 635 PS_ASSERT_PTR_NON_NULL(readout, false); … … 682 675 } 683 676 psImage *image = readout->image; // Image to correct 677 psImage *mask = readout->mask; // Corresponding mask 684 678 685 679 if (exptime <= 0.0) { … … 688 682 for (int y = 0; y < image->numRows; y++) { 689 683 for (int x = 0; x < image->numCols; x++) { 684 if (mask && !isfinite(shutterImage->data.F32[y][x])) { 685 mask->data.PS_TYPE_MASK_DATA[y][x] |= blank; 686 image->data.F32[y][x] = NAN; 687 continue; 688 } 690 689 image->data.F32[y][x] *= 1.0 / (exptime + shutterImage->data.F32[y][x]); 691 690 } … … 699 698 for (int y = 0; y < image->numRows; y++) { 700 699 for (int x = 0; x < image->numCols; x++) { 700 if (mask && !isfinite(shutterImage->data.F32[y][x])) { 701 mask->data.PS_TYPE_MASK_DATA[y][x] |= blank; 702 image->data.F32[y][x] = NAN; 703 continue; 704 } 701 705 image->data.F32[y][x] *= exptime / (exptime + shutterImage->data.F32[y][x]); 702 706 } … … 903 907 PS_ASSERT_ARRAY_NON_NULL(inputs, false); 904 908 PS_ASSERT_INT_EQUAL(data->num, inputs->n, false); 909 PS_ASSERT_INT_NONNEGATIVE(nIter, false); 910 PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false); 905 911 906 912 int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine … … 915 921 if (pattern) { 916 922 pmReadoutUpdateSize(pattern, minInputCols, minInputRows, xSize, ySize, false, false, maskVal); 923 } 924 925 psImage *nums = pmReadoutAnalysisImage(shutter, PM_READOUT_STACK_ANALYSIS_COUNT, xSize, ySize, 926 PS_TYPE_U16, 0); // Image with number fitted per pixel 927 if (!nums) { 928 return false; 929 } 930 psImage *sigma = pmReadoutAnalysisImage(shutter, PM_READOUT_STACK_ANALYSIS_SIGMA, xSize, ySize, 931 PS_TYPE_F32, NAN); // Image with stdev per pixel 932 if (!sigma) { 933 psFree(nums); 934 return false; 917 935 } 918 936 … … 955 973 pmShutterCorrection *corr = pmShutterCorrectionLinFit(data->exptimes, counts, errors, mask, 956 974 reference, nIter, rej, maskVal); 975 if (!corr) { 976 // Nothing we can do about it 977 psErrorClear(); 978 shutterImage->data.F32[yOut][xOut] = NAN; 979 patternImage->data.F32[yOut][xOut] = NAN; 980 nums->data.U16[yOut][xOut] = 0; 981 sigma->data.F32[yOut][xOut] = NAN; 982 continue; 983 } 957 984 shutterImage->data.F32[yOut][xOut] = corr->offset; 958 985 patternImage->data.F32[yOut][xOut] = corr->scale; 986 nums->data.U16[yOut][xOut] = corr->num; 987 sigma->data.F32[yOut][xOut] = corr->stdev; 959 988 psFree(corr); 960 989 } … … 963 992 psFree(errors); 964 993 psFree(counts); 994 psFree(nums); 995 psFree(sigma); 965 996 966 997 // Update the "concepts" -
trunk/psModules/src/detrend/pmShutterCorrection.h
r13870 r17228 5 5 * @author Paul Price, IfA 6 6 * 7 * @version $Revision: 1.1 4$ $Name: not supported by cvs2svn $8 * @date $Date: 200 7-06-19 03:40:48$7 * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2008-03-29 03:10:17 $ 9 9 * Copyright 2006 Institute for Astronomy, University of Hawaii 10 10 */ … … 61 61 62 62 /// Shutter correction parameters, applicable for a single pixel 63 typedef struct 64 { 63 typedef struct { 65 64 double scale; ///< The normalisation for an exposure, A(k) 66 65 double offset; ///< The time offset, dTk 67 66 double offref; ///< The reference time offset, dTo 68 } 69 pmShutterCorrection; 67 int num; ///< Number of points used 68 float stdev; ///< Standard deviation 69 } pmShutterCorrection; 70 70 71 71 /// Allocator for shutter correction parameters … … 76 76 /// This function is used before doing the full non-linear fit, to get parameters close to the true. Assumes 77 77 /// exptime vector is sorted (ascending order; longest is last) prior to input. 78 pmShutterCorrection *pmShutterCorrectionGuess(const psVector *exptime, ///< Exposure times for each exposure 79 const psVector *counts ///< Counts for each exposure 80 ); 78 pmShutterCorrection *pmShutterCorrectionGuess( 79 const psVector *exptime, ///< Exposure times for each exposure 80 const psVector *counts ///< Counts for each exposure 81 ); 81 82 82 83 /// Generate shutter correction based on a linear fit … … 84 85 /// Performs a linear fit to counts as a function of exposure time, with the reference time offset fixed (so 85 86 /// that the system is linear). Performs iterative clipping, if nIter > 1. 86 pmShutterCorrection *pmShutterCorrectionLinFit(const psVector *exptime, ///< Exposure times for each exposure 87 const psVector *counts, ///< Counts for each exposure 88 const psVector *cntError, ///< Error in the counts 89 const psVector *mask, ///< Mask for each exposure 90 float offref, ///< Reference time offset 91 int nIter, ///< Number of iterations 92 float rej, ///< Rejection threshold (sigma) 93 psMaskType maskVal ///< Mask value 94 ); 87 pmShutterCorrection *pmShutterCorrectionLinFit( 88 const psVector *exptime, ///< Exposure times for each exposure 89 const psVector *counts, ///< Counts for each exposure 90 const psVector *cntError, ///< Error in the counts 91 const psVector *mask, ///< Mask for each exposure 92 float offref, ///< Reference time offset 93 int nIter, ///< Number of iterations 94 float rej, ///< Rejection threshold (sigma) 95 psMaskType maskVal ///< Mask value 96 ); 95 97 96 98 /// Generate shutter correction based on a full non-linear fit … … 99 101 /// the reference time offset, so that future fits may be performed using linear fitting with the reference 100 102 /// time offset fixed. 101 pmShutterCorrection *pmShutterCorrectionFullFit(const psVector *exptime, ///< Exposure times for each exposure 102 const psVector *counts, ///< Counts for each exposure 103 const psVector *cntError, ///< Error in the counts 104 const pmShutterCorrection *guess ///< Initial guess 105 ); 103 pmShutterCorrection *pmShutterCorrectionFullFit( 104 const psVector *exptime, ///< Exposure times for each exposure 105 const psVector *counts, ///< Counts for each exposure 106 const psVector *cntError, ///< Error in the counts 107 const pmShutterCorrection *guess ///< Initial guess 108 ); 106 109 107 110 /// Measure a shutter correction image from an array of images … … 111 114 /// measuring the reference time offset using the full non-linear fit for a small number of representative 112 115 /// regions (middle and corners), and then using that to perform a linear fit to each pixel. 113 bool pmShutterCorrectionMeasure(pmReadout *output, ///< Output readout 114 const psArray *readouts, ///< Array of readouts 115 int size, ///< Size of samples for statistics for non-linear fit 116 psStatsOptions meanStat, ///< Statistic to use for mean 117 psStatsOptions stdevStat, ///< Statistic to use for stdev 118 int nIter, ///< Number of iterations 119 float rej, ///< Rejection threshold (sigma) 120 psMaskType maskVal ///< Mask value 116 bool pmShutterCorrectionMeasure( 117 pmReadout *output, ///< Output readout 118 const psArray *readouts, ///< Array of readouts 119 int size, ///< Size of samples for statistics for non-linear fit 120 psStatsOptions meanStat, ///< Statistic to use for mean 121 psStatsOptions stdevStat, ///< Statistic to use for stdev 122 int nIter, ///< Number of iterations 123 float rej, ///< Rejection threshold (sigma) 124 psMaskType maskVal ///< Mask value 121 125 ); 122 126 … … 124 128 /// 125 129 /// Given a shutter correction (with dT for each pixel), applies this correction to an input image. 126 bool pmShutterCorrectionApply(pmReadout *readout, ///< Readout to which to apply shutter correction 127 const pmReadout *shutter ///< Shutter correction readout, with dT for each pixel 128 ); 130 bool pmShutterCorrectionApply( 131 pmReadout *readout, ///< Readout to which to apply shutter correction 132 const pmReadout *shutter, ///< Shutter correction readout, with dT for each pixel 133 psMaskType blank ///< Value to give blank pixels 134 ); 129 135 130 136 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// … … 156 162 /// 157 163 /// Performs statistics on the readout, recording the data 158 bool pmShutterCorrectionAddReadout(pmShutterCorrectionData *data, ///< Correction data 159 const pmReadout *readout, ///< Readout to add 160 psStatsOptions meanStat, ///< Statistic to use for mean 161 psStatsOptions stdevStat, ///< Statistic to use for stdev 162 psMaskType maskVal, ///< Mask value 163 psRandom *rng ///< Random number generator 164 bool pmShutterCorrectionAddReadout( 165 pmShutterCorrectionData *data, ///< Correction data 166 const pmReadout *readout, ///< Readout to add 167 psStatsOptions meanStat, ///< Statistic to use for mean 168 psStatsOptions stdevStat, ///< Statistic to use for stdev 169 psMaskType maskVal, ///< Mask value 170 psRandom *rng ///< Random number generator 164 171 ); 165 172 166 173 /// Calculate the reference shutter time from the correction data 167 float pmShutterCorrectionReference(const pmShutterCorrectionData *data ///< Correction data 174 float pmShutterCorrectionReference( 175 const pmShutterCorrectionData *data ///< Correction data 168 176 ); 169 177 … … 171 179 /// 172 180 /// Performs the linear fit to each pixel in the stack. 173 bool pmShutterCorrectionGenerate(pmReadout *shutter, ///< Shutter correction 174 pmReadout *pattern, ///< Background pattern (or NULL) 175 const psArray *inputs, ///< Stack of input pmReadouts 176 float reference, ///< Reference shutter time (from pmShutterCorrectionRef) 177 const pmShutterCorrectionData *data, ///< Correction data 178 int nIter, ///< Number of iterations 179 float rej, ///< Rejection threshold (sigma) 180 psMaskType maskVal ///< Mask value 181 bool pmShutterCorrectionGenerate( 182 pmReadout *shutter, ///< Shutter correction 183 pmReadout *pattern, ///< Background pattern (or NULL) 184 const psArray *inputs, ///< Stack of input pmReadouts 185 float reference, ///< Reference shutter time (from pmShutterCorrectionRef) 186 const pmShutterCorrectionData *data, ///< Correction data 187 int nIter, ///< Number of iterations 188 float rej, ///< Rejection threshold (sigma) 189 psMaskType maskVal ///< Mask value 181 190 ); 182 191 -
trunk/psModules/src/imcombine/pmReadoutCombine.c
r17137 r17228 18 18 19 19 //#define SHOW_BUSY 1 // Show that the function is busy 20 20 21 21 22 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// … … 94 95 } 95 96 96 psStats *stats = psStatsAlloc(params->combine); // The statistics to use in the combination 97 psStatsOptions combineStdev = 0; // Statistics option for weights 98 switch (params->combine) { 99 case PS_STAT_SAMPLE_MEAN: 100 case PS_STAT_SAMPLE_MEDIAN: 101 combineStdev = PS_STAT_SAMPLE_STDEV; 102 break; 103 case PS_STAT_ROBUST_MEDIAN: 104 combineStdev = PS_STAT_ROBUST_STDEV; 105 break; 106 case PS_STAT_FITTED_MEAN: 107 combineStdev = PS_STAT_FITTED_STDEV; 108 break; 109 case PS_STAT_CLIPPED_MEAN: 110 combineStdev = PS_STAT_CLIPPED_STDEV; 111 break; 112 default: 113 psAbort("Should never get here --- checked params->combine before.\n"); 114 } 115 116 psStats *stats = psStatsAlloc(params->combine | combineStdev); // The statistics to use in the combination 97 117 if (params->combine == PS_STAT_CLIPPED_MEAN) { 98 118 stats->clipSigma = params->rej; … … 107 127 } 108 128 } 129 if (params->weights && first) { 130 psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, 131 "Using input weights to combine images", ""); 132 } 109 133 110 134 int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine … … 120 144 psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0); 121 145 122 psStatsOptions combineStdev = 0; // Statistics option for weights 123 if (params->weights) { 124 if (first) { 125 psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, 126 "Using input weights to combine images", ""); 127 } 128 129 // Get the correct statistics option for weights 130 switch (params->combine) { 131 case PS_STAT_SAMPLE_MEAN: 132 case PS_STAT_SAMPLE_MEDIAN: 133 combineStdev = PS_STAT_SAMPLE_STDEV; 134 break; 135 case PS_STAT_ROBUST_MEDIAN: 136 combineStdev = PS_STAT_ROBUST_STDEV; 137 break; 138 case PS_STAT_FITTED_MEAN: 139 combineStdev = PS_STAT_FITTED_STDEV; 140 break; 141 case PS_STAT_CLIPPED_MEAN: 142 combineStdev = PS_STAT_CLIPPED_STDEV; 143 break; 144 default: 145 psAbort("Should never get here --- checked params->combine before.\n"); 146 } 147 stats->options |= combineStdev; 148 } 146 psImage *counts = pmReadoutAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_COUNT, xSize, ySize, 147 PS_TYPE_U16, 0); 148 if (!counts) { 149 return false; 150 } 151 psImage *sigma = pmReadoutAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_SIGMA, xSize, ySize, 152 PS_TYPE_F32, NAN); 153 if (!sigma) { 154 psFree(counts); 155 return false; 156 } 157 158 stats->options |= combineStdev; 149 159 150 160 // We loop through each pixel in the output image. We loop through each input readout. We determine if … … 273 283 outputMask[yOut][xOut] = params->blank; 274 284 outputImage[yOut][xOut] = NAN; 285 counts->data.U16[yOut][xOut] = 0; 286 sigma->data.F32[yOut][xOut] = NAN; 275 287 continue; 276 288 } … … 288 300 maskData[indexData[k]] = 1; 289 301 numMasked++; 302 numValid--; 290 303 } 291 304 } … … 293 306 for (int k = pixels->n - 1, numMasked = 0; numMasked < numHigh && k >= 0; k--) { 294 307 // Don't count the ones that are already masked 295 if (! 308 if (!maskData[indexData[k]]) { 296 309 maskData[indexData[k]] = 1; 297 310 numMasked++; 311 numValid--; 298 312 } 299 313 } 300 314 } 315 counts->data.U16[yOut][xOut] = numValid; 301 316 302 317 // XXXXX this step probably is very expensive : convert errors to variance everywhere? … … 314 329 outputWeight[yOut][xOut] = NAN; 315 330 } 331 sigma->data.F32[yOut][xOut] = NAN; 316 332 } else { 317 333 outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine); … … 323 339 // also, the weighted mean is not obviously the correct thing here 324 340 } 341 sigma->data.F32[yOut][xOut] = psStatsGetValue(stats, combineStdev); 325 342 } 326 343 } … … 338 355 psFree(stats); 339 356 psFree(invScale); 357 psFree(counts); 358 psFree(sigma); 340 359 341 360 // Update the "concepts" -
trunk/psModules/src/imcombine/pmReadoutCombine.h
r13591 r17228 5 5 * @author Paul Price, IfA 6 6 * 7 * @version $Revision: 1.1 2$ $Name: not supported by cvs2svn $8 * @date $Date: 200 7-06-02 03:51:03$7 * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2008-03-29 03:10:17 $ 9 9 * Copyright 2004-2006 Institute for Astronomy, University of Hawaii 10 10 */ … … 20 20 /// These values define how the combination is performed, and should not vary by detector, so that it can be 21 21 /// re-used for multiple combinations. 22 typedef struct 23 { 22 typedef struct { 24 23 psStatsOptions combine; ///< Statistic to use when performing the combination 25 24 psMaskType maskVal; ///< Mask value … … 31 30 float rej; ///< Rejection threshould for clipping (for CLIPPED_MEAN only) 32 31 bool weights; ///< Use the supplied weights (instead of calculated stdev)? 33 } 34 pmCombineParams; 32 } pmCombineParams; 35 33 36 34 // Allocator for pmCombineParams
Note:
See TracChangeset
for help on using the changeset viewer.
