IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 29170


Ignore:
Timestamp:
Sep 17, 2010, 11:29:31 AM (16 years ago)
Author:
eugene
Message:

rework the dual convolution paradigm which is used to prevent excessively large kernels: replace penalties with an explicit coupled minimization of the second moment of the square of the output psf; the derivatives of this term in the minimization have the form of second moments of the various kernel and image combinations; to ensure the relative scaling is invarient between different flux normalizations and to ensure the second moments are not a function of the source fluxes, it is necessary to change where the normalization is applied

Location:
branches/eam_branches/ipp-20100823/psModules/src/imcombine
Files:
2 added
5 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.c

    r29165 r29170  
    2121//#define USE_WINDOW                      // Include weight (1/variance) in equation?
    2222
     23# define PENALTY false
     24# define MOMENTS (!PENALTY)
    2325
    2426//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    … …  
    170172
    171173// Calculate the least-squares matrix and vector for dual convolution
    172 static bool calculateDualMatrixVector(psImage *matrix,        // Least-squares matrix, updated
    173                                       psVector *vector,       // Least-squares vector, updated
     174// XXX we could avoid calculating these values on successive passes *if* the stamp has not changed.
     175static bool calculateDualMatrixVector(pmSubtractionStamp *stamp,              // stamp of interest
    174176                                      double normValue,       // Normalisation, updated
    175                                       const psKernel *image1, // Image 1
    176                                       const psKernel *image2, // Image 2
     177                                      double normValue2,      // Normalisation, updated
    177178                                      const psKernel *weight,  // Weight image
    178179                                      const psKernel *window,  // Window image
    179                                       const psArray *convolutions1, // Convolutions of image 1 for each kernel
    180                                       const psArray *convolutions2, // Convolutions of image 2 for each kernel
    181180                                      const pmSubtractionKernels *kernels, // Kernels
    182181                                      const psImage *polyValues, // Spatial polynomial values
    … …  
    184183                                      )
    185184{
    186     int numKernels = kernels->num;                      // Number of kernels
    187     int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
     185    int numKernels = kernels->num;                        // Number of kernels
     186    int spatialOrder = kernels->spatialOrder;             // Order of spatial variation
    188187    int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
    189188    double poly[numPoly];                                 // Polynomial terms
    190189    double poly2[numPoly][numPoly];                       // Polynomial-polynomial values
    191190
    192     int numParams = numKernels * numPoly;  // Number of regular parameters
    193     int numParams2 = numKernels * numPoly; // Number of additional parameters for dual
    194     int numDual = numParams + numParams2;  // Total number of parameters for dual
     191    int numParams = numKernels * numPoly;                 // Number of regular parameters
     192    int numParams2 = numKernels * numPoly;         // Number of additional parameters for dual
     193    int numDual = numParams + numParams2;          // Total number of parameters for dual
     194
     195    psImage *matrix = stamp->matrix;               // Least-squares matrix, updated
     196    psVector *vector = stamp->vector;              // Least-squares vector, updated
     197    psKernel *image1 = stamp->image1;              // Image 1
     198    psKernel *image2 = stamp->image2;              // Image 2
     199    psArray *convolutions1 = stamp->convolutions1; // Convolutions of image 1 for each kernel
     200    psArray *convolutions2 = stamp->convolutions2; // Convolutions of image 2 for each kernel
    195201
    196202    psAssert(matrix &&
    … …  
    199205             matrix->numRows == numDual,
    200206             "Least-squares matrix is bad.");
     207
    201208    psAssert(vector &&
    202209             vector->type.type == PS_TYPE_F64 &&
    … …  
    232239    // [kernel 0, x^n y^m][kernel 1 x^n y^m]...[kernel N, x^n y^m]
    233240
     241    // for DUAL convolution analysis, we apply the normalization to I1 as follows:
     242    // norm = I2 / I1
     243    //
     244    // I1c = norm I1 + \Sum_i a_i norm I1 \cross k_i
     245    // I2c =      I2 + \Sum_i b_i      I2 \cross k_i
     246
     247    // we cannot absorb the normalization into a_i until the analysis is complete, or the
     248    // second moment terms are incorrectly calculated.
     249
    234250    for (int i = 0; i < numKernels; i++) {
    235251        psKernel *iConv1 = convolutions1->data[i]; // Convolution 1 for index i
    … …  
    242258            double sumBB = 0.0;         // Sum of convolution products between image 2
    243259            double sumAB = 0.0;         // Sum of convolution products across images 1 and 2
     260
     261            double MxxAA = 0.0;
     262            double MyyAA = 0.0;
     263            double MxxBB = 0.0;
     264            double MyyBB = 0.0;
    244265            for (int y = - footprint; y <= footprint; y++) {
    245266                for (int x = - footprint; x <= footprint; x++) {
    246                     double aa = iConv1->kernel[y][x] * jConv1->kernel[y][x];
     267                    double aa = iConv1->kernel[y][x] * jConv1->kernel[y][x] * PS_SQR(normValue);
    247268                    double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x];
    248                     double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x];
     269                    double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x] * normValue;
    249270                    if (weight) {
    250271                        float wtVal = weight->kernel[y][x];
    … …  
    262283                    sumBB += bb;
    263284                    sumAB += ab;
    264                 }
    265             }
     285
     286                    if (MOMENTS) {
     287                        MxxAA += x*x*aa;
     288                        MyyAA += y*y*aa;
     289                        MxxBB += x*x*bb;
     290                        MyyBB += y*y*bb;
     291                    }
     292                }
     293            }
     294
     295            if (MOMENTS) {
     296                MxxAA /= stamp->normSquare1 * PS_SQR(normValue);
     297                MyyAA /= stamp->normSquare1 * PS_SQR(normValue);
     298                MxxBB /= stamp->normSquare2;
     299                MyyBB /= stamp->normSquare2;
     300            }
     301
     302            // XXX this makes the Chisq portion independent of the normalization and star flux
     303            // but may be mis-scaling between stars of different fluxes
     304            sumAA /= PS_SQR(stamp->normI2);
     305            sumAB /= PS_SQR(stamp->normI2);
     306            sumBB /= PS_SQR(stamp->normI2);
     307
     308            // fprintf (stderr, "i,j : %d %d : M(xx,yy)(AA,BB) : %f %f %f %f\n", i, j, MxxAA, MyyAA, MxxBB, MyyBB);
    266309
    267310            // Spatial variation of kernel coeffs
    … …  
    278321                    matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb;
    279322
     323                    // add in second moments
     324                    if (MOMENTS) {
     325                        matrix->data.F64[iIndex][jIndex] += kernels->penalty * MxxAA;
     326                        matrix->data.F64[iIndex][jIndex] += kernels->penalty * MyyAA;
     327
     328                        matrix->data.F64[jIndex][iIndex] += kernels->penalty * MxxAA;
     329                        matrix->data.F64[jIndex][iIndex] += kernels->penalty * MyyAA;
     330
     331                        matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MxxBB;
     332                        matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MyyBB;
     333
     334                        matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MxxBB;
     335                        matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MyyBB;
     336
     337                        // XXX this uses the single moments, first try
     338                        // matrix->data.F64[iIndex][jIndex] += kernels->penalty * stamp->MxxI1->data.F32[i]*stamp->MxxI1->data.F32[j];
     339                        // matrix->data.F64[iIndex][jIndex] += kernels->penalty * stamp->MyyI1->data.F32[i]*stamp->MyyI1->data.F32[j];
     340                        //
     341                        // matrix->data.F64[jIndex][iIndex] += kernels->penalty * stamp->MxxI1->data.F32[i]*stamp->MxxI1->data.F32[j];
     342                        // matrix->data.F64[jIndex][iIndex] += kernels->penalty * stamp->MyyI1->data.F32[i]*stamp->MyyI1->data.F32[j];
     343                        //
     344                        // matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * stamp->MxxI2->data.F32[i]*stamp->MxxI2->data.F32[j];
     345                        // matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * stamp->MyyI2->data.F32[i]*stamp->MyyI2->data.F32[j];
     346                        //
     347                        // matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * stamp->MxxI2->data.F32[i]*stamp->MxxI2->data.F32[j];
     348                        // matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * stamp->MyyI2->data.F32[i]*stamp->MyyI2->data.F32[j];
     349                    }
    280350                    matrix->data.F64[iIndex][jIndex + numParams] = ab;
    281351                    matrix->data.F64[jIndex + numParams][iIndex] = ab;
    … …  
    290360            for (int y = - footprint; y <= footprint; y++) {
    291361                for (int x = - footprint; x <= footprint; x++) {
    292                     double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x];
     362                    double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x] * normValue;
    293363                    if (weight) {
    294364                        ab *= weight->kernel[y][x];
    … …  
    301371            }
    302372
     373            // XXX this makes the Chisq portion independent of the normalization and star flux
     374            // but may be mis-scaling between stars of different fluxes
     375            sumAB /= PS_SQR(stamp->normI2);
     376
    303377            // Spatial variation of kernel coeffs
    304378            for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
    … …  
    315389        double sumAI1 = 0.0;            // Sum of A.I_1 products (for matrix, normalisation)
    316390        double sumBI1 = 0.0;            // Sum of B.I_1 products (for matrix, normalisation)
     391
     392        double MxxAI1 = 0.0;
     393        double MyyAI1 = 0.0;
     394        double MxxBI2 = 0.0;
     395        double MyyBI2 = 0.0;
    317396        for (int y = - footprint; y <= footprint; y++) {
    318397            for (int x = - footprint; x <= footprint; x++) {
    … …  
    322401                float i2 = image2->kernel[y][x];
    323402
    324                 double ai2 = a * i2;
     403                double ai2 = a * i2 * normValue;
    325404                double bi2 = b * i2;
    326                 double ai1 = a * i1;
    327                 double bi1 = b * i1;
     405                double ai1 = a * i1 * PS_SQR(normValue);
     406                double bi1 = b * i1 * normValue;
    328407
    329408                if (weight) {
    … …  
    345424                sumAI1 += ai1;
    346425                sumBI1 += bi1;
    347             }
    348         }
     426
     427                if (MOMENTS) {
     428                    MxxAI1 += x*x*ai1;
     429                    MyyAI1 += y*y*ai1;
     430                    MxxBI2 += x*x*bi2;
     431                    MyyBI2 += y*y*bi2;
     432                }
     433            }
     434        }
     435
     436        if (MOMENTS) {
     437            MxxAI1 /= stamp->normSquare1 * PS_SQR(normValue);
     438            MyyAI1 /= stamp->normSquare1 * PS_SQR(normValue);
     439            MxxBI2 /= stamp->normSquare2;
     440            MyyBI2 /= stamp->normSquare2;
     441        }
     442
     443        // fprintf (stderr, "i : %d : M(xx,yy)(AI1,BI2) : %f %f %f %f\n", i, MxxAI1, MyyAI1, MxxBI2, MyyBI2);
     444
     445        // XXX this makes the Chisq portion independent of the normalization and star flux
     446        // but may be mis-scaling between stars of different fluxes
     447        sumAI1 /= PS_SQR(stamp->normI2);
     448        sumBI1 /= PS_SQR(stamp->normI2);
     449        sumAI2 /= PS_SQR(stamp->normI2);
     450        sumBI2 /= PS_SQR(stamp->normI2);
     451
    349452        // Spatial variation
    350453        for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
    … …  
    353456            double ai1 = sumAI1 * poly[iTerm];
    354457            double bi1 = sumBI1 * poly[iTerm];
    355             vector->data.F64[iIndex]             = ai2 - normValue * ai1;
    356             vector->data.F64[iIndex + numParams] = bi2 - normValue * bi1;
     458            vector->data.F64[iIndex]             = ai2 - ai1;
     459            vector->data.F64[iIndex + numParams] = bi2 - bi1;
     460
     461            // fprintf (stderr, "i : %d : V(I1,I2) : %f %f\n", i, vector->data.F64[iIndex], vector->data.F64[iIndex + numParams]);
     462
     463            // add in second moments
     464            if (MOMENTS) {
     465                vector->data.F64[iIndex]             -= kernels->penalty * MxxAI1;
     466                vector->data.F64[iIndex]             -= kernels->penalty * MyyAI1;
     467
     468                vector->data.F64[iIndex + numParams] -= kernels->penalty * MxxBI2;
     469                vector->data.F64[iIndex + numParams] -= kernels->penalty * MyyBI2;
     470            }
    357471        }
    358472    }
    … …  
    414528                // Contribution to chi^2: a_i^2 P_i
    415529                psAssert(isfinite(penalties1->data.F32[i]), "Invalid penalty");
    416                 fprintf (stderr, "penalty: %f + %f (%f * %f)\n", matrix->data.F64[index][index], normSquare1 * penalties1->data.F32[i], normSquare1, penalties1->data.F32[i]);
     530                // fprintf (stderr, "penalty: %f + %f (%f * %f)\n", matrix->data.F64[index][index], normSquare1 * penalties1->data.F32[i], normSquare1, penalties1->data.F32[i]);
    417531                matrix->data.F64[index][index] += normSquare1 * penalties1->data.F32[i];
    418532
    419                 fprintf (stderr, "penalty: (x^%d y^%d fwhm %f) : %f + %f (%f * %f)\n", kernels->u->data.S32[index], kernels->v->data.S32[index], kernels->widths->data.F32[index],
    420                         matrix->data.F64[index + numParams][index + numParams], normSquare2 * penalties2->data.F32[i], normSquare2, penalties2->data.F32[i]);
     533                // fprintf (stderr, "penalty: (x^%d y^%d fwhm %f) : %f + %f (%f * %f)\n", kernels->u->data.S32[index], kernels->v->data.S32[index], kernels->widths->data.F32[index],
     534                // matrix->data.F64[index + numParams][index + numParams], normSquare2 * penalties2->data.F32[i], normSquare2, penalties2->data.F32[i]);
    421535                matrix->data.F64[index + numParams][index + numParams] += normSquare2 * penalties2->data.F32[i];                             
    422536                // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1)
    … …  
    484598//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    485599
     600bool pmSubtractionCalculateMoments(
     601    pmSubtractionKernels *kernels, // Kernels
     602    pmSubtractionStampList *stamps)
     603{
     604    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
     605
     606    // XXX skip this, right?
     607    return true;
     608
     609    // these are only used by DUAL mode
     610    if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) return true;
     611
     612    psTimerStart("pmSubtractionCalculateMoments");
     613
     614    int footprint = stamps->footprint;  // Half-size of stamps
     615
     616    // Loop over each stamp and calculate its normalization factor
     617    for (int i = 0; i < stamps->num; i++) {
     618        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
     619        if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) continue;
     620        if (stamp->status == PM_SUBTRACTION_STAMP_NONE) continue;
     621
     622        pmSubtractionCalculateMomentsStamp(kernels, stamp, footprint, stamps->normWindow2, stamps->normWindow1);
     623    }
     624
     625    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate moments: %f sec", psTimerClear("pmSubtractionCalculateMoments"));
     626
     627    return true;
     628}
     629
     630bool pmSubtractionCalculateMomentsStamp(
     631    pmSubtractionKernels *kernels, // Kernels
     632    pmSubtractionStamp *stamp,          // stamp on which to save normalization)
     633    int footprint,                      // (Half-)Size of stamp
     634    int normWindow1,                    // Window (half-)size for normalisation measurement
     635    int normWindow2                     // Window (half-)size for normalisation measurement
     636    )
     637{
     638    double Mxx, Myy;
     639
     640    int numKernels = kernels->num;
     641
     642    // Generate convolutions: these are generated once and saved
     643    if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) {
     644        psError(psErrorCodeLast(), false, "Unable to convolve stamp");
     645        return false;
     646    }
     647
     648    if (!stamp->MxxI1) {
     649        stamp->MxxI1 = psVectorAlloc (numKernels, PS_TYPE_F32);
     650    }
     651    if (!stamp->MyyI1) {
     652        stamp->MyyI1 = psVectorAlloc (numKernels, PS_TYPE_F32);
     653    }
     654    if (!stamp->MxxI2) {
     655        stamp->MxxI2 = psVectorAlloc (numKernels, PS_TYPE_F32);
     656    }
     657    if (!stamp->MyyI2) {
     658        stamp->MyyI2 = psVectorAlloc (numKernels, PS_TYPE_F32);
     659    }
     660
     661    for (int i = 0; i < numKernels; i++) {
     662        pmSubtractionCalculateMomentsKernel(&Mxx, &Myy, stamp->convolutions1->data[i], footprint, normWindow1);
     663        stamp->MxxI1->data.F32[i] = Mxx / stamp->normI1;
     664        stamp->MyyI1->data.F32[i] = Myy / stamp->normI1;
     665        pmSubtractionCalculateMomentsKernel(&Mxx, &Myy, stamp->convolutions2->data[i], footprint, normWindow2);
     666        stamp->MxxI2->data.F32[i] = Mxx / stamp->normI2;
     667        stamp->MyyI2->data.F32[i] = Myy / stamp->normI2;
     668    }
     669
     670    pmSubtractionCalculateMomentsKernel(&Mxx, &Myy, stamp->image1, footprint, normWindow1);
     671    stamp->MxxI1raw = Mxx / stamp->normI1;
     672    stamp->MyyI1raw = Myy / stamp->normI1;
     673
     674    pmSubtractionCalculateMomentsKernel(&Mxx, &Myy, stamp->image2, footprint, normWindow2);
     675    stamp->MxxI2raw = Mxx / stamp->normI2;
     676    stamp->MyyI2raw = Myy / stamp->normI2;
     677
     678    // fprintf (stderr, "Mxx I1: %f, Myy I1: %f, Mxx I2: %f, Myy I2: %f\n", stamp->MxxI1raw, stamp->MyyI1raw, stamp->MxxI2raw, stamp->MyyI2raw);
     679
     680    return true;
     681}
     682
     683bool pmSubtractionCalculateMomentsKernel(double *Mxx, double *Myy, psKernel *image, int footprint, int window) {
     684
     685    double Sxx = 0.0;
     686    double Syy = 0.0;
     687    for (int y = - footprint; y <= footprint; y++) {
     688        for (int x = - footprint; x <= footprint; x++) {
     689            if (PS_SQR(x) + PS_SQR(y) > PS_SQR(window)) continue;
     690
     691            double flux = image->kernel[y][x];
     692
     693            Sxx += PS_SQR(x) * flux;
     694            Syy += PS_SQR(y) * flux;
     695        }
     696    }
     697    *Mxx = Sxx;
     698    *Myy = Syy;
     699    return true;
     700}
     701
     702///---------
     703
    486704bool pmSubtractionCalculateNormalization(
    487705    pmSubtractionStampList *stamps,
    … …  
    493711
    494712    psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations
     713    psVector *norm2 = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations
    495714
    496715    int footprint = stamps->footprint;  // Half-size of stamps
    … …  
    510729        }
    511730        psVectorAppend(norms, stamp->norm);
     731        psVectorAppend(norm2, stamp->normSquare2 / stamp->normSquare1);
    512732    }
    513733
    … …  
    517737        psFree(stats);
    518738        psFree(norms);
     739        psFree(norm2);
    519740        return false;
    520741    }
    521742    stamps->normValue = stats->robustMedian;
    522743
    523     fprintf (stderr, "norm (1): %f\n", stamps->normValue);
     744    psStatsInit(stats);
     745    if (!psVectorStats(stats, norm2, NULL, NULL, 0)) {
     746        psError(PM_ERR_DATA, false, "Unable to determine median normalisation");
     747        psFree(stats);
     748        psFree(norms);
     749        psFree(norm2);
     750        return false;
     751    }
     752    stamps->normValue2 = stats->robustMedian;
     753
     754    fprintf (stderr, "norm (1): %f (%f)\n", stamps->normValue, stamps->normValue2);
    524755
    525756    psFree(stats);
    526757    psFree(norms);
     758    psFree(norm2);
    527759
    528760    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate normalization: %f sec", psTimerClear("pmSubtractionCalculateNormalization"));
    … …  
    560792    }
    561793    stamp->norm = normI2 / normI1;
     794    stamp->normI1 = normI1;
     795    stamp->normI2 = normI2;
    562796    stamp->normSquare1 = normSquare1;
    563797    stamp->normSquare2 = normSquare2;
    … …  
    675909        break;
    676910      case PM_SUBTRACTION_MODE_DUAL:
    677         status = calculateDualMatrixVector(stamp->matrix, stamp->vector, stamps->normValue, stamp->image1, stamp->image2,
    678                                            weight, window, stamp->convolutions1, stamp->convolutions2,
    679                                            kernels, polyValues, footprint);
     911        status = calculateDualMatrixVector(stamp, stamps->normValue, stamps->normValue2, weight, window, kernels, polyValues, footprint);
    680912        break;
    681913      default:
    … …  
    9241156#endif
    9251157
    926         calculatePenalty(sumMatrix, sumVector, kernels, normSquare1, normSquare2);
     1158        if (PENALTY) {
     1159            calculatePenalty(sumMatrix, sumVector, kernels, normSquare1, normSquare2);
     1160        }
    9271161
    9281162        psVector *solution = NULL;                       // Solution to equation!
    … …  
    9301164        psVectorInit(solution, 0);
    9311165
    932         solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-6);
     1166        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
    9331167
    9341168#if (1)
    … …  
    9471181        }
    9481182
     1183        // for DUAL convolution analysis, we apply the normalization to I1 as follows:
     1184        // I1c = norm I1 + \Sum_i a_i norm I1 \cross k_i
     1185        // I2c =      I2 + \Sum_i b_i      I2 \cross k_i
     1186
     1187        // We absorb the normalization into a_i after the analysis is complete to be consistent
     1188        // with the SINGLE definitions of the convolutions
     1189
    9491190        int numKernels = kernels->num;
    9501191        for (int i = 0; i < numKernels * numSpatial; i++) {
    951             // XXX fprintf (stderr, "keep\n");
    952             kernels->solution1->data.F64[i] = solution->data.F64[i];
     1192            // we solve for coefficients
     1193            kernels->solution1->data.F64[i] = solution->data.F64[i] * stamps->normValue;
    9531194            kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1];
    9541195        }

  • branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.h

    r29165 r29170  
    7878  );
    7979
     80bool pmSubtractionCalculateMoments(
     81    pmSubtractionKernels *kernels, // Kernels
     82    pmSubtractionStampList *stamps);
     83
     84bool pmSubtractionCalculateMomentsStamp(
     85    pmSubtractionKernels *kernels, // Kernels
     86    pmSubtractionStamp *stamp,          // stamp on which to save normalization)
     87    int footprint,                      // (Half-)Size of stamp
     88    int normWindow1,                    // Window (half-)size for normalisation measurement
     89    int normWindow2                     // Window (half-)size for normalisation measurement
     90    );
     91
     92bool pmSubtractionCalculateMomentsKernel(double *Mxx, double *Myy, psKernel *image, int footprint, int window);
     93
    8094#endif

  • branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionMatch.c

    r29165 r29170  
    710710                }
    711711
     712                // step 0b : calculate the moments, pass along to the next steps via stamps->normValue
     713                // XXX this step is now skipped -- delete
     714                psTrace("psModules.imcombine", 3, "Calculating normalization...\n");
     715                if (!pmSubtractionCalculateMoments(kernels, stamps)) {
     716                    psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     717                    goto MATCH_ERROR;
     718                }
     719
    712720                // step 1: generate the elements of the matrix equation Ax = B
    713721                psTrace("psModules.imcombine", 3, "Calculating kernel equations...\n");

  • branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.c

    r29165 r29170  
    344344    stamp->vector = NULL;
    345345    stamp->norm = NAN;
     346    stamp->normI1 = NAN;
     347    stamp->normI2 = NAN;
    346348    stamp->normSquare1 = NAN;
    347349    stamp->normSquare2 = NAN;
     350
     351    stamp->MxxI1 = NULL;
     352    stamp->MyyI1 = NULL;
     353    stamp->MxxI2 = NULL;
     354    stamp->MyyI2 = NULL;
     355
     356    stamp->MxxI1raw = NAN;
     357    stamp->MyyI1raw = NAN;
     358    stamp->MxxI2raw = NAN;
     359    stamp->MyyI2raw = NAN;
    348360
    349361    return stamp;

  • branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.h

    r29165 r29170  
    2626    float normFrac;                     ///< Fraction of flux in window for normalisation window
    2727    float normValue;                    ///< calculated normalization
     28    float normValue2;                   ///< calculated normalization
    2829    psKernel *window1;                  ///< window function generated from ensemble of stamps (input 1)
    2930    psKernel *window2;                  ///< window function generated from ensemble of stamps (input 2)
    … …  
    8687    psVector *vector;                   ///< Least-squares vector, or NULL
    8788    double norm;                        ///< Normalisation difference
     89    double normI1;                       ///< Sum(flux) for image 1
     90    double normI2;                       ///< Sum(flux) for image 2
    8891    double normSquare1;                 ///< Sum(flux^2) for image 1 (used for penalty)
    8992    double normSquare2;                 ///< Sum(flux^2) for image 2 (used for penalty)
    9093    pmSubtractionStampStatus status;    ///< Status of stamp
     94    psVector *MxxI1;                    ///< second moments of convolution images
     95    psVector *MyyI1;                    ///< second moments of convolution images
     96    psVector *MxxI2;                    ///< second moments of convolution images
     97    psVector *MyyI2;                    ///< second moments of convolution images
     98    double MxxI1raw;
     99    double MyyI1raw;
     100    double MxxI2raw;
     101    double MyyI2raw;
    91102} pmSubtractionStamp;
    92103
Note: See TracChangeset for help on using the changeset viewer.