Changeset 36218

Timestamp:

Oct 15, 2013, 11:40:28 AM (13 years ago)

Author:

eugene

Message:

galaxy shape grid defined in recipe

Location:

branches/eam_branches/ipp-20130904/psphot/src

Files:

• psphot.h (modified) (2 diffs)
• psphotGalaxyShape.c (modified) (11 diffs)
• psphotSetThreads.c (modified) (1 diff)

branches/eam_branches/ipp-20130904/psphot/src/psphot.h

@@ -421 +421 @@
 } psphotStackOptions;
 
+typedef struct {
+    float fRmajorMin;
+    float fRmajorMax;
+    float fRmajorDel;
+    float fRminorMin;
+    float fRminorMax;
+    float fRminorDel;
+} psphotGalaxyShapeOptions;
+
 /*** psphotStackMatchPSF prototypes ***/
 bool psphotStackMatchPSFs (pmConfig *config, const pmFPAview *view);
@@ -517 +526 @@
 bool psphotGalaxyShapeReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf);
 bool psphotGalaxyShape_Threaded (psThreadJob *job);
-bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, int psfSize);
+bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psphotGalaxyShapeOptions *opt, psImageMaskType maskVal, int psfSize);
 bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, int psfSize);
+psphotGalaxyShapeOptions *psphotGalaxyShapeOptionsAlloc();
 
 bool psphotRadialProfileWings (pmConfig *config, const pmFPAview *view, const char *filerule);

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

@@ -83 +83 @@
     maskVal |= markVal;
 
+    // bit-mask to mark pixels not used in analysis
+    psphotGalaxyShapeOptions *opt = psphotGalaxyShapeOptionsAlloc();
+    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;
+    opt->fRmajorDel = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MAJOR_DEL"); if (!status) opt->fRmajorDel = 0.1;
+    opt->fRminorMin = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MINOR_MIN"); if (!status) opt->fRminorMin = 0.5;
+    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;
+
     // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
+    // NOTE: this is only used if we are NOT smoothing with a 1D Gaussian
     int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
     assert (status);
@@ -119 +129 @@
             psArrayAdd(job->args, 1, cells->data[j]); // sources
             psArrayAdd(job->args, 1, fitOptions);
+            psArrayAdd(job->args, 1, opt);
             PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK);
             PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
@@ -156 +167 @@
     psFree (cellGroups);
     psFree (fitOptions);
+    psFree (opt);
 
     psLogMsg ("psphot.galaxy", PS_LOG_WARN, "measure galaxy shapes : %f sec for %ld objects\n", psTimerMark ("psphot.galaxy"), sources->n);
@@ -166 +178 @@
     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);
+    psphotGalaxyShapeOptions *opt  = job->args->data[3];
+    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);
+    int psfSize                    = PS_SCALAR_VALUE(job->args->data[6],S32);
 
     float fitRadius;
@@ -202 +215 @@
         // this function populates moments->Mrf,GalaxyShape,GalaxyShapeErr
         // do the following for a set of shapes (Ex,Ey)
-        psphotGalaxyShapeGrid (source, fitOptions, maskVal, psfSize);
+        psphotGalaxyShapeGrid (source, fitOptions, opt, maskVal, psfSize);
+
         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
 
@@ -211 +225 @@
 }
 
-bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, int psfSize) {
+bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psphotGalaxyShapeOptions *opt, psImageMaskType maskVal, int psfSize) {
 
     // we use the 0th model to define the initial guess shape
@@ -230 +244 @@
     // I have some source guess (e0, e1, e2)
     psEllipseAxes guessAxes = pmPSF_ModelToAxes (PAR, modelType);
-    psEllipsePol guessPol = psEllipseAxesToPol (guessAxes);
-
-    // 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)
 
     if (!source->galaxyFits) {
-      source->galaxyFits = pmSourceGalaxyFitsAlloc();
-    }
-
-    for (float de1 = -0.2; de1 <= +0.2; de1 += 0.1) {
-        for (float de2 = -0.2; de2 <= +0.2; de2 += 0.1) {
+        source->galaxyFits = pmSourceGalaxyFitsAlloc();
+    }
+
+    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) {
   
-            psEllipsePol testPol = guessPol;
-            testPol.e1 += de1;
-            testPol.e2 += de2;
-
-            psEllipseAxes axes = psEllipsePolToAxes (testPol, 0.1);
-
-            pmPSF_AxesToModel (PAR, axes, modelType);
-
+            psEllipseAxes testAxes = guessAxes;
+            testAxes.major = guessAxes.major * fRmajor;
+            testAxes.minor = guessAxes.minor * fRminor;
+            
+            pmPSF_AxesToModel (PAR, testAxes, modelType);
+            
             // where do the results go??
             psphotGalaxyShapeSource (pcm, source, maskVal, psfSize);
         }
     }
+
     psFree (pcm);
     return true;
@@ -267 +276 @@
 
     // generate the modelFlux
+    // need to reset here each time since we assign below to calculate the flux
     pmPCMMakeModel (source, pcm->modelConv, maskVal, psfSize);
+    pcm->modelConv->params->data.F32[PM_PAR_I0] = 1.0;
         
     int nPix = 0;
@@ -273 +284 @@
     float YM = 0.0;
     float MM = 0.0;
-    bool usePoisson = false;
+    bool usePoisson = true;
 
     for (int iy = 0; iy < source->pixels->numRows; iy++) {
@@ -306 +317 @@
     float Chisq = (YY - 2 * Io * YM + Io * Io * MM) / (float) nPix;
 
-    // fprintf (stderr, "Io: %f, Chisq: %f\n", Io, Chisq);
-
-    psVectorAppend (source->galaxyFits->Io, Io);
-    psVectorAppend (source->galaxyFits->dIo, dIo);
+    pcm->modelConv->params->data.F32[PM_PAR_I0] = Io;
+    float flux = pcm->modelConv->modelFlux (pcm->modelConv->params);
+    float dflux = flux * (dIo / Io);
+
+    psVectorAppend (source->galaxyFits->Flux, flux);
+    psVectorAppend (source->galaxyFits->dFlux, dflux);
     psVectorAppend (source->galaxyFits->chisq, Chisq);
     source->galaxyFits->nPix = nPix;
@@ -315 +328 @@
     return true;
 }
+
+/**** support functions ****/
+
+void psphotGalaxyShapeOptionsFree (psphotGalaxyShapeOptions *opt) {
+    return;
+}
+
+psphotGalaxyShapeOptions *psphotGalaxyShapeOptionsAlloc()
+{
+    psphotGalaxyShapeOptions *opt = (psphotGalaxyShapeOptions *) psAlloc(sizeof(psphotGalaxyShapeOptions));
+    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;
+}
+

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

@@ -57 +57 @@
     psFree(task);
 
-    task = psThreadTaskAlloc("PSPHOT_GALAXY_SHAPES", 6);
+    task = psThreadTaskAlloc("PSPHOT_GALAXY_SHAPES", 7);
     task->function = &psphotGalaxyShape_Threaded;
     psThreadTaskAdd(task);