IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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

Many changes:

  • psphot-related issues

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

  • ppSub-related:

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

Location:
trunk/psModules
Files:
1 deleted
83 edited
11 copied

Legend:

Unmodified
Added
Removed

  • trunk/psModules

    • Property svn:mergeinfo deleted

  • 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
Note: See TracChangeset for help on using the changeset viewer.