Changeset 36828 for branches/eam_branches/ipp-20140423/psphot/src/psphotGalaxyShape.c

Timestamp:

Jun 7, 2014, 5:38:26 AM (12 years ago)

Author:

eugene

Message:

merge changes from trunk

Location:

branches/eam_branches/ipp-20140423/psphot/src

Files:

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

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

@@ -80 +80 @@
     assert (markVal);
 
-    // maskVal is used to test for rejected pixels, and must include markVal
     maskVal |= markVal;
 
-    // bit-mask to mark pixels not used in analysis
     psphotGalaxyShapeOptions *opt = psphotGalaxyShapeOptionsAlloc();
+    opt->Q = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_Q"); 
+    psAssert (status, "missing GALAXY_SHAPES_Q");
+    opt->NSigma = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_NSIGMA"); 
+    psAssert (status, "missing GALAXY_SHAPES_NSIGMA");
+    opt->clampSN = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_CLAMP_SN"); 
+    psAssert (status, "missing GALAXY_SHAPES_NSIGMA");
+    psString modelTypeToSave = psMetadataLookupStr(&status, recipe, "EXT_MODEL_TYPE_FORCE");
+    if (modelTypeToSave && strcmp(modelTypeToSave, "ALL") && strcmp(modelTypeToSave, "BEST")) {
+        opt->extModelType = pmModelClassGetType(modelTypeToSave);
+    } else {
+        opt->extModelType = pmModelClassGetType("PS_MODEL_SERSIC");
+    }
+
+#ifdef notdef
     opt->fRmajorMin = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MAJOR_MIN"); if (!status) opt->fRmajorMin = 0.5;
     opt->fRmajorMax = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MAJOR_MAX"); if (!status) opt->fRmajorMax = 2.0;
@@ -91 +103 @@
     opt->fRminorMax = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MINOR_MAX"); if (!status) opt->fRminorMax = 2.0;
     opt->fRminorDel = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MINOR_DEL"); if (!status) opt->fRminorDel = 0.1;
+#endif
 
     // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
