Changeset 32347 for trunk/psModules/src/objects/models/pmModel_SERSIC.c

Timestamp:

Sep 6, 2011, 1:02:53 PM (15 years ago)

Author:

eugene

Message:

• add concept of saddlePoints to peaks (not actually used in the end)
• add tmp flags to mark sources for analysis or not in psphotStack
• autocode the pmSourceIO_CMF_PS1_* functions
• use 1D gauss approx for convolution in PCM fitting
• added pmSourceExtFitPars (not actually used)
• in model guess, use 1st radial moments to define size (if it exists)
• include PSF_INST_MAG, AP_MAG, KRON_MAG in xfit output
• fix the position for extended source fits (avoid instability)
• Sersic-like models (incl. Exp and Dev) use Reff, not sigma; conversion tools need to respect this
• only use a single pass on the centroid (unwindowed, but limited to 1.5 sigma radius) - this avoids moving the centroid because of nearby neighbors
• use symmetrical averaging (geometric mean) to calculated 1st radial moment (and avoid neighbor biases), do not use symm. averaging for the flux
• fix the integration of the sersic, pgauss, and related model functions.
• fix the central pixel to have the full flux for sersic-like models (interpolated value)

Location:

trunk/psModules/src/objects

Files:

• . (modified) (1 prop)
• models/pmModel_SERSIC.c (modified) (7 diffs)

trunk/psModules/src/objects

@@ -5 +5 @@
 *.la
 *.lo
+pmSourceIO_CMF_PS1_V1.c
+pmSourceIO_CMF_PS1_V2.c
+pmSourceIO_CMF_PS1_V3.c

@@ -5 +5 @@
 *.la
 *.lo
+pmSourceIO_CMF_PS1_V1.c
+pmSourceIO_CMF_PS1_V2.c
+pmSourceIO_CMF_PS1_V3.c

trunk/psModules/src/objects/models/pmModel_SERSIC.c

@@ -13 +13 @@
    * PM_PAR_XPOS 2  - X center of object
    * PM_PAR_YPOS 3  - Y center of object
-   * PM_PAR_SXX 4   - X^2 term of elliptical contour (sqrt(2) / SigmaX)
-   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (sqrt(2) / SigmaY)
+   * PM_PAR_SXX 4   - X^2 term of elliptical contour (SigmaX / sqrt(2))
+   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (SigmaY / sqrt(2))
    * PM_PAR_SXY 6   - X*Y term of elliptical contour
    * PM_PAR_7   7   - normalized sersic parameter
+
+   * note that a Sersic model is usually defined in terms of R_e, the half-light radius.  This
+     construction does not include a factor of 2 in the X^2 term, etc, like for a Gaussian.
+     Conversion from SXX, SYY, SXY to R_major, R_minor, theta can be done by using setting:
+     shape.sx = SXX / sqrt(2), shape.sy = SYY / sqrt(2), shape.sxy = SXY, then calling
+     psEllipseShapeToAxes, and multiplying the values of axes.major, axes.minor by sqrt(2)
 
    note that a standard sersic model uses exp(-K*(z^(1/2n) - 1).  the additional elements (K,
@@ -85 +91 @@
 static bool limitsApply = true;         // Apply limits?
 
+# include "pmModel_SERSIC.CP.h"
+
 psF32 PM_MODEL_FUNC (psVector *deriv,
                      const psVector *params,
@@ -97 +105 @@
     psF32 z  = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
 
-    // XXX if the elliptical contour is defined in valid way, this step should not be required.
-    // other models (like PGAUSS) don't use fractional powers, and thus do not have NaN values
-    // for negative values of z
-    // XXX use an assert here to force the elliptical parameters to be correctly determined
-    // if (z < 0) z = 0;
-    assert (z >= 0);
+    // If the elliptical contour is defined in a valid way, we should never trigger this
+    // assert.  Other models (like PGAUSS) don't use fractional powers, and thus do not have
+    // NaN values for negative values of z
+    psAssert (z >= 0, "do not allow negative z values in model");
 
     float index = 0.5 / PAR[PM_PAR_7];
+    float par7 = PAR[PM_PAR_7];
     float bn = 1.9992*index - 0.3271;
     float Io = exp(bn);
 
-    psF32 f2 = bn*pow(z,PAR[PM_PAR_7]);
+    psF32 f2 = bn*pow(z,par7);
     psF32 f1 = Io*exp(-f2);
+
+    psF32 radius = hypot(X, Y);
+    if (radius < 1.0) {
+
+        // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center
+
+        // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
+        psEllipseShape shape;
+
+        shape.sx  = PAR[PM_PAR_SXX];
+        shape.sy  = PAR[PM_PAR_SYY];
+        shape.sxy = PAR[PM_PAR_SXY];
+
+        // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
+        psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+
+        // get the central pixel flux from the lookup table
+        float xPix = (axes.major - centralPixelXo) / centralPixeldX;
+        xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1);
+        float yPix = (index - centralPixelYo) / centralPixeldY;
+        yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1);
+
+        // the integral of a Sersic has an analytical form as follows:
+        float logGamma = lgamma(2.0*index);
+        float bnFactor = pow(bn, 2.0*index);
+        float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
+
+        // XXX interpolate to get the value
+        // XXX for the moment, just integerize
+        // XXX I need to multiply by the integrated flux to get the flux in the central pixel
+        float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm;
+        
+        float px1 = 1.0 / PAR[PM_PAR_SXX];
+        float py1 = 1.0 / PAR[PM_PAR_SYY];
+        float z10 = PS_SQR(px1);
+        float z01 = PS_SQR(py1);
+
+        // which pixels do we need for this interpolation?
+        // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...)
+        if ((X >= 0) && (Y >= 0)) {
+            float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
+            float V00 = Vcenter;
+            float V10 = Io*exp(-bn*pow(z10,par7));
+            float V01 = Io*exp(-bn*pow(z01,par7));
+            float V11 = Io*exp(-bn*pow(z11,par7));
+            f1 = interpolatePixels(V00, V10, V01, V11, X, Y);
+        }
+        if ((X < 0) && (Y >= 0)) {
+            float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
+            float V00 = Io*exp(-bn*pow(z10,par7));
+            float V10 = Vcenter;
+            float V01 = Io*exp(-bn*pow(z11,par7));
+            float V11 = Io*exp(-bn*pow(z01,par7));
+            f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y);
+        }
+        if ((X >= 0) && (Y < 0)) {
+            float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
+            float V00 = Io*exp(-bn*pow(z01,par7));
+            float V10 = Io*exp(-bn*pow(z11,par7));
+            float V01 = Vcenter;
+            float V11 = Io*exp(-bn*pow(z10,par7));
+            f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y));
+        }
+        if ((X < 0) && (Y < 0)) {
+            float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
+            float V00 = Io*exp(-bn*pow(z11,par7));
+            float V10 = Io*exp(-bn*pow(z10,par7));
+            float V01 = Io*exp(-bn*pow(z01,par7));
+            float V11 = Vcenter;
+            f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y));
+        }
+    }   
+
     psF32 z0 = PAR[PM_PAR_I0]*f1;
     psF32 f0 = PAR[PM_PAR_SKY] + z0;
