IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 29004 for trunk/psModules


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:
1 deleted
104 edited
13 copied

Legend:

Unmodified
Added
Removed

  • trunk/psModules

    • Property svn:mergeinfo deleted

  • trunk/psModules/src/camera/pmFPAfileIO.c

    r28340 r29004  
    2424#include "pmFPAfileFitsIO.h"
    2525#include "pmFPAfileFringeIO.h"
     26
     27#include "pmTrend2D.h"
     28#include "pmResiduals.h"
     29#include "pmGrowthCurve.h"
    2630#include "pmSpan.h"
     31#include "pmFootprintSpans.h"
    2732#include "pmFootprint.h"
    2833#include "pmPeaks.h"
    2934#include "pmMoments.h"
    30 #include "pmResiduals.h"
    31 #include "pmGrowthCurve.h"
    32 #include "pmTrend2D.h"
     35#include "pmModelFuncs.h"
     36#include "pmModel.h"
     37#include "pmSourceMasks.h"
     38#include "pmSourceExtendedPars.h"
     39#include "pmSourceDiffStats.h"
     40#include "pmSource.h"
     41#include "pmSourceFitModel.h"
    3342#include "pmPSF.h"
    34 #include "pmModel.h"
    35 #include "pmSource.h"
     43#include "pmPSFtry.h"
     44
    3645#include "pmSourceIO.h"
    37 #include "pmResiduals.h"
    3846#include "pmPSF_IO.h"
     47
    3948#include "pmAstrometryModel.h"
    4049#include "pmAstrometryRefstars.h"

  • trunk/psModules/src/camera/pmReadoutFake.c

    r28405 r29004  
    1010#include "pmFPA.h"
    1111
    12 #include "pmMoments.h"
     12
     13#include "pmTrend2D.h"
    1314#include "pmResiduals.h"
    1415#include "pmGrowthCurve.h"
    15 #include "pmTrend2D.h"
    16 #include "pmPSF.h"
    17 #include "pmModel.h"
    18 #include "pmModelClass.h"
    1916#include "pmSpan.h"
     17#include "pmFootprintSpans.h"
    2018#include "pmFootprint.h"
    2119#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"
    2228#include "pmSource.h"
    23 #include "pmSourceUtils.h"
    24 #include "pmModelUtils.h"
     29#include "pmSourceFitModel.h"
     30#include "pmPSF.h"
     31#include "pmPSFtry.h"
     32
    2533#include "pmSourceGroups.h"
    26 
    2734#include "pmReadoutFake.h"
    2835

  • trunk/psModules/src/camera/pmReadoutFake.h

    r26450 r29004  
    22#define PM_READOUT_FAKE_H
    33
    4 #include <pslib.h>
    5 #include <pmHDU.h>
    6 #include <pmFPA.h>
    7 
    8 #include <pmMoments.h>
    9 #include <pmResiduals.h>
    10 #include <pmGrowthCurve.h>
    11 #include <pmTrend2D.h>
    12 #include <pmPSF.h>
    13 #include <pmSourceMasks.h>
     4// #include <pslib.h>
     5// #include <pmHDU.h>
     6// #include <pmFPA.h>
     7//
     8// #include <pmMoments.h>
     9// #include <pmResiduals.h>
     10// #include <pmGrowthCurve.h>
     11// #include <pmTrend2D.h>
     12// #include <pmPSF.h>
     13// #include <pmSourceMasks.h>
    1414
    1515/// Set threading

  • trunk/psModules/src/config/pmConfig.c

    r28287 r29004  
    3232
    3333#include "pmConfig.h"
     34#include "pmVisualUtils.h"
    3435
    3536#ifdef HAVE_NEBCLIENT
    … …  
    638639        psArgumentVerbosity(argc, argv);
    639640        // XXX: substitute the string for the default log level for "2".
     641    }
     642
     643    // Set the visualization levels
     644    // argument format is: -visual (facil) (level)
     645    while ((argNum = psArgumentGet(*argc, argv, "-visual"))) {
     646        if ( (*argc < argNum + 3) ) {
     647            psError(PS_ERR_IO, true, "-visual switch specified without facility and level.");
     648            return NULL;
     649        }
     650        psArgumentRemove(argNum, argc, argv);
     651        pmVisualSetLevel(argv[argNum], atoi(argv[argNum+1]));
     652        psArgumentRemove(argNum, argc, argv);
     653        psArgumentRemove(argNum, argc, argv);
     654    }
     655    if ((argNum = psArgumentGet(*argc, argv, "-visual-all"))) {
     656        pmVisualSetLevel(".", 10);
     657        psArgumentRemove(argNum, argc, argv);
     658    }
     659    if ((argNum = psArgumentGet(*argc, argv, "-visual-levels"))) {
     660        pmVisualPrintLevels(stdout);
     661        psArgumentRemove(argNum, argc, argv);
    640662    }
    641663

  • trunk/psModules/src/detrend/pmMaskStats.c

    r28100 r29004  
    77
    88#include <pslib.h>
    9 #include <psmodules.h>
    109
    1110#include "pmHDU.h"
    … …  
    1615#include "pmFPAAstrometry.h"
    1716
     17#include "pmMaskStats.h"
    1818
    1919#define ESCAPE { \

  • trunk/psModules/src/detrend/pmShutterCorrection.c

    r28405 r29004  
    302302    PS_ASSERT_PTR_NON_NULL(guess, NULL);
    303303
    304     psMinimization *minInfo = psMinimizationAlloc(15, 0.1); // Minimization information
     304    psMinimization *minInfo = psMinimizationAlloc(15, 0.1, 1.0); // Minimization information
    305305
    306306    psVector *params = psVectorAlloc (3, PS_TYPE_F32); // Fitting parameters

  • trunk/psModules/src/extras/Makefile.am

    r27750 r29004  
    99        pmKapaPlots.c \
    1010        pmVisual.c \
     11        pmVisualUtils.c \
    1112        ippStages.c
    1213
    … …  
    1718        pmKapaPlots.h \
    1819        pmVisual.h \
     20        pmVisualUtils.h \
    1921        ippDiffMode.h \
    2022        ippStages.h

  • trunk/psModules/src/extras/pmVisual.c

    r28129 r29004  
    2121#include "pmAstrometryObjects.h"
    2222#include "pmSubtractionStamps.h"
     23
    2324#include "pmTrend2D.h"
     25#include "pmResiduals.h"
     26#include "pmGrowthCurve.h"
     27#include "pmSpan.h"
     28#include "pmFootprintSpans.h"
     29#include "pmFootprint.h"
     30#include "pmPeaks.h"
     31#include "pmMoments.h"
     32#include "pmModelFuncs.h"
     33#include "pmModel.h"
     34#include "pmSourceMasks.h"
     35#include "pmSourceExtendedPars.h"
     36#include "pmSourceDiffStats.h"
     37#include "pmSource.h"
     38#include "pmSourceFitModel.h"
    2439#include "pmPSF.h"
    2540#include "pmPSFtry.h"
    26 #include "pmSource.h"
     41
    2742#include "pmFPAExtent.h"
    28 
    29 # if (HAVE_KAPA)
    30 # include <kapa.h>
    31 # include "pmVisual.h"
    32 # include "pmKapaPlots.h"
    33 # define KAPAX 700
    34 # define KAPAY 700
    3543
    3644#include "pmAstrometryVisual.h"
    … …  
    3846#include "pmStackVisual.h"
    3947#include "pmSourceVisual.h"
     48
     49# if (HAVE_KAPA)
     50# include <kapa.h>
     51#include "pmVisual.h"
     52#include "pmKapaPlots.h"
     53
     54# define KAPAX 700
     55# define KAPAY 700
    4056
    4157//#define TESTING
    … …  
    102118
    103119    /* Wait up to 1.0 second for a response, then continue */
    104     timeout.tv_sec = 2;
     120    timeout.tv_sec = 10;
    105121    timeout.tv_usec = 0;
    106122

  • 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

  • trunk/psModules/src/objects/Makefile.am

    r28013 r29004  
    4343        pmSourceIO_CMF_PS1_V1.c \
    4444        pmSourceIO_CMF_PS1_V2.c \
     45        pmSourceIO_CMF_PS1_V3.c \
    4546        pmSourceIO_CMF_PS1_SV1.c \
    4647        pmSourceIO_CMF_PS1_DV1.c \
     48        pmSourceIO_CMF_PS1_DV2.c \
    4749        pmSourceIO_MatchedRefs.c \
    4850        pmSourcePlots.c \
    … …  
    6062        pmPSFtryFitPSF.c \
    6163        pmPSFtryMetric.c \
     64        pmPCMdata.c \
     65        pmPCM_MinimizeChisq.c \
     66        pmSourceFitPCM.c \
    6267        pmTrend2D.c \
    6368        pmGrowthCurveGenerate.c \
    … …  
    7176        models/pmModel_QGAUSS.c \
    7277        models/pmModel_RGAUSS.c \
    73         models/pmModel_SERSIC.c
     78        models/pmModel_SERSIC.c \
     79        models/pmModel_EXP.c \
     80        models/pmModel_DEV.c
    7481
    7582pkginclude_HEADERS = \
    … …  
    8087        pmPeaks.h \
    8188        pmMoments.h \
     89        pmModelFuncs.h \
    8290        pmModel.h \
    8391        pmModelClass.h \
    … …  
    101109        pmPSF_IO.h \
    102110        pmPSFtry.h \
     111        pmPCMdata.h \
    103112        pmTrend2D.h \
    104113        pmGrowthCurve.h \
    … …  
    111120        models/pmModel_QGAUSS.h \
    112121        models/pmModel_RGAUSS.h \
    113         models/pmModel_SERSIC.h
     122        models/pmModel_SERSIC.h \
     123        models/pmModel_EXP.h \
     124        models/pmModel_DEV.h
    114125
    115126CLEANFILES = *~

  • trunk/psModules/src/objects/models/pmModel_GAUSS.c

    r26916 r29004  
    2121#include <stdio.h>
    2222#include <pslib.h>
    23 
     23#include "pmHDU.h"
     24#include "pmFPA.h"
     25
     26#include "pmTrend2D.h"
     27#include "pmResiduals.h"
     28#include "pmGrowthCurve.h"
     29#include "pmSpan.h"
     30#include "pmFootprintSpans.h"
     31#include "pmFootprint.h"
     32#include "pmPeaks.h"
    2433#include "pmMoments.h"
    25 #include "pmPeaks.h"
     34#include "pmModelFuncs.h"
     35#include "pmModel.h"
     36#include "pmModelUtils.h"
     37#include "pmModelClass.h"
     38#include "pmSourceMasks.h"
     39#include "pmSourceExtendedPars.h"
     40#include "pmSourceDiffStats.h"
    2641#include "pmSource.h"
    27 #include "pmModel.h"
     42#include "pmSourceFitModel.h"
     43#include "pmPSF.h"
     44#include "pmPSFtry.h"
     45#include "pmDetections.h"
     46
    2847#include "pmModel_GAUSS.h"
    2948
     49# define PM_MODEL_NPARAM          7
    3050# define PM_MODEL_FUNC            pmModelFunc_GAUSS
    3151# define PM_MODEL_FLUX            pmModelFlux_GAUSS
    … …  
    83103        dPAR[PM_PAR_XPOS] = q*(2*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]);
    84104        dPAR[PM_PAR_YPOS] = q*(2*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]);
     105
    85106        // the extra factor of 2 below is needed to avoid excessive swings
    86107        dPAR[PM_PAR_SXX]  = +4.0*q*px*px/PAR[PM_PAR_SXX];
    … …  
    102123        return true;
    103124    }
    104     psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
     125    psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds");
    105126
    106127    // we need to calculate the limits for SXY specially
    … …  
    347368// this test is invalid if the parameters are derived
    348369// from the PSF model
     370// XXX how is this used?  it prevents forced photometry from ever being 'successful'
    349371bool PM_MODEL_FIT_STATUS (pmModel *model)
    350372{
    … …  
    394416    return;
    395417}
    396 
    397 # undef PM_MODEL_FUNC
    398 # undef PM_MODEL_FLUX
    399 # undef PM_MODEL_GUESS
    400 # undef PM_MODEL_LIMITS
    401 # undef PM_MODEL_RADIUS
    402 # undef PM_MODEL_FROM_PSF
    403 # undef PM_MODEL_PARAMS_FROM_PSF
    404 # undef PM_MODEL_FIT_STATUS
    405 # undef PM_MODEL_SET_LIMITS

  • trunk/psModules/src/objects/models/pmModel_PGAUSS.c

    r27565 r29004  
    2121#include <stdio.h>
    2222#include <pslib.h>
    23 
     23#include "pmHDU.h"
     24#include "pmFPA.h"
     25
     26#include "pmTrend2D.h"
     27#include "pmResiduals.h"
     28#include "pmGrowthCurve.h"
     29#include "pmSpan.h"
     30#include "pmFootprintSpans.h"
     31#include "pmFootprint.h"
     32#include "pmPeaks.h"
    2433#include "pmMoments.h"
    25 #include "pmPeaks.h"
     34#include "pmModelFuncs.h"
     35#include "pmModel.h"
     36#include "pmModelUtils.h"
     37#include "pmModelClass.h"
     38#include "pmSourceMasks.h"
     39#include "pmSourceExtendedPars.h"
     40#include "pmSourceDiffStats.h"
    2641#include "pmSource.h"
    27 #include "pmModel.h"
     42#include "pmSourceFitModel.h"
     43#include "pmPSF.h"
     44#include "pmPSFtry.h"
     45#include "pmDetections.h"
     46
    2847#include "pmModel_PGAUSS.h"
    2948
     49# define PM_MODEL_NPARAM          7
    3050# define PM_MODEL_FUNC            pmModelFunc_PGAUSS
    3151# define PM_MODEL_FLUX            pmModelFlux_PGAUSS
    … …  
    103123        return true;
    104124    }
    105     psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
     125    psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds");
    106126
    107127    // we need to calculate the limits for SXY specially
    … …  
    448468    return;
    449469}
    450 
    451 # undef PM_MODEL_FUNC
    452 # undef PM_MODEL_FLUX
    453 # undef PM_MODEL_GUESS
    454 # undef PM_MODEL_LIMITS
    455 # undef PM_MODEL_RADIUS
    456 # undef PM_MODEL_FROM_PSF
    457 # undef PM_MODEL_PARAMS_FROM_PSF
    458 # undef PM_MODEL_FIT_STATUS
    459 # undef PM_MODEL_SET_LIMITS

  • trunk/psModules/src/objects/models/pmModel_PS1_V1.c

    r27565 r29004  
    2222#include <stdio.h>
    2323#include <pslib.h>
    24 
     24#include "pmHDU.h"
     25#include "pmFPA.h"
     26
     27#include "pmTrend2D.h"
     28#include "pmResiduals.h"
     29#include "pmGrowthCurve.h"
     30#include "pmSpan.h"
     31#include "pmFootprintSpans.h"
     32#include "pmFootprint.h"
     33#include "pmPeaks.h"
    2534#include "pmMoments.h"
    26 #include "pmPeaks.h"
     35#include "pmModelFuncs.h"
     36#include "pmModel.h"
     37#include "pmModelUtils.h"
     38#include "pmModelClass.h"
     39#include "pmSourceMasks.h"
     40#include "pmSourceExtendedPars.h"
     41#include "pmSourceDiffStats.h"
    2742#include "pmSource.h"
    28 #include "pmModel.h"
     43#include "pmSourceFitModel.h"
     44#include "pmPSF.h"
     45#include "pmPSFtry.h"
     46#include "pmDetections.h"
     47
    2948#include "pmModel_PS1_V1.h"
    3049
     50# define PM_MODEL_NPARAM          8
    3151# define PM_MODEL_FUNC            pmModelFunc_PS1_V1
    3252# define PM_MODEL_FLUX            pmModelFlux_PS1_V1
    … …  
    122142        return true;
    123143    }
    124     psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
     144    psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds");
    125145
    126146    // we need to calculate the limits for SXY specially
    … …  
    222242    PAR[PM_PAR_SYY]  = PS_MAX(0.5, M_SQRT2*shape.sy);
    223243    PAR[PM_PAR_SXY]  = shape.sxy;
    224     PAR[PM_PAR_7]    = 1.0;
     244    PAR[PM_PAR_7]    = 0.5;
    225245
    226246    return(true);
    … …  
    468488    return;
    469489}
    470 
    471 # undef PM_MODEL_FUNC
    472 # undef PM_MODEL_FLUX
    473 # undef PM_MODEL_GUESS
    474 # undef PM_MODEL_LIMITS
    475 # undef PM_MODEL_RADIUS
    476 # undef PM_MODEL_FROM_PSF
    477 # undef PM_MODEL_PARAMS_FROM_PSF
    478 # undef PM_MODEL_FIT_STATUS
    479 # undef PM_MODEL_SET_LIMITS
    480 # undef ALPHA
    481 # undef ALPHA_M

  • trunk/psModules/src/objects/models/pmModel_QGAUSS.c

    r27565 r29004  
    2222#include <stdio.h>
    2323#include <pslib.h>
    24 
     24#include "pmHDU.h"
     25#include "pmFPA.h"
     26
     27#include "pmTrend2D.h"
     28#include "pmResiduals.h"
     29#include "pmGrowthCurve.h"
     30#include "pmSpan.h"
     31#include "pmFootprintSpans.h"
     32#include "pmFootprint.h"
     33#include "pmPeaks.h"
    2534#include "pmMoments.h"
    26 #include "pmPeaks.h"
     35#include "pmModelFuncs.h"
     36#include "pmModel.h"
     37#include "pmModelUtils.h"
     38#include "pmModelClass.h"
     39#include "pmSourceMasks.h"
     40#include "pmSourceExtendedPars.h"
     41#include "pmSourceDiffStats.h"
    2742#include "pmSource.h"
    28 #include "pmModel.h"
     43#include "pmSourceFitModel.h"
     44#include "pmPSF.h"
     45#include "pmPSFtry.h"
     46#include "pmDetections.h"
     47
    2948#include "pmModel_QGAUSS.h"
    3049
     50# define PM_MODEL_NPARAM          8
    3151# define PM_MODEL_FUNC            pmModelFunc_QGAUSS
    3252# define PM_MODEL_FLUX            pmModelFlux_QGAUSS
    … …  
    123143        return true;
    124144    }
    125     psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
     145    psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds");
    126146
    127147    // we need to calculate the limits for SXY specially
    … …  
    469489    return;
    470490}
    471 
    472 # undef PM_MODEL_FUNC
    473 # undef PM_MODEL_FLUX
    474 # undef PM_MODEL_GUESS
    475 # undef PM_MODEL_LIMITS
    476 # undef PM_MODEL_RADIUS
    477 # undef PM_MODEL_FROM_PSF
    478 # undef PM_MODEL_PARAMS_FROM_PSF
    479 # undef PM_MODEL_FIT_STATUS
    480 # undef PM_MODEL_SET_LIMITS
    481 # undef ALPHA
    482 # undef ALPHA_M

  • trunk/psModules/src/objects/models/pmModel_RGAUSS.c

    r27565 r29004  
    2222#include <stdio.h>
    2323#include <pslib.h>
    24 
     24#include "pmHDU.h"
     25#include "pmFPA.h"
     26
     27#include "pmTrend2D.h"
     28#include "pmResiduals.h"
     29#include "pmGrowthCurve.h"
     30#include "pmSpan.h"
     31#include "pmFootprintSpans.h"
     32#include "pmFootprint.h"
     33#include "pmPeaks.h"
    2534#include "pmMoments.h"
    26 #include "pmPeaks.h"
     35#include "pmModelFuncs.h"
     36#include "pmModel.h"
     37#include "pmModelUtils.h"
     38#include "pmModelClass.h"
     39#include "pmSourceMasks.h"
     40#include "pmSourceExtendedPars.h"
     41#include "pmSourceDiffStats.h"
    2742#include "pmSource.h"
    28 #include "pmModel.h"
     43#include "pmSourceFitModel.h"
     44#include "pmPSF.h"
     45#include "pmPSFtry.h"
     46#include "pmDetections.h"
     47
    2948#include "pmModel_RGAUSS.h"
    3049
     50# define PM_MODEL_NPARAM          8
    3151# define PM_MODEL_FUNC            pmModelFunc_RGAUSS
    3252# define PM_MODEL_FLUX            pmModelFlux_RGAUSS
    … …  
    112132        return true;
    113133    }
    114     psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
     134    psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds");
    115135
    116136    // we need to calculate the limits for SXY specially
    … …  
    209229    PAR[PM_PAR_XPOS] = peak->xf;
    210230    PAR[PM_PAR_YPOS] = peak->yf;
    211     PAR[PM_PAR_SXX]  = PS_MAX(0.5, M_SQRT2*shape.sx);
    212     PAR[PM_PAR_SYY]  = PS_MAX(0.5, M_SQRT2*shape.sy);
     231    PAR[PM_PAR_SXX]  = PS_MAX(0.5, shape.sx);
     232    PAR[PM_PAR_SYY]  = PS_MAX(0.5, shape.sy);
    213233    PAR[PM_PAR_SXY]  = shape.sxy;
    214     PAR[PM_PAR_7]    = 2.25;
     234    PAR[PM_PAR_7]    = 1.5;
    215235
    216236    return(true);
    … …  
    463483    return;
    464484}
    465 
    466 # undef PM_MODEL_FUNC
    467 # undef PM_MODEL_FLUX
    468 # undef PM_MODEL_GUESS
    469 # undef PM_MODEL_LIMITS
    470 # undef PM_MODEL_RADIUS
    471 # undef PM_MODEL_FROM_PSF
    472 # undef PM_MODEL_PARAMS_FROM_PSF
    473 # undef PM_MODEL_FIT_STATUS
    474 # undef PM_MODEL_SET_LIMITS

  • trunk/psModules/src/objects/models/pmModel_SERSIC.c

    r26916 r29004  
    1818   * PM_PAR_7   7   - normalized sersic parameter
    1919
    20    note that a standard sersic model uses exp(-K*(z^(1/n) - 1).  the additional elements (K,
     20   note that a standard sersic model uses exp(-K*(z^(1/2n) - 1).  the additional elements (K,
    2121   the -1 offset) are absorbed in this model by the normalization, the exponent, and the
    2222   radial scale.  We fit the elements in this form, then re-normalize them on output.
    … …  
    2525#include <stdio.h>
    2626#include <pslib.h>
    27 
     27#include "pmHDU.h"
     28#include "pmFPA.h"
     29
     30#include "pmTrend2D.h"
     31#include "pmResiduals.h"
     32#include "pmGrowthCurve.h"
     33#include "pmSpan.h"
     34#include "pmFootprintSpans.h"
     35#include "pmFootprint.h"
     36#include "pmPeaks.h"
    2837#include "pmMoments.h"
    29 #include "pmPeaks.h"
     38#include "pmModelFuncs.h"
     39#include "pmModel.h"
     40#include "pmModelUtils.h"
     41#include "pmModelClass.h"
     42#include "pmSourceMasks.h"
     43#include "pmSourceExtendedPars.h"
     44#include "pmSourceDiffStats.h"
    3045#include "pmSource.h"
    31 #include "pmModel.h"
     46#include "pmSourceFitModel.h"
     47#include "pmPSF.h"
     48#include "pmPSFtry.h"
     49#include "pmDetections.h"
     50
    3251#include "pmModel_SERSIC.h"
    3352
     53# define PM_MODEL_NPARAM          8
    3454# define PM_MODEL_FUNC            pmModelFunc_SERSIC
    3555# define PM_MODEL_FLUX            pmModelFlux_SERSIC
    … …  
    4767
    4868// Lax parameter limits
    49 static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.05, 0.05, -1.0, 0.05 };
     69static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0, 0.05 };
    5070static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 };
    5171
    … …  
    84104    assert (z >= 0);
    85105
    86     psF32 f2 = pow(z,PAR[PM_PAR_7]);
    87     psF32 f1 = exp(-f2);
     106    float index = 0.5 / PAR[PM_PAR_7];
     107    float bn = 1.9992*index - 0.3271;
     108    float Io = exp(bn);
     109
     110    psF32 f2 = bn*pow(z,PAR[PM_PAR_7]);
     111    psF32 f1 = Io*exp(-f2);
    88112    psF32 z0 = PAR[PM_PAR_I0]*f1;
    89113    psF32 f0 = PAR[PM_PAR_SKY] + z0;
    … …  
    98122
    99123        // gradient is infinite for z = 0; saturate at z = 0.01
    100         psF32 z1 = (z < 0.01) ? z0*PAR[PM_PAR_7]*pow(0.01,PAR[PM_PAR_7] - 1.0) : z0*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7] - 1.0);
     124        psF32 z1 = (z < 0.01) ? z0*bn*PAR[PM_PAR_7]*pow(0.01,PAR[PM_PAR_7] - 1.0) : z0*bn*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7] - 1.0);
    101125
    102126        dPAR[PM_PAR_SKY]  = +1.0;
    103127        dPAR[PM_PAR_I0]   = +f1;
    104         dPAR[PM_PAR_7]    = (z == 0.0) ? 0.0 : -z0*f2*log(z);
     128        dPAR[PM_PAR_7]    = (z < 0.01) ? -z0*pow(0.01,PAR[PM_PAR_7])*log(0.01) : -z0*f2*log(z);
     129        dPAR[PM_PAR_7]   *= 3.0;
    105130
    106131        assert (isfinite(z1));
    … …  
    109134        dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]);
    110135        dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]);
    111         dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX];
     136        dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX]; // XXX : increase drag?
    112137        dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
    113         dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
    114138        dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
    115139    }
    … …  
    127151        return true;
    128152    }
    129     psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
     153    psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds");
    130154
    131155    // we need to calculate the limits for SXY specially
    … …  
    201225    psF32     *PAR  = model->params->data.F32;
    202226
     227    // the other parameters depend on the guess for PAR_7
     228    if (!isfinite(PAR[PM_PAR_7])) {
     229        PAR[PM_PAR_7]    = 0.25;
     230    }   
     231    float index = 0.5 / PAR[PM_PAR_7];
     232
     233    // the scale-length is a function of the moments and the index:
     234    // Rmajor_guess = Rmajor_moments * Scale * Zero
     235    float Scale = 0.70 + 0.053 * PAR[PM_PAR_7];
     236    float Zero  = 1.16 - 0.615 * PAR[PM_PAR_7];
     237
    203238    psEllipseMoments emoments;
    204239    emoments.x2 = moments->Mxx;
    … …  
    213248    if (!isfinite(axes.theta)) return false;
    214249
     250    // set a lower limit to avoid absurd solutions..
     251    float Rmajor = PS_MAX(1.0, Scale * axes.major + Zero);
     252    float Rminor = Rmajor * (axes.minor / axes.major);
     253
     254    fprintf (stderr, "guess index: %f : %f, %f -> %f, %f\n", index, axes.major, axes.minor, Rmajor, Rminor);
     255
     256    axes.major = Rmajor;
     257    axes.minor = Rminor;
    215258    psEllipseShape shape = psEllipseAxesToShape (axes);
    216259
    … …  
    219262    if (!isfinite(shape.sxy)) return false;
    220263
     264    float bn = 1.9992*index - 0.3271;
     265    // float fR = 1.0 / (sqrt(2.0) * pow (bn, index));
     266    float Io = exp(0.5*bn);
     267
     268    // XXX do we need this factor of sqrt(2)?
     269    // float Sxx = PS_MAX(0.5, M_SQRT2*shape.sx);
     270    // float Syy = PS_MAX(0.5, M_SQRT2*shape.sy);
     271
     272    float Sxx = PS_MAX(0.5, shape.sx);
     273    float Syy = PS_MAX(0.5, shape.sy);
     274
    221275    PAR[PM_PAR_SKY]  = 0.0;
    222     PAR[PM_PAR_I0]   = peak->flux;
     276    PAR[PM_PAR_I0]   = peak->flux / Io;
    223277    PAR[PM_PAR_XPOS] = peak->xf;
    224278    PAR[PM_PAR_YPOS] = peak->yf;
    225     PAR[PM_PAR_SXX]  = PS_MAX(0.5, M_SQRT2*shape.sx);
    226     PAR[PM_PAR_SYY]  = PS_MAX(0.5, M_SQRT2*shape.sy);
     279    PAR[PM_PAR_SXX]  = Sxx;
     280    PAR[PM_PAR_SYY]  = Syy;
    227281    PAR[PM_PAR_SXY]  = shape.sxy;
    228     PAR[PM_PAR_7]    = 0.5;
    229282
    230283    return(true);
    … …  
    254307    float f1, f2;
    255308    for (z = DZ; z < 50; z += DZ) {
    256         f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
     309        // f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
     310        f1 = exp(-pow(z,PAR[PM_PAR_7]));
    257311        z += DZ;
    258         f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
     312        f2 = exp(-pow(z,PAR[PM_PAR_7]));
    259313        norm += f0 + 4*f1 + f2;
    260314        f0 = f2;
    … …  
    287341
    288342    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
    289     psF64 sigma = axes.major;
    290 
    291     psF64 limit = flux / PAR[PM_PAR_I0];
    292 
    293     psF64 z = pow (-log(limit), (1.0 / PAR[PM_PAR_7]));
    294     psAssert (isfinite(z), "fix this code: z should not be nan for %f", PAR[PM_PAR_7]);
    295 
    296     psF64 radius = sigma * sqrt (2.0 * z);
    297     psAssert (isfinite(radius), "fix this code: radius should not be nan for %f, %f", PAR[PM_PAR_7], sigma);
    298 
    299     if (isnan(radius))
    300         psAbort("error in code: radius is NaN");
    301 
     343
     344    // f = Io exp(-z^n) -> z^n = ln(Io/f)
     345    psF64 zn = log(PAR[PM_PAR_I0] / flux);
     346    psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / PAR[PM_PAR_7]);
     347
     348    fprintf (stderr, "sersic model %f %f, n %f, radius: %f, zn: %f, f/Io: %f, major: %f\n", PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], PAR[PM_PAR_7], radius, zn, flux/PAR[PM_PAR_I0], axes.major);
     349
     350    psAssert (isfinite(radius), "fix this code: z should not be nan for %f", PAR[PM_PAR_7]);
    302351    return (radius);
    303352}
    … …  
    448497    return;
    449498}
    450 
    451 # undef PM_MODEL_FUNC
    452 # undef PM_MODEL_FLUX
    453 # undef PM_MODEL_GUESS
    454 # undef PM_MODEL_LIMITS
    455 # undef PM_MODEL_RADIUS
    456 # undef PM_MODEL_FROM_PSF
    457 # undef PM_MODEL_PARAMS_FROM_PSF
    458 # undef PM_MODEL_FIT_STATUS
    459 # undef PM_MODEL_SET_LIMITS

  • trunk/psModules/src/objects/pmDetEff.c

    r25477 r29004  
    33#endif
    44
     5#include <string.h>
    56#include <pslib.h>
     7#include "pmHDU.h"
     8#include "pmFPA.h"
    69
     10#include "pmTrend2D.h"
     11#include "pmResiduals.h"
     12#include "pmGrowthCurve.h"
     13#include "pmSpan.h"
     14#include "pmFootprintSpans.h"
     15#include "pmFootprint.h"
     16#include "pmPeaks.h"
     17#include "pmMoments.h"
     18#include "pmModelFuncs.h"
     19#include "pmModel.h"
     20#include "pmModelUtils.h"
     21#include "pmModelClass.h"
     22#include "pmSourceMasks.h"
     23#include "pmSourceExtendedPars.h"
     24#include "pmSourceDiffStats.h"
     25#include "pmSource.h"
     26#include "pmSourceFitModel.h"
     27#include "pmPSF.h"
     28#include "pmPSFtry.h"
     29#include "pmDetections.h"
    730#include "pmDetEff.h"
    8 
    931
    1032static void detEffFree(pmDetEff *de)

  • trunk/psModules/src/objects/pmDetEff.h

    r25383 r29004  
    11#ifndef PM_DET_EFF_H
    22#define PM_DET_EFF_H
    3 
    4 #include <pslib.h>
    5 #include <string.h>
    6 
    7 #include "pmFPA.h"
    83
    94#define PM_DETEFF_ANALYSIS "DETEFF"     // Location of detection efficiency on pmReadout.analysis

  • trunk/psModules/src/objects/pmFootprint.c

    r27672 r29004  
    1919#include <pslib.h>
    2020#include "pmSpan.h"
     21#include "pmFootprintSpans.h"
    2122#include "pmFootprint.h"
    2223#include "pmPeaks.h"

  • trunk/psModules/src/objects/pmFootprint.h

    r27672 r29004  
    1010#ifndef PM_FOOTPRINT_H
    1111#define PM_FOOTPRINT_H
    12 
    13 #include <pslib.h>
    14 #include "pmSpan.h"
    15 #include "pmFootprintSpans.h"
    1612
    1713typedef struct {

  • trunk/psModules/src/objects/pmFootprintArrayGrow.c

    r21183 r29004  
    1818#include <pslib.h>
    1919#include "pmSpan.h"
     20#include "pmFootprintSpans.h"
    2021#include "pmFootprint.h"
    2122#include "pmPeaks.h"

  • trunk/psModules/src/objects/pmFootprintArraysMerge.c

    r24274 r29004  
    1818#include <pslib.h>
    1919#include "pmSpan.h"
     20#include "pmFootprintSpans.h"
    2021#include "pmFootprint.h"
    2122#include "pmPeaks.h"

  • trunk/psModules/src/objects/pmFootprintAssignPeaks.c

    r20937 r29004  
    1818#include <pslib.h>
    1919#include "pmSpan.h"
     20#include "pmFootprintSpans.h"
    2021#include "pmFootprint.h"
    2122#include "pmPeaks.h"

  • trunk/psModules/src/objects/pmFootprintCullPeaks.c

    r27672 r29004  
    1919#include <pslib.h>
    2020#include "pmSpan.h"
     21#include "pmFootprintSpans.h"
    2122#include "pmFootprint.h"
    22 #include "pmFootprintSpans.h"
    2323#include "pmPeaks.h"
    2424

  • trunk/psModules/src/objects/pmFootprintFind.c

    r20937 r29004  
    1919#include <pslib.h>
    2020#include "pmSpan.h"
     21#include "pmFootprintSpans.h"
    2122#include "pmFootprint.h"
    2223#include "pmPeaks.h"

  • trunk/psModules/src/objects/pmFootprintFindAtPoint.c

    r27672 r29004  
    1919#include <pslib.h>
    2020#include "pmSpan.h"
     21#include "pmFootprintSpans.h"
    2122#include "pmFootprint.h"
    2223#include "pmPeaks.h"
    23 #include "pmFootprintSpans.h"
    2424
    2525/************************************************************************************************************/

  • trunk/psModules/src/objects/pmFootprintIDs.c

    r27672 r29004  
    1919#include <pslib.h>
    2020#include "pmSpan.h"
     21#include "pmFootprintSpans.h"
    2122#include "pmFootprint.h"
    2223#include "pmPeaks.h"

  • trunk/psModules/src/objects/pmFootprintSpans.c

    r27672 r29004  
    1919#include <pslib.h>
    2020#include "pmSpan.h"
     21#include "pmFootprintSpans.h"
    2122#include "pmFootprint.h"
    2223#include "pmPeaks.h"
    23 #include "pmFootprintSpans.h"
    2424
    2525static void pmStartSpanFree(pmStartSpan *sspan) {

  • trunk/psModules/src/objects/pmFootprintSpans.h

    r27672 r29004  
    1010#ifndef PM_FOOTPRINT_SPANS_H
    1111#define PM_FOOTPRINT_SPANS_H
    12 
    13 #include <pslib.h>
    14 #include "pmSpan.h"
    1512
    1613/* We define two helper structures used in building the pmFootprints:

  • trunk/psModules/src/objects/pmGrowthCurve.c

    r20937 r29004  
    2020#include "pmHDU.h"
    2121#include "pmFPA.h"
     22
     23#include "pmTrend2D.h"
     24#include "pmResiduals.h"
     25#include "pmGrowthCurve.h"
    2226#include "pmSpan.h"
     27#include "pmFootprintSpans.h"
    2328#include "pmFootprint.h"
    2429#include "pmPeaks.h"
    2530#include "pmMoments.h"
    26 #include "pmResiduals.h"
    27 #include "pmGrowthCurve.h"
    28 #include "pmTrend2D.h"
     31#include "pmModelFuncs.h"
     32#include "pmModel.h"
     33#include "pmModelUtils.h"
     34#include "pmModelClass.h"
     35#include "pmSourceMasks.h"
     36#include "pmSourceExtendedPars.h"
     37#include "pmSourceDiffStats.h"
     38#include "pmSource.h"
     39#include "pmSourceFitModel.h"
    2940#include "pmPSF.h"
    30 #include "pmModel.h"
    31 #include "pmSource.h"
     41#include "pmPSFtry.h"
     42#include "pmDetections.h"
     43
    3244#include "psVectorBracket.h"
    3345

  • trunk/psModules/src/objects/pmGrowthCurveGenerate.c

    r27531 r29004  
    2424#include "pmHDU.h"
    2525#include "pmFPA.h"
     26#include "pmFPAMaskWeight.h"
     27#include "psVectorBracket.h"
     28#include "pmErrorCodes.h"
     29
     30#include "pmTrend2D.h"
     31#include "pmResiduals.h"
     32#include "pmGrowthCurve.h"
    2633#include "pmSpan.h"
     34#include "pmFootprintSpans.h"
    2735#include "pmFootprint.h"
    2836#include "pmPeaks.h"
    2937#include "pmMoments.h"
    30 #include "pmResiduals.h"
    31 #include "pmGrowthCurve.h"
    32 #include "pmTrend2D.h"
     38#include "pmModelFuncs.h"
     39#include "pmModel.h"
     40#include "pmModelUtils.h"
     41#include "pmModelClass.h"
     42#include "pmSourceMasks.h"
     43#include "pmSourceExtendedPars.h"
     44#include "pmSourceDiffStats.h"
     45#include "pmSource.h"
     46#include "pmSourceFitModel.h"
    3347#include "pmPSF.h"
    34 #include "pmModel.h"
    35 #include "pmSource.h"
    36 #include "pmModelClass.h"
    37 #include "pmModelUtils.h"
     48#include "pmPSFtry.h"
     49#include "pmDetections.h"
     50
    3851#include "pmSourcePhotometry.h"
    39 #include "pmFPAMaskWeight.h"
    40 #include "psVectorBracket.h"
    41 #include "pmErrorCodes.h"
    4252
    4353pmGrowthCurve *pmGrowthCurveForPosition (psImage *image, pmPSF *psf, bool ignore, psImageMaskType maskVal, psImageMaskType markVal, float xc, float yc);

  • trunk/psModules/src/objects/pmModel.c

    r25754 r29004  
    2323#include "pmHDU.h"
    2424#include "pmFPA.h"
     25
     26#include "pmTrend2D.h"
    2527#include "pmResiduals.h"
    2628#include "pmGrowthCurve.h"
    27 #include "pmTrend2D.h"
    28 #include "pmPSF.h"
     29#include "pmSpan.h"
     30#include "pmFootprintSpans.h"
     31#include "pmFootprint.h"
     32#include "pmPeaks.h"
     33#include "pmMoments.h"
     34#include "pmModelFuncs.h"
    2935#include "pmModel.h"
    3036#include "pmModelClass.h"
    … …  
    7581
    7682    for (psS32 i = 0; i < tmp->params->n; i++) {
    77         tmp->params->data.F32[i] = 0.0;
    78         tmp->dparams->data.F32[i] = 0.0;
     83        tmp->params->data.F32[i] = NAN;
     84        tmp->dparams->data.F32[i] = NAN;
    7985    }
    8086

  • trunk/psModules/src/objects/pmModel.h

    r26916 r29004  
    1414# define PM_MODEL_H
    1515
    16 #include <pslib.h>
    17 #include "pmPSF.h"
    18 
    1916/// @addtogroup Objects Object Detection / Analysis Functions
    2017/// @{
    … …  
    2219/* pointers for the functions types below are supplied to each pmModel, and can be used by
    2320   the programmer without needing to know the model class */
    24 
    25 typedef enum {
    26     PM_MODEL_STATUS_NONE         = 0x00, ///< model fit not yet attempted, no other info
    27     PM_MODEL_STATUS_FITTED       = 0x01, ///< model fit completed
    28     PM_MODEL_STATUS_NONCONVERGE  = 0x02, ///< model fit did not converge
    29     PM_MODEL_STATUS_OFFIMAGE     = 0x04, ///< model fit drove out of range
    30     PM_MODEL_STATUS_BADARGS      = 0x08, ///< model fit called with invalid args
    31     PM_MODEL_STATUS_LIMITS       = 0x10  ///< model parameters hit limits
    32 } pmModelStatus;
    33 
    34 typedef enum {
    35     PM_MODEL_OP_NONE    = 0x00,
    36     PM_MODEL_OP_FUNC    = 0x01,
    37     PM_MODEL_OP_RES0    = 0x02,
    38     PM_MODEL_OP_RES1    = 0x04,
    39     PM_MODEL_OP_FULL    = 0x07,
    40     PM_MODEL_OP_SKY     = 0x08,
    41     PM_MODEL_OP_CENTER  = 0x10,
    42     PM_MODEL_OP_NORM    = 0x20,
    43     PM_MODEL_OP_NOISE   = 0x40,
    44 } pmModelOpMode;
    45 
    46 /// Parameter limit types
    47 typedef enum {
    48     PM_MODEL_LIMITS_NONE,               ///< Apply no limits: suitable for debugging
    49     PM_MODEL_LIMITS_IGNORE,             ///< Ignore all limits: fit can go to town
    50     PM_MODEL_LIMITS_LAX,                ///< Lax limits: attempting to reproduce even bad data
    51     PM_MODEL_LIMITS_MODERATE,           ///< Moderate limits: cope with mildly bad data
    52     PM_MODEL_LIMITS_STRICT,             ///< Strict limits: insist on good quality data
    53 } pmModelLimitsType;
    54 
    55 typedef struct pmModel pmModel;
    56 typedef struct pmSource pmSource;
    57 
    58 //  This function is the model chi-square minimization function for this model.
    59 typedef psMinimizeLMChi2Func pmModelFunc;
    60 
    61 //  This function sets the parameter limits for this model.
    62 typedef psMinimizeLMLimitFunc pmModelLimits;
    63 
    64 // This function returns the integrated flux for the given model parameters.
    65 typedef psF64 (*pmModelFlux)(const psVector *params);
    66 
    67 // This function returns the radius at which the given model and parameters
    68 // achieves the given flux.
    69 typedef psF64 (*pmModelRadius)(const psVector *params, double flux);
    70 
    71 //  This function provides the model guess parameters based on the details of
    72 //  the given source.
    73 typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source);
    74 
    75 //  This function constructs the PSF model for the given source based on the
    76 //  supplied psf and the EXT model for the object.
    77 typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelEXT, const pmPSF *psf);
    78 
    79 //  This function sets the model parameters based on the PSF for a given coordinate and central
    80 //  intensity
    81 typedef bool (*pmModelParamsFromPSF)(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io);
    82 
    83 //  This function returns the success / failure status of the given model fit
    84 typedef bool (*pmModelFitStatusFunc)(pmModel *model);
    85 
    86 //  This function sets the parameter limits for the given model
    87 typedef bool (*pmModelSetLimitsFunc)(pmModelLimitsType type);
    8821
    8922/** pmModel data structure
    … …  
    12255    pmModelSetLimitsFunc modelSetLimits;
    12356};
    124 
    125 /** Symbolic names for the elements of [d]params
    126  * Note: these are #defines not enums as a given element of [d]params
    127  * may/will correspond to different parameters in different contexts
    128  */
    129 #define PM_PAR_SKY 0   ///< Sky
    130 #define PM_PAR_I0 1   ///< Central intensity
    131 #define PM_PAR_XPOS 2   ///< X centre of object
    132 #define PM_PAR_YPOS 3   ///< Y centre of object
    133 #define PM_PAR_SXX 4   ///< shape X^2 moment
    134 #define PM_PAR_SYY 5   ///< shape Y^2 moment
    135 #define PM_PAR_SXY 6   ///< shape XY moment
    136 #define PM_PAR_7 7   ///< ??? Unknown parameter
    137 #define PM_PAR_8 8   ///< ??? Unknown parameter
    13857
    13958/** pmModelAlloc()

  • trunk/psModules/src/objects/pmModelClass.c

    r25754 r29004  
    2323#include "pmHDU.h"
    2424#include "pmFPA.h"
     25
     26#include "pmTrend2D.h"
     27#include "pmResiduals.h"
     28#include "pmGrowthCurve.h"
    2529#include "pmSpan.h"
     30#include "pmFootprintSpans.h"
    2631#include "pmFootprint.h"
    2732#include "pmPeaks.h"
    2833#include "pmMoments.h"
    29 #include "pmResiduals.h"
    30 #include "pmGrowthCurve.h"
    31 #include "pmTrend2D.h"
    32 #include "pmPSF.h"
     34#include "pmModelFuncs.h"
    3335#include "pmModel.h"
    34 #include "pmSource.h"
     36#include "pmModelUtils.h"
    3537#include "pmModelClass.h"
     38
    3639#include "pmErrorCodes.h"
    3740
    … …  
    4649# include "models/pmModel_RGAUSS.h"
    4750# include "models/pmModel_SERSIC.h"
     51# include "models/pmModel_EXP.h"
     52# include "models/pmModel_DEV.h"
    4853
    4954static pmModelClass defaultModels[] = {
    … …  
    5358    {"PS_MODEL_PS1_V1",       8, (pmModelFunc)pmModelFunc_PS1_V1,  (pmModelFlux)pmModelFlux_PS1_V1,  (pmModelRadius)pmModelRadius_PS1_V1,  (pmModelLimits)pmModelLimits_PS1_V1,  (pmModelGuessFunc)pmModelGuess_PS1_V1, (pmModelFromPSFFunc)pmModelFromPSF_PS1_V1, (pmModelParamsFromPSF)pmModelParamsFromPSF_PS1_V1, (pmModelFitStatusFunc)pmModelFitStatus_PS1_V1, (pmModelSetLimitsFunc)pmModelSetLimits_PS1_V1 },
    5459    {"PS_MODEL_RGAUSS",       8, (pmModelFunc)pmModelFunc_RGAUSS,  (pmModelFlux)pmModelFlux_RGAUSS,  (pmModelRadius)pmModelRadius_RGAUSS,  (pmModelLimits)pmModelLimits_RGAUSS,  (pmModelGuessFunc)pmModelGuess_RGAUSS, (pmModelFromPSFFunc)pmModelFromPSF_RGAUSS, (pmModelParamsFromPSF)pmModelParamsFromPSF_RGAUSS, (pmModelFitStatusFunc)pmModelFitStatus_RGAUSS, (pmModelSetLimitsFunc)pmModelSetLimits_RGAUSS },
    55     {"PS_MODEL_SERSIC",       8, (pmModelFunc)pmModelFunc_SERSIC,  (pmModelFlux)pmModelFlux_SERSIC,  (pmModelRadius)pmModelRadius_SERSIC,  (pmModelLimits)pmModelLimits_SERSIC,  (pmModelGuessFunc)pmModelGuess_SERSIC, (pmModelFromPSFFunc)pmModelFromPSF_SERSIC, (pmModelParamsFromPSF)pmModelParamsFromPSF_SERSIC, (pmModelFitStatusFunc)pmModelFitStatus_SERSIC, (pmModelSetLimitsFunc)pmModelSetLimits_SERSIC }
     60    {"PS_MODEL_SERSIC",       8, (pmModelFunc)pmModelFunc_SERSIC,  (pmModelFlux)pmModelFlux_SERSIC,  (pmModelRadius)pmModelRadius_SERSIC,  (pmModelLimits)pmModelLimits_SERSIC,  (pmModelGuessFunc)pmModelGuess_SERSIC, (pmModelFromPSFFunc)pmModelFromPSF_SERSIC, (pmModelParamsFromPSF)pmModelParamsFromPSF_SERSIC, (pmModelFitStatusFunc)pmModelFitStatus_SERSIC, (pmModelSetLimitsFunc)pmModelSetLimits_SERSIC },
     61    {"PS_MODEL_EXP",          7, (pmModelFunc)pmModelFunc_EXP,     (pmModelFlux)pmModelFlux_EXP,     (pmModelRadius)pmModelRadius_EXP,     (pmModelLimits)pmModelLimits_EXP,     (pmModelGuessFunc)pmModelGuess_EXP,    (pmModelFromPSFFunc)pmModelFromPSF_EXP,    (pmModelParamsFromPSF)pmModelParamsFromPSF_EXP,    (pmModelFitStatusFunc)pmModelFitStatus_EXP,    (pmModelSetLimitsFunc)pmModelSetLimits_EXP },
     62    {"PS_MODEL_DEV",          7, (pmModelFunc)pmModelFunc_DEV,     (pmModelFlux)pmModelFlux_DEV,     (pmModelRadius)pmModelRadius_DEV,     (pmModelLimits)pmModelLimits_DEV,     (pmModelGuessFunc)pmModelGuess_DEV,    (pmModelFromPSFFunc)pmModelFromPSF_DEV,    (pmModelParamsFromPSF)pmModelParamsFromPSF_DEV,    (pmModelFitStatusFunc)pmModelFitStatus_DEV,    (pmModelSetLimitsFunc)pmModelSetLimits_DEV },
    5663};
    5764

  • trunk/psModules/src/objects/pmModelClass.h

    r25738 r29004  
    2828# ifndef PM_MODEL_CLASS_H
    2929# define PM_MODEL_CLASS_H
    30 
    31 #include <pmModel.h>
    3230
    3331/// @addtogroup Objects Object Detection / Analysis Functions

  • trunk/psModules/src/objects/pmModelUtils.c

    r25754 r29004  
    2222#include "pmHDU.h"
    2323#include "pmFPA.h"
     24
     25#include "pmTrend2D.h"
    2426#include "pmResiduals.h"
    2527#include "pmGrowthCurve.h"
    26 #include "pmTrend2D.h"
     28#include "pmSpan.h"
     29#include "pmFootprintSpans.h"
     30#include "pmFootprint.h"
     31#include "pmPeaks.h"
     32#include "pmMoments.h"
     33#include "pmModelFuncs.h"
     34#include "pmModel.h"
     35#include "pmModelUtils.h"
     36#include "pmModelClass.h"
     37#include "pmSourceMasks.h"
     38#include "pmSourceExtendedPars.h"
     39#include "pmSourceDiffStats.h"
     40#include "pmSource.h"
     41#include "pmSourceFitModel.h"
    2742#include "pmPSF.h"
    28 #include "pmModel.h"
     43
    2944#include "pmErrorCodes.h"
    30 #include "pmModelUtils.h"
    3145
    3246/*****************************************************************************

  • trunk/psModules/src/objects/pmMoments.c

    r24498 r29004  
    2929    psTrace("psModules.objects", 10, "---- %s() begin ----\n", __func__);
    3030    pmMoments *tmp = (pmMoments *) psAlloc(sizeof(pmMoments));
     31
     32    tmp->Mrf = 0.0;
     33    tmp->Mrh = 0.0;
     34    tmp->KronFlux = 0.0;
     35    tmp->KronFluxErr = 0.0;
     36
     37    tmp->KronFinner = 0.0;
     38    tmp->KronFouter = 0.0;
     39
    3140    tmp->Mx = 0.0;
    3241    tmp->My = 0.0;

  • trunk/psModules/src/objects/pmMoments.h

    r23487 r29004  
    2424typedef struct
    2525{
     26    float Mrf;    ///< radial first moment
     27    float Mrh;    ///< radial half moment
     28
    2629    float Mx;     ///< X-coord of centroid.
    2730    float My;     ///< Y-coord of centroid.
    … …  
    4750    float SN;     ///< approx signal-to-noise
    4851    int nPixels;  ///< Number of pixels used.
     52
     53    float KronFlux;    ///< Kron flux (flux in 2.5*Mrf)
     54    float KronFluxErr; ///< Kron flux error
     55
     56    float KronFinner;    ///< Kron flux (flux in 2.5*Mrf)
     57    float KronFouter;    ///< Kron flux (flux in 2.5*Mrf)
     58
    4959}
    5060pmMoments;

  • trunk/psModules/src/objects/pmPSF.c

    r26893 r29004  
    2525#include "pmHDU.h"
    2626#include "pmFPA.h"
     27#include "pmFPAMaskWeight.h"
     28#include "psVectorBracket.h"
     29
     30#include "pmTrend2D.h"
     31#include "pmResiduals.h"
     32#include "pmGrowthCurve.h"
    2733#include "pmSpan.h"
     34#include "pmFootprintSpans.h"
    2835#include "pmFootprint.h"
    2936#include "pmPeaks.h"
    3037#include "pmMoments.h"
    31 #include "pmResiduals.h"
    32 #include "pmGrowthCurve.h"
    33 #include "pmTrend2D.h"
     38#include "pmModelFuncs.h"
     39#include "pmModel.h"
     40#include "pmModelUtils.h"
     41#include "pmModelClass.h"
     42#include "pmSourceMasks.h"
     43#include "pmSourceExtendedPars.h"
     44#include "pmSourceDiffStats.h"
     45#include "pmSource.h"
     46#include "pmSourceFitModel.h"
    3447#include "pmPSF.h"
    35 #include "pmModel.h"
    36 #include "pmSource.h"
    37 #include "pmModelClass.h"
    38 #include "pmModelUtils.h"
    39 #include "pmSourcePhotometry.h"
    40 #include "pmFPAMaskWeight.h"
    41 #include "psVectorBracket.h"
     48#include "pmPSFtry.h"
     49#include "pmDetections.h"
     50
    4251#include "pmErrorCodes.h"
    4352
    … …  
    5463
    5564    psFree (options->stats);
     65    psFree (options->fitOptions);
    5666    return;
    5767}
    … …  
    6575
    6676    options->stats         = NULL;
     77    options->fitOptions    = NULL; // XXX this has to be set before calling pmPSF fit functions
    6778
    6879    options->psfTrendMode  = PM_TREND_NONE;

  • trunk/psModules/src/objects/pmPSF.h

    r26893 r29004  
    1414# define PM_PSF_H
    1515
    16 #include <pslib.h>
    17 #include "pmTrend2D.h"
    18 #include "pmGrowthCurve.h"
    19 #include "pmResiduals.h"
    20 #include "pmFPA.h"
    21 
    2216/// @addtogroup Objects Object Detection / Analysis Functions
    2317/// @{
    24 
    25 // type of model carried by the pmModel structure
    26 typedef int pmModelType;
    2718
    2819/** pmPSF data structure
    … …  
    3829 *
    3930 */
    40 typedef struct {
     31struct pmPSF {
    4132    pmModelType type;                   ///< PSF Model in use
    4233    psArray *params;                    ///< Model parameters (psPolynomial2D)
    … …  
    6455    pmGrowthCurve *growth;              ///< apMag vs Radius
    6556    pmResiduals *residuals;             ///< normalized residual image (no spatial variation)
    66 } pmPSF;
     57};
    6758
    6859typedef struct {
    … …  
    8273    float         apRadius;
    8374    bool          chiFluxTrend;         // Fit a trend in Chi2 as a function of flux?
     75    pmSourceFitOptions *fitOptions;
    8476} pmPSFOptions;
    8577

  • trunk/psModules/src/objects/pmPSF_IO.c

    r28449 r29004  
    3838#include "pmFPAfileFitsIO.h"
    3939
     40#include "pmTrend2D.h"
     41#include "pmResiduals.h"
     42#include "pmGrowthCurve.h"
    4043#include "pmSpan.h"
     44#include "pmFootprintSpans.h"
    4145#include "pmFootprint.h"
    4246#include "pmPeaks.h"
    4347#include "pmMoments.h"
    44 #include "pmGrowthCurve.h"
    45 #include "pmResiduals.h"
    46 #include "pmTrend2D.h"
     48#include "pmModelFuncs.h"
     49#include "pmModel.h"
     50#include "pmModelUtils.h"
     51#include "pmModelClass.h"
     52#include "pmSourceMasks.h"
     53#include "pmSourceExtendedPars.h"
     54#include "pmSourceDiffStats.h"
     55#include "pmSource.h"
     56#include "pmSourceFitModel.h"
    4757#include "pmPSF.h"
    48 #include "pmModel.h"
     58#include "pmPSFtry.h"
     59#include "pmDetections.h"
     60
    4961#include "pmPSF_IO.h"
    50 #include "pmSource.h"
    51 #include "pmModelClass.h"
    52 #include "pmModelUtils.h"
    5362#include "pmSourceIO.h"
    5463

  • trunk/psModules/src/objects/pmPSFtry.c

    r25754 r29004  
    1919#include "pmFPA.h"
    2020#include "pmFPAMaskWeight.h"
     21
     22#include "pmTrend2D.h"
     23#include "pmResiduals.h"
     24#include "pmGrowthCurve.h"
    2125#include "pmSpan.h"
     26#include "pmFootprintSpans.h"
    2227#include "pmFootprint.h"
    2328#include "pmPeaks.h"
    2429#include "pmMoments.h"
    25 #include "pmResiduals.h"
    26 #include "pmGrowthCurve.h"
    27 #include "pmTrend2D.h"
     30#include "pmModelFuncs.h"
     31#include "pmModel.h"
     32#include "pmModelUtils.h"
     33#include "pmModelClass.h"
     34#include "pmSourceMasks.h"
     35#include "pmSourceExtendedPars.h"
     36#include "pmSourceDiffStats.h"
     37#include "pmSource.h"
     38#include "pmSourceFitModel.h"
    2839#include "pmPSF.h"
    29 #include "pmModel.h"
    30 #include "pmSource.h"
    31 #include "pmSourceUtils.h"
    3240#include "pmPSFtry.h"
    33 #include "pmModelClass.h"
    34 #include "pmModelUtils.h"
    35 #include "pmSourceFitModel.h"
     41#include "pmDetections.h"
     42
    3643#include "pmSourcePhotometry.h"
    3744#include "pmSourceVisual.h"

  • trunk/psModules/src/objects/pmPSFtryFitEXT.c

    r25754 r29004  
    1919#include "pmFPA.h"
    2020#include "pmFPAMaskWeight.h"
     21
     22#include "pmTrend2D.h"
     23#include "pmResiduals.h"
     24#include "pmGrowthCurve.h"
    2125#include "pmSpan.h"
     26#include "pmFootprintSpans.h"
    2227#include "pmFootprint.h"
    2328#include "pmPeaks.h"
    2429#include "pmMoments.h"
    25 #include "pmResiduals.h"
    26 #include "pmGrowthCurve.h"
    27 #include "pmTrend2D.h"
    28 #include "pmPSF.h"
     30#include "pmModelFuncs.h"
    2931#include "pmModel.h"
     32#include "pmModelUtils.h"
     33#include "pmModelClass.h"
     34#include "pmSourceMasks.h"
     35#include "pmSourceExtendedPars.h"
     36#include "pmSourceDiffStats.h"
    3037#include "pmSource.h"
    3138#include "pmSourceUtils.h"
     39#include "pmSourceFitModel.h"
     40#include "pmPSF.h"
    3241#include "pmPSFtry.h"
    33 #include "pmModelClass.h"
    34 #include "pmModelUtils.h"
    35 #include "pmSourceFitModel.h"
     42#include "pmDetections.h"
     43
    3644#include "pmSourcePhotometry.h"
    3745#include "pmSourceVisual.h"
    … …  
    4452
    4553    psTimerStart ("psf.fit");
     54
     55    // in this segment, we are fitting the full PSF model class (shape unconstrained)
     56    options->fitOptions->mode = PM_SOURCE_FIT_EXT;
    4657
    4758    // maskVal is used to test for rejected pixels, and must include markVal
    … …  
    7384
    7485        // fit model as EXT, not PSF
    75         status = pmSourceFitModel (source, source->modelEXT, PM_SOURCE_FIT_EXT, maskVal);
     86        status = pmSourceFitModel (source, source->modelEXT, options->fitOptions, maskVal);
    7687
    7788        // clear object mask to define valid pixels

  • trunk/psModules/src/objects/pmPSFtryFitPSF.c

    r26027 r29004  
    1717#include "pmFPA.h"
    1818#include "pmFPAMaskWeight.h"
     19
     20#include "pmTrend2D.h"
     21#include "pmResiduals.h"
     22#include "pmGrowthCurve.h"
    1923#include "pmSpan.h"
     24#include "pmFootprintSpans.h"
    2025#include "pmFootprint.h"
    2126#include "pmPeaks.h"
    2227#include "pmMoments.h"
    23 #include "pmResiduals.h"
    24 #include "pmGrowthCurve.h"
    25 #include "pmTrend2D.h"
     28#include "pmModelFuncs.h"
     29#include "pmModel.h"
     30#include "pmModelUtils.h"
     31#include "pmModelClass.h"
     32#include "pmSourceMasks.h"
     33#include "pmSourceExtendedPars.h"
     34#include "pmSourceDiffStats.h"
     35#include "pmSource.h"
     36#include "pmSourceFitModel.h"
    2637#include "pmPSF.h"
    27 #include "pmModel.h"
    28 #include "pmSource.h"
    29 #include "pmSourceUtils.h"
    3038#include "pmPSFtry.h"
    31 #include "pmModelClass.h"
    32 #include "pmModelUtils.h"
    33 #include "pmSourceFitModel.h"
     39#include "pmDetections.h"
     40
    3441#include "pmSourcePhotometry.h"
    3542#include "pmSourceVisual.h"
    … …  
    4249
    4350    psTimerStart ("psf.fit");
     51
     52    // in this segment, we are fitting the fitted PSF model class (shape constrained)
     53    options->fitOptions->mode = PM_SOURCE_FIT_PSF;
    4454
    4555    // maskVal is used to test for rejected pixels, and must include markVal
    … …  
    8191
    8292        // fit the PSF model to the source
    83         status = pmSourceFitModel (source, source->modelPSF, PM_SOURCE_FIT_NORM, maskVal);
     93        status = pmSourceFitModel (source, source->modelPSF, options->fitOptions, maskVal);
    8494
    8595        // skip poor fits
    … …  
    98108
    99109        // This function calculates the psf and aperture magnitudes
    100         status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal); // raw PSF mag, AP mag
     110        status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal, markVal); // raw PSF mag, AP mag
    101111        if (!status || isnan(source->apMag)) {
    102112            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask

  • trunk/psModules/src/objects/pmPSFtryMakePSF.c

    r25754 r29004  
    1818#include "pmFPA.h"
    1919#include "pmFPAMaskWeight.h"
     20
     21#include "pmTrend2D.h"
     22#include "pmResiduals.h"
     23#include "pmGrowthCurve.h"
    2024#include "pmSpan.h"
     25#include "pmFootprintSpans.h"
    2126#include "pmFootprint.h"
    2227#include "pmPeaks.h"
    2328#include "pmMoments.h"
    24 #include "pmResiduals.h"
    25 #include "pmGrowthCurve.h"
    26 #include "pmTrend2D.h"
     29#include "pmModelFuncs.h"
     30#include "pmModel.h"
     31#include "pmModelUtils.h"
     32#include "pmModelClass.h"
     33#include "pmSourceMasks.h"
     34#include "pmSourceExtendedPars.h"
     35#include "pmSourceDiffStats.h"
     36#include "pmSource.h"
     37#include "pmSourceFitModel.h"
    2738#include "pmPSF.h"
    28 #include "pmModel.h"
    29 #include "pmSource.h"
    30 #include "pmSourceUtils.h"
    3139#include "pmPSFtry.h"
    32 #include "pmModelClass.h"
    33 #include "pmModelUtils.h"
    34 #include "pmSourceFitModel.h"
     40#include "pmDetections.h"
     41
    3542#include "pmSourcePhotometry.h"
    3643#include "pmSourceVisual.h"

  • trunk/psModules/src/objects/pmPSFtryMetric.c

    r26260 r29004  
    1818#include "pmFPA.h"
    1919#include "pmFPAMaskWeight.h"
     20
     21#include "pmTrend2D.h"
     22#include "pmResiduals.h"
     23#include "pmGrowthCurve.h"
    2024#include "pmSpan.h"
     25#include "pmFootprintSpans.h"
    2126#include "pmFootprint.h"
    2227#include "pmPeaks.h"
    2328#include "pmMoments.h"
    24 #include "pmResiduals.h"
    25 #include "pmGrowthCurve.h"
    26 #include "pmTrend2D.h"
     29#include "pmModelFuncs.h"
     30#include "pmModel.h"
     31#include "pmModelUtils.h"
     32#include "pmModelClass.h"
     33#include "pmSourceMasks.h"
     34#include "pmSourceExtendedPars.h"
     35#include "pmSourceDiffStats.h"
     36#include "pmSource.h"
     37#include "pmSourceFitModel.h"
    2738#include "pmPSF.h"
    28 #include "pmModel.h"
    29 #include "pmSource.h"
    30 #include "pmSourceUtils.h"
    3139#include "pmPSFtry.h"
    32 #include "pmModelClass.h"
    33 #include "pmModelUtils.h"
    34 #include "pmSourceFitModel.h"
     40#include "pmDetections.h"
     41
    3542#include "pmSourcePhotometry.h"
    3643#include "pmSourceVisual.h"

  • trunk/psModules/src/objects/pmPSFtryModel.c

    r27818 r29004  
    1919#include "pmFPA.h"
    2020#include "pmFPAMaskWeight.h"
     21
     22#include "pmTrend2D.h"
     23#include "pmResiduals.h"
     24#include "pmGrowthCurve.h"
    2125#include "pmSpan.h"
     26#include "pmFootprintSpans.h"
    2227#include "pmFootprint.h"
    2328#include "pmPeaks.h"
    2429#include "pmMoments.h"
    25 #include "pmResiduals.h"
    26 #include "pmGrowthCurve.h"
    27 #include "pmTrend2D.h"
     30#include "pmModelFuncs.h"
     31#include "pmModel.h"
     32#include "pmModelUtils.h"
     33#include "pmModelClass.h"
     34#include "pmSourceMasks.h"
     35#include "pmSourceExtendedPars.h"
     36#include "pmSourceDiffStats.h"
     37#include "pmSource.h"
     38#include "pmSourceFitModel.h"
    2839#include "pmPSF.h"
    29 #include "pmModel.h"
    30 #include "pmSource.h"
    31 #include "pmSourceUtils.h"
    3240#include "pmPSFtry.h"
    33 #include "pmModelClass.h"
    34 #include "pmModelUtils.h"
    35 #include "pmSourceFitModel.h"
    36 #include "pmSourcePhotometry.h"
     41#include "pmDetections.h"
     42
    3743#include "pmSourceVisual.h"
    3844
    … …  
    4955pmPSFtry *pmPSFtryModel (const psArray *sources, const char *modelName, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal)
    5056{
     57    assert (options->fitOptions);
     58
    5159    // validate the requested model name
    5260    options->type = pmModelClassGetType (modelName);

  • trunk/psModules/src/objects/pmPeaks.c

    r26893 r29004  
    2222#include <pslib.h>
    2323#include "pmSpan.h"
     24#include "pmFootprintSpans.h"
    2425#include "pmFootprint.h"
    2526#include "pmPeaks.h"

  • trunk/psModules/src/objects/pmPeaks.h

    r27657 r29004  
    1717# ifndef PM_PEAKS_H
    1818# define PM_PEAKS_H
    19 
    20 #include <pslib.h>
    21 #include "pmFootprint.h"
    2219
    2320/// @addtogroup Objects Object Detection / Analysis Functions

  • trunk/psModules/src/objects/pmPhotObj.c

    r28013 r29004  
    1616#include <string.h>
    1717#include <pslib.h>
     18
     19#include "pmHDU.h"
     20#include "pmFPA.h"
     21
     22#include "pmTrend2D.h"
     23#include "pmResiduals.h"
     24#include "pmGrowthCurve.h"
     25#include "pmSpan.h"
     26#include "pmFootprintSpans.h"
     27#include "pmFootprint.h"
     28#include "pmPeaks.h"
     29#include "pmMoments.h"
     30#include "pmModelFuncs.h"
     31#include "pmModel.h"
     32#include "pmModelUtils.h"
     33#include "pmModelClass.h"
     34#include "pmSourceMasks.h"
     35#include "pmSourceExtendedPars.h"
     36#include "pmSourceDiffStats.h"
     37#include "pmSource.h"
     38
    1839#include "pmPhotObj.h"
    19 #include "pmSource.h"
     40
    2041
    2142static void pmPhotObjFree (pmPhotObj *tmp)
    … …  
    82103    return (0);
    83104}
     105
     106// sort by X (ascending)
     107int pmPhotObjSortByX (const void **a, const void **b)
     108{
     109    pmPhotObj *objA = *(pmPhotObj **)a;
     110    pmPhotObj *objB = *(pmPhotObj **)b;
     111
     112    psF32 fA = objA->x;
     113    psF32 fB = objB->x;
     114
     115    psF32 diff = fA - fB;
     116    if (diff > FLT_EPSILON) return (+1);
     117    if (diff < FLT_EPSILON) return (-1);
     118    return (0);
     119}

  • trunk/psModules/src/objects/pmPhotObj.h

    r28013 r29004  
    1010# ifndef PM_PHOT_OBJ_H
    1111# define PM_PHOT_OBJ_H
    12 
    13 #include <pslib.h>
    14 #include "pmPeaks.h"
    15 #include "pmModel.h"
    16 #include "pmMoments.h"
    1712
    1813/// @addtogroup Objects Object Detection / Analysis Functions
    … …  
    4742
    4843int pmPhotObjSortBySN (const void **a, const void **b);
     44int pmPhotObjSortByX (const void **a, const void **b);
    4945
    5046/// @}

  • trunk/psModules/src/objects/pmSource.c

    r28013 r29004  
    2323#include "pmFPA.h"
    2424#include "pmFPAMaskWeight.h"
     25
     26#include "pmTrend2D.h"
     27#include "pmResiduals.h"
     28#include "pmGrowthCurve.h"
    2529#include "pmSpan.h"
     30#include "pmFootprintSpans.h"
    2631#include "pmFootprint.h"
    2732#include "pmPeaks.h"
    2833#include "pmMoments.h"
    29 #include "pmResiduals.h"
    30 #include "pmGrowthCurve.h"
    31 #include "pmTrend2D.h"
    32 #include "pmPSF.h"
     34#include "pmModelFuncs.h"
    3335#include "pmModel.h"
     36#include "pmModelUtils.h"
     37#include "pmModelClass.h"
     38#include "pmSourceMasks.h"
     39#include "pmSourceExtendedPars.h"
     40#include "pmSourceDiffStats.h"
    3441#include "pmSource.h"
    3542
    … …  
    98105    static int id = 1;
    99106    pmSource *source = (pmSource *) psAlloc(sizeof(pmSource));
    100     *(int *)&source->id = id++;
     107    P_PM_SOURCE_SET_ID(source, id++);
     108
    101109    source->seq = -1;
    102110    source->peak = NULL;
    … …  
    114122    source->type = PM_SOURCE_TYPE_UNKNOWN;
    115123    source->mode = PM_SOURCE_MODE_DEFAULT;
     124    source->mode2 = PM_SOURCE_MODE_DEFAULT;
    116125    source->tmpFlags = 0;
    117126    source->extpars = NULL;
    … …  
    131140    source->sky    = NAN;
    132141    source->skyErr = NAN;
    133     source->pixWeight = NAN;
     142    source->pixWeightNotBad = NAN;
     143    source->pixWeightNotPoor = NAN;
    134144
    135145    source->psfChisq = NAN;
    … …  
    142152
    143153/******************************************************************************
    144 pmSourceCopy(): copy the pmSource structure and contents
    145 XXX : are we OK with incrementing the ID?
     154pmSourceCopy(): copy the pmSource, yielding a copy of the source that can be used without
     155affecting the original.  This Copy can be used to allow multiple fit attempts on the same
     156object.  The pixels, variance, and mask arrays all point to the same original subarrays.  The
     157peak and moments point at the original values.
    146158*****************************************************************************/
    147159pmSource *pmSourceCopy(pmSource *in)
     160{
     161    if (in == NULL) {
     162        return(NULL);
     163    }
     164    pmSource *source = pmSourceAlloc ();
     165
     166    // keep the original ID so we can find map back to the original
     167    P_PM_SOURCE_SET_ID(source, in->id);
     168
     169    // peak has the same values as the original
     170    if (in->peak != NULL) {
     171        source->peak = pmPeakAlloc (in->peak->x, in->peak->y, in->peak->value, in->peak->type);
     172        source->peak->xf = in->peak->xf;
     173        source->peak->yf = in->peak->yf;
     174        source->peak->flux = in->peak->flux;
     175        source->peak->SN = in->peak->SN;
     176    }
     177
     178    // copy the values in the moments structure
     179    if (in->moments != NULL) {
     180        source->moments  =  pmMomentsAlloc();
     181        *source->moments = *in->moments;
     182    }
     183
     184    // These images are all views to the parent.  We want a new view, but pointing at the same
     185    // pixels.  Modifying these pixels (ie, subtracting the model) will affect the pixels seen
     186    // by all copies.
     187    source->pixels   = psImageCopyView(NULL, in->pixels);
     188    source->variance   = psImageCopyView(NULL, in->variance);
     189    source->maskView = in->maskView ? psImageCopyView(NULL, in->maskView) : NULL;
     190
     191    // the maskObj is a unique mask array; create a new mask image
     192    source->maskObj = in->maskObj ? psImageCopy (NULL, in->maskObj, PS_TYPE_IMAGE_MASK) : NULL;
     193
     194    source->type = in->type;
     195    source->mode = in->mode;
     196    source->imageID = in->imageID;
     197
     198    return(source);
     199}
     200
     201/******************************************************************************
     202pmSourceCopyData(): this creates a new, duplicate source with the same parameters as the
     203original (but is actually a new source at the same location)
     204*****************************************************************************/
     205pmSource *pmSourceCopyData(pmSource *in)
    148206{
    149207    if (in == NULL) {
    … …  
    482540        }
    483541        psfClump.X  = stats->clippedMean;
    484         psfClump.dX = stats->clippedStdev;
     542        psfClump.dX = hypot(stats->clippedStdev, PSF_CLUMP_GRID_SCALE);
    485543
    486544        if (!psVectorStats (stats, tmpSy, NULL, NULL, 0)) {
    … …  
    489547        }
    490548        psfClump.Y  = stats->clippedMean;
    491         psfClump.dY = stats->clippedStdev;
     549        psfClump.dY = hypot(stats->clippedStdev, PSF_CLUMP_GRID_SCALE);
    492550
    493551        psTrace ("psModules.objects", 2, "clump  X,  Y: %f, %f\n", psfClump.X, psfClump.Y);
    … …  
    910968    bool addNoise = mode & PM_MODEL_OP_NOISE;
    911969
    912     if (source->modelFlux) {
     970    // require the use of pmModelAddWithOffset if we are adding noise (because the model size and norm are rescaled)
     971    if (!addNoise && source->modelFlux) {
    913972        // add in the pixels from the modelFlux image
    914973        int dX = source->modelFlux->col0 - source->pixels->col0;
    … …  
    931990
    932991        psF32 **target = source->pixels->data.F32;
    933         if (addNoise) {
    934             // when adding noise, we assume the shape and Io have been modified
    935             target = source->variance->data.F32;
    936         }
    937992
    938993        for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
    … …  
    9491004            }
    9501005        }
    951         if (!addNoise) {
    952             if (add) {
    953                 source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
    954             } else {
    955                 source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
    956             }
     1006        if (add) {
     1007            source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
     1008        } else {
     1009            source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
    9571010        }
    9581011        return true;
    … …  
    9731026        }
    9741027    }
     1028
     1029    return true;
     1030}
     1031
     1032// should we call pmSourceCacheModel if it does not exist?
     1033bool pmSourceNoiseOp (pmSource *source, pmModelOpMode mode, float FACTOR, float SIZE, bool add, psImageMaskType maskVal, int dx, int dy)
     1034{
     1035    assert (mode & PM_MODEL_OP_NOISE);
     1036    PS_ASSERT_PTR_NON_NULL(source, false);
     1037    PS_ASSERT_PTR_NON_NULL(source->peak, false);
     1038
     1039    if (add) {
     1040        psTrace ("psphot", 3, "adding noise to object at %f,%f\n", source->peak->xf, source->peak->yf);
     1041    } else {
     1042        psTrace ("psphot", 3, "removing noise from object at %f,%f\n", source->peak->xf, source->peak->yf);
     1043    }
     1044
     1045    pmSourceNoiseOpModel (source->modelPSF, source, mode, FACTOR, SIZE, add, maskVal, dx, dy);
     1046
     1047    if (source->modelEXT) {
     1048        pmSourceNoiseOpModel (source->modelEXT, source, mode, FACTOR, SIZE, add, maskVal, dx, dy);
     1049    }
     1050
     1051    return true;
     1052}
     1053
     1054bool pmSourceNoiseOpModel (pmModel *model, pmSource *source, pmModelOpMode mode, float FACTOR, float SIZE, bool add, psImageMaskType maskVal, int dx, int dy)
     1055{
     1056    bool status;
     1057    psEllipseShape oldshape;
     1058    psEllipseShape newshape;
     1059    psEllipseAxes axes;
     1060
     1061    if (add) {
     1062        psTrace ("psphot", 4, "adding noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
     1063    } else {
     1064        psTrace ("psphot", 4, "remove noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
     1065    }
     1066
     1067    psF32 *PAR = model->params->data.F32;
     1068
     1069    // save original values
     1070    float oldI0  = PAR[PM_PAR_I0];
     1071    oldshape.sx  = PAR[PM_PAR_SXX];
     1072    oldshape.sy  = PAR[PM_PAR_SYY];
     1073    oldshape.sxy = PAR[PM_PAR_SXY];
     1074
     1075    // XXX can this be done more intelligently?
     1076    if (oldI0 == 0.0) return false;
     1077    if (!isfinite(oldI0)) return false;
     1078
     1079    // increase size and height of source
     1080    axes = psEllipseShapeToAxes (oldshape, 20.0);
     1081    axes.major *= SIZE;
     1082    axes.minor *= SIZE;
     1083    newshape = psEllipseAxesToShape (axes);
     1084    PAR[PM_PAR_I0]  = FACTOR*oldI0;
     1085    PAR[PM_PAR_SXX] = newshape.sx;
     1086    PAR[PM_PAR_SYY] = newshape.sy;
     1087    PAR[PM_PAR_SXY] = newshape.sxy;
     1088
     1089    psImage *target = source->variance;
     1090
     1091    if (add) {
     1092        status = pmModelAddWithOffset (target, source->maskObj, model, mode, maskVal, dx, dy);
     1093    } else {
     1094        status = pmModelSubWithOffset (target, source->maskObj, model, mode, maskVal, dx, dy);
     1095    }
     1096
     1097    // restore original values
     1098    PAR[PM_PAR_I0]  = oldI0;
     1099    PAR[PM_PAR_SXX] = oldshape.sx;
     1100    PAR[PM_PAR_SYY] = oldshape.sy;
     1101    PAR[PM_PAR_SXY] = oldshape.sxy;
    9751102
    9761103    return true;

  • trunk/psModules/src/objects/pmSource.h

    r28013 r29004  
    1111# define PM_SOURCE_H
    1212
    13 #include <pslib.h>
    14 #include "pmPeaks.h"
    15 #include "pmModel.h"
    16 #include "pmMoments.h"
    17 #include "pmSourceExtendedPars.h"
    18 #include "pmSourceDiffStats.h"
    19 
    2013/// @addtogroup Objects Object Detection / Analysis Functions
    2114/// @{
    22 
    23 #include <pmSourceMasks.h>
    2415
    2516/** pmSourceType enumeration
    … …  
    7465    pmSourceType type;                  ///< Best identification of object.
    7566    pmSourceMode mode;                  ///< analysis flags set for object.
     67    pmSourceMode2 mode2;                ///< analysis flags set for object.
    7668    pmSourceTmpF tmpFlags;              ///< internal-only flags
    7769    psArray *blends;                    ///< collection of sources thought to be confused with object
    … …  
    8274    float errMag;                       ///< error in psfMag OR extMag (depending on type)
    8375    float apMag;                        ///< apMag corresponding to psfMag or extMag (depending on type)
    84     float pixWeight;                    ///< model-weighted coverage of valid pixels
     76    float apMagRaw;                     ///< raw mag in given aperture
     77    float apRadius;                     ///< radius for aperture magnitude
     78
     79    float pixWeightNotBad;              ///< PSF-weighted coverage of unmasked (not BAD) pixels
     80    float pixWeightNotPoor;             ///< PSF-weighted coverage of unmasked (not POOR) pixels
     81
    8582    float psfChisq;                     ///< probability of PSF
    8683    float crNsigma;                     ///< Nsigma deviation from PSF to CR
    8784    float extNsigma;                    ///< Nsigma deviation from PSF to EXT
    8885    float sky, skyErr;                  ///< The sky and its error at the center of the object
    89     float apRadius;
    9086    psRegion region;                    ///< area on image covered by selected pixels
    9187    pmSourceExtendedPars *extpars;      ///< extended source parameters
    … …  
    113109pmPSFClump;
    114110
     111// private macro to set the source ID (a const)
     112#define P_PM_SOURCE_SET_ID(S,V) { *(int *)&(S)->id = (V); }
    115113
    116114/** pmSourceAlloc()
    … …  
    126124
    127125pmSource  *pmSourceCopy(pmSource *source);
     126pmSource *pmSourceCopyData(pmSource *in);
    128127
    129128// free just the pixels for a source, keeping derived data
    … …  
    242241bool pmSourceSubWithOffset (pmSource *source, pmModelOpMode mode, psImageMaskType maskVal, int dx, int dy);
    243242
     243bool pmSourceNoiseOpModel (pmModel *model, pmSource *source, pmModelOpMode mode, float FACTOR, float SIZE, bool add, psImageMaskType maskVal, int dx, int dy);
     244bool pmSourceNoiseOp (pmSource *source, pmModelOpMode mode, float FACTOR, float SIZE, bool add, psImageMaskType maskVal, int dx, int dy);
     245
    244246bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add, psImageMaskType maskVal, int dx, int dy);
    245247bool pmSourceCacheModel (pmSource *source, psImageMaskType maskVal);

  • trunk/psModules/src/objects/pmSourceContour.c

    r24887 r29004  
    2323#include "pmHDU.h"
    2424#include "pmFPA.h"
     25
     26#include "pmTrend2D.h"
     27#include "pmResiduals.h"
     28#include "pmGrowthCurve.h"
    2529#include "pmSpan.h"
     30#include "pmFootprintSpans.h"
    2631#include "pmFootprint.h"
    2732#include "pmPeaks.h"
    2833#include "pmMoments.h"
    29 #include "pmResiduals.h"
    30 #include "pmGrowthCurve.h"
    31 #include "pmTrend2D.h"
    32 #include "pmPSF.h"
     34#include "pmModelFuncs.h"
    3335#include "pmModel.h"
     36#include "pmModelUtils.h"
     37#include "pmModelClass.h"
     38#include "pmSourceMasks.h"
     39#include "pmSourceExtendedPars.h"
     40#include "pmSourceDiffStats.h"
    3441#include "pmSource.h"
     42
    3543#include "pmSourceContour.h"
    3644

  • trunk/psModules/src/objects/pmSourceDiffStats.c

    r27531 r29004  
    2929    diffStats->nRatioAll = NAN;
    3030    diffStats->nGood = 0;
     31
     32    diffStats->SNp = NAN;
     33    diffStats->SNm = NAN;
     34    diffStats->Rp = NAN;
     35    diffStats->Rm = NAN;
    3136}
    3237

  • trunk/psModules/src/objects/pmSourceDiffStats.h

    r27531 r29004  
    2929    float nRatioAll;                    // = nGood / (nGood + nMask + nBad)
    3030    int   nGood;                        // nGood as defined above
     31    float SNp;                          // S/N of matched source in positive image
     32    float SNm;                          // S/N of matched source in negative image
     33    float Rp;                           // radius of matched source in positive image
     34    float Rm;                           // radius of matched source in negative image
    3135} pmSourceDiffStats;
    3236

  • trunk/psModules/src/objects/pmSourceFitModel.c

    r26070 r29004  
    2323#include "pmHDU.h"
    2424#include "pmFPA.h"
     25
     26#include "pmTrend2D.h"
     27#include "pmResiduals.h"
     28#include "pmGrowthCurve.h"
    2529#include "pmSpan.h"
     30#include "pmFootprintSpans.h"
    2631#include "pmFootprint.h"
    2732#include "pmPeaks.h"
    2833#include "pmMoments.h"
    29 #include "pmGrowthCurve.h"
    30 #include "pmResiduals.h"
    31 #include "pmTrend2D.h"
    32 #include "pmPSF.h"
     34#include "pmModelFuncs.h"
    3335#include "pmModel.h"
     36#include "pmModelUtils.h"
     37#include "pmModelClass.h"
     38#include "pmSourceMasks.h"
     39#include "pmSourceExtendedPars.h"
     40#include "pmSourceDiffStats.h"
    3441#include "pmSource.h"
    35 #include "pmModelClass.h"
    3642#include "pmSourceFitModel.h"
    3743
    38 // save as static values so they may be set externally
    39 static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15;
    40 static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1;
    41 static psF32 PM_SOURCE_FIT_MODEL_WEIGHT = 1.0;
    42 static bool  PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = true;
    43 
    44 bool pmSourceFitModelInit (float nIter, float tol, float weight, bool poissonErrors)
     44void pmSourceFitOptionsFree(pmSourceFitOptions *opt)
    4545{
    46 
    47     PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = nIter;
    48     PM_SOURCE_FIT_MODEL_TOLERANCE = tol;
    49     PM_SOURCE_FIT_MODEL_WEIGHT = weight;
    50     PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = poissonErrors;
    51 
    52     return true;
     46    return;
     47}
     48
     49pmSourceFitOptions *pmSourceFitOptionsAlloc(void) {
     50
     51    pmSourceFitOptions *opt = (pmSourceFitOptions *) psAlloc(sizeof(pmSourceFitOptions));
     52    psMemSetDeallocator(opt, (psFreeFunc) pmSourceFitOptionsFree);
     53
     54    opt->mode = PM_SOURCE_FIT_PSF;
     55    opt->nIter  = 15;
     56    opt->minTol = 0.01;
     57    opt->maxTol = 1.00;
     58    opt->weight = 1.00;
     59    opt->maxChisqDOF = NAN;
     60    opt->poissonErrors = true;
     61
     62    return opt;
    5363}
    5464
    5565bool pmSourceFitModel (pmSource *source,
    5666                       pmModel *model,
    57                        pmSourceFitMode mode,
     67                       pmSourceFitOptions *options,
    5868                       psImageMaskType maskVal)
    5969{
    … …  
    7686    psVector *yErr = psVectorAllocEmpty(nPix, PS_TYPE_F32);
    7787
     88    // XXX for a test, skip the central pixel in the sersic fit
     89    bool skipCenter = false && (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
     90    float Xo = model->params->data.F32[PM_PAR_XPOS];
     91    float Yo = model->params->data.F32[PM_PAR_YPOS];
     92
    7893    // fill in the coordinate and value entries
    7994    nPix = 0;
    … …  
    95110            // skip nan values in image
    96111            if (!isfinite(source->variance->data.F32[i][j])) {
    97               fprintf (stderr, "impossible! %x vs %x\n", source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[i][j], maskVal);
    98               continue;
    99             }
    100 
    101             psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
     112                fprintf (stderr, "impossible! %x vs %x\n", source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[i][j], maskVal);
     113                continue;
     114            }
    102115
    103116            // Convert i/j to image space:
    104117            // 0.5 PIX: the coordinate values must be in pixel coords, not index           
    105             coord->data.F32[0] = (psF32) (j + 0.5 + source->pixels->col0);
    106             coord->data.F32[1] = (psF32) (i + 0.5 + source->pixels->row0);
     118            float Xv = (psF32) (j + 0.5 + source->pixels->col0);
     119            float Yv = (psF32) (i + 0.5 + source->pixels->row0);
     120
     121            // XXX possible skip of center pixel:
     122            if (skipCenter) {
     123                float r = hypot(Xv - Xo, Yv - Yo);
     124                if (r < 0.75) {
     125                    continue;
     126                }
     127            }
     128
     129            psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
     130            coord->data.F32[0] = Xv;
     131            coord->data.F32[1] = Yv;
    107132            x->data[nPix] = (psPtr *) coord;
    108133            y->data.F32[nPix] = source->pixels->data.F32[i][j];
    … …  
    111136            // as variance to avoid the bias from systematic errors here we would just use the
    112137            // source sky variance
    113             if (PM_SOURCE_FIT_MODEL_PIX_WEIGHTS) {
     138            if (options->poissonErrors) {
    114139                yErr->data.F32[nPix] = 1.0 / source->variance->data.F32[i][j];
    115140            } else {
    116                 yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT;
     141                yErr->data.F32[nPix] = 1.0 / options->weight;
    117142            }
    118143            nPix++;
    … …  
    133158    // set parameter mask based on fitting mode
    134159    int nParams = 0;
    135     switch (mode) {
    136     case PM_SOURCE_FIT_NORM:
     160    switch (options->mode) {
     161      case PM_SOURCE_FIT_NORM:
    137162        // NORM-only model fits only source normalization (Io)
    138163        nParams = 1;
    … …  
    140165        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
    141166        break;
    142     case PM_SOURCE_FIT_PSF:
     167      case PM_SOURCE_FIT_PSF:
    143168        // PSF model only fits x,y,Io
    144169        nParams = 3;
    … …  
    148173        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 0;
    149174        break;
    150     case PM_SOURCE_FIT_EXT:
     175      case PM_SOURCE_FIT_EXT:
    151176        // EXT model fits all params (except sky)
    152177        nParams = params->n - 1;
    … …  
    154179        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
    155180        break;
    156     default:
    157         psAbort("invalid fitting mode");
     181      case PM_SOURCE_FIT_INDEX:
     182        // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
     183        psVectorInit (constraint->paramMask, 1);
     184        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
     185        if (params->n == 7) {
     186            nParams = 1;
     187        } else {
     188            nParams = 2;
     189            constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 0;
     190        }
     191        break;
     192      case PM_SOURCE_FIT_NO_INDEX:
     193        // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
     194        psVectorInit (constraint->paramMask, 0);
     195        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
     196        if (params->n == 7) {
     197            nParams = params->n - 1;
     198        } else {
     199            nParams = params->n - 2;
     200            constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1;
     201        }
     202        break;
     203      default:
     204        psAbort("invalid fitting mode");
    158205    }
    159206    // force the floating parameters to fall within the contraint ranges
    160207    for (int i = 0; i < params->n; i++) {
    161         model->modelLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
    162         model->modelLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
     208        model->modelLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
     209        model->modelLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
    163210    }
    164211
    … …  
    173220    }
    174221
    175     psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE);
     222    psMinimization *myMin = psMinimizationAlloc (options->nIter, options->minTol, options->maxTol);
    176223
    177224    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
    … …  
    194241    model->flags |= PM_MODEL_STATUS_FITTED;
    195242    if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE;
     243    if (myMin->lastDelta > myMin->minTol) model->flags |= PM_MODEL_STATUS_WEAK_FIT;
    196244
    197245    // get the Gauss-Newton distance for fixed model parameters

  • trunk/psModules/src/objects/pmSourceFitModel.h

    r21183 r29004  
    1919    PM_SOURCE_FIT_EXT,
    2020    PM_SOURCE_FIT_PSF_AND_SKY,
    21     PM_SOURCE_FIT_EXT_AND_SKY
     21    PM_SOURCE_FIT_EXT_AND_SKY,
     22    PM_SOURCE_FIT_INDEX,
     23    PM_SOURCE_FIT_NO_INDEX,
    2224} pmSourceFitMode;
    2325
    24 bool pmSourceFitModelInit(
    25     float nIter,   ///< max number of allowed iterations
    26     float tol,      ///< convergence criterion
    27     float weight,      ///< use this weight for constant-weight fits
    28     bool poissonErrors   // use poisson errors for fits?
    29 );
     26typedef struct {
     27    pmSourceFitMode mode;               ///< optionally fit all or a subset of parameters
     28    float nIter;                        ///< max number of allowed iterations
     29    float minTol;                       ///< convergence criterion
     30    float maxTol;                       ///< convergence criterion
     31    float maxChisqDOF;                  ///< convergence criterion
     32    float weight;                       ///< use this weight for constant-weight fits
     33    bool poissonErrors;                 ///< use poisson errors for fits?
     34} pmSourceFitOptions;
     35
     36// the pmSourceFitOptions structure is used to control details of the fitting process
     37pmSourceFitOptions *pmSourceFitOptionsAlloc(void);
     38
     39// bool pmSourceFitModelInit(
     40//     pmSourceFitMode mode,            ///< what parameter set should be fitted?
     41//     float nIter,                     ///< max number of allowed iterations
     42//     float tol,                               ///< convergence criterion
     43//     float weight,                    ///< use this weight for constant-weight fits
     44//     bool poissonErrors                       ///< use poisson errors for fits?
     45// );
    3046
    3147/** pmSourceFitModel()
    … …  
    3854 */
    3955bool pmSourceFitModel(
    40     pmSource *source,   ///< The input pmSource
    41     pmModel *model,   ///< model to be fitted
    42     pmSourceFitMode mode,  ///< define parameters to be fitted
     56    pmSource *source,                   ///< The input pmSource
     57    pmModel *model,                     ///< model to be fitted
     58    pmSourceFitOptions *options,        ///< define parameters to be fitted
    4359    psImageMaskType maskVal             ///< Value to mask
    4460);
    4561
    46 
    47 // initialize data for a group of object models
    48 bool pmSourceFitSetInit (pmModelType type);
    49 
    50 // clear data for a group of object models
    51 void pmSourceFitSetClear (void);
    52 
    53 // function used to set limits for a group of models
    54 bool pmSourceFitSet_CheckLimits (psMinConstraintMode mode, int nParam, float *params, float *betas);
    55 
    56 // function used to fit a group of object models
    57 psF32 pmSourceFitSet_Function(psVector *deriv,
    58                               const psVector *params,
    59                               const psVector *x);
    60 
    61 /** pmSourceFitSet()
    62  *
    63  * Fit the requested model to the specified source. The starting guess for the model is given
    64  * by the input source.model parameter values. The pixels of interest are specified by the
    65  * source.pixels and source.mask entries. This function calls psMinimizeLMChi2() on the image
    66  * data. The function returns TRUE on success or FALSE on failure.
    67  *
    68  */
    69 bool pmSourceFitSet(
    70     pmSource *source,   ///< The input pmSource
    71     psArray *modelSet,   ///< model to be fitted
    72     pmSourceFitMode mode,  ///< define parameters to be fitted
    73     psImageMaskType maskVal             ///< Vale to mask
    74 
    75 );
     62// // initialize data for a group of object models
     63// bool pmSourceFitSetInit (pmModelType type);
     64//
     65// // clear data for a group of object models
     66// void pmSourceFitSetClear (void);
     67//
     68// // function used to set limits for a group of models
     69// bool pmSourceFitSet_CheckLimits (psMinConstraintMode mode, int nParam, float *params, float *betas);
     70//
     71// // function used to fit a group of object models
     72// psF32 pmSourceFitSet_Function(psVector *deriv,
     73//                               const psVector *params,
     74//                               const psVector *x);
     75//
     76// /** pmSourceFitSet()
     77//  *
     78//  * Fit the requested model to the specified source. The starting guess for the model is given
     79//  * by the input source.model parameter values. The pixels of interest are specified by the
     80//  * source.pixels and source.mask entries. This function calls psMinimizeLMChi2() on the image
     81//  * data. The function returns TRUE on success or FALSE on failure.
     82//  *
     83//  */
     84// bool pmSourceFitSet(
     85//     pmSource *source,   ///< The input pmSource
     86//     psArray *modelSet,   ///< model to be fitted
     87//     pmSourceFitMode mode,  ///< define parameters to be fitted
     88//     psImageMaskType maskVal          ///< Vale to mask
     89//
     90// );
    7691
    7792/// @}

  • trunk/psModules/src/objects/pmSourceFitSet.c

    r27903 r29004  
    2222#include "pmHDU.h"
    2323#include "pmFPA.h"
     24
     25#include "pmTrend2D.h"
     26#include "pmResiduals.h"
     27#include "pmGrowthCurve.h"
    2428#include "pmSpan.h"
     29#include "pmFootprintSpans.h"
    2530#include "pmFootprint.h"
    2631#include "pmPeaks.h"
    2732#include "pmMoments.h"
    28 #include "pmGrowthCurve.h"
    29 #include "pmResiduals.h"
    30 #include "pmTrend2D.h"
    31 #include "pmPSF.h"
     33#include "pmModelFuncs.h"
    3234#include "pmModel.h"
     35#include "pmModelUtils.h"
     36#include "pmModelClass.h"
     37#include "pmSourceMasks.h"
     38#include "pmSourceExtendedPars.h"
     39#include "pmSourceDiffStats.h"
    3340#include "pmSource.h"
    34 #include "pmModelClass.h"
     41
    3542#include "pmSourceFitModel.h"
    3643#include "pmSourceFitSet.h"
    3744
    3845// save as static values so they may be set externally
    39 static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15;
    40 static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1;
    41 static psF32 PM_SOURCE_FIT_MODEL_WEIGHT = 1.0;
    42 static bool  PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = true;
     46// static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15;
     47// static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1;
     48// static psF32 PM_SOURCE_FIT_MODEL_WEIGHT = 1.0;
     49// static bool  PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = true;
    4350
    4451/********************* Source Model Set Functions ***************************/
    … …  
    429436bool pmSourceFitSet (pmSource *source,
    430437                     psArray *modelSet,
    431                      pmSourceFitMode mode,
     438                     pmSourceFitOptions *options,
    432439                     psImageMaskType maskVal)
    433440{
    … …  
    478485            // as variance to avoid the bias from systematic errors here we would just use the
    479486            // source sky variance
    480             if (PM_SOURCE_FIT_MODEL_PIX_WEIGHTS) {
    481                 yErr->data.F32[nPix] = 1.0 / source->variance->data.F32[i][j];
    482             } else {
    483                 yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT;
    484             }
    485             nPix++;
    486         }
     487            if (options->poissonErrors) {
     488                yErr->data.F32[nPix] = 1.0 / source->variance->data.F32[i][j];
     489            } else {
     490                yErr->data.F32[nPix] = 1.0 / options->weight;
     491            }
     492            nPix++;
     493        }
    487494    }
    488495    x->n = nPix;
    … …  
    490497    yErr->n = nPix;
    491498
    492     // create the FitSet for this thread and set the initial parameter guesses
     499// create the FitSet for this thread and set the initial parameter guesses
    493500    pmSourceFitSetData *thisSet = pmSourceFitSetDataSet(modelSet);
    494501
    495     // define param and deriv vectors for complete set of parameters
     502// define param and deriv vectors for complete set of parameters
    496503    psVector *params = psVectorAlloc (thisSet->nParamSet, PS_TYPE_F32);
    497504
    498     // set the param and deriv vectors based on the curent values
     505// set the param and deriv vectors based on the curent values
    499506    pmSourceFitSetJoin (NULL, params, thisSet);
    500507
    501     // create the minimization constraints
     508// create the minimization constraints
    502509    psMinConstraint *constraint = psMinConstraintAlloc();
    503510    constraint->paramMask = psVectorAlloc (thisSet->nParamSet, PS_TYPE_VECTOR_MASK);
    504511    constraint->checkLimits = pmSourceFitSetCheckLimits;
    505512
    506     pmSourceFitSetMasks (constraint, thisSet, mode);
    507 
    508     // force the floating parameters to fall within the contraint ranges
     513    pmSourceFitSetMasks (constraint, thisSet, options->mode);
     514
     515// force the floating parameters to fall within the contraint ranges
    509516    for (int i = 0; i < params->n; i++) {
    510         pmSourceFitSetCheckLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
    511         pmSourceFitSetCheckLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
     517        pmSourceFitSetCheckLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
     518        pmSourceFitSetCheckLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
    512519    }
    513520
    514521    if (psTraceGetLevel("psModules.objects") >= 5) {
    515         for (int i = 0; i < params->n; i++) {
    516             fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]);
    517         }
     522        for (int i = 0; i < params->n; i++) {
     523            fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]);
     524        }
    518525    }
    519526
    520527    if (nPix <  thisSet->nParamSet + 1) {
    521         psTrace (__func__, 4, "insufficient valid pixels\n");
    522         psTrace("psModules.objects", 10, "---- %s() end : fail pixels ----\n", __func__);
    523         for (int i = 0; i < modelSet->n; i++) {
    524             pmModel *model = modelSet->data[i];
    525             model->flags |= PM_MODEL_STATUS_BADARGS;
    526         }
    527         psFree (x);
    528         psFree (y);
    529         psFree (yErr);
    530         psFree (params);
    531         psFree(constraint);
    532         pmSourceFitSetDataClear(); // frees thisSet and removes if from the array of fitSets
    533         return(false);
    534     }
    535 
    536     psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE);
     528        psTrace (__func__, 4, "insufficient valid pixels\n");
     529        psTrace("psModules.objects", 10, "---- %s() end : fail pixels ----\n", __func__);
     530        for (int i = 0; i < modelSet->n; i++) {
     531            pmModel *model = modelSet->data[i];
     532            model->flags |= PM_MODEL_STATUS_BADARGS;
     533        }
     534        psFree (x);
     535        psFree (y);
     536        psFree (yErr);
     537        psFree (params);
     538        psFree(constraint);
     539        pmSourceFitSetDataClear(); // frees thisSet and removes if from the array of fitSets
     540        return(false);
     541    }
     542
     543    psMinimization *myMin = psMinimizationAlloc (options->nIter, options->minTol, options->maxTol);
    537544
    538545    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
    … …  
    540547    fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, pmSourceFitSetFunction);
    541548    if (!fitStatus) {
    542         psTrace("psModules.objects", 4, "Failed to fit model (%ld components)\n", modelSet->n);
    543     }
    544 
    545     // parameter errors from the covariance matrix
     549        psTrace("psModules.objects", 4, "Failed to fit model (%ld components)\n", modelSet->n);
     550    }
     551
     552// parameter errors from the covariance matrix
    546553    psVector *dparams = psVectorAlloc (thisSet->nParamSet, PS_TYPE_F32);
    547554    for (int i = 0; i < dparams->n; i++) {
    548         if ((constraint->paramMask != NULL) && constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i])
    549             continue;
    550         dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
    551     }
    552 
    553     // get the Gauss-Newton distance for fixed model parameters
     555        if ((constraint->paramMask != NULL) && constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i])
     556            continue;
     557        dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
     558    }
     559
     560// get the Gauss-Newton distance for fixed model parameters
    554561    if (constraint->paramMask != NULL) {
    555         psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32);
    556         psVector *altmask = psVectorAlloc (params->n, PS_TYPE_VECTOR_MASK);
    557         altmask->data.PS_TYPE_VECTOR_MASK_DATA[0] = 1;
    558         for (int i = 1; i < dparams->n; i++) {
    559             altmask->data.PS_TYPE_VECTOR_MASK_DATA[i] = (constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) ? 0 : 1;
    560         }
    561         psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, pmSourceFitSetFunction);
    562 
    563         for (int i = 0; i < dparams->n; i++) {
    564             if (!constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i])
    565                 continue;
    566             // note that delta is the value *subtracted* from the parameter
    567             // to get the new guess.  for dparams to represent the direction
    568             // of motion, we need to take -delta
    569             dparams->data.F32[i] = -delta->data.F32[i];
    570         }
    571         psFree (delta);
    572         psFree (altmask);
     562        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32);
     563        psVector *altmask = psVectorAlloc (params->n, PS_TYPE_VECTOR_MASK);
     564        altmask->data.PS_TYPE_VECTOR_MASK_DATA[0] = 1;
     565        for (int i = 1; i < dparams->n; i++) {
     566            altmask->data.PS_TYPE_VECTOR_MASK_DATA[i] = (constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) ? 0 : 1;
     567        }
     568        psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, pmSourceFitSetFunction);
     569
     570        for (int i = 0; i < dparams->n; i++) {
     571            if (!constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i])
     572                continue;
     573            // note that delta is the value *subtracted* from the parameter
     574            // to get the new guess.  for dparams to represent the direction
     575            // of motion, we need to take -delta
     576            dparams->data.F32[i] = -delta->data.F32[i];
     577        }
     578        psFree (delta);
     579        psFree (altmask);
    573580    }
    574581

  • trunk/psModules/src/objects/pmSourceFitSet.h

    r23487 r29004  
    5656    pmSource *source,                   ///< The input pmSource
    5757    psArray *modelSet,                  ///< model to be fitted
    58     pmSourceFitMode mode,               ///< define parameters to be fitted
     58    pmSourceFitOptions *options,        ///< define options for fitting process
    5959    psImageMaskType maskVal             ///< Vale to mask
    6060

  • trunk/psModules/src/objects/pmSourceGroups.c

    r26450 r29004  
    66#include <pslib.h>
    77
     8#include "pmHDU.h"
    89#include "pmFPA.h"
     10
     11#include "pmTrend2D.h"
     12#include "pmResiduals.h"
     13#include "pmGrowthCurve.h"
     14#include "pmSpan.h"
     15#include "pmFootprintSpans.h"
     16#include "pmFootprint.h"
     17#include "pmPeaks.h"
     18#include "pmMoments.h"
     19#include "pmModelFuncs.h"
     20#include "pmModel.h"
     21#include "pmModelUtils.h"
     22#include "pmModelClass.h"
     23#include "pmSourceMasks.h"
     24#include "pmSourceExtendedPars.h"
     25#include "pmSourceDiffStats.h"
    926#include "pmSource.h"
     27#include "pmSourceFitModel.h"
     28#include "pmPSF.h"
     29#include "pmPSFtry.h"
     30#include "pmDetections.h"
     31
    1032#include "pmSourceGroups.h"
    1133

  • trunk/psModules/src/objects/pmSourceGroups.h

    r26182 r29004  
    11#ifndef PM_SOURCE_GROUPS_H
    22#define PM_SOURCE_GROUPS_H
    3 
    4 #include <pslib.h>
    5 
    6 #include "pmFPA.h"
    73
    84/// Groups of sources

  • trunk/psModules/src/objects/pmSourceIO.c

    r28013 r29004  
    3131#include "pmConceptsRead.h"
    3232
     33#include "pmTrend2D.h"
     34#include "pmResiduals.h"
     35#include "pmGrowthCurve.h"
    3336#include "pmSpan.h"
     37#include "pmFootprintSpans.h"
    3438#include "pmFootprint.h"
    3539#include "pmPeaks.h"
    3640#include "pmMoments.h"
    37 #include "pmGrowthCurve.h"
    38 #include "pmResiduals.h"
    39 #include "pmTrend2D.h"
     41#include "pmModelFuncs.h"
     42#include "pmModel.h"
     43#include "pmModelUtils.h"
     44#include "pmModelClass.h"
     45#include "pmSourceMasks.h"
     46#include "pmSourceExtendedPars.h"
     47#include "pmSourceDiffStats.h"
     48#include "pmSource.h"
     49#include "pmSourceFitModel.h"
    4050#include "pmPSF.h"
    41 #include "pmModel.h"
     51#include "pmPSFtry.h"
     52
    4253#include "pmDetections.h"
    43 #include "pmSource.h"
    44 #include "pmModelClass.h"
    4554#include "pmDetEff.h"
    4655#include "pmSourceIO.h"
    … …  
    326335}
    327336
     337# define PM_SOURCES_WRITE(NAME,TYPE)                                    \
     338    if (!strcmp (exttype, NAME)) {                                      \
     339        status &= pmSourcesWrite_##TYPE(file->fits, readout, sources, file->header, outhead, dataname, recipe); \
     340        if (xsrcname) {                                                 \
     341            status &= pmSourcesWrite_##TYPE##_XSRC(file->fits, readout, sources, file->header, xsrcname, recipe); \
     342        }                                                               \
     343        if (xfitname) {                                                 \
     344            status &= pmSourcesWrite_##TYPE##_XFIT (file->fits, readout, sources, file->header, xfitname); \
     345        }                                                               \
     346    }
     347
    328348// write out all readout-level Objects files for this cell
    329349bool pmReadoutWriteObjects (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
    … …  
    360380    }
    361381
     382    // the older types (RAW, OBJ, SX, CMP) are for backwards compatibility -- deprecate eventually?
    362383    switch (file->type) {
    363384      case PM_FPA_FILE_RAW:
    … …  
    518539                psMetadataAddStr (outhead, PS_LIST_TAIL, "XFITNAME", PS_META_REPLACE, "name of XFIT table extension", xfitname);
    519540            }
    520 
    521             // XXX these are case-sensitive since the EXTYPE is case-sensitive
     541   
     542
     543            // these are case-sensitive since the EXTYPE is case-sensitive
    522544            status = true;
    523             if (!strcmp (exttype, "SMPDATA")) {
    524                 status &= pmSourcesWrite_SMPDATA (file->fits, sources, file->header, outhead, dataname);
    525             }
    526             if (!strcmp (exttype, "PS1_DEV_0")) {
    527                 status &= pmSourcesWrite_PS1_DEV_0 (file->fits, sources, file->header, outhead, dataname);
    528             }
    529             if (!strcmp (exttype, "PS1_DEV_1")) {
    530                 status &= pmSourcesWrite_PS1_DEV_1 (file->fits, sources, file->header, outhead, dataname);
    531             }
    532             if (!strcmp (exttype, "PS1_CAL_0")) {
    533                 status &= pmSourcesWrite_PS1_CAL_0 (file->fits, readout, sources, file->header, outhead, dataname);
    534             }
    535             if (!strcmp (exttype, "PS1_V1")) {
    536                 status &= pmSourcesWrite_CMF_PS1_V1 (file->fits, readout, sources, file->header, outhead, dataname);
    537             }
    538             if (!strcmp (exttype, "PS1_V2")) {
    539                 status &= pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname);
    540             }
    541             if (!strcmp (exttype, "PS1_SV1")) {
    542                 status &= pmSourcesWrite_CMF_PS1_SV1 (file->fits, readout, sources, file->header, outhead, dataname, recipe);
    543             }
    544             if (!strcmp (exttype, "PS1_DV1")) {
    545                 status &= pmSourcesWrite_CMF_PS1_DV1 (file->fits, readout, sources, file->header, outhead, dataname);
    546             }
    547 
    548             if (xsrcname) {
    549                 if (!strcmp (exttype, "PS1_DEV_1")) {
    550                     status &= pmSourcesWrite_PS1_DEV_1_XSRC (file->fits, sources, xsrcname, recipe);
    551                 }
    552                 if (!strcmp (exttype, "PS1_CAL_0")) {
    553                     status &= pmSourcesWrite_PS1_CAL_0_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe);
    554                 }
    555                 if (!strcmp (exttype, "PS1_V1")) {
    556                     status &= pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe);
    557                 }
    558                 if (!strcmp (exttype, "PS1_V2")) {
    559                     status &= pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe);
    560                 }
    561                 if (!strcmp (exttype, "PS1_SV1")) {
    562                     status &= pmSourcesWrite_CMF_PS1_SV1_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe);
    563                 }
    564                 if (!strcmp (exttype, "PS1_DV1")) {
    565                     status &= pmSourcesWrite_CMF_PS1_DV1_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe);
    566                 }
    567             }
    568             if (xfitname) {
    569                 if (!strcmp (exttype, "PS1_DEV_1")) {
    570                     status &= pmSourcesWrite_PS1_DEV_1_XFIT (file->fits, sources, xfitname);
    571                 }
    572                 if (!strcmp (exttype, "PS1_CAL_0")) {
    573                     status &= pmSourcesWrite_PS1_CAL_0_XFIT (file->fits, readout, sources, file->header, xfitname);
    574                 }
    575                 if (!strcmp (exttype, "PS1_V1")) {
    576                     status &= pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, readout, sources, xfitname);
    577                 }
    578                 if (!strcmp (exttype, "PS1_V2")) {
    579                     status &= pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, readout, sources, xfitname);
    580                 }
    581                 if (!strcmp (exttype, "PS1_SV1")) {
    582                     status &= pmSourcesWrite_CMF_PS1_SV1_XFIT (file->fits, readout, sources, xfitname);
    583                 }
    584                 if (!strcmp (exttype, "PS1_DV1")) {
    585                     status &= pmSourcesWrite_CMF_PS1_DV1_XFIT (file->fits, readout, sources, xfitname);
    586                 }
    587             }
     545            PM_SOURCES_WRITE("SMPDATA",   SMPDATA);
     546            PM_SOURCES_WRITE("PS1_DEV_0", PS1_DEV_0);
     547            PM_SOURCES_WRITE("PS1_DEV_1", PS1_DEV_1);
     548            PM_SOURCES_WRITE("PS1_CAL_0", PS1_CAL_0);
     549            PM_SOURCES_WRITE("PS1_V1",    CMF_PS1_V1);
     550            PM_SOURCES_WRITE("PS1_V2",    CMF_PS1_V2);
     551            PM_SOURCES_WRITE("PS1_V3",    CMF_PS1_V3);
     552            PM_SOURCES_WRITE("PS1_SV1",   CMF_PS1_SV1);
     553            PM_SOURCES_WRITE("PS1_DV1",   CMF_PS1_DV1);
     554            PM_SOURCES_WRITE("PS1_DV2",   CMF_PS1_DV2);
     555
    588556            psFree (outhead);
    589557            psFree (exttype);
    … …  
    10391007                sources = pmSourcesRead_CMF_PS1_DV1 (file->fits, hdu->header);
    10401008            }
     1009            if (!strcmp (exttype, "PS1_DV2")) {
     1010                sources = pmSourcesRead_CMF_PS1_DV2 (file->fits, hdu->header);
     1011            }
    10411012
    10421013            if (!pmReadoutReadDetEff(file->fits, readout, deteffname)) {

  • trunk/psModules/src/objects/pmSourceIO.h

    r28013 r29004  
    2424bool pmSourcesWriteCMP (psArray *sources, char *filename, psMetadata *header);
    2525
    26 bool pmSourcesWrite_SMPDATA (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
    27 bool pmSourcesWrite_PS1_DEV_0 (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
     26bool pmSource_CMF_WritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config);
    2827
    29 bool pmSourcesWrite_PS1_DEV_1 (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
    30 bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe);
    31 bool pmSourcesWrite_PS1_DEV_1_XFIT (psFits *fits, psArray *sources, char *extname);
     28// All of these functions need to use the same API, even if not all elements are used in a specific case
     29bool pmSourcesWrite_SMPDATA(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
     30bool pmSourcesWrite_SMPDATA_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
     31bool pmSourcesWrite_SMPDATA_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
    3232
    33 bool pmSourcesWrite_PS1_CAL_0 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
    34 bool pmSourcesWrite_PS1_CAL_0_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    35 bool pmSourcesWrite_PS1_CAL_0_XFIT (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
     33bool pmSourcesWrite_PS1_DEV_0(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
     34bool pmSourcesWrite_PS1_DEV_0_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
     35bool pmSourcesWrite_PS1_DEV_0_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
    3636
    37 bool pmSourcesWrite_CMF_PS1_V1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
    38 bool pmSourcesWrite_CMF_PS1_V1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    39 bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname);
     37bool pmSourcesWrite_PS1_DEV_1(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
     38bool pmSourcesWrite_PS1_DEV_1_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
     39bool pmSourcesWrite_PS1_DEV_1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
    4040
    41 bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
    42 bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    43 bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname);
     41bool pmSourcesWrite_PS1_CAL_0(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
     42bool pmSourcesWrite_PS1_CAL_0_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
     43bool pmSourcesWrite_PS1_CAL_0_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
    4444
    45 bool pmSourcesWrite_CMF_PS1_SV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
    46 bool pmSourcesWrite_CMF_PS1_SV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    47 bool pmSourcesWrite_CMF_PS1_SV1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname);
     45bool pmSourcesWrite_CMF_PS1_V1(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
     46bool pmSourcesWrite_CMF_PS1_V1_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
     47bool pmSourcesWrite_CMF_PS1_V1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
    4848
    49 bool pmSourcesWrite_CMF_PS1_DV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
    50 bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
    51 bool pmSourcesWrite_CMF_PS1_DV1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname);
     49bool pmSourcesWrite_CMF_PS1_V2(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
     50bool pmSourcesWrite_CMF_PS1_V2_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
     51bool pmSourcesWrite_CMF_PS1_V2_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
    5252
    53 bool pmSource_CMF_WritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config);
     53bool pmSourcesWrite_CMF_PS1_V3(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
     54bool pmSourcesWrite_CMF_PS1_V3_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
     55bool pmSourcesWrite_CMF_PS1_V3_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
     56
     57bool pmSourcesWrite_CMF_PS1_SV1(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
     58bool pmSourcesWrite_CMF_PS1_SV1_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
     59bool pmSourcesWrite_CMF_PS1_SV1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
     60
     61bool pmSourcesWrite_CMF_PS1_DV1(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
     62bool pmSourcesWrite_CMF_PS1_DV1_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
     63bool pmSourcesWrite_CMF_PS1_DV1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
     64
     65bool pmSourcesWrite_CMF_PS1_DV2(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe);
     66bool pmSourcesWrite_CMF_PS1_DV2_XSRC(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
     67bool pmSourcesWrite_CMF_PS1_DV2_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
    5468
    5569psArray *pmSourcesReadCMP (char *filename, psMetadata *header);
    … …  
    6377psArray *pmSourcesRead_CMF_PS1_SV1 (psFits *fits, psMetadata *header);
    6478psArray *pmSourcesRead_CMF_PS1_DV1 (psFits *fits, psMetadata *header);
     79psArray *pmSourcesRead_CMF_PS1_DV2 (psFits *fits, psMetadata *header);
    6580
    6681bool pmSourcesWritePSFs (psArray *sources, char *filename);

  • trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c

    r28013 r29004  
    2828#include "pmFPAfile.h"
    2929
     30#include "pmTrend2D.h"
     31#include "pmResiduals.h"
     32#include "pmGrowthCurve.h"
    3033#include "pmSpan.h"
     34#include "pmFootprintSpans.h"
    3135#include "pmFootprint.h"
    3236#include "pmPeaks.h"
    3337#include "pmMoments.h"
    34 #include "pmGrowthCurve.h"
    35 #include "pmResiduals.h"
    36 #include "pmTrend2D.h"
     38#include "pmModelFuncs.h"
     39#include "pmModel.h"
     40#include "pmModelUtils.h"
     41#include "pmModelClass.h"
     42#include "pmSourceMasks.h"
     43#include "pmSourceExtendedPars.h"
     44#include "pmSourceDiffStats.h"
     45#include "pmSource.h"
     46#include "pmSourceFitModel.h"
    3747#include "pmPSF.h"
    38 #include "pmModel.h"
    39 #include "pmSource.h"
    40 #include "pmModelClass.h"
     48#include "pmPSFtry.h"
     49
    4150#include "pmSourceIO.h"
    4251
    … …  
    4756// This version has elements intended for difference images & forced photometry:
    4857// diffStats entries (good for dipoles); flux + flux error (for insignificant detections)
    49 bool pmSourcesWrite_CMF_PS1_DV1 (psFits *fits, pmReadout *readout, psArray *sources,
    50                                 psMetadata *imageHeader, psMetadata *tableHeader, char *extname)
     58bool pmSourcesWrite_CMF_PS1_DV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe)
    5159{
    5260    PS_ASSERT_PTR_NON_NULL(fits, false);
    … …  
    194202        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
    195203        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
    196         psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeight);
     204        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeightNotBad);
    197205        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
    198206        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
    … …  
    324332        source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
    325333        source->peak->flux = peakFlux;
     334        source->peak->xf   = PAR[PM_PAR_XPOS]; // more accurate position
     335        source->peak->yf   = PAR[PM_PAR_YPOS]; // more accurate position
    326336        source->peak->dx   = dPAR[PM_PAR_XPOS];
    327337        source->peak->dy   = dPAR[PM_PAR_YPOS];
    328         source->peak->SN   = sqrt(source->peak->flux); // XXX a proxy: various functions sort by peak S/N
    329 
    330         source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF");
     338        if (isfinite (source->errMag) && (source->errMag > 0.0)) {
     339          source->peak->SN = 1.0 / source->errMag;
     340        } else {
     341          source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N
     342        }
     343
     344        source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
    331345        source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
    332346        source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA");
    … …  
    536550
    537551// XXX this layout is still the same as PS1_DEV_1
    538 bool pmSourcesWrite_CMF_PS1_DV1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname)
     552bool pmSourcesWrite_CMF_PS1_DV1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname)
    539553{
    540554
    … …  
    585599            assert (model);
    586600
     601            // skip models which were not actually fitted
     602            if (model->flags & PM_MODEL_STATUS_BADARGS) continue;
     603
    587604            PAR = model->params->data.F32;
    588605            dPAR = model->dparams->data.F32;

  • trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c

    r28013 r29004  
    2828#include "pmFPAfile.h"
    2929
     30#include "pmTrend2D.h"
     31#include "pmResiduals.h"
     32#include "pmGrowthCurve.h"
    3033#include "pmSpan.h"
     34#include "pmFootprintSpans.h"
    3135#include "pmFootprint.h"
    3236#include "pmPeaks.h"
    3337#include "pmMoments.h"
    34 #include "pmGrowthCurve.h"
    35 #include "pmResiduals.h"
    36 #include "pmTrend2D.h"
     38#include "pmModelFuncs.h"
     39#include "pmModel.h"
     40#include "pmModelUtils.h"
     41#include "pmModelClass.h"
     42#include "pmSourceMasks.h"
     43#include "pmSourceExtendedPars.h"
     44#include "pmSourceDiffStats.h"
     45#include "pmSource.h"
     46#include "pmSourceFitModel.h"
    3747#include "pmPSF.h"
    38 #include "pmModel.h"
    39 #include "pmSource.h"
    40 #include "pmModelClass.h"
     48#include "pmPSFtry.h"
     49
    4150#include "pmSourceIO.h"
    4251
    … …  
    4655
    4756// NOTE: this output function is intended for psphotStack analysis: it includes per-psf radial fluxes
    48 // XXX currently in the 'read' function is NOT consistent with the 'write' function (does not read radial fluxes)
     57// XXX currently, the 'read' function is NOT consistent with the 'write' function (does not read radial fluxes)
    4958
    5059bool pmSourcesWrite_CMF_PS1_SV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe)
    … …  
    6271    psF32 errMag, chisq, apRadius;
    6372    psS32 nPix, nDOF;
     73    char keyword1[80], keyword2[80];
    6474
    6575    pmChip *chip = readout->parent->parent;
    … …  
    95105    // we use this just to define the output vectors (which must be present for all objects)
    96106    bool status = false;
     107    psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
    97108    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
    98109    psAssert (radMax, "this must have been defined and tested earlier!");
    99110    psAssert (radMax->n, "this must have been defined and tested earlier!");
     111    psAssert (radMin->n == radMax->n, "inconsistent annular bins");
     112
     113    // write the radial profile apertures to header
     114    for (int i = 0; i < radMax->n; i++) {
     115      sprintf (keyword1, "RMIN_%02d", i);
     116      sprintf (keyword2, "RMAX_%02d", i);
     117      psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]);
     118      psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]);
     119    }
     120
    100121
    101122    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
    … …  
    193214        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
    194215        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
    195         psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeight);
     216        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor (bad)",          source->pixWeightNotBad);
     217        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
    196218        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
    197219        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
    198220
    199221        // distinguish moments measure from window vs S/N > XX ??
    200         float mxx = source->moments ? source->moments->Mxx : NAN;
    201         float mxy = source->moments ? source->moments->Mxy : NAN;
    202         float myy = source->moments ? source->moments->Myy : NAN;
    203         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      mxx);
    204         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      mxy);
    205         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      myy);
    206 
    207         psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
     222        float Mxx = source->moments ? source->moments->Mxx : NAN;
     223        float Mxy = source->moments ? source->moments->Mxy : NAN;
     224        float Myy = source->moments ? source->moments->Myy : NAN;
     225
     226        float Mrf  = source->moments ? source->moments->Mrf : NAN;
     227        float Mrh  = source->moments ? source->moments->Mrh : NAN;
     228        float Krf  = source->moments ? source->moments->KronFlux : NAN;
     229        float dKrf = source->moments ? source->moments->KronFluxErr : NAN;
     230
     231        float Kinner = source->moments ? source->moments->KronFinner : NAN;
     232        float Kouter = source->moments ? source->moments->KronFouter : NAN;
     233
     234        float M_c3 = source->moments ? 1.0*source->moments->Mxxx - 3.0*source->moments->Mxyy : NAN;
     235        float M_s3 = source->moments ? 3.0*source->moments->Mxxy - 1.0*source->moments->Myyy : NAN;
     236        float M_c4 = source->moments ? 1.0*source->moments->Mxxxx - 6.0*source->moments->Mxxyy + 1.0*source->moments->Myyyy : NAN;
     237        float M_s4 = source->moments ? 4.0*source->moments->Mxxxy - 4.0*source->moments->Mxyyy : NAN;
     238
     239        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      Mxx);
     240        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      Mxy);
     241        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      Myy);
     242
     243        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                    M_c3);
     244        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                    M_s3);
     245        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                   M_c4);
     246        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                   M_s4);
     247
     248        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                       Mrf);
     249        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                        Mrh);
     250        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Krf);
     251        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                          dKrf);
     252
     253        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Kinner);
     254        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Kouter);
     255
     256        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                     source->mode);
     257        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                     source->mode2);
    208258
    209259        psVector *radFlux    = psVectorAlloc(radMax->n, PS_TYPE_F32);
    … …  
    352402        source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
    353403        source->peak->flux = peakFlux;
     404        source->peak->xf   = PAR[PM_PAR_XPOS]; // more accurate position
     405        source->peak->yf   = PAR[PM_PAR_YPOS]; // more accurate position
    354406        source->peak->dx   = dPAR[PM_PAR_XPOS];
    355407        source->peak->dy   = dPAR[PM_PAR_YPOS];
    356         source->peak->SN   = sqrt(source->peak->flux); // XXX a proxy: various functions sort by peak S/N
    357 
    358         source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF");
     408        if (isfinite (source->errMag) && (source->errMag > 0.0)) {
     409          source->peak->SN = 1.0 / source->errMag;
     410        } else {
     411          source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N
     412        }
     413
     414        source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
     415        source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT");
    359416        source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
    360417        source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA");
    … …  
    371428        source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY");
    372429
     430        source->moments->Mrf         = psMetadataLookupF32 (&status, row, "MOMENTS_R1");
     431        source->moments->Mrh         = psMetadataLookupF32 (&status, row, "MOMENTS_RH");
     432        source->moments->KronFlux    = psMetadataLookupF32 (&status, row, "KRON_FLUX");
     433        source->moments->KronFluxErr = psMetadataLookupF32 (&status, row, "KRON_FLUX_ERR");
     434
     435        source->moments->KronFinner  = psMetadataLookupF32 (&status, row, "KRON_FLUX_INNER");
     436        source->moments->KronFouter  = psMetadataLookupF32 (&status, row, "KRON_FLUX_OUTER");
     437
     438        // XXX we do not save all of the 3rd and 4th moment parameters. when we load in data,
     439        // we are storing enough information so the output will be consistent with the input
     440        source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C");
     441        source->moments->Mxxy = 0.0;
     442        source->moments->Mxyy = 0.0;
     443        source->moments->Myyy = -1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3S");
     444
     445        source->moments->Mxxxx = +1.00 * psMetadataLookupF32 (&status, row, "MOMENTS_M4C");
     446        source->moments->Mxxxy = 0.0;
     447        source->moments->Mxxyy = 0.0;
     448        source->moments->Mxyyy = -0.25 * psMetadataLookupF32 (&status, row, "MOMENTS_M4S");
     449        source->moments->Myyyy = 0.0;
     450
    373451        source->mode = psMetadataLookupU32 (&status, row, "FLAGS");
     452        source->mode2 = psMetadataLookupU32 (&status, row, "FLAGS2");
    374453        assert (status);
    375454
    … …  
    391470    psF32 xErr, yErr;
    392471    int nRow = -1;
     472    char keyword1[80], keyword2[80];
    393473
    394474    // create a header to hold the output data
    … …  
    422502    bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
    423503    bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
    424     // bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
    425     // bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
    426504
    427505    psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
    … …  
    429507    psAssert (radMin->n == radMax->n, "inconsistent annular bins");
    430508
    431     // int nRadialBins = 0;
    432     // if (doAnnuli) {
    433     //  // get the max count of radial bins
    434     //  for (int i = 0; i < sources->n; i++) {
    435     //      pmSource *source = sources->data[i];
    436     //      if (!source->extpars) continue;
    437     //      if (!source->extpars->radProfile ) continue;
    438     //         if (!source->extpars->radProfile->binSB) continue;
    439     //      nRadialBins = PS_MAX(nRadialBins, source->extpars->radProfile->binSB->n);
    440     //  }
    441     // }
     509    // write the radial profile apertures to header
     510    for (int i = 0; i < radMax->n; i++) {
     511      sprintf (keyword1, "RMIN_%02d", i);
     512      sprintf (keyword2, "RMAX_%02d", i);
     513      psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]);
     514      psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]);
     515    }
    442516
    443517    // we write out all sources, regardless of quality.  the source flags tell us the state
    444518    for (int i = 0; i < sources->n; i++) {
    445         // skip source if it is not a ext sourc
    446         // XXX we have two places that extended source parameters are measured:
    447         // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models,
    448         // psphotFitEXT, which does the simple extended source model fit (not psf-convolved)
    449         // should we require both?
    450 
    451519        pmSource *source = sources->data[i];
    452520
    … …  
    514582        }
    515583
    516 # if (0)
    517         // Kron measurements
    518         if (doKron) {
    519             pmSourceKronValues *kron = source->extpars->kron;
    520             if (kron) {
    521                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       kron->mag);
    522                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", kron->magErr);
    523                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          kron->rad);
    524                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    kron->radErr);
    525             } else {
    526                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       NAN);
    527                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", NAN);
    528                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          NAN);
    529                 psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    NAN);
    530             }
    531         }
    532 
    533         // Isophot measurements
    534         // XXX insert header data: isophotal level
    535         if (doIsophotal) {
    536             pmSourceIsophotalValues *isophot = source->extpars->isophot;
    537             if (isophot) {
    538                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       isophot->mag);
    539                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr);
    540                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          isophot->rad);
    541                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    isophot->radErr);
    542             } else {
    543                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       NAN);
    544                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", NAN);
    545                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          NAN);
    546                 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    NAN);
    547             }
    548         }
    549 # endif
    550 
    551584        // Flux Annuli (if we have extended source measurements, we have these.  only optionally save them)
    552585        if (doAnnuli) {
    … …  
    616649
    617650// XXX this layout is still the same as PS1_DEV_1
    618 bool pmSourcesWrite_CMF_PS1_SV1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname)
     651bool pmSourcesWrite_CMF_PS1_SV1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname)
    619652{
    620653
    … …  
    665698            assert (model);
    666699
     700            // skip models which were not actually fitted
     701            if (model->flags & PM_MODEL_STATUS_BADARGS) continue;
     702
    667703            PAR = model->params->data.F32;
    668704            dPAR = model->dparams->data.F32;

  • trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c

    r27818 r29004  
    2828#include "pmFPAfile.h"
    2929
     30#include "pmTrend2D.h"
     31#include "pmResiduals.h"
     32#include "pmGrowthCurve.h"
    3033#include "pmSpan.h"
     34#include "pmFootprintSpans.h"
    3135#include "pmFootprint.h"
    3236#include "pmPeaks.h"
    3337#include "pmMoments.h"
    34 #include "pmGrowthCurve.h"
    35 #include "pmResiduals.h"
    36 #include "pmTrend2D.h"
     38#include "pmModelFuncs.h"
     39#include "pmModel.h"
     40#include "pmModelUtils.h"
     41#include "pmModelClass.h"
     42#include "pmSourceMasks.h"
     43#include "pmSourceExtendedPars.h"
     44#include "pmSourceDiffStats.h"
     45#include "pmSource.h"
     46#include "pmSourceFitModel.h"
    3747#include "pmPSF.h"
    38 #include "pmModel.h"
    39 #include "pmSource.h"
    40 #include "pmModelClass.h"
     48#include "pmPSFtry.h"
     49
    4150#include "pmSourceIO.h"
    4251
    … …  
    4554// followed by a zero-size matrix, followed by the table data
    4655
    47 bool pmSourcesWrite_CMF_PS1_V1 (psFits *fits, pmReadout *readout, psArray *sources,
    48                                 psMetadata *imageHeader, psMetadata *tableHeader, char *extname)
     56bool pmSourcesWrite_CMF_PS1_V1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe)
    4957{
    5058    PS_ASSERT_PTR_NON_NULL(fits, false);
    … …  
    184192        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
    185193        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
    186         psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeight);
     194        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeightNotBad);
    187195        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
    188196        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
    … …  
    310318        source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
    311319        source->peak->flux = peakFlux;
     320        source->peak->xf   = PAR[PM_PAR_XPOS]; // more accurate position
     321        source->peak->yf   = PAR[PM_PAR_YPOS]; // more accurate position
    312322        source->peak->dx   = dPAR[PM_PAR_XPOS];
    313323        source->peak->dy   = dPAR[PM_PAR_YPOS];
    314 
    315         source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF");
     324        if (isfinite (source->errMag) && (source->errMag > 0.0)) {
     325          source->peak->SN = 1.0 / source->errMag;
     326        } else {
     327          source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N
     328        }
     329
     330        source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
    316331        source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
    317332        source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA");
    … …  
    516531
    517532// XXX this layout is still the same as PS1_DEV_1
    518 bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname)
     533bool pmSourcesWrite_CMF_PS1_V1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname)
    519534{
    520535
    … …  
    565580            assert (model);
    566581
     582            // skip models which were not actually fitted
     583            if (model->flags & PM_MODEL_STATUS_BADARGS) continue;
     584
    567585            PAR = model->params->data.F32;
    568586            dPAR = model->dparams->data.F32;

  • trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c

    r27818 r29004  
    2828#include "pmFPAfile.h"
    2929
     30#include "pmTrend2D.h"
     31#include "pmResiduals.h"
     32#include "pmGrowthCurve.h"
    3033#include "pmSpan.h"
     34#include "pmFootprintSpans.h"
    3135#include "pmFootprint.h"
    3236#include "pmPeaks.h"
    3337#include "pmMoments.h"
    34 #include "pmGrowthCurve.h"
    35 #include "pmResiduals.h"
    36 #include "pmTrend2D.h"
     38#include "pmModelFuncs.h"
     39#include "pmModel.h"
     40#include "pmModelUtils.h"
     41#include "pmModelClass.h"
     42#include "pmSourceMasks.h"
     43#include "pmSourceExtendedPars.h"
     44#include "pmSourceDiffStats.h"
     45#include "pmSource.h"
     46#include "pmSourceFitModel.h"
    3747#include "pmPSF.h"
    38 #include "pmModel.h"
    39 #include "pmSource.h"
    40 #include "pmModelClass.h"
     48#include "pmPSFtry.h"
     49
    4150#include "pmSourceIO.h"
    4251
    … …  
    4554// followed by a zero-size matrix, followed by the table data
    4655
    47 bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname)
     56bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe)
    4857{
    4958    PS_ASSERT_PTR_NON_NULL(fits, false);
    … …  
    182191        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
    183192        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
    184         psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeight);
     193        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeightNotBad);
    185194        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
    186195        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
    … …  
    314323        source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
    315324        source->peak->flux = peakFlux;
     325        source->peak->xf   = PAR[PM_PAR_XPOS]; // more accurate position
     326        source->peak->yf   = PAR[PM_PAR_YPOS]; // more accurate position
    316327        source->peak->dx   = dPAR[PM_PAR_XPOS];
    317328        source->peak->dy   = dPAR[PM_PAR_YPOS];
    318         source->peak->SN   = sqrt(source->peak->flux); // XXX a proxy: various functions sort by peak S/N
    319 
    320         source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF");
     329        if (isfinite (source->errMag) && (source->errMag > 0.0)) {
     330          source->peak->SN = 1.0 / source->errMag;
     331        } else {
     332          source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N
     333        }
     334
     335        source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
    321336        source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
    322337        source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA");
    … …  
    578593
    579594// XXX this layout is still the same as PS1_DEV_1
    580 bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname)
     595bool pmSourcesWrite_CMF_PS1_V2_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname)
    581596{
    582597
    … …  
    627642            assert (model);
    628643
     644            // skip models which were not actually fitted
     645            if (model->flags & PM_MODEL_STATUS_BADARGS) continue;
     646
    629647            PAR = model->params->data.F32;
    630648            dPAR = model->dparams->data.F32;

  • trunk/psModules/src/objects/pmSourceIO_CMP.c

    r20937 r29004  
    2828#include "pmFPAfile.h"
    2929
     30#include "pmTrend2D.h"
     31#include "pmResiduals.h"
     32#include "pmGrowthCurve.h"
    3033#include "pmSpan.h"
     34#include "pmFootprintSpans.h"
    3135#include "pmFootprint.h"
    3236#include "pmPeaks.h"
    3337#include "pmMoments.h"
    34 #include "pmGrowthCurve.h"
    35 #include "pmResiduals.h"
    36 #include "pmTrend2D.h"
     38#include "pmModelFuncs.h"
     39#include "pmModel.h"
     40#include "pmModelUtils.h"
     41#include "pmModelClass.h"
     42#include "pmSourceMasks.h"
     43#include "pmSourceExtendedPars.h"
     44#include "pmSourceDiffStats.h"
     45#include "pmSource.h"
     46#include "pmSourceFitModel.h"
    3747#include "pmPSF.h"
    38 #include "pmModel.h"
    39 #include "pmSource.h"
    40 #include "pmModelClass.h"
     48#include "pmPSFtry.h"
     49
    4150#include "pmSourceIO.h"
    4251

  • trunk/psModules/src/objects/pmSourceIO_MatchedRefs.c

    r27177 r29004  
    2828#include "pmFPAfile.h"
    2929
     30#include "pmTrend2D.h"
     31#include "pmResiduals.h"
     32#include "pmGrowthCurve.h"
    3033#include "pmSpan.h"
     34#include "pmFootprintSpans.h"
    3135#include "pmFootprint.h"
    3236#include "pmPeaks.h"
    3337#include "pmMoments.h"
    34 #include "pmGrowthCurve.h"
    35 #include "pmResiduals.h"
    36 #include "pmTrend2D.h"
     38#include "pmModelFuncs.h"
     39#include "pmModel.h"
     40#include "pmModelUtils.h"
     41#include "pmModelClass.h"
     42#include "pmSourceMasks.h"
     43#include "pmSourceExtendedPars.h"
     44#include "pmSourceDiffStats.h"
     45#include "pmSource.h"
     46#include "pmSourceFitModel.h"
    3747#include "pmPSF.h"
    38 #include "pmModel.h"
    39 #include "pmSource.h"
    40 #include "pmModelClass.h"
     48#include "pmPSFtry.h"
     49
    4150#include "pmSourceIO.h"
     51
    4252#include "pmAstrometryObjects.h"
    4353#include "pmAstrometryWCS.h"

  • trunk/psModules/src/objects/pmSourceIO_OBJ.c

    r20937 r29004  
    2828#include "pmFPAfile.h"
    2929
     30#include "pmTrend2D.h"
     31#include "pmResiduals.h"
     32#include "pmGrowthCurve.h"
    3033#include "pmSpan.h"
     34#include "pmFootprintSpans.h"
    3135#include "pmFootprint.h"
    3236#include "pmPeaks.h"
    3337#include "pmMoments.h"
    34 #include "pmGrowthCurve.h"
    35 #include "pmResiduals.h"
    36 #include "pmTrend2D.h"
     38#include "pmModelFuncs.h"
     39#include "pmModel.h"
     40#include "pmModelUtils.h"
     41#include "pmModelClass.h"
     42#include "pmSourceMasks.h"
     43#include "pmSourceExtendedPars.h"
     44#include "pmSourceDiffStats.h"
     45#include "pmSource.h"
     46#include "pmSourceFitModel.h"
    3747#include "pmPSF.h"
    38 #include "pmModel.h"
    39 #include "pmSource.h"
    40 #include "pmModelClass.h"
     48#include "pmPSFtry.h"
     49
    4150#include "pmSourceIO.h"
    4251

  • trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c

    r25979 r29004  
    2828#include "pmFPAfile.h"
    2929
     30#include "pmTrend2D.h"
     31#include "pmResiduals.h"
     32#include "pmGrowthCurve.h"
    3033#include "pmSpan.h"
     34#include "pmFootprintSpans.h"
    3135#include "pmFootprint.h"
    3236#include "pmPeaks.h"
    3337#include "pmMoments.h"
    34 #include "pmGrowthCurve.h"
    35 #include "pmResiduals.h"
    36 #include "pmTrend2D.h"
     38#include "pmModelFuncs.h"
     39#include "pmModel.h"
     40#include "pmModelUtils.h"
     41#include "pmModelClass.h"
     42#include "pmSourceMasks.h"
     43#include "pmSourceExtendedPars.h"
     44#include "pmSourceDiffStats.h"
     45#include "pmSource.h"
     46#include "pmSourceFitModel.h"
    3747#include "pmPSF.h"
    38 #include "pmModel.h"
    39 #include "pmSource.h"
    40 #include "pmModelClass.h"
     48#include "pmPSFtry.h"
     49
    4150#include "pmSourceIO.h"
    4251
    … …  
    4958// XXX how do I generate the source tables which I need to send to PSPS?
    5059
    51 bool pmSourcesWrite_PS1_CAL_0 (psFits *fits, pmReadout *readout, psArray *sources,
    52                                psMetadata *imageHeader, psMetadata *tableHeader, char *extname)
     60bool pmSourcesWrite_PS1_CAL_0 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe)
    5361{
    5462    PS_ASSERT_PTR_NON_NULL(fits, false);
    … …  
    180188        psMetadataAdd (row, PS_LIST_TAIL, "PSF_WIDTH_Y",      PS_DATA_F32, "PSF width in y coordinate",                  axes.minor);
    181189        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
    182         psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeight);
     190        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeightNotBad);
    183191
    184192        // XXX not sure how to get this : need to load Nimages with weight?
    … …  
    289297        source->peak->dx   = dPAR[PM_PAR_XPOS];
    290298        source->peak->dy   = dPAR[PM_PAR_YPOS];
    291 
    292         source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF");
     299        source->peak->xf   = PAR[PM_PAR_XPOS]; // more accurate position
     300        source->peak->yf   = PAR[PM_PAR_YPOS]; // more accurate position
     301
     302        source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
    293303
    294304        // note that some older versions used PSF_PROBABILITY: this was not well defined.
    … …  
    597607            assert (model);
    598608
     609            // skip models which were not actually fitted
     610            if (model->flags & PM_MODEL_STATUS_BADARGS) continue;
     611
    599612            PAR = model->params->data.F32;
    600613            dPAR = model->dparams->data.F32;

  • trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c

    r25979 r29004  
    2828#include "pmFPAfile.h"
    2929
     30#include "pmTrend2D.h"
     31#include "pmResiduals.h"
     32#include "pmGrowthCurve.h"
    3033#include "pmSpan.h"
     34#include "pmFootprintSpans.h"
    3135#include "pmFootprint.h"
    3236#include "pmPeaks.h"
    3337#include "pmMoments.h"
    34 #include "pmGrowthCurve.h"
    35 #include "pmResiduals.h"
    36 #include "pmTrend2D.h"
     38#include "pmModelFuncs.h"
     39#include "pmModel.h"
     40#include "pmModelUtils.h"
     41#include "pmModelClass.h"
     42#include "pmSourceMasks.h"
     43#include "pmSourceExtendedPars.h"
     44#include "pmSourceDiffStats.h"
     45#include "pmSource.h"
     46#include "pmSourceFitModel.h"
    3747#include "pmPSF.h"
    38 #include "pmModel.h"
    39 #include "pmSource.h"
    40 #include "pmModelClass.h"
     48#include "pmPSFtry.h"
     49
    4150#include "pmSourceIO.h"
    4251
    … …  
    4958// XXX how do I generate the source tables which I need to send to PSPS?
    5059// XXX: input parameter imageHeader is never used.
    51 bool pmSourcesWrite_PS1_DEV_0 (psFits *fits, psArray *sources, psMetadata *imageHeader,
    52                                psMetadata *tableHeader, char *extname)
     60bool pmSourcesWrite_PS1_DEV_0 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe)
    5361{
    5462    PS_ASSERT_PTR_NON_NULL(fits, false);
    … …  
    117125        psMetadataAdd (row, PS_LIST_TAIL, "PSF_WIDTH_Y",      PS_DATA_F32, "PSF width in y coordinate",                  axes.minor);
    118126        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
    119         psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeight);
     127        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeightNotBad);
    120128        // XXX not sure how to get this : need to load Nimages with weight
    121129        psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", nImageOverlap);
    … …  
    215223        source->peak->dx   = dPAR[PM_PAR_XPOS];
    216224        source->peak->dy   = dPAR[PM_PAR_YPOS];
    217 
    218         source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF");
     225        source->peak->xf   = PAR[PM_PAR_XPOS]; // more accurate position
     226        source->peak->yf   = PAR[PM_PAR_YPOS]; // more accurate position
     227
     228        source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
    219229
    220230        // XXX other values saved but not loaded?
    … …  
    228238    return (sources);
    229239}
     240
     241bool pmSourcesWrite_PS1_DEV_0_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
     242{
     243    return true;
     244}
     245
     246bool pmSourcesWrite_PS1_DEV_0_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname)
     247{
     248    return true;
     249}

  • trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c

    r25979 r29004  
    2828#include "pmFPAfile.h"
    2929
     30#include "pmTrend2D.h"
     31#include "pmResiduals.h"
     32#include "pmGrowthCurve.h"
    3033#include "pmSpan.h"
     34#include "pmFootprintSpans.h"
    3135#include "pmFootprint.h"
    3236#include "pmPeaks.h"
    3337#include "pmMoments.h"
    34 #include "pmGrowthCurve.h"
    35 #include "pmResiduals.h"
    36 #include "pmTrend2D.h"
     38#include "pmModelFuncs.h"
     39#include "pmModel.h"
     40#include "pmModelUtils.h"
     41#include "pmModelClass.h"
     42#include "pmSourceMasks.h"
     43#include "pmSourceExtendedPars.h"
     44#include "pmSourceDiffStats.h"
     45#include "pmSource.h"
     46#include "pmSourceFitModel.h"
    3747#include "pmPSF.h"
    38 #include "pmModel.h"
    39 #include "pmSource.h"
    40 #include "pmModelClass.h"
     48#include "pmPSFtry.h"
     49
    4150#include "pmSourceIO.h"
    4251
    … …  
    4756// this output format is valid for psphot analysis of an image, and does not include calibrated
    4857// values derived in the DVO database.
    49 // XXX how do I generate the source tables which I need to send to PSPS?
    50 
    51 bool pmSourcesWrite_PS1_DEV_1 (psFits *fits, psArray *sources,
    52                                psMetadata *imageHeader, psMetadata *tableHeader, char *extname)
     58
     59bool pmSourcesWrite_PS1_DEV_1(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe)
    5360{
    5461    PS_ASSERT_PTR_NON_NULL(fits, false);
    … …  
    142149        psMetadataAdd (row, PS_LIST_TAIL, "PSF_WIDTH_Y",      PS_DATA_F32, "PSF width in y coordinate",                  axes.minor);
    143150        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
    144         psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeight);
     151        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeightNotBad);
    145152
    146153        // XXX not sure how to get this : need to load Nimages with weight?
    … …  
    259266        source->peak->dx   = dPAR[PM_PAR_XPOS];
    260267        source->peak->dy   = dPAR[PM_PAR_YPOS];
    261 
    262         source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF");
     268        source->peak->xf   = PAR[PM_PAR_XPOS]; // more accurate position
     269        source->peak->yf   = PAR[PM_PAR_YPOS]; // more accurate position
     270
     271        source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
    263272
    264273        // note that some older versions used PSF_PROBABILITY: this was not well defined.
    … …  
    281290}
    282291
    283 bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe)
     292bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
    284293{
    285294
    … …  
    453462}
    454463
    455 bool pmSourcesWrite_PS1_DEV_1_XFIT (psFits *fits, psArray *sources, char *extname)
     464bool pmSourcesWrite_PS1_DEV_1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname)
    456465{
    457466
    … …  
    502511            assert (model);
    503512
     513            // skip models which were not actually fitted
     514            if (model->flags & PM_MODEL_STATUS_BADARGS) continue;
     515
    504516            PAR = model->params->data.F32;
    505517            dPAR = model->dparams->data.F32;

  • trunk/psModules/src/objects/pmSourceIO_RAW.c

    r25754 r29004  
    2828#include "pmFPAfile.h"
    2929
     30#include "pmTrend2D.h"
     31#include "pmResiduals.h"
     32#include "pmGrowthCurve.h"
    3033#include "pmSpan.h"
     34#include "pmFootprintSpans.h"
    3135#include "pmFootprint.h"
    3236#include "pmPeaks.h"
    3337#include "pmMoments.h"
    34 #include "pmGrowthCurve.h"
    35 #include "pmResiduals.h"
    36 #include "pmTrend2D.h"
     38#include "pmModelFuncs.h"
     39#include "pmModel.h"
     40#include "pmModelUtils.h"
     41#include "pmModelClass.h"
     42#include "pmSourceMasks.h"
     43#include "pmSourceExtendedPars.h"
     44#include "pmSourceDiffStats.h"
     45#include "pmSource.h"
     46#include "pmSourceFitModel.h"
    3747#include "pmPSF.h"
    38 #include "pmModel.h"
    39 #include "pmSource.h"
    40 #include "pmModelClass.h"
    41 #include "pmSourcePhotometry.h"
     48#include "pmPSFtry.h"
     49
    4250#include "pmSourceIO.h"
    4351
    … …  
    120128                 source[0].peak->SN,
    121129                 source[0].apRadius,
    122                  source[0].pixWeight,
     130                 source[0].pixWeightNotBad,
    123131                 model[0].nDOF,
    124132                 model[0].nIter);
    … …  
    182190                 source[0].peak->SN,
    183191                 source[0].apRadius,
    184                  source[0].pixWeight,
     192                 source[0].pixWeightNotBad,
    185193                 model[0].nDOF,
    186194                 model[0].nIter);

  • trunk/psModules/src/objects/pmSourceIO_SMPDATA.c

    r25979 r29004  
    2828#include "pmFPAfile.h"
    2929
     30#include "pmTrend2D.h"
     31#include "pmResiduals.h"
     32#include "pmGrowthCurve.h"
    3033#include "pmSpan.h"
     34#include "pmFootprintSpans.h"
    3135#include "pmFootprint.h"
    3236#include "pmPeaks.h"
    3337#include "pmMoments.h"
    34 #include "pmResiduals.h"
    35 #include "pmGrowthCurve.h"
    36 #include "pmTrend2D.h"
     38#include "pmModelFuncs.h"
     39#include "pmModel.h"
     40#include "pmModelUtils.h"
     41#include "pmModelClass.h"
     42#include "pmSourceMasks.h"
     43#include "pmSourceExtendedPars.h"
     44#include "pmSourceDiffStats.h"
     45#include "pmSource.h"
     46#include "pmSourceFitModel.h"
    3747#include "pmPSF.h"
    38 #include "pmModel.h"
    39 #include "pmSource.h"
    40 #include "pmModelClass.h"
     48#include "pmPSFtry.h"
     49
    4150#include "pmSourceIO.h"
    4251
    … …  
    4554// followed by a zero-size matrix, followed by the table data
    4655// XXX: input parameter imageHeader is never used
    47 bool pmSourcesWrite_SMPDATA (psFits *fits, psArray *sources, psMetadata *imageHeader,
    48                              psMetadata *tableHeader, char *extname)
     56bool pmSourcesWrite_SMPDATA (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe)
    4957{
    5058    PS_ASSERT_PTR_NON_NULL(fits, false);
    … …  
    106114        psMetadataAdd (row, PS_LIST_TAIL, "THETA",   PS_DATA_F32, "", axes.theta);
    107115        psMetadataAdd (row, PS_LIST_TAIL, "DOPHOT",  PS_DATA_U8,  "", type);
    108         psMetadataAdd (row, PS_LIST_TAIL, "WEIGHT",  PS_DATA_U8,  "", PS_MIN (255, PS_MAX(0, 255*source->pixWeight)));
     116        psMetadataAdd (row, PS_LIST_TAIL, "WEIGHT",  PS_DATA_U8,  "", PS_MIN (255, PS_MAX(0, 255*source->pixWeightNotBad)));
    109117        psMetadataAdd (row, PS_LIST_TAIL, "DUMMY",   PS_DATA_U16, "", 0);
    110118
    … …  
    189197        source->apMag  = psMetadataLookupF32 (&status, row, "MAG_AP")  - ZERO_POINT;
    190198
    191         source->pixWeight = psMetadataLookupU8 (&status, row, "WEIGHT")/255.0;
     199        source->pixWeightNotBad = psMetadataLookupU8 (&status, row, "WEIGHT")/255.0;
    192200        int dophot = psMetadataLookupU8 (&status, row, "DOPHOT");
    193201        pmSourceSetDophotType (source, dophot);
    … …  
    204212    return (sources);
    205213}
     214
     215bool pmSourcesWrite_SMPDATA_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
     216{
     217    return true;
     218}
     219
     220bool pmSourcesWrite_SMPDATA_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname)
     221{
     222    return true;
     223}

  • trunk/psModules/src/objects/pmSourceIO_SX.c

    r20937 r29004  
    2828#include "pmFPAfile.h"
    2929
     30#include "pmTrend2D.h"
     31#include "pmResiduals.h"
     32#include "pmGrowthCurve.h"
    3033#include "pmSpan.h"
     34#include "pmFootprintSpans.h"
    3135#include "pmFootprint.h"
    3236#include "pmPeaks.h"
    3337#include "pmMoments.h"
    34 #include "pmGrowthCurve.h"
    35 #include "pmResiduals.h"
    36 #include "pmTrend2D.h"
     38#include "pmModelFuncs.h"
     39#include "pmModel.h"
     40#include "pmModelUtils.h"
     41#include "pmModelClass.h"
     42#include "pmSourceMasks.h"
     43#include "pmSourceExtendedPars.h"
     44#include "pmSourceDiffStats.h"
     45#include "pmSource.h"
     46#include "pmSourceFitModel.h"
    3747#include "pmPSF.h"
    38 #include "pmModel.h"
    39 #include "pmSource.h"
    40 #include "pmModelClass.h"
     48#include "pmPSFtry.h"
     49
    4150#include "pmSourceIO.h"
    4251

  • trunk/psModules/src/objects/pmSourceMasks.h

    r24694 r29004  
    3939} pmSourceMode;
    4040
     41// Bit flags to distinguish analysis results
     42// When adding to or subtracting from this list, please also modify pmSourceMaskHeader
     43typedef enum {
     44    PM_SOURCE_MODE2_DEFAULT          = 0x00000000, ///< Initial value: resets all bits
     45    PM_SOURCE_MODE2_DIFF_WITH_SINGLE = 0x00000001, ///< diff source matched to a single positive detection
     46    PM_SOURCE_MODE2_DIFF_WITH_DOUBLE = 0x00000002, ///< diff source matched to positive detections in both images
     47} pmSourceMode2;
    4148
    4249/// Populate header with mask values

  • trunk/psModules/src/objects/pmSourceMatch.c

    r28518 r29004  
    66#include <string.h>
    77#include <pslib.h>
    8 
     8#include "pmHDU.h"
     9#include "pmFPA.h"
     10
     11#include "pmTrend2D.h"
     12#include "pmResiduals.h"
     13#include "pmGrowthCurve.h"
     14#include "pmSpan.h"
     15#include "pmFootprintSpans.h"
     16#include "pmFootprint.h"
     17#include "pmPeaks.h"
     18#include "pmMoments.h"
     19#include "pmModelFuncs.h"
     20#include "pmModel.h"
     21#include "pmModelUtils.h"
     22#include "pmModelClass.h"
     23#include "pmSourceMasks.h"
     24#include "pmSourceExtendedPars.h"
     25#include "pmSourceDiffStats.h"
    926#include "pmSource.h"
     27#include "pmSourceFitModel.h"
     28#include "pmPSF.h"
     29#include "pmPSFtry.h"
     30#include "pmDetections.h"
     31
    1032#include "pmErrorCodes.h"
    1133

  • trunk/psModules/src/objects/pmSourceMatch.h

    r26076 r29004  
    11#ifndef PM_SOURCE_MATCH_H
    22#define PM_SOURCE_MATCH_H
    3 
    4 #include <pslib.h>
    53
    64/// Mask values for matched sources

  • trunk/psModules/src/objects/pmSourceMoments.c

    r27531 r29004  
    2222#include <strings.h>
    2323#include <pslib.h>
     24
    2425#include "pmHDU.h"
    2526#include "pmFPA.h"
    26 #include "pmFPAMaskWeight.h"
     27
     28#include "pmTrend2D.h"
     29#include "pmResiduals.h"
     30#include "pmGrowthCurve.h"
    2731#include "pmSpan.h"
     32#include "pmFootprintSpans.h"
    2833#include "pmFootprint.h"
    2934#include "pmPeaks.h"
    3035#include "pmMoments.h"
    31 #include "pmResiduals.h"
    32 #include "pmGrowthCurve.h"
    33 #include "pmTrend2D.h"
    34 #include "pmPSF.h"
     36#include "pmModelFuncs.h"
    3537#include "pmModel.h"
     38#include "pmModelUtils.h"
     39#include "pmModelClass.h"
     40#include "pmSourceMasks.h"
     41#include "pmSourceExtendedPars.h"
     42#include "pmSourceDiffStats.h"
    3643#include "pmSource.h"
    3744
    … …  
    5461# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)
    5562
     63static bool beVerbose = false;
     64void pmSourceMomentsSetVerbose(bool state){ beVerbose = state; }
     65
    5666bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal)
    5767{
    … …  
    6171    PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
    6272
    63     // use sky from moments if defined, 0.0 otherwise
     73    // use sky from moments if defined, 0.0 otherwise
     74
     75    // XXX this value comes from the sky model at the source center, and tends to over-estimate
     76    // the sky in the vicinity of bright sources.  we are better off assuming the model worked
     77    // well:
    6478    psF32 sky = 0.0;
    65     if (source->moments == NULL) {
    66         source->moments = pmMomentsAlloc();
    67     } else {
    68         sky = source->moments->Sky;
    69     }
     79    // XXX if (source->moments == NULL) {
     80    // XXX     source->moments = pmMomentsAlloc();
     81    // XXX } else {
     82    // XXX     sky = source->moments->Sky;
     83    // XXX }
    7084
    7185    // First Pass: calculate the first moments (these are subtracted from the coordinates below)
    … …  
    131145            psF32 wDiff = *vWgt;
    132146
    133             // skip pixels below specified significance level.  this is allowed, but should be
    134             // avoided -- the over-weights the wings of bright stars compared to those of faint
    135             // stars.
    136             if (PS_SQR(pDiff) < minSN2*wDiff) continue;
    137             // if (pDiff < 0) continue; // XXX : MWV says I should include < 0.0 valued points...
     147            // skip pixels below specified significance level.  for a PSFs, this
     148            // over-weights the wings of bright stars compared to those of faint stars.
     149            // for the estimator used for extended source analysis (where the window
     150            // function is allowed to be arbitrarily large), we need to clip to avoid
     151            // negative second moments.
     152            if (PS_SQR(pDiff) < minSN2*wDiff) continue; //
     153            if ((minSN > 0.0) && (pDiff < 0)) continue; //
    138154
    139155            // Apply a Gaussian window function.  Be careful with the window function.  S/N
    … …  
    197213    // Xn  = SUM (x - xc)^n * (z - sky)
    198214
     215    psF32 RF = 0.0;
     216    psF32 RH = 0.0;
     217    psF32 RS = 0.0;
    199218    psF32 XX = 0.0;
    200219    psF32 XY = 0.0;
    … …  
    244263            if (r > radius) continue;
    245264
    246             psF32 pDiff = *vPix - sky;
     265            psF32 fDiff = *vPix - sky;
     266            psF32 pDiff = fDiff;
    247267            psF32 wDiff = *vWgt;
    248268
    … …  
    257277            if (sigma > 0.0) {
    258278                // XXX a lot of extra flops; can we do pre-calculate?
    259                 psF32 z  = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2;
     279                // XXX we were re-calculating r2 (maybe the compiler caught this?)
     280                // psF32 z  = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2;
     281                psF32 z  = r2 * rsigma2;
    260282                assert (z >= 0.0);
    261283                psF32 weight  = exp(-z);
    … …  
    266288
    267289            Sum += pDiff;
     290
     291# if (1)
     292# if (0)
     293            if (fDiff < 0) continue;
     294# endif
     295            psF32 rf = r * fDiff;
     296            psF32 rh = sqrt(r) * fDiff;
     297            psF32 rs = fDiff;
     298# else
     299            psF32 rf = r * pDiff;
     300            psF32 rh = sqrt(r) * pDiff;
     301            psF32 rs = pDiff;
     302# endif
    268303
    269304            psF32 x = xDiff * pDiff;
    … …  
    284319            psF32 xyyy = xDiff * yyy / r2;
    285320            psF32 yyyy = yDiff * yyy / r2;
     321
     322            RF  += rf;
     323            RH  += rh;
     324            RS  += rs;
    286325
    287326            XX  += xx;
    … …  
    302341    }
    303342
     343    source->moments->Mrf = RF/RS;
     344    source->moments->Mrh = RH/RS;
     345
    304346    source->moments->Mxx = XX/Sum;
    305347    source->moments->Mxy = XY/Sum;
    … …  
    317359    source->moments->Myyyy = YYYY/Sum;
    318360
    319     // if (source->moments->Mxx < 0) {
    320     //  fprintf (stderr, "error: neg second moment??\n");
    321     // }
    322     // if (source->moments->Myy < 0) {
    323     //  fprintf (stderr, "error: neg second moment??\n");
    324     // }
    325 
    326     psTrace ("psModules.objects", 4, "Mxx: %f  Mxy: %f  Myy: %f  Mxxx: %f  Mxxy: %f  Mxyy: %f  Myyy: %f  Mxxxx: %f  Mxxxy: %f  Mxxyy: %f  Mxyyy: %f  Mxyyy: %f\n",
     361    // Calculate the Kron magnitude (make this block optional?)
     362    // float radKron = 2.5*source->moments->Mrf;
     363    float radKinner = 1.0*source->moments->Mrf;
     364    float radKron   = 2.5*source->moments->Mrf;
     365    float radKouter = 4.0*source->moments->Mrf;
     366
     367    int nKronPix = 0;
     368    Sum = Var = 0.0;
     369    float SumInner = 0.0;
     370    float SumOuter = 0.0;
     371
     372    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
     373
     374        psF32 yDiff = row - yCM;
     375        if (fabs(yDiff) > radKron) continue;
     376
     377        psF32 *vPix = source->pixels->data.F32[row];
     378        psF32 *vWgt = source->variance->data.F32[row];
     379        psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
     380
     381        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
     382            if (vMsk) {
     383                if (*vMsk & maskVal) {
     384                    vMsk++;
     385                    continue;
     386                }
     387                vMsk++;
     388            }
     389            if (isnan(*vPix)) continue;
     390
     391            psF32 xDiff = col - xCM;
     392            if (fabs(xDiff) > radKron) continue;
     393
     394            // radKron is just a function of (xDiff, yDiff)
     395            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
     396            psF32 r  = sqrt(r2);
     397
     398            psF32 pDiff = *vPix - sky;
     399            psF32 wDiff = *vWgt;
     400
     401            // skip pixels below specified significance level.  this is allowed, but should be
     402            // avoided -- the over-weights the wings of bright stars compared to those of faint
     403            // stars.
     404            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
     405
     406            if (r < radKron) {
     407                Sum += pDiff;
     408                Var += wDiff;
     409                nKronPix ++;
     410                // if (beVerbose) fprintf (stderr, "mome: %d %d  %f  %f  %f\n", col, row, sky, *vPix, Sum);
     411            }
     412
     413            if ((r > radKinner) && (r < radKron)) {
     414                SumInner += pDiff;
     415            }
     416            if ((r > radKron)  && (r < radKouter)) {
     417                SumOuter += pDiff;
     418            }
     419        }
     420    }
     421    source->moments->KronFlux = Sum;
     422    source->moments->KronFinner = SumInner;
     423    source->moments->KronFouter = SumOuter;
     424    source->moments->KronFluxErr = sqrt(Var);
     425
     426    psTrace ("psModules.objects", 4, "Mrf: %f  KronFlux: %f  Mxx: %f  Mxy: %f  Myy: %f  Mxxx: %f  Mxxy: %f  Mxyy: %f  Myyy: %f  Mxxxx: %f  Mxxxy: %f  Mxxyy: %f  Mxyyy: %f  Mxyyy: %f\n",
     427             source->moments->Mrf,   source->moments->KronFlux,
    327428             source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy,
    328429             source->moments->Mxxx,  source->moments->Mxxy,  source->moments->Mxyy,  source->moments->Myyy,

  • trunk/psModules/src/objects/pmSourcePhotometry.c

    r28517 r29004  
    2222#include "pmFPA.h"
    2323#include "pmFPAMaskWeight.h"
     24
     25#include "pmTrend2D.h"
     26#include "pmResiduals.h"
     27#include "pmGrowthCurve.h"
    2428#include "pmSpan.h"
     29#include "pmFootprintSpans.h"
    2530#include "pmFootprint.h"
    2631#include "pmPeaks.h"
    2732#include "pmMoments.h"
    28 #include "pmGrowthCurve.h"
    29 #include "pmResiduals.h"
    30 #include "pmTrend2D.h"
     33#include "pmModelFuncs.h"
     34#include "pmModel.h"
     35#include "pmModelUtils.h"
     36#include "pmModelClass.h"
     37#include "pmSourceMasks.h"
     38#include "pmSourceExtendedPars.h"
     39#include "pmSourceDiffStats.h"
     40#include "pmSource.h"
     41#include "pmSourceFitModel.h"
    3142#include "pmPSF.h"
    32 #include "pmModel.h"
    33 #include "pmSource.h"
    34 #include "pmModelClass.h"
     43#include "pmPSFtry.h"
     44
    3545#include "pmSourcePhotometry.h"
    3646
    … …  
    6676
    6777// XXX masked region should be (optionally) elliptical
    68 bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal)
     78bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal, psImageMaskType markVal)
    6979{
    7080    PS_ASSERT_PTR_NON_NULL(source, false);
    … …  
    122132        for (int i = 0; i < source->modelFits->n; i++) {
    123133            pmModel *model = source->modelFits->data[i];
     134            if (model->flags & PM_MODEL_STATUS_BADARGS) continue;
    124135            status = pmSourcePhotometryModel (&model->mag, NULL, model);
    125136            if (model == source->modelEXT) foundEXT = true;
    … …  
    145156    // measure the contribution of included pixels
    146157    if (mode & PM_SOURCE_PHOT_WEIGHT) {
    147         pmSourcePixelWeight (&source->pixWeight, model, source->maskObj, maskVal);
     158        pmSourcePixelWeight (&source->pixWeightNotBad, &source->pixWeightNotPoor, model, source->maskObj, maskVal, markVal);
    148159    }
    149160
    150161    // measure the contribution of included pixels
    151162    if (mode & PM_SOURCE_PHOT_DIFFSTATS) {
    152         pmSourceMeasureDiffStats (source, maskVal);
     163        pmSourceMeasureDiffStats (source, maskVal, markVal);
    153164    }
    154165
    … …  
    191202
    192203    // measure object aperture photometry
    193     status = pmSourcePhotometryAper  (&source->apMag, model, flux, mask, maskVal);
     204    status = pmSourcePhotometryAper  (&source->apMagRaw, model, flux, mask, maskVal);
    194205    if (!status) {
    195206        psTrace ("psModules.objects", 3, "fail mag : bad Ap Mag");
    … …  
    199210    // if the aper mag is NAN, the flux < 0.  this can happen for sources near the
    200211    // detection limits (esp near bright neighbors)
     212    source->apMag = source->apMagRaw;
    201213    if (isfinite (source->apMag) && isPSF && psf) {
    202214        if (psf->growth && (mode & PM_SOURCE_PHOT_GROWTH)) {
    203             source->apMag += pmGrowthCurveCorrect (psf->growth, source->apRadius);
     215            source->apMag = source->apMagRaw + pmGrowthCurveCorrect (psf->growth, source->apRadius);
    204216        }
    205217        if (mode & PM_SOURCE_PHOT_APCORR) {
    206218            // XXX this should be removed -- we no longer fit for the 'sky bias'
     219            // XXX is this happening???
    207220            rflux   = pow (10.0, 0.4*source->psfMag);
     221            psAssert (psf->skyBias == 0.0, "sky bias not 0");
     222            psAssert (psf->skySat == 0.0, "sky sat not 0");
    208223            source->apMag -= PS_SQR(source->apRadius)*rflux * psf->skyBias + psf->skySat / rflux;
    209224        }
    … …  
    257272    PS_ASSERT_PTR_NON_NULL(image, false);
    258273    PS_ASSERT_PTR_NON_NULL(mask, false);
    259     PS_ASSERT_PTR_NON_NULL(model, false);
     274
     275    if (DO_SKY) {
     276        PS_ASSERT_PTR_NON_NULL(model, false);
     277    }
    260278
    261279    float apSum = 0;
    … …  
    271289    psF32 **imData = image->data.F32;
    272290    psImageMaskType **mkData = mask->data.PS_TYPE_IMAGE_MASK_DATA;
     291    int nAperPix = 0;
    273292
    274293    // measure apMag
    275     for (int ix = 0; ix < image->numCols; ix++) {
    276         for (int iy = 0; iy < image->numRows; iy++) {
     294    for (int iy = 0; iy < image->numRows; iy++) {
     295        for (int ix = 0; ix < image->numCols; ix++) {
    277296            if (mkData[iy][ix] & maskVal)
    278297                continue;
    279298            apSum += imData[iy][ix] - sky;
     299            nAperPix ++;
     300            // fprintf (stderr, "aper: %d %d  %f  %f  %f\n", ix, iy, sky, imData[iy][ix], apSum);
    280301        }
    281302    }
    … …  
    290311
    291312// return source aperture magnitude
    292 bool pmSourcePixelWeight (float *pixWeight, pmModel *model, psImage *mask, psImageMaskType maskVal)
    293 {
    294     PS_ASSERT_PTR_NON_NULL(pixWeight, false);
     313bool pmSourcePixelWeight (float *pixWeightNotBad, float *pixWeightNotPoor, pmModel *model, psImage *mask, psImageMaskType maskVal, psImageMaskType markVal)
     314{
     315    PS_ASSERT_PTR_NON_NULL(pixWeightNotBad, false);
     316    PS_ASSERT_PTR_NON_NULL(pixWeightNotPoor, false);
    295317    PS_ASSERT_PTR_NON_NULL(mask, false);
    296318    PS_ASSERT_PTR_NON_NULL(model, false);
    297319
    298320    float modelSum = 0;
    299     float validSum = 0;
     321    float notBadSum = 0;
     322    float notPoorSum = 0;
    300323    float sky = 0;
    301324    float value;
    … …  
    305328    int dY, DY, NY;
    306329
    307     *pixWeight = 0.0;
     330    *pixWeightNotBad = 0.0;
     331    *pixWeightNotPoor = 0.0;
    308332
    309333    // we only care about the value of the object model, not the local sky
    … …  
    345369
    346370            // for the full model, add all points
    347             value = model->modelFunc (NULL, params, coord) - sky;
     371            value = fabs(model->modelFunc (NULL, params, coord) - sky);
    348372            modelSum += value;
    349373
    350374            // include count only the unmasked pixels within the image area
    351             if (mx < 0)
    352                 continue;
    353             if (my < 0)
    354                 continue;
    355             if (mx >= NX)
    356                 continue;
    357             if (my >= NY)
    358                 continue;
    359             if (mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskVal)
    360                 continue;
    361 
    362             validSum += value;
     375            if (mx < 0) continue;
     376            if (my < 0) continue;
     377            if (mx >= NX) continue;
     378            if (my >= NY) continue;
     379
     380            if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskVal)) {
     381                notBadSum += value;
     382            }
     383            if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & ~markVal)) {
     384                notPoorSum += value;
     385            }
    363386        }
    364387    }
    365388    psFree (coord);
    366389
    367     if (validSum <= 0)
    368         return false;
    369 
    370     *pixWeight = validSum / modelSum;
     390    *pixWeightNotBad  = notBadSum  / modelSum;
     391    *pixWeightNotPoor = notPoorSum / modelSum;
     392
    371393    return (true);
    372394}
    … …  
    374396# define FLUX_LIMIT 3.0
    375397
    376 // return source aperture magnitude
    377 bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal)
     398// measure stats that may be using in difference images for distinguishing real sources from bad residuals
     399bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal, psImageMaskType markVal)
    378400{
    379401    PS_ASSERT_PTR_NON_NULL(source, false);
    … …  
    399421    for (int iy = 0; iy < flux->numRows; iy++) {
    400422        for (int ix = 0; ix < flux->numCols; ix++) {
     423            // only count up the stats in the unmarked region (ie, the aperture)
     424            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & markVal) {
     425                continue;
     426            }
    401427            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) {
    402428                nMask ++;

  • trunk/psModules/src/objects/pmSourcePhotometry.h

    r28426 r29004  
    5252
    5353bool pmSourceMagnitudesInit (psMetadata *config);
    54 bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal);
    55 bool pmSourcePixelWeight (float *pixWeight, pmModel *model, psImage *mask, psImageMaskType maskVal);
     54bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal, psImageMaskType markVal);
     55
     56bool pmSourcePixelWeight (float *pixWeightNotBad, float *pixWeightNotPoor, pmModel *model, psImage *mask, psImageMaskType maskVal, psImageMaskType markVal);
     57
    5658bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, psImageMaskType maskVal, const float covarFactor);
    5759
    58 bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal);
     60bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
    5961
    6062double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);

  • trunk/psModules/src/objects/pmSourcePlotApResid.c

    r26893 r29004  
    2525#include "pmFPAview.h"
    2626#include "pmFPAfile.h"
     27
     28#include "pmTrend2D.h"
     29#include "pmResiduals.h"
     30#include "pmGrowthCurve.h"
    2731#include "pmSpan.h"
     32#include "pmFootprintSpans.h"
    2833#include "pmFootprint.h"
    2934#include "pmPeaks.h"
    3035#include "pmMoments.h"
    31 #include "pmResiduals.h"
    32 #include "pmGrowthCurve.h"
    33 #include "pmTrend2D.h"
     36#include "pmModelFuncs.h"
     37#include "pmModel.h"
     38#include "pmModelUtils.h"
     39#include "pmModelClass.h"
     40#include "pmSourceMasks.h"
     41#include "pmSourceExtendedPars.h"
     42#include "pmSourceDiffStats.h"
     43#include "pmSource.h"
     44#include "pmSourceFitModel.h"
    3445#include "pmPSF.h"
    35 #include "pmModel.h"
    3646#include "pmDetections.h"
    37 #include "pmSource.h"
     47
    3848#include "pmSourcePlots.h"
    3949#include "pmKapaPlots.h"

  • trunk/psModules/src/objects/pmSourcePlotMoments.c

    r26893 r29004  
    2828#include "pmFPAview.h"
    2929#include "pmFPAfile.h"
     30
     31#include "pmTrend2D.h"
     32#include "pmResiduals.h"
     33#include "pmGrowthCurve.h"
    3034#include "pmSpan.h"
     35#include "pmFootprintSpans.h"
    3136#include "pmFootprint.h"
    3237#include "pmPeaks.h"
    3338#include "pmMoments.h"
    34 #include "pmGrowthCurve.h"
    35 #include "pmResiduals.h"
    36 #include "pmTrend2D.h"
     39#include "pmModelFuncs.h"
     40#include "pmModel.h"
     41#include "pmModelUtils.h"
     42#include "pmModelClass.h"
     43#include "pmSourceMasks.h"
     44#include "pmSourceExtendedPars.h"
     45#include "pmSourceDiffStats.h"
     46#include "pmSource.h"
     47#include "pmSourceFitModel.h"
    3748#include "pmPSF.h"
    38 #include "pmModel.h"
    3949#include "pmDetections.h"
    40 #include "pmSource.h"
    4150#include "pmSourcePlots.h"
    4251#include "pmKapaPlots.h"

  • trunk/psModules/src/objects/pmSourcePlotPSFModel.c

    r26893 r29004  
    2828#include "pmFPAview.h"
    2929#include "pmFPAfile.h"
     30
     31
     32#include "pmTrend2D.h"
     33#include "pmResiduals.h"
     34#include "pmGrowthCurve.h"
    3035#include "pmSpan.h"
     36#include "pmFootprintSpans.h"
    3137#include "pmFootprint.h"
    3238#include "pmPeaks.h"
    3339#include "pmMoments.h"
    34 #include "pmGrowthCurve.h"
    35 #include "pmResiduals.h"
    36 #include "pmTrend2D.h"
     40#include "pmModelFuncs.h"
     41#include "pmModel.h"
     42#include "pmModelUtils.h"
     43#include "pmModelClass.h"
     44#include "pmSourceMasks.h"
     45#include "pmSourceExtendedPars.h"
     46#include "pmSourceDiffStats.h"
     47#include "pmSource.h"
     48#include "pmSourceFitModel.h"
    3749#include "pmPSF.h"
    38 #include "pmModel.h"
    3950#include "pmDetections.h"
    40 #include "pmSource.h"
    4151#include "pmSourcePlots.h"
    4252#include "pmKapaPlots.h"

  • trunk/psModules/src/objects/pmSourceSky.c

    r23187 r29004  
    2121#include <string.h>
    2222#include <pslib.h>
     23
    2324#include "pmHDU.h"
    2425#include "pmFPA.h"
    25 #include "pmFPAMaskWeight.h"
     26
     27#include "pmTrend2D.h"
     28#include "pmResiduals.h"
     29#include "pmGrowthCurve.h"
    2630#include "pmSpan.h"
     31#include "pmFootprintSpans.h"
    2732#include "pmFootprint.h"
    2833#include "pmPeaks.h"
    2934#include "pmMoments.h"
    30 #include "pmResiduals.h"
    31 #include "pmGrowthCurve.h"
    32 #include "pmTrend2D.h"
    33 #include "pmPSF.h"
     35#include "pmModelFuncs.h"
    3436#include "pmModel.h"
     37#include "pmModelUtils.h"
     38#include "pmModelClass.h"
     39#include "pmSourceMasks.h"
     40#include "pmSourceExtendedPars.h"
     41#include "pmSourceDiffStats.h"
    3542#include "pmSource.h"
     43
    3644#include "pmSourceSky.h"
    3745

  • trunk/psModules/src/objects/pmSourceUtils.c

    r24578 r29004  
    2121#include <string.h>
    2222#include <pslib.h>
     23
    2324#include "pmHDU.h"
    2425#include "pmFPA.h"
    25 #include "pmFPAMaskWeight.h"
     26
     27#include "pmTrend2D.h"
     28#include "pmResiduals.h"
     29#include "pmGrowthCurve.h"
    2630#include "pmSpan.h"
     31#include "pmFootprintSpans.h"
    2732#include "pmFootprint.h"
    2833#include "pmPeaks.h"
    2934#include "pmMoments.h"
    30 #include "pmResiduals.h"
    31 #include "pmGrowthCurve.h"
    32 #include "pmTrend2D.h"
    33 #include "pmPSF.h"
     35#include "pmModelFuncs.h"
    3436#include "pmModel.h"
     37#include "pmModelUtils.h"
     38#include "pmModelClass.h"
     39#include "pmSourceMasks.h"
     40#include "pmSourceExtendedPars.h"
     41#include "pmSourceDiffStats.h"
    3542#include "pmSource.h"
     43
    3644#include "pmSourceUtils.h"
    3745

  • trunk/psModules/src/objects/pmSourceVisual.c

    r26260 r29004  
    44
    55#include <pslib.h>
     6#include "pmHDU.h"
     7#include "pmFPA.h"
     8
    69#include "pmTrend2D.h"
     10#include "pmResiduals.h"
     11#include "pmGrowthCurve.h"
     12#include "pmSpan.h"
     13#include "pmFootprintSpans.h"
     14#include "pmFootprint.h"
     15#include "pmPeaks.h"
     16#include "pmMoments.h"
     17#include "pmModelFuncs.h"
     18#include "pmModel.h"
     19#include "pmModelUtils.h"
     20#include "pmModelClass.h"
     21#include "pmSourceMasks.h"
     22#include "pmSourceExtendedPars.h"
     23#include "pmSourceDiffStats.h"
     24#include "pmSource.h"
     25#include "pmSourceFitModel.h"
    726#include "pmPSF.h"
    827#include "pmPSFtry.h"
    9 #include "pmSource.h"
     28#include "pmDetections.h"
     29
    1030#include "pmSourceVisual.h"
    1131
    … …  
    1333#include <kapa.h>
    1434#include "pmVisual.h"
     35#include "pmVisualUtils.h"
    1536
    1637// functions used to visualize the analysis as it goes
    … …  
    3455    Graphdata graphdata;
    3556
    36     if (!pmVisualIsVisual()) return true;
     57    if (!pmVisualTestLevel("psphot.psf.metric", 2)) return true;
    3758
    3859    if (kapa1 == -1) {
    … …  
    118139    Graphdata graphdata;
    119140
    120     if (!pmVisualIsVisual()) return true;
     141    if (!pmVisualTestLevel("psphot.psf.subpix", 3)) return true;
    121142
    122143    if (kapa1 == -1) {
    … …  
    280301bool pmSourceVisualShowModelFits (pmPSF *psf, psArray *sources, psImageMaskType maskVal) {
    281302
    282     if (!pmVisualIsVisual()) return true;
     303    if (!pmVisualTestLevel("psphot.psf.fits", 2)) return true;
    283304
    284305    if (kapa2 == -1) {
    … …  
    360381bool pmSourceVisualShowModelFit (pmSource *source) {
    361382
    362     if (!pmVisualIsVisual()) return true;
     383    if (!pmVisualTestLevel("psphot.psf.fitresid", 2)) return true;
     384
    363385    if (!source->pixels) return false;
    364386    if (!source->modelFlux) return false;
    … …  
    404426    Graphdata graphdata;
    405427
    406     if (!pmVisualIsVisual() || !plotPSF) return true;
     428    if (!plotPSF) return true;
     429    if (!pmVisualTestLevel("psphot.psf.resid", 2)) return true;
    407430
    408431    if (kapa1 == -1) {

  • trunk/psModules/src/objects/pmSpan.h

    r20945 r29004  
    1010# ifndef PM_SPAN_H
    1111# define PM_SPAN_H
    12 
    13 #include <pslib.h>
    1412
    1513/// @addtogroup Objects Object Detection / Analysis Functions

  • trunk/psModules/src/objects/pmTrend2D.h

    r25754 r29004  
    1212# ifndef PM_TREND_2D_H
    1313# define PM_TREND_2D_H
    14 
    15 #include <pslib.h>
    1614
    1715/// @addtogroup Objects Object Detection / Analysis Functions

  • trunk/psModules/src/psmodules.h

    r28043 r29004  
    1010#include <pmKapaPlots.h>
    1111#include <pmVisual.h>
     12#include <pmVisualUtils.h>
    1213#include <ippStages.h>
    1314#include <ippDiffMode.h>
    … …  
    113114
    114115// the following headers are from psModule:objects
     116#include <pmTrend2D.h>
     117#include <pmResiduals.h>
     118#include <pmGrowthCurve.h>
     119
    115120#include <pmSpan.h>
    116121#include <pmFootprintSpans.h>
    … …  
    119124#include <pmDetections.h>
    120125#include <pmMoments.h>
     126
     127#include <pmModelFuncs.h>
     128#include <pmModel.h>
     129
     130#include <pmSourceMasks.h>
    121131#include <pmSourceExtendedPars.h>
    122132#include <pmSourceDiffStats.h>
    123 #include <pmResiduals.h>
    124 #include <pmGrowthCurve.h>
    125 #include <pmTrend2D.h>
     133#include <pmSource.h>
     134#include <pmSourceFitModel.h>
    126135#include <pmPSF.h>
    127 #include <pmModel.h>
    128 #include <pmSourceMasks.h>
    129 #include <pmSource.h>
     136#include <pmPSFtry.h>
    130137#include <pmPhotObj.h>
    131138#include <pmSourceUtils.h>
    132139#include <pmSourceIO.h>
    133140#include <pmSourceSky.h>
    134 #include <pmSourceFitModel.h>
    135141#include <pmSourceFitSet.h>
    136142#include <pmSourceContour.h>
    137143#include <pmSourcePlots.h>
    138144#include <pmPSF_IO.h>
    139 #include <pmPSFtry.h>
    140145#include <pmModelClass.h>
    141146#include <pmModelUtils.h>
    … …  
    144149#include <pmSourceMatch.h>
    145150#include <pmDetEff.h>
     151#include <pmPCMdata.h>
    146152
    147153// The following headers are from random locations, here because they cross bounds
Note: See TracChangeset for help on using the changeset viewer.