Changeset 17153

Timestamp:

Mar 26, 2008, 11:00:19 AM (18 years ago)

Author:

eugene

Message:

adding function to measure 1 sigma contour

File:

• branches/eam_branch_20080324/psphot/src/psphotSourceSize.c (modified) (2 diffs)

branches/eam_branch_20080324/psphot/src/psphotSourceSize.c

@@ -55 +55 @@
         if (!keep) continue;
 
-        // need a reasonably fast way to follow the psf R = 1 annulus
-
-        // save the PSF model from the Ensemble fit
-        pmModel *PSF = source->modelPSF;
-
-        // convert model to polar terms
-        psEllipseShape shape;
-        psF32 *PAR = PSF->params->data.F32;
-        shape.sx   = PAR[PM_PAR_SXX] / M_SQRT2;
-        shape.sy   = PAR[PM_PAR_SYY] / M_SQRT2;
-        shape.sxy  = PAR[PM_PAR_SXY];
-
-        psEllipseShapeToAxes (shape, 20.0);
-        
-
-        for (theta = 0.0; theta < 360.0; theta += 10.0) {
-          
-
+        // measure the flux at the 1 sigma contour
+        psVector *contour = psphotModelContour (source->pixels, source->maskObj, source->modelPSF, 1.0);
 
         // XXX for now, just skip any masked pixels
@@ -198 +182 @@
  */
 
+
+// given the PSF ellipse parameters, navigate around the 1sigma contour
+// XXX return the Nsigma total deviation?  
+// XXX return just the sum around the contour?  
+// this is measure on the residual image - should we ignore negative deviations?
+psVector *psphotModelContour (psImage *image, psImage *mask, pmModel *model, float Ro) {
+
+    psF32 *PAR = model->params->data.F32;
+    psVector *contour = psVectorAllocEmpty (50, PS_TYPE_F32);
+    int nPts = 0;
+
+    // Ro = (x / SXX)^2 + (y / SYY)^2 + x y SXY;
+    // y^2 (1/SYY^2) + y (x SXY) + (x / SXX)^2 - Ro = 0;
+    // y = [-(x SXY) +/- sqrt ((x SXY)^2 - 4 (1/SYY^2) ((x/SXX)^2 - Ro))] * [SYY^2 / 2];
+    // y = [-B +/- sqrt (B^2 - 4 A C)] / [2 A];
+
+    // min/max value of x is where T -> 0
+    // solve this for x2:
+    float Q = Ro * PS_SQR(PAR[PM_PAR_SXX]) / (1.0 - PS_SQR(PAR[PM_PAR_SXX]*PAR[PM_PAR_SYY]*PAR[PM_PAR_SXY]) / 4.0);
+    if (Q < 0.0) return contour; // ellipse is imaginary
+
+    int xMax = sqrt(Q);
+    int xMin = -1.0*xMax;
+
+    for (int x = MIN; x <= MAX; x++) {
+        float A = PS_SQR (1.0 / PAR[PM_PAR_SYY]);
+        float B = x * PAR[PM_PAR_SXY];
+        float C = PS_SQR (x / PAR[PM_PAR_SXX]) - Ro;
+
+        float T = PS_SQR(B) - 4*A*C;
+        if (T < 0.0) continue;
+    
+        float yP = (-B + sqrt (T)) / (2.0 * A);
+        float yM = (-B - sqrt (T)) / (2.0 * A);
+
+        int xPix  = x  + PAR[PM_PAR_X0] - source->pixels->col0 + 0.5;
+        int yPixM = yM + PAR[PM_PAR_Y0] - source->pixels->row0 + 0.5;
+        int yPixP = yP + PAR[PM_PAR_Y0] - source->pixels->row0 + 0.5;
+
+        if (xPix < 0) continue;
+        if (xPix >= source->pixels->numCols) continue;
+
+        if ((yPixM >= 0) && (yPixM < source->pixels->numRows)) {
+            if (!mask || !mask->data.U8[xPix][yPixM]) {
+                contour->data.F32[nPts] = image->data.F32[xPix][yPixM];
+                psVectorExtend (contour, 100, 1);
+                nPts++;
+            }
+        }
+        
+        if (yPixM == yPixP) continue;
+
+        if ((yPixP >= 0) && (yPixP < source->pixels->numRows)) {
+            if (!mask || !mask->data.U8[xPix][yPixP]) {
+                contour->data.F32[nPts] = image->data.F32[xPix][yPixP];
+                psVectorExtend (contour, 100, 1);
+                nPts++;
+            }
+        }
+    }   
+
+    return contour;
+}