IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 20, 2010, 1:14:11 PM (16 years ago)
Author:
eugene
Message:

Many changes:

  • psphot-related issues

1) added a new feature to psLib/src/fft/psImageFFT to allow for a pre-FFTed kernel generated for a specific image dimension which can then be applied to an arbitrary number of images for convolution -- this avoids the kernel FFT for every image convolution.
2) updated the psMinimizeLMM to include 2 levels of tolerance: minTol and maxTol: if the minimization reaches minTol, it stops. if it only reaches maxTol within maxIter, it stops and is considered successful, if it fails to reach maxTol, then it fails. This allows us to accept fits that are actually acceptable (within error) even if they are not as close as ideally possible.
3) add new stat: psfWeightNotPoor (vs psfWeightNotBad) : the first gives the fraction of psf-weighted unmasked pixels considering any mask bits (except the internal 'mark' bit), while the former considers only 'bad' mask bits -- these are written to QF_PSF and QF_PSF_PERFECT in the CMF files.
4) define user-set parameters for max and min valid flux in the linear fit analysis: Note this was the cause of the non-negative fluxes in forced photometry -- the min limit was hard-wired to 0.0.
5) significance image is now constructed as (image + (image/1000)2) / sqrt(variance) so that bright sources have well-defined peaks (other wise, they become somewhat flat-topped as image = sqrt(variance), leading to ill-defined peaks).
6) modification of the visualization functions to accept facility and level values akin to psTrace
7) modification of the source fitting APIs to allow thread- and source- independent fit options (iterations, tolerances, etc)
8) use Kron magnitude as test for source size (CR vs EXT) instead of PSF-based aperture
9) set a min systematic error in the aperture mags when used for size classification significance
10) threaded the psf-convolved model (PCM) fitting process (extended source fits)
11) for extended sources, adjust the radius based on the footprint and re-calculate moments for guess; save the psf-based moments for output in psf table
12) for Sersic fitting, choose the best index with a grid search using only a few iterations; iterate fully on the selected value.
13) same for Sersic fitting using the PCM fitting process
14) fixed a bug in which the PSF candidate stars which failed in the psf fitting were not correctly unmarked as being PSF stars
15) define galaxy fit radius based on sky stdev (global, not local)
16) move the PCM code to psModules
17) save and report the raw aperture mag in addition to the curve-of-growth corrected one (PS1_V3)
18) save and report the Kron parameters and higher-order moments
19) some psModule header reorganization for code clarity
20) fixed a bug in which the diff stats were counting 'marked' pixels (outside area of interest) as 'masked'.
21) save the radial profile aperture sizes in the headers
22) better guess for Sersic parameters based on moments (requires setting the functional form so that the scale length is right)

  • ppSub-related:

1) ensure masked pixels are NANed in output diff image
2) add code to flag detections if they have bright positive neighbors (+ output of these in PS1_DV2)
3) define separate penalties for each image (based on their fwhm values) (this requires the penalties to be calculated later in the code).
4) define separate apertures for each image for flux normalization
5) choose aperture based on curve-of-growth (was based on fixed fraction of full aperture flux, and thus noisy)
6) some fine tuning of the penalty factor (this still seems arbitrary, and results are somewhat sensitive to the right value)

