Changeset 34409

Timestamp:

Sep 6, 2012, 6:40:26 AM (14 years ago)

Author:

eugene

Message:

fixing model id image, adding satstars

File:

• branches/eam_branches/ipp-20120905/psphot/src/psphotBlendFit.c (modified) (5 diffs)

branches/eam_branches/ipp-20120905/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
@@ -385 +387 @@
     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 +434 @@
     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 +446 @@
 
             float z = rMxx * PS_SQR(dX) + rMyy * PS_SQR(dY) + Sxy*dX*dY;
+            if (z > 2.311) continue;
 
             float f = Io*exp(-z);
@@ -453 +457 @@
     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 = 0; i < logFmodel->n; i++) {
+        float logF = logFmodel->data.F32[i];
+        float flux = pow(10.0, logF);
+        if (flux > threshold) continue;
+        if (i == 0) 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;
+
+    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;
+}
+