LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
SdssCentroid.cc
Go to the documentation of this file.
1// -*- lsst-c++ -*-
2/*
3 * LSST Data Management System
4 * Copyright 2008-2013 LSST Corporation.
5 *
6 * This product includes software developed by the
7 * LSST Project (http://www.lsst.org/).
8 *
9 * This program is free software: you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation, either version 3 of the License, or
12 * (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the LSST License Statement and
20 * the GNU General Public License along with this program. If not,
21 * see <http://www.lsstcorp.org/LegalNotices/>.
22 */
23#include <iostream>
24#include <cmath>
25#include <numeric>
26#include "ndarray/eigen.h"
27#include "lsst/geom/Angle.h"
30#include "lsst/pex/exceptions.h"
35
36namespace lsst {
37namespace meas {
38namespace base {
39namespace {
40FlagDefinitionList flagDefinitions;
41} // namespace
42
43FlagDefinition const SdssCentroidAlgorithm::FAILURE = flagDefinitions.addFailureFlag();
44FlagDefinition const SdssCentroidAlgorithm::EDGE =
45 flagDefinitions.add("flag_edge", "Object too close to edge; peak used.");
47 flagDefinitions.add("flag_noSecondDerivative", "Vanishing second derivative");
49 flagDefinitions.add("flag_almostNoSecondDerivative", "Almost vanishing second derivative");
50FlagDefinition const SdssCentroidAlgorithm::NOT_AT_MAXIMUM =
51 flagDefinitions.add("flag_notAtMaximum", "Object is not at a maximum");
52FlagDefinition const SdssCentroidAlgorithm::NEAR_EDGE =
53 flagDefinitions.add("flag_near_edge", "Object close to edge; fallback kernel used.");
54
56
57namespace {
58
59/************************************************************************************************************/
60
61float const AMPAST4 = 1.33; // amplitude of `4th order' corr compared to theory
62
63/*
64 * Do the Gaussian quartic interpolation for the position
65 * of the maximum for three equally spaced values vm,v0,vp, assumed to be at
66 * abscissae -1,0,1; the answer is returned as *cen
67 *
68 * Return 0 is all is well, otherwise 1
69 */
70static int inter4(float vm, float v0, float vp, float *cen, bool negative = false) {
71 float const sp = v0 - vp;
72 float const sm = v0 - vm;
73 float const d2 = sp + sm;
74 float const s = 0.5 * (vp - vm);
75
76 if ((!negative && (d2 <= 0.0f || v0 <= 0.0f)) || (negative && (d2 >= 0.0f || v0 >= 0.0f))) {
77 return (1);
78 }
79
80 *cen = s / d2 * (1.0 + AMPAST4 * sp * sm / (d2 * v0));
81
82 return fabs(*cen) < 1 ? 0 : 1;
83}
84
85/*****************************************************************************/
86/*
87 * Calculate error in centroid
88 */
89float astrom_errors(float skyVar, // variance of pixels at the sky level
90 float sourceVar, // variance in peak due to excess counts over sky
91 float A, // abs(peak value in raw image)
92 float tau2, // Object is N(0,tau2)
93 float As, // abs(peak value in smoothed image)
94 float s, // slope across central pixel
95 float d, // curvature at central pixel
96 float sigma, // width of smoothing filter
97 int quarticBad) { // was quartic estimate bad?
98
99 double const k = quarticBad ? 0 : AMPAST4; /* quartic correction coeff */
100 double const sigma2 = sigma * sigma; /* == sigma^2 */
101 double sVar, dVar; /* variances of s and d */
102 double xVar; /* variance of centroid, x */
103
105 return (1e3);
106 }
107
108 if (sigma <= 0) { /* no smoothing; no covariance */
109 sVar = 0.5 * skyVar; /* due to sky */
110 dVar = 6 * skyVar;
111
112 sVar += 0.5 * sourceVar * exp(-1 / (2 * tau2));
113 dVar += sourceVar * (4 * exp(-1 / (2 * tau2)) + 2 * exp(-1 / (2 * tau2)));
114 } else { /* smoothed */
115 sVar = skyVar / (8 * geom::PI * sigma2) * (1 - exp(-1 / sigma2));
116 dVar = skyVar / (2 * geom::PI * sigma2) * (3 - 4 * exp(-1 / (4 * sigma2)) + exp(-1 / sigma2));
117
118 sVar += sourceVar / (12 * geom::PI * sigma2) * (exp(-1 / (3 * sigma2)) - exp(-1 / sigma2));
119 dVar += sourceVar / (3 * geom::PI * sigma2) * (2 - 3 * exp(-1 / (3 * sigma2)) + exp(-1 / sigma2));
120 }
121
122 xVar = sVar * pow(1 / d + k / (4 * As) * (1 - 12 * s * s / (d * d)), 2) +
123 dVar * pow(s / (d * d) - k / (4 * As) * 8 * s * s / (d * d * d), 2);
124
125 return (xVar >= 0 ? sqrt(xVar) : NAN);
126}
127
128/************************************************************************************************************/
129/*
130 * Estimate the position of an object, assuming we know that it's approximately the size of the PSF
131 */
132
133template <typename ImageXy_locatorT, typename VarImageXy_locatorT>
134int doMeasureCentroidImpl(double *xCenter, // output; x-position of object
135 double *dxc, // output; error in xCenter
136 double *yCenter, // output; y-position of object
137 double *dyc, // output; error in yCenter
138 double *sizeX2, double *sizeY2, // output; object widths^2 in x and y directions
139 ImageXy_locatorT im, // Locator for the pixel values
140 VarImageXy_locatorT vim, // Locator for the image containing the variance
141 double smoothingSigma, // Gaussian sigma of already-applied smoothing filter
142 FlagHandler flagHandler) {
143 /*
144 * find a first quadratic estimate
145 */
146 double const d2x = 2 * im(0, 0) - im(-1, 0) - im(1, 0);
147 double const d2y = 2 * im(0, 0) - im(0, -1) - im(0, 1);
148 double const sx = 0.5 * (im(1, 0) - im(-1, 0));
149 double const sy = 0.5 * (im(0, 1) - im(0, -1));
150
151 if (d2x == 0.0 || d2y == 0.0) {
153 }
154 if (d2x < 0.0 || d2y < 0.0) {
156 }
157
158 double const dx0 = sx / d2x;
159 double const dy0 = sy / d2y; // first guess
160
161 if (fabs(dx0) > 10.0 || fabs(dy0) > 10.0) {
163 }
164
165 double vpk = im(0, 0) + 0.5 * (sx * dx0 + sy * dy0); // height of peak in image
166 if (vpk < 0) {
167 vpk = -vpk;
168 }
169 /*
170 * now evaluate maxima on stripes
171 */
172 float m0x = 0, m1x = 0, m2x = 0;
173 float m0y = 0, m1y = 0, m2y = 0;
174
175 int quarticBad = 0;
176 quarticBad += inter4(im(-1, -1), im(0, -1), im(1, -1), &m0x);
177 quarticBad += inter4(im(-1, 0), im(0, 0), im(1, 0), &m1x);
178 quarticBad += inter4(im(-1, 1), im(0, 1), im(1, 1), &m2x);
179
180 quarticBad += inter4(im(-1, -1), im(-1, 0), im(-1, 1), &m0y);
181 quarticBad += inter4(im(0, -1), im(0, 0), im(0, 1), &m1y);
182 quarticBad += inter4(im(1, -1), im(1, 0), im(1, 1), &m2y);
183
184 double xc, yc; // position of maximum
185 double sigmaX2, sigmaY2; // widths^2 in x and y of smoothed object
186
187 if (quarticBad) { // >= 1 quartic interpolator is bad
188 xc = dx0;
189 yc = dy0;
190 sigmaX2 = vpk / d2x; // widths^2 in x
191 sigmaY2 = vpk / d2y; // and y
192 } else {
193 double const smx = 0.5 * (m2x - m0x);
194 double const smy = 0.5 * (m2y - m0y);
195 double const dm2x = m1x - 0.5 * (m0x + m2x);
196 double const dm2y = m1y - 0.5 * (m0y + m2y);
197 double const dx = m1x + dy0 * (smx - dy0 * dm2x); // first quartic approx
198 double const dy = m1y + dx0 * (smy - dx0 * dm2y);
199 double const dx4 = m1x + dy * (smx - dy * dm2x); // second quartic approx
200 double const dy4 = m1y + dx * (smy - dx * dm2y);
201
202 xc = dx4;
203 yc = dy4;
204 sigmaX2 = vpk / d2x - (1 + 6 * dx0 * dx0) / 4; // widths^2 in x
205 sigmaY2 = vpk / d2y - (1 + 6 * dy0 * dy0) / 4; // and y
206 }
207 /*
208 * Now for the errors.
209 */
210 float tauX2 = sigmaX2; // width^2 of _un_ smoothed object
211 float tauY2 = sigmaY2;
212 tauX2 -= smoothingSigma * smoothingSigma; // correct for smoothing
213 tauY2 -= smoothingSigma * smoothingSigma;
214
215 if (tauX2 <= smoothingSigma * smoothingSigma) { // problem; sigmaX2 must be bad
216 tauX2 = smoothingSigma * smoothingSigma;
217 }
218 if (tauY2 <= smoothingSigma * smoothingSigma) { // sigmaY2 must be bad
219 tauY2 = smoothingSigma * smoothingSigma;
220 }
221
222 float const skyVar = (vim(-1, -1) + vim(0, -1) + vim(1, -1) + vim(-1, 0) + vim(1, 0) + vim(-1, 1) +
223 vim(0, 1) + vim(1, 1)) /
224 8.0; // Variance in sky
225 float const sourceVar = vim(0, 0); // extra variance of peak due to its photons
226 float const A = vpk * sqrt((sigmaX2 / tauX2) * (sigmaY2 / tauY2)); // peak of Unsmoothed object
227
228 *xCenter = xc;
229 *yCenter = yc;
230
231 *dxc = astrom_errors(skyVar, sourceVar, A, tauX2, vpk, sx, d2x, fabs(smoothingSigma), quarticBad);
232 *dyc = astrom_errors(skyVar, sourceVar, A, tauY2, vpk, sy, d2y, fabs(smoothingSigma), quarticBad);
233
234 *sizeX2 = tauX2; // return the estimates of the (object size)^2
235 *sizeY2 = tauY2;
236 return 0;
237}
238
239template <typename MaskedImageXy_locatorT>
240int doMeasureCentroidImpl(double *xCenter, // output; x-position of object
241 double *dxc, // output; error in xCenter
242 double *yCenter, // output; y-position of object
243 double *dyc, // output; error in yCenter
244 double *sizeX2, double *sizeY2, // output; object widths^2 in x and y directions
245 double *peakVal, // output; peak of object
246 MaskedImageXy_locatorT mim, // Locator for the pixel values
247 double smoothingSigma, // Gaussian sigma of already-applied smoothing filter
248 bool negative, FlagHandler flagHandler) {
249 /*
250 * find a first quadratic estimate
251 */
252 double const d2x = 2 * mim.image(0, 0) - mim.image(-1, 0) - mim.image(1, 0);
253 double const d2y = 2 * mim.image(0, 0) - mim.image(0, -1) - mim.image(0, 1);
254 double const sx = 0.5 * (mim.image(1, 0) - mim.image(-1, 0));
255 double const sy = 0.5 * (mim.image(0, 1) - mim.image(0, -1));
256
257 if (d2x == 0.0 || d2y == 0.0) {
259 }
260 if ((!negative && (d2x < 0.0 || d2y < 0.0)) || (negative && (d2x > 0.0 || d2y > 0.0))) {
262 }
263
264 double const dx0 = sx / d2x;
265 double const dy0 = sy / d2y; // first guess
266
267 if (fabs(dx0) > 10.0 || fabs(dy0) > 10.0) {
269 }
270
271 double vpk = mim.image(0, 0) + 0.5 * (sx * dx0 + sy * dy0); // height of peak in image
272 // if (vpk < 0) {
273 // vpk = -vpk;
274 //}
275 /*
276 * now evaluate maxima on stripes
277 */
278 float m0x = 0, m1x = 0, m2x = 0;
279 float m0y = 0, m1y = 0, m2y = 0;
280
281 int quarticBad = 0;
282 quarticBad += inter4(mim.image(-1, -1), mim.image(0, -1), mim.image(1, -1), &m0x, negative);
283 quarticBad += inter4(mim.image(-1, 0), mim.image(0, 0), mim.image(1, 0), &m1x, negative);
284 quarticBad += inter4(mim.image(-1, 1), mim.image(0, 1), mim.image(1, 1), &m2x, negative);
285
286 quarticBad += inter4(mim.image(-1, -1), mim.image(-1, 0), mim.image(-1, 1), &m0y, negative);
287 quarticBad += inter4(mim.image(0, -1), mim.image(0, 0), mim.image(0, 1), &m1y, negative);
288 quarticBad += inter4(mim.image(1, -1), mim.image(1, 0), mim.image(1, 1), &m2y, negative);
289
290 double xc, yc; // position of maximum
291 double sigmaX2, sigmaY2; // widths^2 in x and y of smoothed object
292
293 if (quarticBad) { // >= 1 quartic interpolator is bad
294 xc = dx0;
295 yc = dy0;
296 sigmaX2 = vpk / d2x; // widths^2 in x
297 sigmaY2 = vpk / d2y; // and y
298 } else {
299 double const smx = 0.5 * (m2x - m0x);
300 double const smy = 0.5 * (m2y - m0y);
301 double const dm2x = m1x - 0.5 * (m0x + m2x);
302 double const dm2y = m1y - 0.5 * (m0y + m2y);
303 double const dx = m1x + dy0 * (smx - dy0 * dm2x); // first quartic approx
304 double const dy = m1y + dx0 * (smy - dx0 * dm2y);
305 double const dx4 = m1x + dy * (smx - dy * dm2x); // second quartic approx
306 double const dy4 = m1y + dx * (smy - dx * dm2y);
307
308 xc = dx4;
309 yc = dy4;
310 sigmaX2 = vpk / d2x - (1 + 6 * dx0 * dx0) / 4; // widths^2 in x
311 sigmaY2 = vpk / d2y - (1 + 6 * dy0 * dy0) / 4; // and y
312 }
313 /*
314 * Now for the errors.
315 */
316 float tauX2 = sigmaX2; // width^2 of _un_ smoothed object
317 float tauY2 = sigmaY2;
318 tauX2 -= smoothingSigma * smoothingSigma; // correct for smoothing
319 tauY2 -= smoothingSigma * smoothingSigma;
320
321 if (tauX2 <= smoothingSigma * smoothingSigma) { // problem; sigmaX2 must be bad
322 tauX2 = smoothingSigma * smoothingSigma;
323 }
324 if (tauY2 <= smoothingSigma * smoothingSigma) { // sigmaY2 must be bad
325 tauY2 = smoothingSigma * smoothingSigma;
326 }
327
328 float const skyVar =
329 (mim.variance(-1, -1) + mim.variance(0, -1) + mim.variance(1, -1) + mim.variance(-1, 0) +
330 mim.variance(1, 0) + mim.variance(-1, 1) + mim.variance(0, 1) + mim.variance(1, 1)) /
331 8.0; // Variance in sky
332 float const sourceVar = mim.variance(0, 0); // extra variance of peak due to its photons
333 float const A = vpk * sqrt((sigmaX2 / tauX2) * (sigmaY2 / tauY2)); // peak of Unsmoothed object
334
335 *xCenter = xc;
336 *yCenter = yc;
337
338 *dxc = astrom_errors(skyVar, sourceVar, A, tauX2, vpk, sx, d2x, fabs(smoothingSigma), quarticBad);
339 *dyc = astrom_errors(skyVar, sourceVar, A, tauY2, vpk, sy, d2y, fabs(smoothingSigma), quarticBad);
340
341 *sizeX2 = tauX2; // return the estimates of the (object size)^2
342 *sizeY2 = tauY2;
343
344 *peakVal = vpk;
345 return 0;
346}
347
348template <typename MaskedImageT>
350 const int y, MaskedImageT const &mimage, int binX, int binY,
351 FlagHandler _flagHandler) {
352 geom::Point2D const center(x + mimage.getX0(), y + mimage.getY0());
353 afw::geom::ellipses::Quadrupole const &shape = psf->computeShape(center);
354 double const smoothingSigma = shape.getDeterminantRadius();
355#if 0
356 double const nEffective = psf->computeEffectiveArea(); // not implemented yet (#2821)
357#else
358 double const nEffective = 4 * M_PI * smoothingSigma * smoothingSigma; // correct for a Gaussian
359#endif
360
361 std::shared_ptr<afw::math::Kernel const> kernel = psf->getLocalKernel(center);
362 int const kWidth = kernel->getWidth();
363 int const kHeight = kernel->getHeight();
364
365 geom::BoxI bbox(geom::Point2I(x - binX * (2 + kWidth / 2), y - binY * (2 + kHeight / 2)),
366 geom::ExtentI(binX * (3 + kWidth + 1), binY * (3 + kHeight + 1)));
367
368 // image to smooth, a shallow copy
370 if (!mimage.getBBox(afw::image::LOCAL).contains(bbox)) {
371 return std::make_tuple(mimage, 0, SdssCentroidAlgorithm::EDGE.number);
372 }
373 subImage.reset(new MaskedImageT(mimage, bbox, afw::image::LOCAL));
374 std::shared_ptr<MaskedImageT> binnedImage = afw::math::binImage(*subImage, binX, binY, afw::math::MEAN);
375 binnedImage->setXY0(subImage->getXY0());
376 // image to smooth into, a deep copy.
377 MaskedImageT smoothedImage = MaskedImageT(*binnedImage, true);
378 if(smoothedImage.getWidth() / 2 != kWidth / 2 + 2 || // assumed by the code that uses smoothedImage
379 smoothedImage.getHeight() / 2 != kHeight / 2 + 2) {
380 throw LSST_EXCEPT(lsst::pex::exceptions::LengthError, "invalid image dimensions");
381 }
382
383 afw::math::convolve(smoothedImage, *binnedImage, *kernel, afw::math::ConvolutionControl());
384 *smoothedImage.getVariance() *= binX * binY * nEffective; // We want the per-pixel variance, so undo the
385 // effects of binning and smoothing
386
387 return std::make_tuple(smoothedImage, smoothingSigma, 0);
388}
389
390} // end anonymous namespace
391
393 afw::table::Schema &schema)
394 : _ctrl(ctrl),
395 _centroidKey(CentroidResultKey::addFields(schema, name, "centroid from Sdss Centroid algorithm",
396 SIGMA_ONLY)),
397 _flagHandler(FlagHandler::addFields(schema, name, getFlagDefinitions())),
398 _centroidExtractor(schema, name, true),
399 _centroidChecker(schema, name, ctrl.doFootprintCheck, ctrl.maxDistToPeak) {}
401 afw::image::Exposure<float> const &exposure) const {
402 // get our current best guess about the centroid: either a centroider measurement or peak.
403 geom::Point2D center = _centroidExtractor(measRecord, _flagHandler);
405 result.x = center.getX();
406 result.y = center.getY();
407 measRecord.set(_centroidKey, result); // better than NaN
408
410 typedef MaskedImageT::Image ImageT;
411 typedef MaskedImageT::Variance VarianceT;
412 bool negative = false;
413 try {
414 negative = measRecord.get(measRecord.getSchema().find<afw::table::Flag>("flags_negative").key);
415 } catch (pexExcept::Exception &e) {
416 }
417
418 MaskedImageT const &mimage = exposure.getMaskedImage();
419 ImageT const &image = *mimage.getImage();
420 std::shared_ptr<afw::detection::Psf const> psf = exposure.getPsf();
421
422 int const x = image.positionToIndex(center.getX(), afw::image::X).first;
423 int const y = image.positionToIndex(center.getY(), afw::image::Y).first;
424
425 if (!image.getBBox().contains(geom::Extent2I(x, y) + image.getXY0())) {
426 _flagHandler.setValue(measRecord, EDGE.number, true);
427 _flagHandler.setValue(measRecord, SdssCentroidAlgorithm::FAILURE.number, true);
428 return;
429 }
430
431 // Algorithm uses a least-squares fit (implemented via a convolution) to a symmetrized PSF model.
432 // If you don't have a Psf, you need to use another centroider, such as GaussianCentroider.
433 if (!psf) {
434 throw LSST_EXCEPT(FatalAlgorithmError, "SdssCentroid algorithm requires a Psf with every exposure");
435 }
436
437 int binX = 1;
438 int binY = 1;
439 double xc = 0., yc = 0., dxc = 0., dyc = 0.; // estimated centre and error therein
440 bool stopBinning = false;
441 for (int binsize = 1; binsize <= _ctrl.binmax; binsize *= 2) {
443 smoothAndBinImage(psf, x, y, mimage, binX, binY, _flagHandler);
444 int errorFlag = std::get<2>(smoothResult);
445 if (errorFlag == static_cast<int>(EDGE.number)) {
446 psf = std::make_shared<afw::detection::GaussianPsf>(5, 5, 0.5);
447 smoothResult = smoothAndBinImage(psf, x, y, mimage, binX, binY, _flagHandler);
448 stopBinning = true;
449 errorFlag = std::get<2>(smoothResult);
450 if (errorFlag == 0) {
451 errorFlag = NEAR_EDGE.number;
452 }
453 }
454 if (errorFlag > 0) {
455 _flagHandler.setValue(measRecord, errorFlag, true);
456 _flagHandler.setValue(measRecord, SdssCentroidAlgorithm::FAILURE.number, true);
457 // if NEAR_EDGE is not a fatal error we continue otherwise return
458 if (errorFlag != static_cast<int>(NEAR_EDGE.number)) {
459 return;
460 }
461 }
462 MaskedImageT const smoothedImage = std::get<0>(smoothResult);
463 double const smoothingSigma = std::get<1>(smoothResult);
464
465 MaskedImageT::xy_locator mim =
466 smoothedImage.xy_at(smoothedImage.getWidth() / 2, smoothedImage.getHeight() / 2);
467
468 double sizeX2, sizeY2; // object widths^2 in x and y directions
469 double peakVal; // peak intensity in image
470
471 errorFlag = doMeasureCentroidImpl(&xc, &dxc, &yc, &dyc, &sizeX2, &sizeY2, &peakVal, mim, smoothingSigma, negative,
472 _flagHandler);
473 if (errorFlag > 0) {
474 _flagHandler.setValue(measRecord, errorFlag, true);
475 _flagHandler.setValue(measRecord, SdssCentroidAlgorithm::FAILURE.number, true);
476 return;
477 }
478
479 if (binsize > 1) {
480 // dilate from the lower left corner of central pixel
481 xc = (xc + 0.5) * binX - 0.5;
482 dxc *= binX;
483 sizeX2 *= binX * binX;
484
485 yc = (yc + 0.5) * binY - 0.5;
486 dyc *= binY;
487 sizeY2 *= binY * binY;
488 }
489
490 xc += x; // xc, yc are measured relative to pixel (x, y)
491 yc += y;
492
493 if (stopBinning) {
494 break;
495 }
496
497 double const fac = _ctrl.wfac * (1 + smoothingSigma * smoothingSigma);
498 double const facX2 = fac * binX * binX;
499 double const facY2 = fac * binY * binY;
500
501 if (sizeX2 < facX2 && ::pow(xc - x, 2) < facX2 && sizeY2 < facY2 && ::pow(yc - y, 2) < facY2) {
502 if (binsize > 1 || _ctrl.peakMin < 0.0 || peakVal > _ctrl.peakMin) {
503 break;
504 }
505 }
506
507 if (sizeX2 >= facX2 || ::pow(xc - x, 2) >= facX2) {
508 binX *= 2;
509 }
510 if (sizeY2 >= facY2 || ::pow(yc - y, 2) >= facY2) {
511 binY *= 2;
512 }
513 }
514 result.x = afw::image::indexToPosition(xc + image.getX0());
515 result.y = afw::image::indexToPosition(yc + image.getY0());
516
517 result.xErr = sqrt(dxc * dxc);
518 result.yErr = sqrt(dyc * dyc);
519 measRecord.set(_centroidKey, result);
520 _centroidChecker(measRecord);
521}
522
524 _flagHandler.handleFailure(measRecord, error);
525}
526
529 : CentroidTransform{name, mapper} {
532 if (flag == SdssCentroidAlgorithm::FAILURE) continue;
533 if (mapper.getInputSchema().getNames().count(mapper.getInputSchema().join(name, flag.name)) == 0)
534 continue;
536 mapper.getInputSchema().find<afw::table::Flag>(name + "_" + flag.name).key;
537 mapper.addMapping(key);
538 }
539}
540
541} // namespace base
542} // namespace meas
543} // namespace lsst
AmpInfoBoxKey bbox
Definition Amplifier.cc:117
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
table::Key< double > sigma2
afw::table::Key< double > sigma
#define M_PI
Definition ListMatch.cc:31
SchemaMapper * mapper
This implements the SdssCentroid algorithm within the meas_base measurement framework.
int y
Definition SpanSet.cc:48
A class to contain the data, WCS, and other information needed to describe an image of the sky.
Definition Exposure.h:72
A class to manipulate images, masks, and variance as a single object.
Definition MaskedImage.h:74
Tag types used to declare specialized field types.
Definition misc.h:31
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Definition BaseRecord.h:151
Schema getSchema() const
Return the Schema that holds this record's fields and keys.
Definition BaseRecord.h:80
void set(Key< T > const &key, U const &value)
Set value of a field for the given key.
Definition BaseRecord.h:164
Defines the fields and offsets for a table.
Definition Schema.h:51
SchemaItem< T > find(std::string const &name) const
Find a SchemaItem in the Schema by name.
Definition Schema.cc:467
A mapping between the keys of two Schemas, used to copy data between them.
Record class that contains measurements made on a single exposure.
Definition Source.h:78
An integer coordinate rectangle.
Definition Box.h:55
A FunctorKey for CentroidResult.
Base for centroid measurement transformations.
Exception to be thrown when a measurement algorithm experiences a fatal error.
Definition exceptions.h:76
vector-type utility class to build a collection of FlagDefinitions
Definition FlagHandler.h:60
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.
Definition exceptions.h:48
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
static FlagDefinition const NEAR_EDGE
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" ;
SdssCentroidTransform(Control const &ctrl, std::string const &name, afw::table::SchemaMapper &mapper)
Provides consistent interface for LSST exceptions.
Definition Exception.h:107
Reports attempts to exceed implementation-defined length limits for some classes.
Definition Runtime.h:76
T fabs(T... args)
T make_tuple(T... args)
double indexToPosition(double ind)
Convert image index to image position.
Definition ImageUtils.h:55
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
Definition Statistics.h:57
std::shared_ptr< ImageT > binImage(ImageT const &inImage, int const binX, int const binY, lsst::afw::math::Property const flags=lsst::afw::math::MEAN)
Definition binImage.cc:44
double constexpr PI
The ratio of a circle's circumference to diameter.
Definition Angle.h:40
@ SIGMA_ONLY
Only the diagonal elements of the covariance matrix are provided.
Definition constants.h:45
T pow(T... args)
T reset(T... args)
A reusable struct for centroid measurements.
CentroidElement x
x (column) coordinate of the measured position
Simple class used to define and document flags The name and doc constitute the identity of the FlagDe...
Definition FlagHandler.h:40