Changeset 34772 for branches/czw_branch/20120906/psphot/src/psphotBlendFit.c

Timestamp:

Dec 6, 2012, 6:32:19 PM (14 years ago)

Author:

watersc1

Message:

merge from trunk. Preliminary versions of stacktool/warptool updates and regenerate_background.pl script. The warp code is finished and tested, and I still need to get the stack version resolved.

Location:

branches/czw_branch/20120906/psphot

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/psphotBlendFit.c (modified) (7 diffs)

branches/czw_branch/20120906/psphot/src/psphotBlendFit.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 bool psphotBlendFitSetSource(pmSource *source, psImage *IDimage, float sigSigma);
+bool psphotBlendSetSourceSatstar (psImage *image, pmSource *source);
+float InterpolateValues (float X0, float Y0, float X1, float Y1, float X);
 
 // for now, let's store the detections on the readout->analysis for each readout
@@ -256 +258 @@
 
         int TEST_ON = false;
-# if (1)
+# if (0)
 # define TEST_X 653
 # define TEST_Y 466
         if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
             fprintf (stderr, "test object\n");
-            //  psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
+            psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
             TEST_ON = true;
         }
@@ -377 +379 @@
     if (!source) return false;
     if (!source->peak) return false; // XXX how can we have a peak-less source?
+
+# if (0)
+# define TEST_X 653
+# define TEST_Y 466
+        if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+            fprintf (stderr, "test object\n");
+        }
+# undef TEST_X
+# undef TEST_Y
+# endif
+
     if (!source->moments) return false;
     if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
@@ -385 +398 @@
     if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) {
         // find the radius at which we hit something like 0.1*SATURATION
-        return false;
+        psphotBlendSetSourceSatstar (IDimage, source);
+        return true;
     }
 
@@ -431 +445 @@
     float rMxx = 0.5 / PS_SQR(shape.sx);
     float rMyy = 0.5 / PS_SQR(shape.sy);
-    float Sxy = -1. * shape.sxy;    // factor of -1 is included to match the previous window function
-                                    // implementation. XXX: Is this correct?
+    float Sxy = shape.sxy;    // factor of -1 is included to match the previous window function
+    // implementation. XXX: Is this correct?
 
     int ID = source->id;
@@ -443 +457 @@
 
             float z = rMxx * PS_SQR(dX) + rMyy * PS_SQR(dY) + Sxy*dX*dY;
+            if (z > 2.311) continue;
 
             float f = Io*exp(-z);
@@ -453 +468 @@
     return true;
 }
+
+bool psphotBlendSetSourceSatstar (psImage *IDimage, pmSource *source) {
+
+    float Xo = source->satstar->Xo;
+    float Yo = source->satstar->Yo;
+    psVector *logRmodel = source->satstar->logRmodel;
+    psVector *logFmodel = source->satstar->logFmodel;
+
+    float lPeak = logFmodel->data.F32[0];
+    float fPeak = pow(10.0, lPeak);
+    float threshold = 0.1*fPeak;
+    float logThresh = log10(threshold);
+
+    float radius = NAN;
+
+    // what is the radius for a specific peak fraction?
+    for (int i = 1; i < logFmodel->n; i++) {
+        float logF = logFmodel->data.F32[i];
+        if (!isfinite(logF)) continue;
+
+        float flux = pow(10.0, logF);
+        if (flux > threshold) continue;
+        
+        float logF0 = logFmodel->data.F32[i - 1];
+        float logR0 = logRmodel->data.F32[i - 1];
+
+        float logF1 = logFmodel->data.F32[i];
+        float logR1 = logRmodel->data.F32[i];
+
+        float logR = InterpolateValues (logF0, logR0, logF1, logR1, logThresh);
+        radius = pow(10.0, logR);
+        break;
+    }
+
+    if (!isfinite(radius)) {
+        for (int i = logFmodel->n - 1; i >= 0; i--) {
+            if (!isfinite(logFmodel->data.F32[i])) continue;
+            float logR = logRmodel->data.F32[i];
+            radius = pow(10.0, logR);
+            break;
+        }
+    }
+
+    if (!isfinite(radius)) return false;
+
+    int Nx = IDimage->numCols;
+    int Ny = IDimage->numRows;
+
+    // region to mask
+    int minX = PS_MIN(PS_MAX(Xo - radius, 0), Nx - 1);
+    int maxX = PS_MIN(PS_MAX(Xo + radius, 0), Nx - 1);
+    int minY = PS_MIN(PS_MAX(Yo - radius, 0), Ny - 1);
+    int maxY = PS_MIN(PS_MAX(Yo + radius, 0), Ny - 1);
+
+    int ID = source->id;
+
+    fprintf (stderr, "Xo,Yo: %f,%f; radius: %f\n", Xo, Yo, radius);
+
+    for (int iy = minY; iy < maxY; iy++) {
+        for (int ix = minX; ix < maxX; ix++) {
+
+            float dX = (ix - Xo);
+            float dY = (iy - Yo);
+            float R = hypot (dX, dY) ;
+            if (R > radius) continue;
+
+            IDimage->data.S32[iy][ix] = ID;
+        }
+    }
+    return true;
+}
+