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_DEV.c

    r35876 r35961  
    1616   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (sqrt(2) / SigmaY)
    1717   * PM_PAR_SXY 6   - X*Y term of elliptical contour
    18    * PM_PAR_7   7   - normalized dev parameter
    1918
    2019   note that a standard dev model uses exp(-K*(z^(1/2n) - 1).  the additional elements (K,
    … …  
    8685static float *paramsMinUse = paramsMinLax;
    8786static float *paramsMaxUse = paramsMaxLax;
    88 static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5 };
     87static float betaUse[] = { 2, 3e6, 5, 5, 3.0, 3.0, 0.5 };
    8988
    9089static bool limitsApply = true;         // Apply limits?
    9190
    92 # include "pmModel_SERSIC.CP.h"
     91// # include "pmModel_SERSIC.CP.h"
    9392
    9493psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    109108    psAssert (z >= 0, "do not allow negative z values in model");
    110109
    111     float index = 0.5 / ALPHA;
    112     float par7 = ALPHA;
    113     float bn = 1.9992*index - 0.3271;
    114     float Io = exp(bn);
    115 
    116     psF32 f2 = bn*pow(z,ALPHA);
    117     psF32 f1 = Io*exp(-f2);
    118 
     110    // for DEV, we can hard-wire kappa(4):
     111    // float index = 4.0;
     112    float kappa = 7.670628;
     113
     114    // r = sqrt(z)
     115    float q = kappa*pow(z,ALPHA);
     116    psF32 f0 = exp(-q);
     117
     118    psF32 f1 = PAR[PM_PAR_I0]*f0;
     119    psF32 f = PAR[PM_PAR_SKY] + f1;
     120
     121    assert (isfinite(q));
     122    assert (isfinite(f0));
     123    assert (isfinite(f1));
     124    assert (isfinite(f));
     125
     126    // only worry about the central 4 pixels at most
     127    // If I use DELTA = 0.2, I'm way off for the total flux
     128    // If I use DELTA = 0.02, I'm totally good (but I am under on the total flux for R = 30 by 0.2 mags -- aperture failure)
     129    // For DELTA = 0.02 & Rmin/Rmaj = 0.25, I'm over flux by 0.15 mags (due to the central pixel)
    119130    psF32 radius = hypot(X, Y);
    120131    if (radius < 1.0) {
    121 
    122         // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center
    123 
    124         // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
    125         psEllipseAxes axes;
    126         pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    127 
    128         // get the central pixel flux from the lookup table
    129         float xPix = (axes.major - centralPixelXo) / centralPixeldX;
    130         xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1);
    131         float yPix = (index - centralPixelYo) / centralPixeldY;
    132         yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1);
    133 
    134         // the integral of a Sersic has an analytical form as follows:
    135         float logGamma = lgamma(2.0*index);
    136         float bnFactor = pow(bn, 2.0*index);
    137         float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
    138 
    139         // XXX interpolate to get the value
    140         // XXX for the moment, just integerize
    141         // XXX I need to multiply by the integrated flux to get the flux in the central pixel
    142         float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm;
    143        
    144         float px1 = 1.0 / PAR[PM_PAR_SXX];
    145         float py1 = 1.0 / PAR[PM_PAR_SYY];
    146         float z10 = PS_SQR(px1);
    147         float z01 = PS_SQR(py1);
    148 
    149         // which pixels do we need for this interpolation?
    150         // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...)
    151         if ((X >= 0) && (Y >= 0)) {
    152             float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
    153             float V00 = Vcenter;
    154             float V10 = Io*exp(-bn*pow(z10,par7));
    155             float V01 = Io*exp(-bn*pow(z01,par7));
    156             float V11 = Io*exp(-bn*pow(z11,par7));
    157             f1 = interpolatePixels(V00, V10, V01, V11, X, Y);
     132      // subdivide the central 2,3,4 pixels by Nx,Ny
     133      float Npix = 0.0;
     134      float Fpix = 0.0;
     135      float Xpix = floor(pixcoord->data.F32[0]) - PAR[PM_PAR_XPOS];
     136      float Ypix = floor(pixcoord->data.F32[1]) - PAR[PM_PAR_YPOS];
     137      # define DELTA 0.02
     138      for (float ix = 0.1; ix <= 0.9; ix += DELTA) {
     139        for (float iy = 0.1; iy <= 0.9; iy += DELTA) {
     140          psF32 X  = Xpix + ix;
     141          psF32 Y  = Ypix + iy;
     142          psF32 px = X / PAR[PM_PAR_SXX];
     143          psF32 py = Y / PAR[PM_PAR_SYY];
     144          psF32 z  = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
     145         
     146          // sqrt(z) is r
     147          float q = kappa*pow(z,ALPHA);
     148          psF32 f0 = exp(-q);
     149         
     150          psF32 f1 = PAR[PM_PAR_I0]*f0;
     151          psF32 fx = PAR[PM_PAR_SKY] + f1;
     152          Fpix += fx;
     153          Npix += 1.0;
    158154        }
    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(z10,par7));
    162             float V10 = Vcenter;
    163             float V01 = Io*exp(-bn*pow(z11,par7));
    164             float V11 = Io*exp(-bn*pow(z01,par7));
    165             f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y);
    166         }
    167         if ((X >= 0) && (Y < 0)) {
    168             float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
    169             float V00 = Io*exp(-bn*pow(z01,par7));
    170             float V10 = Io*exp(-bn*pow(z11,par7));
    171             float V01 = Vcenter;
    172             float V11 = Io*exp(-bn*pow(z10,par7));
    173             f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y));
    174         }
    175         if ((X < 0) && (Y < 0)) {
    176             float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
    177             float V00 = Io*exp(-bn*pow(z11,par7));
    178             float V10 = Io*exp(-bn*pow(z10,par7));
    179             float V01 = Io*exp(-bn*pow(z01,par7));
    180             float V11 = Vcenter;
    181             f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y));
    182         }
     155      }
     156      f = Fpix / Npix;
    183157    }   
    184 
    185     psF32 z0 = PAR[PM_PAR_I0]*f1;
    186     psF32 f0 = PAR[PM_PAR_SKY] + z0;
    187 
    188     assert (isfinite(f2));
    189     assert (isfinite(f1));
    190     assert (isfinite(z0));
    191     assert (isfinite(f0));
    192158
    193159    if (deriv != NULL) {
    … …  
    195161
    196162        dPAR[PM_PAR_SKY]  = +1.0;
    197         dPAR[PM_PAR_I0]   = +2.0*f1; // XXX extra damping..
    198 
    199         // gradient is infinite for z = 0; saturate at z = 0.01
    200         psF32 z1 = (z < 0.01) ? z0*bn*ALPHA*pow(0.01,ALPHA - 1.0) : z0*bn*ALPHA*pow(z,ALPHA - 1.0);
    201 
    202         assert (isfinite(z1));
    203 
    204         dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]);
    205         dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]);
    206         dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX];
    207         dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
    208         dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
    209     }
    210     return (f0);
     163        dPAR[PM_PAR_I0]   = +f0;
     164
     165        if (z > 0.01) {
     166          float z1 = f1*kappa*ALPHA*pow(z,ALPHA-1.0);
     167          dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px + Y*PAR[PM_PAR_SXY]);
     168          dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py + X*PAR[PM_PAR_SXY]);
     169          dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX];
     170          dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
     171          dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
     172        } else {
     173          // gradient -> 0 for z -> 0, but has undef form
     174          float z1 = f1*kappa*ALPHA*pow(z,ALPHA);
     175          dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]);
     176          dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]);
     177          dPAR[PM_PAR_SXX]  = +2.0*z1*px/PAR[PM_PAR_SXX]/PAR[PM_PAR_SXX];
     178          dPAR[PM_PAR_SYY]  = +2.0*z1*py/PAR[PM_PAR_SYY]/PAR[PM_PAR_SYY];
     179          dPAR[PM_PAR_SXY]  = -1.0*z1;
     180        }
     181    }
     182    return (f);
    211183}
    212184
    … …  
    302274    }
    303275
    304     // the normalization is modified by the slope
    305     float index = 0.5 / ALPHA;
    306     float bn = 1.9992*index - 0.3271;
    307     float Io = exp(0.5*bn);
    308 
    309276    // set the model normalization
    310277    if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
    311278      return false;
    312279    }
    313     PAR[PM_PAR_I0] /= Io;
    314280
    315281    // set the model position
    … …  
    328294    psEllipseAxes axes;
    329295    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    330     float AspectRatio = axes.minor / axes.major;
    331 
    332     float index = 4.0;
    333     float bn = 1.9992*index - 0.3271;
    334 
    335     // the integral of a Sersic has an analytical form as follows:
    336     float logGamma = lgamma(2.0*index);
    337     float bnFactor = pow(bn, 2.0*index);
    338     float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
    339    
    340     psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio;
    341 
    342     return(Flux);
     296
     297    float norm = 0.00168012;
     298    float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm;
     299
     300    return(flux);
    343301}
    344302
    … …  
    359317    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    360318
    361     // f = Io exp(-z^n) -> z^n = ln(Io/f)
    362     psF64 zn = log(PAR[PM_PAR_I0] / flux);
     319    // static value for DEV:
     320    float kappa = 7.670628;
     321
     322    // f = Io exp(-kappa*z^n) -> z^n = ln(Io/f) / kappa
     323    psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa;
    363324    psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / ALPHA);
    364325
Note: See TracChangeset for help on using the changeset viewer.