IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36021


Ignore:
Timestamp:
Aug 25, 2013, 12:08:56 PM (13 years ago)
Author:
eugene
Message:

convert SERSIC to use pmModelCP_SersicSubpix; clean up DEV and EXP versions; add pmPCMMakeModel, pmSersicKappa, pmSersicNorm functions

Location:
branches/eam_branches/ipp-20130711/psModules/src/objects
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r35966 r36021  
    6363# define PM_MODEL_SET_LIMITS      pmModelSetLimits_DEV
    6464
    65 // f = exp(-z^0.125)
     65// f = exp(-kappa*r^(1/index))
     66// f = exp(-kappa*z^(0.5/index))
     67// index = 4, 0.5/index = 0.125
    6668# define ALPHA 0.125
    67 // # define ALPHA 0.25
    6869
    6970// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
    … …  
    8687static float *paramsMinUse = paramsMinLax;
    8788static float *paramsMaxUse = paramsMaxLax;
    88 static float betaUse[] = { 2, 3e6, 5, 5, 3.0, 3.0, 0.5 };
     89static float betaUse[] = { 2, 3e6, 5, 5, 10.0, 10.0, 0.5 };
    8990
    9091static bool limitsApply = true;         // Apply limits?
    … …  
    113114    // r = sqrt(z)
    114115    float q = kappa*pow(z,ALPHA);
    115     psF32 f0 = exp(-q);
    116 
    117     psF32 f1 = PAR[PM_PAR_I0]*f0;
    118     psF32 f = PAR[PM_PAR_SKY] + f1;
     116    float f0 = exp(-q);
    119117
    120118    assert (isfinite(q));
    121119    assert (isfinite(f0));
    122     assert (isfinite(f1));
    123     assert (isfinite(f));
    124 
    125     // only worry about the central 4 pixels at most
    126     // If I use DELTA = 0.2, I'm way off for the total flux
    127     // 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)
    128     // For DELTA = 0.02 & Rmin/Rmaj = 0.25, I'm over flux by 0.15 mags (due to the central pixel)
    129     psF32 radius = hypot(X, Y);
     120
     121    // only worry about the central pixels at most
     122    float radius = hypot(X, Y);
    130123    if (radius <= 1.5) {
    131124        // Nsub ~ 10*index^2 + 1
    … …  
    135128        Nsub = PS_MAX (Nsub, 11);
    136129        f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], 4.0, Nsub);
    137         f1 = PAR[PM_PAR_I0]*f0;
    138         f = PAR[PM_PAR_SKY] + f1;
    139130    }   
     131
     132    float f1 = PAR[PM_PAR_I0]*f0;
     133    float f = PAR[PM_PAR_SKY] + f1;
     134
     135    assert (isfinite(f1));
     136    assert (isfinite(f));
    140137
    141138    if (deriv != NULL) {
    … …  
    304301    // f = Io exp(-kappa*z^n) -> z^n = ln(Io/f) / kappa
    305302    psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa;
    306     psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / ALPHA);
     303    psF64 radius = axes.major * pow(zn, 0.5 / ALPHA);
    307304
    308305    psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)",

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

    r35966 r36021  
    7979static float *paramsMinUse = paramsMinLax;
    8080static float *paramsMaxUse = paramsMaxLax;
    81 static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5};
     81static float betaUse[] = { 2, 3e6, 5, 5, 10.0, 10.0, 0.5};
    8282
    8383static bool limitsApply = true;         // Apply limits?
    … …  
    115115    psF32 f0 = exp(-q);
    116116
    117     psF32 f1 = PAR[PM_PAR_I0]*f0;
    118     psF32 f = PAR[PM_PAR_SKY] + f1;
    119 
    120117    assert (isfinite(q));
    121     assert (isfinite(f0));
    122     assert (isfinite(f1));
    123     assert (isfinite(f));
    124118
    125119    // only worry about the central 4 pixels at most
    … …  
    127121    if (radius <= 1.5) {
    128122        f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], 1.0, 51);
    129         f1 = PAR[PM_PAR_I0]*f0;
    130         f = PAR[PM_PAR_SKY] + f1;
    131     }
     123    }
     124    assert (isfinite(f0));
     125
     126    psF32 f1 = PAR[PM_PAR_I0]*f0;
     127    psF32 f = PAR[PM_PAR_SKY] + f1;
     128
     129    assert (isfinite(f1));
     130    assert (isfinite(f));
    132131
    133132    if (deriv != NULL) {
    … …  
    298297    // f = Io exp(-kappa*sqrt(z)) -> sqrt(z) = ln(Io/f) / kappa
    299298    psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa;
    300     psF64 radius = axes.major * sqrt (2.0) * zn;
     299    psF64 radius = axes.major * zn;
    301300
    302301    psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)",

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

    r35876 r36021  
    2020   * note that a Sersic model is usually defined in terms of R_e, the half-light radius.  This
    2121     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:
     22     Conversion from SXX, SYY, SXY to R_major, R_minor, theta can be done by using:
    2323     shape.sx = SXX / sqrt(2), shape.sy = SYY / sqrt(2), shape.sxy = SXY, then calling
    2424     psEllipseShapeToAxes, and multiplying the values of axes.major, axes.minor by sqrt(2)
    … …  
    5555#include "pmPSFtry.h"
    5656#include "pmDetections.h"
     57#include "pmModel_CentralPixel.h"
    5758
    5859#include "pmModel_SERSIC.h"
    … …  
    7475
    7576// Lax parameter limits
    76 static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0, 0.05 };
    77 static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 4.0 };
     77static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0, 0.1 };
     78static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 1.0 };
    7879
    7980// Moderate parameter limits
    … …  
    8889static float *paramsMinUse = paramsMinLax;
    8990static float *paramsMaxUse = paramsMaxLax;
    90 static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5, 2.0 };
     91static float betaUse[] = { 2, 3e6, 5, 5, 10.0, 10.0, 0.5, 1.0};
    9192
    9293static bool limitsApply = true;         // Apply limits?
    9394
    94 # include "pmModel_SERSIC.CP.h"
     95// # include "pmModel_SERSIC.CP.h"
    9596
    9697psF32 PM_MODEL_FUNC (psVector *deriv,
    … …  
    111112    psAssert (z >= 0, "do not allow negative z values in model");
    112113
    113     float index = 0.5 / PAR[PM_PAR_7];
    114     float par7 = PAR[PM_PAR_7];
    115     float bn = 1.9992*index - 0.3271;
    116     float Io = exp(bn);
    117 
    118     psF32 f2 = bn*pow(z,par7);
    119     psF32 f1 = Io*exp(-f2);
     114    float Sindex = 0.5 / PAR[PM_PAR_7];
     115    float kappa = pmSersicKappa (Sindex);
     116
     117    float q = kappa*pow(z,PAR[PM_PAR_7]);
     118    psF32 f0 = exp(-q);
     119
     120    assert (isfinite(q));
    120121
    121122    psF32 radius = hypot(X, Y);
    122     if (radius < 1.0) {
    123 
    124         // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center
    125 
    126         // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
    127         psEllipseAxes axes;
    128         pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    129 
    130         // get the central pixel flux from the lookup table
    131         float xPix = (axes.major - centralPixelXo) / centralPixeldX;
    132         xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1);
    133         float yPix = (index - centralPixelYo) / centralPixeldY;
    134         yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1);
    135 
    136         // the integral of a Sersic has an analytical form as follows:
    137         float logGamma = lgamma(2.0*index);
    138         float bnFactor = pow(bn, 2.0*index);
    139         float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
    140 
    141         // XXX interpolate to get the value
    142         // XXX for the moment, just integerize
    143         // XXX I need to multiply by the integrated flux to get the flux in the central pixel
    144         float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm;
    145        
    146         float px1 = 1.0 / PAR[PM_PAR_SXX];
    147         float py1 = 1.0 / PAR[PM_PAR_SYY];
    148         float z10 = PS_SQR(px1);
    149         float z01 = PS_SQR(py1);
    150 
    151         // which pixels do we need for this interpolation?
    152         // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...)
    153         if ((X >= 0) && (Y >= 0)) {
    154             float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
    155             float V00 = Vcenter;
    156             float V10 = Io*exp(-bn*pow(z10,par7));
    157             float V01 = Io*exp(-bn*pow(z01,par7));
    158             float V11 = Io*exp(-bn*pow(z11,par7));
    159             f1 = interpolatePixels(V00, V10, V01, V11, X, Y);
    160         }
    161         if ((X < 0) && (Y >= 0)) {
    162             float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
    163             float V00 = Io*exp(-bn*pow(z10,par7));
    164             float V10 = Vcenter;
    165             float V01 = Io*exp(-bn*pow(z11,par7));
    166             float V11 = Io*exp(-bn*pow(z01,par7));
    167             f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y);
    168         }
    169         if ((X >= 0) && (Y < 0)) {
    170             float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
    171             float V00 = Io*exp(-bn*pow(z01,par7));
    172             float V10 = Io*exp(-bn*pow(z11,par7));
    173             float V01 = Vcenter;
    174             float V11 = Io*exp(-bn*pow(z10,par7));
    175             f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y));
    176         }
    177         if ((X < 0) && (Y < 0)) {
    178             float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
    179             float V00 = Io*exp(-bn*pow(z11,par7));
    180             float V10 = Io*exp(-bn*pow(z10,par7));
    181             float V01 = Io*exp(-bn*pow(z01,par7));
    182             float V11 = Vcenter;
    183             f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y));
    184         }
    185     }   
    186 
    187     psF32 z0 = PAR[PM_PAR_I0]*f1;
    188     psF32 f0 = PAR[PM_PAR_SKY] + z0;
    189 
    190     if (!isfinite(z0)) {
    191         fprintf(stderr, "z0 is not finite for %f %f %f %f %f.  Parameters: \n", X, Y, radius, z, f1);
     123    if (radius <= 1.5) {
     124        // Nsub ~ 10*index^2 + 1
     125        psEllipseAxes axes = pmPSF_ModelToAxes(PAR, pmModelClassGetType ("PS_MODEL_SERSIC"));
     126        int Nsub = 2 * ((int)(6.0*Sindex / axes.minor)) + 1;
     127        Nsub = PS_MIN (Nsub, 121);
     128        Nsub = PS_MAX (Nsub, 11);
     129        f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], Sindex, Nsub);
     130    }
     131    if (!isfinite(f0)) {
     132        fprintf(stderr, "f0 is not finite for %f %f %f %f %f.  Parameters: \n", X, Y, radius, z, q);
    192133        fprintf(stderr, "%f %f %f %f %f %f %f %f\n", PAR[0], PAR[1], PAR[2], PAR[3], PAR[4],
    193134            PAR[5], PAR[6], PAR[7]);
    194135    }
    195 
    196     assert (isfinite(f2));
     136    assert (isfinite(f0));
     137
     138    psF32 f1 = PAR[PM_PAR_I0]*f0;
     139    psF32 f = PAR[PM_PAR_SKY] + f1;
     140
    197141    assert (isfinite(f1));
    198     assert (isfinite(z0));
    199     assert (isfinite(f0));
     142    assert (isfinite(f));
    200143
    201144    if (deriv != NULL) {
    … …  
    203146
    204147        dPAR[PM_PAR_SKY]  = +1.0;
    205         dPAR[PM_PAR_I0]   = +f1;
    206 
    207         // gradient is infinite for z = 0; saturate at z = 0.01
    208         psF32 z1 = (z < 0.01) ? z0*bn*par7*pow(0.01,par7 - 1.0) : z0*bn*par7*pow(z,par7 - 1.0);
    209 
    210         dPAR[PM_PAR_7]    = (z < 0.01) ? -z0*pow(0.01,par7)*log(0.01) : -z0*f2*log(z);
    211         dPAR[PM_PAR_7]   *= 3.0;
    212 
    213         assert (isfinite(z1));
     148        dPAR[PM_PAR_I0]   = +f0;
     149
     150        if (z > 0.01) {
     151          float z1 = f1*kappa*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7]-1.0);
     152          dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px + Y*PAR[PM_PAR_SXY]);
     153          dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py + X*PAR[PM_PAR_SXY]);
     154          dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX];
     155          dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
     156          dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
     157          dPAR[PM_PAR_7]    = -1.0*f1*q*log(z);
     158        } else {
     159          // gradient -> 0 for z -> 0, but has undef form
     160          float z1 = f1*kappa*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7]);
     161          dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]);
     162          dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]);
     163          dPAR[PM_PAR_SXX]  = +2.0*z1*px/PAR[PM_PAR_SXX]/PAR[PM_PAR_SXX];
     164          dPAR[PM_PAR_SYY]  = +2.0*z1*py/PAR[PM_PAR_SYY]/PAR[PM_PAR_SYY];
     165          dPAR[PM_PAR_SXY]  = -1.0*z1;
     166          dPAR[PM_PAR_7]    = -1.0*f1*q*log(z + 0.0001);
     167        }
    214168        assert (isfinite(dPAR[PM_PAR_7]));
    215 
    216         dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]);
    217         dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]);
    218         dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX]; // XXX : increase drag?
    219         dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
    220         dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
    221     }
    222     return (f0);
     169    }
     170    return (f);
    223171}
    224172
    … …  
    370318    psEllipseAxes axes;
    371319    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    372     float AspectRatio = axes.minor / axes.major;
    373 
    374     float index = 0.5 / PAR[PM_PAR_7];
    375     float bn = 1.9992*index - 0.3271;
    376 
    377     // the integral of a Sersic has an analytical form as follows:
    378     float logGamma = lgamma(2.0*index);
    379     float bnFactor = pow(bn, 2.0*index);
    380     float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
    381    
    382     psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio;
    383 
    384     return(Flux);
     320
     321    float Sindex = 0.5 / PAR[PM_PAR_7];
     322    float norm = pmSersicNorm (Sindex);
     323
     324    float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm;
     325
     326    return(flux);
    385327}
    386328
    … …  
    401343    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
    402344
     345    float Sindex = 0.5 / PAR[PM_PAR_7];
     346    float kappa = pmSersicKappa (Sindex);
     347
    403348    // f = Io exp(-z^n) -> z^n = ln(Io/f)
    404     psF64 zn = log(PAR[PM_PAR_I0] / flux);
    405     psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / PAR[PM_PAR_7]);
     349    psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa;
     350    psF64 radius = axes.major * pow(zn, Sindex);
    406351
    407352    psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f), par 7 = %f",

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

    r35967 r36021  
    723723}
    724724
     725float pmSersicKappa (float Sindex) {
     726    // this function is empirically derived from a fit to data for Sindex = 0.5 - 5.5
     727    // constrain Sindex or kappa?
     728    float kappa = -0.275552 + 1.972625*Sindex + 0.003487 * PS_SQR(Sindex);
     729    return kappa;
     730}
     731
     732float pmSersicNorm (float Sindex) {
     733
     734    float C0 = NAN;
     735    float C1 = NAN;
     736    float C2 = NAN;
     737
     738    // y = 0.201545 x^0 -0.950965 x^1 -0.315248 x^2
     739    // y = 0.402084 x^0 -1.357775 x^1 -0.105102 x^2
     740    // y = 0.619093 x^0 -1.591674 x^1 -0.041576 x^2
     741    // y = 0.770263 x^0 -1.696421 x^1 -0.023363 x^2
     742    // y = 0.885891 x^0 -1.755684 x^1 -0.015753 x^2
     743
     744    if ((Sindex >= 0.0) && (Sindex < 1.0)) {
     745        C0 = 0.201545; C1 = -0.950965; C2 = -0.315248;
     746        // y = 0.201545 x^0 -0.950965 x^1 -0.315248 x^2
     747    }
     748    if ((Sindex >= 1.0) && (Sindex < 2.0)) {
     749        C0 = 0.402084; C1 = -1.357775; C2 = -0.105102;
     750        // y = 0.402084 x^0 -1.357775 x^1 -0.105102 x^2
     751    }
     752    if ((Sindex >= 2.0) && (Sindex < 3.0)) {
     753        C0 = 0.619093; C1 = -1.591674; C2 = -0.041576;
     754        // y = 0.619093 x^0 -1.591674 x^1 -0.041576 x^2
     755    }
     756    if ((Sindex >= 3.0) && (Sindex < 4.0)) {
     757        C0 = 0.770263; C1 = -1.696421; C2 = -0.023363;
     758        // y = 0.770263 x^0 -1.696421 x^1 -0.023363 x^2
     759    }
     760    if ((Sindex >= 4.0) && (Sindex < 5.5)) {
     761        C0 = 0.885891; C1 = -1.755684; C2 = -0.015753;
     762        // y = 0.885891 x^0 -1.755684 x^1 -0.015753 x^2
     763    }
     764
     765    if (isnan(C0)) return NAN;
     766
     767    float lnorm = C0 + C1*Sindex + C2*Sindex*Sindex;
     768    float norm = exp(lnorm);
     769    return norm;
     770}
     771
    725772# if (0)
    726773// create a vector containing only the unique entries in the input vector

  • branches/eam_branches/ipp-20130711/psModules/src/objects/pmModel_CentralPixel.h

    r35967 r36021  
    6464float pmModelCP_SersicSubpix (float dx, float dy, float Rxx, float Rxy, float Ryy, float Sindex, int Nsub);
    6565
     66float pmSersicKappa (float Sindex);
     67float pmSersicNorm (float Sindex);
     68
    6669# endif

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

    r35961 r36021  
    261261        // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
    262262        psVectorInit (constraint->paramMask, 0);
    263         constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
     263        // constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
    264264        if (params->n == 7) {
     265            nParams = params->n;
     266        } else {
    265267            nParams = params->n - 1;
    266         } else {
    267             nParams = params->n - 2;
    268268            constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1;
    269269        }
    … …  
    488488}
    489489
     490// construct a realization of the source model
     491bool pmPCMMakeModel (pmSource *source, pmModel *model, psImageMaskType maskVal, int psfSize) {
     492
     493    PS_ASSERT_PTR_NON_NULL(source, false);
     494
     495    // if we already have a cached image, re-use that memory
     496    source->modelFlux = psImageCopy (source->modelFlux, source->pixels, PS_TYPE_F32);
     497    psImageInit (source->modelFlux, 0.0);
     498
     499    // modelFlux always has unity normalization (I0 = 1.0)
     500    pmModelAdd (source->modelFlux, source->maskObj, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
     501
     502    // convolve the model image with the PSF
     503    if (USE_1D_GAUSS) {
     504        // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
     505        // * the model flux is not masked
     506        // * threading takes place above this level
     507       
     508        // define the Gauss parameters from the psf
     509        pmModel *modelPSF = source->modelPSF;
     510        psAssert (modelPSF, "psf model must be defined");
     511   
     512        psEllipseAxes axes;
     513        bool useReff = pmModelUseReff (modelPSF->type);
     514        psF32 *PAR = modelPSF->params->data.F32;
     515        pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
     516   
     517        float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
     518        float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
     519
     520        float sigma = 0.5 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
     521        float nsigma = 2.0;
     522
     523        psImageSmooth (source->modelFlux, sigma, nsigma);
     524    } else {
     525        // make sure we save a cached copy of the psf flux
     526        pmSourceCachePSF (source, maskVal);
     527
     528        // convert the cached cached psf model for this source to a psKernel
     529        psKernel *psf = pmPCMkernelFromPSF (source, psfSize);
     530        if (!psf) {
     531            // NOTE: this only happens if the source is too close to an edge
     532            model->flags |= PM_MODEL_STATUS_BADARGS;
     533            return NULL;
     534        }
     535
     536        // XXX not sure if I can place the output on top of the input
     537        psImageConvolveFFT (source->modelFlux, source->modelFlux, NULL, 0, psf);
     538    }
     539    return true;
     540}
     541

  • branches/eam_branches/ipp-20130711/psModules/src/objects/pmPCMdata.h

    r32725 r36021  
    9898bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize);
    9999
     100bool pmPCMMakeModel (pmSource *source, pmModel *model, psImageMaskType maskVal, int psfSize);
     101
    100102/// @}
    101103# endif /* PM_PCM_DATA_H */
Note: See TracChangeset for help on using the changeset viewer.