Changeset 27494

Timestamp:

Mar 28, 2010, 11:02:04 AM (16 years ago)

Author:

eugene

Message:

measuring diff stats; adding diff stats and flux in output

Location:

branches/eam_branches/20100225/psModules/src

Files:

• objects/Makefile.am (modified) (3 diffs)
• objects/pmGrowthCurveGenerate.c (modified) (1 diff)
• objects/pmSource.c (modified) (8 diffs)
• objects/pmSource.h (modified) (4 diffs)
• objects/pmSourceIO.c (modified) (4 diffs)
• objects/pmSourceIO.h (modified) (2 diffs)
• objects/pmSourcePhotometry.c (modified) (6 diffs)
• objects/pmSourcePhotometry.h (modified) (2 diffs)
• psmodules.h (modified) (1 diff)

branches/eam_branches/20100225/psModules/src/objects/Makefile.am

@@ -22 +22 @@
         pmSourceMasks.c \
         pmSourceMoments.c \
+        pmSourceDiffStats.c \
         pmSourceExtendedPars.c \
         pmSourceUtils.c \
@@ -40 +41 @@
         pmSourceIO_CMF_PS1_V1.c \
         pmSourceIO_CMF_PS1_V2.c \
+        pmSourceIO_CMF_PS1_DV1.c \
         pmSourceIO_MatchedRefs.c \
         pmSourcePlots.c \
@@ -79 +81 @@
         pmSource.h \
         pmSourceMasks.h \
+        pmSourceDiffStats.h \
         pmSourceExtendedPars.h \
         pmSourceUtils.h \

branches/eam_branches/20100225/psModules/src/objects/pmGrowthCurveGenerate.c

@@ -157 +157 @@
 
     // measure the fitMag for this model
-    pmSourcePhotometryModel (&fitMag, model);
+    pmSourcePhotometryModel (&fitMag, NULL, model);
     growth->fitMag = fitMag;
 

branches/eam_branches/20100225/psModules/src/objects/pmSource.c

@@ -46 +46 @@
     psFree(tmp->maskView);
     psFree(tmp->modelFlux);
-    psFree(tmp->psfFlux);
+    psFree(tmp->psfImage);
     psFree(tmp->moments);
     psFree(tmp->modelPSF);
@@ -52 +52 @@
     psFree(tmp->modelFits);
     psFree(tmp->extpars);
-    psFree(tmp->blends);
+    psFree(tmp->moments);
+    psFree(tmp->diffStats);
     psTrace("psModules.objects", 10, "---- end ----\n");
 }
@@ -68 +69 @@
     psFree (source->maskView);
     psFree (source->modelFlux);
-    psFree (source->psfFlux);
+    psFree (source->psfImage);
 
     source->pixels = NULL;
@@ -75 +76 @@
     source->maskView = NULL;
     source->modelFlux = NULL;
-    source->psfFlux = NULL;
+    source->psfImage = NULL;
     return;
 }
@@ -103 +104 @@
     source->maskView = NULL;
     source->modelFlux = NULL;
-    source->psfFlux = NULL;
+    source->psfImage = NULL;
     source->moments = NULL;
     source->blends = NULL;
@@ -113 +114 @@
     source->tmpFlags = 0;
     source->extpars = NULL;
+    source->diffStats = NULL;
+
     source->region = psRegionSet(NAN, NAN, NAN, NAN);
     psMemSetDeallocator(source, (psFreeFunc) sourceFree);
 
     // default values are NAN
-    source->psfMag = NAN;
+    source->psfMag     = NAN;
+    source->psfFlux    = NAN;
+    source->psfFluxErr = NAN;
     source->extMag = NAN;
     source->errMag = NAN;
@@ -259 +264 @@
         mySource->modelFlux = NULL;
 
-        // drop the old psfFlux pixels and force the user to re-create
-        psFree (mySource->psfFlux);
-        mySource->psfFlux = NULL;
+        // drop the old psfImage pixels and force the user to re-create
+        psFree (mySource->psfImage);
+        mySource->psfImage = NULL;
     }
     return extend;
@@ -873 +878 @@
 
     // if we already have a cached image, re-use that memory
-    source->psfFlux = psImageCopy (source->psfFlux, source->pixels, PS_TYPE_F32);
-    psImageInit (source->psfFlux, 0.0);
+    source->psfImage = psImageCopy (source->psfImage, source->pixels, PS_TYPE_F32);
+    psImageInit (source->psfImage, 0.0);
 
     // in some places (psphotEnsemble), we need a normalized version
     // in others, we just want the model.  which is more commonly used?
