IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 15, 2013, 5:56:56 PM (13 years ago)
Author:
eugene
Message:

more work on the central pixel optimizations -- perhaps not needed (not so expensive?); add some interactive support for PCM chisq fitting; EXP and DEV are for the moment using subdivided central pixels, but this is perhaps too slow?; turn on sky fitting for the PCM model fitting

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r35876 r35961  
    8282static bool limitsApply = true;         // Apply limits?
    8383
    84 # include "pmModel_SERSIC.CP.h"
     84// # include "pmModel_SERSIC.CP.h"
     85
     86// the problems I'm having with the SERSIC-like functions are:
     87// 1) making sure I have the right functional form so that PAR[SXX,etc] represent R_eff (half-light radius)
     88// 2) getting the central pixel right
     89// 3) getting the derivaties right.
    8590
    8691psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    101106    psAssert (z >= 0, "do not allow negative z values in model");
    102107
    103     float index = 1.0;
    104     float par7 = 0.5;
    105     float bn = 1.9992*index - 0.3271;
    106     float Io = exp(bn);
    107 
    108     psF32 f2 = bn*sqrt(z);
    109     psF32 f1 = Io*exp(-f2);
    110 
     108    // for EXP, we can hard-wire kappa(1):
     109    // float index = 1.0;
     110    float kappa = 1.70056;
     111
     112    // sqrt(z) is r
     113    float q = kappa*sqrt(z);
     114    psF32 f0 = exp(-q);
     115
     116    psF32 f1 = PAR[PM_PAR_I0]*f0;
     117    psF32 f = PAR[PM_PAR_SKY] + f1;
     118
     119    assert (isfinite(q));
     120    assert (isfinite(f0));
     121    assert (isfinite(f1));
     122    assert (isfinite(f));
     123
     124    // only worry about the central 4 pixels at most
    111125    psF32 radius = hypot(X, Y);
    112126    if (radius < 1.0) {
    113 
    114         // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center
    115 
    116         // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
    117         psEllipseAxes axes;
    118         pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    119 
    120         // get the central pixel flux from the lookup table
    121         float xPix = (axes.major - centralPixelXo) / centralPixeldX;
    122         xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1);
    123         float yPix = (index - centralPixelYo) / centralPixeldY;
    124         yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1);
    125 
    126         // the integral of a Sersic has an analytical form as follows:
    127         float logGamma = lgamma(2.0*index);
    128         float bnFactor = pow(bn, 2.0*index);
    129         float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
    130 
    131         // XXX interpolate to get the value
    132         // XXX for the moment, just integerize
    133         // XXX I need to multiply by the integrated flux to get the flux in the central pixel
    134         float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm;
    135        
    136         float px1 = 1.0 / PAR[PM_PAR_SXX];
    137         float py1 = 1.0 / PAR[PM_PAR_SYY];
    138         float z10 = PS_SQR(px1);
    139         float z01 = PS_SQR(py1);
    140 
    141         // which pixels do we need for this interpolation?
    142         // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...)
    143         if ((X >= 0) && (Y >= 0)) {
    144             float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
    145             float V00 = Vcenter;
    146             float V10 = Io*exp(-bn*pow(z10,par7));
    147             float V01 = Io*exp(-bn*pow(z01,par7));
    148             float V11 = Io*exp(-bn*pow(z11,par7));
    149             f1 = interpolatePixels(V00, V10, V01, V11, X, Y);
     127      // subdivide the central 2,3,4 pixels by Nx,Ny
     128      float Npix = 0.0;
     129      float Fpix = 0.0;
     130      float Xpix = floor(pixcoord->data.F32[0]) - PAR[PM_PAR_XPOS];
     131      float Ypix = floor(pixcoord->data.F32[1]) - PAR[PM_PAR_YPOS];
     132      for (float ix = 0.1; ix < 1.0; ix += 0.2) {
     133        for (float iy = 0.1; iy < 1.0; iy += 0.2) {
     134          psF32 X  = Xpix + ix;
     135          psF32 Y  = Ypix + iy;
     136          psF32 px = X / PAR[PM_PAR_SXX];
     137          psF32 py = Y / PAR[PM_PAR_SYY];
     138          psF32 z  = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
     139         
     140          // sqrt(z) is r
     141          float q = kappa*sqrt(z);
     142          psF32 f0 = exp(-q);
     143         
     144          psF32 f1 = PAR[PM_PAR_I0]*f0;
     145          psF32 fx = PAR[PM_PAR_SKY] + f1;
     146          Fpix += fx;
     147          Npix += 1.0;
    150148        }
    151         if ((X < 0) && (Y >= 0)) {
    152             float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
    153             float V00 = Io*exp(-bn*pow(z10,par7));
    154             float V10 = Vcenter;
    155             float V01 = Io*exp(-bn*pow(z11,par7));
    156             float V11 = Io*exp(-bn*pow(z01,par7));
    157             f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y);
    158         }
    159         if ((X >= 0) && (Y < 0)) {
    160             float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
    161             float V00 = Io*exp(-bn*pow(z01,par7));
    162             float V10 = Io*exp(-bn*pow(z11,par7));
    163             float V01 = Vcenter;
    164             float V11 = Io*exp(-bn*pow(z10,par7));
    165             f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y));
    166         }
    167         if ((X < 0) && (Y < 0)) {
    168             float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
    169             float V00 = Io*exp(-bn*pow(z11,par7));
    170             float V10 = Io*exp(-bn*pow(z10,par7));
    171             float V01 = Io*exp(-bn*pow(z01,par7));
    172             float V11 = Vcenter;
    173             f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y));
    174         }
    175     }   
    176 
    177     psF32 z0 = PAR[PM_PAR_I0]*f1;
    178     psF32 f0 = PAR[PM_PAR_SKY] + z0;
    179 
    180     assert (isfinite(f2));
    181     assert (isfinite(f1));
    182     assert (isfinite(z0));
    183     assert (isfinite(f0));
     149      }
     150      f = Fpix / Npix;
     151    }
    184152
    185153    if (deriv != NULL) {
    … …  
    187155
    188156        dPAR[PM_PAR_SKY]  = +1.0;
    189         dPAR[PM_PAR_I0]   = +f1;
    190 
    191         // gradient is infinite for z = 0; saturate at z = 0.01
    192         // z1 is -df/dz (the negative sign is canceled by most of dz/dPAR[i]
    193         psF32 z1 = (z < 0.01) ? 0.5*bn*z0/sqrt(0.01) : 0.5*bn*z0/sqrt(z);
    194 
    195         // XXX dampen SXX and SYY as in GAUSS?
    196         dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]);
    197         dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]);
    198         dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX];
    199         dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
    200         dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
    201     }
    202     return (f0);
     157        dPAR[PM_PAR_I0]   = +f0;
     158
     159        if (z > 0.01) {
     160          float z1 = 0.5*f1*kappa/sqrt(z);
     161          dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px + Y*PAR[PM_PAR_SXY]);
     162          dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py + X*PAR[PM_PAR_SXY]);
     163          dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX];
     164          dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
     165          dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
     166        } else {
     167          // gradient -> 0 for z -> 0, but has undef form
     168          float z1 = 0.5*f1*kappa;
     169          dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]);
     170          dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]);
     171          dPAR[PM_PAR_SXX]  = +2.0*z1*px/PAR[PM_PAR_SXX]/PAR[PM_PAR_SXX];
     172          dPAR[PM_PAR_SYY]  = +2.0*z1*py/PAR[PM_PAR_SYY]/PAR[PM_PAR_SYY];
     173          dPAR[PM_PAR_SXY]  = -1.0*z1;
     174        }
     175    }
     176    return (f);
    203177}
    204178
    … …  
    314288    psEllipseAxes axes;
    315289    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    316     float AspectRatio = axes.minor / axes.major;
    317 
    318     float index = 1.0;
    319     float bn = 1.9992*index - 0.3271;
    320 
    321     // the integral of a Sersic has an analytical form as follows:
    322     float logGamma = lgamma(2.0*index);
    323     float bnFactor = pow(bn, 2.0*index);
    324     float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
    325    
    326     psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio;
    327 
    328     return(Flux);
     290
     291    // static value for EXP:
     292    float norm = 0.34578; // \int exp(-kappa*sqrt(z)) r dr
     293
     294    float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm;
     295
     296    return(flux);
    329297}
    330298
    … …  
    345313    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    346314
    347     // f = Io exp(-sqrt(z)) -> sqrt(z) = ln(Io/f)
    348     psF64 zn = log(PAR[PM_PAR_I0] / flux);
     315    // static value for EXP:
     316    float kappa = 1.70056;
     317
     318    // f = Io exp(-kappa*sqrt(z)) -> sqrt(z) = ln(Io/f) / kappa
     319    psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa;
    349320    psF64 radius = axes.major * sqrt (2.0) * zn;
    350321
    … …  
    501472    return;
    502473}
     474
     475# if (0)
     476void bilin_inter_function () {
     477        // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
     478        psEllipseAxes axes;
     479        pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
     480
     481        // get the central pixel flux from the lookup table
     482        float xPix = (axes.major - centralPixelXo) / centralPixeldX;
     483        xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1);
     484        float yPix = (index - centralPixelYo) / centralPixeldY;
     485        yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1);
     486
     487        // the integral of a Sersic has an analytical form as follows:
     488        float logGamma = lgamma(2.0*index);
     489        float bnFactor = pow(bn, 2.0*index);
     490        float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
     491
     492        // XXX interpolate to get the value
     493        // XXX for the moment, just integerize
     494        // XXX I need to multiply by the integrated flux to get the flux in the central pixel
     495        float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm;
     496       
     497        float px1 = 1.0 / PAR[PM_PAR_SXX];
     498        float py1 = 1.0 / PAR[PM_PAR_SYY];
     499        float z10 = PS_SQR(px1);
     500        float z01 = PS_SQR(py1);
     501
     502        // which pixels do we need for this interpolation?
     503        // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...)
     504        if ((X >= 0) && (Y >= 0)) {
     505            float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
     506            float V00 = Vcenter;
     507            float V10 = Io*exp(-bn*pow(z10,par7));
     508            float V01 = Io*exp(-bn*pow(z01,par7));
     509            float V11 = Io*exp(-bn*pow(z11,par7));
     510            f1 = interpolatePixels(V00, V10, V01, V11, X, Y);
     511        }
     512        if ((X < 0) && (Y >= 0)) {
     513            float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
     514            float V00 = Io*exp(-bn*pow(z10,par7));
     515            float V10 = Vcenter;
     516            float V01 = Io*exp(-bn*pow(z11,par7));
     517            float V11 = Io*exp(-bn*pow(z01,par7));
     518            f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y);
     519        }
     520        if ((X >= 0) && (Y < 0)) {
     521            float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
     522            float V00 = Io*exp(-bn*pow(z01,par7));
     523            float V10 = Io*exp(-bn*pow(z11,par7));
     524            float V01 = Vcenter;
     525            float V11 = Io*exp(-bn*pow(z10,par7));
     526            f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y));
     527        }
     528        if ((X < 0) && (Y < 0)) {
     529            float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
     530            float V00 = Io*exp(-bn*pow(z11,par7));
     531            float V10 = Io*exp(-bn*pow(z10,par7));
     532            float V01 = Io*exp(-bn*pow(z01,par7));
     533            float V11 = Vcenter;
     534            f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y));
     535        }
     536}
     537# endif
Note: See TracChangeset for help on using the changeset viewer.