Changeset 36375 for trunk/psModules/src/objects/pmSourceIO_CMF.c.in

Timestamp:

Dec 10, 2013, 2:55:11 PM (12 years ago)

Author:

eugene

Message:

merge changes from eam_branches/ipp-20130904

Location:

trunk

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/objects/pmSourceIO_CMF.c.in (modified) (13 diffs)

trunk/psModules/src/objects/pmSourceIO_CMF.c.in

@@ -789 +789 @@
     char name[64];
 
+    pmModelType modelTypeTrail = pmModelClassGetType("PS_MODEL_TRAIL");
+
     // create a header to hold the output data
     psMetadata *outhead = psMetadataAlloc ();
@@ -798 +800 @@
     sources = psArraySort (sources, pmSourceSortByFlux);
 
-    @>PS1_DV2@ float magOffset; 
-    @>PS1_DV2@ float zeroptErr; 
-    @>PS1_DV2@ float fwhmMajor; 
-    @>PS1_DV2@ float fwhmMinor;
-    @>PS1_DV2@ pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
+    float magOffset; 
+    float zeroptErr; 
+    float fwhmMajor; 
+    float fwhmMinor;
+    pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
 
     // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams
@@ -823 +825 @@
     @>PS1_DV2@ pmChip *chip = readout->parent->parent;
 
+    pmModelStatus badModel = PM_MODEL_STATUS_NONE;
+    badModel |= PM_MODEL_STATUS_BADARGS;
+    badModel |= PM_MODEL_STATUS_OFFIMAGE;
+    badModel |= PM_MODEL_STATUS_NAN_CHISQ;
+    badModel |= PM_MODEL_SERSIC_PCM_FAIL_GUESS;
+    badModel |= PM_MODEL_SERSIC_PCM_FAIL_GRID;
+    badModel |= PM_MODEL_PCM_FAIL_GUESS;
+
     table = psArrayAllocEmpty (sources->n);
 
@@ -847 +857 @@
 
             // skip models which were not actually fitted
-            if (model->flags & PM_MODEL_STATUS_BADARGS) continue;
+            // XXX 
+            if (model->flags & badModel) continue;
 
             PAR = model->params->data.F32;
@@ -853 +864 @@
             xPos = PAR[PM_PAR_XPOS];
             yPos = PAR[PM_PAR_YPOS];
-            xErr = dPAR[PM_PAR_XPOS];
-            yErr = dPAR[PM_PAR_YPOS];
+
+            // for the extended source models, we do not always fit the centroid in the non-linear fitting process
+            // current situation (hard-wired into psphotSourceFits.c:psphotFitPCM,
+            // SERSIC, DEV, EXP : X,Y not fitted (PCM and not PCM)
+            // TRAIL : X,Y are fitted
+            // 
+            
+            // XXX this should be based on what happened, not on the model type
+            if (model->type == modelTypeTrail) {
+                xErr = dPAR[PM_PAR_XPOS];
+                yErr = dPAR[PM_PAR_YPOS];
+            } else {
+                // this is definitely an underestimate since it does not
+                // account for the extent of the source
+                xErr = fwhmMajor * model->magErr / 2.35;
+                yErr = fwhmMinor * model->magErr / 2.35;
+            }
 
             @>PS1_DV2@ psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
@@ -870 +896 @@
             row = psMetadataAlloc ();
 
-            // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
             // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation
             // This set of psMetadataAdd Entries marks the "----" "Start of the XFIT segment"
@@ -888 +913 @@
             @>PS1_DV2@ float calMag = isfinite(magOffset) ? model->mag + magOffset : NAN;
             @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CAL_MAG", PS_DATA_F32, "EXT Magnitude using supplied calibration",   calMag);
