Index: trunk/psModules/src/objects/Makefile.am =================================================================== --- trunk/psModules/src/objects/Makefile.am (revision 32046) +++ trunk/psModules/src/objects/Makefile.am (revision 32347) @@ -130,14 +130,14 @@ # pmSourceID_CMF_* functions use a common framework -BUILT_SOURCES = pmSourceIO_CMF_PS1_V1.v1.c pmSourceIO_CMF_PS1_V2.v1.c pmSourceIO_CMF_PS1_V3.v1.c +BUILT_SOURCES = pmSourceIO_CMF_PS1_V1.c pmSourceIO_CMF_PS1_V2.c pmSourceIO_CMF_PS1_V3.c -pmSourceIO_CMF_PS1_V1.v1.c : pmSourceIO_CMF.c.in mksource.pl - mksource.pl pmSourceIO_CMF.c.in PS1_V1 pmSourceIO_CMF_PS1_V1.v1.c +pmSourceIO_CMF_PS1_V1.c : pmSourceIO_CMF.c.in mksource.pl + mksource.pl pmSourceIO_CMF.c.in PS1_V1 pmSourceIO_CMF_PS1_V1.c -pmSourceIO_CMF_PS1_V2.v1.c : pmSourceIO_CMF.c.in mksource.pl - mksource.pl pmSourceIO_CMF.c.in PS1_V2 pmSourceIO_CMF_PS1_V2.v1.c +pmSourceIO_CMF_PS1_V2.c : pmSourceIO_CMF.c.in mksource.pl + mksource.pl pmSourceIO_CMF.c.in PS1_V2 pmSourceIO_CMF_PS1_V2.c -pmSourceIO_CMF_PS1_V3.v1.c : pmSourceIO_CMF.c.in mksource.pl - mksource.pl pmSourceIO_CMF.c.in PS1_V2 pmSourceIO_CMF_PS1_V3.v1.c +pmSourceIO_CMF_PS1_V3.c : pmSourceIO_CMF.c.in mksource.pl + mksource.pl pmSourceIO_CMF.c.in PS1_V3 pmSourceIO_CMF_PS1_V3.c # EXTRA_DIST = pmErrorCodes.h.in pmErrorCodes.dat pmErrorCodes.c.in Index: trunk/psModules/src/objects/mksource.pl =================================================================== --- trunk/psModules/src/objects/mksource.pl (revision 32046) +++ trunk/psModules/src/objects/mksource.pl (revision 32347) @@ -14,7 +14,10 @@ # see if we can add in PS1_DV* and PS1_SV* as well... -@cmfmodes = ("PS1_V1", 1, +%cmfmodes = ("PS1_V1", 1, "PS1_V2", 2, "PS1_V3", 3); + +print "1: $cmfmodes{1}\n"; +print "PS1_V1: $cmfmodes{'PS1_V1'}\n"; open (FILE, "$template") || die "failed to open template $template\n"; @@ -50,15 +53,19 @@ if ($gtMode) { + # print "gtMode : $line\n"; $thisLevel = $cmfmodes{$gtMode}; $myLevel = $cmfmodes{$cmfmode}; - if ($thisLevel <= $myLevel) { next; } - $line =~ s|\@<\S*\@\s*||; + print "gtMode : $gtMode vs $cmfmode, $thisLevel, $myLevel\n"; + if ($myLevel <= $thisLevel) { next; } + $line =~ s|\@>\S*\@\s*||; } if ($ltMode) { + # print "ltMode : $line\n"; $thisLevel = $cmfmodes{$ltMode}; $myLevel = $cmfmodes{$cmfmode}; - if ($thisLevel >= $myLevel) { next; } - $line =~ s|\@>\S*\@\s*||; + print "ltMode : $ltMode vs $cmfmode, $thisLevel, $myLevel\n"; + if ($myLevel >= $thisLevel) { next; } + $line =~ s|\@<\S*\@\s*||; } Index: trunk/psModules/src/objects/models/pmModel_DEV.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_DEV.c (revision 32046) +++ trunk/psModules/src/objects/models/pmModel_DEV.c (revision 32347) @@ -89,4 +89,6 @@ static bool limitsApply = true; // Apply limits? +# include "pmModel_SERSIC.CP.h" + psF32 PM_MODEL_FUNC (psVector *deriv, const psVector *params, @@ -94,8 +96,4 @@ { psF32 *PAR = params->data.F32; - - float index = 0.5 / ALPHA; - float bn = 1.9992*index - 0.3271; - float Io = exp(bn); psF32 X = pixcoord->data.F32[0] - PAR[PM_PAR_XPOS]; @@ -105,8 +103,89 @@ psF32 z = (PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y); - assert (z >= 0); + // If the elliptical contour is defined in a valid way, we should never trigger this + // assert. Other models (like PGAUSS) don't use fractional powers, and thus do not have + // NaN values for negative values of z + 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) + psEllipseShape shape; + + shape.sx = PAR[PM_PAR_SXX]; + shape.sy = PAR[PM_PAR_SYY]; + shape.sxy = PAR[PM_PAR_SXY]; + + // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio + psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); + + // 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; @@ -120,9 +199,9 @@ psF32 *dPAR = deriv->data.F32; + 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); - - dPAR[PM_PAR_SKY] = +1.0; - dPAR[PM_PAR_I0] = +2.0*f1; assert (isfinite(z1)); @@ -223,4 +302,5 @@ // set the shape parameters + // XXX adjust this? if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) { return false; @@ -246,36 +326,28 @@ } +// A DeVaucouleur model is equivalent to a Sersic with index = 4.0 psF64 PM_MODEL_FLUX (const psVector *params) { - float z, norm; psEllipseShape shape; psF32 *PAR = params->data.F32; - shape.sx = PAR[PM_PAR_SXX] / M_SQRT2; - shape.sy = PAR[PM_PAR_SYY] / M_SQRT2; + shape.sx = PAR[PM_PAR_SXX]; + shape.sy = PAR[PM_PAR_SYY]; shape.sxy = PAR[PM_PAR_SXY]; - // Area is equivalent to 2 pi sigma^2 + // for a non-circular DeVaucouleur, the flux of the Rmajor equivalent is scaled by the AspectRatio psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); - psF64 Area = 2.0 * M_PI * axes.major * axes.minor; - - // the area needs to be multiplied by the integral of f(z) - norm = 0.0; - - # define DZ 0.25 - - float f0 = 1.0; - float f1, f2; - for (z = DZ; z < 150; z += DZ) { - f1 = exp(-pow(z,ALPHA)); - z += DZ; - f2 = exp(-pow(z,ALPHA)); - norm += f0 + 4*f1 + f2; - f0 = f2; - } - norm *= DZ / 3.0; - - psF64 Flux = PAR[PM_PAR_I0] * Area * norm; + 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); @@ -297,6 +369,6 @@ return (1.0); - shape.sx = PAR[PM_PAR_SXX] / M_SQRT2; - shape.sy = PAR[PM_PAR_SYY] / M_SQRT2; + shape.sx = PAR[PM_PAR_SXX]; + shape.sy = PAR[PM_PAR_SYY]; shape.sxy = PAR[PM_PAR_SXY]; Index: trunk/psModules/src/objects/models/pmModel_EXP.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_EXP.c (revision 32046) +++ trunk/psModules/src/objects/models/pmModel_EXP.c (revision 32347) @@ -81,4 +81,6 @@ static bool limitsApply = true; // Apply limits? +# include "pmModel_SERSIC.CP.h" + psF32 PM_MODEL_FUNC (psVector *deriv, const psVector *params, @@ -93,13 +95,89 @@ psF32 z = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y; - // XXX if the elliptical contour is defined in valid way, this step should not be required. - // other models (like PGAUSS) don't use fractional powers, and thus do not have NaN values - // for negative values of z - // XXX use an assert here to force the elliptical parameters to be correctly determined - // if (z < 0) z = 0; - assert (z >= 0); - - psF32 f2 = sqrt(z); - psF32 f1 = exp(-f2); + // If the elliptical contour is defined in a valid way, we should never trigger this + // assert. Other models (like PGAUSS) don't use fractional powers, and thus do not have + // NaN values for negative values of z + 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); + + 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) + psEllipseShape shape; + + shape.sx = PAR[PM_PAR_SXX]; + shape.sy = PAR[PM_PAR_SYY]; + shape.sxy = PAR[PM_PAR_SXY]; + + // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio + psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); + + // 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; @@ -118,5 +196,5 @@ // 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*z0/sqrt(0.01) : 0.5*z0/sqrt(z); + 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? @@ -216,4 +294,5 @@ // set the shape parameters + // XXX adjust this? if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) { return false; @@ -233,36 +312,28 @@ } +// An exponential model is equivalent to a Sersic with index = 1.0 psF64 PM_MODEL_FLUX (const psVector *params) { - float z, norm; psEllipseShape shape; psF32 *PAR = params->data.F32; - shape.sx = PAR[PM_PAR_SXX] / M_SQRT2; - shape.sy = PAR[PM_PAR_SYY] / M_SQRT2; + shape.sx = PAR[PM_PAR_SXX]; + shape.sy = PAR[PM_PAR_SYY]; shape.sxy = PAR[PM_PAR_SXY]; - // Area is equivalent to 2 pi sigma^2 + // for a non-circular Exponential, the flux of the Rmajor equivalent is scaled by the AspectRatio psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); - psF64 Area = 2.0 * M_PI * axes.major * axes.minor; - - // the area needs to be multiplied by the integral of f(z) = exp(-sqrt(z)) [0 to infinity] - norm = 0.0; - - # define DZ 0.25 - - float f0 = 1.0; - float f1, f2; - for (z = DZ; z < 150; z += DZ) { - f1 = exp(-sqrt(z)); - z += DZ; - f2 = exp(-sqrt(z)); - norm += f0 + 4*f1 + f2; - f0 = f2; - } - norm *= DZ / 3.0; - - psF64 Flux = PAR[PM_PAR_I0] * Area * norm; + 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); @@ -284,6 +355,6 @@ return (1.0); - shape.sx = PAR[PM_PAR_SXX] / M_SQRT2; - shape.sy = PAR[PM_PAR_SYY] / M_SQRT2; + shape.sx = PAR[PM_PAR_SXX]; + shape.sy = PAR[PM_PAR_SYY]; shape.sxy = PAR[PM_PAR_SXY]; Index: trunk/psModules/src/objects/models/pmModel_PGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 32046) +++ trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 32347) @@ -217,7 +217,8 @@ } -psF64 PM_MODEL_FLUX(const psVector *params) -{ - float norm, z; +// integrate the model to get the full flux +psF64 PM_MODEL_FLUX (const psVector *params) +{ + float z, norm; psEllipseShape shape; @@ -228,25 +229,27 @@ shape.sxy = PAR[PM_PAR_SXY]; - // Area is equivalent to 2 pi sigma^2 psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); - psF64 Area = 2.0 * M_PI * axes.major * axes.minor; - - // the area needs to be multiplied by the integral of f(z) + float AspectRatio = axes.minor / axes.major; + + // flux = 2 \pi \int f(r) r dr norm = 0.0; - # define DZ 0.25 - - float f0 = 1.0; + # define DR 0.25 + + // f = f(r) * r + float f0 = 0.0; float f1, f2; - for (z = DZ; z < 150; z += DZ) { - f1 = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0); - z += DZ; - f2 = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0); + for (float r = DR; r < 150; r += DR) { + z = 0.5 * PS_SQR(r / axes.major); + f1 = r / (1 + z + z*z/2.0 + z*z*z/6.0); + r += DR; + z = 0.5 * PS_SQR(r / axes.major); + f2 = r / (1 + z + z*z/2.0 + z*z*z/6.0); norm += f0 + 4*f1 + f2; f0 = f2; } - norm *= DZ / 3.0; - - psF64 Flux = PAR[PM_PAR_I0] * Area * norm; + norm *= DR / 3.0; + + psF64 Flux = PAR[PM_PAR_I0] * norm * 2.0 * M_PI * AspectRatio; return(Flux); Index: trunk/psModules/src/objects/models/pmModel_PS1_V1.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_PS1_V1.c (revision 32046) +++ trunk/psModules/src/objects/models/pmModel_PS1_V1.c (revision 32347) @@ -14,6 +14,6 @@ * PM_PAR_XPOS 2 - X center of object * PM_PAR_YPOS 3 - Y center of object - * PM_PAR_SXX 4 - X^2 term of elliptical contour (sqrt(2) / SigmaX) - * PM_PAR_SYY 5 - Y^2 term of elliptical contour (sqrt(2) / SigmaY) + * PM_PAR_SXX 4 - X^2 term of elliptical contour (SigmaX / sqrt(2)) + * PM_PAR_SYY 5 - Y^2 term of elliptical contour (SigmaY / sqrt(2)) * PM_PAR_SXY 6 - X*Y term of elliptical contour * PM_PAR_7 7 - amplitude of the linear component (k) @@ -239,4 +239,5 @@ } +// integrate the model to get the full flux psF64 PM_MODEL_FLUX (const psVector *params) { @@ -250,25 +251,27 @@ shape.sxy = PAR[PM_PAR_SXY]; - // Area is equivalent to 2 pi sigma^2 psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); - psF64 Area = 2.0 * M_PI * axes.major * axes.minor; - - // the area needs to be multiplied by the integral of f(z) + float AspectRatio = axes.minor / axes.major; + + // flux = 2 \pi \int f(r) r dr norm = 0.0; - # define DZ 0.25 - - float f0 = 1.0; + # define DR 0.25 + + // f = f(r) * r + float f0 = 0.0; float f1, f2; - for (z = DZ; z < 150; z += DZ) { - f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA)); - z += DZ; - f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA)); + for (float r = DR; r < 150; r += DR) { + z = 0.5 * PS_SQR(r / axes.major); + f1 = r / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA)); + r += DR; + z = 0.5 * PS_SQR(r / axes.major); + f2 = r / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA)); norm += f0 + 4*f1 + f2; f0 = f2; } - norm *= DZ / 3.0; - - psF64 Flux = PAR[PM_PAR_I0] * Area * norm; + norm *= DR / 3.0; + + psF64 Flux = PAR[PM_PAR_I0] * norm * 2.0 * M_PI * AspectRatio; return(Flux); Index: trunk/psModules/src/objects/models/pmModel_QGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 32046) +++ trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 32347) @@ -240,4 +240,5 @@ } +// integrate the model to get the full flux psF64 PM_MODEL_FLUX (const psVector *params) { @@ -251,25 +252,27 @@ shape.sxy = PAR[PM_PAR_SXY]; - // Area is equivalent to 2 pi sigma^2 psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); - psF64 Area = 2.0 * M_PI * axes.major * axes.minor; - - // the area needs to be multiplied by the integral of f(z) + float AspectRatio = axes.minor / axes.major; + + // flux = 2 \pi \int f(r) r dr norm = 0.0; - # define DZ 0.25 - - float f0 = 1.0; + # define DR 0.25 + + // f = f(r) * r + float f0 = 0.0; float f1, f2; - for (z = DZ; z < 150; z += DZ) { - f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA)); - z += DZ; - f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA)); + for (float r = DR; r < 150; r += DR) { + z = 0.5 * PS_SQR(r / axes.major); + f1 = r / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA)); + r += DR; + z = 0.5 * PS_SQR(r / axes.major); + f2 = r / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA)); norm += f0 + 4*f1 + f2; f0 = f2; } - norm *= DZ / 3.0; - - psF64 Flux = PAR[PM_PAR_I0] * Area * norm; + norm *= DR / 3.0; + + psF64 Flux = PAR[PM_PAR_I0] * norm * 2.0 * M_PI * AspectRatio; return(Flux); Index: trunk/psModules/src/objects/models/pmModel_RGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 32046) +++ trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 32347) @@ -229,7 +229,8 @@ } +// integrate the model to get the full flux psF64 PM_MODEL_FLUX (const psVector *params) { - float norm, z; + float z, norm; psEllipseShape shape; @@ -240,25 +241,27 @@ shape.sxy = PAR[PM_PAR_SXY]; - // Area is equivalent to 2 pi sigma^2 psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); - psF64 Area = 2.0 * M_PI * axes.major * axes.minor; - - // the area needs to be multiplied by the integral of f(z) + float AspectRatio = axes.minor / axes.major; + + // flux = 2 \pi \int f(r) r dr norm = 0.0; - # define DZ 0.25 - - float f0 = 1.0; + # define DR 0.25 + + // f = f(r) * r + float f0 = 0.0; float f1, f2; - for (z = DZ; z < 150; z += DZ) { - f1 = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7])); - z += DZ; - f2 = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7])); + for (float r = DR; r < 150; r += DR) { + z = 0.5 * PS_SQR(r / axes.major); + f1 = r / (1 + z + pow(z, PAR[PM_PAR_7])); + r += DR; + z = 0.5 * PS_SQR(r / axes.major); + f2 = r / (1 + z + pow(z, PAR[PM_PAR_7])); norm += f0 + 4*f1 + f2; f0 = f2; } - norm *= DZ / 3.0; - - psF64 Flux = PAR[PM_PAR_I0] * Area * norm; + norm *= DR / 3.0; + + psF64 Flux = PAR[PM_PAR_I0] * norm * 2.0 * M_PI * AspectRatio; return(Flux); Index: trunk/psModules/src/objects/models/pmModel_SERSIC.CP.h =================================================================== --- trunk/psModules/src/objects/models/pmModel_SERSIC.CP.h (revision 32347) +++ trunk/psModules/src/objects/models/pmModel_SERSIC.CP.h (revision 32347) @@ -0,0 +1,30 @@ +static int centralPixelNX = 29; +static int centralPixelNY = 12; +static float centralPixelXo = 2.0; +static float centralPixelYo = 0.5; +static float centralPixeldX = 0.5; +static float centralPixeldY = 0.5; + +static float centralPixel[12][29] = { +{ 0.0521, 0.0336, 0.0235, 0.0173, 0.0133, 0.0105, 0.0085, 0.0070, 0.0059, 0.0051, 0.0044, 0.0038, 0.0033, 0.0030, 0.0026, 0.0024, 0.0021, 0.0019, 0.0018, 0.0016, 0.0015, 0.0014, 0.0013, 0.0012, 0.0011, 0.0010, 0.0010, 0.0009, 0.0008, }, +{ 0.0817, 0.0556, 0.0402, 0.0304, 0.0238, 0.0191, 0.0157, 0.0131, 0.0111, 0.0096, 0.0083, 0.0073, 0.0064, 0.0057, 0.0051, 0.0046, 0.0042, 0.0038, 0.0035, 0.0032, 0.0029, 0.0027, 0.0025, 0.0023, 0.0022, 0.0020, 0.0019, 0.0018, 0.0017, }, +{ 0.1096, 0.0785, 0.0591, 0.0462, 0.0371, 0.0305, 0.0255, 0.0217, 0.0187, 0.0162, 0.0143, 0.0126, 0.0113, 0.0101, 0.0091, 0.0083, 0.0076, 0.0069, 0.0064, 0.0059, 0.0054, 0.0050, 0.0047, 0.0044, 0.0041, 0.0038, 0.0036, 0.0034, 0.0032, }, +{ 0.1339, 0.0998, 0.0777, 0.0624, 0.0513, 0.0431, 0.0367, 0.0317, 0.0277, 0.0244, 0.0217, 0.0194, 0.0175, 0.0158, 0.0144, 0.0132, 0.0121, 0.0112, 0.0103, 0.0096, 0.0089, 0.0083, 0.0078, 0.0073, 0.0069, 0.0065, 0.0061, 0.0058, 0.0055, }, +{ 0.1547, 0.1190, 0.0951, 0.0781, 0.0655, 0.0559, 0.0484, 0.0423, 0.0374, 0.0334, 0.0300, 0.0271, 0.0246, 0.0225, 0.0206, 0.0190, 0.0176, 0.0163, 0.0152, 0.0142, 0.0133, 0.0125, 0.0117, 0.0110, 0.0104, 0.0099, 0.0093, 0.0089, 0.0084, }, +{ 0.1726, 0.1361, 0.1110, 0.0928, 0.0791, 0.0685, 0.0600, 0.0531, 0.0475, 0.0427, 0.0387, 0.0353, 0.0323, 0.0297, 0.0275, 0.0255, 0.0237, 0.0221, 0.0207, 0.0195, 0.0183, 0.0172, 0.0163, 0.0154, 0.0146, 0.0139, 0.0132, 0.0126, 0.0120, }, +{ 0.1881, 0.1513, 0.1255, 0.1065, 0.0919, 0.0805, 0.0713, 0.0637, 0.0574, 0.0521, 0.0476, 0.0437, 0.0403, 0.0373, 0.0347, 0.0323, 0.0303, 0.0284, 0.0267, 0.0252, 0.0238, 0.0225, 0.0214, 0.0203, 0.0193, 0.0184, 0.0176, 0.0168, 0.0161, }, +{ 0.2017, 0.1649, 0.1387, 0.1191, 0.1040, 0.0919, 0.0821, 0.0740, 0.0672, 0.0614, 0.0564, 0.0521, 0.0483, 0.0450, 0.0420, 0.0394, 0.0370, 0.0349, 0.0329, 0.0312, 0.0296, 0.0281, 0.0268, 0.0255, 0.0244, 0.0233, 0.0223, 0.0214, 0.0205, }, +{ 0.2138, 0.1771, 0.1507, 0.1308, 0.1152, 0.1027, 0.0924, 0.0839, 0.0767, 0.0705, 0.0651, 0.0604, 0.0563, 0.0527, 0.0494, 0.0465, 0.0439, 0.0415, 0.0393, 0.0374, 0.0356, 0.0339, 0.0324, 0.0310, 0.0296, 0.0284, 0.0273, 0.0262, 0.0252, }, +{ 0.2246, 0.1883, 0.1618, 0.1416, 0.1257, 0.1128, 0.1022, 0.0933, 0.0857, 0.0792, 0.0735, 0.0686, 0.0642, 0.0603, 0.0568, 0.0536, 0.0508, 0.0482, 0.0458, 0.0436, 0.0416, 0.0398, 0.0381, 0.0365, 0.0351, 0.0337, 0.0324, 0.0313, 0.0301, }, +{ 0.2347, 0.1986, 0.1721, 0.1517, 0.1355, 0.1224, 0.1115, 0.1023, 0.0945, 0.0877, 0.0817, 0.0765, 0.0719, 0.0677, 0.0640, 0.0606, 0.0576, 0.0548, 0.0522, 0.0499, 0.0477, 0.0458, 0.0439, 0.0422, 0.0406, 0.0391, 0.0377, 0.0364, 0.0352, }, +{ 0.2448, 0.2086, 0.1819, 0.1614, 0.1450, 0.1316, 0.1204, 0.1110, 0.1029, 0.0958, 0.0897, 0.0842, 0.0793, 0.0750, 0.0711, 0.0675, 0.0643, 0.0613, 0.0586, 0.0561, 0.0538, 0.0517, 0.0497, 0.0479, 0.0462, 0.0445, 0.0430, 0.0416, 0.0403, }, +}; + +static float interpolatePixels (float V00, float V10, float V01, float V11, float dx, float dy) { + // bilinear interpolation + float rx = 1.0 - dx; + float ry = 1.0 - dy; + float value = V00*rx*ry + V10*dx*ry + V01*rx*dy + V11*dx*dy; + return value; +} + Index: trunk/psModules/src/objects/models/pmModel_SERSIC.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 32046) +++ trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 32347) @@ -13,8 +13,14 @@ * PM_PAR_XPOS 2 - X center of object * PM_PAR_YPOS 3 - Y center of object - * PM_PAR_SXX 4 - X^2 term of elliptical contour (sqrt(2) / SigmaX) - * PM_PAR_SYY 5 - Y^2 term of elliptical contour (sqrt(2) / SigmaY) + * PM_PAR_SXX 4 - X^2 term of elliptical contour (SigmaX / sqrt(2)) + * PM_PAR_SYY 5 - Y^2 term of elliptical contour (SigmaY / sqrt(2)) * PM_PAR_SXY 6 - X*Y term of elliptical contour * PM_PAR_7 7 - normalized sersic parameter + + * 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: + 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) note that a standard sersic model uses exp(-K*(z^(1/2n) - 1). the additional elements (K, @@ -85,4 +91,6 @@ static bool limitsApply = true; // Apply limits? +# include "pmModel_SERSIC.CP.h" + psF32 PM_MODEL_FUNC (psVector *deriv, const psVector *params, @@ -97,17 +105,89 @@ psF32 z = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y; - // XXX if the elliptical contour is defined in valid way, this step should not be required. - // other models (like PGAUSS) don't use fractional powers, and thus do not have NaN values - // for negative values of z - // XXX use an assert here to force the elliptical parameters to be correctly determined - // if (z < 0) z = 0; - assert (z >= 0); + // If the elliptical contour is defined in a valid way, we should never trigger this + // assert. Other models (like PGAUSS) don't use fractional powers, and thus do not have + // NaN values for negative values of z + 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,PAR[PM_PAR_7]); + psF32 f2 = bn*pow(z,par7); 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) + psEllipseShape shape; + + shape.sx = PAR[PM_PAR_SXX]; + shape.sy = PAR[PM_PAR_SYY]; + shape.sxy = PAR[PM_PAR_SXY]; + + // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio + psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); + + // 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; @@ -121,10 +201,11 @@ psF32 *dPAR = deriv->data.F32; - // gradient is infinite for z = 0; saturate at z = 0.01 - psF32 z1 = (z < 0.01) ? z0*bn*PAR[PM_PAR_7]*pow(0.01,PAR[PM_PAR_7] - 1.0) : z0*bn*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7] - 1.0); - dPAR[PM_PAR_SKY] = +1.0; dPAR[PM_PAR_I0] = +f1; - dPAR[PM_PAR_7] = (z < 0.01) ? -z0*pow(0.01,PAR[PM_PAR_7])*log(0.01) : -z0*f2*log(z); + + // 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; @@ -269,8 +350,4 @@ float Io = exp(0.5*bn); - // XXX do we need this factor of sqrt(2)? - // float Sxx = PS_MAX(0.5, M_SQRT2*shape.sx); - // float Syy = PS_MAX(0.5, M_SQRT2*shape.sy); - float Sxx = PS_MAX(0.5, shape.sx); float Syy = PS_MAX(0.5, shape.sy); @@ -294,37 +371,28 @@ } +// A sersic model has an integral flux which can be analytically represented psF64 PM_MODEL_FLUX (const psVector *params) { - float z, norm; psEllipseShape shape; psF32 *PAR = params->data.F32; - shape.sx = PAR[PM_PAR_SXX] / M_SQRT2; - shape.sy = PAR[PM_PAR_SYY] / M_SQRT2; + shape.sx = PAR[PM_PAR_SXX]; + shape.sy = PAR[PM_PAR_SYY]; shape.sxy = PAR[PM_PAR_SXY]; - // Area is equivalent to 2 pi sigma^2 + // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); - psF64 Area = 2.0 * M_PI * axes.major * axes.minor; - - // the area needs to be multiplied by the integral of f(z) - norm = 0.0; - - # define DZ 0.25 - - float f0 = 1.0; - float f1, f2; - for (z = DZ; z < 150; z += DZ) { - // f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25)); - f1 = exp(-pow(z,PAR[PM_PAR_7])); - z += DZ; - f2 = exp(-pow(z,PAR[PM_PAR_7])); - norm += f0 + 4*f1 + f2; - f0 = f2; - } - norm *= DZ / 3.0; - - psF64 Flux = PAR[PM_PAR_I0] * Area * norm; + 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); @@ -346,6 +414,6 @@ return (1.0); - shape.sx = PAR[PM_PAR_SXX] / M_SQRT2; - shape.sy = PAR[PM_PAR_SYY] / M_SQRT2; + shape.sx = PAR[PM_PAR_SXX]; + shape.sy = PAR[PM_PAR_SYY]; shape.sxy = PAR[PM_PAR_SXY]; Index: trunk/psModules/src/objects/pmFootprintCullPeaks.c =================================================================== --- trunk/psModules/src/objects/pmFootprintCullPeaks.c (revision 32046) +++ trunk/psModules/src/objects/pmFootprintCullPeaks.c (revision 32347) @@ -181,4 +181,5 @@ if (!myFP->n) { psWarning ("empty footprint?"); + psFree (myFP); continue; } Index: trunk/psModules/src/objects/pmModelUtils.c =================================================================== --- trunk/psModules/src/objects/pmModelUtils.c (revision 32046) +++ trunk/psModules/src/objects/pmModelUtils.c (revision 32347) @@ -122,4 +122,9 @@ if (!isfinite(axes.theta)) return false; + // Mxx, Mxy, Myy define the elliptical shape, but Mrf defines the width + float scale = (isfinite(moments->Mrf) && (moments->Mrf > 0.0)) ? moments->Mrf / axes.major : 1.0; + axes.major *= scale; + axes.minor *= scale; + psEllipseShape shape = psEllipseAxesToShape (axes); Index: trunk/psModules/src/objects/pmMoments.c =================================================================== --- trunk/psModules/src/objects/pmMoments.c (revision 32046) +++ trunk/psModules/src/objects/pmMoments.c (revision 32347) @@ -32,4 +32,8 @@ tmp->Mrf = 0.0; tmp->Mrh = 0.0; + + tmp->KronCore = 0.0; + tmp->KronCoreErr = 0.0; + tmp->KronFlux = 0.0; tmp->KronFluxErr = 0.0; @@ -37,4 +41,8 @@ tmp->KronFinner = 0.0; tmp->KronFouter = 0.0; + + tmp->KronFluxPSF = 0.0; + tmp->KronFluxPSFErr = 0.0; + tmp->KronRadiusPSF = 0.0; tmp->Mx = 0.0; Index: trunk/psModules/src/objects/pmMoments.h =================================================================== --- trunk/psModules/src/objects/pmMoments.h (revision 32046) +++ trunk/psModules/src/objects/pmMoments.h (revision 32347) @@ -51,4 +51,8 @@ int nPixels; ///< Number of pixels used. + float KronFluxPSF; ///< Kron Flux using PSF-optimized window + float KronFluxPSFErr; ///< Kron Flux Error using PSF-optimized window + float KronRadiusPSF; ///< Kron Radius using PSF-optimized window (Flux in 2.5 Radius) + float KronCore; ///< flux in r < 1.0*Mrf float KronCoreErr; ///< error on flux in r < 1.0*Mrf Index: trunk/psModules/src/objects/pmPCM_MinimizeChisq.c =================================================================== --- trunk/psModules/src/objects/pmPCM_MinimizeChisq.c (revision 32046) +++ trunk/psModules/src/objects/pmPCM_MinimizeChisq.c (revision 32347) @@ -45,4 +45,5 @@ # define USE_FFT 1 # define PRE_CONVOLVE 1 +# define TESTCOPY 0 bool pmPCM_MinimizeChisq ( @@ -93,5 +94,9 @@ psFree (pcm->psfFFT); } - pcm->psfFFT = psImageConvolveKernelInit(pcm->modelFlux, pcm->psf); +# if (!TESTCOPY) + if (!pcm->use1Dgauss) { + pcm->psfFFT = psImageConvolveKernelInit(pcm->modelFlux, pcm->psf); + } +# endif # endif @@ -285,6 +290,6 @@ // Convert i/j to image space: - coord->data.F32[0] = (psF32) (j + source->pixels->col0); - coord->data.F32[1] = (psF32) (i + source->pixels->row0); + coord->data.F32[0] = (psF32) (j + 0.5 + source->pixels->col0); + coord->data.F32[1] = (psF32) (i + 0.5 + source->pixels->row0); pcm->modelFlux->data.F32[i][j] = pcm->modelConv->modelFunc (deriv, params, coord); @@ -309,5 +314,20 @@ # if (PRE_CONVOLVE) // convolve model image and derivative images with pre-convolved kernel - psImageConvolveKernel (pcm->modelConvFlux, pcm->modelFlux, NULL, 0, pcm->psfFFT); + +// XXX for a test, just copy, rather than convolve +# if (TESTCOPY) + psImageCopy (pcm->modelConvFlux, pcm->modelFlux, pcm->modelFlux->type.type); +# else + if (pcm->use1Dgauss) { + // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread): + // * the model flux is not masked + // * threading takes place above this level + pcm->modelConvFlux = psImageCopy (pcm->modelConvFlux, pcm->modelFlux, pcm->modelFlux->type.type); + psImageSmooth (pcm->modelConvFlux, pcm->sigma, pcm->nsigma); + } else { + psImageConvolveKernel (pcm->modelConvFlux, pcm->modelFlux, NULL, 0, pcm->psfFFT); + } +# endif + for (int n = 0; n < pcm->dmodelsFlux->n; n++) { if (pcm->dmodelsFlux->data[n] == NULL) continue; @@ -315,5 +335,17 @@ psImage *dmodel = pcm->dmodelsFlux->data[n]; psImage *dmodelConv = pcm->dmodelsConvFlux->data[n]; - psImageConvolveKernel (dmodelConv, dmodel, NULL, 0, pcm->psfFFT); +# if (TESTCOPY) + psImageCopy (dmodelConv, dmodel, dmodel->type.type); +# else + if (pcm->use1Dgauss) { + // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread): + // * the model flux is not masked + // * threading takes place above this level + dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type); + psImageSmooth (dmodelConv, pcm->sigma, pcm->nsigma); + } else { + psImageConvolveKernel (dmodelConv, dmodel, NULL, 0, pcm->psfFFT); + } +# endif } # else @@ -325,5 +357,14 @@ psImage *dmodel = pcm->dmodelsFlux->data[n]; psImage *dmodelConv = pcm->dmodelsConvFlux->data[n]; - psImageConvolveFFT (dmodelConv, dmodel, NULL, 0, pcm->psf); + + if (pcm->use1Dgauss) { + // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread): + // * the model flux is not masked + // * threading takes place above this level + dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type); + psImageSmooth (dmodelConv, pcm->sigma, pcm->nsigma); + } else { + psImageConvolveFFT (dmodelConv, dmodel, NULL, 0, pcm->psf); + } } # endif // PRE-CONVOLVE Index: trunk/psModules/src/objects/pmPCMdata.c =================================================================== --- trunk/psModules/src/objects/pmPCMdata.c (revision 32046) +++ trunk/psModules/src/objects/pmPCMdata.c (revision 32347) @@ -41,4 +41,7 @@ #include "pmPCMdata.h" +# define USE_DELTA_PSF 0 +# define USE_1D_GAUSS 1 + static void pmPCMdataFree (pmPCMdata *pcm) { @@ -88,4 +91,10 @@ pcm->constraint = NULL; pcm->nDOF = 0; + + // full convolution with the PSF is expensive. if we have to save time, we can do a 1D + // convolution with a Gaussian approximation to the kernel + pcm->use1Dgauss = false; + pcm->nsigma = 3.0; + pcm->sigma = 1.0; // this should be set to something sensible when the psf is known return pcm; @@ -257,4 +266,24 @@ pcm->nDOF = nPix - nParams; +# if (USE_1D_GAUSS) + pmModel *modelPSF = source->modelPSF; + psAssert (modelPSF, "psf model must be defined"); + + psEllipseShape shape; + psEllipseAxes axes; + + shape.sx = modelPSF->params->data.F32[PM_PAR_SXX]; + shape.sy = modelPSF->params->data.F32[PM_PAR_SYY]; + shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY]; + axes = psEllipseShapeToAxes (shape, 20.0); + + float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*modelPSF->params->data.F32[PM_PAR_I0]); + float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); + + pcm->use1Dgauss = true; + pcm->sigma = 0.5 * (FWHM_MAJOR + FWHM_MINOR) / 2.35; + pcm->nsigma = 2.0; +# endif + return pcm; } @@ -368,2 +397,62 @@ return true; } + +// construct a realization of the source model +bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize) { + + PS_ASSERT_PTR_NON_NULL(source, false); + + // select appropriate model + pmModel *model = pmSourceGetModel (NULL, source); + if (model == NULL) return false; // model must be defined + + // 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"); + + psEllipseShape shape; + psEllipseAxes axes; + + shape.sx = modelPSF->params->data.F32[PM_PAR_SXX]; + shape.sy = modelPSF->params->data.F32[PM_PAR_SYY]; + shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY]; + axes = psEllipseShapeToAxes (shape, 20.0); + + float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*modelPSF->params->data.F32[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 32046) +++ trunk/psModules/src/objects/pmPCMdata.h (revision 32347) @@ -35,4 +35,8 @@ int nPar; int nDOF; + + bool use1Dgauss; + float sigma; + float nsigma; } pmPCMdata; @@ -90,4 +94,5 @@ bool pmSourceFitPCM (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize); +bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize); /// @} Index: trunk/psModules/src/objects/pmPSF.c =================================================================== --- trunk/psModules/src/objects/pmPSF.c (revision 32046) +++ trunk/psModules/src/objects/pmPSF.c (revision 32347) @@ -281,4 +281,5 @@ // convert the parameters used in the fitted source model // to the parameters used in the 2D PSF model +// XXX this function may be invalid for SERSIC, DEV, EXP models (SQRT2 not used?) bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis) { @@ -307,4 +308,5 @@ // convert the PSF parameters used in the 2D PSF model fit into the // parameters used in the source model +// XXX this function may be invalid for SERSIC, DEV, EXP models (SQRT2 not used?) psEllipsePol pmPSF_ModelToFit (psF32 *modelPar) { @@ -329,5 +331,5 @@ // convert the parameters used in the fitted source model to the psEllipseAxes representation // (major,minor,theta) -psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, double maxAR) +psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, double maxAR, pmModelType type) { psEllipseShape shape; @@ -336,10 +338,21 @@ axes.minor = NAN; axes.theta = NAN; -// XXX: must assert non-NULL input parameter + // XXX: must assert non-NULL input parameter PS_ASSERT_PTR_NON_NULL(modelPar, axes); - shape.sx = modelPar[PM_PAR_SXX] / M_SQRT2; - shape.sy = modelPar[PM_PAR_SYY] / M_SQRT2; - shape.sxy = modelPar[PM_PAR_SXY]; + bool useReff = true; + useReff |= (type == pmModelClassGetType ("PS_MODEL_SERSIC")); + useReff |= (type == pmModelClassGetType ("PS_MODEL_DEV")); + useReff |= (type == pmModelClassGetType ("PS_MODEL_EXP")); + + if (useReff) { + shape.sx = modelPar[PM_PAR_SXX]; + shape.sy = modelPar[PM_PAR_SYY]; + shape.sxy = modelPar[PM_PAR_SXY]; + } else { + shape.sx = modelPar[PM_PAR_SXX] / M_SQRT2; + shape.sy = modelPar[PM_PAR_SYY] / M_SQRT2; + shape.sxy = modelPar[PM_PAR_SXY]; + } if ((shape.sx == 0) || (shape.sy == 0)) { @@ -357,5 +370,5 @@ // convert the psEllipseAxes representation (major,minor,theta) to the parameters used in the // fitted source model -bool pmPSF_AxesToModel (psF32 *modelPar, psEllipseAxes axes) +bool pmPSF_AxesToModel (psF32 *modelPar, psEllipseAxes axes, pmModelType type) { PS_ASSERT_PTR_NON_NULL(modelPar, false); @@ -370,8 +383,18 @@ psEllipseShape shape = psEllipseAxesToShape (axes); - modelPar[PM_PAR_SXX] = shape.sx * M_SQRT2; - modelPar[PM_PAR_SYY] = shape.sy * M_SQRT2; - modelPar[PM_PAR_SXY] = shape.sxy; - + bool useReff = true; + useReff |= pmModelClassGetType ("PS_MODEL_SERSIC"); + useReff |= pmModelClassGetType ("PS_MODEL_DEV"); + useReff |= pmModelClassGetType ("PS_MODEL_EXP"); + + if (useReff) { + modelPar[PM_PAR_SXX] = shape.sx; + modelPar[PM_PAR_SYY] = shape.sy; + modelPar[PM_PAR_SXY] = shape.sxy; + } else { + modelPar[PM_PAR_SXX] = shape.sx * M_SQRT2; + modelPar[PM_PAR_SYY] = shape.sy * M_SQRT2; + modelPar[PM_PAR_SXY] = shape.sxy; + } return true; } @@ -438,5 +461,5 @@ # if (0) psF32 *params = model->params->data.F32; // Model parameters - psEllipseAxes axes = pmPSF_ModelToAxes(params, MAX_AXIS_RATIO); // Ellipse axes + psEllipseAxes axes = pmPSF_ModelToAxes(params, MAX_AXIS_RATIO, model->type); // Ellipse axes // Curiously, the minor axis can be larger than the major axis, so need to check. Index: trunk/psModules/src/objects/pmPSF.h =================================================================== --- trunk/psModules/src/objects/pmPSF.h (revision 32046) +++ trunk/psModules/src/objects/pmPSF.h (revision 32347) @@ -106,9 +106,9 @@ pmPSF *pmPSFBuildSimple (char *typeName, float sxx, float syy, float sxy, ...); -bool pmPSF_AxesToModel (psF32 *modelPar, psEllipseAxes axes); +bool pmPSF_AxesToModel (psF32 *modelPar, psEllipseAxes axes, pmModelType type); bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis); psEllipsePol pmPSF_ModelToFit (psF32 *modelPar); -psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, double maxAR); +psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, double maxAR, pmModelType type); /// Calculate FWHM value from a PSF Index: trunk/psModules/src/objects/pmPeaks.c =================================================================== --- trunk/psModules/src/objects/pmPeaks.c (revision 32046) +++ trunk/psModules/src/objects/pmPeaks.c (revision 32347) @@ -137,5 +137,9 @@ *****************************************************************************/ static void peakFree(pmPeak *tmp) -{} // +{ + if (!tmp) return; + psFree (tmp->saddlePoints); + return; +} pmPeak *pmPeakAlloc(psS32 x, @@ -162,4 +166,5 @@ tmp->type = type; tmp->footprint = NULL; + tmp->saddlePoints = NULL; psMemSetDeallocator(tmp, (psFreeFunc) peakFree); Index: trunk/psModules/src/objects/pmPeaks.h =================================================================== --- trunk/psModules/src/objects/pmPeaks.h (revision 32046) +++ trunk/psModules/src/objects/pmPeaks.h (revision 32347) @@ -72,4 +72,5 @@ pmPeakType type; ///< Description of peak. pmFootprint *footprint; ///< reference to containing footprint + psArray *saddlePoints; ///< set of saddle points between this peak and near neighbors } pmPeak; Index: trunk/psModules/src/objects/pmSource.c =================================================================== --- trunk/psModules/src/objects/pmSource.c (revision 32046) +++ trunk/psModules/src/objects/pmSource.c (revision 32347) @@ -58,4 +58,5 @@ psFree(tmp->modelEXT); psFree(tmp->modelFits); + psFree(tmp->extFitPars); psFree(tmp->extpars); psFree(tmp->moments); @@ -119,4 +120,5 @@ source->modelEXT = NULL; source->modelFits = NULL; + source->extFitPars = NULL; source->type = PM_SOURCE_TYPE_UNKNOWN; source->mode = PM_SOURCE_MODE_DEFAULT; @@ -196,4 +198,5 @@ // NOTE : because of the const id element, we cannot just assign *source = *in + source->imageID = in->imageID; source->type = in->type; source->mode = in->mode; Index: trunk/psModules/src/objects/pmSource.h =================================================================== --- trunk/psModules/src/objects/pmSource.h (revision 32046) +++ trunk/psModules/src/objects/pmSource.h (revision 32347) @@ -36,4 +36,6 @@ PM_SOURCE_TMPF_MOMENTS_MEASURED = 0x0010, PM_SOURCE_TMPF_CANDIDATE_PSFSTAR = 0x0020, + PM_SOURCE_TMPF_STACK_KEEP = 0x0040, + PM_SOURCE_TMPF_STACK_SKIP = 0x0080, } pmSourceTmpF; @@ -76,4 +78,5 @@ pmModel *modelEXT; ///< EXT Model fit used for subtraction (parameters and type) psArray *modelFits; ///< collection of extended source models (best == modelEXT) + psArray *extFitPars; ///< extra extended fit parameters pmSourceType type; ///< Best identification of object. pmSourceMode mode; ///< analysis flags set for object. Index: trunk/psModules/src/objects/pmSourceExtendedPars.c =================================================================== --- trunk/psModules/src/objects/pmSourceExtendedPars.c (revision 32046) +++ trunk/psModules/src/objects/pmSourceExtendedPars.c (revision 32347) @@ -258,2 +258,27 @@ return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceExtendedFluxFree); } + +// *** pmSourceExtFitPars describes extra metadata related to an extended fit +static void pmSourceExtFitParsFree (pmSourceExtFitPars *pars) { + return; +} + +pmSourceExtFitPars *pmSourceExtFitParsAlloc (void) { + + pmSourceExtFitPars *pars = (pmSourceExtFitPars *) psAlloc(sizeof(pmSourceExtFitPars)); + psMemSetDeallocator(pars, (psFreeFunc) pmSourceExtFitParsFree); + + pars->Mxx = NAN; + pars->Mxy = NAN; + pars->Myy = NAN; + + pars->Mrf = NAN; + pars->Mrh = NAN; + + pars->apMag = NAN; + pars->krMag = NAN; + pars->psfMag = NAN; + pars->peakMag = NAN; + + return pars; +} Index: trunk/psModules/src/objects/pmSourceExtendedPars.h =================================================================== --- trunk/psModules/src/objects/pmSourceExtendedPars.h (revision 32046) +++ trunk/psModules/src/objects/pmSourceExtendedPars.h (revision 32347) @@ -67,4 +67,19 @@ } pmSourceExtendedPars; +// additional measurements related to the model fits +typedef struct { + float Mxx; + float Mxy; + float Myy; + + float Mrf; + float Mrh; + + float apMag; + float krMag; + float psfMag; + float peakMag; +} pmSourceExtFitPars; + pmSourceRadialFlux *pmSourceRadialFluxAlloc(); bool psMemCheckSourceRadialFlux(psPtr ptr); @@ -92,5 +107,5 @@ bool pmSourceRadialProfileSortPair(psVector *index, psVector *extra); - +pmSourceExtFitPars *pmSourceExtFitParsAlloc (void); /// @} Index: trunk/psModules/src/objects/pmSourceFitModel.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitModel.c (revision 32046) +++ trunk/psModules/src/objects/pmSourceFitModel.c (revision 32347) @@ -176,7 +176,9 @@ break; case PM_SOURCE_FIT_EXT: - // EXT model fits all params (except sky) - nParams = params->n - 1; + // EXT model fits all shape params and Io (not Xo, Yo, sky) + nParams = params->n - 3; psVectorInit (constraint->paramMask, 0); + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 1; + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 1; constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1; break; @@ -193,11 +195,13 @@ break; case PM_SOURCE_FIT_NO_INDEX: - // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params + // PSF model only fits Io, Sxx, Sxy, Syy psVectorInit (constraint->paramMask, 0); + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 1; + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 1; constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1; if (params->n == 7) { - nParams = params->n - 1; + nParams = params->n - 3; } else { - nParams = params->n - 2; + nParams = params->n - 4; constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1; } Index: trunk/psModules/src/objects/pmSourceFitPCM.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitPCM.c (revision 32046) +++ trunk/psModules/src/objects/pmSourceFitPCM.c (revision 32347) @@ -48,6 +48,10 @@ // convolved model image. +# define TIMING 0 + bool pmSourceFitPCM (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { + if (TIMING) { psTimerStart ("pmSourceFitPCM"); } + psVector *params = pcm->modelConv->params; psVector *dparams = pcm->modelConv->dparams; @@ -64,5 +68,10 @@ psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); + float t1, t2, t3, t4, t5; + if (TIMING) { t1 = psTimerMark ("pmSourceFitPCM"); } + bool fitStatus = pmPCM_MinimizeChisq (myMin, covar, params, source, pcm); + if (TIMING) { t2 = psTimerMark ("pmSourceFitPCM"); } + for (int i = 0; i < dparams->n; i++) { if ((pcm->constraint->paramMask != NULL) && pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) @@ -76,4 +85,5 @@ } psTrace ("psphot", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value); + if (TIMING) { t3 = psTimerMark ("pmSourceFitPCM"); } // renormalize output model image (generated by fitting process) @@ -97,4 +107,5 @@ pmSourceChisqUnsubtracted (source, pcm->modelConv, maskVal); } + if (TIMING) { t4 = psTimerMark ("pmSourceFitPCM"); } // set the model success or failure status @@ -114,4 +125,9 @@ source->mode |= PM_SOURCE_MODE_FITTED; // XXX is this needed? + if (TIMING) { t5 = psTimerMark ("pmSourceFitPCM"); } + + if (TIMING) { + fprintf (stderr, "nIter: %2d, npix: %5d, t1: %6.4f, t2: %6.4f, t3: %6.4f, t4: %6.4f, t5: %6.4f\n", myMin->iter, pcm->nPix, t1, t2, t3, t4, t5); + } psFree(myMin); @@ -121,4 +137,5 @@ } +// XXX deprecate this function or merge with the empirical model bool pmSourceModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal) { Index: trunk/psModules/src/objects/pmSourceIO_CMF.c.in =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF.c.in (revision 32046) +++ trunk/psModules/src/objects/pmSourceIO_CMF.c.in (revision 32347) @@ -55,5 +55,4 @@ // followed by a zero-size matrix, followed by the table data -// # define MODE @CMFMODE@ bool pmSourcesWrite_CMF_@CMFMODE@ (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) { @@ -171,7 +170,11 @@ @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kouter); + // XXX do not keep this long term, just a TEST: + // @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_PSF", PS_DATA_F32, "Kron Flux", moments.KronPSF); + // @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_PSF_SIG",PS_DATA_F32, "Kron Flux", moments.KronPSFErr); // Do NOT write these : not consistent with the definition of PS1_V3 in Ohana/src/libautocode/dev/cmf-ps1-v3.d // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_FLUX", PS_DATA_F32, "Kron Flux (in 1.0 R1)", moments.KronCore); // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_ERROR", PS_DATA_F32, "Kron Error (in 1.0 R1)", moments.KronCoreErr); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); @@ -222,5 +225,5 @@ // read in a readout from the fits file -psArray *pmSourcesRead_CMF_PS1_V3 (psFits *fits, psMetadata *header) +psArray *pmSourcesRead_CMF_@CMFMODE@ (psFits *fits, psMetadata *header) { PS_ASSERT_PTR_NON_NULL(fits, false); @@ -281,6 +284,11 @@ // XXX use these to determine PAR[PM_PAR_I0]? @ALL@ source->psfMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG"); - @ALL@ source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); + @ALL@ source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); @ALL@ source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); + @=PS1_V3@ source->apMagRaw = psMetadataLookupF32 (&status, row, "AP_MAG_RAW"); + + // XXX use these to determine PAR[PM_PAR_I0] if they exist? + @=PS1_V3@ source->psfFlux = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX"); + @=PS1_V3@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG"); // XXX this scaling is incorrect: does not include the 2 \pi AREA factor @@ -288,5 +296,5 @@ @ALL@ dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; - pmPSF_AxesToModel (PAR, axes); + pmPSF_AxesToModel (PAR, axes, modelType); @ALL@ float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); @@ -353,5 +361,5 @@ } -bool pmSourcesWrite_CMF_PS1_V3_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) +bool pmSourcesWrite_CMF_@CMFMODE@_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) { bool status; @@ -541,5 +549,5 @@ // XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_V3_XFIT (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) +bool pmSourcesWrite_CMF_@CMFMODE@_XFIT (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) { @@ -590,4 +598,7 @@ assert (model); + // pmSourceExtFitPars *extPars = source->extFitPars->data[j]; + // assert (extPars); + // skip models which were not actually fitted if (model->flags & PM_MODEL_STATUS_BADARGS) continue; @@ -600,5 +611,8 @@ yErr = dPAR[PM_PAR_YPOS]; - axes = pmPSF_ModelToAxes (PAR, 20.0); + axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); + + float kronFlux = source->moments ? source->moments->KronFlux : NAN; + float kronMag = isfinite(kronFlux) ? -2.5*log10(kronFlux) : NAN; row = psMetadataAlloc (); @@ -612,4 +626,14 @@ psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG", 0, "EXT fit instrumental magnitude", model->mag); psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude", model->magErr); + + // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_XX", 0, "second moment in x", extPars->Mxx); + // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_XY", 0, "second moment in x,y", extPars->Mxy); + // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_YY", 0, "second moment in y", extPars->Myy); + // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_R1", 0, "first radial moment", extPars->Mrf); + // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_RH", 0, "half radial moment", extPars->Mrh); + + psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_INST_MAG", 0, "PSF fit instrumental magnitude", source->psfMag); + psMetadataAddF32 (row, PS_LIST_TAIL, "AP_MAG", 0, "PSF-sized aperture magnitude", source->apMag); + psMetadataAddF32 (row, PS_LIST_TAIL, "KRON_MAG", 0, "Kron Mag", kronMag); psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS", 0, "number of model parameters", model->params->n); @@ -673,5 +697,5 @@ } -bool pmSourcesWrite_CMF_PS1_V3_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) +bool pmSourcesWrite_CMF_@CMFMODE@_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) { return true; Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c (revision 32046) +++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c (revision 32347) @@ -263,5 +263,5 @@ dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; - pmPSF_AxesToModel (PAR, axes); + pmPSF_AxesToModel (PAR, axes, modelType); float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); @@ -547,5 +547,5 @@ yErr = dPAR[PM_PAR_YPOS]; - axes = pmPSF_ModelToAxes (PAR, 20.0); + axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); row = psMetadataAlloc (); Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c (revision 32046) +++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c (revision 32347) @@ -281,5 +281,5 @@ dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; - pmPSF_AxesToModel (PAR, axes); + pmPSF_AxesToModel (PAR, axes, modelType); float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); @@ -572,5 +572,5 @@ yErr = dPAR[PM_PAR_YPOS]; - axes = pmPSF_ModelToAxes (PAR, 20.0); + axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); row = psMetadataAlloc (); Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 32046) +++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 32347) @@ -280,5 +280,5 @@ dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; - pmPSF_AxesToModel (PAR, axes); + pmPSF_AxesToModel (PAR, axes, modelType); float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); @@ -607,5 +607,5 @@ yErr = dPAR[PM_PAR_YPOS]; - axes = pmPSF_ModelToAxes (PAR, 20.0); + axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); row = psMetadataAlloc (); Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c (revision 32046) +++ (revision ) @@ -1,664 +1,0 @@ -/** @file pmSourceIO.c - * - * @author EAM, IfA - * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-18 02:44:19 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include - -#include "pmConfig.h" -#include "pmDetrendDB.h" - -#include "pmHDU.h" -#include "pmFPA.h" -#include "pmFPALevel.h" -#include "pmFPAview.h" -#include "pmFPAfile.h" - -#include "pmTrend2D.h" -#include "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmSpan.h" -#include "pmFootprintSpans.h" -#include "pmFootprint.h" -#include "pmPeaks.h" -#include "pmMoments.h" -#include "pmModelFuncs.h" -#include "pmModel.h" -#include "pmModelUtils.h" -#include "pmModelClass.h" -#include "pmSourceMasks.h" -#include "pmSourceExtendedPars.h" -#include "pmSourceDiffStats.h" -#include "pmSource.h" -#include "pmSourceFitModel.h" -#include "pmPSF.h" -#include "pmPSFtry.h" - -#include "pmSourceIO.h" -#include "pmSourceOutputs.h" - -// panstarrs-style FITS table output (header + table in 1st extension) -// this format consists of a header derived from the image header -// followed by a zero-size matrix, followed by the table data - -bool pmSourcesWrite_CMF_PS1_V1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) -{ - PS_ASSERT_PTR_NON_NULL(fits, false); - PS_ASSERT_PTR_NON_NULL(sources, false); - PS_ASSERT_PTR_NON_NULL(extname, false); - - psArray *table; - psMetadata *row; - - pmChip *chip = readout->parent->parent; - - // if the sequence is defined, write these in seq order; otherwise - // write them in S/N order: - if (sources->n > 0) { - pmSource *source = (pmSource *) sources->data[0]; - if (source->seq == -1) { - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - } else { - sources = psArraySort (sources, pmSourceSortBySeq); - } - } - - table = psArrayAllocEmpty (sources->n); - - short nImageOverlap; - float magOffset; - float zeroptErr; - float fwhmMajor; - float fwhmMinor; - pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); - - // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state - // by the time we call this function, all values should be assigned. let's use asserts to be sure in some cases. - for (int i = 0; i < sources->n; i++) { - pmSource *source = (pmSource *) sources->data[i]; - - // If source->seq is -1, source was generated in this analysis. If source->seq is - // not -1, source was read from elsewhere: in the latter case, preserve the source - // ID. source.seq is used instead of source.id since the latter is a const - // generated on Alloc, and would thus be wrong for read in sources. - if (source->seq == -1) { - source->seq = i; - } - - // set the 'best' values for various output fields: - pmSourceOutputs outputs; - pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); - - pmSourceOutputsMoments moments; - pmSourceOutputsSetMoments (&moments, source); - - row = psMetadataAlloc (); - psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F32, "PSF RA coordinate (degrees)", outputs.ra); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F32, "PSF DEC coordinate (degrees)", outputs.dec); - psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", outputs.posAngle); - psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfMagErr); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", outputs.calMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr); - psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); - psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); - - psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq); - psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); - - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor", source->pixWeightNotBad); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); - - psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); - - // XXX not sure how to get this : need to load Nimages with weight? - psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", nImageOverlap); - psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", 0); - - psArrayAdd (table, 100, row); - psFree (row); - - // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not - // subtracted - - // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image - // edge; 3) any pixels in the 3x3 peak region are masked; - } - - psMetadata *header = psMetadataCopy(NULL, tableHeader); - pmSourceMasksHeader(header); - - if (table->n == 0) { - if (!psFitsWriteBlank(fits, header, extname)) { - psError(psErrorCodeLast(), false, "Unable to write empty sources file."); - psFree(table); - psFree(header); - return false; - } - psFree(table); - psFree(header); - return true; - } - - psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); - if (!psFitsWriteTable(fits, header, table, extname)) { - psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); - psFree(table); - psFree(header); - return false; - } - psFree(table); - psFree(header); - - return true; -} - -// read in a readout from the fits file -psArray *pmSourcesRead_CMF_PS1_V1 (psFits *fits, psMetadata *header) -{ - PS_ASSERT_PTR_NON_NULL(fits, false); - PS_ASSERT_PTR_NON_NULL(header, false); - - bool status; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - - // define PSF model type - // XXX need to carry the extra model parameters - int modelType = pmModelClassGetType ("PS_MODEL_GAUSS"); - - char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); - if (PSF_NAME != NULL) { - modelType = pmModelClassGetType (PSF_NAME); - } - assert (modelType > -1); - - // We get the size of the table, and allocate the array of sources first because the table - // is large and ephemeral --- when the table gets blown away, whatever is allocated after - // the table is read blocks the free. In fact, it's better to read the table row by row. - long numSources = psFitsTableSize(fits); // Number of sources in table - psArray *sources = psArrayAlloc(numSources); // Array of sources, to return - - // convert the table to the pmSource entriesa - for (int i = 0; i < numSources; i++) { - psMetadata *row = psFitsReadTableRow(fits, i); // Table row - - pmSource *source = pmSourceAlloc (); - pmModel *model = pmModelAlloc (modelType); - source->modelPSF = model; - source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags - - // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - - source->seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); - PAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF"); - PAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF"); - dPAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF_SIG"); - dPAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF_SIG"); - axes.major = psMetadataLookupF32 (&status, row, "PSF_MAJOR"); - axes.minor = psMetadataLookupF32 (&status, row, "PSF_MINOR"); - axes.theta = psMetadataLookupF32 (&status, row, "PSF_THETA"); - - PAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY"); - dPAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY_SIGMA"); - source->sky = PAR[PM_PAR_SKY]; - source->skyErr = dPAR[PM_PAR_SKY]; - - // XXX use these to determine PAR[PM_PAR_I0]? - source->psfMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG"); - source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); - source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); - - // XXX this scaling is incorrect: does not include the 2 \pi AREA factor - PAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN; - dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; - - pmPSF_AxesToModel (PAR, axes); - - float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); - float peakFlux = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN; - - // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) - source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); - source->peak->rawFlux = peakFlux; - source->peak->smoothFlux = peakFlux; - source->peak->xf = PAR[PM_PAR_XPOS]; // more accurate position - source->peak->yf = PAR[PM_PAR_YPOS]; // more accurate position - source->peak->dx = dPAR[PM_PAR_XPOS]; - source->peak->dy = dPAR[PM_PAR_YPOS]; - - source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF"); - source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); - source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); - source->apRadius = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); - - // note that some older versions used PSF_PROBABILITY: this was not well defined. - model->chisq = psMetadataLookupF32 (&status, row, "PSF_CHISQ"); - model->nDOF = psMetadataLookupS32 (&status, row, "PSF_NDOF"); - model->nPix = psMetadataLookupS32 (&status, row, "PSF_NPIX"); - - source->moments = pmMomentsAlloc (); - source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf - source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf - - source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX"); - source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY"); - source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY"); - - source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); - assert (status); - - sources->data[i] = source; - psFree(row); - } - - return sources; -} - -// XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_V1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) -{ - - bool status; - psArray *table; - psMetadata *row; - psF32 *PAR, *dPAR; - psF32 xPos, yPos; - psF32 xErr, yErr; - - // create a header to hold the output data - psMetadata *outhead = psMetadataAlloc (); - - // write the links to the image header - psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); - - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - - table = psArrayAllocEmpty (sources->n); - - // which extended source analyses should we perform? - // bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); - // bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); - // bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); - // bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); - - psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); - psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); - assert (radialBinsLower->n == radialBinsUpper->n); - - // we write out all sources, regardless of quality. the source flags tell us the state - for (int i = 0; i < sources->n; i++) { - // skip source if it is not a ext sourc - // XXX we have two places that extended source parameters are measured: - // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models, - // psphotFitEXT, which does the simple extended source model fit (not psf-convolved) - // should we require both? - - pmSource *source = sources->data[i]; - - // skip sources without measurements - if (source->extpars == NULL) continue; - - // we require a PSF model fit (ignore the real crud) - pmModel *model = source->modelPSF; - if (model == NULL) continue; - - // XXX I need to split the extended models from the extended aperture measurements - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - - row = psMetadataAlloc (); - - // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) - psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_EXT", PS_DATA_F32, "EXT model x coordinate", xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT", PS_DATA_F32, "EXT model y coordinate", yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_F32, "Sigma in EXT x coordinate", xErr); - psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); - -# if (0) - // Petrosian measurements - // XXX insert header data: petrosian ref radius, flux ratio - if (doPetrosian) { - pmSourcePetrosianValues *petrosian = source->extpars->petrosian; - if (petrosian) { - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", petrosian->mag); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", petrosian->rad); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", petrosian->radErr); - } else { - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", NAN); - } - } - - // Kron measurements - if (doKron) { - pmSourceKronValues *kron = source->extpars->kron; - if (kron) { - psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", kron->mag); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", kron->magErr); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", kron->rad); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", kron->radErr); - } else { - psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", NAN); - } - } - - // Isophot measurements - // XXX insert header data: isophotal level - if (doIsophotal) { - pmSourceIsophotalValues *isophot = source->extpars->isophot; - if (isophot) { - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", isophot->mag); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", isophot->rad); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", isophot->radErr); - } else { - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", NAN); - } - } - - // Flux Annuli - if (doAnnuli) { - pmSourceAnnuli *annuli = source->extpars->annuli; - if (annuli) { - psVector *fluxVal = annuli->flux; - psVector *fluxErr = annuli->fluxErr; - psVector *fluxVar = annuli->fluxVar; - - for (int j = 0; j < fluxVal->n; j++) { - char name[32]; - sprintf (name, "FLUX_VAL_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]); - sprintf (name, "FLUX_ERR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]); - sprintf (name, "FLUX_VAR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]); - } - } else { - for (int j = 0; j < radialBinsLower->n; j++) { - char name[32]; - sprintf (name, "FLUX_VAL_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN); - sprintf (name, "FLUX_ERR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN); - sprintf (name, "FLUX_VAR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN); - } - } - } - -# endif - psArrayAdd (table, 100, row); - psFree (row); - } - - if (table->n == 0) { - if (!psFitsWriteBlank (fits, outhead, extname)) { - psError(psErrorCodeLast(), false, "Unable to write empty sources."); - psFree(outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - return true; - } - - psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); - if (!psFitsWriteTable (fits, outhead, table, extname)) { - psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); - psFree (outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - - return true; -} - -// XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_V1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) -{ - - psArray *table; - psMetadata *row; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - psF32 xPos, yPos; - psF32 xErr, yErr; - char name[64]; - - // create a header to hold the output data - psMetadata *outhead = psMetadataAlloc (); - - // write the links to the image header - psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); - - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - - // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams - int nParamMax = 0; - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source->modelFits == NULL) continue; - for (int j = 0; j < source->modelFits->n; j++) { - pmModel *model = source->modelFits->data[j]; - assert (model); - nParamMax = PS_MAX (nParamMax, model->params->n); - } - } - - table = psArrayAllocEmpty (sources->n); - - // we write out all sources, regardless of quality. the source flags tell us the state - for (int i = 0; i < sources->n; i++) { - - pmSource *source = sources->data[i]; - - // XXX if no model fits are saved, write out modelEXT? - if (source->modelFits == NULL) continue; - - // We have multiple sources : need to flag the one used to subtract the light (the 'best' model) - for (int j = 0; j < source->modelFits->n; j++) { - - // choose the convolved EXT model, if available, otherwise the simple one - pmModel *model = source->modelFits->data[j]; - assert (model); - - // skip models which were not actually fitted - if (model->flags & PM_MODEL_STATUS_BADARGS) continue; - - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - - axes = pmPSF_ModelToAxes (PAR, 20.0); - - row = psMetadataAlloc (); - - // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) - psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); - psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT", 0, "EXT model x coordinate", xPos); - psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT", 0, "EXT model y coordinate", yPos); - psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG", 0, "Sigma in EXT x coordinate", xErr); - psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG", 0, "Sigma in EXT y coordinate", yErr); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG", 0, "EXT fit instrumental magnitude", model->mag); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude", model->magErr); - - psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS", 0, "number of model parameters", model->params->n); - psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "name of model", pmModelClassGetName (model->type)); - - // XXX these should be major and minor, not 'x' and 'y' - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width in x coordinate", axes.major); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width in y coordinate", axes.minor); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", axes.theta); - - // write out the other generic parameters - for (int k = 0; k < nParamMax; k++) { - if (k == PM_PAR_I0) continue; - if (k == PM_PAR_SKY) continue; - if (k == PM_PAR_XPOS) continue; - if (k == PM_PAR_YPOS) continue; - if (k == PM_PAR_SXX) continue; - if (k == PM_PAR_SXY) continue; - if (k == PM_PAR_SYY) continue; - - snprintf (name, 64, "EXT_PAR_%02d", k); - - if (k < model->params->n) { - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]); - } else { - psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN); - } - } - - // XXX other parameters which may be set. - // XXX flag / value to define the model - // XXX write out the model type, fit status flags - - psArrayAdd (table, 100, row); - psFree (row); - } - } - - if (table->n == 0) { - if (!psFitsWriteBlank (fits, outhead, extname)) { - psError(psErrorCodeLast(), false, "Unable to write empty sources."); - psFree(outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - return true; - } - - psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); - if (!psFitsWriteTable (fits, outhead, table, extname)) { - psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); - psFree (outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - return true; -} - -bool pmSourceLocalAstrometry (psSphere *ptSky, float *posAngle, float *pltScale, pmChip *chip, float xPos, float yPos) { - - pmFPA *fpa = chip->parent; - - if (!chip->toFPA) goto escape; - if (!fpa->toTPA) goto escape; - if (!fpa->toSky) goto escape; - - // generate RA,DEC - psPlane ptCH, ptFP, ptTP_o, ptTP_x, ptTP_y; - - // calculate the astrometry for the coordinate of interest - ptCH.x = xPos; - ptCH.y = yPos; - psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); - psPlaneTransformApply (&ptTP_o, fpa->toTPA, &ptFP); - psDeproject (ptSky, &ptTP_o, fpa->toSky); - - // calculate the astrometry for the coordinate + 1pix in X - ptCH.x = xPos + 1.0; - ptCH.y = yPos; - psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); - psPlaneTransformApply (&ptTP_x, fpa->toTPA, &ptFP); - - // calculate the astrometry for the coordinate + 1pix in Y - ptCH.x = xPos; - ptCH.y = yPos + 1.0; - psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); - psPlaneTransformApply (&ptTP_y, fpa->toTPA, &ptFP); - - // the resulting Tangent Plane coordinates are in TP pixels; convert to local Tangent Plane - // degrees - - float dTPx_dCHx = fpa->toSky->Xs * (ptTP_x.x - ptTP_o.x); - float dTPy_dCHx = fpa->toSky->Ys * (ptTP_x.y - ptTP_o.y); - - float dTPx_dCHy = fpa->toSky->Xs * (ptTP_y.x - ptTP_o.x); - float dTPy_dCHy = fpa->toSky->Ys * (ptTP_y.y - ptTP_o.y); - - float pltScale_x = hypot(dTPx_dCHx, dTPy_dCHx); - float pltScale_y = hypot(dTPx_dCHy, dTPy_dCHy); - *pltScale = 0.5*(pltScale_x + pltScale_y); - - float posAngle_x = atan2 (+dTPy_dCHx, +dTPx_dCHx); - float posAngle_y = atan2 (-dTPy_dCHy, +dTPx_dCHy); - *posAngle = 0.5*(posAngle_x + posAngle_y); - - return true; - -escape: - // no astrometry calibration, give up - ptSky->r = NAN; - ptSky->d = NAN; - *posAngle = NAN; - *pltScale = NAN; - - return false; -} - -bool pmSourcesWrite_CMF_PS1_V1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) -{ - return true; -} Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 32046) +++ (revision ) @@ -1,669 +1,0 @@ -/** @file pmSourceIO.c - * - * @author EAM, IfA - * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-18 02:44:19 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include - -#include "pmConfig.h" -#include "pmDetrendDB.h" - -#include "pmHDU.h" -#include "pmFPA.h" -#include "pmFPALevel.h" -#include "pmFPAview.h" -#include "pmFPAfile.h" - -#include "pmTrend2D.h" -#include "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmSpan.h" -#include "pmFootprintSpans.h" -#include "pmFootprint.h" -#include "pmPeaks.h" -#include "pmMoments.h" -#include "pmModelFuncs.h" -#include "pmModel.h" -#include "pmModelUtils.h" -#include "pmModelClass.h" -#include "pmSourceMasks.h" -#include "pmSourceExtendedPars.h" -#include "pmSourceDiffStats.h" -#include "pmSource.h" -#include "pmSourceFitModel.h" -#include "pmPSF.h" -#include "pmPSFtry.h" - -#include "pmSourceIO.h" -#include "pmSourceOutputs.h" - -// panstarrs-style FITS table output (header + table in 1st extension) -// this format consists of a header derived from the image header -// followed by a zero-size matrix, followed by the table data - -bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) -{ - PS_ASSERT_PTR_NON_NULL(fits, false); - PS_ASSERT_PTR_NON_NULL(sources, false); - PS_ASSERT_PTR_NON_NULL(extname, false); - - psArray *table; - psMetadata *row; - - pmChip *chip = readout->parent->parent; - - // if the sequence is defined, write these in seq order; otherwise - // write them in S/N order: - if (sources->n > 0) { - pmSource *source = (pmSource *) sources->data[0]; - if (source->seq == -1) { - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - } else { - sources = psArraySort (sources, pmSourceSortBySeq); - } - } - - table = psArrayAllocEmpty (sources->n); - - short nImageOverlap; - float magOffset; - float zeroptErr; - float fwhmMajor; - float fwhmMinor; - pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); - - // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state - // by the time we call this function, all values should be assigned. let's use asserts to be sure in some cases. - for (int i = 0; i < sources->n; i++) { - pmSource *source = (pmSource *) sources->data[i]; - - // If source->seq is -1, source was generated in this analysis. If source->seq is - // not -1, source was read from elsewhere: in the latter case, preserve the source - // ID. source.seq is used instead of source.id since the latter is a const - // generated on Alloc, and would thus be wrong for read in sources. - if (source->seq == -1) { - source->seq = i; - } - - // set the 'best' values for various output fields: - pmSourceOutputs outputs; - pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); - - pmSourceOutputsMoments moments; - pmSourceOutputsSetMoments (&moments, source); - - row = psMetadataAlloc (); - psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr); // XXX this is only measured for non-linear fits - - psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", outputs.posAngle); - psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfMagErr); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", outputs.calMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr); - psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", outputs.ra); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", outputs.dec); - psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); - psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); - - psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq); - psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); - - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor", source->pixWeightNotBad); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); - - psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); - - // XXX not sure how to get this : need to load Nimages with weight? - psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", nImageOverlap); - psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", 0); - - psArrayAdd (table, 100, row); - psFree (row); - - // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not - // subtracted - - // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image - // edge; 3) any pixels in the 3x3 peak region are masked; - } - - // XXX why do we make a copy here to be supplemented with the masks? why not do this in the calling function? - psMetadata *header = psMetadataCopy(NULL, tableHeader); - pmSourceMasksHeader(header); - - if (table->n == 0) { - if (!psFitsWriteBlank(fits, header, extname)) { - psError(psErrorCodeLast(), false, "Unable to write blank sources file."); - psFree(table); - psFree(header); - return false; - } - psFree(table); - psFree(header); - return true; - } - - psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); - if (!psFitsWriteTable(fits, header, table, extname)) { - psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); - psFree(table); - psFree(header); - return false; - } - psFree(table); - psFree(header); - - return true; -} - -// read in a readout from the fits file -psArray *pmSourcesRead_CMF_PS1_V2 (psFits *fits, psMetadata *header) -{ - PS_ASSERT_PTR_NON_NULL(fits, false); - PS_ASSERT_PTR_NON_NULL(header, false); - - bool status; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - - // define PSF model type - // XXX need to carry the extra model parameters - int modelType = pmModelClassGetType ("PS_MODEL_GAUSS"); - - char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); - if (PSF_NAME != NULL) { - modelType = pmModelClassGetType (PSF_NAME); - } - assert (modelType > -1); - - // We get the size of the table, and allocate the array of sources first because the table - // is large and ephemeral --- when the table gets blown away, whatever is allocated after - // the table is read blocks the free. In fact, it's better to read the table row by row. - long numSources = psFitsTableSize(fits); // Number of sources in table - psArray *sources = psArrayAlloc(numSources); // Array of sources, to return - - // convert the table to the pmSource entriesa - for (int i = 0; i < numSources; i++) { - psMetadata *row = psFitsReadTableRow(fits, i); // Table row - if (!row) { - psError(psErrorCodeLast(), false, "Unable to read row %d of sources", i); - psFree(sources); - return NULL; - } - - pmSource *source = pmSourceAlloc (); - pmModel *model = pmModelAlloc (modelType); - source->modelPSF = model; - source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags - - // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - - source->seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); - PAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF"); - PAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF"); - dPAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF_SIG"); - dPAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF_SIG"); - axes.major = psMetadataLookupF32 (&status, row, "PSF_MAJOR"); - axes.minor = psMetadataLookupF32 (&status, row, "PSF_MINOR"); - axes.theta = psMetadataLookupF32 (&status, row, "PSF_THETA"); - - PAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY"); - dPAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY_SIGMA"); - source->sky = PAR[PM_PAR_SKY]; - source->skyErr = dPAR[PM_PAR_SKY]; - - // XXX use these to determine PAR[PM_PAR_I0]? - source->psfMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG"); - source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); - source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); - - // XXX this scaling is incorrect: does not include the 2 \pi AREA factor - PAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN; - dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; - - pmPSF_AxesToModel (PAR, axes); - - float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); - float peakFlux = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN; - - // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) - source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); - source->peak->rawFlux = peakFlux; - source->peak->smoothFlux = peakFlux; - source->peak->xf = PAR[PM_PAR_XPOS]; // more accurate position - source->peak->yf = PAR[PM_PAR_YPOS]; // more accurate position - source->peak->dx = dPAR[PM_PAR_XPOS]; - source->peak->dy = dPAR[PM_PAR_YPOS]; - - source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF"); - source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); - source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); - source->apRadius = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); - - // note that some older versions used PSF_PROBABILITY: this was not well defined. - model->chisq = psMetadataLookupF32 (&status, row, "PSF_CHISQ"); - model->nDOF = psMetadataLookupS32 (&status, row, "PSF_NDOF"); - model->nPix = psMetadataLookupS32 (&status, row, "PSF_NPIX"); - - source->moments = pmMomentsAlloc (); - source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf - source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf - - source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX"); - source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY"); - source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY"); - - source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); - assert (status); - - sources->data[i] = source; - psFree(row); - } - - return sources; -} - -bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) -{ - - bool status; - psArray *table; - psMetadata *row; - psF32 *PAR, *dPAR; - psF32 xPos, yPos; - psF32 xErr, yErr; - int nRow = -1; - - // create a header to hold the output data - psMetadata *outhead = psMetadataAlloc (); - - // write the links to the image header - psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); - - pmChip *chip = readout->parent->parent; - pmFPA *fpa = chip->parent; - - // zero point corrections - bool status1 = false; - bool status2 = false; - float magOffset = 0.0; - float exptime = psMetadataLookupF32(&status1, fpa->concepts, "FPA.EXPOSURE"); - float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); - if (!isfinite(zeropt)) { - zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); - } - if (status1 && status2 && (exptime > 0.0)) { - magOffset = zeropt + 2.5*log10(exptime); - } - - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - - table = psArrayAllocEmpty (sources->n); - - // which extended source analyses should we perform? - bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); - bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); - // bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); - // bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); - - psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); - psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); - psAssert (radMin->n == radMax->n, "inconsistent annular bins"); - - // int nRadialBins = 0; - // if (doAnnuli) { - // // get the max count of radial bins - // for (int i = 0; i < sources->n; i++) { - // pmSource *source = sources->data[i]; - // if (!source->extpars) continue; - // if (!source->extpars->radProfile ) continue; - // if (!source->extpars->radProfile->binSB) continue; - // nRadialBins = PS_MAX(nRadialBins, source->extpars->radProfile->binSB->n); - // } - // } - - // we write out all sources, regardless of quality. the source flags tell us the state - for (int i = 0; i < sources->n; i++) { - // skip source if it is not a ext sourc - // XXX we have two places that extended source parameters are measured: - // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models, - // psphotFitEXT, which does the simple extended source model fit (not psf-convolved) - // should we require both? - - pmSource *source = sources->data[i]; - - // skip sources without measurements - if (source->extpars == NULL) continue; - - // we require a PSF model fit (ignore the real crud) - pmModel *model = source->modelPSF; - if (model == NULL) continue; - - // XXX I need to split the extended models from the extended aperture measurements - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - - row = psMetadataAlloc (); - - // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) - psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_EXT", PS_DATA_F32, "EXT model x coordinate", xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT", PS_DATA_F32, "EXT model y coordinate", yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_F32, "Sigma in EXT x coordinate", xErr); - psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); - - float AxialRatio = NAN; - float AxialTheta = NAN; - pmSourceExtendedPars *extpars = source->extpars; - if (extpars) { - AxialRatio = extpars->axes.minor / extpars->axes.major; - AxialTheta = extpars->axes.theta; - } - psMetadataAdd (row, PS_LIST_TAIL, "F25_ARATIO", PS_DATA_F32, "Axial Ratio of radial profile", AxialRatio); - psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA", PS_DATA_F32, "Angle of radial profile ellipse", AxialTheta); - - // Petrosian measurements - // XXX insert header data: petrosian ref radius, flux ratio - // XXX check flags to see if Pet was measured - if (doPetrosian) { - pmSourceExtendedPars *extpars = source->extpars; - if (extpars) { - // XXX note that this mag is either calibrated or instrumental depending on existence of zero point - float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point - float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", mag); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", magErr); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius (pix)", extpars->petrosianRadius); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error (pix)", extpars->petrosianRadiusErr); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50", PS_DATA_F32, "Petrosian R50 (pix)", extpars->petrosianR50); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)", extpars->petrosianR50Err); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90", PS_DATA_F32, "Petrosian R90 (pix)", extpars->petrosianR90); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)", extpars->petrosianR90Err); - } else { - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50", PS_DATA_F32, "Petrosian R50 (pix)", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)",NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90", PS_DATA_F32, "Petrosian R90 (pix)", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)",NAN); - } - } - -# if (0) - // Kron measurements - if (doKron) { - pmSourceKronValues *kron = source->extpars->kron; - if (kron) { - psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", kron->mag); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", kron->magErr); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", kron->rad); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", kron->radErr); - } else { - psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", NAN); - } - } - - // Isophot measurements - // XXX insert header data: isophotal level - if (doIsophotal) { - pmSourceIsophotalValues *isophot = source->extpars->isophot; - if (isophot) { - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", isophot->mag); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", isophot->rad); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", isophot->radErr); - } else { - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", NAN); - } - } -# endif - - // Flux Annuli (if we have extended source measurements, we have these. only optionally save them) - if (doAnnuli) { - psVector *radSB = psVectorAlloc(radMin->n, PS_TYPE_F32); - psVector *radFlux = psVectorAlloc(radMin->n, PS_TYPE_F32); - psVector *radFill = psVectorAlloc(radMin->n, PS_TYPE_F32); - psVectorInit (radSB, NAN); - psVectorInit (radFlux, NAN); - psVectorInit (radFill, NAN); - if (!source->extpars) goto empty_annuli; - if (!source->extpars->radProfile) goto empty_annuli; - if (!source->extpars->radProfile->binSB) goto empty_annuli; - psAssert (source->extpars->radProfile->binSum, "programming error"); - psAssert (source->extpars->radProfile->binFill, "programming error"); - psAssert (source->extpars->radProfile->binSB->n <= radFlux->n, "inconsistent vector lengths"); - psAssert (source->extpars->radProfile->binSum->n <= radFlux->n, "inconsistent vector lengths"); - psAssert (source->extpars->radProfile->binFill->n <= radFlux->n, "inconsistent vector lengths"); - - // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux - for (int j = 0; j < source->extpars->radProfile->binSB->n; j++) { - radSB->data.F32[j] = source->extpars->radProfile->binSB->data.F32[j]; - radFlux->data.F32[j] = source->extpars->radProfile->binSum->data.F32[j]; - radFill->data.F32[j] = source->extpars->radProfile->binFill->data.F32[j]; - } - - empty_annuli: - psMetadataAdd (row, PS_LIST_TAIL, "PROF_SB", PS_DATA_VECTOR, "mean surface brightness annuli", radSB); - psMetadataAdd (row, PS_LIST_TAIL, "PROF_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux); - psMetadataAdd (row, PS_LIST_TAIL, "PROF_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); - psFree (radSB); - psFree (radFlux); - psFree (radFill); - } - if (nRow < 0) { - nRow = row->list->n; - } else { - psAssert (nRow == row->list->n, "inconsistent row lengths"); - } - psArrayAdd (table, 100, row); - psFree (row); - } - - if (table->n == 0) { - if (!psFitsWriteBlank (fits, outhead, extname)) { - psError(psErrorCodeLast(), false, "Unable to write empty sources file."); - psFree(outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - return true; - } - - psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); - if (!psFitsWriteTable (fits, outhead, table, extname)) { - psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); - psFree (outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - - return true; -} - -// XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_V2_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) -{ - - psArray *table; - psMetadata *row; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - psF32 xPos, yPos; - psF32 xErr, yErr; - char name[64]; - - // create a header to hold the output data - psMetadata *outhead = psMetadataAlloc (); - - // write the links to the image header - psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); - - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - - // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams - int nParamMax = 0; - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source->modelFits == NULL) continue; - for (int j = 0; j < source->modelFits->n; j++) { - pmModel *model = source->modelFits->data[j]; - assert (model); - nParamMax = PS_MAX (nParamMax, model->params->n); - } - } - - table = psArrayAllocEmpty (sources->n); - - // we write out all sources, regardless of quality. the source flags tell us the state - for (int i = 0; i < sources->n; i++) { - - pmSource *source = sources->data[i]; - - // XXX if no model fits are saved, write out modelEXT? - if (source->modelFits == NULL) continue; - - // We have multiple sources : need to flag the one used to subtract the light (the 'best' model) - for (int j = 0; j < source->modelFits->n; j++) { - - // choose the convolved EXT model, if available, otherwise the simple one - pmModel *model = source->modelFits->data[j]; - assert (model); - - // skip models which were not actually fitted - if (model->flags & PM_MODEL_STATUS_BADARGS) continue; - - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - - axes = pmPSF_ModelToAxes (PAR, 20.0); - - row = psMetadataAlloc (); - - // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) - psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); - psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT", 0, "EXT model x coordinate", xPos); - psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT", 0, "EXT model y coordinate", yPos); - psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG", 0, "Sigma in EXT x coordinate", xErr); - psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG", 0, "Sigma in EXT y coordinate", yErr); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG", 0, "EXT fit instrumental magnitude", model->mag); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude", model->magErr); - - psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS", 0, "number of model parameters", model->params->n); - psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "name of model", pmModelClassGetName (model->type)); - - // XXX these should be major and minor, not 'x' and 'y' - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width in x coordinate", axes.major); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width in y coordinate", axes.minor); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", axes.theta); - - // write out the other generic parameters - for (int k = 0; k < nParamMax; k++) { - if (k == PM_PAR_I0) continue; - if (k == PM_PAR_SKY) continue; - if (k == PM_PAR_XPOS) continue; - if (k == PM_PAR_YPOS) continue; - if (k == PM_PAR_SXX) continue; - if (k == PM_PAR_SXY) continue; - if (k == PM_PAR_SYY) continue; - - snprintf (name, 64, "EXT_PAR_%02d", k); - - if (k < model->params->n) { - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]); - } else { - psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN); - } - } - - // XXX other parameters which may be set. - // XXX flag / value to define the model - // XXX write out the model type, fit status flags - - psArrayAdd (table, 100, row); - psFree (row); - } - } - - if (table->n == 0) { - if (!psFitsWriteBlank (fits, outhead, extname)) { - psError(psErrorCodeLast(), false, "Unable to write empty sources file."); - psFree(outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - return true; - } - - psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); - if (!psFitsWriteTable (fits, outhead, table, extname)) { - psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); - psFree (outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - return true; -} - -bool pmSourcesWrite_CMF_PS1_V2_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) -{ - return true; -} Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c (revision 32046) +++ (revision ) @@ -1,680 +1,0 @@ -/** @file pmSourceIO.c - * - * @author EAM, IfA - * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-18 02:44:19 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include - -#include "pmConfig.h" -#include "pmDetrendDB.h" - -#include "pmHDU.h" -#include "pmFPA.h" -#include "pmFPALevel.h" -#include "pmFPAview.h" -#include "pmFPAfile.h" - -#include "pmTrend2D.h" -#include "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmSpan.h" -#include "pmFootprintSpans.h" -#include "pmFootprint.h" -#include "pmPeaks.h" -#include "pmMoments.h" -#include "pmModelFuncs.h" -#include "pmModel.h" -#include "pmModelUtils.h" -#include "pmModelClass.h" -#include "pmSourceMasks.h" -#include "pmSourceExtendedPars.h" -#include "pmSourceDiffStats.h" -#include "pmSource.h" -#include "pmSourceFitModel.h" -#include "pmPSF.h" -#include "pmPSFtry.h" - -#include "pmSourceIO.h" -#include "pmSourceOutputs.h" - -// panstarrs-style FITS table output (header + table in 1st extension) -// this format consists of a header derived from the image header -// followed by a zero-size matrix, followed by the table data - -bool pmSourcesWrite_CMF_PS1_V3 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) -{ - PS_ASSERT_PTR_NON_NULL(fits, false); - PS_ASSERT_PTR_NON_NULL(sources, false); - PS_ASSERT_PTR_NON_NULL(extname, false); - - psArray *table; - psMetadata *row; - - pmChip *chip = readout->parent->parent; - - // if the sequence is defined, write these in seq order; otherwise - // write them in S/N order: - if (sources->n > 0) { - pmSource *source = (pmSource *) sources->data[0]; - if (source->seq == -1) { - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - } else { - sources = psArraySort (sources, pmSourceSortBySeq); - } - } - - table = psArrayAllocEmpty (sources->n); - - short nImageOverlap; - float magOffset; - float zeroptErr; - float fwhmMajor; - float fwhmMinor; - pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); - - // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state - // by the time we call this function, all values should be assigned. let's use asserts to be sure in some cases. - for (int i = 0; i < sources->n; i++) { - pmSource *source = (pmSource *) sources->data[i]; - - // If source->seq is -1, source was generated in this analysis. If source->seq is - // not -1, source was read from elsewhere: in the latter case, preserve the source - // ID. source.seq is used instead of source.id since the latter is a const - // generated on Alloc, and would thus be wrong for read in sources. - if (source->seq == -1) { - source->seq = i; - } - - // set the 'best' values for various output fields: - pmSourceOutputs outputs; - pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); - - pmSourceOutputsMoments moments; - pmSourceOutputsSetMoments (&moments, source); - - row = psMetadataAlloc (); - psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr); - psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", outputs.posAngle); - psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfMagErr); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental flux (counts)", source->psfFlux); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux", source->psfFluxErr); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in reported aperture", source->apMagRaw); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); - - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", outputs.calMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr); - - // NOTE: RA & DEC (both double) need to be on an 8-byte boundary... - psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", outputs.ra); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", outputs.dec); - - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); - psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); - psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); - - psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq); - psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); - - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor (bad)", source->pixWeightNotBad); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT", PS_DATA_F32, "PSF coverage/quality factor (poor)", source->pixWeightNotPoor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", moments.M_c3); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", moments.M_s3); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", moments.M_c4); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", moments.M_s4); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", moments.Mrf); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", moments.Mrh); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Krf); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", moments.dKrf); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kinner); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kouter); - - // Do NOT write these : not consistent with the definition of PS1_V3 in Ohana/src/libautocode/dev/cmf-ps1-v3.d - // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_FLUX", PS_DATA_F32, "Kron Flux (in 1.0 R1)", moments.KronCore); - // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_ERROR", PS_DATA_F32, "Kron Error (in 1.0 R1)", moments.KronCoreErr); - psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); - psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); - psMetadataAdd (row, PS_LIST_TAIL, "PADDING2", PS_DATA_S32, "more padding", 0); - - // XXX not sure how to get this : need to load Nimages with weight? - psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", nImageOverlap); - psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", 0); - - psArrayAdd (table, 100, row); - psFree (row); - - // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not - // subtracted - - // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image - // edge; 3) any pixels in the 3x3 peak region are masked; - } - - // XXX why do we make a copy here to be supplemented with the masks? why not do this in the calling function? - psMetadata *header = psMetadataCopy(NULL, tableHeader); - pmSourceMasksHeader(header); - - if (table->n == 0) { - if (!psFitsWriteBlank(fits, header, extname)) { - psError(psErrorCodeLast(), false, "Unable to write blank sources file."); - psFree(table); - psFree(header); - return false; - } - psFree(table); - psFree(header); - return true; - } - - psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); - if (!psFitsWriteTable(fits, header, table, extname)) { - psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); - psFree(table); - psFree(header); - return false; - } - psFree(table); - psFree(header); - - return true; -} - -// read in a readout from the fits file -psArray *pmSourcesRead_CMF_PS1_V3 (psFits *fits, psMetadata *header) -{ - PS_ASSERT_PTR_NON_NULL(fits, false); - PS_ASSERT_PTR_NON_NULL(header, false); - - bool status; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - - // define PSF model type - // XXX need to carry the extra model parameters - int modelType = pmModelClassGetType ("PS_MODEL_GAUSS"); - - char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); - if (PSF_NAME != NULL) { - modelType = pmModelClassGetType (PSF_NAME); - } - assert (modelType > -1); - - // We get the size of the table, and allocate the array of sources first because the table - // is large and ephemeral --- when the table gets blown away, whatever is allocated after - // the table is read blocks the free. In fact, it's better to read the table row by row. - long numSources = psFitsTableSize(fits); // Number of sources in table - psArray *sources = psArrayAlloc(numSources); // Array of sources, to return - - // convert the table to the pmSource entriesa - for (int i = 0; i < numSources; i++) { - psMetadata *row = psFitsReadTableRow(fits, i); // Table row - if (!row) { - psError(psErrorCodeLast(), false, "Unable to read row %d of sources", i); - psFree(sources); - return NULL; - } - - pmSource *source = pmSourceAlloc (); - pmModel *model = pmModelAlloc (modelType); - source->modelPSF = model; - source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags - - // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - - source->seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); - PAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF"); - PAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF"); - dPAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF_SIG"); - dPAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF_SIG"); - axes.major = psMetadataLookupF32 (&status, row, "PSF_MAJOR"); - axes.minor = psMetadataLookupF32 (&status, row, "PSF_MINOR"); - axes.theta = psMetadataLookupF32 (&status, row, "PSF_THETA"); - - PAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY"); - dPAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY_SIGMA"); - source->sky = PAR[PM_PAR_SKY]; - source->skyErr = dPAR[PM_PAR_SKY]; - - source->psfMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG"); - source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); - source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); - source->apMagRaw = psMetadataLookupF32 (&status, row, "AP_MAG_RAW"); - - // XXX use these to determine PAR[PM_PAR_I0] if they exist? - source->psfFlux = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX"); - source->psfFluxErr = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG"); - - // XXX this scaling is incorrect: does not include the 2 \pi AREA factor - PAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN; - dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; - - pmPSF_AxesToModel (PAR, axes); - - float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); - float peakFlux = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN; - - // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) - source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); - source->peak->rawFlux = peakFlux; - source->peak->smoothFlux = peakFlux; - source->peak->xf = PAR[PM_PAR_XPOS]; // more accurate position - source->peak->yf = PAR[PM_PAR_YPOS]; // more accurate position - source->peak->dx = dPAR[PM_PAR_XPOS]; - source->peak->dy = dPAR[PM_PAR_YPOS]; - - // we no longer sort by S/N, only flux - // if (isfinite (source->psfMagErr) && (source->psfMagErr > 0.0)) { - // source->peak->SN = 1.0 / source->psfMagErr; - // } else { - // source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N - // } - - source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF"); - source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT"); - source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); - source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); - source->apRadius = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); - - // note that some older versions used PSF_PROBABILITY: this was not well defined. - model->chisq = psMetadataLookupF32 (&status, row, "PSF_CHISQ"); - model->nDOF = psMetadataLookupS32 (&status, row, "PSF_NDOF"); - model->nPix = psMetadataLookupS32 (&status, row, "PSF_NPIX"); - - source->moments = pmMomentsAlloc (); - source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf - source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf - - source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX"); - source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY"); - source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY"); - - source->moments->Mrf = psMetadataLookupF32 (&status, row, "MOMENTS_R1"); - source->moments->Mrh = psMetadataLookupF32 (&status, row, "MOMENTS_RH"); - source->moments->KronFlux = psMetadataLookupF32 (&status, row, "KRON_FLUX"); - source->moments->KronFluxErr = psMetadataLookupF32 (&status, row, "KRON_FLUX_ERR"); - - source->moments->KronFinner = psMetadataLookupF32 (&status, row, "KRON_FLUX_INNER"); - source->moments->KronFouter = psMetadataLookupF32 (&status, row, "KRON_FLUX_OUTER"); - - // XXX we do not save all of the 3rd and 4th moment parameters. when we load in data, - // we are storing enough information so the output will be consistent with the input - source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C"); - source->moments->Mxxy = 0.0; - source->moments->Mxyy = 0.0; - source->moments->Myyy = -1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3S"); - - source->moments->Mxxxx = +1.00 * psMetadataLookupF32 (&status, row, "MOMENTS_M4C"); - source->moments->Mxxxy = 0.0; - source->moments->Mxxyy = 0.0; - source->moments->Mxyyy = -0.25 * psMetadataLookupF32 (&status, row, "MOMENTS_M4S"); - source->moments->Myyyy = 0.0; - - source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); - source->mode2 = psMetadataLookupU32 (&status, row, "FLAGS2"); - assert (status); - - sources->data[i] = source; - psFree(row); - } - - return sources; -} - -bool pmSourcesWrite_CMF_PS1_V3_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) -{ - bool status; - psArray *table; - psMetadata *row; - psF32 *PAR, *dPAR; - psF32 xPos, yPos; - psF32 xErr, yErr; - int nRow = -1; - char keyword1[80], keyword2[80]; - - // create a header to hold the output data - psMetadata *outhead = psMetadataAlloc (); - - // write the links to the image header - psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); - - pmChip *chip = readout->parent->parent; - pmFPA *fpa = chip->parent; - - // zero point corrections - bool status1 = false; - bool status2 = false; - float magOffset = 0.0; - float exptime = psMetadataLookupF32(&status1, fpa->concepts, "FPA.EXPOSURE"); - float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); - if (!isfinite(zeropt)) { - zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); - } - if (status1 && status2 && (exptime > 0.0)) { - magOffset = zeropt + 2.5*log10(exptime); - } - - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - - table = psArrayAllocEmpty (sources->n); - - // which extended source analyses should we perform? - bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); - bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); - - psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); - psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); - psAssert (radMin->n == radMax->n, "inconsistent annular bins"); - - // write the radial profile apertures to header - for (int i = 0; i < radMax->n; i++) { - sprintf (keyword1, "RMIN_%02d", i); - sprintf (keyword2, "RMAX_%02d", i); - psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]); - psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]); - } - - // we write out all sources, regardless of quality. the source flags tell us the state - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - // skip sources without measurements - if (source->extpars == NULL) continue; - - // we require a PSF model fit (ignore the real crud) - pmModel *model = source->modelPSF; - if (model == NULL) continue; - - // XXX I need to split the extended models from the extended aperture measurements - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - - row = psMetadataAlloc (); - - // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) - psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_EXT", PS_DATA_F32, "EXT model x coordinate", xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT", PS_DATA_F32, "EXT model y coordinate", yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_F32, "Sigma in EXT x coordinate", xErr); - psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); - - float AxialRatio = NAN; - float AxialTheta = NAN; - pmSourceExtendedPars *extpars = source->extpars; - if (extpars) { - AxialRatio = extpars->axes.minor / extpars->axes.major; - AxialTheta = extpars->axes.theta; - } - psMetadataAdd (row, PS_LIST_TAIL, "F25_ARATIO", PS_DATA_F32, "Axial Ratio of radial profile", AxialRatio); - psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA", PS_DATA_F32, "Angle of radial profile ellipse", AxialTheta); - - // Petrosian measurements - // XXX insert header data: petrosian ref radius, flux ratio - // XXX check flags to see if Pet was measured - if (doPetrosian) { - pmSourceExtendedPars *extpars = source->extpars; - if (extpars) { - // XXX note that this mag is either calibrated or instrumental depending on existence of zero point - float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point - float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", mag); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", magErr); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius (pix)", extpars->petrosianRadius); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error (pix)", extpars->petrosianRadiusErr); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50", PS_DATA_F32, "Petrosian R50 (pix)", extpars->petrosianR50); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)", extpars->petrosianR50Err); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90", PS_DATA_F32, "Petrosian R90 (pix)", extpars->petrosianR90); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)", extpars->petrosianR90Err); - } else { - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50", PS_DATA_F32, "Petrosian R50 (pix)", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)",NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90", PS_DATA_F32, "Petrosian R90 (pix)", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)",NAN); - } - } - - // Flux Annuli (if we have extended source measurements, we have these. only optionally save them) - if (doAnnuli) { - psVector *radSB = psVectorAlloc(radMin->n, PS_TYPE_F32); - psVector *radFlux = psVectorAlloc(radMin->n, PS_TYPE_F32); - psVector *radFill = psVectorAlloc(radMin->n, PS_TYPE_F32); - psVectorInit (radSB, NAN); - psVectorInit (radFlux, NAN); - psVectorInit (radFill, NAN); - if (!source->extpars) goto empty_annuli; - if (!source->extpars->radProfile) goto empty_annuli; - if (!source->extpars->radProfile->binSB) goto empty_annuli; - psAssert (source->extpars->radProfile->binSum, "programming error"); - psAssert (source->extpars->radProfile->binFill, "programming error"); - psAssert (source->extpars->radProfile->binSB->n <= radFlux->n, "inconsistent vector lengths"); - psAssert (source->extpars->radProfile->binSum->n <= radFlux->n, "inconsistent vector lengths"); - psAssert (source->extpars->radProfile->binFill->n <= radFlux->n, "inconsistent vector lengths"); - - // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux - for (int j = 0; j < source->extpars->radProfile->binSB->n; j++) { - radSB->data.F32[j] = source->extpars->radProfile->binSB->data.F32[j]; - radFlux->data.F32[j] = source->extpars->radProfile->binSum->data.F32[j]; - radFill->data.F32[j] = source->extpars->radProfile->binFill->data.F32[j]; - } - - empty_annuli: - psMetadataAdd (row, PS_LIST_TAIL, "PROF_SB", PS_DATA_VECTOR, "mean surface brightness annuli", radSB); - psMetadataAdd (row, PS_LIST_TAIL, "PROF_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux); - psMetadataAdd (row, PS_LIST_TAIL, "PROF_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); - psFree (radSB); - psFree (radFlux); - psFree (radFill); - } - if (nRow < 0) { - nRow = row->list->n; - } else { - psAssert (nRow == row->list->n, "inconsistent row lengths"); - } - psArrayAdd (table, 100, row); - psFree (row); - } - - if (table->n == 0) { - if (!psFitsWriteBlank (fits, outhead, extname)) { - psError(psErrorCodeLast(), false, "Unable to write empty sources file."); - psFree(outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - return true; - } - - psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); - if (!psFitsWriteTable (fits, outhead, table, extname)) { - psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); - psFree (outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - - return true; -} - -// XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_V3_XFIT (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) -{ - - psArray *table; - psMetadata *row; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - psF32 xPos, yPos; - psF32 xErr, yErr; - char name[64]; - - // create a header to hold the output data - psMetadata *outhead = psMetadataAlloc (); - - // write the links to the image header - psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); - - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - - // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams - int nParamMax = 0; - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source->modelFits == NULL) continue; - for (int j = 0; j < source->modelFits->n; j++) { - pmModel *model = source->modelFits->data[j]; - assert (model); - nParamMax = PS_MAX (nParamMax, model->params->n); - } - } - - table = psArrayAllocEmpty (sources->n); - - // we write out all sources, regardless of quality. the source flags tell us the state - for (int i = 0; i < sources->n; i++) { - - pmSource *source = sources->data[i]; - - // XXX if no model fits are saved, write out modelEXT? - if (source->modelFits == NULL) continue; - - // We have multiple sources : need to flag the one used to subtract the light (the 'best' model) - for (int j = 0; j < source->modelFits->n; j++) { - - // choose the convolved EXT model, if available, otherwise the simple one - pmModel *model = source->modelFits->data[j]; - assert (model); - - // skip models which were not actually fitted - if (model->flags & PM_MODEL_STATUS_BADARGS) continue; - - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - - axes = pmPSF_ModelToAxes (PAR, 20.0); - - row = psMetadataAlloc (); - - // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) - psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); - psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT", 0, "EXT model x coordinate", xPos); - psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT", 0, "EXT model y coordinate", yPos); - psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG", 0, "Sigma in EXT x coordinate", xErr); - psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG", 0, "Sigma in EXT y coordinate", yErr); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG", 0, "EXT fit instrumental magnitude", model->mag); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude", model->magErr); - - psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS", 0, "number of model parameters", model->params->n); - psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "name of model", pmModelClassGetName (model->type)); - - // XXX these should be major and minor, not 'x' and 'y' - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width in x coordinate", axes.major); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width in y coordinate", axes.minor); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", axes.theta); - - // write out the other generic parameters - for (int k = 0; k < nParamMax; k++) { - if (k == PM_PAR_I0) continue; - if (k == PM_PAR_SKY) continue; - if (k == PM_PAR_XPOS) continue; - if (k == PM_PAR_YPOS) continue; - if (k == PM_PAR_SXX) continue; - if (k == PM_PAR_SXY) continue; - if (k == PM_PAR_SYY) continue; - - snprintf (name, 64, "EXT_PAR_%02d", k); - - if (k < model->params->n) { - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]); - } else { - psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN); - } - } - - // XXX other parameters which may be set. - // XXX flag / value to define the model - // XXX write out the model type, fit status flags - - psArrayAdd (table, 100, row); - psFree (row); - } - } - - if (table->n == 0) { - if (!psFitsWriteBlank (fits, outhead, extname)) { - psError(psErrorCodeLast(), false, "Unable to write empty sources file."); - psFree(outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - return true; - } - - psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); - if (!psFitsWriteTable (fits, outhead, table, extname)) { - psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); - psFree (outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - return true; -} - -bool pmSourcesWrite_CMF_PS1_V3_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) -{ - return true; -} Index: trunk/psModules/src/objects/pmSourceIO_CMP.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMP.c (revision 32046) +++ trunk/psModules/src/objects/pmSourceIO_CMP.c (revision 32347) @@ -135,5 +135,5 @@ lsky = (source->sky < 1.0) ? 0.0 : log10(source->sky); - axes = pmPSF_ModelToAxes (PAR, 20.0); + axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); float psfMagErr = isfinite(source->psfMagErr) ? source->psfMagErr : 999; @@ -293,5 +293,5 @@ goto skip_source; - pmPSF_AxesToModel (PAR, axes); + pmPSF_AxesToModel (PAR, axes, modelType); psArrayAdd (sources, 100, source); Index: trunk/psModules/src/objects/pmSourceIO_OBJ.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_OBJ.c (revision 32046) +++ trunk/psModules/src/objects/pmSourceIO_OBJ.c (revision 32347) @@ -91,5 +91,5 @@ } - axes = pmPSF_ModelToAxes (PAR, 20.0); + axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); psLineInit (line); Index: trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 32046) +++ trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 32347) @@ -113,5 +113,5 @@ yErr = dPAR[PM_PAR_YPOS]; if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { - axes = pmPSF_ModelToAxes (PAR, 20.0); + axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); } else { axes.major = NAN; @@ -288,5 +288,5 @@ dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; - pmPSF_AxesToModel (PAR, axes); + pmPSF_AxesToModel (PAR, axes, modelType); float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); @@ -618,5 +618,5 @@ yErr = dPAR[PM_PAR_YPOS]; - axes = pmPSF_ModelToAxes (PAR, 20.0); + axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); // generate RA,DEC Index: trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c (revision 32046) +++ trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c (revision 32347) @@ -89,5 +89,5 @@ yErr = dPAR[PM_PAR_YPOS]; - axes = pmPSF_ModelToAxes (PAR, 20.0); + axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); } else { // XXX: This code seg faults if source->peak is NULL. @@ -214,5 +214,5 @@ source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); - pmPSF_AxesToModel (PAR, axes); + pmPSF_AxesToModel (PAR, axes, modelType); float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); Index: trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 32046) +++ trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 32347) @@ -95,5 +95,5 @@ yErr = dPAR[PM_PAR_YPOS]; if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { - axes = pmPSF_ModelToAxes (PAR, 20.0); + axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); } else { axes.major = NAN; @@ -257,5 +257,5 @@ dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; - pmPSF_AxesToModel (PAR, axes); + pmPSF_AxesToModel (PAR, axes, modelType); float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); @@ -522,5 +522,5 @@ yErr = dPAR[PM_PAR_YPOS]; - axes = pmPSF_ModelToAxes (PAR, 20.0); + axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); row = psMetadataAlloc (); Index: trunk/psModules/src/objects/pmSourceIO_SMPDATA.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_SMPDATA.c (revision 32046) +++ trunk/psModules/src/objects/pmSourceIO_SMPDATA.c (revision 32347) @@ -92,5 +92,5 @@ lsky = (source->sky < 1.0) ? 0.0 : log10(source->sky); - axes = pmPSF_ModelToAxes (PAR, 20.0); + axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); } else { @@ -190,5 +190,5 @@ axes.theta = psMetadataLookupF32 (&status, row, "THETA"); - pmPSF_AxesToModel (PAR, axes); + pmPSF_AxesToModel (PAR, axes, modelType); source->psfMag = psMetadataLookupF32 (&status, row, "MAG_RAW") - ZERO_POINT; Index: trunk/psModules/src/objects/pmSourceIO_SX.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_SX.c (revision 32046) +++ trunk/psModules/src/objects/pmSourceIO_SX.c (revision 32347) @@ -81,5 +81,5 @@ // pmSourceSextractType (source, &type, &flags); - axes = pmPSF_ModelToAxes (PAR, 20.0); + axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); psLineInit (line); Index: trunk/psModules/src/objects/pmSourceMoments.c =================================================================== --- trunk/psModules/src/objects/pmSourceMoments.c (revision 32046) +++ trunk/psModules/src/objects/pmSourceMoments.c (revision 32347) @@ -96,11 +96,18 @@ // (int) so they can be used in the image index below. - // do 2 passes : the first pass should use a somewhat smaller radius and no sigma window to - // get an unbiased (but probably noisy) centroid - if (!pmSourceMomentsGetCentroid (source, 0.75*radius, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) { - return false; - } - // second pass applies the Gaussian window and uses the centroid from the first pass - if (!pmSourceMomentsGetCentroid (source, radius, sigma, minSN, maskVal, source->moments->Mx, source->moments->My)) { + // XXX // do 2 passes : the first pass should use a somewhat smaller radius and no sigma window to + // XXX // get an unbiased (but probably noisy) centroid + // XXX if (!pmSourceMomentsGetCentroid (source, 0.75*radius, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) { + // XXX return false; + // XXX } + // XXX // second pass applies the Gaussian window and uses the centroid from the first pass + // XXX if (!pmSourceMomentsGetCentroid (source, radius, sigma, minSN, maskVal, source->moments->Mx, source->moments->My)) { + // XXX return false; + // XXX } + + // If we use a large radius for the centroid, it will be biased by any neighbors. The flux + // of any object drops pretty quickly outside 1-2 sigmas. (The exception is bright + // saturated stars, for which we need to use a very large radius here) + if (!pmSourceMomentsGetCentroid (source, 1.5*sigma, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) { return false; } @@ -108,11 +115,4 @@ // Now calculate higher-order moments, using the above-calculated first moments to adjust coordinates // Xn = SUM (x - xc)^n * (z - sky) - - float RFW = 0.0; - float RHW = 0.0; - - float RF = 0.0; - float RH = 0.0; - float RS = 0.0; float XX = 0.0; float XY = 0.0; @@ -143,4 +143,5 @@ float yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage + // calculate the higher-order moments using Xo,Yo for (psS32 row = 0; row < source->pixels->numRows ; row++) { @@ -194,12 +195,5 @@ Sum += pDiff; - // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?) float r = sqrt(r2); - float rf = r * fDiff; - float rh = sqrt(r) * fDiff; - float rs = fDiff; - - float rfw = r * pDiff; - float rhw = sqrt(r) * pDiff; float x = xDiff * pDiff; @@ -221,11 +215,4 @@ float yyyy = yDiff * yyy / r2; - RF += rf; - RH += rh; - RS += rs; - - RFW += rfw; - RHW += rhw; - XX += xx; XY += xy; @@ -242,10 +229,22 @@ XYYY += xyyy; YYYY += yyyy; + + // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?) + // XXX float r = sqrt(r2); + // XXX float rf = r * fDiff; + // XXX float rh = sqrt(r) * fDiff; + // XXX float rs = fDiff; + // XXX + // XXX float rfw = r * pDiff; + // XXX float rhw = sqrt(r) * pDiff; + // XXX + // XXX RF += rf; + // XXX RH += rh; + // XXX RS += rs; + // XXX + // XXX RFW += rfw; + // XXX RHW += rhw; } } - - source->moments->Mrf = RF/RS; - source->moments->Mrh = RH/RS; - source->moments->Mxx = XX/Sum; source->moments->Mxy = XY/Sum; @@ -262,4 +261,64 @@ source->moments->Mxyyy = XYYY/Sum; source->moments->Myyyy = YYYY/Sum; + + // *** now calculate the 1st radial moment (for kron flux) -- symmetrical averaging + + float **vPix = source->pixels->data.F32; + float **vWgt = source->variance->data.F32; + psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; + + float RF = 0.0; + float RH = 0.0; + float RS = 0.0; + + for (psS32 row = 0; row < source->pixels->numRows ; row++) { + + float yDiff = row - yCM; + if (fabs(yDiff) > radius) continue; + + // coordinate of mirror pixel + int yFlip = yCM - yDiff; + if (yFlip < 0) continue; + if (yFlip >= source->pixels->numRows) continue; + + for (psS32 col = 0; col < source->pixels->numCols ; col++) { + // check mask and value for this pixel + if (vMsk && (vMsk[row][col] & maskVal)) continue; + if (isnan(vPix[row][col])) continue; + + float xDiff = col - xCM; + if (fabs(xDiff) > radius) continue; + + // coordinate of mirror pixel + int xFlip = xCM - xDiff; + if (xFlip < 0) continue; + if (xFlip >= source->pixels->numCols) continue; + + // check mask and value for mirror pixel + if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue; + if (isnan(vPix[yFlip][xFlip])) continue; + + // radius is just a function of (xDiff, yDiff) + float r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + if (r2 > R2) continue; + + float fDiff1 = vPix[row][col] - sky; + float fDiff2 = vPix[yFlip][xFlip] - sky; + float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2)); + + // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?) + float r = sqrt(r2); + float rf = r * pDiff; + float rh = sqrt(r) * pDiff; + float rs = 0.5 * (fDiff1 + fDiff2); + + RF += rf; + RH += rh; + RS += rs; + } + } + + source->moments->Mrf = RF/RS; + source->moments->Mrh = RH/RS; // if Mrf (first radial moment) is very small, we are getting into low-significance @@ -270,4 +329,6 @@ kronRefRadius = MIN(radius, kronRefRadius); } + + // *** now calculate the kron flux values using the 1st radial moment float radKinner = 1.0*kronRefRadius; @@ -283,14 +344,130 @@ float SumOuter = 0.0; + // calculate the Kron flux, and related fluxes (NO symmetrical averaging) for (psS32 row = 0; row < source->pixels->numRows ; row++) { - + float yDiff = row - yCM; if (fabs(yDiff) > radKouter) continue; + + for (psS32 col = 0; col < source->pixels->numCols ; col++) { + // check mask and value for this pixel + if (vMsk && (vMsk[row][col] & maskVal)) continue; + if (isnan(vPix[row][col])) continue; + + float xDiff = col - xCM; + if (fabs(xDiff) > radKouter) continue; + + // radKron is just a function of (xDiff, yDiff) + float r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + + float fDiff1 = vPix[row][col] - sky; + float pDiff = fDiff1; + float wDiff = vWgt[row][col]; + + // skip pixels below specified significance level. this is allowed, but should be + // avoided -- the over-weights the wings of bright stars compared to those of faint + // stars. + if (PS_SQR(pDiff) < minSN2*wDiff) continue; + + float r = sqrt(r2); + if (r < radKron) { + Sum += pDiff; + Var += wDiff; + nKronPix ++; + // if (beVerbose) fprintf (stderr, "mome: %d %d %f %f %f\n", col, row, sky, *vPix, Sum); + } + + // use sigma (fixed by psf) not a radKron based value + if (r < sigma) { + SumCore += pDiff; + VarCore += wDiff; + nCorePix ++; + } + + if ((r > radKinner) && (r < radKron)) { + SumInner += pDiff; + nInner ++; + } + if ((r > radKron) && (r < radKouter)) { + SumOuter += pDiff; + nOuter ++; + } + } + } + // *** should I rescale these fluxes by pi R^2 / nNpix? + // XXX source->moments->KronCore = SumCore * M_PI * PS_SQR(sigma) / nCorePix; + // XXX source->moments->KronCoreErr = sqrt(VarCore) * M_PI * PS_SQR(sigma) / nCorePix; + // XXX source->moments->KronFlux = Sum * M_PI * PS_SQR(radKron) / nKronPix; + // XXX source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix; + // XXX source->moments->KronFinner = SumInner * M_PI * (PS_SQR(radKron) - PS_SQR(radKinner)) / nInner; + // XXX source->moments->KronFouter = SumOuter * M_PI * (PS_SQR(radKouter) - PS_SQR(radKron)) / nOuter; + + source->moments->KronCore = SumCore; + source->moments->KronCoreErr = sqrt(VarCore); + source->moments->KronFlux = Sum; + source->moments->KronFluxErr = sqrt(Var); + source->moments->KronFinner = SumInner; + source->moments->KronFouter = SumOuter; + + // XXX not sure I should save this here... + source->moments->KronFluxPSF = source->moments->KronFlux; + source->moments->KronFluxPSFErr = source->moments->KronFluxErr; + source->moments->KronRadiusPSF = source->moments->Mrf; + + psTrace ("psModules.objects", 4, "Mrf: %f KronFlux: %f Mxx: %f Mxy: %f Myy: %f Mxxx: %f Mxxy: %f Mxyy: %f Myyy: %f Mxxxx: %f Mxxxy: %f Mxxyy: %f Mxyyy: %f Mxyyy: %f\n", + source->moments->Mrf, source->moments->KronFlux, + source->moments->Mxx, source->moments->Mxy, source->moments->Myy, + source->moments->Mxxx, source->moments->Mxxy, source->moments->Mxyy, source->moments->Myyy, + source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy); + + psTrace ("psModules.objects", 3, "peak %f %f (%f = %f) Mx: %f My: %f Sum: %f Mxx: %f Mxy: %f Myy: %f sky: %f Npix: %d\n", + source->peak->xf, source->peak->yf, source->peak->rawFlux, sqrt(source->peak->detValue), source->moments->Mx, source->moments->My, Sum, source->moments->Mxx, source->moments->Mxy, source->moments->Myy, sky, source->moments->nPixels); + + return(true); +} + +bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) { + + // First Pass: calculate the first moments (these are subtracted from the coordinates below) + // Sum = SUM (z - sky) + // X1 = SUM (x - xc)*(z - sky) + // .. etc + + float sky = 0.0; + + float peakPixel = -PS_MAX_F32; + psS32 numPixels = 0; + float Sum = 0.0; + float Var = 0.0; + float X1 = 0.0; + float Y1 = 0.0; + float R2 = PS_SQR(radius); + float minSN2 = PS_SQR(minSN); + float rsigma2 = 0.5 / PS_SQR(sigma); + + float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage + float yPeak = yGuess - source->pixels->row0; // coord of peak in subimage + + // we are guaranteed to have a valid pixel and variance at this location (right? right?) + // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]); + // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel"); + // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel"); + // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel"); + + // the moments [Sum(x*f) / Sum(f)] are calculated in pixel index values, and should + // not depend on the fractional pixel location of the source. However, the aperture + // (radius) and the Gaussian window (sigma) depend subtly on the fractional pixel + // position of the expected centroid + + for (psS32 row = 0; row < source->pixels->numRows ; row++) { + + float yDiff = row + 0.5 - yPeak; + if (fabs(yDiff) > radius) continue; float *vPix = source->pixels->data.F32[row]; float *vWgt = source->variance->data.F32[row]; - psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; - // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row]; + psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; + // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row]; for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { @@ -304,117 +481,4 @@ if (isnan(*vPix)) continue; - float xDiff = col - xCM; - if (fabs(xDiff) > radKouter) continue; - - // radKron is just a function of (xDiff, yDiff) - float r2 = PS_SQR(xDiff) + PS_SQR(yDiff); - - float pDiff = *vPix - sky; - float wDiff = *vWgt; - - // skip pixels below specified significance level. this is allowed, but should be - // avoided -- the over-weights the wings of bright stars compared to those of faint - // stars. - if (PS_SQR(pDiff) < minSN2*wDiff) continue; - - float r = sqrt(r2); - if (r < radKron) { - Sum += pDiff; - Var += wDiff; - nKronPix ++; - // if (beVerbose) fprintf (stderr, "mome: %d %d %f %f %f\n", col, row, sky, *vPix, Sum); - } - - // use sigma (fixed by psf) not a radKron based value - if (r < sigma) { - SumCore += pDiff; - VarCore += wDiff; - nCorePix ++; - } - - if ((r > radKinner) && (r < radKron)) { - SumInner += pDiff; - nInner ++; - } - if ((r > radKron) && (r < radKouter)) { - SumOuter += pDiff; - nOuter ++; - } - } - } - // *** should I rescale these fluxes by pi R^2 / nNpix? - source->moments->KronCore = SumCore * M_PI * PS_SQR(sigma) / nCorePix; - source->moments->KronCoreErr = sqrt(VarCore) * M_PI * PS_SQR(sigma) / nCorePix; - source->moments->KronFlux = Sum * M_PI * PS_SQR(radKron) / nKronPix; - source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix; - source->moments->KronFinner = SumInner * M_PI * (PS_SQR(radKron) - PS_SQR(radKinner)) / nInner; - source->moments->KronFouter = SumOuter * M_PI * (PS_SQR(radKouter) - PS_SQR(radKron)) / nOuter; - - psTrace ("psModules.objects", 4, "Mrf: %f KronFlux: %f Mxx: %f Mxy: %f Myy: %f Mxxx: %f Mxxy: %f Mxyy: %f Myyy: %f Mxxxx: %f Mxxxy: %f Mxxyy: %f Mxyyy: %f Mxyyy: %f\n", - source->moments->Mrf, source->moments->KronFlux, - source->moments->Mxx, source->moments->Mxy, source->moments->Myy, - source->moments->Mxxx, source->moments->Mxxy, source->moments->Mxyy, source->moments->Myyy, - source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy); - - psTrace ("psModules.objects", 3, "peak %f %f (%f = %f) Mx: %f My: %f Sum: %f Mxx: %f Mxy: %f Myy: %f sky: %f Npix: %d\n", - source->peak->xf, source->peak->yf, source->peak->rawFlux, sqrt(source->peak->detValue), source->moments->Mx, source->moments->My, Sum, source->moments->Mxx, source->moments->Mxy, source->moments->Myy, sky, source->moments->nPixels); - - return(true); -} - -bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) { - - // First Pass: calculate the first moments (these are subtracted from the coordinates below) - // Sum = SUM (z - sky) - // X1 = SUM (x - xc)*(z - sky) - // .. etc - - float sky = 0.0; - - float peakPixel = -PS_MAX_F32; - psS32 numPixels = 0; - float Sum = 0.0; - float Var = 0.0; - float X1 = 0.0; - float Y1 = 0.0; - float R2 = PS_SQR(radius); - float minSN2 = PS_SQR(minSN); - float rsigma2 = 0.5 / PS_SQR(sigma); - - float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage - float yPeak = yGuess - source->pixels->row0; // coord of peak in subimage - - // we are guaranteed to have a valid pixel and variance at this location (right? right?) - // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]); - // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel"); - // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel"); - // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel"); - - // the moments [Sum(x*f) / Sum(f)] are calculated in pixel index values, and should - // not depend on the fractional pixel location of the source. However, the aperture - // (radius) and the Gaussian window (sigma) depend subtly on the fractional pixel - // position of the expected centroid - - for (psS32 row = 0; row < source->pixels->numRows ; row++) { - - float yDiff = row + 0.5 - yPeak; - if (fabs(yDiff) > radius) continue; - - float *vPix = source->pixels->data.F32[row]; - float *vWgt = source->variance->data.F32[row]; - - psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; - // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row]; - - for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { - if (vMsk) { - if (*vMsk & maskVal) { - vMsk++; - continue; - } - vMsk++; - } - if (isnan(*vPix)) continue; - float xDiff = col + 0.5 - xPeak; if (fabs(xDiff) > radius) continue; Index: trunk/psModules/src/objects/pmSourceOutputs.c =================================================================== --- trunk/psModules/src/objects/pmSourceOutputs.c (revision 32046) +++ trunk/psModules/src/objects/pmSourceOutputs.c (revision 32347) @@ -107,5 +107,5 @@ } if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXY]) && isfinite(PAR[PM_PAR_SYY])) { - axes = pmPSF_ModelToAxes (PAR, 20.0); + axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); outputs->psfMajor = axes.major; outputs->psfMinor = axes.minor; @@ -178,5 +178,67 @@ moments->KronCore = source->moments ? source->moments->KronCore : NAN; moments->KronCoreErr = source->moments ? source->moments->KronCoreErr : NAN; - - return true; -} + moments->KronPSF = source->moments ? source->moments->KronFluxPSF : NAN; + moments->KronPSFErr = source->moments ? source->moments->KronFluxPSFErr : NAN; + + return true; +} + +bool pmSourceLocalAstrometry (psSphere *ptSky, float *posAngle, float *pltScale, pmChip *chip, float xPos, float yPos) { + + pmFPA *fpa = chip->parent; + + if (!chip->toFPA) goto escape; + if (!fpa->toTPA) goto escape; + if (!fpa->toSky) goto escape; + + // generate RA,DEC + psPlane ptCH, ptFP, ptTP_o, ptTP_x, ptTP_y; + + // calculate the astrometry for the coordinate of interest + ptCH.x = xPos; + ptCH.y = yPos; + psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); + psPlaneTransformApply (&ptTP_o, fpa->toTPA, &ptFP); + psDeproject (ptSky, &ptTP_o, fpa->toSky); + + // calculate the astrometry for the coordinate + 1pix in X + ptCH.x = xPos + 1.0; + ptCH.y = yPos; + psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); + psPlaneTransformApply (&ptTP_x, fpa->toTPA, &ptFP); + + // calculate the astrometry for the coordinate + 1pix in Y + ptCH.x = xPos; + ptCH.y = yPos + 1.0; + psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); + psPlaneTransformApply (&ptTP_y, fpa->toTPA, &ptFP); + + // the resulting Tangent Plane coordinates are in TP pixels; convert to local Tangent Plane + // degrees + + float dTPx_dCHx = fpa->toSky->Xs * (ptTP_x.x - ptTP_o.x); + float dTPy_dCHx = fpa->toSky->Ys * (ptTP_x.y - ptTP_o.y); + + float dTPx_dCHy = fpa->toSky->Xs * (ptTP_y.x - ptTP_o.x); + float dTPy_dCHy = fpa->toSky->Ys * (ptTP_y.y - ptTP_o.y); + + float pltScale_x = hypot(dTPx_dCHx, dTPy_dCHx); + float pltScale_y = hypot(dTPx_dCHy, dTPy_dCHy); + *pltScale = 0.5*(pltScale_x + pltScale_y); + + float posAngle_x = atan2 (+dTPy_dCHx, +dTPx_dCHx); + float posAngle_y = atan2 (-dTPy_dCHy, +dTPx_dCHy); + *posAngle = 0.5*(posAngle_x + posAngle_y); + + return true; + +escape: + // no astrometry calibration, give up + ptSky->r = NAN; + ptSky->d = NAN; + *posAngle = NAN; + *pltScale = NAN; + + return false; +} + Index: trunk/psModules/src/objects/pmSourceOutputs.h =================================================================== --- trunk/psModules/src/objects/pmSourceOutputs.h (revision 32046) +++ trunk/psModules/src/objects/pmSourceOutputs.h (revision 32347) @@ -52,4 +52,6 @@ float KronCore; float KronCoreErr; + float KronPSF; + float KronPSFErr; } pmSourceOutputsMoments; Index: trunk/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- trunk/psModules/src/objects/pmSourcePhotometry.c (revision 32046) +++ trunk/psModules/src/objects/pmSourcePhotometry.c (revision 32347) @@ -94,5 +94,5 @@ { PS_ASSERT_PTR_NON_NULL(source, false); - PS_ASSERT_PTR_NON_NULL(psf, false); + // PS_ASSERT_PTR_NON_NULL(psf, false); int status = false; Index: trunk/psModules/src/objects/pmSourcePlotPSFModel.c =================================================================== --- trunk/psModules/src/objects/pmSourcePlotPSFModel.c (revision 32046) +++ trunk/psModules/src/objects/pmSourcePlotPSFModel.c (revision 32347) @@ -145,5 +145,5 @@ // force the axis ratio to be < 20.0 psEllipseAxes axes_mnt = psEllipseMomentsToAxes (moments, 20.0); - psEllipseAxes axes_psf = pmPSF_ModelToAxes (PAR, 20.0); + psEllipseAxes axes_psf = pmPSF_ModelToAxes (PAR, 20.0, model->type); // moments major axis