IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 29, 2010, 6:02:38 PM (16 years ago)
Author:
Paul Price
Message:

Fairly extensive change to make subtractions more stable on images with limited lit area in the presence of spatial variation. Took the polynomials which were scaled over the entire image and scale them only over the lit area. Removed some fprintf statements. Fixed the problem causing the background term to always come out 1.0.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r26737 r26739  
    213213
    214214    *norm = normI2 / normI1;
    215     fprintf(stderr, "Sums: %f %f %f, normWindow: %d\n", normI1, normI2, *norm, normWindow);
    216215
    217216    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
    … …  
    522521
    523522    *norm = normI2 / normI1;
    524     fprintf(stderr, "Sums: %f %f %f - I1I1: %f\n", normI1, normI2, *norm, sumI1I1);
    525523
    526524    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
    … …  
    590588                // Contribution to chi^2: a_i^2 P_i
    591589                psAssert(isfinite(penalties->data.F32[i]), "Invalid penalty");
    592                 fprintf (stderr, "penalty main: %d %e vs %e x %f : %f\n",
    593                          index,
    594                          matrix->data.F64[index][index],
    595                          norm, penalties->data.F32[i],
    596                          (norm * penalties->data.F32[i]) / matrix->data.F64[index][index] );
    597590                matrix->data.F64[index][index] += norm * penalties->data.F32[i];
    598591                if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
    599                     fprintf (stderr, "penalty dual: %d %e vs %e x %f : %f\n",
    600                              index + numParams + 2,
    601                              matrix->data.F64[index + numParams + 2][index + numParams + 2],
    602                              norm, penalties->data.F32[i],
    603                              (norm * penalties->data.F32[i]) / matrix->data.F64[index + numParams + 2][index + numParams + 2] );
    604592                    matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties->data.F32[i];
    605593                    // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1)
    … …  
    997985            // Solve normalisation and background separately
    998986            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     987            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
    999988
    1000989            psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm
    … …  
    10141003
    10151004            fprintf(stderr, "Norm: %lf\n", normValue);
    1016             // fprintf(stderr, "BG: %lf\n", bgValue);
    10171005
    10181006            // Solve kernel components
    10191007            for (int i = 0; i < numSolution1; i++) {
    1020                 sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; // + bgValue * sumMatrix->data.F64[bgIndex][i];
     1008                sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i];
    10211009
    10221010                sumMatrix->data.F64[i][normIndex] = 0.0;
    10231011                sumMatrix->data.F64[normIndex][i] = 0.0;
    1024 
    1025                 // sumMatrix->data.F64[i][bgIndex] = 0.0;
    1026                 // sumMatrix->data.F64[bgIndex][i] = 0.0;
    1027             }
     1012            }
     1013            sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex];
     1014            sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
     1015            sumMatrix->data.F64[normIndex][bgIndex] = 0.0;
    10281016
    10291017            sumMatrix->data.F64[normIndex][normIndex] = 1.0;
    1030             // sumMatrix->data.F64[bgIndex][bgIndex] = 1.0;
    1031             // sumMatrix->data.F64[normIndex][bgIndex] = sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
    1032 
    10331018            sumVector->data.F64[normIndex] = 0.0;
    1034             // sumVector->data.F64[bgIndex] = 0.0;
    10351019
    10361020            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
    10371021
    10381022            solution->data.F64[normIndex] = normValue;
    1039             // solution->data.F64[bgIndex] = bgValue;
    10401023        }
    10411024# endif
    10421025
    10431026        if (!kernels->solution1) {
    1044             kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);
    1045             psVectorInit (kernels->solution1, 0.0);
     1027            kernels->solution1 = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     1028            psVectorInit(kernels->solution1, 0.0);
    10461029        }
    10471030
    … …  
    11261109            // Solve normalisation and background separately
    11271110            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     1111            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
    11281112
    11291113#if 0
    … …  
    11591143
    11601144            double normValue = stats->robustMedian;
    1161             // double bgValue = 0.0;
    11621145
    11631146            psFree(stats);
    11641147
    11651148            fprintf(stderr, "Norm: %lf\n", normValue);
    1166             // fprintf(stderr, "BG: %lf\n", bgValue);
    11671149
    11681150            // Solve kernel components
    11691151            for (int i = 0; i < numSolution2; i++) {
    1170                 sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; // + bgValue * sumMatrix->data.F64[bgIndex][i];
    1171                 sumVector->data.F64[i + numSolution1] -= normValue * sumMatrix->data.F64[normIndex][i + numSolution1]; // + bgValue * sumMatrix->data.F64[bgIndex][i + numSolution1];
     1152                sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i];
     1153                sumVector->data.F64[i + numSolution1] -= normValue * sumMatrix->data.F64[normIndex][i + numSolution1];
    11721154
    11731155                sumMatrix->data.F64[i][normIndex] = 0.0;
    11741156                sumMatrix->data.F64[normIndex][i] = 0.0;
    11751157
    1176                 // sumMatrix->data.F64[i][bgIndex] = 0.0;
    1177                 // sumMatrix->data.F64[bgIndex][i] = 0.0;
    11781158                sumMatrix->data.F64[i + numSolution1][normIndex] = 0.0;
    11791159                sumMatrix->data.F64[normIndex][i + numSolution1] = 0.0;
    1180                 // sumMatrix->data.F64[i + numSolution1][bgIndex] = 0.0;
    1181                 // sumMatrix->data.F64[bgIndex][i + numSolution1] = 0.0;
    1182             }
     1160            }
     1161            sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex];
     1162            sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
     1163            sumMatrix->data.F64[normIndex][bgIndex] = 0.0;
    11831164
    11841165            sumMatrix->data.F64[normIndex][normIndex] = 1.0;
    1185             // sumMatrix->data.F64[bgIndex][bgIndex] = 1.0;
    1186             // sumMatrix->data.F64[normIndex][bgIndex] = sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
    11871166
    11881167            sumVector->data.F64[normIndex] = 0.0;
    1189             // sumVector->data.F64[bgIndex] = 0.0;
    11901168
    11911169            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
    11921170
    11931171            solution->data.F64[normIndex] = normValue;
    1194             // solution->data.F64[bgIndex] = bgValue;
    11951172        }
    11961173#endif
Note: See TracChangeset for help on using the changeset viewer.