-            @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CHISQ",   PS_DATA_F32, "EXT Magnitude using supplied calibration",   model->chisq);
-            @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NDOF",    PS_DATA_S32, "EXT Magnitude using supplied calibration",   model->nDOF);
+            @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CHISQ",   PS_DATA_F32, "EXT Model Chisq",   model->chisq);
+            @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NDOF",    PS_DATA_S32, "EXT Model num degrees of freedom",   model->nDOF);
+            @>PS1_SV1,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_MODEL_TYPE",    PS_DATA_S32, "type for chosen EXT_MODEL",   source->modelEXT ? source->modelEXT->type : -1);
+
+            // EAM : adding for PV2 outputs:
+            @>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_FLAGS", PS_DATA_S16, "model fit flags (pmModelStatus)", source->modelEXT ? source->modelEXT->flags : 0);
 
             @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "POSANGLE",   0, "position angle at source (degrees)",         posAngle);
@@ -946 +975 @@
 
                 snprintf (name, 64, "EXT_PAR_%02d", k);
-
+               
                 if (k < model->params->n) {
                     psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->params->data.F32[k]);
@@ -956 +985 @@
             // optionally, write out the covariance matrix values
             // XXX do I need to pad this to match the biggest covar matrix?
-            if (model->covar) {
+            if (false && model->covar) {
                 for (int iy = 0; iy < model->covar->numCols; iy++) {
                     for (int ix = iy; ix < model->covar->numCols; ix++) {
@@ -1058 +1087 @@
         model->magErr = psMetadataLookupF32(&status, row, "EXT_INST_MAG_SIG");
 
+        model->chisq = psMetadataLookupF32(&status, row, "EXT_CHISQ");
+        model->nDOF = psMetadataLookupF32(&status, row, "EXT_NDOF");
+
+        // EXT_MODEL_TYPE gives the model chosen by psphot as the best.
+        // Putting this into the XFIT table makes 3 copies of it (one for each model)
+        // but since we have many fewer XFIT rows than psf rows that is cheaper than putting it
+        // in the psf table.
+        psS32 extModelType = psMetadataLookupS32(&status, row, "EXT_MODEL_TYPE");
+        if (!status) {
+            // older cmfs don't have this column
+            extModelType = -1;
+        }
+
         psEllipseAxes axes;
         axes.major = psMetadataLookupF32(&status, row, "EXT_WIDTH_MAJ");
@@ -1072 +1114 @@
         if (model->params->n > 7) {
             PAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07");
-        }
-        // read the covariance matrix
-        int nparams = model->params->n;
-        psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32);
-        for (int y = 0; y < nparams; y++) {
-            for (int x = 0; x < nparams; x++) {
-                char name[64];
-                snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x);
-                covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name);
+            // XXX add an error:
+            // dPAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07_");
+        }
+
+        // XXX : make this depend on what is in the cmf
+        if (0) {
+            // read the covariance matrix
+            int nparams = model->params->n;
+            psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32);
+            for (int y = 0; y < nparams; y++) {
+                for (int x = 0; x < nparams; x++) {
+                    char name[64];
+                    snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x);
+                    covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name);
+                }
+            }
+            model->covar = covar;
+        }
+
+        if (modelType == extModelType) {
+            // The software that created this source picked this model as the best of the fits. 
+            // Set the extModel to point to it.
+            // This is important for programs like psastro (skycal) so that its output cmfs
+            // will have valid EXT_MODEL_TYPE
+            psFree(source->modelEXT);
+            source->modelEXT = psMemIncrRefCounter(model);
+            source->type = PM_SOURCE_TYPE_EXTENDED;
+            if (0) {
+                // since FLAGS were read we don't need to do this
+                source->mode |= PM_SOURCE_MODE_EXTMODEL;
+                source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
             }
         }
-        model->covar = covar;
 
         psArrayAdd(source->modelFits, 1, model);
         psFree(model);
-
         psFree(row);
     }
@@ -1209 +1271 @@
 
         write_annuli:
