Changeset 36148 for branches/eam_branches/ipp-20130904/psphot/src/psphotGalaxyShape.c

Timestamp:

Sep 24, 2013, 12:09:27 PM (13 years ago)

Author:

eugene

Message:

first version of psphotFullForce which builds

Location:

branches/eam_branches/ipp-20130904/psphot/src

Files:

• . (modified) (1 prop)
• psphotGalaxyShape.c (modified) (13 diffs)

branches/eam_branches/ipp-20130904/psphot/src

@@ -24 +24 @@
 psphotModelTest
 psphotMinimal
+psphotFullForce

@@ -24 +24 @@
 psphotModelTest
 psphotMinimal
+psphotFullForce

branches/eam_branches/ipp-20130904/psphot/src/psphotGalaxyShape.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-bool psphotGalaxyShape (pmConfig *config, const pmFPAview *view, const char *filerule, int pass)
+bool psphotGalaxyShape (pmConfig *config, const pmFPAview *view, const char *filerule)
 {
     bool status = true;
@@ -38 +38 @@
         // psAssert (psf, "missing psf?");
 
-        if (!psphotGalaxyShapeReadout (config, recipe, view, filerule, readout, sources, psf, i, pass)) {
+        if (!psphotGalaxyShapeReadout (config, recipe, view, filerule, readout, sources, psf)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i);
             return false;
@@ -46 +46 @@
 }
 
-bool psphotGalaxyShapeReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index, int pass) {
+bool psphotGalaxyShapeReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf) {
 
     bool status = false;
@@ -55 +55 @@
     }
 
-    psTimerStart ("psphot.kron");
+    psTimerStart ("psphot.galaxy");
 
     // determine the number of allowed threads
@@ -63 +63 @@
     }
 
-    bool KRON_APPLY_WEIGHT = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WEIGHT");
-    if (!status) {
-        KRON_APPLY_WEIGHT = true;
-    }
-
-    bool KRON_APPLY_WINDOW = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WINDOW");
-    if (!status) {
-        KRON_APPLY_WINDOW = false;
-    }
-
     // bit-masks to test for good/bad pixels
     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
@@ -83 +73 @@
     // maskVal is used to test for rejected pixels, and must include markVal
     maskVal |= markVal;
+
+    // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
+    int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
+    assert (status);
+
+    // source fitting parameters for extended source fits
+    int fitIter = psMetadataLookupS32(&status, recipe, "EXT_FIT_ITER"); // Max number of fit iterations
+    assert (status && fitIter > 0);
+
+    float fitMinTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MIN_TOL"); // Fit tolerance
+    if (!status || !isfinite(fitMinTol) || fitMinTol <= 0) {
+        fitMinTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_TOL"); // Fit tolerance
+        if (!status || !isfinite(fitMinTol) || fitMinTol <= 0) {
+            psAbort("PSF_FIT_MIN_TOL (and PSF_FIT_TOL) not defined or positive");
+        }
+    }
+
+    float fitMaxTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_TOL"); // Fit tolerance
+    if (!status || !isfinite(fitMaxTol) || fitMaxTol <= 0) {
+        fitMaxTol = 1.0;
+    }
+
+    bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance
+    if (!status) {
+        // default to the old method (chisqConvergence)
+        chisqConvergence = true;
+    }
+
+    int gainFactorMode = psMetadataLookupS32 (&status, recipe, "LMM_FIT_GAIN_FACTOR_MODE"); // Fit tolerance
+    if (!status) {
+        // default to the old method (chisqConvergence)
+        gainFactorMode = 0;
+    }
+
+    // perform full extended source non-linear fits?
+    bool isInteractive = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS_INTERACTIVE");
+    if (!status) isInteractive = false;
+
+    // Define source fitting parameters for extended source fits
+    pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
+    fitOptions->mode           = PM_SOURCE_FIT_EXT_AND_SKY;
+    fitOptions->saveCovariance = true;  // XXX make this a user option?
+    fitOptions->covarFactor    = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
+    fitOptions->nIter          = fitIter;
+    fitOptions->minTol         = fitMinTol;
+    fitOptions->maxTol         = fitMaxTol;
+
+    fitOptions->gainFactorMode   = gainFactorMode;
+    fitOptions->chisqConvergence = chisqConvergence;
+    fitOptions->isInteractive    = isInteractive;
 
     // source analysis is done in S/N order (brightest first)