-    // psfFlux always has unity normalization (I0 = 1.0)
-    pmModelAdd (source->psfFlux, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
+    // psfImage always has unity normalization (I0 = 1.0)
+    pmModelAdd (source->psfImage, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
     return true;
 }

branches/eam_branches/20100225/psModules/src/objects/pmSource.h

@@ -16 +16 @@
 #include "pmMoments.h"
 #include "pmSourceExtendedPars.h"
+#include "pmSourceDiffStats.h"
 
 /// @addtogroup Objects Object Detection / Analysis Functions
@@ -65 +66 @@
     psImage *maskView;                  ///< view into global image mask for this object region
     psImage *modelFlux;                 ///< cached copy of the best model for this source
-    psImage *psfFlux;                   ///< cached copy of the psf model for this source
+    psImage *psfImage;                   ///< cached copy of the psf model for this source
     pmMoments *moments;                 ///< Basic moments measured for the object.
     pmModel *modelPSF;                  ///< PSF Model fit (parameters and type)
@@ -75 +76 @@
     psArray *blends;                    ///< collection of sources thought to be confused with object
     float psfMag;                       ///< calculated from flux in modelPSF
+    float psfFlux;                      ///< calculated from flux in modelPSF
+    float psfFluxErr;                   ///< calculated from flux in modelPSF
     float extMag;                       ///< calculated from flux in modelEXT
     float errMag;                       ///< error in psfMag OR extMag (depending on type)
@@ -86 +89 @@
     psRegion region;                    ///< area on image covered by selected pixels
     pmSourceExtendedPars *extpars;      ///< extended source parameters
+    pmSourceDiffStats *diffStats;       ///< extra parameters for difference detections
 };
 

branches/eam_branches/20100225/psModules/src/objects/pmSourceIO.c

@@ -539 +539 @@
                 status &= pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname);
             }
+            if (!strcmp (exttype, "PS1_DV1")) {
+                status &= pmSourcesWrite_CMF_PS1_DV1 (file->fits, readout, sources, file->header, outhead, dataname);
+            }
 
             if (xsrcname) {
@@ -553 +556 @@
                     status &= pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, sources, xsrcname, recipe);
                 }
+                if (!strcmp (exttype, "PS1_DV1")) {
+                    status &= pmSourcesWrite_CMF_PS1_DV1_XSRC (file->fits, sources, xsrcname, recipe);
+                }
             }
             if (xfitname) {
@@ -567 +573 @@
                     status &= pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, sources, xfitname);
                 }
+                if (!strcmp (exttype, "PS1_DV1")) {
+                    status &= pmSourcesWrite_CMF_PS1_DV1_XFIT (file->fits, sources, xfitname);
+                }
             }
             psFree (outhead);
@@ -1015 +1024 @@
                 sources = pmSourcesRead_CMF_PS1_V2 (file->fits, hdu->header);
             }
+            if (!strcmp (exttype, "PS1_DV1")) {
+                sources = pmSourcesRead_CMF_PS1_DV1 (file->fits, hdu->header);
+            }
 
             if (!pmReadoutReadDetEff(file->fits, readout, deteffname)) {

branches/eam_branches/20100225/psModules/src/objects/pmSourceIO.h

@@ -43 +43 @@
 bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, psArray *sources, char *extname);
 
+bool pmSourcesWrite_CMF_PS1_DV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
+bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe);
+bool pmSourcesWrite_CMF_PS1_DV1_XFIT (psFits *fits, psArray *sources, char *extname);
+
 bool pmSource_CMF_WritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config);
 
@@ -53 +57 @@
 psArray *pmSourcesRead_CMF_PS1_V1 (psFits *fits, psMetadata *header);
 psArray *pmSourcesRead_CMF_PS1_V2 (psFits *fits, psMetadata *header);
+psArray *pmSourcesRead_CMF_PS1_DV1 (psFits *fits, psMetadata *header);
 
 bool pmSourcesWritePSFs (psArray *sources, char *filename);

branches/eam_branches/20100225/psModules/src/objects/pmSourcePhotometry.c

@@ -109 +109 @@
         psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y);
         psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y);
-        source->psfMag = -2.5*log10(fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0]);
+        source->psfFlux = fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0];
+        source->psfFluxErr = fluxScale * source->modelPSF->dparams->data.F32[PM_PAR_I0];
+        source->psfMag = -2.5*log10(source->psfFlux);
     } else {
-        status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF);
+        status = pmSourcePhotometryModel (&source->psfMag, &source->psfFlux, source->modelPSF);
+        source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]);
     }
 
@@ -119 +122 @@
         for (int i = 0; i < source->modelFits->n; i++) {
             pmModel *model = source->modelFits->data[i];
-            status = pmSourcePhotometryModel (&model->mag, model);
+            status = pmSourcePhotometryModel (&model->mag, NULL, model);
             if (model == source->modelEXT) foundEXT = true;
         }
@@ -125 +128 @@
             source->extMag = source->modelEXT->mag;
         } else {
-            status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
+            status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
         }
     } else {
         if (source->modelEXT) {
-            status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
+            status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
         }
     }
