LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
SdssShape.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
2 /*
3  * LSST Data Management System
4  * Copyright 2008-2016 AURA/LSST.
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 <https://www.lsstcorp.org/LegalNotices/>.
22  */
23 
24 #include <cmath>
25 #include <tuple>
26 
27 #include "boost/tuple/tuple.hpp"
28 #include "Eigen/LU"
29 #include "lsst/geom/Angle.h"
30 #include "lsst/geom/Box.h"
32 #include "lsst/afw/image.h"
33 #include "lsst/afw/detection/Psf.h"
34 #include "lsst/afw/geom/ellipses.h"
35 #include "lsst/afw/table/Source.h"
38 
39 namespace lsst {
40 namespace meas {
41 namespace base {
42 namespace {
43 FlagDefinitionList flagDefinitions;
44 } // namespace
45 
46 FlagDefinition const SdssShapeAlgorithm::FAILURE = flagDefinitions.addFailureFlag();
47 FlagDefinition const SdssShapeAlgorithm::UNWEIGHTED_BAD =
48  flagDefinitions.add("flag_unweightedBad", "Both weighted and unweighted moments were invalid");
49 FlagDefinition const SdssShapeAlgorithm::UNWEIGHTED = flagDefinitions.add(
50  "flag_unweighted", "Weighted moments converged to an invalid value; using unweighted moments");
51 FlagDefinition const SdssShapeAlgorithm::SHIFT =
52  flagDefinitions.add("flag_shift", "centroid shifted by more than the maximum allowed amount");
53 FlagDefinition const SdssShapeAlgorithm::MAXITER =
54  flagDefinitions.add("flag_maxIter", "Too many iterations in adaptive moments");
55 FlagDefinition const SdssShapeAlgorithm::PSF_SHAPE_BAD =
56  flagDefinitions.add("flag_psf", "Failure in measuring PSF model shape");
57 
58 FlagDefinitionList const &SdssShapeAlgorithm::getFlagDefinitions() { return flagDefinitions; }
59 
60 namespace { // anonymous
61 
62 typedef Eigen::Matrix<double, 4, 4, Eigen::DontAlign> Matrix4d;
63 
64 /*****************************************************************************/
65 /*
66  * Error analysis, courtesy of David Johnston, University of Chicago
67  */
68 /*
69  * This function takes the 4 Gaussian parameters A, sigmaXXW and the
70  * sky variance and fills in the Fisher matrix from the least squares fit.
71  *
72  * Following "Numerical Recipes in C" section 15.5, it ignores the 2nd
73  * derivative parts and so the fisher matrix is just a function of these
74  * best fit model parameters. The components are calculated analytically.
75  */
76 Matrix4d calc_fisher(SdssShapeResult const &shape, // the Shape that we want the the Fisher matrix for
77  float bkgd_var // background variance level for object
78  ) {
79  float const A = shape.instFlux; // amplitude; will be converted to instFlux later
80  float const sigma11W = shape.xx;
81  float const sigma12W = shape.xy;
82  float const sigma22W = shape.yy;
83 
84  double const D = sigma11W * sigma22W - sigma12W * sigma12W;
85 
87  throw LSST_EXCEPT(pex::exceptions::DomainError, "Determinant is too small calculating Fisher matrix");
88  }
89  /*
90  * a normalization factor
91  */
92  if (bkgd_var <= 0.0) {
94  (boost::format("Background variance must be positive (saw %g)") % bkgd_var).str());
95  }
96  double const F = geom::PI * sqrt(D) / bkgd_var;
97  /*
98  * Calculate the 10 independent elements of the 4x4 Fisher matrix
99  */
100  Matrix4d fisher;
101 
102  double fac = F * A / (4.0 * D);
103  fisher(0, 0) = F;
104  fisher(0, 1) = fac * sigma22W;
105  fisher(1, 0) = fisher(0, 1);
106  fisher(0, 2) = fac * sigma11W;
107  fisher(2, 0) = fisher(0, 2);
108  fisher(0, 3) = -fac * 2 * sigma12W;
109  fisher(3, 0) = fisher(0, 3);
110 
111  fac = 3.0 * F * A * A / (16.0 * D * D);
112  fisher(1, 1) = fac * sigma22W * sigma22W;
113  fisher(2, 2) = fac * sigma11W * sigma11W;
114  fisher(3, 3) = fac * 4.0 * (sigma12W * sigma12W + D / 3.0);
115 
116  fisher(1, 2) = fisher(3, 3) / 4.0;
117  fisher(2, 1) = fisher(1, 2);
118  fisher(1, 3) = fac * (-2 * sigma22W * sigma12W);
119  fisher(3, 1) = fisher(1, 3);
120  fisher(2, 3) = fac * (-2 * sigma11W * sigma12W);
121  fisher(3, 2) = fisher(2, 3);
122 
123  return fisher;
124 }
125 //
126 // Here's a class to allow us to get the Image and variance from an Image or MaskedImage
127 //
128 template <typename ImageT> // general case
129 struct ImageAdaptor {
130  typedef ImageT Image;
131 
132  static bool const hasVariance = false;
133 
134  Image const &getImage(ImageT const &image) const { return image; }
135 
136  double getVariance(ImageT const &, int, int) { return std::numeric_limits<double>::quiet_NaN(); }
137 };
138 
139 template <typename T> // specialise to a MaskedImage
140 struct ImageAdaptor<afw::image::MaskedImage<T> > {
141  typedef typename afw::image::MaskedImage<T>::Image Image;
142 
143  static bool const hasVariance = true;
144 
145  Image const &getImage(afw::image::MaskedImage<T> const &mimage) const { return *mimage.getImage(); }
146 
147  double getVariance(afw::image::MaskedImage<T> const &mimage, int ix, int iy) {
148  return mimage.at(ix, iy).variance();
149  }
150 };
151 
153 std::tuple<std::pair<bool, double>, double, double, double> getWeights(double sigma11, double sigma12,
154  double sigma22) {
155  double const NaN = std::numeric_limits<double>::quiet_NaN();
156  if (std::isnan(sigma11) || std::isnan(sigma12) || std::isnan(sigma22)) {
157  return std::make_tuple(std::make_pair(false, NaN), NaN, NaN, NaN);
158  }
159  double const det = sigma11 * sigma22 - sigma12 * sigma12; // determinant of sigmaXX matrix
160  if (std::isnan(det) || det < std::numeric_limits<float>::epsilon()) { // a suitably small number
161  return std::make_tuple(std::make_pair(false, det), NaN, NaN, NaN);
162  }
163  return std::make_tuple(std::make_pair(true, det), sigma22 / det, -sigma12 / det, sigma11 / det);
164 }
165 
167 bool shouldInterp(double sigma11, double sigma22, double det) {
168  float const xinterp = 0.25; // I.e. 0.5*0.5
169  return (sigma11 < xinterp || sigma22 < xinterp || det < xinterp * xinterp);
170 }
171 
172 // Decide on the bounding box for the region to examine while calculating the adaptive moments
173 // This routine will work in either LOCAL or PARENT coordinates (but of course which you pass
174 // determine which you will get back).
175 geom::Box2I computeAdaptiveMomentsBBox(geom::Box2I const &bbox, // full image bbox
176  geom::Point2D const &center, // centre of object
177  double sigma11_w, // quadratic moments of the
178  double, // weighting function
179  double sigma22_w, // xx, xy, and yy
180  double maxRadius = 1000 // Maximum radius of area to use
181  ) {
182  double radius = std::min(4 * std::sqrt(std::max(sigma11_w, sigma22_w)), maxRadius);
183  geom::Extent2D offset(radius, radius);
184  geom::Box2I result(geom::Box2D(center - offset, center + offset));
185  result.clip(bbox);
186  return result;
187 }
188 
189 /*****************************************************************************/
190 /*
191  * Calculate weighted moments of an object up to 2nd order
192  */
193 template <bool instFluxOnly, typename ImageT>
194 static int calcmom(ImageT const &image, // the image data
195  float xcen, float ycen, // centre of object
196  geom::BoxI bbox, // bounding box to consider
197  float bkgd, // data's background level
198  bool interpflag, // interpolate within pixels?
199  double w11, double w12, double w22, // weights
200  double *pI0, // amplitude of fit (if !NULL)
201  double *psum, // sum w*I (if !NULL)
202  double *psumx, double *psumy, // sum [xy]*w*I (if !instFluxOnly)
203  double *psumxx, double *psumxy, double *psumyy, // sum [xy]^2*w*I (if !instFluxOnly)
204  double *psums4, // sum w*I*weight^2 (if !instFluxOnly && !NULL)
205  bool negative = false) {
206  float tmod, ymod;
207  float X, Y; // sub-pixel interpolated [xy]
208  float weight;
209  float tmp;
210  double sum, sumx, sumy, sumxx, sumyy, sumxy, sums4;
211 #define RECALC_W 0 // estimate sigmaXX_w within BBox?
212 #if RECALC_W
213  double wsum, wsumxx, wsumxy, wsumyy;
214 
215  wsum = wsumxx = wsumxy = wsumyy = 0;
216 #endif
217 
218  assert(w11 >= 0); // i.e. it was set
219  if (fabs(w11) > 1e6 || fabs(w12) > 1e6 || fabs(w22) > 1e6) {
220  return (-1);
221  }
222 
223  sum = sumx = sumy = sumxx = sumxy = sumyy = sums4 = 0;
224 
225  int const ix0 = bbox.getMinX(); // corners of the box being analyzed
226  int const ix1 = bbox.getMaxX();
227  int const iy0 = bbox.getMinY(); // corners of the box being analyzed
228  int const iy1 = bbox.getMaxY();
229 
230  if (ix0 < 0 || ix1 >= image.getWidth() || iy0 < 0 || iy1 >= image.getHeight()) {
231  return -1;
232  }
233 
234  for (int i = iy0; i <= iy1; ++i) {
235  typename ImageT::x_iterator ptr = image.x_at(ix0, i);
236  float const y = i - ycen;
237  float const y2 = y * y;
238  float const yl = y - 0.375;
239  float const yh = y + 0.375;
240  for (int j = ix0; j <= ix1; ++j, ++ptr) {
241  float x = j - xcen;
242  if (interpflag) {
243  float const xl = x - 0.375;
244  float const xh = x + 0.375;
245 
246  float expon = xl * xl * w11 + yl * yl * w22 + 2.0 * xl * yl * w12;
247  tmp = xh * xh * w11 + yh * yh * w22 + 2.0 * xh * yh * w12;
248  expon = (expon > tmp) ? expon : tmp;
249  tmp = xl * xl * w11 + yh * yh * w22 + 2.0 * xl * yh * w12;
250  expon = (expon > tmp) ? expon : tmp;
251  tmp = xh * xh * w11 + yl * yl * w22 + 2.0 * xh * yl * w12;
252  expon = (expon > tmp) ? expon : tmp;
253 
254  if (expon <= 9.0) {
255  tmod = *ptr - bkgd;
256  for (Y = yl; Y <= yh; Y += 0.25) {
257  double const interpY2 = Y * Y;
258  for (X = xl; X <= xh; X += 0.25) {
259  double const interpX2 = X * X;
260  double const interpXy = X * Y;
261  expon = interpX2 * w11 + 2 * interpXy * w12 + interpY2 * w22;
262  weight = std::exp(-0.5 * expon);
263 
264  ymod = tmod * weight;
265  sum += ymod;
266  if (!instFluxOnly) {
267  sumx += ymod * (X + xcen);
268  sumy += ymod * (Y + ycen);
269 #if RECALC_W
270  wsum += weight;
271 
272  tmp = interpX2 * weight;
273  wsumxx += tmp;
274  sumxx += tmod * tmp;
275 
276  tmp = interpXy * weight;
277  wsumxy += tmp;
278  sumxy += tmod * tmp;
279 
280  tmp = interpY2 * weight;
281  wsumyy += tmp;
282  sumyy += tmod * tmp;
283 #else
284  sumxx += interpX2 * ymod;
285  sumxy += interpXy * ymod;
286  sumyy += interpY2 * ymod;
287 #endif
288  sums4 += expon * expon * ymod;
289  }
290  }
291  }
292  }
293  } else {
294  float x2 = x * x;
295  float xy = x * y;
296  float expon = x2 * w11 + 2 * xy * w12 + y2 * w22;
297 
298  if (expon <= 14.0) {
299  weight = std::exp(-0.5 * expon);
300  tmod = *ptr - bkgd;
301  ymod = tmod * weight;
302  sum += ymod;
303  if (!instFluxOnly) {
304  sumx += ymod * j;
305  sumy += ymod * i;
306 #if RECALC_W
307  wsum += weight;
308 
309  tmp = x2 * weight;
310  wsumxx += tmp;
311  sumxx += tmod * tmp;
312 
313  tmp = xy * weight;
314  wsumxy += tmp;
315  sumxy += tmod * tmp;
316 
317  tmp = y2 * weight;
318  wsumyy += tmp;
319  sumyy += tmod * tmp;
320 #else
321  sumxx += x2 * ymod;
322  sumxy += xy * ymod;
323  sumyy += y2 * ymod;
324 #endif
325  sums4 += expon * expon * ymod;
326  }
327  }
328  }
329  }
330  }
331 
332  if (pI0) {
333  std::tuple<std::pair<bool, double>, double, double, double> const weights = getWeights(w11, w12, w22);
334  double const detW = std::get<1>(weights) * std::get<3>(weights) - std::pow(std::get<2>(weights), 2);
335  *pI0 = sum / (geom::PI * sqrt(detW));
336  }
337  if (psum) {
338  *psum = sum;
339  }
340  if (!instFluxOnly) {
341  *psumx = sumx;
342  *psumy = sumy;
343  *psumxx = sumxx;
344  *psumxy = sumxy;
345  *psumyy = sumyy;
346  if (psums4 != NULL) {
347  *psums4 = sums4;
348  }
349  }
350 
351 #if RECALC_W
352  if (wsum > 0 && !instFluxOnly) {
353  double det = w11 * w22 - w12 * w12;
354  wsumxx /= wsum;
355  wsumxy /= wsum;
356  wsumyy /= wsum;
357  printf("%g %g %g %g %g %g\n", w22 / det, -w12 / det, w11 / det, wsumxx, wsumxy, wsumyy);
358  }
359 #endif
360 
361  if (negative) {
362  return (instFluxOnly || (sum < 0 && sumxx < 0 && sumyy < 0)) ? 0 : -1;
363  } else {
364  return (instFluxOnly || (sum > 0 && sumxx > 0 && sumyy > 0)) ? 0 : -1;
365  }
366 }
367 
368 /*
369  * Workhorse for adaptive moments
370  *
371  * All inputs are expected to be in LOCAL image coordinates
372  */
373 template <typename ImageT>
374 bool getAdaptiveMoments(ImageT const &mimage, double bkgd, double xcen, double ycen, double shiftmax,
375  SdssShapeResult *shape, int maxIter, float tol1, float tol2, bool negative) {
376  double I0 = 0; // amplitude of best-fit Gaussian
377  double sum; // sum of intensity*weight
378  double sumx, sumy; // sum ((int)[xy])*intensity*weight
379  double sumxx, sumxy, sumyy; // sum {x^2,xy,y^2}*intensity*weight
380  double sums4; // sum intensity*weight*exponent^2
381  float const xcen0 = xcen; // initial centre
382  float const ycen0 = ycen; // of object
383 
384  double sigma11W = 1.5; // quadratic moments of the
385  double sigma12W = 0.0; // weighting fcn;
386  double sigma22W = 1.5; // xx, xy, and yy
387 
388  double w11 = -1, w12 = -1, w22 = -1; // current weights for moments; always set when iter == 0
389  float e1_old = 1e6, e2_old = 1e6; // old values of shape parameters e1 and e2
390  float sigma11_ow_old = 1e6; // previous version of sigma11_ow
391 
392  typename ImageAdaptor<ImageT>::Image const &image = ImageAdaptor<ImageT>().getImage(mimage);
393 
394  if (std::isnan(xcen) || std::isnan(ycen)) {
395  // Can't do anything
396  shape->flags[SdssShapeAlgorithm::UNWEIGHTED_BAD.number] = true;
397  return false;
398  }
399 
400  bool interpflag = false; // interpolate finer than a pixel?
402  int iter = 0; // iteration number
403  for (; iter < maxIter; iter++) {
404  bbox = computeAdaptiveMomentsBBox(image.getBBox(afw::image::LOCAL), geom::Point2D(xcen, ycen),
405  sigma11W, sigma12W, sigma22W);
406  std::tuple<std::pair<bool, double>, double, double, double> weights =
407  getWeights(sigma11W, sigma12W, sigma22W);
408  if (!std::get<0>(weights).first) {
409  shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
410  break;
411  }
412 
413  double const detW = std::get<0>(weights).second;
414 
415 #if 0 // this form was numerically unstable on my G4 powerbook
416  assert(detW >= 0.0);
417 #else
418  assert(sigma11W * sigma22W >= sigma12W * sigma12W - std::numeric_limits<float>::epsilon());
419 #endif
420 
421  {
422  const double ow11 = w11; // old
423  const double ow12 = w12; // values
424  const double ow22 = w22; // of w11, w12, w22
425 
426  w11 = std::get<1>(weights);
427  w12 = std::get<2>(weights);
428  w22 = std::get<3>(weights);
429 
430  if (shouldInterp(sigma11W, sigma22W, detW)) {
431  if (!interpflag) {
432  interpflag = true; // N.b.: stays set for this object
433  if (iter > 0) {
434  sigma11_ow_old = 1.e6; // force at least one more iteration
435  w11 = ow11;
436  w12 = ow12;
437  w22 = ow22;
438  iter--; // we didn't update wXX
439  }
440  }
441  }
442  }
443 
444  if (calcmom<false>(image, xcen, ycen, bbox, bkgd, interpflag, w11, w12, w22, &I0, &sum, &sumx, &sumy,
445  &sumxx, &sumxy, &sumyy, &sums4, negative) < 0) {
446  shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
447  break;
448  }
449 
450 #if 0
451 /*
452  * Find new centre
453  *
454  * This is only needed if we update the centre; if we use the input position we've already done the work
455  */
456  xcen = sumx/sum;
457  ycen = sumy/sum;
458 #endif
459  shape->x = sumx / sum; // update centroid. N.b. we're not setting errors here
460  shape->y = sumy / sum;
461 
462  if (fabs(shape->x - xcen0) > shiftmax || fabs(shape->y - ycen0) > shiftmax) {
463  shape->flags[SdssShapeAlgorithm::SHIFT.number] = true;
464  }
465  /*
466  * OK, we have the centre. Proceed to find the second moments.
467  */
468  float const sigma11_ow = sumxx / sum; // quadratic moments of
469  float const sigma22_ow = sumyy / sum; // weight*object
470  float const sigma12_ow = sumxy / sum; // xx, xy, and yy
471 
472  if (sigma11_ow <= 0 || sigma22_ow <= 0) {
473  shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
474  break;
475  }
476 
477  float const d = sigma11_ow + sigma22_ow; // current values of shape parameters
478  float const e1 = (sigma11_ow - sigma22_ow) / d;
479  float const e2 = 2.0 * sigma12_ow / d;
480  /*
481  * Did we converge?
482  */
483  if (iter > 0 && fabs(e1 - e1_old) < tol1 && fabs(e2 - e2_old) < tol1 &&
484  fabs(sigma11_ow / sigma11_ow_old - 1.0) < tol2) {
485  break; // yes; we converged
486  }
487 
488  e1_old = e1;
489  e2_old = e2;
490  sigma11_ow_old = sigma11_ow;
491  /*
492  * Didn't converge, calculate new values for weighting function
493  *
494  * The product of two Gaussians is a Gaussian:
495  * <x^2 exp(-a x^2 - 2bxy - cy^2) exp(-Ax^2 - 2Bxy - Cy^2)> =
496  * <x^2 exp(-(a + A) x^2 - 2(b + B)xy - (c + C)y^2)>
497  * i.e. the inverses of the covariances matrices add.
498  *
499  * We know sigmaXX_ow and sigmaXXW, the covariances of the weighted object
500  * and of the weights themselves. We can estimate the object's covariance as
501  * sigmaXX_ow^-1 - sigmaXXW^-1
502  * and, as we want to find a set of weights with the _same_ covariance as the
503  * object we take this to be the an estimate of our correct weights.
504  *
505  * N.b. This assumes that the object is roughly Gaussian.
506  * Consider the object:
507  * O == delta(x + p) + delta(x - p)
508  * the covariance of the weighted object is equal to that of the unweighted
509  * object, and this prescription fails badly. If we detect this, we set
510  * the UNWEIGHTED flag, and calculate the UNweighted moments
511  * instead.
512  */
513  {
514  float n11, n12, n22; // elements of inverse of next guess at weighting function
515  float ow11, ow12, ow22; // elements of inverse of sigmaXX_ow
516 
517  std::tuple<std::pair<bool, double>, double, double, double> weights =
518  getWeights(sigma11_ow, sigma12_ow, sigma22_ow);
519  if (!std::get<0>(weights).first) {
520  shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
521  break;
522  }
523 
524  ow11 = std::get<1>(weights);
525  ow12 = std::get<2>(weights);
526  ow22 = std::get<3>(weights);
527 
528  n11 = ow11 - w11;
529  n12 = ow12 - w12;
530  n22 = ow22 - w22;
531 
532  weights = getWeights(n11, n12, n22);
533  if (!std::get<0>(weights).first) {
534  // product-of-Gaussians assumption failed
535  shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
536  break;
537  }
538 
539  sigma11W = std::get<1>(weights);
540  sigma12W = std::get<2>(weights);
541  sigma22W = std::get<3>(weights);
542  }
543 
544  if (sigma11W <= 0 || sigma22W <= 0) {
545  shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
546  break;
547  }
548  }
549 
550  if (iter == maxIter) {
551  shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
552  shape->flags[SdssShapeAlgorithm::MAXITER.number] = true;
553  }
554 
555  if (sumxx + sumyy == 0.0) {
556  shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
557  }
558  /*
559  * Problems; try calculating the un-weighted moments
560  */
561  if (shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number]) {
562  w11 = w22 = w12 = 0;
563  if (calcmom<false>(image, xcen, ycen, bbox, bkgd, interpflag, w11, w12, w22, &I0, &sum, &sumx, &sumy,
564  &sumxx, &sumxy, &sumyy, NULL, negative) < 0 ||
565  (!negative && sum <= 0) || (negative && sum >= 0)) {
566  shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = false;
567  shape->flags[SdssShapeAlgorithm::UNWEIGHTED_BAD.number] = true;
568 
569  if (sum > 0) {
570  shape->xx = 1 / 12.0; // a single pixel
571  shape->xy = 0.0;
572  shape->yy = 1 / 12.0;
573  }
574 
575  return false;
576  }
577 
578  sigma11W = sumxx / sum; // estimate of object moments
579  sigma12W = sumxy / sum; // usually, object == weight
580  sigma22W = sumyy / sum; // at this point
581  }
582 
583  shape->instFlux = I0;
584  shape->xx = sigma11W;
585  shape->xy = sigma12W;
586  shape->yy = sigma22W;
587 
588  if (shape->xx + shape->yy != 0.0) {
589  int const ix = afw::image::positionToIndex(xcen);
590  int const iy = afw::image::positionToIndex(ycen);
591 
592  if (ix >= 0 && ix < mimage.getWidth() && iy >= 0 && iy < mimage.getHeight()) {
593  float const bkgd_var = ImageAdaptor<ImageT>().getVariance(
594  mimage, ix, iy); // XXX Overestimate as it includes object
595 
596  if (bkgd_var > 0.0) { // NaN is not > 0.0
597  if (!(shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number])) {
598  Matrix4d fisher = calc_fisher(*shape, bkgd_var); // Fisher matrix
599  Matrix4d cov = fisher.inverse();
600  // convention is to order moments (xx, yy, xy)
601  shape->instFluxErr = std::sqrt(cov(0, 0));
602  shape->xxErr = std::sqrt(cov(1, 1));
603  shape->yyErr = std::sqrt(cov(2, 2));
604  shape->xyErr = std::sqrt(cov(3, 3));
605  shape->instFlux_xx_Cov = cov(0, 1);
606  shape->instFlux_yy_Cov = cov(0, 2);
607  shape->instFlux_xy_Cov = cov(0, 3);
608  shape->xx_yy_Cov = cov(1, 2);
609  shape->xx_xy_Cov = cov(1, 3);
610  shape->yy_xy_Cov = cov(2, 3);
611  }
612  }
613  }
614  }
615 
616  return true;
617 }
618 
619 } // namespace
620 
622  : instFlux_xx_Cov(std::numeric_limits<ErrElement>::quiet_NaN()),
623  instFlux_yy_Cov(std::numeric_limits<ErrElement>::quiet_NaN()),
624  instFlux_xy_Cov(std::numeric_limits<ErrElement>::quiet_NaN()) {}
625 
627  bool doMeasurePsf) {
629  r._shapeResult =
630  ShapeResultKey::addFields(schema, name, "elliptical Gaussian adaptive moments", SIGMA_ONLY);
631  r._centroidResult = CentroidResultKey::addFields(schema, name, "elliptical Gaussian adaptive moments",
633  r._instFluxResult = FluxResultKey::addFields(schema, name, "elliptical Gaussian adaptive moments");
634 
635  // Only include PSF shape fields if doMeasurePsf = True
636  if (doMeasurePsf) {
637  r._includePsf = true;
638  r._psfShapeResult = afw::table::QuadrupoleKey::addFields(
639  schema, schema.join(name, "psf"), "adaptive moments of the PSF model at the object position");
640  } else {
641  r._includePsf = false;
642  }
643 
644  r._instFlux_xx_Cov =
645  schema.addField<ErrElement>(schema.join(name, "instFlux", "xx", "Cov"),
646  (boost::format("uncertainty covariance between %s and %s") %
647  schema.join(name, "instFlux") % schema.join(name, "xx"))
648  .str(),
649  "count*pixel^2");
650  r._instFlux_yy_Cov =
651  schema.addField<ErrElement>(schema.join(name, "instFlux", "yy", "Cov"),
652  (boost::format("uncertainty covariance between %s and %s") %
653  schema.join(name, "instFlux") % schema.join(name, "yy"))
654  .str(),
655  "count*pixel^2");
656  r._instFlux_xy_Cov =
657  schema.addField<ErrElement>(schema.join(name, "instFlux", "xy", "Cov"),
658  (boost::format("uncertainty covariance between %s and %s") %
659  schema.join(name, "instFlux") % schema.join(name, "xy"))
660  .str(),
661  "count*pixel^2");
662 
663  // Skip the psf flag if not recording the PSF shape.
664  if (r._includePsf) {
666  } else {
669  }
670  return r;
671 }
672 
674  : _shapeResult(s),
675  _centroidResult(s),
676  _instFluxResult(s),
677  _instFlux_xx_Cov(s["instFlux"]["xx"]["Cov"]),
678  _instFlux_yy_Cov(s["instFlux"]["yy"]["Cov"]),
679  _instFlux_xy_Cov(s["instFlux"]["xy"]["Cov"]) {
680  // The input SubSchema may optionally provide for a PSF.
681  try {
682  _psfShapeResult = afw::table::QuadrupoleKey(s["psf"]);
684  _includePsf = true;
685  } catch (pex::exceptions::NotFoundError &e) {
686  _flagHandler =
688  _includePsf = false;
689  }
690 }
691 
694  static_cast<ShapeResult &>(result) = record.get(_shapeResult);
695  static_cast<CentroidResult &>(result) = record.get(_centroidResult);
696  static_cast<FluxResult &>(result) = record.get(_instFluxResult);
697  result.instFlux_xx_Cov = record.get(_instFlux_xx_Cov);
698  result.instFlux_yy_Cov = record.get(_instFlux_yy_Cov);
699  result.instFlux_xy_Cov = record.get(_instFlux_xy_Cov);
700  for (size_t n = 0; n < SdssShapeAlgorithm::N_FLAGS; ++n) {
701  if (n == SdssShapeAlgorithm::PSF_SHAPE_BAD.number && !_includePsf) continue;
702  result.flags[n] = _flagHandler.getValue(record, n);
703  }
704  return result;
705 }
706 
708  return record.get(_psfShapeResult);
709 }
710 
712  record.set(_shapeResult, value);
713  record.set(_centroidResult, value);
714  record.set(_instFluxResult, value);
715  record.set(_instFlux_xx_Cov, value.instFlux_xx_Cov);
716  record.set(_instFlux_yy_Cov, value.instFlux_yy_Cov);
717  record.set(_instFlux_xy_Cov, value.instFlux_xy_Cov);
718  for (size_t n = 0; n < SdssShapeAlgorithm::N_FLAGS; ++n) {
719  if (n == SdssShapeAlgorithm::PSF_SHAPE_BAD.number && !_includePsf) continue;
720  _flagHandler.setValue(record, n, value.flags[n]);
721  }
722 }
723 
725  afw::geom::ellipses::Quadrupole const &value) const {
726  record.set(_psfShapeResult, value);
727 }
728 
730  return _shapeResult == other._shapeResult && _centroidResult == other._centroidResult &&
731  _instFluxResult == other._instFluxResult && _psfShapeResult == other._psfShapeResult &&
732  _instFlux_xx_Cov == other._instFlux_xx_Cov && _instFlux_yy_Cov == other._instFlux_yy_Cov &&
733  _instFlux_xy_Cov == other._instFlux_xy_Cov;
734  // don't bother with flags - if we've gotten this far, it's basically impossible the flags don't match
735 }
736 
738  return _shapeResult.isValid() && _centroidResult.isValid() && _instFluxResult.isValid() &&
739  _psfShapeResult.isValid() && _instFlux_xx_Cov.isValid() && _instFlux_yy_Cov.isValid() &&
740  _instFlux_xy_Cov.isValid();
741  // don't bother with flags - if we've gotten this far, it's basically impossible the flags are invalid
742 }
743 
746  : _ctrl(ctrl),
747  _resultKey(ResultKey::addFields(schema, name, ctrl.doMeasurePsf)),
748  _centroidExtractor(schema, name) {}
749 
750 template <typename ImageT>
752  bool negative, Control const &control) {
753  double xcen = center.getX(); // object's column position
754  double ycen = center.getY(); // object's row position
755 
756  xcen -= image.getX0(); // work in image Pixel coordinates
757  ycen -= image.getY0();
758 
759  float shiftmax = control.maxShift; // Max allowed centroid shift
760  if (shiftmax < 2) {
761  shiftmax = 2;
762  } else if (shiftmax > 10) {
763  shiftmax = 10;
764  }
765 
767  try {
768  result.flags[FAILURE.number] =
769  !getAdaptiveMoments(image, control.background, xcen, ycen, shiftmax, &result, control.maxIter,
770  control.tol1, control.tol2, negative);
771  } catch (pex::exceptions::Exception &err) {
772  result.flags[FAILURE.number] = true;
773  }
774  if (result.flags[UNWEIGHTED.number] || result.flags[SHIFT.number]) {
775  // These are also considered fatal errors in terms of the quality of the results,
776  // even though they do produce some results.
777  result.flags[FAILURE.number] = true;
778  }
779  double IxxIyy = result.getQuadrupole().getIxx() * result.getQuadrupole().getIyy();
780  double Ixy_sq = result.getQuadrupole().getIxy();
781  Ixy_sq *= Ixy_sq;
782  double epsilon = 1.0e-6;
783  if (IxxIyy < (1.0 + epsilon) * Ixy_sq)
784  // We are checking that Ixx*Iyy > (1 + epsilon)*Ixy*Ixy where epsilon is suitably small. The
785  // value of epsilon used here is a magic number subject to future review (DM-5801 was marked won't fix).
786  {
787  if (!result.flags[FAILURE.number]) {
789  (boost::format("computeAdaptiveMoments IxxIxy %d < (1 + eps=%d)*(Ixy^2=%d);"
790  " implying singular moments without any flag set")
791  % IxxIyy % epsilon % Ixy_sq).str());
792  }
793  }
794 
795  // getAdaptiveMoments() just computes the zeroth moment in result.instFlux (and its error in
796  // result.instFluxErr, result.instFlux_xx_Cov, etc.) That's related to the instFlux by some geometric
797  // factors, which we apply here.
798  // This scale factor is just the inverse of the normalization constant in a bivariate normal distribution:
799  // 2*pi*sigma_x*sigma_y*sqrt(1-rho^2), where rho is the correlation.
800  // This happens to be twice the ellipse area pi*sigma_maj*sigma_min, which is identically equal to:
801  // pi*sqrt(I_xx*I_yy - I_xy^2); i.e. pi*sqrt(determinant(I))
802  double instFluxScale = geom::TWOPI * std::sqrt(IxxIyy - Ixy_sq);
803 
804  result.instFlux *= instFluxScale;
805  result.instFluxErr *= instFluxScale;
806  result.x += image.getX0();
807  result.y += image.getY0();
808 
809  if (ImageAdaptor<ImageT>::hasVariance) {
810  result.instFlux_xx_Cov *= instFluxScale;
811  result.instFlux_yy_Cov *= instFluxScale;
812  result.instFlux_xy_Cov *= instFluxScale;
813  }
814 
815  return result;
816 }
817 
818 template <typename ImageT>
820  afw::geom::ellipses::Quadrupole const &shape,
821  geom::Point2D const &center) {
822  // while arguments to computeFixedMomentsFlux are in PARENT coordinates, the implementation is LOCAL.
823  geom::Point2D localCenter = center - geom::Extent2D(image.getXY0());
824 
825  geom::BoxI const bbox = computeAdaptiveMomentsBBox(image.getBBox(afw::image::LOCAL), localCenter,
826  shape.getIxx(), shape.getIxy(), shape.getIyy());
827 
828  std::tuple<std::pair<bool, double>, double, double, double> weights =
829  getWeights(shape.getIxx(), shape.getIxy(), shape.getIyy());
830 
832 
833  if (!std::get<0>(weights).first) {
834  throw pex::exceptions::InvalidParameterError("Input shape is singular");
835  }
836 
837  double const w11 = std::get<1>(weights);
838  double const w12 = std::get<2>(weights);
839  double const w22 = std::get<3>(weights);
840  bool const interp = shouldInterp(shape.getIxx(), shape.getIyy(), std::get<0>(weights).second);
841 
842  double sum0 = 0; // sum of pixel values weighted by a Gaussian
843  if (calcmom<true>(ImageAdaptor<ImageT>().getImage(image), localCenter.getX(), localCenter.getY(), bbox,
844  0.0, interp, w11, w12, w22, NULL, &sum0, NULL, NULL, NULL, NULL, NULL, NULL) < 0) {
845  throw LSST_EXCEPT(pex::exceptions::RuntimeError, "Error from calcmom");
846  }
847 
848  result.instFlux = sum0 * 2.0;
849 
850  if (ImageAdaptor<ImageT>::hasVariance) {
851  int ix = static_cast<int>(center.getX() - image.getX0());
852  int iy = static_cast<int>(center.getY() - image.getY0());
853  if (!image.getBBox(afw::image::LOCAL).contains(geom::Point2I(ix, iy))) {
855  (boost::format("Center (%d,%d) not in image (%dx%d)") % ix % iy %
856  image.getWidth() % image.getHeight())
857  .str());
858  }
859  double var = ImageAdaptor<ImageT>().getVariance(image, ix, iy);
860  // 0th moment (i0) error = sqrt(var / wArea); instFlux (error) = 2 * wArea * i0 (error)
861  double const wArea = geom::PI * std::sqrt(shape.getDeterminant());
862  result.instFluxErr = 2 * std::sqrt(var * wArea);
863  }
864 
865  return result;
866 }
867 
869  afw::image::Exposure<float> const &exposure) const {
870  bool negative = false;
871 
872  try {
873  negative = measRecord.get(measRecord.getSchema().find<afw::table::Flag>("flags_negative").key);
874  } catch (pexExcept::Exception &e) {
875  }
877  exposure.getMaskedImage(), _centroidExtractor(measRecord, _resultKey.getFlagHandler()), negative,
878  _ctrl);
879 
880  if (_ctrl.doMeasurePsf) {
881  // Compute moments of Psf model. In the interest of implementing this quickly, we're just
882  // calling Psf::computeShape(), which delegates to SdssShapeResult::computeAdaptiveMoments
883  // for all nontrivial Psf classes. But this could in theory save the results of a shape
884  // computed some other way as part of base_SdssShape, which might be confusing. We should
885  // fix this eventually either by making Psf shape measurement not part of base_SdssShape, or
886  // by making the measurements stored with shape.sdss always computed via the
887  // SdssShapeAlgorithm instead of delegating to the Psf class.
888  try {
890  if (!psf) {
891  result.flags[PSF_SHAPE_BAD.number] = true;
892  } else {
893  _resultKey.setPsfShape(measRecord, psf->computeShape(geom::Point2D(result.x, result.y)));
894  }
895  } catch (pex::exceptions::Exception &err) {
896  result.flags[PSF_SHAPE_BAD.number] = true;
897  }
898  }
899 
900  measRecord.set(_resultKey, result);
901 }
902 
904  _resultKey.getFlagHandler().handleFailure(measRecord, error);
905 }
906 
907 #define INSTANTIATE_IMAGE(IMAGE) \
908  template SdssShapeResult SdssShapeAlgorithm::computeAdaptiveMoments( \
909  IMAGE const &, geom::Point2D const &, bool, Control const &); \
910  template FluxResult SdssShapeAlgorithm::computeFixedMomentsFlux( \
911  IMAGE const &, afw::geom::ellipses::Quadrupole const &, geom::Point2D const &)
912 
913 #define INSTANTIATE_PIXEL(PIXEL) \
914  INSTANTIATE_IMAGE(afw::image::Image<PIXEL>); \
915  INSTANTIATE_IMAGE(afw::image::MaskedImage<PIXEL>);
916 
917 INSTANTIATE_PIXEL(int);
918 INSTANTIATE_PIXEL(float);
919 INSTANTIATE_PIXEL(double);
920 
923  : BaseTransform{name}, _instFluxTransform{name, mapper}, _centroidTransform{name, mapper} {
924  // If the input schema has a PSF flag -- it's optional -- assume we are also transforming the PSF.
925  _transformPsf = mapper.getInputSchema().getNames().count("sdssShape_flag_psf") ? true : false;
926 
927  // Skip the last flag if not transforming the PSF shape.
928  for (std::size_t i = 0; i < SdssShapeAlgorithm::getFlagDefinitions().size(); i++) {
930  if (flag == SdssShapeAlgorithm::FAILURE) continue;
931  if (mapper.getInputSchema().getNames().count(mapper.getInputSchema().join(name, flag.name)) == 0)
932  continue;
934  mapper.getInputSchema().find<afw::table::Flag>(name + "_" + flag.name).key;
935  mapper.addMapping(key);
936  }
937 
938  _outShapeKey = ShapeResultKey::addFields(mapper.editOutputSchema(), name, "Shape in celestial moments",
940  if (_transformPsf) {
941  _outPsfShapeKey = afw::table::QuadrupoleKey::addFields(mapper.editOutputSchema(), name + "_psf",
942  "PSF shape in celestial moments",
944  }
945 }
946 
948  afw::table::BaseCatalog &outputCatalog, afw::geom::SkyWcs const &wcs,
949  afw::image::PhotoCalib const &photoCalib) const {
950  // The instFlux and cetroid transforms will check that the catalog lengths
951  // match and throw if not, so we don't repeat the test here.
952  _instFluxTransform(inputCatalog, outputCatalog, wcs, photoCalib);
953  _centroidTransform(inputCatalog, outputCatalog, wcs, photoCalib);
954 
955  CentroidResultKey centroidKey(inputCatalog.getSchema()[_name]);
956  ShapeResultKey inShapeKey(inputCatalog.getSchema()[_name]);
957  afw::table::QuadrupoleKey inPsfShapeKey;
958  if (_transformPsf) {
959  inPsfShapeKey = afw::table::QuadrupoleKey(
960  inputCatalog.getSchema()[inputCatalog.getSchema().join(_name, "psf")]);
961  }
962 
963  afw::table::SourceCatalog::const_iterator inSrc = inputCatalog.begin();
964  afw::table::BaseCatalog::iterator outSrc = outputCatalog.begin();
965  for (; inSrc != inputCatalog.end(); ++inSrc, ++outSrc) {
966  ShapeResult inShape = inShapeKey.get(*inSrc);
967  ShapeResult outShape;
968 
969  // The transformation from the (x, y) to the (Ra, Dec) basis.
970  geom::AffineTransform crdTr =
971  wcs.linearizePixelToSky(centroidKey.get(*inSrc).getCentroid(), geom::radians);
972  outShape.setShape(inShape.getShape().transform(crdTr.getLinear()));
973 
974  // Transformation matrix from pixel to celestial basis.
976  outShape.setShapeErr((m * inShape.getShapeErr().cast<double>() * m.transpose()).cast<ErrElement>());
977 
978  _outShapeKey.set(*outSrc, outShape);
979 
980  if (_transformPsf) {
981  _outPsfShapeKey.set(*outSrc, inPsfShapeKey.get(*inSrc).transform(crdTr.getLinear()));
982  }
983  }
984 }
985 
986 } // namespace base
987 } // namespace meas
988 } // namespace lsst
py::object result
Definition: _schema.cc:429
table::Key< std::string > name
Definition: Amplifier.cc:116
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
double x
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
afw::table::Key< afw::table::Array< ImagePixelT > > image
uint64_t * ptr
Definition: RangeSet.cc:88
SchemaMapper * mapper
Definition: SchemaMapper.cc:71
table::Key< table::Array< std::uint8_t > > wcs
Definition: SkyWcs.cc:66
int y
Definition: SpanSet.cc:48
int m
Definition: SpanSet.cc:48
table::Schema schema
Definition: python.h:134
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Definition: SkyWcs.h:117
Transformer transform(lsst::geom::LinearTransform const &transform)
Return the transform that maps the ellipse to the unit circle.
Definition: Transformer.h:116
An ellipse core with quadrupole moments as parameters.
Definition: Quadrupole.h:47
double getDeterminant() const
Return the determinant of the matrix representation.
Definition: Quadrupole.h:83
MaskedImageT getMaskedImage()
Return the MaskedImage.
Definition: Exposure.h:228
std::shared_ptr< lsst::afw::detection::Psf const > getPsf() const
Return the Exposure's Psf object.
Definition: Exposure.h:327
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:73
lsst::afw::image::Image< ImagePixelT > Image
Definition: MaskedImage.h:85
The photometric calibration of an exposure.
Definition: PhotoCalib.h:114
Base class for all records.
Definition: BaseRecord.h:31
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
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Definition: BaseRecord.h:151
Iterator class for CatalogT.
Definition: Catalog.h:40
iterator begin()
Iterator access.
Definition: Catalog.h:401
bool isValid() const noexcept
Return true if the key was initialized to valid offset.
Definition: Key.h:97
A FunctorKey used to get or set a geom::ellipses::Quadrupole from a tuple of constituent Keys.
Definition: aggregates.h:282
geom::ellipses::Quadrupole get(BaseRecord const &record) const override
Get a Quadrupole from the given record.
Definition: aggregates.cc:110
void set(BaseRecord &record, geom::ellipses::Quadrupole const &value) const override
Set a Quadrupole in the given record.
Definition: aggregates.cc:114
static QuadrupoleKey addFields(Schema &schema, std::string const &name, std::string const &doc, CoordinateType coordType=CoordinateType::PIXEL)
Add a set of quadrupole subfields to a schema and return a QuadrupoleKey that points to them.
Definition: aggregates.cc:100
bool isValid() const noexcept
Return True if all the constituent Keys are valid.
Definition: aggregates.h:342
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.
Definition: SchemaMapper.h:21
typename Base::const_iterator const_iterator
Definition: SortedCatalog.h:50
Record class that contains measurements made on a single exposure.
Definition: Source.h:78
A proxy type for name lookups in a Schema.
Definition: Schema.h:367
An affine coordinate transformation consisting of a linear transformation and an offset.
LinearTransform const & getLinear() const noexcept
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
An integer coordinate rectangle.
Definition: Box.h:55
Abstract base class for all C++ measurement transformations.
Definition: Transform.h:86
A FunctorKey for CentroidResult.
virtual CentroidResult get(afw::table::BaseRecord const &record) const
Get a CentroidResult from the given record.
bool isValid() const
Return True if the centroid key is valid.
static CentroidResultKey addFields(afw::table::Schema &schema, std::string const &name, std::string const &doc, UncertaintyEnum uncertainty)
Add the appropriate fields to a Schema, and return a CentroidResultKey that manages them.
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
Definition: FlagHandler.h:125
Utility class for handling flag fields that indicate the failure modes of an algorithm.
Definition: FlagHandler.h:148
void handleFailure(afw::table::BaseRecord &record, MeasurementError const *error=nullptr) const
Handle an expected or unexpected Exception thrown by a measurement algorithm.
Definition: FlagHandler.cc:76
static FlagHandler addFields(afw::table::Schema &schema, std::string const &prefix, FlagDefinitionList const &flagDefs, FlagDefinitionList const &exclDefs=FlagDefinitionList::getEmptyList())
Add Flag fields to a schema, creating a FlagHandler object to manage them.
Definition: FlagHandler.cc:37
void setValue(afw::table::BaseRecord &record, std::size_t i, bool value) const
Set the flag field corresponding to the given flag index.
Definition: FlagHandler.h:262
bool getValue(afw::table::BaseRecord const &record, std::size_t i) const
Return the value of the flag field corresponding to the given flag index.
Definition: FlagHandler.h:242
bool isValid() const
Return True if both the instFlux and instFluxErr Keys are valid.
static FluxResultKey addFields(afw::table::Schema &schema, std::string const &name, std::string const &doc)
Add a pair of _instFlux, _instFluxErr fields to a Schema, and return a FluxResultKey that points to t...
Exception to be thrown when a measurement algorithm experiences a known failure mode.
Definition: exceptions.h:48
static FlagDefinition const SHIFT
Definition: SdssShape.h:158
static FlagDefinitionList const & getFlagDefinitions()
Definition: SdssShape.cc:58
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: SdssShape.cc:868
static Result computeAdaptiveMoments(ImageT const &image, geom::Point2D const &position, bool negative=false, Control const &ctrl=Control())
Compute the adaptive Gaussian-weighted moments of an image.
static FlagDefinition const FAILURE
Definition: SdssShape.h:155
static unsigned int const N_FLAGS
Definition: SdssShape.h:154
static FluxResult computeFixedMomentsFlux(ImageT const &image, afw::geom::ellipses::Quadrupole const &shape, geom::Point2D const &position)
Compute the instFlux within a fixed Gaussian aperture.
Definition: SdssShape.cc:819
static FlagDefinition const MAXITER
Definition: SdssShape.h:159
static FlagDefinition const UNWEIGHTED
Definition: SdssShape.h:157
static FlagDefinition const PSF_SHAPE_BAD
Definition: SdssShape.h:160
static FlagDefinition const UNWEIGHTED_BAD
Definition: SdssShape.h:156
SdssShapeAlgorithm(Control const &ctrl, std::string const &name, afw::table::Schema &schema)
Definition: SdssShape.cc:744
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: SdssShape.cc:903
A C++ control class to handle SdssShapeAlgorithm's configuration.
Definition: SdssShape.h:52
double maxShift
"Maximum centroid shift, limited to 2-10" ;
Definition: SdssShape.h:56
int maxIter
"Maximum number of iterations" ;
Definition: SdssShape.h:55
float tol2
"Convergence tolerance for FWHM" ;
Definition: SdssShape.h:58
bool doMeasurePsf
"Whether to also compute the shape of the PSF model" ;
Definition: SdssShape.h:59
float tol1
"Convergence tolerance for e1,e2" ;
Definition: SdssShape.h:57
double background
"Additional value to add to background" ;
Definition: SdssShape.h:54
Result object SdssShapeAlgorithm.
Definition: SdssShape.h:224
ErrElement instFlux_yy_Cov
instFlux, yy term in the uncertainty covariance matrix
Definition: SdssShape.h:227
SdssShapeResult()
Constructor; initializes everything to NaN.
Definition: SdssShape.cc:621
ErrElement instFlux_xy_Cov
instFlux, xy term in the uncertainty covariance matrix
Definition: SdssShape.h:228
std::bitset< SdssShapeAlgorithm::N_FLAGS > flags
Status flags (see SdssShapeAlgorithm).
Definition: SdssShape.h:230
ErrElement instFlux_xx_Cov
instFlux, xx term in the uncertainty covariance matrix
Definition: SdssShape.h:226
A FunctorKey that maps SdssShapeResult to afw::table Records.
Definition: SdssShape.h:73
FlagHandler const & getFlagHandler() const
Definition: SdssShape.h:125
SdssShapeResultKey()
Default constructor; instance will not be usuable unless subsequently assigned to.
Definition: SdssShape.h:90
virtual SdssShapeResult get(afw::table::BaseRecord const &record) const
Get an SdssShapeResult from the given record.
Definition: SdssShape.cc:692
static SdssShapeResultKey addFields(afw::table::Schema &schema, std::string const &name, bool doMeasurePsf)
Add the appropriate fields to a Schema, and return a SdssShapeResultKey that manages them.
Definition: SdssShape.cc:626
virtual afw::geom::ellipses::Quadrupole getPsfShape(afw::table::BaseRecord const &record) const
Get a Quadrupole for the Psf from the given record.
Definition: SdssShape.cc:707
bool isValid() const
Return True if the key is valid.
Definition: SdssShape.cc:737
bool operator==(SdssShapeResultKey const &other) const
Compare the FunctorKey for equality with another, using the underlying Keys.
Definition: SdssShape.cc:729
virtual void set(afw::table::BaseRecord &record, SdssShapeResult const &value) const
Set an SdssShapeResult in the given record.
Definition: SdssShape.cc:711
virtual void setPsfShape(afw::table::BaseRecord &record, afw::geom::ellipses::Quadrupole const &value) const
Set a Quadrupole for the Psf at the position of the given record.
Definition: SdssShape.cc:724
virtual void operator()(afw::table::SourceCatalog const &inputCatalog, afw::table::BaseCatalog &outputCatalog, afw::geom::SkyWcs const &wcs, afw::image::PhotoCalib const &photoCalib) const
Definition: SdssShape.cc:947
SdssShapeTransform(Control const &ctrl, std::string const &name, afw::table::SchemaMapper &mapper)
Definition: SdssShape.cc:921
A FunctorKey for ShapeResult.
static ShapeResultKey addFields(afw::table::Schema &schema, std::string const &name, std::string const &doc, UncertaintyEnum uncertainty, afw::table::CoordinateType coordType=afw::table::CoordinateType::PIXEL)
Add the appropriate fields to a Schema, and return a ShapeResultKey that manages them.
virtual ShapeResult get(afw::table::BaseRecord const &record) const
Get a ShapeResult from the given record.
bool isValid() const
Return True if the shape key is valid.
virtual void set(afw::table::BaseRecord &record, ShapeResult const &value) const
Set a ShapeResult in the given record.
Reports arguments outside the domain of an operation.
Definition: Runtime.h:57
Provides consistent interface for LSST exceptions.
Definition: Exception.h:107
Reports invalid arguments.
Definition: Runtime.h:66
Reports errors in the logical structure of the program.
Definition: Runtime.h:46
Reports attempts to access elements using an invalid key.
Definition: Runtime.h:151
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104
T exp(T... args)
T fabs(T... args)
T printf(T... args)
T isnan(T... args)
T make_pair(T... args)
T make_tuple(T... args)
T max(T... args)
T min(T... args)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
int positionToIndex(double pos)
Convert image position to nearest integer index.
Definition: ImageUtils.h:69
constexpr double PI
The ratio of a circle's circumference to diameter.
Definition: Angle.h:39
Extent< double, 2 > Extent2D
Definition: Extent.h:400
constexpr double TWOPI
Definition: Angle.h:40
constexpr AngleUnit radians
constant with units of radians
Definition: Angle.h:108
@ SIGMA_ONLY
Only the diagonal elements of the covariance matrix are provided.
Definition: constants.h:45
@ FULL_COVARIANCE
The full covariance matrix is provided.
Definition: constants.h:46
@ NO_UNCERTAINTY
Algorithm provides no uncertainy information at all.
Definition: constants.h:44
Eigen::Matrix< ShapeElement, 3, 3, Eigen::DontAlign > ShapeTrMatrix
Definition: constants.h:62
float ErrElement
Definition: constants.h:55
ShapeTrMatrix makeShapeTransformMatrix(geom::LinearTransform const &xform)
Construct a matrix suitable for transforming second moments.
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
STL namespace.
T pow(T... args)
T quiet_NaN(T... args)
T sqrt(T... args)
afw::table::Key< double > weight
#define INSTANTIATE_PIXEL(PIXEL)
Definition: SdssShape.cc:913
A reusable struct for centroid measurements.
Centroid const getCentroid() const
Return a Point object containing the measured x and y.
Simple class used to define and document flags The name and doc constitute the identity of the FlagDe...
Definition: FlagHandler.h:40
A reusable result struct for instFlux measurements.
Definition: FluxUtilities.h:41
meas::base::Flux instFlux
Measured instFlux in DN.
Definition: FluxUtilities.h:42
A reusable struct for moments-based shape measurements.
Shape const getShape() const
Return an afw::geom::ellipses object corresponding to xx, yy, xy.
ShapeCov const getShapeErr() const
Return the 3x3 symmetric covariance matrix, with rows and columns ordered (xx, yy,...
void setShape(Shape const &shape)
Set struct elements from the given Quadrupole object.
ShapeElement xy
image or model second moment for xy^2
ShapeElement xx
image or model second moment for x^2
void setShapeErr(ShapeCov const &matrix)
Set the struct standard deviation elements from the given matrix, with rows and columns ordered (xx,...
ShapeElement yy
image or model second moment for y^2
Key< int > psf
Definition: Exposure.cc:65
Key< int > photoCalib
Definition: Exposure.cc:67