- Timestamp:
- Nov 8, 2007, 12:59:04 PM (19 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionParams.c
r15247 r15516 50 50 #endif 51 51 52 /// Select the appropriate convolution, given the kernel basis function and subtraction mode 53 static inline psKernel *selectConvolution(const pmSubtractionStamp *stamp, // Stamp 54 int kernelIndex, // Index for kernel component 55 pmSubtractionMode mode // Mode of subtraction 56 ) 57 { 58 switch (mode) { 59 case PM_SUBTRACTION_MODE_1: 60 return stamp->convolutions1->data[kernelIndex]; 61 case PM_SUBTRACTION_MODE_2: 62 return stamp->convolutions2->data[kernelIndex]; 63 default: 64 psAbort("Unsupported subtraction mode: %x", mode); 65 } 66 return NULL; // Unreached 67 } 68 52 69 // Accumulate cross-term sums for a stamp 53 70 static void accumulateCross(double *sumI, // Sum of I(x)/sigma(x)^2 … … 55 72 double *sumIC, // Sum of I(x)conv(x)/sigma(x)^2 56 73 const pmSubtractionStamp *stamp, // Stamp with weight 57 const psKernel * input, // Input image, I(x)74 const psKernel *target, // Target stamp 58 75 int kernelIndex, // Index for kernel component 59 int footprint // Size of region of interest 76 int footprint, // Size of region of interest 77 pmSubtractionMode mode // Mode of subtraction 60 78 ) 61 79 { 62 80 psKernel *weight = stamp->weight; // Weight, sigma(x)^2 63 psKernel *convolution = s tamp->convolutions->data[kernelIndex]; // Convolution of interest81 psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest 64 82 65 83 for (int y = -footprint; y <= footprint; y++) { 66 psF32 *in = & input->kernel[y][-footprint]; // Dereference input84 psF32 *in = &target->kernel[y][-footprint]; // Dereference input 67 85 psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight 68 86 psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution … … 82 100 const pmSubtractionStamp *stamp, // Stamp with input and weight 83 101 int kernelIndex, // Index for kernel component 84 int footprint // Size of region of interest 102 int footprint, // Size of region of interest 103 pmSubtractionMode mode // Mode of subtraction 85 104 ) 86 105 { 87 106 psKernel *weight = stamp->weight; // Weight, sigma(x)^2 88 psKernel *convolution = s tamp->convolutions->data[kernelIndex]; // Convolution of interest107 psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest 89 108 90 109 for (int y = -footprint; y <= footprint; y++) { … … 100 119 } 101 120 102 static double accumulateChi2( psKernel *input, // Input stamp121 static double accumulateChi2(const psKernel *target, // Target stamp 103 122 pmSubtractionStamp *stamp, // Stamp with weight 104 123 int kernelIndex, // Index for kernel component 105 124 double coeff, // Coefficient of convolution 106 125 double bg, // Background term 107 int footprint // Size of region of interest 126 int footprint, // Size of region of interest 127 pmSubtractionMode mode // Mode of subtraction 108 128 ) 109 129 { 110 130 double chi2 = 0.0; 111 131 psKernel *weight = stamp->weight; // Weight, sigma(x)^2 112 psKernel *convolution = s tamp->convolutions->data[kernelIndex]; // Convolution of interest132 psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest 113 133 114 134 for (int y = -footprint; y <= footprint; y++) { 115 psF32 *in = & input->kernel[y][-footprint]; // Dereference input135 psF32 *in = &target->kernel[y][-footprint]; // Dereference input 116 136 psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight 117 137 psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution … … 125 145 126 146 // Return the initial value of chi^2 127 static double initialChi2( psKernel *input, // Input stamp147 static double initialChi2(const psKernel *target, // Target stamp 128 148 const pmSubtractionStamp *stamp, // Stamp with weight 129 int footprint // Size of convolution 149 int footprint, // Size of convolution 150 pmSubtractionMode mode // Mode of subtraction 130 151 ) 131 152 { 132 153 psKernel *weight = stamp->weight; // Weight map 133 psKernel *reference = stamp->reference; // Reference stamp 154 psKernel *source; // Source stamp 155 switch (mode) { 156 case PM_SUBTRACTION_MODE_1: 157 source = stamp->image1; 158 break; 159 case PM_SUBTRACTION_MODE_2: 160 source = stamp->image2; 161 break; 162 default: 163 psAbort("Unsupported subtraction mode: %x", mode); 164 } 134 165 135 166 double chi2 = 0.0; // Chi^2 136 167 for (int y = -footprint; y <= footprint; y++) { 137 psF32 *in = & input->kernel[y][-footprint]; // Dereference input168 psF32 *in = &target->kernel[y][-footprint]; // Dereference input 138 169 psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight 139 psF32 *ref = & reference->kernel[y][-footprint]; // Derference reference170 psF32 *ref = &source->kernel[y][-footprint]; // Derference reference 140 171 for (int x = -footprint; x <= footprint; x++, in++, wt++, ref++) { 141 172 float diff = *in - *ref; // Temporary value … … 148 179 149 180 // Subtract a convolution from the input 150 static void subtractConvolution(psKernel * input, // Input stamp181 static void subtractConvolution(psKernel *target, // Target stamp 151 182 const pmSubtractionStamp *stamp, // Stamp with weight 152 183 int kernelIndex, // Index for kernel component 153 184 float coeff, // Coefficient of subtraction 154 185 float bg, // Background term 155 int footprint // Size of region of interest156 )157 { 158 psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest 159 186 int footprint, // Size of region of interest 187 pmSubtractionMode mode // Mode of subtraction 188 ) 189 { 190 psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest 160 191 for (int y = -footprint; y <= footprint; y++) { 161 psF32 *in = & input->kernel[y][-footprint]; // Dereference input192 psF32 *in = &target->kernel[y][-footprint]; // Dereference input 162 193 psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution 163 194 for (int x = -footprint; x <= footprint; x++, in++, conv++) { … … 173 204 int spatialOrder, const psVector *fwhms, int maxOrder, 174 205 const pmSubtractionStampList *stamps, int footprint, 175 float tolerance )206 float tolerance, pmSubtractionMode mode) 176 207 { 177 208 if (type != PM_SUBTRACTION_KERNEL_ISIS && type != PM_SUBTRACTION_KERNEL_GUNK) { … … 210 241 // Need to save the stamp inputs --- we're changing the values! 211 242 int numStamps = stamps->num; // Number of stamps 212 psArray * inputs = psArrayAlloc(numStamps); // Deep copies of the inputs243 psArray *targets = psArrayAlloc(numStamps); // Deep copies of the targets 213 244 psVector *badStamps = psVectorAlloc(numStamps, PS_TYPE_U8); // Mark the bad stamps 214 245 psVectorInit(badStamps, 0); … … 219 250 continue; 220 251 } 221 psKernel *input = stamp->input; // Input image of interest 222 psImage *copy = psImageCopy(NULL, input->image, PS_TYPE_F32); // Copy of the image 223 inputs->data[i] = psKernelAllocFromImage(copy, size + footprint, size + footprint); 252 psKernel *target; // Target image of interest 253 switch (mode) { 254 case PM_SUBTRACTION_MODE_1: 255 target = stamp->image2; 256 break; 257 case PM_SUBTRACTION_MODE_2: 258 target = stamp->image1; 259 break; 260 default: 261 psAbort("Unsupported subtraction mode: %x", mode); 262 } 263 psImage *copy = psImageCopy(NULL, target->image, PS_TYPE_F32); // Copy of the image 264 targets->data[i] = psKernelAllocFromImage(copy, size + footprint, size + footprint); 224 265 psFree(copy); // Drop reference 225 266 } … … 238 279 continue; 239 280 } 240 if (!pmSubtractionConvolveStamp(stamp, kernels, footprint )) {281 if (!pmSubtractionConvolveStamp(stamp, kernels, footprint, mode)) { 241 282 psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i); 242 psFree( inputs);283 psFree(targets); 243 284 psFree(kernels); 244 285 psFree(badStamps); … … 258 299 "Sum of 1/sigma^2 is non-finite for stamp %d (%d,%d)\n", 259 300 i, (int)stamp->x, (int)stamp->y); 260 psFree( inputs);301 psFree(targets); 261 302 psFree(kernels); 262 303 psFree(badStamps); … … 265 306 266 307 for (int j = 0; j < numKernels; j++) { 267 accumulateConvolutions(&sumC->data.F64[j], &sumCC->data.F64[j], stamp, j, footprint );268 } 269 270 lastChi2 += initialChi2( inputs->data[i], stamp, footprint);308 accumulateConvolutions(&sumC->data.F64[j], &sumCC->data.F64[j], stamp, j, footprint, mode); 309 } 310 311 lastChi2 += initialChi2(targets->data[i], stamp, footprint, mode); 271 312 numPixels += PS_SQR(2 * footprint + 1); 272 313 } … … 297 338 } 298 339 pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest 299 accumulateCross(&sumI, &sumII, &sumIC, stamp, inputs->data[j], i, footprint);340 accumulateCross(&sumI, &sumII, &sumIC, stamp, targets->data[j], i, footprint, mode); 300 341 } 301 342 … … 310 351 } 311 352 pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest 312 chi2 += accumulateChi2( inputs->data[j], stamp, i, coeff, bg, footprint);353 chi2 += accumulateChi2(targets->data[j], stamp, i, coeff, bg, footprint, mode); 313 354 } 314 355 … … 328 369 if (bestIndex == -1) { 329 370 psError(PS_ERR_UNKNOWN, false, "Unable to find best kernel component in round %d.", iter); 330 psFree( inputs);371 psFree(targets); 331 372 psFree(sumC); 332 373 psFree(sumCC); … … 345 386 } 346 387 pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest 347 subtractConvolution( inputs->data[j], stamp, bestIndex, bestCoeff, bestBG, footprint);388 subtractConvolution(targets->data[j], stamp, bestIndex, bestCoeff, bestBG, footprint, mode); 348 389 } 349 390 … … 361 402 lastChi2 = bestChi2; 362 403 } 363 psFree( inputs);404 psFree(targets); 364 405 psFree(sumC); 365 406 psFree(sumCC); … … 401 442 pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest 402 443 psArray *convolutions = convNew->data[j]; // Convolutions for this stamp 403 convolutions->data[rank] = psMemIncrRefCounter(s tamp->convolutions->data[i]);444 convolutions->data[rank] = psMemIncrRefCounter(selectConvolution(stamp, i, mode)); 404 445 } 405 446 } … … 420 461 } 421 462 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest 422 psFree(stamp->convolutions); 423 stamp->convolutions = convNew->data[i]; 463 psFree(stamp->convolutions1); 464 psFree(stamp->convolutions2); 465 switch (mode) { 466 case PM_SUBTRACTION_MODE_1: 467 stamp->convolutions1 = convNew->data[i]; 468 stamp->convolutions2 = NULL; 469 break; 470 case PM_SUBTRACTION_MODE_2: 471 stamp->convolutions1 = NULL; 472 stamp->convolutions2 = convNew->data[i]; 473 break; 474 default: 475 psAbort("Unsupported subtraction mode: %x", mode); 476 } 424 477 } 425 478
Note:
See TracChangeset
for help on using the changeset viewer.
