Changeset 15756 for trunk/psModules/src/imcombine/pmSubtractionKernels.c
- Timestamp:
- Dec 6, 2007, 3:57:15 PM (18 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/pmSubtractionKernels.c
r15432 r15756 4 4 5 5 #include <stdio.h> 6 #include <string.h> 6 7 #include <strings.h> 7 8 #include <pslib.h> 8 9 10 #include "pmSubtraction.h" 9 11 #include "pmSubtractionKernels.h" 10 12 … … 22 24 psFree(kernels->vStop); 23 25 psFree(kernels->preCalc); 24 } 25 26 // Raise a number to an integer power 26 psFree(kernels->solution1); 27 psFree(kernels->solution2); 28 } 29 30 // Raise an integer to an integer power 27 31 static inline long power(int value, // Value 28 32 int exp // Exponent … … 45 49 bool p_pmSubtractionKernelsAddGrid(pmSubtractionKernels *kernels, int start, int size) 46 50 { 47 PS_ASSERT_PTR_NON_NULL(kernels, false); 48 PS_ASSERT_VECTOR_NON_NULL(kernels->u, false); 49 PS_ASSERT_VECTOR_NON_NULL(kernels->v, false); 50 PS_ASSERT_VECTOR_NON_NULL(kernels->widths, false); 51 PS_ASSERT_VECTOR_TYPE(kernels->u, PS_TYPE_S32, false); 52 PS_ASSERT_VECTOR_TYPE(kernels->v, PS_TYPE_S32, false); 53 PS_ASSERT_VECTOR_TYPE(kernels->widths, PS_TYPE_F32, false); 54 PS_ASSERT_ARRAY_NON_NULL(kernels->preCalc, false); 51 PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); 55 52 PS_ASSERT_INT_NONNEGATIVE(start, false); 56 53 PS_ASSERT_INT_NONNEGATIVE(size, false); 57 54 58 int numNew = PS_SQR(2 * size + 1) ; // Number of new kernel parameters to add55 int numNew = PS_SQR(2 * size + 1) - 1; // Number of new kernel parameters to add 59 56 60 57 // Ensure the sizes match … … 68 65 for (int v = - size, index = start; v <= size; v++) { 69 66 for (int u = - size; u <= size; u++, index++) { 67 if (v == 0 && u == 0) { 68 // Skip normalisation component: added explicitly 69 index--; 70 continue; 71 } 70 72 kernels->widths->data.F32[index] = NAN; 71 73 kernels->u->data.S32[index] = u; … … 81 83 82 84 pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder, 83 const psVector *fwhms, const psVector *orders) 85 const psVector *fwhms, const psVector *orders, 86 pmSubtractionMode mode) 84 87 { 85 88 PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL); … … 97 100 for (int i = 0; i < numGaussians; i++) { 98 101 int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian 99 psStringAppend(¶ms, "(%. 2f,%d)", fwhms->data.F32[i], orders->data.S32[i]);102 psStringAppend(¶ms, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]); 100 103 num += (gaussOrder + 1) * (gaussOrder + 2) / 2; 101 104 } 102 105 103 106 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, 104 size, spatialOrder ); // The kernels107 size, spatialOrder, mode); // The kernels 105 108 psStringAppend(&kernels->description, "ISIS(%d,%s,%d)", size, params, spatialOrder); 106 109 … … 149 152 } 150 153 151 kernels->subIndex = -1;152 153 154 return kernels; 154 155 } … … 159 160 160 161 pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type, 161 int size, int spatialOrder )162 int size, int spatialOrder, pmSubtractionMode mode) 162 163 { 163 164 pmSubtractionKernels *kernels = psAlloc(sizeof(pmSubtractionKernels)); // Kernels, to return … … 173 174 kernels->uStop = NULL; 174 175 kernels->vStop = NULL; 175 kernels->subIndex = 0;176 176 kernels->size = size; 177 177 kernels->inner = 0; 178 178 kernels->spatialOrder = spatialOrder; 179 179 kernels->bgOrder = 0; 180 181 return kernels; 182 } 183 184 pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder) 180 kernels->mode = mode; 181 kernels->solution1 = NULL; 182 kernels->solution2 = NULL; 183 184 return kernels; 185 } 186 187 pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, pmSubtractionMode mode) 185 188 { 186 189 PS_ASSERT_INT_POSITIVE(size, NULL); 187 190 PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); 188 191 189 int num = PS_SQR(2 * size + 1) ; // Number of basis functions192 int num = PS_SQR(2 * size + 1) - 1; // Number of basis functions 190 193 191 194 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS, 192 size, spatialOrder ); // The kernels195 size, spatialOrder, mode); // The kernels 193 196 psStringAppend(&kernels->description, "POIS(%d,%d)", size, spatialOrder); 194 197 psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements", … … 199 202 } 200 203 201 kernels->subIndex = num / 2;202 assert(kernels->u->data.S32[kernels->subIndex] == 0 &&203 kernels->v->data.S32[kernels->subIndex] == 0);204 205 204 return kernels; 206 205 } … … 208 207 209 208 pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder, 210 const psVector *fwhms, const psVector *orders) 209 const psVector *fwhms, const psVector *orders, 210 pmSubtractionMode mode) 211 211 { 212 212 pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, 213 fwhms, orders ); // Kernels213 fwhms, orders, mode); // Kernels 214 214 if (!kernels) { 215 215 return NULL; 216 216 } 217 217 218 // Subtract a reference kernelfrom all the others to maintain flux scaling218 // Subtract normalisation from all the others to maintain flux scaling 219 219 if (spatialOrder > 0) { 220 int subIndex = 0; // Index of kernel to subtract from others221 psKernel *subtract = kernels->preCalc->data[subIndex]; // Kernel to subtract from others222 220 int numGaussians = fwhms->n; // Number of Gaussians 223 221 for (int i = 0, index = 0; i < numGaussians; i++) { 224 222 for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { 225 223 for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { 226 if (uOrder % 2 == 0 && vOrder % 2 == 0 && index != subIndex) {224 if (uOrder % 2 == 0 && vOrder % 2 == 0) { 227 225 psKernel *kernel = kernels->preCalc->data[index]; // Kernel to subtract 228 psBinaryOp(kernel->image, kernel->image, "-", subtract->image);226 kernel->kernel[0][0] -= 1.0; 229 227 } 230 228 } … … 242 240 psFitsWriteImage(kernelFile, NULL, kernel->image, 0, NULL); 243 241 psFitsClose(kernelFile); 244 } 245 } 246 247 return kernels; 248 } 249 250 /// Generate SPAM kernels 251 pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, ///< Half-size of the kernel 252 int spatialOrder, ///< Order of spatial variations 253 int inner, ///< Inner radius to preserve unbinned 254 int binning ///< Kernel binning factor 255 ) 242 double sum = 0.0; 243 for (int y = 0; y < kernel->image->numRows; y++) { 244 for (int x = 0; x < kernel->image->numCols; x++) { 245 sum += kernel->image->data.F32[y][x]; 246 } 247 } 248 psTrace("psModules.imcombine.kernel", 10, "Kernel %d sum: %lf\n", i, sum); 249 } 250 } 251 252 return kernels; 253 } 254 255 pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, int spatialOrder, int inner, int binning, 256 pmSubtractionMode mode) 256 257 { 257 258 PS_ASSERT_INT_POSITIVE(size, NULL); … … 269 270 psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); 270 271 271 int num = PS_SQR(2 * numTotal + 1) ; // Number of basis functions272 int num = PS_SQR(2 * numTotal + 1) - 1; // Number of basis functions 272 273 273 274 psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); 274 275 275 276 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, 276 size, spatialOrder ); // The kernels277 size, spatialOrder, mode); // The kernels 277 278 psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d)", size, inner, binning, spatialOrder); 278 279 … … 280 281 size, inner, binning, spatialOrder, num); 281 282 282 kernels->uStop = psVectorAlloc(num, PS_TYPE_ F32);283 kernels->vStop = psVectorAlloc(num, PS_TYPE_ F32);283 kernels->uStop = psVectorAlloc(num, PS_TYPE_S32); 284 kernels->vStop = psVectorAlloc(num, PS_TYPE_S32); 284 285 285 286 psVector *locations = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Locations for each kernel element … … 314 315 315 316 for (int j = - numTotal; j <= numTotal; j++, index++) { 317 if (i == 0 && j == 0) { 318 // Skip normalisation component: added explicitly 319 index--; 320 continue; 321 } 316 322 int v = locations->data.S32[numTotal + j]; // Location of pixel 317 323 int vStop = v + widths->data.S32[numTotal + j]; // Width of pixel … … 327 333 } 328 334 329 kernels->subIndex = num / 2;330 assert(kernels->u->data.S32[kernels->subIndex] == 0 &&331 kernels->v->data.S32[kernels->subIndex] == 0 &&332 kernels->uStop->data.S32[kernels->subIndex] == 0 &&333 kernels->vStop->data.S32[kernels->subIndex] == 0);334 335 335 psFree(locations); 336 336 psFree(widths); … … 340 340 341 341 342 /// Generate FRIES kernels 343 pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, ///< Half-size of the kernel 344 int spatialOrder, ///< Order of spatial variations 345 int inner ///< Inner radius to preserve unbinned 346 ) 342 pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, pmSubtractionMode mode) 347 343 { 348 344 PS_ASSERT_INT_POSITIVE(size, NULL); … … 366 362 psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); 367 363 368 int num = PS_SQR(2 * numTotal + 1) ; // Number of basis functions364 int num = PS_SQR(2 * numTotal + 1) - 1; // Number of basis functions 369 365 370 366 psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); 371 367 372 368 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, 373 size, spatialOrder ); // The kernels369 size, spatialOrder, mode); // The kernels 374 370 psStringAppend(&kernels->description, "FRIES(%d,%d,%d)", size, inner, spatialOrder); 375 371 … … 377 373 size, inner, spatialOrder, num); 378 374 379 kernels->uStop = psVectorAlloc(num, PS_TYPE_ F32);380 kernels->vStop = psVectorAlloc(num, PS_TYPE_ F32);375 kernels->uStop = psVectorAlloc(num, PS_TYPE_S32); 376 kernels->vStop = psVectorAlloc(num, PS_TYPE_S32); 381 377 382 378 psVector *start = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); … … 409 405 int uStop = stop->data.S32[numTotal + i]; // Width of pixel 410 406 for (int j = - numTotal; j <= numTotal; j++, index++) { 407 if (i == 0 && j == 0) { 408 // Skip normalisation component: added explicitly 409 index--; 410 continue; 411 } 411 412 int v = start->data.S32[numTotal + j]; // Location of pixel 412 413 int vStop = stop->data.S32[numTotal + j]; // Width of pixel … … 422 423 } 423 424 424 kernels->subIndex = num / 2;425 assert(kernels->u->data.S32[kernels->subIndex] == 0 &&426 kernels->v->data.S32[kernels->subIndex] == 0 &&427 kernels->uStop->data.S32[kernels->subIndex] == 0 &&428 kernels->vStop->data.S32[kernels->subIndex] == 0);429 430 425 psFree(start); 431 426 psFree(stop); … … 436 431 // Grid United with Normal Kernel 437 432 pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms, 438 const psVector *orders, int inner )433 const psVector *orders, int inner, pmSubtractionMode mode) 439 434 { 440 435 PS_ASSERT_INT_POSITIVE(size, NULL); … … 449 444 450 445 pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, 451 fwhms, orders ); // Kernels446 fwhms, orders, mode); // Kernels 452 447 psStringPrepend(&kernels->description, "GUNK="); 453 448 psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder); … … 469 464 470 465 int numISIS = kernels->num; // Number of ISIS kernels 471 int numInner = PS_SQR(2 * inner + 1); // Number of inner kernel elements472 466 473 467 if (!p_pmSubtractionKernelsAddGrid(kernels, numISIS, inner)) { … … 475 469 } 476 470 477 kernels->subIndex = numInner/2 + numISIS;478 assert(kernels->u->data.S32[kernels->subIndex] == 0 &&479 kernels->v->data.S32[kernels->subIndex] == 0);480 481 471 return kernels; 482 472 } 483 473 484 474 // RINGS --- just what it says 485 pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder) 475 pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder, 476 pmSubtractionMode mode) 486 477 { 487 478 PS_ASSERT_INT_POSITIVE(size, NULL); … … 505 496 } 506 497 507 int numInner = inner ;// Number of pixels in the inner region498 int numInner = inner - 1; // Number of pixels in the inner region 508 499 int numOuter = fibNum; // Number of summed pixels in the outer region 509 500 510 501 int numRings = numOuter + numInner; // Number of rings (not including the central pixel) 511 int numPoly = (ringsOrder + 1) * (ringsOrder + 2) / 2; // Number of polynomial variants of each ring512 513 int num = numRings * numPoly + 1; // Total number of basis functions502 int numPoly = PM_SUBTRACTION_POLYTERMS(ringsOrder); // Number of polynomial variants of each ring 503 504 int num = numRings * numPoly; // Total number of basis functions 514 505 515 506 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, 516 size, spatialOrder ); // The kernels507 size, spatialOrder, mode); // The kernels 517 508 psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d)", size, inner, ringsOrder, spatialOrder); 518 509 … … 522 513 // Set the Gaussian kernel parameters 523 514 int fibIndex = 1, fibIndexMinus1 = 0; // Fibonnacci parameters 524 int radiusLast = 0; // Last radius 525 for (int i = 0, index = 0; i < numRings + 1; i++) { 526 515 int radiusLast = 1; // Last radius 516 for (int i = 1, index = 0; i < numRings + 1; i++) { 527 517 float lower2; // Lower limit of radius^2 528 518 float upper2; // Upper limit of radius^2 529 if (i == 0) { 530 lower2 = 0; 531 upper2 = PS_SQR(0.5); 532 } else if (i <= numInner) { 519 if (i <= inner) { 533 520 // A ring every pixel width 534 521 float radius = i; … … 572 559 for (int v = -size; v <= size; v++) { 573 560 int v2 = PS_SQR(v); // Square of v 574 // float vPoly = power(v, vOrder); // Value of v^vOrder575 561 float vPoly = powf(v/(float)size, vOrder); // Value of v^vOrder 576 562 … … 579 565 int distance2 = u2 + v2; // Distance from the centre 580 566 if (distance2 > lower2 && distance2 < upper2) { 581 // float uPoly = power(u, uOrder); // Value of u^uOrder582 567 float uPoly = powf(u/(float)size, uOrder); // Value of u^uOrder 583 568 … … 627 612 } 628 613 629 kernels->subIndex = 0;630 631 614 return kernels; 632 615 } … … 634 617 pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder, 635 618 const psVector *fwhms, const psVector *orders, int inner, 636 int binning, int ringsOrder )619 int binning, int ringsOrder, pmSubtractionMode mode) 637 620 { 638 621 switch (type) { 639 622 case PM_SUBTRACTION_KERNEL_POIS: 640 return pmSubtractionKernelsPOIS(size, spatialOrder );623 return pmSubtractionKernelsPOIS(size, spatialOrder, mode); 641 624 case PM_SUBTRACTION_KERNEL_ISIS: 642 return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders );625 return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, mode); 643 626 case PM_SUBTRACTION_KERNEL_SPAM: 644 return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning );627 return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, mode); 645 628 case PM_SUBTRACTION_KERNEL_FRIES: 646 return pmSubtractionKernelsFRIES(size, spatialOrder, inner );629 return pmSubtractionKernelsFRIES(size, spatialOrder, inner, mode); 647 630 case PM_SUBTRACTION_KERNEL_GUNK: 648 return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner );631 return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, mode); 649 632 case PM_SUBTRACTION_KERNEL_RINGS: 650 return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder );633 return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, mode); 651 634 default: 652 635 psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type);
Note:
See TracChangeset
for help on using the changeset viewer.
