IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 35630


Ignore:
Timestamp:
Jun 3, 2013, 6:13:52 PM (13 years ago)
Author:
eugene
Message:

add iterative model window for extended sources (crude and hackish still)

Location:
branches/eam_branches/ipp-20130509/psphot/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r34542 r35630  
    232232bool            psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal);
    233233bool            psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal);
    234 pmModel        *psphotFitEXT (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal);
     234pmModel        *psphotFitEXT (pmModel *guessModel, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal);
    235235psArray        *psphotFitDBL (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal);
    236236

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

    r35559 r35630  
    11# include "psphotInternal.h"
     2bool psphotSetWindowTrail (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float newRadius);
    23
    34// for now, let's store the detections on the readout->analysis for each readout
    … …  
    384385          assert (status);
    385386
    386           // limit selection to some SN limit
    387           // assert (source->peak); // how can a source not have a peak?
    388387          // limit selection to some SN limit
    389388          bool skipSource = false;
    … …  
    424423              }
    425424          } else {
    426               psFree (source->modelFlux);
    427               source->modelFlux = NULL;
    428               modelFit = psphotFitEXT (readout, source, fitOptions, modelType, maskVal, markVal);
    429               if (!modelFit) {
    430                   psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->Mx, source->moments->My);
    431                   Nfail ++;
    432                   continue;
    433               }
    434               psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n",
    435                        source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter);
    436               Nplain ++;
    437               if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
    438                   NplainPass ++;
    439                   source->mode |= PM_SOURCE_MODE_EXTENDED_FIT;
    440               }
     425              bool doneFits = false;
     426              while (!doneFits) {
     427                  psFree (source->modelFlux);
     428                  source->modelFlux = NULL;
     429                  modelFit = psphotFitEXT (modelFit, readout, source, fitOptions, modelType, maskVal, markVal);
     430                  if (!modelFit) {
     431                      psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->Mx, source->moments->My);
     432                      Nfail ++;
     433                      doneFits = true;
     434                      continue;
     435                  }
     436                  psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n", source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter);
     437                  Nplain ++;
     438                  if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
     439                      NplainPass ++;
     440                      source->mode |= PM_SOURCE_MODE_EXTENDED_FIT;
     441                  }
     442                  doneFits = true;
     443
     444                  if (modelType == pmModelClassGetType("PS_MODEL_TRAIL")) {
     445                      // calculate the object end points:
     446                      float *PAR = modelFit->params->data.F32;
     447                      float Xmin = PAR[PM_PAR_XPOS] - 0.5*PAR[PM_PAR_LENGTH]*sin(PAR[PM_PAR_THETA]);
     448                      float Ymin = PAR[PM_PAR_YPOS] - 0.5*PAR[PM_PAR_LENGTH]*cos(PAR[PM_PAR_THETA]);
     449                      float Xmax = PAR[PM_PAR_XPOS] + 0.5*PAR[PM_PAR_LENGTH]*sin(PAR[PM_PAR_THETA]);
     450                      float Ymax = PAR[PM_PAR_YPOS] + 0.5*PAR[PM_PAR_LENGTH]*cos(PAR[PM_PAR_THETA]);
     451
     452                      if ((source->peak->xf > 1100) &&
     453                          (source->peak->xf < 1400) &&
     454                          (source->peak->yf < 245)) {
     455                          fprintf (stderr, "src vs fit : %d %d - %d %d | %f %f - %f %f\n",
     456                                   source->pixels->col0, source->pixels->row0,
     457                                   source->pixels->col0 + source->pixels->numCols, source->pixels->row0 + source->pixels->numRows,
     458                                   Xmin, Ymin, Xmax, Ymax);
     459                      }
     460                      if (PAR[PM_PAR_LENGTH] > 0.9*fitRadius) {
     461                          doneFits = false;
     462                          fprintf (stderr, "update window : %f %f : %f -> %f\n", source->peak->xf, source->peak->yf, fitRadius, 2*fitRadius);
     463                          psphotSetWindowTrail (&fitRadius, &windowRadius, readout, source, markVal, fitRadius*2.0);
     464                      }
     465                  }
     466              }
    441467          }
    442 
    443           // XXX deprecate?
    444           // XXX pmSourceExtFitPars *extFitPars = pmSourceExtFitParsAlloc();
    445           // XXX extFitPars->Mxx = source->moments->Mxx;
    446           // XXX extFitPars->Mxy = source->moments->Mxy;
    447           // XXX extFitPars->Myy = source->moments->Myy;
    448           // XXX extFitPars->Mrf = source->moments->Mrf;
    449           // XXX extFitPars->Mrh = source->moments->Mrh;
    450           // XXX extFitPars->peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN;
    451          
    452           // save kron mag, but assign apMag & psfMag on output (not yet calculated)
    453           // XXX extFitPars->krMag = -2.5*log10(source->moments->KronFlux);
    454 
    455           // psArrayAdd (source->extFitPars, 4, extFitPars);
    456           // psFree (extFitPars);
     468          // XXX really need to do this in a cleaner way:
     469          if (!modelFit) continue;
    457470
    458471          // test for fit quality / result
    … …  
    545558    return true;
    546559}
     560
     561# define PAD_WINDOW 3.0
     562bool psphotSetWindowTrail (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float newRadius) {
     563
     564    psRegion newRegion;
     565
     566    psAssert (source, "source not defined??");
     567    psAssert (source->moments, "moments not defined??");
     568
     569    *fitRadius = newRadius;
     570    *windowRadius = *fitRadius + PAD_WINDOW;
     571
     572    // check to see if new region is completely contained within old region
     573    newRegion = psRegionForSquare (source->peak->xf, source->peak->yf, *windowRadius);
     574    newRegion = psRegionForImage (readout->image, newRegion);
     575
     576    // redefine the pixels to match
     577    pmSourceRedefinePixelsByRegion (source, readout, newRegion);
     578
     579    // clear the circular mask
     580    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     581
     582    // set the mask to flag the excluded pixels
     583    psImageKeepCircle (source->maskObj, source->peak->xf, source->peak->yf, *fitRadius, "OR", markVal);
     584
     585    return true;
     586}

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

    r35559 r35630  
    268268    {
    269269        // XXX need to handle failures better here
    270         EXT = psphotFitEXT (readout, source, fitOptions, modelTypeEXT, maskVal, markVal);
     270        EXT = psphotFitEXT (NULL, readout, source, fitOptions, modelTypeEXT, maskVal, markVal);
    271271        if (!EXT) goto escape;
    272272        if (!isfinite(EXT->params->data.F32[PM_PAR_I0])) goto escape;
    … …  
    443443}
    444444
    445 pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) {
     445pmModel *psphotFitEXT (pmModel *guessModel, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) {
    446446
    447447    if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) {
    … …  
    457457
    458458    // use the source moments, etc to guess basic model parameters
    459     pmModel *model = pmSourceModelGuess (source, modelType, maskVal, markVal);
     459    pmModel *model = guessModel ? guessModel : pmSourceModelGuess (source, modelType, maskVal, markVal);
    460460    if (!model) {
    461461        psTrace ("psphot", 5, "failed to generate a model for source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
Note: See TracChangeset for help on using the changeset viewer.