@@ -238 +251 @@
 bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psphotGalaxyShapeOptions *opt, psImageMaskType maskVal, int psfSize) {
 
-    // we use the 0th model to define the initial guess shape
-    pmModel *model = source->modelFits->data[0];
-    if (!model) return false;
-
-    pmModelType modelType = model->type;
-
-    // we are using fitOptions->mode : be sure this makes sense
-    pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize);
-    if (!pcm) return false;
-
-    // we are fitting only PM_PAR_I0; the shape elements are generated from a grid
-    psF32 *PAR = pcm->modelConv->params->data.F32;
-
-    // double check that the guess is carried along...
-
-    // I have some source guess (e0, e1, e2)
-    psEllipseAxes guessAxes = pmPSF_ModelToAxes (PAR, modelType);
-
-    if (!source->galaxyFits) {
-        source->galaxyFits = pmSourceGalaxyFitsAlloc();
-    }
-
-    float fRmajorBest = NAN;
-    float fRminorBest = NAN;
-    float chisqBest = NAN;
-    for (float fRmajor = opt->fRmajorMin; fRmajor < opt->fRmajorMax + 0.5*opt->fRmajorDel; fRmajor += opt->fRmajorDel) {
-        for (float fRminor = opt->fRminorMin; fRminor < opt->fRminorMax + 0.5*opt->fRminorDel; fRminor += opt->fRminorDel) {
-  
-            psEllipseAxes testAxes = guessAxes;
-            testAxes.major = guessAxes.major * fRmajor;
-            testAxes.minor = guessAxes.minor * fRminor;
-            
-            pmPSF_AxesToModel (PAR, testAxes, modelType);
-            
-            psphotGalaxyShapeSource (pcm, source, maskVal, psfSize, true);
-
-            int i = source->galaxyFits->chisq->n - 1;
-            float flux = source->galaxyFits->Flux->data.F32[i];
-            if (isfinite(flux)) {
-                float thisChisq = source->galaxyFits->chisq->data.F32[i];
-                if (isfinite(thisChisq) && isfinite(flux) && (!isfinite(chisqBest) || thisChisq < chisqBest)) {
-                    chisqBest = thisChisq;
-                    fRmajorBest = fRmajor;
-                    fRminorBest = fRminor;
+
+    for (int iModel = 0 ; iModel < source->modelFits->n; iModel++) {
+        pmModel *model = source->modelFits->data[iModel];
+        if (!model) return false;
+
+        pmModelType modelType = model->type;
+        model->flags = PM_MODEL_STATUS_NONE;
+
+        // we are using fitOptions->mode : be sure this makes sense
+        pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize);
+        if (!pcm) return false;
+
+        // we are fitting only PM_PAR_I0; the shape elements are generated from a grid
+        psF32 *PAR = pcm->modelConv->params->data.F32;
+
+        pmSourceGalaxyFits *galaxyFits = pmSourceGalaxyFitsAlloc();
+
+        // set the options for this source
+        if (!psphotGalaxyShapeOptionsSet(source, galaxyFits, opt)) {
+            // this source won't work
+            psFree(pcm);
+            psFree(galaxyFits);
+            continue;
+        }
+
+        if (!source->galaxyFits) {
+            source->galaxyFits = psArrayAllocEmpty(1);
+        }
+        psArrayAdd(source->galaxyFits, 1, galaxyFits);
+        psFree(galaxyFits);
+
+        // I have some source guess (e0, e1, e2)
+        psEllipseAxes guessAxes = pmPSF_ModelToAxes (PAR, modelType);
+
+        float fRmajorBest = NAN;
+        float fRminorBest = NAN;
+        float chisqBest = NAN;
+        for (float fRmajor = galaxyFits->fRmajorMin; fRmajor < galaxyFits->fRmajorMax + 0.5*galaxyFits->fRmajorDel; fRmajor += galaxyFits->fRmajorDel) {
+            for (float fRminor = galaxyFits->fRminorMin; fRminor < galaxyFits->fRminorMax + 0.5*galaxyFits->fRminorDel; fRminor += galaxyFits->fRminorDel) {
+      
+                psEllipseAxes testAxes = guessAxes;
+                testAxes.major = guessAxes.major * fRmajor;
+                testAxes.minor = guessAxes.minor * fRminor;
+                
+                pmPSF_AxesToModel (PAR, testAxes, modelType);
+                
+                psphotGalaxyShapeSource (pcm, source, galaxyFits, maskVal, psfSize, true);
+
+                int i = galaxyFits->chisq->n - 1;
+                float flux = galaxyFits->Flux->data.F32[i];
+                if (isfinite(flux)) {
+                    float thisChisq = galaxyFits->chisq->data.F32[i];
+                    if (isfinite(thisChisq) && isfinite(flux) && (!isfinite(chisqBest) || thisChisq < chisqBest)) {
+                        chisqBest = thisChisq;
+                        fRmajorBest = fRmajor;
+                        fRminorBest = fRminor;
+                    }
                 }
+                // reset I0 to avoid potential problems on the next iteration
+                PAR[PM_PAR_I0] = 1.0;
             }
-            // reset I0 to avoid potential problems on the next iteration
-            PAR[PM_PAR_I0] = 1.0;
-        }
-    }
-
-#define SAVE_BEST_MODEL
-#ifdef SAVE_BEST_MODEL
-    // Save model with smallest chisq
-    if (isfinite(chisqBest)) {
-#else 
-    // Save model with nominal parameters
-    {
-        fRmajorBest = 1;
-        fRminorBest = 1;
+        }
+
+        if (isfinite(chisqBest)) {
+            // now save the best fitting model as the source's extended model ...
+            psEllipseAxes testAxes = guessAxes;
+
+            // ... unless this macro is defined
+#ifndef SAVE_NOMINAL_MODEL
+            testAxes.major = guessAxes.major * fRmajorBest;
+            testAxes.minor = guessAxes.minor * fRminorBest;
 #endif
-        // now save the best fitting model as the source's extended model
-        psEllipseAxes testAxes = guessAxes;
-        testAxes.major = guessAxes.major * fRmajorBest;
-        testAxes.minor = guessAxes.minor * fRminorBest;
-        
-        pmPSF_AxesToModel (PAR, testAxes, modelType);
-        
-        psphotGalaxyShapeSource (pcm, source, maskVal, psfSize, false);
-
-        // Replace modelEXT with this model, if the model is good
-        if (isfinite(PAR[PM_PAR_I0])) {
-            psFree (source->modelEXT);
-
-            source->modelEXT = psMemIncrRefCounter (pcm->modelConv);
-            source->type = PM_SOURCE_TYPE_EXTENDED;
-            source->mode |= PM_SOURCE_MODE_EXTMODEL;
-            source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
-
-            // cache the model flux
-            pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma);
-        }
-    }
-
-    psFree (pcm);
+            
+            pmPSF_AxesToModel (PAR, testAxes, modelType);
+                
+            psphotGalaxyShapeSource (pcm, source, galaxyFits, maskVal, psfSize, false);
+
+            // Replace modelEXT with the best model from the first of the model fits if one of them is good
+            if (isfinite(PAR[PM_PAR_I0]) && modelType == opt->extModelType) {
+                psFree (source->modelEXT);
+
+                source->modelEXT = psMemIncrRefCounter (pcm->modelConv);
+                source->type = PM_SOURCE_TYPE_EXTENDED;
+                source->mode |= PM_SOURCE_MODE_EXTMODEL;
+                source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
+
+                // cache the model flux
+                pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma);
+            }
+        }
+
+        psFree (pcm);
+    }
     return true;
 }
