Changeset 34258

Timestamp:

Jul 31, 2012, 3:58:10 PM (14 years ago)

Author:

eugene

Message:

add TRAIL models (and some trace test output); create psphotMinimal to test the psphotReadoutMinimal variant

Location:

trunk/psphot

Files:

• . (modified) (1 prop)
• src (modified) (2 props)
• src/Makefile.am (modified) (3 diffs)
• src/psphot.h (modified) (2 diffs)
• src/psphotBlendFit.c (modified) (1 diff)
• src/psphotExtendedSourceFits.c (modified) (2 diffs)
• src/psphotForced.c (modified) (1 diff)
• src/psphotImageLoop.c (modified) (1 diff)
• src/psphotMinimal.c (copied) (copied from branches/eam_branches/ipp-20120627/psphot/src/psphotMinimal.c )
• src/psphotMinimalArguments.c (copied) (copied from branches/eam_branches/ipp-20120627/psphot/src/psphotMinimalArguments.c )
• src/psphotReadoutMinimal.c (modified) (2 diffs)
• src/psphotSourceFits.c (modified) (2 diffs)
• src/psphotSourceSize.c (modified) (2 diffs)
• src/psphotStackImageLoop.c (modified) (1 prop)

trunk/psphot/src

@@ -23 +23 @@
 psphotStack
 psphotModelTest
+psphotMinimal

@@ -23 +23 @@
 psphotStack
 psphotModelTest
+psphotMinimal

trunk/psphot/src/Makefile.am

@@ -25 +25 @@
 libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
 
-bin_PROGRAMS = psphot psphotForced psphotMakePSF psphotStack psphotModelTest
+bin_PROGRAMS = psphot psphotForced psphotMinimal psphotMakePSF psphotStack psphotModelTest
 # bin_PROGRAMS = psphotPetrosianStudy psphotTest psphotMomentsStudy 
 
@@ -35 +35 @@
 psphotForced_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
 psphotForced_LDADD = libpsphot.la
+
+psphotMinimal_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
+psphotMinimal_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
+psphotMinimal_LDADD = libpsphot.la
 
 psphotMakePSF_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
@@ -73 +77 @@
         psphotForced.c             \
         psphotForcedArguments.c    \
+        psphotParseCamera.c        \
+        psphotImageLoop.c          \
+        psphotMosaicChip.c         \
+        psphotCleanup.c
+
+# forced photometry of specified positions given a specified psf
+psphotMinimal_SOURCES = \
+        psphotMinimal.c            \
+        psphotMinimalArguments.c           \
         psphotParseCamera.c        \
         psphotImageLoop.c          \

trunk/psphot/src/psphot.h

@@ -20 +20 @@
     PSPHOT_MAKE_PSF,
     PSPHOT_MODEL_TEST,
+    PSPHOT_MINIMAL,
 } psphotImageLoopMode;
 
@@ -332 +333 @@
 bool psphotForcedReadout(pmConfig *config, const pmFPAview *view, const char *filerule);
 
+pmConfig *psphotMinimalArguments(int argc, char **argv);
+bool psphotReadoutMinimal(pmConfig *config, const pmFPAview *view, const char *filerule);
+
 pmConfig *psphotMakePSFArguments(int argc, char **argv);
 bool psphotMakePSFReadout(pmConfig *config, const pmFPAview *view, const char *filerule);

trunk/psphot/src/psphotBlendFit.c

@@ -61 +61 @@
     if (!status) {
         nThreads = 0;
+    }
+
+    // do not thread if we are trying to study the fitting process
+    if (psTraceGetLevel ("psphot.psphotFitEXT") >= 6) {
+      nThreads = 0;
     }
 

trunk/psphot/src/psphotExtendedSourceFits.c

@@ -77 +77 @@
     if (!status) {
         nThreads = 0;
+    }
+    // do not thread if we are trying to study the fitting process
+    if (psTraceGetLevel ("psphot.psphotFitEXT") >= 6) {
+      nThreads = 0;
     }
 
@@ -468 +472 @@
         psFree (source->modelEXT);
         source->modelEXT = psMemIncrRefCounter (source->modelFits->data[minModel]);
+        source->type = PM_SOURCE_TYPE_EXTENDED;
+        source->mode |= PM_SOURCE_MODE_EXTMODEL;
+        source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
 
         // adjust the window so the subtraction covers the faint wings

trunk/psphot/src/psphotForced.c

@@ -1 +1 @@
 # include "psphotStandAlone.h"
-# define FORCED_PHOTOMETRY 1
 
 int main (int argc, char **argv) {

trunk/psphot/src/psphotImageLoop.c

@@ -118 +118 @@
                     }
                     break;
