IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26430


Ignore:
Timestamp:
Dec 15, 2009, 6:40:02 PM (16 years ago)
Author:
eugene
Message:

attempt to handle degeneracy with svn (minimal success)

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r26373 r26430  
    801801        }
    802802
     803        if ((stamp->x <= 0.0) && (stamp->y <= 0.0)) {
     804            psAbort ("bad stamp");
     805        }
     806        if (!isfinite(stamp->x) && !isfinite(stamp->y)) {
     807            psAbort ("bad stamp");
     808        }
     809
    803810        if (pmSubtractionThreaded()) {
    804811            psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_CALCULATE_EQUATION");
    … …  
    822829    }
    823830
    824     pmSubtractionVisualPlotLeastSquares(stamps);
     831    // pmSubtractionVisualPlotLeastSquares(stamps);
    825832
    826833    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate equation: %f sec",
    … …  
    838845psImage *p_pmSubSolve_VwUt (psImage *V, psImage *wUt);
    839846
    840 psVector *p_pmSubSolve_UtB (psImage *U, psVector *B);
    841 psVector *p_pmSubSolve_wUtB (psVector *w, psVector *UtB);
    842 psVector *p_pmSubSolve_VwUtB (psImage *V, psVector *wUtB);
    843 
    844 psVector *p_pmSubSolve_Ax (psImage *A, psVector *x);
    845 
    846 double p_pmSubSolve_y2 (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps);
    847 double p_pmSubSolve_VdV (psVector *x, psVector *y);
     847bool p_pmSubSolve_SetWeights (psVector *wApply, psVector *w, psVector *wMask);
     848
     849bool p_pmSubSolve_UtB (psVector **UtB, psImage *U, psVector *B);
     850bool p_pmSubSolve_wUtB (psVector **wUtB, psVector *w, psVector *UtB);
     851bool p_pmSubSolve_VwUtB (psVector **VwUtB, psImage *V, psVector *wUtB);
     852
     853bool p_pmSubSolve_Ax (psVector **B, psImage *A, psVector *x);
     854bool p_pmSubSolve_VdV (double *value, psVector *x, psVector *y);
     855bool p_pmSubSolve_y2 (double *y2, pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps);
     856
     857psImage *p_pmSubSolve_Xvar (psImage *V, psVector *w);
     858
    848859double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps);
    849860
    … …  
    976987# define SVD_ANALYSIS 1
    977988# define COEFF_SIG 0.0
    978 # define SVD_TOL 1e-9
     989# define SVD_TOL 0.0
    979990
    980991        psVector *solution = NULL;
    981992        psVector *solutionErr = NULL;
     993
     994#ifdef TESTING
     995        psFitsWriteImageSimple("A.fits", sumMatrix, NULL);
     996        psVectorWriteFile ("B.dat", sumVector);
     997#endif
    982998
    983999        // Use SVD to determine the kernel coeffs (and validate)
    … …  
    10031019            // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0
    10041020
    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 
    10201021            // 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
     1022            // kernel and normalization terms from the background term.
     1023            // XXX is this true?  or was this due to an error in the analysis?
     1024
    10241025            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;
    10301026
    10311027            // now pull out the kernel elements into their own square matrix
    … …  
    10331029            psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64);
    10341030
    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)
    10371031            for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) {
    1038                 // if (ix == normIndex) continue;
    10391032                if (ix == bgIndex) continue;
    10401033                for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) {
    1041                     // if (iy == normIndex) continue;
    10421034                    if (iy == bgIndex) continue;
    10431035                    kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix];
    10441036                    ky++;
    10451037                }
    1046                 // kernelVector->data.F64[kx] = sumVector->data.F64[ix] - norm*sumMatrix->data.F64[normIndex][ix];
    10471038                kernelVector->data.F64[kx] = sumVector->data.F64[ix];
    10481039                kx++;
    10491040            }
    1050 
    1051             psFitsWriteImageSimple("A.fits", kernelMatrix, NULL);
    10521041
    10531042            psImage *U = NULL;
    … …  
    10591048            }
    10601049
    1061             // TEST psFitsWriteImageSimple("U.fits", U, NULL);
    1062             // TEST psFitsWriteImageSimple("V.fits", V, NULL);
    1063             // TEST psVectorWriteFile("w.dat", w);
    1064 
    10651050            // calculate A_inverse:
    10661051            // Ainv = V * w * U^T
    10671052            psImage *wUt  = p_pmSubSolve_wUt (w, U);
    10681053            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 
     1054            psImage *Xvar = NULL;
     1055            psFree (wUt);
     1056
     1057# ifdef TESTING
    10741058            // kernel terms:
    10751059            for (int i = 0; i < w->n; i++) {
    10761060                fprintf (stderr, "w: %f\n", w->data.F64[i]);
    10771061            }
    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);
     1062# endif
     1063            // loop over w adding in more and more of the values until chisquare is no longer
     1064            // dropping significantly.
     1065            // XXX this does not seem to work very well: we seem to need all terms even for
     1066            // simple cases...
     1067
     1068            psVector *Xsvd = NULL;
     1069            {
     1070                psVector *Ax = NULL;
     1071                psVector *UtB = NULL;
     1072                psVector *wUtB = NULL;
     1073
     1074                psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64);
     1075                psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8);
     1076                psVectorInit (wMask, 1); // start by masking everything
     1077
     1078                double chiSquareLast = NAN;
     1079                int maxWeight = 0;
     1080
     1081                double Axx, Bx, y2;
     1082
     1083                // XXX this is an attempt to exclude insignificant modes.
     1084                // it was not successful with the ISIS kernel set: removing even
     1085                // the least significant mode leaves additional ringing / noise
     1086                // because the terms are so coupled.
     1087                for (int k = 0; false && (k < w->n); k++) {
     1088
     1089                    // unmask the k-th weight
     1090                    wMask->data.U8[k] = 0;
     1091                    p_pmSubSolve_SetWeights(wApply, w, wMask);
     1092
     1093                    // solve for x:
     1094                    // x = V * w * (U^T * B)
     1095                    p_pmSubSolve_UtB (&UtB, U, kernelVector);
     1096                    p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
     1097                    p_pmSubSolve_VwUtB (&Xsvd, 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                    p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
     1103                    p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
     1104                    p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
     1105                    p_pmSubSolve_y2 (&y2, kernels, stamps);
     1106
     1107                    // apparently, this works (compare with the brute force value below
     1108                    double chiSquare = Axx - 2.0*Bx + y2;
     1109                    double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare;
     1110                    chiSquareLast = chiSquare;
     1111
     1112                    // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi);
     1113                    if (k && !maxWeight && (deltaChi < 1.0)) {
     1114                        maxWeight = k;
     1115                    }
     1116                }
     1117
     1118                // keep all terms or we get extra ringing
     1119                maxWeight = w->n;
     1120                psVectorInit (wMask, 1);
     1121                for (int k = 0; k < maxWeight; k++) {
     1122                    wMask->data.U8[k] = 0;
     1123                }
     1124                p_pmSubSolve_SetWeights(wApply, w, wMask);
     1125
     1126                // solve for x:
     1127                // x = V * w * (U^T * B)
     1128                p_pmSubSolve_UtB (&UtB, U, kernelVector);
     1129                p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
     1130                p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
     1131
     1132                // chi-square for this system of equations:
     1133                // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
     1134                // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
     1135                p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
     1136                p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
     1137                p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
     1138                p_pmSubSolve_y2 (&y2, kernels, stamps);
     1139
     1140                // apparently, this works (compare with the brute force value below
     1141                double chiSquare = Axx - 2.0*Bx + y2;
     1142                psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare);
     1143
     1144                // re-calculate A^{-1} to get new variances:
     1145                // Ainv = V * w * U^T
     1146                // XXX since we keep all terms, this is identical to Ainv
     1147                psImage *wUt  = p_pmSubSolve_wUt (wApply, U);
     1148                Xvar = p_pmSubSolve_VwUt (V, wUt);
     1149                psFree (wUt);
     1150
     1151                psFree (Ax);
     1152                psFree (UtB);
     1153                psFree (wUtB);
     1154                psFree (wApply);
     1155                psFree (wMask);
     1156            }
    11101157
    11111158            // copy the kernel solutions to the full solution vector:
    … …  
    11141161
    11151162            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 }
    11221163                if (ix == bgIndex) {
    11231164                    solution->data.F64[ix] = 0;
    … …  
    11301171            }
    11311172
    1132             // TEST psVectorWriteFile("X.dat", Xsvd);
    1133             // TEST psVector *Bo = p_pmSubSolve_Ax (sumMatrix, Xsvd);
    1134             // TEST psVectorWriteFile("Bo.dat", Bo);
    1135                
    11361173            psFree (kernelMatrix);
    11371174            psFree (kernelVector);
    … …  
    11411178            psFree (w);
    11421179
    1143             psFree (wUt);
    11441180            psFree (Ainv);
    1145 
    1146             psFree (UtB);
    1147             psFree (wUtB);
    11481181            psFree (Xsvd);
    11491182        } else {
    11501183            psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
    11511184            psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
    1152             psFree(sumMatrix);
    11531185            if (!luMatrix) {
    11541186                psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n");
     1187                psFree(solution);
    11551188                psFree(sumVector);
     1189                psFree(sumMatrix);
    11561190                psFree(luMatrix);
    11571191                psFree(permutation);
    … …  
    11601194
    11611195            solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation);
    1162             psFree(sumVector);
    11631196            psFree(luMatrix);
    11641197            psFree(permutation);
    11651198            if (!solution) {
    11661199                psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
     1200                psFree(solution);
     1201                psFree(sumVector);
     1202                psFree(sumMatrix);
    11671203                return NULL;
     1204            }
     1205
     1206            // XXX LUD does not provide A^{-1}?  fake the error for now
     1207            solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     1208            for (int ix = 0; ix < sumVector->n; ix++) {
     1209                solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix];
    11681210            }
    11691211        }
    … …  
    11881230            int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
    11891231            for (int i = 0; i < numKernels * numPoly; i++) {
    1190                 fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i]));
     1232                // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i]));
    11911233                if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) {
    1192                     fprintf (stderr, "drop\n");
     1234                    // XXX fprintf (stderr, "drop\n");
    11931235                    kernels->solution1->data.F64[i] = 0.0;
    11941236                } else {
    1195                     fprintf (stderr, "keep\n");
     1237                    // XXX fprintf (stderr, "keep\n");
    11961238                    kernels->solution1->data.F64[i] = solution->data.F64[i];
    11971239                }
    … …  
    12011243        // fprintf (stderr, "chi square Brute = %f\n", chiSquare);
    12021244
    1203         psFree (solution);
     1245        psFree(solution);
     1246        psFree(sumVector);
     1247        psFree(sumMatrix);
    12041248
    12051249#ifdef TESTING
    … …  
    16621706
    16631707// we are supplied U, not Ut
    1664 psVector *p_pmSubSolve_UtB (psImage *U, psVector *B) {
     1708bool p_pmSubSolve_UtB (psVector **UtB, psImage *U, psVector *B) {
    16651709
    16661710    psAssert (U->numRows == B->n, "U and B dimensions do not match");
    16671711
    1668     psVector *UtB = psVectorAlloc (U->numCols, PS_TYPE_F64);
     1712    UtB[0] = psVectorRecycle (UtB[0], U->numCols, PS_TYPE_F64);
    16691713
    16701714    for (int i = 0; i < U->numCols; i++) {
    … …  
    16731717            sum += B->data.F64[j] * U->data.F64[j][i];
    16741718        }
    1675         UtB->data.F64[i] = sum;
    1676     }
    1677     return UtB;
     1719        UtB[0]->data.F64[i] = sum;
     1720    }
     1721    return true;
    16781722}
    16791723
    16801724// w is diagonal
    1681 psVector *p_pmSubSolve_wUtB (psVector *w, psVector *UtB) {
     1725bool p_pmSubSolve_wUtB (psVector **wUtB, psVector *w, psVector *UtB) {
    16821726
    16831727    psAssert (w->n == UtB->n, "w and UtB dimensions do not match");
    16841728
    16851729    // wUt has dimensions transposed relative to Ut.
    1686     psVector *wUtB = psVectorAlloc (w->n, PS_TYPE_F64);
    1687     psVectorInit (wUtB, 0.0);
     1730    wUtB[0] = psVectorRecycle (wUtB[0], w->n, PS_TYPE_F64);
     1731    psVectorInit (wUtB[0], 0.0);
    16881732
    16891733    for (int i = 0; i < w->n; i++) {
    16901734        if (!isfinite(w->data.F64[i])) continue;
    16911735        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;
     1736        wUtB[0]->data.F64[i] = UtB->data.F64[i] / w->data.F64[i];
     1737    }
     1738    return true;
    16951739}
    16961740
    16971741// this is basically matrix * vector
    1698 psVector *p_pmSubSolve_VwUtB (psImage *V, psVector *wUtB) {
     1742bool p_pmSubSolve_VwUtB (psVector **VwUtB, psImage *V, psVector *wUtB) {
    16991743
    17001744    psAssert (V->numCols == wUtB->n, "V and wUtB dimensions do not match");
    17011745
    1702     psVector *VwUtB = psVectorAlloc (V->numRows, PS_TYPE_F64);
     1746    VwUtB[0] = psVectorRecycle (*VwUtB, V->numRows, PS_TYPE_F64);
    17031747
    17041748    for (int j = 0; j < V->numRows; j++) {
    … …  
    17071751            sum += V->data.F64[j][i] * wUtB->data.F64[i];
    17081752        }
    1709         VwUtB->data.F64[j] = sum;
    1710     }
    1711     return VwUtB;
     1753        VwUtB[0]->data.F64[j] = sum;
     1754    }
     1755    return true;
    17121756}
    17131757
    17141758// this is basically matrix * vector
    1715 psVector *p_pmSubSolve_Ax (psImage *A, psVector *x) {
     1759bool p_pmSubSolve_Ax (psVector **B, psImage *A, psVector *x) {
    17161760
    17171761    psAssert (A->numCols == x->n, "A and x dimensions do not match");
    17181762
    1719     psVector *B = psVectorAlloc (A->numRows, PS_TYPE_F64);
     1763    B[0] = psVectorRecycle (*B, A->numRows, PS_TYPE_F64);
    17201764
    17211765    for (int j = 0; j < A->numRows; j++) {
    … …  
    17241768            sum += A->data.F64[j][i] * x->data.F64[i];
    17251769        }
    1726         B->data.F64[j] = sum;
    1727     }
    1728     return B;
     1770        B[0]->data.F64[j] = sum;
     1771    }
     1772    return true;
    17291773}
    17301774
    17311775// this is basically Vector * vector
    1732 double p_pmSubSolve_VdV (psVector *x, psVector *y) {
     1776bool p_pmSubSolve_VdV (double *value, psVector *x, psVector *y) {
    17331777
    17341778    psAssert (x->n == y->n, "x and y dimensions do not match");
    … …  
    17381782        sum += x->data.F64[i] * y->data.F64[i];
    17391783    }
    1740     return sum;
    1741 }
    1742 
    1743 double p_pmSubSolve_y2 (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) {
     1784    *value = sum;
     1785    return true;
     1786}
     1787
     1788bool p_pmSubSolve_y2 (double *y2, pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) {
    17441789
    17451790    int footprint = stamps->footprint; // Half-size of stamps
    … …  
    17901835        }
    17911836    }
    1792     return sum;
     1837    *y2 = sum;
     1838    return true;
    17931839}
    17941840
    … …  
    18781924}
    18791925
     1926bool p_pmSubSolve_SetWeights (psVector *wApply, psVector *w, psVector *wMask) {
     1927
     1928    for (int i = 0; i < w->n; i++) {
     1929        wApply->data.F64[i] = wMask->data.U8[i] ? 0.0 : w->data.F64[i];
     1930    }
     1931    return true;
     1932}
     1933
     1934// we are supplied V and w; w represents a diagonal matrix (also, we apply 1/w instead of w)
     1935psImage *p_pmSubSolve_Xvar (psImage *V, psVector *w) {
     1936
     1937    psAssert (w->n == V->numCols, "w and U dimensions do not match");
     1938
     1939    psImage *Vn = psImageAlloc (V->numCols, V->numRows, PS_TYPE_F64);
     1940    psImageInit (Vn, 0.0);
     1941
     1942    // generate Vn = V * w^{-1}
     1943    for (int j = 0; j < Vn->numRows; j++) {
     1944        for (int i = 0; i < Vn->numCols; i++) {
     1945            if (!isfinite(w->data.F64[i])) continue;
     1946            if (w->data.F64[i] == 0.0) continue;
     1947            Vn->data.F64[j][i] = V->data.F64[j][i] / w->data.F64[i];
     1948        }
     1949    }
     1950
     1951    psImage *Xvar = psImageAlloc (V->numCols, V->numRows, PS_TYPE_F64);
     1952    psImageInit (Xvar, 0.0);
     1953
     1954    // generate Xvar = Vn * Vn^T
     1955    for (int j = 0; j < Vn->numRows; j++) {
     1956        for (int i = 0; i < Vn->numCols; i++) {
     1957            double sum = 0.0;
     1958            for (int k = 0; k < Vn->numCols; k++) {
     1959                sum += Vn->data.F64[k][i]*Vn->data.F64[k][j];
     1960            }
     1961            Xvar->data.F64[j][i] = sum;
     1962        }
     1963    }
     1964    return Xvar;
     1965}
     1966
    18801967// I get confused by the index values between the image vs matrix usage:  In terms
    18811968// of the elements of an image A(x,y) = A->data.F64[y][x] = A_x,y, a matrix
Note: See TracChangeset for help on using the changeset viewer.