Changeset 36154

Timestamp:

Sep 27, 2013, 3:47:04 PM (13 years ago)

Author:

eugene

Message:

updates for psphotFullForce I/O

Location:

branches/eam_branches/ipp-20130904

Files:

• psModules/src/objects/pmSourceIO_CFF.c (modified) (6 diffs)
• psphot/src/psphotGalaxyShape.c (modified) (3 diffs)
• psphot/src/psphotMergeSources.c (modified) (2 diffs)
• psphot/src/psphotPetroFlux.c (modified) (3 diffs)

branches/eam_branches/ipp-20130904/psModules/src/objects/pmSourceIO_CFF.c

@@ -72 +72 @@
     assert (modelType > -1);
 
+    // how to deal with this?
+    int galaxyModelType = pmModelClassGetType ("PS_MODEL_EXP");
+
     // We get the size of the table, and allocate the array of sources first because the table
     // is large and ephemeral --- when the table gets blown away, whatever is allocated after
@@ -87 +90 @@
         }
 
-        // here are the things we read from the FITS table
-
-        unsigned int ID = psMetadataLookupU32 (&status, row, "IPP_IDET");
-        float X = psMetadataLookupF32 (&status, row, "X_PSF");
-        float Y = psMetadataLookupF32 (&status, row, "Y_PSF");
+        // here are the things we read from the FITS table (XXX modify names if needed)
+
+        unsigned int ID  = psMetadataLookupU32 (&status, row, "IPP_IDET");
+        float X          = psMetadataLookupF32 (&status, row, "X_PSF");
+        float Y          = psMetadataLookupF32 (&status, row, "Y_PSF");
 
         float apRadius   = psMetadataLookupS32 (&status, row, "AP_RADIUS");
@@ -97 +100 @@
         float petRadius  = psMetadataLookupU32 (&status, row, "PETRO_RADIUS");
         bool fitGalaxy   = psMetadataLookupU32 (&status, row, "FIT_GALAXY");
+        bool psfStar     = psMetadataLookupU32 (&status, row, "PSF_STAR");
 
         float Rmajor     = psMetadataLookupU32 (&status, row, "R_MAJOR");
@@ -108 +112 @@
         source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags
 
+        // XXX we can set this in general, but for a specific image, we need to weed out SATSTARS
+        if (psfStar) {
+            source->tmpFlags |= PM_SOURCE_TMPF_CANDIDATE_PSFSTAR;
+        }
+
         // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models.
         psF32 *PAR = model->params->data.F32;
         psF32 *dPAR = model->dparams->data.F32;
 
-        source->seq          = ID;
-        PAR[PM_PAR_XPOS]     = X;
-        PAR[PM_PAR_YPOS]     = Y;
+        source->seq       = ID;
+        PAR[PM_PAR_XPOS]  = X;
+        PAR[PM_PAR_YPOS]  = Y;
 
         dPAR[PM_PAR_XPOS] = 0.0;
@@ -129 +138 @@
 
         source->psfMag    = 0.0;
-        source->psfMagErr    = 0.0;
+        source->psfMagErr = 0.0;
         source->apMag     = 0.0;
-
+        source->apRadius  = apRadius;
+
+        // we generate a somewhat fake PSF model here -- 
+        // in most (all?) contexts, we will replace this with a measured psf model
+        // elsewhere
         axes.major        = 1.0;
         axes.minor        = 1.0;
@@ -138 +151 @@
 
         float peakFlux    = 1.0;
-
-        source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
+        source->peak      = pmPeakAlloc(X, Y, peakFlux, PM_PEAK_LONE);
+        source->peak->xf  = X; // more accurate position
+        source->peak->yf  = Y; // more accurate position
+        source->peak->dx  = 0.0;
+        source->peak->dy  = 0.0;
         source->peak->rawFlux = peakFlux;
         source->peak->smoothFlux = peakFlux;
-        source->peak->xf   = PAR[PM_PAR_XPOS]; // more accurate position
-        source->peak->yf   = PAR[PM_PAR_YPOS]; // more accurate position
-        source->peak->dx   = dPAR[PM_PAR_XPOS];
-        source->peak->dy   = dPAR[PM_PAR_YPOS];
-
-        assert (status);
+
+        source->moments = pmMomentsAlloc ();
+        source->moments->Mx = X;
+        source->moments->My = Y;
+        source->moments->Mrf = kronRadius;
+        source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED;
+        // XXX does the above cause a problem with pmSourceMoments (do we run it?)
+
+        if (isfinite(petRadius)) {
+            source->extpars = pmSourceExtendedParsAlloc ();
+            source->extpars->petrosianRadius = petRadius;
+        }
+
+        if (fitGalaxy) {
+            source->modelFits = psArrayAllocEmpty (models->list->n);
+            pmModel *model = pmModelAlloc(modelType);
+            psF32 *xPAR = model->params->data.F32;
+
+            xPAR[PM_PAR_SKY]  = 0.0;
+            xPAR[PM_PAR_I0]   = 1.0;
+            xPAR[PM_PAR_XPOS] = X;
+            xPAR[PM_PAR_YPOS] = Y;
+            
+            psEllipseAxes galAxes;
+            galAxes.major = Rmajor;
+            galAxes.minor = Rminor;
+            galAxes.theta = theta; // XXX degrees or radians?
+
+            pmPSF_AxesToModel (PAR, guessAxes, galModelType);
+            xPAR[PM_PAR_7] = Sindex;
+
+            psArrayAdd (source->modelFits, 1, model);
+            psFree (model);
+            
+        }
 
         sources->data[i] = source;

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

