Index: trunk/psModules/src/objects/Makefile.am =================================================================== --- trunk/psModules/src/objects/Makefile.am (revision 35925) +++ trunk/psModules/src/objects/Makefile.am (revision 36085) @@ -20,4 +20,5 @@ pmModelClass.c \ pmModelUtils.c \ + pmModel_CentralPixel.c \ pmSource.c \ pmPhotObj.c \ @@ -97,4 +98,5 @@ pmModelClass.h \ pmModelUtils.h \ + pmModel_CentralPixel.h \ pmSource.h \ pmPhotObj.h \ @@ -110,5 +112,5 @@ pmSourceOutputs.h \ pmSourceIO.h \ - pmSourceSatstar.h \ + pmSourceSatstar.h \ pmSourcePlots.h \ pmSourceVisual.h \ Index: trunk/psModules/src/objects/models/pmModel_DEV.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_DEV.c (revision 35925) +++ trunk/psModules/src/objects/models/pmModel_DEV.c (revision 36085) @@ -16,5 +16,4 @@ * PM_PAR_SYY 5 - Y^2 term of elliptical contour (sqrt(2) / SigmaY) * PM_PAR_SXY 6 - X*Y term of elliptical contour - * PM_PAR_7 7 - normalized dev parameter note that a standard dev model uses exp(-K*(z^(1/2n) - 1). the additional elements (K, @@ -49,4 +48,5 @@ #include "pmPSFtry.h" #include "pmDetections.h" +#include "pmModel_CentralPixel.h" #include "pmModel_DEV.h" @@ -63,7 +63,8 @@ # define PM_MODEL_SET_LIMITS pmModelSetLimits_DEV -// f = exp(-z^0.125) +// f = exp(-kappa*r^(1/index)) +// f = exp(-kappa*z^(0.5/index)) +// index = 4, 0.5/index = 0.125 # define ALPHA 0.125 -// # define ALPHA 0.25 // the model is a function of the pixel coordinate (pixcoord[0,1] = x,y) @@ -73,5 +74,5 @@ // Lax parameter limits static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0 }; -static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0 }; // Moderate parameter limits @@ -86,9 +87,7 @@ static float *paramsMinUse = paramsMinLax; static float *paramsMaxUse = paramsMaxLax; -static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5 }; +static float betaUse[] = { 2, 3e6, 5, 5, 10.0, 10.0, 0.5 }; static bool limitsApply = true; // Apply limits? - -# include "pmModel_SERSIC.CP.h" psF32 PM_MODEL_FUNC (psVector *deriv, @@ -109,85 +108,31 @@ psAssert (z >= 0, "do not allow negative z values in model"); - float index = 0.5 / ALPHA; - float par7 = ALPHA; - float bn = 1.9992*index - 0.3271; - float Io = exp(bn); - - psF32 f2 = bn*pow(z,ALPHA); - psF32 f1 = Io*exp(-f2); - - psF32 radius = hypot(X, Y); - if (radius < 1.0) { - - // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center - - // first, use Rmajor and index to find the central pixel flux (fraction of total flux) - psEllipseAxes axes; - pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); - - // get the central pixel flux from the lookup table - float xPix = (axes.major - centralPixelXo) / centralPixeldX; - xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1); - float yPix = (index - centralPixelYo) / centralPixeldY; - yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1); - - // the integral of a Sersic has an analytical form as follows: - float logGamma = lgamma(2.0*index); - float bnFactor = pow(bn, 2.0*index); - float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor; - - // XXX interpolate to get the value - // XXX for the moment, just integerize - // XXX I need to multiply by the integrated flux to get the flux in the central pixel - float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm; - - float px1 = 1.0 / PAR[PM_PAR_SXX]; - float py1 = 1.0 / PAR[PM_PAR_SYY]; - float z10 = PS_SQR(px1); - float z01 = PS_SQR(py1); - - // which pixels do we need for this interpolation? - // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...) - if ((X >= 0) && (Y >= 0)) { - float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive - float V00 = Vcenter; - float V10 = Io*exp(-bn*pow(z10,par7)); - float V01 = Io*exp(-bn*pow(z01,par7)); - float V11 = Io*exp(-bn*pow(z11,par7)); - f1 = interpolatePixels(V00, V10, V01, V11, X, Y); - } - if ((X < 0) && (Y >= 0)) { - float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative - float V00 = Io*exp(-bn*pow(z10,par7)); - float V10 = Vcenter; - float V01 = Io*exp(-bn*pow(z11,par7)); - float V11 = Io*exp(-bn*pow(z01,par7)); - f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y); - } - if ((X >= 0) && (Y < 0)) { - float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative - float V00 = Io*exp(-bn*pow(z01,par7)); - float V10 = Io*exp(-bn*pow(z11,par7)); - float V01 = Vcenter; - float V11 = Io*exp(-bn*pow(z10,par7)); - f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y)); - } - if ((X < 0) && (Y < 0)) { - float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive - float V00 = Io*exp(-bn*pow(z11,par7)); - float V10 = Io*exp(-bn*pow(z10,par7)); - float V01 = Io*exp(-bn*pow(z01,par7)); - float V11 = Vcenter; - f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y)); - } + // for DEV, we can hard-wire kappa(4): + // float index = 4.0; + float kappa = 7.670628; + + // r = sqrt(z) + float q = kappa*pow(z,ALPHA); + float f0 = exp(-q); + + assert (isfinite(q)); + assert (isfinite(f0)); + + // only worry about the central pixels at most + float radius = hypot(X, Y); + if (radius <= 1.5) { + // Nsub ~ 10*index^2 + 1 + psEllipseAxes axes = pmPSF_ModelToAxes(PAR, pmModelClassGetType ("PS_MODEL_DEV")); + int Nsub = 2 * ((int)(25 / axes.minor)) + 1; + Nsub = PS_MIN (Nsub, 121); + Nsub = PS_MAX (Nsub, 11); + f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], 4.0, Nsub); } - psF32 z0 = PAR[PM_PAR_I0]*f1; - psF32 f0 = PAR[PM_PAR_SKY] + z0; - - assert (isfinite(f2)); + float f1 = PAR[PM_PAR_I0]*f0; + float f = PAR[PM_PAR_SKY] + f1; + assert (isfinite(f1)); - assert (isfinite(z0)); - assert (isfinite(f0)); + assert (isfinite(f)); if (deriv != NULL) { @@ -195,18 +140,24 @@ dPAR[PM_PAR_SKY] = +1.0; - dPAR[PM_PAR_I0] = +2.0*f1; // XXX extra damping.. - - // gradient is infinite for z = 0; saturate at z = 0.01 - psF32 z1 = (z < 0.01) ? z0*bn*ALPHA*pow(0.01,ALPHA - 1.0) : z0*bn*ALPHA*pow(z,ALPHA - 1.0); - - assert (isfinite(z1)); - - dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]); - dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]); - dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; - dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; - dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; - } - return (f0); + dPAR[PM_PAR_I0] = +f0; + + if (z > 0.01) { + float z1 = f1*kappa*ALPHA*pow(z,ALPHA-1.0); + dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px + Y*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py + X*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; + } else { + // gradient -> 0 for z -> 0, but has undef form + float z1 = f1*kappa*ALPHA*pow(z,ALPHA); + dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*z1*px/PAR[PM_PAR_SXX]/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*z1*py/PAR[PM_PAR_SYY]/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -1.0*z1; + } + } + return (f); } @@ -302,14 +253,8 @@ } - // the normalization is modified by the slope - float index = 0.5 / ALPHA; - float bn = 1.9992*index - 0.3271; - float Io = exp(0.5*bn); - // set the model normalization if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) { return false; } - PAR[PM_PAR_I0] /= Io; // set the model position @@ -328,17 +273,9 @@ psEllipseAxes axes; pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); - float AspectRatio = axes.minor / axes.major; - - float index = 4.0; - float bn = 1.9992*index - 0.3271; - - // the integral of a Sersic has an analytical form as follows: - float logGamma = lgamma(2.0*index); - float bnFactor = pow(bn, 2.0*index); - float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor; - - psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio; - - return(Flux); + + float norm = 0.00168012; + float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm; + + return(flux); } @@ -359,7 +296,10 @@ pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); - // f = Io exp(-z^n) -> z^n = ln(Io/f) - psF64 zn = log(PAR[PM_PAR_I0] / flux); - psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / ALPHA); + // static value for DEV: + float kappa = 7.670628; + + // f = Io exp(-kappa*z^n) -> z^n = ln(Io/f) / kappa + psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa; + psF64 radius = axes.major * pow(zn, 0.5 / ALPHA); psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", Index: trunk/psModules/src/objects/models/pmModel_EXP.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_EXP.c (revision 35925) +++ trunk/psModules/src/objects/models/pmModel_EXP.c (revision 36085) @@ -45,4 +45,5 @@ #include "pmPSFtry.h" #include "pmDetections.h" +#include "pmModel_CentralPixel.h" #include "pmModel_EXP.h" @@ -65,5 +66,5 @@ // Lax parameter limits static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.05, 0.05, -1.0 }; -static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0 }; // Moderate parameter limits @@ -78,9 +79,14 @@ static float *paramsMinUse = paramsMinLax; static float *paramsMaxUse = paramsMaxLax; -static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5}; +static float betaUse[] = { 2, 3e6, 5, 5, 10.0, 10.0, 0.5}; static bool limitsApply = true; // Apply limits? -# include "pmModel_SERSIC.CP.h" +// # include "pmModel_SERSIC.CP.h" + +// the problems I'm having with the SERSIC-like functions are: +// 1) making sure I have the right functional form so that PAR[SXX,etc] represent R_eff (half-light radius) +// 2) getting the central pixel right +// 3) getting the derivaties right. psF32 PM_MODEL_FUNC (psVector *deriv, @@ -101,85 +107,26 @@ psAssert (z >= 0, "do not allow negative z values in model"); - float index = 1.0; - float par7 = 0.5; - float bn = 1.9992*index - 0.3271; - float Io = exp(bn); - - psF32 f2 = bn*sqrt(z); - psF32 f1 = Io*exp(-f2); - + // for EXP, we can hard-wire kappa(1): + // float index = 1.0; + float kappa = 1.70056; + + // sqrt(z) is r + float q = kappa*sqrt(z); + psF32 f0 = exp(-q); + + assert (isfinite(q)); + + // only worry about the central 4 pixels at most psF32 radius = hypot(X, Y); - if (radius < 1.0) { - - // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center - - // first, use Rmajor and index to find the central pixel flux (fraction of total flux) - psEllipseAxes axes; - pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); - - // get the central pixel flux from the lookup table - float xPix = (axes.major - centralPixelXo) / centralPixeldX; - xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1); - float yPix = (index - centralPixelYo) / centralPixeldY; - yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1); - - // the integral of a Sersic has an analytical form as follows: - float logGamma = lgamma(2.0*index); - float bnFactor = pow(bn, 2.0*index); - float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor; - - // XXX interpolate to get the value - // XXX for the moment, just integerize - // XXX I need to multiply by the integrated flux to get the flux in the central pixel - float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm; - - float px1 = 1.0 / PAR[PM_PAR_SXX]; - float py1 = 1.0 / PAR[PM_PAR_SYY]; - float z10 = PS_SQR(px1); - float z01 = PS_SQR(py1); - - // which pixels do we need for this interpolation? - // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...) - if ((X >= 0) && (Y >= 0)) { - float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive - float V00 = Vcenter; - float V10 = Io*exp(-bn*pow(z10,par7)); - float V01 = Io*exp(-bn*pow(z01,par7)); - float V11 = Io*exp(-bn*pow(z11,par7)); - f1 = interpolatePixels(V00, V10, V01, V11, X, Y); - } - if ((X < 0) && (Y >= 0)) { - float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative - float V00 = Io*exp(-bn*pow(z10,par7)); - float V10 = Vcenter; - float V01 = Io*exp(-bn*pow(z11,par7)); - float V11 = Io*exp(-bn*pow(z01,par7)); - f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y); - } - if ((X >= 0) && (Y < 0)) { - float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative - float V00 = Io*exp(-bn*pow(z01,par7)); - float V10 = Io*exp(-bn*pow(z11,par7)); - float V01 = Vcenter; - float V11 = Io*exp(-bn*pow(z10,par7)); - f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y)); - } - if ((X < 0) && (Y < 0)) { - float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive - float V00 = Io*exp(-bn*pow(z11,par7)); - float V10 = Io*exp(-bn*pow(z10,par7)); - float V01 = Io*exp(-bn*pow(z01,par7)); - float V11 = Vcenter; - f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y)); - } - } - - psF32 z0 = PAR[PM_PAR_I0]*f1; - psF32 f0 = PAR[PM_PAR_SKY] + z0; - - assert (isfinite(f2)); + if (radius <= 1.5) { + f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], 1.0, 51); + } + assert (isfinite(f0)); + + psF32 f1 = PAR[PM_PAR_I0]*f0; + psF32 f = PAR[PM_PAR_SKY] + f1; + assert (isfinite(f1)); - assert (isfinite(z0)); - assert (isfinite(f0)); + assert (isfinite(f)); if (deriv != NULL) { @@ -187,18 +134,24 @@ dPAR[PM_PAR_SKY] = +1.0; - dPAR[PM_PAR_I0] = +f1; - - // gradient is infinite for z = 0; saturate at z = 0.01 - // z1 is -df/dz (the negative sign is canceled by most of dz/dPAR[i] - psF32 z1 = (z < 0.01) ? 0.5*bn*z0/sqrt(0.01) : 0.5*bn*z0/sqrt(z); - - // XXX dampen SXX and SYY as in GAUSS? - dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]); - dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]); - dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; - dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; - dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; - } - return (f0); + dPAR[PM_PAR_I0] = +f0; + + if (z > 0.01) { + float z1 = 0.5*f1*kappa/sqrt(z); + dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px + Y*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py + X*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; + } else { + // gradient -> 0 for z -> 0, but has undef form + float z1 = 0.5*f1*kappa; + dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*z1*px/PAR[PM_PAR_SXX]/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*z1*py/PAR[PM_PAR_SYY]/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -1.0*z1; + } + } + return (f); } @@ -314,17 +267,11 @@ psEllipseAxes axes; pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); - float AspectRatio = axes.minor / axes.major; - - float index = 1.0; - float bn = 1.9992*index - 0.3271; - - // the integral of a Sersic has an analytical form as follows: - float logGamma = lgamma(2.0*index); - float bnFactor = pow(bn, 2.0*index); - float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor; - - psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio; - - return(Flux); + + // static value for EXP: + float norm = 0.34578; // \int exp(-kappa*sqrt(z)) r dr + + float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm; + + return(flux); } @@ -345,7 +292,10 @@ pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); - // f = Io exp(-sqrt(z)) -> sqrt(z) = ln(Io/f) - psF64 zn = log(PAR[PM_PAR_I0] / flux); - psF64 radius = axes.major * sqrt (2.0) * zn; + // static value for EXP: + float kappa = 1.70056; + + // f = Io exp(-kappa*sqrt(z)) -> sqrt(z) = ln(Io/f) / kappa + psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa; + psF64 radius = axes.major * zn; psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", @@ -501,2 +451,66 @@ return; } + +# if (0) +void bilin_inter_function () { + // first, use Rmajor and index to find the central pixel flux (fraction of total flux) + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); + + // get the central pixel flux from the lookup table + float xPix = (axes.major - centralPixelXo) / centralPixeldX; + xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1); + float yPix = (index - centralPixelYo) / centralPixeldY; + yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1); + + // the integral of a Sersic has an analytical form as follows: + float logGamma = lgamma(2.0*index); + float bnFactor = pow(bn, 2.0*index); + float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor; + + // XXX interpolate to get the value + // XXX for the moment, just integerize + // XXX I need to multiply by the integrated flux to get the flux in the central pixel + float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm; + + float px1 = 1.0 / PAR[PM_PAR_SXX]; + float py1 = 1.0 / PAR[PM_PAR_SYY]; + float z10 = PS_SQR(px1); + float z01 = PS_SQR(py1); + + // which pixels do we need for this interpolation? + // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...) + if ((X >= 0) && (Y >= 0)) { + float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive + float V00 = Vcenter; + float V10 = Io*exp(-bn*pow(z10,par7)); + float V01 = Io*exp(-bn*pow(z01,par7)); + float V11 = Io*exp(-bn*pow(z11,par7)); + f1 = interpolatePixels(V00, V10, V01, V11, X, Y); + } + if ((X < 0) && (Y >= 0)) { + float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative + float V00 = Io*exp(-bn*pow(z10,par7)); + float V10 = Vcenter; + float V01 = Io*exp(-bn*pow(z11,par7)); + float V11 = Io*exp(-bn*pow(z01,par7)); + f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y); + } + if ((X >= 0) && (Y < 0)) { + float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative + float V00 = Io*exp(-bn*pow(z01,par7)); + float V10 = Io*exp(-bn*pow(z11,par7)); + float V01 = Vcenter; + float V11 = Io*exp(-bn*pow(z10,par7)); + f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y)); + } + if ((X < 0) && (Y < 0)) { + float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive + float V00 = Io*exp(-bn*pow(z11,par7)); + float V10 = Io*exp(-bn*pow(z10,par7)); + float V01 = Io*exp(-bn*pow(z01,par7)); + float V11 = Vcenter; + f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y)); + } +} +# endif Index: trunk/psModules/src/objects/models/pmModel_GAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 35925) +++ trunk/psModules/src/objects/models/pmModel_GAUSS.c (revision 36085) @@ -61,5 +61,5 @@ // Lax parameter limits static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0 }; -static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0 }; // Moderate parameter limits Index: trunk/psModules/src/objects/models/pmModel_PGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 35925) +++ trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 36085) @@ -61,5 +61,5 @@ // Lax parameter limits static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0 }; -static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0 }; // Moderate parameter limits Index: trunk/psModules/src/objects/models/pmModel_PS1_V1.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_PS1_V1.c (revision 35925) +++ trunk/psModules/src/objects/models/pmModel_PS1_V1.c (revision 36085) @@ -70,15 +70,15 @@ // Lax parameter limits static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -1.0 }; -static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 }; // Moderate parameter limits // Tolerate a small divot (k < 0) static float paramsMinModerate[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -0.05 }; -static float paramsMaxModerate[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 }; +static float paramsMaxModerate[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 }; // Strict parameter limits // k = PAR_7 < 0 is very undesirable (big divot in the middle) static float paramsMinStrict[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.0 }; -static float paramsMaxStrict[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 }; +static float paramsMaxStrict[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 }; // Parameter limits to use Index: trunk/psModules/src/objects/models/pmModel_QGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 35925) +++ trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 36085) @@ -70,15 +70,15 @@ // Lax parameter limits static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -1.0 }; -static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 }; // Moderate parameter limits // Tolerate a small divot (k < 0) static float paramsMinModerate[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -0.05 }; -static float paramsMaxModerate[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 }; +static float paramsMaxModerate[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 }; // Strict parameter limits // k = PAR_7 < 0 is very undesirable (big divot in the middle) static float paramsMinStrict[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.0 }; -static float paramsMaxStrict[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 }; +static float paramsMaxStrict[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 }; // Parameter limits to use Index: trunk/psModules/src/objects/models/pmModel_RGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 35925) +++ trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 36085) @@ -66,5 +66,5 @@ // Lax parameter limits static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 1.25 }; -static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 4.0 }; // Moderate parameter limits Index: trunk/psModules/src/objects/models/pmModel_SERSIC.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 35925) +++ trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 36085) @@ -20,5 +20,5 @@ * note that a Sersic model is usually defined in terms of R_e, the half-light radius. This construction does not include a factor of 2 in the X^2 term, etc, like for a Gaussian. - Conversion from SXX, SYY, SXY to R_major, R_minor, theta can be done by using setting: + Conversion from SXX, SYY, SXY to R_major, R_minor, theta can be done by using: shape.sx = SXX / sqrt(2), shape.sy = SYY / sqrt(2), shape.sxy = SXY, then calling psEllipseShapeToAxes, and multiplying the values of axes.major, axes.minor by sqrt(2) @@ -55,4 +55,5 @@ #include "pmPSFtry.h" #include "pmDetections.h" +#include "pmModel_CentralPixel.h" #include "pmModel_SERSIC.h" @@ -74,6 +75,6 @@ // Lax parameter limits -static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0, 0.05 }; -static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 }; +static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0, 0.1 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 1.0 }; // Moderate parameter limits @@ -88,9 +89,9 @@ static float *paramsMinUse = paramsMinLax; static float *paramsMaxUse = paramsMaxLax; -static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5, 2.0 }; +static float betaUse[] = { 2, 3e6, 5, 5, 10.0, 10.0, 0.5, 1.0}; static bool limitsApply = true; // Apply limits? -# include "pmModel_SERSIC.CP.h" +// # include "pmModel_SERSIC.CP.h" psF32 PM_MODEL_FUNC (psVector *deriv, @@ -111,91 +112,33 @@ psAssert (z >= 0, "do not allow negative z values in model"); - float index = 0.5 / PAR[PM_PAR_7]; - float par7 = PAR[PM_PAR_7]; - float bn = 1.9992*index - 0.3271; - float Io = exp(bn); - - psF32 f2 = bn*pow(z,par7); - psF32 f1 = Io*exp(-f2); + float Sindex = 0.5 / PAR[PM_PAR_7]; + float kappa = pmSersicKappa (Sindex); + + float q = kappa*pow(z,PAR[PM_PAR_7]); + psF32 f0 = exp(-q); + + assert (isfinite(q)); psF32 radius = hypot(X, Y); - if (radius < 1.0) { - - // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center - - // first, use Rmajor and index to find the central pixel flux (fraction of total flux) - psEllipseAxes axes; - pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); - - // get the central pixel flux from the lookup table - float xPix = (axes.major - centralPixelXo) / centralPixeldX; - xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1); - float yPix = (index - centralPixelYo) / centralPixeldY; - yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1); - - // the integral of a Sersic has an analytical form as follows: - float logGamma = lgamma(2.0*index); - float bnFactor = pow(bn, 2.0*index); - float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor; - - // XXX interpolate to get the value - // XXX for the moment, just integerize - // XXX I need to multiply by the integrated flux to get the flux in the central pixel - float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm; - - float px1 = 1.0 / PAR[PM_PAR_SXX]; - float py1 = 1.0 / PAR[PM_PAR_SYY]; - float z10 = PS_SQR(px1); - float z01 = PS_SQR(py1); - - // which pixels do we need for this interpolation? - // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...) - if ((X >= 0) && (Y >= 0)) { - float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive - float V00 = Vcenter; - float V10 = Io*exp(-bn*pow(z10,par7)); - float V01 = Io*exp(-bn*pow(z01,par7)); - float V11 = Io*exp(-bn*pow(z11,par7)); - f1 = interpolatePixels(V00, V10, V01, V11, X, Y); - } - if ((X < 0) && (Y >= 0)) { - float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative - float V00 = Io*exp(-bn*pow(z10,par7)); - float V10 = Vcenter; - float V01 = Io*exp(-bn*pow(z11,par7)); - float V11 = Io*exp(-bn*pow(z01,par7)); - f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y); - } - if ((X >= 0) && (Y < 0)) { - float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative - float V00 = Io*exp(-bn*pow(z01,par7)); - float V10 = Io*exp(-bn*pow(z11,par7)); - float V01 = Vcenter; - float V11 = Io*exp(-bn*pow(z10,par7)); - f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y)); - } - if ((X < 0) && (Y < 0)) { - float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive - float V00 = Io*exp(-bn*pow(z11,par7)); - float V10 = Io*exp(-bn*pow(z10,par7)); - float V01 = Io*exp(-bn*pow(z01,par7)); - float V11 = Vcenter; - f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y)); - } - } - - psF32 z0 = PAR[PM_PAR_I0]*f1; - psF32 f0 = PAR[PM_PAR_SKY] + z0; - - if (!isfinite(z0)) { - fprintf(stderr, "z0 is not finite for %f %f %f %f %f. Parameters: \n", X, Y, radius, z, f1); + if (radius <= 1.5) { + // Nsub ~ 10*index^2 + 1 + psEllipseAxes axes = pmPSF_ModelToAxes(PAR, pmModelClassGetType ("PS_MODEL_SERSIC")); + int Nsub = 2 * ((int)(6.0*Sindex / axes.minor)) + 1; + Nsub = PS_MIN (Nsub, 121); + Nsub = PS_MAX (Nsub, 11); + f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], Sindex, Nsub); + } + if (!isfinite(f0)) { + fprintf(stderr, "f0 is not finite for %f %f %f %f %f. Parameters: \n", X, Y, radius, z, q); fprintf(stderr, "%f %f %f %f %f %f %f %f\n", PAR[0], PAR[1], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[7]); } - - assert (isfinite(f2)); + assert (isfinite(f0)); + + psF32 f1 = PAR[PM_PAR_I0]*f0; + psF32 f = PAR[PM_PAR_SKY] + f1; + assert (isfinite(f1)); - assert (isfinite(z0)); - assert (isfinite(f0)); + assert (isfinite(f)); if (deriv != NULL) { @@ -203,22 +146,28 @@ dPAR[PM_PAR_SKY] = +1.0; - dPAR[PM_PAR_I0] = +f1; - - // gradient is infinite for z = 0; saturate at z = 0.01 - psF32 z1 = (z < 0.01) ? z0*bn*par7*pow(0.01,par7 - 1.0) : z0*bn*par7*pow(z,par7 - 1.0); - - dPAR[PM_PAR_7] = (z < 0.01) ? -z0*pow(0.01,par7)*log(0.01) : -z0*f2*log(z); - dPAR[PM_PAR_7] *= 3.0; - - assert (isfinite(z1)); + dPAR[PM_PAR_I0] = +f0; + + if (z > 0.01) { + float z1 = f1*kappa*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7]-1.0); + dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px + Y*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py + X*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; + dPAR[PM_PAR_7] = -1.0*f1*q*log(z); + } else { + // gradient -> 0 for z -> 0, but has undef form + float z1 = f1*kappa*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7]); + dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]); + dPAR[PM_PAR_SXX] = +2.0*z1*px/PAR[PM_PAR_SXX]/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +2.0*z1*py/PAR[PM_PAR_SYY]/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -1.0*z1; + // dPAR[PM_PAR_7] = -1.0*f1*q*log(z + 0.0001); + dPAR[PM_PAR_7] = -1.0*f1*q*log(z + 0.0001); // factor of 16 to reduce the gain + } assert (isfinite(dPAR[PM_PAR_7])); - - dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]); - dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]); - dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; // XXX : increase drag? - dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; - dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; - } - return (f0); + } + return (f); } @@ -370,17 +319,11 @@ psEllipseAxes axes; pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); - float AspectRatio = axes.minor / axes.major; - - float index = 0.5 / PAR[PM_PAR_7]; - float bn = 1.9992*index - 0.3271; - - // the integral of a Sersic has an analytical form as follows: - float logGamma = lgamma(2.0*index); - float bnFactor = pow(bn, 2.0*index); - float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor; - - psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio; - - return(Flux); + + float Sindex = 0.5 / PAR[PM_PAR_7]; + float norm = pmSersicNorm (Sindex); + + float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm; + + return(flux); } @@ -401,7 +344,10 @@ pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); + float Sindex = 0.5 / PAR[PM_PAR_7]; + float kappa = pmSersicKappa (Sindex); + // f = Io exp(-z^n) -> z^n = ln(Io/f) - psF64 zn = log(PAR[PM_PAR_I0] / flux); - psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / PAR[PM_PAR_7]); + psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa; + psF64 radius = axes.major * pow(zn, Sindex); psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f), par 7 = %f", Index: trunk/psModules/src/objects/models/pmModel_TRAIL.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_TRAIL.c (revision 35925) +++ trunk/psModules/src/objects/models/pmModel_TRAIL.c (revision 36085) @@ -61,5 +61,5 @@ // Lax parameter limits static float paramsMinLax[] = { -1.0e3, 1.0e-2, -1.0e2, -1.0e2, 0.5, -3.3, -0.5 }; -static float paramsMaxLax[] = { 1.0e5, 1.0e+8, +1.0e4, +1.0e4, 150.0, +3.3 , 5.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.00+9, +1.0e5, +1.0e5, 150.0, +3.3 , 5.0 }; // Moderate parameter limits Index: trunk/psModules/src/objects/pmFootprintCullPeaks.c =================================================================== --- trunk/psModules/src/objects/pmFootprintCullPeaks.c (revision 35925) +++ trunk/psModules/src/objects/pmFootprintCullPeaks.c (revision 36085) @@ -25,10 +25,10 @@ bool dumpfootprints (pmFootprint *fp, pmFootprintSpans *fpSp); - /* - * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently - * isolated. More precisely, for each peak find the highest coll that you'd have to traverse - * to reach a still higher peak --- and if that coll's more (less?) than nsigma DN below your - * starting point, discard the peak. - */ +/* + * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently + * isolated. More precisely, for each peak find the highest coll that you'd have to traverse + * to reach a still higher peak --- and if that coll's more (less?) than nsigma DN below your + * starting point, discard the peak. + */ # define IN_PEAK 1 @@ -48,5 +48,5 @@ if (fp->peaks == NULL || fp->peaks->n < 2) { // nothing to do - return PS_ERR_NONE; + return PS_ERR_NONE; } @@ -91,16 +91,16 @@ // max flux is above threshold for brightest peak - pmPeak *maxPeak = NULL; - for (int i = 0; i < fp->peaks->n; i++) { - pmPeak *testPeak = fp->peaks->data[i]; - float this_peak = useSmoothedImage ? testPeak->smoothFlux : testPeak->rawFlux; + pmPeak *maxPeak = NULL; + for (int i = 0; i < fp->peaks->n; i++) { + pmPeak *testPeak = fp->peaks->data[i]; + float this_peak = useSmoothedImage ? testPeak->smoothFlux : testPeak->rawFlux; - if (isfinite(this_peak)) { - maxPeak = fp->peaks->data[i]; - break; - } - } - psAssert(maxPeak,"maxPeak was not set in these peaks"); - // = fp->peaks->data[0]; + if (isfinite(this_peak)) { + maxPeak = fp->peaks->data[i]; + break; + } + } + psAssert(maxPeak,"maxPeak was not set in these peaks"); + // = fp->peaks->data[0]; float maxFlux = useSmoothedImage ? maxPeak->smoothFlux : maxPeak->rawFlux; @@ -130,8 +130,8 @@ } #if (0) - if (threshbounds->data.F32[threshbounds->n-1] > maxFlux) { - psWarning ("upper limit: %f does not include max flux: %f", - threshbounds->data.F32[threshbounds->n-1], maxFlux); - } + if (threshbounds->data.F32[threshbounds->n-1] > maxFlux) { + psWarning ("upper limit: %f does not include max flux: %f", + threshbounds->data.F32[threshbounds->n-1], maxFlux); + } #endif psHistogram *threshist = psHistogramAllocGeneric(threshbounds); Index: trunk/psModules/src/objects/pmModel_CentralPixel.c =================================================================== --- trunk/psModules/src/objects/pmModel_CentralPixel.c (revision 36085) +++ trunk/psModules/src/objects/pmModel_CentralPixel.c (revision 36085) @@ -0,0 +1,813 @@ +/* @file pmModel_CentralPixel.c + * @brief Functions to manage the central pixel for sersic-like models + * @author EAM, IfA + * + * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-16 22:30:50 $ + * + * Copyright 2013 Institute for Astronomy, University of Hawaii + */ + +/****************************************************************************** + * this file contains functions to determine the flux of the central pixel(s) for an + * exponential / devaucouleur / sersic style galaxy model. The problem is that (a) these + * models are so centrally-peaked that it is necessary to determine the true mean flux in the + * central pixel by integration of fractional pixels (at least 0.02 pixels in the case of a + * DEV model) and (b) the process of integrating the central pixel is too slow to be used for + * real processing. + + * we bypass this problem by defining a set of pre-calculated central pixel images, with + * subpixel resolution > 1 pixel (maybe 11 subpixels per real pixel). These pre-calculated + * images are generated for a series of values for the following parameters: sersic index, + * effective radius, axial ratio. We then select the closest image to our specific case, and + * integrate over the true sub-pixels relevant for our position and model. We have thus turned + * problem from 2500 evaluations of the full sersic model to ~100 straight additions (possibly + * x 6 if we need to interpolate in each of the dimensions). + + * we need a number of support functions: + + * pmModelCP_Load : load CP model data from the specified file. + * pmModelCP_GetImage : choose an appropriate CP model image for a given set of parameters + * pmModelCP_GetValue : calculate the true CP value for the given image and parameters + + *****************************************************************************/ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include + +#include "pmModel_CentralPixel.h" + +static void pmModelCP_Free(pmModelCP *cp) { + psFree (cp->flux); +} + +pmModelCP *pmModelCP_Alloc(void) +{ + pmModelCP *tmp = (pmModelCP *) psAlloc(sizeof(pmModelCP)); + psMemSetDeallocator(tmp, (psFreeFunc) pmModelCP_Free); + + tmp->flux = NULL; + tmp->Rmajor = NAN; + tmp->Aratio = NAN; + tmp->Sindex = NAN; + + return tmp; +} + +static void pmModelCPset_Free(pmModelCPset *tmp) { + + for (int i = 0; i < tmp->RmajorNitem; i++) { + for (int j = 0; j < tmp->AratioNitem; j++) { + psFree (tmp->lookupCube[i][j]); + } + psFree (tmp->lookupCube[i]); + } + psFree (tmp->lookupCube); + psFree (tmp->images); +} + +pmModelCPset *pmModelCPset_Alloc(void) +{ + pmModelCPset *tmp = (pmModelCPset *) psAlloc(sizeof(pmModelCPset)); + psMemSetDeallocator(tmp, (psFreeFunc) pmModelCPset_Free); + + tmp->RmajorMin = NAN; + tmp->RmajorMax = NAN; + tmp->RmajorDel = NAN; + + tmp->AratioMin = NAN; + tmp->AratioMax = NAN; + tmp->AratioDel = NAN; + + tmp->SindexMin = NAN; + tmp->SindexMax = NAN; + tmp->SindexDel = NAN; + + tmp->RmajorNitem = 0; + tmp->AratioNitem = 0; + tmp->SindexNitem = 0; + + tmp->lookupCube = NULL; + tmp->images = NULL; + + return tmp; +} + +// load the central-pixel maps as an array of pmModelCP +pmModelCPset *pmModelCP_Load (char *filename) { + + bool status; + + psFits *fits = psFitsOpen (filename, "r"); + if (!fits) { + return false; + } + + // read the PHU -- it defines descriptive metadata + psMetadata *PHU = psFitsReadHeader (NULL, fits); + if (!PHU) { + psFitsClose (fits); + return false; + } + + pmModelCPset *CPset = pmModelCPset_Alloc(); + + // NOTE : RMAJOR refers to the LOG_10 of the major axis + CPset->RmajorMin = psMetadataLookupF32 (&status, PHU, "RMAJ_MIN"); psAssert (status, "missing keyword RMAJ_MIN"); + CPset->RmajorMax = psMetadataLookupF32 (&status, PHU, "RMAJ_MAX"); psAssert (status, "missing keyword RMAJ_MAX"); + CPset->RmajorDel = psMetadataLookupF32 (&status, PHU, "RMAJ_DEL"); psAssert (status, "missing keyword RMAJ_DEL"); + CPset->AratioMin = psMetadataLookupF32 (&status, PHU, "ARAT_MIN"); psAssert (status, "missing keyword ARAT_MIN"); + CPset->AratioMax = psMetadataLookupF32 (&status, PHU, "ARAT_MAX"); psAssert (status, "missing keyword ARAT_MAX"); + CPset->AratioDel = psMetadataLookupF32 (&status, PHU, "ARAT_DEL"); psAssert (status, "missing keyword ARAT_DEL"); + CPset->SindexMin = psMetadataLookupF32 (&status, PHU, "SIDX_MIN"); psAssert (status, "missing keyword SIDX_MIN"); + CPset->SindexMax = psMetadataLookupF32 (&status, PHU, "SIDX_MAX"); psAssert (status, "missing keyword SIDX_MAX"); + CPset->SindexDel = psMetadataLookupF32 (&status, PHU, "SIDX_DEL"); psAssert (status, "missing keyword SIDX_DEL"); + + CPset->RmajorNitem = 1 + (int)(0.5 + (CPset->RmajorMax - CPset->RmajorMin) / CPset->RmajorDel); + CPset->AratioNitem = 1 + (int)(0.5 + (CPset->AratioMax - CPset->AratioMin) / CPset->AratioDel); + CPset->SindexNitem = 1 + (int)(0.5 + (CPset->SindexMax - CPset->SindexMin) / CPset->SindexDel); + + // array entry = lookupCube[RmajorBin][AratioBin][SindexBin] + + CPset->lookupCube = (int ***) psAlloc (sizeof(int **)*CPset->RmajorNitem); + for (int i = 0; i < CPset->RmajorNitem; i++) { + CPset->lookupCube[i] = (int **) psAlloc (sizeof(int *)*CPset->AratioNitem); + for (int j = 0; j < CPset->AratioNitem; j++) { + CPset->lookupCube[i][j] = (int *) psAlloc (sizeof(int)*CPset->SindexNitem); + for (int k = 0; k < CPset->SindexNitem; k++) { + CPset->lookupCube[i][j][k] = -1; + } + } + } + + CPset->images = psArrayAllocEmpty (CPset->RmajorNitem*CPset->AratioNitem*CPset->SindexNitem); + + // the CP file contains a set of 2D images; load them all + + psRegion fullImage = psRegionSet (0, 0, 0, 0); + while (true) { + bool status = psFitsMoveExtNum (fits, 1, true); + if (!status) break; + + psMetadata *header = psFitsReadHeader (NULL, fits); + if (!header) { + fprintf (stderr, "error reading header\n"); + return false; + } + + pmModelCP *cp = pmModelCP_Alloc (); + + cp->flux = psFitsReadImageBuffer (NULL, fits, fullImage, 0); + if (!cp->flux) { + fprintf (stderr, "error reading image\n"); + return false; + } + + cp->Rmajor = psMetadataLookupF32 (&status, header, "R_MAJOR"); + cp->Aratio = psMetadataLookupF32 (&status, header, "A_RATIO"); + cp->Sindex = psMetadataLookupF32 (&status, header, "S_INDEX"); + + int RmajorBin = (int)((cp->Rmajor - CPset->RmajorMin) / CPset->RmajorDel); psAssert ((RmajorBin < CPset->RmajorNitem) && (RmajorBin >= 0), "bad bin"); + int AratioBin = (int)((cp->Aratio - CPset->AratioMin) / CPset->AratioDel); psAssert ((AratioBin < CPset->AratioNitem) && (AratioBin >= 0), "bad bin"); + int SindexBin = (int)((cp->Sindex - CPset->SindexMin) / CPset->SindexDel); psAssert ((SindexBin < CPset->SindexNitem) && (SindexBin >= 0), "bad bin"); + + CPset->lookupCube[RmajorBin][AratioBin][SindexBin] = CPset->images->n; + + psArrayAdd (CPset->images, 121, cp); + psFree (cp); + psFree (header); + } + + psFree (PHU); + psFitsClose (fits); + + return CPset; +} + +// choose the closest image to the given coords +pmModelCP *pmModelCP_GetImage (pmModelCPset *CPset, float Rmajor, float Aratio, float Sindex) { + + // the pmModelCP set is defined for a grid of Rmajor, Aratio, Sindex values + + int RmajorBin = (int)((Rmajor - CPset->RmajorMin) / CPset->RmajorDel); psAssert ((RmajorBin < CPset->RmajorNitem) && (RmajorBin >= 0), "bad bin"); + int AratioBin = (int)((Aratio - CPset->AratioMin) / CPset->AratioDel); psAssert ((AratioBin < CPset->AratioNitem) && (AratioBin >= 0), "bad bin"); + int SindexBin = (int)((Sindex - CPset->SindexMin) / CPset->SindexDel); psAssert ((SindexBin < CPset->SindexNitem) && (SindexBin >= 0), "bad bin"); + + int entry = CPset->lookupCube[RmajorBin][AratioBin][SindexBin]; + + pmModelCP *cp = CPset->images->data[entry]; + + return (cp); +} + +// XXX for test purposes only: +# define TEST_IMAGE 0 +# if (TEST_IMAGE) +static psImage *map = NULL; +# endif + +float pmModelCP_GetFlux_RotSquare (pmModelCP *cp, float dx, float dy, float theta); + +float pmModelCP_GetFlux (pmModelCP *cp, float dx, float dy, float theta) { + +# if (TEST_IMAGE) + map = psImageCopy (map, cp->flux, PS_TYPE_S32); + psImageInit (map, 0.0); +# endif + + // float flux = pmModelCP_GetFlux_Bresen (cp, dx, dy, theta); + // float flux = pmModelCP_GetFlux_Old (cp, dx, dy, theta); + float flux = pmModelCP_GetFlux_RotSquare (cp, dx, dy, theta); + + // RotSquare for theta = 0.0 & Bresen give the same answer + // if I count from x[0] <= ix < x[1] + +# if (TEST_IMAGE) + psFits *fits = psFitsOpen ("map.fits", "w"); + psFitsWriteImage (fits, NULL, map, 0, NULL); + psFitsClose (fits); + psFree (map); +# endif + + return flux; +} + +float pmModelCP_GetFlux_Old (pmModelCP *cp, float dx, float dy, float theta) { + + // the cp data is defined for the central 3x3 pixels. we allow dx,dy to have values of + // -1.0 <= dx,dy <= +1.0 + + // Xsub = (Xim * cos(theta) - Yim * sin(theta) + 1.5) * Nsub + // Ysub = (Yim * cos(theta) + Xim * sin(theta) + 1.5) * Nsub + + // integrate from (dx - 0.5 to dx + 0.5), (dy - 0.5 to dy + 0.5), + + // get the Xsub,Ysub values for the 4 corners, find the Xmin,Xmax, Ymin,Ymax in the + // subrastered image + + float cs = cos(theta*PS_RAD_DEG); + float sn = sin(theta*PS_RAD_DEG); + + float Nsub = 11.0; + int Xsub00 = ((dx - 0.5)*cs - (dy - 0.5)*sn + 1.5)*Nsub; + int Ysub00 = ((dx - 0.5)*sn + (dy - 0.5)*cs + 1.5)*Nsub; + int Xsub01 = ((dx - 0.5)*cs - (dy + 0.5)*sn + 1.5)*Nsub; + int Ysub01 = ((dx - 0.5)*sn + (dy + 0.5)*cs + 1.5)*Nsub; + int Xsub10 = ((dx + 0.5)*cs - (dy - 0.5)*sn + 1.5)*Nsub; + int Ysub10 = ((dx + 0.5)*sn + (dy - 0.5)*cs + 1.5)*Nsub; + int Xsub11 = ((dx + 0.5)*cs - (dy + 0.5)*sn + 1.5)*Nsub; + int Ysub11 = ((dx + 0.5)*sn + (dy + 0.5)*cs + 1.5)*Nsub; + + int Xmin, Xmax, Ymin, Ymax; + + Xmin = PS_MIN(Xsub00,Xsub01); + Xmin = PS_MIN(Xsub10,Xmin); + Xmin = PS_MIN(Xsub11,Xmin); + Xmin = PS_MIN(Xmin, cp->flux->numCols - 1); + Xmin = PS_MAX(Xmin, 0); + Xmax = PS_MAX(Xsub00,Xsub01); + Xmax = PS_MAX(Xsub10,Xmax); + Xmax = PS_MAX(Xsub11,Xmax); + Xmax = PS_MIN(Xmax, cp->flux->numCols - 1); + Xmax = PS_MAX(Xmax, 0); + Ymin = PS_MIN(Ysub00,Ysub01); + Ymin = PS_MIN(Ysub10,Ymin); + Ymin = PS_MIN(Ysub11,Ymin); + Ymin = PS_MIN(Ymin, cp->flux->numRows - 1); + Ymin = PS_MAX(Ymin, 0); + Ymax = PS_MAX(Ysub00,Ysub01); + Ymax = PS_MAX(Ysub10,Ymax); + Ymax = PS_MAX(Ysub11,Ymax); + Ymax = PS_MIN(Ymax, cp->flux->numRows - 1); + Ymax = PS_MAX(Ymax, 0); + + // integrate pixels from Xmin,Ymin to Xmax,Ymax, only include pixels contained in the + // target pixel + + float flux = 0.0; + int npix = 0; + for (int i = Xmin; i < Xmax; i++) { + float dX = i / Nsub - 1.5; + for (int j = Ymin; j < Ymax; j++) { + float dY = j / Nsub - 1.5; + + float Xim = dX*cs + dY*sn; + if (Xim < (dx - 0.5)) continue; + if (Xim > (dx + 0.5)) continue; + + float Yim = -dX*sn + dY*cs; + if (Yim < (dy - 0.5)) continue; + if (Yim > (dy + 0.5)) continue; + + flux += cp->flux->data.F32[j][i]; + npix ++; + } + } + + float normFlux = flux / npix; + return normFlux; +} + +// *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors +# define COMPARE_INDEX(A,B) (y[A] < y[B]) +# define SWAP_INDEX(TYPE,A,B) { \ + int tmp; \ + if (A != B) { \ + tmp = x[A]; \ + x[A] = x[B]; \ + x[B] = tmp; \ + tmp = y[A]; \ + y[A] = y[B]; \ + y[B] = tmp; \ + } \ + } + +bool pmModelCP_SortCorners (int *x, int *y, int Npar) { + + if (Npar < 2) return true; + + // sort the vector set by the radius + PSSORT (Npar, COMPARE_INDEX, SWAP_INDEX, NONE); + return true; +} + +float pmModelCP_GetFlux_RotSquare (pmModelCP *cp, float dx, float dy, float theta) { + + // the cp data is defined for the central 3x3 pixels. we allow dx,dy to have values of + // -1.0 <= dx,dy <= +1.0 + + // Xsub = (Xim * cos(theta) - Yim * sin(theta) + 1.5) * Nsub + // Ysub = (Yim * cos(theta) + Xim * sin(theta) + 1.5) * Nsub + + // integrate from (dx - 0.5 to dx + 0.5), (dy - 0.5 to dy + 0.5), + + // get the Xsub,Ysub values for the 4 corners, find the Xmin,Xmax, Ymin,Ymax in the + // subrastered image + + float cs = cos(theta*PS_RAD_DEG); + float sn = sin(theta*PS_RAD_DEG); + float Nsub = 11.0; + + int Xsub[4], Ysub[4]; + + Xsub[0] = ((dx - 0.5)*cs - (dy - 0.5)*sn + 1.5)*Nsub; + Ysub[0] = ((dx - 0.5)*sn + (dy - 0.5)*cs + 1.5)*Nsub; + Xsub[1] = ((dx - 0.5)*cs - (dy + 0.5)*sn + 1.5)*Nsub; + Ysub[1] = ((dx - 0.5)*sn + (dy + 0.5)*cs + 1.5)*Nsub; + Xsub[2] = ((dx + 0.5)*cs - (dy - 0.5)*sn + 1.5)*Nsub; + Ysub[2] = ((dx + 0.5)*sn + (dy - 0.5)*cs + 1.5)*Nsub; + Xsub[3] = ((dx + 0.5)*cs - (dy + 0.5)*sn + 1.5)*Nsub; + Ysub[3] = ((dx + 0.5)*sn + (dy + 0.5)*cs + 1.5)*Nsub; + + // first, sort the corners in order of the Y coordinate: + pmModelCP_SortCorners (Xsub, Ysub, 4); + + float flux = 0.0; + float npix = 0.0; + + // if (Ysub[0] == Ysub[1]), we have a simple square + if (Ysub[0] == Ysub[1]) { + psAssert (Ysub[2] == Ysub[3], "not square?"); + int Xmin = PS_MIN(Xsub[0], Xsub[1]); + int Xmax = PS_MAX(Xsub[0], Xsub[1]); + for (int iy = Ysub[0]; iy < Ysub[3]; iy++) { + for (int ix = Xmin; ix < Xmax; ix++) { + flux += cp->flux->data.F32[iy][ix]; + npix += 1.0; +# if (TEST_IMAGE) + fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]); + map->data.S32[iy][ix] ++; +# endif + } + } + float normFlux = flux / npix; + return normFlux; + } + + // second case: Xsub[1] > Xsub[2]: + if (Xsub[1] > Xsub[2]) { + float dYdXp, dYdXm; + // first segment, Ysub[0] to Ysub[1]: + dYdXp = (Ysub[1] - Ysub[0]) / (float) (Xsub[1] - Xsub[0]); + dYdXm = (Ysub[2] - Ysub[0]) / (float) (Xsub[2] - Xsub[0]); + for (int iy = Ysub[0]; iy < Ysub[1]; iy++) { + int Xs = (iy - Ysub[0]) / dYdXm + Xsub[0]; + int Xe = (iy - Ysub[0]) / dYdXp + Xsub[0]; + for (int ix = Xs; ix < Xe; ix ++) { + flux += cp->flux->data.F32[iy][ix]; + npix += 1.0; +# if (TEST_IMAGE) + fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]); + map->data.S32[iy][ix] ++; +# endif + } + } + // 2nd segment, Ysub[1] to Ysub[2]: + dYdXp = (Ysub[3] - Ysub[1]) / (float) (Xsub[3] - Xsub[1]); + dYdXm = (Ysub[2] - Ysub[0]) / (float) (Xsub[2] - Xsub[0]); + for (int iy = Ysub[1]; iy < Ysub[2]; iy++) { + int Xs = (iy - Ysub[0]) / dYdXm + Xsub[0]; + int Xe = (iy - Ysub[1]) / dYdXp + Xsub[1]; + for (int ix = Xs; ix < Xe; ix ++) { + flux += cp->flux->data.F32[iy][ix]; + npix += 1.0; +# if (TEST_IMAGE) + fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]); + map->data.S32[iy][ix] ++; +# endif + } + } + // first segment, Ysub[0] to Ysub[1]: + dYdXp = (Ysub[3] - Ysub[1]) / (float) (Xsub[3] - Xsub[1]); + dYdXm = (Ysub[3] - Ysub[2]) / (float) (Xsub[3] - Xsub[2]); + for (int iy = Ysub[2]; iy < Ysub[3]; iy++) { + int Xs = (iy - Ysub[2]) / dYdXm + Xsub[2]; + int Xe = (iy - Ysub[1]) / dYdXp + Xsub[1]; + for (int ix = Xs; ix < Xe; ix ++) { + flux += cp->flux->data.F32[iy][ix]; + npix += 1.0; +# if (TEST_IMAGE) + fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]); + map->data.S32[iy][ix] ++; +# endif + } + } + float normFlux = flux / npix; + return normFlux; + } + + // third case: Xsub[1] < Xsub[2]: + if (Xsub[2] > Xsub[1]) { + // first segment, Ysub[0] to Ysub[1]: + float dYdXp, dYdXm; + dYdXp = (Ysub[2] - Ysub[0]) / (float) (Xsub[2] - Xsub[0]); + dYdXm = (Ysub[1] - Ysub[0]) / (float) (Xsub[1] - Xsub[0]); + for (int iy = Ysub[0]; iy < Ysub[1]; iy++) { + int Xs = (iy - Ysub[0]) / dYdXm + Xsub[0]; + int Xe = (iy - Ysub[0]) / dYdXp + Xsub[0]; + for (int ix = Xs; ix < Xe; ix ++) { + flux += cp->flux->data.F32[iy][ix]; + npix += 1.0; +# if (TEST_IMAGE) + fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]); + map->data.S32[iy][ix] ++; +# endif + } + } + // 2nd segment, Ysub[1] to Ysub[2]: + dYdXp = (Ysub[2] - Ysub[0]) / (float) (Xsub[2] - Xsub[0]); + dYdXm = (Ysub[3] - Ysub[1]) / (float) (Xsub[3] - Xsub[1]); + for (int iy = Ysub[1]; iy < Ysub[2]; iy++) { + int Xs = (iy - Ysub[1]) / dYdXm + Xsub[1]; + int Xe = (iy - Ysub[0]) / dYdXp + Xsub[0]; + for (int ix = Xs; ix < Xe; ix ++) { + flux += cp->flux->data.F32[iy][ix]; + npix += 1.0; +# if (TEST_IMAGE) + fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]); + map->data.S32[iy][ix] ++; +# endif + } + } + // first segment, Ysub[0] to Ysub[1]: + dYdXp = (Ysub[3] - Ysub[2]) / (float) (Xsub[3] - Xsub[2]); + dYdXm = (Ysub[3] - Ysub[1]) / (float) (Xsub[3] - Xsub[1]); + for (int iy = Ysub[2]; iy < Ysub[3]; iy++) { + int Xs = (iy - Ysub[1]) / dYdXm + Xsub[1]; + int Xe = (iy - Ysub[2]) / dYdXp + Xsub[2]; + for (int ix = Xs; ix < Xe; ix ++) { + flux += cp->flux->data.F32[iy][ix]; + npix += 1.0; +# if (TEST_IMAGE) + fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]); + map->data.S32[iy][ix] ++; +# endif + } + } + float normFlux = flux / npix; + return normFlux; + } + myAbort ("impossible case?"); +} + +float pmModelCP_GetFlux_Bresen (pmModelCP *cp, float dx, float dy, float theta) { + + // the cp data is defined for the central 3x3 pixels. we allow dx,dy to have values of + // -1.0 <= dx,dy <= +1.0 + + // Xsub = ( Xim * cos(theta) + Yim * sin(theta) + 1.5) * Nsub + // Ysub = (-Yim * cos(theta) + Xim * sin(theta) + 1.5) * Nsub + + // integrate from (dx - 0.5 to dx + 0.5), (dy - 0.5 to dy + 0.5), + + // get the Xsub,Ysub values for the 4 corners, find the Xmin,Xmax, Ymin,Ymax in the + // subrastered image + + float cs = cos(theta*PS_RAD_DEG); + float sn = sin(theta*PS_RAD_DEG); + + float Nsub = 11.0; + int Xsub00 = 0.5 + ((dx - 0.5)*cs + (dy - 0.5)*sn + 1.5)*Nsub; + int Ysub00 = 0.5 + ((dy - 0.5)*cs - (dx - 0.5)*sn + 1.5)*Nsub; + int Xsub01 = 0.5 + ((dx - 0.5)*cs + (dy + 0.5)*sn + 1.5)*Nsub; + int Ysub01 = 0.5 + ((dy + 0.5)*cs - (dx - 0.5)*sn + 1.5)*Nsub; + int Xsub10 = 0.5 + ((dx + 0.5)*cs + (dy - 0.5)*sn + 1.5)*Nsub; + int Ysub10 = 0.5 + ((dy - 0.5)*cs - (dx + 0.5)*sn + 1.5)*Nsub; + int Xsub11 = 0.5 + ((dx + 0.5)*cs + (dy + 0.5)*sn + 1.5)*Nsub; + int Ysub11 = 0.5 + ((dy + 0.5)*cs - (dx + 0.5)*sn + 1.5)*Nsub; + + float flux = pmModelCP_GetFlux_BresenSquare (cp, Xsub00, Ysub00, Xsub10, Ysub10, Xsub01, Ysub01, Xsub11, Ysub11); + return flux; +} + +// first line is (X00,Y00) - (X10,Y1) : last line is (X01,Y01) - (X11,Y11) +float pmModelCP_GetFlux_BresenSquare (pmModelCP *cp, int X00, int Y00, int X10, int Y10, int X01, int Y01, int X11, int Y11) { + + int dX0 = X01 - X00; + int dY0 = Y01 - Y00; + + // int dX1 = X11 - X10; + // int dY1 = Y11 - Y10; + + // myAssert ((dX0 == dX1) && (dY0 == dY1), "pixel is not square?"); + + bool FlipCoords = (abs(dX0) < abs(dY0)); + bool FlipDirect = FlipCoords ? (Y00 > Y10) : (X00 > X10); + + float flux = 0.0; + if (!FlipDirect && !FlipCoords) flux = pmModelCP_GetFlux_BresenSquareBase (cp, X00, Y00, X10, Y10, X01, Y01, X11, Y11, false); + if ( FlipDirect && !FlipCoords) flux = pmModelCP_GetFlux_BresenSquareBase (cp, X10, Y10, X00, Y00, X11, Y11, X01, Y01, false); + if (!FlipDirect && FlipCoords) flux = pmModelCP_GetFlux_BresenSquareBase (cp, Y00, X00, Y10, X10, Y01, X01, Y11, X11, true); + if ( FlipDirect && FlipCoords) flux = pmModelCP_GetFlux_BresenSquareBase (cp, Y10, X10, Y00, X00, Y11, X11, Y01, X01, true); + return flux; +} + +// draw a line between (X00,Y00) & (X01,Y01) and increment to the next line segment until endpoints (X10,Y10) & (X11,Y11) +float pmModelCP_GetFlux_BresenSquareBase (pmModelCP *cp, int X00, int Y00, int X10, int Y10, int X01, int Y01, int X11, int Y11, bool swapcoords) { + + int dX0 = X01 - X00; + int dY0 = Y01 - Y00; + + // int dX1 = X11 - X10; + // int dY1 = Y11 - Y10; + + // myAssert ((dX0 == dX1) && (dY0 == dY1), "pixel is not square?"); + + float flux = 0.0; + float npix = 0.0; + + int Ys = Y00; + int Ye = Y10; + int e = 0; + for (int Xs = X00, Xe = X10; Xs < X01; Xs++, Xe++) { + if (swapcoords) { + pmModelCP_GetFlux_BresenLine (&flux, &npix, cp, Ys, Xs, Ye, Xe); + } else { + pmModelCP_GetFlux_BresenLine (&flux, &npix, cp, Xs, Ys, Xe, Ye); + } + e += dY0; + float e2 = 2 * e; + if (e2 > dX0) { + Ys++; + Ye++; + e -= dX0; + } + if (e2 < -dX0) { + Ys--; + Ye--; + e += dX0; + } + } + float normFlux = flux / npix; + // fprintf (stderr, "bres: %f %f %f\n", flux, (float) npix, normFlux); + return normFlux; +} + +// get the sequence right: +// if abs(dY) > abs(dX) : we will run in the Y direction not the X direction (we swap X and Y going in) +// if the direction (dX or dY) is negative, go the opposite direction +bool pmModelCP_GetFlux_BresenLine (float *flux, float *npix, pmModelCP *cp, int X0, int Y0, int X1, int Y1) { + + int dX = X1 - X0; + int dY = Y1 - Y0; + + bool FlipCoords = (abs(dX) < abs(dY)); + bool FlipDirect = FlipCoords ? (Y0 > Y1) : (X0 > X1); + + if (!FlipDirect && !FlipCoords) pmModelCP_GetFlux_BresenLineBase (flux, npix, cp, X0, Y0, X1, Y1, false); + if ( FlipDirect && !FlipCoords) pmModelCP_GetFlux_BresenLineBase (flux, npix, cp, X1, Y1, X0, Y0, false); + if (!FlipDirect && FlipCoords) pmModelCP_GetFlux_BresenLineBase (flux, npix, cp, Y0, X0, Y1, X1, true); + if ( FlipDirect && FlipCoords) pmModelCP_GetFlux_BresenLineBase (flux, npix, cp, Y1, X1, Y0, X0, true); + return true; +} + +bool pmModelCP_GetFlux_BresenLineBase (float *flux, float *npix, pmModelCP *cp, int X0, int Y0, int X1, int Y1, bool swapcoords) { + + int dX = X1 - X0; + int dY = Y1 - Y0; + + int Y = Y0; + int e = 0; + for (int X = X0; X < X1; X++) { + if (swapcoords) { + *flux += cp->flux->data.F32[X][Y]; + *npix += 1.0; +# if (TEST_IMAGE) + fprintf (stderr, "%d %d | %f %f | %f\n", X, Y, *flux, *npix, cp->flux->data.F32[X][Y]); + map->data.S32[X][Y] ++; +# endif + } else { + *flux += cp->flux->data.F32[Y][X]; + *npix += 1.0; +# if (TEST_IMAGE) + fprintf (stderr, "%d %d | %f %f | %f\n", X, Y, *flux, *npix, cp->flux->data.F32[Y][X]); + map->data.S32[Y][X] ++; +# endif + } + e += dY; + float e2 = 2 * e; + if (e2 > dX) { + Y++; + e -= dX; + } + if (e2 < -dX) { + Y--; + e += dX; + } + } + return true; +} + +// this is a test function which generates a full sersic model evaluation with sub-pixel sampling +float pmModelCP_FullSersic (float dx, float dy, float theta, float Rmajor, float Aratio, float Sindex) { + + float flux = 0.0; + int npix = 0; + + float Rminor = Aratio * Rmajor; + float f1 = 1.0 / PS_SQR(Rminor) + 1.0 / PS_SQR(Rmajor); + float f2 = 1.0 / PS_SQR(Rminor) - 1.0 / PS_SQR(Rmajor); + + float sxr = 0.5*f1 - 0.5*f2*cos(2.0*theta*PS_RAD_DEG); + float syr = 0.5*f1 + 0.5*f2*cos(2.0*theta*PS_RAD_DEG); + + float Rxx = +1.0 / sqrt(sxr); + float Ryy = +1.0 / sqrt(syr); + float Rxy = -f2*sin(2.0*theta*PS_RAD_DEG); + + float kappa = -0.275552 + 1.972625*Sindex + 0.003487 * PS_SQR(Sindex); + float rindex = 0.5 / Sindex; + + float off = -60.0/(11*11); + float delta = 1.0 / (11*11); + for (float ix = off; ix < 0.5; ix += delta) { + for (float iy = off; iy < 0.5; iy += delta) { + + float dX = dx + ix; + float dY = dy + iy; + float z = PS_SQR(dX / Rxx) + PS_SQR(dY / Ryy) + dX * dY * Rxy; + + float q = pow (z, rindex); + float f = exp(-kappa*q); + + flux += f; + npix ++; + } + } + float normFlux = flux / npix; + // fprintf (stderr, "full : %f %f %f\n", flux, (float) npix, normFlux); + return normFlux; +} + +// this is a test function which generates a full sersic model evaluation with sub-pixel sampling +// Nsub is the number of sub-pixel samplings and must be odd +// dx,dy are the centroid offset +float pmModelCP_SersicSubpix (float dx, float dy, float Rxx, float Rxy, float Ryy, float Sindex, int Nsub) { + + float flux = 0.0; + int npix = 0; + + float kappa = -0.275552 + 1.972625*Sindex + 0.003487 * PS_SQR(Sindex); + float rindex = 0.5 / Sindex; + + // have the resolution be a user-parameter? + psAssert (Nsub % 2 == 1, "Nsub is not odd"); + int Nsub2 = (Nsub - 1) / 2; + + float delta = 1.0 / (float) Nsub; + float off = -Nsub2 * delta; + for (float ix = off; ix < 0.5; ix += delta) { + for (float iy = off; iy < 0.5; iy += delta) { + + float dX = dx + ix; + float dY = dy + iy; + float z = PS_SQR(dX / Rxx) + PS_SQR(dY / Ryy) + dX * dY * Rxy; + + float q = pow (z, rindex); + float f = exp(-kappa*q); + + flux += f; + npix ++; + } + } + float normFlux = flux / npix; + // fprintf (stderr, "full : %f %f %f\n", flux, (float) npix, normFlux); + return normFlux; +} + +float pmSersicKappa (float Sindex) { + // this function is empirically derived from a fit to data for Sindex = 0.5 - 5.5 + // constrain Sindex or kappa? + float kappa = -0.275552 + 1.972625*Sindex + 0.003487 * PS_SQR(Sindex); + return kappa; +} + +float pmSersicNorm (float Sindex) { + + float C0 = NAN; + float C1 = NAN; + float C2 = NAN; + + // y = 0.201545 x^0 -0.950965 x^1 -0.315248 x^2 + // y = 0.402084 x^0 -1.357775 x^1 -0.105102 x^2 + // y = 0.619093 x^0 -1.591674 x^1 -0.041576 x^2 + // y = 0.770263 x^0 -1.696421 x^1 -0.023363 x^2 + // y = 0.885891 x^0 -1.755684 x^1 -0.015753 x^2 + + if ((Sindex >= 0.0) && (Sindex < 1.0)) { + C0 = 0.201545; C1 = -0.950965; C2 = -0.315248; + // y = 0.201545 x^0 -0.950965 x^1 -0.315248 x^2 + } + if ((Sindex >= 1.0) && (Sindex < 2.0)) { + C0 = 0.402084; C1 = -1.357775; C2 = -0.105102; + // y = 0.402084 x^0 -1.357775 x^1 -0.105102 x^2 + } + if ((Sindex >= 2.0) && (Sindex < 3.0)) { + C0 = 0.619093; C1 = -1.591674; C2 = -0.041576; + // y = 0.619093 x^0 -1.591674 x^1 -0.041576 x^2 + } + if ((Sindex >= 3.0) && (Sindex < 4.0)) { + C0 = 0.770263; C1 = -1.696421; C2 = -0.023363; + // y = 0.770263 x^0 -1.696421 x^1 -0.023363 x^2 + } + if ((Sindex >= 4.0) && (Sindex < 5.5)) { + C0 = 0.885891; C1 = -1.755684; C2 = -0.015753; + // y = 0.885891 x^0 -1.755684 x^1 -0.015753 x^2 + } + + if (isnan(C0)) return NAN; + + float lnorm = C0 + C1*Sindex + C2*Sindex*Sindex; + float norm = exp(lnorm); + return norm; +} + +# if (0) +// create a vector containing only the unique entries in the input vector +psVector *psVectorUniqueSubset (psVector *input) { + + // sort the input vector (to new temp vector) + // run through the sorted vector, copying to a new output vector if the current value is + // new + + psVector *temp = psVectorSort (input); + + psVector *output = psVectorAllocEmpty (0.5*input->n, PS_TYPE_F32); + + psVectorAppend (output, temp->data.F32[0]); + float lastValue = temp->data.F32[0]; + for (int i = 0; i < temp->n; i++) { + if (temp->data.F32[i] == lastValue) continue; + psVectorAppend (output, temp->data.F32[i]); + float lastValue = temp->data.F32[i]; + } + psFree (temp); + return output; +} + +getUnique() { + // we need to convert the collection of Rmajor, Aratio, Sindex values to a cube + // such that entry[RmajorBin][AratioBin][SindexBin] is the CPset array element + + // create full vectors will all Rmajor, Aratio, Sindex values: + psVector *RmajorAll = psVectorAllocEmpty (CPset->n, PS_TYPE_F32); + psVector *AratioAll = psVectorAllocEmpty (CPset->n, PS_TYPE_F32); + psVector *SindexAll = psVectorAllocEmpty (CPset->n, PS_TYPE_F32); + for (int i = 0; i < CPset->n; i++) { + pmModelCP *cp = CPset->data[i]; + psVectorAppend (RmajorAll, cp->Rmajor); + psVectorAppend (AratioAll, cp->Aratio); + psVectorAppend (SindexAll, cp->Sindex); + } + psVector *RmajorUniq = psVectorUniqueSubset (RmajorAll); + psVector *AratioUniq = psVectorUniqueSubset (AratioAll); + psVector *SindexUniq = psVectorUniqueSubset (SindexAll); +} +# endif Index: trunk/psModules/src/objects/pmModel_CentralPixel.h =================================================================== --- trunk/psModules/src/objects/pmModel_CentralPixel.h (revision 36085) +++ trunk/psModules/src/objects/pmModel_CentralPixel.h (revision 36085) @@ -0,0 +1,69 @@ +/* @file pmModel_CentralPixel.h + * @brief Functions to manage the central pixel for sersic-like models + * @author EAM, IfA + * + * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-16 22:30:50 $ + * + * Copyright 2013 Institute for Astronomy, University of Hawaii + */ + +# ifndef PM_MODEL_CENTRAL_PIXEL_H +# define PM_MODEL_CENTRAL_PIXEL_H + +/// @addtogroup Objects Object Detection / Analysis Functions +/// @{ + +typedef struct { + psImage *flux; + float Rmajor; + float Aratio; + float Sindex; +} pmModelCP; + +typedef struct { + psArray *images; + + float RmajorMin; + float RmajorMax; + float RmajorDel; + + float AratioMin; + float AratioMax; + float AratioDel; + + float SindexMin; + float SindexMax; + float SindexDel; + + int RmajorNitem; + int AratioNitem; + int SindexNitem; + + int ***lookupCube; + +} pmModelCPset; + +pmModelCP *pmModelCP_Alloc(void); +pmModelCPset *pmModelCPset_Alloc(void); + +pmModelCPset *pmModelCP_Load (char *filename); + +pmModelCP *pmModelCP_GetImage (pmModelCPset *CPset, float Rmajor, float Aratio, float Sindex); + +float pmModelCP_GetFlux (pmModelCP *cp, float dx, float dy, float theta); +float pmModelCP_FullSersic (float dx, float dy, float theta, float Rmajor, float Aratio, float Sindex); + +bool pmModelCP_GetFlux_BresenLineBase (float *flux, float *npix, pmModelCP *cp, int X0, int Y0, int X1, int Y1, bool swapcoords); +bool pmModelCP_GetFlux_BresenLine (float *flux, float *npix, pmModelCP *cp, int X0, int Y0, int X1, int Y1); +float pmModelCP_GetFlux_BresenSquareBase (pmModelCP *cp, int X00, int Y00, int X01, int Y01, int X10, int Y10, int X11, int Y11, bool swapcoords); +float pmModelCP_GetFlux_BresenSquare (pmModelCP *cp, int X00, int Y00, int X01, int Y01, int X10, int Y10, int X11, int Y11); +float pmModelCP_GetFlux_Bresen (pmModelCP *cp, float dx, float dy, float theta); +float pmModelCP_GetFlux_Old (pmModelCP *cp, float dx, float dy, float theta); + +float pmModelCP_SersicSubpix (float dx, float dy, float Rxx, float Rxy, float Ryy, float Sindex, int Nsub); + +float pmSersicKappa (float Sindex); +float pmSersicNorm (float Sindex); + +# endif Index: trunk/psModules/src/objects/pmPCM_MinimizeChisq.c =================================================================== --- trunk/psModules/src/objects/pmPCM_MinimizeChisq.c (revision 35925) +++ trunk/psModules/src/objects/pmPCM_MinimizeChisq.c (revision 36085) @@ -42,4 +42,9 @@ #include "pmPCMdata.h" +# define SAVE_IMAGES 0 +# if (SAVE_IMAGES) +int psphotSaveImage (psMetadata *header, psImage *image, char *filename); +# endif + # define FACILITY "psModules.objects" @@ -91,5 +96,5 @@ psF32 lambda = 0.001; psF32 dLinear = 0.0; - psF32 nu = 2.0; + psF32 nu = 3.0; # if (USE_FFT && PRE_CONVOLVE) @@ -130,6 +135,30 @@ } + char key[10]; // used for interactive responses + bool testValue = false; + // set a new guess for Alpha, Beta, Params if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda, &dLinear)) { + if (min->isInteractive) { + fprintf (stdout, "guess failed (singular matrix or NaN values), continue? [Y,n] "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + switch (key[0]) { + case 'n': + case 'N': + done = true; + break; + case 'y': + case 'Y': + case '\n': + lambda *= 10.0; + continue; + default: + lambda *= 10.0; + continue; + } + if (done) break; + } min->iter ++; if (min->iter >= min->maxIter) break; @@ -138,4 +167,40 @@ } + if (min->isInteractive) { + p_psVectorPrint(psTraceGetDestination(), Params, "current parameters: "); + fprintf (stdout, "last chisq : %f\n", min->value); + bool getOptions = true; + while (getOptions) { + fprintf (stdout, "options: (m)odify, (g)o, (q)uit: "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + switch (key[0]) { + case 'm': + case 'M': + testValue = TRUE; + fprintf (stdout, "enter (Npar) (value): "); + int Npar = 0; + float value= 0; + int Nscan = fscanf (stdin, "%d %f", &Npar, &value); + if (Nscan != 2) { + fprintf (stderr, "scan failure\n"); + } + Params->data.F32[Npar] = value; + break; + case 'g': + case 'G': + case '\n': + getOptions = false; + break; + default: + done = true; + break; + } + fprintf (stderr, "foo\n"); + } + if (done) break; + } + // dump some useful info if trace is defined if (psTraceGetLevel(FACILITY) >= 6) { @@ -202,5 +267,5 @@ // XXX : Madsen gives suggestion for better use of rho // rho is positive if the new chisq is smaller - if (rho >= -1e-6) { + if (testValue || (rho >= -1e-6)) { min->value = Chisq; alpha = psImageCopy(alpha, Alpha, PS_TYPE_F32); @@ -215,5 +280,5 @@ case 0: if (rho >= -1e-6) { - lambda *= 0.25; + lambda *= 0.1; } else { lambda *= 10.0; @@ -234,8 +299,8 @@ if (rho > 0.0) { lambda *= PS_MAX(0.33, (1.0 - pow(2.0*rho - 1.0, 3.0))); - nu = 2.0; + nu = 3.0; } else { lambda *= nu; - nu *= 2.0; + nu *= 3.0; } break; @@ -474,10 +539,18 @@ // XXX TEST : SAVE IMAGES # if (SAVE_IMAGES) - psphotSaveImage (NULL, pcm->psf->image, "psf.fits"); - psphotSaveImage (NULL, pcm->modelFlux, "model.fits"); - psphotSaveImage (NULL, pcm->modelConvFlux, "modelConv.fits"); - psphotSaveImage (NULL, source->pixels, "obj.fits"); - psphotSaveImage (NULL, source->maskObj, "mask.fits"); - psphotSaveImage (NULL, source->variance, "variance.fits"); + static int Npass = 0; + char name[128]; + snprintf (name, 128, "psf.%03d.fits", Npass); psphotSaveImage (NULL, pcm->psf->image, name); + snprintf (name, 128, "mod.%03d.fits", Npass); psphotSaveImage (NULL, pcm->modelFlux, name); + snprintf (name, 128, "cnv.%03d.fits", Npass); psphotSaveImage (NULL, pcm->modelConvFlux, name); + snprintf (name, 128, "obj.%03d.fits", Npass); psphotSaveImage (NULL, source->pixels, name); + snprintf (name, 128, "msk.%03d.fits", Npass); psphotSaveImage (NULL, source->maskObj, name); + snprintf (name, 128, "var.%03d.fits", Npass); psphotSaveImage (NULL, source->variance, name); + for (int n = 0; n < pcm->dmodelsFlux->n; n++) { + psImage *dmodelConv = pcm->dmodelsConvFlux->data[n]; + if (!dmodelConv) continue; + snprintf (name, 128, "dpar.%01d.%03d.fits", n, Npass); psphotSaveImage (NULL, dmodelConv, name); + } + Npass ++; # endif Index: trunk/psModules/src/objects/pmPCMdata.c =================================================================== --- trunk/psModules/src/objects/pmPCMdata.c (revision 35925) +++ trunk/psModules/src/objects/pmPCMdata.c (revision 36085) @@ -173,4 +173,78 @@ } +int pmPCMsetParams (psMinConstraint *constraint, pmSourceFitMode mode) { + + // set parameter mask based on fitting mode + int nParams = 0; + + switch (mode) { + case PM_SOURCE_FIT_NORM: + // fits only source normalization (Io) + nParams = 1; + psVectorInit (constraint->paramMask, 1); + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0; + break; + + case PM_SOURCE_FIT_PSF: + // fits only x,y,Io + nParams = 3; + psVectorInit (constraint->paramMask, 1); + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0; + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 0; + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 0; + break; + + case PM_SOURCE_FIT_EXT: + // fits all params except sky + nParams = params->n - 1; + psVectorInit (constraint->paramMask, 0); + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1; + break; + + case PM_SOURCE_FIT_EXT_AND_SKY: + // fits all params including sky + nParams = params->n; + psVectorInit (constraint->paramMask, 0); + break; + + case PM_SOURCE_FIT_SHAPE: + // fits shape (Sxx, Sxy, Syy) and Io + nParams = 5; + psVectorInit (constraint->paramMask, 1); + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 0; + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0; + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SXX] = 0; + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SXY] = 0; + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SYY] = 0; + break; + + case PM_SOURCE_FIT_INDEX: + // fits only Io, index (PAR7) -- only Io for models with < 8 params + psVectorInit (constraint->paramMask, 1); + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0; + if (params->n == 7) { + nParams = 1; + } else { + nParams = 2; + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 0; + } + break; + + case PM_SOURCE_FIT_NO_INDEX: + // fits all but index (PAR7) including sky + psVectorInit (constraint->paramMask, 0); + if (params->n == 7) { + nParams = params->n; + } else { + nParams = params->n - 1; + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1; + } + break; + default: + psAbort("invalid fitting mode"); + } + return nParams; +} + pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model, psImageMaskType maskVal, float psfSize) { @@ -219,52 +293,5 @@ constraint->checkLimits = model->modelLimits; - // set parameter mask based on fitting mode - int nParams = 0; - switch (fitOptions->mode) { - case PM_SOURCE_FIT_NORM: - // NORM-only model fits only source normalization (Io) - nParams = 1; - psVectorInit (constraint->paramMask, 1); - constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0; - break; - case PM_SOURCE_FIT_PSF: - // PSF model only fits x,y,Io - nParams = 3; - psVectorInit (constraint->paramMask, 1); - constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0; - constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 0; - constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 0; - break; - case PM_SOURCE_FIT_EXT: - // EXT model fits all params (except sky) - nParams = params->n - 1; - psVectorInit (constraint->paramMask, 0); - constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1; - break; - case PM_SOURCE_FIT_INDEX: - // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params - psVectorInit (constraint->paramMask, 1); - constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0; - if (params->n == 7) { - nParams = 1; - } else { - nParams = 2; - constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 0; - } - break; - case PM_SOURCE_FIT_NO_INDEX: - // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params - psVectorInit (constraint->paramMask, 0); - constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1; - if (params->n == 7) { - nParams = params->n - 1; - } else { - nParams = params->n - 2; - constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1; - } - break; - default: - psAbort("invalid fitting mode"); - } + int nParams = pmPCMsetParams (constraint, fitOptions->mode); if (nPix < nParams + 1) { @@ -341,53 +368,5 @@ } - // if we changed the fit mode, we need to update nDOF - int nParams = 0; - // set parameter mask based on fitting mode - switch (fitOptions->mode) { - case PM_SOURCE_FIT_NORM: - // NORM-only model fits only source normalization (Io) - nParams = 1; - psVectorInit (pcm->constraint->paramMask, 1); - pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0; - break; - case PM_SOURCE_FIT_PSF: - // PSF model only fits x,y,Io - nParams = 3; - psVectorInit (pcm->constraint->paramMask, 1); - pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0; - pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 0; - pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 0; - break; - case PM_SOURCE_FIT_EXT: - // EXT model fits all params (except sky) - nParams = model->params->n - 1; - psVectorInit (pcm->constraint->paramMask, 0); - pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1; - break; - case PM_SOURCE_FIT_INDEX: - // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params - psVectorInit (pcm->constraint->paramMask, 1); - pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0; - if (model->params->n == 7) { - nParams = 1; - } else { - nParams = 2; - pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 0; - } - break; - case PM_SOURCE_FIT_NO_INDEX: - // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params - psVectorInit (pcm->constraint->paramMask, 0); - pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1; - if (model->params->n == 7) { - nParams = model->params->n - 1; - } else { - nParams = model->params->n - 2; - pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1; - } - break; - default: - psAbort("invalid fitting mode"); - } + int nParams = pmPCMsetParams (pcm->constraint, fitOptions->mode); if (pcm->nPix < nParams + 1) { @@ -478,2 +457,54 @@ } +// construct a realization of the source model +bool pmPCMMakeModel (pmSource *source, pmModel *model, psImageMaskType maskVal, int psfSize) { + + PS_ASSERT_PTR_NON_NULL(source, false); + + // if we already have a cached image, re-use that memory + source->modelFlux = psImageCopy (source->modelFlux, source->pixels, PS_TYPE_F32); + psImageInit (source->modelFlux, 0.0); + + // modelFlux always has unity normalization (I0 = 1.0) + pmModelAdd (source->modelFlux, source->maskObj, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal); + + // convolve the model image with the PSF + if (USE_1D_GAUSS) { + // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread): + // * the model flux is not masked + // * threading takes place above this level + + // define the Gauss parameters from the psf + pmModel *modelPSF = source->modelPSF; + psAssert (modelPSF, "psf model must be defined"); + + psEllipseAxes axes; + bool useReff = pmModelUseReff (modelPSF->type); + psF32 *PAR = modelPSF->params->data.F32; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff); + + float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]); + float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); + + float sigma = 0.5 * (FWHM_MAJOR + FWHM_MINOR) / 2.35; + float nsigma = 2.0; + + psImageSmooth (source->modelFlux, sigma, nsigma); + } else { + // make sure we save a cached copy of the psf flux + pmSourceCachePSF (source, maskVal); + + // convert the cached cached psf model for this source to a psKernel + psKernel *psf = pmPCMkernelFromPSF (source, psfSize); + if (!psf) { + // NOTE: this only happens if the source is too close to an edge + model->flags |= PM_MODEL_STATUS_BADARGS; + return NULL; + } + + // XXX not sure if I can place the output on top of the input + psImageConvolveFFT (source->modelFlux, source->modelFlux, NULL, 0, psf); + } + return true; +} + Index: trunk/psModules/src/objects/pmPCMdata.h =================================================================== --- trunk/psModules/src/objects/pmPCMdata.h (revision 35925) +++ trunk/psModules/src/objects/pmPCMdata.h (revision 36085) @@ -98,4 +98,6 @@ bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize); +bool pmPCMMakeModel (pmSource *source, pmModel *model, psImageMaskType maskVal, int psfSize); + /// @} # endif /* PM_PCM_DATA_H */ Index: trunk/psModules/src/objects/pmSourceFitModel.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitModel.c (revision 35925) +++ trunk/psModules/src/objects/pmSourceFitModel.c (revision 36085) @@ -66,4 +66,5 @@ opt->gainFactorMode = 0; opt->chisqConvergence = true; + opt->isInteractive = false; return opt; @@ -247,4 +248,5 @@ myMin->gainFactorMode = options->gainFactorMode; myMin->chisqConvergence = options->chisqConvergence; + myMin->isInteractive = options->isInteractive; psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); Index: trunk/psModules/src/objects/pmSourceFitModel.h =================================================================== --- trunk/psModules/src/objects/pmSourceFitModel.h (revision 35925) +++ trunk/psModules/src/objects/pmSourceFitModel.h (revision 36085) @@ -21,4 +21,5 @@ PM_SOURCE_FIT_EXT_AND_SKY, PM_SOURCE_FIT_INDEX, + PM_SOURCE_FIT_SHAPE, PM_SOURCE_FIT_NO_INDEX, PM_SOURCE_FIT_TRAIL, @@ -37,4 +38,5 @@ int gainFactorMode; bool chisqConvergence; + bool isInteractive; } pmSourceFitOptions; Index: trunk/psModules/src/objects/pmSourceFitPCM.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitPCM.c (revision 35925) +++ trunk/psModules/src/objects/pmSourceFitPCM.c (revision 36085) @@ -68,4 +68,5 @@ myMin->chisqConvergence = fitOptions->chisqConvergence; myMin->gainFactorMode = fitOptions->gainFactorMode; + myMin->isInteractive = fitOptions->isInteractive; psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); Index: trunk/psModules/src/objects/pmSourceFitSet.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitSet.c (revision 35925) +++ trunk/psModules/src/objects/pmSourceFitSet.c (revision 36085) @@ -570,4 +570,5 @@ myMin->gainFactorMode = options->gainFactorMode; myMin->chisqConvergence = options->chisqConvergence; + myMin->isInteractive = options->isInteractive; psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);