Changeset 27839 for branches/simmosaic_branches/psModules/src/objects/pmSource.c

Timestamp:

May 3, 2010, 8:45:22 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simmosaic_branches

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/objects/pmSource.c (modified) (26 diffs)

branches/simmosaic_branches/psModules/src/objects/pmSource.c

@@ -3 +3 @@
  *  Functions to define and manipulate sources on images
  *
- *  @author GLG, MHPCC
- *  @author EAM, IfA: significant modifications.
+ *  @author EAM, IfA
+ *  @author GLG, MHPCC (initial code base)
  *
  *  @version $Revision: 1.70 $ $Name: not supported by cvs2svn $
  *  @date $Date: 2009-02-16 22:29:59 $
- *
- *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
- *
+ *  Copyright 2009 Institute for Astronomy, University of Hawaii
  */
 
@@ -48 +46 @@
     psFree(tmp->maskView);
     psFree(tmp->modelFlux);
-    psFree(tmp->psfFlux);
+    psFree(tmp->psfImage);
     psFree(tmp->moments);
     psFree(tmp->modelPSF);
@@ -54 +52 @@
     psFree(tmp->modelFits);
     psFree(tmp->extpars);
+    psFree(tmp->moments);
+    psFree(tmp->diffStats);
     psFree(tmp->blends);
     psTrace("psModules.objects", 10, "---- end ----\n");
@@ -70 +70 @@
     psFree (source->maskView);
     psFree (source->modelFlux);
-    psFree (source->psfFlux);
+    psFree (source->psfImage);
 
     source->pixels = NULL;
@@ -77 +77 @@
     source->maskView = NULL;
     source->modelFlux = NULL;
-    source->psfFlux = NULL;
+    source->psfImage = NULL;
     return;
 }
@@ -105 +105 @@
     source->maskView = NULL;
     source->modelFlux = NULL;
-    source->psfFlux = NULL;
+    source->psfImage = NULL;
     source->moments = NULL;
     source->blends = NULL;
@@ -115 +115 @@
     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;
@@ -176 +180 @@
     source->type = in->type;
     source->mode = in->mode;
+    source->imageID = in->imageID;
 
     return(source);
@@ -242 +247 @@
     extend |= (mySource->maskView == NULL);
 
