Changeset 14709 for trunk/psModules/src/imcombine/pmSubtractionParams.c
- Timestamp:
- Aug 30, 2007, 11:45:26 AM (19 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/pmSubtractionParams.c
r14671 r14709 9 9 10 10 #include "pmSubtractionStamps.h" 11 #include "pmSubtraction.h" 12 #include "pmSubtractionParams.h" 11 13 12 14 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 13 15 16 #if 0 14 17 // Convolve the reference stamp by the kernel 15 18 static psKernel *convolveStamp(const pmSubtractionStamp *stamp, // Stamp to be convolved … … 45 48 return convolution; 46 49 } 47 50 #endif 48 51 49 52 // Accumulate cross-term sums for a stamp … … 53 56 const pmSubtractionStamp *stamp, // Stamp with weight 54 57 const psKernel *input, // Input image, I(x) 55 const psKernel *convolution // Convolution, conv(x) 56 ) 57 { 58 // Range of convolution 59 int xMin = convolution->xMin; 60 int xMax = convolution->xMax; 61 int yMin = convolution->yMin; 62 int yMax = convolution->yMax; 63 58 int kernelIndex, // Index for kernel component 59 int footprint // Size of region of interest 60 ) 61 { 64 62 psKernel *weight = stamp->weight; // Weight, sigma(x)^2 65 66 for (int y = yMin; y <= yMax; y++) { 67 psF32 *in = &input->kernel[y][xMin]; // Dereference input 68 psF32 *wt = &weight->kernel[y][xMin]; // Dereference weight 69 psF32 *conv = &convolution->kernel[y][xMin]; // Dereference convolution 70 for (int x = xMin; x <= xMax; x++, in++, wt++, conv++) { 63 psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest 64 65 for (int y = -footprint; y <= footprint; y++) { 66 psF32 *in = &input->kernel[y][-footprint]; // Dereference input 67 psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight 68 psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution 69 for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) { 71 70 double temp = *in / *wt; // Temporary product 72 71 *sumI += temp; … … 82 81 double *sumCC, // Sum of conv(x)^2/sigma(x)^2 83 82 const pmSubtractionStamp *stamp, // Stamp with input and weight 84 const psKernel *convolution // Convolution, conv(x) 85 ) 86 { 87 // Range of convolution 88 int xMin = convolution->xMin; 89 int xMax = convolution->xMax; 90 int yMin = convolution->yMin; 91 int yMax = convolution->yMax; 92 83 int kernelIndex, // Index for kernel component 84 int footprint // Size of region of interest 85 ) 86 { 93 87 psKernel *weight = stamp->weight; // Weight, sigma(x)^2 94 95 for (int y = yMin; y <= yMax; y++) { 96 psF32 *wt = &weight->kernel[y][xMin]; // Dereference weight 97 psF32 *conv = &convolution->kernel[y][xMin]; // Dereference convolution 98 for (int x = xMin; x <= xMax; x++, wt++, conv++) { 88 psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest 89 90 for (int y = -footprint; y <= footprint; y++) { 91 psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight 92 psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution 93 for (int x = -footprint; x <= footprint; x++, wt++, conv++) { 99 94 double convNoise = *conv / *wt; // Temporary product 100 95 *sumC += convNoise; … … 106 101 107 102 static double accumulateChi2(psKernel *input, // Input stamp 108 psKernel *convolution, // Convolution, R(x)*k(u)109 103 pmSubtractionStamp *stamp, // Stamp with weight 104 int kernelIndex, // Index for kernel component 110 105 double coeff, // Coefficient of convolution 111 double bg // Background term 112 ) 113 { 114 // Range of convolution 115 int xMin = convolution->xMin; 116 int xMax = convolution->xMax; 117 int yMin = convolution->yMin; 118 int yMax = convolution->yMax; 119 106 double bg, // Background term 107 int footprint // Size of region of interest 108 ) 109 { 120 110 double chi2 = 0.0; 121 psKernel *weight = stamp->weight; 122 123 for (int y = yMin; y <= yMax; y++) { 124 psF32 *in = &input->kernel[y][xMin]; // Dereference input 125 psF32 *wt = &weight->kernel[y][xMin]; // Dereference weight 126 psF32 *conv = &convolution->kernel[y][xMin]; // Dereference convolution 127 for (int x = xMin; x <= xMax; x++, in++, wt++, conv++) { 111 psKernel *weight = stamp->weight; // Weight, sigma(x)^2 112 psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest 113 114 for (int y = -footprint; y <= footprint; y++) { 115 psF32 *in = &input->kernel[y][-footprint]; // Dereference input 116 psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight 117 psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution 118 for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) { 128 119 chi2 += PS_SQR(*in - bg - coeff * *conv) / *wt; 129 120 } … … 159 150 static void subtractConvolution(psKernel *input, // Input stamp 160 151 const pmSubtractionStamp *stamp, // Stamp with weight 161 const psKernel *convolution, // Convolution, R(x)*k(u)152 int kernelIndex, // Index for kernel component 162 153 float coeff, // Coefficient of subtraction 163 float bg // Background term 164 ) 165 { 166 // Range of convolution 167 int xMin = convolution->xMin; 168 int xMax = convolution->xMax; 169 int yMin = convolution->yMin; 170 int yMax = convolution->yMax; 171 172 psKernel *weight = stamp->weight; // Weight map 173 174 for (int y = yMin; y <= yMax; y++) { 175 psF32 *in = &input->kernel[y][xMin]; // Dereference input 176 psF32 *conv = &convolution->kernel[y][xMin]; // Dereference convolution 177 psF32 *wt = &weight->kernel[y][xMin]; // Dereference weight 178 for (int x = xMin; x <= xMax; x++, in++, conv++, wt++) { 154 float bg, // Background term 155 int footprint // Size of region of interest 156 ) 157 { 158 psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest 159 160 for (int y = -footprint; y <= footprint; y++) { 161 psF32 *in = &input->kernel[y][-footprint]; // Dereference input 162 psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution 163 for (int x = -footprint; x <= footprint; x++, in++, conv++) { 179 164 *in -= *conv * coeff + bg; 180 165 } … … 229 214 psKernel *input = stamp->input; // Input image of interest 230 215 psImage *copy = psImageCopy(NULL, input->image, PS_TYPE_F32); // Copy of the image 231 inputs->data[i] = psKernelAllocFromImage(copy, size + footprint + 1, size + footprint + 1);216 inputs->data[i] = psKernelAllocFromImage(copy, size + footprint, size + footprint); 232 217 psFree(copy); // Drop reference 233 218 } 234 219 235 220 // Generate the convolutions 236 psArray *convolutions = psArrayAlloc(numKernels); // Convolutions per kernel component and stamp237 for (int i = 0; i < numKernels; i++) {238 psKernel *kernel = kernels->preCalc->data[i]; // Precalculated kernel239 240 psArray *stampConv = psArrayAlloc(numStamps); // Convolutions for each stamp241 convolutions->data[i] = stampConv; 242 243 for (int j = 0; j < numStamps; j++) {244 stamp Conv->data[j] = convolveStamp(stamps->data[j], kernel, footprint);221 for (int i = 0; i < numStamps; i++) { 222 pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest 223 if (!stamp->convolutions) { 224 stamp->convolutions = psArrayAlloc(numKernels); 225 } 226 227 for (int j = 0; j < numKernels; j++) { 228 psKernel *kernel = kernels->preCalc->data[j]; // Precalculated kernel 229 stamp->convolutions->data[j] = p_pmSubtractionConvolveStampPrecalc(stamp->reference, kernel); 245 230 } 246 231 } … … 264 249 psVectorInit(sumCC, 0.0); 265 250 for (int i = 0; i < numKernels; i++) { 266 psArray *stampsConv = convolutions->data[i]; // Convolutions per stamp267 251 for (int j = 0; j < numStamps; j++) { 268 accumulateConvolutions(&sumC->data.F64[i], &sumCC->data.F64[i], stamps->data[j], 269 stampsConv->data[j]); 252 accumulateConvolutions(&sumC->data.F64[i], &sumCC->data.F64[i], stamps->data[j], i, footprint); 270 253 } 271 254 } … … 298 281 double sumIC = 0.0; // sum of I(x)C(x)/sigma(x)^2 299 282 300 psArray *stampsConv = convolutions->data[i]; // Convolutions per stamp301 302 283 for (int j = 0; j < numStamps; j++) { 303 accumulateCross(&sumI, &sumII, &sumIC, stamps->data[j], inputs->data[j], stampsConv->data[j]);284 accumulateCross(&sumI, &sumII, &sumIC, stamps->data[j], inputs->data[j], i, footprint); 304 285 } 305 286 … … 310 291 double chi2 = 0.0; // Chi^2 311 292 for (int j = 0; j < numStamps; j++) { 312 chi2 += accumulateChi2(inputs->data[j], stamps Conv->data[j], stamps->data[j], coeff, bg);293 chi2 += accumulateChi2(inputs->data[j], stamps->data[j], i, coeff, bg, footprint); 313 294 } 314 295 … … 329 310 ranking->data.S32[bestIndex] = iter; 330 311 // Remove its contribution, and don't include it in the future. 331 psArray *stampsConv = convolutions->data[bestIndex]; // Convolutions per stamp332 312 for (int j = 0; j < numStamps; j++) { 333 subtractConvolution(inputs->data[j], stamps->data[j], stampsConv->data[j], 334 bestCoeff, bestBG); 313 subtractConvolution(inputs->data[j], stamps->data[j], bestIndex, bestCoeff, bestBG, footprint); 335 314 } 336 315 … … 349 328 } 350 329 psFree(inputs); 351 psFree(convolutions);352 330 psFree(sumC); 353 331 psFree(sumCC); … … 392 370 psKernel *subtract = kernels->preCalc->data[0]; // Kernel to subtract from the rest 393 371 for (int i = 1; i < newSize; i++) { 394 psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest 395 psBinaryOp(kernel->image, kernel->image, "-", subtract->image); 372 if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) { 373 psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest 374 psBinaryOp(kernel->image, kernel->image, "-", subtract->image); 375 } 396 376 } 397 377 } else if (type == PM_SUBTRACTION_KERNEL_GUNK) { … … 400 380 401 381 for (int i = 0; i < newSize; i++) { 402 psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest 403 kernel->kernel[0][0] -= 1.0; 382 if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) { 383 psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest 384 kernel->kernel[0][0] -= 1.0; 385 } 404 386 } 405 387
Note:
See TracChangeset
for help on using the changeset viewer.
