IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26373


Ignore:
Timestamp:
Dec 8, 2009, 7:11:00 PM (16 years ago)
Author:
eugene
Message:

implement SVD analysis of the chi-square equation

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

Legend:

Unmodified
Added
Removed

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

    r26332 r26373  
    2121#include "pmSubtractionStamps.h"
    2222#include "pmSubtractionEquation.h"
     23#include "pmSubtractionVisual.h"
    2324#include "pmSubtractionThreads.h"
    2425
    … …  
    525526    int x0 = - image->xMin, y0 = - image->yMin; // Position of centre of convolved image
    526527    psKernel *convolved = psKernelAllocFromImage(conv, x0, y0); // Kernel version
     528
     529    // pmSubtractionVisualShowSubtraction(image->image, kernel->image, conv);
     530
    527531    psFree(conv);
    528532    return convolved;

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

    r26347 r26373  
    831831}
    832832
     833// private functions used on pmSubtractionSolveEquation
     834bool psVectorWriteFile (char *filename, const psVector *vector);
     835bool psFitsWriteImageSimple (char *filename, psImage *image, psMetadata *header);
     836
     837psImage *p_pmSubSolve_wUt (psVector *w, psImage *U);
     838psImage *p_pmSubSolve_VwUt (psImage *V, psImage *wUt);
     839
     840psVector *p_pmSubSolve_UtB (psImage *U, psVector *B);
     841psVector *p_pmSubSolve_wUtB (psVector *w, psVector *UtB);
     842psVector *p_pmSubSolve_VwUtB (psImage *V, psVector *wUtB);
     843
     844psVector *p_pmSubSolve_Ax (psImage *A, psVector *x);
     845
     846double p_pmSubSolve_y2 (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps);
     847double p_pmSubSolve_VdV (psVector *x, psVector *y);
     848double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps);
     849
    833850bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels,
    834851                                const pmSubtractionStampList *stamps,
    … …  
    943960        {
    944961            psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL);
    945             psFits *fits = psFitsOpen("matrixInv.fits", "w");
    946             psFitsWriteImage(fits, NULL, inverse, 0, NULL);
    947             psFitsClose(fits);
     962            psFitsWriteImageSimple("matrixInv.fits", inverse, NULL);
    948963            psFree(inverse);
    949964        }
    … …  
    952967            psImage *Xt = psMatrixTranspose(NULL, X);
    953968            psImage *XtX = psMatrixMultiply(NULL, Xt, X);
    954             psFits *fits = psFitsOpen("matrixErr.fits", "w");
    955             psFitsWriteImage(fits, NULL, XtX, 0, NULL);
    956             psFitsClose(fits);
     969            psFitsWriteImageSimple("matrixErr.fits", XtX, NULL);
    957970            psFree(X);
    958971            psFree(Xt);
    … …  
    961974#endif
    962975
    963         // XXX test: save the matrix A and vector b:
    964         {
    965             psFits *fits = psFitsOpen("matrix.fits", "w");
    966             psFitsWriteImage(fits, NULL, sumMatrix, 0, NULL);
    967             psFitsClose(fits);
    968 
    969             FILE *f = fopen ("vector.dat", "w");
    970             int fd = fileno(f);
    971             p_psVectorPrint (fd, sumVector, "B");
    972             fclose (f);
    973         }
    974 
    975976# define SVD_ANALYSIS 1
    976         // re-do this section with SVD:
     977# define COEFF_SIG 0.0
     978# define SVD_TOL 1e-9
     979
     980        psVector *solution = NULL;
     981        psVector *solutionErr = NULL;
     982
     983        // Use SVD to determine the kernel coeffs (and validate)
    977984        if (SVD_ANALYSIS) {
    978             // We have sumVector and sumMatrix.  we are trying to solve the equation: sumMatrix * x = sumVector.
     985
     986            // We have sumVector and sumMatrix.  we are trying to solve the following equation:
     987            // sumMatrix * x = sumVector.
     988
    979989            // we can use any standard matrix inversion to solve this.  However, the basis
    980990            // functions in general have substantial correlation, so that the solution may be
    981             // somewhat poorly determiend or unstable.  If not numerically ill-conditioned, the
     991            // somewhat poorly determined or unstable.  If not numerically ill-conditioned, the
    982992            // system of equations may be statistically ill-conditioned.  Noise in the image
    983             // will drive insignificant but correlated terms in the solution.  To avoid these
     993            // will drive insignificant, but correlated, terms in the solution.  To avoid these
    984994            // problems, we can use SVD to identify numerically unconstrained values and to
    985995            // avoid statistically badly determined value.
    … …  
    9931003            // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0
    9941004
    995             // I get confused by the index values between the image vs matrix usage:  In terms
    996             // of the elements of an image A(x,y) = A->data.F32[y][x] = A_x,y, a matrix
    997             // multiplication is: A_k,j * B_i,k = C_i,j
     1005            // XXX TEST with known matrix:
     1006            // psFree (sumMatrix);
     1007            // sumMatrix = psImageAlloc (2, 2, PS_TYPE_F64);
     1008            // sumMatrix->data.F64[0][0] = 1.0;
     1009            // sumMatrix->data.F64[0][1] = 0.5;
     1010            //
     1011            // sumMatrix->data.F64[1][0] = 0.5;
     1012            // sumMatrix->data.F64[1][1] = 2.0;
     1013
     1014            // sumMatrix->data.F64[2][0] = 0.0;
     1015            // sumMatrix->data.F64[2][1] = 1.0;
     1016
     1017            // TEST psFitsWriteImageSimple("A.fits", sumMatrix, NULL);
     1018            // TEST psVectorWriteFile("B.dat", sumVector);;
     1019
     1020            // If I use the SVD trick to re-condition the matrix, I need to break out the
     1021            // kernel terms from the normalization and background terms.
     1022
     1023            // XX int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     1024            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
     1025
     1026            // solve for the normalization first (independent of all other terms)
     1027            // XX double sumRR = sumMatrix->data.F64[normIndex][normIndex];
     1028            // XX double sumIR = sumVector->data.F64[normIndex];
     1029            // XX double norm = sumIR / sumRR;
     1030
     1031            // now pull out the kernel elements into their own square matrix
     1032            psImage  *kernelMatrix = psImageAlloc  (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64);
     1033            psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64);
     1034
     1035            // note that we need to subtract norm*sumRC from the sumVector elements (since we
     1036            // are fitting to the residual image of I - norm*Ref)
     1037            for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) {
     1038                // if (ix == normIndex) continue;
     1039                if (ix == bgIndex) continue;
     1040                for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) {
     1041                    // if (iy == normIndex) continue;
     1042                    if (iy == bgIndex) continue;
     1043                    kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix];
     1044                    ky++;
     1045                }
     1046                // kernelVector->data.F64[kx] = sumVector->data.F64[ix] - norm*sumMatrix->data.F64[normIndex][ix];
     1047                kernelVector->data.F64[kx] = sumVector->data.F64[ix];
     1048                kx++;
     1049            }
     1050
     1051            psFitsWriteImageSimple("A.fits", kernelMatrix, NULL);
    9981052
    9991053            psImage *U = NULL;
    10001054            psImage *V = NULL;
    10011055            psVector *w = NULL;
    1002             if (!psMatrixSVD (&U, &w, &V, sumMatrix)) {
     1056            if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) {
    10031057                psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n");
    10041058                return NULL;
    10051059            }
    10061060
    1007             // calculate A_inverse and x_svd
    1008             // A'_i,j = \sum_k (V_k,j * w_k * Ut_i,k)
    1009             // but U^T_i,j = U_j,i, so:
    1010             // A'_i,j = \sum_k (V_k,j * w_k * U_k,i)
    1011 
    1012             // v1_j = \sum_k (UT_k,j * B_k)
    1013             // v1_j = w_j * \sum_k (U_j,k * B_k)
    1014             // x_i = \sum_j V_j,i * w_j * \sum_k (U_j,k * B_k)
    1015            
    1016             // A'_i,j = \sum_k (V_k,j * w_k * U_k,i)
    1017 
    1018             psImage  *Ainv = psImageAlloc(A->numCols, A->numRows, PS_TYPE_F32);
    1019             psVector *Xsvd = psVectorAlloc(A->numCols, PS_TYPE_F32);
    1020             for (int iy = 0; iy < A->numRows; iy++) {
    1021                 double vsum = 0.0;
    1022                 for (int ix = 0; ix < A->numCols; ix++) {
    1023                     double msum = 0.0;
    1024                     double vsum1 = 0.0;
    1025                     for (int k = 0; k < A->numCols; k++) {
    1026                         if (fabs(w->data.F32[k]) < FLT_EPSILON) continue;
    1027                         msum += V->data.F32[iy][k] * U->data.F32[ix][k] / w->data.F32[k];
    1028                         vsum1 += U->data.F32[k][iy] * B->data.F32[k];
    1029                     }
    1030                     Ainv->data.F32[iy][ix] = sum;
    1031                     if (fabs(w->data.F32[ix]) < FLT_EPSILON) continue;
    1032                     vsum += V->data.F32[iy][ix] * vsum1 / w->data.F32[ix];
    1033                 }
    1034                 Xsvd->data.F32[iy] = vsum;
    1035             }
    1036 
    1037             // compare Xsvd[k] to dXsvd and zero out w entries that are not significant
    1038             # define COEFF_SIG 1.0
    1039             for (int k = 0; k < A->numRows; k++) {
    1040                 float dXsvd = sqrt(Ainv->data.F32[k][k]);
    1041                 if (fabs(Xsvd->data.F32[k] / dXsvd) < COEFF_SIG) {
    1042                     w->data.F32[k] = 0.0;
    1043                 }
    1044             }
    1045            
    1046             //
    1047 
    1048 
    1049         }
    1050 
    1051         psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
    1052         psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
    1053         psFree(sumMatrix);
    1054         if (!luMatrix) {
    1055             psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n");
    1056             psFree(sumVector);
    1057             psFree(luMatrix);
    1058             psFree(permutation);
    1059             return NULL;
    1060         }
    1061 
    1062         psVector *solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation);
    1063         psFree(sumVector);
    1064         psFree(luMatrix);
    1065         psFree(permutation);
    1066         if (!solution) {
    1067             psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
    1068             return NULL;
    1069         }
     1061            // TEST psFitsWriteImageSimple("U.fits", U, NULL);
     1062            // TEST psFitsWriteImageSimple("V.fits", V, NULL);
     1063            // TEST psVectorWriteFile("w.dat", w);
     1064
     1065            // calculate A_inverse:
     1066            // Ainv = V * w * U^T
     1067            psImage *wUt  = p_pmSubSolve_wUt (w, U);
     1068            psImage *Ainv = p_pmSubSolve_VwUt (V, wUt);
     1069
     1070            // TEST psFitsWriteImageSimple("AInv.fits", Ainv, NULL);
     1071            // TEST psImage *Io = p_pmSubSolve_VwUt (Ainv, sumMatrix);
     1072            // TEST psFitsWriteImageSimple("I.fits", Io, NULL);
     1073
     1074            // kernel terms:
     1075            for (int i = 0; i < w->n; i++) {
     1076                fprintf (stderr, "w: %f\n", w->data.F64[i]);
     1077            }
     1078
     1079            // XXX loop over w adding in more and more of the values until chisquare is no
     1080            // longer dropping significantly
     1081
     1082            // zero-out ill-conditioned SVD components:
     1083            // largest w value is in w->data.F64[0]:
     1084            double Wmax = w->data.F64[0];
     1085            for (int k = 1; k < w->n; k++) {
     1086                if (fabs(w->data.F64[k] / Wmax) < SVD_TOL) {
     1087                    fprintf (stderr, "masking basis function %d : %f\n", k, w->data.F64[k] / Wmax);
     1088                    w->data.F64[k] = NAN;
     1089                    psTrace("psModules.imcombine", 5, "masking basis function %d : %f", k, w->data.F64[k] / Wmax);
     1090                }
     1091            }
     1092
     1093            // solve for x:
     1094            // x = V * w * (U^T * B)
     1095            psVector *UtB  = p_pmSubSolve_UtB (U, kernelVector);
     1096            psVector *wUtB = p_pmSubSolve_wUtB (w, UtB);
     1097            psVector *Xsvd = p_pmSubSolve_VwUtB (V, wUtB);
     1098
     1099            // chi-square for this system of equations:
     1100            // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
     1101            // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
     1102            psVector *Ax = p_pmSubSolve_Ax (kernelMatrix, Xsvd);
     1103            double Axx   = p_pmSubSolve_VdV (Ax, Xsvd);
     1104            double Bx    = p_pmSubSolve_VdV (kernelVector, Xsvd);
     1105            double y2    = p_pmSubSolve_y2 (kernels, stamps);
     1106
     1107            // apparently, this works (compare with the brute force value below
     1108            double chiSquare = Axx - 2.0*Bx + y2;
     1109            fprintf (stderr, "chi square = %f\n", chiSquare);
     1110
     1111            // copy the kernel solutions to the full solution vector:
     1112            solution = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     1113            solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     1114
     1115            for (int ix = 0, kx = 0; ix < sumVector->n; ix++) {
     1116                // XX if (ix == normIndex) {
     1117                // XX     // solution->data.F64[ix] = norm; //
     1118                // XX     solution->data.F64[ix] = 1.0; // XXX a CHEAT
     1119                // XX     solutionErr->data.F64[ix] = 0.001;
     1120                // XX     continue;
     1121                // XX }
     1122                if (ix == bgIndex) {
     1123                    solution->data.F64[ix] = 0;
     1124                    solutionErr->data.F64[ix] = 0.001;
     1125                    continue;
     1126                }
     1127                solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]);
     1128                solution->data.F64[ix] = Xsvd->data.F64[kx];
     1129                kx++;
     1130            }
     1131
     1132            // TEST psVectorWriteFile("X.dat", Xsvd);
     1133            // TEST psVector *Bo = p_pmSubSolve_Ax (sumMatrix, Xsvd);
     1134            // TEST psVectorWriteFile("Bo.dat", Bo);
     1135               
     1136            psFree (kernelMatrix);
     1137            psFree (kernelVector);
     1138
     1139            psFree (U);
     1140            psFree (V);
     1141            psFree (w);
     1142
     1143            psFree (wUt);
     1144            psFree (Ainv);
     1145
     1146            psFree (UtB);
     1147            psFree (wUtB);
     1148            psFree (Xsvd);
     1149        } else {
     1150            psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
     1151            psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
     1152            psFree(sumMatrix);
     1153            if (!luMatrix) {
     1154                psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n");
     1155                psFree(sumVector);
     1156                psFree(luMatrix);
     1157                psFree(permutation);
     1158                return NULL;
     1159            }
     1160
     1161            solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation);
     1162            psFree(sumVector);
     1163            psFree(luMatrix);
     1164            psFree(permutation);
     1165            if (!solution) {
     1166                psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
     1167                return NULL;
     1168            }
     1169        }
    10701170
    10711171        if (!kernels->solution1) {
    … …  
    10881188            int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
    10891189            for (int i = 0; i < numKernels * numPoly; i++) {
    1090                 kernels->solution1->data.F64[i] = solution->data.F64[i];
    1091             }
    1092         }
     1190                fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i]));
     1191                if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) {
     1192                    fprintf (stderr, "drop\n");
     1193                    kernels->solution1->data.F64[i] = 0.0;
     1194                } else {
     1195                    fprintf (stderr, "keep\n");
     1196                    kernels->solution1->data.F64[i] = solution->data.F64[i];
     1197                }
     1198            }
     1199        }
     1200        // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps);
     1201        // fprintf (stderr, "chi square Brute = %f\n", chiSquare);
     1202
    10931203        psFree (solution);
    10941204
    … …  
    11411251
    11421252#ifdef TESTING
    1143         {
    1144             psFits *fits = psFitsOpen("sumMatrix.fits", "w");
    1145             psFitsWriteImage(fits, NULL, sumMatrix, 0, NULL);
    1146             psFitsClose(fits);
    1147         }
     1253        psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, NULL);
    11481254        {
    11491255            psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64);
    … …  
    13451451
    13461452    for (int i = 0; i < stamps->num; i++) {
    1347         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
    1348         if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
    1349             deviations->data.F32[i] = NAN;
    1350             continue;
    1351         }
    1352 
    1353         // Calculate coefficients of the kernel basis functions
    1354         polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
    1355         double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
    1356         double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
    1357 
    1358         // Calculate residuals
    1359         psKernel *weight = stamp->weight; // Weight postage stamp
    1360         psImageInit(residual->image, 0.0);
    1361         if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
    1362             psKernel *target;           // Target postage stamp
    1363             psKernel *source;           // Source postage stamp
    1364             psArray *convolutions;      // Convolution postage stamps for each kernel basis function
    1365             switch (kernels->mode) {
    1366               case PM_SUBTRACTION_MODE_1:
    1367                 target = stamp->image2;
    1368                 source = stamp->image1;
    1369                 convolutions = stamp->convolutions1;
    1370 
    1371                 // Having convolved image1 and changed its normalisation, we need to renormalise the residual
    1372                 // so that it is on the scale of image1.
    1373                 psImage *image = pmSubtractionKernelImage(kernels, stamp->xNorm, stamp->yNorm,
    1374                                                           false); // Kernel image
    1375                 if (!image) {
    1376                     psError(PS_ERR_UNKNOWN, false, "Unable to generate image of kernel.");
    1377                     return false;
    1378                 }
    1379                 double sumKernel = 0;   // Sum of kernel, for normalising residual
    1380                 int size = kernels->size; // Half-size of kernel
    1381                 int fullSize = 2 * size + 1; // Full size of kernel
    1382                 for (int y = 0; y < fullSize; y++) {
    1383                     for (int x = 0; x < fullSize; x++) {
    1384                         sumKernel += image->data.F32[y][x];
    1385                     }
    1386                 }
    1387                 psFree(image);
    1388                 devNorm = 1.0 / sumKernel / numPixels;
    1389                 break;
    1390               case PM_SUBTRACTION_MODE_2:
    1391                 target = stamp->image1;
    1392                 source = stamp->image2;
    1393                 convolutions = stamp->convolutions2;
    1394                 break;
    1395               default:
    1396                 psAbort("Unsupported subtraction mode: %x", kernels->mode);
    1397             }
    1398 
    1399             for (int j = 0; j < numKernels; j++) {
    1400                 psKernel *convolution = convolutions->data[j]; // Convolution
    1401                 double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j,
    1402                                                                   false); // Coefficient
    1403                 for (int y = - footprint; y <= footprint; y++) {
    1404                     for (int x = - footprint; x <= footprint; x++) {
    1405                         residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient;
    1406                     }
    1407                 }
    1408             }
     1453        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
     1454        if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
     1455            deviations->data.F32[i] = NAN;
     1456            continue;
     1457        }
     1458
     1459        // Calculate coefficients of the kernel basis functions
     1460        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
     1461        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
     1462        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
     1463
     1464        // Calculate residuals
     1465        psKernel *weight = stamp->weight; // Weight postage stamp
     1466        psImageInit(residual->image, 0.0);
     1467        if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
     1468            psKernel *target;           // Target postage stamp
     1469            psKernel *source;           // Source postage stamp
     1470            psArray *convolutions;      // Convolution postage stamps for each kernel basis function
     1471            switch (kernels->mode) {
     1472              case PM_SUBTRACTION_MODE_1:
     1473                target = stamp->image2;
     1474                source = stamp->image1;
     1475                convolutions = stamp->convolutions1;
     1476
     1477                // Having convolved image1 and changed its normalisation, we need to renormalise the residual
     1478                // so that it is on the scale of image1.
     1479                psImage *image = pmSubtractionKernelImage(kernels, stamp->xNorm, stamp->yNorm,
     1480                                                          false); // Kernel image
     1481                if (!image) {
     1482                    psError(PS_ERR_UNKNOWN, false, "Unable to generate image of kernel.");
     1483                    return false;
     1484                }
     1485                double sumKernel = 0;   // Sum of kernel, for normalising residual
     1486                int size = kernels->size; // Half-size of kernel
     1487                int fullSize = 2 * size + 1; // Full size of kernel
     1488                for (int y = 0; y < fullSize; y++) {
     1489                    for (int x = 0; x < fullSize; x++) {
     1490                        sumKernel += image->data.F32[y][x];
     1491                    }
     1492                }
     1493                psFree(image);
     1494                devNorm = 1.0 / sumKernel / numPixels;
     1495                break;
     1496              case PM_SUBTRACTION_MODE_2:
     1497                target = stamp->image1;
     1498                source = stamp->image2;
     1499                convolutions = stamp->convolutions2;
     1500                break;
     1501              default:
     1502                psAbort("Unsupported subtraction mode: %x", kernels->mode);
     1503            }
     1504
     1505            for (int j = 0; j < numKernels; j++) {
     1506                psKernel *convolution = convolutions->data[j]; // Convolution
     1507                double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j,
     1508                                                                  false); // Coefficient
     1509                for (int y = - footprint; y <= footprint; y++) {
     1510                    for (int x = - footprint; x <= footprint; x++) {
     1511                        residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient;
     1512                    }
     1513                }
     1514            }
    14091515
    14101516            // XXX visualize the target, source, convolution and residual
    14111517            pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i);
    14121518
    1413             for (int y = - footprint; y <= footprint; y++) {
    1414                 for (int x = - footprint; x <= footprint; x++) {
    1415                     residual->kernel[y][x] += target->kernel[y][x] - background - source->kernel[y][x] * norm;
    1416                 }
    1417             }
    1418         } else {
    1419             // Dual convolution
    1420             psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image
    1421             psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image
    1422             psKernel *image1 = stamp->image1; // The first image
    1423             psKernel *image2 = stamp->image2; // The second image
    1424 
    1425             for (int j = 0; j < numKernels; j++) {
    1426                 psKernel *conv1 = convolutions1->data[j]; // Convolution of first image
    1427                 psKernel *conv2 = convolutions2->data[j]; // Convolution of second image
    1428                 double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1
    1429                 double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2
    1430 
    1431                 for (int y = - footprint; y <= footprint; y++) {
    1432                     for (int x = - footprint; x <= footprint; x++) {
    1433                         residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 - conv1->kernel[y][x] * coeff1;
    1434                     }
    1435                 }
    1436             }
    1437             for (int y = - footprint; y <= footprint; y++) {
    1438                 for (int x = - footprint; x <= footprint; x++) {
    1439                     residual->kernel[y][x] += image2->kernel[y][x] - background - image1->kernel[y][x] * norm;
    1440                 }
    1441             }
    1442         }
    1443 
    1444         double deviation = 0.0;         // Sum of differences
    1445         for (int y = - footprint; y <= footprint; y++) {
    1446             for (int x = - footprint; x <= footprint; x++) {
    1447                 double dev = PS_SQR(residual->kernel[y][x]) * weight->kernel[y][x];
    1448                 deviation += dev;
     1519            for (int y = - footprint; y <= footprint; y++) {
     1520                for (int x = - footprint; x <= footprint; x++) {
     1521                    residual->kernel[y][x] += target->kernel[y][x] - background - source->kernel[y][x] * norm;
     1522                }
     1523            }
     1524        } else {
     1525            // Dual convolution
     1526            psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image
     1527            psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image
     1528            psKernel *image1 = stamp->image1; // The first image
     1529            psKernel *image2 = stamp->image2; // The second image
     1530
     1531            for (int j = 0; j < numKernels; j++) {
     1532                psKernel *conv1 = convolutions1->data[j]; // Convolution of first image
     1533                psKernel *conv2 = convolutions2->data[j]; // Convolution of second image
     1534                double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1
     1535                double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2
     1536
     1537                for (int y = - footprint; y <= footprint; y++) {
     1538                    for (int x = - footprint; x <= footprint; x++) {
     1539                        residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 - conv1->kernel[y][x] * coeff1;
     1540                    }
     1541                }
     1542            }
     1543            for (int y = - footprint; y <= footprint; y++) {
     1544                for (int x = - footprint; x <= footprint; x++) {
     1545                    residual->kernel[y][x] += image2->kernel[y][x] - background - image1->kernel[y][x] * norm;
     1546                }
     1547            }
     1548        }
     1549
     1550        double deviation = 0.0;         // Sum of differences
     1551        for (int y = - footprint; y <= footprint; y++) {
     1552            for (int x = - footprint; x <= footprint; x++) {
     1553                double dev = PS_SQR(residual->kernel[y][x]) * weight->kernel[y][x];
     1554                deviation += dev;
    14491555#ifdef TESTING
    1450                 residual->kernel[y][x] = dev;
    1451 #endif
    1452             }
    1453         }
    1454         deviations->data.F32[i] = devNorm * deviation;
    1455         psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n",
    1456                 i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]);
    1457         if (!isfinite(deviations->data.F32[i])) {
    1458             stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
    1459             psTrace("psModules.imcombine", 5,
    1460                     "Rejecting stamp %d (%d,%d) because of non-finite deviation\n",
    1461                     i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
    1462             continue;
    1463         }
     1556                residual->kernel[y][x] = dev;
     1557#endif
     1558            }
     1559        }
     1560        deviations->data.F32[i] = devNorm * deviation;
     1561        psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n",
     1562                i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]);
     1563        if (!isfinite(deviations->data.F32[i])) {
     1564            stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
     1565            psTrace("psModules.imcombine", 5,
     1566                    "Rejecting stamp %d (%d,%d) because of non-finite deviation\n",
     1567                    i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
     1568            continue;
     1569        }
    14641570
    14651571#ifdef TESTING
    1466         {
    1467             psString filename = NULL;
    1468             psStringAppend(&filename, "resid_%03d.fits", i);
    1469             psFits *fits = psFitsOpen(filename, "w");
    1470             psFree(filename);
    1471             psFitsWriteImage(fits, NULL, residual->image, 0, NULL);
    1472             psFitsClose(fits);
    1473         }
    1474         if (stamp->image1) {
    1475             psString filename = NULL;
    1476             psStringAppend(&filename, "stamp_image1_%03d.fits", i);
    1477             psFits *fits = psFitsOpen(filename, "w");
    1478             psFree(filename);
    1479             psFitsWriteImage(fits, NULL, stamp->image1->image, 0, NULL);
    1480             psFitsClose(fits);
    1481         }
    1482         if (stamp->image2) {
    1483             psString filename = NULL;
    1484             psStringAppend(&filename, "stamp_image2_%03d.fits", i);
    1485             psFits *fits = psFitsOpen(filename, "w");
    1486             psFree(filename);
    1487             psFitsWriteImage(fits, NULL, stamp->image2->image, 0, NULL);
    1488             psFitsClose(fits);
    1489         }
    1490         if (stamp->weight) {
    1491             psString filename = NULL;
    1492             psStringAppend(&filename, "stamp_weight_%03d.fits", i);
    1493             psFits *fits = psFitsOpen(filename, "w");
    1494             psFree(filename);
    1495             psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL);
    1496             psFitsClose(fits);
    1497         }
     1572        {
     1573            psString filename = NULL;
     1574            psStringAppend(&filename, "resid_%03d.fits", i);
     1575            psFits *fits = psFitsOpen(filename, "w");
     1576            psFree(filename);
     1577            psFitsWriteImage(fits, NULL, residual->image, 0, NULL);
     1578            psFitsClose(fits);
     1579        }
     1580        if (stamp->image1) {
     1581            psString filename = NULL;
     1582            psStringAppend(&filename, "stamp_image1_%03d.fits", i);
     1583            psFits *fits = psFitsOpen(filename, "w");
     1584            psFree(filename);
     1585            psFitsWriteImage(fits, NULL, stamp->image1->image, 0, NULL);
     1586            psFitsClose(fits);
     1587        }
     1588        if (stamp->image2) {
     1589            psString filename = NULL;
     1590            psStringAppend(&filename, "stamp_image2_%03d.fits", i);
     1591            psFits *fits = psFitsOpen(filename, "w");
     1592            psFree(filename);
     1593            psFitsWriteImage(fits, NULL, stamp->image2->image, 0, NULL);
     1594            psFitsClose(fits);
     1595        }
     1596        if (stamp->weight) {
     1597            psString filename = NULL;
     1598            psStringAppend(&filename, "stamp_weight_%03d.fits", i);
     1599            psFits *fits = psFitsOpen(filename, "w");
     1600            psFree(filename);
     1601            psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL);
     1602            psFitsClose(fits);
     1603        }
    14981604#endif
    14991605   
    … …  
    15161622    return deviations;
    15171623}
     1624
     1625// we are supplied U, not Ut; w represents a diagonal matrix (also, we apply 1/w instead of w)
     1626psImage *p_pmSubSolve_wUt (psVector *w, psImage *U) {
     1627
     1628    psAssert (w->n == U->numCols, "w and U dimensions do not match");
     1629
     1630    // wUt has dimensions transposed relative to Ut.
     1631    psImage *wUt = psImageAlloc (U->numRows, U->numCols, PS_TYPE_F64);
     1632    psImageInit (wUt, 0.0);
     1633
     1634    for (int i = 0; i < wUt->numCols; i++) {
     1635        for (int j = 0; j < wUt->numRows; j++) {
     1636            if (!isfinite(w->data.F64[j])) continue;
     1637            if (w->data.F64[j] == 0.0) continue;
     1638            wUt->data.F64[j][i] = U->data.F64[i][j] / w->data.F64[j];
     1639        }
     1640    }
     1641    return wUt;
     1642}
     1643
     1644// XXX this is just standard matrix multiplication: use psMatrixMultiply?
     1645psImage *p_pmSubSolve_VwUt (psImage *V, psImage *wUt) {
     1646
     1647    psAssert (V->numCols == wUt->numRows, "matrix dimensions do not match");
     1648
     1649    psImage *Ainv = psImageAlloc (wUt->numCols, V->numRows, PS_TYPE_F64);
     1650
     1651    for (int i = 0; i < Ainv->numCols; i++) {
     1652        for (int j = 0; j < Ainv->numRows; j++) {
     1653            double sum = 0.0;
     1654            for (int k = 0; k < V->numCols; k++) {
     1655                sum += V->data.F64[j][k] * wUt->data.F64[k][i];
     1656            }
     1657            Ainv->data.F64[j][i] = sum;
     1658        }
     1659    }
     1660    return Ainv;
     1661}
     1662
     1663// we are supplied U, not Ut
     1664psVector *p_pmSubSolve_UtB (psImage *U, psVector *B) {
     1665
     1666    psAssert (U->numRows == B->n, "U and B dimensions do not match");
     1667
     1668    psVector *UtB = psVectorAlloc (U->numCols, PS_TYPE_F64);
     1669
     1670    for (int i = 0; i < U->numCols; i++) {
     1671        double sum = 0.0;
     1672        for (int j = 0; j < U->numRows; j++) {
     1673            sum += B->data.F64[j] * U->data.F64[j][i];
     1674        }
     1675        UtB->data.F64[i] = sum;
     1676    }
     1677    return UtB;
     1678}
     1679
     1680// w is diagonal
     1681psVector *p_pmSubSolve_wUtB (psVector *w, psVector *UtB) {
     1682
     1683    psAssert (w->n == UtB->n, "w and UtB dimensions do not match");
     1684
     1685    // wUt has dimensions transposed relative to Ut.
     1686    psVector *wUtB = psVectorAlloc (w->n, PS_TYPE_F64);
     1687    psVectorInit (wUtB, 0.0);
     1688
     1689    for (int i = 0; i < w->n; i++) {
     1690        if (!isfinite(w->data.F64[i])) continue;
     1691        if (w->data.F64[i] == 0.0) continue;
     1692        wUtB->data.F64[i] = UtB->data.F64[i] / w->data.F64[i];
     1693    }
     1694    return wUtB;
     1695}
     1696
     1697// this is basically matrix * vector
     1698psVector *p_pmSubSolve_VwUtB (psImage *V, psVector *wUtB) {
     1699
     1700    psAssert (V->numCols == wUtB->n, "V and wUtB dimensions do not match");
     1701
     1702    psVector *VwUtB = psVectorAlloc (V->numRows, PS_TYPE_F64);
     1703
     1704    for (int j = 0; j < V->numRows; j++) {
     1705        double sum = 0.0;
     1706        for (int i = 0; i < V->numCols; i++) {
     1707            sum += V->data.F64[j][i] * wUtB->data.F64[i];
     1708        }
     1709        VwUtB->data.F64[j] = sum;
     1710    }
     1711    return VwUtB;
     1712}
     1713
     1714// this is basically matrix * vector
     1715psVector *p_pmSubSolve_Ax (psImage *A, psVector *x) {
     1716
     1717    psAssert (A->numCols == x->n, "A and x dimensions do not match");
     1718
     1719    psVector *B = psVectorAlloc (A->numRows, PS_TYPE_F64);
     1720
     1721    for (int j = 0; j < A->numRows; j++) {
     1722        double sum = 0.0;
     1723        for (int i = 0; i < A->numCols; i++) {
     1724            sum += A->data.F64[j][i] * x->data.F64[i];
     1725        }
     1726        B->data.F64[j] = sum;
     1727    }
     1728    return B;
     1729}
     1730
     1731// this is basically Vector * vector
     1732double p_pmSubSolve_VdV (psVector *x, psVector *y) {
     1733
     1734    psAssert (x->n == y->n, "x and y dimensions do not match");
     1735
     1736    double sum = 0.0;
     1737    for (int i = 0; i < x->n; i++) {
     1738        sum += x->data.F64[i] * y->data.F64[i];
     1739    }
     1740    return sum;
     1741}
     1742
     1743double p_pmSubSolve_y2 (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) {
     1744
     1745    int footprint = stamps->footprint; // Half-size of stamps
     1746
     1747    double sum = 0.0;
     1748    for (int i = 0; i < stamps->num; i++) {
     1749
     1750        pmSubtractionStamp *stamp = stamps->stamps->data[i];
     1751        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
     1752
     1753        psKernel *weight = NULL;
     1754        psKernel *window = NULL;
     1755        psKernel *input = NULL;
     1756
     1757#ifdef USE_WEIGHT
     1758        weight = stamp->weight;
     1759#endif
     1760#ifdef USE_WINDOW
     1761        window = stamps->window;
     1762#endif
     1763
     1764        switch (kernels->mode) {
     1765            // MODE_1 : convolve image 1 to match image 2 (and vice versa)
     1766          case PM_SUBTRACTION_MODE_1:
     1767            input = stamp->image2;
     1768            break;
     1769          case PM_SUBTRACTION_MODE_2:
     1770            input = stamp->image1;
     1771            break;
     1772          default:
     1773            psAbort ("programming error");
     1774        }
     1775
     1776        for (int y = - footprint; y <= footprint; y++) {
     1777            for (int x = - footprint; x <= footprint; x++) {
     1778                double in = input->kernel[y][x];
     1779                double value = in*in;
     1780                if (weight) {
     1781                    float wtVal = weight->kernel[y][x];
     1782                    value *= wtVal;
     1783                }
     1784                if (window) {
     1785                    float  winVal = window->kernel[y][x];
     1786                    value *= winVal;
     1787                }
     1788                sum += value;
     1789            }
     1790        }
     1791    }
     1792    return sum;
     1793}
     1794
     1795double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) {
     1796
     1797    int footprint = stamps->footprint; // Half-size of stamps
     1798    int numKernels = kernels->num;      // Number of kernels
     1799
     1800    double sum = 0.0;
     1801
     1802    psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
     1803    psImageInit(residual->image, 0.0);
     1804
     1805    psImage *polyValues = NULL;         // Polynomial values
     1806
     1807    for (int i = 0; i < stamps->num; i++) {
     1808
     1809        pmSubtractionStamp *stamp = stamps->stamps->data[i];
     1810        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
     1811
     1812        psKernel *weight = NULL;
     1813        psKernel *window = NULL;
     1814        psKernel *target = NULL;
     1815        psKernel *source = NULL;
     1816
     1817        psArray *convolutions = NULL;
     1818
     1819#ifdef USE_WEIGHT
     1820        weight = stamp->weight;
     1821#endif
     1822#ifdef USE_WINDOW
     1823        window = stamps->window;
     1824#endif
     1825
     1826        switch (kernels->mode) {
     1827            // MODE_1 : convolve image 1 to match image 2 (and vice versa)
     1828          case PM_SUBTRACTION_MODE_1:
     1829            target = stamp->image2;
     1830            source = stamp->image1;
     1831            convolutions = stamp->convolutions1;
     1832            break;
     1833          case PM_SUBTRACTION_MODE_2:
     1834            target = stamp->image1;
     1835            source = stamp->image2;
     1836            convolutions = stamp->convolutions2;
     1837            break;
     1838          default:
     1839            psAbort ("programming error");
     1840        }
     1841
     1842        // Calculate coefficients of the kernel basis functions
     1843        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
     1844        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
     1845        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
     1846
     1847        psImageInit(residual->image, 0.0);
     1848        for (int j = 0; j < numKernels; j++) {
     1849            psKernel *convolution = convolutions->data[j]; // Convolution
     1850            double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient
     1851            for (int y = - footprint; y <= footprint; y++) {
     1852                for (int x = - footprint; x <= footprint; x++) {
     1853                    residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient;
     1854                }
     1855            }
     1856        }
     1857
     1858        for (int y = - footprint; y <= footprint; y++) {
     1859            for (int x = - footprint; x <= footprint; x++) {
     1860                double resid = target->kernel[y][x] - background - source->kernel[y][x] * norm + residual->kernel[y][x];
     1861                double value = PS_SQR(resid);
     1862                if (weight) {
     1863                    float wtVal = weight->kernel[y][x];
     1864                    value *= wtVal;
     1865                }
     1866                if (window) {
     1867                    float  winVal = window->kernel[y][x];
     1868                    value *= winVal;
     1869                }
     1870                sum += value;
     1871            }
     1872        }
     1873    }
     1874    psFree (polyValues);
     1875    psFree (residual);
     1876
     1877    return sum;
     1878}
     1879
     1880// I get confused by the index values between the image vs matrix usage:  In terms
     1881// of the elements of an image A(x,y) = A->data.F64[y][x] = A_x,y, a matrix
     1882// multiplication is: A_k,j * B_i,k = C_i,j
     1883
     1884
     1885bool psFitsWriteImageSimple (char *filename, psImage *image, psMetadata *header) {
     1886
     1887    psFits *fits = psFitsOpen(filename, "w");
     1888    psFitsWriteImage(fits, header, image, 0, NULL);
     1889    psFitsClose(fits);
     1890   
     1891    return true;
     1892}
     1893
     1894bool psVectorWriteFile (char *filename, const psVector *vector) {
     1895
     1896    FILE *f = fopen (filename, "w");
     1897    int fd = fileno(f);
     1898    p_psVectorPrint (fd, vector, "unnamed");
     1899    fclose (f);
     1900
     1901    return true;
     1902}
     1903

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

    r26342 r26373  
    646646    psVector *flux = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32);
    647647
     648    double maxValue = 0.0;
     649
    648650    // generate the window pixels
    649651    for (int y = -size; y <= size; y++) {
    … …  
    666668            }
    667669            stamps->window->kernel[y][x] = stats->robustMedian;
     670            if (maxValue < stats->robustMedian) {
     671                maxValue = stats->robustMedian;
     672            }
     673        }
     674    }
     675
     676    // re-normalize so chisquare values are sensible
     677    for (int y = -size; y <= size; y++) {
     678        for (int x = -size; x <= size; x++) {
     679            stamps->window->kernel[y][x] /= maxValue;
    668680        }
    669681    }

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

    r26341 r26373  
    308308    // if (!pmVisualIsVisual()) return true;
    309309
     310    double sum;
     311
    310312    int NXoff = index % NX;
    311313    int NYoff = index / NX;
    … …  
    315317
    316318    // insert the (target) kernel into the (target) image:
     319    sum = 0.0;
    317320    for (int y = -footprint; y <= footprint; y++) {
    318321        for (int x = -footprint; x <= footprint; x++) {
    319322            targetImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x];
    320         }
    321     }
     323            sum += targetImage->data.F32[y + NYpix][x + NXpix];
     324        }
     325    }
     326    targetImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
    322327
    323328    // insert the (source) kernel into the (source) image:
     329    sum = 0.0;
    324330    for (int y = -footprint; y <= footprint; y++) {
    325331        for (int x = -footprint; x <= footprint; x++) {
    326332            sourceImage->data.F32[y + NYpix][x + NXpix] = source->kernel[y][x];
    327         }
    328     }
    329    
     333            sum += sourceImage->data.F32[y + NYpix][x + NXpix];
     334        }
     335    }
     336    sourceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
     337
    330338    // insert the (convolution) kernel into the (convolution) image:
     339    sum = 0.0;
    331340    for (int y = -footprint; y <= footprint; y++) {
    332341        for (int x = -footprint; x <= footprint; x++) {
    333342            convolutionImage->data.F32[y + NYpix][x + NXpix] = -convolution->kernel[y][x];
    334         }
    335     }
     343            sum += convolutionImage->data.F32[y + NYpix][x + NXpix];
     344        }
     345    }
     346    convolutionImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
    336347   
    337348    // insert the (difference) kernel into the (difference) image:
     349    sum = 0.0;
    338350    for (int y = -footprint; y <= footprint; y++) {
    339351        for (int x = -footprint; x <= footprint; x++) {
    340352            differenceImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm;
    341         }
    342     }
     353            sum += differenceImage->data.F32[y + NYpix][x + NXpix];
     354        }
     355    }
     356    differenceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
    343357
    344358    // insert the (residual) kernel into the (residual) image:
     359    sum = 0.0;
    345360    for (int y = -footprint; y <= footprint; y++) {
    346361        for (int x = -footprint; x <= footprint; x++) {
    347362            residualImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm + convolution->kernel[y][x];
    348         }
    349     }
     363            sum += residualImage->data.F32[y + NYpix][x + NXpix];
     364        }
     365    }
     366    residualImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
    350367
    351368    // insert the (fresidual) kernel into the (fresidual) image:
    352369    for (int y = -footprint; y <= footprint; y++) {
    353370        for (int x = -footprint; x <= footprint; x++) {
    354             fresidualImage->data.F32[y + NYpix][x + NXpix] = PS_SQR(residualImage->data.F32[y + NYpix][x + NXpix]) / (target->kernel[y][x] + 25.0);
     371            fresidualImage->data.F32[y + NYpix][x + NXpix] = residualImage->data.F32[y + NYpix][x + NXpix] / sqrt(PS_MAX(target->kernel[y][x], 100.0));
    355372        }
    356373    }
Note: See TracChangeset for help on using the changeset viewer.