@@ -105 +145 @@
 
             // allocate a job -- if threads are not defined, this just runs the job
-            psThreadJob *job = psThreadJobAlloc ("PSPHOT_KRON_ITERATE");
+            psThreadJob *job = psThreadJobAlloc ("PSPHOT_GALAXY_SHAPES");
 
             psArrayAdd(job->args, 1, readout);
             psArrayAdd(job->args, 1, cells->data[j]); // sources
-            PS_ARRAY_ADD_SCALAR(job->args, markVal,            PS_TYPE_IMAGE_MASK);
-            PS_ARRAY_ADD_SCALAR(job->args, maskVal,            PS_TYPE_IMAGE_MASK);
-            PS_ARRAY_ADD_SCALAR(job->args, pass,               PS_TYPE_S32);
+            psArrayAdd(job->args, 1, fitOptions);
+            PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK);
+            PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
+            PS_ARRAY_ADD_SCALAR(job->args, psfSize, PS_TYPE_S32);
 
 // set this to 0 to run without threading
@@ -146 +187 @@
     psFree (cellGroups);
 
-    psLogMsg ("psphot.kron", PS_LOG_WARN, "measure kron fluxes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n);
+    psLogMsg ("psphot.galaxy", PS_LOG_WARN, "measure galaxy shapes : %f sec for %ld objects\n", psTimerMark ("psphot.galaxy"), sources->n);
     return true;
 }
@@ -152 +193 @@
 bool psphotGalaxyShape_Threaded (psThreadJob *job) {
 
-    pmReadout *readout              = job->args->data[0];
-    psArray *sources                = job->args->data[2];
-    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
-    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
+    pmReadout *readout             = job->args->data[0];
+    psArray *sources               = job->args->data[1];
+    pmSourceFitOptions *fitOptions = job->args->data[2];
+    psImageMaskType markVal        = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType maskVal        = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
+    int psfSize                    = PS_SCALAR_VALUE(job->args->data[5],S32);
 
     for (int i = 0; i < sources->n; i++) {
-
-        // XXX what do i need for the kron flux (given a supplied kron radius?)
 
         pmSource *source = sources->data[i];
@@ -184 +225 @@
 
         // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
-        bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
+        pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
         psAssert (source->pixels, "WTF?");
 
         // this function populates moments->Mrf,GalaxyShape,GalaxyShapeErr
         // do the following for a set of shapes (Ex,Ey)
-        psphotGalaxyShapeSource (source, maskVal);
+        psphotGalaxyShapeGrid (source, fitOptions, maskVal, psfSize);
         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
 
@@ -198 +239 @@
 }
 
-bool psphotGalaxyShapeGrid (pmSource *source, psImageMaskType maskVal, float psfSize) {
-
-    pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize);
-
-  // I have some source guess (e0, e1, e2)
-  psEllipsePol guessPol;
-
-  // I am going to retain the value of e0, and grid search around e1,e2
-  // I need to study a valid range (and distribution?) for e1,e2)
-
-  for (float de1 = -0.2; de1 <= +0.2; de1 += 0.1) {
-    for (float de2 = -0.2; de2 <= +0.2; de2 += 0.1) {
+bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, int psfSize) {
+
+    // allocate the model
+    // XXX need to figure out how to get the model into the program
+    pmModelType modelType = pmModelClassGetType ("PS_MODEL_EXP");
+
+    pmModel *model = pmModelAlloc(modelType);
+    if (!model) {
+        return NULL;
+    }
+
+    pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize);
+
+    psF32 *PAR = pcm->modelConv->params->data.F32;
+
+    // I have some source guess (e0, e1, e2)
+    psEllipsePol guessPol;
+
+    // I am going to retain the value of e0, and grid search around e1,e2
+    // I need to study a valid range (and distribution?) for e1,e2)
+
+    for (float de1 = -0.2; de1 <= +0.2; de1 += 0.1) {
+        for (float de2 = -0.2; de2 <= +0.2; de2 += 0.1) {
   
-      psEllipsePol testPol = guessPol;
-      testPol.e1 += de1;
-      testPol.e2 += de2;
-
-      psEllipseShape shape;
-      psEllipsePolToShape (&shape, testPol);
-      
-      PAR[PM_PAR_SXX] = shape.sxx;
-      PAR[PM_PAR_SXY] = shape.sxy;
-      PAR[PM_PAR_SYY] = shape.syy;
-
-      // where do the results go??
-      psphotGalaxyShapeSource (pcm, source, maskVal);
-    }
-  }
-
-  return true;
+            psEllipsePol testPol = guessPol;
+            testPol.e1 += de1;
+            testPol.e2 += de2;
+
+            psEllipseAxes axes = psEllipsePolToAxes (testPol, 0.1);
+
+            pmPSF_AxesToModel (PAR, axes, modelType);
+
+            // where do the results go??
+            psphotGalaxyShapeSource (pcm, source, maskVal, psfSize);
+        }
+    }
+
+    return true;
 }
 
 // fit the given model to the source and find chisq & normalization
 // XXX is this a single-component model? sersic with a supplied index, Reff, axis ratio, and theta?
-bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal) {
+bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, int psfSize) {
 
     PS_ASSERT_PTR_NON_NULL(source, false);
     PS_ASSERT_PTR_NON_NULL(source->peak, false);
     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-
-    psF32 *PAR = pcm->modelConv->params->data.F32;
 
     // generate the modelFlux
@@ -249 +296 @@
 
     for (int iy = 0; iy < source->pixels->numRows; iy++) {
-      for (int ix = 0; ix < source->pixels->numCols; ix++) {
-        // skip masked points
-        if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
-          continue;
+        for (int ix = 0; ix < source->pixels->numCols; ix++) {
+            // skip masked points
+            if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
+                continue;
+            }
+            // skip zero-variance points
+            if (source->variance->data.F32[iy][ix] == 0) {
+                continue;
+            }
+            // skip nan value points
+            if (!isfinite(source->pixels->data.F32[iy][ix])) {
+                continue;
+            }
+
+            float fy = source->pixels->data.F32[iy][ix];
+            float fm = source->modelFlux->data.F32[iy][ix];
+            float wt = (usePoisson) ? 1.0 / source->variance->data.F32[iy][ix] : 1.0;
+
+            YY += PS_SQR(fy) * wt;
+            YM += fm * fy * wt;
+            MM += PS_SQR(fm) * wt;
         }
-        // skip zero-variance points
-        if (source->variance->data.F32[iy][ix] == 0) {
-          continue;
-        }
-        // skip nan value points
-        if (!isfinite(source->pixels->data.F32[iy][ix])) {
-          continue;
-        }
-
-        float fy = source->pixels->data.F32[iy][ix];
-        float fm = source->modelFlux->data.F32[iy][ix];
-        float wt = (usePoisson) ? 1.0 / source->variance->data.F32[iy][ix] : 1.0;
-
-        YY += PS_SQR(fy) * wt;
-        YM += fm * fy * wt;
-        MM += PS_SQR(fm) * wt;
-      }
     }
 
@@ -276 +323 @@
     float Chisq = YY - 2 * Io * YM + Io * Io * MM;
 
-    return true;
-}
+    fprintf (stderr, "Io: %f, Chisq: %f\n", Io, Chisq);
+
+    return true;
+}

@@ -24 +24 @@
 psphotModelTest
 psphotMinimal
+psphotFullForce