IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26586


Ignore:
Timestamp:
Jan 13, 2010, 2:37:39 PM (16 years ago)
Author:
Paul Price
Message:

Removing unnecessary parts.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r26583 r26586  
    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            }
    … …  
    595595    int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms
    596596
    597     // 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).
    598598    // = \sum_i^N_Gaussians [(order + 1) * (order + 2) / 2], eg for 1 Gauss and 1st order, numKernels = 3
    599599
    … …  
    670670        status = calculateMatrixVector(stamp->matrix, stamp->vector, stamp->image2, stamp->image1,
    671671                                       weight, window, stamp->convolutions1, kernels,
    672                                        polyValues, footprint, mode);
     672                                       polyValues, footprint, mode);
    673673        break;
    674674      case PM_SUBTRACTION_MODE_2:
    675675        status = calculateMatrixVector(stamp->matrix, stamp->vector, stamp->image1, stamp->image2,
    676676                                       weight, window, stamp->convolutions2, kernels,
    677                                        polyValues, footprint, mode);
     677                                       polyValues, footprint, mode);
    678678        break;
    679679      case PM_SUBTRACTION_MODE_DUAL:
    … …  
    736736        }
    737737
    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         }
     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        }
    744744
    745745        if (pmSubtractionThreaded()) {
    … …  
    796796double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps);
    797797
    798 bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, 
    799                                 const pmSubtractionStampList *stamps,
    800                                 const pmSubtractionEquationCalculationMode mode)
     798bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels,
     799                                const pmSubtractionStampList *stamps,
     800                                const pmSubtractionEquationCalculationMode mode)
    801801{
    802802    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
    … …  
    913913# define SVD_TOL 0.0
    914914
    915         psVector *solution = NULL;
    916         psVector *solutionErr = NULL;
     915        psVector *solution = NULL;
     916        psVector *solutionErr = NULL;
    917917
    918918#ifdef TESTING
    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);
     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);
    980980
    981981# ifdef TESTING
    982             // kernel terms:
    983             for (int i = 0; i < w->n; i++) {
    984                 fprintf (stderr, "w: %f\n", w->data.F64[i]);
    985             }
     982            // kernel terms:
     983            for (int i = 0; i < w->n; i++) {
     984                fprintf (stderr, "w: %f\n", w->data.F64[i]);
     985            }
    986986# endif
    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);
     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);
    11721172
    11731173#ifdef TESTING
    … …  
    12421242                    sumMatrix->data.F64[index][index] = 1.0;            \
    12431243                    sumVector->data.F64[index] = 0.0;                   \
    1244                 }                                                       \
    1245             }
    1246 
    1247             // Remove a kernel basis for image 1 from the equation
    1248 #define MASK_BASIS_1(INDEX)                                             \
    1249             {                                                           \
    1250                 for (int j = 0, index = INDEX; j < numSpatial; j++, index += numKernels) { \
    1251                     for (int k = 0; k < numParams2; k++) {              \
    1252                         sumMatrix1->data.F64[k][index] = 0.0;           \
    1253                         sumMatrix1->data.F64[index][k] = 0.0;           \
    1254                         sumMatrixX->data.F64[k][index] = 0.0;           \
    1255                     }                                                   \
    1256                     sumMatrix1->data.F64[bgIndex][index] = 0.0;         \
    1257                     sumMatrix1->data.F64[index][bgIndex] = 0.0;         \
    1258                     sumMatrix1->data.F64[normIndex][index] = 0.0;       \
    1259                     sumMatrix1->data.F64[index][normIndex] = 0.0;       \
    1260                     sumMatrix1->data.F64[index][index] = 1.0;           \
    1261                     sumVector1->data.F64[index] = 0.0;                  \
    1262                 }                                                       \
    1263             }
    1264 
    1265             // Remove a kernel basis for image 2 from the equation
    1266 #define MASK_BASIS_2(INDEX)                                             \
    1267             {                                                           \
    1268                 for (int j = 0, index = INDEX; j < numSpatial; j++, index += numKernels) { \
    1269                     for (int k = 0; k < numParams2; k++) {              \
    1270                         sumMatrix2->data.F64[k][index] = 0.0;           \
    1271                         sumMatrix2->data.F64[index][k] = 0.0;           \
    1272                         sumMatrixX->data.F64[index][k] = 0.0;           \
    1273                     }                                                   \
    1274                     sumMatrix2->data.F64[index][index] = 1.0;           \
    1275                     sumMatrixX->data.F64[index][normIndex] = 0.0;       \
    1276                     sumMatrixX->data.F64[index][bgIndex] = 0.0;         \
    1277                     sumVector2->data.F64[index] = 0.0;                  \
    12781244                }                                                       \
    12791245            }
    … …  
    13171283            psVectorWriteFile("sumVectorFix.dat", sumVector);
    13181284#endif
    1319         }
     1285        }
    13201286#endif /*** kernel-masking block ***/
    13211287        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
    13221288
    1323 
     1289#ifdef TESTING
    13241290        for (int i = 0; i < solution->n; i++) {
    13251291            fprintf(stderr, "Dual solution %d: %lf\n", i, solution->data.F64[i]);
    13261292        }
     1293#endif
    13271294
    13281295        psFree(sumMatrix);
    … …  
    13831350    pmSubtractionVisualShowFitInit (stamps);
    13841351
    1385     psVector *fSigRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); 
    1386     psVector *fMinRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); 
    1387     psVector *fMaxRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); 
     1352    psVector *fSigRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
     1353    psVector *fMinRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
     1354    psVector *fMaxRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
    13881355
    13891356    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 
    1462             // XXX measure some useful stats on the residuals
    1463             float sum = 0.0;
    1464             float peak = 0.0;
    1465             for (int y = - footprint; y <= footprint; y++) {
    1466                 for (int x = - footprint; x <= footprint; x++) {
    1467                     sum += 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);
    1468                     peak = PS_MAX(peak, 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm));
    1469                 }
    1470             }
    1471 
    1472             // only count pixels with more than X% of the source flux
    1473             // calculate stdev(dflux)
    1474             float dflux1 = 0.0;
    1475             float dflux2 = 0.0;
    1476             int npix = 0;
    1477 
    1478             float dmax = 0.0;
    1479             float dmin = 0.0;
    1480 
    1481             for (int y = - footprint; y <= footprint; y++) {
    1482                 for (int x = - footprint; x <= footprint; x++) {
    1483                     float dflux = 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);
    1484                     if (dflux < 0.02*sum) continue;
    1485                     dflux1 += residual->kernel[y][x];
    1486                     dflux2 += PS_SQR(residual->kernel[y][x]);
    1487                     dmax = PS_MAX(residual->kernel[y][x], dmax);
    1488                     dmin = PS_MIN(residual->kernel[y][x], dmin);
    1489                     npix ++;
    1490                 }
    1491             }
    1492             float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix));
    1493             // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dmax/peak, dmin/peak);
    1494             psVectorAppend(fSigRes, sigma/sum);
    1495             psVectorAppend(fMaxRes, dmax/peak);
    1496             psVectorAppend(fMinRes, dmin/peak);
    1497         } else {
    1498             // Dual convolution
    1499             psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image
    1500             psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image
    1501             psKernel *image1 = stamp->image1; // The first image
    1502             psKernel *image2 = stamp->image2; // The second image
    1503 
    1504             for (int j = 0; j < numKernels; j++) {
    1505                 psKernel *conv1 = convolutions1->data[j]; // Convolution of first image
    1506                 psKernel *conv2 = convolutions2->data[j]; // Convolution of second image
    1507                 double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1
    1508                 double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2
    1509 
    1510                 for (int y = - footprint; y <= footprint; y++) {
    1511                     for (int x = - footprint; x <= footprint; x++) {
    1512                         residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 + conv1->kernel[y][x] * coeff1;
    1513                     }
    1514                 }
    1515             }
    1516 
    1517             // XXX visualize the target, source, convolution and residual
    1518             pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i);
    1519 
    1520             for (int y = - footprint; y <= footprint; y++) {
    1521                 for (int x = - footprint; x <= footprint; x++) {
    1522                     residual->kernel[y][x] += background + image1->kernel[y][x] * norm - image2->kernel[y][x];
    1523                 }
    1524             }
    1525         }
    1526 
    1527         double deviation = 0.0;         // Sum of differences
    1528         for (int y = - footprint; y <= footprint; y++) {
    1529             for (int x = - footprint; x <= footprint; x++) {
    1530                 double dev = PS_SQR(residual->kernel[y][x]) * weight->kernel[y][x];
    1531                 deviation += dev;
     1357        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
     1358        if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
     1359            deviations->data.F32[i] = NAN;
     1360            continue;
     1361        }
     1362
     1363        // Calculate coefficients of the kernel basis functions
     1364        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
     1365        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
     1366        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
     1367
     1368        // Calculate residuals
     1369        psKernel *weight = stamp->weight; // Weight postage stamp
     1370        psImageInit(residual->image, 0.0);
     1371        if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
     1372            psKernel *target;           // Target postage stamp
     1373            psKernel *source;           // Source postage stamp
     1374            psArray *convolutions;      // Convolution postage stamps for each kernel basis function
     1375            switch (kernels->mode) {
     1376              case PM_SUBTRACTION_MODE_1:
     1377                target = stamp->image2;
     1378                source = stamp->image1;
     1379                convolutions = stamp->convolutions1;
     1380
     1381                // Having convolved image1 and changed its normalisation, we need to renormalise the residual
     1382                // so that it is on the scale of image1.
     1383                psImage *image = pmSubtractionKernelImage(kernels, stamp->xNorm, stamp->yNorm,
     1384                                                          false); // Kernel image
     1385                if (!image) {
     1386                    psError(PS_ERR_UNKNOWN, false, "Unable to generate image of kernel.");
     1387                    return false;
     1388                }
     1389                double sumKernel = 0;   // Sum of kernel, for normalising residual
     1390                int size = kernels->size; // Half-size of kernel
     1391                int fullSize = 2 * size + 1; // Full size of kernel
     1392                for (int y = 0; y < fullSize; y++) {
     1393                    for (int x = 0; x < fullSize; x++) {
     1394                        sumKernel += image->data.F32[y][x];
     1395                    }
     1396                }
     1397                psFree(image);
     1398                devNorm = 1.0 / sumKernel / numPixels;
     1399                break;
     1400              case PM_SUBTRACTION_MODE_2:
     1401                target = stamp->image1;
     1402                source = stamp->image2;
     1403                convolutions = stamp->convolutions2;
     1404                break;
     1405              default:
     1406                psAbort("Unsupported subtraction mode: %x", kernels->mode);
     1407            }
     1408
     1409            for (int j = 0; j < numKernels; j++) {
     1410                psKernel *convolution = convolutions->data[j]; // Convolution
     1411                double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j,
     1412                                                                  false); // Coefficient
     1413                for (int y = - footprint; y <= footprint; y++) {
     1414                    for (int x = - footprint; x <= footprint; x++) {
     1415                        residual->kernel[y][x] += convolution->kernel[y][x] * coefficient;
     1416                    }
     1417                }
     1418            }
     1419
     1420            // XXX visualize the target, source, convolution and residual
     1421            pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i);
     1422
     1423            for (int y = - footprint; y <= footprint; y++) {
     1424                for (int x = - footprint; x <= footprint; x++) {
     1425                    residual->kernel[y][x] += background + source->kernel[y][x] * norm - target->kernel[y][x];
     1426                }
     1427            }
     1428
     1429            // XXX measure some useful stats on the residuals
     1430            float sum = 0.0;
     1431            float peak = 0.0;
     1432            for (int y = - footprint; y <= footprint; y++) {
     1433                for (int x = - footprint; x <= footprint; x++) {
     1434                    sum += 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);
     1435                    peak = PS_MAX(peak, 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm));
     1436                }
     1437            }
     1438
     1439            // only count pixels with more than X% of the source flux
     1440            // calculate stdev(dflux)
     1441            float dflux1 = 0.0;
     1442            float dflux2 = 0.0;
     1443            int npix = 0;
     1444
     1445            float dmax = 0.0;
     1446            float dmin = 0.0;
     1447
     1448            for (int y = - footprint; y <= footprint; y++) {
     1449                for (int x = - footprint; x <= footprint; x++) {
     1450                    float dflux = 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);
     1451                    if (dflux < 0.02*sum) continue;
     1452                    dflux1 += residual->kernel[y][x];
     1453                    dflux2 += PS_SQR(residual->kernel[y][x]);
     1454                    dmax = PS_MAX(residual->kernel[y][x], dmax);
     1455                    dmin = PS_MIN(residual->kernel[y][x], dmin);
     1456                    npix ++;
     1457                }
     1458            }
     1459            float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix));
     1460            // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dmax/peak, dmin/peak);
     1461            psVectorAppend(fSigRes, sigma/sum);
     1462            psVectorAppend(fMaxRes, dmax/peak);
     1463            psVectorAppend(fMinRes, dmin/peak);
     1464        } else {
     1465            // Dual convolution
     1466            psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image
     1467            psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image
     1468            psKernel *image1 = stamp->image1; // The first image
     1469            psKernel *image2 = stamp->image2; // The second image
     1470
     1471            for (int j = 0; j < numKernels; j++) {
     1472                psKernel *conv1 = convolutions1->data[j]; // Convolution of first image
     1473                psKernel *conv2 = convolutions2->data[j]; // Convolution of second image
     1474                double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1
     1475                double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2
     1476
     1477                for (int y = - footprint; y <= footprint; y++) {
     1478                    for (int x = - footprint; x <= footprint; x++) {
     1479                        residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 + conv1->kernel[y][x] * coeff1;
     1480                    }
     1481                }
     1482            }
     1483
     1484            // XXX visualize the target, source, convolution and residual
     1485            pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i);
     1486
     1487            for (int y = - footprint; y <= footprint; y++) {
     1488                for (int x = - footprint; x <= footprint; x++) {
     1489                    residual->kernel[y][x] += background + image1->kernel[y][x] * norm - image2->kernel[y][x];
     1490                }
     1491            }
     1492        }
     1493
     1494        double deviation = 0.0;         // Sum of differences
     1495        for (int y = - footprint; y <= footprint; y++) {
     1496            for (int x = - footprint; x <= footprint; x++) {
     1497                double dev = PS_SQR(residual->kernel[y][x]) * weight->kernel[y][x];
     1498                deviation += dev;
    15321499#ifdef TESTING
    1533                 residual->kernel[y][x] = dev;
    1534 #endif
    1535             }
    1536         }
    1537         deviations->data.F32[i] = devNorm * deviation;
    1538         psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n",
    1539                 i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]);
    1540         if (!isfinite(deviations->data.F32[i])) {
    1541             stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
    1542             psTrace("psModules.imcombine", 5,
    1543                     "Rejecting stamp %d (%d,%d) because of non-finite deviation\n",
    1544                     i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
    1545             continue;
    1546         }
     1500                residual->kernel[y][x] = dev;
     1501#endif
     1502            }
     1503        }
     1504        deviations->data.F32[i] = devNorm * deviation;
     1505        psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n",
     1506                i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]);
     1507        if (!isfinite(deviations->data.F32[i])) {
     1508            stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
     1509            psTrace("psModules.imcombine", 5,
     1510                    "Rejecting stamp %d (%d,%d) because of non-finite deviation\n",
     1511                    i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
     1512            continue;
     1513        }
    15471514
    15481515#ifdef TESTING
    1549         {
    1550             psString filename = NULL;
    1551             psStringAppend(&filename, "resid_%03d.fits", i);
    1552             psFits *fits = psFitsOpen(filename, "w");
    1553             psFree(filename);
    1554             psFitsWriteImage(fits, NULL, residual->image, 0, NULL);
    1555             psFitsClose(fits);
    1556         }
    1557         if (stamp->image1) {
    1558             psString filename = NULL;
    1559             psStringAppend(&filename, "stamp_image1_%03d.fits", i);
    1560             psFits *fits = psFitsOpen(filename, "w");
    1561             psFree(filename);
    1562             psFitsWriteImage(fits, NULL, stamp->image1->image, 0, NULL);
    1563             psFitsClose(fits);
    1564         }
    1565         if (stamp->image2) {
    1566             psString filename = NULL;
    1567             psStringAppend(&filename, "stamp_image2_%03d.fits", i);
    1568             psFits *fits = psFitsOpen(filename, "w");
    1569             psFree(filename);
    1570             psFitsWriteImage(fits, NULL, stamp->image2->image, 0, NULL);
    1571             psFitsClose(fits);
    1572         }
    1573         if (stamp->weight) {
    1574             psString filename = NULL;
    1575             psStringAppend(&filename, "stamp_weight_%03d.fits", i);
    1576             psFits *fits = psFitsOpen(filename, "w");
    1577             psFree(filename);
    1578             psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL);
    1579             psFitsClose(fits);
    1580         }
    1581 #endif
    1582    
     1516        {
     1517            psString filename = NULL;
     1518            psStringAppend(&filename, "resid_%03d.fits", i);
     1519            psFits *fits = psFitsOpen(filename, "w");
     1520            psFree(filename);
     1521            psFitsWriteImage(fits, NULL, residual->image, 0, NULL);
     1522            psFitsClose(fits);
     1523        }
     1524        if (stamp->image1) {
     1525            psString filename = NULL;
     1526            psStringAppend(&filename, "stamp_image1_%03d.fits", i);
     1527            psFits *fits = psFitsOpen(filename, "w");
     1528            psFree(filename);
     1529            psFitsWriteImage(fits, NULL, stamp->image1->image, 0, NULL);
     1530            psFitsClose(fits);
     1531        }
     1532        if (stamp->image2) {
     1533            psString filename = NULL;
     1534            psStringAppend(&filename, "stamp_image2_%03d.fits", i);
     1535            psFits *fits = psFitsOpen(filename, "w");
     1536            psFree(filename);
     1537            psFitsWriteImage(fits, NULL, stamp->image2->image, 0, NULL);
     1538            psFitsClose(fits);
     1539        }
     1540        if (stamp->weight) {
     1541            psString filename = NULL;
     1542            psStringAppend(&filename, "stamp_weight_%03d.fits", i);
     1543            psFits *fits = psFitsOpen(filename, "w");
     1544            psFree(filename);
     1545            psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL);
     1546            psFitsClose(fits);
     1547        }
     1548#endif
     1549
    15831550    }
    15841551
    15851552    // calculate and report the normalization and background for the image center
    1586     { 
    1587         polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, 0.0, 0.0);
    1588         double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
    1589         double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
    1590         psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f, background: %f", norm, background);
    1591 
    1592         pmSubtractionVisualShowFit(norm);
    1593         pmSubtractionVisualPlotFit(kernels);
    1594 
    1595         psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    1596         psVectorStats (stats, fSigRes, NULL, NULL, 0);
    1597         kernels->fSigResMean = stats->robustMedian;
    1598         kernels->fSigResStdev = stats->robustStdev;
    1599 
    1600         psStatsInit (stats);
    1601         psVectorStats (stats, fMaxRes, NULL, NULL, 0);
    1602         kernels->fMaxResMean = stats->robustMedian;
    1603         kernels->fMaxResStdev = stats->robustStdev;
    1604 
    1605         psStatsInit (stats);
    1606         psVectorStats (stats, fMinRes, NULL, NULL, 0);
    1607         kernels->fMinResMean = stats->robustMedian;
    1608         kernels->fMinResStdev = stats->robustStdev;
    1609 
    1610         // XXX save these values somewhere
    1611         psLogMsg("psModules.imcombine", PS_LOG_INFO, "fSigma: %f +/- %f, fMaxRes: %f +/- %f, fMinRes: %f +/- %f",
    1612                  kernels->fSigResMean, kernels->fSigResStdev,
    1613                  kernels->fMaxResMean, kernels->fMaxResStdev,
    1614                 kernels->fMinResMean, kernels->fMinResStdev);
    1615 
    1616         psFree (fSigRes);
    1617         psFree (fMaxRes);
    1618         psFree (fMinRes);
    1619         psFree (stats);
     1553    {
     1554        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, 0.0, 0.0);
     1555        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
     1556        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
     1557        psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f, background: %f", norm, background);
     1558
     1559        pmSubtractionVisualShowFit(norm);
     1560        pmSubtractionVisualPlotFit(kernels);
     1561
     1562        psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     1563        psVectorStats (stats, fSigRes, NULL, NULL, 0);
     1564        kernels->fSigResMean = stats->robustMedian;
     1565        kernels->fSigResStdev = stats->robustStdev;
     1566
     1567        psStatsInit (stats);
     1568        psVectorStats (stats, fMaxRes, NULL, NULL, 0);
     1569        kernels->fMaxResMean = stats->robustMedian;
     1570        kernels->fMaxResStdev = stats->robustStdev;
     1571
     1572        psStatsInit (stats);
     1573        psVectorStats (stats, fMinRes, NULL, NULL, 0);
     1574        kernels->fMinResMean = stats->robustMedian;
     1575        kernels->fMinResStdev = stats->robustStdev;
     1576
     1577        // XXX save these values somewhere
     1578        psLogMsg("psModules.imcombine", PS_LOG_INFO, "fSigma: %f +/- %f, fMaxRes: %f +/- %f, fMinRes: %f +/- %f",
     1579                 kernels->fSigResMean, kernels->fSigResStdev,
     1580                 kernels->fMaxResMean, kernels->fMaxResStdev,
     1581                kernels->fMinResMean, kernels->fMinResStdev);
     1582
     1583        psFree (fSigRes);
     1584        psFree (fMaxRes);
     1585        psFree (fMinRes);
     1586        psFree (stats);
    16201587    }
    16211588
    … …  
    16371604
    16381605    for (int i = 0; i < wUt->numCols; i++) {
    1639         for (int j = 0; j < wUt->numRows; j++) {
    1640             if (!isfinite(w->data.F64[j])) continue;
    1641             if (w->data.F64[j] == 0.0) continue;
    1642             wUt->data.F64[j][i] = U->data.F64[i][j] / w->data.F64[j];
    1643         }
     1606        for (int j = 0; j < wUt->numRows; j++) {
     1607            if (!isfinite(w->data.F64[j])) continue;
     1608            if (w->data.F64[j] == 0.0) continue;
     1609            wUt->data.F64[j][i] = U->data.F64[i][j] / w->data.F64[j];
     1610        }
    16441611    }
    16451612    return wUt;
    … …  
    16541621
    16551622    for (int i = 0; i < Ainv->numCols; i++) {
    1656         for (int j = 0; j < Ainv->numRows; j++) {
    1657             double sum = 0.0;
    1658             for (int k = 0; k < V->numCols; k++) {
    1659                 sum += V->data.F64[j][k] * wUt->data.F64[k][i];
    1660             }
    1661             Ainv->data.F64[j][i] = sum;
    1662         }
     1623        for (int j = 0; j < Ainv->numRows; j++) {
     1624            double sum = 0.0;
     1625            for (int k = 0; k < V->numCols; k++) {
     1626                sum += V->data.F64[j][k] * wUt->data.F64[k][i];
     1627            }
     1628            Ainv->data.F64[j][i] = sum;
     1629        }
    16631630    }
    16641631    return Ainv;
    … …  
    16731640
    16741641    for (int i = 0; i < U->numCols; i++) {
    1675         double sum = 0.0;
    1676         for (int j = 0; j < U->numRows; j++) {
    1677             sum += B->data.F64[j] * U->data.F64[j][i];
    1678         }
    1679         UtB[0]->data.F64[i] = sum;
     1642        double sum = 0.0;
     1643        for (int j = 0; j < U->numRows; j++) {
     1644            sum += B->data.F64[j] * U->data.F64[j][i];
     1645        }
     1646        UtB[0]->data.F64[i] = sum;
    16801647    }
    16811648    return true;
    … …  
    16921659
    16931660    for (int i = 0; i < w->n; i++) {
    1694         if (!isfinite(w->data.F64[i])) continue;
    1695         if (w->data.F64[i] == 0.0) continue;
    1696         wUtB[0]->data.F64[i] = UtB->data.F64[i] / w->data.F64[i];
     1661        if (!isfinite(w->data.F64[i])) continue;
     1662        if (w->data.F64[i] == 0.0) continue;
     1663        wUtB[0]->data.F64[i] = UtB->data.F64[i] / w->data.F64[i];
    16971664    }
    16981665    return true;
    … …  
    17071674
    17081675    for (int j = 0; j < V->numRows; j++) {
    1709         double sum = 0.0;
    1710         for (int i = 0; i < V->numCols; i++) {
    1711             sum += V->data.F64[j][i] * wUtB->data.F64[i];
    1712         }
    1713         VwUtB[0]->data.F64[j] = sum;
     1676        double sum = 0.0;
     1677        for (int i = 0; i < V->numCols; i++) {
     1678            sum += V->data.F64[j][i] * wUtB->data.F64[i];
     1679        }
     1680        VwUtB[0]->data.F64[j] = sum;
    17141681    }
    17151682    return true;
    … …  
    17241691
    17251692    for (int j = 0; j < A->numRows; j++) {
    1726         double sum = 0.0;
    1727         for (int i = 0; i < A->numCols; i++) {
    1728             sum += A->data.F64[j][i] * x->data.F64[i];
    1729         }
    1730         B[0]->data.F64[j] = sum;
     1693        double sum = 0.0;
     1694        for (int i = 0; i < A->numCols; i++) {
     1695            sum += A->data.F64[j][i] * x->data.F64[i];
     1696        }
     1697        B[0]->data.F64[j] = sum;
    17311698    }
    17321699    return true;
    … …  
    17401707    double sum = 0.0;
    17411708    for (int i = 0; i < x->n; i++) {
    1742         sum += x->data.F64[i] * y->data.F64[i];
     1709        sum += x->data.F64[i] * y->data.F64[i];
    17431710    }
    17441711    *value = sum;
    … …  
    17531720    for (int i = 0; i < stamps->num; i++) {
    17541721
    1755         pmSubtractionStamp *stamp = stamps->stamps->data[i];
    1756         if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
    1757 
    1758         psKernel *weight = NULL;
    1759         psKernel *window = NULL;
    1760         psKernel *input = NULL;
     1722        pmSubtractionStamp *stamp = stamps->stamps->data[i];
     1723        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
     1724
     1725        psKernel *weight = NULL;
     1726        psKernel *window = NULL;
     1727        psKernel *input = NULL;
    17611728
    17621729#ifdef USE_WEIGHT
    1763         weight = stamp->weight;
     1730        weight = stamp->weight;
    17641731#endif
    17651732#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             input = stamp->image2;
    1773             break;
    1774           case PM_SUBTRACTION_MODE_2:
    1775             input = stamp->image1;
    1776             break;
    1777           default:
    1778             psAbort ("programming error");
    1779         }
    1780 
    1781         for (int y = - footprint; y <= footprint; y++) {
    1782             for (int x = - footprint; x <= footprint; x++) {
    1783                 double in = input->kernel[y][x];
    1784                 double value = in*in;
    1785                 if (weight) {
    1786                     float wtVal = weight->kernel[y][x];
    1787                     value *= wtVal;
    1788                 }
    1789                 if (window) {
    1790                     float  winVal = window->kernel[y][x];
    1791                     value *= winVal;
    1792                 }
    1793                 sum += value;
    1794             }
    1795         }
     1733        window = stamps->window;
     1734#endif
     1735
     1736        switch (kernels->mode) {
     1737            // MODE_1 : convolve image 1 to match image 2 (and vice versa)
     1738          case PM_SUBTRACTION_MODE_1:
     1739            input = stamp->image2;
     1740            break;
     1741          case PM_SUBTRACTION_MODE_2:
     1742            input = stamp->image1;
     1743            break;
     1744          default:
     1745            psAbort ("programming error");
     1746        }
     1747
     1748        for (int y = - footprint; y <= footprint; y++) {
     1749            for (int x = - footprint; x <= footprint; x++) {
     1750                double in = input->kernel[y][x];
     1751                double value = in*in;
     1752                if (weight) {
     1753                    float wtVal = weight->kernel[y][x];
     1754                    value *= wtVal;
     1755                }
     1756                if (window) {
     1757                    float  winVal = window->kernel[y][x];
     1758                    value *= winVal;
     1759                }
     1760                sum += value;
     1761            }
     1762        }
    17961763    }
    17971764    *y2 = sum;
    … …  
    18131780    for (int i = 0; i < stamps->num; i++) {
    18141781
    1815         pmSubtractionStamp *stamp = stamps->stamps->data[i];
    1816         if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
    1817 
    1818         psKernel *weight = NULL;
    1819         psKernel *window = NULL;
    1820         psKernel *target = NULL;
    1821         psKernel *source = NULL;
    1822 
    1823         psArray *convolutions = NULL;
     1782        pmSubtractionStamp *stamp = stamps->stamps->data[i];
     1783        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
     1784
     1785        psKernel *weight = NULL;
     1786        psKernel *window = NULL;
     1787        psKernel *target = NULL;
     1788        psKernel *source = NULL;
     1789
     1790        psArray *convolutions = NULL;
    18241791
    18251792#ifdef USE_WEIGHT
    1826         weight = stamp->weight;
     1793        weight = stamp->weight;
    18271794#endif
    18281795#ifdef USE_WINDOW
    1829         window = stamps->window;
    1830 #endif
    1831 
    1832         switch (kernels->mode) {
    1833             // MODE_1 : convolve image 1 to match image 2 (and vice versa)
    1834           case PM_SUBTRACTION_MODE_1:
    1835             target = stamp->image2;
    1836             source = stamp->image1;
    1837             convolutions = stamp->convolutions1;
    1838             break;
    1839           case PM_SUBTRACTION_MODE_2:
    1840             target = stamp->image1;
    1841             source = stamp->image2;
    1842             convolutions = stamp->convolutions2;
    1843             break;
    1844           default:
    1845             psAbort ("programming error");
    1846         }
    1847 
    1848         // Calculate coefficients of the kernel basis functions
    1849         polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
    1850         double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
    1851         double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
    1852 
    1853         psImageInit(residual->image, 0.0);
    1854         for (int j = 0; j < numKernels; j++) {
    1855             psKernel *convolution = convolutions->data[j]; // Convolution
    1856             double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient
    1857             for (int y = - footprint; y <= footprint; y++) {
    1858                 for (int x = - footprint; x <= footprint; x++) {
    1859                     residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient;
    1860                 }
    1861             }
    1862         }
    1863 
    1864         for (int y = - footprint; y <= footprint; y++) {
    1865             for (int x = - footprint; x <= footprint; x++) {
    1866                 double resid = target->kernel[y][x] - background - source->kernel[y][x] * norm + residual->kernel[y][x];
    1867                 double value = PS_SQR(resid);
    1868                 if (weight) {
    1869                     float wtVal = weight->kernel[y][x];
    1870                     value *= wtVal;
    1871                 }
    1872                 if (window) {
    1873                     float  winVal = window->kernel[y][x];
    1874                     value *= winVal;
    1875                 }
    1876                 sum += value;
    1877             }
    1878         }
     1796        window = stamps->window;
     1797#endif
     1798
     1799        switch (kernels->mode) {
     1800            // MODE_1 : convolve image 1 to match image 2 (and vice versa)
     1801          case PM_SUBTRACTION_MODE_1:
     1802            target = stamp->image2;
     1803            source = stamp->image1;
     1804            convolutions = stamp->convolutions1;
     1805            break;
     1806          case PM_SUBTRACTION_MODE_2:
     1807            target = stamp->image1;
     1808            source = stamp->image2;
     1809            convolutions = stamp->convolutions2;
     1810            break;
     1811          default:
     1812            psAbort ("programming error");
     1813        }
     1814
     1815        // Calculate coefficients of the kernel basis functions
     1816        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
     1817        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
     1818        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
     1819
     1820        psImageInit(residual->image, 0.0);
     1821        for (int j = 0; j < numKernels; j++) {
     1822            psKernel *convolution = convolutions->data[j]; // Convolution
     1823            double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient
     1824            for (int y = - footprint; y <= footprint; y++) {
     1825                for (int x = - footprint; x <= footprint; x++) {
     1826                    residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient;
     1827                }
     1828            }
     1829        }
     1830
     1831        for (int y = - footprint; y <= footprint; y++) {
     1832            for (int x = - footprint; x <= footprint; x++) {
     1833                double resid = target->kernel[y][x] - background - source->kernel[y][x] * norm + residual->kernel[y][x];
     1834                double value = PS_SQR(resid);
     1835                if (weight) {
     1836                    float wtVal = weight->kernel[y][x];
     1837                    value *= wtVal;
     1838                }
     1839                if (window) {
     1840                    float  winVal = window->kernel[y][x];
     1841                    value *= winVal;
     1842                }
     1843                sum += value;
     1844            }
     1845        }
    18791846    }
    18801847    psFree (polyValues);
    … …  
    18871854
    18881855    for (int i = 0; i < w->n; i++) {
    1889         wApply->data.F64[i] = wMask->data.U8[i] ? 0.0 : w->data.F64[i];
     1856        wApply->data.F64[i] = wMask->data.U8[i] ? 0.0 : w->data.F64[i];
    18901857    }
    18911858    return true;
    … …  
    19021869    // generate Vn = V * w^{-1}
    19031870    for (int j = 0; j < Vn->numRows; j++) {
    1904         for (int i = 0; i < Vn->numCols; i++) {
    1905             if (!isfinite(w->data.F64[i])) continue;
    1906             if (w->data.F64[i] == 0.0) continue;
    1907             Vn->data.F64[j][i] = V->data.F64[j][i] / w->data.F64[i];
    1908         }
     1871        for (int i = 0; i < Vn->numCols; i++) {
     1872            if (!isfinite(w->data.F64[i])) continue;
     1873            if (w->data.F64[i] == 0.0) continue;
     1874            Vn->data.F64[j][i] = V->data.F64[j][i] / w->data.F64[i];
     1875        }
    19091876    }
    19101877
    … …  
    19141881    // generate Xvar = Vn * Vn^T
    19151882    for (int j = 0; j < Vn->numRows; j++) {
    1916         for (int i = 0; i < Vn->numCols; i++) {
    1917             double sum = 0.0;
    1918             for (int k = 0; k < Vn->numCols; k++) {
    1919                 sum += Vn->data.F64[k][i]*Vn->data.F64[k][j];
    1920             }
    1921             Xvar->data.F64[j][i] = sum;
    1922         }
     1883        for (int i = 0; i < Vn->numCols; i++) {
     1884            double sum = 0.0;
     1885            for (int k = 0; k < Vn->numCols; k++) {
     1886                sum += Vn->data.F64[k][i]*Vn->data.F64[k][j];
     1887            }
     1888            Xvar->data.F64[j][i] = sum;
     1889        }
    19231890    }
    19241891    return Xvar;
    … …  
    19351902    psFitsWriteImage(fits, header, image, 0, NULL);
    19361903    psFitsClose(fits);
    1937    
     1904
    19381905    return true;
    19391906}
Note: See TracChangeset for help on using the changeset viewer.