@@ -210 +210 @@
         if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
 
-        // XXX where are we storing the supplied kron radius?
+        // psphotSetRadiusMomentsExact sets the radius based on Mrf
         if (!isfinite(source->moments->Mrf)) continue;
 
+        // modelFits is allocated if a galaxy fit is requested
+        if (!source->modelFits) continue;
+
+        psphotSetRadiusMomentsExact(&fitRadius, &windowRadius, readout, source, markVal); // NOTE : 6 allocs
+
         // skip saturated stars modeled with a radial profile 
+        // XXX worry about this at some point..
         if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
 
@@ -221 +227 @@
         }
 
-        // XXX this is not right (use same concept as general source fit)
-        float windowRadius = 6.0*source->moments->Mrf ;
-
-        // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
-        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)
@@ -245 +244 @@
     pmModelType modelType = pmModelClassGetType ("PS_MODEL_EXP");
 
-    pmModel *model = pmModelAlloc(modelType);
-    if (!model) {
-        return NULL;
-    }
-
+    pmModel *model = source->modelFits->data[0];
+    if (!model) return false;
+
+    // we are using fitOptions->mode : be sure this makes sense
     pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize);
 
+    // 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)
-    psEllipsePol guessPol;
+    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

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

@@ -105 +105 @@
     {
         pmReadout *readoutCMF = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT.CMF");
-        if (!readoutCMF) goto loadTXT;
+        if (!readoutCMF) goto loadCFF;
 
         extCMF = psMetadataLookupPtr (NULL, readoutCMF->analysis, "PSPHOT.DETECTIONS");
@@ -111 +111 @@
             for (int i = 0; i < extCMF->allSources->n; i++) {
                 pmSource *source = extCMF->allSources->data[i];
+                source->mode |= PM_SOURCE_MODE_EXTERNAL;
+
+                // the supplied peak flux needs to be re-normalized
+                source->peak->rawFlux = 1.0;
+                source->peak->smoothFlux = 1.0;
+                source->peak->detValue = 1.0;
+
+                // drop the loaded source modelPSF
+                psFree (source->modelPSF);
+                source->modelPSF = NULL;
+                source->imageID = index;
+
+                psArrayAdd (detections->newSources, 100, source);
+            }
+        }
+    }
+
+loadCFF:
+    // load data from input CFF file:
+    {
+        pmReadout *readoutCFF = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT.CFF");
+        if (!readoutCFF) goto loadTXT;
+
+        extCFF = psMetadataLookupPtr (NULL, readoutCFF->analysis, "PSPHOT.DETECTIONS");
+        if (extCFF) {
+            for (int i = 0; i < extCFF->allSources->n; i++) {
+                pmSource *source = extCFF->allSources->data[i];
                 source->mode |= PM_SOURCE_MODE_EXTERNAL;
 

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

@@ -187 +187 @@
     PS_ASSERT_PTR_NON_NULL(source->peak, false);
     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
+    PS_ASSERT_PTR_NON_NULL(source->variance, false);
+    PS_ASSERT_PTR_NON_NULL(source->moments, false);
+    PS_ASSERT_PTR_NON_NULL(source->extpars, false);
 
     // the peak position is less accurate but less subject to extreme deviations
@@ -202 +205 @@
     // Calculate the Petro magnitude (make this block optional?)
     // XXX set the aperture here
-    float radPetro  = 2.0*source->moments->Mrf;
+    float radPetro  = PETROSIAN_RADII*source->extpars->petrosianRadius;
     float radPetro2 = radPetro*radPetro;
 
@@ -242 +245 @@
 
     // return the flux and error to parameters?
-    // source->moments->PetroFlux    = Sum;
-    // source->moments->PetroFluxErr = sqrt(Var);
+    source->extpars->petrosianFlux    = Sum;
+    source->extpars->petrosianFluxErr = sqrt(Var);
+    source->extpars->petrosianFill = nPetroPix / (M_PI * radPetro2);
 
     fprintf (stderr, "petro flux: %f +/- %f\n", Sum, sqrt(Var));