@@ -121 +201 @@
         psF32 *dPAR = deriv->data.F32;
 
-        // gradient is infinite for z = 0; saturate at z = 0.01
-        psF32 z1 = (z < 0.01) ? z0*bn*PAR[PM_PAR_7]*pow(0.01,PAR[PM_PAR_7] - 1.0) : z0*bn*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7] - 1.0);
-
         dPAR[PM_PAR_SKY]  = +1.0;
         dPAR[PM_PAR_I0]   = +f1;
-        dPAR[PM_PAR_7]    = (z < 0.01) ? -z0*pow(0.01,PAR[PM_PAR_7])*log(0.01) : -z0*f2*log(z);
+
+        // gradient is infinite for z = 0; saturate at z = 0.01
+        psF32 z1 = (z < 0.01) ? z0*bn*par7*pow(0.01,par7 - 1.0) : z0*bn*par7*pow(z,par7 - 1.0);
+
+        dPAR[PM_PAR_7]    = (z < 0.01) ? -z0*pow(0.01,par7)*log(0.01) : -z0*f2*log(z);
         dPAR[PM_PAR_7]   *= 3.0;
 
@@ -269 +350 @@
     float Io = exp(0.5*bn);
 
-    // XXX do we need this factor of sqrt(2)?
-    // float Sxx = PS_MAX(0.5, M_SQRT2*shape.sx);
-    // float Syy = PS_MAX(0.5, M_SQRT2*shape.sy);
-
     float Sxx = PS_MAX(0.5, shape.sx);
     float Syy = PS_MAX(0.5, shape.sy);
@@ -294 +371 @@
 }
 
+// A sersic model has an integral flux which can be analytically represented
 psF64 PM_MODEL_FLUX (const psVector *params)
 {
-    float z, norm;
     psEllipseShape shape;
 
     psF32 *PAR = params->data.F32;
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
+    shape.sx  = PAR[PM_PAR_SXX];
+    shape.sy  = PAR[PM_PAR_SYY];
     shape.sxy = PAR[PM_PAR_SXY];
 
-    // Area is equivalent to 2 pi sigma^2
+    // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
-    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
-
-    // the area needs to be multiplied by the integral of f(z)
-    norm = 0.0;
-
-    # define DZ 0.25
-
-    float f0 = 1.0;
-    float f1, f2;
-    for (z = DZ; z < 150; z += DZ) {
-        // f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
-        f1 = exp(-pow(z,PAR[PM_PAR_7]));
-        z += DZ;
-        f2 = exp(-pow(z,PAR[PM_PAR_7]));
-        norm += f0 + 4*f1 + f2;
-        f0 = f2;
-    }
-    norm *= DZ / 3.0;
-
-    psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
+    float AspectRatio = axes.minor / axes.major;
+
+    float index = 0.5 / PAR[PM_PAR_7];
+    float bn = 1.9992*index - 0.3271;
+
+    // the integral of a Sersic has an analytical form as follows:
+    float logGamma = lgamma(2.0*index);
+    float bnFactor = pow(bn, 2.0*index);
+    float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
+    
+    psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio;
 
     return(Flux);
@@ -346 +414 @@
         return (1.0);
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
+    shape.sx  = PAR[PM_PAR_SXX];
+    shape.sy  = PAR[PM_PAR_SYY];
     shape.sxy = PAR[PM_PAR_SXY];
 

@@ -5 +5 @@
 *.la
 *.lo
+pmSourceIO_CMF_PS1_V1.c
+pmSourceIO_CMF_PS1_V2.c
+pmSourceIO_CMF_PS1_V3.c