LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
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/afw/detection/Psf.h"
28 #include "lsst/pex/exceptions.h"
29 #include "lsst/pex/logging/Trace.h"
32 #include "lsst/afw/geom/Angle.h"
33 #include "lsst/afw/table/Source.h"
35 
36 
37 namespace lsst { namespace meas { namespace base {
38 
39 namespace {
40 
41 /************************************************************************************************************/
42 
43 float const AMPAST4 = 1.33; // amplitude of `4th order' corr compared to theory
44 
45 /*
46  * Do the Gaussian quartic interpolation for the position
47  * of the maximum for three equally spaced values vm,v0,vp, assumed to be at
48  * abscissae -1,0,1; the answer is returned as *cen
49  *
50  * Return 0 is all is well, otherwise 1
51  */
52 static int inter4(float vm, float v0, float vp, float *cen) {
53  float const sp = v0 - vp;
54  float const sm = v0 - vm;
55  float const d2 = sp + sm;
56  float const s = 0.5*(vp - vm);
57 
58  if (d2 <= 0.0f || v0 <= 0.0f) {
59  return(1);
60  }
61 
62  *cen = s/d2*(1.0 + AMPAST4*sp*sm/(d2*v0));
63 
64  return fabs(*cen) < 1 ? 0 : 1;
65 }
66 
67 /*****************************************************************************/
68 /*
69  * Calculate error in centroid
70  */
71 float astrom_errors(float skyVar, // variance of pixels at the sky level
72  float sourceVar, // variance in peak due to excess counts over sky
73  float A, // abs(peak value in raw image)
74  float tau2, // Object is N(0,tau2)
75  float As, // abs(peak value in smoothed image)
76  float s, // slope across central pixel
77  float d, // curvature at central pixel
78  float sigma, // width of smoothing filter
79  int quarticBad) { // was quartic estimate bad?
80 
81  float const k = quarticBad ? 0 : AMPAST4; /* quartic correction coeff */
82  float const sigma2 = sigma*sigma; /* == sigma^2 */
83  float sVar, dVar; /* variances of s and d */
84  float xVar; /* variance of centroid, x */
85 
86  if (fabs(As) < std::numeric_limits<float>::min() ||
87  fabs(d) < std::numeric_limits<float>::min()) {
88  return(1e3);
89  }
90 
91  if (sigma <= 0) { /* no smoothing; no covariance */
92  sVar = 0.5*skyVar; /* due to sky */
93  dVar = 6*skyVar;
94 
95  sVar += 0.5*sourceVar*exp(-1/(2*tau2));
96  dVar += sourceVar*(4*exp(-1/(2*tau2)) + 2*exp(-1/(2*tau2)));
97  } else { /* smoothed */
98  sVar = skyVar/(8*lsst::afw::geom::PI*sigma2)*(1 - exp(-1/sigma2));
99  dVar = skyVar/(2*lsst::afw::geom::PI*sigma2)*(3 - 4*exp(-1/(4*sigma2)) + exp(-1/sigma2));
100 
101  sVar += sourceVar/(12*lsst::afw::geom::PI*sigma2)*(exp(-1/(3*sigma2)) - exp(-1/sigma2));
102  dVar += sourceVar/(3*lsst::afw::geom::PI*sigma2)*(2 - 3*exp(-1/(3*sigma2)) + exp(-1/sigma2));
103  }
104 
105  xVar = sVar*pow(1/d + k/(4*As)*(1 - 12*s*s/(d*d)), 2) +
106  dVar*pow(s/(d*d) - k/(4*As)*8*s*s/(d*d*d), 2);
107 
108  return(xVar >= 0 ? sqrt(xVar) : NAN);
109 }
110 
111 /************************************************************************************************************/
112 /*
113  * Estimate the position of an object, assuming we know that it's approximately the size of the PSF
114  */
115 
116 template<typename ImageXy_locatorT, typename VarImageXy_locatorT>
117 void doMeasureCentroidImpl(double *xCenter, // output; x-position of object
118  double *dxc, // output; error in xCenter
119  double *yCenter, // output; y-position of object
120  double *dyc, // output; error in yCenter
121  double *sizeX2, double *sizeY2, // output; object widths^2 in x and y directions
122  ImageXy_locatorT im, // Locator for the pixel values
123  VarImageXy_locatorT vim, // Locator for the image containing the variance
124  double smoothingSigma, // Gaussian sigma of already-applied smoothing filter
125  FlagHandler flagHandler
126  )
127 {
128  /*
129  * find a first quadratic estimate
130  */
131  double const d2x = 2*im(0, 0) - im(-1, 0) - im(1, 0);
132  double const d2y = 2*im(0, 0) - im( 0, -1) - im(0, 1);
133  double const sx = 0.5*(im(1, 0) - im(-1, 0));
134  double const sy = 0.5*(im(0, 1) - im( 0, -1));
135 
136  if (d2x == 0.0 || d2y == 0.0) {
137  throw LSST_EXCEPT(
138  MeasurementError,
139  flagHandler.getDefinition(SdssCentroidAlgorithm::NO_SECOND_DERIVATIVE).doc,
141  );
142  }
143  if (d2x < 0.0 || d2y < 0.0) {
144  throw LSST_EXCEPT(
145  MeasurementError,
146  flagHandler.getDefinition(SdssCentroidAlgorithm::NOT_AT_MAXIMUM).doc +
147  (boost::format(": d2I/dx2, d2I/dy2 = %g %g") % d2x % d2y).str(),
149  );
150  }
151 
152  double const dx0 = sx/d2x;
153  double const dy0 = sy/d2y; // first guess
154 
155  if (fabs(dx0) > 10.0 || fabs(dy0) > 10.0) {
156  throw LSST_EXCEPT(
157  MeasurementError,
158  flagHandler.getDefinition(SdssCentroidAlgorithm::ALMOST_NO_SECOND_DERIVATIVE).doc +
159  (boost::format(": sx, d2x, sy, d2y = %f %f %f %f") % sx % d2x % sy % d2y).str(),
161  );
162  }
163 
164  double vpk = im(0, 0) + 0.5*(sx*dx0 + sy*dy0); // height of peak in image
165  if (vpk < 0) {
166  vpk = -vpk;
167  }
168 /*
169  * now evaluate maxima on stripes
170  */
171  float m0x = 0, m1x = 0, m2x = 0;
172  float m0y = 0, m1y = 0, m2y = 0;
173 
174  int quarticBad = 0;
175  quarticBad += inter4(im(-1, -1), im( 0, -1), im( 1, -1), &m0x);
176  quarticBad += inter4(im(-1, 0), im( 0, 0), im( 1, 0), &m1x);
177  quarticBad += inter4(im(-1, 1), im( 0, 1), im( 1, 1), &m2x);
178 
179  quarticBad += inter4(im(-1, -1), im(-1, 0), im(-1, 1), &m0y);
180  quarticBad += inter4(im( 0, -1), im( 0, 0), im( 0, 1), &m1y);
181  quarticBad += inter4(im( 1, -1), im( 1, 0), im( 1, 1), &m2y);
182 
183  double xc, yc; // position of maximum
184  double sigmaX2, sigmaY2; // widths^2 in x and y of smoothed object
185 
186  if (quarticBad) { // >= 1 quartic interpolator is bad
187  xc = dx0;
188  yc = dy0;
189  sigmaX2 = vpk/d2x; // widths^2 in x
190  sigmaY2 = vpk/d2y; // and y
191  } else {
192  double const smx = 0.5*(m2x - m0x);
193  double const smy = 0.5*(m2y - m0y);
194  double const dm2x = m1x - 0.5*(m0x + m2x);
195  double const dm2y = m1y - 0.5*(m0y + m2y);
196  double const dx = m1x + dy0*(smx - dy0*dm2x); // first quartic approx
197  double const dy = m1y + dx0*(smy - dx0*dm2y);
198  double const dx4 = m1x + dy*(smx - dy*dm2x); // second quartic approx
199  double const dy4 = m1y + dx*(smy - dx*dm2y);
200 
201  xc = dx4;
202  yc = dy4;
203  sigmaX2 = vpk/d2x - (1 + 6*dx0*dx0)/4; // widths^2 in x
204  sigmaY2 = vpk/d2y - (1 + 6*dy0*dy0)/4; // and y
205  }
206  /*
207  * Now for the errors.
208  */
209  float tauX2 = sigmaX2; // width^2 of _un_ smoothed object
210  float tauY2 = sigmaY2;
211  tauX2 -= smoothingSigma*smoothingSigma; // correct for smoothing
212  tauY2 -= smoothingSigma*smoothingSigma;
213 
214  if (tauX2 <= smoothingSigma*smoothingSigma) { // problem; sigmaX2 must be bad
215  tauX2 = smoothingSigma*smoothingSigma;
216  }
217  if (tauY2 <= smoothingSigma*smoothingSigma) { // sigmaY2 must be bad
218  tauY2 = smoothingSigma*smoothingSigma;
219  }
220 
221  float const skyVar = (vim(-1, -1) + vim( 0, -1) + vim( 1, -1) +
222  vim(-1, 0) + vim( 1, 0) +
223  vim(-1, 1) + vim( 0, 1) + vim( 1, 1))/8.0; // Variance in sky
224  float const sourceVar = vim(0, 0); // extra variance of peak due to its photons
225  float const A = vpk*sqrt((sigmaX2/tauX2)*(sigmaY2/tauY2)); // peak of Unsmoothed object
226 
227  *xCenter = xc;
228  *yCenter = yc;
229 
230  *dxc = astrom_errors(skyVar, sourceVar, A, tauX2, vpk, sx, d2x, fabs(smoothingSigma), quarticBad);
231  *dyc = astrom_errors(skyVar, sourceVar, A, tauY2, vpk, sy, d2y, fabs(smoothingSigma), quarticBad);
232 
233  *sizeX2 = tauX2; // return the estimates of the (object size)^2
234  *sizeY2 = tauY2;
235 }
236 
237 template<typename MaskedImageXy_locatorT>
238 void doMeasureCentroidImpl(double *xCenter, // output; x-position of object
239  double *dxc, // output; error in xCenter
240  double *yCenter, // output; y-position of object
241  double *dyc, // output; error in yCenter
242  double *sizeX2, double *sizeY2, // output; object widths^2 in x and y directions
243  double *peakVal, // output; peak of object
244  MaskedImageXy_locatorT mim, // Locator for the pixel values
245  double smoothingSigma, // Gaussian sigma of already-applied smoothing filter
246  FlagHandler flagHandler
247  )
248 {
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) {
258  throw LSST_EXCEPT(
259  MeasurementError,
260  flagHandler.getDefinition(SdssCentroidAlgorithm::NO_SECOND_DERIVATIVE).doc,
262  );
263  }
264  if (d2x < 0.0 || d2y < 0.0) {
265  throw LSST_EXCEPT(
266  MeasurementError,
267  flagHandler.getDefinition(SdssCentroidAlgorithm::NOT_AT_MAXIMUM).doc +
268  (boost::format(": d2I/dx2, d2I/dy2 = %g %g") % d2x % d2y).str(),
270  );
271  }
272 
273  double const dx0 = sx/d2x;
274  double const dy0 = sy/d2y; // first guess
275 
276  if (fabs(dx0) > 10.0 || fabs(dy0) > 10.0) {
277  throw LSST_EXCEPT(
278  MeasurementError,
279  flagHandler.getDefinition(SdssCentroidAlgorithm::NO_SECOND_DERIVATIVE).doc +
280  (boost::format(": sx, d2x, sy, d2y = %f %f %f %f") % sx % d2x % sy % d2y).str(),
282  );
283  }
284 
285  double vpk = mim.image(0, 0) + 0.5*(sx*dx0 + sy*dy0); // height of peak in image
286  if (vpk < 0) {
287  vpk = -vpk;
288  }
289 /*
290  * now evaluate maxima on stripes
291  */
292  float m0x = 0, m1x = 0, m2x = 0;
293  float m0y = 0, m1y = 0, m2y = 0;
294 
295  int quarticBad = 0;
296  quarticBad += inter4(mim.image(-1, -1), mim.image( 0, -1), mim.image( 1, -1), &m0x);
297  quarticBad += inter4(mim.image(-1, 0), mim.image( 0, 0), mim.image( 1, 0), &m1x);
298  quarticBad += inter4(mim.image(-1, 1), mim.image( 0, 1), mim.image( 1, 1), &m2x);
299 
300  quarticBad += inter4(mim.image(-1, -1), mim.image(-1, 0), mim.image(-1, 1), &m0y);
301  quarticBad += inter4(mim.image( 0, -1), mim.image( 0, 0), mim.image( 0, 1), &m1y);
302  quarticBad += inter4(mim.image( 1, -1), mim.image( 1, 0), mim.image( 1, 1), &m2y);
303 
304  double xc, yc; // position of maximum
305  double sigmaX2, sigmaY2; // widths^2 in x and y of smoothed object
306 
307  if (quarticBad) { // >= 1 quartic interpolator is bad
308  xc = dx0;
309  yc = dy0;
310  sigmaX2 = vpk/d2x; // widths^2 in x
311  sigmaY2 = vpk/d2y; // and y
312  } else {
313  double const smx = 0.5*(m2x - m0x);
314  double const smy = 0.5*(m2y - m0y);
315  double const dm2x = m1x - 0.5*(m0x + m2x);
316  double const dm2y = m1y - 0.5*(m0y + m2y);
317  double const dx = m1x + dy0*(smx - dy0*dm2x); // first quartic approx
318  double const dy = m1y + dx0*(smy - dx0*dm2y);
319  double const dx4 = m1x + dy*(smx - dy*dm2x); // second quartic approx
320  double const dy4 = m1y + dx*(smy - dx*dm2y);
321 
322  xc = dx4;
323  yc = dy4;
324  sigmaX2 = vpk/d2x - (1 + 6*dx0*dx0)/4; // widths^2 in x
325  sigmaY2 = vpk/d2y - (1 + 6*dy0*dy0)/4; // and y
326  }
327  /*
328  * Now for the errors.
329  */
330  float tauX2 = sigmaX2; // width^2 of _un_ smoothed object
331  float tauY2 = sigmaY2;
332  tauX2 -= smoothingSigma*smoothingSigma; // correct for smoothing
333  tauY2 -= smoothingSigma*smoothingSigma;
334 
335  if (tauX2 <= smoothingSigma*smoothingSigma) { // problem; sigmaX2 must be bad
336  tauX2 = smoothingSigma*smoothingSigma;
337  }
338  if (tauY2 <= smoothingSigma*smoothingSigma) { // sigmaY2 must be bad
339  tauY2 = smoothingSigma*smoothingSigma;
340  }
341 
342  float const skyVar = (mim.variance(-1, -1) + mim.variance( 0, -1) + mim.variance( 1, -1) +
343  mim.variance(-1, 0) + mim.variance( 1, 0) +
344  mim.variance(-1, 1) + mim.variance( 0, 1) + mim.variance( 1, 1)
345  )/8.0; // Variance in sky
346  float const sourceVar = mim.variance(0, 0); // extra variance of peak due to its photons
347  float const A = vpk*sqrt((sigmaX2/tauX2)*(sigmaY2/tauY2)); // peak of Unsmoothed object
348 
349  *xCenter = xc;
350  *yCenter = yc;
351 
352  *dxc = astrom_errors(skyVar, sourceVar, A, tauX2, vpk, sx, d2x, fabs(smoothingSigma), quarticBad);
353  *dyc = astrom_errors(skyVar, sourceVar, A, tauY2, vpk, sy, d2y, fabs(smoothingSigma), quarticBad);
354 
355  *sizeX2 = tauX2; // return the estimates of the (object size)^2
356  *sizeY2 = tauY2;
357 
358  *peakVal = vpk;
359 }
360 
361 template<typename MaskedImageT>
362 std::pair<MaskedImageT, double>
363 smoothAndBinImage(CONST_PTR(lsst::afw::detection::Psf) psf,
364  int const x, const int y,
365  MaskedImageT const& mimage,
366  int binX, int binY,
367  FlagHandler _flagHandler)
368 {
369  lsst::afw::geom::Point2D const center(x + mimage.getX0(), y + mimage.getY0());
370  lsst::afw::geom::ellipses::Quadrupole const& shape = psf->computeShape(center);
371  double const smoothingSigma = shape.getDeterminantRadius();
372 #if 0
373  double const nEffective = psf->computeEffectiveArea(); // not implemented yet (#2821)
374 #else
375  double const nEffective = 4*M_PI*smoothingSigma*smoothingSigma; // correct for a Gaussian
376 #endif
377 
379  int const kWidth = kernel->getWidth();
380  int const kHeight = kernel->getHeight();
381 
382  lsst::afw::geom::BoxI bbox(lsst::afw::geom::Point2I(x - binX*(2 + kWidth/2), y - binY*(2 + kHeight/2)),
383  lsst::afw::geom::ExtentI(binX*(3 + kWidth + 1), binY*(3 + kHeight + 1)));
384 
385  // image to smooth, a shallow copy
386  PTR(MaskedImageT) subImage;
387  try {
388  subImage.reset(new MaskedImageT(mimage, bbox, lsst::afw::image::LOCAL));
389  } catch (pex::exceptions::LengthError & err) {
390  throw LSST_EXCEPT(
391  MeasurementError,
392  _flagHandler.getDefinition(SdssCentroidAlgorithm::EDGE).doc,
394  );
395  }
396  PTR(MaskedImageT) binnedImage = lsst::afw::math::binImage(*subImage, binX, binY, lsst::afw::math::MEAN);
397  binnedImage->setXY0(subImage->getXY0());
398  // image to smooth into, a deep copy.
399  MaskedImageT smoothedImage = MaskedImageT(*binnedImage, true);
400  assert(smoothedImage.getWidth()/2 == kWidth/2 + 2); // assumed by the code that uses smoothedImage
401  assert(smoothedImage.getHeight()/2 == kHeight/2 + 2);
402 
403  lsst::afw::math::convolve(smoothedImage, *binnedImage, *kernel, lsst::afw::math::ConvolutionControl());
404  *smoothedImage.getVariance() *= binX*binY*nEffective; // We want the per-pixel variance, so undo the
405  // effects of binning and smoothing
406 
407  return std::make_pair(smoothedImage, smoothingSigma);
408 }
409 
410 boost::array<FlagDefinition,SdssCentroidAlgorithm::N_FLAGS> const & getFlagDefinitions() {
411  static boost::array<FlagDefinition,SdssCentroidAlgorithm::N_FLAGS> const flagDefs = {{
412  {"flag", "general failure flag, set if anything went wrong"},
413  {"flag_edge", "Object too close to edge"},
414  {"flag_noSecondDerivative", "Vanishing second derivative"},
415  {"flag_almostNoSecondDerivative", "Almost vanishing second derivative"},
416  {"flag_notAtMaximum", "Object is not at a maximum"}
417  }};
418  return flagDefs;
419 }
420 
421 } // end anonymous namespace
422 
424  Control const & ctrl,
425  std::string const & name,
427 ) : _ctrl(ctrl),
428  _centroidKey(
429  CentroidResultKey::addFields(schema, name, "centroid from Sdss Centroid algorithm", SIGMA_ONLY)
430  ),
431  _centroidExtractor(schema, name, true)
432 {
433  _flagHandler = FlagHandler::addFields(schema, name,
434  getFlagDefinitions().begin(), getFlagDefinitions().end());
435 }
437  afw::table::SourceRecord & measRecord,
438  afw::image::Exposure<float> const & exposure
439 ) const {
440 
441  // get our current best guess about the centroid: either a centroider measurement or peak.
442  afw::geom::Point2D center = _centroidExtractor(measRecord, _flagHandler);
443  CentroidResult result;
444  result.x = center.getX();
445  result.y = center.getY();
446  measRecord.set(_centroidKey, result); // better than NaN
447 
448  typedef afw::image::Exposure<float>::MaskedImageT MaskedImageT;
449  typedef MaskedImageT::Image ImageT;
450  typedef MaskedImageT::Variance VarianceT;
451 
452  MaskedImageT const& mimage = exposure.getMaskedImage();
453  ImageT const& image = *mimage.getImage();
454  CONST_PTR(lsst::afw::detection::Psf) psf = exposure.getPsf();
455 
456  int const x = image.positionToIndex(center.getX(), lsst::afw::image::X).first;
457  int const y = image.positionToIndex(center.getY(), lsst::afw::image::Y).first;
458 
459  if (!image.getBBox().contains(lsst::afw::geom::Extent2I(x,y) + image.getXY0())) {
460  throw LSST_EXCEPT(
462  _flagHandler.getDefinition(EDGE).doc,
463  EDGE
464  );
465  }
466 
467  // Algorithm uses a least-squares fit (implemented via a convolution) to a symmetrized PSF model.
468  // If you don't have a Psf, you need to use another centroider, such as GaussianCentroider.
469  if (!psf) {
470  throw LSST_EXCEPT(
471  FatalAlgorithmError,
472  "SdssCentroid algorithm requires a Psf with every exposure"
473  );
474  }
475 
476  int binX = 1;
477  int binY = 1;
478  double xc=0., yc=0., dxc=0., dyc=0.; // estimated centre and error therein
479  for(int binsize = 1; binsize <= _ctrl.binmax; binsize *= 2) {
480  std::pair<MaskedImageT, double> result = smoothAndBinImage(psf, x, y, mimage, binX, binY, _flagHandler);
481  MaskedImageT const smoothedImage = result.first;
482  double const smoothingSigma = result.second;
483 
484  MaskedImageT::xy_locator mim = smoothedImage.xy_at(smoothedImage.getWidth()/2,
485  smoothedImage.getHeight()/2);
486 
487  double sizeX2, sizeY2; // object widths^2 in x and y directions
488  double peakVal; // peak intensity in image
489 
490  doMeasureCentroidImpl(&xc, &dxc, &yc, &dyc, &sizeX2, &sizeY2, &peakVal, mim,
491  smoothingSigma, _flagHandler);
492 
493  if(binsize > 1) {
494  // dilate from the lower left corner of central pixel
495  xc = (xc + 0.5)*binX - 0.5;
496  dxc *= binX;
497  sizeX2 *= binX*binX;
498 
499  yc = (yc + 0.5)*binY - 0.5;
500  dyc *= binY;
501  sizeY2 *= binY*binY;
502  }
503 
504  xc += x; // xc, yc are measured relative to pixel (x, y)
505  yc += y;
506 
507  double const fac = _ctrl.wfac*(1 + smoothingSigma*smoothingSigma);
508  double const facX2 = fac*binX*binX;
509  double const facY2 = fac*binY*binY;
510 
511  if (sizeX2 < facX2 && ::pow(xc - x, 2) < facX2 &&
512  sizeY2 < facY2 && ::pow(yc - y, 2) < facY2) {
513  if (binsize > 1 || _ctrl.peakMin < 0.0 || peakVal > _ctrl.peakMin) {
514  break;
515  }
516  }
517 
518  if (sizeX2 >= facX2 || ::pow(xc - x, 2) >= facX2) {
519  binX *= 2;
520  }
521  if (sizeY2 >= facY2 || ::pow(yc - y, 2) >= facY2) {
522  binY *= 2;
523  }
524  }
525  result.x = lsst::afw::image::indexToPosition(xc + image.getX0());
526  result.y = lsst::afw::image::indexToPosition(yc + image.getY0());
527 
528  result.xSigma = sqrt(dxc*dxc);
529  result.ySigma = sqrt(dyc*dyc);
530  measRecord.set(_centroidKey, result);
531  _flagHandler.setValue(measRecord, FAILURE, false);
532 
533 }
534 
535 
537  _flagHandler.handleFailure(measRecord, error);
538 }
539 
541  Control const & ctrl,
542  std::string const & name,
543  afw::table::SchemaMapper & mapper
544 ) :
545  CentroidTransform{name, mapper}
546 {
547  for (auto flag = getFlagDefinitions().begin() + 1; flag < getFlagDefinitions().end(); ++flag) {
548  mapper.addMapping(mapper.getInputSchema().find<afw::table::Flag>(
549  mapper.getInputSchema().join(name, flag->name)).key);
550  }
551 }
552 
553 }}} // end namespace lsst::meas::base
554 
int y
An ellipse core with quadrupole moments as parameters.
Definition: Quadrupole.h:45
Defines the fields and offsets for a table.
Definition: Schema.h:46
geom::ellipses::Quadrupole computeShape(geom::Point2D position=makeNullPoint(), image::Color color=image::Color()) const
Compute the ellipse corresponding to the second moments of the Psf.
MaskedImageT getMaskedImage()
Return the MaskedImage.
Definition: Exposure.h:150
SdssCentroidAlgorithm(Control const &ctrl, std::string const &name, afw::table::Schema &schema)
boost::shared_ptr< math::Kernel const > getLocalKernel(geom::Point2D position=makeNullPoint(), image::Color color=image::Color()) const
Return a FixedKernel corresponding to the Psf image at the given point.
ErrElement ySigma
1-Sigma uncertainty on y (sqrt of variance)
estimate sample mean
Definition: Statistics.h:67
#define PTR(...)
Definition: base.h:41
double indexToPosition(double ind)
Convert image index to image position.
Definition: ImageUtils.h:54
table::Key< std::string > name
Definition: ApCorrMap.cc:71
Eigen matrix objects that present a view into an ndarray::Array.
boost::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:48
boost::shared_ptr< Kernel const > ConstPtr
Definition: Kernel.h:142
A mapping between the keys of two Schemas, used to copy data between them.
Definition: SchemaMapper.h:19
Only the diagonal elements of the covariance matrix are provided.
Definition: constants.h:43
A reusable struct for centroid measurements.
#define CONST_PTR(...)
Definition: base.h:47
definition of the Trace messaging facilities
ErrElement xSigma
1-Sigma uncertainty on x (sqrt of variance)
CentroidElement x
x (column) coordinate of the measured position
A C++ control class to handle SdssCentroidAlgorithm&#39;s configuration.
Definition: SdssCentroid.h:46
int binmax
&quot;maximum allowed binning&quot; ;
Definition: SdssCentroid.h:49
Exception to be thrown when a measurement algorithm experiences a known failure mode.
Definition: exceptions.h:48
An integer coordinate rectangle.
Definition: Box.h:53
Key< T > addMapping(Key< T > const &inputKey, bool doReplace=false)
Add a new field to the output Schema that is a copy of a field in the input Schema.
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:43
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
def error
Definition: log.py:108
SafeCentroidExtractor _centroidExtractor
Definition: SdssCentroid.h:98
double getDeterminantRadius() const
Return the radius defined as the 4th root of the determinant of the quadrupole matrix.
tbl::Schema schema
boost::shared_ptr< lsst::afw::detection::Psf > getPsf()
Return the Exposure&#39;s Psf object.
Definition: Exposure.h:224
int x
Convolve and convolveAtAPoint functions for Image and Kernel.
double peakMin
&quot;if the peak&#39;s less than this insist on binning at least once&quot; ;
Definition: SdssCentroid.h:50
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
double wfac
&quot;fiddle factor for adjusting the binning&quot; ;
Definition: SdssCentroid.h:51
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...
virtual void measure(afw::table::SourceRecord &measRecord, afw::image::Exposure< float > const &exposure) const
void set(Key< T > const &key, U const &value)
Set value of a field for the given key.
Definition: BaseRecord.h:136
A FunctorKey for CentroidResult.
CentroidElement y
y (row) coordinate of the measured position
Record class that contains measurements made on a single exposure.
Definition: Source.h:81
double const PI
The ratio of a circle&#39;s circumference to diameter.
Definition: Angle.h:18
A polymorphic base class for representing an image&#39;s Point Spread Function.
Definition: Psf.h:68
afw::table::Key< double > sigma2
virtual void fail(afw::table::SourceRecord &measRecord, MeasurementError *error=NULL) const
static FlagHandler addFields(afw::table::Schema &schema, std::string const &prefix, FlagDefinition const *begin, FlagDefinition const *end)
Definition: FlagHandler.cc:28
Include files required for standard LSST Exception handling.
SdssCentroidTransform(Control const &ctrl, std::string const &name, afw::table::SchemaMapper &mapper)