Location:
trunk/psModules
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules

    • Property svn:mergeinfo deleted

  • trunk/psModules/src/imcombine/pmPSFEnvelope.c

    r28332 r29004  
    1010#include "pmHDU.h"
    1111#include "pmFPA.h"
    12 #include "pmReadoutFake.h"
    13 
    14 #include "pmMoments.h"
     12
     13#include "pmTrend2D.h"
    1514#include "pmResiduals.h"
    1615#include "pmGrowthCurve.h"
    17 #include "pmTrend2D.h"
    18 #include "pmPSF.h"
    19 #include "pmModel.h"
    20 #include "pmModelClass.h"
    21 #include "pmModelUtils.h"
    2216#include "pmSpan.h"
     17#include "pmFootprintSpans.h"
    2318#include "pmFootprint.h"
    2419#include "pmPeaks.h"
     20#include "pmMoments.h"
     21#include "pmModelFuncs.h"
     22#include "pmModel.h"
     23#include "pmModelUtils.h"
     24#include "pmModelClass.h"
     25#include "pmSourceMasks.h"
     26#include "pmSourceExtendedPars.h"
     27#include "pmSourceDiffStats.h"
    2528#include "pmSource.h"
    26 #include "pmSourceUtils.h"
    2729#include "pmSourceFitModel.h"
     30#include "pmPSF.h"
    2831#include "pmPSFtry.h"
    2932
    30 
     33#include "pmReadoutFake.h"
    3134#include "pmPSFEnvelope.h"
    3235#include "pmStackVisual.h"
    … …  
    410413    options->chiFluxTrend = false;      // All sources have similar flux, so fitting a trend often fails
    411414
    412     pmSourceFitModelInit(SOURCE_FIT_ITERATIONS, 0.01, VARIANCE_VAL, true);
     415    // options which modify the behavior of the model fitting
     416    options->fitOptions                = pmSourceFitOptionsAlloc();
     417    options->fitOptions->nIter         = SOURCE_FIT_ITERATIONS;
     418    options->fitOptions->minTol        = 0.01;
     419    options->fitOptions->maxTol        = 1.00;
     420    options->fitOptions->poissonErrors = true;
     421    options->fitOptions->weight        = VARIANCE_VAL;
     422    options->fitOptions->mode          = PM_SOURCE_FIT_PSF;
     423
    413424    pmModelClassSetLimits(PM_MODEL_LIMITS_STRICT); // Important for getting a good stack target PSF
    414425

  • trunk/psModules/src/imcombine/pmSubtraction.c

    r28667 r29004  
    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:
    … …  
    12151223    bool threaded = pmSubtractionThreaded(); // Running threaded?
    12161224
     1225    // XXX This is no longer used
    12171226    psImage *convMask = NULL;           // Convolved mask image (common to inputs 1 and 2)
    12181227    if (subMask) {
    … …  
    14131422    return true;
    14141423}
     1424
     1425bool pmSubtractionGetFWHMs(float *fwhm1, float *fwhm2) {
     1426
     1427  *fwhm1 = FWHM1;
     1428  *fwhm2 = FWHM2;
     1429  return true;
     1430}
     1431
     1432bool pmSubtractionSetFWHMs(float fwhm1, float fwhm2) {
     1433
     1434  FWHM1 = fwhm1;
     1435  FWHM2 = fwhm2;
     1436  return true;
     1437}

  • trunk/psModules/src/imcombine/pmSubtraction.h

    r26893 r29004  
    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

  • trunk/psModules/src/imcombine/pmSubtractionEquation.c

    r28405 r29004  
    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]);

  • trunk/psModules/src/imcombine/pmSubtractionEquation.h

    r26893 r29004  
    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    );

  • trunk/psModules/src/imcombine/pmSubtractionKernels.c

    r27335 r29004  
    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
    … …  
    207218        if (uOrder % 2 == 0 && vOrder % 2 == 0) {
    208219            // Even functions: normalise to unit sum and subtract null pixel so that sum is zero
    209             scale2D = 1.0 / fabs(sum);
     220            // Re-normalize
     221            // scale2D  = 1.0 / fabs(sum);
     222            scale2D  = 1.0 / sqrt(sum2);
    210223            zeroNull = true;
    211224        } else {
    … …  
    239252
    240253    psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "*", psScalarAlloc(scale2D, PS_TYPE_F32));
    241     penalty *= 1.0 / sum2;
    242254
    243255    if (zeroNull) {
    244         preCalc->kernel->kernel[0][0] -= 1.0;
    245     }
    246 
    247 #if 0
     256        // preCalc->kernel->kernel[0][0] -= 1.0;
     257        preCalc->kernel->kernel[0][0] -= sum / sqrt (sum2);
     258    }
     259
     260#if 1
    248261    {
    249         double sum = 0.0;   // Sum of kernel component
     262        double Sum = 0.0;   // Sum of kernel component
     263        double Sum2 = 0.0;   // Sum of kernel component
    250264        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++) {
    253                 sum += preCalc->kernel->kernel[v][u];
     265        for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
     266            for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) {
     267                double value = preCalc->kernel->kernel[v][u];
     268                Sum += value;
     269                Sum2 += PS_SQR(value);
    254270                min = PS_MIN(preCalc->kernel->kernel[v][u], min);
    255271                max = PS_MAX(preCalc->kernel->kernel[v][u], max);
    256272            }
    257273        }
    258         fprintf(stderr, "%d mod: %lf, null: %f, min: %lf, max: %lf, scale: %f\n", index, sum, preCalc->kernel->kernel[0][0], min, max, scale2D);
     274        fprintf(stderr, "%d sum: %lf, sum2: %lf, null: %f, min: %lf, max: %lf, scale: %f\n", index, Sum, Sum2, preCalc->kernel->kernel[0][0], min, max, scale2D);
    259275    }
    260276#endif
    … …  
    267283    }
    268284    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);
     285    psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, fwhm, uOrder, vOrder);
    272286
    273287    return true;
    … …  
    321335
    322336                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);
     337                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
     338                // pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], false, false);
    325339            }
    326340        }
    … …  
    379393            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    380394                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);
     395                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
    382396            }
    383397        }
    … …  
    385399            // XXX modify size for hermitians to account for sqrt(2) in Hermitian definition (relative to ISIS Gaussian)
    386400            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);
     401            pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, order, order, fwhms->data.F32[i], true, true);
    388402        }
    389403    }
    … …  
    437451            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    438452                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);
     453                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
    440454            }
    441455        }
    … …  
    506520
    507521                // 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);
     522                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
    509523
    510524                // XXXX test demo that deconvolved kernel is valid
    … …  
    572586    kernels->preCalc = psArrayAlloc(numBasisFunctions);
    573587    kernels->penalty = penalty;
    574     kernels->penalties = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
     588    kernels->penalties1 = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
     589    psVectorInit(kernels->penalties1, NAN);
     590    kernels->penalties2 = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
     591    psVectorInit(kernels->penalties2, NAN);
     592    kernels->havePenalties = false;
    575593    kernels->size = size;
    576594    kernels->inner = 0;
    … …  
    771789
    772790    psWarning("Kernel penalty for dual-convolution is not configured for SPAM kernels.");
    773     psVectorInit(kernels->penalties, 0.0);
     791    psVectorInit(kernels->penalties1, 0.0);
     792    psVectorInit(kernels->penalties2, 0.0);
    774793
    775794    return kernels;
    … …  
    866885
    867886    psWarning("Kernel penalty for dual-convolution is not configured for FRIES kernels.");
    868     psVectorInit(kernels->penalties, 0.0);
     887    psVectorInit(kernels->penalties1, 0.0);
     888    psVectorInit(kernels->penalties2, 0.0);
    869889
    870890    return kernels;
    … …  
    10401060                kernels->u->data.S32[index] = uOrder;
    10411061                kernels->v->data.S32[index] = vOrder;
    1042                 kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
    1043                 if (!isfinite(kernels->penalties->data.F32[index])) {
     1062
     1063                // XXX convert to use the convolved 2nd moment
     1064                kernels->penalties1->data.F32[index] = kernels->penalty * fabsf(moment);
     1065                if (!isfinite(kernels->penalties1->data.F32[index])) {
     1066                    psAbort ("invalid penalty");
     1067                }
     1068                kernels->penalties2->data.F32[index] = kernels->penalty * fabsf(moment);
     1069                if (!isfinite(kernels->penalties2->data.F32[index])) {
    10441070                    psAbort ("invalid penalty");
    10451071                }
    … …  
    12471273    out->preCalc = psMemIncrRefCounter(in->preCalc);
    12481274    out->penalty = in->penalty;
    1249     out->penalties = psMemIncrRefCounter(in->penalties);
     1275    out->penalties1 = psMemIncrRefCounter(in->penalties1);
     1276    out->penalties2 = psMemIncrRefCounter(in->penalties2);
    12501277    out->uStop = psMemIncrRefCounter(in->uStop);
    12511278    out->vStop = psMemIncrRefCounter(in->vStop);

  • trunk/psModules/src/imcombine/pmSubtractionKernels.h

    r26893 r29004  
    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

  • trunk/psModules/src/imcombine/pmSubtractionMatch.c

    r28405 r29004  
    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;

  • trunk/psModules/src/imcombine/pmSubtractionMatch.h

    r26893 r29004  
    110110    );
    111111
    112 
    113112#endif

  • trunk/psModules/src/imcombine/pmSubtractionStamps.c

    r27402 r29004  
    77#include <pslib.h>
    88
     9#include "pmErrorCodes.h"
     10#include "pmHDU.h"
     11#include "pmFPA.h"
     12
    913// All these includes required to get stamps out of an array of pmSources
    10 #include "pmMoments.h"
     14#include "pmTrend2D.h"
     15#include "pmResiduals.h"
     16#include "pmGrowthCurve.h"
    1117#include "pmSpan.h"
     18#include "pmFootprintSpans.h"
    1219#include "pmFootprint.h"
    1320#include "pmPeaks.h"
    14 #include "pmResiduals.h"
    15 #include "pmHDU.h"
    16 #include "pmFPA.h"
    17 #include "pmGrowthCurve.h"
    18 #include "pmTrend2D.h"
    19 #include "pmPSF.h"
     21#include "pmMoments.h"
     22#include "pmModelFuncs.h"
    2023#include "pmModel.h"
     24#include "pmSourceMasks.h"
     25#include "pmSourceExtendedPars.h"
     26#include "pmSourceDiffStats.h"
    2127#include "pmSource.h"
    22 #include "pmErrorCodes.h"
    2328
    2429#include "pmSubtraction.h"
    … …  
    4651    psFree(list->y);
    4752    psFree(list->flux);
    48     psFree(list->window);
     53    psFree(list->window1);
     54    psFree(list->window2);
    4955}
    5056
    … …  
    225231    list->y = NULL;
    226232    list->flux = NULL;
    227     list->window = NULL;
    228233    list->normFrac = normFrac;
    229     list->normWindow = 0;
     234    list->window1 = NULL;
     235    list->window2 = NULL;
     236    list->normWindow1 = 0;
     237    list->normWindow2 = 0;
    230238    list->footprint = footprint;
    231239    list->sysErr = sysErr;
    … …  
    248256    out->y = NULL;
    249257    out->flux = NULL;
    250     out->window = psMemIncrRefCounter(in->window);
     258    out->window1 = psMemIncrRefCounter(in->window1);
     259    out->window2 = psMemIncrRefCounter(in->window2);
    251260    out->footprint = in->footprint;
    252     out->normWindow = in->normWindow;
     261    out->normWindow1 = in->normWindow1;
     262    out->normWindow2 = in->normWindow2;
    253263
    254264    for (int i = 0; i < num; i++) {
    … …  
    638648    int size = stamps->footprint; // Size of postage stamps
    639649
    640     psFree (stamps->window);
    641     stamps->window = psKernelAlloc(-size, size, -size, size);
    642     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);
    643657
    644658    // generate normalizations for each stamp
    … …  
    669683
    670684    // generate the window pixels
    671     double sum = 0.0;                   // Sum inside the window
    672     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
    673689    for (int y = -size; y <= size; y++) {
    674690        for (int x = -size; x <= size; x++) {
    … …  
    691707            }
    692708            float f1 = stats->robustMedian;
     709
    693710            psStatsInit (stats);
    694711            if (!psVectorStats (stats, flux2, NULL, NULL, 0)) {
    … …  
    697714            float f2 = stats->robustMedian;
    698715
    699             stamps->window->kernel[y][x] = f1 + f2;
     716            stamps->window1->kernel[y][x] = f1;
    700717            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(size)) {
    701                 sum += stamps->window->kernel[y][x];
    702             }
    703             maxValue = PS_MAX(maxValue, stamps->window->kernel[y][x]);
    704         }
    705     }
    706 
    707     psTrace("psModules.imcombine", 3, "Window total: %f, threshold: %f\n",
    708             sum, (1.0 - stamps->normFrac) * sum);
    709     bool done = false;
    710     for (int radius = 1; radius <= size && !done; radius++) {
    711         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;
    712766        for (int y = -radius; y <= radius; y++) {
    713767            for (int x = -radius; x <= radius; x++) {
    714768                if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) {
    715                     within += stamps->window->kernel[y][x];
     769                    within1 += stamps->window1->kernel[y][x];
    716770                }
    717             }
    718         }
    719         psTrace("psModules.imcombine", 5, "Radius %d: %f\n", radius, within);
    720         if (within > (1.0 - stamps->normFrac) * sum) {
    721             stamps->normWindow = radius;
    722             done = true;
    723         }
    724     }
    725 
    726     psTrace("psModules.imcombine", 3, "Normalisation window radius set to %d\n", stamps->normWindow);
    727     if (stamps->normWindow == 0 || stamps->normWindow >= size) {
    728         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).");
    729816        return false;
    730817    }
    … …  
    733820    for (int y = -size; y <= size; y++) {
    734821        for (int x = -size; x <= size; x++) {
    735             stamps->window->kernel[y][x] /= maxValue;
    736         }
    737     }
    738 
    739 #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
    740833    {
    741         psFits *fits = psFitsOpen ("window.fits", "w");
    742         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);
    743840        psFitsClose (fits);
    744841    }
    … …  
    747844    psFree (stats);
    748845    psFree (flux1);
    749     psFree(flux2);
     846    psFree (flux2);
    750847    psFree (norm1);
    751848    psFree (norm2);
    752849    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    bool zeroNull = false;
     894    int uOrder = kernels->u->data.S32[index];
     895    int vOrder = kernels->v->data.S32[index];
     896    if (uOrder % 2 == 0 && vOrder % 2 == 0) zeroNull = true;
     897
     898    if (EMPIRICAL) {
     899        psKernel *convolution1 = stamp->convolutions1->data[index];
     900        penalty1 = pmSubtractionKernelPenaltySingle(convolution1, zeroNull);
     901
     902        psKernel *convolution2 = stamp->convolutions2->data[index];
     903        penalty2 = pmSubtractionKernelPenaltySingle(convolution2, zeroNull);
     904    } else {
     905        pmSubtractionKernelPreCalc *kernel = kernels->preCalc->data[index];
     906        float M2 = pmSubtractionKernelPenaltySingle(kernel->kernel, zeroNull);
     907
     908        penalty1 = M2 + PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment
     909        penalty2 = M2 + PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment
     910    }
     911    kernels->penalties1->data.F32[index] = kernels->penalty * penalty1;
     912    psAssert (isfinite(kernels->penalties1->data.F32[index]), "invalid penalty");
     913
     914    kernels->penalties2->data.F32[index] = kernels->penalty * penalty2;
     915    psAssert (isfinite(kernels->penalties2->data.F32[index]), "invalid penalty");
     916
     917    fprintf(stderr, "penalty1: %f, penalty2: %f\n", penalty1, penalty2);
     918
     919    return true;
     920}
     921
     922float pmSubtractionKernelPenaltySingle(psKernel *kernel, bool zeroNull)
     923{
     924    // Calculate moments
     925    double penalty = 0.0;                   // Moment, for penalty
     926    double sum = 0.0, sum2 = 0.0;           // Sum of kernel component
     927    float min = INFINITY, max = -INFINITY;  // Minimum and maximum kernel value
     928    for (int v = kernel->yMin; v <= kernel->yMax; v++) {
     929        for (int u = kernel->xMin; u <= kernel->xMax; u++) {
     930            double value = kernel->kernel[v][u];
     931            if (false && zeroNull && (u == 0) && (v == 0)) {
     932                value += 1.0;
     933            }
     934            double value2 = PS_SQR(value);
     935            sum += value;
     936            sum2 += value2;
     937            penalty += value2 * PS_SQR((PS_SQR(u) + PS_SQR(v)));
     938            min = PS_MIN(value, min);
     939            max = PS_MAX(value, max);
     940        }
     941    }
     942    penalty *= 1.0 / sum2;
     943
     944    if (1) {
     945        fprintf(stderr, "min: %lf, max: %lf, moment: %lf, flux^2: %lf\n", min, max, penalty, sum2);
     946        // psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, fwhm, uOrder, vOrder, penalty);
     947    }
     948
     949    return penalty;
    753950}
    754951

  • trunk/psModules/src/imcombine/pmSubtractionStamps.h

    r26893 r29004  
    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, bool zeroNull);
     203
    197204#endif

  • trunk/psModules/src/imcombine/pmSubtractionVisual.c

    r26893 r29004  
    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.kernels.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.