IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 21, 2011, 10:27:18 AM (15 years ago)
Author:
eugene
Message:

0.5 pixel error in pmPCM fitting; update sersic-like models to use a more-accurate central pixel flux value; update sersic-like models to have consistent Sxx / axes.major relationship (axes.major is effective radius, not sigma); fix flux integration (sersic-like models have an analytical form, other models were not correctly integrating f(r) rdr)

Location:
branches/eam_branches/ipp-20110710/psModules/src/objects
Files:
1 added
10 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20110710/psModules/src/objects

    • Property svn:ignore

      •  

        old new  
        55*.la
        66*.lo
         7pmSourceIO_CMF_PS1_V1.c
         8pmSourceIO_CMF_PS1_V2.c
         9pmSourceIO_CMF_PS1_V3.c

  • branches/eam_branches/ipp-20110710/psModules/src/objects/models/pmModel_DEV.c

    r31451 r32153  
    8989static bool limitsApply = true;         // Apply limits?
    9090
     91# include "pmModel_SERSIC.CP.h"
     92
    9193psF32 PM_MODEL_FUNC (psVector *deriv,
    9294                     const psVector *params,
    … …  
    9496{
    9597    psF32 *PAR = params->data.F32;
    96 
    97     float index = 0.5 / ALPHA;
    98     float bn = 1.9992*index - 0.3271;
    99     float Io = exp(bn);
    10098
    10199    psF32 X  = pixcoord->data.F32[0] - PAR[PM_PAR_XPOS];
    … …  
    105103    psF32 z  = (PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y);
    106104
    107     assert (z >= 0);
     105    // If the elliptical contour is defined in a valid way, we should never trigger this
     106    // assert.  Other models (like PGAUSS) don't use fractional powers, and thus do not have
     107    // NaN values for negative values of z
     108    psAssert (z >= 0, "do not allow negative z values in model");
     109
     110    float index = 0.5 / ALPHA;
     111    float par7 = ALPHA;
     112    float bn = 1.9992*index - 0.3271;
     113    float Io = exp(bn);
    108114
    109115    psF32 f2 = bn*pow(z,ALPHA);
    110116    psF32 f1 = Io*exp(-f2);
     117
     118    psF32 radius = hypot(X, Y);
     119    if (radius < 1.0) {
     120
     121        // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center
     122
     123        // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
     124        psEllipseShape shape;
     125
     126        shape.sx  = PAR[PM_PAR_SXX];
     127        shape.sy  = PAR[PM_PAR_SYY];
     128        shape.sxy = PAR[PM_PAR_SXY];
     129
     130        // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
     131        psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     132
     133        // get the central pixel flux from the lookup table
     134        float xPix = (axes.major - centralPixelXo) / centralPixeldX;
     135        xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1);
     136        float yPix = (index - centralPixelYo) / centralPixeldY;
     137        yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1);
     138
     139        // the integral of a Sersic has an analytical form as follows:
     140        float logGamma = lgamma(2.0*index);
     141        float bnFactor = pow(bn, 2.0*index);
     142        float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
     143
     144        // XXX interpolate to get the value
     145        // XXX for the moment, just integerize
     146        // XXX I need to multiply by the integrated flux to get the flux in the central pixel
     147        float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm;
     148       
     149        float px1 = 1.0 / PAR[PM_PAR_SXX];
     150        float py1 = 1.0 / PAR[PM_PAR_SYY];
     151        float z10 = PS_SQR(px1);
     152        float z01 = PS_SQR(py1);
     153
     154        // which pixels do we need for this interpolation?
     155        // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...)
     156        if ((X >= 0) && (Y >= 0)) {
     157            float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
     158            float V00 = Vcenter;
     159            float V10 = Io*exp(-bn*pow(z10,par7));
     160            float V01 = Io*exp(-bn*pow(z01,par7));
     161            float V11 = Io*exp(-bn*pow(z11,par7));
     162            f1 = interpolatePixels(V00, V10, V01, V11, X, Y);
     163        }
     164        if ((X < 0) && (Y >= 0)) {
     165            float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
     166            float V00 = Io*exp(-bn*pow(z10,par7));
     167            float V10 = Vcenter;
     168            float V01 = Io*exp(-bn*pow(z11,par7));
     169            float V11 = Io*exp(-bn*pow(z01,par7));
     170            f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y);
     171        }
     172        if ((X >= 0) && (Y < 0)) {
     173            float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
     174            float V00 = Io*exp(-bn*pow(z01,par7));
     175            float V10 = Io*exp(-bn*pow(z11,par7));
     176            float V01 = Vcenter;
     177            float V11 = Io*exp(-bn*pow(z10,par7));
     178            f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y));
     179        }
     180        if ((X < 0) && (Y < 0)) {
     181            float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
     182            float V00 = Io*exp(-bn*pow(z11,par7));
     183            float V10 = Io*exp(-bn*pow(z10,par7));
     184            float V01 = Io*exp(-bn*pow(z01,par7));
     185            float V11 = Vcenter;
     186            f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y));
     187        }
     188    }   
     189
    111190    psF32 z0 = PAR[PM_PAR_I0]*f1;
    112191    psF32 f0 = PAR[PM_PAR_SKY] + z0;
    … …  
    120199        psF32 *dPAR = deriv->data.F32;
    121200
     201        dPAR[PM_PAR_SKY]  = +1.0;
     202        dPAR[PM_PAR_I0]   = +2.0*f1; // XXX extra damping..
     203
    122204        // gradient is infinite for z = 0; saturate at z = 0.01
    123205        psF32 z1 = (z < 0.01) ? z0*bn*ALPHA*pow(0.01,ALPHA - 1.0) : z0*bn*ALPHA*pow(z,ALPHA - 1.0);
    124 
    125         dPAR[PM_PAR_SKY]  = +1.0;
    126         dPAR[PM_PAR_I0]   = +2.0*f1;
    127206
    128207        assert (isfinite(z1));
    … …  
    223302
    224303    // set the shape parameters
     304    // XXX adjust this?
    225305    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
    226306      return false;
    … …  
    246326}
    247327
     328// A DeVaucouleur model is equivalent to a Sersic with index = 4.0
    248329psF64 PM_MODEL_FLUX (const psVector *params)
    249330{
    250     float z, norm;
    251331    psEllipseShape shape;
    252332
    253333    psF32 *PAR = params->data.F32;
    254334
    255     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    256     shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
     335    shape.sx  = PAR[PM_PAR_SXX];
     336    shape.sy  = PAR[PM_PAR_SYY];
    257337    shape.sxy = PAR[PM_PAR_SXY];
    258338
    259     // Area is equivalent to 2 pi sigma^2
     339    // for a non-circular DeVaucouleur, the flux of the Rmajor equivalent is scaled by the AspectRatio
    260340    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
    261     psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
    262 
    263     // the area needs to be multiplied by the integral of f(z)
    264     norm = 0.0;
    265 
    266     # define DZ 0.25
    267 
    268     float f0 = 1.0;
    269     float f1, f2;
    270     for (z = DZ; z < 150; z += DZ) {
    271         f1 = exp(-pow(z,ALPHA));
    272         z += DZ;
    273         f2 = exp(-pow(z,ALPHA));
    274         norm += f0 + 4*f1 + f2;
    275         f0 = f2;
    276     }
    277     norm *= DZ / 3.0;
    278 
    279     psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
     341    float AspectRatio = axes.minor / axes.major;
     342
     343    float index = 4.0;
     344    float bn = 1.9992*index - 0.3271;
     345
     346    // the integral of a Sersic has an analytical form as follows:
     347    float logGamma = lgamma(2.0*index);
     348    float bnFactor = pow(bn, 2.0*index);
     349    float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
     350   
     351    psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio;
    280352
    281353    return(Flux);
    … …  
    297369        return (1.0);
    298370
    299     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    300     shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
     371    shape.sx  = PAR[PM_PAR_SXX];
     372    shape.sy  = PAR[PM_PAR_SYY];
    301373    shape.sxy = PAR[PM_PAR_SXY];
    302374

  • branches/eam_branches/ipp-20110710/psModules/src/objects/models/pmModel_EXP.c

    r31451 r32153  
    8181static bool limitsApply = true;         // Apply limits?
    8282
     83# include "pmModel_SERSIC.CP.h"
     84
    8385psF32 PM_MODEL_FUNC (psVector *deriv,
    8486                     const psVector *params,
    … …  
    9395    psF32 z  = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
    9496
    95     // XXX if the elliptical contour is defined in valid way, this step should not be required.
    96     // other models (like PGAUSS) don't use fractional powers, and thus do not have NaN values
    97     // for negative values of z
    98     // XXX use an assert here to force the elliptical parameters to be correctly determined
    99     // if (z < 0) z = 0;
    100     assert (z >= 0);
    101 
    102     psF32 f2 = sqrt(z);
    103     psF32 f1 = exp(-f2);
     97    // If the elliptical contour is defined in a valid way, we should never trigger this
     98    // assert.  Other models (like PGAUSS) don't use fractional powers, and thus do not have
     99    // NaN values for negative values of z
     100    psAssert (z >= 0, "do not allow negative z values in model");
     101
     102    float index = 1.0;
     103    float par7 = 0.5;
     104    float bn = 1.9992*index - 0.3271;
     105    float Io = exp(bn);
     106
     107    psF32 f2 = bn*sqrt(z);
     108    psF32 f1 = Io*exp(-f2);
     109
     110    psF32 radius = hypot(X, Y);
     111    if (radius < 1.0) {
     112
     113        // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center
     114
     115        // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
     116        psEllipseShape shape;
     117
     118        shape.sx  = PAR[PM_PAR_SXX];
     119        shape.sy  = PAR[PM_PAR_SYY];
     120        shape.sxy = PAR[PM_PAR_SXY];
     121
     122        // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
     123        psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     124
     125        // get the central pixel flux from the lookup table
     126        float xPix = (axes.major - centralPixelXo) / centralPixeldX;
     127        xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1);
     128        float yPix = (index - centralPixelYo) / centralPixeldY;
     129        yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1);
     130
     131        // the integral of a Sersic has an analytical form as follows:
     132        float logGamma = lgamma(2.0*index);
     133        float bnFactor = pow(bn, 2.0*index);
     134        float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
     135
     136        // XXX interpolate to get the value
     137        // XXX for the moment, just integerize
     138        // XXX I need to multiply by the integrated flux to get the flux in the central pixel
     139        float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm;
     140       
     141        float px1 = 1.0 / PAR[PM_PAR_SXX];
     142        float py1 = 1.0 / PAR[PM_PAR_SYY];
     143        float z10 = PS_SQR(px1);
     144        float z01 = PS_SQR(py1);
     145
     146        // which pixels do we need for this interpolation?
     147        // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...)
     148        if ((X >= 0) && (Y >= 0)) {
     149            float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
     150            float V00 = Vcenter;
     151            float V10 = Io*exp(-bn*pow(z10,par7));
     152            float V01 = Io*exp(-bn*pow(z01,par7));
     153            float V11 = Io*exp(-bn*pow(z11,par7));
     154            f1 = interpolatePixels(V00, V10, V01, V11, X, Y);
     155        }
     156        if ((X < 0) && (Y >= 0)) {
     157            float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
     158            float V00 = Io*exp(-bn*pow(z10,par7));
     159            float V10 = Vcenter;
     160            float V01 = Io*exp(-bn*pow(z11,par7));
     161            float V11 = Io*exp(-bn*pow(z01,par7));
     162            f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y);
     163        }
     164        if ((X >= 0) && (Y < 0)) {
     165            float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
     166            float V00 = Io*exp(-bn*pow(z01,par7));
     167            float V10 = Io*exp(-bn*pow(z11,par7));
     168            float V01 = Vcenter;
     169            float V11 = Io*exp(-bn*pow(z10,par7));
     170            f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y));
     171        }
     172        if ((X < 0) && (Y < 0)) {
     173            float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
     174            float V00 = Io*exp(-bn*pow(z11,par7));
     175            float V10 = Io*exp(-bn*pow(z10,par7));
     176            float V01 = Io*exp(-bn*pow(z01,par7));
     177            float V11 = Vcenter;
     178            f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y));
     179        }
     180    }   
     181
    104182    psF32 z0 = PAR[PM_PAR_I0]*f1;
    105183    psF32 f0 = PAR[PM_PAR_SKY] + z0;
    … …  
    118196        // gradient is infinite for z = 0; saturate at z = 0.01
    119197        // z1 is -df/dz (the negative sign is canceled by most of dz/dPAR[i]
    120         psF32 z1 = (z < 0.01) ? 0.5*z0/sqrt(0.01) : 0.5*z0/sqrt(z);
     198        psF32 z1 = (z < 0.01) ? 0.5*bn*z0/sqrt(0.01) : 0.5*bn*z0/sqrt(z);
    121199
    122200        // XXX dampen SXX and SYY as in GAUSS?
    … …  
    216294
    217295    // set the shape parameters
     296    // XXX adjust this?
    218297    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
    219298      return false;
    … …  
    233312}
    234313
     314// An exponential model is equivalent to a Sersic with index = 1.0
    235315psF64 PM_MODEL_FLUX (const psVector *params)
    236316{
    237     float z, norm;
    238317    psEllipseShape shape;
    239318
    240319    psF32 *PAR = params->data.F32;
    241320
    242     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    243     shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
     321    shape.sx  = PAR[PM_PAR_SXX];
     322    shape.sy  = PAR[PM_PAR_SYY];
    244323    shape.sxy = PAR[PM_PAR_SXY];
    245324
    246     // Area is equivalent to 2 pi sigma^2
     325    // for a non-circular Exponential, the flux of the Rmajor equivalent is scaled by the AspectRatio
    247326    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
    248     psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
    249 
    250     // the area needs to be multiplied by the integral of f(z) = exp(-sqrt(z)) [0 to infinity]
    251     norm = 0.0;
    252 
    253     # define DZ 0.25
    254 
    255     float f0 = 1.0;
    256     float f1, f2;
    257     for (z = DZ; z < 150; z += DZ) {
    258         f1 = exp(-sqrt(z));
    259         z += DZ;
    260         f2 = exp(-sqrt(z));
    261         norm += f0 + 4*f1 + f2;
    262         f0 = f2;
    263     }
    264     norm *= DZ / 3.0;
    265 
    266     psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
     327    float AspectRatio = axes.minor / axes.major;
     328
     329    float index = 1.0;
     330    float bn = 1.9992*index - 0.3271;
     331
     332    // the integral of a Sersic has an analytical form as follows:
     333    float logGamma = lgamma(2.0*index);
     334    float bnFactor = pow(bn, 2.0*index);
     335    float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
     336   
     337    psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio;
    267338
    268339    return(Flux);
    … …  
    284355        return (1.0);
    285356
    286     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    287     shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
     357    shape.sx  = PAR[PM_PAR_SXX];
     358    shape.sy  = PAR[PM_PAR_SYY];
    288359    shape.sxy = PAR[PM_PAR_SXY];
    289360

  • branches/eam_branches/ipp-20110710/psModules/src/objects/models/pmModel_PGAUSS.c

    r31451 r32153  
    217217}
    218218
    219 psF64 PM_MODEL_FLUX(const psVector *params)
    220 {
    221     float norm, z;
     219// integrate the model to get the full flux
     220psF64 PM_MODEL_FLUX (const psVector *params)
     221{
     222    float z, norm;
    222223    psEllipseShape shape;
    223224
    … …  
    228229    shape.sxy = PAR[PM_PAR_SXY];
    229230
    230     // Area is equivalent to 2 pi sigma^2
    231231    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
    232     psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
    233 
    234     // the area needs to be multiplied by the integral of f(z)
     232    float AspectRatio = axes.minor / axes.major;
     233
     234    // flux = 2 \pi \int f(r) r dr
    235235    norm = 0.0;
    236236
    237     # define DZ 0.25
    238 
    239     float f0 = 1.0;
     237    # define DR 0.25
     238
     239    // f = f(r) * r
     240    float f0 = 0.0;
    240241    float f1, f2;
    241     for (z = DZ; z < 150; z += DZ) {
    242         f1 = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0);
    243         z += DZ;
    244         f2 = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0);
     242    for (float r = DR; r < 150; r += DR) {
     243        z = 0.5 * PS_SQR(r / axes.major);
     244        f1 = r / (1 + z + z*z/2.0 + z*z*z/6.0);
     245        r += DR;
     246        z = 0.5 * PS_SQR(r / axes.major);
     247        f2 = r / (1 + z + z*z/2.0 + z*z*z/6.0);
    245248        norm += f0 + 4*f1 + f2;
    246249        f0 = f2;
    247250    }
    248     norm *= DZ / 3.0;
    249 
    250     psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
     251    norm *= DR / 3.0;
     252
     253    psF64 Flux = PAR[PM_PAR_I0] * norm * 2.0 * M_PI * AspectRatio;
    251254
    252255    return(Flux);

  • branches/eam_branches/ipp-20110710/psModules/src/objects/models/pmModel_PS1_V1.c

    r31670 r32153  
    1414   * PM_PAR_XPOS 2  - X center of object
    1515   * PM_PAR_YPOS 3  - Y center of object
    16    * PM_PAR_SXX 4   - X^2 term of elliptical contour (sqrt(2) / SigmaX)
    17    * PM_PAR_SYY 5   - Y^2 term of elliptical contour (sqrt(2) / SigmaY)
     16   * PM_PAR_SXX 4   - X^2 term of elliptical contour (SigmaX / sqrt(2))
     17   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (SigmaY / sqrt(2))
    1818   * PM_PAR_SXY 6   - X*Y term of elliptical contour
    1919   * PM_PAR_7   7   - amplitude of the linear component (k)
    … …  
    239239}
    240240
     241// integrate the model to get the full flux
    241242psF64 PM_MODEL_FLUX (const psVector *params)
    242243{
    … …  
    250251    shape.sxy = PAR[PM_PAR_SXY];
    251252
    252     // Area is equivalent to 2 pi sigma^2
    253253    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
    254     psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
    255 
    256     // the area needs to be multiplied by the integral of f(z)
     254    float AspectRatio = axes.minor / axes.major;
     255
     256    // flux = 2 \pi \int f(r) r dr
    257257    norm = 0.0;
    258258
    259     # define DZ 0.25
    260 
    261     float f0 = 1.0;
     259    # define DR 0.25
     260
     261    // f = f(r) * r
     262    float f0 = 0.0;
    262263    float f1, f2;
    263     for (z = DZ; z < 150; z += DZ) {
    264         f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
    265         z += DZ;
    266         f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
     264    for (float r = DR; r < 150; r += DR) {
     265        z = 0.5 * PS_SQR(r / axes.major);
     266        f1 = r / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
     267        r += DR;
     268        z = 0.5 * PS_SQR(r / axes.major);
     269        f2 = r / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
    267270        norm += f0 + 4*f1 + f2;
    268271        f0 = f2;
    269272    }
    270     norm *= DZ / 3.0;
    271 
    272     psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
     273    norm *= DR / 3.0;
     274
     275    psF64 Flux = PAR[PM_PAR_I0] * norm * 2.0 * M_PI * AspectRatio;
    273276
    274277    return(Flux);

  • branches/eam_branches/ipp-20110710/psModules/src/objects/models/pmModel_QGAUSS.c

    r31670 r32153  
    240240}
    241241
     242// integrate the model to get the full flux
    242243psF64 PM_MODEL_FLUX (const psVector *params)
    243244{
    … …  
    251252    shape.sxy = PAR[PM_PAR_SXY];
    252253
    253     // Area is equivalent to 2 pi sigma^2
    254254    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
    255     psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
    256 
    257     // the area needs to be multiplied by the integral of f(z)
     255    float AspectRatio = axes.minor / axes.major;
     256
     257    // flux = 2 \pi \int f(r) r dr
    258258    norm = 0.0;
    259259
    260     # define DZ 0.25
    261 
    262     float f0 = 1.0;
     260    # define DR 0.25
     261
     262    // f = f(r) * r
     263    float f0 = 0.0;
    263264    float f1, f2;
    264     for (z = DZ; z < 150; z += DZ) {
    265         f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
    266         z += DZ;
    267         f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
     265    for (float r = DR; r < 150; r += DR) {
     266        z = 0.5 * PS_SQR(r / axes.major);
     267        f1 = r / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
     268        r += DR;
     269        z = 0.5 * PS_SQR(r / axes.major);
     270        f2 = r / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
    268271        norm += f0 + 4*f1 + f2;
    269272        f0 = f2;
    270273    }
    271     norm *= DZ / 3.0;
    272 
    273     psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
     274    norm *= DR / 3.0;
     275
     276    psF64 Flux = PAR[PM_PAR_I0] * norm * 2.0 * M_PI * AspectRatio;
    274277
    275278    return(Flux);

  • branches/eam_branches/ipp-20110710/psModules/src/objects/models/pmModel_RGAUSS.c

    r31451 r32153  
    229229}
    230230
     231// integrate the model to get the full flux
    231232psF64 PM_MODEL_FLUX (const psVector *params)
    232233{
    233     float norm, z;
     234    float z, norm;
    234235    psEllipseShape shape;
    235236
    … …  
    240241    shape.sxy = PAR[PM_PAR_SXY];
    241242
    242     // Area is equivalent to 2 pi sigma^2
    243243    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
    244     psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
    245 
    246     // the area needs to be multiplied by the integral of f(z)
     244    float AspectRatio = axes.minor / axes.major;
     245
     246    // flux = 2 \pi \int f(r) r dr
    247247    norm = 0.0;
    248248
    249     # define DZ 0.25
    250 
    251     float f0 = 1.0;
     249    # define DR 0.25
     250
     251    // f = f(r) * r
     252    float f0 = 0.0;
    252253    float f1, f2;
    253     for (z = DZ; z < 150; z += DZ) {
    254         f1 = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7]));
    255         z += DZ;
    256         f2 = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7]));
     254    for (float r = DR; r < 150; r += DR) {
     255        z = 0.5 * PS_SQR(r / axes.major);
     256        f1 = r / (1 + z + pow(z, PAR[PM_PAR_7]));
     257        r += DR;
     258        z = 0.5 * PS_SQR(r / axes.major);
     259        f2 = r / (1 + z + pow(z, PAR[PM_PAR_7]));
    257260        norm += f0 + 4*f1 + f2;
    258261        f0 = f2;
    259262    }
    260     norm *= DZ / 3.0;
    261 
    262     psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
     263    norm *= DR / 3.0;
     264
     265    psF64 Flux = PAR[PM_PAR_I0] * norm * 2.0 * M_PI * AspectRatio;
    263266
    264267    return(Flux);

  • branches/eam_branches/ipp-20110710/psModules/src/objects/models/pmModel_SERSIC.c

    r31451 r32153  
    1313   * PM_PAR_XPOS 2  - X center of object
    1414   * PM_PAR_YPOS 3  - Y center of object
    15    * PM_PAR_SXX 4   - X^2 term of elliptical contour (sqrt(2) / SigmaX)
    16    * PM_PAR_SYY 5   - Y^2 term of elliptical contour (sqrt(2) / SigmaY)
     15   * PM_PAR_SXX 4   - X^2 term of elliptical contour (SigmaX / sqrt(2))
     16   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (SigmaY / sqrt(2))
    1717   * PM_PAR_SXY 6   - X*Y term of elliptical contour
    1818   * PM_PAR_7   7   - normalized sersic parameter
     19
     20   * note that a Sersic model is usually defined in terms of R_e, the half-light radius.  This
     21     construction does not include a factor of 2 in the X^2 term, etc, like for a Gaussian.
     22     Conversion from SXX, SYY, SXY to R_major, R_minor, theta can be done by using setting:
     23     shape.sx = SXX / sqrt(2), shape.sy = SYY / sqrt(2), shape.sxy = SXY, then calling
     24     psEllipseShapeToAxes, and multiplying the values of axes.major, axes.minor by sqrt(2)
    1925
    2026   note that a standard sersic model uses exp(-K*(z^(1/2n) - 1).  the additional elements (K,
    … …  
    8591static bool limitsApply = true;         // Apply limits?
    8692
     93# include "pmModel_SERSIC.CP.h"
     94
    8795psF32 PM_MODEL_FUNC (psVector *deriv,
    8896                     const psVector *params,
    … …  
    97105    psF32 z  = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
    98106
    99     // XXX if the elliptical contour is defined in valid way, this step should not be required.
    100     // other models (like PGAUSS) don't use fractional powers, and thus do not have NaN values
    101     // for negative values of z
    102     // XXX use an assert here to force the elliptical parameters to be correctly determined
    103     // if (z < 0) z = 0;
    104     assert (z >= 0);
     107    // If the elliptical contour is defined in a valid way, we should never trigger this
     108    // assert.  Other models (like PGAUSS) don't use fractional powers, and thus do not have
     109    // NaN values for negative values of z
     110    psAssert (z >= 0, "do not allow negative z values in model");
    105111
    106112    float index = 0.5 / PAR[PM_PAR_7];
     113    float par7 = PAR[PM_PAR_7];
    107114    float bn = 1.9992*index - 0.3271;
    108115    float Io = exp(bn);
    109116
    110     psF32 f2 = bn*pow(z,PAR[PM_PAR_7]);
     117    psF32 f2 = bn*pow(z,par7);
    111118    psF32 f1 = Io*exp(-f2);
     119
     120    psF32 radius = hypot(X, Y);
     121    if (radius < 1.0) {
     122
     123        // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center
     124
     125        // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
     126        psEllipseShape shape;
     127
     128        shape.sx  = PAR[PM_PAR_SXX];
     129        shape.sy  = PAR[PM_PAR_SYY];
     130        shape.sxy = PAR[PM_PAR_SXY];
     131
     132        // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
     133        psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     134
     135        // get the central pixel flux from the lookup table
     136        float xPix = (axes.major - centralPixelXo) / centralPixeldX;
     137        xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1);
     138        float yPix = (index - centralPixelYo) / centralPixeldY;
     139        yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1);
     140
     141        // the integral of a Sersic has an analytical form as follows:
     142        float logGamma = lgamma(2.0*index);
     143        float bnFactor = pow(bn, 2.0*index);
     144        float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
     145
     146        // XXX interpolate to get the value
     147        // XXX for the moment, just integerize
     148        // XXX I need to multiply by the integrated flux to get the flux in the central pixel
     149        float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm;
     150       
     151        float px1 = 1.0 / PAR[PM_PAR_SXX];
     152        float py1 = 1.0 / PAR[PM_PAR_SYY];
     153        float z10 = PS_SQR(px1);
     154        float z01 = PS_SQR(py1);
     155
     156        // which pixels do we need for this interpolation?
     157        // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...)
     158        if ((X >= 0) && (Y >= 0)) {
     159            float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
     160            float V00 = Vcenter;
     161            float V10 = Io*exp(-bn*pow(z10,par7));
     162            float V01 = Io*exp(-bn*pow(z01,par7));
     163            float V11 = Io*exp(-bn*pow(z11,par7));
     164            f1 = interpolatePixels(V00, V10, V01, V11, X, Y);
     165        }
     166        if ((X < 0) && (Y >= 0)) {
     167            float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
     168            float V00 = Io*exp(-bn*pow(z10,par7));
     169            float V10 = Vcenter;
     170            float V01 = Io*exp(-bn*pow(z11,par7));
     171            float V11 = Io*exp(-bn*pow(z01,par7));
     172            f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y);
     173        }
     174        if ((X >= 0) && (Y < 0)) {
     175            float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
     176            float V00 = Io*exp(-bn*pow(z01,par7));
     177            float V10 = Io*exp(-bn*pow(z11,par7));
     178            float V01 = Vcenter;
     179            float V11 = Io*exp(-bn*pow(z10,par7));
     180            f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y));
     181        }
     182        if ((X < 0) && (Y < 0)) {
     183            float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
     184            float V00 = Io*exp(-bn*pow(z11,par7));
     185            float V10 = Io*exp(-bn*pow(z10,par7));
     186            float V01 = Io*exp(-bn*pow(z01,par7));
     187            float V11 = Vcenter;
     188            f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y));
     189        }
     190    }   
     191
    112192    psF32 z0 = PAR[PM_PAR_I0]*f1;
    113193    psF32 f0 = PAR[PM_PAR_SKY] + z0;
    … …  
    121201        psF32 *dPAR = deriv->data.F32;
    122202
    123         // gradient is infinite for z = 0; saturate at z = 0.01
    124         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);
    125 
    126203        dPAR[PM_PAR_SKY]  = +1.0;
    127204        dPAR[PM_PAR_I0]   = +f1;
    128         dPAR[PM_PAR_7]    = (z < 0.01) ? -z0*pow(0.01,PAR[PM_PAR_7])*log(0.01) : -z0*f2*log(z);
     205
     206        // gradient is infinite for z = 0; saturate at z = 0.01
     207        psF32 z1 = (z < 0.01) ? z0*bn*par7*pow(0.01,par7 - 1.0) : z0*bn*par7*pow(z,par7 - 1.0);
     208
     209        dPAR[PM_PAR_7]    = (z < 0.01) ? -z0*pow(0.01,par7)*log(0.01) : -z0*f2*log(z);
    129210        dPAR[PM_PAR_7]   *= 3.0;
    130211
    … …  
    269350    float Io = exp(0.5*bn);
    270351
    271     // XXX do we need this factor of sqrt(2)?
    272     // float Sxx = PS_MAX(0.5, M_SQRT2*shape.sx);
    273     // float Syy = PS_MAX(0.5, M_SQRT2*shape.sy);
    274 
    275352    float Sxx = PS_MAX(0.5, shape.sx);
    276353    float Syy = PS_MAX(0.5, shape.sy);
    … …  
    294371}
    295372
     373// A sersic model has an integral flux which can be analytically represented
    296374psF64 PM_MODEL_FLUX (const psVector *params)
    297375{
    298     float z, norm;
    299376    psEllipseShape shape;
    300377
    301378    psF32 *PAR = params->data.F32;
    302379
    303     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    304     shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
     380    shape.sx  = PAR[PM_PAR_SXX];
     381    shape.sy  = PAR[PM_PAR_SYY];
    305382    shape.sxy = PAR[PM_PAR_SXY];
    306383
    307     // Area is equivalent to 2 pi sigma^2
     384    // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
    308385    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
    309     psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
    310 
    311     // the area needs to be multiplied by the integral of f(z)
    312     norm = 0.0;
    313 
    314     # define DZ 0.25
    315 
    316     float f0 = 1.0;
    317     float f1, f2;
    318     for (z = DZ; z < 150; z += DZ) {
    319         // f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
    320         f1 = exp(-pow(z,PAR[PM_PAR_7]));
    321         z += DZ;
    322         f2 = exp(-pow(z,PAR[PM_PAR_7]));
    323         norm += f0 + 4*f1 + f2;
    324         f0 = f2;
    325     }
    326     norm *= DZ / 3.0;
    327 
    328     psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
     386    float AspectRatio = axes.minor / axes.major;
     387
     388    float index = 0.5 / PAR[PM_PAR_7];
     389    float bn = 1.9992*index - 0.3271;
     390
     391    // the integral of a Sersic has an analytical form as follows:
     392    float logGamma = lgamma(2.0*index);
     393    float bnFactor = pow(bn, 2.0*index);
     394    float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
     395   
     396    psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio;
    329397
    330398    return(Flux);
    … …  
    346414        return (1.0);
    347415
    348     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
    349     shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
     416    shape.sx  = PAR[PM_PAR_SXX];
     417    shape.sy  = PAR[PM_PAR_SYY];
    350418    shape.sxy = PAR[PM_PAR_SXY];
    351419

  • branches/eam_branches/ipp-20110710/psModules/src/objects/pmPCM_MinimizeChisq.c

    r31926 r32153  
    290290
    291291            // Convert i/j to image space:
    292             coord->data.F32[0] = (psF32) (j + source->pixels->col0);
    293             coord->data.F32[1] = (psF32) (i + source->pixels->row0);
     292            coord->data.F32[0] = (psF32) (j + 0.5 + source->pixels->col0);
     293            coord->data.F32[1] = (psF32) (i + 0.5 + source->pixels->row0);
    294294
    295295            pcm->modelFlux->data.F32[i][j] = pcm->modelConv->modelFunc (deriv, params, coord);

  • branches/eam_branches/ipp-20110710/psModules/src/objects/pmSourceFitPCM.c

    r31992 r32153  
    137137}
    138138
     139// XXX deprecate this function or merge with the empirical model
    139140bool pmSourceModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal) {
    140141
Note: See TracChangeset for help on using the changeset viewer.