Changeset 14801
- Timestamp:
- Sep 10, 2007, 10:10:05 AM (19 years ago)
- Location:
- trunk/psModules/src/imcombine
- Files:
-
- 7 edited
-
pmSubtraction.c (modified) (9 diffs)
-
pmSubtraction.h (modified) (3 diffs)
-
pmSubtractionMatch.c (modified) (3 diffs)
-
pmSubtractionParams.c (modified) (14 diffs)
-
pmSubtractionParams.h (modified) (2 diffs)
-
pmSubtractionStamps.c (modified) (13 diffs)
-
pmSubtractionStamps.h (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/pmSubtraction.c
r14753 r14801 4 4 * @author GLG, MHPCC 5 5 * 6 * @version $Revision: 1.5 8$ $Name: not supported by cvs2svn $7 * @date $Date: 2007-09- 05 21:20:20$6 * @version $Revision: 1.59 $ $Name: not supported by cvs2svn $ 7 * @date $Date: 2007-09-10 20:10:05 $ 8 8 * 9 9 * Copyright 2004-2007 Institute for Astronomy, University of Hawaii … … 582 582 583 583 584 bool pmSubtractionCalculateEquation(psArray *stamps, const pmSubtractionKernels *kernels, int footprint) 585 { 586 PS_ASSERT_ARRAY_NON_NULL(stamps, false); 584 bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, 585 int footprint) 586 { 587 PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); 587 588 PS_ASSERT_PTR_NON_NULL(kernels, false); 588 589 PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false); … … 601 602 // We iterate over each stamp, allocate the matrix and vectors if 602 603 // necessary, and then calculate those matrix/vectors. 603 for (int i = 0; i < stamps->n ; i++) {604 pmSubtractionStamp *stamp = stamps-> data[i]; // Stamp of interest604 for (int i = 0; i < stamps->num; i++) { 605 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest 605 606 if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) { 606 607 continue; … … 736 737 737 738 738 psVector *pmSubtractionSolveEquation(psVector *solution, const p sArray*stamps)739 { 740 P S_ASSERT_ARRAY_NON_NULL(stamps, NULL);739 psVector *pmSubtractionSolveEquation(psVector *solution, const pmSubtractionStampList *stamps) 740 { 741 PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL); 741 742 742 743 // Check inputs, while summing the stamp matrices and vectors 743 744 long numParams = -1; // Number of parameters 744 for (int i = 0; i < stamps->n ; i++) {745 pmSubtractionStamp *stamp = stamps-> data[i]; // Stamp of interest745 for (int i = 0; i < stamps->num; i++) { 746 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest 746 747 PS_ASSERT_PTR_NON_NULL(stamp, NULL); 747 748 if (stamp->status != PM_SUBTRACTION_STAMP_USED) { … … 772 773 psVectorInit(sumVector, 0.0); 773 774 psImageInit(sumMatrix, 0.0); 774 for (int i = 0; i < stamps->n ; i++) {775 pmSubtractionStamp *stamp = stamps-> data[i]; // Stamp of interest775 for (int i = 0; i < stamps->num; i++) { 776 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest 776 777 if (stamp->status == PM_SUBTRACTION_STAMP_USED) { 777 778 (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); … … 810 811 811 812 812 int pmSubtractionRejectStamps(p sArray*stamps, psImage *subMask, const psVector *solution,813 int pmSubtractionRejectStamps(pmSubtractionStampList *stamps, psImage *subMask, const psVector *solution, 813 814 int footprint, float sigmaRej, const pmSubtractionKernels *kernels) 814 815 { 815 P S_ASSERT_ARRAY_NON_NULL(stamps, -1);816 PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, -1); 816 817 PS_ASSERT_IMAGE_NON_EMPTY(subMask, -1); 817 818 PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, -1); … … 825 826 826 827 // Measure deviations 827 psVector *deviations = psVectorAlloc(stamps->n , PS_TYPE_F32); // Mean deviation for stamps828 psVector *deviations = psVectorAlloc(stamps->num, PS_TYPE_F32); // Mean deviation for stamps 828 829 double totalSquareDev = 0.0; // Total square deviation from zero 829 830 int numStamps = 0; // Number of used stamps … … 834 835 float background = solution->data.F64[solution->n-1]; // The difference in background 835 836 836 for (int i = 0; i < stamps->n ; i++) {837 pmSubtractionStamp *stamp = stamps-> data[i]; // The stamp of interest837 for (int i = 0; i < stamps->num; i++) { 838 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest 838 839 if (stamp->status != PM_SUBTRACTION_STAMP_USED) { 839 840 continue; … … 899 900 900 901 int numRejected = 0; 901 for (int i = 0; i < stamps->n ; i++) {902 pmSubtractionStamp *stamp = stamps-> data[i]; // Stamp of interest902 for (int i = 0; i < stamps->num; i++) { 903 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest 903 904 if (stamp->status == PM_SUBTRACTION_STAMP_USED && deviations->data.F32[i] > limit) { 904 905 // Mask out the stamp in the image so you it's not found again -
trunk/psModules/src/imcombine/pmSubtraction.h
r14739 r14801 6 6 * @author GLG, MHPCC 7 7 * 8 * @version $Revision: 1.1 6$ $Name: not supported by cvs2svn $9 * @date $Date: 2007-09- 05 00:15:28$8 * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2007-09-10 20:10:05 $ 10 10 * 11 11 * Copyright 2004-207 Institute for Astronomy, University of Hawaii … … 47 47 48 48 /// Calculate the least-squares equation to match the image quality 49 bool pmSubtractionCalculateEquation(p sArray *stamps, ///< The stamps for which to calculate the equation,49 bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps 50 50 const pmSubtractionKernels *kernels, ///< Kernel parameters 51 51 int footprint ///< Half-size of region over which to calculate equation … … 54 54 /// Solve the least-squares equation to match the image quality 55 55 psVector *pmSubtractionSolveEquation(psVector *solution, ///< Solution vector, or NULL 56 const p sArray *stamps ///< Array of stamps56 const pmSubtractionStampList *stamps ///< Stamps 57 57 ); 58 58 59 59 /// Reject stamps 60 int pmSubtractionRejectStamps(p sArray *stamps, ///< Array of stamps to check for rejection60 int pmSubtractionRejectStamps(pmSubtractionStampList *stamps, ///< Stamps 61 61 psImage *subMask, ///< Subtraction mask 62 62 const psVector *solution, ///< Solution vector -
trunk/psModules/src/imcombine/pmSubtractionMatch.c
r14766 r14801 46 46 47 47 48 static bool getStamps(p sArray**stamps, // Stamps to read48 static bool getStamps(pmSubtractionStampList **stamps, // Stamps to read 49 49 const psArray *stampsData, // Stamp data from a file 50 50 const pmReadout *reference, // Reference readout … … 84 84 85 85 psFree(*stamps); 86 *stamps = pmSubtractionSetStamps(xStamp, yStamp, fluxStamp, subMask, region); 87 } else { 88 psTrace("psModules.imcombine", 3, "Finding stamps...\n"); 89 *stamps = pmSubtractionFindStamps(*stamps, reference->image, subMask, region, 90 threshold, stampSpacing); 91 } 92 if (!stamps) { 86 // Apply exclusion zone if we're matching to a nominated PSF; otherwise don't care 87 *stamps = pmSubtractionSetStamps(xStamp, yStamp, fluxStamp, reference->image, subMask, 88 region, stampSpacing, input ? 0 : footprint); 89 } 90 psTrace("psModules.imcombine", 3, "Finding stamps...\n"); 91 *stamps = pmSubtractionFindStamps(*stamps, reference->image, subMask, region, 92 threshold, stampSpacing); 93 if (!*stamps) { 93 94 psError(PS_ERR_UNKNOWN, false, "Unable to find stamps."); 94 95 return false; … … 220 221 psRegion *region = NULL; // Iso-kernel region 221 222 psString regionString = NULL; // String for region 222 p sArray *stamps = NULL;// Stamps for matching PSF223 pmSubtractionStampList *stamps = NULL; // Stamps for matching PSF 223 224 psVector *solution = NULL; // Solution to match PSF 224 225 pmSubtractionKernels *kernels = NULL; // Kernel basis functions -
trunk/psModules/src/imcombine/pmSubtractionParams.c
r14738 r14801 171 171 172 172 pmSubtractionKernels *pmSubtractionKernelsOptimumISIS(pmSubtractionKernelsType type, int size, int inner, 173 int spatialOrder, psVector *fwhms, int maxOrder, 174 const psArray *stamps, int footprint, float tolerance) 173 int spatialOrder, psVector *fwhms, int maxOrder, 174 const pmSubtractionStampList *stamps, int footprint, 175 float tolerance) 175 176 { 176 177 if (type != PM_SUBTRACTION_KERNEL_ISIS && type != PM_SUBTRACTION_KERNEL_GUNK) { … … 184 185 PS_ASSERT_VECTOR_TYPE(fwhms, PS_TYPE_F32, NULL); 185 186 PS_ASSERT_INT_NONNEGATIVE(maxOrder, NULL); 186 P S_ASSERT_ARRAY_NON_NULL(stamps, NULL);187 PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL); 187 188 PS_ASSERT_INT_NONNEGATIVE(footprint, NULL); 188 189 PS_ASSERT_FLOAT_LARGER_THAN(tolerance, 0.0, NULL); … … 208 209 209 210 // Need to save the stamp inputs --- we're changing the values! 210 int numStamps = stamps->n ;// Number of stamps211 int numStamps = stamps->num; // Number of stamps 211 212 psArray *inputs = psArrayAlloc(numStamps); // Deep copies of the inputs 213 psVector *badStamps = psVectorAlloc(numStamps, PS_TYPE_U8); // Mark the bad stamps 214 psVectorInit(badStamps, 0); 212 215 for (int i = 0; i < numStamps; i++) { 213 pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest 214 if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED || stamp->status == PM_SUBTRACTION_STAMP_NONE) { 216 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest 217 if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE && stamp->status != PM_SUBTRACTION_STAMP_USED) { 218 badStamps->data.U8[i] = 0xff; 215 219 continue; 216 220 } … … 230 234 int numPixels = 0; // Number of pixels contributing to chi^2 231 235 for (int i = 0; i < numStamps; i++) { 232 pmSubtractionStamp *stamp = stamps-> data[i]; // Stamp of interest233 if ( stamp->status == PM_SUBTRACTION_STAMP_REJECTED || stamp->status == PM_SUBTRACTION_STAMP_NONE) {236 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest 237 if (badStamps->data.U8[i]) { 234 238 continue; 235 239 } … … 238 242 psFree(inputs); 239 243 psFree(kernels); 244 psFree(badStamps); 240 245 return NULL; 241 246 } … … 255 260 psFree(inputs); 256 261 psFree(kernels); 262 psFree(badStamps); 257 263 return NULL; 258 264 } … … 287 293 288 294 for (int j = 0; j < numStamps; j++) { 289 pmSubtractionStamp *stamp = stamps->data[j]; // Stamp of interest 290 if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED || 291 stamp->status == PM_SUBTRACTION_STAMP_NONE) { 295 if (badStamps->data.U8[j]) { 292 296 continue; 293 297 } 298 pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest 294 299 accumulateCross(&sumI, &sumII, &sumIC, stamp, inputs->data[j], i, footprint); 295 300 } … … 301 306 double chi2 = 0.0; // Chi^2 302 307 for (int j = 0; j < numStamps; j++) { 303 pmSubtractionStamp *stamp = stamps->data[j]; // Stamp of interest 304 if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED || 305 stamp->status == PM_SUBTRACTION_STAMP_NONE) { 308 if (badStamps->data.U8[j]) { 306 309 continue; 307 310 } 311 pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest 308 312 chi2 += accumulateChi2(inputs->data[j], stamp, i, coeff, bg, footprint); 309 313 } … … 329 333 psFree(ranking); 330 334 psFree(kernels); 335 psFree(badStamps); 331 336 return NULL; 332 337 } … … 336 341 // Remove its contribution, and don't include it in the future. 337 342 for (int j = 0; j < numStamps; j++) { 338 pmSubtractionStamp *stamp = stamps->data[j]; // Stamp of interest 339 if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED || 340 stamp->status == PM_SUBTRACTION_STAMP_NONE) { 343 if (badStamps->data.U8[j]) { 341 344 continue; 342 345 } 346 pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest 343 347 subtractConvolution(inputs->data[j], stamp, bestIndex, bestCoeff, bestBG, footprint); 344 348 } … … 365 369 psFree(ranking); 366 370 psFree(kernels); 371 psFree(badStamps); 367 372 return NULL; 368 373 } … … 376 381 psArray *convNew = psArrayAlloc(numStamps); 377 382 for (int i = 0; i < numStamps; i++) { 383 if (badStamps->data.U8[i]) { 384 continue; 385 } 378 386 convNew->data[i] = psArrayAlloc(newSize); 379 387 } … … 388 396 389 397 for (int j = 0; j < numStamps; j++) { 398 if (badStamps->data.U8[j]) { 399 continue; 400 } 401 pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest 390 402 psArray *convolutions = convNew->data[j]; // Convolutions for this stamp 391 pmSubtractionStamp *stamp = stamps->data[j]; // Stamp of interest392 403 convolutions->data[rank] = psMemIncrRefCounter(stamp->convolutions->data[i]); 393 404 } … … 405 416 406 417 for (int i = 0; i < numStamps; i++) { 407 pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest 418 if (badStamps->data.U8[i]) { 419 continue; 420 } 421 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest 408 422 psFree(stamp->convolutions); 409 423 stamp->convolutions = convNew->data[i]; 410 424 } 411 425 426 psFree(badStamps); 412 427 psFree(ranking); 413 428 -
trunk/psModules/src/imcombine/pmSubtractionParams.h
r14671 r14801 3 3 4 4 #include <pslib.h> 5 #include "pmSubtractionKernels.h" 5 #include <pmSubtractionKernels.h> 6 #include <pmSubtractionStamps.h> 6 7 7 8 /// Generate a set of optimum kernels for ISIS (or GUNK) … … 12 13 psVector *fwhms, ///< Gaussian FWHMs to try 13 14 int maxOrder, ///< Maximum polynomial order 14 const p sArray*stamps, ///< Stamps15 const pmSubtractionStampList *stamps, ///< Stamps 15 16 int footprint, ///< Convolution footprint for stamps 16 17 float tolerance ///< Maximum difference in chi^2 -
trunk/psModules/src/imcombine/pmSubtractionStamps.c
r14701 r14801 9 9 #include "pmSubtractionStamps.h" 10 10 11 #define STAMP_LIST_BUFFER 20 // Number of stamps to add to list at a time 12 11 13 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 12 14 // Private (file-static) functions 13 15 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 14 16 17 // Free function for pmSubtractionStampList 18 static void subtractionStampListFree(pmSubtractionStampList *list // Stamp list to free 19 ) 20 { 21 psFree(list->stamps); 22 psFree(list->regions); 23 psFree(list->x); 24 psFree(list->y); 25 psFree(list->flux); 26 } 27 15 28 // Free function for pmSubtractionStamp 16 void subtractionStampFree(pmSubtractionStamp *stamp // Stamp to free17 )29 static void subtractionStampFree(pmSubtractionStamp *stamp // Stamp to free 30 ) 18 31 { 19 32 psFree(stamp->matrix); … … 28 41 static bool checkStampRegion(int x, int y, // Coordinates of stamp 29 42 const psRegion *region // Region of interest 30 )43 ) 31 44 { 32 45 if (!region) { … … 40 53 static bool checkStampMask(int x, int y, // Coordinates of stamp 41 54 const psImage *mask 42 )55 ) 43 56 { 44 57 if (!mask) { … … 54 67 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 55 68 56 pmSubtractionStamp *pmSubtractionStampAlloc(pmSubtractionStampStatus status) 69 pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, int numRows, const psRegion *region, 70 float spacing) 71 { 72 pmSubtractionStampList *list = psAlloc(sizeof(pmSubtractionStampList)); // Stamp list to return 73 psMemSetDeallocator(list, (psFreeFunc)subtractionStampListFree); 74 75 // Get region in which to find stamps: [xMin:xMax,yMin:yMax] 76 int xMin = 0, xMax = numCols, yMin = 0, yMax = numRows; 77 if (region) { 78 xMin = PS_MAX(region->x0, xMin); 79 xMax = PS_MIN(region->x1, xMax); 80 yMin = PS_MAX(region->y0, yMin); 81 yMax = PS_MIN(region->y1, yMax); 82 } 83 int xSize = xMax - xMin, ySize = yMax - yMin; // Size of region of interest 84 int xStamps = xSize / spacing + 1, yStamps = ySize / spacing + 1; // Number of stamps in x and y 85 86 list->num = xStamps * yStamps; 87 list->stamps = psArrayAlloc(list->num); 88 list->regions = psArrayAlloc(list->num); 89 90 for (int y = 0, index = 0; y < yStamps; y++) { 91 int yStart = yMin + y * ySize / (yStamps + 1); // Subregion starts here 92 int yStop = yMin + (y + 1) * ySize / (yStamps + 1) - 1; // Subregion stops here 93 assert(yStart >= yMin && yStop < yMax); 94 95 for (int x = 0; x < xStamps; x++, index++) { 96 int xStart = xMin + x * xSize / (xStamps + 1); // Subregion starts here 97 int xStop = xMin + (x + 1) * xSize / (xStamps + 1) - 1; // Subregion stops here 98 assert(xStart >= xMin && xStop < xMax); 99 100 list->stamps->data[index] = pmSubtractionStampAlloc(); 101 list->regions->data[index] = psRegionAlloc(xStart, xStop, yStart, yStop); 102 } 103 } 104 105 list->x = NULL; 106 list->y = NULL; 107 list->flux = NULL; 108 109 return list; 110 } 111 112 pmSubtractionStamp *pmSubtractionStampAlloc(void) 57 113 { 58 114 pmSubtractionStamp *stamp = psAlloc(sizeof(pmSubtractionStamp)); // Stamp to return … … 66 122 stamp->matrix = NULL; 67 123 stamp->vector = NULL; 68 stamp->status = status;124 stamp->status = PM_SUBTRACTION_STAMP_INIT; 69 125 70 126 stamp->reference = NULL; … … 77 133 78 134 79 psArray *pmSubtractionFindStamps(psArray *stamps, const psImage *image, const psImage *subMask, 80 const psRegion *region, float threshold, float spacing) 135 pmSubtractionStampList *pmSubtractionFindStamps(pmSubtractionStampList *stamps, const psImage *image, 136 const psImage *subMask, const psRegion *region, 137 float threshold, float spacing) 81 138 { 82 139 PS_ASSERT_IMAGE_NON_NULL(image, NULL); … … 107 164 int numRows = image->numRows, numCols = image->numCols; // Size of image 108 165 109 // Get region in which to find stamps: [xMin:xMax,yMin:yMax] 110 int xMin = 0, xMax = numCols, yMin = 0, yMax = numRows; 111 if (region) { 112 xMin = PS_MAX(region->x0, xMin); 113 xMax = PS_MIN(region->x1, xMax); 114 yMin = PS_MAX(region->y0, yMin); 115 yMax = PS_MIN(region->y1, yMax); 116 } 117 int xSize = xMax - xMin, ySize = yMax - yMin; // Size of region of interest 118 119 // Number of stamps 120 int xNumStamps = xSize / spacing + 1; 121 int yNumStamps = ySize / spacing + 1; 122 166 if (!stamps) { 167 stamps = pmSubtractionStampListAlloc(numCols, numRows, region, spacing); 168 } 169 170 int numStamps = stamps->num; // Number of stamp regions 123 171 int numFound = 0; // Number of stamps found 124 125 if (stamps) { 126 PS_ASSERT_INT_EQUAL(stamps->n, xNumStamps * yNumStamps, NULL); 127 // Just in case, check for NULL stamps 128 for (int i = 0; i < xNumStamps * yNumStamps; i++) { 129 if (!stamps->data[i]) { 130 stamps->data[i] = pmSubtractionStampAlloc(PM_SUBTRACTION_STAMP_REJECTED); 172 for (int i = 0; i < numStamps; i++) { 173 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest 174 175 // Only find a new stamp if we need to 176 if (stamp->status != PM_SUBTRACTION_STAMP_REJECTED && stamp->status != PM_SUBTRACTION_STAMP_INIT) { 177 continue; 178 } 179 180 float xStamp = 0, yStamp = 0; // Coordinates of stamp 181 float fluxStamp = NAN; // Flux of stamp 182 bool goodStamp = false; // Found a good stamp? 183 184 // A couple different ways of finding stamps: 185 if (stamps->x && stamps->y) { 186 // Get the next stamp from the list 187 psVector *xList = stamps->x->data[i], *yList = stamps->y->data[i]; // Coordinate lists 188 psVector *fluxList = stamps->flux->data[i]; // List of stamp fluxes 189 190 // Take stamp off the top of the (sorted) list 191 if (xList->n > 0) { 192 int index = xList->n - 1; // Index of new stamp 193 xStamp = xList->data.F32[index]; 194 yStamp = yList->data.F32[index]; 195 fluxStamp = fluxList->data.F32[index]; 196 197 // Chop off the top of the list 198 xList->n = index; 199 yList->n = index; 200 fluxList->n = index; 201 202 goodStamp = true; 131 203 } 132 } 133 } else { 134 stamps = psArrayAlloc(xNumStamps * yNumStamps); 135 for (int i = 0; i < xNumStamps * yNumStamps ; i++) { 136 stamps->data[i] = pmSubtractionStampAlloc(PM_SUBTRACTION_STAMP_INIT); 137 } 138 } 139 140 for (int j = 0, index = 0; j < yNumStamps; j++) { 141 for (int i = 0; i < xNumStamps; i++, index++) { 142 pmSubtractionStamp *stamp = stamps->data[index]; // Stamp of interest 143 144 // Only find a new stamp if we need to 145 if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED || 146 stamp->status == PM_SUBTRACTION_STAMP_INIT) { 147 148 // Find maximum non-masked value in the image section, 149 // but don't include a footprint around the edge 150 float fluxBest = threshold; // Flux of best stamp pixel 151 int xBest = 0, yBest = 0; // Coordinates of best stamp 152 153 // Bounds of region to search for stamp 154 int yStart = yMin + j * ySize / (yNumStamps + 1); 155 int yStop = yMin + (j + 1) * ySize / (yNumStamps + 1) - 1; 156 int xStart = xMin + i * xSize / (xNumStamps + 1); 157 int xStop = xMin + (i + 1) * xSize / (xNumStamps + 1) - 1; 158 assert(yStop < numRows && xStop < numCols && yStart >= 0 && xStart >= 0); 159 160 for (int y = yStart; y <= yStop ; y++) { 161 for (int x = xStart; x <= xStop ; x++) { 162 if (checkStampMask(x, y, subMask) && image->data.F32[y][x] > fluxBest) { 163 fluxBest = image->data.F32[y][x]; 164 xBest = x; 165 yBest = y; 166 } 204 } else { 205 // Use a simple method of automatically finding stamps --- take the highest pixel in the subregion 206 fluxStamp = threshold; 207 psRegion *subRegion = stamps->regions->data[i]; // Sub-region of interest 208 for (int y = subRegion->y0; y <= subRegion->y1 ; y++) { 209 for (int x = subRegion->x0; x <= subRegion->y1 ; x++) { 210 if (checkStampMask(x, y, subMask) && image->data.F32[y][x] > fluxStamp) { 211 fluxStamp = image->data.F32[y][x]; 212 xStamp = x; 213 yStamp = y; 214 goodStamp = true; 167 215 } 168 216 } 169 170 stamp->x = xBest;171 stamp->y = yBest;172 stamp->flux = fluxBest;173 174 // Reset the postage stamps since we're making a new stamp175 psFree(stamp->reference);176 psFree(stamp->input);177 psFree(stamp->weight);178 stamp->reference = stamp->input = stamp->weight = NULL;179 180 if (fluxBest > threshold) {181 stamp->status = PM_SUBTRACTION_STAMP_CALCULATE;182 numFound++;183 psTrace("psModules.imcombine", 5, "Found stamp in region (%d,%d): %d,%d\n",184 i, j, (int)stamp->x, (int)stamp->y);185 } else {186 stamp->status = PM_SUBTRACTION_STAMP_NONE;187 }188 217 } 189 218 } 190 } 191 192 psLogMsg("psModules.imcombine", PS_LOG_INFO, "Found %d stamps in %dx%d regions", 193 numFound, xNumStamps, yNumStamps); 219 220 if (goodStamp) { 221 stamp->x = xStamp; 222 stamp->y = yStamp; 223 stamp->flux = fluxStamp; 224 225 // Reset the postage stamps since we're making a new stamp 226 psFree(stamp->reference); 227 psFree(stamp->input); 228 psFree(stamp->weight); 229 stamp->reference = stamp->input = stamp->weight = NULL; 230 231 stamp->status = PM_SUBTRACTION_STAMP_CALCULATE; 232 numFound++; 233 psTrace("psModules.imcombine", 5, "Found stamp in subregion %d: %d,%d\n", 234 i, (int)stamp->x, (int)stamp->y); 235 } else { 236 stamp->status = PM_SUBTRACTION_STAMP_NONE; 237 } 238 } 239 240 psLogMsg("psModules.imcombine", PS_LOG_INFO, "Found %d stamps", numFound); 194 241 195 242 return stamps; … … 197 244 198 245 199 psArray *pmSubtractionSetStamps(const psVector *x, const psVector *y, const psVector *flux, 200 const psImage *subMask, const psRegion *region) 246 pmSubtractionStampList *pmSubtractionSetStamps(const psVector *x, const psVector *y, const psVector *flux, 247 const psImage *image, const psImage *subMask, 248 const psRegion *region, float spacing, int exclusionZone) 201 249 { 202 250 PS_ASSERT_VECTOR_NON_NULL(x, false); … … 209 257 PS_ASSERT_VECTOR_TYPE(flux, PS_TYPE_F32, false); 210 258 PS_ASSERT_VECTORS_SIZE_EQUAL(flux, x, false); 259 } else { 260 PS_ASSERT_IMAGE_NON_NULL(image, false); 211 261 } 212 262 if (subMask) { 213 263 PS_ASSERT_IMAGE_NON_NULL(subMask, false); 214 264 PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, false); 215 } 216 217 int num = x->n; // Number of stamps 218 psArray *stamps = psArrayAllocEmpty(num); // Stamps, to return 219 220 for (int i = 0; i < num; i++) { 221 int xStamp = x->data.F32[i] + 0.5, yStamp = y->data.F32[i] + 0.5; // Pixel coords of stamp 222 if (!checkStampRegion(xStamp, yStamp, region)) { 223 continue; 224 } 225 if (!checkStampMask(xStamp, yStamp, subMask)) { 226 continue; 227 } 228 229 pmSubtractionStamp *stamp = pmSubtractionStampAlloc(PM_SUBTRACTION_STAMP_CALCULATE); // Stamp 230 stamp->x = x->data.F32[i]; 231 stamp->y = y->data.F32[i]; 232 if (flux) { 233 stamp->flux = flux->data.F32[i]; 234 } 235 stamps->data[stamps->n++] = stamp; 236 } 265 if (image) { 266 PS_ASSERT_IMAGE_NON_NULL(image, false); 267 PS_ASSERT_IMAGES_SIZE_EQUAL(image, subMask, false); 268 } 269 } 270 271 int numStars = x->n; // Number of stars 272 psVector *exclude = psVectorAlloc(numStars, PS_TYPE_U8); // Exclude a star? 273 psVectorInit(exclude, 0); 274 275 // Apply exclusion zone around each star --- important when we're convolving to a specified PSF; less so 276 // when we're trying to get a good subtraction. 277 if (exclusionZone > 0) { 278 psTrace("psModules.imcombine", 2, "Applying exclusion zone of %d pixels for stamps\n", exclusionZone); 279 // We use something resembling a binary search --- coordinates are sorted in the x dimension, and then 280 // to exclude stars within a nominated distance, we need only examine (i.e., calculate the 281 // 2-dimensional distance to) other stars in the list that are within that distance of the x 282 // coordinate of the star of interest. By marking both stars when we find a violation of the 283 // exclusion zone, we need only search upwards from the x coordinate of the star of interest. 284 285 int minDistance2 = PS_SQR(exclusionZone); // Minimum square distance for other stars 286 psVector *indexes = psVectorSortIndex(NULL, x); // Indices for sorting in x 287 for (int i = 0; i < numStars; i++) { 288 int iIndex = indexes->data.S32[i]; // Index for star i 289 if (exclude->data.U8[iIndex]) { 290 continue; 291 } 292 float ix = x->data.F32[iIndex], iy = y->data.F32[iIndex]; // Coordinates for star i 293 float jx, jy; // Coordinates for star j 294 for (int j = i + 1, jIndex = indexes->data.S32[j]; 295 (jx = x->data.F32[jIndex]) < ix + exclusionZone && j < numStars; 296 j++, jIndex = indexes->data.S32[j]) { 297 jy = y->data.F32[jIndex]; 298 if (PS_SQR(ix - jx) + PS_SQR(iy - jy) < minDistance2) { 299 exclude->data.U8[iIndex] = 0xff; 300 exclude->data.U8[jIndex] = 0xff; 301 psTrace("psModules.imcombine", 5, "Excluding stamps %d,%d and %d,%d\n", 302 (int)x->data.F32[iIndex], (int)y->data.F32[iIndex], 303 (int)x->data.F32[jIndex], (int)y->data.F32[jIndex]); 304 // Would 'break' here, but there might be more stamps within the exclusion zone. 305 } 306 } 307 } 308 psFree(indexes); 309 } 310 311 pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows, 312 region, spacing); // Stamp list 313 int numStamps = stamps->num; // Number of stamps 314 315 // Initialise the lists 316 stamps->x = psArrayAlloc(numStamps); 317 stamps->y = psArrayAlloc(numStamps); 318 stamps->flux = psArrayAlloc(numStamps); 319 for (int i = 0; i < numStamps; i++) { 320 stamps->x->data[i] = psVectorAllocEmpty(STAMP_LIST_BUFFER, PS_TYPE_F32); 321 stamps->y->data[i] = psVectorAllocEmpty(STAMP_LIST_BUFFER, PS_TYPE_F32); 322 stamps->flux->data[i] = psVectorAllocEmpty(STAMP_LIST_BUFFER, PS_TYPE_F32); 323 } 324 325 // Put the stars into their appropriate subregions 326 for (int i = 0; i < numStars; i++) { 327 if (exclude->data.U8[i]) { 328 // Star has been excluded 329 continue; 330 } 331 float xStamp = x->data.F32[i], yStamp = y->data.F32[i]; // Coordinates of stamp 332 int xPix = xStamp + 0.5, yPix = yStamp + 0.5; // Pixel coordinate of stamp 333 if (!checkStampRegion(xPix, yPix, region)) { 334 // It's not in the big region 335 continue; 336 } 337 if (!checkStampMask(xPix, yPix, subMask)) { 338 // Not a good stamp 339 continue; 340 } 341 342 for (int j = 0; j < numStamps; j++) { 343 psRegion *subRegion = stamps->regions->data[j]; // Subregion of interest 344 if (checkStampRegion(xPix, yPix, subRegion)) { 345 psVector *xList = stamps->x->data[j], *yList = stamps->y->data[j]; // Pixel lists 346 psVector *fluxList = stamps->flux->data[j]; // Flux list 347 348 int index = xList->n; // Index of new stamp candidate 349 xList->data.F32[index] = xStamp; 350 yList->data.F32[index] = yStamp; 351 352 if (flux) { 353 fluxList->data.F32[index] = flux->data.F32[i]; 354 } else { 355 fluxList->data.F32[index] = image->data.F32[yPix][xPix]; 356 } 357 358 psVectorExtend(xList, STAMP_LIST_BUFFER, 1); 359 psVectorExtend(yList, STAMP_LIST_BUFFER, 1); 360 psVectorExtend(fluxList, STAMP_LIST_BUFFER, 1); 361 362 break; 363 } 364 } 365 } 366 psFree(exclude); 367 368 // Sort the list by flux, with the brightest last 369 for (int i = 0; i < numStamps; i++) { 370 psVector *xList = stamps->x->data[i], *yList = stamps->y->data[i]; // Pixel lists 371 psVector *fluxList = stamps->flux->data[i]; // Flux list 372 373 psVector *indexes = psVectorSortIndex(NULL, fluxList); // Indices to sort flux 374 int num = indexes->n; // Number of candidate stamps in this subregion 375 376 psVector *xSorted = psVectorAlloc(num, PS_TYPE_F32); // Sorted version of x list 377 psVector *ySorted = psVectorAlloc(num, PS_TYPE_F32); // Sorted version of y list 378 psVector *fluxSorted = psVectorAlloc(num, PS_TYPE_F32); // Sorted version of flux list 379 for (int j = 0; j < num; j++) { 380 int k = indexes->data.S32[j]; // Sorted index 381 xSorted->data.F32[j] = xList->data.F32[k]; 382 ySorted->data.F32[j] = yList->data.F32[k]; 383 fluxSorted->data.F32[j] = fluxList->data.F32[k]; 384 } 385 psFree(indexes); 386 387 psFree(stamps->x->data[i]); 388 psFree(stamps->y->data[i]); 389 psFree(stamps->flux->data[i]); 390 391 stamps->x->data[i] = xSorted; 392 stamps->y->data[i] = ySorted; 393 stamps->flux->data[i] = fluxSorted; 394 } 395 237 396 238 397 return stamps; … … 240 399 241 400 242 bool pmSubtractionExtractStamps(p sArray *stamps, psImage *reference, psImage *input, psImage *weight,243 int footprint, int kernelSize)244 { 245 P S_ASSERT_ARRAY_NON_NULL(stamps, false);401 bool pmSubtractionExtractStamps(pmSubtractionStampList *stamps, psImage *reference, psImage *input, 402 psImage *weight, int footprint, int kernelSize) 403 { 404 PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); 246 405 PS_ASSERT_IMAGE_NON_NULL(reference, false); 247 406 PS_ASSERT_IMAGE_TYPE(reference, PS_TYPE_F32, false); … … 263 422 264 423 if (!weight) { 265 // Use the input as a rough approximation to the variance map, and HOPE that it's positive. 424 // Use the input (or if I must, the reference) as a rough approximation to the variance map, and HOPE 425 // that it's positive. 266 426 weight = input ? input : reference; 267 427 } 268 428 269 for (int i = 0; i < stamps->n ; i++) {270 pmSubtractionStamp *stamp = stamps-> data[i]; // Stamp of interest429 for (int i = 0; i < stamps->num; i++) { 430 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest 271 431 if (!stamp || stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) { 272 432 continue; … … 307 467 308 468 309 bool pmSubtractionGenerateStamps(p sArray*stamps, float fwhm, int footprint, int kernelSize)310 { 311 P S_ASSERT_ARRAY_NON_NULL(stamps, false);469 bool pmSubtractionGenerateStamps(pmSubtractionStampList *stamps, float fwhm, int footprint, int kernelSize) 470 { 471 PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); 312 472 PS_ASSERT_FLOAT_LARGER_THAN(fwhm, 0.0, false); 313 473 PS_ASSERT_INT_NONNEGATIVE(footprint, false); 474 PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); 314 475 315 476 int size = kernelSize + footprint; // Size of postage stamps 316 int num = stamps->n ;// Number of stamps477 int num = stamps->num; // Number of stamps 317 478 float sigma = fwhm / (2.0 * sqrtf(2.0 * log(2.0))); // Gaussian sigma 318 479 319 480 for (int i = 0; i < num; i++) { 320 pmSubtractionStamp *stamp = stamps-> data[i]; // Stamp of interest481 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest 321 482 if (!(stamp->status & PM_SUBTRACTION_STAMP_CALCULATE)) { 322 483 continue; … … 347 508 } 348 509 349 return stamps;350 } 510 return true; 511 } -
trunk/psModules/src/imcombine/pmSubtractionStamps.h
r14701 r14801 15 15 } pmSubtractionStampStatus; 16 16 17 /// A list of stamps 18 typedef struct { 19 long num; ///< Number of stamps 20 psArray *stamps; ///< The stamps 21 psArray *regions; ///< Regions for each stamp 22 psArray *x, *y; ///< Coordinates for possible stamps (or NULL) 23 psArray *flux; ///< Fluxes for possible stamps (or NULL) 24 } pmSubtractionStampList; 25 26 /// Allocate a list of stamps 27 pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, // Number of columns in image 28 int numRows, // Number of rows in image 29 const psRegion *region, // Region for stamps, or NULL 30 float spacing // Rough average spacing between stamps 31 ); 32 33 /// Assertion for stamp list to be valid 34 #define PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(LIST, RETURNVALUE) { \ 35 PS_ASSERT_PTR_NON_NULL(LIST, RETURNVALUE); \ 36 PS_ASSERT_ARRAY_NON_NULL((LIST)->stamps, RETURNVALUE); \ 37 PS_ASSERT_ARRAY_NON_NULL((LIST)->regions, RETURNVALUE); \ 38 PS_ASSERT_INT_POSITIVE((LIST)->num, RETURNVALUE); \ 39 PS_ASSERT_ARRAY_SIZE((LIST)->stamps, (LIST)->num, RETURNVALUE); \ 40 PS_ASSERT_ARRAY_SIZE((LIST)->regions, (LIST)->num, RETURNVALUE); \ 41 if ((LIST)->x || (LIST)->y || (LIST)->flux) { \ 42 PS_ASSERT_ARRAY_NON_NULL((LIST)->x, RETURNVALUE); \ 43 PS_ASSERT_ARRAY_NON_NULL((LIST)->y, RETURNVALUE); \ 44 PS_ASSERT_ARRAY_NON_NULL((LIST)->flux, RETURNVALUE); \ 45 PS_ASSERT_ARRAY_SIZE((LIST)->x, (LIST)->num, RETURNVALUE); \ 46 PS_ASSERT_ARRAY_SIZE((LIST)->y, (LIST)->num, RETURNVALUE); \ 47 PS_ASSERT_ARRAY_SIZE((LIST)->flux, (LIST)->num, RETURNVALUE); \ 48 } \ 49 } 50 17 51 /// A stamp for image subtraction 18 52 typedef struct { … … 26 60 psImage *matrix; ///< Associated matrix 27 61 psVector *vector; ///< Assoicated vector 28 pmSubtractionStampStatus status; ///< Status of stamp62 pmSubtractionStampStatus status; ///< Status of stamp 29 63 } pmSubtractionStamp; 30 64 65 /// Allocate a stamp 66 pmSubtractionStamp *pmSubtractionStampAlloc(void); 67 31 68 /// Find stamps on an image 32 p sArray *pmSubtractionFindStamps(psArray*stamps, ///< Output stamps, or NULL33 const psImage *image, ///< Image for which to find stamps34 const psImage *mask, ///< Mask, or NULL35 const psRegion *region, ///< Region in which to search for stamps, or NULL36 float threshold, ///< Threshold for stamps in the image37 float spacing ///< Rough spacing for stamps69 pmSubtractionStampList *pmSubtractionFindStamps(pmSubtractionStampList *stamps, ///< Output stamps, or NULL 70 const psImage *image, ///< Image for which to find stamps 71 const psImage *mask, ///< Mask, or NULL 72 const psRegion *region, ///< Region to search, or NULL 73 float threshold, ///< Threshold for stamps in the image 74 float spacing ///< Rough spacing for stamps 38 75 ); 39 76 40 77 /// Set stamps based on a list of x,y 41 psArray *pmSubtractionSetStamps(const psVector *x, ///< x coordinates for each stamp 42 const psVector *y, ///< y coordinates for each stamp 43 const psVector *flux, ///< Flux for each stamp, or NULL 44 const psImage *mask, ///< Mask, or NULL 45 const psRegion *region ///< Region in which to search for stamps, or NULL 78 /// 79 /// We may optionally apply an exclusion zone around each star --- this is important when we're convolving to 80 /// a specified PSF; less so when we're only trying to get a good subtraction. 81 pmSubtractionStampList *pmSubtractionSetStamps(const psVector *x, ///< x coordinates for each stamp 82 const psVector *y, ///< y coordinates for each stamp 83 const psVector *flux, ///< Flux for each stamp, or NULL 84 const psImage *image, ///< Image for flux of stamp 85 const psImage *mask, ///< Mask, or NULL 86 const psRegion *region, ///< Region to search, or NULL 87 float spacing, ///< Rough spacing for stamps 88 int exclusionZone // Exclusion zone around each star 46 89 ); 47 90 48 91 /// Extract stamps from the images 49 bool pmSubtractionExtractStamps(p sArray*stamps, ///< Stamps92 bool pmSubtractionExtractStamps(pmSubtractionStampList *stamps, ///< Stamps 50 93 psImage *reference, ///< Reference image 51 94 psImage *input, ///< Input image (or NULL) … … 56 99 57 100 /// Generate target for stamps based on list 58 bool pmSubtractionGenerateStamps(p sArray*stamps, ///< Stamps101 bool pmSubtractionGenerateStamps(pmSubtractionStampList *stamps, ///< Stamps 59 102 float fwhm, ///< FWHM for each stamp 60 103 int footprint, ///< Stamp footprint size
Note:
See TracChangeset
for help on using the changeset viewer.