-    // extend = true;
     if (extend) {
         // re-create the subimage
@@ -250 +254 @@
 
         mySource->pixels   = psImageSubset(readout->image,  newRegion);
-        mySource->variance   = psImageSubset(readout->variance, newRegion);
+        mySource->variance = psImageSubset(readout->variance, newRegion);
         mySource->maskView = psImageSubset(readout->mask,   newRegion);
         mySource->region   = newRegion;
@@ -262 +266 @@
         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;
@@ -278 +282 @@
 // psphot-specific function which applies the recipe values
 // only apply selection to sources within specified region
-pmPSFClump pmSourcePSFClump(psRegion *region, psArray *sources, psMetadata *recipe)
+pmPSFClump pmSourcePSFClump(psImage **savedImage, psRegion *region, psArray *sources, float PSF_SN_LIM, float PSF_CLUMP_GRID_SCALE, psF32 SX_MAX, psF32 SY_MAX, psF32 AR_MAX)
 {
     psTrace("psModules.objects", 10, "---- begin ----\n");
 
     psArray *peaks  = NULL;
-    pmPSFClump errorClump = {-1.0, -1.0, 0.0, 0.0};
-    pmPSFClump emptyClump = {+0.0, +0.0, 0.0, 0.0};
+    pmPSFClump errorClump = {-1.0, -1.0, 0.0, 0.0, 0, 0.0};
+    pmPSFClump emptyClump = {+0.0, +0.0, 0.0, 0.0, 0, 0.0};
     pmPSFClump psfClump;
 
     PS_ASSERT_PTR_NON_NULL(sources, errorClump);
-    PS_ASSERT_PTR_NON_NULL(recipe, errorClump);
-
-    bool status = true;                 // Status of MD lookup
-    float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM");
-    if (!status) {
-        PSF_SN_LIM = 0;
-    }
-    float PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_GRID_SCALE");
-    if (!status) {
-        PSF_CLUMP_GRID_SCALE = 0.1;
-    }
 
     // find the sigmaX, sigmaY clump
     {
-        psF32 SX_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SX_MAX");
-        if (!status) {
-            psWarning("MOMENTS_SX_MAX not set in recipe");
-            SX_MAX = 10.0;
-        }
-        psF32 SY_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SY_MAX");
-        if (!status) {
-            psWarning("MOMENTS_SY_MAX not set in recipe");
-            SY_MAX = 10.0;
-        }
-        psF32 AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX");
-        if (!status) {
-            psWarning("MOMENTS_AR_MAX not set in recipe");
-            AR_MAX =  3.0;
-        }
         psF32 AR_MIN = 1.0 / AR_MAX;
 
         // construct a sigma-plane image
-        int numCols = SX_MAX / PSF_CLUMP_GRID_SCALE, numRows = SY_MAX / PSF_CLUMP_GRID_SCALE; // Size of sigma-plane image
+        int numCols = 1 + SX_MAX / PSF_CLUMP_GRID_SCALE; // Size of sigma-plane image
+        int numRows = 1 + SY_MAX / PSF_CLUMP_GRID_SCALE; // Size of sigma-plane image
         psTrace("psModules.objects", 10, "sigma-plane dimensions: %dx%d\n", numCols, numRows);
         psImage *splane = psImageAlloc(numCols, numRows, PS_TYPE_F32); // sigma-plane image
@@ -333 +312 @@
             }
 
-            int x = src->peak->x, y = src->peak->y; // Coordinates of peak
-            if (x < region->x0 || x > region->x1 || y < region->y0 || y > region->y1) {
-                continue;
-            }
+            if (region) {
+                int x = src->peak->x, y = src->peak->y; // Coordinates of peak
+                if (x < region->x0 || x > region->x1 || y < region->y0 || y > region->y1) {
+                    continue;
+                }
+            }
 
             if (src->mode & PM_SOURCE_MODE_BLEND) {
                 continue;
             }
+
+            if (!src->moments->nPixels) continue;
 
             if (src->moments->SN < PSF_SN_LIM) {
@@ -385 +368 @@
 
         // find the peak in this image
-        psStats *stats = psStatsAlloc (PS_STAT_MAX);
+        psStats *stats = psStatsAlloc (PS_STAT_MAX | PS_STAT_SAMPLE_STDEV);
         if (!psImageStats (stats, splane, NULL, 0)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n");
@@ -395 +378 @@
         psTrace ("psModules.objects", 2, "clump threshold is %f\n", stats[0].max/2);
 
-        const bool keep_psf_clump = psMetadataLookupBool(NULL, recipe, "KEEP_PSF_CLUMP");
-        if (keep_psf_clump)
-        {
-            psMetadataAdd(recipe, PS_LIST_TAIL,
-                          "PSF_CLUMP", PS_DATA_IMAGE, "Image of PSF coefficients", splane);
+        psfClump.nSigma = stats->sampleStdev;
+
+        if (savedImage) {
+            *savedImage = psMemIncrRefCounter(splane);
         }
         psFree (splane);
@@ -405 +387 @@
 
         // if we failed to find a valid peak, return the empty clump (failure signal)
-        if (peaks == NULL)
+        if (peaks == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "failure in peak analysis for PSF clump.\n");
+            psFree (peaks);
+            return emptyClump;
+        }
+
+        if (peaks->n == 0)
         {
             psLogMsg ("psphot", 3, "failed to find a peak in the PSF clump image\n");
@@ -413 +401 @@
                 psLogMsg ("psphot", 3, "no significant peak\n");
             }
+            psFree (peaks);
             return (emptyClump);
         }
@@ -431 +420 @@
         psTrace ("psModules.objects", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->value);
 
+        // XXX store the mean sigma?
+        float meanSigma = psfClump.nSigma;
+        psfClump.nStars = clump->value;
+        psfClump.nSigma = clump->value / meanSigma;
+
         // define section window for clump
         minSx = clump->x * PSF_CLUMP_GRID_SCALE - 2.0*PSF_CLUMP_GRID_SCALE;
@@ -453 +447 @@
                 continue;
 
-            if (tmpSrc->peak->x < region->x0) continue;
-            if (tmpSrc->peak->x > region->x1) continue;
-            if (tmpSrc->peak->y < region->y0) continue;
-            if (tmpSrc->peak->y > region->y1) continue;
+            if (region) {
+                if (tmpSrc->peak->x < region->x0) continue;
+                if (tmpSrc->peak->x > region->x1) continue;
+                if (tmpSrc->peak->y < region->y0) continue;
+                if (tmpSrc->peak->y > region->y1) continue;
+            }
 
             if (tmpSrc->moments->Mxx < minSx)
@@ -512 +508 @@
 *****************************************************************************/
 
-bool pmSourceRoughClass(psRegion *region, psArray *sources, psMetadata *recipe, pmPSFClump clump, psImageMaskType maskSat)
+bool pmSourceRoughClass(psRegion *region, psArray *sources, float PSF_SN_LIM, float PSF_CLUMP_NSIGMA, pmPSFClump clump, psImageMaskType maskSat)
 {
     psTrace("psModules.objects", 10, "---- begin ----");
 
     PS_ASSERT_PTR_NON_NULL(sources, false);
-    PS_ASSERT_PTR_NON_NULL(recipe, false);
 
     int Nsat     = 0;
@@ -530 +525 @@
     psVector *starsn_peaks = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     psVector *starsn_moments = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
-
-    // get basic parameters, or set defaults
-    bool status;
-    float PSF_SN_LIM = psMetadataLookupF32 (&status, recipe, "PSF_SN_LIM");
-    if (!status) PSF_SN_LIM = 20.0;
-    float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
-    if (!status) PSF_CLUMP_NSIGMA = 1.5;
-
-    // float INNER_RADIUS = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
 
     pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN;
@@ -894 +880 @@
 
     // 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;
 }
@@ -919 +905 @@
     pmModel *model = pmSourceGetModel (NULL, source);
     if (model == NULL) return false;  // model must be defined
+
+    bool addNoise = mode & PM_MODEL_OP_NOISE;
 
     if (source->modelFlux) {
@@ -941 +929 @@
 
         psF32 **target = source->pixels->data.F32;
-        if (mode & PM_MODEL_OP_NOISE) {
-            // XXX need to scale by the gain...
+        if (addNoise) {
+            // when adding noise, we assume the shape and Io have been modified
             target = source->variance->data.F32;
         }
 
-        // XXX need to respect the source and model masks?
         for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
             int oy = iy + dY;
@@ -960 +947 @@
             }
         }
+        if (!addNoise) {
+            if (add) {
+                source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+            } else {
+                source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
+            }
+        }
         return true;
     }
 
     psImage *target = source->pixels;
-    if (mode & PM_MODEL_OP_NOISE) {
+    if (addNoise) {
         target = source->variance;
     }
 
-    if (add) {
-        status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
-    } else {
-        status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
+    if (!addNoise) {
+        if (add) {
+            status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
+            source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+        } else {
+            status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
+            source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
+        }
     }
 
@@ -1061 +1059 @@
     psF32 fA = (A->peak == NULL) ? 0 : A->peak->y;
     psF32 fB = (B->peak == NULL) ? 0 : B->peak->y;
+
+    psF32 diff = fA - fB;
+    if (diff > FLT_EPSILON) return (+1);
+    if (diff < FLT_EPSILON) return (-1);
+    return (0);
+}
+
+// sort by X (ascending)
+int pmSourceSortByX (const void **a, const void **b)
+{
+    pmSource *A = *(pmSource **)a;
+    pmSource *B = *(pmSource **)b;
+
+    psF32 fA = (A->peak == NULL) ? 0 : A->peak->x;
+    psF32 fB = (B->peak == NULL) ? 0 : B->peak->x;
 
     psF32 diff = fA - fB;