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