+                  case PSPHOT_MINIMAL:
+                    if (!psphotReadoutMinimal (config, view, "PSPHOT.INPUT")) {
+                        psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
+                        psFree (view);
+                        return false;
+                    }
+                    break;
                   case PSPHOT_FORCED:
                     if (!psphotForcedReadout (config, view, "PSPHOT.INPUT")) {

trunk/psphot/src/psphotReadoutMinimal.c

@@ -25 +25 @@
     // Generate the mask and weight images, including the user-defined analysis region of interest
     psphotSetMaskAndVariance (config, view, filerule);
+
+    // only subtract background if needed?
+    // activate this for a clean test with psphotMinimal. (add to recipe!)
+    if (0) {
+      // generate a background model (median, smoothed image)
+      if (!psphotModelBackground (config, view, filerule)) {
+        return psphotReadoutCleanup (config, view, filerule);
+      }
+      if (!psphotSubtractBackground (config, view, filerule)) {
+        return psphotReadoutCleanup (config, view, filerule);
+      }
+    }
 
     // load the psf model, if suppled.  FWHM_X,FWHM_Y,etc are saved on readout->analysis
@@ -72 +84 @@
     psphotFitSourcesLinear (config, view, filerule, false, false);
 
-// XXX eventually, add the extended source fits here
-# if (0)
+    // measure the radial profiles to the sky
+    psphotRadialProfileWings (config, view, filerule);
+
+    // re-measure the kron mags with models subtracted and more appropriate windows
+    psphotKronIterate(config, view, filerule);
+
     // measure source size for the remaining sources
-    psphotSourceSize (config, view, filerule);
+    psphotSourceSize (config, view, filerule, false);
 
-    psphotExtendedSourceAnalysis (config, view, filerule);
+    // NOTE: Petrosian and Isophotal mags are not relevant at this time
+    // psphotExtendedSourceAnalysis (config, view, filerule);
 
+    // in ppSub context, this is used to fit TRAILs (and maybe EXP for comets)
     psphotExtendedSourceFits (config, view, filerule);
-# endif
 
     // calculate source magnitudes

trunk/psphot/src/psphotSourceFits.c

@@ -463 +463 @@
     }
 
+# if (PS_TRACE_ON) 
+    if (psTraceGetLevel ("psphot") >= 6) {
+
+      // Moments-based shapes parameters
+      psEllipseMoments moments;
+      moments.x2 = source->moments->Mxx;
+      moments.xy = source->moments->Mxy;
+      moments.y2 = source->moments->Myy;
+      // force the axis ratio to be < 20.0
+      psEllipseAxes axes = psEllipseMomentsToAxes (moments, 20.0);
+
+      // PSF shape parameters
+      psEllipseShape psfShape;
+      psfShape.sx  = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2;
+      psfShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY];
+      psfShape.sy  = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2;
+      psEllipseAxes psfAxes = psEllipseShapeToAxes (psfShape, 20.0);
+
+      fprintf (stderr, "--- guess values ---\n");
+      fprintf (stderr, "(x,y): %f, %f  Mxx: %f, Myy: %f, Mxy: %f -> major: %f, minor: %f, theta: %f (%f deg)\n", source->peak->xf, source->peak->yf, source->moments->Mxx, source->moments->Myy, source->moments->Mxy, axes.major, axes.minor, axes.theta, axes.theta*PS_DEG_RAD);
+      fprintf (stderr, "psf: major: %f, minor: %f, theta: %f (%f deg)\n", psfAxes.major, psfAxes.minor, psfAxes.theta, psfAxes.theta*PS_DEG_RAD);
+      for (int i = 0; i < model->params->n; i++) {
+        fprintf (stderr, "par %d: %f\n", i, model->params->data.F32[i]);
+      }
+    }
+    if (psTraceGetLevel ("psphot") >= 7) {
+      psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
+    }
+# endif
+
     // for sersic models, use a grid search to choose an index, then float the params there
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
@@ -468 +498 @@
     }
 
+    options.mode = PM_SOURCE_FIT_EXT;
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
         options.mode = PM_SOURCE_FIT_NO_INDEX;
-    } else {
-        options.mode = PM_SOURCE_FIT_EXT;
-    }
-
-    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
+    }
+    if (modelType == pmModelClassGetType("PS_MODEL_TRAIL")) {
+        options.mode = PM_SOURCE_FIT_TRAIL;
+    }
+
     pmSourceFitModel (source, model, &options, maskVal);
-    // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
-    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
+
+# if (PS_TRACE_ON) 
+    if (psTraceGetLevel ("psphot") >= 5) {
+      if (psTraceGetLevel ("psphot") >= 6) {
+        fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
+      }
+      fprintf (stderr, "--- fitted values ---\n");
+      for (int i = 0; i < model->params->n; i++) {
+        fprintf (stderr, "par %d: %f\n", i, model->params->data.F32[i]);
+      }
+      psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
+    }
+# endif
 
     return (model);

trunk/psphot/src/psphotSourceSize.c

@@ -462 +462 @@
     psImageMaskType maskVal = options->maskVal | options->markVal;
 
+    // in the ppSub context, do we get sensible values for ApResid?
+    float ApResidPSF = options->ApResid;
+    if (!isfinite(ApResidPSF)) {
+      ApResidPSF = 0.0;
+    }
+    float ApSysErrPSF = options->ApSysErr;
+    if (!isfinite(ApSysErrPSF)) {
+      ApSysErrPSF = 0.0;
+    }
+
     for (psS32 i = 0 ; i < sources->n ; i++) {
 
@@ -529 +539 @@
         // set nSigmaMAG to include both systematic and poisson error terms.  we include a hard
         // floor on the Ap Sys Err (to be a bit generous).  XXX put the floor in the recipe...
-        float nSigmaMAG = (dMag - options->ApResid) / hypot(source->psfMagErr, hypot(options->ApSysErr, 0.02));
+        float nSigmaMAG = (dMag - ApResidPSF) / hypot(source->psfMagErr, hypot(ApSysErrPSF, 0.02));
         source->extNsigma = nSigmaMAG;
 

@@ -23 +23 @@
 psphotStack
 psphotModelTest
+psphotMinimal