-            psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX",       PS_DATA_VECTOR, "flux within annuli",       radFlux);
-            psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR",   PS_DATA_VECTOR, "flux error in annuli",     radFluxErr);
-            psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_DATA_VECTOR, "flux standard deviation",  radFluxStdev);
-            psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL",       PS_DATA_VECTOR, "fill factor of annuli",    radFill);
+            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX",       PS_META_REPLACE, "flux within annuli",       radFlux);
+            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_ERR",   PS_META_REPLACE, "flux error in annuli",     radFluxErr);
+            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_META_REPLACE, "flux standard deviation",  radFluxStdev);
+            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FILL",       PS_META_REPLACE, "fill factor of annuli",    radFill);
             psFree (radFlux);
             psFree (radFluxErr);
@@ -1343 +1405 @@
     return true;
 }
+
+// XXX where should I record the number of columns??
+bool pmSourcesWrite_CMF_@CMFMODE@_XGAL (psFits *fits, psArray *sources, char *extname, psMetadata *recipe)
+{
+    bool status = false;
+
+    // perform full non-linear fits / extended source analysis?
+    if (!psMetadataLookupBool (&status, recipe, "GALAXY_SHAPES")) {
+        psLogMsg ("psphot", PS_LOG_INFO, "galaxy shapes were not measured, skipping\n");
+        return true;
+    }
+
+    // create a header to hold the output data
+    psMetadata *outhead = psMetadataAlloc ();
+
+    // write the links to the image header
+    psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "galaxy table extension", extname);
+
+    // let's write these out in S/N order
+    sources = psArraySort (sources, pmSourceSortByFlux);
+
+    psArray *table = psArrayAllocEmpty (sources->n);
+
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *thisSource = sources->data[i];
+
+        // this is the "real" version of this source 
+        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
+
+        // if we did not fit the galaxy model, modelFits will be NULL
+        if (source->modelFits == NULL) continue;
+
+        // if we did not fit the galaxy model, galaxyFits will also be NULL
+        if (source->galaxyFits == NULL) continue;
+
+        pmModel *model = source->modelFits->data[0];
+        if (!model) return false;
+
+        // X,Y coordinates are stored with the model parameters
+        psF32 *PAR = model->params->data.F32;
+
+        psMetadata *row = psMetadataAlloc ();
+
+        // we write out the x,y positions so people can link to the psf either way (position or ID)
+        psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index", source->seq);
+        psMetadataAddF32 (row, PS_LIST_TAIL, "X_FIT",            0, "model x coordinate",             PAR[PM_PAR_XPOS]);
+        psMetadataAddF32 (row, PS_LIST_TAIL, "Y_FIT",            0, "model y coordinate",             PAR[PM_PAR_YPOS]);
+        psMetadataAddF32 (row, PS_LIST_TAIL, "NPIX",             0, "number of pixels for fits",      source->galaxyFits->nPix);
+
+        psVector *Flux = source->galaxyFits->Flux;
+        psVector *dFlux = source->galaxyFits->dFlux;
+        psVector *chisq = source->galaxyFits->chisq;
+
+        psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX",     PS_META_REPLACE, "normalization for galaxy flux", Flux);
+        psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX_ERR", PS_META_REPLACE, "error on normalization", dFlux);
+        psMetadataAddVector (row, PS_LIST_TAIL, "GAL_CHISQ",    PS_META_REPLACE, "galaxy fit chisq", chisq);
+
+        psArrayAdd (table, 100, row);
+        psFree (row);
+    }
+
+    if (table->n == 0) {
+        if (!psFitsWriteBlank (fits, outhead, extname)) {
+            psError(psErrorCodeLast(), false, "Unable to write empty sources file.");
+            psFree(outhead);
+            psFree(table);
+            return false;
+        }
+        psFree (outhead);
+        psFree (table);
+        return true;
+    }
+
+    psTrace ("pmFPAfile", 5, "writing galaxy data %s\n", extname);
+    if (!psFitsWriteTable (fits, outhead, table, extname)) {
+        psError(psErrorCodeLast(), false, "writing galaxy data %s\n", extname);
+        psFree (outhead);
+        psFree(table);
+        return false;
+    }
+    psFree (outhead);
+    psFree (table);
+    return true;
+}
+