@@ -143 +146 @@
     if (mode & PM_SOURCE_PHOT_WEIGHT) {
         pmSourcePixelWeight (&source->pixWeight, model, source->maskObj, maskVal);
+    }
+
+    // measure the contribution of included pixels
+    if (mode & PM_SOURCE_PHOT_DIFFSTATS) {
+        pmSourceMeasureDiffStats (source, maskVal);
     }
 
@@ -217 +225 @@
 
 // return source model magnitude
-bool pmSourcePhotometryModel (float *fitMag, pmModel *model)
-{
-    PS_ASSERT_PTR_NON_NULL(fitMag, false);
-    if (model == NULL) {
-        return false;
-    }
-
-    float fitSum = 0;
-    *fitMag = NAN;
+bool pmSourcePhotometryModel (float *fitMag, float *fitFlux, pmModel *model)
+{
+    psAssert (fitMag || fitFlux, "at least one of magnitude or flux must be requested (not NULL)");
+    if (model == NULL) return false;
+
+    float mag  = NAN;
+    float flux = NAN;
 
     // measure fitMag
-    fitSum = model->modelFlux (model->params);
-    if (fitSum <= 0)
-        return false;
-    if (!isfinite(fitSum))
-        return false;
-    *fitMag = -2.5*log10(fitSum);
+    flux = model->modelFlux (model->params);
+    if (flux > 0) {
+        mag = -2.5*log10(flux);
+    }
+    if (fitMag) {
+        *fitMag = mag;
+    }
+    if (fitFlux) {
+        *fitFlux = flux;
+    }
+
+    if (flux <= 0) return false;
+    if (!isfinite(flux)) return false;
 
     return (true);
@@ -356 +369 @@
 
     *pixWeight = validSum / modelSum;
+    return (true);
+}
+
+# define FLUX_LIMIT 3.0
+
+// return source aperture magnitude
+bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal)
+{
+    PS_ASSERT_PTR_NON_NULL(source, false);
+
+    if (source->diffStats == NULL) {
+        source->diffStats = pmSourceDiffStatsAlloc();
+    }
+
+    float fGood = 0.0;
+    float fBad  = 0.0;
+    int   nGood = 0;
+    int   nMask = 0;
+    int   nBad  = 0;
+    
+    psImage *flux     = source->pixels;
+    psImage *variance = source->variance;
+    psImage *mask     = source->maskObj;
+
+    for (int iy = 0; iy < flux->numRows; iy++) {
+        for (int ix = 0; ix < flux->numCols; ix++) {
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) {
+                nMask ++;
+                continue;
+            }
+
+            float SN = flux->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
+
+            if (SN > +FLUX_LIMIT) { 
+                nGood ++;
+                fGood += flux->data.F32[iy][ix];
+            }
+
+            if (SN < -FLUX_LIMIT) { 
+                nBad ++;
+                fBad += flux->data.F32[iy][ix];
+            }
+        }
+    }
+
+    source->diffStats->fRatio     = fGood / (fGood + fBad);        
+    source->diffStats->nRatioBad  = nGood / (nGood + nBad);        
+    source->diffStats->nRatioMask = nGood / (nGood + nMask);       
+    source->diffStats->nRatioAll  = nGood / (nGood + nMask + nBad);
+    source->diffStats->nGood      = nGood;
+
     return (true);
 }

branches/eam_branches/20100225/psModules/src/objects/pmSourcePhotometry.h

@@ -29 +29 @@
 
 typedef enum {
-    PM_SOURCE_PHOT_NONE   = 0x0000,
-    PM_SOURCE_PHOT_GROWTH = 0x0001,
-    PM_SOURCE_PHOT_APCORR = 0x0002,
-    PM_SOURCE_PHOT_WEIGHT = 0x0004,
-    PM_SOURCE_PHOT_INTERP = 0x0008,
+    PM_SOURCE_PHOT_NONE      = 0x0000,
+    PM_SOURCE_PHOT_GROWTH    = 0x0001,
+    PM_SOURCE_PHOT_APCORR    = 0x0002,
+    PM_SOURCE_PHOT_WEIGHT    = 0x0004,
+    PM_SOURCE_PHOT_INTERP    = 0x0008,
+    PM_SOURCE_PHOT_DIFFSTATS = 0x0010,
 } pmSourcePhotometryMode;
 
 bool pmSourcePhotometryModel(
     float *fitMag,                      ///< integrated fit magnitude
+    float *fitFlux,                     ///< integrated fit magnitude
     pmModel *model                      ///< model used for photometry
 );
@@ -54 +56 @@
 bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, psImageMaskType maskVal, const float covarFactor);
 
+bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal);
 
 double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor);

branches/eam_branches/20100225/psModules/src/psmodules.h

@@ -115 +115 @@
 #include <pmDetections.h>
 #include <pmMoments.h>
+#include <pmSourceExtendedPars.h>
+#include <pmSourceDiffPars.h>
 #include <pmResiduals.h>
 #include <pmGrowthCurve.h>