IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26562


Ignore:
Timestamp:
Jan 12, 2010, 12:11:55 PM (16 years ago)
Author:
Paul Price
Message:

Dual convolution working properly. Reworked the math to be more simple: carrying around a single matrix for the stamps now. Masking terms in the solution when doing dual convolution doesn't work as well as using a penalty function, so disabled, soon to be removed completely. Resolved conflicts.

Location:
branches/eam_branches/20091201/psModules/src/imcombine
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r26547 r26562  
    907907                psFree(stamp->weight);
    908908                stamp->image1 = stamp->image2 = stamp->weight = NULL;
    909                 psFree(stamp->matrix1);
    910                 psFree(stamp->matrix2);
    911                 psFree(stamp->matrixX);
    912                 stamp->matrix1 = stamp->matrix2 = stamp->matrixX = NULL;
    913                 psFree(stamp->vector1);
    914                 psFree(stamp->vector2);
    915                 stamp->vector1 = stamp->vector2 = NULL;
     909                psFree(stamp->matrix);
     910                stamp->matrix = NULL;
     911                psFree(stamp->vector);
     912                stamp->vector = NULL;
    916913            } else {
    917914                numGood++;

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

    r26552 r26562  
    1515#include "pmSubtractionVisual.h"
    1616
    17 #define TESTING                         // TESTING output for debugging; may not work with threads!
     17// #define TESTING                         // TESTING output for debugging; may not work with threads!
    1818
    1919// #define USE_WEIGHT                      // Include weight (1/variance) in equation?
    … …  
    2626
    2727// Calculate the least-squares matrix and vector
    28 // XXX why are we calculating these values on each iteration?  aren't they identical (no change in pixel values...)
    2928static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated
    3029                                  psVector *vector, // Least-squares vector, updated
    … …  
    221220}
    222221
     222
    223223// Calculate the least-squares matrix and vector for dual convolution
    224 static bool calculateDualMatrixVector(psImage *matrix1, // Least-squares matrix, updated
    225                                       psVector *vector1, // Least-squares vector, updated
    226                                       psImage *matrix2,  // Least-squares matrix, updated
    227                                       psVector *vector2, // Least-squares vector, updated
    228                                       psImage *matrixX,  // Cross-matrix
     224static bool calculateDualMatrixVector(psImage *matrix, // Least-squares matrix, updated
     225                                      psVector *vector, // Least-squares vector, updated
    229226                                      const psKernel *image1, // Image 1
    230227                                      const psKernel *image2, // Image 2
    … …  
    238235                                      )
    239236{
    240     // A_ij = A_i A_j
    241     // B_ij = B_i B_j
    242     // C_ij = A_i B_j
    243     // d_i = A_i I_2
    244     // e_i = B_i I_2
    245 
    246     // A_i = I_1 * k_i
    247     // B_i = I_2 * k_i
    248 
    249     // Background: A_i = 1.0
    250     // Normalisation: A_i = I_1
    251 
    252237    int numKernels = kernels->num;                      // Number of kernels
    253238    int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
    … …  
    257242    double poly[numPoly];                                 // Polynomial terms
    258243    double poly2[numPoly][numPoly];                       // Polynomial-polynomial values
     244
     245    int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms
     246    int numParams = numKernels * numPoly + 1 + numBackground;       // Number of regular parameters
     247    int numParams2 = numKernels * numPoly;                          // Number of additional parameters for dual
     248    int numDual = numParams + numParams2;                           // Total number of parameters for dual
     249
     250    psAssert(matrix &&
     251             matrix->type.type == PS_TYPE_F64 &&
     252             matrix->numCols == numDual &&
     253             matrix->numRows == numDual,
     254             "Least-squares matrix is bad.");
     255    psAssert(vector &&
     256             vector->type.type == PS_TYPE_F64 &&
     257             vector->n == numDual,
     258             "Least-squares vector is bad.");
    259259
    260260    // Evaluate polynomial-polynomial terms
    … …  
    280280            psKernel *jConv2 = convolutions2->data[j]; // Convolution 2 for index j
    281281
    282             double sumAA = 0.0;         // Sum of convolution products for matrix A
    283             double sumBB = 0.0;         // Sum of convolution products for matrix B
    284             double sumAB = 0.0;         // Sum of convolution products for matrix C
     282            double sumAA = 0.0;         // Sum of convolution products between image 1
     283            double sumBB = 0.0;         // Sum of convolution products between image 2
     284            double sumAB = 0.0;         // Sum of convolution products across images 1 and 2
    285285            for (int y = - footprint; y <= footprint; y++) {
    286286                for (int x = - footprint; x <= footprint; x++) {
    … …  
    312312                    double bb = sumBB * poly2[iTerm][jTerm];
    313313                    double ab = sumAB * poly2[iTerm][jTerm];
    314                     matrix1->data.F64[iIndex][jIndex] = aa;
    315                     matrix1->data.F64[jIndex][iIndex] = aa;
    316                     matrix2->data.F64[iIndex][jIndex] = bb;
    317                     matrix2->data.F64[jIndex][iIndex] = bb;
    318                     matrixX->data.F64[iIndex][jIndex] = ab;
     314
     315                    matrix->data.F64[iIndex][jIndex] = aa;
     316                    matrix->data.F64[jIndex][iIndex] = aa;
     317
     318                    matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb;
     319                    matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb;
     320
     321                    matrix->data.F64[iIndex][jIndex + numParams] = ab;
     322                    matrix->data.F64[jIndex + numParams][iIndex] = ab;
    319323                }
    320324            }
    … …  
    340344                for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
    341345                    double ab = sumAB * poly2[iTerm][jTerm];
    342                     matrixX->data.F64[iIndex][jIndex] = ab;
    343                 }
    344             }
    345         }
    346 
    347         double sumAI2 = 0.0;            // Sum of A.I_2 products (for vector 1)
    348         double sumBI2 = 0.0;            // Sum of B.I_2 products (for vector 2)
    349         double sumAI1 = 0.0;            // Sum of A.I_1 products (for matrix 1, normalisation)
    350         double sumA = 0.0;              // Sum of A (for matrix 1, background)
    351         double sumBI1 = 0.0;            // Sum of B.I_1 products (for matrix X, normalisation)
    352         double sumB = 0.0;              // Sum of B products (for matrix X, background)
    353         double sumI2 = 0.0;             // Sum of I_2 (for vector 1, background)
     346                    matrix->data.F64[iIndex][jIndex + numParams] = ab;
     347                    matrix->data.F64[jIndex + numParams][iIndex] = ab;
     348                }
     349            }
     350        }
     351
     352        double sumAI2 = 0.0;            // Sum of A.I_2 products (for vector)
     353        double sumBI2 = 0.0;            // Sum of B.I_2 products (for vector)
     354        double sumAI1 = 0.0;            // Sum of A.I_1 products (for matrix, normalisation)
     355        double sumA = 0.0;              // Sum of A (for matrix, background)
     356        double sumBI1 = 0.0;            // Sum of B.I_1 products (for matrix, normalisation)
     357        double sumB = 0.0;              // Sum of B products (for matrix, background)
     358        double sumI2 = 0.0;             // Sum of I_2 (for vector, background)
    354359        for (int y = - footprint; y <= footprint; y++) {
    355360            for (int x = - footprint; x <= footprint; x++) {
    … …  
    402407            double b = sumB * poly[iTerm];
    403408
    404             matrix1->data.F64[iIndex][normIndex] = ai1;
    405             matrix1->data.F64[normIndex][iIndex] = ai1;
    406             matrix1->data.F64[iIndex][bgIndex] = a;
    407             matrix1->data.F64[bgIndex][iIndex] = a;
    408             vector1->data.F64[iIndex] = ai2;
    409             vector2->data.F64[iIndex] = bi2;
    410             matrixX->data.F64[iIndex][normIndex] = bi1;
    411             matrixX->data.F64[iIndex][bgIndex] = b;
    412         }
    413     }
    414 
    415     double sumI1 = 0.0;                 // Sum of I_1 (for matrix 1, background-normalisation)
    416     double sumI1I1 = 0.0;               // Sum of I_1^2 (for matrix 1, normalisation-normalisation)
    417     double sum1 = 0.0;                  // Sum of 1 (for matrix 1, background-background)
    418     double sumI2 = 0.0;                 // Sum of I_2 (for vector 1, background)
    419     double sumI1I2 = 0.0;               // Sum of I_1.I_2 (for vector 1, normalisation)
     409            matrix->data.F64[iIndex][normIndex] = ai1;
     410            matrix->data.F64[normIndex][iIndex] = ai1;
     411            matrix->data.F64[iIndex][bgIndex] = a;
     412            matrix->data.F64[bgIndex][iIndex] = a;
     413            matrix->data.F64[iIndex + numParams][normIndex] = bi1;
     414            matrix->data.F64[normIndex][iIndex + numParams] = bi1;
     415            matrix->data.F64[iIndex + numParams][bgIndex] = b;
     416            matrix->data.F64[bgIndex][iIndex + numParams] = b;
     417            vector->data.F64[iIndex] = ai2;
     418            vector->data.F64[iIndex + numParams] = bi2;
     419        }
     420    }
     421
     422    double sumI1 = 0.0;                 // Sum of I_1 (for matrix, background-normalisation)
     423    double sumI1I1 = 0.0;               // Sum of I_1^2 (for matrix, normalisation-normalisation)
     424    double sum1 = 0.0;                  // Sum of 1 (for matrix, background-background)
     425    double sumI2 = 0.0;                 // Sum of I_2 (for vector, background)
     426    double sumI1I2 = 0.0;               // Sum of I_1.I_2 (for vector, normalisation)
    420427    for (int y = - footprint; y <= footprint; y++) {
    421428        for (int x = - footprint; x <= footprint; x++) {
    … …  
    450457        }
    451458    }
    452     matrix1->data.F64[bgIndex][normIndex] = sumI1;
    453     matrix1->data.F64[normIndex][bgIndex] = sumI1;
    454     matrix1->data.F64[normIndex][normIndex] = sumI1I1;
    455     matrix1->data.F64[bgIndex][bgIndex] = sum1;
    456     vector1->data.F64[bgIndex] = sumI2;
    457     vector1->data.F64[normIndex] = sumI1I2;
     459    matrix->data.F64[bgIndex][normIndex] = sumI1;
     460    matrix->data.F64[normIndex][bgIndex] = sumI1;
     461    matrix->data.F64[normIndex][normIndex] = sumI1I1;
     462    matrix->data.F64[bgIndex][bgIndex] = sum1;
     463    vector->data.F64[bgIndex] = sumI2;
     464    vector->data.F64[normIndex] = sumI1I2;
    458465
    459466    return true;
    460467}
    461 
    462 // Merge dual matrices and vectors into single matrix equation
    463 // Have: Aa = Ct.b + d
    464 // Have: Ca = Bb + e
    465 // Set: F = ( A -Ct ;  C -B )
    466 // Set: g = ( a ; b )
    467 // Set: h = ( d ; e )
    468 // So that we combine the above two equations: Fg = h
    469 static bool calculateEquationDual(psImage **outMatrix,
    470                                   psVector **outVector,
    471                                   const psImage *sumMatrix1,
    472                                   const psImage *sumMatrix2,
    473                                   const psImage *sumMatrixX,
    474                                   const psVector *sumVector1,
    475                                   const psVector *sumVector2
    476                                   )
    477 {
    478     psAssert(sumMatrix1 && sumMatrix2 && sumMatrixX, "Require input matrices");
    479     psAssert(sumVector1 && sumVector2, "Require input vectors");
    480     int num1 = sumVector1->n;   // Number of parameters in first set
    481     int num2 = sumVector2->n;   // Number of parameters in second set
    482     int num = num1 + num2;      // Number of parameters in new set
    483 
    484     psAssert(sumMatrix1->type.type == PS_TYPE_F64 &&
    485              sumMatrix2->type.type == PS_TYPE_F64 &&
    486              sumMatrixX->type.type == PS_TYPE_F64 &&
    487              sumVector1->type.type == PS_TYPE_F64 &&
    488              sumVector2->type.type == PS_TYPE_F64,
    489              "Require input type is F64");
    490 
    491     psAssert(outMatrix, "Require output matrix");
    492     psAssert(outVector, "Require output vector");
    493     if (!*outMatrix) {
    494         *outMatrix = psImageAlloc(num, num, PS_TYPE_F64);
    495     }
    496     if (!*outVector) {
    497         *outVector = psVectorAlloc(num, PS_TYPE_F64);
    498     }
    499     psImage *matrix = *outMatrix;
    500     psVector *vector = *outVector;
    501 
    502     psAssert(sumMatrix1->numCols == num1 && sumMatrix1->numRows == num1, "Require size NxN");
    503     psAssert(sumMatrix2->numCols == num2 && sumMatrix2->numRows == num2, "Require size MxM");
    504     psAssert(sumMatrixX->numCols == num1 && sumMatrixX->numRows == num2, "Require size MxN");
    505 
    506     memcpy(vector->data.F64, sumVector1->data.F64, num1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
    507     memcpy(&vector->data.F64[num1], sumVector2->data.F64, num2 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
    508 
    509     for (int i = 0; i < num1; i++) {
    510         memcpy(matrix->data.F64[i], sumMatrix1->data.F64[i], num1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
    511         for (int j = 0, k = num1; j < num2; j++, k++) {
    512             matrix->data.F64[i][k] = sumMatrixX->data.F64[j][i];
    513         }
    514     }
    515     for (int i1 = 0, i2 = num1; i1 < num2; i1++, i2++) {
    516         memcpy(matrix->data.F64[i2], sumMatrixX->data.F64[i1], num1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
    517         for (int j = 0, k = num1; j < num2; j++, k++) {
    518             matrix->data.F64[i2][k] = sumMatrix2->data.F64[i1][j];
    519         }
    520     }
    521 
    522     return true;
    523 }
    524 
    525468
    526469#if 1
    … …  
    562505// Calculate the value of a polynomial, specified by coefficients and polynomial values
    563506double p_pmSubtractionCalculatePolynomial(const psVector *coeff, // Coefficients
    564                                                  const psImage *polyValues, // Polynomial values
    565                                                  int order, // Order of polynomials
    566                                                  int index, // Index at which to begin
    567                                                  int step // Step between subsequent indices
    568                                                  )
     507                                          const psImage *polyValues, // Polynomial values
     508                                          int order, // Order of polynomials
     509                                          int index, // Index at which to begin
     510                                          int step // Step between subsequent indices
     511                                          )
    569512{
    570513    double sum = 0.0;                   // Value of the polynomial sum
    … …  
    581524
    582525double p_pmSubtractionSolutionCoeff(const pmSubtractionKernels *kernels, const psImage *polyValues,
    583                                            int index, bool wantDual)
     526                                    int index, bool wantDual)
    584527{
    585528#if 0
    … …  
    659602    // number of coefficients for the spatial polynomial, normalisation and a constant background offset.
    660603    int numParams = numKernels * numSpatial + 1 + numBackground;
     604    if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
     605        // An additional image is convolved
     606        numParams += numKernels * numSpatial;
     607    }
    661608
    662609    pmSubtractionStamp *stamp = stamps->stamps->data[index]; // Stamp of interest
    … …  
    698645                                                    stamp->xNorm, stamp->yNorm); // Polynomial terms
    699646
    700     bool new = stamp->vector1 ? false : true; // Is this a new run?
     647    bool new = stamp->vector ? false : true; // Is this a new run?
    701648    if (new) {
    702         stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64);
    703         stamp->vector1 = psVectorAlloc(numParams, PS_TYPE_F64);
     649        stamp->matrix = psImageAlloc(numParams, numParams, PS_TYPE_F64);
     650        stamp->vector = psVectorAlloc(numParams, PS_TYPE_F64);
    704651    }
    705652#ifdef TESTING
    706     psImageInit(stamp->matrix1, NAN);
    707     psVectorInit(stamp->vector1, NAN);
     653    psImageInit(stamp->matrix, NAN);
     654    psVectorInit(stamp->vector, NAN);
    708655#endif
    709656
    … …  
    722669    switch (kernels->mode) {
    723670      case PM_SUBTRACTION_MODE_1:
    724         status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image2, stamp->image1,
     671        status = calculateMatrixVector(stamp->matrix, stamp->vector, stamp->image2, stamp->image1,
    725672                                       weight, window, stamp->convolutions1, kernels,
    726673                                       polyValues, footprint, mode);
    727674        break;
    728675      case PM_SUBTRACTION_MODE_2:
    729         status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image1, stamp->image2,
     676        status = calculateMatrixVector(stamp->matrix, stamp->vector, stamp->image1, stamp->image2,
    730677                                       weight, window, stamp->convolutions2, kernels,
    731678                                       polyValues, footprint, mode);
    732679        break;
    733680      case PM_SUBTRACTION_MODE_DUAL:
    734         if (new) {
    735             stamp->matrix2 = psImageAlloc(numKernels * numSpatial, numKernels * numSpatial, PS_TYPE_F64);
    736             stamp->matrixX = psImageAlloc(numParams, numKernels * numSpatial, PS_TYPE_F64);
    737             stamp->vector2 = psVectorAlloc(numKernels * numSpatial, PS_TYPE_F64);
    738         }
    739 #ifdef TESTING
    740         psImageInit(stamp->matrix2, NAN);
    741         psImageInit(stamp->matrixX, NAN);
    742         psVectorInit(stamp->vector2, NAN);
    743 #endif
    744         status = calculateDualMatrixVector(stamp->matrix1, stamp->vector1, stamp->matrix2, stamp->vector2,
    745                                            stamp->matrixX, stamp->image1, stamp->image2, weight, window,
    746                                            stamp->convolutions1, stamp->convolutions2, kernels, polyValues,
    747                                            footprint);
     681        status = calculateDualMatrixVector(stamp->matrix, stamp->vector,stamp->image1, stamp->image2,
     682                                           weight, window, stamp->convolutions1, stamp->convolutions2,
     683                                           kernels, polyValues, footprint);
    748684        break;
    749685      default:
    … …  
    762698    {
    763699        psString matrixName = NULL;
    764         psStringAppend(&matrixName, "matrix1_%d.fits", index);
     700        psStringAppend(&matrixName, "matrix_%d.fits", index);
    765701        psFits *matrixFile = psFitsOpen(matrixName, "w");
    766702        psFree(matrixName);
    767         psFitsWriteImage(matrixFile, NULL, stamp->matrix1, 0, NULL);
     703        psFitsWriteImage(matrixFile, NULL, stamp->matrix, 0, NULL);
    768704        psFitsClose(matrixFile);
    769705
    770706        matrixName = NULL;
    771         psStringAppend(&matrixName, "vector1_%d.fits", index);
    772         psImage *dummy = psImageAlloc(stamp->vector1->n, 1, PS_TYPE_F64);
    773         memcpy(dummy->data.F64[0], stamp->vector1->data.F64,
    774                PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector1->n);
     707        psStringAppend(&matrixName, "vector_%d.fits", index);
     708        psImage *dummy = psImageAlloc(stamp->vector->n, 1, PS_TYPE_F64);
     709        memcpy(dummy->data.F64[0], stamp->vector->data.F64,
     710               PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector->n);
    775711        matrixFile = psFitsOpen(matrixName, "w");
    776712        psFree(matrixName);
    … …  
    778714        psFree(dummy);
    779715        psFitsClose(matrixFile);
    780 
    781         if (stamp->vector2) {
    782             matrixName = NULL;
    783             psStringAppend(&matrixName, "vector2_%d.fits", index);
    784             dummy = psImageAlloc(stamp->vector2->n, 1, PS_TYPE_F64);
    785             memcpy(dummy->data.F64[0], stamp->vector2->data.F64,
    786                    PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector2->n);
    787             matrixFile = psFitsOpen(matrixName, "w");
    788             psFree(matrixName);
    789             psFitsWriteImage(matrixFile, NULL, dummy, 0, NULL);
    790             psFree(dummy);
    791             psFitsClose(matrixFile);
    792         }
    793 
    794         if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
    795             matrixName = NULL;
    796             psStringAppend(&matrixName, "matrix2_%d.fits", index);
    797             matrixFile = psFitsOpen(matrixName, "w");
    798             psFree(matrixName);
    799             psFitsWriteImage(matrixFile, NULL, stamp->matrix2, 0, NULL);
    800             psFitsClose(matrixFile);
    801 
    802             matrixName = NULL;
    803             psStringAppend(&matrixName, "matrixX_%d.fits", index);
    804             matrixFile = psFitsOpen(matrixName, "w");
    805             psFree(matrixName);
    806             psFitsWriteImage(matrixFile, NULL, stamp->matrixX, 0, NULL);
    807             psFitsClose(matrixFile);
    808         }
    809716    }
    810717#endif
    … …  
    898805
    899806    // Check inputs
    900     int numParams = -1;                // Number of parameters
    901     int numParams2 = 0;                // Number of parameters for part solution (DUAL mode)
     807    int numKernels = kernels->num;      // Number of kernel basis functions
     808    int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations
     809    int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms
     810    int numParams = numKernels * numSpatial + 1 + numBackground;    // Number of parameters being solved for
     811    int numSolution1 = numParams, numSolution2 = 0;                 // Number of parameters for each solution
     812    if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
     813        // An additional image is convolved
     814        numSolution2 = numKernels * numSpatial;
     815        numParams += numSolution2;
     816    }
     817
    902818    for (int i = 0; i < stamps->num; i++) {
    903819        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    … …  
    907823        }
    908824
    909         PS_ASSERT_VECTOR_NON_NULL(stamp->vector1, false);
    910         if (numParams == -1) {
    911             numParams = stamp->vector1->n;
    912         }
    913         PS_ASSERT_VECTOR_SIZE(stamp->vector1, (long)numParams, false);
    914         PS_ASSERT_VECTOR_TYPE(stamp->vector1, PS_TYPE_F64, false);
    915         PS_ASSERT_IMAGE_NON_NULL(stamp->matrix1, false);
    916         PS_ASSERT_IMAGE_SIZE(stamp->matrix1, numParams, numParams, false);
    917         PS_ASSERT_IMAGE_TYPE(stamp->matrix1, PS_TYPE_F64, false);
    918 
    919         if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
    920             PS_ASSERT_IMAGE_NON_NULL(stamp->matrix2, false);
    921             PS_ASSERT_IMAGE_NON_NULL(stamp->matrixX, false);
    922             if (numParams2 == 0) {
    923                 numParams2 = stamp->matrix2->numCols;
    924             }
    925             PS_ASSERT_IMAGE_SIZE(stamp->matrix2, numParams2, numParams2, false);
    926             PS_ASSERT_IMAGE_SIZE(stamp->matrixX, numParams, numParams2, false);
    927             PS_ASSERT_IMAGE_TYPE(stamp->matrix2, PS_TYPE_F64, false);
    928             PS_ASSERT_IMAGE_TYPE(stamp->matrixX, PS_TYPE_F64, false);
    929             PS_ASSERT_VECTOR_NON_NULL(stamp->vector2, false);
    930             PS_ASSERT_VECTOR_SIZE(stamp->vector2, (long)numParams2, false);
    931             PS_ASSERT_VECTOR_TYPE(stamp->vector2, PS_TYPE_F64, false);
    932         }
    933     }
    934     if (numParams == -1) {
    935         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "No suitable stamps found.");
    936         return NULL;
     825        PS_ASSERT_VECTOR_NON_NULL(stamp->vector, false);
     826        PS_ASSERT_VECTOR_SIZE(stamp->vector, (long)numParams, false);
     827        PS_ASSERT_VECTOR_TYPE(stamp->vector, PS_TYPE_F64, false);
     828        PS_ASSERT_IMAGE_NON_NULL(stamp->matrix, false);
     829        PS_ASSERT_IMAGE_SIZE(stamp->matrix, numParams, numParams, false);
     830        PS_ASSERT_IMAGE_TYPE(stamp->matrix, PS_TYPE_F64, false);
    937831    }
    938832
    … …  
    958852#ifdef TESTING
    959853              // XXX double-check for NAN in data:
    960                 for (int iy = 0; iy < stamp->matrix1->numRows; iy++) {
    961                     for (int ix = 0; ix < stamp->matrix1->numCols; ix++) {
    962                         if (!isfinite(stamp->matrix1->data.F64[iy][ix])) {
    963                             fprintf (stderr, "WARNING: NAN in matrix1\n");
     854                for (int iy = 0; iy < stamp->matrix->numRows; iy++) {
     855                    for (int ix = 0; ix < stamp->matrix->numCols; ix++) {
     856                        if (!isfinite(stamp->matrix->data.F64[iy][ix])) {
     857                            fprintf (stderr, "WARNING: NAN in matrix\n");
    964858                        }
    965859                    }
    966860                }
    967                 for (int ix = 0; ix < stamp->vector1->n; ix++) {
    968                     if (!isfinite(stamp->vector1->data.F64[ix])) {
    969                         fprintf (stderr, "WARNING: NAN in vector1\n");
     861                for (int ix = 0; ix < stamp->vector->n; ix++) {
     862                    if (!isfinite(stamp->vector->data.F64[ix])) {
     863                        fprintf (stderr, "WARNING: NAN in vector\n");
    970864                    }
    971865                }
    972866#endif
    973867
    974                 (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix1);
    975                 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector1);
     868                (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix);
     869                (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
    976870                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
    977871                numStamps++;
    … …  
    984878        for (int ix = 0; ix < sumVector->n; ix++) {
    985879            if (!isfinite(sumVector->data.F64[ix])) {
    986                 fprintf (stderr, "WARNING: NAN in vector1\n");
     880                fprintf (stderr, "WARNING: NAN in vector\n");
    987881            }
    988882        }
    … …  
    997891        for (int ix = 0; ix < sumVector->n; ix++) {
    998892            if (!isfinite(sumVector->data.F64[ix])) {
    999                 fprintf (stderr, "WARNING: NAN in vector1\n");
     893                fprintf (stderr, "WARNING: NAN in vector\n");
    1000894            }
    1001895        }
    … …  
    1016910#endif
    1017911
    1018 # define SVD_ANALYSIS 1
     912# define SVD_ANALYSIS 0
    1019913# define COEFF_SIG 0.0
    1020914# define SVD_TOL 0.0
    … …  
    12821176        for (int ix = 0; ix < kernels->solution1->n; ix++) {
    12831177            if (!isfinite(kernels->solution1->data.F64[ix])) {
    1284                 fprintf (stderr, "WARNING: NAN in vector1\n");
     1178                fprintf (stderr, "WARNING: NAN in vector\n");
    12851179            }
    12861180        }
    … …  
    12911185
    12921186        // Accumulation of stamp matrices/vectors
    1293         psImage *sumMatrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64);
    1294         psImage *sumMatrix2 = psImageAlloc(numParams2, numParams2, PS_TYPE_F64);
    1295         psImage *sumMatrixX = psImageAlloc(numParams, numParams2, PS_TYPE_F64);
    1296         psVector *sumVector1 = psVectorAlloc(numParams, PS_TYPE_F64);
    1297         psVector *sumVector2 = psVectorAlloc(numParams2, PS_TYPE_F64);
    1298         psImageInit(sumMatrix1, 0.0);
    1299         psImageInit(sumMatrix2, 0.0);
    1300         psImageInit(sumMatrixX, 0.0);
    1301         psVectorInit(sumVector1, 0.0);
    1302         psVectorInit(sumVector2, 0.0);
     1187        psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64);
     1188        psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64);
     1189        psImageInit(sumMatrix, 0.0);
     1190        psVectorInit(sumVector, 0.0);
    13031191
    13041192        int numStamps = 0;              // Number of good stamps
    … …  
    13061194            pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    13071195            if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
    1308                 (void)psBinaryOp(sumMatrix1, sumMatrix1, "+", stamp->matrix1);
    1309                 (void)psBinaryOp(sumMatrix2, sumMatrix2, "+", stamp->matrix2);
    1310                 (void)psBinaryOp(sumMatrixX, sumMatrixX, "+", stamp->matrixX);
    1311                 (void)psBinaryOp(sumVector1, sumVector1, "+", stamp->vector1);
    1312                 (void)psBinaryOp(sumVector2, sumVector2, "+", stamp->vector2);
     1196                (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix);
     1197                (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
     1198
    13131199                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
    13141200                numStamps++;
    … …  
    13181204
    13191205#ifdef TESTING
    1320         psFitsWriteImageSimple ("sumMatrix1.fits", sumMatrix1, NULL);
    1321         psFitsWriteImageSimple ("sumMatrix2.fits", sumMatrix2, NULL);
    1322         psFitsWriteImageSimple ("sumMatrixX.fits", sumMatrixX, NULL);
    1323         psVectorWriteFile("sumVector1.dat", sumVector1);
    1324         psVectorWriteFile("sumVector2.dat", sumVector2);
    1325 #endif
    1326 
     1206        psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, NULL);
     1207        psVectorWriteFile("sumVector.dat", sumVector);
     1208#endif
    13271209
    13281210#if 1
    13291211        int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
    1330         calculatePenalty(sumMatrix1, sumVector1, kernels, sumMatrix1->data.F64[bgIndex][bgIndex]);
    1331         calculatePenalty(sumMatrix2, sumVector2, kernels, sumMatrix1->data.F64[bgIndex][bgIndex]);
    1332 #endif
    1333 
    1334         psImage *sumMatrix = NULL;      // Combined matrix
    1335         psVector *sumVector = NULL;     // Combined vector
    1336         calculateEquationDual(&sumMatrix, &sumVector, sumMatrix1, sumMatrix2,
    1337                               sumMatrixX, sumVector1, sumVector2);
    1338 
    1339 #ifdef TESTING
    1340         psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, NULL);
    1341         {
    1342             psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64);
    1343             psFits *fits = psFitsOpen("sumVector.fits", "w");
    1344             for (int i = 0; i < numParams + numParams2; i++) {
    1345                 image->data.F64[0][i] = sumVector->data.F64[i];
    1346             }
    1347             psFitsWriteImage(fits, NULL, image, 0, NULL);
    1348             psFree(image);
    1349             psFitsClose(fits);
    1350         }
    1351         psVectorWriteFile("sumVector.dat", sumVector);
     1212        calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]);
    13521213#endif
    13531214
    13541215        psVector *solution = NULL;                       // Solution to equation!
    1355         solution = psVectorAlloc(numParams + numParams2, PS_TYPE_F64);
     1216        solution = psVectorAlloc(numParams, PS_TYPE_F64);
    13561217        psVectorInit(solution, 0);
    13571218
    13581219#if 0
    13591220        {
    1360             solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-5);
     1221            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
    13611222            if (!solution) {
    13621223                psError(PS_ERR_UNKNOWN, false, "SVD solution of least-squares equation failed.\n");
    … …  
    13671228
    13681229#ifdef TESTING
    1369             {
    1370                 psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64);
    1371                 psFits *fits = psFitsOpen("solnVector.fits", "w");
    1372                 for (int i = 0; i < numParams + numParams2; i++) {
    1373                     fprintf(stderr, "Solution %d: %lf\n", i, solution->data.F64[i]);
    1374                     image->data.F64[0][i] = solution->data.F64[i];
    1375                 }
    1376                 psFitsWriteImage(fits, NULL, image, 0, NULL);
    1377                 psFree(image);
    1378                 psFitsClose(fits);
     1230            for (int i = 0; i < numParams; i++) {
     1231                fprintf(stderr, "Solution %d: %lf\n", i, solution->data.F64[i]);
    13791232            }
    13801233#endif
    … …  
    13821235            int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations
    13831236            int numKernels = kernels->num; // Number of kernel basis functions
     1237
     1238#define MASK_BASIS(INDEX)                                               \
     1239            {                                                           \
     1240                for (int j = 0, index = INDEX; j < numSpatial; j++, index += numKernels) { \
     1241                    for (int k = 0; k < numParams; k++) {               \
     1242                        sumMatrix->data.F64[index][k] = sumMatrix->data.F64[k][index] = 0.0; \
     1243                    }                                                   \
     1244                    sumMatrix->data.F64[index][index] = 1.0;            \
     1245                    sumVector->data.F64[index] = 0.0;                   \
     1246                }                                                       \
     1247            }
    13841248
    13851249            // Remove a kernel basis for image 1 from the equation
    … …  
    14271291                // have to change this to properly determine the value of the parameter at the centre.
    14281292                double param1 = solution->data.F64[i],
    1429                     param2 = solution->data.F64[numParams + i]; // Parameters of interest
     1293                    param2 = solution->data.F64[numSolution1 + i]; // Parameters of interest
    14301294                bool mask1 = false, mask2 = false;              // Masked the parameter?
    14311295                if (fabs(param1) < thresh) {
    14321296                    psTrace("psModules.imcombine", 7, "Parameter %d: 1 below threshold\n", i);
    1433                     MASK_BASIS_1(i);
     1297                    MASK_BASIS(i);
    14341298                    mask1 = true;
    14351299                }
    14361300                if (fabs(param2) < thresh) {
    14371301                    psTrace("psModules.imcombine", 7, "Parameter %d: 2 below threshold\n", i);
    1438                     MASK_BASIS_2(i);
     1302                    MASK_BASIS(numSolution1 + i);
    14391303                    mask2 = true;
    14401304                }
    … …  
    14431307                    if (fabs(param1) > fabs(param2)) {
    14441308                        psTrace("psModules.imcombine", 7, "Parameter %d: 1 > 2\n", i);
    1445                         MASK_BASIS_2(i);
     1309                        MASK_BASIS(numSolution1 + i);
    14461310                    } else {
    14471311                        psTrace("psModules.imcombine", 7, "Parameter %d: 2 > 1\n", i);
    1448                         MASK_BASIS_1(i);
     1312                        MASK_BASIS(i);
    14491313                    }
    14501314                }
    14511315            }
    14521316
    1453             calculateEquationDual(&sumMatrix, &sumVector, sumMatrix1, sumMatrix2,
    1454                                   sumMatrixX, sumVector1, sumVector2);
    1455 
    14561317#ifdef TESTING
    1457             {
    1458                 psFits *fits = psFitsOpen("sumMatrixFix.fits", "w");
    1459                 psFitsWriteImage(fits, NULL, sumMatrix, 0, NULL);
    1460                 psFitsClose(fits);
    1461             }
    1462             {
    1463                 psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64);
    1464                 psFits *fits = psFitsOpen("sumVectorFix.fits", "w");
    1465                 for (int i = 0; i < numParams + numParams2; i++) {
    1466                     image->data.F64[0][i] = sumVector->data.F64[i];
    1467                 }
    1468                 psFitsWriteImage(fits, NULL, image, 0, NULL);
    1469                 psFree(image);
    1470                 psFitsClose(fits);
    1471             }
    1472 #endif
    1473         }
    1474 #endif
    1475 
    1476         solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-7);
    1477         if (!solution) {
    1478             psError(PS_ERR_UNKNOWN, false, "SVD solution of least-squares equation failed.\n");
    1479             psFree(sumMatrix);
    1480             psFree(sumVector);
    1481             return NULL;
    1482         }
    1483 
    1484         psFree(sumMatrix1);
    1485         psFree(sumMatrix2);
    1486         psFree(sumMatrixX);
    1487         psFree(sumVector1);
    1488         psFree(sumVector2);
     1318            psFitsWriteImageSimple ("sumMatrixFix.fits", sumMatrix, NULL);
     1319            psVectorWriteFile("sumVectorFix.dat", sumVector);
     1320#endif
     1321
     1322        }
     1323#endif
     1324
     1325        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
     1326
     1327
     1328        for (int i = 0; i < solution->n; i++) {
     1329            fprintf(stderr, "Dual solution %d: %lf\n", i, solution->data.F64[i]);
     1330        }
    14891331
    14901332        psFree(sumMatrix);
    14911333        psFree(sumVector);
    14921334
    1493 
    1494 #ifdef TESTING
    1495         {
    1496             psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64);
    1497             psFits *fits = psFitsOpen("solnVectorFix.fits", "w");
    1498             for (int i = 0; i < numParams + numParams2; i++) {
    1499                 fprintf(stderr, "Solution %d: %lf\n", i, solution->data.F64[i]);
    1500                 image->data.F64[0][i] = solution->data.F64[i];
    1501             }
    1502             psFitsWriteImage(fits, NULL, image, 0, NULL);
    1503             psFree(image);
    1504             psFitsClose(fits);
    1505         }
    1506 #endif
    1507 
    15081335        if (!kernels->solution1) {
    1509             kernels->solution1 = psVectorAlloc(numParams, PS_TYPE_F64);
     1336            kernels->solution1 = psVectorAlloc(numSolution1, PS_TYPE_F64);
    15101337        }
    15111338        if (!kernels->solution2) {
    1512             kernels->solution2 = psVectorAlloc(numParams2, PS_TYPE_F64);
    1513         }
    1514 
    1515         memcpy(kernels->solution1->data.F64, solution->data.F64, numParams * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
    1516         memcpy(kernels->solution2->data.F64, &solution->data.F64[numParams],
    1517                numParams2 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
     1339            kernels->solution2 = psVectorAlloc(numSolution2, PS_TYPE_F64);
     1340        }
     1341
     1342        memcpy(kernels->solution1->data.F64, solution->data.F64,
     1343               numSolution1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
     1344        memcpy(kernels->solution2->data.F64, &solution->data.F64[numSolution1],
     1345               numSolution2 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
    15181346
    15191347        psFree(solution);

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

    r26505 r26562  
    5858    psFree(stamp->convolutions1);
    5959    psFree(stamp->convolutions2);
    60 
    61     psFree(stamp->matrix1);
    62     psFree(stamp->matrix2);
    63     psFree(stamp->matrixX);
    64     psFree(stamp->vector1);
    65     psFree(stamp->vector2);
    66 
     60    psFree(stamp->matrix);
     61    psFree(stamp->vector);
    6762}
    6863
    … …  
    284279        }
    285280
    286         outStamp->matrix1 = inStamp->matrix1 ? psImageCopy(NULL, inStamp->matrix1, PS_TYPE_F64) : NULL;
    287         outStamp->matrix2 = inStamp->matrix2 ? psImageCopy(NULL, inStamp->matrix2, PS_TYPE_F64) : NULL;
    288         outStamp->matrixX = inStamp->matrixX ? psImageCopy(NULL, inStamp->matrixX, PS_TYPE_F64) : NULL;
    289         outStamp->vector1 = inStamp->vector1 ? psVectorCopy(NULL, inStamp->vector1, PS_TYPE_F64) : NULL;
    290         outStamp->vector2 = inStamp->vector2 ? psVectorCopy(NULL, inStamp->vector2, PS_TYPE_F64) : NULL;
     281        outStamp->matrix = inStamp->matrix ? psImageCopy(NULL, inStamp->matrix, PS_TYPE_F64) : NULL;
     282        outStamp->vector = inStamp->vector ? psVectorCopy(NULL, inStamp->vector, PS_TYPE_F64) : NULL;
    291283
    292284        out->stamps->data[i] = outStamp;
    … …  
    314306    stamp->convolutions2 = NULL;
    315307
    316     stamp->matrix1 = NULL;
    317     stamp->matrix2 = NULL;
    318     stamp->matrixX = NULL;
    319     stamp->vector1 = NULL;
    320     stamp->vector2 = NULL;
     308    stamp->matrix = NULL;
     309    stamp->vector = NULL;
    321310
    322311    return stamp;
    … …  
    384373
    385374    if (!stamps) {
    386         stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, sysErr, skyErr);
     375        stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, sysErr, skyErr);
    387376    }
    388377
    … …  
    439428                                            subMask, xMin, xMax, yMin, yMax, numCols, numRows, border);
    440429
    441                     // XXX reset for a test:
    442                     xStamp = xList->data.F32[j];
    443                     yStamp = yList->data.F32[j];
    444                     // fprintf (stderr, "find: %d %d ==> %5.1f %5.1f\n", xCentre, yCentre, xStamp, yStamp);
     430                    // XXX reset for a test:
     431                    xStamp = xList->data.F32[j];
     432                    yStamp = yList->data.F32[j];
     433                    // fprintf (stderr, "find: %d %d ==> %5.1f %5.1f\n", xCentre, yCentre, xStamp, yStamp);
    445434                }
    446435            } else {
    … …  
    549538        }
    550539
    551         // fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", xStamp, yStamp, xPix, yPix);
     540        // fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", xStamp, yStamp, xPix, yPix);
    552541
    553542        bool found = false;
    … …  
    635624        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    636625        if (!stamp) continue;
    637         if (!stamp->image1) continue;
    638         if (!stamp->image2) continue;
    639 
    640         float sum1 = 0.0;
    641         float sum2 = 0.0;
    642         for (int y = -size; y <= size; y++) {
    643             for (int x = -size; x <= size; x++) {
    644                 sum1 += stamp->image1->kernel[y][x];
    645                 sum2 += stamp->image2->kernel[y][x];
    646             }
    647         }
    648         norm1->data.F32[i] = sum1;
    649         norm2->data.F32[i] = sum2;
     626        if (!stamp->image1) continue;
     627        if (!stamp->image2) continue;
     628
     629        float sum1 = 0.0;
     630        float sum2 = 0.0;
     631        for (int y = -size; y <= size; y++) {
     632            for (int x = -size; x <= size; x++) {
     633                sum1 += stamp->image1->kernel[y][x];
     634                sum2 += stamp->image2->kernel[y][x];
     635            }
     636        }
     637        norm1->data.F32[i] = sum1;
     638        norm2->data.F32[i] = sum2;
    650639    }
    651640
    … …  
    658647    // generate the window pixels
    659648    for (int y = -size; y <= size; y++) {
    660         for (int x = -size; x <= size; x++) {
    661 
    662             flux->n = 0;
    663             for (int i = 0; i < stamps->num; i++) {
    664                 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    665                 if (!stamp) continue;
    666                 if (!stamp->image1) continue;
    667                 if (!stamp->image2) continue;
    668                
    669                 psVectorAppend (flux, stamp->image1->kernel[y][x] / norm1->data.F32[i]);
    670                 psVectorAppend (flux, stamp->image2->kernel[y][x] / norm2->data.F32[i]);
    671             }
    672 
    673             psStatsInit (stats);
    674             if (!psVectorStats (stats, flux, NULL, NULL, 0)) {
    675                 psAbort ("failed to generate stats");
    676             }
    677             stamps->window->kernel[y][x] = stats->robustMedian;
    678             if (maxValue < stats->robustMedian) {
    679                 maxValue = stats->robustMedian;
    680             }
    681         }
     649        for (int x = -size; x <= size; x++) {
     650
     651            flux->n = 0;
     652            for (int i = 0; i < stamps->num; i++) {
     653                pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
     654                if (!stamp) continue;
     655                if (!stamp->image1) continue;
     656                if (!stamp->image2) continue;
     657
     658                psVectorAppend (flux, stamp->image1->kernel[y][x] / norm1->data.F32[i]);
     659                psVectorAppend (flux, stamp->image2->kernel[y][x] / norm2->data.F32[i]);
     660            }
     661
     662            psStatsInit (stats);
     663            if (!psVectorStats (stats, flux, NULL, NULL, 0)) {
     664                psAbort ("failed to generate stats");
     665            }
     666            stamps->window->kernel[y][x] = stats->robustMedian;
     667            if (maxValue < stats->robustMedian) {
     668                maxValue = stats->robustMedian;
     669            }
     670        }
    682671    }
    683672
    684673    // re-normalize so chisquare values are sensible
    685674    for (int y = -size; y <= size; y++) {
    686         for (int x = -size; x <= size; x++) {
    687             stamps->window->kernel[y][x] /= maxValue;
    688         }
     675        for (int x = -size; x <= size; x++) {
     676            stamps->window->kernel[y][x] /= maxValue;
     677        }
    689678    }
    690679
    691680# ifdef TESTING
    692681    {
    693         psFits *fits = psFitsOpen ("window.fits", "w");
    694         psFitsWriteImage (fits, NULL, stamps->window->image, 0, NULL);
    695         psFitsClose (fits);
     682        psFits *fits = psFitsOpen ("window.fits", "w");
     683        psFitsWriteImage (fits, NULL, stamps->window->image, 0, NULL);
     684        psFitsClose (fits);
    696685    }
    697686# endif
    … …  
    743732            return false;
    744733        }
    745         // fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", stamp->x, stamp->y, x, y);
     734        // fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", stamp->x, stamp->y, x, y);
    746735
    747736        // Catch memory leaks --- these should have been freed and NULLed before
    … …  
    768757
    769758            if ((isfinite(stamps->skyErr) && (stamps->skyErr > 0)) ||
    770                 (isfinite(stamps->sysErr) && (stamps->sysErr > 0))) {
     759                (isfinite(stamps->sysErr) && (stamps->sysErr > 0))) {
    771760                float sysErr = 0.25 * PS_SQR(stamps->sysErr); // Systematic error
    772                 float skyErr = stamps->skyErr;
     761                float skyErr = stamps->skyErr;
    773762                psKernel *image1 = stamp->image1, *image2 = stamp->image2; // Input images
    774763                for (int y = -size; y <= size; y++) {
    … …  
    823812            y->data.F32[numOut] = source->peak->yf;
    824813        }
    825         // fprintf (stderr, "stamp: %5.1f %5.1f\n", x->data.F32[numOut], y->data.F32[numOut]);
     814        // fprintf (stderr, "stamp: %5.1f %5.1f\n", x->data.F32[numOut], y->data.F32[numOut]);
    826815        numOut++;
    827816    }
    … …  
    873862
    874863    for (int i = 0; i < stamps->num; i++) {
    875         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    876         if (!stamp) continue;
    877         if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
    878         stamp->status = PM_SUBTRACTION_STAMP_CALCULATE;
     864        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
     865        if (!stamp) continue;
     866        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
     867        stamp->status = PM_SUBTRACTION_STAMP_CALCULATE;
    879868    }
    880869    return true;

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.h

    r26505 r26562  
    2424    psArray *flux;                      ///< Fluxes for possible stamps (or NULL)
    2525    int footprint;                      ///< Half-size of stamps
    26     psKernel *window;                   ///< window function generated from ensemble of stamps
     26    psKernel *window;                   ///< window function generated from ensemble of stamps
    2727    float sysErr;                       ///< Systematic error
    2828    float skyErr;                       ///< Systematic error
    … …  
    7676    psArray *convolutions1;             ///< Convolutions of image 1 for each kernel component, or NULL
    7777    psArray *convolutions2;             ///< Convolutions of image 2 for each kernel component, or NULL
    78     psImage *matrix1, *matrix2;         ///< Least-squares matrices for each image, or NULL
    79     psImage *matrixX;                   ///< Cross-matrix (for mode DUAL), or NULL
    80     psVector *vector1, *vector2;        ///< Least-squares vectors for each image, or NULL
     78    psImage *matrix;                    ///< Least-squares matrix, or NULL
     79    psVector *vector;                   ///< Least-squares vector, or NULL
    8180    pmSubtractionStampStatus status;    ///< Status of stamp
    8281} pmSubtractionStamp;

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

    r26547 r26562  
    159159        if (stamp == NULL) continue;
    160160
    161         psImage *im = stamp->matrix1;
    162         if (im == NULL) im = stamp->matrix2;
    163         if (im == NULL) im = stamp->matrixX;
     161        psImage *im = stamp->matrix;
    164162        if (im == NULL) continue;
    165163
    … …  
    189187        if (stamp == NULL) continue;
    190188
    191         psImage *im = stamp->matrix1;
    192         if (im == NULL) im = stamp->matrix2;
    193         if (im == NULL) im = stamp->matrixX;
     189        psImage *im = stamp->matrix;
    194190        if (im == NULL) continue;
    195191
Note: See TracChangeset for help on using the changeset viewer.