26 #include "ndarray/eigen.h"
39 FlagDefinitionList flagDefinitions;
44 flagDefinitions.add(
"flag_edge",
"Object too close to edge");
46 flagDefinitions.add(
"flag_noSecondDerivative",
"Vanishing second derivative");
48 flagDefinitions.add(
"flag_almostNoSecondDerivative",
"Almost vanishing second derivative");
50 flagDefinitions.add(
"flag_notAtMaximum",
"Object is not at a maximum");
58 float const AMPAST4 = 1.33;
67 static int inter4(
float vm,
float v0,
float vp,
float *cen,
bool negative =
false) {
68 float const sp = v0 - vp;
69 float const sm = v0 - vm;
70 float const d2 = sp + sm;
71 float const s = 0.5 * (vp - vm);
73 if ((!negative && (d2 <= 0.0f || v0 <= 0.0f)) || (negative && (d2 >= 0.0f || v0 >= 0.0f))) {
77 *cen = s / d2 * (1.0 + AMPAST4 * sp * sm / (d2 * v0));
79 return fabs(*cen) < 1 ? 0 : 1;
86 float astrom_errors(
float skyVar,
96 double const k = quarticBad ? 0 : AMPAST4;
109 sVar += 0.5 * sourceVar * exp(-1 / (2 * tau2));
110 dVar += sourceVar * (4 * exp(-1 / (2 * tau2)) + 2 * exp(-1 / (2 * tau2)));
119 xVar = sVar *
pow(1 / d + k / (4 * As) * (1 - 12 * s * s / (d * d)), 2) +
120 dVar *
pow(s / (d * d) - k / (4 * As) * 8 * s * s / (d * d * d), 2);
122 return (xVar >= 0 ?
sqrt(xVar) : NAN);
130 template <
typename ImageXy_locatorT,
typename VarImageXy_locatorT>
131 void doMeasureCentroidImpl(
double *xCenter,
135 double *sizeX2,
double *sizeY2,
137 VarImageXy_locatorT vim,
138 double smoothingSigma,
139 FlagHandler flagHandler) {
143 double const d2x = 2 * im(0, 0) - im(-1, 0) - im(1, 0);
144 double const d2y = 2 * im(0, 0) - im(0, -1) - im(0, 1);
145 double const sx = 0.5 * (im(1, 0) - im(-1, 0));
146 double const sy = 0.5 * (im(0, 1) - im(0, -1));
148 if (d2x == 0.0 || d2y == 0.0) {
152 if (d2x < 0.0 || d2y < 0.0) {
155 (
boost::format(
": d2I/dx2, d2I/dy2 = %g %g") % d2x % d2y).str(),
159 double const dx0 = sx / d2x;
160 double const dy0 = sy / d2y;
162 if (
fabs(dx0) > 10.0 ||
fabs(dy0) > 10.0) {
166 (
boost::format(
": sx, d2x, sy, d2y = %f %f %f %f") % sx % d2x % sy % d2y).str(),
170 double vpk = im(0, 0) + 0.5 * (sx * dx0 + sy * dy0);
177 float m0x = 0, m1x = 0, m2x = 0;
178 float m0y = 0, m1y = 0, m2y = 0;
181 quarticBad += inter4(im(-1, -1), im(0, -1), im(1, -1), &m0x);
182 quarticBad += inter4(im(-1, 0), im(0, 0), im(1, 0), &m1x);
183 quarticBad += inter4(im(-1, 1), im(0, 1), im(1, 1), &m2x);
185 quarticBad += inter4(im(-1, -1), im(-1, 0), im(-1, 1), &m0y);
186 quarticBad += inter4(im(0, -1), im(0, 0), im(0, 1), &m1y);
187 quarticBad += inter4(im(1, -1), im(1, 0), im(1, 1), &m2y);
190 double sigmaX2, sigmaY2;
198 double const smx = 0.5 * (m2x - m0x);
199 double const smy = 0.5 * (m2y - m0y);
200 double const dm2x = m1x - 0.5 * (m0x + m2x);
201 double const dm2y = m1y - 0.5 * (m0y + m2y);
202 double const dx = m1x + dy0 * (smx - dy0 * dm2x);
203 double const dy = m1y + dx0 * (smy - dx0 * dm2y);
204 double const dx4 = m1x + dy * (smx - dy * dm2x);
205 double const dy4 = m1y + dx * (smy - dx * dm2y);
209 sigmaX2 = vpk / d2x - (1 + 6 * dx0 * dx0) / 4;
210 sigmaY2 = vpk / d2y - (1 + 6 * dy0 * dy0) / 4;
215 float tauX2 = sigmaX2;
216 float tauY2 = sigmaY2;
217 tauX2 -= smoothingSigma * smoothingSigma;
218 tauY2 -= smoothingSigma * smoothingSigma;
220 if (tauX2 <= smoothingSigma * smoothingSigma) {
221 tauX2 = smoothingSigma * smoothingSigma;
223 if (tauY2 <= smoothingSigma * smoothingSigma) {
224 tauY2 = smoothingSigma * smoothingSigma;
227 float const skyVar = (vim(-1, -1) + vim(0, -1) + vim(1, -1) + vim(-1, 0) + vim(1, 0) + vim(-1, 1) +
228 vim(0, 1) + vim(1, 1)) /
230 float const sourceVar = vim(0, 0);
231 float const A = vpk *
sqrt((sigmaX2 / tauX2) * (sigmaY2 / tauY2));
236 *dxc = astrom_errors(skyVar, sourceVar, A, tauX2, vpk, sx, d2x,
fabs(smoothingSigma), quarticBad);
237 *dyc = astrom_errors(skyVar, sourceVar, A, tauY2, vpk, sy, d2y,
fabs(smoothingSigma), quarticBad);
243 template <
typename MaskedImageXy_locatorT>
244 void doMeasureCentroidImpl(
double *xCenter,
248 double *sizeX2,
double *sizeY2,
250 MaskedImageXy_locatorT mim,
251 double smoothingSigma,
252 bool negative, FlagHandler flagHandler) {
256 double const d2x = 2 * mim.image(0, 0) - mim.image(-1, 0) - mim.image(1, 0);
257 double const d2y = 2 * mim.image(0, 0) - mim.image(0, -1) - mim.image(0, 1);
258 double const sx = 0.5 * (mim.image(1, 0) - mim.image(-1, 0));
259 double const sy = 0.5 * (mim.image(0, 1) - mim.image(0, -1));
261 if (d2x == 0.0 || d2y == 0.0) {
265 if ((!negative && (d2x < 0.0 || d2y < 0.0)) || (negative && (d2x > 0.0 || d2y > 0.0))) {
268 (
boost::format(
": d2I/dx2, d2I/dy2 = %g %g") % d2x % d2y).str(),
272 double const dx0 = sx / d2x;
273 double const dy0 = sy / d2y;
275 if (
fabs(dx0) > 10.0 ||
fabs(dy0) > 10.0) {
279 (
boost::format(
": sx, d2x, sy, d2y = %f %f %f %f") % sx % d2x % sy % d2y).str(),
283 double vpk = mim.image(0, 0) + 0.5 * (sx * dx0 + sy * dy0);
290 float m0x = 0, m1x = 0, m2x = 0;
291 float m0y = 0, m1y = 0, m2y = 0;
294 quarticBad += inter4(mim.image(-1, -1), mim.image(0, -1), mim.image(1, -1), &m0x, negative);
295 quarticBad += inter4(mim.image(-1, 0), mim.image(0, 0), mim.image(1, 0), &m1x, negative);
296 quarticBad += inter4(mim.image(-1, 1), mim.image(0, 1), mim.image(1, 1), &m2x, negative);
298 quarticBad += inter4(mim.image(-1, -1), mim.image(-1, 0), mim.image(-1, 1), &m0y, negative);
299 quarticBad += inter4(mim.image(0, -1), mim.image(0, 0), mim.image(0, 1), &m1y, negative);
300 quarticBad += inter4(mim.image(1, -1), mim.image(1, 0), mim.image(1, 1), &m2y, negative);
303 double sigmaX2, sigmaY2;
311 double const smx = 0.5 * (m2x - m0x);
312 double const smy = 0.5 * (m2y - m0y);
313 double const dm2x = m1x - 0.5 * (m0x + m2x);
314 double const dm2y = m1y - 0.5 * (m0y + m2y);
315 double const dx = m1x + dy0 * (smx - dy0 * dm2x);
316 double const dy = m1y + dx0 * (smy - dx0 * dm2y);
317 double const dx4 = m1x + dy * (smx - dy * dm2x);
318 double const dy4 = m1y + dx * (smy - dx * dm2y);
322 sigmaX2 = vpk / d2x - (1 + 6 * dx0 * dx0) / 4;
323 sigmaY2 = vpk / d2y - (1 + 6 * dy0 * dy0) / 4;
328 float tauX2 = sigmaX2;
329 float tauY2 = sigmaY2;
330 tauX2 -= smoothingSigma * smoothingSigma;
331 tauY2 -= smoothingSigma * smoothingSigma;
333 if (tauX2 <= smoothingSigma * smoothingSigma) {
334 tauX2 = smoothingSigma * smoothingSigma;
336 if (tauY2 <= smoothingSigma * smoothingSigma) {
337 tauY2 = smoothingSigma * smoothingSigma;
341 (mim.variance(-1, -1) + mim.variance(0, -1) + mim.variance(1, -1) + mim.variance(-1, 0) +
342 mim.variance(1, 0) + mim.variance(-1, 1) + mim.variance(0, 1) + mim.variance(1, 1)) /
344 float const sourceVar = mim.variance(0, 0);
345 float const A = vpk *
sqrt((sigmaX2 / tauX2) * (sigmaY2 / tauY2));
350 *dxc = astrom_errors(skyVar, sourceVar, A, tauX2, vpk, sx, d2x,
fabs(smoothingSigma), quarticBad);
351 *dyc = astrom_errors(skyVar, sourceVar, A, tauY2, vpk, sy, d2y,
fabs(smoothingSigma), quarticBad);
359 template <
typename MaskedImageT>
361 const int y, MaskedImageT
const &mimage,
int binX,
int binY,
362 FlagHandler _flagHandler) {
364 afw::geom::ellipses::Quadrupole
const &shape =
psf->computeShape(center);
365 double const smoothingSigma = shape.getDeterminantRadius();
367 double const nEffective =
psf->computeEffectiveArea();
369 double const nEffective = 4 *
M_PI * smoothingSigma * smoothingSigma;
373 int const kWidth = kernel->getWidth();
374 int const kHeight = kernel->getHeight();
377 geom::ExtentI(binX * (3 + kWidth + 1), binY * (3 + kHeight + 1)));
383 }
catch (pex::exceptions::LengthError &err) {
388 binnedImage->setXY0(subImage->getXY0());
390 MaskedImageT smoothedImage = MaskedImageT(*binnedImage,
true);
391 assert(smoothedImage.getWidth() / 2 == kWidth / 2 + 2);
392 assert(smoothedImage.getHeight() / 2 == kHeight / 2 + 2);
394 afw::math::convolve(smoothedImage, *binnedImage, *kernel, afw::math::ConvolutionControl());
395 *smoothedImage.getVariance() *= binX * binY * nEffective;
410 _centroidChecker(
schema,
name, ctrl.doFootprintCheck, ctrl.maxDistToPeak) {}
414 geom::Point2D center = _centroidExtractor(measRecord, _flagHandler);
421 typedef MaskedImageT::Image ImageT;
422 typedef MaskedImageT::Variance VarianceT;
423 bool negative =
false;
425 negative = measRecord.
get(measRecord.
getSchema().
find<afw::table::Flag>(
"flags_negative").key);
430 ImageT
const &
image = *mimage.getImage();
448 double xc = 0., yc = 0., dxc = 0., dyc = 0.;
449 for (
int binsize = 1; binsize <= _ctrl.
binmax; binsize *= 2) {
451 smoothAndBinImage(
psf,
x,
y, mimage, binX, binY, _flagHandler);
452 MaskedImageT
const smoothedImage =
result.first;
453 double const smoothingSigma =
result.second;
455 MaskedImageT::xy_locator mim =
456 smoothedImage.xy_at(smoothedImage.getWidth() / 2, smoothedImage.getHeight() / 2);
458 double sizeX2, sizeY2;
461 doMeasureCentroidImpl(&xc, &dxc, &yc, &dyc, &sizeX2, &sizeY2, &peakVal, mim, smoothingSigma, negative,
466 xc = (xc + 0.5) * binX - 0.5;
468 sizeX2 *= binX * binX;
470 yc = (yc + 0.5) * binY - 0.5;
472 sizeY2 *= binY * binY;
478 double const fac = _ctrl.
wfac * (1 + smoothingSigma * smoothingSigma);
479 double const facX2 = fac * binX * binX;
480 double const facY2 = fac * binY * binY;
482 if (sizeX2 < facX2 && ::pow(xc -
x, 2) < facX2 && sizeY2 < facY2 && ::pow(yc -
y, 2) < facY2) {
483 if (binsize > 1 || _ctrl.
peakMin < 0.0 || peakVal > _ctrl.
peakMin) {
488 if (sizeX2 >= facX2 || ::pow(xc -
x, 2) >= facX2) {
491 if (sizeY2 >= facY2 || ::pow(yc -
y, 2) >= facY2) {
498 result.xErr = sqrt(dxc * dxc);
499 result.yErr = sqrt(dyc * dyc);
501 _centroidChecker(measRecord);
505 _flagHandler.handleFailure(measRecord,
error);
514 if (
mapper.getInputSchema().getNames().count(
mapper.getInputSchema().join(
name, flag.
name)) == 0)
517 mapper.getInputSchema().find<afw::table::Flag>(
name +
"_" + flag.
name).key;
table::Key< std::string > name
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
table::Key< double > sigma2
afw::table::Key< double > sigma
This implements the SdssCentroid algorithm within the meas_base measurement framework.
MaskedImageT getMaskedImage()
Return the MaskedImage.
std::shared_ptr< lsst::afw::detection::Psf const > getPsf() const
Return the Exposure's Psf object.
A class to manipulate images, masks, and variance as a single object.
Schema getSchema() const
Return the Schema that holds this record's fields and keys.
void set(Key< T > const &key, U const &value)
Set value of a field for the given key.
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Defines the fields and offsets for a table.
SchemaItem< T > find(std::string const &name) const
Find a SchemaItem in the Schema by name.
A mapping between the keys of two Schemas, used to copy data between them.
Record class that contains measurements made on a single exposure.
An integer coordinate rectangle.
A FunctorKey for CentroidResult.
Exception to be thrown when a measurement algorithm experiences a fatal error.
vector-type utility class to build a collection of FlagDefinitions
std::size_t size() const
return the current size (number of defined elements) of the collection
Utility class for handling flag fields that indicate the failure modes of an algorithm.
Exception to be thrown when a measurement algorithm experiences a known failure mode.
static FlagDefinition const NOT_AT_MAXIMUM
virtual void fail(afw::table::SourceRecord &measRecord, MeasurementError *error=nullptr) const
Handle an exception thrown by the current algorithm by setting flags in the given record.
SdssCentroidAlgorithm(Control const &ctrl, std::string const &name, afw::table::Schema &schema)
static FlagDefinition const NO_SECOND_DERIVATIVE
static FlagDefinition const FAILURE
virtual void measure(afw::table::SourceRecord &measRecord, afw::image::Exposure< float > const &exposure) const
Called to measure a single child source in an image.
static FlagDefinitionList const & getFlagDefinitions()
static FlagDefinition const ALMOST_NO_SECOND_DERIVATIVE
static FlagDefinition const EDGE
A C++ control class to handle SdssCentroidAlgorithm's configuration.
double wfac
"fiddle factor for adjusting the binning" ;
int binmax
"maximum allowed binning" ;
double peakMin
"if the peak's less than this insist on binning at least once" ;
Provides consistent interface for LSST exceptions.
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
double indexToPosition(double ind)
Convert image index to image position.
void convolve(OutImageT &convolvedImage, InImageT const &inImage, KernelT const &kernel, ConvolutionControl const &convolutionControl=ConvolutionControl())
Convolve an Image or MaskedImage with a Kernel, setting pixels of an existing output image.
@ MEAN
estimate sample mean
std::shared_ptr< ImageT > binImage(ImageT const &inImage, int const binX, int const binY, lsst::afw::math::Property const flags=lsst::afw::math::MEAN)
constexpr double PI
The ratio of a circle's circumference to diameter.
@ SIGMA_ONLY
Only the diagonal elements of the covariance matrix are provided.
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
A base class for image defects.
A reusable struct for centroid measurements.
Simple class used to define and document flags The name and doc constitute the identity of the FlagDe...