IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 28782


Ignore:
Timestamp:
Jul 29, 2010, 2:36:33 PM (16 years ago)
Author:
eugene
Message:

attempting to get a better starting guess; attempting to downweight neighbor detections

Location:
branches/eam_branches/ipp-20100621/psphot/src
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20100621/psphot/src/psphot.h

    r28702 r28782  
    180180
    181181// functions to set the correct source pixels
    182 bool            psphotInitRadiusPSF(const psMetadata *recipe, const psMetadata *analysis, const pmModelType type);
     182bool            psphotInitRadiusPSF (psMetadata *recipe, pmReadout *readout);
     183bool            psphotInitRadiusEXT (psMetadata *recipe, pmReadout *readout);
    183184
    184185bool            psphotCheckRadiusPSF (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal);
    185186bool            psphotCheckRadiusPSFBlend (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal, float dR);
    186 bool            psphotInitRadiusEXT (psMetadata *recipe, pmModelType type);
    187 bool            psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal);
    188 float           psphotSetRadiusEXT (pmReadout *readout, pmSource *source, psImageMaskType markVal);
     187bool            psphotSetRadiusFootprint (float *radius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float factor);
     188bool            psphotSetRadiusModel (pmModel *model, pmReadout *readout, pmSource *source, psImageMaskType markVal, bool deep);
    189189
    190190bool            psphotDumpMoments (psMetadata *recipe, psArray *sources);
    … …  
    202202//  functions to support the source fitting process
    203203bool            psphotInitLimitsPSF (psMetadata *recipe, pmReadout *readout);
    204 bool            psphotInitLimitsEXT (psMetadata *recipe);
     204bool            psphotInitLimitsEXT (psMetadata *recipe, pmReadout *readout);
    205205bool            psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal);
    206206bool            psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal);
    … …  
    426426bool psphotStackObjectsUnifyPosition (psArray *objects);
    427427
    428 bool psphotFitSersicIndexPCM (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize);
     428bool psphotFitSersicIndex (pmModel *model, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal);
     429
     430bool psphotFitSersicIndexPCM (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize);
    429431pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize);
    430432

  • branches/eam_branches/ipp-20100621/psphot/src/psphotAddNoise.c

    r28013 r28782  
    3434
    3535    bool status = false;
    36     psEllipseShape oldshape;
    37     psEllipseShape newshape;
    38     psEllipseAxes axes;
    3936
    4037    // find the currently selected readout
    … …  
    7875    }
    7976
     77    psphotSaveImage (NULL, readout->variance, "test.var0.fits");
     78
    8079    // loop over all source
    8180    for (int i = 0; i < sources->n; i++) {
    … …  
    8685        if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue;
    8786
    88         // select appropriate model
    89         pmModel *model = pmSourceGetModel (NULL, source);
    90         if (model == NULL) continue;  // model must be defined
    91 
    92         if (add) {
    93             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]);
    94         } else {
    95             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]);
    96         }
    97 
    98         psF32 *PAR = model->params->data.F32;
    99 
    100         // save original values
    101         float oldI0  = PAR[PM_PAR_I0];
    102         oldshape.sx  = PAR[PM_PAR_SXX];
    103         oldshape.sy  = PAR[PM_PAR_SYY];
    104         oldshape.sxy = PAR[PM_PAR_SXY];
    105 
    106         // XXX can this be done more intelligently?
    107         if (oldI0 == 0.0) continue;
    108         if (!isfinite(oldI0)) continue;
    109 
    110         // increase size and height of source
    111         axes = psEllipseShapeToAxes (oldshape, 20.0);
    112         axes.major *= SIZE;
    113         axes.minor *= SIZE;
    114         newshape = psEllipseAxesToShape (axes);
    115         PAR[PM_PAR_I0]  = FACTOR*oldI0;
    116         PAR[PM_PAR_SXX] = newshape.sx;
    117         PAR[PM_PAR_SYY] = newshape.sy;
    118         PAR[PM_PAR_SXY] = newshape.sxy;
    119 
    120         // XXX if we use pmSourceOp, the size (and possibly Io) will not be respected
    121         pmSourceOp (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, add, maskVal, 0, 0);
    122 
    123         // restore original values
    124         PAR[PM_PAR_I0]  = oldI0;
    125         PAR[PM_PAR_SXX] = oldshape.sx;
    126         PAR[PM_PAR_SYY] = oldshape.sy;
    127         PAR[PM_PAR_SXY] = oldshape.sxy;
     87        pmSourceNoiseOp (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, FACTOR, SIZE, add, maskVal, 0, 0);
    12888    }
    12989    if (add) {
    … …  
    13292        psLogMsg ("psphot.noise", PS_LOG_INFO, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
    13393    }
     94
     95    psphotSaveImage (NULL, readout->variance, "test.var1.fits");
    13496    return true;
    13597}

  • branches/eam_branches/ipp-20100621/psphot/src/psphotBlendFit.c

    r28657 r28782  
    9090
    9191    psphotInitLimitsPSF (recipe, readout);
    92     psphotInitLimitsEXT (recipe);
    93     psphotInitRadiusPSF (recipe, readout->analysis, psf->type);
     92    psphotInitLimitsEXT (recipe, readout);
     93    psphotInitRadiusPSF (recipe, readout);
    9494
    9595    // starts the timer, sets up the array of fitSets

  • branches/eam_branches/ipp-20100621/psphot/src/psphotExtendedSourceFits.c

    r28687 r28782  
    3939    int NplainPass = 0;
    4040    bool savePics = false;
     41    float radius;
    4142
    4243    // find the currently selected readout
    … …  
    164165        Next ++;
    165166
     167        // set the radius based on the footprint (also sets the mask pixels)
     168        if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) return false;
     169
     170        // XXX note that this changes the source moments that are published...
     171        // recalculate the source moments using the larger extended-source moments radius
     172        // at this stage, skip Gaussian windowing, and do not clip pixels by S/N
     173        // this uses the footprint to judge both radius and aperture?
     174        // XXX save the psf-based moments for output
     175        if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) {
     176            // subtract the best fit from the object, leave local sky
     177            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     178            // XXX raise an error of some kind
     179            continue;
     180        }
     181
    166182        // save the modelFlux here in case we need to subtract it (for failure)
    167183        psImage *modelFluxStart = psMemIncrRefCounter (source->modelFlux);
    … …  
    242258
    243259          // test for fit quality / result
     260          modelFit->fitRadius = radius;
    244261          psArrayAdd (source->modelFits, 4, modelFit);
    245262

  • branches/eam_branches/ipp-20100621/psphot/src/psphotGuessModels.c

    r28405 r28782  
    8080
    8181    // setup the PSF fit radius details
    82     psphotInitRadiusPSF (recipe, readout->analysis, psf->type);
     82    psphotInitRadiusPSF (recipe, readout);
    8383
    8484    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)

  • branches/eam_branches/ipp-20100621/psphot/src/psphotRadiusChecks.c

    r28440 r28782  
    88                                        // and a per-object radius is calculated)
    99
    10 bool psphotInitRadiusPSF(const psMetadata *recipe, const psMetadata *analysis, const pmModelType type) {
     10bool psphotInitRadiusPSF(psMetadata *recipe, pmReadout *readout) {
    1111
    1212    bool status = true;
    … …  
    1515    PSF_FIT_PADDING = psMetadataLookupF32(&status, recipe, "PSF_FIT_PADDING");
    1616
    17     PSF_FIT_RADIUS =  psMetadataLookupF32(&status, analysis, "PSF_FIT_RADIUS");
     17    PSF_FIT_RADIUS =  psMetadataLookupF32(&status, readout->analysis, "PSF_FIT_RADIUS");
    1818    if (!status) {
    1919        PSF_FIT_RADIUS = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS");
    2020    }
    2121
    22     PSF_APERTURE =  psMetadataLookupF32(&status, analysis, "PSF_APERTURE");
     22    PSF_APERTURE =  psMetadataLookupF32(&status, readout->analysis, "PSF_APERTURE");
    2323    if (!status) {
    2424        PSF_APERTURE =  psMetadataLookupF32(&status, recipe, "PSF_APERTURE");
    … …  
    2828
    2929    if (PSF_FIT_RADIUS == 0.0) {
    30         float gaussSigma = psMetadataLookupF32(&status, analysis, "MOMENTS_GAUSS_SIGMA");
     30        float gaussSigma = psMetadataLookupF32(&status, readout->analysis, "MOMENTS_GAUSS_SIGMA");
    3131        if (!status) {
    3232            gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA");
    … …  
    3737
    3838    if (PSF_APERTURE == 0.0) {
    39         float gaussSigma = psMetadataLookupF32(&status, analysis, "MOMENTS_GAUSS_SIGMA");
     39        float gaussSigma = psMetadataLookupF32(&status, readout->analysis, "MOMENTS_GAUSS_SIGMA");
    4040        if (!status) {
    4141            gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA");
    … …  
    122122}
    123123
     124static float EXT_FIT_SKY_SIG;
    124125static float EXT_FIT_NSIGMA;
    125126static float EXT_FIT_PADDING;
    126127static float EXT_FIT_MAX_RADIUS;
    127128
    128 bool psphotInitRadiusEXT (psMetadata *recipe, pmModelType type) {
     129bool psphotInitRadiusEXT (psMetadata *recipe, pmReadout *readout) {
    129130
    130131    bool status;
    … …  
    134135    EXT_FIT_MAX_RADIUS = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_RADIUS");
    135136
     137    float skyMean  = psMetadataLookupF32 (&status, readout->analysis, "SKY_MEAN");
     138    float skyStdev = psMetadataLookupF32 (&status, readout->analysis, "SKY_STDEV");
     139
     140    fprintf (stderr, "sky: %f +/- %f\n", skyMean, skyStdev);
     141
     142    EXT_FIT_SKY_SIG = skyStdev;
     143
    136144    return true;
    137145}
    138146
    139147// call this function whenever you (re)-define the EXT model
    140 float psphotSetRadiusEXT (pmReadout *readout, pmSource *source, psImageMaskType markVal) {
     148bool psphotSetRadiusFootprint (float *radius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float factor) {
    141149
    142150    psAssert (source, "source not defined??");
    … …  
    146154
    147155    // set the radius based on the footprint:
    148     if (!peak->footprint) goto escape;
     156    if (!peak->footprint) return false;
    149157    pmFootprint *footprint = peak->footprint;
    150     if (!footprint->spans) goto escape;
    151     if (footprint->spans->n < 1) goto escape;
     158    if (!footprint->spans) return false;
     159    if (footprint->spans->n < 1) return false;
    152160
    153161    // find the max radius
    154     float radius = 0.0;
     162    float rawRadius = 0.0;
    155163    for (int j = 0; j < footprint->spans->n; j++) {
    156164        pmSpan *span = footprint->spans->data[j];
    … …  
    160168        float dX1 = span->x1 - peak->xf;
    161169
    162         radius = PS_MAX (radius, hypot(dY, dX0));
    163         radius = PS_MAX (radius, hypot(dY, dX1));
    164     }
    165 
    166     radius += EXT_FIT_PADDING;
    167     if (isnan(radius)) psAbort("error in radius");
    168 
    169     radius = PS_MIN (radius, EXT_FIT_MAX_RADIUS);
     170        rawRadius = PS_MAX (rawRadius, hypot(dY, dX0));
     171        rawRadius = PS_MAX (rawRadius, hypot(dY, dX1));
     172    }
     173    if (isnan(rawRadius)) return false;
     174    rawRadius = PS_MIN (factor*rawRadius + EXT_FIT_PADDING, EXT_FIT_MAX_RADIUS);
    170175
    171176    // redefine the pixels if needed
    172     pmSourceRedefinePixels (source, readout, peak->xf, peak->yf, radius);
    173 
    174     // set the mask to flag the excluded pixels
    175     psImageKeepCircle (source->maskObj, peak->xf, peak->yf, radius, "OR", markVal);
    176     return radius;
    177 
    178 escape:
    179     return NAN;
    180     // bool result = psphotCheckRadiusEXT (readout, source, model, markVal);
    181     // return result;
     177    pmSourceRedefinePixels (source, readout, peak->xf, peak->yf, rawRadius);
     178
     179    // set the mask to flag the excluded pixels
     180    psImageKeepCircle (source->maskObj, peak->xf, peak->yf, rawRadius, "OR", markVal);
     181
     182    *radius = rawRadius;
     183    return true;
    182184}
    183185
    184186// alternative EXT radius based on model guess (for use without footprints)
    185 bool psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal) {
    186 
    187     psAbort ("do not use this function");
     187bool psphotSetRadiusModel (pmModel *model, pmReadout *readout, pmSource *source, psImageMaskType markVal, bool deep) {
    188188
    189189    psF32 *PAR = model->params->data.F32;
    … …  
    193193
    194194    // set the fit radius based on the object flux limit and the model
    195     float rawRadius = model->modelRadius (model->params, EXT_FIT_NSIGMA*moments->dSky);
    196 
    197     model->fitRadius = rawRadius + EXT_FIT_PADDING;
    198     if (isnan(model->fitRadius)) psAbort("error in radius");
     195    float flux = deep ? EXT_FIT_NSIGMA*EXT_FIT_SKY_SIG : 0.1 * model->params->data.F32[PM_PAR_I0];
     196
     197    float rawRadius = model->modelRadius (model->params, flux);
     198    if (isnan(rawRadius)) return false;
     199
     200    rawRadius = PS_MIN (rawRadius + EXT_FIT_PADDING, EXT_FIT_MAX_RADIUS);
     201    model->fitRadius = rawRadius;
    199202
    200203    // redefine the pixels if needed
    201     bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius);
     204    pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius);
    202205
    203206    // set the mask to flag the excluded pixels
    204207    psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius, "OR", markVal);
    205     return status;
    206 }
     208    return true;
     209}

  • branches/eam_branches/ipp-20100621/psphot/src/psphotSourceFits.c

    r28702 r28782  
    11# include "psphotInternal.h"
    2 bool psphotFitSersicIndex (pmSource *source, pmModel *model, pmSourceFitOptions *fitOptions, psImageMaskType maskVal);
    32
    43// given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful
    … …  
    1918
    2019    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "fitted %5d psf, %5d blend, %5d ext, %5d dbl : %6.2f sec\n",
    21              NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits"));
     20              NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits"));
    2221    return true;
    2322}
    … …  
    205204}
    206205
     206// save a local, static copy of the EXT model type so we don't have to lookup for each object
    207207static pmModelType modelTypeEXT;
    208208
    209 bool psphotInitLimitsEXT (psMetadata *recipe) {
     209bool psphotInitLimitsEXT (psMetadata *recipe, pmReadout *readout) {
    210210
    211211    bool status;
    … …  
    214214    char *modelNameEXT = psMetadataLookupStr (&status, recipe, "EXT_MODEL");
    215215    modelTypeEXT = pmModelClassGetType (modelNameEXT);
    216     psphotInitRadiusEXT (recipe, modelTypeEXT);
     216
     217    psphotInitRadiusEXT (recipe, readout);
    217218
    218219    return true;
    … …  
    221222bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *newSources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
    222223
     224    float radius;
    223225    bool okEXT, okDBL;
    224226    float chiEXT, chiDBL;
    … …  
    230232
    231233    // skip the source if we don't think it is extended
     234    // XXX are these robust, or do we have better info from the source-size analysis??
    232235    if (source->type == PM_SOURCE_TYPE_UNKNOWN) return false;
    233236    if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
    … …  
    235238
    236239    // set the radius based on the footprint (also sets the mask pixels)
    237     float radius = psphotSetRadiusEXT (readout, source, markVal);
     240    if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) return false;
    238241
    239242    // XXX note that this changes the source moments that are published...
     243    // XXX all published moments should use the same measurement
    240244    // recalculate the source moments using the larger extended-source moments radius
    241245    // at this stage, skip Gaussian windowing, and do not clip pixels by S/N
    242246    // this uses the footprint to judge both radius and aperture?
     247    // XXX save the psf-based moments for output
    243248    if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) return false;
    244249
    … …  
    250255    // this temporary source is used as a place-holder by the psphotEval functions below
    251256    tmpSrc = pmSourceAlloc ();
    252 
    253     // XXX need to handle failures better here
    254     EXT = psphotFitEXT (readout, source, fitOptions, modelTypeEXT, maskVal, markVal);
    255     if (!EXT) goto escape;
    256     if (!isfinite(EXT->params->data.F32[PM_PAR_I0])) goto escape;
    257 
    258     okEXT = psphotEvalEXT (tmpSrc, EXT);
    259     chiEXT = EXT ? EXT->chisq / EXT->nDOF : NAN;
    260 
    261     // DBL will always be defined, but DBL->data[n] might not
    262     DBL = psphotFitDBL (readout, source, fitOptions, maskVal, markVal);
    263     if (!DBL) goto escape;
    264     if (!DBL->n) goto escape;
    265 
    266     okDBL  = psphotEvalDBL (tmpSrc, DBL->data[0]);
    267     okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]);
    268     // XXX should I keep / save the flags set in the eval functions?
     257    {
     258        // DBL will always be defined, but DBL->data[n] might not
     259        DBL = psphotFitDBL (readout, source, fitOptions, maskVal, markVal);
     260        if (!DBL) goto escape;
     261        if (!DBL->n) goto escape;
     262
     263        okDBL  = psphotEvalDBL (tmpSrc, DBL->data[0]);
     264        okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]);
     265        // XXX should I keep / save the flags set in the eval functions?
     266
     267        // correct first model chisqs for flux trend
     268        chiDBL = NAN;
     269        ONE = DBL->data[0];
     270        if (ONE) {
     271            if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
     272            chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
     273            ONE->chisqNorm = ONE->chisq / chiTrend;
     274            chiDBL = ONE->chisq / ONE->nDOF; // save chisq for double-star/galaxy comparison
     275            ONE->fitRadius = radius;
     276        }
     277
     278        // correct second model chisqs for flux trend
     279        ONE = DBL->data[1];
     280        if (ONE) {
     281            if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
     282            chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
     283            ONE->chisqNorm = ONE->chisq / chiTrend;
     284            ONE->fitRadius = radius;
     285        }
     286    }
     287
     288    {
     289        // XXX need to handle failures better here
     290        EXT = psphotFitEXT (readout, source, fitOptions, modelTypeEXT, maskVal, markVal);
     291        if (!EXT) goto escape;
     292        if (!isfinite(EXT->params->data.F32[PM_PAR_I0])) goto escape;
     293
     294        okEXT = psphotEvalEXT (tmpSrc, EXT);
     295        chiEXT = EXT ? EXT->chisq / EXT->nDOF : NAN;
     296    }
    269297
    270298    // clear the circular mask
    271299    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    272 
    273     // correct first model chisqs for flux trend
    274     chiDBL = NAN;
    275     ONE = DBL->data[0];
    276     if (ONE) {
    277         if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
    278       chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
    279       ONE->chisqNorm = ONE->chisq / chiTrend;
    280       chiDBL = ONE->chisq / ONE->nDOF; // save chisq for double-star/galaxy comparison
    281     }
    282 
    283     // correct second model chisqs for flux trend
    284     ONE = DBL->data[1];
    285     if (ONE) {
    286         if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
    287       chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
    288       ONE->chisqNorm = ONE->chisq / chiTrend;
    289     }
    290300
    291301    psFree (tmpSrc);
    … …  
    314324
    315325    // save new model
     326    // XXX save the correct radius...
    316327    source->modelEXT = EXT;
    317     source->modelEXT->fitRadius = radius;
    318328    source->type = PM_SOURCE_TYPE_EXTENDED;
    319329    source->mode |= PM_SOURCE_MODE_EXTMODEL;
    … …  
    343353    source->modelPSF = psMemIncrRefCounter (DBL->data[0]);
    344354    source->mode     |= PM_SOURCE_MODE_PAIR;
    345     source->modelPSF->fitRadius = radius;
    346355
    347356    // copy most data from the primary source (modelEXT, blends stay NULL)
    348357    pmSource *newSrc = pmSourceCopyData (source);
    349358    newSrc->modelPSF = psMemIncrRefCounter (DBL->data[1]);
    350     newSrc->modelPSF->fitRadius = radius;
    351359
    352360    // build cached models and subtract
    … …  
    452460    maskVal |= markVal;
    453461
    454     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
    455 
    456462    // use the source moments, etc to guess basic model parameters
    457     pmModel *EXT = pmSourceModelGuess (source, modelType);
    458     if (!EXT) {
     463    pmModel *model = pmSourceModelGuess (source, modelType);
     464    if (!model) {
    459465        psTrace ("psphot", 5, "failed to generate a model for source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
    460466        return NULL;
    … …  
    463469    // for sersic models, use a grid search to choose an index, then float the params there
    464470    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
    465         psphotFitSersicIndex (source, EXT, fitOptions, maskVal);
     471        // for the test fits, use a somewhat smaller radius
     472        psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
     473        psphotFitSersicIndex (model, readout, source, fitOptions, maskVal, markVal);
     474    }
     475
     476    if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
     477        psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0);
    466478    }
    467479
    … …  
    471483        options.mode = PM_SOURCE_FIT_EXT;
    472484    }
    473     pmSourceFitModel (source, EXT, &options, maskVal);
     485
     486    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
     487    pmSourceFitModel (source, model, &options, maskVal);
     488    fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
    474489
    475490    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
    476     return (EXT);
     491    return (model);
    477492}
    478493
    479494pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
     495
     496    if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) {
     497        psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
     498    }
     499
     500    pmSourceFitOptions options = *fitOptions;
     501
     502    NfitPCM ++;
     503
     504    // maskVal is used to test for rejected pixels, and must include markVal
     505    maskVal |= markVal;
    480506
    481507    // allocate the model
    … …  
    484510        return NULL;
    485511    }
    486 
    487     pmSourceFitOptions options = *fitOptions;
    488 
    489     if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) {
    490         psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
    491     }
    492 
    493     float radius = psphotSetRadiusEXT (readout, source, markVal);
    494 
    495     // XXX note that this changes the source moments that are published...
    496     // recalculate the source moments using the larger extended-source moments radius
    497     // at this stage, skip Gaussian windowing, and do not clip pixels by S/N
    498     // this uses the footprint to judge both radius and aperture?
    499     if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) {
    500         // XXX set some mask bit/
    501         model->flags |= PM_MODEL_STATUS_BADARGS;
    502         return model;
    503     }
    504 
    505     NfitPCM ++;
    506 
    507     // maskVal is used to test for rejected pixels, and must include markVal
    508     maskVal |= markVal;
    509 
    510     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
    511512
    512513    pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize);
    … …  
    518519
    519520    // use the source moments, etc to guess basic model parameters
    520     pmSourceModelGuessPCM (pcm, source, maskVal, markVal);
     521    if (!pmSourceModelGuessPCM (pcm, source, maskVal, markVal)) {
     522        psFree (pcm);
     523        model->flags |= PM_MODEL_STATUS_BADARGS;
     524        return model;
     525    }
    521526
    522527    // for sersic models, use a grid search to choose an index, then float the params there
    523528    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
    524         psphotFitSersicIndexPCM (pcm, source, fitOptions, maskVal, markVal, psfSize);
     529        // for the test fits, use a somewhat smaller radius
     530        psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
     531        if (!psphotFitSersicIndexPCM (pcm, readout, source, fitOptions, maskVal, markVal, psfSize)) {
     532            model->flags |= PM_MODEL_STATUS_BADARGS;
     533            return model;
     534        }
     535    }
     536
     537    if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
     538        psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0);
    525539    }
    526540
    … …  
    530544        options.mode = PM_SOURCE_FIT_EXT;
    531545    }
     546
     547    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
    532548    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
     549    fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
    533550
    534551    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
    … …  
    544561// A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
    545562// for a sersic model, attempt to fit just the index and normalization with a modest number of iterations
    546 bool psphotFitSersicIndex (pmSource *source, pmModel *model, pmSourceFitOptions *fitOptions, psImageMaskType maskVal) {
     563bool psphotFitSersicIndex (pmModel *model, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
    547564
    548565    assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
    … …  
    552569    // fit EXT (not PSF) model (set/unset the pixel mask)
    553570    options.mode = PM_SOURCE_FIT_NO_INDEX;
    554     options.nIter = 3;
     571    options.nIter = 4;
    555572
    556573    int iMin = -1;
    … …  
    561578        model->params->data.F32[PM_PAR_7] = indexGuess[i];
    562579
    563         model->modelGuess(model, source);
     580        if (!model->modelGuess(model, source)) {
     581            model->flags |= PM_MODEL_STATUS_BADARGS;
     582            return false;
     583        }
     584
     585        // each time we change the model guess, we need to adjust the radius
     586        // XXX this did not work : we do not need such a large radius -- just uses moments-based radius
     587        if (false && !psphotSetRadiusModel (model, readout, source, markVal, false)) {
     588            psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
     589        }
     590       
    564591        pmSourceFitModel (source, model, &options, maskVal);
     592        fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
     593
     594        chiSquare[i] = model->chisqNorm;
     595        if (i == 0) {
     596            xMin = chiSquare[i];
     597            iMin = i;
     598        } else {
     599            if (chiSquare[i] < xMin) {
     600                xMin = chiSquare[i];
     601                iMin = i;
     602            }
     603        }
     604    }
     605    assert (iMin >= 0);
     606
     607    model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
     608    model->params->data.F32[PM_PAR_7] = indexGuess[iMin];
     609    model->modelGuess(model, source);
     610
     611    // each time we change the model guess, we need to adjust the radius
     612    // if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
     613    //  psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal);
     614    // }
     615
     616    return true;
     617}
     618
     619// A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
     620// for a sersic model, attempt to fit just the index and normalization with a modest number of iterations
     621bool psphotFitSersicIndexPCM (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
     622
     623    pmModel *model = pcm->modelConv;
     624
     625    assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
     626
     627    pmSourceFitOptions options = *fitOptions;
     628   
     629    // fit EXT (not PSF) model (set/unset the pixel mask)
     630    options.mode = PM_SOURCE_FIT_NO_INDEX;
     631    options.nIter = 4;
     632
     633    int iMin = -1;
     634    float xMin = NAN;
     635    float chiSquare[N_INDEX_GUESS];
     636
     637    for (int i = 0; i < N_INDEX_GUESS; i++) {
     638        model->params->data.F32[PM_PAR_7] = indexGuess[i];
     639       
     640        if (!model->modelGuess(model, source)) {
     641            model->flags |= PM_MODEL_STATUS_BADARGS;
     642            return false;
     643        }
     644
     645        // each time we change the model guess, we need to adjust the radius
     646        // XXX this did not work : we do not need such a large radius -- just uses moments-based radius
     647        if (false && !psphotSetRadiusModel (model, readout, source, markVal, false)) {
     648            psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
     649        }
     650       
     651        pmSourceFitModel (source, model, &options, maskVal);
     652        fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
     653
     654        // pmSourceModelGuessPCM(pcm, source, maskVal, markVal);
     655        // pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
    565656
    566657        chiSquare[i] = model->chisq;
    … …  
    575666        }
    576667    }
    577 
    578     model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
    579     model->params->data.F32[PM_PAR_7] = indexGuess[iMin];
    580     model->modelGuess(model, source);
    581 
    582     return true;
    583 }
    584 
    585 // A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
    586 // for a sersic model, attempt to fit just the index and normalization with a modest number of iterations
    587 bool psphotFitSersicIndexPCM (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
    588 
    589     pmModel *model = pcm->modelConv;
    590 
    591     assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
    592 
    593     pmSourceFitOptions options = *fitOptions;
    594    
    595     // fit EXT (not PSF) model (set/unset the pixel mask)
    596     options.mode = PM_SOURCE_FIT_NO_INDEX;
    597     options.nIter = 3;
    598 
    599     int iMin = -1;
    600     float xMin = NAN;
    601     float chiSquare[N_INDEX_GUESS];
    602 
    603     for (int i = 0; i < N_INDEX_GUESS; i++) {
    604         model->params->data.F32[PM_PAR_7] = indexGuess[i];
    605        
    606         model->modelGuess(model, source);
    607         pmSourceFitModel (source, model, &options, maskVal);
    608 
    609         // pmSourceModelGuessPCM(pcm, source, maskVal, markVal);
    610         // pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
    611 
    612         chiSquare[i] = model->chisq;
    613         if (i == 0) {
    614             xMin = chiSquare[i];
    615             iMin = i;
    616         } else {
    617             if (chiSquare[i] < xMin) {
    618                 xMin = chiSquare[i];
    619                 iMin = i;
    620             }
    621         }
    622     }
    623 
    624668    assert (iMin >= 0);
    625669   
Note: See TracChangeset for help on using the changeset viewer.