IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 3, 2010, 8:41:49 AM (16 years ago)
Author:
eugene
Message:

updates from trunk

Location:
branches/tap_branches
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/tap_branches

  • branches/tap_branches/psModules

  • branches/tap_branches/psModules/src/imcombine/pmSubtractionParams.c

    r21363 r27838  
    88#include <pslib.h>
    99
     10#include "pmErrorCodes.h"
    1011#include "pmSubtractionStamps.h"
    1112#include "pmSubtraction.h"
    … …  
    7172                            double *sumII, // Sum of I(x)^2/sigma(x)^2
    7273                            double *sumIC, // Sum of I(x)conv(x)/sigma(x)^2
    73                             const pmSubtractionStamp *stamp, // Stamp with variance
     74                            const pmSubtractionStamp *stamp, // Stamp
    7475                            const psKernel *target, // Target stamp
    7576                            int kernelIndex, // Index for kernel component
    … …  
    7879    )
    7980{
    80     psKernel *variance = stamp->variance;   // Variance, sigma(x)^2
     81    psKernel *weight = stamp->weight;   // Weight image
    8182    psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest
    8283
    8384    for (int y = -footprint; y <= footprint; y++) {
    8485        psF32 *in = &target->kernel[y][-footprint]; // Dereference input
    85         psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance
     86        psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
    8687        psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
    8788        for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) {
    88             double temp = *in / *wt; // Temporary product
     89            double temp = *in * *wt; // Temporary product
    8990            *sumI += temp;
    9091            *sumII += *in * temp;
    … …  
    9899static void accumulateConvolutions(double *sumC, // Sum of conv(x)/sigma(x)^2
    99100                                   double *sumCC, // Sum of conv(x)^2/sigma(x)^2
    100                                    const pmSubtractionStamp *stamp, // Stamp with input and variance
     101                                   const pmSubtractionStamp *stamp, // Stamp with input and weight
    101102                                   int kernelIndex, // Index for kernel component
    102103                                   int footprint, // Size of region of interest
    … …  
    104105    )
    105106{
    106     psKernel *variance = stamp->variance;   // Variance, sigma(x)^2
     107    psKernel *weight = stamp->weight;   // Weight image
    107108    psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest
    108109
    109110    for (int y = -footprint; y <= footprint; y++) {
    110         psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance
     111        psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
    111112        psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
    112113        for (int x = -footprint; x <= footprint; x++, wt++, conv++) {
    113             double convNoise = *conv / *wt; // Temporary product
     114            double convNoise = *conv * *wt; // Temporary product
    114115            *sumC += convNoise;
    115116            *sumCC += *conv * convNoise;
    … …  
    120121
    121122static double accumulateChi2(const psKernel *target, // Target stamp
    122                              pmSubtractionStamp *stamp, // Stamp with variance
     123                             pmSubtractionStamp *stamp, // Stamp with weight
    123124                             int kernelIndex, // Index for kernel component
    124125                             double coeff, // Coefficient of convolution
    … …  
    129130{
    130131    double chi2 = 0.0;
    131     psKernel *variance = stamp->variance;   // Variance, sigma(x)^2
     132    psKernel *weight = stamp->weight;   // Weight image
    132133    psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest
    133134
    134135    for (int y = -footprint; y <= footprint; y++) {
    135136        psF32 *in = &target->kernel[y][-footprint]; // Dereference input
    136         psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance
     137        psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
    137138        psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
    138139        for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) {
    139             chi2 += PS_SQR(*in - bg - coeff * *conv) / *wt;
     140            chi2 += PS_SQR(*in - bg - coeff * *conv) * *wt;
    140141        }
    141142    }
    … …  
    146147// Return the initial value of chi^2
    147148static double initialChi2(const psKernel *target, // Target stamp
    148                           const pmSubtractionStamp *stamp, // Stamp with variance
     149                          const pmSubtractionStamp *stamp, // Stamp
    149150                          int footprint, // Size of convolution
    150151                          pmSubtractionMode mode // Mode of subtraction
    151152    )
    152153{
    153     psKernel *variance = stamp->variance;   // Variance map
     154    psKernel *weight = stamp->weight;   // Weight image
    154155    psKernel *source;                   // Source stamp
    155156    switch (mode) {
    … …  
    167168    for (int y = -footprint; y <= footprint; y++) {
    168169        psF32 *in = &target->kernel[y][-footprint]; // Dereference input
    169         psF32 *wt = &variance->kernel[y][-footprint]; // Dereference variance
     170        psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
    170171        psF32 *ref = &source->kernel[y][-footprint]; // Derference reference
    171172        for (int x = -footprint; x <= footprint; x++, in++, wt++, ref++) {
    172173            float diff = *in - *ref;    // Temporary value
    173             chi2 += PS_SQR(diff) / *wt;
     174            chi2 += PS_SQR(diff) * *wt;
    174175        }
    175176    }
    … …  
    180181// Subtract a convolution from the input
    181182static void subtractConvolution(psKernel *target, // Target stamp
    182                                 const pmSubtractionStamp *stamp, // Stamp with variance
     183                                const pmSubtractionStamp *stamp, // Stamp
    183184                                int kernelIndex, // Index for kernel component
    184185                                float coeff, // Coefficient of subtraction
    … …  
    204205                                                      int spatialOrder, const psVector *fwhms, int maxOrder,
    205206                                                      const pmSubtractionStampList *stamps, int footprint,
    206                                                       float tolerance, float penalty, pmSubtractionMode mode)
     207                                                      float tolerance, float penalty, psRegion bounds,
     208                                                      pmSubtractionMode mode)
    207209{
    208210    if (type != PM_SUBTRACTION_KERNEL_ISIS && type != PM_SUBTRACTION_KERNEL_GUNK) {
    … …  
    232234    psVectorInit(orders, maxOrder);
    233235    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
    234                                                                   penalty, mode); // Kernels
     236                                                                  penalty, bounds, mode); // Kernels
    235237    psFree(orders);
    236238    psFree(kernels->description);
    … …  
    280282        }
    281283        if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) {
    282             psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i);
     284            psError(psErrorCodeLast(), false, "Unable to convolve stamp %d.", i);
    283285            psFree(targets);
    284286            psFree(kernels);
    … …  
    288290
    289291        // This sum is invariant to the kernel
    290         psKernel *variance = stamp->variance; // Variance map for stamp
     292        psKernel *weight = stamp->weight; // Weight image
    291293        for (int v = -footprint; v <= footprint; v++) {
    292             psF32 *wt = &variance->kernel[v][-footprint]; // Dereference variance map
     294            psF32 *wt = &weight->kernel[v][-footprint]; // Dereference weight
    293295            for (int u = -footprint; u <= footprint; u++, wt++) {
    294                 sum1 += 1.0 / *wt;
     296                sum1 += 1.0 * *wt;
    295297            }
    296298        }
    297299        if (!isfinite(sum1)) {
    298             psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     300            psError(PM_ERR_DATA, true,
    299301                    "Sum of 1/sigma^2 is non-finite for stamp %d (%d,%d)\n",
    300302                    i, (int)stamp->x, (int)stamp->y);
    … …  
    368370
    369371        if (bestIndex == -1) {
    370             psError(PS_ERR_UNKNOWN, false, "Unable to find best kernel component in round %d.", iter);
     372            psError(PM_ERR_DATA, true, "Unable to find best kernel component in round %d.", iter);
    371373            psFree(targets);
    372374            psFree(sumC);
    … …  
    407409
    408410    if (cutIndex < 0) {
    409         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to converge to tolerance %g\n", tolerance);
     411        psError(PM_ERR_DATA, true, "Unable to converge to tolerance %g\n", tolerance);
    410412        psFree(ranking);
    411413        psFree(kernels);
    … …  
    482484    // Maintain photometric scaling
    483485    if (type == PM_SUBTRACTION_KERNEL_ISIS) {
    484         psKernel *subtract = kernels->preCalc->data[0]; // Kernel to subtract from the rest
     486
     487        // XXX in r26035, this code was just wrong.  we had:
     488
     489        // psKernel *subtract = kernels->preCalc->data[0]
     490
     491        // but, kernels->preCalc was an array of psArray, not an array of kernels.  It is now
     492        // an array of pmSubtractionKernelPreCalc.
     493
     494        pmSubtractionKernelPreCalc *subtract = kernels->preCalc->data[0]; // Kernel to subtract from the rest
     495
    485496        for (int i = 1; i < newSize; i++) {
    486497            if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) {
    487                 psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
    488                 psBinaryOp(kernel->image, kernel->image, "-", subtract->image);
     498                pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i]; // Kernel of interest
     499                psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "-", subtract->kernel->image);
    489500            }
    490501        }
    … …  
    495506        for (int i = 0; i < newSize; i++) {
    496507            if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) {
    497                 psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
    498                 kernel->kernel[0][0] -= 1.0;
     508                pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i]; // Kernel of interest
     509                preCalc->kernel->kernel[0][0] -= 1.0;
    499510            }
    500511        }
Note: See TracChangeset for help on using the changeset viewer.