Changeset 6556 for branches/rel10_ifa/psModules/src/objects/pmMoments.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/pmMoments.c (modified) (3 diffs)

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

@@ -1 @@
-/** @file  pmObjects.c
+/** @file  pmMoments.c
  *
- *  This file will ...
+ *  Functions defining the pmMoments structure
  *
  *  @author GLG, MHPCC
- *  @author EAM, IfA: significant modifications.
+ *  @author EAM, IfA
  *
- *  @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
@@ -13 +13 @@
  */
 
-#include <stdio.h>
-#include <math.h>
-#include <string.h>
 #include "pslib.h"
-#include "pmObjects.h"
-#include "pmModelGroup.h"
+#include "pmMoments.h"
 
 /******************************************************************************
@@ -42 +38 @@
 }
 
-/******************************************************************************
-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);
-}
-