IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26587


Ignore:
Timestamp:
Jan 13, 2010, 4:06:30 PM (16 years ago)
Author:
Paul Price
Message:

Cleaning up penalty functions.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c

    r26574 r26587  
    9090
    9191    for (int i = 0, x = -size; x <= size; i++, x++) {
    92         float xf = x / sigma;
    93         float z = -0.25*xf*xf;
     92        float xf = x / sigma;
     93        float z = -0.25*xf*xf;
    9494        kernel->data.F32[i] = norm * p_pmSubtractionHermitianPolynomial(xf, order) * exp(z);
    9595    }
    … …  
    100100// Generate 1D convolution kernel for HERM (normalized for 2D)
    101101psKernel *pmSubtractionKernelHERM_RADIAL(float sigma, // Gaussian width
    102                                         int order, // Polynomial order
    103                                         int size // Kernel half-size
     102                                        int order, // Polynomial order
     103                                        int size // Kernel half-size
    104104    )
    105105{
    … …  
    112112    // generate 2D radial hermitian
    113113    for (int v = -size; v <= size; v++) {
    114         for (int u = -size; u <= size; u++) {
    115             float r = hypot(u, v) / sigma;
    116             float z = -0.25*r*r;
    117             kernel->kernel[v][u] = norm * p_pmSubtractionHermitianPolynomial(r, order) * exp(z);
    118         }
     114        for (int u = -size; u <= size; u++) {
     115            float r = hypot(u, v) / sigma;
     116            float z = -0.25*r*r;
     117            kernel->kernel[v][u] = norm * p_pmSubtractionHermitianPolynomial(r, order) * exp(z);
     118        }
    119119    }
    120120
    … …  
    156156            kernels->preCalc->data[index] = NULL;
    157157            kernels->penalties->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v));
    158             psAssert (isfinite(kernels->penalties->data.F32[index]), "invalid penalty");
     158            psAssert (isfinite(kernels->penalties->data.F32[index]), "invalid penalty");
    159159            psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v);
    160160        }
    … …  
    170170}
    171171
    172 bool pmSubtractionKernelPreCalcNormalize (pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc, int index, int size, int uOrder, int vOrder, float fwhm, bool AlardLuptonStyle, bool forceZeroNull) {
    173 
    174     // Calculate moments
    175     double moment = 0.0;    // Moment, for penalty
    176     for (int v = -size; v <= size; v++) {
    177         for (int u = -size; u <= size; u++) {
    178             double value = preCalc->kernel->kernel[v][u];
    179             moment += PS_SQR(value) * PS_SQR((PS_SQR(u) + PS_SQR(v)));
    180         }
    181     }
    182 
     172bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc,
     173                                         int index, int size, int uOrder, int vOrder, float fwhm,
     174                                         bool AlardLuptonStyle, bool forceZeroNull)
     175{
    183176    // we have 4 cases here:
    184177    // 1) for odd functions, normalize the kernel by the maximum swing / Npix
    … …  
    187180    // 4) for deconvolved hermitians, subtract 1 from the 0,0 pixel for the 0,0 function(s)
    188181
    189     double sum = 0.0;   // Sum of kernel component
    190     double min = FLT_MAX;
    191     double max = FLT_MIN;
    192 
     182    // Calculate moments
     183    double penalty = 0.0;                   // Moment, for penalty
     184    double sum = 0.0, sum2 = 0.0;           // Sum of kernel component
     185    float min = INFINITY, max = -INFINITY;  // Minimum and maximum kernel value
    193186    for (int v = -size; v <= size; v++) {
    194         for (int u = -size; u <= size; u++) {
    195             sum += preCalc->kernel->kernel[v][u];
    196             min = PS_MIN(preCalc->kernel->kernel[v][u], min);
    197             max = PS_MAX(preCalc->kernel->kernel[v][u], max);
    198         }
    199     }
     187        for (int u = -size; u <= size; u++) {
     188            double value = preCalc->kernel->kernel[v][u];
     189            double value2 = PS_SQR(value);
     190            sum += value;
     191            sum2 += value2;
     192            penalty += value2 * PS_SQR((PS_SQR(u) + PS_SQR(v)));
     193            min = PS_MIN(value, min);
     194            max = PS_MAX(value, max);
     195        }
     196    }
     197
    200198#if 0
    201     fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf, moment: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max, moment);
     199    fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf, moment: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max, penalty);
    202200#endif
    203201
    204     // only even terms have non-zero sums
    205     if ((uOrder % 2 == 0) && (vOrder % 2 == 0)) {
    206         moment /= PS_SQR(sum);
    207     } else {
    208         // XXX keep this?
    209         moment = 0.0;
    210     }
    211 
    212     bool zeroNull = false;
    213     float scale1D = 1.0 / sqrt(sum);
    214     float scale2D = 1.0 / sum;
    215 
    216     if (AlardLuptonStyle && (uOrder % 2 == 0 && vOrder % 2 == 0)) {
    217         zeroNull = true;
    218     }
     202    bool zeroNull = false;              // Zero out using the null position?
     203    float scale2D = NAN;                // Scaling for 2-D kernels
     204
     205    if (AlardLuptonStyle) {
     206        if (uOrder % 2 == 0 && vOrder % 2 == 0) {
     207            // Even functions: normalise to unit sum and subtract null pixel so that sum is zero
     208            scale2D = 1.0 / fabs(sum);
     209            zeroNull = true;
     210        } else {
     211            // Odd functions: choose normalisation so that parameters have about the same strength as for even
     212            // functions, no subtraction of null pixel because the sum is already (near) zero
     213            scale2D = 1.0 / max;
     214            zeroNull = false;
     215        }
     216    }
     217
    219218    if (!AlardLuptonStyle && (uOrder == 0 && vOrder == 0)) {
    220         zeroNull = true;
     219        zeroNull = true;
    221220    }
    222221    if (forceZeroNull) {
    223         zeroNull = true;
     222        zeroNull = true;
    224223    }
    225224    if (!forceZeroNull && ((uOrder % 2) || (vOrder % 2))) {
    226         // scale2D = 1.0 / (preCalc->kernel->image->numCols * preCalc->kernel->image->numRows * max);
    227         scale2D = 1.0 / max;
    228         scale1D = sqrt(scale2D);
    229     }
    230 
    231     if (preCalc->xKernel) { 
    232         psBinaryOp(preCalc->xKernel, preCalc->xKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32));
     225        // Odd function
     226        scale2D = 1.0 / max;
     227    }
     228
     229    float scale1D = sqrtf(scale2D);     // Scaling for 1-D kernels
     230    if (preCalc->xKernel) {
     231        psBinaryOp(preCalc->xKernel, preCalc->xKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32));
    233232    }
    234233    if (preCalc->yKernel) {
    235         psBinaryOp(preCalc->yKernel, preCalc->yKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32));
    236     }
     234        psBinaryOp(preCalc->yKernel, preCalc->yKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32));
     235    }
     236
    237237    psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "*", psScalarAlloc(scale2D, PS_TYPE_F32));
     238    penalty *= PS_SQR(scale2D);
     239
    238240    if (zeroNull) {
    239         preCalc->kernel->kernel[0][0] -= 1.0;
     241        preCalc->kernel->kernel[0][0] -= 1.0;
    240242    }
    241243
    … …  
    245247    max = FLT_MIN;
    246248    for (int v = -size; v <= size; v++) {
    247         for (int u = -size; u <= size; u++) {
    248             sum += preCalc->kernel->kernel[v][u];
    249             min = PS_MIN(preCalc->kernel->kernel[v][u], min);
    250             max = PS_MAX(preCalc->kernel->kernel[v][u], max);
    251         }
     249        for (int u = -size; u <= size; u++) {
     250            sum += preCalc->kernel->kernel[v][u];
     251            min = PS_MIN(preCalc->kernel->kernel[v][u], min);
     252            max = PS_MAX(preCalc->kernel->kernel[v][u], max);
     253        }
    252254    }
    253255    fprintf(stderr, "%d mod: %lf, null: %f, min: %lf, max: %lf, scale: %f\n", index, sum, preCalc->kernel->kernel[0][0], min, max, scale2D);
    … …  
    258260    kernels->v->data.S32[index] = vOrder;
    259261    if (kernels->preCalc->data[index]) {
    260         psFree(kernels->preCalc->data[index]);
     262        psFree(kernels->preCalc->data[index]);
    261263    }
    262264    kernels->preCalc->data[index] = preCalc;
    263     kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
     265    kernels->penalties->data.F32[index] = kernels->penalty * penalty;
    264266    psAssert (isfinite(kernels->penalties->data.F32[index]), "invalid penalty");
    265     psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, fwhm, uOrder, vOrder, fabsf(moment));
     267    psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, fwhm, uOrder, vOrder, penalty);
    266268
    267269    return true;
    … …  
    284286    psVector *orders = psVectorAllocEmpty (ordersIN->n, PS_TYPE_S32);
    285287    for (int i = 0; i < fwhmsIN->n; i++) {
    286         if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
    287         psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
    288         psVectorAppend(orders, ordersIN->data.S32[i]);
     288        if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
     289        psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
     290        psVectorAppend(orders, ordersIN->data.S32[i]);
    289291    }
    290292
    … …  
    315317
    316318                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values
    317                 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
    318                 // pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], false, false);
     319                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
     320                // pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], false, false);
    319321            }
    320322        }
    … …  
    325327
    326328pmSubtractionKernels *pmSubtractionKernelsISIS_RADIAL(int size, int spatialOrder,
    327                                                       const psVector *fwhmsIN, const psVector *ordersIN,
    328                                                       float penalty, pmSubtractionMode mode)
     329                                                      const psVector *fwhmsIN, const psVector *ordersIN,
     330                                                      float penalty, pmSubtractionMode mode)
    329331{
    330332    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
    … …  
    340342    psVector *orders = psVectorAllocEmpty (ordersIN->n, PS_TYPE_S32);
    341343    for (int i = 0; i < fwhmsIN->n; i++) {
    342         if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
    343         psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
    344         psVectorAppend(orders, ordersIN->data.S32[i]);
     344        if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
     345        psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
     346        psVectorAppend(orders, ordersIN->data.S32[i]);
    345347    }
    346348
    … …  
    353355        psStringAppend(&params, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]);
    354356        num += (gaussOrder + 1) * (gaussOrder + 2) / 2;
    355         num += (11 - gaussOrder - 1);   // include all higher order radial terms
     357        num += (11 - gaussOrder - 1);   // include all higher order radial terms
    356358    }
    357359
    … …  
    369371            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    370372                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values
    371                 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
     373                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
    372374            }
    373375        }
    374         for (int order = orders->data.S32[i] + 1; order < 11; order ++, index ++) {
    375             // XXX modify size for hermitians to account for sqrt(2) in Hermitian definition (relative to ISIS Gaussian)
    376             pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS_RADIAL, order, order, size, sigma / sqrt(2.0)); // structure to hold precalculated values
    377             pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, order, order, fwhms->data.F32[i], true, true);
    378         }
     376        for (int order = orders->data.S32[i] + 1; order < 11; order ++, index ++) {
     377            // XXX modify size for hermitians to account for sqrt(2) in Hermitian definition (relative to ISIS Gaussian)
     378            pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS_RADIAL, order, order, size, sigma / sqrt(2.0)); // structure to hold precalculated values
     379            pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, order, order, fwhms->data.F32[i], true, true);
     380        }
    379381    }
    380382    return kernels;
    … …  
    382384
    383385pmSubtractionKernels *pmSubtractionKernelsHERM(int size, int spatialOrder,
    384                                                const psVector *fwhmsIN, const psVector *ordersIN,
    385                                                float penalty, pmSubtractionMode mode)
     386                                               const psVector *fwhmsIN, const psVector *ordersIN,
     387                                               float penalty, pmSubtractionMode mode)
    386388{
    387389    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
    … …  
    397399    psVector *orders = psVectorAllocEmpty (ordersIN->n, PS_TYPE_S32);
    398400    for (int i = 0; i < fwhmsIN->n; i++) {
    399         if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
    400         psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
    401         psVectorAppend(orders, ordersIN->data.S32[i]);
     401        if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
     402        psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
     403        psVectorAppend(orders, ordersIN->data.S32[i]);
    402404    }
    403405
    … …  
    425427            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    426428                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values
    427                 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
     429                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
    428430            }
    429431        }
    … …  
    434436
    435437pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, int spatialOrder,
    436                                                      const psVector *fwhmsIN, const psVector *ordersIN,
    437                                                      float penalty, pmSubtractionMode mode)
     438                                                     const psVector *fwhmsIN, const psVector *ordersIN,
     439                                                     float penalty, pmSubtractionMode mode)
    438440{
    439441    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
    … …  
    449451    psVector *orders = psVectorAllocEmpty (ordersIN->n, PS_TYPE_S32);
    450452    for (int i = 0; i < fwhmsIN->n; i++) {
    451         if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
    452         psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
    453         psVectorAppend(orders, ordersIN->data.S32[i]);
     453        if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
     454        psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
     455        psVectorAppend(orders, ordersIN->data.S32[i]);
    454456    }
    455457
    … …  
    488490                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values
    489491
    490                 // save the generated 2D kernel as the target, deconvolve it by Gaussian, replacing the generated 2D kernel
    491                 psKernel *kernelTarget = preCalc->kernel;
     492                // save the generated 2D kernel as the target, deconvolve it by Gaussian, replacing the generated 2D kernel
     493                psKernel *kernelTarget = preCalc->kernel;
    492494                preCalc->kernel = pmSubtractionDeconvolveKernel(kernelTarget, kernelGauss); // Kernel
    493495
    494                 // XXX do we use Alard-Lupton normalization (last param true) or not?
    495                 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
    496 
    497                 // XXXX test demo that deconvolved kernel is valid
     496                // XXX do we use Alard-Lupton normalization (last param true) or not?
     497                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
     498
     499                // XXXX test demo that deconvolved kernel is valid
    498500# if 1
    499                 psImage *kernelConv = psImageConvolveFFT(NULL, preCalc->kernel->image, NULL, 0, kernelGauss);
    500                 psArrayAdd (deconKernels, 100, kernelConv);
    501                 psFree (kernelConv);
    502 
    503                 if (!uOrder && !vOrder){
    504                     pmSubtractionVisualShowSubtraction (kernelTarget->image, preCalc->kernel->image, kernelConv);
    505                 }
     501                psImage *kernelConv = psImageConvolveFFT(NULL, preCalc->kernel->image, NULL, 0, kernelGauss);
     502                psArrayAdd (deconKernels, 100, kernelConv);
     503                psFree (kernelConv);
     504
     505                if (!uOrder && !vOrder){
     506                    pmSubtractionVisualShowSubtraction (kernelTarget->image, preCalc->kernel->image, kernelConv);
     507                }
    506508# endif
    507509            }
    … …  
    512514    psImage *dot = psImageAlloc(deconKernels->n, deconKernels->n, PS_TYPE_F32);
    513515    for (int i = 0; i < deconKernels->n; i++) {
    514         for (int j = 0; j <= i; j++) {
    515             psImage *t1 = deconKernels->data[i];
    516             psImage *t2 = deconKernels->data[j];
    517 
    518             double sum = 0.0;
    519             for (int iy = 0; iy < t1->numRows; iy++) {
    520                 for (int ix = 0; ix < t1->numCols; ix++) {
    521                     sum += t1->data.F32[iy][ix] * t2->data.F32[iy][ix];
    522                 }
    523             }
    524             dot->data.F32[j][i] = sum;
    525             dot->data.F32[i][j] = sum;
    526         }
     516        for (int j = 0; j <= i; j++) {
     517            psImage *t1 = deconKernels->data[i];
     518            psImage *t2 = deconKernels->data[j];
     519
     520            double sum = 0.0;
     521            for (int iy = 0; iy < t1->numRows; iy++) {
     522                for (int ix = 0; ix < t1->numCols; ix++) {
     523                    sum += t1->data.F32[iy][ix] * t2->data.F32[iy][ix];
     524                }
     525            }
     526            dot->data.F32[j][i] = sum;
     527            dot->data.F32[i][j] = sum;
     528        }
    527529    }
    528530    pmSubtractionVisualShowSubtraction (dot, NULL, NULL);
    … …  
    587589    switch (type) {
    588590      case PM_SUBTRACTION_KERNEL_ISIS:
    589         preCalc->xKernel = pmSubtractionKernelISIS(sigma, uOrder, size);
    590         preCalc->yKernel = pmSubtractionKernelISIS(sigma, vOrder, size);
    591         preCalc->uCoords = NULL;
    592         preCalc->vCoords = NULL;
    593         preCalc->poly    = NULL;
    594         break;
     591        preCalc->xKernel = pmSubtractionKernelISIS(sigma, uOrder, size);
     592        preCalc->yKernel = pmSubtractionKernelISIS(sigma, vOrder, size);
     593        preCalc->uCoords = NULL;
     594        preCalc->vCoords = NULL;
     595        preCalc->poly    = NULL;
     596        break;
    595597      case PM_SUBTRACTION_KERNEL_HERM:
    596         preCalc->xKernel = pmSubtractionKernelHERM(sigma, uOrder, size);
    597         preCalc->yKernel = pmSubtractionKernelHERM(sigma, vOrder, size);
    598         preCalc->uCoords = NULL;
    599         preCalc->vCoords = NULL;
    600         preCalc->poly    = NULL;
    601         break;
     598        preCalc->xKernel = pmSubtractionKernelHERM(sigma, uOrder, size);
     599        preCalc->yKernel = pmSubtractionKernelHERM(sigma, vOrder, size);
     600        preCalc->uCoords = NULL;
     601        preCalc->vCoords = NULL;
     602        preCalc->poly    = NULL;
     603        break;
    602604      case PM_SUBTRACTION_KERNEL_RINGS:
    603         // the RINGS kernel uses the uCoords, vCoords, and poly elements of the structure
    604         // we allocate these vectors here, but leave the kernel generation to the main function
    605         preCalc->xKernel = NULL;
    606         preCalc->yKernel = NULL;
    607         preCalc->kernel  = NULL;
    608         preCalc->uCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // u coords
    609         preCalc->vCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // v coords
    610         preCalc->poly    = psVectorAllocEmpty(size, PS_TYPE_F32); // Polynomial
    611         return preCalc;
     605        // the RINGS kernel uses the uCoords, vCoords, and poly elements of the structure
     606        // we allocate these vectors here, but leave the kernel generation to the main function
     607        preCalc->xKernel = NULL;
     608        preCalc->yKernel = NULL;
     609        preCalc->kernel  = NULL;
     610        preCalc->uCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // u coords
     611        preCalc->vCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // v coords
     612        preCalc->poly    = psVectorAllocEmpty(size, PS_TYPE_F32); // Polynomial
     613        return preCalc;
    612614      case PM_SUBTRACTION_KERNEL_ISIS_RADIAL:
    613         preCalc->kernel  = pmSubtractionKernelHERM_RADIAL(sigma, uOrder, size);
    614         preCalc->xKernel = NULL;
    615         preCalc->yKernel = NULL;
    616         preCalc->uCoords = NULL;
    617         preCalc->vCoords = NULL;
    618         preCalc->poly    = NULL;
    619         return preCalc;
     615        preCalc->kernel  = pmSubtractionKernelHERM_RADIAL(sigma, uOrder, size);
     616        preCalc->xKernel = NULL;
     617        preCalc->yKernel = NULL;
     618        preCalc->uCoords = NULL;
     619        preCalc->vCoords = NULL;
     620        preCalc->poly    = NULL;
     621        return preCalc;
    620622      default:
    621         psAbort("programming error: invalid type for PreCalc kernel");
     623        psAbort("programming error: invalid type for PreCalc kernel");
    622624    }
    623625
    … …  
    626628    // generate 2D kernel from 1D realizations
    627629    for (int v = -size, y = 0; v <= size; v++, y++) {
    628         for (int u = -size, x = 0; u <= size; u++, x++) {
    629             preCalc->kernel->kernel[v][u] = preCalc->xKernel->data.F32[x] * preCalc->yKernel->data.F32[y]; // Value of kernel
    630         }
    631     }
    632    
     630        for (int u = -size, x = 0; u <= size; u++, x++) {
     631            preCalc->kernel->kernel[v][u] = preCalc->xKernel->data.F32[x] * preCalc->yKernel->data.F32[y]; // Value of kernel
     632        }
     633    }
     634
    633635    return preCalc;
    634636}
    … …  
    958960            for (int vOrder = 0; vOrder <= (i == 0 ? 0 : ringsOrder - uOrder); vOrder++, index++) {
    959961
    960                 pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc (PM_SUBTRACTION_KERNEL_RINGS, 0, 0, RINGS_BUFFER, 0.0);
     962                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc (PM_SUBTRACTION_KERNEL_RINGS, 0, 0, RINGS_BUFFER, 0.0);
    961963                double moment = 0.0;    // Moment, for penalty
    962964
    … …  
    964966                    // Central pixel is easy
    965967                    preCalc->uCoords->data.S32[0] = 0;
    966                     preCalc->vCoords->data.S32[0] = 0;
     968                    preCalc->vCoords->data.S32[0] = 0;
    967969                    preCalc->poly->data.F32[0] = 1.0;
    968970                    preCalc->uCoords->n = 1;
    969                     preCalc->vCoords->n = 1;
    970                     preCalc->poly->n = 1;
     971                    preCalc->vCoords->n = 1;
     972                    preCalc->poly->n = 1;
    971973                    radiusLast = 0;
    972974                    moment = 0.0;
    … …  
    10251027                kernels->v->data.S32[index] = vOrder;
    10261028                kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
    1027                 if (!isfinite(kernels->penalties->data.F32[index])) {
    1028                     psAbort ("invalid penalty");
    1029                 }
     1029                if (!isfinite(kernels->penalties->data.F32[index])) {
     1030                    psAbort ("invalid penalty");
     1031                }
    10301032
    10311033                psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d\n", index,
    … …  
    11171119
    11181120    // currently known descriptions:
    1119     // ISIS(...), ISIS_RADIAL(...), HERM(...), DECONV_HERM(...), POIS(...), SPAM(...), 
    1120     // FRIES(...), GUNK=ISIS(...)+POIS(...), RINGS(...), 
     1121    // ISIS(...), ISIS_RADIAL(...), HERM(...), DECONV_HERM(...), POIS(...), SPAM(...),
     1122    // FRIES(...), GUNK=ISIS(...)+POIS(...), RINGS(...),
    11211123    // the descriptive name is the set of characters before the (
    1122    
     1124
    11231125    type = pmSubtractionKernelsTypeFromString (description);
    11241126    char *ptr = strchr(description, '(');
    … …  
    11301132      case PM_SUBTRACTION_KERNEL_HERM:
    11311133      case PM_SUBTRACTION_KERNEL_DECONV_HERM:
    1132         PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt);
    1133 
    1134         // Count the number of Gaussians
    1135         int numGauss = 0;
    1136         for (char *string = ptr; string; string = strchr(string + 1, '(')) {
    1137             numGauss++;
    1138         }
    1139 
    1140         fwhms = psVectorAlloc(numGauss, PS_TYPE_F32);
    1141         orders = psVectorAlloc(numGauss, PS_TYPE_S32);
    1142 
    1143         for (int i = 0; i < numGauss; i++) {
    1144             ptr++;                                                              // Eat the '('
    1145             PARSE_STRING_NUMBER(fwhms->data.F32[i], ptr, ',', parseStringFloat); // Eat "1.234,"
    1146             PARSE_STRING_NUMBER(orders->data.S32[i], ptr, ')', parseStringInt); // Eat "3)"
    1147         }
    1148 
    1149         ptr++;                      // Eat ','
    1150         PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
    1151         penalty = parseStringFloat(ptr);
    1152 
    1153         return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
     1134        PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt);
     1135
     1136        // Count the number of Gaussians
     1137        int numGauss = 0;
     1138        for (char *string = ptr; string; string = strchr(string + 1, '(')) {
     1139            numGauss++;
     1140        }
     1141
     1142        fwhms = psVectorAlloc(numGauss, PS_TYPE_F32);
     1143        orders = psVectorAlloc(numGauss, PS_TYPE_S32);
     1144
     1145        for (int i = 0; i < numGauss; i++) {
     1146            ptr++;                                                              // Eat the '('
     1147            PARSE_STRING_NUMBER(fwhms->data.F32[i], ptr, ',', parseStringFloat); // Eat "1.234,"
     1148            PARSE_STRING_NUMBER(orders->data.S32[i], ptr, ')', parseStringInt); // Eat "3)"
     1149        }
     1150
     1151        ptr++;                      // Eat ','
     1152        PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
     1153        penalty = parseStringFloat(ptr);
     1154
     1155        return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
    11541156
    11551157      case PM_SUBTRACTION_KERNEL_RINGS:
    … …  
    11591161        PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
    11601162        PARSE_STRING_NUMBER(penalty, ptr, ')', parseStringInt);
    1161         return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
     1163        return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
    11621164      default:
    1163         psAbort("Deciphering kernels other than ISIS, HERM, DECONV_HERM or RINGS is not currently supported.");
    1164     } 
     1165        psAbort("Deciphering kernels other than ISIS, HERM, DECONV_HERM or RINGS is not currently supported.");
     1166    }
    11651167    return NULL;
    11661168}
    … …  
    11771179    char *ptr = strchr(type, '(');
    11781180    if (ptr) {
    1179         nameLength = ptr - type;
     1181        nameLength = ptr - type;
    11801182    }
    11811183
Note: See TracChangeset for help on using the changeset viewer.