Changeset 11115
- Timestamp:
- Jan 16, 2007, 1:51:51 PM (19 years ago)
- Location:
- trunk/psModules/src
- Files:
-
- 4 edited
-
imcombine/Makefile.am (modified) (1 diff)
-
imcombine/pmImageCombine.c (modified) (5 diffs)
-
imcombine/pmImageCombine.h (modified) (2 diffs)
-
psmodules.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/Makefile.am
r9610 r11115 5 5 6 6 libpsmodulesimcombine_la_SOURCES = \ 7 pmReadoutCombine.c 8 #pmImageCombine.c7 pmReadoutCombine.c \ 8 pmImageCombine.c 9 9 # pmImageSubtract.c 10 10 11 11 pkginclude_HEADERS = \ 12 pmReadoutCombine.h 13 #pmImageCombine.h12 pmReadoutCombine.h \ 13 pmImageCombine.h 14 14 # pmImageSubtract.h 15 15 -
trunk/psModules/src/imcombine/pmImageCombine.c
r10551 r11115 8 8 * @author GLG, MHPCC 9 9 * 10 * @version $Revision: 1.1 0$ $Name: not supported by cvs2svn $11 * @date $Date: 200 6-12-08 11:39:30$10 * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ 11 * @date $Date: 2007-01-16 23:51:51 $ 12 12 * 13 13 * XXX: pmRejectPixels() has a known bug with the pmImageTransform() call. … … 26 26 #include "pslib.h" 27 27 28 // 29 // The following macros define how big the initial pixel list will be, and 30 // how much it should be incremented when realloc'ed. 31 // 32 #define PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH 100 33 #define PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH_INC 100 34 #define PS_COMBINE_IMAGE_MAX_QUESTIONABLE_PIXELS 1000 35 /****************************************************************************** 36 pmCombineImages(combine, questionablePixels, images, errors, masks, maskVal, 37 pixels, numIter, sigmaClip, stats) 38 39 XXX: Allocate a dummy psStats structure so that we don't destroy away its data. 40 *****************************************************************************/ 28 #define PIXEL_LIST_BUFFER 100 // Size of the pixel list buffer 29 30 // Data structure for use as a buffer in combining pixels 31 typedef struct 32 { 33 psVector *pixels; // Pixel values 34 psVector *masks; // Pixel masks 35 psVector *errors; // Pixel errors 36 psStats *stats; // Statistics to use with combination 37 } 38 combineBuffer; 39 40 void combineBufferFree(combineBuffer *buffer) 41 { 42 psFree(buffer->pixels); 43 psFree(buffer->masks); 44 psFree(buffer->errors); 45 psFree(buffer->stats); 46 } 47 48 combineBuffer *combineBufferAlloc(long numImages // Number of images that will be combined 49 ) 50 { 51 combineBuffer *buffer = psAlloc(sizeof(combineBuffer)); 52 psMemSetDeallocator(buffer, (psFreeFunc)combineBufferFree); 53 54 buffer->pixels = psVectorAlloc(numImages, PS_TYPE_F32); 55 buffer->masks = psVectorAlloc(numImages, PS_TYPE_MASK); 56 buffer->errors = psVectorAlloc(numImages, PS_TYPE_F32); 57 buffer->stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); 58 59 return buffer; 60 } 61 62 63 static bool combinePixels(psImage *combine, // Combined image, for output 64 psArray *questionablePixels, // Array of rejection masks 65 int x, int y, // Position in the images 66 const psArray *images, // Array of input images 67 const psArray *errors, // Array of input error images 68 const psArray *masks, // Array of input masks 69 psU32 maskVal, // Mask value 70 psS32 numIter, // Number of rejection iterations 71 psF32 sigmaClip, // Number of standard deviations at which to reject 72 combineBuffer *buffer // Buffer for combination; to avoid multiple allocations 73 ) 74 { 75 assert(combine); 76 assert(x >= 0 && x < combine->numCols); 77 assert(y >= 0 && y < combine->numRows); 78 assert(images); 79 int numImages = images->n; // Number of images to combine 80 if (masks) { 81 assert(masks->n == numImages); 82 } 83 if (errors) { 84 assert(errors->n == numImages); 85 } 86 assert(numIter >= 0); 87 assert(sigmaClip > 0); 88 assert(!questionablePixels || (questionablePixels && questionablePixels->n == numImages)); 89 90 if (buffer) { 91 psMemIncrRefCounter(buffer); 92 } else { 93 buffer = combineBufferAlloc(numImages); 94 } 95 96 psVector *pixelData = buffer->pixels; // Values for the pixel of interest 97 psVector *pixelErrors = buffer->errors; // Errors for the pixel of interest 98 psVector *pixelMasks = buffer->masks; // Masks for the pixel of interest 99 psStats *stats = buffer->stats; // Statistics for combination 100 101 // 102 // Loop through each image, extract the pixel/mask/error data into psVectors. 103 // 104 if (!masks) { 105 psVectorInit(pixelMasks, 0); 106 } 107 if (!errors) { 108 pixelErrors = NULL; 109 } 110 for (int i = 0; i < numImages; i++) { 111 // Set the pixel data 112 psImage *image = images->data[i]; // Image of interest 113 pixelData->data.F32[i] = image->data.F32[y][x]; 114 // Set the pixel mask data, if necessary 115 if (masks) { 116 psImage *mask = masks->data[i]; // Mask of interest 117 pixelMasks->data.U8[i] = mask->data.U8[y][x]; 118 } // Set the pixel error data, if necessary 119 if (errors) { 120 psImage *error = errors->data[i]; // Error image of interest 121 pixelErrors->data.F32[i] = error->data.F32[y][x]; 122 } 123 } 124 125 // 126 // Iterate on the pixels, rejecting outliers 127 // 128 for (int iter = 0; iter < numIter; iter++) { 129 // Combine all the pixels, using the specified stat. 130 if (!psVectorStats(stats, pixelData, pixelErrors, pixelMasks, maskVal)) { 131 combine->data.F32[y][x] = NAN; 132 psFree(buffer); 133 return false; 134 } 135 float combinedPixel = stats->sampleMean; // Value of the combination 136 137 if (iter == 0) { 138 // Save the value produced with no rejection, since it may be useful later 139 // (if the rejection turns out to be unnecessary) 140 combine->data.F32[y][x] = combinedPixel; 141 } 142 143 // 144 // Reject all pixels that lie more that sigmaClip standard deviations from 145 // the combined pixel value. 146 // 147 int numRejects = 0; // Number of rejections 148 float stdev = stats->sampleStdev; 149 for (int i = 0; i < numImages; i++) { 150 if (!(pixelMasks->data.U8[i] & maskVal) && 151 fabs(pixelData->data.F32[i] - combinedPixel) > sigmaClip * stdev) { 152 // Reject pixel as questionable 153 numRejects++; 154 pixelMasks->data.U8[i] = maskVal; 155 if (questionablePixels) { 156 // Mark the pixel as questionable 157 psPixels *qp = questionablePixels->data[i]; // Questionable pixels for this image 158 int qpNum = qp->n; // Number of QPs in the image of interest 159 if (qpNum >= qp->nalloc) { 160 // Grow dynamically, if required 161 qp = psPixelsRealloc(qp, qp->nalloc + PIXEL_LIST_BUFFER); 162 questionablePixels->data[i] = qp; 163 } 164 qp->data[qpNum].x = x; 165 qp->data[qpNum].y = y; 166 qp->n++; 167 } 168 } 169 } 170 171 // 172 // If the number of rejected pixels is zero, then there's no point continuing the loop. 173 // 174 if (numRejects == 0) { 175 break; 176 } 177 178 // 179 // XXX: Is it possible to have all pixels rejected? If so, we should exit the loop. 180 // 181 } 182 183 psFree(buffer); 184 return true; 185 } 186 41 187 psImage *pmCombineImages( 42 188 psImage *combine, ///< Combined image (output) … … 48 194 const psPixels *pixels, ///< Pixels to combine 49 195 psS32 numIter, ///< Number of rejection iterations 50 psF32 sigmaClip, ///< Number of standard deviations at which to reject 51 psStats *stats) ///< Statistics to use in the combination 52 { 53 54 PS_ASSERT_PTR_NON_NULL(images, combine); 55 psU32 numImages = images->n; 196 psF32 sigmaClip ///< Number of standard deviations at which to reject 197 ) 198 { 56 199 psTrace("psModules.imcombine", 3, "Calling pmCombineImages(%ld)\n", images->n); 57 200 58 if (errors != NULL) { 59 if (images->n != errors->n) { 60 psError(PS_ERR_UNKNOWN, true, "images and errors args must have same length (%ld != %ld)\n", 61 images->n, errors->n); 62 return(combine); 63 } 64 } 65 if (masks != NULL) { 66 if (images->n != masks->n) { 67 psError(PS_ERR_UNKNOWN, true, "images and masks args must have same length (%ld != %ld)\n", 68 images->n, masks->n); 69 return(combine); 70 } 71 } 72 73 psImage *tmpImg = (psImage *) images->data[0]; 74 psU32 numRows = tmpImg->numRows; 75 psU32 numCols = tmpImg->numCols; 76 77 // 78 // Check that all images have the appropriate size and type. 79 // 80 for (psS32 i = 0 ; i < numImages ; i++) { 81 psImage *tmpDataImg = (psImage *) images->data[i]; 82 PS_ASSERT_IMAGE_NON_NULL(tmpDataImg, combine); 83 PS_ASSERT_IMAGE_TYPE(tmpDataImg, PS_TYPE_F32, combine); 84 if ((tmpDataImg->numRows != numRows) || (tmpDataImg->numCols != numCols)) { 85 psError(PS_ERR_UNKNOWN, true, "image %d has size (%d, %d); should be (%d, %d).\n", 86 i, tmpDataImg->numRows, tmpDataImg->numCols, numRows, numCols); 87 } 88 89 if (errors != NULL) { 90 psImage *tmpErrorImg = (psImage *) errors->data[i]; 91 PS_ASSERT_IMAGE_NON_NULL(tmpErrorImg, combine); 92 PS_ASSERT_IMAGE_TYPE(tmpErrorImg, PS_TYPE_F32, combine); 93 PS_ASSERT_IMAGES_SIZE_EQUAL(tmpDataImg, tmpErrorImg, NULL); 94 } 95 96 if (masks != NULL) { 97 psImage *tmpMaskImg = (psImage *) masks->data[i]; 98 PS_ASSERT_IMAGE_NON_NULL(tmpMaskImg, combine); 99 PS_ASSERT_IMAGE_TYPE(tmpMaskImg, PS_TYPE_U8, combine); 100 PS_ASSERT_IMAGES_SIZE_EQUAL(tmpDataImg, tmpMaskImg, NULL); 101 } 102 } 103 PS_ASSERT_PTR_NON_NULL(stats, combine); 104 psStatsOptions statistic = psStatsSingleOption(stats->options); 105 if (statistic == 0) { 106 psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics options set: %x\n", statistic); 107 return combine; 108 } 201 PS_ASSERT_ARRAY_NON_NULL(images, NULL); 202 PS_ASSERT_INT_POSITIVE(images->n, NULL); 203 long numImages = images->n; // Number of images 204 int numCols = ((psImage*)images->data[0])->numCols; // Size in x 205 int numRows = ((psImage*)images->data[0])->numRows; // Size in y 206 207 if (combine) { 208 PS_ASSERT_IMAGE_NON_NULL(combine, NULL); 209 PS_ASSERT_IMAGE_SIZE(combine, numCols, numRows, NULL); 210 } 211 if (questionablePixels && !*questionablePixels) { 212 PS_ASSERT_ARRAY_NON_NULL(*questionablePixels, NULL); 213 PS_ASSERT_ARRAY_SIZE(*questionablePixels, numImages, 0); 214 } 215 for (int i = 1; i < images->n; i++) { 216 psImage *image = images->data[i]; // Image of interest 217 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); 218 PS_ASSERT_IMAGE_SIZE(image, numCols, numRows, NULL); 219 } 220 if (errors) { 221 PS_ASSERT_ARRAYS_SIZE_EQUAL(images, errors, NULL); 222 for (int i = 0; i < images->n; i++) { 223 psImage *error = errors->data[i]; 224 PS_ASSERT_IMAGE_SIZE(error, numCols, numRows, NULL); 225 PS_ASSERT_IMAGE_TYPE(error, PS_TYPE_F32, NULL); 226 } 227 } 228 if (masks) { 229 PS_ASSERT_ARRAYS_SIZE_EQUAL(images, masks, NULL); 230 for (int i = 0; i < images->n; i++) { 231 psImage *mask = masks->data[i]; 232 PS_ASSERT_IMAGE_SIZE(mask, numCols, numRows, NULL); 233 PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, NULL); 234 } 235 } 236 PS_ASSERT_INT_POSITIVE(numIter, NULL); 237 PS_ASSERT_FLOAT_LARGER_THAN(sigmaClip, 0.0, NULL); 109 238 110 239 // Allocate and initialize the combined image, if necessary. 111 if ( combine == NULL) {240 if (!combine) { 112 241 combine = psImageAlloc(numCols, numRows, PS_TYPE_F32); 113 if (pixels != NULL) { 114 psImageInit(combine, 0.0); 242 if (pixels) { 243 // Set everything we're not working on to NAN 244 psImageInit(combine, NAN); 115 245 } 116 246 } … … 120 250 // struct for each image. 121 251 // 122 if (*questionablePixels == NULL) { 123 *questionablePixels = psArrayAlloc(numImages); 124 } else if ((*questionablePixels)->n != numImages) { 125 *questionablePixels = psArrayRealloc(*questionablePixels, numImages); 126 } 127 for (psS32 im = 0 ; im < numImages ; im++) { 128 psFree((*questionablePixels)->data[im]); 129 ((*questionablePixels)->data[im]) = 130 (psPtr *) psPixelsAlloc(PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH); 131 ((psPixels *) ((*questionablePixels)->data[im]))->n = 0; 132 } 133 // 134 // qpPtr is used to maintain a count of the questionable pixels for each image. 135 // 136 psVector *qpPtr = psVectorAlloc(numImages, PS_TYPE_S32); 137 psVectorInit(qpPtr, 0); 138 // 139 // Allocate the necessary psVectors for the call to psVectorStats(). 140 // These vectors will be used whether we are combining a list of pixels, 141 // or every pixel in the input images. 142 // 143 psVector *pixelData = psVectorAlloc(numImages, PS_TYPE_F32); 144 psVector *pixelMask = NULL; 145 if (masks != NULL) { 146 pixelMask = psVectorAlloc(numImages, PS_TYPE_U8); 147 psVectorInit(pixelMask, 0); 148 } 149 150 psVector *pixelErrors = NULL; 151 if (errors != NULL) { 152 pixelErrors = psVectorAlloc(numImages, PS_TYPE_F32); 153 psVectorInit(pixelErrors, 1.0); 154 } 155 156 if (pixels != NULL) { 252 if (questionablePixels) { 253 if (*questionablePixels == NULL) { 254 *questionablePixels = psArrayAlloc(numImages); 255 } else if ((*questionablePixels)->n != numImages) { 256 *questionablePixels = psArrayRealloc(*questionablePixels, numImages); 257 } 258 for (int i = 0; i < numImages; i++) { 259 psFree((*questionablePixels)->data[i]); 260 (*questionablePixels)->data[i] = psPixelsAlloc(PIXEL_LIST_BUFFER); 261 } 262 } 263 264 combineBuffer *buffer = combineBufferAlloc(numImages); // Buffer for combination 265 266 if (pixels) { 157 267 // Only those specified pixels should be combined. 158 268 159 psStats *stdevStats = psStatsAlloc(PS_STAT_SAMPLE_STDEV); 160 161 for (psS32 p = 0 ; p < pixels->n ; p++) { 162 // Must initialize the mask to 0 for every pixel. 163 psVectorInit(pixelMask, 0); 164 psS32 col = (pixels->data[p]).x; 165 psS32 row = (pixels->data[p]).y; 166 167 // 168 // Loop through each image, extract the pixel/mask/error data 169 // into psVectors. 170 // 171 for (psS32 im = 0 ; im < numImages ; im++) { 172 // Set the pixel data 173 pixelData->data.F32[im] = ((psImage *) images->data[im])->data.F32[row][col]; 174 // Set the pixel mask data, if necessary 175 if (masks != NULL) { 176 pixelMask->data.U8[im] = ((psImage *) masks->data[im])->data.F32[row][col]; 177 } 178 179 // Set the pixel error data, if necessary 180 if (errors != NULL) { 181 pixelErrors->data.F32[im] = ((psImage *) errors->data[im])->data.F32[row][col]; 182 } 269 for (int p = 0; p < pixels->n; p++) { 270 int x = pixels->data[p].x; // Column of interest 271 int y = pixels->data[p].y; // Row of interest 272 273 if (!combinePixels(combine, questionablePixels ? *questionablePixels : NULL, x, y, 274 images, errors, masks, maskVal, numIter, sigmaClip, buffer)) { 275 // Bad pixel --- no big deal 276 psErrorClear(); 183 277 } 184 185 // 186 // Iterate on the pixels, rejecting outliers 187 // 188 for (psS32 iter = 0 ; iter < numIter ; iter++) { 189 // Combine all the pixels, using the specified stat. 190 191 psVectorStats(stats, pixelData, pixelErrors, pixelMask, maskVal); 192 psF64 combinedPixel = psStatsGetValue(stats, statistic); 193 194 if (iter == 0) { 195 combine->data.F32[row][col] = (psF32) combinedPixel; 196 } 197 198 // 199 // Reject all pixels that lie more that sigmaClip standard deviations from 200 // the combined pixel value. 201 // 202 psS32 numRejects = 0; 203 for (psS32 im = 0 ; im < numImages ; im++) { 204 psVectorStats(stdevStats, pixelData, pixelErrors, pixelMask, maskVal); 205 psF64 stdev = stdevStats->sampleStdev; 206 207 if (!(pixelMask->data.U8[im] & maskVal)) { 208 if (fabs(pixelData->data.F32[im]) > 209 (fabs(sigmaClip * stdev) + fabs(combinedPixel))) { 210 numRejects++; 211 pixelMask->data.U8[im] = maskVal; 212 // 213 // XXX: These data structures indirections are getting complicated. 214 // 215 psS32 ptr = qpPtr->data.S32[im]; 216 psPixels *pixelListPtr = ((psPixels *) ((*questionablePixels)->data[im])); 217 218 if (ptr >= pixelListPtr->nalloc) { 219 (*questionablePixels)->data[im] = 220 (psPtr *) psPixelsRealloc(((psPixels *) ((*questionablePixels)->data[im])), 221 ((((psPixels *) ((*questionablePixels)->data[im]))->nalloc) + 222 PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH_INC)); 223 // XXX: Can the realloc() fail? Must we check for NULL? 224 } 225 ((psPixels *) ((*questionablePixels)->data[im]))->data[ptr].x = col; 226 ((psPixels *) ((*questionablePixels)->data[im]))->data[ptr].y = row; 227 (qpPtr->data.S32[im])++; 228 // XXX: this pixel ->n increment is wierd 229 ((psPixels *) ((*questionablePixels)->data[im]))->n = qpPtr->data.S32[im]; 230 } 231 } 232 } 233 234 // 235 // If the number of rejected pixels is zero, then there's no point 236 // continuing the loop. 237 // 238 if (numRejects == 0) { 239 break; 240 } 241 psS32 totalRejects = 0; 242 for (psS32 im = 0 ; im < numImages ; im++) { 243 if (pixelMask->data.U8[im] & maskVal) { 244 totalRejects++; 245 } 246 } 247 248 // 249 // XXX: Is it possible to have all pixels rejected? If so, we should 250 // exit the loop. 251 // 252 if (totalRejects == numImages) { 253 break; 254 } 255 } 256 } 257 psFree(stdevStats); 278 } 258 279 } else { 259 280 // … … 267 288 // combine image. 268 289 // 269 for (psS32 row = 0 ; row < numRows ; row++) { 270 for (psS32 col = 0 ; col < numCols ; col++) { 271 for (psS32 im = 0 ; im < numImages ; im++) { 272 // Set the pixel data 273 pixelData->data.F32[im] = ((psImage *) images->data[im])->data.F32[row][col]; 274 275 // Set the pixel mask data, if necessary 276 if (masks != NULL) { 277 pixelMask->data.U8[im] = ((psImage *) masks->data[im])->data.F32[row][col]; 278 } 279 280 // Set the pixel error data, if necessary 281 if (errors != NULL) { 282 pixelErrors->data.F32[im] = ((psImage *) errors->data[im])->data.F32[row][col]; 283 } 284 290 for (int y = 0; y < numRows; y++) { 291 for (int x = 0; x < numCols; x++) { 292 if (!combinePixels(combine, questionablePixels ? *questionablePixels : NULL, x, y, 293 images, errors, masks, maskVal, numIter, sigmaClip, buffer)) { 294 // Bad pixel --- no big deal 295 psErrorClear(); 285 296 } 286 // Combine all the pixels, using the specified stat.287 psVectorStats(stats, pixelData, pixelErrors, pixelMask, maskVal);288 combine->data.F32[row][col] = psStatsGetValue(stats, statistic);289 297 } 290 298 } 291 299 } 292 300 293 psFree(pixelData); 294 psFree(pixelMask); 295 psFree(pixelErrors); 296 psFree(qpPtr); 301 psFree(buffer); 297 302 298 303 psTrace("psModules.imcombine", 3, "Exiting pmCombineImages(%ld)\n", images->n); 299 return (combine);304 return combine; 300 305 } 301 306 -
trunk/psModules/src/imcombine/pmImageCombine.h
r9621 r11115 8 8 * @author GLG, MHPCC 9 9 * 10 * @version $Revision: 1. 3$ $Name: not supported by cvs2svn $11 * @date $Date: 200 6-10-18 02:29:15$10 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ 11 * @date $Date: 2007-01-16 23:51:51 $ 12 12 * 13 13 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 31 31 const psPixels *pixels, ///< Pixels to combine 32 32 psS32 numIter, ///< Number of rejection iterations 33 psF32 sigmaClip, ///< Number of standard deviations at which to reject 34 psStats *stats ///< Statistics to use in the combination 33 psF32 sigmaClip ///< Number of standard deviations at which to reject 35 34 ); 36 35 -
trunk/psModules/src/psmodules.h
r10827 r11115 61 61 62 62 // the following headers are from psModule:imcombine 63 //#include <pmImageCombine.h>63 #include <pmImageCombine.h> 64 64 //#include <pmImageSubtract.h> 65 65 #include <pmReadoutCombine.h>
Note:
See TracChangeset
for help on using the changeset viewer.
