Changeset 18287 for trunk/psModules/src/imcombine/pmSubtractionKernels.c
- Timestamp:
- Jun 23, 2008, 12:41:08 PM (18 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/pmSubtractionKernels.c
r18146 r18287 24 24 psFree(kernels->vStop); 25 25 psFree(kernels->preCalc); 26 psFree(kernels->penalties); 26 27 psFree(kernels->solution1); 27 28 psFree(kernels->solution2); … … 60 61 kernels->v = psVectorRealloc(kernels->v, start + numNew); 61 62 kernels->preCalc = psArrayRealloc(kernels->preCalc, start + numNew); 63 kernels->penalties = psVectorRealloc(kernels->penalties, start + numNew); 62 64 kernels->inner = start; 63 65 … … 74 76 kernels->v->data.S32[index] = v; 75 77 kernels->preCalc->data[index] = NULL; 78 kernels->penalties->data.F32[index] = kernels->penalty * (PS_SQR(u) + PS_SQR(v)); 76 79 77 80 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v); … … 84 87 pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder, 85 88 const psVector *fwhms, const psVector *orders, 86 pmSubtractionMode mode)89 float penalty, pmSubtractionMode mode) 87 90 { 88 91 PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL); … … 104 107 } 105 108 106 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, 107 s ize, spatialOrder, mode); // The kernels108 psStringAppend(&kernels->description, "ISIS(%d,%s,%d )", size, params, spatialOrder);109 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, size, 110 spatialOrder, penalty, mode); // The kernels 111 psStringAppend(&kernels->description, "ISIS(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty); 109 112 110 113 psLogMsg("psModules.imcombine", PS_LOG_INFO, "ISIS kernel: %s,%d --> %d elements", … … 122 125 psKernel *preCalc = psKernelAlloc(-size, size, -size, size); 123 126 double sum = 0.0; // Normalisation 127 double moment = 0.0; // Moment, for penalty 124 128 for (int v = -size; v <= size; v++) { 125 129 for (int u = -size; u <= size; u++) { 126 130 sum += preCalc->kernel[v][u] = norm * power(u, uOrder) * power(v, vOrder) * 127 131 expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigma)); 132 moment += preCalc->kernel[v][u] * (PS_SQR(u) + PS_SQR(v)); 128 133 } 129 134 } … … 146 151 } 147 152 kernels->preCalc->data[index] = preCalc; 148 149 psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, 150 fwhms->data.F32[i], uOrder, vOrder); 153 kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment); 154 155 psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, 156 fwhms->data.F32[i], uOrder, vOrder, fabsf(moment)); 151 157 } 152 158 } … … 161 167 162 168 pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type, 163 int size, int spatialOrder, pmSubtractionMode mode) 169 int size, int spatialOrder, float penalty, 170 pmSubtractionMode mode) 164 171 { 165 172 pmSubtractionKernels *kernels = psAlloc(sizeof(pmSubtractionKernels)); // Kernels, to return … … 173 180 kernels->widths = psVectorAlloc(numBasisFunctions, PS_TYPE_F32); 174 181 kernels->preCalc = psArrayAlloc(numBasisFunctions); 182 kernels->penalty = penalty; 183 kernels->penalties = psVectorAlloc(numBasisFunctions, PS_TYPE_F32); 175 184 kernels->uStop = NULL; 176 185 kernels->vStop = NULL; … … 188 197 } 189 198 190 pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, pmSubtractionMode mode) 199 pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty, 200 pmSubtractionMode mode) 191 201 { 192 202 PS_ASSERT_INT_POSITIVE(size, NULL); … … 195 205 int num = PS_SQR(2 * size + 1) - 1; // Number of basis functions 196 206 197 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS, 198 s ize, spatialOrder, mode); // The kernels199 psStringAppend(&kernels->description, "POIS(%d,%d )", size, spatialOrder);207 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS, size, 208 spatialOrder, penalty, mode); // The kernels 209 psStringAppend(&kernels->description, "POIS(%d,%d,%.2e)", size, spatialOrder, penalty); 200 210 psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements", 201 211 size, spatialOrder, num); … … 211 221 pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder, 212 222 const psVector *fwhms, const psVector *orders, 213 pmSubtractionMode mode)214 { 215 pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, 216 fwhms, orders, mode); // Kernels223 float penalty, pmSubtractionMode mode) 224 { 225 pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders, 226 penalty, mode); // Kernels 217 227 if (!kernels) { 218 228 return NULL; … … 242 252 243 253 pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, int spatialOrder, int inner, int binning, 244 pmSubtractionMode mode)254 float penalty, pmSubtractionMode mode) 245 255 { 246 256 PS_ASSERT_INT_POSITIVE(size, NULL); … … 262 272 psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); 263 273 264 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, 265 size, spatialOrder, mode); // The kernels 266 psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d)", size, inner, binning, spatialOrder); 274 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, size, 275 spatialOrder, penalty, mode); // The kernels 276 psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d,%.2e)", size, inner, binning, spatialOrder, 277 penalty); 267 278 268 279 psLogMsg("psModules.imcombine", PS_LOG_INFO, "SPAM kernel: %d,%d,%d,%d --> %d elements", … … 324 335 psFree(widths); 325 336 337 psWarning("Kernel penalty for dual-convolution is not configured for SPAM kernels."); 338 psVectorInit(kernels->penalties, 0.0); 339 326 340 return kernels; 327 341 } 328 342 329 343 330 pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, pmSubtractionMode mode) 344 pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, float penalty, 345 pmSubtractionMode mode) 331 346 { 332 347 PS_ASSERT_INT_POSITIVE(size, NULL); … … 354 369 psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); 355 370 356 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, 357 s ize, spatialOrder, mode); // The kernels358 psStringAppend(&kernels->description, "FRIES(%d,%d,%d )", size, inner, spatialOrder);371 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, size, 372 spatialOrder, penalty, mode); // The kernels 373 psStringAppend(&kernels->description, "FRIES(%d,%d,%d,%.2e)", size, inner, spatialOrder, penalty); 359 374 360 375 psLogMsg("psModules.imcombine", PS_LOG_INFO, "FRIES kernel: %d,%d,%d --> %d elements", … … 414 429 psFree(stop); 415 430 431 psWarning("Kernel penalty for dual-convolution is not configured for FRIES kernels."); 432 psVectorInit(kernels->penalties, 0.0); 433 416 434 return kernels; 417 435 } … … 419 437 // Grid United with Normal Kernel 420 438 pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms, 421 const psVector *orders, int inner, pmSubtractionMode mode) 439 const psVector *orders, int inner, float penalty, 440 pmSubtractionMode mode) 422 441 { 423 442 PS_ASSERT_INT_POSITIVE(size, NULL); … … 431 450 PS_ASSERT_INT_LESS_THAN(inner, size, NULL); 432 451 433 pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, 434 fwhms, orders, mode); // Kernels 452 // XXX GUNK doesn't seem to work --- doesn't add the POIS components, or at least, they're not noticed 453 454 pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders, 455 penalty, mode); // Kernels 435 456 psStringPrepend(&kernels->description, "GUNK="); 436 457 psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder); … … 447 468 // RINGS --- just what it says 448 469 pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder, 449 pmSubtractionMode mode)470 float penalty, pmSubtractionMode mode) 450 471 { 451 472 PS_ASSERT_INT_POSITIVE(size, NULL); … … 477 498 int num = numRings * numPoly; // Total number of basis functions 478 499 479 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, 480 size, spatialOrder, mode); // The kernels 481 psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d)", size, inner, ringsOrder, spatialOrder); 500 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, size, 501 spatialOrder, penalty, mode); // The kernels 502 psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d,%.2e)", size, inner, ringsOrder, spatialOrder, 503 penalty); 482 504 483 505 psLogMsg("psModules.imcombine", PS_LOG_INFO, "RINGS kernel: %d,%d,%d,%d --> %d elements", … … 520 542 psVector *vCoords = data->data[1] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // v coords 521 543 psVector *poly = data->data[2] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_F32); // Polynomial 544 double moment = 0.0; // Moment, for penalty 522 545 523 546 if (i == 0) { … … 527 550 uCoords->n = vCoords->n = poly->n = 1; 528 551 radiusLast = 0; 552 moment = 0.0; 529 553 } else { 530 554 int j = 0; // Index for data … … 546 570 poly->data.F32[j] = polyVal; 547 571 norm += polyVal; 572 moment += polyVal * (PS_SQR(u) + PS_SQR(v)); 548 573 549 574 psVectorExtend(uCoords, RINGS_BUFFER, 1); … … 571 596 psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32)); 572 597 } 598 // moment /= norm; 573 599 } 574 600 … … 578 604 kernels->u->data.S32[index] = uOrder; 579 605 kernels->v->data.S32[index] = vOrder; 606 kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment); 580 607 581 608 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d\n", index, … … 590 617 pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder, 591 618 const psVector *fwhms, const psVector *orders, int inner, 592 int binning, int ringsOrder, pmSubtractionMode mode) 619 int binning, int ringsOrder, float penalty, 620 pmSubtractionMode mode) 593 621 { 594 622 switch (type) { 595 623 case PM_SUBTRACTION_KERNEL_POIS: 596 return pmSubtractionKernelsPOIS(size, spatialOrder, mode);624 return pmSubtractionKernelsPOIS(size, spatialOrder, penalty, mode); 597 625 case PM_SUBTRACTION_KERNEL_ISIS: 598 return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, mode);626 return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, mode); 599 627 case PM_SUBTRACTION_KERNEL_SPAM: 600 return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, mode);628 return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, mode); 601 629 case PM_SUBTRACTION_KERNEL_FRIES: 602 return pmSubtractionKernelsFRIES(size, spatialOrder, inner, mode);630 return pmSubtractionKernelsFRIES(size, spatialOrder, inner, penalty, mode); 603 631 case PM_SUBTRACTION_KERNEL_GUNK: 604 return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, mode);632 return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, penalty, mode); 605 633 case PM_SUBTRACTION_KERNEL_RINGS: 606 return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, mode);634 return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, penalty, mode); 607 635 default: 608 636 psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type); … … 657 685 int binning = 0; // Binning to use 658 686 int ringsOrder = 0; // Polynomial order for rings 687 float penalty = 0.0; // Penalty for wideness 659 688 660 689 if (strncmp(description, "ISIS", 4) == 0) { … … 684 713 685 714 ptr++; // Eat ',' 686 spatialOrder = parseStringInt(ptr); 715 PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt); 716 penalty = parseStringFloat(ptr); 687 717 } 688 718 } else if (strncmp(description, "RINGS", 5) == 0) { … … 692 722 PARSE_STRING_NUMBER(inner, ptr, ',', parseStringInt); 693 723 PARSE_STRING_NUMBER(ringsOrder, ptr, ',', parseStringInt); 694 PARSE_STRING_NUMBER(spatialOrder, ptr, ')', parseStringInt); 724 PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt); 725 PARSE_STRING_NUMBER(penalty, ptr, ')', parseStringInt); 695 726 } else { 696 727 psAbort("Deciphering kernels other than ISIS and RINGS is not currently supported."); … … 699 730 700 731 return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, 701 inner, binning, ringsOrder, mode);732 inner, binning, ringsOrder, penalty, mode); 702 733 } 703 734
Note:
See TracChangeset
for help on using the changeset viewer.
