Changeset 35630 for branches/eam_branches/ipp-20130509/psphot/src/psphotExtendedSourceFits.c

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)

File:

• branches/eam_branches/ipp-20130509/psphot/src/psphotExtendedSourceFits.c (modified) (4 diffs)

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

@@ -1 +1 @@
 # include "psphotInternal.h"
+bool psphotSetWindowTrail (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float newRadius);
 
 // for now, let's store the detections on the readout->analysis for each readout
@@ -384 +385 @@
           assert (status);
 
-          // limit selection to some SN limit
-          // assert (source->peak); // how can a source not have a peak?
           // limit selection to some SN limit
           bool skipSource = false;
@@ -424 +423 @@
               }
           } else {
-              psFree (source->modelFlux);
-              source->modelFlux = NULL;
-              modelFit = psphotFitEXT (readout, source, fitOptions, modelType, maskVal, markVal);
-              if (!modelFit) {
-                  psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->Mx, source->moments->My);
-                  Nfail ++;
-                  continue;
-              }
-              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);
-              Nplain ++;
-              if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
-                  NplainPass ++;
-                  source->mode |= PM_SOURCE_MODE_EXTENDED_FIT;
-              }
+              bool doneFits = false;
+              while (!doneFits) {
+                  psFree (source->modelFlux);
+                  source->modelFlux = NULL;
+                  modelFit = psphotFitEXT (modelFit, readout, source, fitOptions, modelType, maskVal, markVal);
+                  if (!modelFit) {
+                      psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->Mx, source->moments->My);
+                      Nfail ++;
+                      doneFits = true;
+                      continue;
+                  }
+                  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);
+                  Nplain ++;
+                  if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
+                      NplainPass ++;
+                      source->mode |= PM_SOURCE_MODE_EXTENDED_FIT;
+                  }
+                  doneFits = true;
+
+                  if (modelType == pmModelClassGetType("PS_MODEL_TRAIL")) {
+                      // calculate the object end points:
+                      float *PAR = modelFit->params->data.F32;
+                      float Xmin = PAR[PM_PAR_XPOS] - 0.5*PAR[PM_PAR_LENGTH]*sin(PAR[PM_PAR_THETA]);
+                      float Ymin = PAR[PM_PAR_YPOS] - 0.5*PAR[PM_PAR_LENGTH]*cos(PAR[PM_PAR_THETA]);
+                      float Xmax = PAR[PM_PAR_XPOS] + 0.5*PAR[PM_PAR_LENGTH]*sin(PAR[PM_PAR_THETA]);
+                      float Ymax = PAR[PM_PAR_YPOS] + 0.5*PAR[PM_PAR_LENGTH]*cos(PAR[PM_PAR_THETA]);
+
+                      if ((source->peak->xf > 1100) && 
+                          (source->peak->xf < 1400) && 
+                          (source->peak->yf < 245)) {
+                          fprintf (stderr, "src vs fit : %d %d - %d %d | %f %f - %f %f\n", 
+                                   source->pixels->col0, source->pixels->row0, 
+                                   source->pixels->col0 + source->pixels->numCols, source->pixels->row0 + source->pixels->numRows, 
+                                   Xmin, Ymin, Xmax, Ymax);
+                      }
+                      if (PAR[PM_PAR_LENGTH] > 0.9*fitRadius) {
+                          doneFits = false;
+                          fprintf (stderr, "update window : %f %f : %f -> %f\n", source->peak->xf, source->peak->yf, fitRadius, 2*fitRadius);
+                          psphotSetWindowTrail (&fitRadius, &windowRadius, readout, source, markVal, fitRadius*2.0);
+                      }
+                  }
+              }
           }
-
-          // XXX deprecate?
-          // XXX pmSourceExtFitPars *extFitPars = pmSourceExtFitParsAlloc();
-          // XXX extFitPars->Mxx = source->moments->Mxx;
-          // XXX extFitPars->Mxy = source->moments->Mxy;
-          // XXX extFitPars->Myy = source->moments->Myy;
-          // XXX extFitPars->Mrf = source->moments->Mrf;
-          // XXX extFitPars->Mrh = source->moments->Mrh;
-          // XXX extFitPars->peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN;
-          
-          // save kron mag, but assign apMag & psfMag on output (not yet calculated)
-          // XXX extFitPars->krMag = -2.5*log10(source->moments->KronFlux);
-
-          // psArrayAdd (source->extFitPars, 4, extFitPars);
-          // psFree (extFitPars);
+          // XXX really need to do this in a cleaner way:
+          if (!modelFit) continue;
 
           // test for fit quality / result
@@ -545 +558 @@
     return true;
 }
+
+# define PAD_WINDOW 3.0
+bool psphotSetWindowTrail (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float newRadius) {
+
+    psRegion newRegion;
+
+    psAssert (source, "source not defined??");
+    psAssert (source->moments, "moments not defined??");
+
+    *fitRadius = newRadius;
+    *windowRadius = *fitRadius + PAD_WINDOW;
+
+    // check to see if new region is completely contained within old region
+    newRegion = psRegionForSquare (source->peak->xf, source->peak->yf, *windowRadius);
+    newRegion = psRegionForImage (readout->image, newRegion);
+
+    // redefine the pixels to match
+    pmSourceRedefinePixelsByRegion (source, readout, newRegion);
+
+    // clear the circular mask
+    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 
+
+    // set the mask to flag the excluded pixels
+    psImageKeepCircle (source->maskObj, source->peak->xf, source->peak->yf, *fitRadius, "OR", markVal);
+
+    return true;
+}