IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5772


Ignore:
Timestamp:
Dec 13, 2005, 6:43:44 AM (20 years ago)
Author:
eugene
Message:

substantial work on ensemble PSF fitting, better functions for evaluation, overall code cleanups

Location:
trunk/psphot
Files:
4 added
17 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/Makefile

    r5718 r5772  
    3434$(SRC)/psphotOutput.$(ARCH).o      \
    3535$(SRC)/psphotMarkPSF.$(ARCH).o     \
    36 $(SRC)/psphotSubtractPSF.$(ARCH).o \
    3736$(SRC)/psphotSortBySN.$(ARCH).o    \
    3837$(SRC)/psphotDefinePixels.$(ARCH).o\
    … …  
    4140$(SRC)/psphotBasicDeblend.$(ARCH).o  \
    4241$(SRC)/psphotModelGroupInit.$(ARCH).o  \
     42$(SRC)/psphotRadiusChecks.$(ARCH).o  \
     43$(SRC)/psphotReplaceUnfit.$(ARCH).o  \
     44$(SRC)/psphotEvalPSF.$(ARCH).o  \
     45$(SRC)/psphotEvalFLT.$(ARCH).o  \
     46$(SRC)/psphotFullFit.$(ARCH).o  \
    4347$(SRC)/pmSourceContour.$(ARCH).o \
    4448$(SRC)/psLine.$(ARCH).o            \

  • trunk/psphot/src/models/pmModel_SGAUSS.c

    r5593 r5772  
    8181    params_min[0][0].data.F32[6] = -5.0;
    8282    params_min[0][0].data.F32[7] = 0.5;
    83     params_min[0][0].data.F32[8] = 0.001;
     83    params_min[0][0].data.F32[8] = 0.05;
    8484
    8585    params_max[0][0].data.F32[0] = 1e5;

  • trunk/psphot/src/pmPeaksSigmaLimit.c

    r5672 r5772  
    2424    NSIGMA    = psMetadataLookupF32 (&status, config, "PEAKS_NSIGMA_LIMIT");
    2525   
    26     // threshold = NSIGMA*sky->sampleStdev + sky->sampleMean;
    27     threshold = NSIGMA*sky->sampleStdev;
     26    threshold = NSIGMA*sky->sampleStdev + sky->sampleMean;
     27    // threshold = NSIGMA*sky->sampleStdev;
    2828    psLogMsg ("psphot", 3, "threshold: %f DN\n", threshold);
    2929

  • trunk/psphot/src/psModulesUtils.c

    r5350 r5772  
    101101    return true;
    102102}
     103
     104pmModel *pmModelCopy (pmModel *model) {
     105
     106    pmModel *new = pmModelAlloc (model->type);
     107   
     108    new->chisq = model->chisq;
     109    new->nIter = model->nIter;
     110
     111    for (int i = 0; i < new->params->n; i++) {
     112        new->params->data.F32[i]  = model->params->data.F32[i];
     113        new->dparams->data.F32[i] = model->dparams->data.F32[i];
     114    }
     115   
     116    return (new);
     117}

  • trunk/psphot/src/psphot.c

    r5756 r5772  
    4343    psfClump = pmSourcePSFClump (sources, config);
    4444
    45     // group into STAR, COSMIC, GALAXY, SATURATED
     45    // group into STAR, COSMIC, GALAXY, SATURATED, etc.
    4646    pmSourceRoughClass (sources, config, psfClump);
     47
     48    // optional dump of all rough source data
     49    char *output = psMetadataLookupPtr (&status, config, "MOMENTS_OUTPUT_FILE");
     50    if (status && (output != NULL) && (output[0]))
     51    {
     52      pmMomentsWriteText (sources, output);
     53      psFree (output);
     54    }
    4755    if (!strcasecmp (breakPt, "CLASS")) exit (0);
    4856
     57    // mark blended peaks PS_SOURCE_BLEND
    4958    psphotBasicDeblend (sources, config, sky);
    5059    if (!strcasecmp (breakPt, "DEBLEND")) exit (0);
    … …  
    5766    if (!strcasecmp (breakPt, "PSFFIT")) exit (0);
    5867
     68    // XXX add this as conditional output
     69    psphotSamplePSFs (psf, imdata->image);
     70
    5971    int FITMODE = psMetadataLookupS32 (&status, config, "FIT_MODE");
    6072    if (!status) FITMODE = 2;
    6173
    6274    switch (FITMODE) {
    63       case -2:
     75      case 0:
    6476        psphotEnsemblePSF (imdata, config, sources, psf, sky);
    65         psphotReapplyPSF (imdata, config, sources, psf, sky);
    66         break;
    67 
    68       case -1:
    69         psphotEnsemblePSF (imdata, config, sources, psf, sky);
    70         break;
    71 
    72       case 0:
    73         psphotFixedPSF (imdata, config, sources, psf, sky);
     77        sources = psArraySort (sources, psphotSortBySN);
     78        psphotFullFit (imdata, config, sources, psf, sky);
     79//      psphotReplaceUnfit (sources);
    7480        break;
    7581
    7682      case 1:
    77         // test PSF on all except SATURATE and DEFECT (artifacts)
     83        // fit extended objects with galaxy models
    7884        psphotApplyPSF (imdata, config, sources, psf, sky);
     85        psphotFitGalaxies (imdata, config, sources, sky);
    7986        break;
    8087
    … …  
    8289        // fit extended objects with galaxy models
    8390        psphotApplyPSF (imdata, config, sources, psf, sky);
    84         psphotFitGalaxies (imdata, config, sources, sky);
    8591        break;
    8692

  • trunk/psphot/src/psphot.h

    r5718 r5772  
    7474psMetadata     *psphotTestArguments (int *argc, char **argv);
    7575bool            psphotBasicDeblend (psArray *sources, psMetadata *config, psStats *sky);
     76
     77bool psphotFullFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
     78bool psphotInitLimitsPSF (psMetadata *config);
     79bool psphotEvalPSF (pmSource *source);
     80bool psphotEvalFLT (pmSource *source);
     81bool psphotInitRadiusPSF (psMetadata *config, psStats *sky, pmModelType type);
     82bool psphotCheckRadiusPSF (eamReadout *imdata, pmSource *source);
     83bool psphotInitRadiusFLT (psMetadata *config, psStats *sky, pmModelType type);
     84bool psphotCheckRadiusFLT (eamReadout *imdata, pmSource *source);
     85bool psphotSamplePSFs (pmPSF *psf, psImage *image);
     86bool psphotReplaceUnfit (psArray *sources);

  • trunk/psphot/src/psphotApplyPSF.c

    r5131 r5772  
    22
    33// fit psf model to all objects
    4 // PSFSTAR objects will be refitted
    5 // run this function to a specific flux limit?
    6 
    74bool psphotApplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky)
    85{
    9     bool  status;
    106    float x;
    117    float y;
    12     int   Nfit = 0;
    13     int   Nsub = 0;
     8    int   Nfit  = 0;
     9    int   Nsub  = 0;
    1410    int   Niter = 0;
    1511
    1612    psTimerStart ("psphot");
    1713
    18     // we may set this differently here from the value used to mark likely saturated stars
    19     float SATURATION   = psMetadataLookupF32 (&status, config, "SATURATION");
    20     float OUTER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
    21 
    22     float PSF_MIN_SN       = psMetadataLookupF32 (&status, config, "PSF_MIN_SN");
    23     float PSF_MAX_CHI      = psMetadataLookupF32 (&status, config, "PSF_MAX_CHI");
    24     float PSF_FIT_NSIGMA   = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA");
    25     float PSF_FIT_PADDING  = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING");
    26     float PSF_SHAPE_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA");
    27 
    28     // set the object surface-brightness limit for fitted pixels
    29     float FLUX_LIMIT  = PSF_FIT_NSIGMA * sky->sampleStdev;
    30     psLogMsg ("psphot.apply_psf_model", 4, "fitting pixels with at least %f object counts\n", FLUX_LIMIT);
    31 
    32     // this function specifies the radius at this the model hits the given flux
    33     pmModelRadius modelRadius = pmModelRadius_GetFunction (psf->type);
     14    psphotInitLimitsPSF (config);
     15    psphotInitRadiusPSF (config, sky, psf->type);
    3416
    3517    for (int i = 0; i < sources->n; i++) {
    … …  
    3820
    3921        // skip non-astronomical objects (very likely defects)
    40         // XXX EAM : should we try these anyway?
     22        if (source->mode  & PM_SOURCE_BLEND) continue;
    4123        if (source->type == PM_SOURCE_DEFECT) continue;
    4224        if (source->type == PM_SOURCE_SATURATED) continue;
    … …  
    4729        // set PSF parameters for this model
    4830        pmModel *model = pmModelFromPSF (modelFLT, psf);
     31        psFree (modelFLT);
     32
     33        source->modelPSF = model;
     34
     35        // sets the model radius (via source->model) and adjusts pixels as needed
     36        psphotCheckRadiusPSF (imdata, source);
     37
    4938        x = model->params->data.F32[2];
    5039        y = model->params->data.F32[3];
    51         psFree (modelFLT);
    5240
    53         // set the fit radius based on the object flux limit and the model
    54         // XXX EAM : FLUX_LIMIT should be set based on local sky model (not global median)
    55         model->radius = modelRadius (model->params, FLUX_LIMIT) + PSF_FIT_PADDING;
    56         if (isnan(model->radius)) {
    57           psAbort ("apply_psf_model", "error in radius");
    58         }
    59        
    60         // check if we need to redefine the pixels
    61         // XXX EAM : a better test would examine the source pixels
    62         if (model->radius > OUTER_RADIUS) {
    63           // (re)-allocate image, weight, mask arrays for each peak (square of radius OUTER)
    64             psphotDefinePixels (source, imdata, x, y, model->radius);
    65         }
    66 
    67         // if (i > 66) psTraceSetLevel (".psLib.dataManip.psMinimizeLMChi2", 5);
    68            
    6941        // fit PSF model (set/unset the pixel mask)
    7042        psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
    71         status = pmSourceFitModel (source, model, true);
     43        pmSourceFitModel (source, model, true);
    7244        psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED);
    7345
    74         if (!status || (model->params->data.F32[1] < 0)) {
    75           psLogMsg ("psphot", 5, "PSF fit failed for %f, %f (%d iterations)\n", x, y, model->nIter);
    76           source->type = PM_SOURCE_FAIL_FIT_PSF;  // better choice?
    77           psFree (model);
    78           continue;
    79         }
    80         // XXX EAM : this was an attempt to correct for fit-sky biases
    81         // pmModelSkyOffset (model, x, y, model->radius);
    82 
    83         // model succeeded : tentatively keep it
    84         source->modelPSF = model;
    8546        Niter += model[0].nIter;
    8647        Nfit ++;
    8748
    88         // is it a good model?
    89         psphotMarkPSF (source, PSF_SHAPE_NSIGMA, PSF_MIN_SN, PSF_MAX_CHI, SATURATION);
    90         if (psphotSubtractPSF (source)) {
    91           Nsub ++;
     49        // check if model fit is acceptable
     50        if (psphotEvalPSF (source)) {
     51            pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false);
     52            source->mode |= PM_SOURCE_SUBTRACTED;
     53            Nsub ++;
    9254        }
    9355    }
    9456
    9557    psLogMsg ("psphot", 3, "fit PSF models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter);
    96     psLogMsg ("psphot", 4, "subtracted %d PSF objects\n", Nsub);
     58    psLogMsg ("psphot", 4, "subtracted %d PSF models\n", Nsub);
    9759    return (true);
    9860}

  • trunk/psphot/src/psphotBasicDeblend.c

    r5718 r5772  
    1414    psTimerStart ("psphot");
    1515
    16     // this should be added to the PM_SOURCE flags:
    17     int PM_SOURCE_BLEND = PM_SOURCE_OTHER + 1;
     16    float FRACTION = psMetadataLookupF32 (&status, config, "DEBLEND_PEAK_FRACTION");
     17    if (!status) FRACTION = 0.25;
    1818
    19     float FRACTION = psMetadataLookupF32 (&status, config, "DEBLEND_PEAK_FRACTION");
    2019    float NSIGMA   = psMetadataLookupF32 (&status, config, "DEBLEND_SKY_NSIGMA");
     20    if (!status) NSIGMA = 5.0;
     21
    2122    float minThreshold = NSIGMA*sky->sampleStdev;
    22     fprintf (stderr, "min threshold: %f\n", minThreshold);
    2323
    2424    // we need sources spatially-sorted to find overlaps
    … …  
    4343        source = sources->data[N];
    4444
    45         if (source->type == PM_SOURCE_BLEND) continue;
     45        if (source->mode & PM_SOURCE_BLEND) continue;
    4646
    4747        overlap->n = 0;
    … …  
    119119                             );
    120120
    121                     testSource->type = PM_SOURCE_BLEND;
     121                    testSource->mode |= PM_SOURCE_BLEND;
    122122                    Nblend ++;
    123123                    j = xv->n;

  • trunk/psphot/src/psphotChoosePSF.c

    r5058 r5772  
    1313    // array to store candidate PSF stars
    1414    int NSTARS = psMetadataLookupS32 (&status, config, "PSF_MAX_NSTARS");
    15     if (!status) NSTARS = sources->n;
     15    if (!status) NSTARS = PS_MIN (sources->n, 200);
     16
    1617    stars = psArrayAlloc (sources->n);
    1718    stars->n = 0;
    … …  
    2021    for (int i = 0; (i < sources->n) && (stars->n < NSTARS); i++) {
    2122        pmSource *source = sources->data[i];
    22         if (source->type != PM_SOURCE_PSFSTAR) continue;
    23         psArrayAdd (stars, 200, source);
     23        if (source->mode & PM_SOURCE_PSFSTAR) psArrayAdd (stars, 200, source);
    2424    }
    25     psTrace (".psphot.pspsf", 3, "selected candidate %d PSF objects\n", stars->n);
     25    psLogMsg ("psphot.pspsf", 3, "selected candidate %d PSF objects\n", stars->n);
    2626
    2727    // get the fixed PSF fit radius
    … …  
    6666        }
    6767    }
     68    // XXX EAM : need to unflag the PSF stars which are not used to build the PSF
     69    // XXX EAM : each pmPSFtry needs to have its own mask array
    6870
    6971    // keep only the selected model:

  • trunk/psphot/src/psphotEnsemblePSF.c

    r5718 r5772  
    1212    // source analysis is done in spatial order
    1313    sources = psArraySort (sources, psphotSortByY);
    14 
    15     // this should be added to the PM_SOURCE flags:
    16     int PM_SOURCE_BLEND = PM_SOURCE_OTHER + 1;
    1714
    1815    float OUTER_RADIUS     = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
    … …  
    3532        // skip non-astronomical objects (very likely defects)
    3633        // XXX EAM : should we try these anyway?
    37         if (inSource->type == PM_SOURCE_BLEND) continue;
     34        if (inSource->mode & PM_SOURCE_BLEND) continue;
    3835        if (inSource->type == PM_SOURCE_DEFECT) continue;
    3936        if (inSource->type == PM_SOURCE_SATURATED) continue;
    … …  
    7370        psImage *mask = otSource->mask;
    7471
    75         // XXX EAM : set model to unit peak, zero sky, maybe use peak (x,y)
     72        // XXX EAM : use these lines to fit to the peak
    7673        // model->params->data.F32[2] = inSource->peak->x;
    7774        // model->params->data.F32[3] = inSource->peak->y;
     75        // XXX EAM : better option: improve the peak with 2D poly fit 3x3
     76
     77        // set model to unit peak, zero sky (we assume sky is constant)
    7878        model->params->data.F32[0] = 0.0;
    7979        model->params->data.F32[1] = 1.0;

  • trunk/psphot/src/psphotFitGalaxies.c

    r5672 r5772  
    33bool psphotFitGalaxies (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats)
    44{
    5     bool  status, goodfit;
     5    bool  status;
    66    float x;
    77    float y;
    8     int   Nfit = 0;
    9     int   Nfail = 0;
     8    int   Nfit  = 0;
     9    int   Nsub = 0;
    1010    int   Niter = 0;
    1111
    12     float OUTER_RADIUS     = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
     12    psTimerStart ("psphot");
    1313
    1414    float GAL_MIN_SN       = psMetadataLookupF32 (&status, config, "GAL_MIN_SN");
    15     float GAL_FIT_NSIGMA   = psMetadataLookupF32 (&status, config, "GAL_FIT_NSIGMA");
    16     float GAL_FIT_PADDING  = psMetadataLookupF32 (&status, config, "GAL_FIT_PADDING");
    1715    float GAL_MOMENTS_RAD  = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS");
     16    char       *modelName  = psMetadataLookupPtr (&status, config, "GAL_MODEL");
     17    pmModelType modelType  = pmModelSetType (modelName);
    1818
    19     float FLUX_LIMIT       = GAL_FIT_NSIGMA * skyStats->sampleStdev;
    20 
    21     char         *modelName   = psMetadataLookupPtr (&status, config, "GAL_MODEL");
    22     pmModelType   modelType   = pmModelSetType (modelName);
    23     pmModelRadius modelRadius = pmModelRadius_GetFunction (modelType);
    24 
    25     psTimerStart ("psphot");
     19    psphotInitRadiusFLT (config, skyStats, modelType);
    2620
    2721    for (int i = 0; i < sources->n; i++) {
    2822        pmSource *source = sources->data[i];
    2923
    30         // sources which should not be fitted
    31         // skip all valid stars
    32         if (source->type == PM_SOURCE_PSFSTAR) continue;
    33         if (source->type == PM_SOURCE_SATSTAR) continue;
    34         if (source->type == PM_SOURCE_GOODSTAR) continue;
     24        // only fit suspected extended sources
     25        if (source->type != PM_SOURCE_GALAXY) continue;
     26        if (source->mode  & PM_SOURCE_BLEND) continue;
    3527
    36         // skip all likely defects
    37         if (source->type == PM_SOURCE_DEFECT) continue;
    38         if (source->type == PM_SOURCE_SATURATED) continue;
    39 
    40         // skip poorly fitted stars
    41         if (source->type == PM_SOURCE_FAINTSTAR) continue;
    42         if (source->type == PM_SOURCE_POOR_FIT_PSF) continue;
    43 
    44         // XXX EAM when do we pick these up again?
    45         if (source->moments->SN < GAL_MIN_SN) {
    46           source->type = PM_SOURCE_FAINT_GALAXY;  // better choice?
    47           continue;
    48         }
     28        // skip possible extended sources below fit threshold
     29        if (source->moments->SN < GAL_MIN_SN) continue;
    4930
    5031        // recalculate the source moments using the larger galaxy moments radius
    5132        status = pmSourceMoments (source, GAL_MOMENTS_RAD);
    5233        if (!status) {
    53           source->type = PM_SOURCE_DROP_GALAXY;  // better choice?
     34          source->mode |= PM_SOURCE_FAIL;  // better choice?
    5435          continue;
    5536        }
    … …  
    5738        // use the source moments, etc to guess basic model parameters
    5839        pmModel  *model  = pmSourceModelGuess (source, modelType);
     40        source->modelFLT = model;
    5941        x = model->params->data.F32[2];
    6042        y = model->params->data.F32[3];
    6143
    62         // set the fit radius based on the object flux limit and the model
    63         // XXX EAM : FLUX_LIMIT should be set based on local sky model (not global median)
    64         model->radius = modelRadius (model->params, FLUX_LIMIT) + GAL_FIT_PADDING;
    65         if (isnan(model->radius)) psAbort ("fit_galaxies", "error in radius");
    66 
    67         // check if we need to redefine the pixels
    68         // XXX EAM : a better test would examine the source pixels
    69         if (model->radius > OUTER_RADIUS) {
    70           // (re)-allocate image, weight, mask arrays for each peak (square of radius OUTER)
    71           psphotDefinePixels (source, imdata, x, y, model->radius);
    72         }
     44        // sets the model radius (via source->model) and adjusts pixels as needed
     45        psphotCheckRadiusFLT (imdata, source);
    7346
    7447        // fit FLT (not PSF) model (set/unset the pixel mask)
    7548        psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
    76         status = pmSourceFitModel (source, model, false);
     49        pmSourceFitModel (source, model, false);
    7750        psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED);
    78         if (!status) {
    79           // if the fit fails, we need to change the classification
    80           psLogMsg ("psphot", 5, "GAL fit failed for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius);
    81           source->type = PM_SOURCE_FAIL_FIT_GAL;  // better choice?
    82           source->modelFLT = model;
    83           Nfail ++;
    84           continue;
    85         }
    8651
    87         // check if model fit is acceptable
    88         goodfit = pmModelFitStatus (model);
    89         if (!goodfit) {
    90           // if the fit fails, we need to change the classification
    91           psLogMsg ("psphot", 5, "GAL fit poor for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius);
    92           source->type = PM_SOURCE_POOR_FIT_GAL;  // better choice?
    93           source->modelFLT = model;
    94           Nfail ++;
    95           continue;
    96         }
    97 
    98         // accept the model
    99         source->type = PM_SOURCE_GALAXY;
    100         source->modelFLT = model;
    10152        Niter += model[0].nIter;
    10253        Nfit++;
    10354
    104         // subtract object, leave local sky
    105         pmSourceSubModel (source->pixels, source->mask, model, false, false);
     55        // check if model fit is acceptable
     56        if (pmModelFitStatus (model)) {
     57            pmSourceSubModel (source->pixels, source->mask, model, false, false);
     58            source->mode |= PM_SOURCE_SUBTRACTED;
     59            Nsub ++;
     60        }
    10661    }
    107     psLogMsg ("psphot", 3, "fit galaxies: %f sec for %d galaxies (%d failures, %d total iterations)\n", psTimerMark ("psphot"), Nfit, Nfail, Niter);
     62    psLogMsg ("psphot", 3, "fit GAL models: %f sec for %d galaxies (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter);
     63    psLogMsg ("psphot", 4, "subtracted %d GAL models\n", Nsub);
    10864    return (true);
    10965}

  • trunk/psphot/src/psphotFixedPSF.c

    r5672 r5772  
    6666        if (!status || (model->params->data.F32[1] < 0)) {
    6767          psLogMsg ("psphot", 5, "PSF fit failed for %f, %f (peak %f\n", x, y, model->params->data.F32[1]);
    68           source->type = PM_SOURCE_FAIL_FIT_PSF;  // better choice?
     68          source->mode |= PM_SOURCE_FAIL;  // better choice?
    6969          psFree (model);
    7070          continue;

  • trunk/psphot/src/psphotImageBackground.c

    r5718 r5772  
    6464        }
    6565    }
     66    sky->sampleMean = 0.0;
    6667
    6768    psLogMsg ("psphot", 5, "back: %f sec (fit model)\n", psTimerMark ("psphot"));
    6869
    6970    psphotSaveImage (imdata->header, imdata->image, "backsub.fits");
    70 
    7171    return (skyModel);
    7272}

  • trunk/psphot/src/psphotMagnitudes.c

    r5672 r5772  
    4949   
    5050        // use PSF model with stars
    51       case PM_SOURCE_PSFSTAR:
    52       case PM_SOURCE_SATSTAR:
    53       case PM_SOURCE_GOODSTAR:
     51      case PM_SOURCE_STAR:
    5452        model = source->modelPSF;
    5553        break;
    … …  
    6159           
    6260        // skip defects and poorly fitted stars or galaxies
    63       case PM_SOURCE_DEFECT:
    64       case PM_SOURCE_SATURATED:
    65       case PM_SOURCE_FAINTSTAR:
    66       case PM_SOURCE_POOR_FIT_PSF:
    67       case PM_SOURCE_POOR_FIT_GAL:
    68       case PM_SOURCE_FAIL_FIT_GAL:
    6961      default:
    7062        model = NULL;

  • trunk/psphot/src/psphotMarkPSF.c

    r5672 r5772  
    11# include "psphot.h"
     2
     3// given a pmSource which has been fitted using modelPSF, evaluate the
     4// resulting fit: did the fit succeed? is this object PSF-like? is this object
     5// extended?  return status is TRUE for a valid PSF, FALSE otherwise.
    26
    37// identify objects consistent with PSF shape/magnitude distribution
    … …  
    1418// this includes a minimum buffer (DS) for the brighter objects
    1519
    16 // saturated stars should fall outside (but are already IDed)
    17 // galaxies should be larger, cosmic rays smaller, but need to test?
     20// saturated stars should fall outside (larger), but have peaks above SATURATE
     21// galaxies should be larger, cosmic rays smaller
    1822// we also reject objects with S/N too low or ChiSquare to high
    19 
    20 // any object which meets the criterion is marked as PM_SOURCE_BRIGHTSTAR
    2123
    2224// floor for DS value
    … …  
    3335
    3436    // do we actually have a valid PSF model?
    35     if (model == NULL) return (false);
     37    if (model == NULL) {
     38        source->mode &= ~PM_SOURCE_FITTED;
     39        return false;
     40    }
    3641
    3742    // did the model fit fail for one or another reason?
    3843    switch (model->status) {
     44      case PM_MODEL_UNTRIED:
     45        source->mode &= ~PM_SOURCE_FITTED;
     46        return false;
    3947      case PM_MODEL_BADARGS:
    40       case PM_MODEL_UNTRIED:
    41         source->type = PM_SOURCE_FAIL_FIT_PSF;  // better choice?
     48        source->mode |= PM_SOURCE_FAIL;
    4249        return false;
    4350      case PM_MODEL_SUCCESS:
    … …  
    4653      case PM_MODEL_OFFIMAGE:
    4754        psLogMsg ("psphot", 5, "PSF fit failed for %f, %f (%d iterations)\n", model->params->data.F32[2], model->params->data.F32[3], model->nIter);
    48         source->type = PM_SOURCE_FAIL_FIT_PSF;  // better choice?
     55        source->mode |= PM_SOURCE_FAIL;
    4956        return false;
    5057      default:
    … …  
    5259    }
    5360
     61    // unless we prove otherwise, this object is a star.
     62    source->type = PM_SOURCE_STAR;
     63
    5464    // if the object has fitted peak above saturation, label as SATSTAR
    5565    // this is a valid PSF object, but ignore the other quality tests
    5666    // remember: fit does not use saturated pixels (masked)
     67    // XXX no extended object can saturate and stay extended...
    5768    if (model->params->data.F32[1] >= SATURATE) {
    58         if (source->type == PM_SOURCE_PSFSTAR) {
     69        if (source->mode & PM_SOURCE_PSFSTAR) {
    5970            psLogMsg ("psphot", 5, "PSFSTAR marked SATSTAR\n");
    6071        }
    61         source->type = PM_SOURCE_SATSTAR;
    62         return (true);
     72        source->mode &= ~PM_SOURCE_PSFSTAR;
     73        source->mode |=  PM_SOURCE_SATSTAR;
     74        return true;
    6375    }
    6476
    6577    // if the object has a fitted peak below 0, the fit did not converge cleanly
    6678    if (model->params->data.F32[1] < 0) {
    67         source->type = PM_SOURCE_FAIL_FIT_PSF;
    68         return (false);
     79        source->mode |= PM_SOURCE_FAIL;
     80        return false;
    6981    }
    7082
    7183    // if the source was predicted to be a SATSTAR, but it fitted below saturation,
    7284    // make a note to the user
    73     if (source->type == PM_SOURCE_SATSTAR) {
    74         psLogMsg ("psphot", 5, "SATSTAR marked GOODSTAR (fitted peak below saturation)\n");
    75         source->type = PM_SOURCE_GOODSTAR;
     85    if (source->mode & PM_SOURCE_SATSTAR) {
     86        psLogMsg ("psphot", 5, "SATSTAR marked normal (fitted peak below saturation)\n");
     87        source->mode &= ~PM_SOURCE_SATSTAR;
    7688    }
    7789
    … …  
    93105    keep &= (SN > minSN);
    94106    keep &= (Chi < maxChi);
    95     if (keep) {
    96         if (source->type == PM_SOURCE_PSFSTAR) return (true);
    97         source->type = PM_SOURCE_GOODSTAR;
    98         return (true);
    99     }
    100    
    101     if (source->type == PM_SOURCE_PSFSTAR) {
     107    if (keep) return true;
     108
     109    // this source is not a star, unflag as PSFSTAR
     110    // XXX : if this object was used to build the PSF, this flag should
     111    //       be set even if the object is not a star...
     112    if (source->mode & PM_SOURCE_PSFSTAR) {
     113        source->mode &= ~PM_SOURCE_PSFSTAR;
    102114        psLogMsg ("psphot", 5, "PSFSTAR demoted based on fit quality\n");
    103115    }
    … …  
    106118    if ((nSx <= -shapeNsigma) || (nSy <= -shapeNsigma)) {
    107119        source->type = PM_SOURCE_DEFECT;
    108         return (false);
     120        return false;
    109121    }
    110122
    … …  
    112124    if ((nSx >= shapeNsigma) || (nSy >= shapeNsigma)) {
    113125        source->type = PM_SOURCE_GALAXY;
    114         return (false);
     126        return false;
    115127    }
    116128
    117     // object appears to be extremely faint: what is this?
    118     if (SN < minSN) {
    119         source->type = PM_SOURCE_FAINTSTAR;
    120         return (false);
    121     }
    122 
    123     // these are pooly fitted, probable stars near other stars?
    124     source->type = PM_SOURCE_POOR_FIT_PSF;
    125     return (false);
     129    // poor-quality fit; only keep if nothing else works...
     130    source->mode |= PM_SOURCE_POOR;
     131    return false;
    126132}       

  • trunk/psphot/src/psphotOutput.c

    r5754 r5772  
    389389        pmSource *source = (pmSource *) sources->data[i];
    390390
    391         // valid source types for this function
    392         if (source->type == PM_SOURCE_SATSTAR) goto valid;
    393         if (source->type == PM_SOURCE_PSFSTAR) goto valid;
    394         if (source->type == PM_SOURCE_GOODSTAR) goto valid;
    395         continue;
    396 
    397     valid:
     391        // write out sources fitted as PSFs
     392        if (source->type != PM_SOURCE_STAR) continue;
     393        if (source->mode & PM_SOURCE_FAIL) continue;
     394        if (source->mode & PM_SOURCE_POOR) continue;
     395
    398396        model = pmSourceMagnitudes (source, psf, RADIUS);
    399397        if (model == NULL) continue;
    … …  
    452450        model = (pmModel  *) source->modelFLT;
    453451        if (model == NULL) continue;
    454         if (source->type == PM_SOURCE_GALAXY) goto valid;
    455         continue;
    456 
    457     valid:
     452        if (source->type != PM_SOURCE_GALAXY) continue;
     453        if (source->mode & PM_SOURCE_FAIL) continue;
     454        if (source->mode & PM_SOURCE_POOR) continue;
     455       
    458456        params = model->params;
    459457        dparams = model->dparams;
    … …  
    513511
    514512        // skip these sources (in PSF or FLT)
    515         if (source->type == PM_SOURCE_GALAXY) continue;
    516         if (source->type == PM_SOURCE_PSFSTAR) continue;
    517         if (source->type == PM_SOURCE_SATSTAR) continue;
    518         if (source->type == PM_SOURCE_GOODSTAR) continue;
    519 
     513        if (source->type == PM_SOURCE_STAR) {
     514            if (source->mode & PM_SOURCE_FAIL) goto isNULL;
     515            if (source->mode & PM_SOURCE_POOR) goto isNULL;
     516            continue;
     517        }
     518        if (source->type == PM_SOURCE_GALAXY) {
     519            if (source->mode & PM_SOURCE_FAIL) goto isNULL;
     520            if (source->mode & PM_SOURCE_POOR) goto isNULL;
     521            continue;
     522        }
     523           
     524    isNULL:
     525       
    520526        if (source->moments == NULL) {
    521527          fprintf (f, "%5d %5d  %7.1f  %7.1f %7.1f  %6.3f %6.3f  %8.1f %7.1f %7.1f %7.1f  %4d %2d\n",
    … …  
    579585        return (8);
    580586
    581       case PM_SOURCE_SATSTAR:
    582         return (10);
    583 
    584       case PM_SOURCE_PSFSTAR:
    585       case PM_SOURCE_GOODSTAR:
     587      case PM_SOURCE_STAR:
     588        if (source->mode & PM_SOURCE_SATSTAR) return (10);
     589        if (source->mode & PM_SOURCE_POOR) return (7);
     590        if (source->mode & PM_SOURCE_FAIL) return (4);
    586591        return (1);
    587592
    588       case PM_SOURCE_POOR_FIT_PSF:
    589         return (7);
    590 
    591       case PM_SOURCE_FAIL_FIT_PSF:
    592       case PM_SOURCE_FAINTSTAR:
    593         return (4);
    594 
    595593      case PM_SOURCE_GALAXY:
    596       case PM_SOURCE_FAINT_GALAXY:
    597       case PM_SOURCE_DROP_GALAXY:
    598       case PM_SOURCE_FAIL_FIT_GAL:
    599       case PM_SOURCE_POOR_FIT_GAL:
    600594        return (2);
    601595
    602       case PM_SOURCE_OTHER:
    603596      default:
    604597        return (0);

  • trunk/psphot/src/psphotReapplyPSF.c

    r5718 r5772  
    1414    sources = psArraySort (sources, psphotSortBySN);
    1515   
    16      // this should be added to the PM_SOURCE flags:
    17     int PM_SOURCE_BLEND = PM_SOURCE_OTHER + 1;
    18 
    1916   // we may set this differently here from the value used to mark likely saturated stars
    2017    float SATURATION       = psMetadataLookupF32 (&status, config, "SATURATION");
    … …  
    2926        // skip non-astronomical objects (very likely defects)
    3027        // XXX EAM : should we try these anyway?
    31         if (source->type == PM_SOURCE_BLEND) continue;
     28        if (source->mode |= PM_SOURCE_BLEND) continue;
    3229        if (source->type == PM_SOURCE_DEFECT) continue;
    3330        if (source->type == PM_SOURCE_SATURATED) continue;
    … …  
    7067    return (true);
    7168}
    72 
    73 pmModel *pmModelCopy (pmModel *model) {
    74 
    75     pmModel *new = pmModelAlloc (model->type);
    76    
    77     new->chisq = model->chisq;
    78     new->nIter = model->nIter;
    79 
    80     for (int i = 0; i < new->params->n; i++) {
    81         new->params->data.F32[i]  = model->params->data.F32[i];
    82         new->dparams->data.F32[i] = model->dparams->data.F32[i];
    83     }
    84    
    85     return (new);
    86 }
Note: See TracChangeset for help on using the changeset viewer.