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

Timestamp:

Nov 8, 2013, 10:41:29 AM (13 years ago)

Author:

bills

Message:

implement creation of CFF files from sources. some tweaks to reading of cffs

File:

• branches/eam_branches/ipp-20130904/psModules/src/objects/pmSourceIO_CFF.c (modified) (5 diffs)

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

@@ -65 +65 @@
     int modelType = pmModelClassGetType ("PS_MODEL_GAUSS");
 
-    char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME");
+    char *PSF_NAME = psMetadataLookupStr (&status, header, "PSFMODEL");
     if (PSF_NAME != NULL) {
         modelType = pmModelClassGetType (PSF_NAME);
@@ -107 +107 @@
         // XXX: we need to put a lookup table in the cff header to define the correspondence of the
         // model type values in the cff with our models. (We want to use an interger for efficiency
-        // but the values depend on the array in pmModelClass.c
+        // but the value for each model type is set on the organization of the the array in pmModelClass.c
+        // For now use the input values verbatim and trust the user that this is valid value
         int   galaxyModelType = psMetadataLookupS32(&status, row, "MODEL_TYPE");
         float Sindex     = psMetadataLookupF32 (&status, row, "INDEX"); // Should this be PAR_07 not sersic index
@@ -173 +174 @@
         }
 
-        if (fitGalaxy) {
+        // XXX: should use < 0 as invalid galaxyModelType
+
+        if (fitGalaxy && galaxyModelType > 0) {
             source->modelFits = psArrayAllocEmpty (1);
             pmModel *model = pmModelAlloc(galaxyModelType);
@@ -192 +195 @@
                 xPAR[PM_PAR_7] = 0.5 / Sindex;
             }
+            // XXX: set source->modelEXT to this model and flag as extended so that we have a better
+            // shot at subtracting extended sources?
 
             psArrayAdd (source->modelFits, 1, model);
@@ -204 +209 @@
 }
 
+bool pmSourcesWrite_CFF(psFits *fits, psArray *sources, psMetadata *header, psMetadata *recipe) {
+    
+    char *extname = "SkyChip.cff";
+    
+    psArray *table = psArrayAllocEmpty(sources->n);
+
+    pmModelType sersicModelType = pmModelClassGetType("PS_MODEL_SERSCIC");
+
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *thisSource = sources->data[i];
+        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
+
+        psF32 xPos, yPos, flux, rMajor, rMinor, theta;
+        psS32 modelType = 0;
+        bool fitGalaxy = false;
+        bool psfStar = false;
+        psS32 sersicIndex = 0;
+        if (source->modelFits == NULL) {
+            pmModel *model = source->modelPSF;
+            if (model == NULL) continue;
+            psF32 *PAR = model->params->data.F32;
+            if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY])  || !isfinite(PAR[PM_PAR_SXY]) ||
+                !isfinite(source->psfFlux)) {
+                continue;
+            }
+
+            xPos = model->params->data.F32[PM_PAR_XPOS];
+            yPos = model->params->data.F32[PM_PAR_YPOS];
+            flux = source->psfFlux;
+            psfStar = (source->mode & PM_SOURCE_MODE_PSFSTAR) ? true : false;
+            rMajor = 0;
+            rMinor = 0;
+            theta = 0;
+        } else {
+            // Choose the extended model to use for this source.
+            // For now we use the one set as the modelEXT for this source which had the best fit.
+            // XXX: make this controllable by recipe
+            //   use best (as implemented here)
+            //   choose specific type
+            int iModel = -1;
+            if (source->modelEXT) {
+                pmModelType ext_model_type = source->modelEXT->type;
+                for (int j=0; j<source->modelFits->n; j++) {
+                    pmModel *aModel = source->modelFits->data[j];
+                    if (aModel->type == ext_model_type) {
+                        iModel = j;
+                        break;
+                    }
+                }
+            }
+            if (iModel == -1) {
+                // Can this happen? perhaps if model type is qgauss (extended but below s/n for model fits?
+                // XXX: Should this be an assert?
+                continue;
+            }
+            pmModel *model = source->modelFits->data[iModel];
+            psF32 *PAR = model->params->data.F32;
+            xPos = PAR[PM_PAR_XPOS];
+            yPos = PAR[PM_PAR_YPOS];
+            if (model->type == pmModelClassGetType("PS_MODEL_TRAIL")) {
+                // XXX: do we need to handle this type
+                continue;
+            } else {
+                if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY])  || !isfinite(PAR[PM_PAR_SXY]) ||
+                    !isfinite(model->mag)) { 
+                    // bad model
+                    continue;
+                }
+                psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type);
+                rMajor = axes.major;
+                rMinor = axes.minor;
+                theta  = axes.theta*PS_DEG_RAD;
+                flux = pow(10.0, -0.4*model->mag);
+                fitGalaxy = true;
+            }
+            modelType = model->type;
+            if (modelType == sersicModelType) {
+                sersicIndex = 0.5 / PAR[PM_PAR_7];
+            }
+        }
+        psMetadata *row = psMetadataAlloc();
+        psMetadataAddU32 (row, PS_LIST_TAIL, "ID",         0,   "IPP detection identifier",  source->seq);
+        psMetadataAddF32 (row, PS_LIST_TAIL, "X",                0, "x coordinate",          xPos);
+        psMetadataAddF32 (row, PS_LIST_TAIL, "Y",                0, "y coordinate",          yPos);
+        psMetadataAddF32 (row, PS_LIST_TAIL, "FLUX",             0, "flux",                  flux);
+        psMetadataAddF32 (row, PS_LIST_TAIL, "AP_RADIUS",        0, "aperture radius",       source->apRadius);
+        psMetadataAddF32 (row, PS_LIST_TAIL, "KRON_RADIUS",      0, "Kron radius",           source->moments->Mrf * 2.5);
+        psMetadataAddF32 (row, PS_LIST_TAIL, "PETRO_RADIUS",     0, "Petrosian Radius",      source->extpars ? source->extpars->petrosianRadius : NAN);
+        psMetadataAddBool (row, PS_LIST_TAIL, "FIT_GALAXY",      0, "source has xfit",       fitGalaxy); 
+        psMetadataAddBool (row, PS_LIST_TAIL, "PSF_STAR",        0, "source was psf star",   psfStar);
+        psMetadataAddF32 (row, PS_LIST_TAIL, "R_MAJOR",          0, "radius of major axis",  rMajor);
+        psMetadataAddF32 (row, PS_LIST_TAIL, "R_MINOR",          0, "radius of minor axis",  rMinor);
+        psMetadataAddF32 (row, PS_LIST_TAIL, "THETA",            0, "theta",                 theta);
+        psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TYPE",       0, "model type",            modelType);
+        psMetadataAddS32 (row, PS_LIST_TAIL, "INDEX",            0, "sersic index",          sersicIndex);
+
+        psArrayAdd(table, 100, row);
+    }
+
+    if (!psFitsWriteTable(fits, NULL, table, extname)) {
+        psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
+        psFree(table);
+        psFree(header);
+        return false;
+    }
+    psFree(table);
+    // psFree(header);
+
+    return true;
+}