Changeset 36707

Timestamp:

Apr 30, 2014, 3:48:08 PM (12 years ago)

Author:

bills

Message:

read and write CFF files with more than one row per source

File:

• trunk/psModules/src/objects/pmSourceIO_CFF.c (modified) (6 diffs)

trunk/psModules/src/objects/pmSourceIO_CFF.c

@@ -53 +53 @@
 #include "pmSourceOutputs.h"
 
-// read in a readout from the fits file
+// read in sources readout from a cff fits file
 psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header)
 {
@@ -78 +78 @@
     // is large and ephemeral --- when the table gets blown away, whatever is allocated after
     // the table is read blocks the free.  In fact, it's better to read the table row by row.
-    long numSources = psFitsTableSize(fits); // Number of sources in table
-    psArray *sources = psArrayAlloc(numSources); // Array of sources, to return
-
-    // convert the table to the pmSource entriesa
-    for (int i = 0; i < numSources; i++) {
+    long numRows = psFitsTableSize(fits); // Number of rows in table
+    psArray *sources = psArrayAllocEmpty(numRows); // Array of sources, to return
+
+    // convert the table to the pmSource entries
+    unsigned int ID_last = 0;
+    pmSource *source = NULL;
+    for (int i = 0; i < numRows; i++) {
         psMetadata *row = psFitsReadTableRow(fits, i); // Table row
         if (!row) {
@@ -123 +125 @@
         float Sindex     = psMetadataLookupF32 (&status, row, "INDEX"); // Should this be PAR_07 not sersic index
 
-        pmSource *source = pmSourceAlloc ();
-        pmModel *model = pmModelAlloc (modelType);
-        source->modelPSF  = model;
-//        RoughClass wants source type to be unknown
-//        source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags
-        source->type = PM_SOURCE_TYPE_UNKNOWN;
-
-        // XXX we can set this in general, but for a specific image, we need to weed out SATSTARS and
-        // stars that are masked
-        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;
-
-        dPAR[PM_PAR_XPOS] = 0.0;
-        dPAR[PM_PAR_YPOS] = 0.0;
-
-        PAR[PM_PAR_SKY]   = 0.0;
-        dPAR[PM_PAR_SKY]  = 0.0;
-
-        PAR[PM_PAR_I0]    = 1.0;
-        dPAR[PM_PAR_I0]   = 0.0;
-
-        source->sky       = PAR[PM_PAR_SKY];
-        source->skyErr    = dPAR[PM_PAR_SKY];
-
-        source->psfMag    = 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;
-        axes.theta        = 0.0;
-        pmPSF_AxesToModel (PAR, axes, modelType);
-
-        // peak->detValue, rawFlux, smoothFlux are all set to the flux argument which is counts per second
-        source->peak      = pmPeakAlloc(X, Y, flux, PM_PEAK_LONE);
-        source->peak->xf  = X; // pmPeakAlloc converts X,Y to int, so reset here
-        source->peak->yf  = Y; // pmPeakAlloc converts X,Y to int, so reset here
-        source->peak->dx  = 0.0;
-        source->peak->dy  = 0.0;
-
-        source->moments = pmMomentsAlloc ();
-        source->moments->Mx = X;
-        source->moments->My = Y;
-        source->moments->Mrf = kronRadius * 0.4; // kronRadius is 2.5 * first radial moment
-
-        // Don't mark the moments as measured because that causes many fields to be left blank.
-        // The moments code knows not to change the position or the Mrf for external sources
-        // source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED;
-
-        if (isfinite(petRadius)) {
-            source->extpars = pmSourceExtendedParsAlloc ();
-            source->extpars->petrosianRadius = petRadius;
-        }
-
+        if (!source || ID != source->seq) {
+            if (source) {
+                psArrayAdd (sources, 1, source);
+                psFree(source);
+            }
+            source = pmSourceAlloc ();
+            pmModel *model = pmModelAlloc (modelType);
+            source->modelPSF  = model;
+            //        RoughClass wants source type to be unknown
+            //        source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags
+            source->type = PM_SOURCE_TYPE_UNKNOWN;
+
+            // XXX we can set this in general, but for a specific image, we need to weed out SATSTARS and
+            // stars that are masked
+            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;
+            ID_last = ID;
+
+            PAR[PM_PAR_XPOS]  = X;
+            PAR[PM_PAR_YPOS]  = Y;
+
+            dPAR[PM_PAR_XPOS] = 0.0;
+            dPAR[PM_PAR_YPOS] = 0.0;
+
+            PAR[PM_PAR_SKY]   = 0.0;
+            dPAR[PM_PAR_SKY]  = 0.0;
+
+            PAR[PM_PAR_I0]    = 1.0;
+            dPAR[PM_PAR_I0]   = 0.0;
+
+            source->sky       = PAR[PM_PAR_SKY];
+            source->skyErr    = dPAR[PM_PAR_SKY];
+
+            source->psfMag    = 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;
+            axes.theta        = 0.0;
+            pmPSF_AxesToModel (PAR, axes, modelType);
+
+            // peak->detValue, rawFlux, smoothFlux are all set to the flux argument which is counts per second
+            source->peak      = pmPeakAlloc(X, Y, flux, PM_PEAK_LONE);
+            source->peak->xf  = X; // pmPeakAlloc converts X,Y to int, so reset here
+            source->peak->yf  = Y; // pmPeakAlloc converts X,Y to int, so reset here
+            source->peak->dx  = 0.0;
+            source->peak->dy  = 0.0;
+
+            source->moments = pmMomentsAlloc ();
+            source->moments->Mx = X;
+            source->moments->My = Y;
+            source->moments->Mrf = kronRadius * 0.4; // kronRadius is 2.5 * first radial moment
+
+            // Don't mark the moments as measured because that causes many fields to be left blank.
+            // The moments code knows not to change the position or the Mrf for external sources
+            // source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED;
+
+            if (isfinite(petRadius)) {
+                source->extpars = pmSourceExtendedParsAlloc ();
+                source->extpars->petrosianRadius = petRadius;
+            }
+
+        }
         if (fitGalaxy && galaxyModelType >= 0) {
-            source->modelFits = psArrayAllocEmpty (1);
-            pmModel *model = pmModelAlloc(galaxyModelType);
-            psF32 *xPAR = model->params->data.F32;
-
-            xPAR[PM_PAR_SKY]  = 0.0;
-            xPAR[PM_PAR_I0]   = 1.0;
-            xPAR[PM_PAR_XPOS] = X;
+            if (!source->modelFits) {
+                source->modelFits = psArrayAllocEmpty (1);
+            }
+            pmModel *model = pmModelAlloc(galaxyModelType);
+            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;
             
@@ -212 +224 @@
             psArrayAdd (source->modelFits, 1, model);
 
-#ifdef notyet
-            // XXX: set source->modelEXT to this model and flag as extended so that we have a better
-            // shot at subtracting extended sources?
-
-            // This doesn't work right. Need to do some more work to flesh out the model before it can be
-            // used in subtaction. Better idea might be to do 2 passes in psphotFullForceReadout to
-            // First find the best extended models and then do everything else.
-            source->modelEXT = psMemIncrRefCounter(model);
-            // is this safe? (no see above)
-            source->type = PM_SOURCE_TYPE_EXTENDED;
-            source->mode |= PM_SOURCE_MODE_EXTMODEL; 
-#endif
-
             psFree (model);
         }
 
-        sources->data[i] = source;
         psFree(row);
+    }
+    if (source) {
+        // close out last source
+        psArrayAdd (sources, 1, source);
+        psFree(source);
     }
 
@@ -250 +253 @@
 
     pmModelType sersicModelType = pmModelClassGetType("PS_MODEL_SERSIC");
+    pmModelType devModelType    = pmModelClassGetType("PS_MODEL_DEV");
+    pmModelType selectedModelType = -1;
+    bool chooseBest = false;
+    bool chooseAll = false;
 
     psString modelToChoose = psMetadataLookupStr(&mdok, recipe, "EXT_MODEL_TYPE_FOR_CFF");
-    pmModelType selectedModelType = -1;
+
     if (mdok && modelToChoose != NULL) {
-        if (strcmp(modelToChoose, "BEST")) {
+        if (!strcmp(modelToChoose, "BEST")) {
+            chooseBest = true;
+        } else if (!strcmp(modelToChoose, "ALL")) {
+            chooseAll = true;
+        } else if (strcmp(modelToChoose, "PS_MODEL_SERSIC")) {
+            // We have selected a model type other than Sersic. 
+            // Save it's type for use below.  Sersic is handled specially
             selectedModelType = pmModelClassGetType(modelToChoose);
         }
     }
+
+    // minimum sersic index to force devModel
+    psF32 sersicMinDev = psMetadataLookupF32(&mdok, recipe, "EXT_MODEL_FORCE_DEV_SERSIC_MIN");
+    if (!mdok) {
+        sersicMinDev = NAN;
+    }
+
+    sources = psArraySort (sources, pmSourceSortBySeq);
 
     for (int i = 0; i < sources->n; i++) {
@@ -263 +284 @@
         pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
 
-        psF32 xPos, yPos, flux, rMajor, rMinor, theta;
-        psS32 modelType = 0;
+        #define MAX_ROWS_PER_SRC 10
+        psF32 xPos[MAX_ROWS_PER_SRC], yPos[MAX_ROWS_PER_SRC], flux[MAX_ROWS_PER_SRC];
+        psF32 rMajor[MAX_ROWS_PER_SRC], rMinor[MAX_ROWS_PER_SRC], theta[MAX_ROWS_PER_SRC];
+        psS32 modelType[MAX_ROWS_PER_SRC];
+        psF32 sersicIndex[MAX_ROWS_PER_SRC];
         bool fitGalaxy = false;
         bool psfStar = (source->mode & PM_SOURCE_MODE_PSFSTAR) ? true : false;
-        psF32 sersicIndex = 0;
-        // For now only perform galaxy fits on extended objects
-        if (source->modelEXT == 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;
-            rMajor = 0;
-            rMinor = 0;
-            theta = 0;
-        } else {
-            //   Find the model with the selected type. If selected type is -1 choose the one selected ad
-            //   modelEXT which was the best
-            int iModel = -1;
-            if (source->modelEXT) {
-                pmModelType ext_model_type =  selectedModelType != -1  ? selectedModelType : source->modelEXT->type;
+        int n_rows = 0;
+
+        // start with psf model
+        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[0] = model->params->data.F32[PM_PAR_XPOS];
+        yPos[0] = model->params->data.F32[PM_PAR_YPOS];
+        flux[0] = source->psfFlux;
+        rMajor[0] = 0;
+        rMinor[0] = 0;
+        theta[0] = 0;
+        modelType[0] = -1;
+        sersicIndex[0] = NAN;
+
+        if (source->modelFits != NULL) {
+            // figure out which models to use based on recipe paramters
+            if (chooseAll) {
+                // Save parameters for all valid extended models
+
+                // but make sure we aren't going to overflow our arrays
+                assert (source->modelFits->n < MAX_ROWS_PER_SRC);
+
                 for (int j=0; j<source->modelFits->n; j++) {
-                    pmModel *aModel = source->modelFits->data[j];
-                    if (aModel->type == ext_model_type) {
-                        iModel = j;
-                        break;
+                    pmModel *model = source->modelFits->data[j];
+                    psF32 *PAR = model->params->data.F32;
+
+                    if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SYY])  && isfinite(PAR[PM_PAR_SXY]) &&
+                        isfinite(model->mag)) {
+
+                        xPos[n_rows] = PAR[PM_PAR_XPOS];
+                        yPos[n_rows] = PAR[PM_PAR_YPOS];
+
+                        psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type);
+                        rMajor[n_rows] = axes.major;
+                        rMinor[n_rows] = axes.minor;
+                        theta[n_rows]  = axes.theta*PS_DEG_RAD;
+                        flux[n_rows] = pow(10.0, -0.4*model->mag);
+                        modelType[n_rows] = model->type;
+                        fitGalaxy = true;
+                        if (model->params->n == 8) {
+                            // this will save the sersic index in all models but that might
+                            // be useful
+                            sersicIndex[n_rows] = 0.5 / PAR[PM_PAR_7];
+                        } else {
+                            sersicIndex[n_rows] = NAN;
+                        }
+
+                        n_rows++;
+                    }
+                }
+            } else {
+                int jModelSersic = -1;
+                int jModelDev = -1;
+                int jModelSelected = -1;
+                psF32 minChisq = NAN;
+                for (int j=0; j<source->modelFits->n; j++) {
+                    pmModel *model = source->modelFits->data[j];
+                    if (chooseBest) {
+                        // choose the model with lowest chisq
+                        if (isfinite(model->chisq) && (!isfinite(minChisq) || model->chisq < minChisq)) {
+                            jModelSelected = j;
+                            minChisq = model->chisq;
+                        }
+                    } else {
+                        // find the index of models of interest
+                        if (model->type == selectedModelType) {
+                            jModelSelected = j;
+                        } else if (model->type == sersicModelType) {
+                            jModelSersic = j;
+                        } else if (model->type == devModelType) {
+                            jModelDev = j;
+                        }
+                    }
+                }
+                if (jModelSelected >= 0 || jModelSersic >= 0) {
+                    // If a specific non-sersic model we take paramers from that one.
+                    // Otherwise we do the sersic model.
+                    pmModel *model = jModelSelected >= 0 ? source->modelFits->data[jModelSelected] :
+                                                           source->modelFits->data[jModelSersic];
+                    psF32 *PAR = model->params->data.F32;
+
+                    if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SYY])  && isfinite(PAR[PM_PAR_SXY]) &&
+                        isfinite(model->mag)) {
+
+                        xPos[0] = PAR[PM_PAR_XPOS];
+                        yPos[0] = PAR[PM_PAR_YPOS];
+
+                        psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type);
+                        rMajor[0] = axes.major;
+                        rMinor[0] = axes.minor;
+                        theta[0]  = axes.theta*PS_DEG_RAD;
+                        flux[0] = pow(10.0, -0.4*model->mag);
+                        modelType[0] = model->type;
+                        fitGalaxy = true;
+                        if (model->type == sersicModelType) {
+                            PS_ASSERT_FLOAT_LARGER_THAN(PAR[PM_PAR_7], 0.0, false);
+                            sersicIndex[0] = 0.5 / PAR[PM_PAR_7];
+                        }
+
+                        n_rows = 1;
+
+                        // Unless a specific non-sersic model type was selected do dev model for sources with
+                        // sersic index above the recipe limit.
+                        if (jModelSelected == -1 && jModelDev >= 0 && isfinite(sersicMinDev) &&
+                            sersicIndex[0] > sersicMinDev) {
+
+                            model = source->modelFits->data[jModelDev];
+                            if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SYY])  && isfinite(PAR[PM_PAR_SXY]) &&
+                                isfinite(model->mag)) { 
+
+                                xPos[1] = PAR[PM_PAR_XPOS];
+                                yPos[1] = PAR[PM_PAR_YPOS];
+
+                                psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type);
+                                rMajor[1] = axes.major;
+                                rMinor[1] = axes.minor;
+                                theta[1]  = axes.theta*PS_DEG_RAD;
+                                flux[1] = pow(10.0, -0.4*model->mag);
+                                modelType[1] = model->type;
+                                sersicIndex[1] = NAN;
+                                n_rows = 2;
+                            }
+                        }
                     }
                 }
             }
-            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) {
-                PS_ASSERT_FLOAT_LARGER_THAN(PAR[PM_PAR_7], 0.0, false);
-                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 per second",       flux/exptime);
-        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);
-        psMetadataAddF32 (row, PS_LIST_TAIL, "INDEX",            0, "sersic index",          sersicIndex);
-
-        psArrayAdd(table, 100, row);
-        psFree(row);
+        }
+        if (n_rows == 0) {
+            // didn't find a suitable extended model we have one row with the psf parameters
+            n_rows = 1;
+        }
+
+        for (int j = 0; j < n_rows; j++) {
+            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[j]);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "Y",                0, "y coordinate",          yPos[j]);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "FLUX",             0, "flux per second",       flux[j]/exptime);
+            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[j]);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "R_MINOR",          0, "radius of minor axis",  rMinor[j]);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "THETA",            0, "theta",                 theta[j]);
+            psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TYPE",       0, "model type",            modelType[j]);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "INDEX",            0, "sersic index",          sersicIndex[j]);
+
+            psArrayAdd(table, 100, row);
+            psFree(row);
+        }
     }