Changeset 4217
- Timestamp:
- Jun 13, 2005, 9:36:12 AM (21 years ago)
- Location:
- trunk/psModules/src
- Files:
-
- 2 edited
-
pmImageCombine.c (modified) (25 diffs)
-
pmImageCombine.h (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/pmImageCombine.c
r4185 r4217 1 1 /** @file pmImageCombine.c 2 2 * 3 * This file will ... 3 * This file will perform image combination of several images of the 4 * same field, produce a list of questionable pixels, then tag some 5 * of those pixels as cosmic rays. 4 6 * 5 * @author Paul Price, IfA 7 * @author Paul Price, IfA (original prototype) 6 8 * @author GLG, MHPCC 7 9 * 8 * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2005-06-09 06:16:43 $ 10 * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ 11 * @date $Date: 2005-06-13 19:36:12 $ 12 * 13 * XXX: pmRejectPixels() has a known bug with the pmImageTransform() call. 10 14 * 11 15 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 17 21 #include <math.h> 18 22 #include "pslib.h" 19 #include "psConstants.h"20 // XXX: Remove this, once psPixels.h is handled properly.21 #include "psPixels.h"22 23 23 24 // … … 27 28 #define PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH 100 28 29 #define PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH_INC 100 30 #define PS_COMBINE_IMAGE_MAX_QUESTIONABLE_PIXELS 1000 29 31 /****************************************************************************** 30 32 pmCombineImages(combine, questionablePixels, images, errors, masks, maskVal, 31 33 pixels, numIter, sigmaClip, stats) 34 35 XXX: Allocate a dummy psStats structure so that we don't destroy away its data. 32 36 *****************************************************************************/ 33 37 psImage *pmCombineImages(psImage *combine, ///< Combined image (output) … … 70 74 for (psS32 i = 0 ; i < numImages ; i++) { 71 75 psImage *tmpDataImg = (psImage *) images->data[i]; 72 psImage *tmpErrorImg = (psImage *) errors->data[i];73 psImage *tmpMaskImg = (psImage *) masks->data[i];74 75 76 PS_ASSERT_IMAGE_NON_NULL(tmpDataImg, combine); 76 PS_ASSERT_IMAGE_NON_NULL(tmpErrorImg, combine);77 PS_ASSERT_IMAGE_NON_NULL(tmpMaskImg, combine);78 79 77 PS_ASSERT_IMAGE_TYPE(tmpDataImg, PS_TYPE_F32, combine); 80 PS_ASSERT_IMAGE_TYPE(tmpErrorImg, PS_TYPE_F32, combine);81 PS_ASSERT_IMAGE_TYPE(tmpMaskImg, PS_TYPE_U8, combine);82 83 78 if ((tmpDataImg->numRows != numRows) || (tmpDataImg->numCols != numCols)) { 84 psError(PS_ERR_UNKNOWN, true, " errorimage %d has size (%d, %d); should be (%d, %d).\n",79 psError(PS_ERR_UNKNOWN, true, "image %d has size (%d, %d); should be (%d, %d).\n", 85 80 i, tmpDataImg->numRows, tmpDataImg->numCols, numRows, numCols); 86 81 } 87 82 88 PS_ASSERT_IMAGES_SIZE_EQUAL(tmpDataImg, tmpErrorImg, NULL); 89 PS_ASSERT_IMAGES_SIZE_EQUAL(tmpDataImg, tmpMaskImg, NULL); 90 } 91 PS_ASSERT_PTR_NON_NULL(pixels, combine); 83 if (errors != NULL) { 84 psImage *tmpErrorImg = (psImage *) errors->data[i]; 85 PS_ASSERT_IMAGE_NON_NULL(tmpErrorImg, combine); 86 PS_ASSERT_IMAGE_TYPE(tmpErrorImg, PS_TYPE_F32, combine); 87 PS_ASSERT_IMAGES_SIZE_EQUAL(tmpDataImg, tmpErrorImg, NULL); 88 } 89 90 if (masks != NULL) { 91 psImage *tmpMaskImg = (psImage *) masks->data[i]; 92 PS_ASSERT_IMAGE_NON_NULL(tmpMaskImg, combine); 93 PS_ASSERT_IMAGE_TYPE(tmpMaskImg, PS_TYPE_U8, combine); 94 PS_ASSERT_IMAGES_SIZE_EQUAL(tmpDataImg, tmpMaskImg, NULL); 95 } 96 } 92 97 PS_ASSERT_PTR_NON_NULL(stats, combine); 93 98 … … 112 117 psFree((*questionablePixels)->data[im]); 113 118 ((*questionablePixels)->data[im]) = (psPtr *) psPixelsAlloc(PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH); 119 ((psPixels *) ((*questionablePixels)->data[im]))->n = 0; 114 120 } 115 121 // … … 128 134 psVector *pixelMask = NULL; 129 135 if (masks != NULL) { 130 pixelMask = psVectorAlloc(numImages, PS_TYPE_ F32);136 pixelMask = psVectorAlloc(numImages, PS_TYPE_U8); 131 137 PS_VECTOR_SET_U8(pixelMask, 0); 132 138 } … … 140 146 if (pixels != NULL) { 141 147 // Only those specified pixels should be combined. 148 142 149 psStats *stdevStats = psStatsAlloc(PS_STAT_SAMPLE_STDEV); 150 143 151 for (psS32 p = 0 ; p < pixels->n ; p++) { 144 152 // Must initialize the mask to 0 for every pixel. … … 154 162 // Set the pixel data 155 163 pixelData->data.F32[im] = ((psImage *) images->data[im])->data.F32[row][col]; 156 157 164 // Set the pixel mask data, if necessary 158 165 if (masks != NULL) { 159 pixelMask->data. F32[im] = ((psImage *) masks->data[im])->data.F32[row][col];166 pixelMask->data.U8[im] = ((psImage *) masks->data[im])->data.F32[row][col]; 160 167 } 161 168 … … 171 178 for (psS32 iter = 0 ; iter < numIter ; iter++) { 172 179 // Combine all the pixels, using the specified stat. 180 173 181 stats = psVectorStats((psStats *) stats, pixelData, pixelErrors, pixelMask, maskVal); 174 182 psF64 combinedPixel; 175 183 psBool rc = p_psGetStatValue(stats, &combinedPixel); 176 184 if (rc != true) { 177 // XXX: Warning.185 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not combine pixels (%d, %d) with the specified stat.\n", row, col); 178 186 } 179 187 if (iter == 0) { … … 185 193 // the combined pixel value. 186 194 // 195 psS32 numRejects = 0; 187 196 for (psS32 im = 0 ; im < numImages ; im++) { 188 st ats = psVectorStats(stdevStats, pixelData, pixelErrors, pixelMask, maskVal);197 stdevStats = psVectorStats(stdevStats, pixelData, pixelErrors, pixelMask, maskVal); 189 198 psF64 stdev; 190 199 psBool rc = p_psGetStatValue(stdevStats, &stdev); 191 200 if (rc != true) { 192 // XXX: Warning. 201 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not compute the standard deviation of pixel (%d, %d).\n", row, col); 202 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not combine pixels (%d, %d) with the specified stat.\n", row, col); 193 203 } 194 204 195 if (fabs(pixelData->data.F32[im] - combinedPixel) > fabs(sigmaClip * stdev)) { 196 pixelMask->data.F32[im] = maskVal; 197 // 198 // XXX: These data structures indirections are getting complicated. 199 // 200 psS32 ptr = qpPtr->data.S32[im]; 201 psPixels *pixelListPtr = ((psPixels *) ((*questionablePixels)->data[im])); 202 if (ptr >= pixelListPtr->n) { 203 (*questionablePixels)->data[im] = (psPtr *) psPixelsRealloc(((psPixels *) ((*questionablePixels)->data[im])), ((((psPixels *) ((*questionablePixels)->data[im]))->n) + PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH_INC)); 205 if (!(pixelMask->data.U8[im] & maskVal)) { 206 if (fabs(pixelData->data.F32[im]) > 207 (fabs(sigmaClip * stdev) + fabs(combinedPixel))) { 208 numRejects++; 209 //printf("Rejecting pixel. Image %d. (row, col) (%d, %d)\n", im, row, col); 210 pixelMask->data.U8[im] = maskVal; 211 // 212 // XXX: These data structures indirections are getting complicated. 213 // 214 psS32 ptr = qpPtr->data.S32[im]; 215 psPixels *pixelListPtr = ((psPixels *) ((*questionablePixels)->data[im])); 216 217 if (ptr >= pixelListPtr->nalloc) { 218 (*questionablePixels)->data[im] = 219 (psPtr *) psPixelsRealloc(((psPixels *) ((*questionablePixels)->data[im])), 220 ((((psPixels *) ((*questionablePixels)->data[im]))->nalloc) + 221 PS_COMBINE_IMAGE_INITIAL_PIXEL_LIST_LENGTH_INC)); 222 // XXX: Can the realloc() fail? Must we check for NULL? 223 } 224 ((psPixels *) ((*questionablePixels)->data[im]))->data[ptr].x = col; 225 ((psPixels *) ((*questionablePixels)->data[im]))->data[ptr].y = row; 226 (qpPtr->data.S32[im])++; 227 ((psPixels *) ((*questionablePixels)->data[im]))->n = qpPtr->data.S32[im]; 204 228 } 205 // XXX: Can the realloc() fail? Must we check for NULL?206 ((psPixels *) ((*questionablePixels)->data[im]))->data[ptr].x = col;207 ((psPixels *) ((*questionablePixels)->data[im]))->data[ptr].y = row;208 (qpPtr->data.S32[im])++;209 229 } 230 } 231 232 // 233 // If the number of rejected pixels is zero, then there's no point 234 // continuing the loop. 235 // 236 if (numRejects == 0) { 237 break; 238 } 239 psS32 totalRejects = 0; 240 for (psS32 im = 0 ; im < numImages ; im++) { 241 if (pixelMask->data.U8[im] & maskVal) { 242 totalRejects++; 243 } 244 } 245 246 // 247 // XXX: Is it possible to have all pixels rejected? If so, we should 248 // exit the loop. 249 // 250 if (totalRejects == numImages) { 251 break; 210 252 } 211 253 } … … 231 273 // Set the pixel mask data, if necessary 232 274 if (masks != NULL) { 233 pixelMask->data. F32[im] = ((psImage *) masks->data[im])->data.F32[row][col];275 pixelMask->data.U8[im] = ((psImage *) masks->data[im])->data.F32[row][col]; 234 276 } 235 277 … … 246 288 combine->data.F32[row][col] = (psF32) tmpF64; 247 289 if (rc != true) { 248 // XXX: Warning.290 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not combine pixels (%d, %d) with the specified stat.\n", row, col); 249 291 } 250 292 } … … 299 341 return median / image->data.F32[y][x]; 300 342 } 343 344 345 /****************************************************************************** 346 DetermineRegion(image, myOutToIn): for a psImage and a psPlaneTransform to that 347 image, this routine determines the size of the input image which maps to that 348 image, and returns the result in a psRegion struct. 349 350 XXX: Basically, this routine is only guaranteed to work if the transform is 351 linear. 352 353 XXX: Shouldn't this functionality be part of psImageTransform()? 354 *****************************************************************************/ 355 static psRegion DetermineRegion(psImage *image, 356 psPlaneTransform *myOutToIn) 357 { 358 psRegion myRegion; 359 myRegion.x0 = PS_MAX_F32; 360 myRegion.x1 = PS_MIN_F32; 361 myRegion.y0 = PS_MAX_F32; 362 myRegion.y1 = PS_MIN_F32; 363 364 psPlane in; 365 psPlane out; 366 367 in.x = 0.0; 368 in.y = 0.0; 369 370 psPlaneTransformApply(&out, myOutToIn, &in); 371 if (out.x < myRegion.x0) { 372 myRegion.x0 = out.x; 373 } 374 if (out.x > myRegion.x1) { 375 myRegion.x1 = out.x; 376 } 377 if (out.y < myRegion.y0) { 378 myRegion.y0 = out.y; 379 } 380 if (out.y > myRegion.y1) { 381 myRegion.y1 = out.y; 382 } 383 384 in.x = (psF32) (image->numCols); 385 in.y = 0.0; 386 psPlaneTransformApply(&out, myOutToIn, &in); 387 if (out.x < myRegion.x0) { 388 myRegion.x0 = out.x; 389 } 390 if (out.x > myRegion.x1) { 391 myRegion.x1 = out.x; 392 } 393 if (out.y < myRegion.y0) { 394 myRegion.y0 = out.y; 395 } 396 if (out.y > myRegion.y1) { 397 myRegion.y1 = out.y; 398 } 399 400 in.x = (psF32) (image->numCols); 401 ; 402 in.y = 0.0; 403 psPlaneTransformApply(&out, myOutToIn, &in); 404 if (out.x < myRegion.x0) { 405 myRegion.x0 = out.x; 406 } 407 if (out.x > myRegion.x1) { 408 myRegion.x1 = out.x; 409 } 410 if (out.y < myRegion.y0) { 411 myRegion.y0 = out.y; 412 } 413 if (out.y > myRegion.y1) { 414 myRegion.y1 = out.y; 415 } 416 417 in.x = (psF32) (image->numCols); 418 in.y = (psF32) (image->numRows); 419 psPlaneTransformApply(&out, myOutToIn, &in); 420 if (out.x < myRegion.x0) { 421 myRegion.x0 = out.x; 422 } 423 if (out.x > myRegion.x1) { 424 myRegion.x1 = out.x; 425 } 426 if (out.y < myRegion.y0) { 427 myRegion.y0 = out.y; 428 } 429 if (out.y > myRegion.y1) { 430 myRegion.y1 = out.y; 431 } 432 433 return(myRegion); 434 } 435 436 /****************************************************************************** 437 XXX: Don't we have a psLib function for this? 438 *****************************************************************************/ 439 static psImage *ImageConvertF32(psImage *image) 440 { 441 psImage *imgF32 = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32); 442 443 for (psS32 i = 0 ; i < image->numRows ; i++) { 444 for (psS32 j = 0 ; j < image->numCols ; j++) { 445 imgF32->data.F32[i][j] = (psF32) image->data.U8[i][j]; 446 } 447 } 448 449 return(imgF32); 450 } 451 301 452 302 453 // … … 314 465 I think they mean syncs with PWP. 315 466 XXX: Must add mask parameter, use it in gradient calculation. 467 XXX: This function does not work. 316 468 *****************************************************************************/ 317 469 psArray *pmRejectPixels(const psArray *images, ///< Array of input images … … 323 475 ) 324 476 { 477 psLogMsg(__func__, PS_LOG_WARN, "WARNING: pmRejectPixels() has known bugs. Specifically, in the psImageTransform() call.\n"); 325 478 PS_ASSERT_PTR_NON_NULL(images, NULL); 326 479 PS_ASSERT_PTR_NON_NULL(errors, NULL); … … 332 485 PS_ASSERT_INT_EQUAL(numImages, inToOut->n, NULL); 333 486 PS_ASSERT_INT_EQUAL(numImages, outToIn->n, NULL); 334 // XXX: Loop through all images, ensure their sizes are equal ?487 // XXX: Loop through all images, ensure their sizes are equal 335 488 336 489 // … … 340 493 for (psS32 im = 0 ; im < numImages ; im++) { 341 494 rejects->data[im] = (psPtr *) psPixelsAlloc(PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH); 495 psPixels *pixels = (psPixels *) rejects->data[im]; 496 pixels->n = 0; 342 497 } 343 498 // … … 346 501 psVector *rPtr = psVectorAlloc(numImages, PS_TYPE_S32); 347 502 PS_VECTOR_SET_S32(rPtr, 0); 348 349 503 350 504 psS32 numCols = ((psImage *) images->data[0])->numCols; 351 505 psS32 numRows = ((psImage *) images->data[0])->numRows; 352 506 psRegion myRegion = psRegionSet(0, numCols-1, 0, numRows-1); 353 psU32 maskVal = 1; 507 psU32 maskVal = 1; // XXX: Is this appropriate? 354 508 355 509 psPlane *inCoords = psAlloc(sizeof(psPlane)); … … 361 515 // 362 516 psPixels *pixelList = (psPixels *) errors->data[im]; 517 363 518 psImage *currImage = (psImage *) images->data[im]; 364 519 myRegion.x0 = 0; 365 myRegion.x1 = currImage->numCols - 1;520 myRegion.x1 = currImage->numCols; 366 521 myRegion.y0 = 0; 367 myRegion.y1 = currImage->numRows - 1;522 myRegion.y1 = currImage->numRows; 368 523 psPlaneTransform *myInToOut = (psPlaneTransform *) inToOut->data[im]; 369 524 psPlaneTransform *myOutToIn = (psPlaneTransform *) outToIn->data[im]; … … 374 529 psImage *maskImage = NULL; 375 530 maskImage = psPixelsToMask(maskImage, pixelList, myRegion, maskVal); 531 psImage *maskImageF32 = ImageConvertF32(maskImage); 376 532 377 533 // 378 534 // Transform that mask image into detector coordinate space 379 535 // 380 // XXX: conform to new psImageTransform() 381 // psImage *transformedImage = psImageTransform(NULL, maskImage, NULL, 0, myOutToIn, NULL, 0); 382 psImage *transformedImage = NULL; 536 psRegion myRegion = DetermineRegion(maskImageF32, myOutToIn); 537 psImage *transformedImage = psImageTransform(NULL, NULL, maskImageF32, NULL, 538 0, myOutToIn, myRegion, NULL, 539 PS_INTERPOLATE_BILINEAR, 0); 540 // 541 // XXX: Currently, a possibly buggy psImageTransform() corrupts the data in 542 // in outToIn, and possibly other places. The following printf() demonstrates 543 // this in conjunction with the test code. 544 // 545 // psPlaneTransform *tmpOutToIn2 = (psPlaneTransform *) outToIn->data[0]; 546 // printf("HMMM3: (%d, %d) (%d %d)\n", tmpOutToIn2->x->nX, tmpOutToIn2->x->nY, tmpOutToIn2->y->nX, tmpOutToIn2->y->nY); 547 // 383 548 384 549 // … … 392 557 inCoords->y = 0.5 + (psF32) (pixelList->data[p]).y; 393 558 psPlaneTransformApply(outCoords, myInToOut, inCoords); 394 395 559 psF32 maskVal = (psF32) psImagePixelInterpolate(transformedImage, outCoords->x, outCoords->y, 396 560 NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); … … 408 572 if (im != otherImg) { 409 573 // Map the outCoords to inCoords that for otherImg space. 410 psPlaneTransformApply(inCoords, (psPlaneTransform * )outToIn->data[otherImg], outCoords); 411 574 psPlaneTransformApply(inCoords, 575 (psPlaneTransform * )outToIn->data[otherImg], 576 outCoords); 412 577 psS32 xPix = (int)(inCoords->x + 0.5); 413 578 psS32 yPix = (int)(inCoords->y + 0.5); … … 417 582 numGrads++; 418 583 } 419 420 584 } 421 585 } … … 434 598 // pixelList is large enough; if not, we realloc() 435 599 // 600 436 601 psS32 ptr = rPtr->data.S32[im]; 437 602 psPixels *pixelListPtr = (psPixels *) rejects->data[im]; 438 if (ptr >= pixelListPtr->n ) {603 if (ptr >= pixelListPtr->nalloc) { 439 604 rejects->data[im] = (psPtr *) psPixelsRealloc(((psPixels *) rejects->data[im]), 440 ((((psPixels *) rejects->data[im])->n) + PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH_INC)); 605 ((((psPixels *) rejects->data[im])->nalloc) + PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH_INC)); 606 // XXX: Can the realloc() fail? Must we check for NULL? 441 607 } 442 // XXX: Can the realloc() fail? Must we check for NULL? 608 443 609 ((psPixels *) rejects->data[im])->data[ptr].x = (pixelList->data[p]).x; 444 610 ((psPixels *) rejects->data[im])->data[ptr].y = (pixelList->data[p]).y; … … 450 616 psFree(myOutToIn); 451 617 psFree(maskImage); 618 psFree(maskImageF32); 452 619 psFree(transformedImage); 453 620 } 454 621 455 // psFree(myRegion);456 622 psFree(inCoords); 457 623 psFree(outCoords); 458 return( NULL);624 return(rejects); 459 625 } 460 626 -
trunk/psModules/src/pmImageCombine.h
r4060 r4217 1 1 /** @file pmImageCombine.h 2 2 * 3 * This file will ... 3 * This file will perform image combination of several images of the 4 * same field, produce a list of questionable pixels, then tag some 5 * of those pixels as cosmic rays. 4 6 * 7 * @author Paul Price, IfA (original prototype) 5 8 * @author GLG, MHPCC 6 9 * 7 * @version $Revision: 1. 1$ $Name: not supported by cvs2svn $8 * @date $Date: 2005-0 5-31 22:32:42 $10 * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ 11 * @date $Date: 2005-06-13 19:36:12 $ 9 12 * 10 13 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 22 25 #include<math.h> 23 26 #include "pslib.h" 24 #include "psConstants.h"25 27 26 28 psImage *pmCombineImages(psImage *combine, ///< Combined image (output)
Note:
See TracChangeset
for help on using the changeset viewer.
