IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26547


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

Dual convolution seems to be working again!

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

Legend:

Unmodified
Added
Removed

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

    r26539 r26547  
    117117    for (int i = 0; i < numKernels; i++) {
    118118        double value = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, wantDual); // Polynomial value
     119        if (wantDual) {
     120            // The model is built with the dual convolution terms added, so to produce zero residual the
     121            // equation results in negative coefficients which we must undo.
     122            value *= -1.0;
     123        }
    119124
    120125        switch (kernels->type) {

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

    r25999 r26547  
    1616#define KERNEL_MOSAIC 2                 // Half-number of kernel instances in the mosaic image
    1717
    18 //#define TESTING
     18#define TESTING
    1919
    2020bool pmSubtractionAnalysis(psMetadata *analysis, psMetadata *header,

  • 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}

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

    r26505 r26547  
    9090
    9191    for (int i = 0, x = -size; x <= size; i++, x++) {
    92         float xf = x / sigma;
    93         float z = -0.25*xf*xf;
     92        float xf = x / sigma;
     93        float z = -0.25*xf*xf;
    9494        kernel->data.F32[i] = norm * p_pmSubtractionHermitianPolynomial(xf, order) * exp(z);
    9595    }
    … …  
    132132            kernels->preCalc->data[index] = NULL;
    133133            kernels->penalties->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v));
    134             if (!isfinite(kernels->penalties->data.F32[index])) {
    135                 psAbort ("invalid penalty");
    136             }
     134            psAssert(isfinite(kernels->penalties->data.F32[index]), "Invalid penalty");
    137135
    138136            psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v);
    … …  
    154152    double moment = 0.0;    // Moment, for penalty
    155153    for (int v = -size; v <= size; v++) {
    156         for (int u = -size; u <= size; u++) {
    157             double value = preCalc->kernel->kernel[v][u];
    158             moment += value * PS_SQR((PS_SQR(u) + PS_SQR(v)));
    159         }
     154        for (int u = -size; u <= size; u++) {
     155            double value = preCalc->kernel->kernel[v][u];
     156            moment += value * PS_SQR((PS_SQR(u) + PS_SQR(v)));
     157        }
    160158    }
    161159
    … …  
    171169
    172170    for (int v = -size; v <= size; v++) {
    173         for (int u = -size; u <= size; u++) {
    174             sum += preCalc->kernel->kernel[v][u];
    175             min = PS_MIN(preCalc->kernel->kernel[v][u], min);
    176             max = PS_MAX(preCalc->kernel->kernel[v][u], max);
    177         }
     171        for (int u = -size; u <= size; u++) {
     172            sum += preCalc->kernel->kernel[v][u];
     173            min = PS_MIN(preCalc->kernel->kernel[v][u], min);
     174            max = PS_MAX(preCalc->kernel->kernel[v][u], max);
     175        }
    178176    }
    179177#if 0
    … …  
    183181    // only even terms have non-zero sums
    184182    if ((uOrder % 2 == 0) && (vOrder % 2 == 0)) {
    185         moment /= sum;
     183        moment /= sum;
    186184    } else {
    187         moment = 0.0;
     185        moment = 0.0;
    188186    }
    189187
    … …  
    193191
    194192    if (AlardLuptonStyle && (uOrder % 2 == 0 && vOrder % 2 == 0)) {
    195         zeroNull = true;
     193        zeroNull = true;
    196194    }
    197195    if (!AlardLuptonStyle && (uOrder == 0 && vOrder == 0)) {
    198         zeroNull = true;
     196        zeroNull = true;
    199197    }
    200198    if ((uOrder % 2) || (vOrder % 2)) {
    201         // scale2D = 1.0 / (preCalc->kernel->image->numCols * preCalc->kernel->image->numRows * max);
    202         scale2D = 1.0 / max;
    203         scale1D = sqrt(scale2D);
     199        // scale2D = 1.0 / (preCalc->kernel->image->numCols * preCalc->kernel->image->numRows * max);
     200        scale2D = 1.0 / max;
     201        scale1D = sqrt(scale2D);
    204202    }
    205203
    … …  
    208206    psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "*", psScalarAlloc(scale2D, PS_TYPE_F32));
    209207    if (zeroNull) {
    210         preCalc->kernel->kernel[0][0] -= 1.0;
     208        preCalc->kernel->kernel[0][0] -= 1.0;
    211209    }
    212210
    … …  
    216214    max = FLT_MIN;
    217215    for (int v = -size; v <= size; v++) {
    218         for (int u = -size; u <= size; u++) {
    219             sum += preCalc->kernel->kernel[v][u];
    220             min = PS_MIN(preCalc->kernel->kernel[v][u], min);
    221             max = PS_MAX(preCalc->kernel->kernel[v][u], max);
    222         }
     216        for (int u = -size; u <= size; u++) {
     217            sum += preCalc->kernel->kernel[v][u];
     218            min = PS_MIN(preCalc->kernel->kernel[v][u], min);
     219            max = PS_MAX(preCalc->kernel->kernel[v][u], max);
     220        }
    223221    }
    224222    fprintf(stderr, "%d mod: %lf, null: %f, min: %lf, max: %lf, scale: %f\n", index, sum, preCalc->kernel->kernel[0][0], min, max, scale2D);
    … …  
    229227    kernels->v->data.S32[index] = vOrder;
    230228    if (kernels->preCalc->data[index]) {
    231         psFree(kernels->preCalc->data[index]);
     229        psFree(kernels->preCalc->data[index]);
    232230    }
    233231    kernels->preCalc->data[index] = preCalc;
    234232    kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
    235233    if (!isfinite(kernels->penalties->data.F32[index])) {
    236         psAbort ("invalid penalty");
     234        psAbort ("invalid penalty");
    237235    }
    238236
    239237    psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index,
    240             fwhm, uOrder, vOrder, fabsf(moment));
     238            fwhm, uOrder, vOrder, fabsf(moment));
    241239
    242240    return true;
    … …  
    259257    psVector *orders = psVectorAllocEmpty (ordersIN->n, PS_TYPE_S32);
    260258    for (int i = 0; i < fwhmsIN->n; i++) {
    261         if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
    262         psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
    263         psVectorAppend(orders, ordersIN->data.S32[i]);
     259        if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
     260        psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
     261        psVectorAppend(orders, ordersIN->data.S32[i]);
    264262    }
    265263
    … …  
    288286        for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
    289287            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    290 
    291288                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values
    292                 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true);
     289                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true);
    293290            }
    294291        }
    … …  
    299296
    300297pmSubtractionKernels *pmSubtractionKernelsHERM(int size, int spatialOrder,
    301                                                const psVector *fwhmsIN, const psVector *ordersIN,
    302                                                float penalty, pmSubtractionMode mode)
     298                                               const psVector *fwhmsIN, const psVector *ordersIN,
     299                                               float penalty, pmSubtractionMode mode)
    303300{
    304301    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
    … …  
    314311    psVector *orders = psVectorAllocEmpty (ordersIN->n, PS_TYPE_S32);
    315312    for (int i = 0; i < fwhmsIN->n; i++) {
    316         if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
    317         psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
    318         psVectorAppend(orders, ordersIN->data.S32[i]);
     313        if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
     314        psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
     315        psVectorAppend(orders, ordersIN->data.S32[i]);
    319316    }
    320317
    … …  
    342339            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    343340                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values
    344                 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true);
     341                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true);
    345342            }
    346343        }
    … …  
    351348
    352349pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, int spatialOrder,
    353                                                      const psVector *fwhmsIN, const psVector *ordersIN,
    354                                                      float penalty, pmSubtractionMode mode)
     350                                                     const psVector *fwhmsIN, const psVector *ordersIN,
     351                                                     float penalty, pmSubtractionMode mode)
    355352{
    356353    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
    … …  
    366363    psVector *orders = psVectorAllocEmpty (ordersIN->n, PS_TYPE_S32);
    367364    for (int i = 0; i < fwhmsIN->n; i++) {
    368         if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
    369         psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
    370         psVectorAppend(orders, ordersIN->data.S32[i]);
     365        if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
     366        psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
     367        psVectorAppend(orders, ordersIN->data.S32[i]);
    371368    }
    372369
    … …  
    405402                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values
    406403
    407                 // save the generated 2D kernel as the target, deconvolve it by Gaussian, replacing the generated 2D kernel
    408                 psKernel *kernelTarget = preCalc->kernel;
     404                // save the generated 2D kernel as the target, deconvolve it by Gaussian, replacing the generated 2D kernel
     405                psKernel *kernelTarget = preCalc->kernel;
    409406                preCalc->kernel = pmSubtractionDeconvolveKernel(kernelTarget, kernelGauss); // Kernel
    410407
    411                 // XXX do we use Alard-Lupton normalization (last param true) or not?
    412                 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true);
    413 
    414                 // XXXX test demo that deconvolved kernel is valid
     408                // XXX do we use Alard-Lupton normalization (last param true) or not?
     409                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true);
     410
     411                // XXXX test demo that deconvolved kernel is valid
    415412# if 1
    416                 psImage *kernelConv = psImageConvolveFFT(NULL, preCalc->kernel->image, NULL, 0, kernelGauss);
    417                 psArrayAdd (deconKernels, 100, kernelConv);
    418                 psFree (kernelConv);
    419 
    420                 if (!uOrder && !vOrder){
    421                     pmSubtractionVisualShowSubtraction (kernelTarget->image, preCalc->kernel->image, kernelConv);
    422                 }
     413                psImage *kernelConv = psImageConvolveFFT(NULL, preCalc->kernel->image, NULL, 0, kernelGauss);
     414                psArrayAdd (deconKernels, 100, kernelConv);
     415                psFree (kernelConv);
     416
     417                if (!uOrder && !vOrder){
     418                    pmSubtractionVisualShowSubtraction (kernelTarget->image, preCalc->kernel->image, kernelConv);
     419                }
    423420# endif
    424421            }
    … …  
    429426    psImage *dot = psImageAlloc(deconKernels->n, deconKernels->n, PS_TYPE_F32);
    430427    for (int i = 0; i < deconKernels->n; i++) {
    431         for (int j = 0; j <= i; j++) {
    432             psImage *t1 = deconKernels->data[i];
    433             psImage *t2 = deconKernels->data[j];
    434 
    435             double sum = 0.0;
    436             for (int iy = 0; iy < t1->numRows; iy++) {
    437                 for (int ix = 0; ix < t1->numCols; ix++) {
    438                     sum += t1->data.F32[iy][ix] * t2->data.F32[iy][ix];
    439                 }
    440             }
    441             dot->data.F32[j][i] = sum;
    442             dot->data.F32[i][j] = sum;
    443         }
     428        for (int j = 0; j <= i; j++) {
     429            psImage *t1 = deconKernels->data[i];
     430            psImage *t2 = deconKernels->data[j];
     431
     432            double sum = 0.0;
     433            for (int iy = 0; iy < t1->numRows; iy++) {
     434                for (int ix = 0; ix < t1->numCols; ix++) {
     435                    sum += t1->data.F32[iy][ix] * t2->data.F32[iy][ix];
     436                }
     437            }
     438            dot->data.F32[j][i] = sum;
     439            dot->data.F32[i][j] = sum;
     440        }
    444441    }
    445442    pmSubtractionVisualShowSubtraction (dot, NULL, NULL);
    … …  
    494491    switch (type) {
    495492      case PM_SUBTRACTION_KERNEL_ISIS:
    496         preCalc->xKernel = pmSubtractionKernelISIS(sigma, uOrder, size);
    497         preCalc->yKernel = pmSubtractionKernelISIS(sigma, vOrder, size);
    498         preCalc->uCoords = NULL;
    499         preCalc->vCoords = NULL;
    500         preCalc->poly    = NULL;
    501         break;
     493        preCalc->xKernel = pmSubtractionKernelISIS(sigma, uOrder, size);
     494        preCalc->yKernel = pmSubtractionKernelISIS(sigma, vOrder, size);
     495        preCalc->uCoords = NULL;
     496        preCalc->vCoords = NULL;
     497        preCalc->poly    = NULL;
     498        break;
    502499      case PM_SUBTRACTION_KERNEL_HERM:
    503         preCalc->xKernel = pmSubtractionKernelHERM(sigma, uOrder, size);
    504         preCalc->yKernel = pmSubtractionKernelHERM(sigma, vOrder, size);
    505         preCalc->uCoords = NULL;
    506         preCalc->vCoords = NULL;
    507         preCalc->poly    = NULL;
    508         break;
     500        preCalc->xKernel = pmSubtractionKernelHERM(sigma, uOrder, size);
     501        preCalc->yKernel = pmSubtractionKernelHERM(sigma, vOrder, size);
     502        preCalc->uCoords = NULL;
     503        preCalc->vCoords = NULL;
     504        preCalc->poly    = NULL;
     505        break;
    509506      case PM_SUBTRACTION_KERNEL_DECONV_HERM:
    510         preCalc->xKernel = pmSubtractionKernelHERM(sigma, uOrder, size);
    511         preCalc->yKernel = pmSubtractionKernelHERM(sigma, vOrder, size);
    512         preCalc->uCoords = NULL;
    513         preCalc->vCoords = NULL;
    514         preCalc->poly    = NULL;
    515         break;
     507        preCalc->xKernel = pmSubtractionKernelHERM(sigma, uOrder, size);
     508        preCalc->yKernel = pmSubtractionKernelHERM(sigma, vOrder, size);
     509        preCalc->uCoords = NULL;
     510        preCalc->vCoords = NULL;
     511        preCalc->poly    = NULL;
     512        break;
    516513      case PM_SUBTRACTION_KERNEL_RINGS:
    517         // the RINGS kernel uses the uCoords, vCoords, and poly elements of the structure
    518         // we allocate these vectors here, but leave the kernel generation to the main function
    519         preCalc->xKernel = NULL;
    520         preCalc->yKernel = NULL;
    521         preCalc->kernel  = NULL;
    522         preCalc->uCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // u coords
    523         preCalc->vCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // v coords
    524         preCalc->poly    = psVectorAllocEmpty(size, PS_TYPE_F32); // Polynomial
    525         return preCalc;
     514        // the RINGS kernel uses the uCoords, vCoords, and poly elements of the structure
     515        // we allocate these vectors here, but leave the kernel generation to the main function
     516        preCalc->xKernel = NULL;
     517        preCalc->yKernel = NULL;
     518        preCalc->kernel  = NULL;
     519        preCalc->uCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // u coords
     520        preCalc->vCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // v coords
     521        preCalc->poly    = psVectorAllocEmpty(size, PS_TYPE_F32); // Polynomial
     522        return preCalc;
    526523      default:
    527         psAbort("programming error: invalid type for PreCalc kernel");
     524        psAbort("programming error: invalid type for PreCalc kernel");
    528525    }
    529526
    … …  
    532529    // generate 2D kernel from 1D realizations
    533530    for (int v = -size, y = 0; v <= size; v++, y++) {
    534         for (int u = -size, x = 0; u <= size; u++, x++) {
    535             preCalc->kernel->kernel[v][u] = preCalc->xKernel->data.F32[x] * preCalc->yKernel->data.F32[y]; // Value of kernel
    536         }
    537     }
    538    
     531        for (int u = -size, x = 0; u <= size; u++, x++) {
     532            preCalc->kernel->kernel[v][u] = preCalc->xKernel->data.F32[x] * preCalc->yKernel->data.F32[y]; // Value of kernel
     533        }
     534    }
     535
    539536    return preCalc;
    540537}
    … …  
    864861            for (int vOrder = 0; vOrder <= (i == 0 ? 0 : ringsOrder - uOrder); vOrder++, index++) {
    865862
    866                 pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc (PM_SUBTRACTION_KERNEL_RINGS, 0, 0, RINGS_BUFFER, 0.0);
     863                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc (PM_SUBTRACTION_KERNEL_RINGS, 0, 0, RINGS_BUFFER, 0.0);
    867864                double moment = 0.0;    // Moment, for penalty
    868865
    … …  
    870867                    // Central pixel is easy
    871868                    preCalc->uCoords->data.S32[0] = 0;
    872                     preCalc->vCoords->data.S32[0] = 0;
     869                    preCalc->vCoords->data.S32[0] = 0;
    873870                    preCalc->poly->data.F32[0] = 1.0;
    874871                    preCalc->uCoords->n = 1;
    875                     preCalc->vCoords->n = 1;
    876                     preCalc->poly->n = 1;
     872                    preCalc->vCoords->n = 1;
     873                    preCalc->poly->n = 1;
    877874                    radiusLast = 0;
    878875                    moment = 0.0;
    … …  
    931928                kernels->v->data.S32[index] = vOrder;
    932929                kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
    933                 if (!isfinite(kernels->penalties->data.F32[index])) {
    934                     psAbort ("invalid penalty");
    935                 }
     930                if (!isfinite(kernels->penalties->data.F32[index])) {
     931                    psAbort ("invalid penalty");
     932                }
    936933
    937934                psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d\n", index,
    … …  
    10261023            type = PM_SUBTRACTION_KERNEL_GUNK;
    10271024            psAbort("Deciphering GUNK kernels (%s) is not currently supported.", description);
    1028         } 
    1029 
    1030         type = pmSubtractionKernelsTypeFromString (description);
    1031         psAssert (type != PM_SUBTRACTION_KERNEL_NONE, "must  be ISIS, HERM or DECONV_HERM");
    1032 
    1033         char *ptr = NULL;
    1034         switch (type) {
    1035           case PM_SUBTRACTION_KERNEL_ISIS:
    1036           case PM_SUBTRACTION_KERNEL_HERM:
    1037             ptr = (char*) description + 5;    // Eat "ISIS(" or "HERM("
    1038             break;
    1039           case PM_SUBTRACTION_KERNEL_DECONV_HERM:
    1040             ptr = (char*) description + 12;    // Eat "DECONV_HERM("
    1041             break;
    1042           default:
    1043             psAbort("programming error: invalid kernel type");
    1044         }
    1045         PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt);
    1046 
    1047         // Count the number of Gaussians
    1048         int numGauss = 0;
    1049         for (char *string = ptr; string; string = strchr(string + 1, '(')) {
    1050             numGauss++;
    1051         }
    1052 
    1053         fwhms = psVectorAlloc(numGauss, PS_TYPE_F32);
    1054         orders = psVectorAlloc(numGauss, PS_TYPE_S32);
    1055 
    1056         for (int i = 0; i < numGauss; i++) {
    1057             ptr++;                  // Eat the '('
    1058             PARSE_STRING_NUMBER(fwhms->data.F32[i], ptr, ',', parseStringFloat); // Eat "1.234,"
    1059             PARSE_STRING_NUMBER(orders->data.S32[i], ptr, ')', parseStringInt); // Eat "3)"
    1060         }
    1061 
    1062         ptr++;                      // Eat ','
    1063         PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
    1064         penalty = parseStringFloat(ptr);
    1065 
    1066         return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
    1067     } 
     1025        }
     1026
     1027        type = pmSubtractionKernelsTypeFromString (description);
     1028        psAssert (type != PM_SUBTRACTION_KERNEL_NONE, "must  be ISIS, HERM or DECONV_HERM");
     1029
     1030        char *ptr = NULL;
     1031        switch (type) {
     1032          case PM_SUBTRACTION_KERNEL_ISIS:
     1033          case PM_SUBTRACTION_KERNEL_HERM:
     1034            ptr = (char*) description + 5;    // Eat "ISIS(" or "HERM("
     1035            break;
     1036          case PM_SUBTRACTION_KERNEL_DECONV_HERM:
     1037            ptr = (char*) description + 12;    // Eat "DECONV_HERM("
     1038            break;
     1039          default:
     1040            psAbort("programming error: invalid kernel type");
     1041        }
     1042        PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt);
     1043
     1044        // Count the number of Gaussians
     1045        int numGauss = 0;
     1046        for (char *string = ptr; string; string = strchr(string + 1, '(')) {
     1047            numGauss++;
     1048        }
     1049
     1050        fwhms = psVectorAlloc(numGauss, PS_TYPE_F32);
     1051        orders = psVectorAlloc(numGauss, PS_TYPE_S32);
     1052
     1053        for (int i = 0; i < numGauss; i++) {
     1054            ptr++;                  // Eat the '('
     1055            PARSE_STRING_NUMBER(fwhms->data.F32[i], ptr, ',', parseStringFloat); // Eat "1.234,"
     1056            PARSE_STRING_NUMBER(orders->data.S32[i], ptr, ')', parseStringInt); // Eat "3)"
     1057        }
     1058
     1059        ptr++;                      // Eat ','
     1060        PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
     1061        penalty = parseStringFloat(ptr);
     1062
     1063        return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
     1064    }
    10681065
    10691066    if (strncmp(description, "RINGS", 5) == 0) {
    … …  
    10751072        PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
    10761073        PARSE_STRING_NUMBER(penalty, ptr, ')', parseStringInt);
    1077         return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
    1078     } 
     1074        return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
     1075    }
    10791076
    10801077    psAbort("Deciphering kernels other than ISIS, HERM, DECONV_HERM or RINGS is not currently supported.");

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

    r26505 r26547  
    242242
    243243    for (int i = 0; i < kernels->num; i++) {
    244         pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i];
    245         psKernel *kernel = preCalc->kernel;
    246 
    247         int xSub = i % NXsub;
    248         int ySub = i / NXsub;
    249 
    250         int xPix = xSub * (2*footprint + 1 + 3) + footprint;
    251         int yPix = ySub * (2*footprint + 1 + 3) + footprint;
    252 
    253         double sum = 0.0;
    254         for (int y = -footprint; y <= footprint; y++) {
    255             for (int x = -footprint; x <= footprint; x++) {
    256                 output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
    257                 sum += kernel->kernel[y][x];
    258             }
    259         }
    260         fprintf (stderr, "kernel %d, sum %f\n", i, sum);
    261     }                                                   
    262        
     244        pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i];
     245        psKernel *kernel = preCalc->kernel;
     246
     247        int xSub = i % NXsub;
     248        int ySub = i / NXsub;
     249
     250        int xPix = xSub * (2*footprint + 1 + 3) + footprint;
     251        int yPix = ySub * (2*footprint + 1 + 3) + footprint;
     252
     253        double sum = 0.0;
     254        for (int y = -footprint; y <= footprint; y++) {
     255            for (int x = -footprint; x <= footprint; x++) {
     256                output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
     257                sum += kernel->kernel[y][x];
     258            }
     259        }
     260        fprintf (stderr, "kernel %d, sum %f\n", i, sum);
     261    }
     262
    263263    pmVisualScaleImage(kapa1, output, "Image", 0, true);
    264264    pmVisualAskUser(&plotImage);
    … …  
    280280    float maxFlux = NAN;
    281281    for (int i = 0; i < stamps->num; i++) {
    282         pmSubtractionStamp *stamp = stamps->stamps->data[i];
    283         if (!isfinite(stamp->flux)) continue;
    284         if (!stamp->convolutions1 && !stamp->convolutions2) continue;
    285         if (!maxStamp) {
    286             maxFlux = stamp->flux;
    287             maxStamp = stamp;
    288             continue;
    289         }
    290         if (stamp->flux > maxFlux) {
    291             maxFlux = stamp->flux;
    292             maxStamp = stamp;
    293         }
     282        pmSubtractionStamp *stamp = stamps->stamps->data[i];
     283        if (!isfinite(stamp->flux)) continue;
     284        if (!stamp->convolutions1 && !stamp->convolutions2) continue;
     285        if (!maxStamp) {
     286            maxFlux = stamp->flux;
     287            maxStamp = stamp;
     288            continue;
     289        }
     290        if (stamp->flux > maxFlux) {
     291            maxFlux = stamp->flux;
     292            maxStamp = stamp;
     293        }
    294294    }
    295295
    296296    if (!isfinite(maxStamp->flux)) {
    297         fprintf (stderr, "no valid stamps?\n");
     297        fprintf (stderr, "no valid stamps?\n");
    298298    }
    299299
    … …  
    301301
    302302    if (maxStamp->convolutions1) {
    303         // output image is a grid of NXsub by NYsub sub-images
    304         nKernels = maxStamp->convolutions1->n;
    305         int NXsub = sqrt(nKernels);
    306         int NYsub = nKernels / NXsub;
    307         if (nKernels % NXsub) NYsub++;
    308 
    309         int NXpix = NXsub * (2*footprint + 1 + 3);
    310         int NYpix = NYsub * (2*footprint + 1 + 3);
    311 
    312         psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32);
    313         psImageInit (output, 0.0);
    314 
    315         for (int i = 0; i < nKernels; i++) {
     303        // output image is a grid of NXsub by NYsub sub-images
     304        nKernels = maxStamp->convolutions1->n;
     305        int NXsub = sqrt(nKernels);
     306        int NYsub = nKernels / NXsub;
     307        if (nKernels % NXsub) NYsub++;
     308
     309        int NXpix = NXsub * (2*footprint + 1 + 3);
     310        int NYpix = NYsub * (2*footprint + 1 + 3);
     311
     312        psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32);
     313        psImageInit (output, 0.0);
     314
     315        for (int i = 0; i < nKernels; i++) {
    316316            psKernel *kernel = maxStamp->convolutions1->data[i];
    317            
    318             int xSub = i % NXsub;
    319             int ySub = i / NXsub;
    320            
    321             int xPix = xSub * (2*footprint + 1 + 3) + footprint;
    322             int yPix = ySub * (2*footprint + 1 + 3) + footprint;
    323            
    324             double sum = 0.0;
    325             for (int y = -footprint; y <= footprint; y++) {
    326                 for (int x = -footprint; x <= footprint; x++) {
    327                     output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
    328                     sum += kernel->kernel[y][x];
    329                 }
    330             }
    331             fprintf (stderr, "kernel %d, sum %f\n", i, sum);
    332         }               
    333         pmVisualScaleImage(kapa2, output, "Image", 0, true);
    334     }                                   
    335        
     317
     318            int xSub = i % NXsub;
     319            int ySub = i / NXsub;
     320
     321            int xPix = xSub * (2*footprint + 1 + 3) + footprint;
     322            int yPix = ySub * (2*footprint + 1 + 3) + footprint;
     323
     324            double sum = 0.0;
     325            for (int y = -footprint; y <= footprint; y++) {
     326                for (int x = -footprint; x <= footprint; x++) {
     327                    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
     328                    sum += kernel->kernel[y][x];
     329                }
     330            }
     331            fprintf (stderr, "kernel %d, sum %f\n", i, sum);
     332        }
     333        pmVisualScaleImage(kapa2, output, "Image", 0, true);
     334    }
     335
    336336    if (maxStamp->convolutions2) {
    337         // output image is a grid of NXsub by NYsub sub-images
    338         nKernels = maxStamp->convolutions2->n;
    339         int NXsub = sqrt(nKernels);
    340         int NYsub = nKernels / NXsub;
    341         if (nKernels % NXsub) NYsub++;
    342 
    343         int NXpix = NXsub * (2*footprint + 1 + 3);
    344         int NYpix = NYsub * (2*footprint + 1 + 3);
    345 
    346         psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32);
    347         psImageInit (output, 0.0);
    348 
    349         for (int i = 0; i < nKernels; i++) {
     337        // output image is a grid of NXsub by NYsub sub-images
     338        nKernels = maxStamp->convolutions2->n;
     339        int NXsub = sqrt(nKernels);
     340        int NYsub = nKernels / NXsub;
     341        if (nKernels % NXsub) NYsub++;
     342
     343        int NXpix = NXsub * (2*footprint + 1 + 3);
     344        int NYpix = NYsub * (2*footprint + 1 + 3);
     345
     346        psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32);
     347        psImageInit (output, 0.0);
     348
     349        for (int i = 0; i < nKernels; i++) {
    350350            psKernel *kernel = maxStamp->convolutions2->data[i];
    351            
    352             int xSub = i % NXsub;
    353             int ySub = i / NXsub;
    354            
    355             int xPix = xSub * (2*footprint + 1 + 3) + footprint;
    356             int yPix = ySub * (2*footprint + 1 + 3) + footprint;
    357            
    358             double sum = 0.0;
    359             for (int y = -footprint; y <= footprint; y++) {
    360                 for (int x = -footprint; x <= footprint; x++) {
    361                     output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
    362                     sum += kernel->kernel[y][x];
    363                 }
    364             }
    365             fprintf (stderr, "kernel %d, sum %f\n", i, sum);
    366         }               
    367         pmVisualScaleImage(kapa2, output, "Image", 1, true);
    368     }                                   
    369        
     351
     352            int xSub = i % NXsub;
     353            int ySub = i / NXsub;
     354
     355            int xPix = xSub * (2*footprint + 1 + 3) + footprint;
     356            int yPix = ySub * (2*footprint + 1 + 3) + footprint;
     357
     358            double sum = 0.0;
     359            for (int y = -footprint; y <= footprint; y++) {
     360                for (int x = -footprint; x <= footprint; x++) {
     361                    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
     362                    sum += kernel->kernel[y][x];
     363                }
     364            }
     365            fprintf (stderr, "kernel %d, sum %f\n", i, sum);
     366        }
     367        pmVisualScaleImage(kapa2, output, "Image", 1, true);
     368    }
     369
    370370    pmVisualAskUser(&plotImage);
    371371    return true;
    … …  
    394394
    395395        overlay[Noverlay].type = KII_OVERLAY_BOX;
    396         if ((stamp->x < 1.0) && (stamp->y < 1.0)) {
    397             // fprintf (stderr, "stamp zero: %f %f\n", stamp->x, stamp->y);
    398             continue;
    399         }
    400         if (!isfinite(stamp->x) && !isfinite(stamp->y)) {
    401             // fprintf (stderr, "stamp nan: %f %f\n", stamp->x, stamp->y);
    402             continue;
    403         }
     396        if ((stamp->x < 1.0) && (stamp->y < 1.0)) {
     397            // fprintf (stderr, "stamp zero: %f %f\n", stamp->x, stamp->y);
     398            continue;
     399        }
     400        if (!isfinite(stamp->x) && !isfinite(stamp->y)) {
     401            // fprintf (stderr, "stamp nan: %f %f\n", stamp->x, stamp->y);
     402            continue;
     403        }
    404404        overlay[Noverlay].x = stamp->x;
    405405        overlay[Noverlay].y = stamp->y;
    … …  
    437437    float NXf = sqrt(stamps->num);
    438438    NX = (int) NXf == NXf ? NXf : NXf + 1.0;
    439    
     439
    440440    float NYf = stamps->num / NX;
    441441    NY = (int) NYf == NY ? NYf : NYf + 1.0;
    … …  
    453453    differenceImage  = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
    454454    convolutionImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
    455    
     455
    456456    psImageInit (sourceImage,      0.0);
    457457    psImageInit (targetImage,      0.0);
    … …  
    479479    sum = 0.0;
    480480    for (int y = -footprint; y <= footprint; y++) {
    481         for (int x = -footprint; x <= footprint; x++) {
    482             targetImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x];
    483             sum += targetImage->data.F32[y + NYpix][x + NXpix];
    484         }
     481        for (int x = -footprint; x <= footprint; x++) {
     482            targetImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x];
     483            sum += targetImage->data.F32[y + NYpix][x + NXpix];
     484        }
    485485    }
    486486    targetImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
    … …  
    489489    sum = 0.0;
    490490    for (int y = -footprint; y <= footprint; y++) {
    491         for (int x = -footprint; x <= footprint; x++) {
    492             sourceImage->data.F32[y + NYpix][x + NXpix] = source->kernel[y][x];
    493             sum += sourceImage->data.F32[y + NYpix][x + NXpix];
    494         }
     491        for (int x = -footprint; x <= footprint; x++) {
     492            sourceImage->data.F32[y + NYpix][x + NXpix] = source->kernel[y][x];
     493            sum += sourceImage->data.F32[y + NYpix][x + NXpix];
     494        }
    495495    }
    496496    sourceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
    … …  
    499499    sum = 0.0;
    500500    for (int y = -footprint; y <= footprint; y++) {
    501         for (int x = -footprint; x <= footprint; x++) {
    502             convolutionImage->data.F32[y + NYpix][x + NXpix] = -convolution->kernel[y][x];
    503             sum += convolutionImage->data.F32[y + NYpix][x + NXpix];
    504         }
     501        for (int x = -footprint; x <= footprint; x++) {
     502            convolutionImage->data.F32[y + NYpix][x + NXpix] = -convolution->kernel[y][x];
     503            sum += convolutionImage->data.F32[y + NYpix][x + NXpix];
     504        }
    505505    }
    506506    convolutionImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
    507    
     507
    508508    // insert the (difference) kernel into the (difference) image:
    509509    sum = 0.0;
    510510    for (int y = -footprint; y <= footprint; y++) {
    511         for (int x = -footprint; x <= footprint; x++) {
    512             differenceImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm;
    513             sum += differenceImage->data.F32[y + NYpix][x + NXpix];
    514         }
     511        for (int x = -footprint; x <= footprint; x++) {
     512            differenceImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm;
     513            sum += differenceImage->data.F32[y + NYpix][x + NXpix];
     514        }
    515515    }
    516516    differenceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
    … …  
    519519    sum = 0.0;
    520520    for (int y = -footprint; y <= footprint; y++) {
    521         for (int x = -footprint; x <= footprint; x++) {
    522             residualImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm + convolution->kernel[y][x];
    523             sum += residualImage->data.F32[y + NYpix][x + NXpix];
    524         }
     521        for (int x = -footprint; x <= footprint; x++) {
     522            residualImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm + convolution->kernel[y][x];
     523            sum += residualImage->data.F32[y + NYpix][x + NXpix];
     524        }
    525525    }
    526526    residualImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
    … …  
    528528    // insert the (fresidual) kernel into the (fresidual) image:
    529529    for (int y = -footprint; y <= footprint; y++) {
    530         for (int x = -footprint; x <= footprint; x++) {
    531             fresidualImage->data.F32[y + NYpix][x + NXpix] = residualImage->data.F32[y + NYpix][x + NXpix] / sqrt(PS_MAX(target->kernel[y][x], 100.0));
    532         }
     530        for (int x = -footprint; x <= footprint; x++) {
     531            fresidualImage->data.F32[y + NYpix][x + NXpix] = residualImage->data.F32[y + NYpix][x + NXpix] / sqrt(PS_MAX(target->kernel[y][x], 100.0));
     532        }
    533533    }
    534534    return true;
    … …  
    536536
    537537bool pmSubtractionVisualShowFit(double norm) {
     538
     539    if (!pmVisualIsVisual()) return true;
    538540
    539541    // XXX a dumb test : dump the residual image and exit
    540542    {
    541         psMetadata *header = psMetadataAlloc();
     543        psMetadata *header = psMetadataAlloc();
    542544        psMetadataAddF32 (header, PS_LIST_TAIL, "NORM", 0, "Normalization", norm);
    543         psFits *fits = psFitsOpen("resid.fits", "w");
    544         psFitsWriteImage(fits, header, residualImage, 0, NULL);
    545         psFitsClose(fits);
    546         exit (0);
    547     }
    548 
    549     // if (!pmVisualIsVisual()) return true;
     545        psFits *fits = psFitsOpen("resid.fits", "w");
     546        psFitsWriteImage(fits, header, residualImage, 0, NULL);
     547        psFitsClose(fits);
     548        exit (0);
     549    }
     550
    550551    if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) return false;
    551552    if (!pmVisualInitWindow(&kapa2, "ppSub:Misc")) return false;
    … …  
    606607    for (int i = 0; i < kernels->num; i++) {
    607608        x->data.F32[i] = i;
    608         y->data.F32[i] = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, false);
     609        y->data.F32[i] = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, false);
    609610        graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[i]);
    610611        graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[i]);
Note: See TracChangeset for help on using the changeset viewer.