IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 28866


Ignore:
Timestamp:
Aug 9, 2010, 8:27:57 AM (16 years ago)
Author:
eugene
Message:

separate penalties for kernels for image 1 and 2, convolve second moments with fwhm to get penalties; separate normalization windows for images 1 and 2; update visual to have facilities like psphot visual; curve of growth analysis for norm window

Location:
branches/eam_branches/ipp-20100621/psModules/src/imcombine
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtraction.c

    r28796 r28866  
    3333#define USE_KERNEL_ERR                  // Use kernel error image?
    3434#define NUM_COVAR_POS 5                 // Number of positions for covariance calculation
     35
     36// XXX we need to pass these fwhm values elsewhere.  These should go on one of the structure, but
     37// things are too confusing to do that now.  just save them here.
     38static float FWHM1 = NAN;
     39static float FWHM2 = NAN;
    3540
    3641//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    … …  
    752757
    753758
    754 bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, const pmSubtractionKernels *kernels, int footprint)
     759bool pmSubtractionConvolveStamp (pmSubtractionStamp *stamp, pmSubtractionKernels *kernels, int footprint)
    755760{
    756761    PS_ASSERT_PTR_NON_NULL(stamp, false);
    … …  
    774779        stamp->convolutions1 = convolveStamp(stamp->convolutions1, stamp->image1, kernels, footprint);
    775780        stamp->convolutions2 = convolveStamp(stamp->convolutions2, stamp->image2, kernels, footprint);
     781        if (!pmSubtractionKernelPenaltiesStamp(stamp, kernels)) {
     782            psAbort("failure in penalties");
     783        }
    776784        break;
    777785      default:
    … …  
    14131421    return true;
    14141422}
     1423
     1424bool pmSubtractionGetFWHMs(float *fwhm1, float *fwhm2) {
     1425
     1426  *fwhm1 = FWHM1;
     1427  *fwhm2 = FWHM2;
     1428  return true;
     1429}
     1430
     1431bool pmSubtractionSetFWHMs(float fwhm1, float fwhm2) {
     1432
     1433  FWHM1 = fwhm1;
     1434  FWHM2 = fwhm2;
     1435  return true;
     1436}

  • branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtraction.h

    r26893 r28866  
    5959/// Convolve the reference stamp with the kernel components
    6060bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, ///< Stamp to convolve
    61                                 const pmSubtractionKernels *kernels, ///< Kernel parameters
     61                                pmSubtractionKernels *kernels, ///< Kernel parameters
    6262                                int footprint ///< Half-size of region over which to calculate equation
    6363    );
    … …  
    157157    );
    158158
     159bool pmSubtractionGetFWHMs(float *fwhm1, float *fwhm2);
     160bool pmSubtractionSetFWHMs(float fwhm1, float fwhm2);
     161
    159162/// @}
    160163#endif

  • branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionEquation.c

    r28405 r28866  
    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: %f + %f (%f * %f)\n", matrix->data.F64[index + numParams + 2][index + numParams + 2], norm * penalties2->data.F32[i], norm, penalties2->data.F32[i]);
     606                    matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties2->data.F32[i];                       
    594607                    // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1)
    595608                    // penalties scale with second moments
    … …  
    682695
    683696    pmSubtractionStampList *stamps = job->args->data[0]; // List of stamps
    684     const pmSubtractionKernels *kernels = job->args->data[1]; // Kernels
     697    pmSubtractionKernels *kernels = job->args->data[1]; // Kernels
    685698    int index = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index
    686699    pmSubtractionEquationCalculationMode mode  = PS_SCALAR_VALUE(job->args->data[3], S32); // calculation model
    … …  
    689702}
    690703
    691 bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels,
     704bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels,
    692705                                         int index, const pmSubtractionEquationCalculationMode mode)
    693706{
    … …  
    778791        status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image2, stamp->image1,
    779792                                       weight, window, stamp->convolutions1, kernels,
    780                                        polyValues, footprint, stamps->normWindow, mode);
     793                                       polyValues, footprint, stamps->normWindow1, stamps->normWindow2, mode);
    781794        break;
    782795      case PM_SUBTRACTION_MODE_2:
    783796        status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image1, stamp->image2,
    784797                                       weight, window, stamp->convolutions2, kernels,
    785                                        polyValues, footprint, stamps->normWindow, mode);
     798                                       polyValues, footprint, stamps->normWindow2, stamps->normWindow1, mode);
    786799        break;
    787800      case PM_SUBTRACTION_MODE_DUAL:
    … …  
    789802                                           stamp->image1, stamp->image2,
    790803                                           weight, window, stamp->convolutions1, stamp->convolutions2,
    791                                            kernels, polyValues, footprint, stamps->normWindow, mode);
     804                                           kernels, polyValues, footprint, stamps->normWindow1, stamps->normWindow2, mode);
    792805        break;
    793806      default:
    … …  
    830843}
    831844
    832 bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels,
     845bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels,
    833846                                    const pmSubtractionEquationCalculationMode mode)
    834847{
    … …  
    9961009            }
    9971010
     1011            // double normValue = 1.0;
    9981012            double normValue = stats->robustMedian;
    9991013            // double bgValue = 0.0;
    … …  
    10231037        }
    10241038# endif
     1039
     1040#if (1)
     1041        for (int i = 0; i < solution->n; i++) {
     1042            fprintf(stderr, "Single solution %d: %lf\n", i, solution->data.F64[i]);
     1043        }
     1044#endif
    10251045
    10261046        if (!kernels->solution1) {
    … …  
    10961116
    10971117        int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
    1098         calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 1000.0);
     1118        calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 10000.0);
    10991119#endif
    11001120
    … …  
    11771197
    11781198
    1179 #ifdef TESTING
     1199#if (1)
    11801200        for (int i = 0; i < solution->n; i++) {
    11811201            fprintf(stderr, "Dual solution %d: %lf\n", i, solution->data.F64[i]);

  • branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionEquation.h

    r26893 r28866  
    1919/// Calculate the least-squares equation to match the image quality for a single stamp
    2020bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, ///< Stamps
    21                                          const pmSubtractionKernels *kernels, ///< Kernel parameters
     21                                         pmSubtractionKernels *kernels, ///< Kernel parameters
    2222                                         int index, ///< Index of stamp
    2323                                         const pmSubtractionEquationCalculationMode mode
    … …  
    2626/// Calculate the least-squares equation to match the image quality
    2727bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps
    28                                     const pmSubtractionKernels *kernels, ///< Kernel parameters
     28                                    pmSubtractionKernels *kernels, ///< Kernel parameters
    2929                                    const pmSubtractionEquationCalculationMode mode
    3030    );

  • branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionKernels.c

    r27335 r28866  
    2626    psFree(kernels->vStop);
    2727    psFree(kernels->preCalc);
    28     psFree(kernels->penalties);
     28    psFree(kernels->penalties1);
     29    psFree(kernels->penalties2);
    2930    psFree(kernels->solution1);
    3031    psFree(kernels->solution2);
    … …  
    140141    kernels->v = psVectorRealloc(kernels->v, start + numNew);
    141142    kernels->preCalc = psArrayRealloc(kernels->preCalc, start + numNew);
    142     kernels->penalties = psVectorRealloc(kernels->penalties, start + numNew);
     143
     144    kernels->penalties1 = psVectorRealloc(kernels->penalties1, start + numNew);
     145    kernels->penalties2 = psVectorRealloc(kernels->penalties2, start + numNew);
     146
    143147    kernels->inner = start;
    144148    kernels->num += numNew;
    … …  
    156160            kernels->v->data.S32[index] = v;
    157161            kernels->preCalc->data[index] = NULL;
    158             kernels->penalties->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v));
    159             psAssert (isfinite(kernels->penalties->data.F32[index]), "invalid penalty");
     162
     163            // XXX this needs to be changed to use the *convolved* second moment
     164            kernels->penalties1->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v));
     165            psAssert (isfinite(kernels->penalties1->data.F32[index]), "invalid penalty");
     166
     167            kernels->penalties2->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v));
     168            psAssert (isfinite(kernels->penalties2->data.F32[index]), "invalid penalty");
     169
    160170            psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v);
    161171        }
    … …  
    166176    kernels->v->n = start + numNew;
    167177    kernels->preCalc->n = start + numNew;
    168     kernels->penalties->n = start + numNew;
     178
     179    kernels->penalties1->n = start + numNew;
     180    kernels->penalties2->n = start + numNew;
    169181
    170182    return true;
    171183}
    172184
    173 bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc,
    174                                          int index, int size, int uOrder, int vOrder, float fwhm,
    175                                          bool AlardLuptonStyle, bool forceZeroNull)
     185static bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc,
     186                                                int index, int uOrder, int vOrder, float fwhm,
     187                                                bool AlardLuptonStyle, bool forceZeroNull)
    176188{
    177189    // we have 4 cases here:
    … …  
    182194
    183195    // Calculate moments
    184     double penalty = 0.0;                   // Moment, for penalty
    185196    double sum = 0.0, sum2 = 0.0;           // Sum of kernel component
    186197    float min = INFINITY, max = -INFINITY;  // Minimum and maximum kernel value
    187     for (int v = -size; v <= size; v++) {
    188         for (int u = -size; u <= size; u++) {
     198
     199    for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
     200        for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) {
    189201            double value = preCalc->kernel->kernel[v][u];
    190202            double value2 = PS_SQR(value);
    191203            sum += value;
    192204            sum2 += value2;
    193             penalty += value2 * PS_SQR((PS_SQR(u) + PS_SQR(v)));
    194205            min = PS_MIN(value, min);
    195206            max = PS_MAX(value, max);
    … …  
    198209
    199210#if 0
    200     fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf, moment: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max, penalty);
     211    fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max);
    201212#endif
    202213
    … …  
    239250
    240251    psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "*", psScalarAlloc(scale2D, PS_TYPE_F32));
    241     penalty *= 1.0 / sum2;
    242252
    243253    if (zeroNull) {
    … …  
    249259        double sum = 0.0;   // Sum of kernel component
    250260        float min = INFINITY, max = -INFINITY;  // Minimum and maximum kernel value
    251         for (int v = -size; v <= size; v++) {
    252             for (int u = -size; u <= size; u++) {
     261        for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
     262            for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) {
    253263                sum += preCalc->kernel->kernel[v][u];
    254264                min = PS_MIN(preCalc->kernel->kernel[v][u], min);
    … …  
    267277    }
    268278    kernels->preCalc->data[index] = preCalc;
    269     kernels->penalties->data.F32[index] = kernels->penalty * penalty;
    270     psAssert (isfinite(kernels->penalties->data.F32[index]), "invalid penalty");
    271     psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, fwhm, uOrder, vOrder, penalty);
     279    psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, fwhm, uOrder, vOrder);
    272280
    273281    return true;
    … …  
    321329
    322330                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values
    323                 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
    324                 // pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], false, false);
     331                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
     332                // pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], false, false);
    325333            }
    326334        }
    … …  
    379387            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    380388                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values
    381                 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
     389                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
    382390            }
    383391        }
    … …  
    385393            // XXX modify size for hermitians to account for sqrt(2) in Hermitian definition (relative to ISIS Gaussian)
    386394            pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS_RADIAL, order, order, size, sigma / sqrt(2.0)); // structure to hold precalculated values
    387             pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, order, order, fwhms->data.F32[i], true, true);
     395            pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, order, order, fwhms->data.F32[i], true, true);
    388396        }
    389397    }
    … …  
    437445            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    438446                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values
    439                 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
     447                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
    440448            }
    441449        }
    … …  
    506514
    507515                // XXX do we use Alard-Lupton normalization (last param true) or not?
    508                 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
     516                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
    509517
    510518                // XXXX test demo that deconvolved kernel is valid
    … …  
    572580    kernels->preCalc = psArrayAlloc(numBasisFunctions);
    573581    kernels->penalty = penalty;
    574     kernels->penalties = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
     582    kernels->penalties1 = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
     583    psVectorInit(kernels->penalties1, NAN);
     584    kernels->penalties2 = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
     585    psVectorInit(kernels->penalties2, NAN);
     586    kernels->havePenalties = false;
    575587    kernels->size = size;
    576588    kernels->inner = 0;
    … …  
    771783
    772784    psWarning("Kernel penalty for dual-convolution is not configured for SPAM kernels.");
    773     psVectorInit(kernels->penalties, 0.0);
     785    psVectorInit(kernels->penalties1, 0.0);
     786    psVectorInit(kernels->penalties2, 0.0);
    774787
    775788    return kernels;
    … …  
    866879
    867880    psWarning("Kernel penalty for dual-convolution is not configured for FRIES kernels.");
    868     psVectorInit(kernels->penalties, 0.0);
     881    psVectorInit(kernels->penalties1, 0.0);
     882    psVectorInit(kernels->penalties2, 0.0);
    869883
    870884    return kernels;
    … …  
    10401054                kernels->u->data.S32[index] = uOrder;
    10411055                kernels->v->data.S32[index] = vOrder;
    1042                 kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
    1043                 if (!isfinite(kernels->penalties->data.F32[index])) {
     1056
     1057                // XXX convert to use the convolved 2nd moment
     1058                kernels->penalties1->data.F32[index] = kernels->penalty * fabsf(moment);
     1059                if (!isfinite(kernels->penalties1->data.F32[index])) {
     1060                    psAbort ("invalid penalty");
     1061                }
     1062                kernels->penalties2->data.F32[index] = kernels->penalty * fabsf(moment);
     1063                if (!isfinite(kernels->penalties2->data.F32[index])) {
    10441064                    psAbort ("invalid penalty");
    10451065                }
    … …  
    12471267    out->preCalc = psMemIncrRefCounter(in->preCalc);
    12481268    out->penalty = in->penalty;
    1249     out->penalties = psMemIncrRefCounter(in->penalties);
     1269    out->penalties1 = psMemIncrRefCounter(in->penalties1);
     1270    out->penalties2 = psMemIncrRefCounter(in->penalties2);
    12501271    out->uStop = psMemIncrRefCounter(in->uStop);
    12511272    out->vStop = psMemIncrRefCounter(in->vStop);

  • branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionKernels.h

    r26893 r28866  
    3939    psArray *preCalc;                   ///< Array of images containing pre-calculated kernel (for ISIS, HERM or DECONV_HERM)
    4040    float penalty;                      ///< Penalty for wideness
    41     psVector *penalties;                ///< Penalty for each kernel component
     41    psVector *penalties1;               ///< Penalty for each kernel component
     42    psVector *penalties2;               ///< Penalty for each kernel component
     43    bool havePenalties;                 ///< flag to test if we have already calculated the penalties or not.
    4244    int size;                           ///< The half-size of the kernel
    4345    int inner;                          ///< The size of an inner region
    … …  
    308310    );
    309311
    310 
    311312#endif

  • branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionMatch.c

    r28405 r28866  
    13041304    float scale = PS_MAX(fwhm1, fwhm2) / scaleRef;      // Scaling factor
    13051305
     1306    // XXX save these values in a static for later use
     1307    pmSubtractionSetFWHMs(fwhm1, fwhm2);
     1308
    13061309    if (isfinite(scaleMin) && scale < scaleMin) {
    13071310        scale = scaleMin;

  • branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionMatch.h

    r26893 r28866  
    110110    );
    111111
    112 
    113112#endif

  • branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionStamps.c

    r28643 r28866  
    5151    psFree(list->y);
    5252    psFree(list->flux);
    53     psFree(list->window);
     53    psFree(list->window1);
     54    psFree(list->window2);
    5455}
    5556
    … …  
    230231    list->y = NULL;
    231232    list->flux = NULL;
    232     list->window = NULL;
    233233    list->normFrac = normFrac;
    234     list->normWindow = 0;
     234    list->window1 = NULL;
     235    list->window2 = NULL;
     236    list->normWindow1 = 0;
     237    list->normWindow2 = 0;
    235238    list->footprint = footprint;
    236239    list->sysErr = sysErr;
    … …  
    253256    out->y = NULL;
    254257    out->flux = NULL;
    255     out->window = psMemIncrRefCounter(in->window);
     258    out->window1 = psMemIncrRefCounter(in->window1);
     259    out->window2 = psMemIncrRefCounter(in->window2);
    256260    out->footprint = in->footprint;
    257     out->normWindow = in->normWindow;
     261    out->normWindow1 = in->normWindow1;
     262    out->normWindow2 = in->normWindow2;
    258263
    259264    for (int i = 0; i < num; i++) {
    … …  
    643648    int size = stamps->footprint; // Size of postage stamps
    644649
    645     psFree (stamps->window);
    646     stamps->window = psKernelAlloc(-size, size, -size, size);
    647     psImageInit(stamps->window->image, 0.0);
     650    psFree (stamps->window1);
     651    stamps->window1 = psKernelAlloc(-size, size, -size, size);
     652    psImageInit(stamps->window1->image, 0.0);
     653
     654    psFree (stamps->window2);
     655    stamps->window2 = psKernelAlloc(-size, size, -size, size);
     656    psImageInit(stamps->window2->image, 0.0);
    648657
    649658    // generate normalizations for each stamp
    … …  
    674683
    675684    // generate the window pixels
    676     double sum = 0.0;                   // Sum inside the window
    677     float maxValue = 0.0;               // Maximum value, for normalisation
     685    double sum1 = 0.0;                   // Sum inside the window
     686    double sum2 = 0.0;                   // Sum inside the window
     687    float maxValue1 = 0.0;               // Maximum value, for normalisation
     688    float maxValue2 = 0.0;               // Maximum value, for normalisation
    678689    for (int y = -size; y <= size; y++) {
    679690        for (int x = -size; x <= size; x++) {
    … …  
    696707            }
    697708            float f1 = stats->robustMedian;
     709
    698710            psStatsInit (stats);
    699711            if (!psVectorStats (stats, flux2, NULL, NULL, 0)) {
    … …  
    702714            float f2 = stats->robustMedian;
    703715
    704             stamps->window->kernel[y][x] = f1 + f2;
     716            stamps->window1->kernel[y][x] = f1;
    705717            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(size)) {
    706                 sum += stamps->window->kernel[y][x];
    707             }
    708             maxValue = PS_MAX(maxValue, stamps->window->kernel[y][x]);
    709         }
    710     }
    711 
    712     psTrace("psModules.imcombine", 3, "Window total: %f, threshold: %f\n",
    713             sum, (1.0 - stamps->normFrac) * sum);
    714     bool done = false;
    715     for (int radius = 1; radius <= size && !done; radius++) {
    716         double within = 0.0;
     718                sum1 += stamps->window1->kernel[y][x];
     719            }
     720            maxValue1 = PS_MAX(maxValue1, stamps->window1->kernel[y][x]);
     721
     722            stamps->window2->kernel[y][x] = f2;
     723            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(size)) {
     724                sum2 += stamps->window2->kernel[y][x];
     725            }
     726            maxValue2 = PS_MAX(maxValue2, stamps->window2->kernel[y][x]);
     727        }
     728    }
     729
     730    psTrace("psModules.imcombine", 3, "Window total (1): %f, threshold: %f\n", sum1, (1.0 - stamps->normFrac) * sum1);
     731    psTrace("psModules.imcombine", 3, "Window total (2): %f, threshold: %f\n", sum2, (1.0 - stamps->normFrac) * sum2);
     732
     733# if (0)
     734    // this block attempts to calculate the radius based on the first radial moment
     735    bool done1 = false;
     736    bool done2 = false;
     737    double prior1 = 0.0;
     738    double prior2 = 0.0;
     739    for (int y = -size; y <= size; y++) {
     740        for (int x = -size; x <= size; x++) {
     741            float r = hypot(x, y);
     742            Sr1 += r * stamps->window1->kernel[y][x];
     743            Sr2 += r * stamps->window2->kernel[y][x];
     744            Sf1 += stamps->window1->kernel[y][x];
     745            Sf2 += stamps->window2->kernel[y][x];
     746        }
     747    }
     748   
     749    float R1 = Sr1 / Sf1;
     750    float R2 = Sr2 / Sf2;
     751
     752    stamps->normWindow1 = 2.5*R1;
     753    stamps->normWindow1 = 2.5*R2;
     754# else
     755    // XXX : this block attempts to calculate the radius by looking at the curve of growth (or something vaguely equivalent).
     756    // It did not do very well (though a true curve-of-growth analysis might be better...)
     757    bool done1 = false;
     758    bool done2 = false;
     759    double prior1 = 0.0;
     760    double prior2 = 0.0;
     761    double delta1o = 1.0;
     762    double delta2o = 1.0;
     763    for (int radius = 1; radius <= size && !(done1 && done2); radius++) {
     764        double within1 = 0.0;
     765        double within2 = 0.0;
    717766        for (int y = -radius; y <= radius; y++) {
    718767            for (int x = -radius; x <= radius; x++) {
    719768                if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) {
    720                     within += stamps->window->kernel[y][x];
     769                    within1 += stamps->window1->kernel[y][x];
    721770                }
    722             }
    723         }
    724         psTrace("psModules.imcombine", 5, "Radius %d: %f\n", radius, within);
    725         if (within > (1.0 - stamps->normFrac) * sum) {
    726             stamps->normWindow = radius;
    727             done = true;
    728         }
    729     }
    730 
    731     psTrace("psModules.imcombine", 3, "Normalisation window radius set to %d\n", stamps->normWindow);
    732     if (stamps->normWindow == 0 || stamps->normWindow >= size) {
    733         psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size.");
     771                if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) {
     772                    within2 += stamps->window2->kernel[y][x];
     773                }
     774            }
     775        }
     776        double delta1 = (within1 - prior1) / within1;
     777        if (!done1 && (fabs(delta1) < stamps->normFrac)) {
     778            // interpolate to the radius at which delta2 is normFrac:
     779            stamps->normWindow1 = radius - (stamps->normFrac - delta1) / (delta1o - delta1);
     780            fprintf (stderr, "choosing %f (%d) for 1 (%f : %f)\n", stamps->normWindow1, radius, within1, delta1);
     781            done1 = true;
     782        }
     783        double delta2 = (within2 - prior2) / within2;
     784        if (!done2 && (fabs(delta2) < stamps->normFrac)) {
     785            // interpolate to the radius at which delta2 is normFrac:
     786            stamps->normWindow2 = radius - (stamps->normFrac - delta2) / (delta2o - delta2);
     787            fprintf (stderr, "choosing %f (%d) for 2 (%f : %f)\n", stamps->normWindow2, radius, within2, delta2);
     788            done2 = true;
     789        }
     790        psTrace("psModules.imcombine", 5, "Radius %d: %f (%f) and %f (%f)\n", radius, within1, delta1, within2, delta2);
     791
     792        prior1 = within1;
     793        prior2 = within2;
     794        delta1o = delta1;
     795        delta2o = delta2;
     796
     797        // if (!done1 && (within1 > (1.0 - stamps->normFrac) * sum1)) {
     798        //     stamps->normWindow1 = radius;
     799        //     done1 = true;
     800        // }
     801        // if (!done2 && (within2 > (1.0 - stamps->normFrac) * sum2)) {
     802        //     stamps->normWindow2 = radius;
     803        //     done2 = true;
     804        // }
     805
     806    }
     807# endif
     808
     809    psTrace("psModules.imcombine", 3, "Normalisation window radii set to %f and %f\n", stamps->normWindow1, stamps->normWindow2);
     810    if (stamps->normWindow1 == 0 || stamps->normWindow1 >= size) {
     811        psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (1).");
     812        return false;
     813    }
     814    if (stamps->normWindow2 == 0 || stamps->normWindow2 >= size) {
     815        psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (2).");
    734816        return false;
    735817    }
    … …  
    738820    for (int y = -size; y <= size; y++) {
    739821        for (int x = -size; x <= size; x++) {
    740             stamps->window->kernel[y][x] /= maxValue;
    741         }
    742     }
    743 
    744 #if 0
     822            stamps->window1->kernel[y][x] /= maxValue1;
     823        }
     824    }
     825    // re-normalize so chisquare values are sensible
     826    for (int y = -size; y <= size; y++) {
     827        for (int x = -size; x <= size; x++) {
     828            stamps->window2->kernel[y][x] /= maxValue2;
     829        }
     830    }
     831
     832#if 1
    745833    {
    746         psFits *fits = psFitsOpen ("window.fits", "w");
    747         psFitsWriteImage (fits, NULL, stamps->window->image, 0, NULL);
     834        psFits *fits = NULL;
     835        fits = psFitsOpen ("window1.fits", "w");
     836        psFitsWriteImage (fits, NULL, stamps->window1->image, 0, NULL);
     837        psFitsClose (fits);
     838        fits = psFitsOpen ("window2.fits", "w");
     839        psFitsWriteImage (fits, NULL, stamps->window2->image, 0, NULL);
    748840        psFitsClose (fits);
    749841    }
    … …  
    752844    psFree (stats);
    753845    psFree (flux1);
    754     psFree(flux2);
     846    psFree (flux2);
    755847    psFree (norm1);
    756848    psFree (norm2);
    757849    return true;
     850}
     851
     852static pthread_mutex_t getPenaltiesMutex = PTHREAD_MUTEX_INITIALIZER;
     853
     854// kernels->penalty is an overall scaling factor (user-supplied)
     855bool pmSubtractionKernelPenaltiesStamp(pmSubtractionStamp *stamp, pmSubtractionKernels *kernels)
     856{
     857    // we only need the penalties if we are doing dual convolution
     858    if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) return true;
     859
     860    // we only calculate the penalties once.
     861    if (kernels->havePenalties) return true;
     862
     863    // in a threaded context, only one thread can calculate the penalties.  attempt to grab a
     864    // mutex before continuing
     865    pthread_mutex_lock(&getPenaltiesMutex);
     866
     867    // did someone else already get the mutex and do this?
     868    if (kernels->havePenalties) {
     869        pthread_mutex_unlock(&getPenaltiesMutex);
     870        return true;
     871    }
     872
     873    for (int i = 0; i < kernels->num; i++) {
     874        pmSubtractionKernelPenalties(stamp, kernels, i);
     875    }
     876
     877    kernels->havePenalties = true;
     878    pthread_mutex_unlock(&getPenaltiesMutex);
     879    return true;
     880}
     881
     882# define EMPIRICAL 0
     883
     884// kernels->penalty is an overall scaling factor (user-supplied)
     885bool pmSubtractionKernelPenalties(pmSubtractionStamp *stamp, pmSubtractionKernels *kernels, int index)
     886{
     887    float penalty1, penalty2;
     888    float fwhm1, fwhm2;
     889
     890    // XXX this is annoyingly hack-ish
     891    pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
     892
     893    if (EMPIRICAL) {
     894        psKernel *convolution1 = stamp->convolutions1->data[index];
     895        penalty1 = pmSubtractionKernelPenaltySingle(convolution1);
     896
     897        psKernel *convolution2 = stamp->convolutions2->data[index];
     898        penalty2 = pmSubtractionKernelPenaltySingle(convolution2);
     899    } else {
     900        pmSubtractionKernelPreCalc *kernel = kernels->preCalc->data[index];
     901        float M2 = pmSubtractionKernelPenaltySingle(kernel->kernel);
     902
     903        penalty1 = M2 * PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment
     904        penalty2 = M2 * PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment
     905    }
     906    kernels->penalties1->data.F32[index] = kernels->penalty * penalty1;
     907    psAssert (isfinite(kernels->penalties1->data.F32[index]), "invalid penalty");
     908
     909    kernels->penalties2->data.F32[index] = kernels->penalty * penalty2;
     910    psAssert (isfinite(kernels->penalties2->data.F32[index]), "invalid penalty");
     911
     912    fprintf(stderr, "penalty1: %f, penalty2: %f\n", penalty1, penalty2);
     913
     914    return true;
     915}
     916
     917float pmSubtractionKernelPenaltySingle(psKernel *kernel)
     918{
     919    // Calculate moments
     920    double penalty = 0.0;                   // Moment, for penalty
     921    double sum = 0.0, sum2 = 0.0;           // Sum of kernel component
     922    float min = INFINITY, max = -INFINITY;  // Minimum and maximum kernel value
     923    for (int v = kernel->yMin; v <= kernel->yMax; v++) {
     924        for (int u = kernel->xMin; u <= kernel->xMax; u++) {
     925            double value = kernel->kernel[v][u];
     926            double value2 = PS_SQR(value);
     927            sum += value;
     928            sum2 += value2;
     929            penalty += value2 * PS_SQR((PS_SQR(u) + PS_SQR(v)));
     930            min = PS_MIN(value, min);
     931            max = PS_MAX(value, max);
     932        }
     933    }
     934    penalty *= 1.0 / sum2;
     935
     936    if (0) {
     937        fprintf(stderr, "min: %lf, max: %lf, moment: %lf\n", min, max, penalty);
     938        // psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, fwhm, uOrder, vOrder, penalty);
     939    }
     940
     941    return penalty;
    758942}
    759943

  • branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionStamps.h

    r26893 r28866  
    2424    psArray *flux;                      ///< Fluxes for possible stamps (or NULL)
    2525    int footprint;                      ///< Half-size of stamps
    26     psKernel *window;                   ///< window function generated from ensemble of stamps
    2726    float normFrac;                     ///< Fraction of flux in window for normalisation window
    28     int normWindow;                     ///< Size of window for measuring normalisation
     27    psKernel *window1;                  ///< window function generated from ensemble of stamps (input 1)
     28    psKernel *window2;                  ///< window function generated from ensemble of stamps (input 2)
     29    float normWindow1;                    ///< Size of window for measuring normalisation
     30    float normWindow2;                    ///< Size of window for measuring normalisation
    2931    float sysErr;                       ///< Systematic error
    3032    float skyErr;                       ///< increase effective readnoise
    … …  
    195197bool pmSubtractionStampsResetStatus (pmSubtractionStampList *stamps);
    196198
     199
     200bool pmSubtractionKernelPenaltiesStamp(pmSubtractionStamp *stamp, pmSubtractionKernels *kernels);
     201bool pmSubtractionKernelPenalties(pmSubtractionStamp *stamp, pmSubtractionKernels *kernels, int index);
     202float pmSubtractionKernelPenaltySingle(psKernel *kernel);
     203
    197204#endif

  • branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionVisual.c

    r26893 r28866  
    2121
    2222#include "pmVisual.h"
     23#include "pmVisualUtils.h"
    2324
    2425#include "pmHDU.h"
    … …  
    6162 *    @return true for success */
    6263bool pmSubtractionVisualPlotConvKernels(psImage *convKernels) {
    63     if (!pmVisualIsVisual() || !plotConvKernels) return true;
     64
     65    if (!pmVisualTestLevel("ppsub.kernels", 1)) return true;
     66
     67    if (!plotConvKernels) return true;
     68
    6469    if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) {
    6570        return false;
    … …  
    7580    @return true for success */
    7681bool pmSubtractionVisualPlotStamps(pmSubtractionStampList *stamps, pmReadout *ro) {
    77     if (!pmVisualIsVisual() || !plotStamps) return true;
     82
     83    if (!pmVisualTestLevel("ppsub.stamps", 1)) return true;
     84
     85    if (!plotStamps) return true;
     86
    7887    if (!pmVisualInitWindow (&kapa1, "ppSub:Images")) {
    7988        return false;
    … …  
    145154bool pmSubtractionVisualPlotLeastSquares (pmSubtractionStampList *stamps, bool dual) {
    146155
    147     if (!pmVisualIsVisual() || !plotLeastSquares) return true;
     156    if (!pmVisualTestLevel("ppsub.chisq", 1)) return true;
     157
     158    if (!plotLeastSquares) return true;
     159
    148160    if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {
    149161        return false;
    … …  
    204216
    205217bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub) {
    206     if (!pmVisualIsVisual() || !plotImage) return true;
     218
     219    if (!pmVisualTestLevel("ppsub.images.sub", 1)) return true;
     220
     221    if (!plotImage) return true;
     222
    207223    if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {
    208224        return false;
    … …  
    218234bool pmSubtractionVisualShowKernels(pmSubtractionKernels *kernels) {
    219235
    220     if (!pmVisualIsVisual()) return true;
     236    if (!pmVisualTestLevel("ppsub.kern.final", 1)) return true;
     237
    221238    if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {
    222239        return false;
    … …  
    264281bool pmSubtractionVisualShowBasis(pmSubtractionStampList *stamps) {
    265282
    266     if (!pmVisualIsVisual()) return true;
     283    if (!pmVisualTestLevel("ppsub.basis", 1)) return true;
     284
    267285    if (!pmVisualInitWindow (&kapa2, "ppSub:StampMasterImage")) {
    268286        return false;
    … …  
    425443bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) {
    426444
    427     if (!pmVisualIsVisual()) return true;
     445    if (!pmVisualTestLevel("ppsub.fit", 1)) return true;
    428446
    429447    // generate 4 storage images large enough to hold the N stamps:
    … …  
    462480bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) {
    463481
    464     if (!pmVisualIsVisual()) return true;
     482    if (!pmVisualTestLevel("ppsub.stamp", 1)) return true;
    465483
    466484    double sum;
    … …  
    543561    }
    544562
    545     if (!pmVisualIsVisual()) return true;
     563    if (!pmVisualTestLevel("ppsub.fit", 1)) return true;
     564
    546565    if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) return false;
    547566    if (!pmVisualInitWindow(&kapa2, "ppSub:Misc")) return false;
    … …  
    605624    Graphdata graphdata;
    606625
    607     if (!pmVisualIsVisual()) return true;
     626    if (!pmVisualTestLevel("ppsub.fit", 1)) return true;
     627
    608628    if (!pmVisualInitWindow(&kapa3, "ppSub:plots")) return false;
    609629
Note: See TracChangeset for help on using the changeset viewer.