Changeset 21363 for trunk/psModules/src/imcombine/pmSubtractionEquation.c
- Timestamp:
- Feb 5, 2009, 4:31:25 PM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/pmSubtractionEquation.c
r20561 r21363 24 24 static inline double calculateSumProduct(const psKernel *image1, // First image in multiplication 25 25 const psKernel *image2, // Second image in multiplication 26 const psKernel * weight, // Weightimage26 const psKernel *variance, // Variance image 27 27 int footprint // (Half-)Size of stamp 28 28 ) … … 31 31 for (int y = - footprint; y <= footprint; y++) { 32 32 for (int x = - footprint; x <= footprint; x++) { 33 sum += image1->kernel[y][x] * image2->kernel[y][x] / 1.0; // weight->kernel[y][x];33 sum += image1->kernel[y][x] * image2->kernel[y][x] / 1.0; // variance->kernel[y][x]; 34 34 } 35 35 } … … 42 42 const psKernel *image1, // First image in multiplication 43 43 const psKernel *image2, // Second image in multiplication 44 const psKernel * weight, // Weightimage44 const psKernel *variance, // Variance image 45 45 const psImage *polyValues, // Spatial polynomial values 46 46 int numKernels, // Number of kernel basis functions … … 50 50 ) 51 51 { 52 double sum = calculateSumProduct(image1, image2, weight, footprint); // Sum of the image products52 double sum = calculateSumProduct(image1, image2, variance, footprint); // Sum of the image products 53 53 if (!isfinite(sum)) { 54 54 return false; … … 79 79 const psKernel *image1, // First image in multiplication 80 80 const psKernel *image2, // Second image in multiplication 81 const psKernel * weight, // Weightimage81 const psKernel *variance, // Variance image 82 82 const psImage *polyValues, // Spatial polynomial values 83 83 int numKernels, // Number of kernel basis functions … … 87 87 ) 88 88 { 89 double sum = calculateSumProduct(image1, image2, weight, footprint); // Sum of the image products89 double sum = calculateSumProduct(image1, image2, variance, footprint); // Sum of the image products 90 90 if (!isfinite(sum)) { 91 91 return false; … … 120 120 const psArray *convolutions1, // Convolutions for element 1 121 121 const psArray *convolutions2, // Convolutions for element 2 122 const psKernel * weight, // Weightimage122 const psKernel *variance, // Variance image 123 123 const psImage *polyValues, // Polynomial values 124 124 int numKernels, // Number of kernel basis functions … … 135 135 psKernel *jConv = convolutions2->data[j]; // Convolution for j-th element 136 136 137 if (!calculateMatrixElement2(matrix, i, j, iConv, jConv, weight, polyValues, numKernels,137 if (!calculateMatrixElement2(matrix, i, j, iConv, jConv, variance, polyValues, numKernels, 138 138 footprint, spatialOrder, symmetric)) { 139 139 psTrace("psModules.imcombine", 2, "Bad sumCC at %d, %d", i, j); … … 151 151 const psArray *convolutions, // Convolutions of source with kernels 152 152 const psKernel *input, // Input stamp, or NULL 153 const psKernel * weight, // Weightstamp153 const psKernel *variance, // Variance stamp 154 154 const psImage *polyValues, // Spatial polynomial values 155 155 int footprint, // (Half-)Size of stamp … … 171 171 172 172 // Square part of the matrix (convolution-convolution products) 173 if (!calculateMatrixSquare(matrix, convolutions, convolutions, weight, polyValues, numKernels,173 if (!calculateMatrixSquare(matrix, convolutions, convolutions, variance, polyValues, numKernels, 174 174 spatialOrder, footprint)) { 175 175 return false; … … 185 185 186 186 // Normalisation-convolution terms 187 if (!calculateMatrixElement1(matrix, i, normIndex, conv, input, weight, polyValues, numKernels,187 if (!calculateMatrixElement1(matrix, i, normIndex, conv, input, variance, polyValues, numKernels, 188 188 footprint, spatialOrder, true)) { 189 189 psTrace("psModules.imcombine", 2, "Bad sumIC at %d", i); … … 195 195 for (int y = - footprint; y <= footprint; y++) { 196 196 for (int x = - footprint; x <= footprint; x++) { 197 sumC += conv->kernel[y][x] / 1.0; // weight->kernel[y][x];197 sumC += conv->kernel[y][x] / 1.0; // variance->kernel[y][x]; 198 198 } 199 199 } … … 218 218 for (int y = - footprint; y <= footprint; y++) { 219 219 for (int x = - footprint; x <= footprint; x++) { 220 double invNoise2 = 1.0 / 1.0; // weight->kernel[y][x];220 double invNoise2 = 1.0 / 1.0; // variance->kernel[y][x]; 221 221 double value = input->kernel[y][x] * invNoise2; 222 222 sumI += value; … … 253 253 const psKernel *input, // Input stamp, or NULL if !normAndBG 254 254 const psKernel *target, // Target stamp 255 const psKernel * weight, // Weightstamp255 const psKernel *variance, // Variance stamp 256 256 const psImage *polyValues, // Spatial polynomial values 257 257 int footprint, // (Half-)Size of stamp … … 277 277 for (int y = - footprint; y <= footprint; y++) { 278 278 for (int x = - footprint; x <= footprint; x++) { 279 sumTC += target->kernel[y][x] * conv->kernel[y][x] / 1.0; // weight->kernel[y][x];279 sumTC += target->kernel[y][x] * conv->kernel[y][x] / 1.0; // variance->kernel[y][x]; 280 280 } 281 281 } … … 297 297 for (int y = - footprint; y <= footprint; y++) { 298 298 for (int x = - footprint; x <= footprint; x++) { 299 float value = target->kernel[y][x] / 1.0; // weight->kernel[y][x];299 float value = target->kernel[y][x] / 1.0; // variance->kernel[y][x]; 300 300 sumIT += value * input->kernel[y][x]; 301 301 sumT += value; … … 329 329 const psArray *convolutions2, // Convolutions of image 2 330 330 const psKernel *image1, // Image 1 stamp 331 const psKernel * weight, // Weightstamp331 const psKernel *variance, // Variance stamp 332 332 const psImage *polyValues, // Spatial polynomial values 333 333 int footprint // (Half-)Size of stamp … … 348 348 PM_SUBTRACTION_INDICES(normIndex, bgIndex, kernels); 349 349 350 if (!calculateMatrixSquare(matrix, convolutions1, convolutions2, weight, polyValues, numKernels,350 if (!calculateMatrixSquare(matrix, convolutions1, convolutions2, variance, polyValues, numKernels, 351 351 spatialOrder, footprint)) { 352 352 return false; … … 356 356 // Normalisation 357 357 psKernel *conv = convolutions2->data[i]; // Convolution 358 if (!calculateMatrixElement1(matrix, i, normIndex, conv, image1, weight, polyValues, numKernels,358 if (!calculateMatrixElement1(matrix, i, normIndex, conv, image1, variance, polyValues, numKernels, 359 359 footprint, spatialOrder, false)) { 360 360 psTrace("psModules.imcombine", 2, "Bad sumIC at %d", i); … … 366 366 for (int y = - footprint; y <= footprint; y++) { 367 367 for (int x = - footprint; x <= footprint; x++) { 368 sumC += conv->kernel[y][x] / 1.0; // weight->kernel[y][x];368 sumC += conv->kernel[y][x] / 1.0; // variance->kernel[y][x]; 369 369 } 370 370 } … … 559 559 case PM_SUBTRACTION_MODE_1: 560 560 status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, 561 stamp-> weight, polyValues, footprint, true);561 stamp->variance, polyValues, footprint, true); 562 562 status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1, 563 stamp->image2, stamp-> weight, polyValues, footprint, true);563 stamp->image2, stamp->variance, polyValues, footprint, true); 564 564 break; 565 565 case PM_SUBTRACTION_MODE_2: 566 566 status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions2, stamp->image2, 567 stamp-> weight, polyValues, footprint, true);567 stamp->variance, polyValues, footprint, true); 568 568 status &= calculateVector(stamp->vector1, kernels, stamp->convolutions2, stamp->image2, 569 stamp->image1, stamp-> weight, polyValues, footprint, true);569 stamp->image1, stamp->variance, polyValues, footprint, true); 570 570 break; 571 571 case PM_SUBTRACTION_MODE_DUAL: … … 581 581 #endif 582 582 status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, 583 stamp-> weight, polyValues, footprint, true);583 stamp->variance, polyValues, footprint, true); 584 584 status &= calculateMatrix(stamp->matrix2, kernels, stamp->convolutions2, NULL, 585 stamp-> weight, polyValues, footprint, false);585 stamp->variance, polyValues, footprint, false); 586 586 status &= calculateMatrixCross(stamp->matrixX, kernels, stamp->convolutions1, 587 stamp->convolutions2, stamp->image1, stamp-> weight, polyValues,587 stamp->convolutions2, stamp->image1, stamp->variance, polyValues, 588 588 footprint); 589 589 status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1, 590 stamp->image2, stamp-> weight, polyValues, footprint, true);590 stamp->image2, stamp->variance, polyValues, footprint, true); 591 591 status &= calculateVector(stamp->vector2, kernels, stamp->convolutions2, NULL, 592 stamp->image2, stamp-> weight, polyValues, footprint, false);592 stamp->image2, stamp->variance, polyValues, footprint, false); 593 593 break; 594 594 default: … … 1033 1033 1034 1034 // Calculate residuals 1035 psKernel * weight = stamp->weight; // Weightpostage stamp1035 psKernel *variance = stamp->variance; // Variance postage stamp 1036 1036 psImageInit(residual->image, 0.0); 1037 1037 if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) { … … 1098 1098 for (int y = - footprint; y <= footprint; y++) { 1099 1099 for (int x = - footprint; x <= footprint; x++) { 1100 double dev = PS_SQR(residual->kernel[y][x]) / weight->kernel[y][x];1100 double dev = PS_SQR(residual->kernel[y][x]) / variance->kernel[y][x]; 1101 1101 deviation += dev; 1102 1102 #ifdef TESTING … … 1141 1141 psFitsClose(fits); 1142 1142 } 1143 if (stamp-> weight) {1143 if (stamp->variance) { 1144 1144 psString filename = NULL; 1145 psStringAppend(&filename, "stamp_ weight_%03d.fits", i);1145 psStringAppend(&filename, "stamp_variance_%03d.fits", i); 1146 1146 psFits *fits = psFitsOpen(filename, "w"); 1147 1147 psFree(filename); 1148 psFitsWriteImage(fits, NULL, stamp-> weight->image, 0, NULL);1148 psFitsWriteImage(fits, NULL, stamp->variance->image, 0, NULL); 1149 1149 psFitsClose(fits); 1150 1150 }
Note:
See TracChangeset
for help on using the changeset viewer.
