Changeset 31451 for trunk/psModules/src/objects/pmSourceMoments.c

Timestamp:

May 5, 2011, 11:02:53 AM (15 years ago)

Author:

eugene

Message:

merge updates from eam branch 20110404

File:

• trunk/psModules/src/objects/pmSourceMoments.c (modified) (19 diffs)

trunk/psModules/src/objects/pmSourceMoments.c

@@ -66 +66 @@
 // if mode & EXTERNAL or mode2 & MATCHED, do not re-calculate the centroid (use peak as centroid)
 
-bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal)
+bool pmSourceMoments(pmSource *source, float radius, float sigma, float minSN, float minKronRadius, psImageMaskType maskVal)
 {
     PS_ASSERT_PTR_NON_NULL(source, false);
@@ -74 +74 @@
 
     // this function assumes the sky has been well-subtracted for the image
-    psF32 sky = 0.0;
+    float sky = 0.0;
 
     if (source->moments == NULL) {
@@ -80 +80 @@
     }
 
-    psF32 Sum = 0.0;
-    psF32 Var = 0.0;
-    psF32 SumCore = 0.0;
-    psF32 VarCore = 0.0;
-    psF32 R2 = PS_SQR(radius);
-    psF32 minSN2 = PS_SQR(minSN);
-    psF32 rsigma2 = 0.5 / PS_SQR(sigma);
+    float Sum = 0.0;
+    float Var = 0.0;
+    float SumCore = 0.0;
+    float VarCore = 0.0;
+    float R2 = PS_SQR(radius);
+    float minSN2 = PS_SQR(minSN);
+    float rsigma2 = 0.5 / PS_SQR(sigma);
 
     // a note about coordinates: coordinates of objects throughout psphot refer to the primary
@@ -109 +109 @@
     // Xn  = SUM (x - xc)^n * (z - sky)
 
-    psF32 RF = 0.0;
-    psF32 RH = 0.0;
-    psF32 RS = 0.0;
-    psF32 XX = 0.0;
-    psF32 XY = 0.0;
-    psF32 YY = 0.0;
-    psF32 XXX = 0.0;
-    psF32 XXY = 0.0;
-    psF32 XYY = 0.0;
-    psF32 YYY = 0.0;
-    psF32 XXXX = 0.0;
-    psF32 XXXY = 0.0;
-    psF32 XXYY = 0.0;
-    psF32 XYYY = 0.0;
-    psF32 YYYY = 0.0;
+    float RFW = 0.0;
+    float RHW = 0.0;
+
+    float RF = 0.0;
+    float RH = 0.0;
+    float RS = 0.0;
+    float XX = 0.0;
+    float XY = 0.0;
+    float YY = 0.0;
+    float XXX = 0.0;
+    float XXY = 0.0;
+    float XYY = 0.0;
+    float YYY = 0.0;
+    float XXXX = 0.0;
+    float XXXY = 0.0;
+    float XXYY = 0.0;
+    float XYYY = 0.0;
+    float YYYY = 0.0;
 
     Sum = 0.0;  // the second pass may include slightly different pixels, re-determine Sum
+
+    // float dX = source->moments->Mx - source->peak->xf;
+    // float dY = source->moments->My - source->peak->yf;
+    // float dR = hypot(dX, dY);
+    // float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
+    // float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
+    float Xo = source->moments->Mx;
+    float Yo = source->moments->My;
 
     // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
     // xCM, yCM from pixel coords to pixel index here.
-    psF32 xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage
-    psF32 yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage
+    float xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
+    float yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
 
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
 
-        psF32 yDiff = row - yCM;
+        float yDiff = row - yCM;
         if (fabs(yDiff) > radius) continue;
 
-        psF32 *vPix = source->pixels->data.F32[row];
-        psF32 *vWgt = source->variance->data.F32[row];
+        float *vPix = source->pixels->data.F32[row];
+        float *vWgt = source->variance->data.F32[row];
+
         psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
+        // psImageMaskType  *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
 
         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
@@ -151 +164 @@
             if (isnan(*vPix)) continue;
 
-            psF32 xDiff = col - xCM;
+            float xDiff = col - xCM;
             if (fabs(xDiff) > radius) continue;
 
             // radius is just a function of (xDiff, yDiff)
-            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
             if (r2 > R2) continue;
 
-            psF32 fDiff = *vPix - sky;
-            psF32 pDiff = fDiff;
-            psF32 wDiff = *vWgt;
+            float fDiff = *vPix - sky;
+            float pDiff = fDiff;
+            float wDiff = *vWgt;
 
             // skip pixels below specified significance level.  this is allowed, but should be
@@ -171 +184 @@
             // weighting over weights the sky for faint sources
             if (sigma > 0.0) {
-                psF32 z = r2 * rsigma2;
+                float z = r2 * rsigma2;
                 assert (z >= 0.0);
-                psF32 weight  = exp(-z);
+                float weight  = exp(-z);
 
                 wDiff *= weight;
@@ -182 +195 @@
 
             // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
-            psF32 r = sqrt(r2);
-            psF32 rf = r * fDiff;
-            psF32 rh = sqrt(r) * fDiff;
-            psF32 rs = fDiff;
-
-            psF32 x = xDiff * pDiff;
-            psF32 y = yDiff * pDiff;
-
-            psF32 xx = xDiff * x;
-            psF32 xy = xDiff * y;
-            psF32 yy = yDiff * y;
-
-            psF32 xxx = xDiff * xx / r;
-            psF32 xxy = xDiff * xy / r;
-            psF32 xyy = xDiff * yy / r;
-            psF32 yyy = yDiff * yy / r;
-
-            psF32 xxxx = xDiff * xxx / r2;
-            psF32 xxxy = xDiff * xxy / r2;
-            psF32 xxyy = xDiff * xyy / r2;
-            psF32 xyyy = xDiff * yyy / r2;
-            psF32 yyyy = yDiff * yyy / r2;
+            float r = sqrt(r2);
+            float rf = r * fDiff;
+            float rh = sqrt(r) * fDiff;
+            float rs = fDiff;
+
+            float rfw = r * pDiff;
+            float rhw = sqrt(r) * pDiff;
+
+            float x = xDiff * pDiff;
+            float y = yDiff * pDiff;
+
+            float xx = xDiff * x;
+            float xy = xDiff * y;
+            float yy = yDiff * y;
+
+            float xxx = xDiff * xx / r;
+            float xxy = xDiff * xy / r;
+            float xyy = xDiff * yy / r;
+            float yyy = yDiff * yy / r;
+
+            float xxxx = xDiff * xxx / r2;
+            float xxxy = xDiff * xxy / r2;
+            float xxyy = xDiff * xyy / r2;
+            float xyyy = xDiff * yyy / r2;
+            float yyyy = yDiff * yyy / r2;
 
             RF  += rf;
             RH  += rh;
             RS  += rs;
+
+            RFW  += rfw;
+            RHW  += rhw;
 
             XX  += xx;
@@ -244 +263 @@
     source->moments->Myyyy = YYYY/Sum;
 
-    // Calculate the Kron magnitude (make this block optional?)
-    float radKinner = 1.0*source->moments->Mrf;
-    float radKron   = 2.5*source->moments->Mrf;
-    float radKouter = 4.0*source->moments->Mrf;
+    // if Mrf (first radial moment) is very small, we are getting into low-significance
+    // territory.  saturate at minKronRadius.  conversely, if Mrf is > radius, we are clearly
+    // making an error.  saturate at radius.
+    float kronRefRadius = MIN(radius, MAX(minKronRadius, source->moments->Mrf));
+
+    float radKinner = 1.0*kronRefRadius;
+    float radKron   = 2.5*kronRefRadius;
+    float radKouter = 4.0*kronRefRadius;
 
     int nKronPix = 0;
@@ -259 +282 @@
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
 
-        psF32 yDiff = row - yCM;
+        float yDiff = row - yCM;
         if (fabs(yDiff) > radKouter) continue;
 
-        psF32 *vPix = source->pixels->data.F32[row];
-        psF32 *vWgt = source->variance->data.F32[row];
+        float *vPix = source->pixels->data.F32[row];
+        float *vWgt = source->variance->data.F32[row];
+
         psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
+        // psImageMaskType  *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
 
         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
@@ -276 +301 @@
             if (isnan(*vPix)) continue;
 
-            psF32 xDiff = col - xCM;
+            float xDiff = col - xCM;
             if (fabs(xDiff) > radKouter) continue;
 
             // radKron is just a function of (xDiff, yDiff)
-            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
-
-            psF32 pDiff = *vPix - sky;
-            psF32 wDiff = *vWgt;
+            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+
+            float pDiff = *vPix - sky;
+            float wDiff = *vWgt;
 
             // skip pixels below specified significance level.  this is allowed, but should be
@@ -290 +315 @@
             if (PS_SQR(pDiff) < minSN2*wDiff) continue;
 
-            psF32 r  = sqrt(r2);
+            float r  = sqrt(r2);
             if (r < radKron) {
                 Sum += pDiff;
@@ -335 +360 @@
 }
 
-bool pmSourceMomentsGetCentroid(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal, float xGuess, float yGuess) { 
+bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) { 
 
     // First Pass: calculate the first moments (these are subtracted from the coordinates below)
@@ -342 +367 @@
     // .. etc
 
-    psF32 sky = 0.0;
-
-    psF32 peakPixel = -PS_MAX_F32;
+    float sky = 0.0;
+
+    float peakPixel = -PS_MAX_F32;
     psS32 numPixels = 0;
-    psF32 Sum = 0.0;
-    psF32 Var = 0.0;
-    psF32 X1 = 0.0;
-    psF32 Y1 = 0.0;
-    psF32 R2 = PS_SQR(radius);
-    psF32 minSN2 = PS_SQR(minSN);
-    psF32 rsigma2 = 0.5 / PS_SQR(sigma);
+    float Sum = 0.0;
+    float Var = 0.0;
+    float X1 = 0.0;
+    float Y1 = 0.0;
+    float R2 = PS_SQR(radius);
+    float minSN2 = PS_SQR(minSN);
+    float rsigma2 = 0.5 / PS_SQR(sigma);
 
     float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage
@@ -358 +383 @@
 
     // we are guaranteed to have a valid pixel and variance at this location (right? right?)
-    // psF32 weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
+    // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
     // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
     // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
@@ -370 +395 @@
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
 
-        psF32 yDiff = row + 0.5 - yPeak;
+        float yDiff = row + 0.5 - yPeak;
         if (fabs(yDiff) > radius) continue;
 
-        psF32 *vPix = source->pixels->data.F32[row];
-        psF32 *vWgt = source->variance->data.F32[row];
+        float *vPix = source->pixels->data.F32[row];
+        float *vWgt = source->variance->data.F32[row];
+
         psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
+        // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
 
         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
@@ -387 +414 @@
             if (isnan(*vPix)) continue;
 
-            psF32 xDiff = col + 0.5 - xPeak;
+            float xDiff = col + 0.5 - xPeak;
             if (fabs(xDiff) > radius) continue;
 
             // radius is just a function of (xDiff, yDiff)
-            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
             if (r2 > R2) continue;
 
-            psF32 pDiff = *vPix - sky;
-            psF32 wDiff = *vWgt;
+            float pDiff = *vPix - sky;
+            float wDiff = *vWgt;
 
             // skip pixels below specified significance level.  for a PSFs, this
@@ -408 +435 @@
             // weighting over weights the sky for faint sources
             if (sigma > 0.0) {
-                psF32 z  = r2*rsigma2;
+                float z  = r2*rsigma2;
                 assert (z >= 0.0);
-                psF32 weight  = exp(-z);
+                float weight  = exp(-z);
 
                 wDiff *= weight;
@@ -419 +446 @@
             Sum += pDiff;
 
-            psF32 xWght = xDiff * pDiff;
-            psF32 yWght = yDiff * pDiff;
+            float xWght = xDiff * pDiff;
+            float yWght = yDiff * pDiff;
 
             X1  += xWght;
@@ -477 +504 @@
     return true;
 }
+
+float pmSourceMinKronRadius(psArray *sources, float PSF_SN_LIM) {
+
+    psVector *radii = psVectorAllocEmpty(100, PS_TYPE_F32);
+
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *src = sources->data[i]; // Source of interest
+        if (!src || !src->moments) {
+            continue;
+        }
+
+        if (src->mode & PM_SOURCE_MODE_BLEND) {
+            continue;
+        }
+
+        if (!src->moments->nPixels) continue;
+
+        if (src->moments->SN < PSF_SN_LIM) continue;
+
+        // XXX put in Mxx,Myy cut based on clump location
+
+        psVectorAppend(radii, src->moments->Mrf);
+    }
+
+    // find the peak in this image
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+    if (!psVectorStats (stats, radii, NULL, NULL, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n");
+        psFree(stats);
+        return NAN;
+    }
+
+    float minRadius = stats->sampleMedian;
+
+    psFree(radii);
+    psFree(stats);
+    return minRadius;
+}
+