IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 9, 2010, 8:27:57 AM (16 years ago)
Author:
eugene
Message:

separate penalties for kernels for image 1 and 2, convolve second moments with fwhm to get penalties; separate normalization windows for images 1 and 2; update visual to have facilities like psphot visual; curve of growth analysis for norm window

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionKernels.c

    r27335 r28866  
    2626    psFree(kernels->vStop);
    2727    psFree(kernels->preCalc);
    28     psFree(kernels->penalties);
     28    psFree(kernels->penalties1);
     29    psFree(kernels->penalties2);
    2930    psFree(kernels->solution1);
    3031    psFree(kernels->solution2);
    … …  
    140141    kernels->v = psVectorRealloc(kernels->v, start + numNew);
    141142    kernels->preCalc = psArrayRealloc(kernels->preCalc, start + numNew);
    142     kernels->penalties = psVectorRealloc(kernels->penalties, start + numNew);
     143
     144    kernels->penalties1 = psVectorRealloc(kernels->penalties1, start + numNew);
     145    kernels->penalties2 = psVectorRealloc(kernels->penalties2, start + numNew);
     146
    143147    kernels->inner = start;
    144148    kernels->num += numNew;
    … …  
    156160            kernels->v->data.S32[index] = v;
    157161            kernels->preCalc->data[index] = NULL;
    158             kernels->penalties->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v));
    159             psAssert (isfinite(kernels->penalties->data.F32[index]), "invalid penalty");
     162
     163            // XXX this needs to be changed to use the *convolved* second moment
     164            kernels->penalties1->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v));
     165            psAssert (isfinite(kernels->penalties1->data.F32[index]), "invalid penalty");
     166
     167            kernels->penalties2->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v));
     168            psAssert (isfinite(kernels->penalties2->data.F32[index]), "invalid penalty");
     169
    160170            psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v);
    161171        }
    … …  
    166176    kernels->v->n = start + numNew;
    167177    kernels->preCalc->n = start + numNew;
    168     kernels->penalties->n = start + numNew;
     178
     179    kernels->penalties1->n = start + numNew;
     180    kernels->penalties2->n = start + numNew;
    169181
    170182    return true;
    171183}
    172184
    173 bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc,
    174                                          int index, int size, int uOrder, int vOrder, float fwhm,
    175                                          bool AlardLuptonStyle, bool forceZeroNull)
     185static bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc,
     186                                                int index, int uOrder, int vOrder, float fwhm,
     187                                                bool AlardLuptonStyle, bool forceZeroNull)
    176188{
    177189    // we have 4 cases here:
    … …  
    182194
    183195    // Calculate moments
    184     double penalty = 0.0;                   // Moment, for penalty
    185196    double sum = 0.0, sum2 = 0.0;           // Sum of kernel component
    186197    float min = INFINITY, max = -INFINITY;  // Minimum and maximum kernel value
    187     for (int v = -size; v <= size; v++) {
    188         for (int u = -size; u <= size; u++) {
     198
     199    for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
     200        for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) {
    189201            double value = preCalc->kernel->kernel[v][u];
    190202            double value2 = PS_SQR(value);
    191203            sum += value;
    192204            sum2 += value2;
    193             penalty += value2 * PS_SQR((PS_SQR(u) + PS_SQR(v)));
    194205            min = PS_MIN(value, min);
    195206            max = PS_MAX(value, max);
    … …  
    198209
    199210#if 0
    200     fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf, moment: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max, penalty);
     211    fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max);
    201212#endif
    202213
    … …  
    239250
    240251    psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "*", psScalarAlloc(scale2D, PS_TYPE_F32));
    241     penalty *= 1.0 / sum2;
    242252
    243253    if (zeroNull) {
    … …  
    249259        double sum = 0.0;   // Sum of kernel component
    250260        float min = INFINITY, max = -INFINITY;  // Minimum and maximum kernel value
    251         for (int v = -size; v <= size; v++) {
    252             for (int u = -size; u <= size; u++) {
     261        for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
     262            for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) {
    253263                sum += preCalc->kernel->kernel[v][u];
    254264                min = PS_MIN(preCalc->kernel->kernel[v][u], min);
    … …  
    267277    }
    268278    kernels->preCalc->data[index] = preCalc;
    269     kernels->penalties->data.F32[index] = kernels->penalty * penalty;
    270     psAssert (isfinite(kernels->penalties->data.F32[index]), "invalid penalty");
    271     psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, fwhm, uOrder, vOrder, penalty);
     279    psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, fwhm, uOrder, vOrder);
    272280
    273281    return true;
    … …  
    321329
    322330                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values
    323                 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
    324                 // pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], false, false);
     331                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
     332                // pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], false, false);
    325333            }
    326334        }
    … …  
    379387            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    380388                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values
    381                 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
     389                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
    382390            }
    383391        }
    … …  
    385393            // XXX modify size for hermitians to account for sqrt(2) in Hermitian definition (relative to ISIS Gaussian)
    386394            pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS_RADIAL, order, order, size, sigma / sqrt(2.0)); // structure to hold precalculated values
    387             pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, order, order, fwhms->data.F32[i], true, true);
     395            pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, order, order, fwhms->data.F32[i], true, true);
    388396        }
    389397    }
    … …  
    437445            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    438446                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values
    439                 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
     447                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
    440448            }
    441449        }
    … …  
    506514
    507515                // XXX do we use Alard-Lupton normalization (last param true) or not?
    508                 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
     516                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
    509517
    510518                // XXXX test demo that deconvolved kernel is valid
    … …  
    572580    kernels->preCalc = psArrayAlloc(numBasisFunctions);
    573581    kernels->penalty = penalty;
    574     kernels->penalties = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
     582    kernels->penalties1 = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
     583    psVectorInit(kernels->penalties1, NAN);
     584    kernels->penalties2 = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
     585    psVectorInit(kernels->penalties2, NAN);
     586    kernels->havePenalties = false;
    575587    kernels->size = size;
    576588    kernels->inner = 0;
    … …  
    771783
    772784    psWarning("Kernel penalty for dual-convolution is not configured for SPAM kernels.");
    773     psVectorInit(kernels->penalties, 0.0);
     785    psVectorInit(kernels->penalties1, 0.0);
     786    psVectorInit(kernels->penalties2, 0.0);
    774787
    775788    return kernels;
    … …  
    866879
    867880    psWarning("Kernel penalty for dual-convolution is not configured for FRIES kernels.");
    868     psVectorInit(kernels->penalties, 0.0);
     881    psVectorInit(kernels->penalties1, 0.0);
     882    psVectorInit(kernels->penalties2, 0.0);
    869883
    870884    return kernels;
    … …  
    10401054                kernels->u->data.S32[index] = uOrder;
    10411055                kernels->v->data.S32[index] = vOrder;
    1042                 kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
    1043                 if (!isfinite(kernels->penalties->data.F32[index])) {
     1056
     1057                // XXX convert to use the convolved 2nd moment
     1058                kernels->penalties1->data.F32[index] = kernels->penalty * fabsf(moment);
     1059                if (!isfinite(kernels->penalties1->data.F32[index])) {
     1060                    psAbort ("invalid penalty");
     1061                }
     1062                kernels->penalties2->data.F32[index] = kernels->penalty * fabsf(moment);
     1063                if (!isfinite(kernels->penalties2->data.F32[index])) {
    10441064                    psAbort ("invalid penalty");
    10451065                }
    … …  
    12471267    out->preCalc = psMemIncrRefCounter(in->preCalc);
    12481268    out->penalty = in->penalty;
    1249     out->penalties = psMemIncrRefCounter(in->penalties);
     1269    out->penalties1 = psMemIncrRefCounter(in->penalties1);
     1270    out->penalties2 = psMemIncrRefCounter(in->penalties2);
    12501271    out->uStop = psMemIncrRefCounter(in->uStop);
    12511272    out->vStop = psMemIncrRefCounter(in->vStop);
Note: See TracChangeset for help on using the changeset viewer.