Changeset 5643

Timestamp:

Nov 30, 2005, 2:32:16 PM (21 years ago)

Author:

eugene

Message:

adding fixed-psf fitting

Location:

trunk/psphot

Files:

• Makefile (modified) (2 diffs)
• src/pmSourceFitFixed.c (modified) (5 diffs)
• src/psphot.c (modified) (1 diff)
• src/psphot.h (modified) (2 diffs)
• src/psphotFixedPSF.c (modified) (4 diffs)

trunk/psphot/Makefile

@@ -28 +28 @@
 $(SRC)/psphotChoosePSF.$(ARCH).o   \
 $(SRC)/psphotApplyPSF.$(ARCH).o    \
+$(SRC)/psphotFixedPSF.$(ARCH).o    \
 $(SRC)/psphotFitGalaxies.$(ARCH).o \
 $(SRC)/psphotOutput.$(ARCH).o      \
@@ -38 +39 @@
 $(SRC)/psModulesUtils.$(ARCH).o    \
 $(SRC)/pmPeaksSigmaLimit.$(ARCH).o \
+$(SRC)/pmSourceFitFixed.$(ARCH).o  \
 $(SRC)/psImageData.$(ARCH).o
 

trunk/psphot/src/pmSourceFitFixed.c

@@ -1 +1 @@
 # include "psphot.h"
+# define USECLIP 0
 
-bool pmSourceFitFixed (pmSource *source,
-                       pmModel *model)
-{
+bool pmSourceFitFixed (pmSource *source, pmModel *model) {
+
     PS_ASSERT_PTR_NON_NULL(source, false);
     PS_ASSERT_PTR_NON_NULL(source->peak, false);
@@ -11 +11 @@
     PS_ASSERT_PTR_NON_NULL(source->weight, false);
 
-    // XXX EAM : is it necessary for the mask & weight to exist?  the
-    //           tests below could be conditions (!NULL)
-
-    psBool fitStatus = true;
-    psBool onPic     = true;
-    psBool rc        = true;
-    psF32  Ro, ymodel;
+    psF32 Ro, ymodel;
 
     psVector *params = model->params;
     psVector *dparams = model->dparams;
-    psVector *paramMask = NULL;
 
     pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
@@ -52 +45 @@
             yErr->data.F32[count] = sqrt (source->weight->data.F32[i][j] * Ro);
 
-            psVectorExtend (x, 1, 100);
-            psVectorExtend (y, 1, 100);
-            psVectorExtend (yErr, 1, 100);
+            psVectorExtend (x,    100, 1);
+            psVectorExtend (y,    100, 1);
+            psVectorExtend (yErr, 100, 1);
             count++;
         }
@@ -61 +54 @@
     psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
 
-    // psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
-    // stats->clipIter = 3;
-    // poly = psVectorClipFitPolynomial1D (poly, stats, NULL, 0, y, yErr, x);
+# if (USECLIP) 
+    psVector *mask = psVectorAlloc (count, PS_TYPE_U8);
+    psVectorInit (mask, 0);
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    stats->clipIter = 3;
+    poly = psVectorClipFitPolynomial1D (poly, stats, mask, 0, y, yErr, x);
+# else
+    poly = psVectorFitPolynomial1D (poly, NULL, 0, y, yErr, x);
+# endif
 
-    poly = psVectorFitPolynomial1D (poly, stats, NULL, 0, y, yErr, x);
-
-    params->data.F32[0] = poly->coeff[0];
-    params->data.F32[1] = poly->coeff[1];
+    // fitted values are variations from current values
+    params->data.F32[0] += poly->coeff[0];
+    params->data.F32[1] *= poly->coeff[1];
 
     // XXX need to get errors back from psVectorClipFit...
@@ -82 +80 @@
     psFree(y);
     psFree(yErr);
+
+# if (USECLIP)
     psFree(stats);
+    psFree(mask);
+# endif
 
     return true;

trunk/psphot/src/psphot.c

@@ -49 +49 @@
     if (!strcasecmp (breakPt, "PSFFIT")) exit (0);
 
-    // test PSF on all except SATURATE and DEFECT (artifacts)
-    psphotApplyPSF (imdata, config, sources, psf, sky);
-    if (!strcasecmp (breakPt, "PSFSUB")) exit (0);
+    int FITMODE = psMetadataLookupS32 (&status, config, "FIT_MODE");
+    if (!status) FITMODE = 0;
 
-    // fit extended objects with galaxy models
-    psphotFitGalaxies (imdata, config, sources, sky);
-    if (!strcasecmp (breakPt, "GALSUB")) exit (0);
+    switch (FITMODE) {
+      case 0:
+        psphotFixedPSF (imdata, config, sources, psf, sky);
+        break;
+
+      case 1:
+        // test PSF on all except SATURATE and DEFECT (artifacts)
+        psphotApplyPSF (imdata, config, sources, psf, sky);
+        break;
+
+      case 2:
+      // fit extended objects with galaxy models
+        psphotApplyPSF (imdata, config, sources, psf, sky);
+        psphotFitGalaxies (imdata, config, sources, sky);
+        break;
+
+      case 3:
+      // fit extended objects with galaxy models
+        psphotFixedPSF (imdata, config, sources, psf, sky);
+        psphotFitGalaxies (imdata, config, sources, sky);
+        break;
+    }
 
     // write out data in appropriate format

trunk/psphot/src/psphot.h

@@ -24 +24 @@
 pmPSF       *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *sky);
 bool         psphotApplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
