IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 4574


Ignore:
Timestamp:
Jul 16, 2005, 6:02:29 AM (21 years ago)
Author:
eugene
Message:

improved galaxy model guesses, better handling of errors

Location:
trunk/psphot
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/Makefile

    r4375 r4574  
    3939$(SRC)/subtract_galaxies.$(ARCH).o \
    4040$(SRC)/LocalSky.$(ARCH).o \
     41$(SRC)/basic_classes.$(ARCH).o \
    4142$(SRC)/find_defects.$(ARCH).o
    4243

  • trunk/psphot/src/apply_psf_model.c

    r4421 r4574  
    6565        status = pmSourceFitModel (source, model, true);
    6666        psImageKeepCircle (source->mask, x, y, model->radius, AND, 0x7f);
    67         if (!status) continue;
     67        if (!status || (model->params->data.F32[1] < 0)) {
     68          // if the fit fails, we need to change the classification
     69          psLogMsg ("psphot", 3, "PSF fit failed for %f, %f (%d iterations)\n", x, y, model->nIter);
     70          source->type = PS_SOURCE_OTHER;  // better choice?
     71          continue;
     72        }
    6873
    6974        source->modelPSF = model;
    … …  
    7176        Nfit ++;
    7277
    73         mark_psf_source (source, shapeNsigma);
     78        mark_psf_source (source, shapeNsigma, SATURATE);
    7479        if (subtract_psf_source (source)) {
    7580          Nsub ++;

  • trunk/psphot/src/basic_classes.c

    r4398 r4574  
    22
    33bool basic_classes (psArray *sources, psMetadata *config) {
     4
     5    pmPSFClump psfClump;
    46
    57    // group into STAR, COSMIC, GALAXY, SATURATED

  • trunk/psphot/src/fit_galaxies.c

    r4251 r4574  
    33// fit selected galaxy model (GAUSS) to all bright objects of type GALAXY
    44
    5 bool fit_galaxies (psImage *image, psMetadata *config, psArray *sources)
     5bool fit_galaxies (psImageData *imdata, psMetadata *config, psArray *sources, psStats *skyStats)
    66{
    77    bool  status;
     8    float x;
     9    float y;
     10    float sky;
    811    int   Nfit = 0;
     12    int   Nfail = 0;
    913    int   Niter = 0;
    1014
     15    float MOMENT_R = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS");
    1116    float RADIUS   = psMetadataLookupF32 (&status, config, "FIT_RADIUS");
    1217    float snFaint  = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM");
    13     psModelType modelType = psModelSetType ("PS_MODEL_PGAUSS");
     18    float OUTER    = psMetadataLookupF32 (&status, config, "OUTER_RADIUS");
     19    float FIT_NSIGMA  = psMetadataLookupF32 (&status, config, "FIT_NSIGMA");
     20    float FIT_PADDING = psMetadataLookupF32 (&status, config, "FIT_PADDING");
     21
     22    float FLUX_LIMIT  = FIT_NSIGMA * skyStats->sampleStdev;
     23
     24    psModelType   modelType   = psModelSetType ("PS_MODEL_RGAUSS");
     25    psModelRadius modelRadius = psModelRadius_GetFunction (modelType);
     26
     27    psTraceSetLevel (".psModules.pmSourceMoments", 5);
    1428
    1529    psTimerStart ("psphot");
    … …  
    1933        if (source->moments->SN < snFaint) continue;
    2034
     35        // recalculate the source moments using the galaxy radius (larger)
     36        status = pmSourceMoments (source, MOMENT_R);
     37
     38        // use the source moments, etc to guess basic model parameters
     39        psModel  *model  = pmSourceModelGuess (source, modelType);
     40
     41        x = model->params->data.F32[2];
     42        y = model->params->data.F32[3];
     43
    2144        // need a better model guess and a better radius choice
    2245        // when radius is not fixed, we will need to check if new radius fits on image
    23 # if (0)
     46
     47        // set the fit radius based on the object flux limit and the model
     48        // FLUX_LIMIT should be set based on local sky model (not global median)
     49        // model->radius = 25.0;
     50        model->radius = modelRadius (model->params, FLUX_LIMIT) + FIT_PADDING;
     51        if (isnan(model->radius)) {
     52          fprintf (stderr, "error in radius\n");
     53          continue;
     54        }
    2455        if (model->radius > OUTER) {
    2556          // allocate image, noise, mask arrays for each peak (square of radius OUTER)
    26           pmSourceDefinePixels (source, imdata, x, y, OUTER);
     57          pmSourceDefinePixels (source, imdata, x, y, model->radius);
    2758        }
    28 # endif
    29         psModel  *model  = pmSourceModelGuess (source, modelType);
    3059
    3160        // fit as FLT, not PSF (skip poor fits)
    32         if (!pmSourceFitModel (source, model, false)) continue;
     61        psImageKeepCircle (source->mask, x, y, model->radius, OR, 0x80);
     62        status = pmSourceFitModel (source, model, false);
     63        psImageKeepCircle (source->mask, x, y, model->radius, AND, 0x7f);
     64        if (!status || (model->params->data.F32[1] < 0)) {
     65          // if the fit fails, we need to change the classification
     66          psLogMsg ("psphot", 3, "GAL fit failed for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius);
     67          source->type = PS_SOURCE_OTHER;  // better choice?
     68          Nfail ++;
     69          continue;
     70        }
     71
    3372        source->modelFLT = model;
    34 
    3573        Niter += model[0].nIter;
    3674        Nfit++;
     75
     76        // subtract object, leave local sky
     77        sky = model->params->data.F32[0];
     78        model->params->data.F32[0] = 0;
     79        pmSourceSubModel (source->pixels, source->mask, model, false);
     80        model->params->data.F32[0] = sky;
     81
    3782    }
    38     psLogMsg ("psphot", 3, "fit galaxies: %f sec for %d galaxies (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter);
     83    psLogMsg ("psphot", 3, "fit galaxies: %f sec for %d galaxies (%d failures, %d total iterations)\n", psTimerMark ("psphot"), Nfit, Nfail, Niter);
    3984    return (true);
    4085}

  • trunk/psphot/src/mark_psf_source.c

    r4375 r4574  
    1515// PS_SOURCE_BRIGHTSTAR
    1616# define MIN_DS 0.01
    17 bool mark_psf_source (psSource *source, float shapeNsigma)
     17bool mark_psf_source (psSource *source, float shapeNsigma, float SATURATE)
    1818{
    1919    float dSX, dSY, SX, SY, SN;
    … …  
    2222    if (source->modelPSF == NULL) return (false);
    2323
    24     // SATSTAR objects are fitted with PSF regardless
    25     if (source->type == PS_SOURCE_SATSTAR) return (true);
     24    // if object has fitted peak above saturation, label as SATSTAR
     25    // remember: fit does not use saturated pixels (masked)
     26    if (source->modelPSF->params->data.F32[1] >= SATURATE) {
     27        if (source->type == PS_SOURCE_PSFSTAR) {
     28            psLogMsg ("psphot", 3, "PSFSTAR marked saturated\n");
     29        }
     30        source->type = PS_SOURCE_SATSTAR;
     31        return (true);
     32    }
     33    if (source->type == PS_SOURCE_SATSTAR) {
     34        psLogMsg ("psphot", 4, "SATSTAR marked bright (fitted peak below saturation)\n");
     35        source->type = PS_SOURCE_BRIGHTSTAR;
     36    }
    2637
    2738    SN  = source->modelPSF->params->data.F32[1]/source->modelPSF->dparams->data.F32[1];
    … …  
    3546    // sigma = 1 / SX
    3647    // dsx = 1 / (SX * SN)
    37     // dsx_o = hypot (1/SX*SN + MIN_DSX)
     48    // dsx_o = hypot (1/(SX*SN), MIN_DSX)
    3849
    3950    // assign PS_SOURCE_BRIGHTSTAR to bright objects within PSF region of dparams[]
    4051    if ((fabs(nSx) < shapeNsigma) && (fabs(nSy) < shapeNsigma)) {
    41       if (source->type == PS_SOURCE_PSFSTAR) return (true);
    42       source->type = PS_SOURCE_BRIGHTSTAR;
    43       return (true);
     52        if (source->type == PS_SOURCE_PSFSTAR) return (true);
     53        source->type = PS_SOURCE_BRIGHTSTAR;
     54        return (true);
    4455    }
    4556   
    4657    if (source->type == PS_SOURCE_PSFSTAR) {
    47       psLogMsg ("psphot", 3, "PSFSTAR demoted based on dSx, dSy\n");
     58        psLogMsg ("psphot", 3, "PSFSTAR demoted based on dSx, dSy\n");
    4859    }
    4960
    5061    if ((nSx >= shapeNsigma) || (nSy >= shapeNsigma)) {
    51       source->type = PS_SOURCE_GALAXY;
    52       return (false);
     62        source->type = PS_SOURCE_GALAXY;
     63        return (false);
    5364    }
    5465    // replace DEFECT with COSMIC?
    5566    if ((nSx <= -shapeNsigma) || (nSy <= -shapeNsigma)) {
    56       source->type = PS_SOURCE_DEFECT;
    57       return (false);
     67        source->type = PS_SOURCE_DEFECT;
     68        return (false);
    5869    }
    5970    psTrace (".psphot.mark_psf_source", 2, "unexpected result: unmarked object\n");

  • trunk/psphot/src/psphot-utils.c

    r4421 r4574  
    9696        model = (psModel  *) source->modelPSF;
    9797        if (model == NULL) continue;
     98        if (source->type == PS_SOURCE_GALAXY) continue;
     99        if (source->type == PS_SOURCE_DEFECT) continue;
     100        if (source->type == PS_SOURCE_SATURATED) continue;
    98101        params = model->params;
    99102        dparams = model->dparams;
    … …  
    107110            fprintf (f, "%9.6f ", dparams[0].data.F32[j+4]);
    108111        }
    109         fprintf (f, ": %2d %7.4f %7.4f %7.2f (%d %d)\n", source[0].type, log10(model[0].chisq), source[0].moments->SN, model[0].radius, model[0].nDOF, model[0].nIter);
     112        fprintf (f, ": %2d %7.4f %7.4f %7.2f (%d %d)\n",
     113                 source[0].type, log10(model[0].chisq), source[0].moments->SN,
     114                 model[0].radius, model[0].nDOF, model[0].nIter);
    110115    }
    111116    fclose (f);
    … …  
    134139        model = (psModel  *) source->modelFLT;
    135140        if (model == NULL) continue;
     141        if (source->type == PS_SOURCE_PSFSTAR) continue;
     142        if (source->type == PS_SOURCE_SATSTAR) continue;
     143        if (source->type == PS_SOURCE_BRIGHTSTAR) continue;
     144        if (source->type == PS_SOURCE_FAINTSTAR) continue;
     145        if (source->type == PS_SOURCE_OTHER) continue;
     146        if (source->type == PS_SOURCE_DEFECT) continue;
     147        if (source->type == PS_SOURCE_SATURATED) continue;
     148
    136149        params = model->params;
    137150        dparams = model->dparams;
    … …  
    287300}
    288301
     302/* this code is in pmObjects.h / psModules:psEllipse */
     303# if (0)
     304typedef struct {
     305  double major;
     306  double minor;
     307  double theta;
     308} EllipseAxes;
     309
     310typedef struct {
     311  double x2;
     312  double y2;
     313  double xy;
     314} EllipseMoments;
     315
     316typedef struct {
     317  double sx;
     318  double sy;
     319  double sxy;
     320} EllipseShape;
     321
     322EllipseAxes EllipseMomentsToAxes (EllipseMoments moments) {
     323
     324  EllipseAxes axes;
     325
     326  double f = sqrt (0.25*PS_SQR(moments.x2 - moments.y2) + PS_SQR(moments.sy));
     327
     328  axes.major = sqrt (0.5*(moments.x2 + moments.y2) + f);
     329  axes.minor = sqrt (0.5*(moments.x2 + moments.y2) - f)
     330  axes.theta = atan2 (2*moments.xy, moments.x2 - moments.y2) / 2;
     331  // theta in radians
     332
     333  return (axes);
     334}
     335
     336EllipseShape EllipseMomentsToShape (EllipseMoments moments) {
     337
     338  EllipseShape shape;
     339
     340  double f = sqrt(0.25*PS_SQR(axes.x2 - axes.y2) + PS_SQR(axes.sxy));
     341
     342  shape.sx = moments.y2 / f;
     343  shape.sy = moments.x2 / f;
     344  shape.sxy = -2.0*moments.xy / f;
     345  // XXX - these are probably wrong
     346
     347  return (shape);
     348}
     349
     350EllipseMoments EllipseShapeToMoments (EllipseShape shape) {
     351
     352  EllipseMoments moments;
     353
     354  return (moments);
     355}
     356
     357EllipseAxes EllipseShapeToAxes (EllipseShape shape) {
     358
     359  EllipseAxes axes;
     360
     361
     362  return (axes);
     363}
     364
     365EllipseShape EllipseAxesToShape (EllipseAxes axes) {
     366
     367  EllipseShape shape;
     368
     369  shape.sx = PS_SQR(cos(axes.theta) / axes.major) + PS_SQR(sin(axes.theta) / axes.minor);
     370  shape.sy = PS_SQR(sin(axes.theta) / axes.major) + PS_SQR(cos(axes.theta) / axes.minor);
     371  shape.sxy = cos(axes.theta)*sin(axes.theta)*(1.0/PS_SQR(axes.major) - 1.0/PS_SQR(axes.minor));
     372
     373  return (shape);
     374}
     375
     376EllipseMoments EllipseAxesToMoments (EllipseAxes axes) {
     377
     378  EllipseMoments moments;
     379
     380  return (moments);
     381}
     382
     383# endif

  • trunk/psphot/src/psphot.c

    r4421 r4574  
    4141
    4242    // fit extended objects with galaxy models
    43     // fit_galaxies (imdata, config, sources);
     43    fit_galaxies (imdata, config, sources, sky);
    4444
    4545    // subtract_galaxies (sources, config);

  • trunk/psphot/src/psphot.h

    r4375 r4574  
    4646bool mark_psf_sources (psArray *sources, psMetadata *config);
    4747bool subtract_psf_sources (psArray *sources);
    48 bool fit_galaxies (psImage *image, psMetadata *config, psArray *sources);
     48bool fit_galaxies (psImageData *imdata, psMetadata *config, psArray *sources, psStats *skyStats);
    4949bool subtract_galaxies (psArray *sources, psMetadata *config);
    5050bool subtract_psf_source (psSource *source);
    51 bool mark_psf_source (psSource *source, float shapeNsigma);
     51bool mark_psf_source (psSource *source, float shapeNsigma, float SATURATE);
    5252bool find_defects (psImage *zapmask, psArray *sources, psMetadata *config, pmPSF *psf);
     53bool basic_classes (psArray *sources, psMetadata *config);
    5354
    5455// psf utilities

  • trunk/psphot/src/source_moments.c

    r4398 r4574  
    55    bool     status  = false;
    66    psArray *sources = NULL;
    7     pmPSFClump psfClump;
    87
    98    psTimerStart ("psphot");

  • trunk/psphot/src/subtract_psf_source.c

    r4251 r4574  
    99
    1010  // non-stellar sources are ignored
     11  if (source->type == PS_SOURCE_OTHER) return (false);
    1112  if (source->type == PS_SOURCE_GALAXY) return (false);
    1213  if (source->type == PS_SOURCE_DEFECT) return (false);
     14  if (source->type == PS_SOURCE_SATURATED) return (false);
    1315
    1416  psModel *model = source->modelPSF;
    … …  
    1820  }         
    1921
    20   psImage  *pixels = source->pixels;
     22  psImage *pixels = source->pixels;
    2123
    2224  // subtract object, leave local sky
Note: See TracChangeset for help on using the changeset viewer.