IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 3, 2009, 12:20:59 PM (17 years ago)
Author:
eugene
Message:

add code to calculate norm separately from kernel coeffs

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r26330 r26332  
    2626
    2727// Calculate the least-squares matrix and vector
     28// XXX why are we calculating these values on each iteration?  aren't they identical (no change in pixel values...)
    2829static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated
    2930                                  psVector *vector, // Least-squares vector, updated
    … …  
    3637                                  const psImage *polyValues, // Spatial polynomial values
    3738                                  int footprint, // (Half-)Size of stamp
    38                                   const bool normOnly
     39                                  const pmSubtractionEquationCalculationMode mode
    3940                                  )
    4041{
    … …  
    5455
    5556    // Evaluate polynomial-polynomial terms
    56     // XXX we can skip this for normOnly
     57    // XXX we can skip this if we are not calculating kernel coeffs
    5758    for (int iyOrder = 0, iIndex = 0; iyOrder <= spatialOrder; iyOrder++) {
    5859        for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex++) {
    … …  
    6869    }
    6970
    70 
    71     // if we only want to calculate the normalization (and background model?)
    72     // then we should assign these values in the matrix and vector to 1 (diagonal) or 0 (off-diagonal)
    73     // XXX need to disable normalization calculation if !normOnly
    74     if (normOnly) {
    75         psImageInit(matrix, 0.0);
    76         psVectorInit(vector, 1.0);
    77         for (int i = 0; i < matrix->numCols; i++) {
    78             matrix->data.F64[i][i] = 1.0;
    79         }
    80     } else {
    81         for (int i = 0; i < numKernels; i++) {
    82             psKernel *iConv = convolutions->data[i]; // Convolution for index i
    83             for (int j = i; j < numKernels; j++) {
    84                 psKernel *jConv = convolutions->data[j]; // Convolution for index j
    85 
    86                 double sumCC = 0.0;         // Sum of convolution products
    87                 for (int y = - footprint; y <= footprint; y++) {
    88                     for (int x = - footprint; x <= footprint; x++) {
    89                         double cc = iConv->kernel[y][x] * jConv->kernel[y][x];
    90 #ifdef USE_WEIGHT
     71    // initialize the matrix and vector for NOP on all coeffs.  we only fill in the coeffs we
     72    // choose to calculate
     73    psImageInit(matrix, 0.0);
     74    psVectorInit(vector, 1.0);
     75    for (int i = 0; i < matrix->numCols; i++) {
     76        matrix->data.F64[i][i] = 1.0;
     77    }
     78
     79    // the order of the elements in the matrix and vector is:
     80    // [kernel 0, x^0 y^0][kernel 1 x^0 y^0]...[kernel N, x^0 y^0]
     81    // [kernel 0, x^1 y^0][kernel 1 x^1 y^0]...[kernel N, x^1 y^0]
     82    // [kernel 0, x^n y^m][kernel 1 x^n y^m]...[kernel N, x^n y^m]
     83    // normalization
     84    // bg 0, bg 1, bg 2 (only 0 is currently used?)
     85
     86    for (int i = 0; i < numKernels; i++) {
     87        psKernel *iConv = convolutions->data[i]; // Convolution for index i
     88        for (int j = i; j < numKernels; j++) {
     89            psKernel *jConv = convolutions->data[j]; // Convolution for index j
     90
     91            double sumCC = 0.0;         // Sum of convolution products
     92            for (int y = - footprint; y <= footprint; y++) {
     93                for (int x = - footprint; x <= footprint; x++) {
     94                    double cc = iConv->kernel[y][x] * jConv->kernel[y][x];
     95                    if (weight) {
    9196                        cc *= weight->kernel[y][x];
    92 #endif
    93 #ifdef USE_WINDOW
    94                         if (window) {
    95                             cc *= window->kernel[y][x];
    96                         }
    97 #endif
    98                         sumCC += cc;
    9997                    }
     98                    if (window) {
     99                        cc *= window->kernel[y][x];
     100                    }
     101                    sumCC += cc;
    100102                }
    101 
    102                 // Spatial variation
     103            }
     104
     105            // Spatial variation of kernel coeffs
     106            if (mode | PM_SUBTRACTION_EQUATION_KERNELS) {
    103107                for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
    104108                    for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
    … …  
    109113                }
    110114            }
    111 
    112             double sumRC = 0.0;             // Sum of the reference-convolution products
    113             double sumIC = 0.0;             // Sum of the input-convolution products
    114             double sumC = 0.0;              // Sum of the convolution
    115             for (int y = - footprint; y <= footprint; y++) {
    116                 for (int x = - footprint; x <= footprint; x++) {
    117                     float conv = iConv->kernel[y][x];
    118                     float in = input->kernel[y][x];
    119                     float ref = reference->kernel[y][x];
    120                     double ic = in * conv;
    121                     double rc = ref * conv;
    122                     double c = conv;
    123 #ifdef USE_WEIGHT
     115        }
     116
     117        double sumRC = 0.0;             // Sum of the reference-convolution products
     118        double sumIC = 0.0;             // Sum of the input-convolution products
     119        double sumC = 0.0;              // Sum of the convolution
     120        for (int y = - footprint; y <= footprint; y++) {
     121            for (int x = - footprint; x <= footprint; x++) {
     122                float conv = iConv->kernel[y][x];
     123                float in = input->kernel[y][x];
     124                float ref = reference->kernel[y][x];
     125                double ic = in * conv;
     126                double rc = ref * conv;
     127                double c = conv;
     128                if (weight) {
    124129                    float wtVal = weight->kernel[y][x];
    125130                    ic *= wtVal;
    126131                    rc *= wtVal;
    127132                    c *= wtVal;
    128 #endif
    129 #ifdef USE_WINDOW
    130                     if (window) {
    131                         float winVal = window->kernel[y][x];
    132                         ic *= winVal;
    133                         rc *= winVal;
    134                         c  *= winVal;
    135                     }
    136 #endif
    137                     sumIC += ic;
    138                     sumRC += rc;
    139                     sumC += c;
    140133                }
     134                if (window) {
     135                    float winVal = window->kernel[y][x];
     136                    ic *= winVal;
     137                    rc *= winVal;
     138                    c  *= winVal;
     139                }
     140                sumIC += ic;
     141                sumRC += rc;
     142                sumC += c;
    141143            }
    142             // Spatial variation
    143             for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
    144                 double normTerm = sumRC * poly[iTerm];
    145                 double bgTerm = sumC * poly[iTerm];
     144        }
     145        // Spatial variation
     146        for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
     147            double normTerm = sumRC * poly[iTerm];
     148            double bgTerm = sumC * poly[iTerm];
     149            if ((mode | PM_SUBTRACTION_EQUATION_NORM) && (mode | PM_SUBTRACTION_EQUATION_KERNELS)) {
    146150                matrix->data.F64[iIndex][normIndex] = normTerm;
    147151                matrix->data.F64[normIndex][iIndex] = normTerm;
     152            }
     153            if ((mode | PM_SUBTRACTION_EQUATION_BG) && (mode | PM_SUBTRACTION_EQUATION_KERNELS)) {
    148154                matrix->data.F64[iIndex][bgIndex] = bgTerm;
    149155                matrix->data.F64[bgIndex][iIndex] = bgTerm;
     156            }
     157            if (mode | PM_SUBTRACTION_EQUATION_KERNELS) {
    150158                vector->data.F64[iIndex] = sumIC * poly[iTerm];
    151159            }
    152160        }
    153161    }
    154 
    155     // XXX need to disable normalization calculation if !normOnly
    156162
    157163    double sumRR = 0.0;                 // Sum of the reference product
    … …  
    167173            double rr = PS_SQR(ref);
    168174            double one = 1.0;
    169 #ifdef USE_WEIGHT
    170             float wtVal = weight->kernel[y][x];
    171             rr *= wtVal;
    172             ir *= wtVal;
    173             in *= wtVal;
    174             ref *= wtVal;
    175             one *= wtVal;
    176 #endif
    177 #ifdef USE_WINDOW
     175            if (weight) {
     176                float wtVal = weight->kernel[y][x];
     177                rr *= wtVal;
     178                ir *= wtVal;
     179                in *= wtVal;
     180                ref *= wtVal;
     181                one *= wtVal;
     182            }
    178183            if (window) {
    179184                float  winVal = window->kernel[y][x];
    … …  
    184189                one *= winVal;
    185190            }
    186 #endif
    187191            sumRR += rr;
    188192            sumIR += ir;
    … …  
    192196        }
    193197    }
    194     matrix->data.F64[normIndex][normIndex] = sumRR;
    195     matrix->data.F64[bgIndex][bgIndex] = sum1;
    196     matrix->data.F64[normIndex][bgIndex] = matrix->data.F64[bgIndex][normIndex] = sumR;
    197     vector->data.F64[normIndex] = sumIR;
    198     vector->data.F64[bgIndex] = sumI;
    199 
     198    if (mode | PM_SUBTRACTION_EQUATION_NORM) {
     199        matrix->data.F64[normIndex][normIndex] = sumRR;
     200        vector->data.F64[normIndex] = sumIR;
     201    }
     202    if (mode | PM_SUBTRACTION_EQUATION_BG) {
     203        matrix->data.F64[bgIndex][bgIndex] = sum1;
     204        vector->data.F64[bgIndex] = sumI;
     205    }
     206    if ((mode | PM_SUBTRACTION_EQUATION_NORM) && (mode | PM_SUBTRACTION_EQUATION_BG)) {
     207        matrix->data.F64[normIndex][bgIndex] = sumR;
     208        matrix->data.F64[bgIndex][normIndex] = sumR;
     209    }
    200210    return true;
    201211}
    … …  
    267277                    double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x];
    268278                    double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x];
    269 #ifdef USE_WEIGHT
    270                     float wtVal = weight->kernel[y][x];
    271                     aa *= wtVal;
    272                     bb *= wtVal;
    273                     ab *= wtVal;
    274 #endif
     279                    if (weight) {
     280                        float wtVal = weight->kernel[y][x];
     281                        aa *= wtVal;
     282                        bb *= wtVal;
     283                        ab *= wtVal;
     284                    }
    275285                    sumAA += aa;
    276286                    sumBB += bb;
    … …  
    299309                for (int x = - footprint; x <= footprint; x++) {
    300310                    double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x];
    301 #ifdef USE_WEIGHT
    302                     ab *= weight->kernel[y][x];
    303 #endif
     311                    if (weight) {
     312                        ab *= weight->kernel[y][x];
     313                    }
    304314                    sumAB += ab;
    305315                }
    … …  
    336346                double i1i2 = i1 * i2;
    337347
    338 #ifdef USE_WEIGHT
    339                 float wtVal = weight->kernel[y][x];
    340                 ai2 *= wtVal;
    341                 bi2 *= wtVal;
    342                 ai1 *= wtVal;
    343                 bi1 *= wtVal;
    344                 i1i2 *= wtVal;
    345                 a *= wtVal;
    346                 b *= wtVal;
    347                 i2 *= wtVal;
    348 #endif
    349 
     348                if (weight) {
     349                    float wtVal = weight->kernel[y][x];
     350                    ai2 *= wtVal;
     351                    bi2 *= wtVal;
     352                    ai1 *= wtVal;
     353                    bi1 *= wtVal;
     354                    i1i2 *= wtVal;
     355                    a *= wtVal;
     356                    b *= wtVal;
     357                    i2 *= wtVal;
     358                }
    350359                sumAI2 += ai2;
    351360                sumBI2 += bi2;
    … …  
    392401            double i1i2 = i1 * i2;
    393402
    394 #ifdef USE_WEIGHT
    395             float wtVal = weight->kernel[y][x];
    396             i1 *= wtVal;
    397             i1i1 *= wtVal;
    398             one *= wtVal;
    399             i2 *= wtVal;
    400             i1i2 *= wtVal;
    401 #endif
    402 
     403            if (weight) {
     404                float wtVal = weight->kernel[y][x];
     405                i1 *= wtVal;
     406                i1i1 *= wtVal;
     407                one *= wtVal;
     408                i2 *= wtVal;
     409                i1i2 *= wtVal;
     410            }
    403411            sumI1 += i1;
    404412            sumI1I1 += i1i1;
    … …  
    588596    const pmSubtractionKernels *kernels = job->args->data[1]; // Kernels
    589597    int index = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index
    590 
    591     return pmSubtractionCalculateEquationStamp(stamps, kernels, index, normOnly);
     598    pmSubtractionEquationCalculationMode mode  = PS_SCALAR_VALUE(job->args->data[3], S32); // calculation model
     599
     600    return pmSubtractionCalculateEquationStamp(stamps, kernels, index, mode);
    592601}
    593602
    594603bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels,
    595                                          int index, bool normOnly)
     604                                         int index, const pmSubtractionEquationCalculationMode mode)
    596605{
    597606    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
    … …  
    662671
    663672    bool status;                    // Status of least-squares matrix/vector calculation
     673
     674    psKernel *weight = NULL;
     675    psKernel *window = NULL;
     676
     677#ifdef USE_WEIGHT
     678    weight = stamp->weight;
     679#endif
     680#ifdef USE_WINDOW
     681    window = stamps->window;
     682#endif
     683
    664684    switch (kernels->mode) {
    665685      case PM_SUBTRACTION_MODE_1:
    666686        status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image2, stamp->image1,
    667                                        stamp->weight, stamps->window, stamp->convolutions1, kernels,
    668                                        polyValues, footprint, normOnly);
     687                                       weight, window, stamp->convolutions1, kernels,
     688                                       polyValues, footprint, mode);
    669689        break;
    670690      case PM_SUBTRACTION_MODE_2:
    671691        status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image1, stamp->image2,
    672                                        stamp->weight, stamps->window, stamp->convolutions2, kernels,
    673                                        polyValues, footprint, normOnly);
     692                                       weight, window, stamp->convolutions2, kernels,
     693                                       polyValues, footprint, mode);
    674694        break;
    675695      case PM_SUBTRACTION_MODE_DUAL:
    676         psAbort ("dual is disabled: need to add window and normOnly");
     696        psAbort ("dual is disabled: need to add window and calculation mode");
    677697        if (new) {
    678698            stamp->matrix2 = psImageAlloc(numKernels * numSpatial, numKernels * numSpatial, PS_TYPE_F64);
    … …  
    686706#endif
    687707        status = calculateDualMatrixVector(stamp->matrix1, stamp->vector1, stamp->matrix2, stamp->vector2,
    688                                            stamp->matrixX, stamp->image1, stamp->image2, stamp->weight,
     708                                           stamp->matrixX, stamp->image1, stamp->image2, weight,
    689709                                           stamp->convolutions1, stamp->convolutions2, kernels, polyValues,
    690710                                           footprint);
    … …  
    758778}
    759779
    760 bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, const bool normOnly)
     780bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, const pmSubtractionEquationCalculationMode mode)
    761781{
    762782    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
    … …  
    778798            psArrayAdd(job->args, 1, (pmSubtractionKernels*)kernels); // Casting away const to put on array
    779799            PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
    780             PS_ARRAY_ADD_SCALAR(job->args, normOnly, PS_TYPE_BOOL);
     800            PS_ARRAY_ADD_SCALAR(job->args, mode, PS_TYPE_S32);
    781801            if (!psThreadJobAddPending(job)) {
    782802                psFree(job);
    … …  
    785805            psFree(job);
    786806        } else {
    787             pmSubtractionCalculateEquationStamp(stamps, kernels, i, normOnly);
     807            pmSubtractionCalculateEquationStamp(stamps, kernels, i, mode);
    788808        }
    789809    }
    … …  
    803823}
    804824
    805 bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps)
     825bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels,
     826                                const pmSubtractionStampList *stamps,
     827                                const pmSubtractionEquationCalculationMode mode)
    806828{
    807829    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
    … …  
    941963            return NULL;
    942964        }
    943         kernels->solution1 = psMatrixLUSolution(kernels->solution1, luMatrix, sumVector, permutation);
     965
     966        psVector *solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation);
     967        psFree(sumVector);
     968        psFree(luMatrix);
     969        psFree(permutation);
     970        if (!solution) {
     971            psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
     972            return NULL;
     973        }
     974
     975        if (!kernels->solution1) {
     976            kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);
     977            psVectorInit (kernels->solution1, 0.0);
     978        }
     979
     980        // only update the solutions that we chose to calculate:
     981        if (mode | PM_SUBTRACTION_EQUATION_NORM) {
     982            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     983            kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
     984        }
     985        if (mode | PM_SUBTRACTION_EQUATION_BG) {
     986            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
     987            kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
     988        }
     989        if (mode | PM_SUBTRACTION_EQUATION_KERNELS) {
     990            int numKernels = kernels->num;
     991            int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
     992            int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
     993            for (int i = 0; i < numKernels * numPoly; i++) {
     994                kernels->solution1->data.F64[i] = solution->data.F64[i];
     995            }
     996        }
     997        psFree (solution);
    944998
    945999#ifdef TESTING
    … …  
    9521006#endif
    9531007
    954         psFree(sumVector);
    955         psFree(luMatrix);
    956         psFree(permutation);
    957         if (!kernels->solution1) {
    958             psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
    959             return NULL;
    960         }
    9611008    } else {
    9621009        // Dual convolution solution
Note: See TracChangeset for help on using the changeset viewer.