IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5802


Ignore:
Timestamp:
Dec 17, 2005, 10:26:59 AM (20 years ago)
Author:
eugene
Message:

clean up of small accounting bugs, adding the post-subtraction ApResid

Location:
trunk/psphot
Files:
1 added
17 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/Makefile

    r5772 r5802  
    5151$(SRC)/pmSourceFitFixed.$(ARCH).o  \
    5252$(SRC)/psSparse.$(ARCH).o            \
    53 $(SRC)/psImageData.$(ARCH).o
     53$(SRC)/psImageData.$(ARCH).o        \
     54$(SRC)/psphotApResid.$(ARCH).o
    5455
    5556PSMODULES = \
    … …  
    7273$(SRC)/modelTestFitSource.$(ARCH).o \
    7374$(SRC)/modelTestArguments.$(ARCH).o \
     75$(SRC)/psphotModelGroupInit.$(ARCH).o  \
    7476$(SRC)/psModulesUtils.$(ARCH).o     \
    7577$(SRC)/psImageData.$(ARCH).o        \
    … …  
    111113
    112114$(BIN)/%.$(ARCH) : $(SRC)/%.$(ARCH).o
    113         @if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi
    114         $(CC) $^ -o $@ $(LFLAGS)
     115        @if [ ! -d $(BIN) ]; then mkdir -p $(BIN) || exit; fi
     116        $(CC) $^ -o $@ $(LFLAGS) || exit
    115117        @echo "done with $@"
    116118
    117119$(DESTBIN)/%: $(BIN)/%.$(ARCH)
    118         @if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN); fi
    119         rm -f $(DESTBIN)/$*
    120         cp $(BIN)/$*.$(ARCH) $(DESTBIN)/$*
     120        @if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN) || exit; fi
     121        rm -f $(DESTBIN)/$* || exit
     122        cp $(BIN)/$*.$(ARCH) $(DESTBIN)/$* || exit
    121123
    122124$(INSTALL): % : $(BIN)/%.$(ARCH)
    … …  
    127129
    128130%.install:
    129         make $(DESTBIN)/$*
     131        make $(DESTBIN)/$* || exit
    130132
    131133# utilities #################################################
    132134
    133135install:
    134         for i in $(INSTALL); do make $$i.install; done
     136        for i in $(INSTALL); do make $$i.install || exit; done
    135137
    136138clean: 

  • trunk/psphot/src/modelTestFitSource.c

    r5755 r5802  
    6767    source->peak->counts = source->moments->Peak;
    6868
     69    fprintf (stderr, "sum: %f\n", source->moments->Sum);
     70
    6971    // get the initial model parameter guess
    7072    pmModel *model = pmSourceModelGuess (source, modelType);
    … …  
    7375    psF32 *params = model->params->data.F32;
    7476    for (int i = 0; i < nParams; i++) {
     77        if (i == 2) {
     78            params[i] = xObj;
     79            continue;
     80        }
     81        if (i == 3) {
     82            params[i] = yObj;
     83            continue;
     84        }
    7585        sprintf (name, "TEST_FIT_PAR%d", i);
    7686        value = psMetadataLookupF32 (&status, config, name);
    7787        if (status) {
     88            fprintf (stderr, "using supplied value %f for PAR %d\n", value, i);
    7889            params[i] = value;
     90        } else {
     91            fprintf (stderr, "using guessed  value %f for PAR %d\n", params[i], i);
    7992        }
    8093    }
     94
     95    float area = params[4]*params[5];
     96    fprintf (stderr, "peak: %f\n", source->moments->Sum*area);
    8197
    8298    // what fitting mode to use?
    … …  
    98114    pmSourceSubModel (source->pixels, source->mask, model, false, false);
    99115   
     116    for (int i = 0; i < nParams; i++) {
     117        fprintf (stderr, "result value %f for PAR %d\n", params[i], i);
     118    }
     119
    100120    // write out
    101121    DumpImage (source->pixels, "resid.fits");

  • trunk/psphot/src/psModulesUtils.c

    r5772 r5802  
    6565    float sky = model->params->data.F32[0];
    6666
     67    *fitMag = 99.0;
     68    *obsMag = 99.0;
     69
    6770    pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type);
    6871    fitSum = modelFluxFunc (model->params);
     72    if (fitSum <= 0) return false;
     73    if (!isfinite(fitSum)) return false;
     74    *fitMag = -2.5*log10(fitSum);
    6975
    7076    for (int ix = 0; ix < image->numCols; ix++) {
    … …  
    7581    }
    7682    if (obsSum <= 0) return false;
    77     if (fitSum <= 0) return false;
     83    *obsMag = -2.5*log10(obsSum);
    7884
    79     *fitMag = -2.5*log10(fitSum);
    80     *obsMag = -2.5*log10(obsSum);
    8185    return (true);
    8286}

  • trunk/psphot/src/psphot.c

    r5773 r5802  
    4646    pmSourceRoughClass (sources, config, psfClump);
    4747
    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     }
    55     if (!strcasecmp (breakPt, "CLASS")) exit (0);
    56 
    5748    // mark blended peaks PS_SOURCE_BLEND
    5849    psphotBasicDeblend (sources, config, sky);
    … …  
    6152    // source analysis is done in S/N order (brightest first)
    6253    sources = psArraySort (sources, psphotSortBySN);
     54
     55    psphotDumpMoments (config, sources);
     56    if (!strcasecmp (breakPt, "CLASS")) exit (0);
    6357
    6458    // use bright stellar objects to measure PSF
    … …  
    7569      case 0:
    7670        psphotEnsemblePSF (imdata, config, sources, psf, sky);
    77         psphotFullFit (imdata, config, sources, psf, sky);
    78 //      psphotReplaceUnfit (sources);
    7971        break;
    8072
    8173      case 1:
     74        psphotEnsemblePSF (imdata, config, sources, psf, sky);
     75        psphotFullFit (imdata, config, sources, psf, sky);
     76        psphotReplaceUnfit (sources);
     77        psphotApResid (sources, config, psf);
     78        break;
     79
     80      case 2:
    8281        // fit extended objects with galaxy models
    8382        psphotApplyPSF (imdata, config, sources, psf, sky);
    … …  
    8584        break;
    8685
    87       case 2:
     86      case 3:
    8887        // fit extended objects with galaxy models
    8988        psphotApplyPSF (imdata, config, sources, psf, sky);
    9089        break;
    9190
    92       case 3:
     91      case 4:
    9392        // fit extended objects with galaxy models
    9493        psphotFixedPSF (imdata, config, sources, psf, sky);

  • trunk/psphot/src/psphot.h

    r5772 r5802  
    8585bool psphotSamplePSFs (pmPSF *psf, psImage *image);
    8686bool psphotReplaceUnfit (psArray *sources);
     87bool psphotDumpMoments (psMetadata *config, psArray *sources);
     88bool psphotApResid (psArray *sources, psMetadata *config, pmPSF *psf);

  • trunk/psphot/src/psphotChoosePSF.c

    r5772 r5802  
    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
     68   
     69    // use the best model:
     70    try = models->data[bestN];
    7071
    71     // keep only the selected model:
    72     try = models->data[bestN];
     72    // unset the PSFSTAR flag for stars not used for PSF model
     73    for (int i = 0; i < try->sources->n; i++) {
     74        pmSource *source = try->sources->data[i];
     75        if (try->mask->data.U8[i] & PSFTRY_MASK_FLT_FAIL) {
     76            source->mode &= ~PM_SOURCE_PSFSTAR;
     77        }
     78    }
     79
     80    // save only the best model;
    7381    pmPSF *psf = psMemCopy(try->psf);
    74     psFree (models);                             // keep only the pmPSF resulting from this analysis
     82    psFree (models);
    7583
    7684    modelName = pmModelGetType (psf->type);

  • trunk/psphot/src/psphotEnsemblePSF.c

    r5772 r5802  
    66    float x;
    77    float y;
    8     float Sky;
    98
    109    psTimerStart ("psphot");
    … …  
    1413
    1514    float OUTER_RADIUS     = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
     15    float INNER_RADIUS     = psMetadataLookupF32 (&status, config, "SKY_INNER_RADIUS");
    1616    float PSF_FIT_NSIGMA   = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA");
    1717    float PSF_FIT_PADDING  = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING");
    … …  
    3838        pmSource *otSource = pmSourceAlloc ();
    3939
     40        // really saturated stars should be re-measured for a better centroid
     41        // XXX EAM : move this to a 'clear satstar function'
     42        // XXX EAM : extend size of fit box around SATSTAR
     43        if (inSource->mode &  PM_SOURCE_SATSTAR) {
     44            status = pmSourceMoments (inSource, INNER_RADIUS);
     45        }
     46
    4047        // XXX EAM : add option to use FLT or PSF form
    4148        // use the source moments, etc to guess basic model parameters
    4249        pmModel *modelFLT = pmSourceModelGuess (inSource, psf->type);
     50        if (inSource->mode &  PM_SOURCE_SATSTAR) {
     51            modelFLT->params->data.F32[2] = inSource->moments->x;
     52            modelFLT->params->data.F32[3] = inSource->moments->y;
     53        }
     54        // XXX EAM : add option to peak-up on peak (for non-sat objects)
    4355
    4456        // set PSF parameters for this model
    … …  
    6981        psImage *flux = otSource->pixels;
    7082        psImage *mask = otSource->mask;
    71 
    72         // XXX EAM : use these lines to fit to the peak
    73         // model->params->data.F32[2] = inSource->peak->x;
    74         // model->params->data.F32[3] = inSource->peak->y;
    75         // XXX EAM : better option: improve the peak with 2D poly fit 3x3
    7683
    7784        // set model to unit peak, zero sky (we assume sky is constant)
    … …  
    136143
    137144        Fi->modelPSF = Mi->modelPSF;
     145
     146        // assign linearly-fitted normalization
    138147        Fi->modelPSF->params->data.F32[1] = norm->data.F32[i];
    139148
    140149        // subtract object
    141150        pmSourceSubModel (Fi->pixels, Fi->mask, Fi->modelPSF, false, false);
    142 
    143         Fi->modelPSF->params->data.F32[0] = Sky;
    144         // need to set this!
    145151    }
    146152

  • trunk/psphot/src/psphotEvalFLT.c

    r5772 r5802  
    3737
    3838    // if the object has a fitted peak below 0, the fit did not converge cleanly
    39     if (model->params->data.F32[1] < 0) {
     39    if (model->params->data.F32[1] <= 0) {
    4040        source->mode |= PM_SOURCE_FAIL;
    4141        return false;

  • trunk/psphot/src/psphotEvalPSF.c

    r5772 r5802  
    9494
    9595    // if the object has a fitted peak below 0, the fit did not converge cleanly
    96     if (model->params->data.F32[1] < 0) {
     96    if (model->params->data.F32[1] <= 0) {
    9797        source->mode |= PM_SOURCE_FAIL;
    9898        return false;
    … …  
    125125    if (keep) return true;
    126126
    127     // this source is not a star, unflag as PSFSTAR
    128     // XXX : if this object was used to build the PSF, this flag should
    129     //       be set even if the object is not a star...
     127    // this source is not a star, warn if it was a PSFSTAR
    130128    if (source->mode & PM_SOURCE_PSFSTAR) {
    131         source->mode &= ~PM_SOURCE_PSFSTAR;
    132         psLogMsg ("psphot", 5, "PSFSTAR demoted based on fit quality\n");
     129        psphotSaveImage (NULL, source->pixels, "failpx.fits");
     130        psphotSaveImage (NULL, source->mask, "failmk.fits");
     131        psphotSaveImage (NULL, source->weight, "failwt.fits");
     132        psLogMsg ("psphot", 5, "PSFSTAR demoted based on fit quality   (%f, %f  :  %f %f %f %f)\n",
     133                  model->params->data.F32[2], model->params->data.F32[3], nSx, nSy, SN, Chi);
    133134    }
    134135

  • trunk/psphot/src/psphotFullFit.c

    r5773 r5802  
    11# include "psphot.h"
     2
     3# define CLEAR_TRACE \
     4  if (TEST_SRC_ON) { \
     5    psTraceSetLevel (".pmObjects.pmSourceFitModel", 0); \
     6    psTraceSetLevel ("psMinimizeLMChi2", 0); \
     7    TEST_SRC_ON = false; } \
    28
    39bool psphotFullFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky)
    410{
    5     bool  skip, status;
     11    bool  status;
    612    float x, y;
    713    int   Nfit = 0;
    … …  
    2632    psphotInitRadiusFLT (config, sky, modelTypeFLT);
    2733
     34    float TEST_SRC_X = psMetadataLookupF32 (&status, config, "TEST_SRC_X");
     35    float TEST_SRC_Y = psMetadataLookupF32 (&status, config, "TEST_SRC_Y");
     36    bool TEST_SRC_ON = false;
     37    bool TEST_SRC    = status;
     38    fprintf (stderr, "test src: %f %f\n", TEST_SRC_X, TEST_SRC_Y);
     39
    2840    for (int i = 0; i < sources->n; i++) {
    2941
    … …  
    3749
    3850        // save the PSF model from the Ensemble fit
    39         pmModel *save  = pmModelCopy (source->modelPSF);
     51        pmModel *LIN  = pmModelCopy (source->modelPSF);
    4052
    4153        // attempt to fit the PSF model
    … …  
    4961        x = modelPSF->params->data.F32[2];
    5062        y = modelPSF->params->data.F32[3];
     63
     64        if (TEST_SRC && (fabs(TEST_SRC_X - x) < 5) && (fabs(TEST_SRC_Y - y) < 5)) {
     65            psTraceSetLevel (".pmObjects.pmSourceFitModel", 6);
     66            psTraceSetLevel ("psMinimizeLMChi2", 6);
     67            TEST_SRC_ON = true;
     68        }
    5169
    5270        // fit PSF model (set/unset the pixel mask)
    … …  
    6381            pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false);
    6482            source->mode |= PM_SOURCE_SUBTRACTED;
    65             psFree (save);
     83            psFree (LIN);
    6684            Nsub ++;
     85            CLEAR_TRACE;
    6786            continue;
    6887        }
    6988
    7089        // skip the source if we don't think it is extended
    71         skip  = (source->type != PM_SOURCE_GALAXY);
    72         skip |= (source->moments->SN < GAL_MIN_SN);
    73         if (skip) {
    74             source->modelPSF = save;
    75             pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false);
    76             source->mode |= PM_SOURCE_SUBTRACTED;
    77             psFree (modelPSF);
    78             Nsub ++;
    79             continue;
    80         }
     90        if (source->type != PM_SOURCE_GALAXY) goto subLINEAR;
     91        if (source->moments->SN < GAL_MIN_SN) goto subLINEAR;
    8192
    8293        // add the double-PSF fit mode
    … …  
    8495       
    8596        // recalculate the source moments using the larger galaxy moments radius
    86         if (!pmSourceMoments (source, GAL_MOMENTS_RAD)) {
    87             source->mode |= PM_SOURCE_FAIL;  // better choice?
    88             continue;
    89         }
     97        if (!pmSourceMoments (source, GAL_MOMENTS_RAD)) goto subLINEAR;
    9098
    9199        // use the source moments, etc to guess basic model parameters
    … …  
    111119            pmSourceSubModel (source->pixels, source->mask, source->modelFLT, false, false);
    112120            source->mode |= PM_SOURCE_SUBTRACTED;
    113             psFree (save);
     121            psFree (LIN);
    114122            Nsub ++;
     123            CLEAR_TRACE;
    115124            continue;
    116125        }
    117126
    118127        // subtract PSF for object, leave local sky
    119         // XXX : we should keep the non-linear PSF fit if object was marked extended but
    120         //       PSF fit is OK otherwise.
     128    subLINEAR:
    121129        psFree (source->modelPSF);
    122         source->modelPSF = save;
     130        source->modelPSF = LIN;
    123131        pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false);
    124132        source->mode |= PM_SOURCE_SUBTRACTED;
     133        source->mode |= PM_SOURCE_LINEAR;
     134        Nsub ++;
     135        CLEAR_TRACE;
    125136    }
    126137

  • trunk/psphot/src/psphotMagnitudes.c

    r5773 r5802  
    88    int status;
    99    float x, y;
    10     float rflux, apMag, fitMag;
     10    float rflux;
    1111    pmModel *model;
    1212

  • trunk/psphot/src/psphotModelTest.c

    r5622 r5802  
    22
    33int main (int argc, char **argv) {
     4
     5    psphotModelGroupInit ();
    46
    57    psMetadata *config = modelTestArguments (&argc, argv);

  • trunk/psphot/src/psphotOutput.c

    r5773 r5802  
    370370bool pmModelWritePSFs (psArray *sources, psMetadata *config, char *filename, pmPSF *psf) {
    371371
    372     double dP, dmag;
     372    double dPos, dMag;
    373373    int i, j;
    374374    FILE *f;
    … …  
    393393        // if (source->mode & PM_SOURCE_POOR) continue;
    394394
    395         model = pmSourceMagnitudes (source, psf, RADIUS, true);
     395        model = pmSourceMagnitudes (source, psf, RADIUS);
    396396        if (model == NULL) continue;
    397397
    … …  
    399399        dPAR = model->dparams->data.F32;
    400400       
    401         // dP is positional error
    402         dP = 0;
    403         dP += PS_SQR(dPAR[2]);
    404         dP += PS_SQR(dPAR[3]);
    405         dP = sqrt (dP);
    406        
    407         dmag = dPAR[1] / PAR[1];
    408 
    409         fprintf (f, "%7.1f %7.1f  %5.1f %7.4f  %7.4f %7.4f  ",
    410                  PAR[2], PAR[3], PAR[0], source->fitMag, dmag, dP);
     401        // dPos is positional error, dMag is mag error
     402        dPos = hypot (dPAR[2], dPAR[3]);
     403        dMag = dPAR[1] / PAR[1];
     404
     405        fprintf (f, "%7.1f %7.1f  %5.1f %8.4f  %7.4f %7.4f  ",
     406                 PAR[2], PAR[3], PAR[0], source->fitMag, dMag, dPos);
    411407
    412408        for (j = 4; j < model->params->n; j++) {
    … …  
    417413            fprintf (f, "%9.6f ", dPAR[j]);
    418414        }
    419         fprintf (f, ": %2d %2d %7.3f %7.3f %7.2f %4d %2d\n",
     415        fprintf (f, ": %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n",
    420416                 source[0].type, source[0].mode,
    421417                 log10(model[0].chisq/model[0].nDOF),
    … …  
    432428bool pmModelWriteFLTs (psArray *sources, char *filename) {
    433429
    434     double dP;
     430    double dPos, dMag;
    435431    int i, j;
    436432    FILE *f;
    437     psVector *params;
    438     psVector *dparams;
     433    psF32 *PAR, *dPAR;
    439434    pmModel  *model;
    440435
    … …  
    450445
    451446        if (source->type != PM_SOURCE_GALAXY) continue;
    452         // if (source->mode & PM_SOURCE_FAIL) continue;
    453         // if (source->mode & PM_SOURCE_POOR) continue;
    454447       
    455         model = pmSourceMagnitudes (source->modelFLT, NULL, model->radius);
    456         if (model == NULL) continue;
     448        if (source->modelFLT == NULL) continue;
     449        model = pmSourceMagnitudes (source, NULL, source->modelFLT->radius);
    457450
    458451        PAR  = model->params->data.F32;
    459452        dPAR = model->dparams->data.F32;
    460453       
    461         // dP is shape error
     454        // dPos is shape error
    462455        // XXX these are hardwired for SGAUSS
    463         dP = 0;
    464         dP += PS_SQR(dPAR[4] / PAR[4]);
    465         dP += PS_SQR(dPAR[5] / PAR[5]);
    466         dP += PS_SQR(dPAR[7] / PAR[7]);
    467         dP = sqrt (dP);
    468 
    469         dmag = dPAR[1] / PAR[1];
    470 
    471         fprintf (f, "%7.1f %7.1f  %5.1f %7.3f  %7.4f %7.4f  ",
    472                  PAR[2], PAR[3], PAR[0], source->fitMag, dmag, dP);
     456        dPos = 0;
     457        dPos += PS_SQR(dPAR[4] / PAR[4]);
     458        dPos += PS_SQR(dPAR[5] / PAR[5]);
     459        // dPos += PS_SQR(dPAR[7] / PAR[7]);
     460        dPos = sqrt (dPos);
     461        dMag = dPAR[1] / PAR[1];
     462
     463        fprintf (f, "%7.1f %7.1f  %7.1f %7.3f  %7.4f %7.4f  ",
     464                 PAR[2], PAR[3], PAR[0], source->fitMag, dMag, dPos);
    473465
    474466        for (j = 4; j < model->params->n; j++) {
    … …  
    479471            fprintf (f, "%9.6f ", dPAR[j]);
    480472        }
    481         fprintf (f, ": %2d %2d %7.3f %7.3f %7.2f %4d %2d\n",
     473        fprintf (f, ": %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n",
    482474                 source[0].type, source[0].mode,
    483475                 log10(model[0].chisq/model[0].nDOF),
    … …  
    548540    }
    549541
     542    fprintf (stderr, "writing out moments\n");
    550543    for (i = 0; i < sources->n; i++) {
    551544        source = sources->data[i];
    552545        if (source->moments == NULL)
    553546            continue;
    554         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",
     547        fprintf (f, "%5d %5d  %7.1f  %7.1f %7.1f  %6.3f %6.3f  %10.1f %7.1f %7.1f %7.1f  %4d %2d %#5x\n",
    555548                 source->peak->x, source->peak->y, source->peak->counts,
    556549                 source->moments->x, source->moments->y,
    … …  
    558551                 source->moments->Sum, source->moments->Peak,
    559552                 source->moments->Sky, source->moments->SN,
    560                  source->moments->nPixels, source->type);
     553                 source->moments->nPixels, source->type, source->mode);
    561554    }
    562555    fclose (f);
    … …  
    599592
    600593    psRegion region = {x - dx, x + dx, y - dy, y + dy};
    601     psImage *sample = psImageSubset (image, region);
     594    psImage *view = psImageSubset (image, region);
     595    psImage *sample = psImageCopy (NULL, view, PS_TYPE_F32);
    602596    psImageInit (sample, 0);
    603597    modelFLT->params->data.F32[2] = x;
    … …  
    635629    return (TRUE);
    636630}
     631
     632bool psphotDumpMoments (psMetadata *config, psArray *sources) {
     633
     634    bool status;
     635
     636    // optional dump of all rough source data
     637    char *output = psMetadataLookupPtr (&status, config, "MOMENTS_OUTPUT_FILE");
     638    if (!status) return false;
     639    if (output == NULL) return false;
     640    if (output[0] == 0) return false;
     641
     642    pmMomentsWriteText (sources, output);
     643    psFree (output);
     644    return true;
     645}

  • trunk/psphot/src/psphotRadiusChecks.c

    r5772 r5802  
    2929    if (isnan(model->radius)) psAbort ("apply_psf_model", "error in radius");
    3030       
     31    if (source->mode &  PM_SOURCE_SATSTAR) {
     32        model->radius *= 2;
     33    }
     34
    3135    // check if we need to redefine the pixels
    3236    // XXX EAM : a better test would examine the source pixels

  • trunk/psphot/src/psphotReplaceUnfit.c

    r5774 r5802  
    33bool psphotReplaceUnfit (psArray *sources) {
    44
    5     bool status;
    6     float x;
    7     float y;
     5    pmModel *model;
     6    pmSource *source;
    87
    98    psTimerStart ("psphot");
    109
    1110    for (int i = 0; i < sources->n; i++) {
    12         pmSource *source = sources->data[i];
     11      source = sources->data[i];
    1312
    1413        // skip non-astronomical objects (very likely defects)
    1514        // these were never fitted and subtracted
    1615        if (source->mode &  PM_SOURCE_BLEND) continue;
    17         if (source->type == PM_SOURCE_DEFECT) continue;
    1816        if (source->type == PM_SOURCE_SATURATED) continue;
     17
     18        if (source->type == PM_SOURCE_DEFECT) {
     19            if (source->mode &  PM_SOURCE_SUBTRACTED) goto addPSF;
     20            continue;
     21        }
     22
     23        if (source->type == PM_SOURCE_GALAXY) {
     24            if (source->mode & PM_SOURCE_FAIL) goto addPSF;
     25            if (source->mode & PM_SOURCE_POOR) goto addPSF;
     26            continue;
     27        }
    1928
    2029        // need to skip the successful fits
    2130        if (source->type == PM_SOURCE_STAR) {
    22             if (source->mode & PM_SOURCE_FAIL) goto subPSF;
    23             if (source->mode & PM_SOURCE_POOR) goto subPSF;
    24             continue;
    25 
    26         subPSF:
    27             pmModel  *model = source->modelPSF;
    28 
    29             x = model->params->data.F32[2];
    30             y = model->params->data.F32[3];
    31 
    32             psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
    33             pmSourceAddModel (source->pixels, source->mask, model, false, false);
    34             psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED);
     31            if (source->mode & PM_SOURCE_FAIL) goto addPSF;
     32            if (source->mode & PM_SOURCE_POOR) goto addPSF;
    3533            continue;
    3634        }
    3735
    38         if (source->type == PM_SOURCE_GALAXY) {
    39             if (source->mode & PM_SOURCE_FAIL) goto subFLT;
    40             if (source->mode & PM_SOURCE_POOR) goto subFLT;
    41             continue;
     36    addPSF:
     37        model = source->modelPSF;
     38        if (model == NULL) continue;
    4239
    43         subPSF:
    44             pmModel  *model = source->modelFLT;
    45 
    46             x = model->params->data.F32[2];
    47             y = model->params->data.F32[3];
    48 
    49             psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
    50             pmSourceAddModel (source->pixels, source->mask, model, false, false);
    51             psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED);
    52             continue;
    53         }
     40        pmSourceAddModel (source->pixels, source->mask, model, false, false);
     41        source->mode &= ~PM_SOURCE_SUBTRACTED;
    5442    }
    5543    psLogMsg ("psphot.replace", 4, "replace models: %f (%d objects)\n", psTimerMark ("psphot"), sources->n);

  • trunk/psphot/src/psphotSetup.c

    r5754 r5802  
    5555        for (int iy = 0; iy < image->numRows; iy++) {
    5656            for (int ix = 0; ix < image->numCols; ix++) {
    57                 weight->data.F32[iy][ix] = image->data.F32[iy][ix] / GAIN + PS_SQR(RDNOISE/GAIN);
     57                weight->data.F32[iy][ix] = PS_MAX (image->data.F32[iy][ix] / GAIN + PS_SQR(RDNOISE/GAIN), 0.0);
    5858            }
    5959        }
    … …  
    104104    psLogMsg ("psphot", 4, "load data: %f sec\n", psTimerMark ("psphot"));
    105105
     106    psphotSaveImage (NULL, weight, "weight.fits");
     107    psphotSaveImage (NULL, mask, "mask.fits");
     108
    106109    // save the image data & return it
    107110    eamReadout *imdata = eamReadoutAlloc(image, weight, mask, header);

  • trunk/psphot/src/psphotSourceStats.c

    r5617 r5802  
    55    bool     status  = false;
    66    psArray *sources = NULL;
     7    float BIG_RADIUS;
    78
    89    psTimerStart ("psphot");
    … …  
    3536
    3637        // measure basic source moments
    37         // XXX EAM : choose between these two versions
    3838        status = pmSourceMoments (source, RADIUS);
    39         if (!status) {
    40           psFree (source);
    41           continue;
     39        if (status) {
     40            // add to the source array
     41            psArrayAdd (sources, 100, source);
     42            psFree (source);
     43            continue;
    4244        }
    43        
    44         // add to the source array
    45         psArrayAdd (sources, 100, source);
     45
     46        // if no valid pixels, or massive swing, likely saturated source,
     47        // try a much larger box
     48        BIG_RADIUS = PS_MIN (INNER, 3*RADIUS);
     49        psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
     50        status = pmSourceMoments (source, BIG_RADIUS);
     51        if (status) {
     52            // add to the source array
     53            psArrayAdd (sources, 100, source);
     54            psFree (source);
     55            continue;
     56        }
     57
    4658        psFree (source);
     59        continue;
    4760    }
    4861     
Note: See TracChangeset for help on using the changeset viewer.