Changeset 14360 for trunk/psModules/src/imcombine/pmSubtractionKernels.c
- Timestamp:
- Jul 20, 2007, 6:30:57 PM (19 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/pmSubtractionKernels.c
r14341 r14360 8 8 #include "pmSubtractionKernels.h" 9 9 10 #define RINGS_BUFFER 10 // Buffer size for RINGS data 11 10 12 11 13 // Free function for pmSubtractionKernels … … 14 16 psFree(kernels->u); 15 17 psFree(kernels->v); 16 psFree(kernels-> sigma);18 psFree(kernels->widths); 17 19 psFree(kernels->uStop); 18 20 psFree(kernels->vStop); … … 50 52 kernels->u = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); 51 53 kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); 52 kernels-> sigma= NULL;54 kernels->widths = NULL; 53 55 kernels->uStop = NULL; 54 56 kernels->vStop = NULL; … … 134 136 psFree(params); 135 137 136 kernels-> sigma= psVectorAlloc(num, PS_TYPE_F32);138 kernels->widths = psVectorAlloc(num, PS_TYPE_F32); 137 139 kernels->preCalc = psArrayAlloc(num); 138 140 … … 156 158 for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { 157 159 for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { 158 kernels-> sigma->data.F32[index] = sigmas->data.F32[i];160 kernels->widths->data.F32[index] = sigmas->data.F32[i]; 159 161 kernels->u->data.S32[index] = uOrder; 160 162 kernels->v->data.S32[index] = vOrder; … … 435 437 psFree(params); 436 438 437 kernels-> sigma= psVectorAlloc(numGaussianVars * numSpatial, PS_TYPE_F32);439 kernels->widths = psVectorAlloc(numGaussianVars * numSpatial, PS_TYPE_F32); 438 440 kernels->preCalc = psArrayAlloc(numGaussianVars * numSpatial); 439 441 kernels->inner = numGaussianVars * numSpatial; … … 459 461 for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { 460 462 for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { 461 kernels-> sigma->data.F32[index] = sigmas->data.F32[i];463 kernels->widths->data.F32[index] = sigmas->data.F32[i]; 462 464 kernels->u->data.S32[index] = uOrder; 463 465 kernels->v->data.S32[index] = vOrder; … … 505 507 } 506 508 509 // RINGS --- just what it says 510 pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int ringsOrder) 511 { 512 PS_ASSERT_INT_POSITIVE(size, NULL); 513 PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); 514 PS_ASSERT_INT_NONNEGATIVE(ringsOrder, NULL); 515 516 int numRings = size + 1; // Number of rings; 0 --> size, inclusive 517 int numPoly = (ringsOrder + 1) * (ringsOrder + 2) / 2; // Number of polynomial variants of each ring 518 int numSpatial = (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of spatial variations of a kernel 519 520 int num = numRings * numPoly * numSpatial; // Total number of basis functions 521 522 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_GUNK, 523 size, spatialOrder); // The kernels 524 525 psLogMsg("psModules.imcombine", PS_LOG_INFO, "RINGS kernel: %d,%d,%d --> %d elements", 526 size, ringsOrder, spatialOrder, num); 527 528 kernels->widths = psVectorAlloc(num, PS_TYPE_S32); 529 kernels->preCalc = psArrayAlloc(num); 530 531 // Set the Gaussian kernel parameters 532 for (int i = 0, index = 0; i < size; i++) { 533 // Iterate over (u,v) order 534 for (int uOrder = 0; uOrder <= (i == 0 ? 0 : ringsOrder); uOrder++) { 535 for (int vOrder = 0; vOrder <= (i == 0 ? 0 : ringsOrder) - uOrder; vOrder++) { 536 537 psArray *data = psArrayAlloc(3); // Container for data 538 psVector *uCoords = data->data[0] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // u coords 539 psVector *vCoords = data->data[1] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // v coords 540 psVector *poly = data->data[2] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_F32); // Polynomial 541 if (i == 0) { 542 uCoords->data.S32[0] = vCoords->data.S32[0] = 0; 543 poly->data.F32[0] = 0; 544 uCoords->n = vCoords->n = poly->n = 1; 545 } else { 546 float radius = i; // Radius of ring 547 float lower2 = PS_SQR(radius - 0.5); // Lower limit of radius^2 548 float upper2 = PS_SQR(radius + 0.5); // Upper limit of radius^2 549 550 int j = 0; // Index for data 551 for (int v = 1; v <= size; v++) { 552 int v2 = PS_SQR(v); // Square of v 553 float vPolyPlus = power(v, vOrder); // Value of (+v)^vOrder 554 float vPolyMinus = power(-v, vOrder); // Value of (-v)^vOrder 555 556 // u = 0 557 uCoords->data.S32[j] = 0; 558 vCoords->data.S32[j] = v; 559 poly->data.F32[j] = (uOrder == 0 ? vPolyPlus : 0.0); 560 j++; 561 562 uCoords->data.S32[j] = 0; 563 vCoords->data.S32[j] = -v; 564 poly->data.F32[j] = (uOrder == 0 ? vPolyMinus : 0.0); 565 j++; 566 567 psVectorExtend(uCoords, RINGS_BUFFER, 2); 568 psVectorExtend(vCoords, RINGS_BUFFER, 2); 569 psVectorExtend(poly, RINGS_BUFFER, 2); 570 571 for (int u = 1; u <= v; u++) { 572 int u2 = PS_SQR(u); // Square of u 573 int distance2 = u2 + v2; // Distance from the centre 574 if (distance2 > lower2 && distance2 < upper2) { 575 float uPolyPlus = power(u, uOrder); // Value of (+u)^uOrder 576 float uPolyMinus = power(-u, uOrder); // Value of (-u)^uOrder 577 578 uCoords->data.S32[j] = u; 579 vCoords->data.S32[j] = v; 580 poly->data.F32[j] = uPolyPlus * vPolyPlus; 581 j++; 582 583 uCoords->data.S32[j] = u; 584 vCoords->data.S32[j] = -v; 585 poly->data.F32[j] = uPolyPlus * vPolyMinus; 586 j++; 587 588 uCoords->data.S32[j] = -u; 589 vCoords->data.S32[j] = v; 590 poly->data.F32[j] = uPolyMinus * vPolyPlus; 591 j++; 592 593 uCoords->data.S32[j] = -u; 594 vCoords->data.S32[j] = -v; 595 poly->data.F32[j] = uPolyMinus * vPolyMinus; 596 j++; 597 598 psVectorExtend(uCoords, RINGS_BUFFER, 4); 599 psVectorExtend(vCoords, RINGS_BUFFER, 4); 600 psVectorExtend(poly, RINGS_BUFFER, 4); 601 } 602 } 603 } 604 } 605 606 // Iterate over spatial order. This loop creates the terms for 607 // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. 608 for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { 609 for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { 610 kernels->preCalc->data[index] = data; 611 kernels->widths->data.S32[index] = PS_SQR(i); 612 kernels->u->data.S32[index] = uOrder; 613 kernels->v->data.S32[index] = vOrder; 614 kernels->xOrder->data.S32[index] = xOrder; 615 kernels->yOrder->data.S32[index] = yOrder; 616 617 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d %d\n", index, 618 i, uOrder, vOrder, xOrder, yOrder); 619 } 620 } 621 } 622 } 623 } 624 625 kernels->subIndex = 0; 626 assert(kernels->xOrder->data.S32[kernels->subIndex] == 0 && 627 kernels->yOrder->data.S32[kernels->subIndex] == 0); 628 629 return kernels; 630 } 631 507 632 pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder, 508 633 const psVector *sigmas, const psVector *orders, int inner, 509 int binning )634 int binning, int ringsOrder) 510 635 { 511 636 switch (type) { … … 520 645 case PM_SUBTRACTION_KERNEL_GUNK: 521 646 return pmSubtractionKernelsGUNK(size, spatialOrder, sigmas, orders, inner); 647 case PM_SUBTRACTION_KERNEL_RINGS: 648 return pmSubtractionKernelsRINGS(size, ringsOrder, spatialOrder); 522 649 default: 523 650 psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type);
Note:
See TracChangeset
for help on using the changeset viewer.
