LSSTApplications  20.0.0
LSSTDataManagementBasePackage
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"
28 #include "lsst/afw/detection/Psf.h"
29 #include "lsst/pex/exceptions.h"
30 #include "lsst/afw/math/ConvolveImage.h"
31 #include "lsst/afw/math/offsetImage.h"
32 #include "lsst/afw/table/Source.h"
33 #include "lsst/meas/base/SdssCentroid.h"
34 
35 namespace lsst {
36 namespace meas {
37 namespace base {
38 namespace {
39 FlagDefinitionList flagDefinitions;
40 } // namespace
41 
42 FlagDefinition const SdssCentroidAlgorithm::FAILURE = flagDefinitions.addFailureFlag();
43 FlagDefinition const SdssCentroidAlgorithm::EDGE =
44  flagDefinitions.add("flag_edge", "Object too close to edge");
45 FlagDefinition const SdssCentroidAlgorithm::NO_SECOND_DERIVATIVE =
46  flagDefinitions.add("flag_noSecondDerivative", "Vanishing second derivative");
47 FlagDefinition const SdssCentroidAlgorithm::ALMOST_NO_SECOND_DERIVATIVE =
48  flagDefinitions.add("flag_almostNoSecondDerivative", "Almost vanishing second derivative");
49 FlagDefinition const SdssCentroidAlgorithm::NOT_AT_MAXIMUM =
50  flagDefinitions.add("flag_notAtMaximum", "Object is not at a maximum");
51 
52 FlagDefinitionList const &SdssCentroidAlgorithm::getFlagDefinitions() { return flagDefinitions; }
53 
54 namespace {
55 
56 /************************************************************************************************************/
57 
58 float const AMPAST4 = 1.33; // amplitude of `4th order' corr compared to theory
59 
60 /*
61  * Do the Gaussian quartic interpolation for the position
62  * of the maximum for three equally spaced values vm,v0,vp, assumed to be at
63  * abscissae -1,0,1; the answer is returned as *cen
64  *
65  * Return 0 is all is well, otherwise 1
66  */
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);
72 
73  if ((!negative && (d2 <= 0.0f || v0 <= 0.0f)) || (negative && (d2 >= 0.0f || v0 >= 0.0f))) {
74  return (1);
75  }
76 
77  *cen = s / d2 * (1.0 + AMPAST4 * sp * sm / (d2 * v0));
78 
79  return fabs(*cen) < 1 ? 0 : 1;
80 }
81 
82 /*****************************************************************************/
83 /*
84  * Calculate error in centroid
85  */
86 float astrom_errors(float skyVar, // variance of pixels at the sky level
87  float sourceVar, // variance in peak due to excess counts over sky
88  float A, // abs(peak value in raw image)
89  float tau2, // Object is N(0,tau2)
90  float As, // abs(peak value in smoothed image)
91  float s, // slope across central pixel
92  float d, // curvature at central pixel
93  float sigma, // width of smoothing filter
94  int quarticBad) { // was quartic estimate bad?
95 
96  double const k = quarticBad ? 0 : AMPAST4; /* quartic correction coeff */
97  double const sigma2 = sigma * sigma; /* == sigma^2 */
98  double sVar, dVar; /* variances of s and d */
99  double xVar; /* variance of centroid, x */
100 
101  if (fabs(As) < std::numeric_limits<float>::min() || fabs(d) < std::numeric_limits<float>::min()) {
102  return (1e3);
103  }
104 
105  if (sigma <= 0) { /* no smoothing; no covariance */
106  sVar = 0.5 * skyVar; /* due to sky */
107  dVar = 6 * skyVar;
108 
109  sVar += 0.5 * sourceVar * exp(-1 / (2 * tau2));
110  dVar += sourceVar * (4 * exp(-1 / (2 * tau2)) + 2 * exp(-1 / (2 * tau2)));
111  } else { /* smoothed */
112  sVar = skyVar / (8 * geom::PI * sigma2) * (1 - exp(-1 / sigma2));
113  dVar = skyVar / (2 * geom::PI * sigma2) * (3 - 4 * exp(-1 / (4 * sigma2)) + exp(-1 / sigma2));
114 
115  sVar += sourceVar / (12 * geom::PI * sigma2) * (exp(-1 / (3 * sigma2)) - exp(-1 / sigma2));
116  dVar += sourceVar / (3 * geom::PI * sigma2) * (2 - 3 * exp(-1 / (3 * sigma2)) + exp(-1 / sigma2));
117  }
118 
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);
121 
122  return (xVar >= 0 ? sqrt(xVar) : NAN);
123 }
124 
125 /************************************************************************************************************/
126 /*
127  * Estimate the position of an object, assuming we know that it's approximately the size of the PSF
128  */
129 
130 template <typename ImageXy_locatorT, typename VarImageXy_locatorT>
131 void doMeasureCentroidImpl(double *xCenter, // output; x-position of object
132  double *dxc, // output; error in xCenter
133  double *yCenter, // output; y-position of object
134  double *dyc, // output; error in yCenter
135  double *sizeX2, double *sizeY2, // output; object widths^2 in x and y directions
136  ImageXy_locatorT im, // Locator for the pixel values
137  VarImageXy_locatorT vim, // Locator for the image containing the variance
138  double smoothingSigma, // Gaussian sigma of already-applied smoothing filter
139  FlagHandler flagHandler) {
140  /*
141  * find a first quadratic estimate
142  */
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));
147 
148  if (d2x == 0.0 || d2y == 0.0) {
149  throw LSST_EXCEPT(MeasurementError, SdssCentroidAlgorithm::NO_SECOND_DERIVATIVE.doc,
150  SdssCentroidAlgorithm::NO_SECOND_DERIVATIVE.number);
151  }
152  if (d2x < 0.0 || d2y < 0.0) {
153  throw LSST_EXCEPT(MeasurementError,
154  SdssCentroidAlgorithm::NOT_AT_MAXIMUM.doc +
155  (boost::format(": d2I/dx2, d2I/dy2 = %g %g") % d2x % d2y).str(),
156  SdssCentroidAlgorithm::NOT_AT_MAXIMUM.number);
157  }
158 
159  double const dx0 = sx / d2x;
160  double const dy0 = sy / d2y; // first guess
161 
162  if (fabs(dx0) > 10.0 || fabs(dy0) > 10.0) {
163  throw LSST_EXCEPT(
164  MeasurementError,
165  SdssCentroidAlgorithm::ALMOST_NO_SECOND_DERIVATIVE.doc +
166  (boost::format(": sx, d2x, sy, d2y = %f %f %f %f") % sx % d2x % sy % d2y).str(),
167  SdssCentroidAlgorithm::ALMOST_NO_SECOND_DERIVATIVE.number);
168  }
169 
170  double vpk = im(0, 0) + 0.5 * (sx * dx0 + sy * dy0); // height of peak in image
171  if (vpk < 0) {
172  vpk = -vpk;
173  }
174  /*
175  * now evaluate maxima on stripes
176  */
177  float m0x = 0, m1x = 0, m2x = 0;
178  float m0y = 0, m1y = 0, m2y = 0;
179 
180  int quarticBad = 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);
184 
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);
188 
189  double xc, yc; // position of maximum
190  double sigmaX2, sigmaY2; // widths^2 in x and y of smoothed object
191 
192  if (quarticBad) { // >= 1 quartic interpolator is bad
193  xc = dx0;
194  yc = dy0;
195  sigmaX2 = vpk / d2x; // widths^2 in x
196  sigmaY2 = vpk / d2y; // and y
197  } else {
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); // first quartic approx
203  double const dy = m1y + dx0 * (smy - dx0 * dm2y);
204  double const dx4 = m1x + dy * (smx - dy * dm2x); // second quartic approx
205  double const dy4 = m1y + dx * (smy - dx * dm2y);
206 
207  xc = dx4;
208  yc = dy4;
209  sigmaX2 = vpk / d2x - (1 + 6 * dx0 * dx0) / 4; // widths^2 in x
210  sigmaY2 = vpk / d2y - (1 + 6 * dy0 * dy0) / 4; // and y
211  }
212  /*
213  * Now for the errors.
214  */
215  float tauX2 = sigmaX2; // width^2 of _un_ smoothed object
216  float tauY2 = sigmaY2;
217  tauX2 -= smoothingSigma * smoothingSigma; // correct for smoothing
218  tauY2 -= smoothingSigma * smoothingSigma;
219 
220  if (tauX2 <= smoothingSigma * smoothingSigma) { // problem; sigmaX2 must be bad
221  tauX2 = smoothingSigma * smoothingSigma;
222  }
223  if (tauY2 <= smoothingSigma * smoothingSigma) { // sigmaY2 must be bad
224  tauY2 = smoothingSigma * smoothingSigma;
225  }
226 
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)) /
229  8.0; // Variance in sky
230  float const sourceVar = vim(0, 0); // extra variance of peak due to its photons
231  float const A = vpk * sqrt((sigmaX2 / tauX2) * (sigmaY2 / tauY2)); // peak of Unsmoothed object
232 
233  *xCenter = xc;
234  *yCenter = yc;
235 
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);
238 
239  *sizeX2 = tauX2; // return the estimates of the (object size)^2
240  *sizeY2 = tauY2;
241 }
242 
243 template <typename MaskedImageXy_locatorT>
244 void doMeasureCentroidImpl(double *xCenter, // output; x-position of object
245  double *dxc, // output; error in xCenter
246  double *yCenter, // output; y-position of object
247  double *dyc, // output; error in yCenter
248  double *sizeX2, double *sizeY2, // output; object widths^2 in x and y directions
249  double *peakVal, // output; peak of object
250  MaskedImageXy_locatorT mim, // Locator for the pixel values
251  double smoothingSigma, // Gaussian sigma of already-applied smoothing filter
252  bool negative, FlagHandler flagHandler) {
253  /*
254  * find a first quadratic estimate
255  */
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));
260 
261  if (d2x == 0.0 || d2y == 0.0) {
262  throw LSST_EXCEPT(MeasurementError, SdssCentroidAlgorithm::NO_SECOND_DERIVATIVE.doc,
263  SdssCentroidAlgorithm::NO_SECOND_DERIVATIVE.number);
264  }
265  if ((!negative && (d2x < 0.0 || d2y < 0.0)) || (negative && (d2x > 0.0 || d2y > 0.0))) {
266  throw LSST_EXCEPT(MeasurementError,
267  SdssCentroidAlgorithm::NOT_AT_MAXIMUM.doc +
268  (boost::format(": d2I/dx2, d2I/dy2 = %g %g") % d2x % d2y).str(),
269  SdssCentroidAlgorithm::NOT_AT_MAXIMUM.number);
270  }
271 
272  double const dx0 = sx / d2x;
273  double const dy0 = sy / d2y; // first guess
274 
275  if (fabs(dx0) > 10.0 || fabs(dy0) > 10.0) {
276  throw LSST_EXCEPT(
277  MeasurementError,
278  SdssCentroidAlgorithm::NO_SECOND_DERIVATIVE.doc +
279  (boost::format(": sx, d2x, sy, d2y = %f %f %f %f") % sx % d2x % sy % d2y).str(),
280  SdssCentroidAlgorithm::ALMOST_NO_SECOND_DERIVATIVE.number);
281  }
282 
283  double vpk = mim.image(0, 0) + 0.5 * (sx * dx0 + sy * dy0); // height of peak in image
284  // if (vpk < 0) {
285  // vpk = -vpk;
286  //}
287  /*
288  * now evaluate maxima on stripes
289  */
290  float m0x = 0, m1x = 0, m2x = 0;
291  float m0y = 0, m1y = 0, m2y = 0;
292 
293  int quarticBad = 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);
297 
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);
301 
302  double xc, yc; // position of maximum
303  double sigmaX2, sigmaY2; // widths^2 in x and y of smoothed object
304 
305  if (quarticBad) { // >= 1 quartic interpolator is bad
306  xc = dx0;
307  yc = dy0;
308  sigmaX2 = vpk / d2x; // widths^2 in x
309  sigmaY2 = vpk / d2y; // and y
310  } else {
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); // first quartic approx
316  double const dy = m1y + dx0 * (smy - dx0 * dm2y);
317  double const dx4 = m1x + dy * (smx - dy * dm2x); // second quartic approx
318  double const dy4 = m1y + dx * (smy - dx * dm2y);
319 
320  xc = dx4;
321  yc = dy4;
322  sigmaX2 = vpk / d2x - (1 + 6 * dx0 * dx0) / 4; // widths^2 in x
323  sigmaY2 = vpk / d2y - (1 + 6 * dy0 * dy0) / 4; // and y
324  }
325  /*
326  * Now for the errors.
327  */
328  float tauX2 = sigmaX2; // width^2 of _un_ smoothed object
329  float tauY2 = sigmaY2;
330  tauX2 -= smoothingSigma * smoothingSigma; // correct for smoothing
331  tauY2 -= smoothingSigma * smoothingSigma;
332 
333  if (tauX2 <= smoothingSigma * smoothingSigma) { // problem; sigmaX2 must be bad
334  tauX2 = smoothingSigma * smoothingSigma;
335  }
336  if (tauY2 <= smoothingSigma * smoothingSigma) { // sigmaY2 must be bad
337  tauY2 = smoothingSigma * smoothingSigma;
338  }
339 
340  float const skyVar =
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)) /
343  8.0; // Variance in sky
344  float const sourceVar = mim.variance(0, 0); // extra variance of peak due to its photons
345  float const A = vpk * sqrt((sigmaX2 / tauX2) * (sigmaY2 / tauY2)); // peak of Unsmoothed object
346 
347  *xCenter = xc;
348  *yCenter = yc;
349 
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);
352 
353  *sizeX2 = tauX2; // return the estimates of the (object size)^2
354  *sizeY2 = tauY2;
355 
356  *peakVal = vpk;
357 }
358 
359 template <typename MaskedImageT>
360 std::pair<MaskedImageT, double> smoothAndBinImage(CONST_PTR(afw::detection::Psf) psf, int const x,
361  const int y, MaskedImageT const &mimage, int binX, int binY,
362  FlagHandler _flagHandler) {
363  geom::Point2D const center(x + mimage.getX0(), y + mimage.getY0());
364  afw::geom::ellipses::Quadrupole const &shape = psf->computeShape(center);
365  double const smoothingSigma = shape.getDeterminantRadius();
366 #if 0
367  double const nEffective = psf->computeEffectiveArea(); // not implemented yet (#2821)
368 #else
369  double const nEffective = 4 * M_PI * smoothingSigma * smoothingSigma; // correct for a Gaussian
370 #endif
371 
372  std::shared_ptr<afw::math::Kernel const> kernel = psf->getLocalKernel(center);
373  int const kWidth = kernel->getWidth();
374  int const kHeight = kernel->getHeight();
375 
376  geom::BoxI bbox(geom::Point2I(x - binX * (2 + kWidth / 2), y - binY * (2 + kHeight / 2)),
377  geom::ExtentI(binX * (3 + kWidth + 1), binY * (3 + kHeight + 1)));
378 
379  // image to smooth, a shallow copy
380  PTR(MaskedImageT) subImage;
381  try {
382  subImage.reset(new MaskedImageT(mimage, bbox, afw::image::LOCAL));
383  } catch (pex::exceptions::LengthError &err) {
384  throw LSST_EXCEPT(MeasurementError, SdssCentroidAlgorithm::EDGE.doc,
385  SdssCentroidAlgorithm::EDGE.number);
386  }
387  PTR(MaskedImageT) binnedImage = afw::math::binImage(*subImage, binX, binY, afw::math::MEAN);
388  binnedImage->setXY0(subImage->getXY0());
389  // image to smooth into, a deep copy.
390  MaskedImageT smoothedImage = MaskedImageT(*binnedImage, true);
391  assert(smoothedImage.getWidth() / 2 == kWidth / 2 + 2); // assumed by the code that uses smoothedImage
392  assert(smoothedImage.getHeight() / 2 == kHeight / 2 + 2);
393 
394  afw::math::convolve(smoothedImage, *binnedImage, *kernel, afw::math::ConvolutionControl());
395  *smoothedImage.getVariance() *= binX * binY * nEffective; // We want the per-pixel variance, so undo the
396  // effects of binning and smoothing
397 
398  return std::make_pair(smoothedImage, smoothingSigma);
399 }
400 
401 } // end anonymous namespace
402 
403 SdssCentroidAlgorithm::SdssCentroidAlgorithm(Control const &ctrl, std::string const &name,
404  afw::table::Schema &schema)
405  : _ctrl(ctrl),
406  _centroidKey(CentroidResultKey::addFields(schema, name, "centroid from Sdss Centroid algorithm",
407  SIGMA_ONLY)),
408  _flagHandler(FlagHandler::addFields(schema, name, getFlagDefinitions())),
409  _centroidExtractor(schema, name, true),
410  _centroidChecker(schema, name, ctrl.doFootprintCheck, ctrl.maxDistToPeak) {}
411 void SdssCentroidAlgorithm::measure(afw::table::SourceRecord &measRecord,
412  afw::image::Exposure<float> const &exposure) const {
413  // get our current best guess about the centroid: either a centroider measurement or peak.
414  geom::Point2D center = _centroidExtractor(measRecord, _flagHandler);
415  CentroidResult result;
416  result.x = center.getX();
417  result.y = center.getY();
418  measRecord.set(_centroidKey, result); // better than NaN
419 
420  typedef afw::image::Exposure<float>::MaskedImageT MaskedImageT;
421  typedef MaskedImageT::Image ImageT;
422  typedef MaskedImageT::Variance VarianceT;
423  bool negative = false;
424  try {
425  negative = measRecord.get(measRecord.getSchema().find<afw::table::Flag>("flags_negative").key);
426  } catch (pexExcept::Exception &e) {
427  }
428 
429  MaskedImageT const &mimage = exposure.getMaskedImage();
430  ImageT const &image = *mimage.getImage();
431  CONST_PTR(afw::detection::Psf) psf = exposure.getPsf();
432 
433  int const x = image.positionToIndex(center.getX(), afw::image::X).first;
434  int const y = image.positionToIndex(center.getY(), afw::image::Y).first;
435 
436  if (!image.getBBox().contains(geom::Extent2I(x, y) + image.getXY0())) {
437  throw LSST_EXCEPT(meas::base::MeasurementError, EDGE.doc, EDGE.number);
438  }
439 
440  // Algorithm uses a least-squares fit (implemented via a convolution) to a symmetrized PSF model.
441  // If you don't have a Psf, you need to use another centroider, such as GaussianCentroider.
442  if (!psf) {
443  throw LSST_EXCEPT(FatalAlgorithmError, "SdssCentroid algorithm requires a Psf with every exposure");
444  }
445 
446  int binX = 1;
447  int binY = 1;
448  double xc = 0., yc = 0., dxc = 0., dyc = 0.; // estimated centre and error therein
449  for (int binsize = 1; binsize <= _ctrl.binmax; binsize *= 2) {
450  std::pair<MaskedImageT, double> result =
451  smoothAndBinImage(psf, x, y, mimage, binX, binY, _flagHandler);
452  MaskedImageT const smoothedImage = result.first;
453  double const smoothingSigma = result.second;
454 
455  MaskedImageT::xy_locator mim =
456  smoothedImage.xy_at(smoothedImage.getWidth() / 2, smoothedImage.getHeight() / 2);
457 
458  double sizeX2, sizeY2; // object widths^2 in x and y directions
459  double peakVal; // peak intensity in image
460 
461  doMeasureCentroidImpl(&xc, &dxc, &yc, &dyc, &sizeX2, &sizeY2, &peakVal, mim, smoothingSigma, negative,
462  _flagHandler);
463 
464  if (binsize > 1) {
465  // dilate from the lower left corner of central pixel
466  xc = (xc + 0.5) * binX - 0.5;
467  dxc *= binX;
468  sizeX2 *= binX * binX;
469 
470  yc = (yc + 0.5) * binY - 0.5;
471  dyc *= binY;
472  sizeY2 *= binY * binY;
473  }
474 
475  xc += x; // xc, yc are measured relative to pixel (x, y)
476  yc += y;
477 
478  double const fac = _ctrl.wfac * (1 + smoothingSigma * smoothingSigma);
479  double const facX2 = fac * binX * binX;
480  double const facY2 = fac * binY * binY;
481 
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) {
484  break;
485  }
486  }
487 
488  if (sizeX2 >= facX2 || ::pow(xc - x, 2) >= facX2) {
489  binX *= 2;
490  }
491  if (sizeY2 >= facY2 || ::pow(yc - y, 2) >= facY2) {
492  binY *= 2;
493  }
494  }
495  result.x = afw::image::indexToPosition(xc + image.getX0());
496  result.y = afw::image::indexToPosition(yc + image.getY0());
497 
498  result.xErr = sqrt(dxc * dxc);
499  result.yErr = sqrt(dyc * dyc);
500  measRecord.set(_centroidKey, result);
501  _centroidChecker(measRecord);
502 }
503 
504 void SdssCentroidAlgorithm::fail(afw::table::SourceRecord &measRecord, MeasurementError *error) const {
505  _flagHandler.handleFailure(measRecord, error);
506 }
507 
508 SdssCentroidTransform::SdssCentroidTransform(Control const &ctrl, std::string const &name,
509  afw::table::SchemaMapper &mapper)
510  : CentroidTransform{name, mapper} {
511  for (std::size_t i = 0; i < SdssCentroidAlgorithm::getFlagDefinitions().size(); i++) {
512  FlagDefinition const &flag = SdssCentroidAlgorithm::getFlagDefinitions()[i];
513  if (flag == SdssCentroidAlgorithm::FAILURE) continue;
514  if (mapper.getInputSchema().getNames().count(mapper.getInputSchema().join(name, flag.name)) == 0)
515  continue;
516  afw::table::Key<afw::table::Flag> key =
517  mapper.getInputSchema().find<afw::table::Flag>(name + "_" + flag.name).key;
518  mapper.addMapping(key);
519  }
520 }
521 
522 } // namespace base
523 } // namespace meas
524 } // namespace lsst
y
int y
Definition: SpanSet.cc:49
schema
table::Schema schema
Definition: Amplifier.cc:115
lsst::afw::image
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Definition: imageAlgorithm.dox:1
lsst::afw::image::Exposure::getPsf
std::shared_ptr< lsst::afw::detection::Psf const > getPsf() const
Return the Exposure's Psf object.
Definition: Exposure.h:307
offsetImage.h
CONST_PTR
#define CONST_PTR(...)
A shared pointer to a const object.
Definition: base.h:47
lsst::afw::image::LOCAL
@ LOCAL
Definition: ImageBase.h:94
lsst::geom::PI
constexpr double PI
The ratio of a circle's circumference to diameter.
Definition: Angle.h:39
lsst::meas::base::FlagDefinition::number
std::size_t number
Definition: FlagHandler.h:54
std::string
STL class.
std::shared_ptr
STL class.
Psf.h
lsst::log.log.logContinued.error
def error(fmt, *args)
Definition: logContinued.py:210
std::fabs
T fabs(T... args)
lsst::afw::table::SourceRecord
Record class that contains measurements made on a single exposure.
Definition: Source.h:80
lsst::afw::table::BaseRecord::get
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Definition: BaseRecord.h:151
lsst::afw::image::Exposure< float >
base
Definition: __init__.py:1
std::pair
lsst::meas::base::SIGMA_ONLY
@ SIGMA_ONLY
Only the diagonal elements of the covariance matrix are provided.
Definition: constants.h:45
psf
Key< int > psf
Definition: Exposure.cc:65
lsst::meas::base::FlagDefinition::doc
std::string doc
Definition: FlagHandler.h:53
pex.config.history.format
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
sigma
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:50
SdssCentroid.h
lsst::meas::base::MeasurementError
Exception to be thrown when a measurement algorithm experiences a known failure mode.
Definition: exceptions.h:48
lsst::meas::base::SdssCentroidAlgorithm::FAILURE
static FlagDefinition const FAILURE
Definition: SdssCentroid.h:73
lsst::meas::base::FatalAlgorithmError
Exception to be thrown when a measurement algorithm experiences a fatal error.
Definition: exceptions.h:76
lsst::afw::table::Schema
Defines the fields and offsets for a table.
Definition: Schema.h:50
lsst::afw::image::indexToPosition
double indexToPosition(double ind)
Convert image index to image position.
Definition: ImageUtils.h:55
lsst::meas::base::CentroidResultKey
A FunctorKey for CentroidResult.
Definition: CentroidUtilities.h:88
lsst::afw::geom.transform.transformContinued.name
string name
Definition: transformContinued.py:32
lsst::afw::image::X
@ X
Definition: ImageUtils.h:36
lsst::meas::base::CentroidTransform
Base for centroid measurement transformations.
Definition: CentroidUtilities.h:168
lsst::meas::base::FlagDefinition::name
std::string name
Definition: FlagHandler.h:52
lsst::meas::base::CentroidResult
A reusable struct for centroid measurements.
Definition: CentroidUtilities.h:41
lsst::afw::image::Exposure::getMaskedImage
MaskedImageT getMaskedImage()
Return the MaskedImage.
Definition: Exposure.h:228
std::sqrt
T sqrt(T... args)
lsst::afw::table::Schema::find
SchemaItem< T > find(std::string const &name) const
Find a SchemaItem in the Schema by name.
Definition: Schema.cc:656
lsst::meas::base::FlagDefinitionList
vector-type utility class to build a collection of FlagDefinitions
Definition: FlagHandler.h:60
sigma2
table::Key< double > sigma2
Definition: FunctionLibrary.cc:21
lsst::meas::base::SdssCentroidAlgorithm::measure
virtual void measure(afw::table::SourceRecord &measRecord, afw::image::Exposure< float > const &exposure) const
Called to measure a single child source in an image.
Definition: SdssCentroid.cc:411
lsst::afw::image::MaskedImage
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:73
lsst::meas::base::SdssCentroidControl::wfac
double wfac
"fiddle factor for adjusting the binning" ;
Definition: SdssCentroid.h:52
lsst::afw::math::MEAN
@ MEAN
estimate sample mean
Definition: Statistics.h:67
Angle.h
lsst::meas::base::SdssCentroidAlgorithm::SdssCentroidAlgorithm
SdssCentroidAlgorithm(Control const &ctrl, std::string const &name, afw::table::Schema &schema)
Definition: SdssCentroid.cc:403
lsst::meas::base::FlagHandler
Utility class for handling flag fields that indicate the failure modes of an algorithm.
Definition: FlagHandler.h:148
x
double x
Definition: ChebyshevBoundedField.cc:277
lsst::afw::table::SchemaMapper
A mapping between the keys of two Schemas, used to copy data between them.
Definition: SchemaMapper.h:21
lsst::afw::table::Key< afw::table::Flag >
lsst::meas::base::SdssCentroidAlgorithm::NOT_AT_MAXIMUM
static FlagDefinition const NOT_AT_MAXIMUM
Definition: SdssCentroid.h:77
lsst::afw::image::Y
@ Y
Definition: ImageUtils.h:36
result
py::object result
Definition: _schema.cc:429
Source.h
lsst
A base class for image defects.
Definition: imageAlgorithm.dox:1
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
lsst::afw::math::convolve
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.
Definition: ConvolveImage.cc:190
lsst::meas::base::FlagDefinitionList::size
std::size_t size() const
return the current size (number of defined elements) of the collection
Definition: FlagHandler.h:125
lsst::meas::base::FlagDefinition
Simple class used to define and document flags The name and doc constitute the identity of the FlagDe...
Definition: FlagHandler.h:40
lsst::meas::base::SdssCentroidAlgorithm::ALMOST_NO_SECOND_DERIVATIVE
static FlagDefinition const ALMOST_NO_SECOND_DERIVATIVE
Definition: SdssCentroid.h:76
std::exp
T exp(T... args)
lsst::meas::base::SdssCentroidControl::binmax
int binmax
"maximum allowed binning" ;
Definition: SdssCentroid.h:50
PTR
#define PTR(...)
Definition: base.h:41
key
Key< U > key
Definition: Schema.cc:281
lsst::geom::Point< double, 2 >
lsst::meas::base::SdssCentroidControl
A C++ control class to handle SdssCentroidAlgorithm's configuration.
Definition: SdssCentroid.h:48
lsst::geom::Box2I
An integer coordinate rectangle.
Definition: Box.h:55
mapper
SchemaMapper * mapper
Definition: SchemaMapper.cc:78
lsst::meas::base::SdssCentroidControl::peakMin
double peakMin
"if the peak's less than this insist on binning at least once" ;
Definition: SdssCentroid.h:51
lsst::afw::table::BaseRecord::getSchema
Schema getSchema() const
Return the Schema that holds this record's fields and keys.
Definition: BaseRecord.h:80
std::size_t
M_PI
#define M_PI
Definition: ListMatch.cc:31
std::make_pair
T make_pair(T... args)
lsst::meas::base::SdssCentroidTransform::SdssCentroidTransform
SdssCentroidTransform(Control const &ctrl, std::string const &name, afw::table::SchemaMapper &mapper)
Definition: SdssCentroid.cc:508
lsst::pex::exceptions::Exception
Provides consistent interface for LSST exceptions.
Definition: Exception.h:107
lsst::afw::math::binImage
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:43
lsst::afw::detection::Psf
A polymorphic base class for representing an image's Point Spread Function.
Definition: Psf.h:76
lsst::afw::table::BaseRecord::set
void set(Key< T > const &key, U const &value)
Set value of a field for the given key.
Definition: BaseRecord.h:164
lsst::meas::base::SdssCentroidAlgorithm::EDGE
static FlagDefinition const EDGE
Definition: SdssCentroid.h:74
lsst::meas::base::SdssCentroidAlgorithm::NO_SECOND_DERIVATIVE
static FlagDefinition const NO_SECOND_DERIVATIVE
Definition: SdssCentroid.h:75
ConvolveImage.h
lsst::geom::Extent< int, 2 >
std::numeric_limits
lsst::meas::base::SdssCentroidAlgorithm::fail
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.
Definition: SdssCentroid.cc:504
lsst::meas::base::SdssCentroidAlgorithm::getFlagDefinitions
static FlagDefinitionList const & getFlagDefinitions()
Definition: SdssCentroid.cc:52
exceptions.h
bbox
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
std::pow
T pow(T... args)
Generated on Wed Jun 24 2020 18:10:10 for LSSTApplications by   doxygen 1.8.18