IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 17, 2014, 12:30:45 PM (12 years ago)
Author:
eugene
Message:

merge changes (from past YEAR) into trunk

Location:
branches/eam_branches/ipp-ops-20130712/psModules
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-ops-20130712/psModules

    • Property svn:mergeinfo deleted

  • branches/eam_branches/ipp-ops-20130712/psModules/src/objects

    • Property svn:ignore

      •  

        old new  
        1212pmSourceIO_CMF_PS1_V1.v1.c
        1313pmSourceIO_CMF_PS1_V4.c
         14pmSourceIO_CMF_PS1_V5.c
        1415pmSourceIO_CMF_PS1_SV1.c
        1516pmSourceIO_CMF_PS1_SV2.c
         17pmSourceIO_CMF_PS1_SV3.c
        1618pmSourceIO_CMF_PS1_DV1.c
        1719pmSourceIO_CMF_PS1_DV2.c
        1820pmSourceIO_CMF_PS1_DV3.c
        19 
         21pmSourceIO_CMF_PS1_DV4.c

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

    r35768 r37066  
    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,
    … …  
    3736#include "pmMoments.h"
    3837#include "pmModelFuncs.h"
     38#include "pmModelClass.h"
    3939#include "pmModel.h"
    4040#include "pmModelUtils.h"
    41 #include "pmModelClass.h"
    4241#include "pmSourceMasks.h"
    4342#include "pmSourceExtendedPars.h"
    4443#include "pmSourceDiffStats.h"
    4544#include "pmSourceSatstar.h"
     45#include "pmSourceLensing.h"
    4646#include "pmSource.h"
    4747#include "pmSourceFitModel.h"
    … …  
    4949#include "pmPSFtry.h"
    5050#include "pmDetections.h"
     51#include "pmModel_CentralPixel.h"
    5152
    5253#include "pmModel_DEV.h"
    … …  
    5859# define PM_MODEL_LIMITS          pmModelLimits_DEV
    5960# define PM_MODEL_RADIUS          pmModelRadius_DEV
     61# define PM_MODEL_SET_FWHM        pmModelSetFWHM_DEV
    6062# define PM_MODEL_FROM_PSF        pmModelFromPSF_DEV
    6163# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_DEV
    … …  
    6365# define PM_MODEL_SET_LIMITS      pmModelSetLimits_DEV
    6466
    65 // f = exp(-z^0.125)
     67// f = exp(-kappa*r^(1/index))
     68// f = exp(-kappa*z^(0.5/index))
     69// index = 4, 0.5/index = 0.125
    6670# define ALPHA 0.125
    67 // # define ALPHA 0.25
    6871
    6972// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
    … …  
    7376// Lax parameter limits
    7477static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0 };
    75 static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 };
     78static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0 };
    7679
    7780// Moderate parameter limits
    … …  
    8689static float *paramsMinUse = paramsMinLax;
    8790static float *paramsMaxUse = paramsMaxLax;
    88 static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5 };
     91static float betaUse[] = { 2, 3e6, 5, 5, 10.0, 10.0, 0.5 };
    8992
    9093static bool limitsApply = true;         // Apply limits?
    91 
    92 # include "pmModel_SERSIC.CP.h"
    9394
    9495psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    109110    psAssert (z >= 0, "do not allow negative z values in model");
    110111
    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 
    119     psF32 radius = hypot(X, Y);
    120     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);
    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(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         }
     112    // for DEV, we can hard-wire kappa(4):
     113    // float index = 4.0;
     114    float kappa = 7.670628;
     115
     116    // r = sqrt(z)
     117    float q = kappa*pow(z,ALPHA);
     118    float f0 = exp(-q);
     119
     120    assert (isfinite(q));
     121
     122    // only worry about the central pixels at most
     123    float radius = hypot(X, Y);
     124    if (radius <= 1.5) {
     125        // Nsub ~ 10*index^2 + 1
     126        psEllipseAxes axes = pmPSF_ModelToAxes(PAR, true); // DEV uses Reff
     127        int Nsub = 2 * ((int)(25 / axes.minor)) + 1;
     128        Nsub = PS_MIN (Nsub, 121);
     129        Nsub = PS_MAX (Nsub, 11);
     130        f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], 4.0, Nsub);
    183131    }   
    184 
    185     psF32 z0 = PAR[PM_PAR_I0]*f1;
    186     psF32 f0 = PAR[PM_PAR_SKY] + z0;
    187 
    188     assert (isfinite(f2));
     132    assert (isfinite(f0));
     133
     134    float f1 = PAR[PM_PAR_I0]*f0;
     135    float f = PAR[PM_PAR_SKY] + f1;
     136
    189137    assert (isfinite(f1));
    190     assert (isfinite(z0));
    191     assert (isfinite(f0));
     138    assert (isfinite(f));
    192139
    193140    if (deriv != NULL) {
    … …  
    195142
    196143        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);
     144        dPAR[PM_PAR_I0]   = +f0;
     145
     146        if (z > 0.01) {
     147          float z1 = f1*kappa*ALPHA*pow(z,ALPHA-1.0);
     148          dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px + Y*PAR[PM_PAR_SXY]);
     149          dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py + X*PAR[PM_PAR_SXY]);
     150          dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX];
     151          dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
     152          dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
     153        } else {
     154          // gradient -> 0 for z -> 0, but has undef form
     155          float z1 = f1*kappa*ALPHA*pow(z,ALPHA);
     156          dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]);
     157          dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]);
     158          dPAR[PM_PAR_SXX]  = +2.0*z1*px/PAR[PM_PAR_SXX]/PAR[PM_PAR_SXX];
     159          dPAR[PM_PAR_SYY]  = +2.0*z1*py/PAR[PM_PAR_SYY]/PAR[PM_PAR_SYY];
     160          dPAR[PM_PAR_SXY]  = -1.0*z1;
     161        }
     162    }
     163    return (f);
    211164}
    212165
    … …  
    292245bool PM_MODEL_GUESS (pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal)
    293246{
     247    // for the moment, we are going to require moments and KronFlux
     248    if (!source->moments) return false;
     249    pmMoments *moments = source->moments;
     250
     251    if (!isfinite(moments->KronFlux)) return false;
     252    if (!isfinite(moments->Mrf)) return false;
     253    if (moments->Mrf < 0.0) return false;
     254
    294255    psF32 *PAR  = model->params->data.F32;
    295256
    … …  
    297258    PAR[PM_PAR_SKY]  = 0.0;
    298259
    299     // set the shape parameters
    300     if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, true)) {
    301       return false;
    302     }
    303 
    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 
    309     // set the model normalization
    310     if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
    311       return false;
    312     }
    313     PAR[PM_PAR_I0] /= Io;
     260    psEllipseMoments emoments;
     261    emoments.x2 = moments->Mxx;
     262    emoments.xy = moments->Mxy;
     263    emoments.y2 = moments->Myy;
     264
     265    // force the axis ratio to be < 20.0
     266    psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0);
     267
     268    if (!isfinite(axes.major)) return false;
     269    if (!isfinite(axes.minor)) return false;
     270    if (!isfinite(axes.theta)) return false;
     271
     272    // Mxx, Mxy, Myy define the elliptical shape, but Mrf defines the width
     273    // the factor of 2.3 comes from Table 1 of Graham and Driver (2005)
     274    float scale = moments->Mrf / axes.major / 2.3;
     275    axes.major *= scale;
     276    axes.minor *= scale;
     277
     278    pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], axes, true);
     279
     280    // psEllipseAxes axes;
     281    // use the code in SetShape here to avoid doing this 2x
     282    // pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
     283
     284    // float norm = pmSersicNorm (4);  // hardwire
     285    float norm = 0.00168012;
     286    float normFlux = 2.0 * M_PI * axes.major * axes.minor * norm;
     287    PAR[PM_PAR_I0] = moments->KronFlux / normFlux;
    314288
    315289    // set the model position
    … …  
    328302    psEllipseAxes axes;
    329303    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);
     304
     305    float norm = 0.00168012;
     306    float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm;
     307
     308    return(flux);
    343309}
    344310
    … …  
    359325    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    360326
    361     // f = Io exp(-z^n) -> z^n = ln(Io/f)
    362     psF64 zn = log(PAR[PM_PAR_I0] / flux);
    363     psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / ALPHA);
     327    // static value for DEV:
     328    float kappa = 7.670628;
     329
     330    // f = Io exp(-kappa*z^n) -> z^n = ln(Io/f) / kappa
     331    psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa;
     332    psF64 radius = axes.major * pow(zn, 0.5 / ALPHA);
    364333
    365334    psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)",
    366335              PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]);
    367336    return (radius);
     337}
     338
     339psF64 PM_MODEL_SET_FWHM (const psVector *params, psF64 sigma) {
     340  return (NAN);
    368341}
    369342
    … …  
    389362    // the 2D PSF model fits polarization terms (E0,E1,E2)
    390363    // convert to shape terms (SXX,SYY,SXY)
    391     bool useReff = pmModelUseReff (modelPSF->type);
     364    bool useReff = modelPSF->class->useReff;
    392365    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
    393366        psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
    … …  
    443416    // convert to shape terms (SXX,SYY,SXY)
    444417    // XXX user-defined value for limit?
    445     bool useReff = pmModelUseReff (model->type);
     418    bool useReff = model->class->useReff;
    446419    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
    447420        psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
Note: See TracChangeset for help on using the changeset viewer.