IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 14360


Ignore:
Timestamp:
Jul 20, 2007, 6:30:57 PM (19 years ago)
Author:
Paul Price
Message:

Adding RINGS kernel --- a series of rings, modified by polynomials.

Location:
trunk/psModules/src/imcombine
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/imcombine/pmSubtraction.c

    r14340 r14360  
    44 *  @author GLG, MHPCC
    55 *
    6  *  @version $Revision: 1.31 $ $Name: not supported by cvs2svn $
    7  *  @date $Date: 2007-07-20 21:18:50 $
     6 *  @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
     7 *  @date $Date: 2007-07-21 04:30:57 $
    88 *
    99 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    … …  
    170170              break; \
    171171          } \
     172          case PM_SUBTRACTION_KERNEL_RINGS: { \
     173              if (i == (KERNELS)->subIndex) { \
     174                  (TARGET)->kernel[0][0] += value; \
     175                  break; \
     176              } \
     177              psArray *preCalc = (KERNELS)->preCalc->data[i]; /* Precalculated data */ \
     178              psVector *uCoords = preCalc->data[0]; /* u coordinates */ \
     179              psVector *vCoords = preCalc->data[1]; /* v coordinates */ \
     180              psVector *poly = preCalc->data[2]; /* Polynomial values */ \
     181              int num = uCoords->n;     /* Number of pixels */ \
     182              for (int j = 0; j < num; j++) { \
     183                  int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; /* Kernel coordinates */ \
     184                  (TARGET)->kernel[v][u] += value * FUNC(poly->data.F32[j]); \
     185              } \
     186              /* The (0,0) kernel is subtracted from other kernels to preserve photometric scaling */ \
     187              if (kernels->spatialOrder > 0) { \
     188                  kernel->kernel[0][0] += subValue; \
     189              } \
     190              break; \
     191          } \
    172192          default: \
    173193            psAbort("Should never get here."); \
    … …  
    274294              break;
    275295          }
     296          case PM_SUBTRACTION_KERNEL_RINGS: {
     297              if (i == kernels->subIndex) {
     298                  kernel->kernel[0][0] += value;
     299                  break;
     300              }
     301              psArray *preCalc = kernels->preCalc->data[i]; // Precalculated data
     302              psVector *uCoords = preCalc->data[0]; // u coordinates
     303              psVector *vCoords = preCalc->data[1]; // v coordinates
     304              psVector *poly = preCalc->data[2]; // Polynomial values
     305              int num = uCoords->n;     // Number of pixels
     306
     307              for (int j = 0; j < num; j++) {
     308                  int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates
     309                  kernel->kernel[v][u] += value * weightFunc(poly->data.F32[j]);
     310              }
     311
     312              // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling
     313              if (kernels->spatialOrder > 0) {
     314                  kernel->kernel[0][0] += subValue;
     315              }
     316              break;
     317          }
    276318          default:
    277319            psAbort("Should never get here.");
    … …  
    367409          }
    368410          return polyValue * sum - sub;
     411      }
     412      case PM_SUBTRACTION_KERNEL_RINGS: {
     413          if (index == kernels->subIndex) {
     414              return image->data.F32[0][0];
     415          }
     416          psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data
     417          psVector *uCoords = preCalc->data[0]; // u coordinates
     418          psVector *vCoords = preCalc->data[1]; // v coordinates
     419          psVector *poly = preCalc->data[2]; // Polynomial values
     420          int num = uCoords->n;         // Number of pixels
     421          double sum = 0.0;             // Accumulated sum from convolution
     422          for (int j = 0; j < num; j++) {
     423              int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates
     424              sum += image->data.F32[y + v][x + u] * poly->data.F32[j];
     425          }
     426          // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling
     427          return polyValue * sum - image->data.F32[0][0];
    369428      }
    370429      default:

  • trunk/psModules/src/imcombine/pmSubtractionKernels.c

    r14341 r14360  
    88#include "pmSubtractionKernels.h"
    99
     10#define RINGS_BUFFER 10                 // Buffer size for RINGS data
     11
    1012
    1113// Free function for pmSubtractionKernels
    … …  
    1416    psFree(kernels->u);
    1517    psFree(kernels->v);
    16     psFree(kernels->sigma);
     18    psFree(kernels->widths);
    1719    psFree(kernels->uStop);
    1820    psFree(kernels->vStop);
    … …  
    5052    kernels->u = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);
    5153    kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);
    52     kernels->sigma = NULL;
     54    kernels->widths = NULL;
    5355    kernels->uStop = NULL;
    5456    kernels->vStop = NULL;
    … …  
    134136    psFree(params);
    135137
    136     kernels->sigma = psVectorAlloc(num, PS_TYPE_F32);
     138    kernels->widths = psVectorAlloc(num, PS_TYPE_F32);
    137139    kernels->preCalc = psArrayAlloc(num);
    138140
    … …  
    156158                for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) {
    157159                    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];
    159161                        kernels->u->data.S32[index] = uOrder;
    160162                        kernels->v->data.S32[index] = vOrder;
    … …  
    435437    psFree(params);
    436438
    437     kernels->sigma = psVectorAlloc(numGaussianVars * numSpatial, PS_TYPE_F32);
     439    kernels->widths = psVectorAlloc(numGaussianVars * numSpatial, PS_TYPE_F32);
    438440    kernels->preCalc = psArrayAlloc(numGaussianVars * numSpatial);
    439441    kernels->inner = numGaussianVars * numSpatial;
    … …  
    459461                for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) {
    460462                    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];
    462464                        kernels->u->data.S32[index] = uOrder;
    463465                        kernels->v->data.S32[index] = vOrder;
    … …  
    505507}
    506508
     509// RINGS --- just what it says
     510pmSubtractionKernels *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
    507632pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder,
    508633                                                   const psVector *sigmas, const psVector *orders, int inner,
    509                                                    int binning)
     634                                                   int binning, int ringsOrder)
    510635{
    511636    switch (type) {
    … …  
    520645      case PM_SUBTRACTION_KERNEL_GUNK:
    521646        return pmSubtractionKernelsGUNK(size, spatialOrder, sigmas, orders, inner);
     647      case PM_SUBTRACTION_KERNEL_RINGS:
     648        return pmSubtractionKernelsRINGS(size, ringsOrder, spatialOrder);
    522649      default:
    523650        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type);

  • trunk/psModules/src/imcombine/pmSubtractionKernels.h

    r14305 r14360  
    1111    PM_SUBTRACTION_KERNEL_FRIES,        ///< Fibonacci Radius Increases Excellence of Subtraction
    1212    PM_SUBTRACTION_KERNEL_GUNK,         ///< Grid United with Normal Kernel --- POIS and ISIS hybrid
     13    PM_SUBTRACTION_KERNEL_RINGS,        ///< Rings Instead of the Normal Gaussian Subtraction
    1314} pmSubtractionKernelsType;
    1415
    … …  
    1718    pmSubtractionKernelsType type;      ///< Type of kernels --- allowing the use of multiple kernels
    1819    psVector *u, *v;                    ///< Offset (for POIS) or polynomial order (for ISIS)
    19     psVector *sigma;                    ///< Gaussian widths (ISIS only)
     20    psVector *widths;                   ///< Gaussian widths (ISIS) or ring radius (RINGS)
    2021    psVector *uStop, *vStop;            ///< Width of kernel element (SPAM,FRIES only)
    2122    psVector *xOrder, *yOrder;          ///< Spatial Polynomial order (for all)
    … …  
    7071    );
    7172
     73/// Generate RINGS kernels
     74pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, ///< Half-size of the kernel
     75                                                int spatialOrder, ///< Order of spatial variations
     76                                                int ringsOrder ///< Polynomial order
     77    );
     78
     79
    7280/// Generate a kernel of a specified type
    7381pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, ///< Kernel type
    … …  
    7785                                                   const psVector *orders, ///< Polynomial order of gaussians
    7886                                                   int inner, ///< Inner radius to preserve unbinned
    79                                                    int binning ///< Kernel binning factor
     87                                                   int binning, ///< Kernel binning factor
     88                                                   int ringsOrder ///< Polynomial order for RINGS
    8089    );
    8190
Note: See TracChangeset for help on using the changeset viewer.