@@ -328 +349 @@
 // 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, int psfSize, bool saveResults) {
+bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, pmSourceGalaxyFits *galaxyFits, psImageMaskType maskVal, int psfSize, bool saveResults) {
 
     PS_ASSERT_PTR_NON_NULL(source, false);
@@ -384 +405 @@
 
     if (saveResults) {
-        psVectorAppend (source->galaxyFits->Flux, flux);
-        psVectorAppend (source->galaxyFits->dFlux, dflux);
-        psVectorAppend (source->galaxyFits->chisq, Chisq);
-        source->galaxyFits->nPix = nPix;
+        psVectorAppend (galaxyFits->Flux, flux);
+        psVectorAppend (galaxyFits->dFlux, dflux);
+        psVectorAppend (galaxyFits->chisq, Chisq);
+        galaxyFits->nPix = nPix;
+        galaxyFits->modelType = pcm->modelConv->type;
     }
 
@@ -404 +426 @@
     psMemSetDeallocator(opt, (psFreeFunc) psphotGalaxyShapeOptionsFree);
     
-    opt->fRmajorMin = 0.5;
-    opt->fRmajorMax = 2.0;
-    opt->fRmajorDel = 0.1;
-    opt->fRminorMin = 0.5;
-    opt->fRminorMax = 2.0;
-    opt->fRminorDel = 0.1;
-
     return opt;
 }
 
+bool psphotGalaxyShapeOptionsSet(pmSource *source, pmSourceGalaxyFits *galaxyFits, psphotGalaxyShapeOptions *options)
+{
+    // XXX: put these in recipe
+    // doesn't make sense to use too small of a spacing
+    float clampedSN = source->extSN < options->clampSN ? source->extSN : options->clampSN;
+
+    float f_del = options->Q / clampedSN;
+    float f_min = 1 - options->NSigma * f_del;
+    float f_max = 1 + options->NSigma * f_del;
+
+    // if f_min goes negative skip this object
+    if (f_min < 0) return false;
+
+    galaxyFits->fRmajorMin = f_min;
+    galaxyFits->fRmajorMax = f_max;
+    galaxyFits->fRmajorDel = f_del;
+    galaxyFits->fRminorMin = f_min;
+    galaxyFits->fRminorMax = f_max;
+    galaxyFits->fRminorDel = f_del;
+
+    return true;
+}
+