Changeset 26035
- Timestamp:
- Nov 4, 2009, 4:52:05 PM (17 years ago)
- Location:
- trunk/psModules/src/imcombine
- Files:
-
- 7 edited
-
pmSubtraction.c (modified) (17 diffs)
-
pmSubtractionEquation.c (modified) (15 diffs)
-
pmSubtractionMatch.c (modified) (12 diffs)
-
pmSubtractionMatch.h (modified) (2 diffs)
-
pmSubtractionParams.c (modified) (10 diffs)
-
pmSubtractionStamps.c (modified) (17 diffs)
-
pmSubtractionStamps.h (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/pmSubtraction.c
r25999 r26035 29 29 #define PIXEL_LIST_BUFFER 100 // Number of entries to add to pixel list at a time 30 30 #define MIN_SAMPLE_STATS 7 // Minimum number to use sample statistics; otherwise use quartiles 31 #define USE_ SYS_ERR // Use systematicerror image?31 #define USE_KERNEL_ERR // Use kernel error image? 32 32 33 33 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// … … 266 266 static void convolveVarianceFFT(psImage *target,// Place the result in here 267 267 psImage *variance, // Variance map to convolve 268 psImage * sys, // Systematicerror image268 psImage *kernelErr, // Kernel error image 269 269 psImage *mask, // Mask image 270 270 psImageMaskType maskVal, // Value to mask … … 278 278 279 279 psImage *subVariance = variance ? psImageSubset(variance, border) : NULL; // Variance map 280 psImage *sub Sys = sys ? psImageSubset(sys, border) : NULL; // Systematicerror image280 psImage *subKE = kernelErr ? psImageSubset(kernelErr, border) : NULL; // Kernel error image 281 281 psImage *subMask = mask ? psImageSubset(mask, border) : NULL; // Mask 282 282 283 283 // XXX Can trim this a little by combining the convolution: only have to take the FFT of the kernel once 284 284 psImage *convVariance = psImageConvolveFFT(NULL, subVariance, subMask, maskVal, kernel); // Convolved variance 285 psImage *conv Sys = subSys ? psImageConvolveFFT(NULL, subSys, subMask, maskVal, kernel) : NULL; // Conv sys285 psImage *convKE = subKE ? psImageConvolveFFT(NULL, subKE, subMask, maskVal, kernel) : NULL; // Conv KE 286 286 287 287 psFree(subVariance); 288 psFree(sub Sys);288 psFree(subKE); 289 289 psFree(subMask); 290 290 291 291 // Now, we have to stick it in where it belongs 292 292 int xMin = region.x0, xMax = region.x1, yMin = region.y0, yMax = region.y1; // Bounds of region 293 if (conv Sys) {293 if (convKE) { 294 294 for (int yTarget = yMin, ySource = size; yTarget < yMax; yTarget++, ySource++) { 295 295 for (int xTarget = xMin, xSource = size; xTarget < xMax; xTarget++, xSource++) { 296 296 target->data.F32[yTarget][xTarget] = convVariance->data.F32[ySource][xSource] + 297 conv Sys->data.F32[ySource][xSource];297 convKE->data.F32[ySource][xSource]; 298 298 } 299 299 } … … 306 306 307 307 psFree(convVariance); 308 psFree(conv Sys);308 psFree(convKE); 309 309 310 310 return; … … 342 342 psImage *image, // Image to convolve 343 343 psImage *variance, // Variance map to convolve, or NULL 344 psImage * sys, // Systematicerror image, or NULL344 psImage *kernelErr, // Kernel error image, or NULL 345 345 psImage *subMask, // Subtraction mask 346 346 const pmSubtractionKernels *kernels, // Kernels … … 379 379 convolveFFT(convImage, image, subMask, subBad, *kernelImage, region, background, kernels->size); 380 380 if (variance) { 381 convolveVarianceFFT(convVariance, variance, sys, subMask, subBad, *kernelVariance, region, kernels->size); 381 convolveVarianceFFT(convVariance, variance, kernelErr, subMask, subBad, *kernelVariance, 382 region, kernels->size); 382 383 } 383 384 } else { … … 429 430 } 430 431 431 #ifdef USE_ SYS_ERR432 // Generate an image that can be used to track systematic errors 433 static psImage *subtraction SysErrImage(const psImage *image, // Image from which to make sys err image434 float sysError // Relative systematic error435 436 { 437 if (!isfinite( sysError) || sysError == 0.0) {432 #ifdef USE_KERNEL_ERR 433 // Generate an image that can be used to track systematic errors in the kernel 434 static psImage *subtractionKernelErrImage(const psImage *image, // Image from which to make kernel error image 435 float kernelError // Relative systematic error in kernel 436 ) 437 { 438 if (!isfinite(kernelError) || kernelError == 0.0) { 438 439 return NULL; 439 440 } 440 441 441 442 int numCols = image->numCols, numRows = image->numRows; // Size of image 442 psImage * sys = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Systematicerror image443 444 float sysError2 = PS_SQR(sysError); // Square of the systematicerror443 psImage *kernelErr = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Kernel error image 444 445 float kernelError2 = PS_SQR(kernelError); // Square of the kernel error 445 446 for (int y = 0; y < numRows; y++) { 446 447 for (int x = 0; x < numCols; x++) { 447 sys->data.F32[y][x] = PS_SQR(image->data.F32[y][x]) * sysError2;448 } 449 } 450 451 return sys;448 kernelErr->data.F32[y][x] = PS_SQR(image->data.F32[y][x]) * kernelError2; 449 } 450 } 451 452 return kernelErr; 452 453 } 453 454 #endif … … 892 893 psFree(stamp->image1); 893 894 psFree(stamp->image2); 894 psFree(stamp-> variance);895 stamp->image1 = stamp->image2 = stamp-> variance= NULL;895 psFree(stamp->weight); 896 stamp->image1 = stamp->image2 = stamp->weight = NULL; 896 897 psFree(stamp->matrix1); 897 898 psFree(stamp->matrix2); … … 1040 1041 psImage *convMask, // Output convolved mask 1041 1042 const pmReadout *ro1, const pmReadout *ro2, // Input readouts 1042 psImage * sys1, psImage *sys2, // Systematicerror images1043 psImage *kernelErr1, psImage *kernelErr2, // Kernel error images 1043 1044 psImage *subMask, // Input subtraction mask 1044 1045 psImageMaskType maskBad, // Mask value to give bad pixels … … 1066 1067 if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { 1067 1068 convolveRegion(out1->image, out1->variance, convMask, &kernelImage, &kernelVariance, 1068 ro1->image, ro1->variance, sys1, subMask, kernels, polyValues, background, *region,1069 maskBad, maskPoor, poorFrac, useFFT, false);1069 ro1->image, ro1->variance, kernelErr1, subMask, kernels, polyValues, background, 1070 *region, maskBad, maskPoor, poorFrac, useFFT, false); 1070 1071 } 1071 1072 if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { 1072 1073 convolveRegion(out2->image, out2->variance, convMask, &kernelImage, &kernelVariance, 1073 ro2->image, ro2->variance, sys2, subMask, kernels, polyValues, background, *region, 1074 maskBad, maskPoor, poorFrac, useFFT, kernels->mode == PM_SUBTRACTION_MODE_DUAL); 1074 ro2->image, ro2->variance, kernelErr2, subMask, kernels, polyValues, background, 1075 *region, maskBad, maskPoor, poorFrac, useFFT, 1076 kernels->mode == PM_SUBTRACTION_MODE_DUAL); 1075 1077 } 1076 1078 … … 1117 1119 const pmReadout *ro1 = args->data[7]; // Input readout 1 1118 1120 const pmReadout *ro2 = args->data[8]; // Input readout 2 1119 psImage * sys1 = args->data[9]; // Systematicerror image 11120 psImage * sys2 = args->data[10]; // Systematicerror image 21121 psImage *kernelErr1 = args->data[9]; // Kernel error image 1 1122 psImage *kernelErr2 = args->data[10]; // Kernel error image 2 1121 1123 psImage *subMask = args->data[11]; // Subtraction mask 1122 1124 psImageMaskType maskBad = PS_SCALAR_VALUE(args->data[12], PS_TYPE_IMAGE_MASK_DATA); // Output mask value for bad pixels … … 1128 1130 bool useFFT = PS_SCALAR_VALUE(args->data[18], U8); // Use FFT for convolution? 1129 1131 1130 return subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, sys1, sys2, 1131 subMask, maskBad, maskPoor, poorFrac, region, kernels, doBG, useFFT); 1132 return subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, kernelErr1, 1133 kernelErr2, subMask, maskBad, maskPoor, poorFrac, region, kernels, 1134 doBG, useFFT); 1132 1135 } 1133 1136 1134 1137 bool pmSubtractionConvolve(pmReadout *out1, pmReadout *out2, const pmReadout *ro1, const pmReadout *ro2, 1135 1138 psImage *subMask, int stride, psImageMaskType maskBad, psImageMaskType maskPoor, 1136 float poorFrac, float sysError, const psRegion *region,1139 float poorFrac, float kernelError, const psRegion *region, 1137 1140 const pmSubtractionKernels *kernels, bool doBG, bool useFFT) 1138 1141 { … … 1172 1175 PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(poorFrac, 0.0, false); 1173 1176 PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(poorFrac, 1.0, false); 1174 PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL( sysError, 0.0, false);1175 PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL( sysError, 1.0, false);1177 PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(kernelError, 0.0, false); 1178 PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(kernelError, 1.0, false); 1176 1179 if (region && psRegionIsNaN(*region)) { 1177 1180 psString string = psRegionToString(*region); … … 1232 1235 } 1233 1236 1234 psImage * sys1 = NULL, *sys2 = NULL; // Systematicerror images1235 #ifdef USE_ SYS_ERR1237 psImage *kernelErr1 = NULL, *kernelErr2 = NULL; // Kernel error images 1238 #ifdef USE_KERNEL_ERR 1236 1239 if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { 1237 sys1 = subtractionSysErrImage(ro1->image, sysError);1240 kernelErr1 = subtractionKernelErrImage(ro1->image, kernelError); 1238 1241 } 1239 1242 if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { 1240 sys2 = subtractionSysErrImage(ro2->image, sysError);1243 kernelErr2 = subtractionKernelErrImage(ro2->image, kernelError); 1241 1244 } 1242 1245 #endif … … 1289 1292 psArrayAdd(args, 1, (pmReadout*)ro1); // Casting away const 1290 1293 psArrayAdd(args, 1, (pmReadout*)ro2); // Casting away const 1291 psArrayAdd(args, 1, sys1);1292 psArrayAdd(args, 1, sys2);1294 psArrayAdd(args, 1, kernelErr1); 1295 psArrayAdd(args, 1, kernelErr2); 1293 1296 psArrayAdd(args, 1, subMask); 1294 1297 PS_ARRAY_ADD_SCALAR(args, maskBad, PS_TYPE_IMAGE_MASK); … … 1306 1309 psFree(job); 1307 1310 } else { 1308 subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, sys1, sys2,1309 subMask, maskBad, maskPoor, poorFrac, subRegion, kernels, doBG,1310 useFFT);1311 subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, 1312 kernelErr1, kernelErr2, subMask, maskBad, maskPoor, poorFrac, 1313 subRegion, kernels, doBG, useFFT); 1311 1314 } 1312 1315 psFree(subRegion); … … 1333 1336 psImageConvolveSetThreads(oldThreads); 1334 1337 1335 psFree( sys1);1336 psFree( sys2);1338 psFree(kernelErr1); 1339 psFree(kernelErr2); 1337 1340 1338 1341 // Calculate covariances -
trunk/psModules/src/imcombine/pmSubtractionEquation.c
r26031 r26035 17 17 //#define TESTING // TESTING output for debugging; may not work with threads! 18 18 19 #define USE_ VARIANCE // Include variancein equation?19 #define USE_WEIGHT // Include weight (1/variance) in equation? 20 20 21 21 … … 29 29 const psKernel *input, // Input image (target) 30 30 const psKernel *reference, // Reference image (convolution source) 31 const psKernel * variance, // Varianceimage31 const psKernel *weight, // Weight image 32 32 const psArray *convolutions, // Convolutions for each kernel 33 33 const pmSubtractionKernels *kernels, // Kernels … … 74 74 for (int x = - footprint; x <= footprint; x++) { 75 75 double cc = iConv->kernel[y][x] * jConv->kernel[y][x]; 76 #ifdef USE_ VARIANCE77 cc /= variance->kernel[y][x];76 #ifdef USE_WEIGHT 77 cc *= weight->kernel[y][x]; 78 78 #endif 79 79 sumCC += cc; … … 102 102 double rc = ref * conv; 103 103 double c = conv; 104 #ifdef USE_ VARIANCE105 float varVal = 1.0 / variance->kernel[y][x];106 ic *= varVal;107 rc *= varVal;108 c *= varVal;104 #ifdef USE_WEIGHT 105 float wtVal = weight->kernel[y][x]; 106 ic *= wtVal; 107 rc *= wtVal; 108 c *= wtVal; 109 109 #endif 110 110 sumIC += ic; … … 137 137 double rr = PS_SQR(ref); 138 138 double one = 1.0; 139 #ifdef USE_ VARIANCE140 float varVal = 1.0 / variance->kernel[y][x];141 rr *= varVal;142 ir *= varVal;143 in *= varVal;144 ref *= varVal;145 one *= varVal;139 #ifdef USE_WEIGHT 140 float wtVal = weight->kernel[y][x]; 141 rr *= wtVal; 142 ir *= wtVal; 143 in *= wtVal; 144 ref *= wtVal; 145 one *= wtVal; 146 146 #endif 147 147 sumRR += rr; … … 169 169 const psKernel *image1, // Image 1 170 170 const psKernel *image2, // Image 2 171 const psKernel * variance, // Varianceimage171 const psKernel *weight, // Weight image 172 172 const psArray *convolutions1, // Convolutions of image 1 for each kernel 173 173 const psArray *convolutions2, // Convolutions of image 2 for each kernel … … 227 227 double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x]; 228 228 double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x]; 229 #ifdef USE_VARIANCE 230 aa /= variance->kernel[y][x]; 231 bb /= variance->kernel[y][x]; 232 ab /= variance->kernel[y][x]; 229 #ifdef USE_WEIGHT 230 float wtVal = weight->kernel[y][x]; 231 aa *= wtVal; 232 bb *= wtVal; 233 ab *= wtVal; 233 234 #endif 234 235 sumAA += aa; … … 258 259 for (int x = - footprint; x <= footprint; x++) { 259 260 double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x]; 260 #ifdef USE_ VARIANCE261 ab /= variance->kernel[y][x];261 #ifdef USE_WEIGHT 262 ab *= weight->kernel[y][x]; 262 263 #endif 263 264 sumAB += ab; … … 295 296 double i1i2 = i1 * i2; 296 297 297 #ifdef USE_ VARIANCE298 float varVal = 1.0 / variance->kernel[y][x];299 ai2 *= varVal;300 bi2 *= varVal;301 ai1 *= varVal;302 bi1 *= varVal;303 i1i2 *= varVal;304 a *= varVal;305 b *= varVal;306 i2 *= varVal;298 #ifdef USE_WEIGHT 299 float wtVal = weight->kernel[y][x]; 300 ai2 *= wtVal; 301 bi2 *= wtVal; 302 ai1 *= wtVal; 303 bi1 *= wtVal; 304 i1i2 *= wtVal; 305 a *= wtVal; 306 b *= wtVal; 307 i2 *= wtVal; 307 308 #endif 308 309 … … 351 352 double i1i2 = i1 * i2; 352 353 353 #ifdef USE_ VARIANCE354 float varVal = 1.0 / variance->kernel[y][x];355 i1 *= varVal;356 i1i1 *= varVal;357 one *= varVal;358 i2 *= varVal;359 i1i2 *= varVal;354 #ifdef USE_WEIGHT 355 float wtVal = weight->kernel[y][x]; 356 i1 *= wtVal; 357 i1i1 *= wtVal; 358 one *= wtVal; 359 i2 *= wtVal; 360 i1i2 *= wtVal; 360 361 #endif 361 362 … … 619 620 case PM_SUBTRACTION_MODE_1: 620 621 status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image2, stamp->image1, 621 stamp-> variance, stamp->convolutions1, kernels, polyValues,622 stamp->weight, stamp->convolutions1, kernels, polyValues, 622 623 footprint); 623 624 break; 624 625 case PM_SUBTRACTION_MODE_2: 625 626 status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image1, stamp->image2, 626 stamp-> variance, stamp->convolutions2, kernels, polyValues,627 stamp->weight, stamp->convolutions2, kernels, polyValues, 627 628 footprint); 628 629 break; … … 639 640 #endif 640 641 status = calculateDualMatrixVector(stamp->matrix1, stamp->vector1, stamp->matrix2, stamp->vector2, 641 stamp->matrixX, stamp->image1, stamp->image2, stamp-> variance,642 stamp->matrixX, stamp->image1, stamp->image2, stamp->weight, 642 643 stamp->convolutions1, stamp->convolutions2, kernels, polyValues, 643 644 footprint); … … 1163 1164 1164 1165 // Calculate residuals 1165 psKernel * variance = stamp->variance; // Variancepostage stamp1166 psKernel *weight = stamp->weight; // Weight postage stamp 1166 1167 psImageInit(residual->image, 0.0); 1167 1168 if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) { … … 1247 1248 for (int y = - footprint; y <= footprint; y++) { 1248 1249 for (int x = - footprint; x <= footprint; x++) { 1249 double dev = PS_SQR(residual->kernel[y][x]) / variance->kernel[y][x];1250 double dev = PS_SQR(residual->kernel[y][x]) * weight->kernel[y][x]; 1250 1251 deviation += dev; 1251 1252 #ifdef TESTING … … 1290 1291 psFitsClose(fits); 1291 1292 } 1292 if (stamp-> variance) {1293 if (stamp->weight) { 1293 1294 psString filename = NULL; 1294 psStringAppend(&filename, "stamp_ variance_%03d.fits", i);1295 psStringAppend(&filename, "stamp_weight_%03d.fits", i); 1295 1296 psFits *fits = psFitsOpen(filename, "w"); 1296 1297 psFree(filename); 1297 psFitsWriteImage(fits, NULL, stamp-> variance->image, 0, NULL);1298 psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL); 1298 1299 psFitsClose(fits); 1299 1300 } -
trunk/psModules/src/imcombine/pmSubtractionMatch.c
r25464 r26035 68 68 float thresh2, // Threshold for stamp finding on readout 2 69 69 float stampSpacing, // Spacing between stamps 70 float sysError, // Relative systematic error in images 70 71 int size, // Kernel half-size 71 72 int footprint, // Convolution footprint for stamps … … 78 79 79 80 *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2, 80 size, footprint, stampSpacing, mode);81 size, footprint, stampSpacing, sysError, mode); 81 82 if (!*stamps) { 82 83 psError(psErrorCodeLast(), false, "Unable to find stamps."); … … 101 102 const pmReadout *ro1, const pmReadout *ro2, // Input images 102 103 int stride, // Size for convolution patches 103 float sysError, // Relative systematic error 104 float sysError, // Systematic error in images 105 float kernelError, // Systematic error in kernel 104 106 psImageMaskType maskVal, // Value to mask for input 105 107 psImageMaskType maskBad, // Mask for output bad pixels … … 151 153 PS_ASSERT_FLOAT_LESS_THAN(sysError, 1.0, false); 152 154 } 155 if (isfinite(kernelError)) { 156 PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(kernelError, 0.0, false); 157 PS_ASSERT_FLOAT_LESS_THAN(kernelError, 1.0, false); 158 } 153 159 // Don't care about maskVal 154 160 // Don't care about maskBad … … 194 200 195 201 bool pmSubtractionMatchPrecalc(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2, 196 psMetadata *analysis, int stride, float sysError,202 psMetadata *analysis, int stride, float kernelError, 197 203 psImageMaskType maskVal, psImageMaskType maskBad, psImageMaskType maskPoor, 198 204 float poorFrac, float badFrac) … … 264 270 } 265 271 266 if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, maskVal, maskBad, maskPoor,267 poorFrac, badFrac, mode)) {272 if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, kernelError, 273 maskVal, maskBad, maskPoor, poorFrac, badFrac, mode)) { 268 274 psFree(kernels); 269 275 psFree(regions); … … 300 306 301 307 if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac, 302 sysError, region, kernel, true, useFFT)) {308 kernelError, region, kernel, true, useFFT)) { 303 309 psError(PS_ERR_UNKNOWN, false, "Unable to convolve image."); 304 310 psFree(outAnalysis); … … 330 336 int inner, int ringsOrder, int binning, float penalty, 331 337 bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold, 332 int iter, float rej, float sysError, psImageMaskType maskVal, psImageMaskType maskBad, 333 psImageMaskType maskPoor, float poorFrac, float badFrac, pmSubtractionMode subMode) 334 { 335 if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, maskVal, maskBad, maskPoor, 336 poorFrac, badFrac, subMode)) { 338 int iter, float rej, float sysError, float kernelError, psImageMaskType maskVal, 339 psImageMaskType maskBad, psImageMaskType maskPoor, float poorFrac, 340 float badFrac, pmSubtractionMode subMode) 341 { 342 if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, kernelError, 343 maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) { 337 344 return false; 338 345 } … … 463 470 if (stampsName && strlen(stampsName) > 0) { 464 471 stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, size, 465 footprint, stampSpacing, s ubMode);472 footprint, stampSpacing, sysError, subMode); 466 473 } else if (sources) { 467 474 stamps = pmSubtractionStampsSetFromSources(sources, ro1->image, subMask, region, size, 468 footprint, stampSpacing, s ubMode);475 footprint, stampSpacing, sysError, subMode); 469 476 } 470 477 … … 472 479 // doesn't matter. 473 480 if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2, 474 stampSpacing, size, footprint, subMode)) {481 stampSpacing, sysError, size, footprint, subMode)) { 475 482 goto MATCH_ERROR; 476 483 } … … 538 545 539 546 if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, 540 stampThresh1, stampThresh2, stampSpacing, 547 stampThresh1, stampThresh2, stampSpacing, sysError, 541 548 size, footprint, subMode)) { 542 549 goto MATCH_ERROR; … … 608 615 psTrace("psModules.imcombine", 2, "Convolving...\n"); 609 616 if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac, 610 sysError, region, kernels, true, useFFT)) {617 kernelError, region, kernels, true, useFFT)) { 611 618 psError(PS_ERR_UNKNOWN, false, "Unable to convolve image."); 612 619 goto MATCH_ERROR; -
trunk/psModules/src/imcombine/pmSubtractionMatch.h
r25120 r26035 39 39 int iter, ///< Rejection iterations 40 40 float rej, ///< Rejection threshold 41 float sysError, ///< Relative systematic error 41 float sysError, ///< Relative systematic error in images 42 float kernelError, ///< Relative systematic error in kernel 42 43 psImageMaskType maskVal, ///< Value to mask for input 43 44 psImageMaskType maskBad, ///< Mask for output bad pixels … … 55 56 psMetadata *analysis, ///< Analysis metadata with pre-calculated kernel, region 56 57 int stride, ///< Size for convolution patches 57 float sysError, ///< Relative systematic error58 float kernelError, ///< Relative systematic error in kernel 58 59 psImageMaskType maskVal, ///< Value to mask for input 59 60 psImageMaskType maskBad, ///< Mask for output bad pixels -
trunk/psModules/src/imcombine/pmSubtractionParams.c
r21363 r26035 71 71 double *sumII, // Sum of I(x)^2/sigma(x)^2 72 72 double *sumIC, // Sum of I(x)conv(x)/sigma(x)^2 73 const pmSubtractionStamp *stamp, // Stamp with variance73 const pmSubtractionStamp *stamp, // Stamp 74 74 const psKernel *target, // Target stamp 75 75 int kernelIndex, // Index for kernel component … … 78 78 ) 79 79 { 80 psKernel * variance = stamp->variance; // Variance, sigma(x)^280 psKernel *weight = stamp->weight; // Weight image 81 81 psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest 82 82 83 83 for (int y = -footprint; y <= footprint; y++) { 84 84 psF32 *in = &target->kernel[y][-footprint]; // Dereference input 85 psF32 *wt = & variance->kernel[y][-footprint]; // Dereference variance85 psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight 86 86 psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution 87 87 for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) { 88 double temp = *in /*wt; // Temporary product88 double temp = *in * *wt; // Temporary product 89 89 *sumI += temp; 90 90 *sumII += *in * temp; … … 98 98 static void accumulateConvolutions(double *sumC, // Sum of conv(x)/sigma(x)^2 99 99 double *sumCC, // Sum of conv(x)^2/sigma(x)^2 100 const pmSubtractionStamp *stamp, // Stamp with input and variance100 const pmSubtractionStamp *stamp, // Stamp with input and weight 101 101 int kernelIndex, // Index for kernel component 102 102 int footprint, // Size of region of interest … … 104 104 ) 105 105 { 106 psKernel * variance = stamp->variance; // Variance, sigma(x)^2106 psKernel *weight = stamp->weight; // Weight image 107 107 psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest 108 108 109 109 for (int y = -footprint; y <= footprint; y++) { 110 psF32 *wt = & variance->kernel[y][-footprint]; // Dereference variance110 psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight 111 111 psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution 112 112 for (int x = -footprint; x <= footprint; x++, wt++, conv++) { 113 double convNoise = *conv /*wt; // Temporary product113 double convNoise = *conv * *wt; // Temporary product 114 114 *sumC += convNoise; 115 115 *sumCC += *conv * convNoise; … … 120 120 121 121 static double accumulateChi2(const psKernel *target, // Target stamp 122 pmSubtractionStamp *stamp, // Stamp with variance122 pmSubtractionStamp *stamp, // Stamp with weight 123 123 int kernelIndex, // Index for kernel component 124 124 double coeff, // Coefficient of convolution … … 129 129 { 130 130 double chi2 = 0.0; 131 psKernel * variance = stamp->variance; // Variance, sigma(x)^2131 psKernel *weight = stamp->weight; // Weight image 132 132 psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest 133 133 134 134 for (int y = -footprint; y <= footprint; y++) { 135 135 psF32 *in = &target->kernel[y][-footprint]; // Dereference input 136 psF32 *wt = & variance->kernel[y][-footprint]; // Dereference variance136 psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight 137 137 psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution 138 138 for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) { 139 chi2 += PS_SQR(*in - bg - coeff * *conv) /*wt;139 chi2 += PS_SQR(*in - bg - coeff * *conv) * *wt; 140 140 } 141 141 } … … 146 146 // Return the initial value of chi^2 147 147 static double initialChi2(const psKernel *target, // Target stamp 148 const pmSubtractionStamp *stamp, // Stamp with variance148 const pmSubtractionStamp *stamp, // Stamp 149 149 int footprint, // Size of convolution 150 150 pmSubtractionMode mode // Mode of subtraction 151 151 ) 152 152 { 153 psKernel * variance = stamp->variance; // Variance map153 psKernel *weight = stamp->weight; // Weight image 154 154 psKernel *source; // Source stamp 155 155 switch (mode) { … … 167 167 for (int y = -footprint; y <= footprint; y++) { 168 168 psF32 *in = &target->kernel[y][-footprint]; // Dereference input 169 psF32 *wt = & variance->kernel[y][-footprint]; // Dereference variance169 psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight 170 170 psF32 *ref = &source->kernel[y][-footprint]; // Derference reference 171 171 for (int x = -footprint; x <= footprint; x++, in++, wt++, ref++) { 172 172 float diff = *in - *ref; // Temporary value 173 chi2 += PS_SQR(diff) /*wt;173 chi2 += PS_SQR(diff) * *wt; 174 174 } 175 175 } … … 180 180 // Subtract a convolution from the input 181 181 static void subtractConvolution(psKernel *target, // Target stamp 182 const pmSubtractionStamp *stamp, // Stamp with variance182 const pmSubtractionStamp *stamp, // Stamp 183 183 int kernelIndex, // Index for kernel component 184 184 float coeff, // Coefficient of subtraction … … 288 288 289 289 // This sum is invariant to the kernel 290 psKernel * variance = stamp->variance; // Variance map for stamp290 psKernel *weight = stamp->weight; // Weight image 291 291 for (int v = -footprint; v <= footprint; v++) { 292 psF32 *wt = & variance->kernel[v][-footprint]; // Dereference variance map292 psF32 *wt = &weight->kernel[v][-footprint]; // Dereference weight 293 293 for (int u = -footprint; u <= footprint; u++, wt++) { 294 sum1 += 1.0 /*wt;294 sum1 += 1.0 * *wt; 295 295 } 296 296 } -
trunk/psModules/src/imcombine/pmSubtractionStamps.c
r25120 r26035 54 54 psFree(stamp->image1); 55 55 psFree(stamp->image2); 56 psFree(stamp-> variance);56 psFree(stamp->weight); 57 57 psFree(stamp->convolutions1); 58 58 psFree(stamp->convolutions2); … … 180 180 181 181 pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, int numRows, const psRegion *region, 182 int footprint, float spacing )182 int footprint, float spacing, float sysErr) 183 183 { 184 184 pmSubtractionStampList *list = psAlloc(sizeof(pmSubtractionStampList)); // Stamp list to return … … 221 221 list->flux = NULL; 222 222 list->footprint = footprint; 223 list->sysErr = sysErr; 223 224 224 225 return list; … … 256 257 outStamp->image1 = inStamp->image1 ? psKernelCopy(inStamp->image1) : NULL; 257 258 outStamp->image2 = inStamp->image2 ? psKernelCopy(inStamp->image2) : NULL; 258 outStamp-> variance = inStamp->variance ? psKernelCopy(inStamp->variance) : NULL;259 outStamp->weight = inStamp->weight ? psKernelCopy(inStamp->weight) : NULL; 259 260 260 261 if (inStamp->convolutions1) { … … 305 306 stamp->image1 = NULL; 306 307 stamp->image2 = NULL; 307 stamp-> variance= NULL;308 stamp->weight = NULL; 308 309 stamp->convolutions1 = NULL; 309 310 stamp->convolutions2 = NULL; … … 322 323 const psImage *image2, const psImage *subMask, 323 324 const psRegion *region, float thresh1, float thresh2, 324 int size, int footprint, float spacing, 325 int size, int footprint, float spacing, float sysErr, 325 326 pmSubtractionMode mode) 326 327 { … … 379 380 380 381 if (!stamps) { 381 stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing );382 stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, sysErr); 382 383 } 383 384 … … 452 453 psFree(stamp->image1); 453 454 psFree(stamp->image2); 454 psFree(stamp-> variance);455 psFree(stamp->weight); 455 456 psFree(stamp->convolutions1); 456 457 psFree(stamp->convolutions2); 457 stamp->image1 = stamp->image2 = stamp-> variance= NULL;458 stamp->image1 = stamp->image2 = stamp->weight = NULL; 458 459 stamp->convolutions1 = stamp->convolutions2 = NULL; 459 460 … … 487 488 const psImage *image, const psImage *subMask, 488 489 const psRegion *region, int size, int footprint, 489 float spacing, pmSubtractionMode mode)490 float spacing, float sysErr, pmSubtractionMode mode) 490 491 491 492 { … … 507 508 int numStars = x->n; // Number of stars 508 509 pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows, 509 region, footprint, spacing); // Stamp list 510 region, footprint, sysErr, 511 spacing); // Stamp list 510 512 int numStamps = stamps->num; // Number of stamps 511 513 … … 616 618 PS_ASSERT_IMAGE_TYPE(image2, PS_TYPE_F32, false); 617 619 } 618 PS_ASSERT_IMAGE_NON_NULL(variance, false); 619 PS_ASSERT_IMAGES_SIZE_EQUAL(variance, image1, false); 620 PS_ASSERT_IMAGE_TYPE(variance, PS_TYPE_F32, false); 621 PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); 620 if (variance) { 621 PS_ASSERT_IMAGE_NON_NULL(variance, false); 622 PS_ASSERT_IMAGES_SIZE_EQUAL(variance, image1, false); 623 PS_ASSERT_IMAGE_TYPE(variance, PS_TYPE_F32, false); 624 PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); 625 } 622 626 623 627 int numCols = image1->numCols, numRows = image1->numRows; // Size of images … … 646 650 assert(stamp->image1 == NULL); 647 651 assert(stamp->image2 == NULL); 648 assert(stamp-> variance== NULL);652 assert(stamp->weight == NULL); 649 653 650 654 psRegion region = psRegionSet(x - size, x + size + 1, y - size, y + size + 1); // Region of interest … … 660 664 } 661 665 662 psImage *wtSub = psImageSubset(variance, region); // Subimage with stamp 663 stamp->variance = psKernelAllocFromImage(wtSub, size, size); 664 psFree(wtSub); // Drop reference 666 psKernel *weight = stamp->weight = psKernelAlloc(-size, size, -size, size); // Weight = 1/variance 667 if (variance) { 668 psImage *varSub = psImageSubset(variance, region); // Subimage with stamp 669 psKernel *var = psKernelAllocFromImage(varSub, size, size); // Variance postage stamp 670 if (isfinite(stamps->sysErr) && stamps->sysErr > 0) { 671 float sysErr = 0.5 * stamps->sysErr; // Systematic error 672 psKernel *image1 = stamp->image1, *image2 = stamp->image2; // Input images 673 for (int y = -size; y <= size; y++) { 674 for (int x = -size; x <= size; x++) { 675 weight->kernel[y][x] = 1.0 / (var->kernel[y][x] + 676 sysErr * (image1->kernel[y][x] + image2->kernel[y][x])); 677 } 678 } 679 } else { 680 for (int y = -size; y <= size; y++) { 681 for (int x = -size; x <= size; x++) { 682 weight->kernel[y][x] = 1.0 / var->kernel[y][x]; 683 } 684 } 685 } 686 psFree(var); 687 psFree(varSub); 688 } else { 689 psImageInit(weight->image, 1.0); 690 } 665 691 666 692 stamp->status = PM_SUBTRACTION_STAMP_CALCULATE; … … 672 698 pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *image, 673 699 const psImage *subMask, const psRegion *region, 674 int size, int footprint, float spacing, 700 int size, int footprint, float spacing, float sysErr, 675 701 pmSubtractionMode mode) 676 702 { … … 702 728 703 729 pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, 704 footprint, spacing, mode); // Stamps to return730 footprint, spacing, sysErr, mode); // Stamps 705 731 psFree(x); 706 732 psFree(y); … … 716 742 pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *image, 717 743 const psImage *subMask, const psRegion *region, 718 int size, int footprint, float spacing, 744 int size, int footprint, float spacing, float sysErr, 719 745 pmSubtractionMode mode) 720 746 { … … 734 760 735 761 pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, footprint, 736 spacing, mode);762 spacing, sysErr, mode); 737 763 psFree(data); 738 764 -
trunk/psModules/src/imcombine/pmSubtractionStamps.h
r25101 r26035 24 24 psArray *flux; ///< Fluxes for possible stamps (or NULL) 25 25 int footprint; ///< Half-size of stamps 26 float sysErr; ///< Systematic error 26 27 } pmSubtractionStampList; 27 28 … … 31 32 const psRegion *region, // Region for stamps, or NULL 32 33 int footprint, // Half-size of stamps 33 float spacing // Rough average spacing between stamps 34 float spacing, // Rough average spacing between stamps 35 float sysErr // Relative systematic error or NAN 34 36 ); 35 37 … … 68 70 psKernel *image1; ///< Reference image postage stamp 69 71 psKernel *image2; ///< Input image postage stamp 70 psKernel * variance; ///< Variance imagepostage stamp, or NULL72 psKernel *weight; ///< Weight image (1/variance) postage stamp, or NULL 71 73 psArray *convolutions1; ///< Convolutions of image 1 for each kernel component, or NULL 72 74 psArray *convolutions2; ///< Convolutions of image 2 for each kernel component, or NULL … … 91 93 int footprint, ///< Half-size for stamps 92 94 float spacing, ///< Rough spacing for stamps 95 float sysErr, ///< Relative systematic error in images 93 96 pmSubtractionMode mode ///< Mode for subtraction 94 97 ); … … 103 106 int footprint, ///< Half-size for stamps 104 107 float spacing, ///< Rough spacing for stamps 108 float sysErr, ///< Systematic error in images 105 109 pmSubtractionMode mode ///< Mode for subtraction 106 110 ); … … 115 119 int footprint, ///< Half-size for stamps 116 120 float spacing, ///< Rough spacing for stamps 121 float sysErr, ///< Systematic error in images 117 122 pmSubtractionMode mode ///< Mode for subtraction 118 123 ); … … 127 132 int footprint, ///< Half-size for stamps 128 133 float spacing, ///< Rough spacing for stamps 134 float sysErr, ///< Systematic error in images 129 135 pmSubtractionMode mode ///< Mode for subtraction 130 136 );
Note:
See TracChangeset
for help on using the changeset viewer.
