IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 26, 2010, 9:18:39 AM (16 years ago)
Author:
Serge CHASTEL
Message:

Merging trunk in branch

Location:
branches/sc_branches/trunkTest
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/sc_branches/trunkTest

  • branches/sc_branches/trunkTest/psModules

    • Property svn:mergeinfo deleted

  • branches/sc_branches/trunkTest/psModules/src/imcombine/pmSubtractionEquation.c

    r28405 r29060  
    3838                                  const psImage *polyValues, // Spatial polynomial values
    3939                                  int footprint, // (Half-)Size of stamp
    40                                   int normWindow, // Window (half-)size for normalisation measurement
     40                                  int normWindow1, // Window (half-)size for normalisation measurement
     41                                  int normWindow2, // Window (half-)size for normalisation measurement
    4142                                  const pmSubtractionEquationCalculationMode mode
    4243                                  )
    … …  
    184185            double one = 1.0;
    185186
    186             if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) {
     187            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) {
    187188                normI1 += ref;
     189            }
     190            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) {
    188191                normI2 += in;
    189192            }
    … …  
    214217
    215218    *norm = normI2 / normI1;
     219
     220    fprintf (stderr, "normValue: %f %f %f\n", normI1, normI2, *norm);
    216221
    217222    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
    … …  
    262267                                      const psImage *polyValues, // Spatial polynomial values
    263268                                      int footprint, // (Half-)Size of stamp
    264                                       int normWindow, // Window (half-)size for normalisation measurement
     269                                      int normWindow1, // Window (half-)size for normalisation measurement
     270                                      int normWindow2, // Window (half-)size for normalisation measurement
    265271                                      const pmSubtractionEquationCalculationMode mode
    266272                                      )
    … …  
    492498            double i1i2 = i1 * i2;
    493499
    494             if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) {
     500            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) {
    495501                normI1 += i1;
     502            }
     503            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) {
    496504                normI2 += i2;
    497505            }
    … …  
    522530
    523531    *norm = normI2 / normI1;
     532    fprintf (stderr, "normValue: %f %f %f\n", normI1, normI2, *norm);
    524533
    525534    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
    … …  
    559568// Add in penalty term to least-squares vector
    560569bool calculatePenalty(psImage *matrix,                     // Matrix to which to add in penalty term
    561                              psVector *vector,                    // Vector to which to add in penalty term
    562                              const pmSubtractionKernels *kernels, // Kernel parameters
    563                              float norm                           // Normalisation
    564     )
     570                      psVector *vector,                    // Vector to which to add in penalty term
     571                      const pmSubtractionKernels *kernels, // Kernel parameters
     572                      float norm                           // Normalisation
     573  )
    565574{
    566575    if (kernels->penalty == 0.0) {
    … …  
    568577    }
    569578
    570     psVector *penalties = kernels->penalties; // Penalties for each kernel component
     579    psVector *penalties1 = kernels->penalties1; // Penalties for each kernel component (input)
     580    psVector *penalties2 = kernels->penalties2; // Penalties for each kernel component (ref)
     581
    571582    int spatialOrder = kernels->spatialOrder; // Order of spatial variations
    572583    int numKernels = kernels->num; // Number of kernel components
    … …  
    588599            for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
    589600                // Contribution to chi^2: a_i^2 P_i
    590                 psAssert(isfinite(penalties->data.F32[i]), "Invalid penalty");
    591                 matrix->data.F64[index][index] += norm * penalties->data.F32[i];
     601                psAssert(isfinite(penalties1->data.F32[i]), "Invalid penalty");
     602                fprintf (stderr, "penalty: %f + %f (%f * %f)\n", matrix->data.F64[index][index], norm * penalties1->data.F32[i], norm, penalties1->data.F32[i]);
     603                matrix->data.F64[index][index] += norm * penalties1->data.F32[i];
    592604                if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
    593                     matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties->data.F32[i];
     605                    fprintf (stderr, "penalty: (x^%d y^%d fwhm %f) : %f + %f (%f * %f)\n", kernels->u->data.S32[index], kernels->v->data.S32[index], kernels->widths->data.F32[index],
     606                             matrix->data.F64[index + numParams + 2][index + numParams + 2], norm * penalties2->data.F32[i], norm, penalties2->data.F32[i]);
     607                    matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties2->data.F32[i];                       
    594608                    // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1)
    595609                    // penalties scale with second moments
    … …  
    682696
    683697    pmSubtractionStampList *stamps = job->args->data[0]; // List of stamps
    684     const pmSubtractionKernels *kernels = job->args->data[1]; // Kernels
     698    pmSubtractionKernels *kernels = job->args->data[1]; // Kernels
    685699    int index = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index
    686700    pmSubtractionEquationCalculationMode mode  = PS_SCALAR_VALUE(job->args->data[3], S32); // calculation model
    … …  
    689703}
    690704
    691 bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels,
     705bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels,
    692706                                         int index, const pmSubtractionEquationCalculationMode mode)
    693707{
    … …  
    778792        status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image2, stamp->image1,
    779793                                       weight, window, stamp->convolutions1, kernels,
    780                                        polyValues, footprint, stamps->normWindow, mode);
     794                                       polyValues, footprint, stamps->normWindow1, stamps->normWindow2, mode);
    781795        break;
    782796      case PM_SUBTRACTION_MODE_2:
    783797        status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image1, stamp->image2,
    784798                                       weight, window, stamp->convolutions2, kernels,
    785                                        polyValues, footprint, stamps->normWindow, mode);
     799                                       polyValues, footprint, stamps->normWindow2, stamps->normWindow1, mode);
    786800        break;
    787801      case PM_SUBTRACTION_MODE_DUAL:
    … …  
    789803                                           stamp->image1, stamp->image2,
    790804                                           weight, window, stamp->convolutions1, stamp->convolutions2,
    791                                            kernels, polyValues, footprint, stamps->normWindow, mode);
     805                                           kernels, polyValues, footprint, stamps->normWindow1, stamps->normWindow2, mode);
    792806        break;
    793807      default:
    … …  
    830844}
    831845
    832 bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels,
     846bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels,
    833847                                    const pmSubtractionEquationCalculationMode mode)
    834848{
    … …  
    9961010            }
    9971011
     1012            // double normValue = 1.0;
    9981013            double normValue = stats->robustMedian;
    9991014            // double bgValue = 0.0;
    … …  
    10231038        }
    10241039# endif
     1040
     1041#if (1)
     1042        for (int i = 0; i < solution->n; i++) {
     1043            fprintf(stderr, "Single solution %d: %lf\n", i, solution->data.F64[i]);
     1044        }
     1045#endif
    10251046
    10261047        if (!kernels->solution1) {
    … …  
    10961117
    10971118        int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
    1098         calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 1000.0);
     1119        calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 100.0);
    10991120#endif
    11001121
    … …  
    11771198
    11781199
    1179 #ifdef TESTING
     1200#if (1)
    11801201        for (int i = 0; i < solution->n; i++) {
    11811202            fprintf(stderr, "Dual solution %d: %lf\n", i, solution->data.F64[i]);
Note: See TracChangeset for help on using the changeset viewer.