IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26576


Ignore:
Timestamp:
Jan 12, 2010, 4:30:30 PM (16 years ago)
Author:
eugene
Message:

measure new residual stats; use WINDOW and WEIGHT

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r26562 r26576  
    1717// #define TESTING                         // TESTING output for debugging; may not work with threads!
    1818
    19 // #define USE_WEIGHT                      // Include weight (1/variance) in equation?
    20 // #define USE_WINDOW                      // Include weight (1/variance) in equation?
     19#define USE_WEIGHT                      // Include weight (1/variance) in equation?
     20#define USE_WINDOW                      // Include weight (1/variance) in equation?
    2121
    2222
    … …  
    3636                                  const psImage *polyValues, // Spatial polynomial values
    3737                                  int footprint, // (Half-)Size of stamp
    38                                   const pmSubtractionEquationCalculationMode mode
     38                                  const pmSubtractionEquationCalculationMode mode
    3939                                  )
    4040{
    … …  
    6868    }
    6969
    70     // initialize the matrix and vector for NOP on all coeffs.  we only fill in the coeffs we
     70    // initialize the matrix and vector for NOP on all coeffs.  we only fill in the coeffs we 
    7171    // choose to calculate
    7272    psImageInit(matrix, 0.0);
    7373    psVectorInit(vector, 1.0);
    7474    for (int i = 0; i < matrix->numCols; i++) {
    75         matrix->data.F64[i][i] = 1.0;
     75        matrix->data.F64[i][i] = 1.0;
    7676    }
    7777
    … …  
    8484
    8585    for (int i = 0; i < numKernels; i++) {
    86         psKernel *iConv = convolutions->data[i]; // Convolution for index i
    87         for (int j = i; j < numKernels; j++) {
    88             psKernel *jConv = convolutions->data[j]; // Convolution for index j
    89 
    90             double sumCC = 0.0;         // Sum of convolution products
    91             for (int y = - footprint; y <= footprint; y++) {
    92                 for (int x = - footprint; x <= footprint; x++) {
    93                     double cc = iConv->kernel[y][x] * jConv->kernel[y][x];
    94                     if (weight) {
    95                         cc *= weight->kernel[y][x];
    96                     }
    97                     if (window) {
    98                         cc *= window->kernel[y][x];
    99                     }
    100                     sumCC += cc;
    101                 }
    102             }
    103 
    104             // Spatial variation of kernel coeffs
    105             if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
    106                 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
    107                     for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
    108                         double value = sumCC * poly2[iTerm][jTerm];
    109                         matrix->data.F64[iIndex][jIndex] = value;
    110                         matrix->data.F64[jIndex][iIndex] = value;
    111                     }
    112                 }
    113             }
    114         }
    115 
    116         double sumRC = 0.0;             // Sum of the reference-convolution products
    117         double sumIC = 0.0;             // Sum of the input-convolution products
    118         double sumC = 0.0;              // Sum of the convolution
    119         for (int y = - footprint; y <= footprint; y++) {
    120             for (int x = - footprint; x <= footprint; x++) {
    121                 float conv = iConv->kernel[y][x];
    122                 float in = input->kernel[y][x];
    123                 float ref = reference->kernel[y][x];
    124                 double ic = in * conv;
    125                 double rc = ref * conv;
    126                 double c = conv;
    127                 if (weight) {
    128                     float wtVal = weight->kernel[y][x];
    129                     ic *= wtVal;
    130                     rc *= wtVal;
    131                     c *= wtVal;
    132                 }
    133                 if (window) {
    134                     float winVal = window->kernel[y][x];
    135                     ic *= winVal;
    136                     rc *= winVal;
    137                     c  *= winVal;
    138                 }
    139                 sumIC += ic;
    140                 sumRC += rc;
    141                 sumC += c;
    142             }
    143         }
    144         // Spatial variation
    145         for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
    146             double normTerm = sumRC * poly[iTerm];
    147             double bgTerm = sumC * poly[iTerm];
    148             if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
    149                 matrix->data.F64[iIndex][normIndex] = normTerm;
    150                 matrix->data.F64[normIndex][iIndex] = normTerm;
    151             }
    152             if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
    153                 matrix->data.F64[iIndex][bgIndex] = bgTerm;
    154                 matrix->data.F64[bgIndex][iIndex] = bgTerm;
    155             }
    156             if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
    157                 vector->data.F64[iIndex] = sumIC * poly[iTerm];
    158                 // XXX TEST for Hermitians: do not calculate A - norm*B - \sum(k x B),
    159                 // instead, calculate A - \sum(k x B), with full hermitians
    160                 if (0 && !(mode & PM_SUBTRACTION_EQUATION_NORM)) {
    161                     // subtract norm * sumRC * poly[iTerm]
    162                     psAssert (kernels->solution1, "programming error: define solution first!");
    163                     int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
    164                     double norm = fabs(kernels->solution1->data.F64[normIndex]);  // Normalisation
    165                     vector->data.F64[iIndex] -= norm * normTerm;
    166                 }
    167             }
    168         }
     86        psKernel *iConv = convolutions->data[i]; // Convolution for index i
     87        for (int j = i; j < numKernels; j++) {
     88            psKernel *jConv = convolutions->data[j]; // Convolution for index j
     89
     90            double sumCC = 0.0;         // Sum of convolution products
     91            for (int y = - footprint; y <= footprint; y++) {
     92                for (int x = - footprint; x <= footprint; x++) {
     93                    double cc = iConv->kernel[y][x] * jConv->kernel[y][x];
     94                    if (weight) {
     95                        cc *= weight->kernel[y][x];
     96                    }
     97                    if (window) {
     98                        cc *= window->kernel[y][x];
     99                    }
     100                    sumCC += cc;
     101                }
     102            }
     103
     104            // Spatial variation of kernel coeffs
     105            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
     106                for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
     107                    for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
     108                        double value = sumCC * poly2[iTerm][jTerm];
     109                        matrix->data.F64[iIndex][jIndex] = value;
     110                        matrix->data.F64[jIndex][iIndex] = value;
     111                    }
     112                }
     113            }
     114        }
     115
     116        double sumRC = 0.0;             // Sum of the reference-convolution products
     117        double sumIC = 0.0;             // Sum of the input-convolution products
     118        double sumC = 0.0;              // Sum of the convolution
     119        for (int y = - footprint; y <= footprint; y++) {
     120            for (int x = - footprint; x <= footprint; x++) {
     121                float conv = iConv->kernel[y][x];
     122                float in = input->kernel[y][x];
     123                float ref = reference->kernel[y][x];
     124                double ic = in * conv;
     125                double rc = ref * conv;
     126                double c = conv;
     127                if (weight) {
     128                    float wtVal = weight->kernel[y][x];
     129                    ic *= wtVal;
     130                    rc *= wtVal;
     131                    c *= wtVal;
     132                }
     133                if (window) {
     134                    float winVal = window->kernel[y][x];
     135                    ic *= winVal;
     136                    rc *= winVal;
     137                    c  *= winVal;
     138                }
     139                sumIC += ic;
     140                sumRC += rc;
     141                sumC += c;
     142            }
     143        }
     144        // Spatial variation
     145        for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
     146            double normTerm = sumRC * poly[iTerm];
     147            double bgTerm = sumC * poly[iTerm];
     148            if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
     149                matrix->data.F64[iIndex][normIndex] = normTerm;
     150                matrix->data.F64[normIndex][iIndex] = normTerm;
     151            }
     152            if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
     153                matrix->data.F64[iIndex][bgIndex] = bgTerm;
     154                matrix->data.F64[bgIndex][iIndex] = bgTerm;
     155            }
     156            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
     157                vector->data.F64[iIndex] = sumIC * poly[iTerm];
     158                // XXX TEST for Hermitians: do not calculate A - norm*B - \sum(k x B),
     159                // instead, calculate A - \sum(k x B), with full hermitians
     160                if (0 && !(mode & PM_SUBTRACTION_EQUATION_NORM)) {
     161                    // subtract norm * sumRC * poly[iTerm]
     162                    psAssert (kernels->solution1, "programming error: define solution first!");
     163                    int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     164                    double norm = fabs(kernels->solution1->data.F64[normIndex]);  // Normalisation
     165                    vector->data.F64[iIndex] -= norm * normTerm;
     166                }
     167            }
     168        }
    169169    }
    170170
    … …  
    181181            double rr = PS_SQR(ref);
    182182            double one = 1.0;
    183             if (weight) {
    184                 float wtVal = weight->kernel[y][x];
    185                 rr *= wtVal;
    186                 ir *= wtVal;
    187                 in *= wtVal;
    188                 ref *= wtVal;
    189                 one *= wtVal;
    190             }
    191             if (window) {
    192                 float  winVal = window->kernel[y][x];
    193                 rr      *= winVal;
    194                 ir      *= winVal;
    195                 in      *= winVal;
    196                 ref *= winVal;
    197                 one *= winVal;
    198             }
     183            if (weight) {
     184                float wtVal = weight->kernel[y][x];
     185                rr *= wtVal;
     186                ir *= wtVal;
     187                in *= wtVal;
     188                ref *= wtVal;
     189                one *= wtVal;
     190            }
     191            if (window) {
     192                float  winVal = window->kernel[y][x];
     193                rr      *= winVal;
     194                ir      *= winVal;
     195                in      *= winVal;
     196                ref *= winVal;
     197                one *= winVal;
     198            }
    199199            sumRR += rr;
    200200            sumIR += ir;
    … …  
    205205    }
    206206    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
    207         matrix->data.F64[normIndex][normIndex] = sumRR;
    208         vector->data.F64[normIndex] = sumIR;
    209         // subtract sum over kernels * kernel solution
     207        matrix->data.F64[normIndex][normIndex] = sumRR;
     208        vector->data.F64[normIndex] = sumIR;
     209        // subtract sum over kernels * kernel solution
    210210    }
    211211    if (mode & PM_SUBTRACTION_EQUATION_BG) {
    212         matrix->data.F64[bgIndex][bgIndex] = sum1;
    213         vector->data.F64[bgIndex] = sumI;
     212        matrix->data.F64[bgIndex][bgIndex] = sum1;
     213        vector->data.F64[bgIndex] = sumI;
    214214    }
    215215    if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) {
    216         matrix->data.F64[normIndex][bgIndex] = sumR;
    217         matrix->data.F64[bgIndex][normIndex] = sumR;
     216        matrix->data.F64[normIndex][bgIndex] = sumR;
     217        matrix->data.F64[bgIndex][normIndex] = sumR;
    218218    }
    219219    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;
    … …  
    330330                for (int x = - footprint; x <= footprint; x++) {
    331331                    double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x];
    332                     if (weight) {
    333                         ab *= weight->kernel[y][x];
    334                     }
    335                     if (window) {
    336                         ab *= window->kernel[y][x];
    337                     }
     332                    if (weight) {
     333                        ab *= weight->kernel[y][x];
     334                    }
     335                    if (window) {
     336                        ab *= window->kernel[y][x];
     337                    }
    338338                    sumAB += ab;
    339339                }
    … …  
    369369                double bi1 = b * i1;
    370370
    371                 if (weight) {
    372                     float wtVal = weight->kernel[y][x];
    373                     ai2 *= wtVal;
    374                     bi2 *= wtVal;
    375                     ai1 *= wtVal;
    376                     bi1 *= wtVal;
    377                     a *= wtVal;
    378                     b *= wtVal;
    379                     i2 *= wtVal;
    380                 }
    381                 if (window) {
    382                     float wtVal = window->kernel[y][x];
    383                     ai2 *= wtVal;
    384                     bi2 *= wtVal;
    385                     ai1 *= wtVal;
    386                     bi1 *= wtVal;
    387                     a *= wtVal;
    388                     b *= wtVal;
    389                     i2 *= wtVal;
    390                 }
     371                if (weight) {
     372                    float wtVal = weight->kernel[y][x];
     373                    ai2 *= wtVal;
     374                    bi2 *= wtVal;
     375                    ai1 *= wtVal;
     376                    bi1 *= wtVal;
     377                    a *= wtVal;
     378                    b *= wtVal;
     379                    i2 *= wtVal;
     380                }
     381                if (window) {
     382                    float wtVal = window->kernel[y][x];
     383                    ai2 *= wtVal;
     384                    bi2 *= wtVal;
     385                    ai1 *= wtVal;
     386                    bi1 *= wtVal;
     387                    a *= wtVal;
     388                    b *= wtVal;
     389                    i2 *= wtVal;
     390                }
    391391                sumAI2 += ai2;
    392392                sumBI2 += bi2;
    … …  
    434434            double i1i2 = i1 * i2;
    435435
    436             if (weight) {
    437                 float wtVal = weight->kernel[y][x];
    438                 i1 *= wtVal;
    439                 i1i1 *= wtVal;
    440                 one *= wtVal;
    441                 i2 *= wtVal;
    442                 i1i2 *= wtVal;
    443             }
    444             if (window) {
    445                 float wtVal = window->kernel[y][x];
    446                 i1 *= wtVal;
    447                 i1i1 *= wtVal;
    448                 one *= wtVal;
    449                 i2 *= wtVal;
    450                 i1i2 *= wtVal;
    451             }
     436            if (weight) {
     437                float wtVal = weight->kernel[y][x];
     438                i1 *= wtVal;
     439                i1i1 *= wtVal;
     440                one *= wtVal;
     441                i2 *= wtVal;
     442                i1i2 *= wtVal;
     443            }
     444            if (window) {
     445                float wtVal = window->kernel[y][x];
     446                i1 *= wtVal;
     447                i1i1 *= wtVal;
     448                one *= wtVal;
     449                i2 *= wtVal;
     450                i1i2 *= wtVal;
     451            }
    452452            sumI1 += i1;
    453453            sumI1I1 += i1i1;
    … …  
    486486            for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
    487487                // Contribution to chi^2: a_i^2 P_i
    488                 psAssert(isfinite(penalties->data.F32[i]), "Invalid penalty");
     488                psAssert(isfinite(penalties->data.F32[i]), "Invalid penalty");
    489489                matrix->data.F64[index][index] -= norm * penalties->data.F32[i];
    490490            }
    … …  
    494494    return true;
    495495}
    496 #endif
    497 
     496# endif
    498497
    499498//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    … …  
    596595    int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms
    597596
    598     // numKernels is the number of unique kernel images (one for each Gaussian modified by a specific polynomial).
     597    // numKernels is the number of unique kernel images (one for each Gaussian modified by a specific polynomial). 
    599598    // = \sum_i^N_Gaussians [(order + 1) * (order + 2) / 2], eg for 1 Gauss and 1st order, numKernels = 3
    600599
    … …  
    671670        status = calculateMatrixVector(stamp->matrix, stamp->vector, stamp->image2, stamp->image1,
    672671                                       weight, window, stamp->convolutions1, kernels,
    673                                        polyValues, footprint, mode);
     672                                       polyValues, footprint, mode);
    674673        break;
    675674      case PM_SUBTRACTION_MODE_2:
    676675        status = calculateMatrixVector(stamp->matrix, stamp->vector, stamp->image1, stamp->image2,
    677676                                       weight, window, stamp->convolutions2, kernels,
    678                                        polyValues, footprint, mode);
     677                                       polyValues, footprint, mode);
    679678        break;
    680679      case PM_SUBTRACTION_MODE_DUAL:
    … …  
    737736        }
    738737
    739         if ((stamp->x <= 0.0) && (stamp->y <= 0.0)) {
    740             psAbort ("bad stamp");
    741         }
    742         if (!isfinite(stamp->x) && !isfinite(stamp->y)) {
    743             psAbort ("bad stamp");
    744         }
     738        if ((stamp->x <= 0.0) && (stamp->y <= 0.0)) {
     739            psAbort ("bad stamp");
     740        }
     741        if (!isfinite(stamp->x) && !isfinite(stamp->y)) {
     742            psAbort ("bad stamp");
     743        }
    745744
    746745        if (pmSubtractionThreaded()) {
    … …  
    797796double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps);
    798797
    799 bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels,
    800                                 const pmSubtractionStampList *stamps,
    801                                 const pmSubtractionEquationCalculationMode mode)
     798bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, 
     799                                const pmSubtractionStampList *stamps,
     800                                const pmSubtractionEquationCalculationMode mode)
    802801{
    803802    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
    … …  
    914913# define SVD_TOL 0.0
    915914
    916         psVector *solution = NULL;
    917         psVector *solutionErr = NULL;
     915        psVector *solution = NULL;
     916        psVector *solutionErr = NULL;
    918917
    919918#ifdef TESTING
    920         psFitsWriteImageSimple("A.fits", sumMatrix, NULL);
    921         psVectorWriteFile("B.dat", sumVector);
    922 #endif
    923 
    924         // Use SVD to determine the kernel coeffs (and validate)
    925         if (SVD_ANALYSIS) {
    926 
    927             // We have sumVector and sumMatrix.  we are trying to solve the following equation:
    928             // sumMatrix * x = sumVector.
    929 
    930             // we can use any standard matrix inversion to solve this.  However, the basis
    931             // functions in general have substantial correlation, so that the solution may be
    932             // somewhat poorly determined or unstable.  If not numerically ill-conditioned, the
    933             // system of equations may be statistically ill-conditioned.  Noise in the image
    934             // will drive insignificant, but correlated, terms in the solution.  To avoid these
    935             // problems, we can use SVD to identify numerically unconstrained values and to
    936             // avoid statistically badly determined value.
    937 
    938             // A = sumMatrix, B = sumVector
    939             // SVD: A = U w V^T  -> A^{-1} = V (1/w) U^T
    940             // x = V (1/w) (U^T B)
    941             // \sigma_x = sqrt(diag(A^{-1}))
    942             // solve for x and A^{-1} to get x & dx
    943             // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0
    944             // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0
    945 
    946             // If I use the SVD trick to re-condition the matrix, I need to break out the
    947             // kernel and normalization terms from the background term.
    948             // XXX is this true?  or was this due to an error in the analysis?
    949 
    950             int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
    951 
    952             // now pull out the kernel elements into their own square matrix
    953             psImage  *kernelMatrix = psImageAlloc  (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64);
    954             psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64);
    955 
    956             for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) {
    957                 if (ix == bgIndex) continue;
    958                 for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) {
    959                     if (iy == bgIndex) continue;
    960                     kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix];
    961                     ky++;
    962                 }
    963                 kernelVector->data.F64[kx] = sumVector->data.F64[ix];
    964                 kx++;
    965             }
    966 
    967             psImage *U = NULL;
    968             psImage *V = NULL;
    969             psVector *w = NULL;
    970             if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) {
    971                 psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n");
    972                 return NULL;
    973             }
    974 
    975             // calculate A_inverse:
    976             // Ainv = V * w * U^T
    977             psImage *wUt  = p_pmSubSolve_wUt (w, U);
    978             psImage *Ainv = p_pmSubSolve_VwUt (V, wUt);
    979             psImage *Xvar = NULL;
    980             psFree (wUt);
     919        psFitsWriteImageSimple("A.fits", sumMatrix, NULL);
     920        psVectorWriteFile ("B.dat", sumVector);
     921#endif
     922
     923        // Use SVD to determine the kernel coeffs (and validate)
     924        if (SVD_ANALYSIS) {
     925
     926            // We have sumVector and sumMatrix.  we are trying to solve the following equation:
     927            // sumMatrix * x = sumVector.
     928
     929            // we can use any standard matrix inversion to solve this.  However, the basis
     930            // functions in general have substantial correlation, so that the solution may be
     931            // somewhat poorly determined or unstable.  If not numerically ill-conditioned, the
     932            // system of equations may be statistically ill-conditioned.  Noise in the image
     933            // will drive insignificant, but correlated, terms in the solution.  To avoid these
     934            // problems, we can use SVD to identify numerically unconstrained values and to
     935            // avoid statistically badly determined value.
     936
     937            // A = sumMatrix, B = sumVector
     938            // SVD: A = U w V^T  -> A^{-1} = V (1/w) U^T
     939            // x = V (1/w) (U^T B)
     940            // \sigma_x = sqrt(diag(A^{-1}))
     941            // solve for x and A^{-1} to get x & dx
     942            // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0
     943            // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0
     944
     945            // If I use the SVD trick to re-condition the matrix, I need to break out the
     946            // kernel and normalization terms from the background term.
     947            // XXX is this true?  or was this due to an error in the analysis?
     948
     949            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
     950
     951            // now pull out the kernel elements into their own square matrix
     952            psImage  *kernelMatrix = psImageAlloc  (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64);
     953            psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64);
     954
     955            for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) {
     956                if (ix == bgIndex) continue;
     957                for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) {
     958                    if (iy == bgIndex) continue;
     959                    kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix];
     960                    ky++;
     961                }
     962                kernelVector->data.F64[kx] = sumVector->data.F64[ix];
     963                kx++;
     964            }
     965
     966            psImage *U = NULL;
     967            psImage *V = NULL;
     968            psVector *w = NULL;
     969            if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) {
     970                psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n");
     971                return NULL;
     972            }
     973
     974            // calculate A_inverse:
     975            // Ainv = V * w * U^T
     976            psImage *wUt  = p_pmSubSolve_wUt (w, U);
     977            psImage *Ainv = p_pmSubSolve_VwUt (V, wUt);
     978            psImage *Xvar = NULL;
     979            psFree (wUt);
    981980
    982981# ifdef TESTING
    983             // kernel terms:
    984             for (int i = 0; i < w->n; i++) {
    985                 fprintf (stderr, "w: %f\n", w->data.F64[i]);
    986             }
     982            // kernel terms:
     983            for (int i = 0; i < w->n; i++) {
     984                fprintf (stderr, "w: %f\n", w->data.F64[i]);
     985            }
    987986# endif
    988             // loop over w adding in more and more of the values until chisquare is no longer
    989             // dropping significantly.
    990             // XXX this does not seem to work very well: we seem to need all terms even for
    991             // simple cases...
    992 
    993             psVector *Xsvd = NULL;
    994             {
    995                 psVector *Ax = NULL;
    996                 psVector *UtB = NULL;
    997                 psVector *wUtB = NULL;
    998 
    999                 psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64);
    1000                 psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8);
    1001                 psVectorInit (wMask, 1); // start by masking everything
    1002 
    1003                 double chiSquareLast = NAN;
    1004                 int maxWeight = 0;
    1005 
    1006                 double Axx, Bx, y2;
    1007 
    1008                 // XXX this is an attempt to exclude insignificant modes.
    1009                 // it was not successful with the ISIS kernel set: removing even
    1010                 // the least significant mode leaves additional ringing / noise
    1011                 // because the terms are so coupled.
    1012                 for (int k = 0; false && (k < w->n); k++) {
    1013 
    1014                     // unmask the k-th weight
    1015                     wMask->data.U8[k] = 0;
    1016                     p_pmSubSolve_SetWeights(wApply, w, wMask);
    1017 
    1018                     // solve for x:
    1019                     // x = V * w * (U^T * B)
    1020                     p_pmSubSolve_UtB (&UtB, U, kernelVector);
    1021                     p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
    1022                     p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
    1023 
    1024                     // chi-square for this system of equations:
    1025                     // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
    1026                     // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
    1027                     p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
    1028                     p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
    1029                     p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
    1030                     p_pmSubSolve_y2 (&y2, kernels, stamps);
    1031 
    1032                     // apparently, this works (compare with the brute force value below
    1033                     double chiSquare = Axx - 2.0*Bx + y2;
    1034                     double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare;
    1035                     chiSquareLast = chiSquare;
    1036 
    1037                     // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi);
    1038                     if (k && !maxWeight && (deltaChi < 1.0)) {
    1039                         maxWeight = k;
    1040                     }
    1041                 }
    1042 
    1043                 // keep all terms or we get extra ringing
    1044                 maxWeight = w->n;
    1045                 psVectorInit (wMask, 1);
    1046                 for (int k = 0; k < maxWeight; k++) {
    1047                     wMask->data.U8[k] = 0;
    1048                 }
    1049                 p_pmSubSolve_SetWeights(wApply, w, wMask);
    1050 
    1051                 // solve for x:
    1052                 // x = V * w * (U^T * B)
    1053                 p_pmSubSolve_UtB (&UtB, U, kernelVector);
    1054                 p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
    1055                 p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
    1056 
    1057                 // chi-square for this system of equations:
    1058                 // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
    1059                 // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
    1060                 p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
    1061                 p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
    1062                 p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
    1063                 p_pmSubSolve_y2 (&y2, kernels, stamps);
    1064 
    1065                 // apparently, this works (compare with the brute force value below
    1066                 double chiSquare = Axx - 2.0*Bx + y2;
    1067                 psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare);
    1068 
    1069                 // re-calculate A^{-1} to get new variances:
    1070                 // Ainv = V * w * U^T
    1071                 // XXX since we keep all terms, this is identical to Ainv
    1072                 psImage *wUt  = p_pmSubSolve_wUt (wApply, U);
    1073                 Xvar = p_pmSubSolve_VwUt (V, wUt);
    1074                 psFree (wUt);
    1075 
    1076                 psFree (Ax);
    1077                 psFree (UtB);
    1078                 psFree (wUtB);
    1079                 psFree (wApply);
    1080                 psFree (wMask);
    1081             }
    1082 
    1083             // copy the kernel solutions to the full solution vector:
    1084             solution = psVectorAlloc(sumVector->n, PS_TYPE_F64);
    1085             solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
    1086 
    1087             for (int ix = 0, kx = 0; ix < sumVector->n; ix++) {
    1088                 if (ix == bgIndex) {
    1089                     solution->data.F64[ix] = 0;
    1090                     solutionErr->data.F64[ix] = 0.001;
    1091                     continue;
    1092                 }
    1093                 solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]);
    1094                 solution->data.F64[ix] = Xsvd->data.F64[kx];
    1095                 kx++;
    1096             }
    1097 
    1098             psFree (kernelMatrix);
    1099             psFree (kernelVector);
    1100 
    1101             psFree (U);
    1102             psFree (V);
    1103             psFree (w);
    1104 
    1105             psFree (Ainv);
    1106             psFree (Xsvd);
    1107         } else {
    1108             psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
    1109             psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
    1110             if (!luMatrix) {
    1111                 psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n");
    1112                 psFree(solution);
    1113                 psFree(sumVector);
    1114                 psFree(sumMatrix);
    1115                 psFree(luMatrix);
    1116                 psFree(permutation);
    1117                 return NULL;
    1118             }
    1119 
    1120             solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation);
    1121             psFree(luMatrix);
    1122             psFree(permutation);
    1123             if (!solution) {
    1124                 psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
    1125                 psFree(solution);
    1126                 psFree(sumVector);
    1127                 psFree(sumMatrix);
    1128                 return NULL;
    1129             }
    1130 
    1131             // XXX LUD does not provide A^{-1}?  fake the error for now
    1132             solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
    1133             for (int ix = 0; ix < sumVector->n; ix++) {
    1134                 solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix];
    1135             }
    1136         }
    1137 
    1138         if (!kernels->solution1) {
    1139             kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);
    1140             psVectorInit (kernels->solution1, 0.0);
    1141         }
    1142 
    1143         // only update the solutions that we chose to calculate:
    1144         if (mode & PM_SUBTRACTION_EQUATION_NORM) {
    1145             int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
    1146             kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
    1147         }
    1148         if (mode & PM_SUBTRACTION_EQUATION_BG) {
    1149             int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
    1150             kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
    1151         }
    1152         if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
    1153             int numKernels = kernels->num;
    1154             int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
    1155             int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
    1156             for (int i = 0; i < numKernels * numPoly; i++) {
    1157                 // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i]));
    1158                 if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) {
    1159                     // XXX fprintf (stderr, "drop\n");
    1160                     kernels->solution1->data.F64[i] = 0.0;
    1161                 } else {
    1162                     // XXX fprintf (stderr, "keep\n");
    1163                     kernels->solution1->data.F64[i] = solution->data.F64[i];
    1164                 }
    1165             }
    1166         }
    1167         // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps);
    1168         // fprintf (stderr, "chi square Brute = %f\n", chiSquare);
    1169 
    1170         psFree(solution);
    1171         psFree(sumVector);
    1172         psFree(sumMatrix);
     987            // loop over w adding in more and more of the values until chisquare is no longer
     988            // dropping significantly.
     989            // XXX this does not seem to work very well: we seem to need all terms even for
     990            // simple cases...
     991
     992            psVector *Xsvd = NULL;
     993            {
     994                psVector *Ax = NULL;
     995                psVector *UtB = NULL;
     996                psVector *wUtB = NULL;
     997
     998                psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64);
     999                psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8);
     1000                psVectorInit (wMask, 1); // start by masking everything
     1001
     1002                double chiSquareLast = NAN;
     1003                int maxWeight = 0;
     1004
     1005                double Axx, Bx, y2;
     1006
     1007                // XXX this is an attempt to exclude insignificant modes.
     1008                // it was not successful with the ISIS kernel set: removing even
     1009                // the least significant mode leaves additional ringing / noise
     1010                // because the terms are so coupled.
     1011                for (int k = 0; false && (k < w->n); k++) {
     1012
     1013                    // unmask the k-th weight
     1014                    wMask->data.U8[k] = 0;
     1015                    p_pmSubSolve_SetWeights(wApply, w, wMask);
     1016
     1017                    // solve for x:
     1018                    // x = V * w * (U^T * B)
     1019                    p_pmSubSolve_UtB (&UtB, U, kernelVector);
     1020                    p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
     1021                    p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
     1022
     1023                    // chi-square for this system of equations:
     1024                    // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
     1025                    // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
     1026                    p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
     1027                    p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
     1028                    p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
     1029                    p_pmSubSolve_y2 (&y2, kernels, stamps);
     1030
     1031                    // apparently, this works (compare with the brute force value below
     1032                    double chiSquare = Axx - 2.0*Bx + y2;
     1033                    double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare;
     1034                    chiSquareLast = chiSquare;
     1035
     1036                    // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi);
     1037                    if (k && !maxWeight && (deltaChi < 1.0)) {
     1038                        maxWeight = k;
     1039                    }
     1040                }
     1041
     1042                // keep all terms or we get extra ringing
     1043                maxWeight = w->n;
     1044                psVectorInit (wMask, 1);
     1045                for (int k = 0; k < maxWeight; k++) {
     1046                    wMask->data.U8[k] = 0;
     1047                }
     1048                p_pmSubSolve_SetWeights(wApply, w, wMask);
     1049
     1050                // solve for x:
     1051                // x = V * w * (U^T * B)
     1052                p_pmSubSolve_UtB (&UtB, U, kernelVector);
     1053                p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
     1054                p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
     1055
     1056                // chi-square for this system of equations:
     1057                // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
     1058                // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
     1059                p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
     1060                p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
     1061                p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
     1062                p_pmSubSolve_y2 (&y2, kernels, stamps);
     1063
     1064                // apparently, this works (compare with the brute force value below
     1065                double chiSquare = Axx - 2.0*Bx + y2;
     1066                psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare);
     1067
     1068                // re-calculate A^{-1} to get new variances:
     1069                // Ainv = V * w * U^T
     1070                // XXX since we keep all terms, this is identical to Ainv
     1071                psImage *wUt  = p_pmSubSolve_wUt (wApply, U);
     1072                Xvar = p_pmSubSolve_VwUt (V, wUt);
     1073                psFree (wUt);
     1074
     1075                psFree (Ax);
     1076                psFree (UtB);
     1077                psFree (wUtB);
     1078                psFree (wApply);
     1079                psFree (wMask);
     1080            }
     1081
     1082            // copy the kernel solutions to the full solution vector:
     1083            solution = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     1084            solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     1085
     1086            for (int ix = 0, kx = 0; ix < sumVector->n; ix++) {
     1087                if (ix == bgIndex) {
     1088                    solution->data.F64[ix] = 0;
     1089                    solutionErr->data.F64[ix] = 0.001;
     1090                    continue;
     1091                }
     1092                solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]);
     1093                solution->data.F64[ix] = Xsvd->data.F64[kx];
     1094                kx++;
     1095            }
     1096
     1097            psFree (kernelMatrix);
     1098            psFree (kernelVector);
     1099
     1100            psFree (U);
     1101            psFree (V);
     1102            psFree (w);
     1103
     1104            psFree (Ainv);
     1105            psFree (Xsvd);
     1106        } else {
     1107            psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
     1108            psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
     1109            if (!luMatrix) {
     1110                psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n");
     1111                psFree(solution);
     1112                psFree(sumVector);
     1113                psFree(sumMatrix);
     1114                psFree(luMatrix);
     1115                psFree(permutation);
     1116                return NULL;
     1117            }
     1118
     1119            solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation);
     1120            psFree(luMatrix);
     1121            psFree(permutation);
     1122            if (!solution) {
     1123                psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
     1124                psFree(solution);
     1125                psFree(sumVector);
     1126                psFree(sumMatrix);
     1127                return NULL;
     1128            }
     1129
     1130            // XXX LUD does not provide A^{-1}?  fake the error for now
     1131            solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     1132            for (int ix = 0; ix < sumVector->n; ix++) {
     1133                solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix];
     1134            }
     1135        }
     1136
     1137        if (!kernels->solution1) {
     1138            kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);
     1139            psVectorInit (kernels->solution1, 0.0);
     1140        }
     1141
     1142        // only update the solutions that we chose to calculate:
     1143        if (mode & PM_SUBTRACTION_EQUATION_NORM) {
     1144            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     1145            kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
     1146        }
     1147        if (mode & PM_SUBTRACTION_EQUATION_BG) {
     1148            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
     1149            kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
     1150        }
     1151        if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
     1152            int numKernels = kernels->num;
     1153            int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
     1154            int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
     1155            for (int i = 0; i < numKernels * numPoly; i++) {
     1156                // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i]));
     1157                if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) {
     1158                    // XXX fprintf (stderr, "drop\n");
     1159                    kernels->solution1->data.F64[i] = 0.0;
     1160                } else {
     1161                    // XXX fprintf (stderr, "keep\n");
     1162                    kernels->solution1->data.F64[i] = solution->data.F64[i];
     1163                }
     1164            }
     1165        }
     1166        // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps);
     1167        // fprintf (stderr, "chi square Brute = %f\n", chiSquare);
     1168
     1169        psFree(solution);
     1170        psFree(sumVector);
     1171        psFree(sumMatrix);
    11731172
    11741173#ifdef TESTING
    … …  
    12011200            }
    12021201        }
    1203 
    12041202
    12051203#ifdef TESTING
    … …  
    13191317            psVectorWriteFile("sumVectorFix.dat", sumVector);
    13201318#endif
    1321 
    1322         }
    1323 #endif
    1324 
     1319        }
     1320#endif
     1321    }
     1322#endif
    13251323        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
    13261324
    … …  
    13691367
    13701368psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps,
    1371                                            const pmSubtractionKernels *kernels)
     1369                                           pmSubtractionKernels *kernels)
    13721370{
    13731371    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL);
    … …  
    13871385    pmSubtractionVisualShowFitInit (stamps);
    13881386
     1387    psVector *fSigRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
     1388    psVector *fMinRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
     1389    psVector *fMaxRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
     1390
    13891391    for (int i = 0; i < stamps->num; i++) {
    1390         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
    1391         if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
    1392             deviations->data.F32[i] = NAN;
    1393             continue;
    1394         }
    1395 
    1396         // Calculate coefficients of the kernel basis functions
    1397         polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
    1398         double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
    1399         double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
    1400 
    1401         // Calculate residuals
    1402         psKernel *weight = stamp->weight; // Weight postage stamp
    1403         psImageInit(residual->image, 0.0);
    1404         if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
    1405             psKernel *target;           // Target postage stamp
    1406             psKernel *source;           // Source postage stamp
    1407             psArray *convolutions;      // Convolution postage stamps for each kernel basis function
    1408             switch (kernels->mode) {
    1409               case PM_SUBTRACTION_MODE_1:
    1410                 target = stamp->image2;
    1411                 source = stamp->image1;
    1412                 convolutions = stamp->convolutions1;
    1413 
    1414                 // Having convolved image1 and changed its normalisation, we need to renormalise the residual
    1415                 // so that it is on the scale of image1.
    1416                 psImage *image = pmSubtractionKernelImage(kernels, stamp->xNorm, stamp->yNorm,
    1417                                                           false); // Kernel image
    1418                 if (!image) {
    1419                     psError(PS_ERR_UNKNOWN, false, "Unable to generate image of kernel.");
    1420                     return false;
    1421                 }
    1422                 double sumKernel = 0;   // Sum of kernel, for normalising residual
    1423                 int size = kernels->size; // Half-size of kernel
    1424                 int fullSize = 2 * size + 1; // Full size of kernel
    1425                 for (int y = 0; y < fullSize; y++) {
    1426                     for (int x = 0; x < fullSize; x++) {
    1427                         sumKernel += image->data.F32[y][x];
    1428                     }
    1429                 }
    1430                 psFree(image);
    1431                 devNorm = 1.0 / sumKernel / numPixels;
    1432                 break;
    1433               case PM_SUBTRACTION_MODE_2:
    1434                 target = stamp->image1;
    1435                 source = stamp->image2;
    1436                 convolutions = stamp->convolutions2;
    1437                 break;
    1438               default:
    1439                 psAbort("Unsupported subtraction mode: %x", kernels->mode);
    1440             }
    1441 
    1442             for (int j = 0; j < numKernels; j++) {
    1443                 psKernel *convolution = convolutions->data[j]; // Convolution
    1444                 double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j,
    1445                                                                   false); // Coefficient
    1446                 for (int y = - footprint; y <= footprint; y++) {
    1447                     for (int x = - footprint; x <= footprint; x++) {
    1448                         residual->kernel[y][x] += convolution->kernel[y][x] * coefficient;
    1449                     }
    1450                 }
    1451             }
    1452 
    1453             // XXX visualize the target, source, convolution and residual
    1454             pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i);
    1455 
    1456             for (int y = - footprint; y <= footprint; y++) {
    1457                 for (int x = - footprint; x <= footprint; x++) {
    1458                     residual->kernel[y][x] += background + source->kernel[y][x] * norm - target->kernel[y][x];
    1459                 }
    1460             }
    1461         } else {
    1462             // Dual convolution
    1463             psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image
    1464             psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image
    1465             psKernel *image1 = stamp->image1; // The first image
    1466             psKernel *image2 = stamp->image2; // The second image
    1467 
    1468             for (int j = 0; j < numKernels; j++) {
    1469                 psKernel *conv1 = convolutions1->data[j]; // Convolution of first image
    1470                 psKernel *conv2 = convolutions2->data[j]; // Convolution of second image
    1471                 double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1
    1472                 double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2
    1473 
    1474                 for (int y = - footprint; y <= footprint; y++) {
    1475                     for (int x = - footprint; x <= footprint; x++) {
    1476                         residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 + conv1->kernel[y][x] * coeff1;
    1477                     }
    1478                 }
    1479             }
    1480 
    1481             // XXX visualize the target, source, convolution and residual
    1482             pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i);
    1483 
    1484             for (int y = - footprint; y <= footprint; y++) {
    1485                 for (int x = - footprint; x <= footprint; x++) {
    1486                     residual->kernel[y][x] += background + image1->kernel[y][x] * norm - image2->kernel[y][x];
    1487                 }
    1488             }
    1489         }
    1490 
    1491         double deviation = 0.0;         // Sum of differences
    1492         for (int y = - footprint; y <= footprint; y++) {
    1493             for (int x = - footprint; x <= footprint; x++) {
    1494                 double dev = PS_SQR(residual->kernel[y][x]) * weight->kernel[y][x];
    1495                 deviation += dev;
     1392        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
     1393        if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
     1394            deviations->data.F32[i] = NAN;
     1395            continue;
     1396        }
     1397
     1398        // Calculate coefficients of the kernel basis functions
     1399        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
     1400        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
     1401        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
     1402
     1403        // Calculate residuals
     1404        psKernel *weight = stamp->weight; // Weight postage stamp
     1405        psImageInit(residual->image, 0.0);
     1406        if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
     1407            psKernel *target;           // Target postage stamp
     1408            psKernel *source;           // Source postage stamp
     1409            psArray *convolutions;      // Convolution postage stamps for each kernel basis function
     1410            switch (kernels->mode) {
     1411              case PM_SUBTRACTION_MODE_1:
     1412                target = stamp->image2;
     1413                source = stamp->image1;
     1414                convolutions = stamp->convolutions1;
     1415
     1416                // Having convolved image1 and changed its normalisation, we need to renormalise the residual
     1417                // so that it is on the scale of image1.
     1418                psImage *image = pmSubtractionKernelImage(kernels, stamp->xNorm, stamp->yNorm,
     1419                                                          false); // Kernel image
     1420                if (!image) {
     1421                    psError(PS_ERR_UNKNOWN, false, "Unable to generate image of kernel.");
     1422                    return false;
     1423                }
     1424                double sumKernel = 0;   // Sum of kernel, for normalising residual
     1425                int size = kernels->size; // Half-size of kernel
     1426                int fullSize = 2 * size + 1; // Full size of kernel
     1427                for (int y = 0; y < fullSize; y++) {
     1428                    for (int x = 0; x < fullSize; x++) {
     1429                        sumKernel += image->data.F32[y][x];
     1430                    }
     1431                }
     1432                psFree(image);
     1433                devNorm = 1.0 / sumKernel / numPixels;
     1434                break;
     1435              case PM_SUBTRACTION_MODE_2:
     1436                target = stamp->image1;
     1437                source = stamp->image2;
     1438                convolutions = stamp->convolutions2;
     1439                break;
     1440              default:
     1441                psAbort("Unsupported subtraction mode: %x", kernels->mode);
     1442            }
     1443
     1444            for (int j = 0; j < numKernels; j++) {
     1445                psKernel *convolution = convolutions->data[j]; // Convolution
     1446                double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j,
     1447                                                                  false); // Coefficient
     1448                for (int y = - footprint; y <= footprint; y++) {
     1449                    for (int x = - footprint; x <= footprint; x++) {
     1450                        residual->kernel[y][x] += convolution->kernel[y][x] * coefficient;
     1451                    }
     1452                }
     1453            }
     1454
     1455            // XXX visualize the target, source, convolution and residual
     1456            pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i);
     1457
     1458            for (int y = - footprint; y <= footprint; y++) {
     1459                for (int x = - footprint; x <= footprint; x++) {
     1460                    residual->kernel[y][x] += background + source->kernel[y][x] * norm - target->kernel[y][x];
     1461                }
     1462            }
     1463
     1464            // XXX measure some useful stats on the residuals
     1465            float sum = 0.0;
     1466            float peak = 0.0;
     1467            for (int y = - footprint; y <= footprint; y++) {
     1468                for (int x = - footprint; x <= footprint; x++) {
     1469                    sum += 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);
     1470                    peak = PS_MAX(peak, 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm));
     1471                }
     1472            }
     1473
     1474            // only count pixels with more than X% of the source flux
     1475            // calculate stdev(dflux)
     1476            float dflux1 = 0.0;
     1477            float dflux2 = 0.0;
     1478            int npix = 0;
     1479
     1480            float dmax = 0.0;
     1481            float dmin = 0.0;
     1482
     1483            for (int y = - footprint; y <= footprint; y++) {
     1484                for (int x = - footprint; x <= footprint; x++) {
     1485                    float dflux = 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);
     1486                    if (dflux < 0.02*sum) continue;
     1487                    dflux1 += residual->kernel[y][x];
     1488                    dflux2 += PS_SQR(residual->kernel[y][x]);
     1489                    dmax = PS_MAX(residual->kernel[y][x], dmax);
     1490                    dmin = PS_MIN(residual->kernel[y][x], dmin);
     1491                    npix ++;
     1492                }
     1493            }
     1494            float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix));
     1495            // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dmax/peak, dmin/peak);
     1496            psVectorAppend(fSigRes, sigma/sum);
     1497            psVectorAppend(fMaxRes, dmax/peak);
     1498            psVectorAppend(fMinRes, dmin/peak);
     1499        } else {
     1500            // Dual convolution
     1501            psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image
     1502            psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image
     1503            psKernel *image1 = stamp->image1; // The first image
     1504            psKernel *image2 = stamp->image2; // The second image
     1505
     1506            for (int j = 0; j < numKernels; j++) {
     1507                psKernel *conv1 = convolutions1->data[j]; // Convolution of first image
     1508                psKernel *conv2 = convolutions2->data[j]; // Convolution of second image
     1509                double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1
     1510                double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2
     1511
     1512                for (int y = - footprint; y <= footprint; y++) {
     1513                    for (int x = - footprint; x <= footprint; x++) {
     1514                        residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 + conv1->kernel[y][x] * coeff1;
     1515                    }
     1516                }
     1517            }
     1518
     1519            // XXX visualize the target, source, convolution and residual
     1520            pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i);
     1521
     1522            for (int y = - footprint; y <= footprint; y++) {
     1523                for (int x = - footprint; x <= footprint; x++) {
     1524                    residual->kernel[y][x] += background + image1->kernel[y][x] * norm - image2->kernel[y][x];
     1525                }
     1526            }
     1527        }
     1528
     1529        double deviation = 0.0;         // Sum of differences
     1530        for (int y = - footprint; y <= footprint; y++) {
     1531            for (int x = - footprint; x <= footprint; x++) {
     1532                double dev = PS_SQR(residual->kernel[y][x]) * weight->kernel[y][x];
     1533                deviation += dev;
    14961534#ifdef TESTING
    1497                 residual->kernel[y][x] = dev;
    1498 #endif
    1499             }
    1500         }
    1501         deviations->data.F32[i] = devNorm * deviation;
    1502         psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n",
    1503                 i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]);
    1504         if (!isfinite(deviations->data.F32[i])) {
    1505             stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
    1506             psTrace("psModules.imcombine", 5,
    1507                     "Rejecting stamp %d (%d,%d) because of non-finite deviation\n",
    1508                     i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
    1509             continue;
    1510         }
     1535                residual->kernel[y][x] = dev;
     1536#endif
     1537            }
     1538        }
     1539        deviations->data.F32[i] = devNorm * deviation;
     1540        psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n",
     1541                i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]);
     1542        if (!isfinite(deviations->data.F32[i])) {
     1543            stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
     1544            psTrace("psModules.imcombine", 5,
     1545                    "Rejecting stamp %d (%d,%d) because of non-finite deviation\n",
     1546                    i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
     1547            continue;
     1548        }
    15111549
    15121550#ifdef TESTING
    1513         {
    1514             psString filename = NULL;
    1515             psStringAppend(&filename, "resid_%03d.fits", i);
    1516             psFits *fits = psFitsOpen(filename, "w");
    1517             psFree(filename);
    1518             psFitsWriteImage(fits, NULL, residual->image, 0, NULL);
    1519             psFitsClose(fits);
    1520         }
    1521         if (stamp->image1) {
    1522             psString filename = NULL;
    1523             psStringAppend(&filename, "stamp_image1_%03d.fits", i);
    1524             psFits *fits = psFitsOpen(filename, "w");
    1525             psFree(filename);
    1526             psFitsWriteImage(fits, NULL, stamp->image1->image, 0, NULL);
    1527             psFitsClose(fits);
    1528         }
    1529         if (stamp->image2) {
    1530             psString filename = NULL;
    1531             psStringAppend(&filename, "stamp_image2_%03d.fits", i);
    1532             psFits *fits = psFitsOpen(filename, "w");
    1533             psFree(filename);
    1534             psFitsWriteImage(fits, NULL, stamp->image2->image, 0, NULL);
    1535             psFitsClose(fits);
    1536         }
    1537         if (stamp->weight) {
    1538             psString filename = NULL;
    1539             psStringAppend(&filename, "stamp_weight_%03d.fits", i);
    1540             psFits *fits = psFitsOpen(filename, "w");
    1541             psFree(filename);
    1542             psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL);
    1543             psFitsClose(fits);
    1544         }
    1545 #endif
    1546 
     1551        {
     1552            psString filename = NULL;
     1553            psStringAppend(&filename, "resid_%03d.fits", i);
     1554            psFits *fits = psFitsOpen(filename, "w");
     1555            psFree(filename);
     1556            psFitsWriteImage(fits, NULL, residual->image, 0, NULL);
     1557            psFitsClose(fits);
     1558        }
     1559        if (stamp->image1) {
     1560            psString filename = NULL;
     1561            psStringAppend(&filename, "stamp_image1_%03d.fits", i);
     1562            psFits *fits = psFitsOpen(filename, "w");
     1563            psFree(filename);
     1564            psFitsWriteImage(fits, NULL, stamp->image1->image, 0, NULL);
     1565            psFitsClose(fits);
     1566        }
     1567        if (stamp->image2) {
     1568            psString filename = NULL;
     1569            psStringAppend(&filename, "stamp_image2_%03d.fits", i);
     1570            psFits *fits = psFitsOpen(filename, "w");
     1571            psFree(filename);
     1572            psFitsWriteImage(fits, NULL, stamp->image2->image, 0, NULL);
     1573            psFitsClose(fits);
     1574        }
     1575        if (stamp->weight) {
     1576            psString filename = NULL;
     1577            psStringAppend(&filename, "stamp_weight_%03d.fits", i);
     1578            psFits *fits = psFitsOpen(filename, "w");
     1579            psFree(filename);
     1580            psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL);
     1581            psFitsClose(fits);
     1582        }
     1583#endif
     1584   
    15471585    }
    15481586
    15491587    // calculate and report the normalization and background for the image center
    1550     {
    1551         polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, 0.0, 0.0);
    1552         double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
    1553         double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
    1554         psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f, background: %f", norm, background);
    1555 
    1556         pmSubtractionVisualShowFit(norm);
    1557         pmSubtractionVisualPlotFit(kernels);
     1588    {
     1589        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, 0.0, 0.0);
     1590        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
     1591        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
     1592        psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f, background: %f", norm, background);
     1593
     1594        pmSubtractionVisualShowFit(norm);
     1595        pmSubtractionVisualPlotFit(kernels);
     1596
     1597        psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     1598        psVectorStats (stats, fSigRes, NULL, NULL, 0);
     1599        kernels->fSigResMean = stats->robustMedian;
     1600        kernels->fSigResStdev = stats->robustStdev;
     1601
     1602        psStatsInit (stats);
     1603        psVectorStats (stats, fMaxRes, NULL, NULL, 0);
     1604        kernels->fMaxResMean = stats->robustMedian;
     1605        kernels->fMaxResStdev = stats->robustStdev;
     1606
     1607        psStatsInit (stats);
     1608        psVectorStats (stats, fMinRes, NULL, NULL, 0);
     1609        kernels->fMinResMean = stats->robustMedian;
     1610        kernels->fMinResStdev = stats->robustStdev;
     1611
     1612        // XXX save these values somewhere
     1613        psLogMsg("psModules.imcombine", PS_LOG_INFO, "fSigma: %f +/- %f, fMaxRes: %f +/- %f, fMinRes: %f +/- %f",
     1614                 kernels->fSigResMean, kernels->fSigResStdev,
     1615                 kernels->fMaxResMean, kernels->fMaxResStdev,
     1616                 kernels->fMinResMean, kernels->fMinResStdev);
     1617
     1618        psFree (fSigRes);
     1619        psFree (fMaxRes);
     1620        psFree (fMinRes);
     1621        psFree (stats);
    15581622    }
    15591623
    15601624    psFree(residual);
    15611625    psFree(polyValues);
     1626
    15621627
    15631628    return deviations;
    … …  
    15741639
    15751640    for (int i = 0; i < wUt->numCols; i++) {
    1576         for (int j = 0; j < wUt->numRows; j++) {
    1577             if (!isfinite(w->data.F64[j])) continue;
    1578             if (w->data.F64[j] == 0.0) continue;
    1579             wUt->data.F64[j][i] = U->data.F64[i][j] / w->data.F64[j];
    1580         }
     1641        for (int j = 0; j < wUt->numRows; j++) {
     1642            if (!isfinite(w->data.F64[j])) continue;
     1643            if (w->data.F64[j] == 0.0) continue;
     1644            wUt->data.F64[j][i] = U->data.F64[i][j] / w->data.F64[j];
     1645        }
    15811646    }
    15821647    return wUt;
    … …  
    15911656
    15921657    for (int i = 0; i < Ainv->numCols; i++) {
    1593         for (int j = 0; j < Ainv->numRows; j++) {
    1594             double sum = 0.0;
    1595             for (int k = 0; k < V->numCols; k++) {
    1596                 sum += V->data.F64[j][k] * wUt->data.F64[k][i];
    1597             }
    1598             Ainv->data.F64[j][i] = sum;
    1599         }
     1658        for (int j = 0; j < Ainv->numRows; j++) {
     1659            double sum = 0.0;
     1660            for (int k = 0; k < V->numCols; k++) {
     1661                sum += V->data.F64[j][k] * wUt->data.F64[k][i];
     1662            }
     1663            Ainv->data.F64[j][i] = sum;
     1664        }
    16001665    }
    16011666    return Ainv;
    … …  
    16101675
    16111676    for (int i = 0; i < U->numCols; i++) {
    1612         double sum = 0.0;
    1613         for (int j = 0; j < U->numRows; j++) {
    1614             sum += B->data.F64[j] * U->data.F64[j][i];
    1615         }
    1616         UtB[0]->data.F64[i] = sum;
     1677        double sum = 0.0;
     1678        for (int j = 0; j < U->numRows; j++) {
     1679            sum += B->data.F64[j] * U->data.F64[j][i];
     1680        }
     1681        UtB[0]->data.F64[i] = sum;
    16171682    }
    16181683    return true;
    … …  
    16291694
    16301695    for (int i = 0; i < w->n; i++) {
    1631         if (!isfinite(w->data.F64[i])) continue;
    1632         if (w->data.F64[i] == 0.0) continue;
    1633         wUtB[0]->data.F64[i] = UtB->data.F64[i] / w->data.F64[i];
     1696        if (!isfinite(w->data.F64[i])) continue;
     1697        if (w->data.F64[i] == 0.0) continue;
     1698        wUtB[0]->data.F64[i] = UtB->data.F64[i] / w->data.F64[i];
    16341699    }
    16351700    return true;
    … …  
    16441709
    16451710    for (int j = 0; j < V->numRows; j++) {
    1646         double sum = 0.0;
    1647         for (int i = 0; i < V->numCols; i++) {
    1648             sum += V->data.F64[j][i] * wUtB->data.F64[i];
    1649         }
    1650         VwUtB[0]->data.F64[j] = sum;
     1711        double sum = 0.0;
     1712        for (int i = 0; i < V->numCols; i++) {
     1713            sum += V->data.F64[j][i] * wUtB->data.F64[i];
     1714        }
     1715        VwUtB[0]->data.F64[j] = sum;
    16511716    }
    16521717    return true;
    … …  
    16611726
    16621727    for (int j = 0; j < A->numRows; j++) {
    1663         double sum = 0.0;
    1664         for (int i = 0; i < A->numCols; i++) {
    1665             sum += A->data.F64[j][i] * x->data.F64[i];
    1666         }
    1667         B[0]->data.F64[j] = sum;
     1728        double sum = 0.0;
     1729        for (int i = 0; i < A->numCols; i++) {
     1730            sum += A->data.F64[j][i] * x->data.F64[i];
     1731        }
     1732        B[0]->data.F64[j] = sum;
    16681733    }
    16691734    return true;
    … …  
    16771742    double sum = 0.0;
    16781743    for (int i = 0; i < x->n; i++) {
    1679         sum += x->data.F64[i] * y->data.F64[i];
     1744        sum += x->data.F64[i] * y->data.F64[i];
    16801745    }
    16811746    *value = sum;
    … …  
    16901755    for (int i = 0; i < stamps->num; i++) {
    16911756
    1692         pmSubtractionStamp *stamp = stamps->stamps->data[i];
    1693         if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
    1694 
    1695         psKernel *weight = NULL;
    1696         psKernel *window = NULL;
    1697         psKernel *input = NULL;
     1757        pmSubtractionStamp *stamp = stamps->stamps->data[i];
     1758        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
     1759
     1760        psKernel *weight = NULL;
     1761        psKernel *window = NULL;
     1762        psKernel *input = NULL;
    16981763
    16991764#ifdef USE_WEIGHT
    1700         weight = stamp->weight;
     1765        weight = stamp->weight;
    17011766#endif
    17021767#ifdef USE_WINDOW
    1703         window = stamps->window;
    1704 #endif
    1705 
    1706         switch (kernels->mode) {
    1707             // MODE_1 : convolve image 1 to match image 2 (and vice versa)
    1708           case PM_SUBTRACTION_MODE_1:
    1709             input = stamp->image2;
    1710             break;
    1711           case PM_SUBTRACTION_MODE_2:
    1712             input = stamp->image1;
    1713             break;
    1714           default:
    1715             psAbort ("programming error");
    1716         }
    1717 
    1718         for (int y = - footprint; y <= footprint; y++) {
    1719             for (int x = - footprint; x <= footprint; x++) {
    1720                 double in = input->kernel[y][x];
    1721                 double value = in*in;
    1722                 if (weight) {
    1723                     float wtVal = weight->kernel[y][x];
    1724                     value *= wtVal;
    1725                 }
    1726                 if (window) {
    1727                     float  winVal = window->kernel[y][x];
    1728                     value *= winVal;
    1729                 }
    1730                 sum += value;
    1731             }
    1732         }
     1768        window = stamps->window;
     1769#endif
     1770
     1771        switch (kernels->mode) {
     1772            // MODE_1 : convolve image 1 to match image 2 (and vice versa)
     1773          case PM_SUBTRACTION_MODE_1:
     1774            input = stamp->image2;
     1775            break;
     1776          case PM_SUBTRACTION_MODE_2:
     1777            input = stamp->image1;
     1778            break;
     1779          default:
     1780            psAbort ("programming error");
     1781        }
     1782
     1783        for (int y = - footprint; y <= footprint; y++) {
     1784            for (int x = - footprint; x <= footprint; x++) {
     1785                double in = input->kernel[y][x];
     1786                double value = in*in;
     1787                if (weight) {
     1788                    float wtVal = weight->kernel[y][x];
     1789                    value *= wtVal;
     1790                }
     1791                if (window) {
     1792                    float  winVal = window->kernel[y][x];
     1793                    value *= winVal;
     1794                }
     1795                sum += value;
     1796            }
     1797        }
    17331798    }
    17341799    *y2 = sum;
    … …  
    17501815    for (int i = 0; i < stamps->num; i++) {
    17511816
    1752         pmSubtractionStamp *stamp = stamps->stamps->data[i];
    1753         if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
    1754 
    1755         psKernel *weight = NULL;
    1756         psKernel *window = NULL;
    1757         psKernel *target = NULL;
    1758         psKernel *source = NULL;
    1759 
    1760         psArray *convolutions = NULL;
     1817        pmSubtractionStamp *stamp = stamps->stamps->data[i];
     1818        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
     1819
     1820        psKernel *weight = NULL;
     1821        psKernel *window = NULL;
     1822        psKernel *target = NULL;
     1823        psKernel *source = NULL;
     1824
     1825        psArray *convolutions = NULL;
    17611826
    17621827#ifdef USE_WEIGHT
    1763         weight = stamp->weight;
     1828        weight = stamp->weight;
    17641829#endif
    17651830#ifdef USE_WINDOW
    1766         window = stamps->window;
    1767 #endif
    1768 
    1769         switch (kernels->mode) {
    1770             // MODE_1 : convolve image 1 to match image 2 (and vice versa)
    1771           case PM_SUBTRACTION_MODE_1:
    1772             target = stamp->image2;
    1773             source = stamp->image1;
    1774             convolutions = stamp->convolutions1;
    1775             break;
    1776           case PM_SUBTRACTION_MODE_2:
    1777             target = stamp->image1;
    1778             source = stamp->image2;
    1779             convolutions = stamp->convolutions2;
    1780             break;
    1781           default:
    1782             psAbort ("programming error");
    1783         }
    1784 
    1785         // Calculate coefficients of the kernel basis functions
    1786         polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
    1787         double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
    1788         double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
    1789 
    1790         psImageInit(residual->image, 0.0);
    1791         for (int j = 0; j < numKernels; j++) {
    1792             psKernel *convolution = convolutions->data[j]; // Convolution
    1793             double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient
    1794             for (int y = - footprint; y <= footprint; y++) {
    1795                 for (int x = - footprint; x <= footprint; x++) {
    1796                     residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient;
    1797                 }
    1798             }
    1799         }
    1800 
    1801         for (int y = - footprint; y <= footprint; y++) {
    1802             for (int x = - footprint; x <= footprint; x++) {
    1803                 double resid = target->kernel[y][x] - background - source->kernel[y][x] * norm + residual->kernel[y][x];
    1804                 double value = PS_SQR(resid);
    1805                 if (weight) {
    1806                     float wtVal = weight->kernel[y][x];
    1807                     value *= wtVal;
    1808                 }
    1809                 if (window) {
    1810                     float  winVal = window->kernel[y][x];
    1811                     value *= winVal;
    1812                 }
    1813                 sum += value;
    1814             }
    1815         }
     1831        window = stamps->window;
     1832#endif
     1833
     1834        switch (kernels->mode) {
     1835            // MODE_1 : convolve image 1 to match image 2 (and vice versa)
     1836          case PM_SUBTRACTION_MODE_1:
     1837            target = stamp->image2;
     1838            source = stamp->image1;
     1839            convolutions = stamp->convolutions1;
     1840            break;
     1841          case PM_SUBTRACTION_MODE_2:
     1842            target = stamp->image1;
     1843            source = stamp->image2;
     1844            convolutions = stamp->convolutions2;
     1845            break;
     1846          default:
     1847            psAbort ("programming error");
     1848        }
     1849
     1850        // Calculate coefficients of the kernel basis functions
     1851        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
     1852        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
     1853        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
     1854
     1855        psImageInit(residual->image, 0.0);
     1856        for (int j = 0; j < numKernels; j++) {
     1857            psKernel *convolution = convolutions->data[j]; // Convolution
     1858            double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient
     1859            for (int y = - footprint; y <= footprint; y++) {
     1860                for (int x = - footprint; x <= footprint; x++) {
     1861                    residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient;
     1862                }
     1863            }
     1864        }
     1865
     1866        for (int y = - footprint; y <= footprint; y++) {
     1867            for (int x = - footprint; x <= footprint; x++) {
     1868                double resid = target->kernel[y][x] - background - source->kernel[y][x] * norm + residual->kernel[y][x];
     1869                double value = PS_SQR(resid);
     1870                if (weight) {
     1871                    float wtVal = weight->kernel[y][x];
     1872                    value *= wtVal;
     1873                }
     1874                if (window) {
     1875                    float  winVal = window->kernel[y][x];
     1876                    value *= winVal;
     1877                }
     1878                sum += value;
     1879            }
     1880        }
    18161881    }
    18171882    psFree (polyValues);
    … …  
    18241889
    18251890    for (int i = 0; i < w->n; i++) {
    1826         wApply->data.F64[i] = wMask->data.U8[i] ? 0.0 : w->data.F64[i];
     1891        wApply->data.F64[i] = wMask->data.U8[i] ? 0.0 : w->data.F64[i];
    18271892    }
    18281893    return true;
    … …  
    18391904    // generate Vn = V * w^{-1}
    18401905    for (int j = 0; j < Vn->numRows; j++) {
    1841         for (int i = 0; i < Vn->numCols; i++) {
    1842             if (!isfinite(w->data.F64[i])) continue;
    1843             if (w->data.F64[i] == 0.0) continue;
    1844             Vn->data.F64[j][i] = V->data.F64[j][i] / w->data.F64[i];
    1845         }
     1906        for (int i = 0; i < Vn->numCols; i++) {
     1907            if (!isfinite(w->data.F64[i])) continue;
     1908            if (w->data.F64[i] == 0.0) continue;
     1909            Vn->data.F64[j][i] = V->data.F64[j][i] / w->data.F64[i];
     1910        }
    18461911    }
    18471912
    … …  
    18511916    // generate Xvar = Vn * Vn^T
    18521917    for (int j = 0; j < Vn->numRows; j++) {
    1853         for (int i = 0; i < Vn->numCols; i++) {
    1854             double sum = 0.0;
    1855             for (int k = 0; k < Vn->numCols; k++) {
    1856                 sum += Vn->data.F64[k][i]*Vn->data.F64[k][j];
    1857             }
    1858             Xvar->data.F64[j][i] = sum;
    1859         }
     1918        for (int i = 0; i < Vn->numCols; i++) {
     1919            double sum = 0.0;
     1920            for (int k = 0; k < Vn->numCols; k++) {
     1921                sum += Vn->data.F64[k][i]*Vn->data.F64[k][j];
     1922            }
     1923            Xvar->data.F64[j][i] = sum;
     1924        }
    18601925    }
    18611926    return Xvar;
    … …  
    18721937    psFitsWriteImage(fits, header, image, 0, NULL);
    18731938    psFitsClose(fits);
    1874 
     1939   
    18751940    return true;
    18761941}
Note: See TracChangeset for help on using the changeset viewer.