IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 26, 2010, 9:18:39 AM (16 years ago)
Author:
Serge CHASTEL
Message:

Merging trunk in branch

Location:
branches/sc_branches/trunkTest
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/sc_branches/trunkTest

  • branches/sc_branches/trunkTest/psModules

    • Property svn:mergeinfo deleted

  • branches/sc_branches/trunkTest/psModules/src/imcombine/pmSubtractionKernels.c

    r27335 r29060  
    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
    … …  
    207218        if (uOrder % 2 == 0 && vOrder % 2 == 0) {
    208219            // Even functions: normalise to unit sum and subtract null pixel so that sum is zero
    209             scale2D = 1.0 / fabs(sum);
     220            // Re-normalize
     221            // scale2D  = 1.0 / fabs(sum);
     222            scale2D  = 1.0 / sqrt(sum2);
    210223            zeroNull = true;
    211224        } else {
    … …  
    239252
    240253    psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "*", psScalarAlloc(scale2D, PS_TYPE_F32));
    241     penalty *= 1.0 / sum2;
    242254
    243255    if (zeroNull) {
    244         preCalc->kernel->kernel[0][0] -= 1.0;
    245     }
    246 
    247 #if 0
     256        // preCalc->kernel->kernel[0][0] -= 1.0;
     257        preCalc->kernel->kernel[0][0] -= sum / sqrt (sum2);
     258    }
     259
     260#if 1
    248261    {
    249         double sum = 0.0;   // Sum of kernel component
     262        double Sum = 0.0;   // Sum of kernel component
     263        double Sum2 = 0.0;   // Sum of kernel component
    250264        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++) {
    253                 sum += preCalc->kernel->kernel[v][u];
     265        for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
     266            for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) {
     267                double value = preCalc->kernel->kernel[v][u];
     268                Sum += value;
     269                Sum2 += PS_SQR(value);
    254270                min = PS_MIN(preCalc->kernel->kernel[v][u], min);
    255271                max = PS_MAX(preCalc->kernel->kernel[v][u], max);
    256272            }
    257273        }
    258         fprintf(stderr, "%d mod: %lf, null: %f, min: %lf, max: %lf, scale: %f\n", index, sum, preCalc->kernel->kernel[0][0], min, max, scale2D);
     274        fprintf(stderr, "%d sum: %lf, sum2: %lf, null: %f, min: %lf, max: %lf, scale: %f\n", index, Sum, Sum2, preCalc->kernel->kernel[0][0], min, max, scale2D);
    259275    }
    260276#endif
    … …  
    267283    }
    268284    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);
     285    psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, fwhm, uOrder, vOrder);
    272286
    273287    return true;
    … …  
    321335
    322336                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);
     337                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
     338                // pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], false, false);
    325339            }
    326340        }
    … …  
    379393            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    380394                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);
     395                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
    382396            }
    383397        }
    … …  
    385399            // XXX modify size for hermitians to account for sqrt(2) in Hermitian definition (relative to ISIS Gaussian)
    386400            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);
     401            pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, order, order, fwhms->data.F32[i], true, true);
    388402        }
    389403    }
    … …  
    437451            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    438452                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);
     453                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
    440454            }
    441455        }
    … …  
    506520
    507521                // 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);
     522                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
    509523
    510524                // XXXX test demo that deconvolved kernel is valid
    … …  
    572586    kernels->preCalc = psArrayAlloc(numBasisFunctions);
    573587    kernels->penalty = penalty;
    574     kernels->penalties = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
     588    kernels->penalties1 = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
     589    psVectorInit(kernels->penalties1, NAN);
     590    kernels->penalties2 = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
     591    psVectorInit(kernels->penalties2, NAN);
     592    kernels->havePenalties = false;
    575593    kernels->size = size;
    576594    kernels->inner = 0;
    … …  
    771789
    772790    psWarning("Kernel penalty for dual-convolution is not configured for SPAM kernels.");
    773     psVectorInit(kernels->penalties, 0.0);
     791    psVectorInit(kernels->penalties1, 0.0);
     792    psVectorInit(kernels->penalties2, 0.0);
    774793
    775794    return kernels;
    … …  
    866885
    867886    psWarning("Kernel penalty for dual-convolution is not configured for FRIES kernels.");
    868     psVectorInit(kernels->penalties, 0.0);
     887    psVectorInit(kernels->penalties1, 0.0);
     888    psVectorInit(kernels->penalties2, 0.0);
    869889
    870890    return kernels;
    … …  
    10401060                kernels->u->data.S32[index] = uOrder;
    10411061                kernels->v->data.S32[index] = vOrder;
    1042                 kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
    1043                 if (!isfinite(kernels->penalties->data.F32[index])) {
     1062
     1063                // XXX convert to use the convolved 2nd moment
     1064                kernels->penalties1->data.F32[index] = kernels->penalty * fabsf(moment);
     1065                if (!isfinite(kernels->penalties1->data.F32[index])) {
     1066                    psAbort ("invalid penalty");
     1067                }
     1068                kernels->penalties2->data.F32[index] = kernels->penalty * fabsf(moment);
     1069                if (!isfinite(kernels->penalties2->data.F32[index])) {
    10441070                    psAbort ("invalid penalty");
    10451071                }
    … …  
    12471273    out->preCalc = psMemIncrRefCounter(in->preCalc);
    12481274    out->penalty = in->penalty;
    1249     out->penalties = psMemIncrRefCounter(in->penalties);
     1275    out->penalties1 = psMemIncrRefCounter(in->penalties1);
     1276    out->penalties2 = psMemIncrRefCounter(in->penalties2);
    12501277    out->uStop = psMemIncrRefCounter(in->uStop);
    12511278    out->vStop = psMemIncrRefCounter(in->vStop);
Note: See TracChangeset for help on using the changeset viewer.