+bool         psphotFixedPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
 bool         psphotFitGalaxies (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats);
 void         psphotOutput (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
@@ -33 +34 @@
 bool         psphotDefinePixels (pmSource *mySource, const eamReadout *imdata, psF32 x, psF32 y, psF32 Radius);
 
+bool         pmSourceFitFixed (pmSource *source, pmModel *model);
 psArray     *pmPeaksSigmaLimit (eamReadout *imdata, psMetadata *config, psStats *sky);
 pmModel     *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius);

trunk/psphot/src/psphotFixedPSF.c

@@ -3 +3 @@
 // fit psf model to all objects with centroid fixed 
 
-bool psphotFitFixedPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) 
+bool psphotFixedPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) 
 { 
     bool  status;
     float x;
     float y;
+    float Sky;
     int   Nfit = 0;
-    int   Nsub = 0;
-    int   Niter = 0;
 
     psTimerStart ("psphot");
 
     // we may set this differently here from the value used to mark likely saturated stars
-    float SATURATION   = psMetadataLookupF32 (&status, config, "SATURATION");
     float OUTER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
 
-    float PSF_MIN_SN       = psMetadataLookupF32 (&status, config, "PSF_MIN_SN");
-    float PSF_MAX_CHI      = psMetadataLookupF32 (&status, config, "PSF_MAX_CHI");
     float PSF_FIT_NSIGMA   = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA");
     float PSF_FIT_PADDING  = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING");
-    float PSF_SHAPE_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA");
 
     // set the object surface-brightness limit for fitted pixels
@@ -45 +40 @@
         // set PSF parameters for this model
         pmModel *model = pmModelFromPSF (modelFLT, psf);
+        // model->params->data.F32[2] = source->peak->x;
+        // model->params->data.F32[3] = source->peak->y;
         x = model->params->data.F32[2];
         y = model->params->data.F32[3];
@@ -69 +66 @@
 
         if (!status || (model->params->data.F32[1] < 0)) {
-          psLogMsg ("psphot", 5, "PSF fit failed for %f, %f (%d iterations)\n", x, y, model->nIter);
+          psLogMsg ("psphot", 5, "PSF fit failed for %f, %f (peak %f\n", x, y, model->params->data.F32[1]);
           source->type = PM_SOURCE_FAIL_FIT_PSF;  // better choice?
           psFree (model);
@@ -77 +74 @@
         // model succeeded : tentatively keep it
         source->modelPSF = model;
-        Niter += model[0].nIter;
         Nfit ++;
 
         // subtract object, leave local sky
-        sky = model->params->data.F32[0];
+        Sky = model->params->data.F32[0];
         model->params->data.F32[0] = 0;
         pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false);
-        model->params->data.F32[0] = sky;
+        model->params->data.F32[0] = Sky;
     }
 
-    psLogMsg ("psphot", 3, "fit fixed PSF models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter);
+    psLogMsg ("psphot", 3, "fit fixed PSF models: %f sec for %d objects\n", psTimerMark ("psphot"), Nfit);
     return (true);
 }