IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 8, 2010, 2:24:01 PM (17 years ago)
Author:
Paul Price
Message:

Dual convolution seems to be working again!

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r26505 r26547  
    1515#include "pmSubtractionVisual.h"
    1616
    17 // #define TESTING                         // TESTING output for debugging; may not work with threads!
    18 
    19 #define USE_WEIGHT                      // Include weight (1/variance) in equation?
     17#define TESTING                         // TESTING output for debugging; may not work with threads!
     18
     19// #define USE_WEIGHT                      // Include weight (1/variance) in equation?
    2020// #define USE_WINDOW                      // Include weight (1/variance) in equation?
    2121
    … …  
    3737                                  const psImage *polyValues, // Spatial polynomial values
    3838                                  int footprint, // (Half-)Size of stamp
    39                                   const pmSubtractionEquationCalculationMode mode
     39                                  const pmSubtractionEquationCalculationMode mode
    4040                                  )
    4141{
    … …  
    6969    }
    7070
    71     // initialize the matrix and vector for NOP on all coeffs.  we only fill in the coeffs we 
     71    // initialize the matrix and vector for NOP on all coeffs.  we only fill in the coeffs we
    7272    // choose to calculate
    7373    psImageInit(matrix, 0.0);
    7474    psVectorInit(vector, 1.0);
    7575    for (int i = 0; i < matrix->numCols; i++) {
    76         matrix->data.F64[i][i] = 1.0;
     76        matrix->data.F64[i][i] = 1.0;
    7777    }
    7878
    … …  
    8585
    8686    for (int i = 0; i < numKernels; i++) {
    87         psKernel *iConv = convolutions->data[i]; // Convolution for index i
    88         for (int j = i; j < numKernels; j++) {
    89             psKernel *jConv = convolutions->data[j]; // Convolution for index j
    90 
    91             double sumCC = 0.0;         // Sum of convolution products
    92             for (int y = - footprint; y <= footprint; y++) {
    93                 for (int x = - footprint; x <= footprint; x++) {
    94                     double cc = iConv->kernel[y][x] * jConv->kernel[y][x];
    95                     if (weight) {
    96                         cc *= weight->kernel[y][x];
    97                     }
    98                     if (window) {
    99                         cc *= window->kernel[y][x];
    100                     }
    101                     sumCC += cc;
    102                 }
    103             }
    104 
    105             // Spatial variation of kernel coeffs
    106             if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
    107                 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
    108                     for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
    109                         double value = sumCC * poly2[iTerm][jTerm];
    110                         matrix->data.F64[iIndex][jIndex] = value;
    111                         matrix->data.F64[jIndex][iIndex] = value;
    112                     }
    113                 }
    114             }
    115         }
    116 
    117         double sumRC = 0.0;             // Sum of the reference-convolution products
    118         double sumIC = 0.0;             // Sum of the input-convolution products
    119         double sumC = 0.0;              // Sum of the convolution
    120         for (int y = - footprint; y <= footprint; y++) {
    121             for (int x = - footprint; x <= footprint; x++) {
    122                 float conv = iConv->kernel[y][x];
    123                 float in = input->kernel[y][x];
    124                 float ref = reference->kernel[y][x];
    125                 double ic = in * conv;
    126                 double rc = ref * conv;
    127                 double c = conv;
    128                 if (weight) {
    129                     float wtVal = weight->kernel[y][x];
    130                     ic *= wtVal;
    131                     rc *= wtVal;
    132                     c *= wtVal;
    133                 }
    134                 if (window) {
    135                     float winVal = window->kernel[y][x];
    136                     ic *= winVal;
    137                     rc *= winVal;
    138                     c  *= winVal;
    139                 }
    140                 sumIC += ic;
    141                 sumRC += rc;
    142                 sumC += c;
    143             }
    144         }
    145         // Spatial variation
    146         for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
    147             double normTerm = sumRC * poly[iTerm];
    148             double bgTerm = sumC * poly[iTerm];
    149             if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
    150                 matrix->data.F64[iIndex][normIndex] = normTerm;
    151                 matrix->data.F64[normIndex][iIndex] = normTerm;
    152             }
    153             if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
    154                 matrix->data.F64[iIndex][bgIndex] = bgTerm;
    155                 matrix->data.F64[bgIndex][iIndex] = bgTerm;
    156             }
    157             if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
    158                 vector->data.F64[iIndex] = sumIC * poly[iTerm];
    159                 // XXX TEST for Hermitians: do not calculate A - norm*B - \sum(k x B),
    160                 // instead, calculate A - \sum(k x B), with full hermitians
    161                 if (0 && !(mode & PM_SUBTRACTION_EQUATION_NORM)) {
    162                     // subtract norm * sumRC * poly[iTerm]
    163                     psAssert (kernels->solution1, "programming error: define solution first!");
    164                     int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
    165                     double norm = fabs(kernels->solution1->data.F64[normIndex]);  // Normalisation
    166                     vector->data.F64[iIndex] -= norm * normTerm;
    167                 }
    168             }
    169         }
     87        psKernel *iConv = convolutions->data[i]; // Convolution for index i
     88        for (int j = i; j < numKernels; j++) {
     89            psKernel *jConv = convolutions->data[j]; // Convolution for index j
     90
     91            double sumCC = 0.0;         // Sum of convolution products
     92            for (int y = - footprint; y <= footprint; y++) {
     93                for (int x = - footprint; x <= footprint; x++) {
     94                    double cc = iConv->kernel[y][x] * jConv->kernel[y][x];
     95                    if (weight) {
     96                        cc *= weight->kernel[y][x];
     97                    }
     98                    if (window) {
     99                        cc *= window->kernel[y][x];
     100                    }
     101                    sumCC += cc;
     102                }
     103            }
     104
     105            // Spatial variation of kernel coeffs
     106            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
     107                for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
     108                    for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
     109                        double value = sumCC * poly2[iTerm][jTerm];
     110                        matrix->data.F64[iIndex][jIndex] = value;
     111                        matrix->data.F64[jIndex][iIndex] = value;
     112                    }
     113                }
     114            }
     115        }
     116
     117        double sumRC = 0.0;             // Sum of the reference-convolution products
     118        double sumIC = 0.0;             // Sum of the input-convolution products
     119        double sumC = 0.0;              // Sum of the convolution
     120        for (int y = - footprint; y <= footprint; y++) {
     121            for (int x = - footprint; x <= footprint; x++) {
     122                float conv = iConv->kernel[y][x];
     123                float in = input->kernel[y][x];
     124                float ref = reference->kernel[y][x];
     125                double ic = in * conv;
     126                double rc = ref * conv;
     127                double c = conv;
     128                if (weight) {
     129                    float wtVal = weight->kernel[y][x];
     130                    ic *= wtVal;
     131                    rc *= wtVal;
     132                    c *= wtVal;
     133                }
     134                if (window) {
     135                    float winVal = window->kernel[y][x];
     136                    ic *= winVal;
     137                    rc *= winVal;
     138                    c  *= winVal;
     139                }
     140                sumIC += ic;
     141                sumRC += rc;
     142                sumC += c;
     143            }
     144        }
     145        // Spatial variation
     146        for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
     147            double normTerm = sumRC * poly[iTerm];
     148            double bgTerm = sumC * poly[iTerm];
     149            if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
     150                matrix->data.F64[iIndex][normIndex] = normTerm;
     151                matrix->data.F64[normIndex][iIndex] = normTerm;
     152            }
     153            if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
     154                matrix->data.F64[iIndex][bgIndex] = bgTerm;
     155                matrix->data.F64[bgIndex][iIndex] = bgTerm;
     156            }
     157            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
     158                vector->data.F64[iIndex] = sumIC * poly[iTerm];
     159                // XXX TEST for Hermitians: do not calculate A - norm*B - \sum(k x B),
     160                // instead, calculate A - \sum(k x B), with full hermitians
     161                if (0 && !(mode & PM_SUBTRACTION_EQUATION_NORM)) {
     162                    // subtract norm * sumRC * poly[iTerm]
     163                    psAssert (kernels->solution1, "programming error: define solution first!");
     164                    int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     165                    double norm = fabs(kernels->solution1->data.F64[normIndex]);  // Normalisation
     166                    vector->data.F64[iIndex] -= norm * normTerm;
     167                }
     168            }
     169        }
    170170    }
    171171
    … …  
    182182            double rr = PS_SQR(ref);
    183183            double one = 1.0;
    184             if (weight) {
    185                 float wtVal = weight->kernel[y][x];
    186                 rr *= wtVal;
    187                 ir *= wtVal;
    188                 in *= wtVal;
    189                 ref *= wtVal;
    190                 one *= wtVal;
    191             }
    192             if (window) {
    193                 float  winVal = window->kernel[y][x];
    194                 rr      *= winVal;
    195                 ir      *= winVal;
    196                 in      *= winVal;
    197                 ref *= winVal;
    198                 one *= winVal;
    199             }
     184            if (weight) {
     185                float wtVal = weight->kernel[y][x];
     186                rr *= wtVal;
     187                ir *= wtVal;
     188                in *= wtVal;
     189                ref *= wtVal;
     190                one *= wtVal;
     191            }
     192            if (window) {
     193                float  winVal = window->kernel[y][x];
     194                rr      *= winVal;
     195                ir      *= winVal;
     196                in      *= winVal;
     197                ref *= winVal;
     198                one *= winVal;
     199            }
    200200            sumRR += rr;
    201201            sumIR += ir;
    … …  
    206206    }
    207207    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
    208         matrix->data.F64[normIndex][normIndex] = sumRR;
    209         vector->data.F64[normIndex] = sumIR;
    210         // subtract sum over kernels * kernel solution
     208        matrix->data.F64[normIndex][normIndex] = sumRR;
     209        vector->data.F64[normIndex] = sumIR;
     210        // subtract sum over kernels * kernel solution
    211211    }
    212212    if (mode & PM_SUBTRACTION_EQUATION_BG) {
    213         matrix->data.F64[bgIndex][bgIndex] = sum1;
    214         vector->data.F64[bgIndex] = sumI;
     213        matrix->data.F64[bgIndex][bgIndex] = sum1;
     214        vector->data.F64[bgIndex] = sumI;
    215215    }
    216216    if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) {
    217         matrix->data.F64[normIndex][bgIndex] = sumR;
    218         matrix->data.F64[bgIndex][normIndex] = sumR;
     217        matrix->data.F64[normIndex][bgIndex] = sumR;
     218        matrix->data.F64[bgIndex][normIndex] = sumR;
    219219    }
    220220    return true;
    … …  
    288288                    double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x];
    289289                    double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x];
    290                     if (weight) {
    291                         float wtVal = weight->kernel[y][x];
    292                         aa *= wtVal;
    293                         bb *= wtVal;
    294                         ab *= wtVal;
    295                     }
    296                     if (window) {
    297                         float wtVal = window->kernel[y][x];
    298                         aa *= wtVal;
    299                         bb *= wtVal;
    300                         ab *= wtVal;
    301                     }
     290                    if (weight) {
     291                        float wtVal = weight->kernel[y][x];
     292                        aa *= wtVal;
     293                        bb *= wtVal;
     294                        ab *= wtVal;
     295                    }
     296                    if (window) {
     297                        float wtVal = window->kernel[y][x];
     298                        aa *= wtVal;
     299                        bb *= wtVal;
     300                        ab *= wtVal;
     301                    }
    302302                    sumAA += aa;
    303303                    sumBB += bb;
    … …  
    326326                for (int x = - footprint; x <= footprint; x++) {
    327327                    double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x];
    328                     if (weight) {
    329                         ab *= weight->kernel[y][x];
    330                     }
    331                     if (window) {
    332                         ab *= window->kernel[y][x];
    333                     }
     328                    if (weight) {
     329                        ab *= weight->kernel[y][x];
     330                    }
     331                    if (window) {
     332                        ab *= window->kernel[y][x];
     333                    }
    334334                    sumAB += ab;
    335335                }
    … …  
    352352        double sumB = 0.0;              // Sum of B products (for matrix X, background)
    353353        double sumI2 = 0.0;             // Sum of I_2 (for vector 1, background)
    354         double sumI1I2 = 0.0;           // Sum of I_1.I_2 (for vector 1, normalisation)
    355354        for (int y = - footprint; y <= footprint; y++) {
    356355            for (int x = - footprint; x <= footprint; x++) {
    357                 float a = iConv1->kernel[y][x];
    358                 float b = iConv2->kernel[y][x];
     356                double a = iConv1->kernel[y][x];
     357                double b = iConv2->kernel[y][x];
    359358                float i1 = image1->kernel[y][x];
    360359                float i2 = image2->kernel[y][x];
    … …  
    364363                double ai1 = a * i1;
    365364                double bi1 = b * i1;
    366                 double i1i2 = i1 * i2;
    367 
    368                 if (weight) {
    369                     float wtVal = weight->kernel[y][x];
    370                     ai2 *= wtVal;
    371                     bi2 *= wtVal;
    372                     ai1 *= wtVal;
    373                     bi1 *= wtVal;
    374                     i1i2 *= wtVal;
    375                     a *= wtVal;
    376                     b *= wtVal;
    377                     i2 *= wtVal;
    378                 }
    379                 if (window) {
    380                     float wtVal = window->kernel[y][x];
    381                     ai2 *= wtVal;
    382                     bi2 *= wtVal;
    383                     ai1 *= wtVal;
    384                     bi1 *= wtVal;
    385                     i1i2 *= wtVal;
    386                     a *= wtVal;
    387                     b *= wtVal;
    388                     i2 *= wtVal;
    389                 }
     365
     366                if (weight) {
     367                    float wtVal = weight->kernel[y][x];
     368                    ai2 *= wtVal;
     369                    bi2 *= wtVal;
     370                    ai1 *= wtVal;
     371                    bi1 *= wtVal;
     372                    a *= wtVal;
     373                    b *= wtVal;
     374                    i2 *= wtVal;
     375                }
     376                if (window) {
     377                    float wtVal = window->kernel[y][x];
     378                    ai2 *= wtVal;
     379                    bi2 *= wtVal;
     380                    ai1 *= wtVal;
     381                    bi1 *= wtVal;
     382                    a *= wtVal;
     383                    b *= wtVal;
     384                    i2 *= wtVal;
     385                }
    390386                sumAI2 += ai2;
    391387                sumBI2 += bi2;
    … …  
    395391                sumB += b;
    396392                sumI2 += i2;
    397                 sumI1I2 += i1i2;
    398393            }
    399394        }
    … …  
    425420    for (int y = - footprint; y <= footprint; y++) {
    426421        for (int x = - footprint; x <= footprint; x++) {
    427             float i1 = image1->kernel[y][x];
    428             float i2 = image2->kernel[y][x];
     422            double i1 = image1->kernel[y][x];
     423            double i2 = image2->kernel[y][x];
    429424
    430425            double i1i1 = i1 * i1;
    … …  
    432427            double i1i2 = i1 * i2;
    433428
    434             if (weight) {
    435                 float wtVal = weight->kernel[y][x];
    436                 i1 *= wtVal;
    437                 i1i1 *= wtVal;
    438                 one *= wtVal;
    439                 i2 *= wtVal;
    440                 i1i2 *= wtVal;
    441             }
    442             if (window) {
    443                 float wtVal = window->kernel[y][x];
    444                 i1 *= wtVal;
    445                 i1i1 *= wtVal;
    446                 one *= wtVal;
    447                 i2 *= wtVal;
    448                 i1i2 *= wtVal;
    449             }
     429            if (weight) {
     430                float wtVal = weight->kernel[y][x];
     431                i1 *= wtVal;
     432                i1i1 *= wtVal;
     433                one *= wtVal;
     434                i2 *= wtVal;
     435                i1i2 *= wtVal;
     436            }
     437            if (window) {
     438                float wtVal = window->kernel[y][x];
     439                i1 *= wtVal;
     440                i1i1 *= wtVal;
     441                one *= wtVal;
     442                i2 *= wtVal;
     443                i1i2 *= wtVal;
     444            }
    450445            sumI1 += i1;
    451446            sumI1I1 += i1i1;
    … …  
    515510        memcpy(matrix->data.F64[i], sumMatrix1->data.F64[i], num1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
    516511        for (int j = 0, k = num1; j < num2; j++, k++) {
    517             matrix->data.F64[i][k] = - sumMatrixX->data.F64[j][i];
     512            matrix->data.F64[i][k] = sumMatrixX->data.F64[j][i];
    518513        }
    519514    }
    … …  
    521516        memcpy(matrix->data.F64[i2], sumMatrixX->data.F64[i1], num1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
    522517        for (int j = 0, k = num1; j < num2; j++, k++) {
    523             matrix->data.F64[i2][k] = - sumMatrix2->data.F64[i1][j];
     518            matrix->data.F64[i2][k] = sumMatrix2->data.F64[i1][j];
    524519        }
    525520    }
    … …  
    529524
    530525
     526#if 1
    531527// Add in penalty term to least-squares vector
    532528static bool calculatePenalty(psImage *matrix,                     // Matrix to which to add in penalty term
    … …  
    547543            for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
    548544                // Contribution to chi^2: a_i^2 P_i
    549                 if (!isfinite(penalties->data.F32[i])) {
    550                     psAbort ("invalid penalty");
    551                 }
     545                if (!isfinite(penalties->data.F32[i])) {
     546                    psAbort ("invalid penalty");
     547                }
    552548                matrix->data.F64[index][index] -= norm * penalties->data.F32[i];
    553549            }
    … …  
    557553    return true;
    558554}
     555#endif
     556
    559557
    560558//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    … …  
    657655    int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms
    658656
    659     // numKernels is the number of unique kernel images (one for each Gaussian modified by a specific polynomial). 
     657    // numKernels is the number of unique kernel images (one for each Gaussian modified by a specific polynomial).
    660658    // = \sum_i^N_Gaussians [(order + 1) * (order + 2) / 2], eg for 1 Gauss and 1st order, numKernels = 3
    661659
    … …  
    728726        status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image2, stamp->image1,
    729727                                       weight, window, stamp->convolutions1, kernels,
    730                                        polyValues, footprint, mode);
     728                                       polyValues, footprint, mode);
    731729        break;
    732730      case PM_SUBTRACTION_MODE_2:
    733731        status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image1, stamp->image2,
    734732                                       weight, window, stamp->convolutions2, kernels,
    735                                        polyValues, footprint, mode);
     733                                       polyValues, footprint, mode);
    736734        break;
    737735      case PM_SUBTRACTION_MODE_DUAL:
    … …  
    834832        }
    835833
    836         if ((stamp->x <= 0.0) && (stamp->y <= 0.0)) {
    837             psAbort ("bad stamp");
    838         }
    839         if (!isfinite(stamp->x) && !isfinite(stamp->y)) {
    840             psAbort ("bad stamp");
    841         }
     834        if ((stamp->x <= 0.0) && (stamp->y <= 0.0)) {
     835            psAbort ("bad stamp");
     836        }
     837        if (!isfinite(stamp->x) && !isfinite(stamp->y)) {
     838            psAbort ("bad stamp");
     839        }
    842840
    843841        if (pmSubtractionThreaded()) {
    … …  
    894892double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps);
    895893
    896 bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, 
    897                                 const pmSubtractionStampList *stamps,
    898                                 const pmSubtractionEquationCalculationMode mode)
     894bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels,
     895                                const pmSubtractionStampList *stamps,
     896                                const pmSubtractionEquationCalculationMode mode)
    899897{
    900898    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
    … …  
    10241022# define SVD_TOL 0.0
    10251023
    1026         psVector *solution = NULL;
    1027         psVector *solutionErr = NULL;
     1024        psVector *solution = NULL;
     1025        psVector *solutionErr = NULL;
    10281026
    10291027#ifdef TESTING
    1030         psFitsWriteImageSimple("A.fits", sumMatrix, NULL);
    1031         psVectorWriteFile ("B.dat", sumVector);
    1032 #endif
    1033 
    1034         // Use SVD to determine the kernel coeffs (and validate)
    1035         if (SVD_ANALYSIS) {
    1036 
    1037             // We have sumVector and sumMatrix.  we are trying to solve the following equation:
    1038             // sumMatrix * x = sumVector.
    1039 
    1040             // we can use any standard matrix inversion to solve this.  However, the basis
    1041             // functions in general have substantial correlation, so that the solution may be
    1042             // somewhat poorly determined or unstable.  If not numerically ill-conditioned, the
    1043             // system of equations may be statistically ill-conditioned.  Noise in the image
    1044             // will drive insignificant, but correlated, terms in the solution.  To avoid these
    1045             // problems, we can use SVD to identify numerically unconstrained values and to
    1046             // avoid statistically badly determined value.
    1047 
    1048             // A = sumMatrix, B = sumVector
    1049             // SVD: A = U w V^T  -> A^{-1} = V (1/w) U^T
    1050             // x = V (1/w) (U^T B)
    1051             // \sigma_x = sqrt(diag(A^{-1}))
    1052             // solve for x and A^{-1} to get x & dx
    1053             // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0
    1054             // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0
    1055 
    1056             // If I use the SVD trick to re-condition the matrix, I need to break out the
    1057             // kernel and normalization terms from the background term.
    1058             // XXX is this true?  or was this due to an error in the analysis?
    1059 
    1060             int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
    1061 
    1062             // now pull out the kernel elements into their own square matrix
    1063             psImage  *kernelMatrix = psImageAlloc  (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64);
    1064             psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64);
    1065 
    1066             for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) {
    1067                 if (ix == bgIndex) continue;
    1068                 for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) {
    1069                     if (iy == bgIndex) continue;
    1070                     kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix];
    1071                     ky++;
    1072                 }
    1073                 kernelVector->data.F64[kx] = sumVector->data.F64[ix];
    1074                 kx++;
    1075             }
    1076 
    1077             psImage *U = NULL;
    1078             psImage *V = NULL;
    1079             psVector *w = NULL;
    1080             if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) {
    1081                 psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n");
    1082                 return NULL;
    1083             }
    1084 
    1085             // calculate A_inverse:
    1086             // Ainv = V * w * U^T
    1087             psImage *wUt  = p_pmSubSolve_wUt (w, U);
    1088             psImage *Ainv = p_pmSubSolve_VwUt (V, wUt);
    1089             psImage *Xvar = NULL;
    1090             psFree (wUt);
     1028        psFitsWriteImageSimple("A.fits", sumMatrix, NULL);
     1029        psVectorWriteFile("B.dat", sumVector);
     1030#endif
     1031
     1032        // Use SVD to determine the kernel coeffs (and validate)
     1033        if (SVD_ANALYSIS) {
     1034
     1035            // We have sumVector and sumMatrix.  we are trying to solve the following equation:
     1036            // sumMatrix * x = sumVector.
     1037
     1038            // we can use any standard matrix inversion to solve this.  However, the basis
     1039            // functions in general have substantial correlation, so that the solution may be
     1040            // somewhat poorly determined or unstable.  If not numerically ill-conditioned, the
     1041            // system of equations may be statistically ill-conditioned.  Noise in the image
     1042            // will drive insignificant, but correlated, terms in the solution.  To avoid these
     1043            // problems, we can use SVD to identify numerically unconstrained values and to
     1044            // avoid statistically badly determined value.
     1045
     1046            // A = sumMatrix, B = sumVector
     1047            // SVD: A = U w V^T  -> A^{-1} = V (1/w) U^T
     1048            // x = V (1/w) (U^T B)
     1049            // \sigma_x = sqrt(diag(A^{-1}))
     1050            // solve for x and A^{-1} to get x & dx
     1051            // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0
     1052            // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0
     1053
     1054            // If I use the SVD trick to re-condition the matrix, I need to break out the
     1055            // kernel and normalization terms from the background term.
     1056            // XXX is this true?  or was this due to an error in the analysis?
     1057
     1058            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
     1059
     1060            // now pull out the kernel elements into their own square matrix
     1061            psImage  *kernelMatrix = psImageAlloc  (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64);
     1062            psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64);
     1063
     1064            for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) {
     1065                if (ix == bgIndex) continue;
     1066                for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) {
     1067                    if (iy == bgIndex) continue;
     1068                    kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix];
     1069                    ky++;
     1070                }
     1071                kernelVector->data.F64[kx] = sumVector->data.F64[ix];
     1072                kx++;
     1073            }
     1074
     1075            psImage *U = NULL;
     1076            psImage *V = NULL;
     1077            psVector *w = NULL;
     1078            if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) {
     1079                psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n");
     1080                return NULL;
     1081            }
     1082
     1083            // calculate A_inverse:
     1084            // Ainv = V * w * U^T
     1085            psImage *wUt  = p_pmSubSolve_wUt (w, U);
     1086            psImage *Ainv = p_pmSubSolve_VwUt (V, wUt);
     1087            psImage *Xvar = NULL;
     1088            psFree (wUt);
    10911089
    10921090# ifdef TESTING
    1093             // kernel terms:
    1094             for (int i = 0; i < w->n; i++) {
    1095                 fprintf (stderr, "w: %f\n", w->data.F64[i]);
    1096             }
     1091            // kernel terms:
     1092            for (int i = 0; i < w->n; i++) {
     1093                fprintf (stderr, "w: %f\n", w->data.F64[i]);
     1094            }
    10971095# endif
    1098             // loop over w adding in more and more of the values until chisquare is no longer
    1099             // dropping significantly.
    1100             // XXX this does not seem to work very well: we seem to need all terms even for
    1101             // simple cases...
    1102 
    1103             psVector *Xsvd = NULL;
    1104             {
    1105                 psVector *Ax = NULL;
    1106                 psVector *UtB = NULL;
    1107                 psVector *wUtB = NULL;
    1108 
    1109                 psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64);
    1110                 psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8);
    1111                 psVectorInit (wMask, 1); // start by masking everything
    1112 
    1113                 double chiSquareLast = NAN;
    1114                 int maxWeight = 0;
    1115 
    1116                 double Axx, Bx, y2;
    1117 
    1118                 // XXX this is an attempt to exclude insignificant modes.
    1119                 // it was not successful with the ISIS kernel set: removing even
    1120                 // the least significant mode leaves additional ringing / noise
    1121                 // because the terms are so coupled.
    1122                 for (int k = 0; false && (k < w->n); k++) {
    1123 
    1124                     // unmask the k-th weight
    1125                     wMask->data.U8[k] = 0;
    1126                     p_pmSubSolve_SetWeights(wApply, w, wMask);
    1127 
    1128                     // solve for x:
    1129                     // x = V * w * (U^T * B)
    1130                     p_pmSubSolve_UtB (&UtB, U, kernelVector);
    1131                     p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
    1132                     p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
    1133 
    1134                     // chi-square for this system of equations:
    1135                     // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
    1136                     // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
    1137                     p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
    1138                     p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
    1139                     p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
    1140                     p_pmSubSolve_y2 (&y2, kernels, stamps);
    1141 
    1142                     // apparently, this works (compare with the brute force value below
    1143                     double chiSquare = Axx - 2.0*Bx + y2;
    1144                     double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare;
    1145                     chiSquareLast = chiSquare;
    1146 
    1147                     // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi);
    1148                     if (k && !maxWeight && (deltaChi < 1.0)) {
    1149                         maxWeight = k;
    1150                     }
    1151                 }
    1152 
    1153                 // keep all terms or we get extra ringing
    1154                 maxWeight = w->n;
    1155                 psVectorInit (wMask, 1);
    1156                 for (int k = 0; k < maxWeight; k++) {
    1157                     wMask->data.U8[k] = 0;
    1158                 }
    1159                 p_pmSubSolve_SetWeights(wApply, w, wMask);
    1160 
    1161                 // solve for x:
    1162                 // x = V * w * (U^T * B)
    1163                 p_pmSubSolve_UtB (&UtB, U, kernelVector);
    1164                 p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
    1165                 p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
    1166 
    1167                 // chi-square for this system of equations:
    1168                 // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
    1169                 // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
    1170                 p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
    1171                 p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
    1172                 p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
    1173                 p_pmSubSolve_y2 (&y2, kernels, stamps);
    1174 
    1175                 // apparently, this works (compare with the brute force value below
    1176                 double chiSquare = Axx - 2.0*Bx + y2;
    1177                 psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare);
    1178 
    1179                 // re-calculate A^{-1} to get new variances:
    1180                 // Ainv = V * w * U^T
    1181                 // XXX since we keep all terms, this is identical to Ainv
    1182                 psImage *wUt  = p_pmSubSolve_wUt (wApply, U);
    1183                 Xvar = p_pmSubSolve_VwUt (V, wUt);
    1184                 psFree (wUt);
    1185 
    1186                 psFree (Ax);
    1187                 psFree (UtB);
    1188                 psFree (wUtB);
    1189                 psFree (wApply);
    1190                 psFree (wMask);
    1191             }
    1192 
    1193             // copy the kernel solutions to the full solution vector:
    1194             solution = psVectorAlloc(sumVector->n, PS_TYPE_F64);
    1195             solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
    1196 
    1197             for (int ix = 0, kx = 0; ix < sumVector->n; ix++) {
    1198                 if (ix == bgIndex) {
    1199                     solution->data.F64[ix] = 0;
    1200                     solutionErr->data.F64[ix] = 0.001;
    1201                     continue;
    1202                 }
    1203                 solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]);
    1204                 solution->data.F64[ix] = Xsvd->data.F64[kx];
    1205                 kx++;
    1206             }
    1207 
    1208             psFree (kernelMatrix);
    1209             psFree (kernelVector);
    1210 
    1211             psFree (U);
    1212             psFree (V);
    1213             psFree (w);
    1214 
    1215             psFree (Ainv);
    1216             psFree (Xsvd);
    1217         } else {
    1218             psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
    1219             psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
    1220             if (!luMatrix) {
    1221                 psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n");
    1222                 psFree(solution);
    1223                 psFree(sumVector);
    1224                 psFree(sumMatrix);
    1225                 psFree(luMatrix);
    1226                 psFree(permutation);
    1227                 return NULL;
    1228             }
    1229 
    1230             solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation);
    1231             psFree(luMatrix);
    1232             psFree(permutation);
    1233             if (!solution) {
    1234                 psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
    1235                 psFree(solution);
    1236                 psFree(sumVector);
    1237                 psFree(sumMatrix);
    1238                 return NULL;
    1239             }
    1240 
    1241             // XXX LUD does not provide A^{-1}?  fake the error for now
    1242             solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
    1243             for (int ix = 0; ix < sumVector->n; ix++) {
    1244                 solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix];
    1245             }
    1246         }
    1247 
    1248         if (!kernels->solution1) {
    1249             kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);
    1250             psVectorInit (kernels->solution1, 0.0);
    1251         }
    1252 
    1253         // only update the solutions that we chose to calculate:
    1254         if (mode & PM_SUBTRACTION_EQUATION_NORM) {
    1255             int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
    1256             kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
    1257         }
    1258         if (mode & PM_SUBTRACTION_EQUATION_BG) {
    1259             int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
    1260             kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
    1261         }
    1262         if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
    1263             int numKernels = kernels->num;
    1264             int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
    1265             int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
    1266             for (int i = 0; i < numKernels * numPoly; i++) {
    1267                 // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i]));
    1268                 if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) {
    1269                     // XXX fprintf (stderr, "drop\n");
    1270                     kernels->solution1->data.F64[i] = 0.0;
    1271                 } else {
    1272                     // XXX fprintf (stderr, "keep\n");
    1273                     kernels->solution1->data.F64[i] = solution->data.F64[i];
    1274                 }
    1275             }
    1276         }
    1277         // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps);
    1278         // fprintf (stderr, "chi square Brute = %f\n", chiSquare);
    1279 
    1280         psFree(solution);
    1281         psFree(sumVector);
    1282         psFree(sumMatrix);
     1096            // loop over w adding in more and more of the values until chisquare is no longer
     1097            // dropping significantly.
     1098            // XXX this does not seem to work very well: we seem to need all terms even for
     1099            // simple cases...
     1100
     1101            psVector *Xsvd = NULL;
     1102            {
     1103                psVector *Ax = NULL;
     1104                psVector *UtB = NULL;
     1105                psVector *wUtB = NULL;
     1106
     1107                psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64);
     1108                psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8);
     1109                psVectorInit (wMask, 1); // start by masking everything
     1110
     1111                double chiSquareLast = NAN;
     1112                int maxWeight = 0;
     1113
     1114                double Axx, Bx, y2;
     1115
     1116                // XXX this is an attempt to exclude insignificant modes.
     1117                // it was not successful with the ISIS kernel set: removing even
     1118                // the least significant mode leaves additional ringing / noise
     1119                // because the terms are so coupled.
     1120                for (int k = 0; false && (k < w->n); k++) {
     1121
     1122                    // unmask the k-th weight
     1123                    wMask->data.U8[k] = 0;
     1124                    p_pmSubSolve_SetWeights(wApply, w, wMask);
     1125
     1126                    // solve for x:
     1127                    // x = V * w * (U^T * B)
     1128                    p_pmSubSolve_UtB (&UtB, U, kernelVector);
     1129                    p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
     1130                    p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
     1131
     1132                    // chi-square for this system of equations:
     1133                    // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
     1134                    // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
     1135                    p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
     1136                    p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
     1137                    p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
     1138                    p_pmSubSolve_y2 (&y2, kernels, stamps);
     1139
     1140                    // apparently, this works (compare with the brute force value below
     1141                    double chiSquare = Axx - 2.0*Bx + y2;
     1142                    double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare;
     1143                    chiSquareLast = chiSquare;
     1144
     1145                    // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi);
     1146                    if (k && !maxWeight && (deltaChi < 1.0)) {
     1147                        maxWeight = k;
     1148                    }
     1149                }
     1150
     1151                // keep all terms or we get extra ringing
     1152                maxWeight = w->n;
     1153                psVectorInit (wMask, 1);
     1154                for (int k = 0; k < maxWeight; k++) {
     1155                    wMask->data.U8[k] = 0;
     1156                }
     1157                p_pmSubSolve_SetWeights(wApply, w, wMask);
     1158
     1159                // solve for x:
     1160                // x = V * w * (U^T * B)
     1161                p_pmSubSolve_UtB (&UtB, U, kernelVector);
     1162                p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
     1163                p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
     1164
     1165                // chi-square for this system of equations:
     1166                // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
     1167                // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
     1168                p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
     1169                p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
     1170                p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
     1171                p_pmSubSolve_y2 (&y2, kernels, stamps);
     1172
     1173                // apparently, this works (compare with the brute force value below
     1174                double chiSquare = Axx - 2.0*Bx + y2;
     1175                psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare);
     1176
     1177                // re-calculate A^{-1} to get new variances:
     1178                // Ainv = V * w * U^T
     1179                // XXX since we keep all terms, this is identical to Ainv
     1180                psImage *wUt  = p_pmSubSolve_wUt (wApply, U);
     1181                Xvar = p_pmSubSolve_VwUt (V, wUt);
     1182                psFree (wUt);
     1183
     1184                psFree (Ax);
     1185                psFree (UtB);
     1186                psFree (wUtB);
     1187                psFree (wApply);
     1188                psFree (wMask);
     1189            }
     1190
     1191            // copy the kernel solutions to the full solution vector:
     1192            solution = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     1193            solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     1194
     1195            for (int ix = 0, kx = 0; ix < sumVector->n; ix++) {
     1196                if (ix == bgIndex) {
     1197                    solution->data.F64[ix] = 0;
     1198                    solutionErr->data.F64[ix] = 0.001;
     1199                    continue;
     1200                }
     1201                solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]);
     1202                solution->data.F64[ix] = Xsvd->data.F64[kx];
     1203                kx++;
     1204            }
     1205
     1206            psFree (kernelMatrix);
     1207            psFree (kernelVector);
     1208
     1209            psFree (U);
     1210            psFree (V);
     1211            psFree (w);
     1212
     1213            psFree (Ainv);
     1214            psFree (Xsvd);
     1215        } else {
     1216            psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
     1217            psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
     1218            if (!luMatrix) {
     1219                psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n");
     1220                psFree(solution);
     1221                psFree(sumVector);
     1222                psFree(sumMatrix);
     1223                psFree(luMatrix);
     1224                psFree(permutation);
     1225                return NULL;
     1226            }
     1227
     1228            solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation);
     1229            psFree(luMatrix);
     1230            psFree(permutation);
     1231            if (!solution) {
     1232                psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
     1233                psFree(solution);
     1234                psFree(sumVector);
     1235                psFree(sumMatrix);
     1236                return NULL;
     1237            }
     1238
     1239            // XXX LUD does not provide A^{-1}?  fake the error for now
     1240            solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     1241            for (int ix = 0; ix < sumVector->n; ix++) {
     1242                solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix];
     1243            }
     1244        }
     1245
     1246        if (!kernels->solution1) {
     1247            kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);
     1248            psVectorInit (kernels->solution1, 0.0);
     1249        }
     1250
     1251        // only update the solutions that we chose to calculate:
     1252        if (mode & PM_SUBTRACTION_EQUATION_NORM) {
     1253            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     1254            kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
     1255        }
     1256        if (mode & PM_SUBTRACTION_EQUATION_BG) {
     1257            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
     1258            kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
     1259        }
     1260        if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
     1261            int numKernels = kernels->num;
     1262            int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
     1263            int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
     1264            for (int i = 0; i < numKernels * numPoly; i++) {
     1265                // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i]));
     1266                if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) {
     1267                    // XXX fprintf (stderr, "drop\n");
     1268                    kernels->solution1->data.F64[i] = 0.0;
     1269                } else {
     1270                    // XXX fprintf (stderr, "keep\n");
     1271                    kernels->solution1->data.F64[i] = solution->data.F64[i];
     1272                }
     1273            }
     1274        }
     1275        // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps);
     1276        // fprintf (stderr, "chi square Brute = %f\n", chiSquare);
     1277
     1278        psFree(solution);
     1279        psFree(sumVector);
     1280        psFree(sumMatrix);
    12831281
    12841282#ifdef TESTING
    … …  
    12991297        psImage *sumMatrixX = psImageAlloc(numParams, numParams2, PS_TYPE_F64);
    13001298        psVector *sumVector1 = psVectorAlloc(numParams, PS_TYPE_F64);
    1301         psVector *sumVector2 = psVectorAlloc(numParams, PS_TYPE_F64);
     1299        psVector *sumVector2 = psVectorAlloc(numParams2, PS_TYPE_F64);
    13021300        psImageInit(sumMatrix1, 0.0);
    13031301        psImageInit(sumMatrix2, 0.0);
    … …  
    13201318        }
    13211319
     1320
     1321#ifdef TESTING
     1322        psFitsWriteImageSimple ("sumMatrix1.fits", sumMatrix1, NULL);
     1323        psFitsWriteImageSimple ("sumMatrix2.fits", sumMatrix2, NULL);
     1324        psFitsWriteImageSimple ("sumMatrixX.fits", sumMatrixX, NULL);
     1325        psVectorWriteFile("sumVector1.dat", sumVector1);
     1326        psVectorWriteFile("sumVector2.dat", sumVector2);
     1327#endif
     1328
     1329
     1330#if 1
    13221331        int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
    13231332        calculatePenalty(sumMatrix1, sumVector1, kernels, sumMatrix1->data.F64[bgIndex][bgIndex]);
    1324         calculatePenalty(sumMatrix2, sumVector2, kernels, -sumMatrix1->data.F64[bgIndex][bgIndex]);
     1333        calculatePenalty(sumMatrix2, sumVector2, kernels, sumMatrix1->data.F64[bgIndex][bgIndex]);
     1334#endif
    13251335
    13261336        psImage *sumMatrix = NULL;      // Combined matrix
    … …  
    13301340
    13311341#ifdef TESTING
    1332         psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, NULL);
     1342        psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, NULL);
    13331343        {
    13341344            psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64);
    … …  
    13411351            psFitsClose(fits);
    13421352        }
     1353        psVectorWriteFile("sumVector.dat", sumVector);
    13431354#endif
    13441355
    13451356        psVector *solution = NULL;                       // Solution to equation!
     1357        solution = psVectorAlloc(numParams + numParams2, PS_TYPE_F64);
     1358        psVectorInit(solution, 0);
     1359
     1360#if 0
    13461361        {
    13471362            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector);
    … …  
    14231438                }
    14241439            }
    1425         }
    1426 
    1427         calculateEquationDual(&sumMatrix, &sumVector, sumMatrix1, sumMatrix2,
    1428                               sumMatrixX, sumVector1, sumVector2);
    14291440
    14301441#ifdef TESTING
    1431         {
    1432             psFits *fits = psFitsOpen("sumMatrixFix.fits", "w");
    1433             psFitsWriteImage(fits, NULL, sumMatrix, 0, NULL);
    1434             psFitsClose(fits);
    1435         }
    1436         {
    1437             psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64);
    1438             psFits *fits = psFitsOpen("sumVectorFix.fits", "w");
    1439             for (int i = 0; i < numParams + numParams2; i++) {
    1440                 image->data.F64[0][i] = sumVector->data.F64[i];
    1441             }
    1442             psFitsWriteImage(fits, NULL, image, 0, NULL);
    1443             psFree(image);
    1444             psFitsClose(fits);
    1445         }
    1446 #endif
     1442            {
     1443                psFits *fits = psFitsOpen("sumMatrixFix.fits", "w");
     1444                psFitsWriteImage(fits, NULL, sumMatrix, 0, NULL);
     1445                psFitsClose(fits);
     1446            }
     1447            {
     1448                psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64);
     1449                psFits *fits = psFitsOpen("sumVectorFix.fits", "w");
     1450                for (int i = 0; i < numParams + numParams2; i++) {
     1451                    image->data.F64[0][i] = sumVector->data.F64[i];
     1452                }
     1453                psFitsWriteImage(fits, NULL, image, 0, NULL);
     1454                psFree(image);
     1455                psFitsClose(fits);
     1456            }
     1457#endif
     1458
     1459            calculateEquationDual(&sumMatrix, &sumVector, sumMatrix1, sumMatrix2,
     1460                                  sumMatrixX, sumVector1, sumVector2);
     1461        }
     1462#endif
     1463
    14471464
    14481465        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector);
    … …  
    14631480        psFree(sumVector);
    14641481
     1482
    14651483#ifdef TESTING
    14661484        {
    … …  
    14681486            psFits *fits = psFitsOpen("solnVector.fits", "w");
    14691487            for (int i = 0; i < numParams + numParams2; i++) {
     1488                fprintf(stderr, "Solution %d: %lf\n", i, solution->data.F64[i]);
    14701489                image->data.F64[0][i] = solution->data.F64[i];
    14711490            }
    … …  
    15301549
    15311550    for (int i = 0; i < stamps->num; i++) {
    1532         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
    1533         if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
    1534             deviations->data.F32[i] = NAN;
    1535             continue;
    1536         }
    1537 
    1538         // Calculate coefficients of the kernel basis functions
    1539         polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
    1540         double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
    1541         double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
    1542 
    1543         // Calculate residuals
    1544         psKernel *weight = stamp->weight; // Weight postage stamp
    1545         psImageInit(residual->image, 0.0);
    1546         if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
    1547             psKernel *target;           // Target postage stamp
    1548             psKernel *source;           // Source postage stamp
    1549             psArray *convolutions;      // Convolution postage stamps for each kernel basis function
    1550             switch (kernels->mode) {
    1551               case PM_SUBTRACTION_MODE_1:
    1552                 target = stamp->image2;
    1553                 source = stamp->image1;
    1554                 convolutions = stamp->convolutions1;
    1555 
    1556                 // Having convolved image1 and changed its normalisation, we need to renormalise the residual
    1557                 // so that it is on the scale of image1.
    1558                 psImage *image = pmSubtractionKernelImage(kernels, stamp->xNorm, stamp->yNorm,
    1559                                                           false); // Kernel image
    1560                 if (!image) {
    1561                     psError(PS_ERR_UNKNOWN, false, "Unable to generate image of kernel.");
    1562                     return false;
    1563                 }
    1564                 double sumKernel = 0;   // Sum of kernel, for normalising residual
    1565                 int size = kernels->size; // Half-size of kernel
    1566                 int fullSize = 2 * size + 1; // Full size of kernel
    1567                 for (int y = 0; y < fullSize; y++) {
    1568                     for (int x = 0; x < fullSize; x++) {
    1569                         sumKernel += image->data.F32[y][x];
    1570                     }
    1571                 }
    1572                 psFree(image);
    1573                 devNorm = 1.0 / sumKernel / numPixels;
    1574                 break;
    1575               case PM_SUBTRACTION_MODE_2:
    1576                 target = stamp->image1;
    1577                 source = stamp->image2;
    1578                 convolutions = stamp->convolutions2;
    1579                 break;
    1580               default:
    1581                 psAbort("Unsupported subtraction mode: %x", kernels->mode);
    1582             }
    1583 
    1584             for (int j = 0; j < numKernels; j++) {
    1585                 psKernel *convolution = convolutions->data[j]; // Convolution
    1586                 double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j,
    1587                                                                   false); // Coefficient
    1588                 for (int y = - footprint; y <= footprint; y++) {
    1589                     for (int x = - footprint; x <= footprint; x++) {
    1590                         residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient;
    1591                     }
    1592                 }
    1593             }
    1594 
    1595             // XXX visualize the target, source, convolution and residual
    1596             pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i);
    1597 
    1598             for (int y = - footprint; y <= footprint; y++) {
    1599                 for (int x = - footprint; x <= footprint; x++) {
    1600                     residual->kernel[y][x] += target->kernel[y][x] - background - source->kernel[y][x] * norm;
    1601                 }
    1602             }
    1603         } else {
    1604             // Dual convolution
    1605             psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image
    1606             psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image
    1607             psKernel *image1 = stamp->image1; // The first image
    1608             psKernel *image2 = stamp->image2; // The second image
    1609 
    1610             for (int j = 0; j < numKernels; j++) {
    1611                 psKernel *conv1 = convolutions1->data[j]; // Convolution of first image
    1612                 psKernel *conv2 = convolutions2->data[j]; // Convolution of second image
    1613                 double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1
    1614                 double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2
    1615 
    1616                 for (int y = - footprint; y <= footprint; y++) {
    1617                     for (int x = - footprint; x <= footprint; x++) {
    1618                         residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 - conv1->kernel[y][x] * coeff1;
    1619                     }
    1620                 }
    1621             }
    1622 
    1623             // XXX visualize the target, source, convolution and residual
    1624             pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i);
    1625 
    1626             for (int y = - footprint; y <= footprint; y++) {
    1627                 for (int x = - footprint; x <= footprint; x++) {
    1628                     residual->kernel[y][x] += image2->kernel[y][x] - background - image1->kernel[y][x] * norm;
    1629                 }
    1630             }
    1631         }
    1632 
    1633         double deviation = 0.0;         // Sum of differences
    1634         for (int y = - footprint; y <= footprint; y++) {
    1635             for (int x = - footprint; x <= footprint; x++) {
    1636                 double dev = PS_SQR(residual->kernel[y][x]) * weight->kernel[y][x];
    1637                 deviation += dev;
     1551        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
     1552        if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
     1553            deviations->data.F32[i] = NAN;
     1554            continue;
     1555        }
     1556
     1557        // Calculate coefficients of the kernel basis functions
     1558        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
     1559        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
     1560        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
     1561
     1562        // Calculate residuals
     1563        psKernel *weight = stamp->weight; // Weight postage stamp
     1564        psImageInit(residual->image, 0.0);
     1565        if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
     1566            psKernel *target;           // Target postage stamp
     1567            psKernel *source;           // Source postage stamp
     1568            psArray *convolutions;      // Convolution postage stamps for each kernel basis function
     1569            switch (kernels->mode) {
     1570              case PM_SUBTRACTION_MODE_1:
     1571                target = stamp->image2;
     1572                source = stamp->image1;
     1573                convolutions = stamp->convolutions1;
     1574
     1575                // Having convolved image1 and changed its normalisation, we need to renormalise the residual
     1576                // so that it is on the scale of image1.
     1577                psImage *image = pmSubtractionKernelImage(kernels, stamp->xNorm, stamp->yNorm,
     1578                                                          false); // Kernel image
     1579                if (!image) {
     1580                    psError(PS_ERR_UNKNOWN, false, "Unable to generate image of kernel.");
     1581                    return false;
     1582                }
     1583                double sumKernel = 0;   // Sum of kernel, for normalising residual
     1584                int size = kernels->size; // Half-size of kernel
     1585                int fullSize = 2 * size + 1; // Full size of kernel
     1586                for (int y = 0; y < fullSize; y++) {
     1587                    for (int x = 0; x < fullSize; x++) {
     1588                        sumKernel += image->data.F32[y][x];
     1589                    }
     1590                }
     1591                psFree(image);
     1592                devNorm = 1.0 / sumKernel / numPixels;
     1593                break;
     1594              case PM_SUBTRACTION_MODE_2:
     1595                target = stamp->image1;
     1596                source = stamp->image2;
     1597                convolutions = stamp->convolutions2;
     1598                break;
     1599              default:
     1600                psAbort("Unsupported subtraction mode: %x", kernels->mode);
     1601            }
     1602
     1603            for (int j = 0; j < numKernels; j++) {
     1604                psKernel *convolution = convolutions->data[j]; // Convolution
     1605                double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j,
     1606                                                                  false); // Coefficient
     1607                for (int y = - footprint; y <= footprint; y++) {
     1608                    for (int x = - footprint; x <= footprint; x++) {
     1609                        residual->kernel[y][x] += convolution->kernel[y][x] * coefficient;
     1610                    }
     1611                }
     1612            }
     1613
     1614            // XXX visualize the target, source, convolution and residual
     1615            pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i);
     1616
     1617            for (int y = - footprint; y <= footprint; y++) {
     1618                for (int x = - footprint; x <= footprint; x++) {
     1619                    residual->kernel[y][x] += background + source->kernel[y][x] * norm - target->kernel[y][x];
     1620                }
     1621            }
     1622        } else {
     1623            // Dual convolution
     1624            psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image
     1625            psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image
     1626            psKernel *image1 = stamp->image1; // The first image
     1627            psKernel *image2 = stamp->image2; // The second image
     1628
     1629            for (int j = 0; j < numKernels; j++) {
     1630                psKernel *conv1 = convolutions1->data[j]; // Convolution of first image
     1631                psKernel *conv2 = convolutions2->data[j]; // Convolution of second image
     1632                double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1
     1633                double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2
     1634
     1635                for (int y = - footprint; y <= footprint; y++) {
     1636                    for (int x = - footprint; x <= footprint; x++) {
     1637                        residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 + conv1->kernel[y][x] * coeff1;
     1638                    }
     1639                }
     1640            }
     1641
     1642            // XXX visualize the target, source, convolution and residual
     1643            pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i);
     1644
     1645            for (int y = - footprint; y <= footprint; y++) {
     1646                for (int x = - footprint; x <= footprint; x++) {
     1647                    residual->kernel[y][x] += background + image1->kernel[y][x] * norm - image2->kernel[y][x];
     1648                }
     1649            }
     1650        }
     1651
     1652        double deviation = 0.0;         // Sum of differences
     1653        for (int y = - footprint; y <= footprint; y++) {
     1654            for (int x = - footprint; x <= footprint; x++) {
     1655                double dev = PS_SQR(residual->kernel[y][x]) * weight->kernel[y][x];
     1656                deviation += dev;
    16381657#ifdef TESTING
    1639                 residual->kernel[y][x] = dev;
    1640 #endif
    1641             }
    1642         }
    1643         deviations->data.F32[i] = devNorm * deviation;
    1644         psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n",
    1645                 i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]);
    1646         if (!isfinite(deviations->data.F32[i])) {
    1647             stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
    1648             psTrace("psModules.imcombine", 5,
    1649                     "Rejecting stamp %d (%d,%d) because of non-finite deviation\n",
    1650                     i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
    1651             continue;
    1652         }
     1658                residual->kernel[y][x] = dev;
     1659#endif
     1660            }
     1661        }
     1662        deviations->data.F32[i] = devNorm * deviation;
     1663        psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n",
     1664                i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]);
     1665        if (!isfinite(deviations->data.F32[i])) {
     1666            stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
     1667            psTrace("psModules.imcombine", 5,
     1668                    "Rejecting stamp %d (%d,%d) because of non-finite deviation\n",
     1669                    i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
     1670            continue;
     1671        }
    16531672
    16541673#ifdef TESTING
    1655         {
    1656             psString filename = NULL;
    1657             psStringAppend(&filename, "resid_%03d.fits", i);
    1658             psFits *fits = psFitsOpen(filename, "w");
    1659             psFree(filename);
    1660             psFitsWriteImage(fits, NULL, residual->image, 0, NULL);
    1661             psFitsClose(fits);
    1662         }
    1663         if (stamp->image1) {
    1664             psString filename = NULL;
    1665             psStringAppend(&filename, "stamp_image1_%03d.fits", i);
    1666             psFits *fits = psFitsOpen(filename, "w");
    1667             psFree(filename);
    1668             psFitsWriteImage(fits, NULL, stamp->image1->image, 0, NULL);
    1669             psFitsClose(fits);
    1670         }
    1671         if (stamp->image2) {
    1672             psString filename = NULL;
    1673             psStringAppend(&filename, "stamp_image2_%03d.fits", i);
    1674             psFits *fits = psFitsOpen(filename, "w");
    1675             psFree(filename);
    1676             psFitsWriteImage(fits, NULL, stamp->image2->image, 0, NULL);
    1677             psFitsClose(fits);
    1678         }
    1679         if (stamp->weight) {
    1680             psString filename = NULL;
    1681             psStringAppend(&filename, "stamp_weight_%03d.fits", i);
    1682             psFits *fits = psFitsOpen(filename, "w");
    1683             psFree(filename);
    1684             psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL);
    1685             psFitsClose(fits);
    1686         }
    1687 #endif
    1688    
     1674        {
     1675            psString filename = NULL;
     1676            psStringAppend(&filename, "resid_%03d.fits", i);
     1677            psFits *fits = psFitsOpen(filename, "w");
     1678            psFree(filename);
     1679            psFitsWriteImage(fits, NULL, residual->image, 0, NULL);
     1680            psFitsClose(fits);
     1681        }
     1682        if (stamp->image1) {
     1683            psString filename = NULL;
     1684            psStringAppend(&filename, "stamp_image1_%03d.fits", i);
     1685            psFits *fits = psFitsOpen(filename, "w");
     1686            psFree(filename);
     1687            psFitsWriteImage(fits, NULL, stamp->image1->image, 0, NULL);
     1688            psFitsClose(fits);
     1689        }
     1690        if (stamp->image2) {
     1691            psString filename = NULL;
     1692            psStringAppend(&filename, "stamp_image2_%03d.fits", i);
     1693            psFits *fits = psFitsOpen(filename, "w");
     1694            psFree(filename);
     1695            psFitsWriteImage(fits, NULL, stamp->image2->image, 0, NULL);
     1696            psFitsClose(fits);
     1697        }
     1698        if (stamp->weight) {
     1699            psString filename = NULL;
     1700            psStringAppend(&filename, "stamp_weight_%03d.fits", i);
     1701            psFits *fits = psFitsOpen(filename, "w");
     1702            psFree(filename);
     1703            psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL);
     1704            psFitsClose(fits);
     1705        }
     1706#endif
     1707
    16891708    }
    16901709
    16911710    // calculate and report the normalization and background for the image center
    1692     { 
    1693         polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, 0.0, 0.0);
    1694         double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
    1695         double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
    1696         psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f, background: %f", norm, background);
    1697 
    1698         pmSubtractionVisualShowFit(norm);
    1699         pmSubtractionVisualPlotFit(kernels);
     1711    {
     1712        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, 0.0, 0.0);
     1713        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
     1714        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
     1715        psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f, background: %f", norm, background);
     1716
     1717        pmSubtractionVisualShowFit(norm);
     1718        pmSubtractionVisualPlotFit(kernels);
    17001719    }
    17011720
    … …  
    17161735
    17171736    for (int i = 0; i < wUt->numCols; i++) {
    1718         for (int j = 0; j < wUt->numRows; j++) {
    1719             if (!isfinite(w->data.F64[j])) continue;
    1720             if (w->data.F64[j] == 0.0) continue;
    1721             wUt->data.F64[j][i] = U->data.F64[i][j] / w->data.F64[j];
    1722         }
     1737        for (int j = 0; j < wUt->numRows; j++) {
     1738            if (!isfinite(w->data.F64[j])) continue;
     1739            if (w->data.F64[j] == 0.0) continue;
     1740            wUt->data.F64[j][i] = U->data.F64[i][j] / w->data.F64[j];
     1741        }
    17231742    }
    17241743    return wUt;
    … …  
    17331752
    17341753    for (int i = 0; i < Ainv->numCols; i++) {
    1735         for (int j = 0; j < Ainv->numRows; j++) {
    1736             double sum = 0.0;
    1737             for (int k = 0; k < V->numCols; k++) {
    1738                 sum += V->data.F64[j][k] * wUt->data.F64[k][i];
    1739             }
    1740             Ainv->data.F64[j][i] = sum;
    1741         }
     1754        for (int j = 0; j < Ainv->numRows; j++) {
     1755            double sum = 0.0;
     1756            for (int k = 0; k < V->numCols; k++) {
     1757                sum += V->data.F64[j][k] * wUt->data.F64[k][i];
     1758            }
     1759            Ainv->data.F64[j][i] = sum;
     1760        }
    17421761    }
    17431762    return Ainv;
    … …  
    17521771
    17531772    for (int i = 0; i < U->numCols; i++) {
    1754         double sum = 0.0;
    1755         for (int j = 0; j < U->numRows; j++) {
    1756             sum += B->data.F64[j] * U->data.F64[j][i];
    1757         }
    1758         UtB[0]->data.F64[i] = sum;
     1773        double sum = 0.0;
     1774        for (int j = 0; j < U->numRows; j++) {
     1775            sum += B->data.F64[j] * U->data.F64[j][i];
     1776        }
     1777        UtB[0]->data.F64[i] = sum;
    17591778    }
    17601779    return true;
    … …  
    17711790
    17721791    for (int i = 0; i < w->n; i++) {
    1773         if (!isfinite(w->data.F64[i])) continue;
    1774         if (w->data.F64[i] == 0.0) continue;
    1775         wUtB[0]->data.F64[i] = UtB->data.F64[i] / w->data.F64[i];
     1792        if (!isfinite(w->data.F64[i])) continue;
     1793        if (w->data.F64[i] == 0.0) continue;
     1794        wUtB[0]->data.F64[i] = UtB->data.F64[i] / w->data.F64[i];
    17761795    }
    17771796    return true;
    … …  
    17861805
    17871806    for (int j = 0; j < V->numRows; j++) {
    1788         double sum = 0.0;
    1789         for (int i = 0; i < V->numCols; i++) {
    1790             sum += V->data.F64[j][i] * wUtB->data.F64[i];
    1791         }
    1792         VwUtB[0]->data.F64[j] = sum;
     1807        double sum = 0.0;
     1808        for (int i = 0; i < V->numCols; i++) {
     1809            sum += V->data.F64[j][i] * wUtB->data.F64[i];
     1810        }
     1811        VwUtB[0]->data.F64[j] = sum;
    17931812    }
    17941813    return true;
    … …  
    18031822
    18041823    for (int j = 0; j < A->numRows; j++) {
    1805         double sum = 0.0;
    1806         for (int i = 0; i < A->numCols; i++) {
    1807             sum += A->data.F64[j][i] * x->data.F64[i];
    1808         }
    1809         B[0]->data.F64[j] = sum;
     1824        double sum = 0.0;
     1825        for (int i = 0; i < A->numCols; i++) {
     1826            sum += A->data.F64[j][i] * x->data.F64[i];
     1827        }
     1828        B[0]->data.F64[j] = sum;
    18101829    }
    18111830    return true;
    … …  
    18191838    double sum = 0.0;
    18201839    for (int i = 0; i < x->n; i++) {
    1821         sum += x->data.F64[i] * y->data.F64[i];
     1840        sum += x->data.F64[i] * y->data.F64[i];
    18221841    }
    18231842    *value = sum;
    … …  
    18321851    for (int i = 0; i < stamps->num; i++) {
    18331852
    1834         pmSubtractionStamp *stamp = stamps->stamps->data[i];
    1835         if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
    1836 
    1837         psKernel *weight = NULL;
    1838         psKernel *window = NULL;
    1839         psKernel *input = NULL;
     1853        pmSubtractionStamp *stamp = stamps->stamps->data[i];
     1854        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
     1855
     1856        psKernel *weight = NULL;
     1857        psKernel *window = NULL;
     1858        psKernel *input = NULL;
    18401859
    18411860#ifdef USE_WEIGHT
    1842         weight = stamp->weight;
     1861        weight = stamp->weight;
    18431862#endif
    18441863#ifdef USE_WINDOW
    1845         window = stamps->window;
    1846 #endif
    1847 
    1848         switch (kernels->mode) {
    1849             // MODE_1 : convolve image 1 to match image 2 (and vice versa)
    1850           case PM_SUBTRACTION_MODE_1:
    1851             input = stamp->image2;
    1852             break;
    1853           case PM_SUBTRACTION_MODE_2:
    1854             input = stamp->image1;
    1855             break;
    1856           default:
    1857             psAbort ("programming error");
    1858         }
    1859 
    1860         for (int y = - footprint; y <= footprint; y++) {
    1861             for (int x = - footprint; x <= footprint; x++) {
    1862                 double in = input->kernel[y][x];
    1863                 double value = in*in;
    1864                 if (weight) {
    1865                     float wtVal = weight->kernel[y][x];
    1866                     value *= wtVal;
    1867                 }
    1868                 if (window) {
    1869                     float  winVal = window->kernel[y][x];
    1870                     value *= winVal;
    1871                 }
    1872                 sum += value;
    1873             }
    1874         }
     1864        window = stamps->window;
     1865#endif
     1866
     1867        switch (kernels->mode) {
     1868            // MODE_1 : convolve image 1 to match image 2 (and vice versa)
     1869          case PM_SUBTRACTION_MODE_1:
     1870            input = stamp->image2;
     1871            break;
     1872          case PM_SUBTRACTION_MODE_2:
     1873            input = stamp->image1;
     1874            break;
     1875          default:
     1876            psAbort ("programming error");
     1877        }
     1878
     1879        for (int y = - footprint; y <= footprint; y++) {
     1880            for (int x = - footprint; x <= footprint; x++) {
     1881                double in = input->kernel[y][x];
     1882                double value = in*in;
     1883                if (weight) {
     1884                    float wtVal = weight->kernel[y][x];
     1885                    value *= wtVal;
     1886                }
     1887                if (window) {
     1888                    float  winVal = window->kernel[y][x];
     1889                    value *= winVal;
     1890                }
     1891                sum += value;
     1892            }
     1893        }
    18751894    }
    18761895    *y2 = sum;
    … …  
    18921911    for (int i = 0; i < stamps->num; i++) {
    18931912
    1894         pmSubtractionStamp *stamp = stamps->stamps->data[i];
    1895         if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
    1896 
    1897         psKernel *weight = NULL;
    1898         psKernel *window = NULL;
    1899         psKernel *target = NULL;
    1900         psKernel *source = NULL;
    1901 
    1902         psArray *convolutions = NULL;
     1913        pmSubtractionStamp *stamp = stamps->stamps->data[i];
     1914        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
     1915
     1916        psKernel *weight = NULL;
     1917        psKernel *window = NULL;
     1918        psKernel *target = NULL;
     1919        psKernel *source = NULL;
     1920
     1921        psArray *convolutions = NULL;
    19031922
    19041923#ifdef USE_WEIGHT
    1905         weight = stamp->weight;
     1924        weight = stamp->weight;
    19061925#endif
    19071926#ifdef USE_WINDOW
    1908         window = stamps->window;
    1909 #endif
    1910 
    1911         switch (kernels->mode) {
    1912             // MODE_1 : convolve image 1 to match image 2 (and vice versa)
    1913           case PM_SUBTRACTION_MODE_1:
    1914             target = stamp->image2;
    1915             source = stamp->image1;
    1916             convolutions = stamp->convolutions1;
    1917             break;
    1918           case PM_SUBTRACTION_MODE_2:
    1919             target = stamp->image1;
    1920             source = stamp->image2;
    1921             convolutions = stamp->convolutions2;
    1922             break;
    1923           default:
    1924             psAbort ("programming error");
    1925         }
    1926 
    1927         // Calculate coefficients of the kernel basis functions
    1928         polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
    1929         double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
    1930         double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
    1931 
    1932         psImageInit(residual->image, 0.0);
    1933         for (int j = 0; j < numKernels; j++) {
    1934             psKernel *convolution = convolutions->data[j]; // Convolution
    1935             double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient
    1936             for (int y = - footprint; y <= footprint; y++) {
    1937                 for (int x = - footprint; x <= footprint; x++) {
    1938                     residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient;
    1939                 }
    1940             }
    1941         }
    1942 
    1943         for (int y = - footprint; y <= footprint; y++) {
    1944             for (int x = - footprint; x <= footprint; x++) {
    1945                 double resid = target->kernel[y][x] - background - source->kernel[y][x] * norm + residual->kernel[y][x];
    1946                 double value = PS_SQR(resid);
    1947                 if (weight) {
    1948                     float wtVal = weight->kernel[y][x];
    1949                     value *= wtVal;
    1950                 }
    1951                 if (window) {
    1952                     float  winVal = window->kernel[y][x];
    1953                     value *= winVal;
    1954                 }
    1955                 sum += value;
    1956             }
    1957         }
     1927        window = stamps->window;
     1928#endif
     1929
     1930        switch (kernels->mode) {
     1931            // MODE_1 : convolve image 1 to match image 2 (and vice versa)
     1932          case PM_SUBTRACTION_MODE_1:
     1933            target = stamp->image2;
     1934            source = stamp->image1;
     1935            convolutions = stamp->convolutions1;
     1936            break;
     1937          case PM_SUBTRACTION_MODE_2:
     1938            target = stamp->image1;
     1939            source = stamp->image2;
     1940            convolutions = stamp->convolutions2;
     1941            break;
     1942          default:
     1943            psAbort ("programming error");
     1944        }
     1945
     1946        // Calculate coefficients of the kernel basis functions
     1947        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
     1948        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
     1949        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
     1950
     1951        psImageInit(residual->image, 0.0);
     1952        for (int j = 0; j < numKernels; j++) {
     1953            psKernel *convolution = convolutions->data[j]; // Convolution
     1954            double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient
     1955            for (int y = - footprint; y <= footprint; y++) {
     1956                for (int x = - footprint; x <= footprint; x++) {
     1957                    residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient;
     1958                }
     1959            }
     1960        }
     1961
     1962        for (int y = - footprint; y <= footprint; y++) {
     1963            for (int x = - footprint; x <= footprint; x++) {
     1964                double resid = target->kernel[y][x] - background - source->kernel[y][x] * norm + residual->kernel[y][x];
     1965                double value = PS_SQR(resid);
     1966                if (weight) {
     1967                    float wtVal = weight->kernel[y][x];
     1968                    value *= wtVal;
     1969                }
     1970                if (window) {
     1971                    float  winVal = window->kernel[y][x];
     1972                    value *= winVal;
     1973                }
     1974                sum += value;
     1975            }
     1976        }
    19581977    }
    19591978    psFree (polyValues);
    … …  
    19661985
    19671986    for (int i = 0; i < w->n; i++) {
    1968         wApply->data.F64[i] = wMask->data.U8[i] ? 0.0 : w->data.F64[i];
     1987        wApply->data.F64[i] = wMask->data.U8[i] ? 0.0 : w->data.F64[i];
    19691988    }
    19701989    return true;
    … …  
    19812000    // generate Vn = V * w^{-1}
    19822001    for (int j = 0; j < Vn->numRows; j++) {
    1983         for (int i = 0; i < Vn->numCols; i++) {
    1984             if (!isfinite(w->data.F64[i])) continue;
    1985             if (w->data.F64[i] == 0.0) continue;
    1986             Vn->data.F64[j][i] = V->data.F64[j][i] / w->data.F64[i];
    1987         }
     2002        for (int i = 0; i < Vn->numCols; i++) {
     2003            if (!isfinite(w->data.F64[i])) continue;
     2004            if (w->data.F64[i] == 0.0) continue;
     2005            Vn->data.F64[j][i] = V->data.F64[j][i] / w->data.F64[i];
     2006        }
    19882007    }
    19892008
    … …  
    19932012    // generate Xvar = Vn * Vn^T
    19942013    for (int j = 0; j < Vn->numRows; j++) {
    1995         for (int i = 0; i < Vn->numCols; i++) {
    1996             double sum = 0.0;
    1997             for (int k = 0; k < Vn->numCols; k++) {
    1998                 sum += Vn->data.F64[k][i]*Vn->data.F64[k][j];
    1999             }
    2000             Xvar->data.F64[j][i] = sum;
    2001         }
     2014        for (int i = 0; i < Vn->numCols; i++) {
     2015            double sum = 0.0;
     2016            for (int k = 0; k < Vn->numCols; k++) {
     2017                sum += Vn->data.F64[k][i]*Vn->data.F64[k][j];
     2018            }
     2019            Xvar->data.F64[j][i] = sum;
     2020        }
    20022021    }
    20032022    return Xvar;
    … …  
    20142033    psFitsWriteImage(fits, header, image, 0, NULL);
    20152034    psFitsClose(fits);
    2016    
     2035
    20172036    return true;
    20182037}
Note: See TracChangeset for help on using the changeset viewer.