IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 25895


Ignore:
Timestamp:
Oct 19, 2009, 11:56:07 AM (17 years ago)
Author:
Paul Price
Message:

For dual convolution, throw out low-significance solutions and the
non-dominant solution when compared across the two kernels. This
seems to work well. It's not yet clear how to do this when spatial
variation is permitted.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/pap/psModules/src/imcombine/pmSubtractionEquation.c

    r25869 r25895  
    1515#include "pmSubtractionVisual.h"
    1616
    17 #define TESTING                         // TESTING output for debugging; may not work with threads!
     17//#define TESTING                         // TESTING output for debugging; may not work with threads!
    1818
    1919#define USE_VARIANCE                    // Include variance in equation?
    … …  
    377377}
    378378
     379// Merge dual matrices and vectors into single matrix equation
     380// Have: Aa = Ct.b + d
     381// Have: Ca = Bb + e
     382// Set: F = ( A -Ct ;  C -B )
     383// Set: g = ( a ; b )
     384// Set: h = ( d ; e )
     385// So that we combine the above two equations: Fg = h
     386static bool calculateEquationDual(psImage **outMatrix,
     387                                  psVector **outVector,
     388                                  const psImage *sumMatrix1,
     389                                  const psImage *sumMatrix2,
     390                                  const psImage *sumMatrixX,
     391                                  const psVector *sumVector1,
     392                                  const psVector *sumVector2
     393                                  )
     394{
     395    psAssert(sumMatrix1 && sumMatrix2 && sumMatrixX, "Require input matrices");
     396    psAssert(sumVector1 && sumVector2, "Require input vectors");
     397    int num1 = sumVector1->n;   // Number of parameters in first set
     398    int num2 = sumVector2->n;   // Number of parameters in second set
     399    int num = num1 + num2;      // Number of parameters in new set
     400
     401    psAssert(sumMatrix1->type.type == PS_TYPE_F64 &&
     402             sumMatrix2->type.type == PS_TYPE_F64 &&
     403             sumMatrixX->type.type == PS_TYPE_F64 &&
     404             sumVector1->type.type == PS_TYPE_F64 &&
     405             sumVector2->type.type == PS_TYPE_F64,
     406             "Require input type is F64");
     407
     408    psAssert(outMatrix, "Require output matrix");
     409    psAssert(outVector, "Require output vector");
     410    if (!*outMatrix) {
     411        *outMatrix = psImageAlloc(num, num, PS_TYPE_F64);
     412    }
     413    if (!*outVector) {
     414        *outVector = psVectorAlloc(num, PS_TYPE_F64);
     415    }
     416    psImage *matrix = *outMatrix;
     417    psVector *vector = *outVector;
     418
     419    psAssert(sumMatrix1->numCols == num1 && sumMatrix1->numRows == num1, "Require size NxN");
     420    psAssert(sumMatrix2->numCols == num2 && sumMatrix2->numRows == num2, "Require size MxM");
     421    psAssert(sumMatrixX->numCols == num1 && sumMatrixX->numRows == num2, "Require size MxN");
     422
     423    memcpy(vector->data.F64, sumVector1->data.F64, num1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
     424    memcpy(&vector->data.F64[num1], sumVector2->data.F64, num2 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
     425
     426    for (int i = 0; i < num1; i++) {
     427        memcpy(matrix->data.F64[i], sumMatrix1->data.F64[i], num1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
     428        for (int j = 0, k = num1; j < num2; j++, k++) {
     429            matrix->data.F64[i][k] = - sumMatrixX->data.F64[j][i];
     430        }
     431    }
     432    for (int i1 = 0, i2 = num1; i1 < num2; i1++, i2++) {
     433        memcpy(matrix->data.F64[i2], sumMatrixX->data.F64[i1], num1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
     434        for (int j = 0, k = num1; j < num2; j++, k++) {
     435            matrix->data.F64[i2][k] = - sumMatrix2->data.F64[i1][j];
     436        }
     437    }
     438
     439    return true;
     440}
     441
    379442
    380443// Add in penalty term to least-squares vector
    … …  
    881944        calculatePenalty(sumMatrix2, sumVector2, kernels, -sumMatrix1->data.F64[bgIndex][bgIndex]);
    882945
    883         // Pure matrix operations
    884         // A * a = Ct * b + d
    885         // C * a = B  * b + e
    886         //
    887         // Set F = ( A -Ct ;  C -B )
    888         // Set g = ( a ; b )
    889         // Set h = ( d ; e )
    890         // So that we combine the above two equations: Fg = h
    891 
    892         int num = numParams + numParams2; // Number of params for new set
    893         psImage *F = psImageAlloc(num, num, PS_TYPE_F64);
    894         psVector *h = psVectorAlloc(num, PS_TYPE_F64);
    895 
    896         for (int i = 0; i < numParams; i++) {
    897             h->data.F64[i] = sumVector1->data.F64[i];
    898             for (int j = 0; j < numParams; j++) {
    899                 F->data.F64[i][j] = sumMatrix1->data.F64[i][j];
    900             }
     946        psImage *sumMatrix = NULL;      // Combined matrix
     947        psVector *sumVector = NULL;     // Combined vector
     948        calculateEquationDual(&sumMatrix, &sumVector, sumMatrix1, sumMatrix2,
     949                              sumMatrixX, sumVector1, sumVector2);
     950
     951#ifdef TESTING
     952        {
     953            psFits *fits = psFitsOpen("sumMatrix.fits", "w");
     954            psFitsWriteImage(fits, NULL, sumMatrix, 0, NULL);
     955            psFitsClose(fits);
     956        }
     957        {
     958            psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64);
     959            psFits *fits = psFitsOpen("sumVector.fits", "w");
     960            for (int i = 0; i < numParams + numParams2; i++) {
     961                image->data.F64[0][i] = sumVector->data.F64[i];
     962            }
     963            psFitsWriteImage(fits, NULL, image, 0, NULL);
     964            psFree(image);
     965            psFitsClose(fits);
     966        }
     967#endif
     968
     969        psVector *solution = NULL;                       // Solution to equation!
     970        psVector *mask = psVectorAlloc(numParams + numParams2, PS_TYPE_U8); // Mask of parameters
     971        psVectorInit(mask, 0);
     972        {
     973            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector);
     974            if (!solution) {
     975                psError(PS_ERR_UNKNOWN, false, "SVD solution of least-squares equation failed.\n");
     976                psFree(sumMatrix);
     977                psFree(sumVector);
     978                psFree(mask);
     979                return NULL;
     980            }
     981
     982            // Remove a kernel basis for image 1 from the equation
     983#define MASK_BASIS_1(INDEX) \
     984            { \
     985                for (int k = 0; k < numParams2; k++) { \
     986                    sumMatrix1->data.F64[k][INDEX] = 0.0; \
     987                    sumMatrix1->data.F64[INDEX][k] = 0.0; \
     988                    sumMatrixX->data.F64[k][INDEX] = 0.0; \
     989                } \
     990                sumMatrix1->data.F64[bgIndex][INDEX] = 0.0; \
     991                sumMatrix1->data.F64[INDEX][bgIndex] = 0.0; \
     992                sumMatrix1->data.F64[normIndex][INDEX] = 0.0; \
     993                sumMatrix1->data.F64[INDEX][normIndex] = 0.0; \
     994                sumMatrix1->data.F64[INDEX][INDEX] = 1.0; \
     995                sumVector1->data.F64[INDEX] = 0.0; \
     996                mask->data.U8[INDEX] = 0xFF; \
     997            }
     998
     999            // Remove a kernel basis for image 2 from the equation
     1000#define MASK_BASIS_2(INDEX) \
     1001            { \
     1002                for (int k = 0; k < numParams2; k++) { \
     1003                    sumMatrix2->data.F64[k][j] = 0.0; \
     1004                    sumMatrix2->data.F64[j][k] = 0.0; \
     1005                    sumMatrixX->data.F64[j][k] = 0.0; \
     1006                } \
     1007                sumMatrix2->data.F64[INDEX][INDEX] = 1.0; \
     1008                sumMatrixX->data.F64[j][normIndex] = 0.0; \
     1009                sumMatrixX->data.F64[j][bgIndex] = 0.0; \
     1010                sumVector2->data.F64[j] = 0.0; \
     1011                mask->data.U8[numParams + j] = 0xFF; \
     1012            }
     1013
     1014            #define TOL 1.0e-5
     1015            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     1016            double norm = solution->data.F64[normIndex];        // Normalisation
     1017            double thresh = norm * TOL;                         // Threshold for low parameters
    9011018            for (int j = 0; j < numParams2; j++) {
    902                 F->data.F64[i][numParams + j] = - sumMatrixX->data.F64[j][i];
    903             }
    904         }
    905         for (int i = 0; i < numParams2; i++) {
    906             h->data.F64[numParams + i] = sumVector2->data.F64[i];
    907             for (int j = 0; j < numParams; j++) {
    908                 F->data.F64[numParams + i][j] = sumMatrixX->data.F64[i][j];
    909             }
    910             for (int j = 0; j < numParams2; j++) {
    911                 F->data.F64[numParams + i][numParams + j] = - sumMatrix2->data.F64[i][j];
    912             }
    913         }
     1019                double param1 = solution->data.F64[j],
     1020                    param2 = solution->data.F64[numParams + j]; // Parameters of interest
     1021                fprintf(stderr, "%lf %lf   ", param1, param2);
     1022                if (fabs(param1) < thresh) {
     1023                    fprintf(stderr, "Parameter %d: 1 below threshold\n", j);
     1024                    MASK_BASIS_1(j);
     1025                }
     1026                if (fabs(param2) < thresh) {
     1027                    fprintf(stderr, "Parameter %d: 2 below threshold\n", j);
     1028                    MASK_BASIS_2(j);
     1029                }
     1030
     1031                if (!mask->data.U8[j] && !mask->data.U8[numParams + j]) {
     1032                    if (fabs(param1) < fabs(param2)) {
     1033                        fprintf(stderr, "Parameter %d: 1 < 2\n", j);
     1034                        MASK_BASIS_1(j);
     1035                    } else {
     1036                        fprintf(stderr, "Parameter %d: 2 < 1\n", j);
     1037                        MASK_BASIS_2(j);
     1038                    }
     1039                }
     1040            }
     1041        }
     1042
     1043        calculateEquationDual(&sumMatrix, &sumVector, sumMatrix1, sumMatrix2,
     1044                              sumMatrixX, sumVector1, sumVector2);
     1045
     1046#ifdef TESTING
     1047        {
     1048            psFits *fits = psFitsOpen("sumMatrixFix.fits", "w");
     1049            psFitsWriteImage(fits, NULL, sumMatrix, 0, NULL);
     1050            psFitsClose(fits);
     1051        }
     1052        {
     1053            psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64);
     1054            psFits *fits = psFitsOpen("sumVectorFix.fits", "w");
     1055            for (int i = 0; i < numParams + numParams2; i++) {
     1056                image->data.F64[0][i] = sumVector->data.F64[i];
     1057            }
     1058            psFitsWriteImage(fits, NULL, image, 0, NULL);
     1059            psFree(image);
     1060            psFitsClose(fits);
     1061        }
     1062#endif
     1063
     1064        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector);
     1065        if (!solution) {
     1066            psError(PS_ERR_UNKNOWN, false, "SVD solution of least-squares equation failed.\n");
     1067            psFree(sumMatrix);
     1068            psFree(sumVector);
     1069            psFree(mask);
     1070            return NULL;
     1071        }
     1072
     1073#if 0
     1074        for (int i = 0; i < num; i++) {
     1075            if (mask->data.U8[i]) {
     1076                solution->data.F64[i] = 0.0;
     1077            }
     1078        }
     1079#endif
     1080        psFree(mask);
     1081
    9141082        psFree(sumMatrix1);
    9151083        psFree(sumMatrix2);
    … …  
    9181086        psFree(sumVector2);
    9191087
     1088        psFree(sumMatrix);
     1089        psFree(sumVector);
     1090
    9201091#ifdef TESTING
    9211092        {
    922             psFits *fits = psFitsOpen("sumMatrix.fits", "w");
    923             psFitsWriteImage(fits, NULL, F, 0, NULL);
    924             psFitsClose(fits);
    925         }
    926         {
    927             psImage *image = psImageAlloc(1, num, PS_TYPE_F64);
    928             psFits *fits = psFitsOpen("sumVector.fits", "w");
    929             for (int i = 0; i < num; i++) {
    930                 image->data.F64[0][i] = h->data.F64[i];
     1093            psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64);
     1094            psFits *fits = psFitsOpen("solnVector.fits", "w");
     1095            for (int i = 0; i < numParams + numParams2; i++) {
     1096                image->data.F64[0][i] = solution->data.F64[i];
    9311097            }
    9321098            psFitsWriteImage(fits, NULL, image, 0, NULL);
    … …  
    9361102#endif
    9371103
    938         psVector *g = NULL;             // Solution!
    939         psVector *mask = psVectorAlloc(num, PS_TYPE_U8); // Mask of parameters
    940         psVectorInit(mask, 0);
    941         {
    942             g = psMatrixSolveSVD(g, F, h);
    943             if (!g) {
    944                 psError(PS_ERR_UNKNOWN, false, "SVD solution of least-squares equation failed.\n");
    945                 psFree(F);
    946                 psFree(h);
    947                 psFree(mask);
    948                 return NULL;
    949             }
    950 
    951             #define TOL 1.0e-5
    952             int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
    953             double norm = g->data.F64[normIndex];               // Normalisation
    954             double thresh = norm * TOL;                         // Threshold for low parameters
    955             for (int j = 0; j < numParams2; j++) {
    956                 double param1 = g->data.F64[j], param2 = g->data.F64[numParams + j]; // Parameters of interest
    957                 fprintf(stderr, "%lf %lf   ", param1, param2);
    958                 if (fabs(param1) < thresh) {
    959                     fprintf(stderr, "Parameter %d 1 below threshold\n", j);
    960                     for (int k = 0; k < numParams2; k++) {
    961                         // Matrix 1
    962                         F->data.F64[k][j] = 0.0;
    963                         F->data.F64[j][k] = 0.0;
    964 #if 1
    965                         // Cross matrices
    966                         F->data.F64[numParams + j][k] = 0.0;
    967                         F->data.F64[j][numParams + k] = 0.0;
    968                         F->data.F64[k][numParams + j] = 0.0;
    969                         F->data.F64[numParams + k][j] = 0.0;
    970 #endif
    971                     }
    972                     F->data.F64[bgIndex][j] = 0.0;
    973                     F->data.F64[j][bgIndex] = 0.0;
    974                     F->data.F64[normIndex][j] = 0.0;
    975                     F->data.F64[j][normIndex] = 0.0;
    976                     F->data.F64[j][j] = 1.0;
    977                     h->data.F64[j] = 0.0;
    978                     mask->data.U8[j] = 0xFF;
    979                 }
    980                 if (fabs(param2) < thresh) {
    981                     fprintf(stderr, "Parameter %d 2 below threshold\n", j);
    982                     for (int k = 0; k < numParams2; k++) {
    983                         // Matrix 2
    984                         F->data.F64[k][numParams + j] = 0.0;
    985                         F->data.F64[numParams + j][k] = 0.0;
    986 #if 1
    987                         // Cross matrices
    988                         F->data.F64[numParams + j][k] = 0.0;
    989                         F->data.F64[j][numParams + k] = 0.0;
    990                         F->data.F64[k][numParams + j] = 0.0;
    991                         F->data.F64[numParams + k][j] = 0.0;
    992 #endif
    993                     }
    994                     F->data.F64[bgIndex][numParams + j] = 0.0;
    995                     F->data.F64[numParams + j][bgIndex] = 0.0;
    996                     F->data.F64[normIndex][numParams + j] = 0.0;
    997                     F->data.F64[numParams + j][normIndex] = 0.0;
    998                     F->data.F64[numParams + j][numParams + j] = 1.0;
    999                     h->data.F64[numParams + j] = 0.0;
    1000                     mask->data.U8[numParams + j] = 0xFF;
    1001                 }
    1002 
    1003                 if (!mask->data.U8[j] && !mask->data.U8[numParams + j]) {
    1004                     if (fabs(param1) < fabs(param2)) {
    1005                         fprintf(stderr, "Parameter %d 1 < 2\n", j);
    1006                         for (int k = 0; k < numParams2; k++) {
    1007                             // Matrix 1
    1008                             F->data.F64[k][j] = 0.0;
    1009                             F->data.F64[j][k] = 0.0;
    1010 #if 0
    1011                             // Cross matrices
    1012                             F->data.F64[numParams + j][k] = 0.0;
    1013                             F->data.F64[j][numParams + k] = 0.0;
    1014                             F->data.F64[k][numParams + j] = 0.0;
    1015                             F->data.F64[numParams + k][j] = 0.0;
    1016 #endif
    1017                         }
    1018                         F->data.F64[bgIndex][j] = 0.0;
    1019                         F->data.F64[j][bgIndex] = 0.0;
    1020                         F->data.F64[normIndex][j] = 0.0;
    1021                         F->data.F64[j][normIndex] = 0.0;
    1022                         F->data.F64[j][j] = 1.0;
    1023                         h->data.F64[j] = 0.0;
    1024                         mask->data.U8[j] = 0xFF;
    1025                     } else {
    1026                         fprintf(stderr, "Parameter %d 2 < 1\n", j);
    1027                         for (int k = 0; k < numParams2; k++) {
    1028                             // Matrix 1
    1029                             F->data.F64[k][numParams + j] = 0.0;
    1030                             F->data.F64[numParams + j][k] = 0.0;
    1031 #if 0
    1032                             // Cross matrices
    1033                             F->data.F64[numParams + j][k] = 0.0;
    1034                             F->data.F64[j][numParams + k] = 0.0;
    1035                             F->data.F64[k][numParams + j] = 0.0;
    1036                             F->data.F64[numParams + k][j] = 0.0;
    1037 #endif
    1038                         }
    1039                         F->data.F64[bgIndex][numParams + j] = 0.0;
    1040                         F->data.F64[numParams + j][bgIndex] = 0.0;
    1041                         F->data.F64[normIndex][numParams + j] = 0.0;
    1042                         F->data.F64[numParams + j][normIndex] = 0.0;
    1043                         F->data.F64[numParams + j][numParams + j] = 1.0;
    1044                         h->data.F64[numParams + j] = 0.0;
    1045                         mask->data.U8[numParams + j] = 0xFF;
    1046                     }
    1047                 }
    1048             }
    1049         }
    1050 
    1051         g = psMatrixSolveSVD(g, F, h);
    1052         if (!g) {
    1053             psError(PS_ERR_UNKNOWN, false, "SVD solution of least-squares equation failed.\n");
    1054             psFree(F);
    1055             psFree(h);
    1056             psFree(mask);
    1057             return NULL;
    1058         }
    1059         for (int i = 0; i < num; i++) {
    1060             if (mask->data.U8[i]) {
    1061                 g->data.F64[i] = 0.0;
    1062             }
    1063         }
    1064         psFree(mask);
    1065 
    1066         psFree(F);
    1067         psFree(h);
    1068 
    1069 #ifdef TESTING
    1070         {
    1071             psImage *image = psImageAlloc(1, num, PS_TYPE_F64);
    1072             psFits *fits = psFitsOpen("solnVector.fits", "w");
    1073             for (int i = 0; i < num; i++) {
    1074                 image->data.F64[0][i] = g->data.F64[i];
    1075             }
    1076             psFitsWriteImage(fits, NULL, image, 0, NULL);
    1077             psFree(image);
    1078             psFitsClose(fits);
    1079         }
    1080 #endif
    1081 
    10821104        if (!kernels->solution1) {
    10831105            kernels->solution1 = psVectorAlloc(numParams, PS_TYPE_F64);
    … …  
    10871109        }
    10881110
    1089         for (int i = 0; i < numParams; i++) {
    1090             kernels->solution1->data.F64[i] = g->data.F64[i];
    1091         }
    1092         for (int i = 0; i < numParams2; i++) {
    1093             kernels->solution2->data.F64[i] = g->data.F64[numParams + i];
    1094         }
    1095         psFree(g);
     1111        memcpy(kernels->solution1->data.F64, solution->data.F64, numParams * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
     1112        memcpy(kernels->solution2->data.F64, &solution->data.F64[numParams],
     1113               numParams2 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
     1114
     1115        psFree(solution);
    10961116
    10971117    }
Note: See TracChangeset for help on using the changeset viewer.