IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 27, 2007, 10:55:23 AM (19 years ago)
Author:
Paul Price
Message:

Adding function pmSubtractionKernelsOptimumISIS, which derives an optimum set of kernels for ISIS (and GUNK). Its use necessitated a few API changes to the other parts of the image subtraction.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r14542 r14671  
    4040
    4141//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    42 // Public functions
     42// Semi-public functions
    4343//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    4444
    45 pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type,
    46                                                 int size, int spatialOrder)
    47 {
    48     pmSubtractionKernels *kernels = psAlloc(sizeof(pmSubtractionKernels)); // Kernels, to return
    49     psMemSetDeallocator(kernels, (psFreeFunc)subtractionKernelsFree);
    50 
    51     kernels->type = type;
    52     kernels->description = NULL;
    53     kernels->num = numBasisFunctions;
    54     kernels->u = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);
    55     kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);
    56     kernels->widths = NULL;
    57     kernels->uStop = NULL;
    58     kernels->vStop = NULL;
    59     kernels->subIndex = 0;
    60     kernels->preCalc = NULL;
    61     kernels->size = size;
    62     kernels->inner = 0;
    63     kernels->spatialOrder = spatialOrder;
    64     kernels->bgOrder = 0;
    65 
    66     return kernels;
    67 }
    68 
    69 pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder)
    70 {
    71     PS_ASSERT_INT_POSITIVE(size, NULL);
    72     PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
    73 
    74     int num = PS_SQR(2 * size + 1); // Number of basis functions
    75 
    76     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS,
    77                                                               size, spatialOrder); // The kernels
    78     psStringAppend(&kernels->description, "POIS(%d,%d)", size, spatialOrder);
    79 
    80     psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements",
    81              size, spatialOrder, num);
     45bool p_pmSubtractionKernelsAddGrid(pmSubtractionKernels *kernels, int start, int size)
     46{
     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);
     55    PS_ASSERT_INT_NONNEGATIVE(start, false);
     56    PS_ASSERT_INT_NONNEGATIVE(size, false);
     57
     58    int numNew = PS_SQR(2 * size + 1);  // Number of new kernel parameters to add
     59
     60    // Ensure the sizes match
     61    kernels->widths = psVectorRealloc(kernels->widths, start + numNew);
     62    kernels->u = psVectorRealloc(kernels->u, start + numNew);
     63    kernels->v = psVectorRealloc(kernels->v, start + numNew);
     64    kernels->preCalc = psArrayRealloc(kernels->preCalc, start + numNew);
    8265
    8366    // Generate a set of kernels for each (u,v)
    84     for (int v = - size, index = 0; v <= size; v++) {
     67    for (int v = - size, index = start; v <= size; v++) {
    8568        for (int u = - size; u <= size; u++, index++) {
     69            kernels->widths->data.F32[index] = NAN;
    8670            kernels->u->data.S32[index] = u;
    8771            kernels->v->data.S32[index] = v;
     72            kernels->preCalc->data[index] = NULL;
    8873
    8974            psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v);
    … …  
    9176    }
    9277
    93     kernels->subIndex = num / 2;
    94     assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
    95            kernels->v->data.S32[kernels->subIndex] == 0);
    96 
    97     return kernels;
    98 }
    99 
    100 pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder,
    101                                                const psVector *fwhms, const psVector *orders)
     78    return true;
     79}
     80
     81pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder,
     82                                                  const psVector *fwhms, const psVector *orders)
    10283{
    10384    PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL);
    … …  
    127108    psFree(params);
    128109
    129     kernels->widths = psVectorAlloc(num, PS_TYPE_F32);
    130     kernels->preCalc = psArrayAlloc(num);
    131     psKernel *subtract = NULL;          // Kernel to subtract to maintain flux scaling
    132 
    133110    // Set the kernel parameters
    134     for (int i = 0, index = 0; i < numGaussians; i++) {
    135         float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma
    136         float norm = 1.0 / (M_2_PI * sqrtf(sigma)); // Normalisation for Gaussian
    137 
    138         // Iterate over (u,v) order
    139         for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
    140             for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    141 
    142                 // Set the pre-calculated kernel
    143                 psKernel *preCalc = psKernelAlloc(-size, size, -size, size);
    144                 double sum = 0.0;       // Normalisation
    145                 for (int v = -size; v <= size; v++) {
    146                     for (int u = -size; u <= size; u++) {
    147                         sum += preCalc->kernel[v][u] = norm * power(u, uOrder) * power(v, vOrder) *
    148                             expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigma));
    149                     }
    150                 }
    151                 if (index == 0) {
    152                     subtract = preCalc;
    153                     for (int v = -size; v <= size; v++) {
    154                         for (int u = -size; u <= size; u++) {
    155                             preCalc->kernel[v][u] /= sum;
    156                         }
    157                     }
    158                 } else if (uOrder % 2 == 0 && vOrder % 2 == 0) {
    159                     // Normalise sum of kernel component to unity for even functions
    160                     for (int v = -size; v <= size; v++) {
    161                         for (int u = -size; u <= size; u++) {
    162                             preCalc->kernel[v][u] = preCalc->kernel[v][u] / sum - subtract->kernel[v][u];
    163                         }
    164                     }
    165                 }
    166 
    167                 kernels->widths->data.F32[index] = fwhms->data.F32[i];
    168                 kernels->u->data.S32[index] = uOrder;
    169                 kernels->v->data.S32[index] = vOrder;
    170                 kernels->preCalc->data[index] = preCalc;
    171 
    172                 psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index,
    173                         fwhms->data.F32[i], uOrder, vOrder);
    174             }
    175         }
    176     }
    177 
    178     kernels->subIndex = -1;
    179 
    180     if (psTraceGetLevel("psModules.imcombine.kernel") >= 10) {
    181         for (int i = 0; i < num; i++) {
    182             psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
    183             psString kernelName = NULL;
    184             psStringAppend(&kernelName, "kernel%03d.fits", i);
    185             psFits *kernelFile = psFitsOpen(kernelName, "w");
    186             psFree(kernelName);
    187             psFitsWriteImage(kernelFile, NULL, kernel->image, 0, NULL);
    188             psFitsClose(kernelFile);
    189         }
    190     }
    191 
    192     return kernels;
    193 }
    194 
    195 /// Generate SPAM kernels
    196 pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, ///< Half-size of the kernel
    197                                                int spatialOrder, ///< Order of spatial variations
    198                                                int inner, ///< Inner radius to preserve unbinned
    199                                                int binning ///< Kernel binning factor
    200     )
    201 {
    202     PS_ASSERT_INT_POSITIVE(size, NULL);
    203     PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
    204     PS_ASSERT_INT_NONNEGATIVE(inner, NULL);
    205     PS_ASSERT_INT_LARGER_THAN(size, inner, NULL);
    206     PS_ASSERT_INT_POSITIVE(binning, NULL);
    207 
    208     // The outer region should be divisible by the "binning"; otherwise allocate remainder to the inner region
    209     int numOuter = (size - inner) / binning; // Number of summed pixels in the outer region
    210     int numInner = inner + (size - inner) % binning; // Number of pixels in the inner region
    211     assert(numOuter * binning + numInner == size);
    212     int numTotal = numOuter + numInner; // Total number of summed pixels
    213 
    214     psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter);
    215 
    216     int num = PS_SQR(2 * numTotal + 1); // Number of basis functions
    217 
    218     psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num);
    219 
    220     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM,
    221                                                               size, spatialOrder); // The kernels
    222     psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d)", size, inner, binning, spatialOrder);
    223 
    224     psLogMsg("psModules.imcombine", PS_LOG_INFO, "SPAM kernel: %d,%d,%d,%d --> %d elements",
    225              size, inner, binning, spatialOrder, num);
    226 
    227     kernels->uStop = psVectorAlloc(num, PS_TYPE_F32);
    228     kernels->vStop = psVectorAlloc(num, PS_TYPE_F32);
    229 
    230     psVector *locations = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Locations for each kernel element
    231     psVector *widths = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Widths for each kernel element
    232     locations->data.S32[numTotal] = 0;
    233     widths->data.S32[numTotal] = 0;
    234     for (int i = 1; i <= numInner; i++) {
    235         locations->data.S32[numTotal + i] = i;
    236         widths->data.S32[numTotal + i] = 0;
    237         locations->data.S32[numTotal - i] = - i;
    238         widths->data.S32[numTotal - i] = 0;
    239     }
    240     for (int i = numInner + 1; i <= numTotal; i++) {
    241         locations->data.S32[numTotal + i] = locations->data.S32[numTotal + i - 1] +
    242             widths->data.S32[numTotal + i - 1] + 1;
    243         widths->data.S32[numTotal + i] = binning - 1;
    244         locations->data.S32[numTotal - i] = locations->data.S32[numTotal - i + 1] - binning;
    245         widths->data.S32[numTotal - i] = binning - 1;
    246     }
    247 
    248     if (psTraceGetLevel("psModules.imcombine") >= 10) {
    249         for (int i = 0; i < 2 * numTotal + 1; i++) {
    250             psTrace("psModules.imcombine", 10, "%d: %d -> %d\n", i, locations->data.S32[i],
    251                     locations->data.S32[i] + widths->data.S32[i]);
    252         }
    253     }
    254 
    255     // Set the kernel parameters
    256     for (int i = - numTotal, index = 0; i <= numTotal; i++) {
    257         int u = locations->data.S32[numTotal + i]; // Location of pixel
    258         int uStop = u + widths->data.S32[numTotal + i]; // Width of pixel
    259 
    260         for (int j = - numTotal; j <= numTotal; j++, index++) {
    261             int v = locations->data.S32[numTotal + j]; // Location of pixel
    262             int vStop = v + widths->data.S32[numTotal + j]; // Width of pixel
    263 
    264             kernels->u->data.S32[index] = u;
    265             kernels->v->data.S32[index] = v;
    266             kernels->uStop->data.S32[index] = uStop;
    267             kernels->vStop->data.S32[index] = vStop;
    268 
    269             psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index,
    270                     u, uStop, v, vStop);
    271         }
    272     }
    273 
    274     kernels->subIndex = num / 2;
    275     assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
    276            kernels->v->data.S32[kernels->subIndex] == 0 &&
    277            kernels->uStop->data.S32[kernels->subIndex] == 0 &&
    278            kernels->vStop->data.S32[kernels->subIndex] == 0);
    279 
    280     psFree(locations);
    281     psFree(widths);
    282 
    283     return kernels;
    284 }
    285 
    286 
    287 /// Generate FRIES kernels
    288 pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, ///< Half-size of the kernel
    289                                                 int spatialOrder, ///< Order of spatial variations
    290                                                 int inner ///< Inner radius to preserve unbinned
    291     )
    292 {
    293     PS_ASSERT_INT_POSITIVE(size, NULL);
    294     PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
    295     PS_ASSERT_INT_NONNEGATIVE(inner, NULL);
    296     PS_ASSERT_INT_LARGER_THAN(size, inner, NULL);
    297 
    298     int fibNum = 0;                     // Number of Fibonacci values
    299     int fibLast = 1, fibTotal = 2;      // Fibonacci sequence
    300     while (fibTotal < size - inner) {
    301         int temp = fibTotal;
    302         fibTotal += fibLast;
    303         fibLast = temp;
    304         fibNum++;
    305     }
    306 
    307     int numInner = inner;               // Number of pixels in the inner region
    308     int numOuter = fibNum;              // Number of summed pixels in the outer region
    309     int numTotal = numOuter + numInner; // Total number of summed pixels
    310 
    311     psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter);
    312 
    313     int num = PS_SQR(2 * numTotal + 1); // Number of basis functions
    314 
    315     psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num);
    316 
    317     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES,
    318                                                               size, spatialOrder); // The kernels
    319     psStringAppend(&kernels->description, "FRIES(%d,%d,%d)", size, inner, spatialOrder);
    320 
    321     psLogMsg("psModules.imcombine", PS_LOG_INFO, "FRIES kernel: %d,%d,%d --> %d elements",
    322              size, inner, spatialOrder, num);
    323 
    324     kernels->uStop = psVectorAlloc(num, PS_TYPE_F32);
    325     kernels->vStop = psVectorAlloc(num, PS_TYPE_F32);
    326 
    327     psVector *start = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32);
    328     psVector *stop = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32);
    329     start->data.S32[numTotal] = 0;
    330     stop->data.S32[numTotal] = 0;
    331     for (int i = 1; i <= numInner; i++) {
    332         start->data.S32[numTotal + i] = i;
    333         stop->data.S32[numTotal + i] = i;
    334         start->data.S32[numTotal - i] = -i;
    335         stop->data.S32[numTotal - i] = -i;
    336     }
    337     for (int i = numInner + 1, fibLast = 1, fib = 2, temp; i <= numTotal;
    338          i++, fib = (temp = fib) + fibLast, fibLast = temp) {
    339         start->data.S32[numTotal + i] = stop->data.S32[numTotal + i - 1] + 1;
    340         stop->data.S32[numTotal + i] = PS_MIN(start->data.S32[numTotal + i] + fib - 1, size);
    341         start->data.S32[numTotal - i] = - stop->data.S32[numTotal + i];
    342         stop->data.S32[numTotal - i] = - start->data.S32[numTotal + i];
    343     }
    344 
    345     if (psTraceGetLevel("psModules.imcombine") >= 10) {
    346         for (int i = 0; i < 2 * numTotal + 1; i++) {
    347             psTrace("psModules.imcombine", 10, "%d: %d -> %d\n", i, start->data.S32[i], stop->data.S32[i]);
    348         }
    349     }
    350 
    351     // Set the kernel parameters
    352     for (int i = - numTotal, index = 0; i <= numTotal; i++) {
    353         int u = start->data.S32[numTotal + i]; // Location of pixel
    354         int uStop = stop->data.S32[numTotal + i]; // Width of pixel
    355         for (int j = - numTotal; j <= numTotal; j++, index++) {
    356             int v = start->data.S32[numTotal + j]; // Location of pixel
    357             int vStop = stop->data.S32[numTotal + j]; // Width of pixel
    358 
    359             kernels->u->data.S32[index] = u;
    360             kernels->v->data.S32[index] = v;
    361             kernels->uStop->data.S32[index] = uStop;
    362             kernels->vStop->data.S32[index] = vStop;
    363 
    364             psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index,
    365                     u, uStop, v, vStop);
    366         }
    367     }
    368 
    369     kernels->subIndex = num / 2;
    370     assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
    371            kernels->v->data.S32[kernels->subIndex] == 0 &&
    372            kernels->uStop->data.S32[kernels->subIndex] == 0 &&
    373            kernels->vStop->data.S32[kernels->subIndex] == 0);
    374 
    375     psFree(start);
    376     psFree(stop);
    377 
    378     return kernels;
    379 }
    380 
    381 // Grid United with Normal Kernel
    382 pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms,
    383                                                const psVector *orders, int inner)
    384 {
    385     PS_ASSERT_INT_POSITIVE(size, NULL);
    386     PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
    387     PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL);
    388     PS_ASSERT_VECTOR_TYPE(fwhms, PS_TYPE_F32, NULL);
    389     PS_ASSERT_VECTOR_NON_NULL(orders, NULL);
    390     PS_ASSERT_VECTOR_TYPE(orders, PS_TYPE_S32, NULL);
    391     PS_ASSERT_VECTORS_SIZE_EQUAL(fwhms, orders, NULL);
    392     PS_ASSERT_INT_NONNEGATIVE(inner, NULL);
    393     PS_ASSERT_INT_LESS_THAN(inner, size, NULL);
    394 
    395     int numGaussians = fwhms->n;       // Number of Gaussians
    396     int numGaussianVars = 0;            // Number of Gaussian variant functions in the kernel
    397     psString params = NULL;             // List of params
    398     for (int i = 0; i < numGaussians; i++) {
    399         int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian
    400         numGaussianVars += (gaussOrder + 1) * (gaussOrder + 2) / 2;
    401         psStringAppend(&params, "(%.2f,%d)", fwhms->data.F32[i], orders->data.S32[i]);
    402     }
    403 
    404     int numInner = PS_SQR(2 * inner + 1); // Number of inner kernel elements
    405 
    406     int num = numGaussianVars + numInner; // Total number of basis functions
    407 
    408     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_GUNK,
    409                                                               size, spatialOrder); // The kernels
    410     psStringAppend(&kernels->description, "GUNK(%d,%d,%s,%d)", size, inner, params, spatialOrder);
    411 
    412     psLogMsg("psModules.imcombine", PS_LOG_INFO, "GUNK kernel: %d,%s,%d --> %d elements",
    413              inner, params, spatialOrder, num);
    414     psFree(params);
    415 
    416     kernels->widths = psVectorAlloc(numGaussianVars, PS_TYPE_F32);
    417     kernels->preCalc = psArrayAlloc(numGaussianVars);
    418     kernels->inner = numGaussianVars;
    419     psKernel *subtract = NULL;          // Kernel to subtract to maintain flux scaling
    420 
    421     // Set the Gaussian kernel parameters
    422111    for (int i = 0, index = 0; i < numGaussians; i++) {
    423112        float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma
    … …  
    432121                    for (int u = -size; u <= size; u++) {
    433122                        sum += preCalc->kernel[v][u] = norm * power(u, uOrder) * power(v, vOrder) *
    434                             expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(fwhms->data.F32[i]));
     123                            expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigma));
    435124                    }
    436125                }
    437                 if (index == 0) {
    438                     subtract = preCalc;
     126
     127                // Normalise sum of kernel component to unity for even functions
     128                if (uOrder % 2 == 0 && vOrder % 2 == 0) {
    439129                    for (int v = -size; v <= size; v++) {
    440130                        for (int u = -size; u <= size; u++) {
    441                             preCalc->kernel[v][u] /= sum;
    442                         }
    443                     }
    444                 } else if (uOrder % 2 == 0 && vOrder % 2 == 0) {
    445                     // Normalise sum of kernel component to unity for even functions
    446                     for (int v = -size; v <= size; v++) {
    447                         for (int u = -size; u <= size; u++) {
    448                             preCalc->kernel[v][u] = preCalc->kernel[v][u] / sum - subtract->kernel[v][u];
     131                            preCalc->kernel[v][u] = preCalc->kernel[v][u] / sum;
    449132                        }
    450133                    }
    … …  
    454137                kernels->u->data.S32[index] = uOrder;
    455138                kernels->v->data.S32[index] = vOrder;
     139                if (kernels->preCalc->data[index]) {
     140                    psFree(kernels->preCalc->data[index]);
     141                }
    456142                kernels->preCalc->data[index] = preCalc;
    457143
    … …  
    462148    }
    463149
    464     // Generate a grid set of kernels for each (u,v)
    465     for (int v = - inner, index = kernels->inner; v <= inner; v++) {
    466         for (int u = - inner; u <= inner; u++, index++) {
     150    kernels->subIndex = -1;
     151
     152    if (psTraceGetLevel("psModules.imcombine.kernel") >= 10) {
     153        for (int i = 0; i < num; i++) {
     154            psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
     155            psString kernelName = NULL;
     156            psStringAppend(&kernelName, "kernel%03d.fits", i);
     157            psFits *kernelFile = psFitsOpen(kernelName, "w");
     158            psFree(kernelName);
     159            psFitsWriteImage(kernelFile, NULL, kernel->image, 0, NULL);
     160            psFitsClose(kernelFile);
     161        }
     162    }
     163
     164    return kernels;
     165}
     166
     167//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     168// Public functions
     169//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     170
     171pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type,
     172                                                int size, int spatialOrder)
     173{
     174    pmSubtractionKernels *kernels = psAlloc(sizeof(pmSubtractionKernels)); // Kernels, to return
     175    psMemSetDeallocator(kernels, (psFreeFunc)subtractionKernelsFree);
     176
     177    kernels->type = type;
     178    kernels->description = NULL;
     179    kernels->num = numBasisFunctions;
     180    kernels->u = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);
     181    kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);
     182    kernels->widths = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
     183    kernels->preCalc = psArrayAlloc(numBasisFunctions);
     184    kernels->uStop = NULL;
     185    kernels->vStop = NULL;
     186    kernels->subIndex = 0;
     187    kernels->size = size;
     188    kernels->inner = 0;
     189    kernels->spatialOrder = spatialOrder;
     190    kernels->bgOrder = 0;
     191
     192    return kernels;
     193}
     194
     195pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder)
     196{
     197    PS_ASSERT_INT_POSITIVE(size, NULL);
     198    PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
     199
     200    int num = PS_SQR(2 * size + 1); // Number of basis functions
     201
     202    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS,
     203                                                              size, spatialOrder); // The kernels
     204    psStringAppend(&kernels->description, "POIS(%d,%d)", size, spatialOrder);
     205    psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements",
     206             size, spatialOrder, num);
     207
     208    if (!p_pmSubtractionKernelsAddGrid(kernels, 0, size)) {
     209        psAbort("Should never get here.");
     210    }
     211
     212    kernels->subIndex = num / 2;
     213    assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
     214           kernels->v->data.S32[kernels->subIndex] == 0);
     215
     216    return kernels;
     217}
     218
     219
     220pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder,
     221                                               const psVector *fwhms, const psVector *orders)
     222{
     223    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder,
     224                                                                  fwhms, orders); // Kernels
     225    if (!kernels) {
     226        return NULL;
     227    }
     228
     229    // Subtract a reference kernel from all the others to maintain flux scaling
     230    psKernel *subtract = kernels->preCalc->data[0]; // Kernel to subtract from all others
     231    for (int i = 1; i < kernels->num; i++) {
     232        psKernel *kernel = kernels->preCalc->data[i]; // Kernel to subtract
     233        psBinaryOp(kernel->image, kernel->image, "-", subtract->image);
     234    }
     235
     236    return kernels;
     237}
     238
     239/// Generate SPAM kernels
     240pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, ///< Half-size of the kernel
     241                                               int spatialOrder, ///< Order of spatial variations
     242                                               int inner, ///< Inner radius to preserve unbinned
     243                                               int binning ///< Kernel binning factor
     244    )
     245{
     246    PS_ASSERT_INT_POSITIVE(size, NULL);
     247    PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
     248    PS_ASSERT_INT_NONNEGATIVE(inner, NULL);
     249    PS_ASSERT_INT_LARGER_THAN(size, inner, NULL);
     250    PS_ASSERT_INT_POSITIVE(binning, NULL);
     251
     252    // The outer region should be divisible by the "binning"; otherwise allocate remainder to the inner region
     253    int numOuter = (size - inner) / binning; // Number of summed pixels in the outer region
     254    int numInner = inner + (size - inner) % binning; // Number of pixels in the inner region
     255    assert(numOuter * binning + numInner == size);
     256    int numTotal = numOuter + numInner; // Total number of summed pixels
     257
     258    psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter);
     259
     260    int num = PS_SQR(2 * numTotal + 1); // Number of basis functions
     261
     262    psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num);
     263
     264    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM,
     265                                                              size, spatialOrder); // The kernels
     266    psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d)", size, inner, binning, spatialOrder);
     267
     268    psLogMsg("psModules.imcombine", PS_LOG_INFO, "SPAM kernel: %d,%d,%d,%d --> %d elements",
     269             size, inner, binning, spatialOrder, num);
     270
     271    kernels->uStop = psVectorAlloc(num, PS_TYPE_F32);
     272    kernels->vStop = psVectorAlloc(num, PS_TYPE_F32);
     273
     274    psVector *locations = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Locations for each kernel element
     275    psVector *widths = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Widths for each kernel element
     276    locations->data.S32[numTotal] = 0;
     277    widths->data.S32[numTotal] = 0;
     278    for (int i = 1; i <= numInner; i++) {
     279        locations->data.S32[numTotal + i] = i;
     280        widths->data.S32[numTotal + i] = 0;
     281        locations->data.S32[numTotal - i] = - i;
     282        widths->data.S32[numTotal - i] = 0;
     283    }
     284    for (int i = numInner + 1; i <= numTotal; i++) {
     285        locations->data.S32[numTotal + i] = locations->data.S32[numTotal + i - 1] +
     286            widths->data.S32[numTotal + i - 1] + 1;
     287        widths->data.S32[numTotal + i] = binning - 1;
     288        locations->data.S32[numTotal - i] = locations->data.S32[numTotal - i + 1] - binning;
     289        widths->data.S32[numTotal - i] = binning - 1;
     290    }
     291
     292    if (psTraceGetLevel("psModules.imcombine") >= 10) {
     293        for (int i = 0; i < 2 * numTotal + 1; i++) {
     294            psTrace("psModules.imcombine", 10, "%d: %d -> %d\n", i, locations->data.S32[i],
     295                    locations->data.S32[i] + widths->data.S32[i]);
     296        }
     297    }
     298
     299    // Set the kernel parameters
     300    for (int i = - numTotal, index = 0; i <= numTotal; i++) {
     301        int u = locations->data.S32[numTotal + i]; // Location of pixel
     302        int uStop = u + widths->data.S32[numTotal + i]; // Width of pixel
     303
     304        for (int j = - numTotal; j <= numTotal; j++, index++) {
     305            int v = locations->data.S32[numTotal + j]; // Location of pixel
     306            int vStop = v + widths->data.S32[numTotal + j]; // Width of pixel
     307
    467308            kernels->u->data.S32[index] = u;
    468309            kernels->v->data.S32[index] = v;
    469 
    470             psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v);
    471         }
    472     }
    473 
    474     kernels->subIndex = numInner/2 + numGaussianVars;
     310            kernels->uStop->data.S32[index] = uStop;
     311            kernels->vStop->data.S32[index] = vStop;
     312
     313            psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index,
     314                    u, uStop, v, vStop);
     315        }
     316    }
     317
     318    kernels->subIndex = num / 2;
     319    assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
     320           kernels->v->data.S32[kernels->subIndex] == 0 &&
     321           kernels->uStop->data.S32[kernels->subIndex] == 0 &&
     322           kernels->vStop->data.S32[kernels->subIndex] == 0);
     323
     324    psFree(locations);
     325    psFree(widths);
     326
     327    return kernels;
     328}
     329
     330
     331/// Generate FRIES kernels
     332pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, ///< Half-size of the kernel
     333                                                int spatialOrder, ///< Order of spatial variations
     334                                                int inner ///< Inner radius to preserve unbinned
     335    )
     336{
     337    PS_ASSERT_INT_POSITIVE(size, NULL);
     338    PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
     339    PS_ASSERT_INT_NONNEGATIVE(inner, NULL);
     340    PS_ASSERT_INT_LARGER_THAN(size, inner, NULL);
     341
     342    int fibNum = 0;                     // Number of Fibonacci values
     343    int fibLast = 1, fibTotal = 2;      // Fibonacci sequence
     344    while (fibTotal < size - inner) {
     345        int temp = fibTotal;
     346        fibTotal += fibLast;
     347        fibLast = temp;
     348        fibNum++;
     349    }
     350
     351    int numInner = inner;               // Number of pixels in the inner region
     352    int numOuter = fibNum;              // Number of summed pixels in the outer region
     353    int numTotal = numOuter + numInner; // Total number of summed pixels
     354
     355    psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter);
     356
     357    int num = PS_SQR(2 * numTotal + 1); // Number of basis functions
     358
     359    psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num);
     360
     361    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES,
     362                                                              size, spatialOrder); // The kernels
     363    psStringAppend(&kernels->description, "FRIES(%d,%d,%d)", size, inner, spatialOrder);
     364
     365    psLogMsg("psModules.imcombine", PS_LOG_INFO, "FRIES kernel: %d,%d,%d --> %d elements",
     366             size, inner, spatialOrder, num);
     367
     368    kernels->uStop = psVectorAlloc(num, PS_TYPE_F32);
     369    kernels->vStop = psVectorAlloc(num, PS_TYPE_F32);
     370
     371    psVector *start = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32);
     372    psVector *stop = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32);
     373    start->data.S32[numTotal] = 0;
     374    stop->data.S32[numTotal] = 0;
     375    for (int i = 1; i <= numInner; i++) {
     376        start->data.S32[numTotal + i] = i;
     377        stop->data.S32[numTotal + i] = i;
     378        start->data.S32[numTotal - i] = -i;
     379        stop->data.S32[numTotal - i] = -i;
     380    }
     381    for (int i = numInner + 1, fibLast = 1, fib = 2, temp; i <= numTotal;
     382         i++, fib = (temp = fib) + fibLast, fibLast = temp) {
     383        start->data.S32[numTotal + i] = stop->data.S32[numTotal + i - 1] + 1;
     384        stop->data.S32[numTotal + i] = PS_MIN(start->data.S32[numTotal + i] + fib - 1, size);
     385        start->data.S32[numTotal - i] = - stop->data.S32[numTotal + i];
     386        stop->data.S32[numTotal - i] = - start->data.S32[numTotal + i];
     387    }
     388
     389    if (psTraceGetLevel("psModules.imcombine") >= 10) {
     390        for (int i = 0; i < 2 * numTotal + 1; i++) {
     391            psTrace("psModules.imcombine", 10, "%d: %d -> %d\n", i, start->data.S32[i], stop->data.S32[i]);
     392        }
     393    }
     394
     395    // Set the kernel parameters
     396    for (int i = - numTotal, index = 0; i <= numTotal; i++) {
     397        int u = start->data.S32[numTotal + i]; // Location of pixel
     398        int uStop = stop->data.S32[numTotal + i]; // Width of pixel
     399        for (int j = - numTotal; j <= numTotal; j++, index++) {
     400            int v = start->data.S32[numTotal + j]; // Location of pixel
     401            int vStop = stop->data.S32[numTotal + j]; // Width of pixel
     402
     403            kernels->u->data.S32[index] = u;
     404            kernels->v->data.S32[index] = v;
     405            kernels->uStop->data.S32[index] = uStop;
     406            kernels->vStop->data.S32[index] = vStop;
     407
     408            psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index,
     409                    u, uStop, v, vStop);
     410        }
     411    }
     412
     413    kernels->subIndex = num / 2;
     414    assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
     415           kernels->v->data.S32[kernels->subIndex] == 0 &&
     416           kernels->uStop->data.S32[kernels->subIndex] == 0 &&
     417           kernels->vStop->data.S32[kernels->subIndex] == 0);
     418
     419    psFree(start);
     420    psFree(stop);
     421
     422    return kernels;
     423}
     424
     425// Grid United with Normal Kernel
     426pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms,
     427                                               const psVector *orders, int inner)
     428{
     429    PS_ASSERT_INT_POSITIVE(size, NULL);
     430    PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
     431    PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL);
     432    PS_ASSERT_VECTOR_TYPE(fwhms, PS_TYPE_F32, NULL);
     433    PS_ASSERT_VECTOR_NON_NULL(orders, NULL);
     434    PS_ASSERT_VECTOR_TYPE(orders, PS_TYPE_S32, NULL);
     435    PS_ASSERT_VECTORS_SIZE_EQUAL(fwhms, orders, NULL);
     436    PS_ASSERT_INT_NONNEGATIVE(inner, NULL);
     437    PS_ASSERT_INT_LESS_THAN(inner, size, NULL);
     438
     439    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder,
     440                                                                  fwhms, orders); // Kernels
     441    psStringPrepend(&kernels->description, "GUNK=");
     442    psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder);
     443
     444    // Subtract unity from the kernels to maintain photometric flux scaling
     445    for (int i = 0; i < kernels->num; i++) {
     446        psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
     447        kernel->kernel[0][0] -= 1.0;
     448    }
     449
     450    int numGaussians = kernels->num;    // Number of ISIS kernels
     451    int numInner = PS_SQR(2 * inner + 1); // Number of inner kernel elements
     452
     453    if (!p_pmSubtractionKernelsAddGrid(kernels, numGaussians, inner)) {
     454        psAbort("Should never get here.");
     455    }
     456
     457    kernels->subIndex = numInner/2 + numGaussians;
    475458    assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
    476459           kernels->v->data.S32[kernels->subIndex] == 0);
    … …  
    518501    psLogMsg("psModules.imcombine", PS_LOG_INFO, "RINGS kernel: %d,%d,%d,%d --> %d elements",
    519502             size, inner, ringsOrder, spatialOrder, num);
    520 
    521     kernels->preCalc = psArrayAlloc(num);
    522503
    523504    // Set the Gaussian kernel parameters
Note: See TracChangeset for help on using the changeset viewer.