IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 14671


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.

Location:
trunk/psModules/src
Files:
2 added
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/imcombine/Makefile.am

    r14632 r14671  
    1111        pmSubtractionKernels.c  \
    1212        pmSubtractionStamps.c   \
    13         pmSubtractionMatch.c
     13        pmSubtractionMatch.c    \
     14        pmSubtractionParams.c
    1415
    1516pkginclude_HEADERS = \
    … …  
    2021        pmSubtractionKernels.h  \
    2122        pmSubtractionStamps.h   \
    22         pmSubtractionMatch.h
     23        pmSubtractionMatch.h    \
     24        pmSubtractionParams.h
    2325
    2426CLEANFILES = *~

  • 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

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

    r14533 r14671  
    3131} pmSubtractionKernels;
    3232
     33/// Generate a delta-function grid for subtraction kernels (like the POIS kernel)
     34bool p_pmSubtractionKernelsAddGrid(pmSubtractionKernels *kernels, ///< The subtraction kernels to append to
     35                                   int start, ///< Index at which to start appending
     36                                   int size ///< Half-size of the grid
     37    );
     38
    3339/// General allocator for pmSubtractionKernels
    3440///
    … …  
    4450pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, ///< Half-size of the kernel (in both dims)
    4551                                               int spatialOrder ///< Order of spatial variations
     52    );
     53
     54/// Generate ISIS kernels without the flux scaling built in
     55pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, ///< Half-size of the kernel
     56                                                    int spatialOrder, ///< Order of spatial variations
     57                                                    const psVector *fwhms, ///< Gaussian FWHMs
     58                                                    const psVector *orders ///< Polynomial order of gaussians
    4659    );
    4760

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

    r14606 r14671  
    99#include "pmHDU.h"
    1010#include "pmFPA.h"
     11#include "pmSubtractionParams.h"
    1112#include "pmSubtractionKernels.h"
    1213#include "pmSubtractionStamps.h"
    … …  
    4445
    4546
     47static bool getStamps(psArray **stamps, // Stamps to read
     48                      const psArray *stampsData, // Stamp data from a file
     49                      const pmReadout *reference, // Reference readout
     50                      const pmReadout *input, // Input readout, or NULL to generate fake stamps
     51                      const psImage *subMask, // Mask for subtraction, or NULL
     52                      psImage *weight,  // Weight map
     53                      const psRegion *region, // Region of interest, or NULL
     54                      float threshold,  // Threshold for stamp finding
     55                      float stampSpacing, // Spacing between stamps
     56                      float targetWidth,// Target width for fake stamps
     57                      int size,         // Kernel half-size
     58                      int footprint     // Convolution footprint for stamps
     59    )
     60{
     61
     62    psImage *inImage = NULL; // Input image
     63    if (input) {
     64        inImage = input->image;
     65    }
     66
     67    if (stampsData) {
     68        psVector *xStamp = NULL, *yStamp = NULL, *fluxStamp = NULL; // Stamp positions and fluxes
     69        if (input) {
     70            // We have x, y because the target is provided by the input image
     71            xStamp = stampsData->data[0];
     72            yStamp = stampsData->data[1];
     73        } else {
     74            // We have x, y and flux in order to generate a target
     75            xStamp = stampsData->data[0];
     76            yStamp = stampsData->data[1];
     77            fluxStamp = stampsData->data[2];
     78        }
     79
     80        psFree(*stamps);
     81        *stamps = pmSubtractionSetStamps(xStamp, yStamp, fluxStamp, subMask, region);
     82    } else {
     83        psTrace("psModules.imcombine", 3, "Finding stamps...\n");
     84        *stamps = pmSubtractionFindStamps(*stamps, reference->image, subMask, region,
     85                                          threshold, stampSpacing);
     86    }
     87    if (!stamps) {
     88        psError(PS_ERR_UNKNOWN, false, "Unable to find stamps.");
     89        return false;
     90    }
     91
     92    memCheck("  find stamps");
     93
     94    if (!input && !pmSubtractionGenerateStamps(*stamps, targetWidth, footprint, size)) {
     95        psError(PS_ERR_UNKNOWN, false, "Unable to generate target stamps.");
     96        return false;
     97    }
     98
     99    psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
     100    if (!pmSubtractionExtractStamps(*stamps, reference->image, inImage, weight, footprint, size)) {
     101        psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps.");
     102        return false;
     103    }
     104
     105    memCheck("   extract stamps");
     106
     107    return true;
     108}
     109
     110
     111
     112
     113
    46114bool pmSubtractionMatch(pmReadout *convolved, const pmReadout *reference, const pmReadout *input,
    47115                        int footprint, float regionSize, float stampSpacing, float threshold,
    48116                        const char *stampsName, float targetWidth, pmSubtractionKernelsType type,
    49                         int size, int order, const psVector *isisWidths, const psVector *isisOrders,
    50                         int inner, int ringsOrder, int binning, int iter, float rej, psMaskType maskBad,
     117                        int size, int spatialOrder, const psVector *isisWidths, const psVector *isisOrders,
     118                        int inner, int ringsOrder, int binning, bool optimum, psVector *optFWHMs,
     119                        int optOrder, float optThreshold, int iter, float rej, psMaskType maskBad,
    51120                        psMaskType maskBlank)
    52121{
    … …  
    88157    // We'll check kernel type when we allocate the kernels
    89158    PS_ASSERT_INT_POSITIVE(size, false);
    90     PS_ASSERT_INT_NONNEGATIVE(order, false);
     159    PS_ASSERT_INT_NONNEGATIVE(spatialOrder, false);
    91160    if (isisWidths || isisOrders) {
    92161        PS_ASSERT_VECTOR_NON_NULL(isisWidths, false);
    … …  
    99168    PS_ASSERT_INT_NONNEGATIVE(ringsOrder, false);
    100169    PS_ASSERT_INT_POSITIVE(binning, false);
     170    if (optimum) {
     171        PS_ASSERT_VECTOR_NON_NULL(optFWHMs, false);
     172        PS_ASSERT_INT_NONNEGATIVE(optOrder, false);
     173        PS_ASSERT_FLOAT_LARGER_THAN(optThreshold, 0.0, false);
     174        PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(optThreshold, 1.0, false);
     175    }
    101176    PS_ASSERT_INT_POSITIVE(iter, false);
    102177    PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false);
    … …  
    136211    }
    137212
     213    // Read stamps from file
     214    psArray *stampsData = NULL;         // Stamps data read from file
     215    if (stampsName && strlen(stampsName) > 0) {
     216        psTrace("psModules.imcombine", 3, "Reading stamps from %s...\n", stampsName);
     217        psVector *xStamp = NULL, *yStamp = NULL; // Stamp positions and fluxes
     218        if (input) {
     219            // We have x, y because the target is provided by the input image
     220            stampsData = psVectorsReadFromFile(stampsName, "%f %f"); // Stamp positions
     221        } else {
     222            // We have x, y and flux in order to generate a target
     223            stampsData = psVectorsReadFromFile(stampsName, "%f %f %f"); // Stamp positions
     224        }
     225
     226        // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions
     227        xStamp = stampsData->data[0];
     228        yStamp = stampsData->data[1];
     229        psBinaryOp(xStamp, xStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
     230        psBinaryOp(yStamp, yStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
     231    }
     232
    138233    int numCols = reference->image->numCols, numRows = reference->image->numRows; // Image dimensions
    139234
    … …  
    145240    memCheck("mask");
    146241
    147     // Kernel basis functions
    148     pmSubtractionKernels *kernels = pmSubtractionKernelsGenerate(type, size, order, isisWidths, isisOrders,
    149                                                                  inner, binning, ringsOrder);
     242
     243    psArray *stamps = NULL;             // Stamps for matching PSF
     244
     245    pmSubtractionKernels *kernels = NULL; // Kernel basis functions
     246    if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
     247        if (!getStamps(&stamps, stampsData, reference, input, subMask, weight, NULL,
     248                       threshold, stampSpacing, targetWidth, size, footprint)) {
     249            goto ERROR;
     250        }
     251        kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
     252                                                  stamps, footprint, optThreshold);
     253        // XXX This is not quite the best thing to do here--- we'll find the same set of stamps again later,
     254        // but since we may not be interested in the entire image on the first pass through, we have to blow
     255        // the whole lot away.  The alternative is to mark stamps that are outside the region of interest, and
     256        // ignore them when generating sums and rejecting.
     257        psFree(stamps);
     258        stamps = NULL;
     259
     260        if (!kernels) {
     261            psErrorClear();
     262            psWarning("Unable to derive optimum ISIS kernel --- switching to default.");
     263        }
     264    }
     265
     266    if (kernels == NULL) {
     267        // Not an ISIS/GUNK kernel, or the optimum kernel search failed
     268        kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
     269                                               inner, binning, ringsOrder);
     270    }
     271
    150272    psMetadataAddPtr(convolved->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL", PS_DATA_UNKNOWN,
    151273                     "Subtraction kernels", kernels);
    152     psArray *stamps = NULL;             // Stamps for matching PSF
    153274    psVector *solution = NULL;          // Solution to match PSF
    154275
    … …  
    186307                psTrace("psModules.imcombine", 2, "Iteration %d...\n", k);
    187308
    188                 if (stampsName && strlen(stampsName) > 0) {
    189                     psTrace("psModules.imcombine", 3, "Reading stamps from %s...\n", stampsName);
    190                     iter = 1;           // There is no iterating because we use all the stamps we have
    191                     psVector *xStamp = NULL, *yStamp = NULL, *fluxStamp = NULL; // Stamp positions and fluxes
    192                     if (input) {
    193                         // We have x, y because the target is provided by the input image
    194                         psArray *stampData = psVectorsReadFromFile(stampsName, "%f %f"); // Stamp positions
    195                         xStamp = stampData->data[0];
    196                         yStamp = stampData->data[1];
    197                     } else {
    198                         // We have x, y and flux in order to generate a target
    199                         psArray *stampData = psVectorsReadFromFile(stampsName, "%f %f %f"); // Stamp positions
    200                         xStamp = stampData->data[0];
    201                         yStamp = stampData->data[1];
    202                         fluxStamp = stampData->data[2];
    203                     }
    204 
    205                     // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions
    206                     psBinaryOp(xStamp, xStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
    207                     psBinaryOp(yStamp, yStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
    208 
    209                     stamps = pmSubtractionSetStamps(xStamp, yStamp, fluxStamp, subMask, region);
    210                 } else {
    211                     psTrace("psModules.imcombine", 3, "Finding stamps...\n");
    212                     stamps = pmSubtractionFindStamps(stamps, reference->image, subMask, region,
    213                                                      threshold, stampSpacing);
    214                 }
    215                 if (!stamps) {
    216                     psError(PS_ERR_UNKNOWN, false, "Unable to find stamps.");
     309                if (!getStamps(&stamps, stampsData, reference, input, subMask, weight, region,
     310                               threshold, stampSpacing, targetWidth, size, footprint)) {
    217311                    goto ERROR;
    218312                }
    219 
    220                 memCheck("  find stamps");
    221 
    222                 if (!input && !pmSubtractionGenerateStamps(stamps, targetWidth, footprint, kernels)) {
    223                     psError(PS_ERR_UNKNOWN, false, "Unable to generate target stamps.");
    224                     goto ERROR;
    225                 }
    226 
    227                 psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
    228                 if (!pmSubtractionExtractStamps(stamps, reference->image, inImage, weight,
    229                                                 footprint, kernels)) {
    230                     psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps.");
    231                     goto ERROR;
    232                 }
    233 
    234                 memCheck("   extract stamps");
    235313
    236314                psTrace("psModules.imcombine", 3, "Calculating equation...\n");
    … …  
    376454    psFree(subMask);
    377455    subMask = NULL;
     456    psFree(stampsData);
     457    stampsData = NULL;
    378458
    379459    if (!pmSubtractionBorder(convolved->image, convolved->weight, convolved->mask, kernels, maskBlank)) {
    … …  
    395475    psFree(stamps);
    396476    psFree(solution);
     477    psFree(stampsData);
    397478    return false;
    398479}

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

    r14584 r14671  
    2828                        int ringsOrder, ///< RINGS polynomial order
    2929                        int binning,    ///< SPAM kernel binning
     30                        bool optimum,   ///< Search for optimum ISIS kernel?
     31                        psVector *optFWHMs, ///< FWHMs for optimum search
     32                        int optOrder,   ///< Maximum order for optimum search
     33                        float optThreshold, ///< Threshold for optimum search (0..1)
    3034                        // Operational parameters
    3135                        int iter,       ///< Rejection iterations

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

    r14534 r14671  
    239239
    240240bool pmSubtractionExtractStamps(psArray *stamps, psImage *reference, psImage *input, psImage *weight,
    241                                 int footprint, const pmSubtractionKernels *kernels)
     241                                int footprint, int kernelSize)
    242242{
    243243    PS_ASSERT_ARRAY_NON_NULL(stamps, false);
    … …  
    254254        PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false);
    255255    }
     256    PS_ASSERT_INT_NONNEGATIVE(footprint, false);
     257    PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
    256258
    257259    int numCols = reference->numCols, numRows = reference->numRows; // Size of images
    258     int size = kernels->size + footprint; // Size of postage stamps
     260    int size = kernelSize + footprint; // Size of postage stamps
    259261
    260262    if (!weight) {
    … …  
    303305
    304306
    305 bool pmSubtractionGenerateStamps(psArray *stamps, float fwhm, int footprint,
    306                                  const pmSubtractionKernels *kernels)
     307bool pmSubtractionGenerateStamps(psArray *stamps, float fwhm, int footprint, int kernelSize)
    307308{
    308309    PS_ASSERT_ARRAY_NON_NULL(stamps, false);
    … …  
    310311    PS_ASSERT_INT_NONNEGATIVE(footprint, false);
    311312
    312     int size = kernels->size + footprint; // Size of postage stamps
     313    int size = kernelSize + footprint; // Size of postage stamps
    313314    int num = stamps->n;                // Number of stamps
    314315    float sigma = fwhm / (2.0 * sqrtf(2.0 * log(2.0))); // Gaussian sigma

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

    r14533 r14671  
    5151                                psImage *weight, ///< Weight (variance) map
    5252                                int footprint, ///< Stamp footprint size
    53                                 const pmSubtractionKernels *kernels ///< Kernel basis functions (for size)
     53                                int kernelSize ///< Kernel half-size
    5454    );
    5555
    … …  
    5858                                 float fwhm, ///< FWHM for each stamp
    5959                                 int footprint, ///< Stamp footprint size
    60                                  const pmSubtractionKernels *kernels ///< Kernel basis functions
     60                                 int size ///< Kernel half-size
    6161    );
    6262

  • trunk/psModules/src/psmodules.h

    r14653 r14671  
    7979#include <pmSubtractionKernels.h>
    8080#include <pmSubtractionMatch.h>
     81#include <pmSubtractionParams.h>
    8182#include <pmReadoutCombine.h>
    8283
Note: See TracChangeset for help on using the changeset viewer.