Changeset 6556 for branches/rel10_ifa/psModules/src/objects/pmSource.c

Timestamp:

Mar 8, 2006, 5:14:23 PM (20 years ago)

Author:

magnier

Message:

major rework of objects code

File:

• branches/rel10_ifa/psModules/src/objects/pmSource.c (modified) (12 diffs)

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

@@ -1 @@
-/** @file  pmObjects.c
+/** @file  pmSource.c
  *
- *  This file will ...
+ *  Functions to define and manipulate sources on images
  *
  *  @author GLG, MHPCC
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-03-07 06:33:35 $
+ *  @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-03-09 03:14:23 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -17 +17 @@
 #include <string.h>
 #include "pslib.h"
-#include "pmObjects.h"
-#include "pmModelGroup.h"
+#include "pmHDU.h"
+#include "pmFPA.h"
+#include "pmPeaks.h"
+#include "pmMoments.h"
+#include "pmModel.h"
+#include "pmSource.h"
 
 static void sourceFree(pmSource *tmp)
@@ -50 +54 @@
     tmp->modelPSF = NULL;
     tmp->modelEXT = NULL;
-    tmp->type = PM_SOURCE_UNKNOWN;
-    tmp->mode = PM_SOURCE_DEFAULT;
+    tmp->type = PM_SOURCE_TYPE_UNKNOWN;
+    tmp->mode = PM_SOURCE_MODE_DEFAULT;
     psMemSetDeallocator(tmp, (psFreeFunc) sourceFree);
 
@@ -219 +223 @@
 
         // XXX EAM : this lets us takes the single highest peak
-        psArraySort (peaks, pmComparePeakDescend);
+        psArraySort (peaks, pmPeaksCompareDescend);
         clump = peaks->data[0];
         psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts);
@@ -331 +335 @@
         // inner = psRegionForSquare (tmpSrc->peak->x - tmpSrc->mask->col0, tmpSrc->peak->y - tmpSrc->mask->row0, 2);
         inner = psRegionForSquare (tmpSrc->peak->x, tmpSrc->peak->y, 2);
-        int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PSPHOT_MASK_SATURATED);
+        int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PM_SOURCE_MASK_SATURATED);
 
         // saturated star (size consistent with PSF or larger)
@@ -338 +342 @@
         big = true;
         if ((Nsatpix > 1) && big) {
-            tmpSrc->type = PM_SOURCE_STAR;
-            tmpSrc->mode = PM_SOURCE_SATSTAR;
+            tmpSrc->type = PM_SOURCE_TYPE_STAR;
+            tmpSrc->mode = PM_SOURCE_MODE_SATSTAR;
             Nsatstar ++;
             continue;
@@ -346 +350 @@
         // saturated object (not a star, eg bleed trails, hot pixels)
         if (Nsatpix > 1) {
-            tmpSrc->type = PM_SOURCE_SATURATED;
-            tmpSrc->mode = PM_SOURCE_DEFAULT;
+            tmpSrc->type = PM_SOURCE_TYPE_SATURATED;
+            tmpSrc->mode = PM_SOURCE_MODE_DEFAULT;
             Nsat ++;
             continue;
@@ -356 +360 @@
         // only set candidate defects if
         if ((sigX < 0.05) || (sigY < 0.05)) {
-            tmpSrc->type = PM_SOURCE_DEFECT;
-            tmpSrc->mode = PM_SOURCE_DEFAULT;
+            tmpSrc->type = PM_SOURCE_TYPE_DEFECT;
+            tmpSrc->mode = PM_SOURCE_MODE_DEFAULT;
             Ncr ++;
             continue;
@@ -364 +368 @@
         // likely unsaturated extended source (too large to be stellar)
         if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
-            tmpSrc->type = PM_SOURCE_EXTENDED;
-            tmpSrc->mode = PM_SOURCE_DEFAULT;
+            tmpSrc->type = PM_SOURCE_TYPE_EXTENDED;
+            tmpSrc->mode = PM_SOURCE_MODE_DEFAULT;
             Next ++;
             continue;
@@ -378 +382 @@
         psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
         if ((tmpSrc->moments->SN > PSF_SN_LIM) && (radius < 1.5)) {
-            tmpSrc->type = PM_SOURCE_STAR;
-            tmpSrc->mode = PM_SOURCE_PSFSTAR;
+            tmpSrc->type = PM_SOURCE_TYPE_STAR;
+            tmpSrc->mode = PM_SOURCE_MODE_PSFSTAR;
             Npsf ++;
             continue;
@@ -385 +389 @@
 
         // random type of star
-        tmpSrc->type = PM_SOURCE_STAR;
-        tmpSrc->mode = PM_SOURCE_DEFAULT;
+        tmpSrc->type = PM_SOURCE_TYPE_STAR;
+        tmpSrc->mode = PM_SOURCE_MODE_DEFAULT;
     }
 
@@ -409 +413 @@
 }
 
+/******************************************************************************
+pmSourceMoments(source, radius): this function takes a subImage defined in the
+pmSource data structure, along with the peak location, and determines the
+various moments associated with that peak.
+ 
+Requires the following to have been created:
+    pmSource
+    pmSource->peak
+    pmSource->pixels
+    pmSource->weight
+    pmSource->mask
+ 
+XXX: The peak calculations are done in image coords, not subImage coords.
+ 
+XXX EAM : this version clips input pixels on S/N
+XXX EAM : this version returns false for several reasons
+*****************************************************************************/
+# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)
+
+bool pmSourceMoments(pmSource *source,
+                     psF32 radius)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
+    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
+    PS_ASSERT_PTR_NON_NULL(source->mask, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
+
+    //
+    // XXX: Verify the setting for sky if source->moments == NULL.
+    //
+    psF32 sky = 0.0;
+    if (source->moments == NULL) {
+        source->moments = pmMomentsAlloc();
+    } else {
+        sky = source->moments->Sky;
+    }
+
+    //
+    // Sum = SUM (z - sky)
+    // X1  = SUM (x - xc)*(z - sky)
+    // X2  = SUM (x - xc)^2 * (z - sky)
+    // XY  = SUM (x - xc)*(y - yc)*(z - sky)
+    //
+    psF32 peakPixel = -PS_MAX_F32;
+    psS32 numPixels = 0;
+    psF32 Sum = 0.0;
+    psF32 Var = 0.0;
+    psF32 X1 = 0.0;
+    psF32 Y1 = 0.0;
+    psF32 X2 = 0.0;
+    psF32 Y2 = 0.0;
+    psF32 XY = 0.0;
+    psF32 x  = 0;
+    psF32 y  = 0;
+    psF32 R2 = PS_SQR(radius);
+
+    psF32 xPeak = source->peak->x;
+    psF32 yPeak = source->peak->y;
+    psF32 xOff = source->pixels->col0 - source->peak->x;
+    psF32 yOff = source->pixels->row0 - source->peak->y;
+
+    // XXX why do I get different results for these two methods of finding Sx?
+    // XXX Sx, Sy would be better measured if we clip pixels close to sky
+    // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2?
+    // We loop through all pixels in this subimage (source->pixels), and for each
+    // pixel that is not masked, AND within the radius of the peak pixel, we
+    // proceed with the moments calculation.  need to do two loops for a
+    // numerically stable result.  first loop: get the sums.
+    // XXX EAM : mask == 0 is valid
+
+    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
+
+        psF32 *vPix = source->pixels->data.F32[row];
+        psF32 *vWgt = source->weight->data.F32[row];
+        psU8  *vMsk = (source->mask == NULL) ? NULL : source->mask->data.U8[row];
+
+        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
+            if ((vMsk != NULL) && *vMsk) {
+                vMsk++;
+                continue;
+            }
+
+            psF32 xDiff = col + xOff;
+            psF32 yDiff = row + yOff;
+
+            // radius is just a function of (xDiff, yDiff)
+            if (!VALID_RADIUS(xDiff, yDiff, R2)) {
+                if (vMsk != NULL)
+                    vMsk++;
+                continue;
+            }
+
+            psF32 pDiff = *vPix - sky;
+            psF32 wDiff = *vWgt;
+
+            // XXX EAM : check for valid S/N in pixel
+            // XXX EAM : should this limit be user-defined?
+            if (PS_SQR(pDiff) < wDiff) {
+                if (vMsk != NULL)
+                    vMsk++;
+                continue;
+            }
+
+            Var += wDiff;
+            Sum += pDiff;
+
+            psF32 xWght = xDiff * pDiff;
+            psF32 yWght = yDiff * pDiff;
+
+            X1  += xWght;
+            Y1  += yWght;
+
+            XY  += xDiff * yWght;
+            X2  += xDiff * xWght;
+            Y2  += yDiff * yWght;
+
+            peakPixel = PS_MAX (*vPix, peakPixel);
+            numPixels++;
+            if (vMsk != NULL)
+                vMsk++;
+        }
+    }
+
+    // if we have less than (1/4) of the possible pixels, force a retry
+    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
+    if ((numPixels < 0.75*R2) || (Sum <= 0)) {
+        psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n");
+        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        return (false);
+    }
+
+    psTrace (".psModules.pmSourceMoments", 5,
+             "sky: %f  Sum: %f  X1: %f  Y1: %f  X2: %f  Y2: %f  XY: %f  Npix: %d\n",
+             sky, Sum, X1, Y1, X2, Y2, XY, numPixels);
+
+    //
+    // first moment X  = X1/Sum + xc
+    // second moment X = sqrt (X2/Sum - (X1/Sum)^2)
+    // Sxy             = XY / Sum
+    //
+    x = X1/Sum;
+    y = Y1/Sum;
+    if ((fabs(x) > radius) || (fabs(y) > radius)) {
+        psTrace (".psModules.pmSourceMoments", 3,
+                 "large centroid swing; invalid peak %d, %d\n",
+                 source->peak->x, source->peak->y);
+        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        return (false);
+    }
+
+    source->moments->x = x + xPeak;
+    source->moments->y = y + yPeak;
+
+    // XXX EAM : Sxy needs to have x*y subtracted
+    source->moments->Sxy = XY/Sum - x*y;
+    source->moments->Sum = Sum;
+    source->moments->SN  = Sum / sqrt(Var);
+    source->moments->Peak = peakPixel;
+    source->moments->nPixels = numPixels;
+
+    // XXX EAM : these values can be negative, so we need to limit the range
+    source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0));
+    source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0));
+
+    psTrace (".psModules.pmSourceMoments", 4,
+             "sky: %f  Sum: %f  x: %f  y: %f  Sx: %f  Sy: %f  Sxy: %f\n",
+             sky, Sum, source->moments->x, source->moments->y,
+             source->moments->Sx, source->moments->Sy, source->moments->Sxy);
+
+    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
+    return(true);
+}
+
+pmModel *pmSourceSelectModel (pmSource *source)
+{
+    switch (source->type) {
+    case PM_SOURCE_TYPE_STAR:
+        return source->modelPSF;
+
+    case PM_SOURCE_TYPE_EXTENDED:
+        return source->modelEXT;
+
+    default:
+        return NULL;
+    }
+    return NULL;
+}