LSSTApplications  21.0.0+75b29a8a7f,21.0.0+e70536a077,21.0.0-1-ga51b5d4+62c747d40b,21.0.0-11-ga6ea59e8e+47cba9fc36,21.0.0-2-g103fe59+914993bf7c,21.0.0-2-g1367e85+e2614ded12,21.0.0-2-g45278ab+e70536a077,21.0.0-2-g4bc9b9f+7b2b5f8678,21.0.0-2-g5242d73+e2614ded12,21.0.0-2-g54e2caa+6403186824,21.0.0-2-g7f82c8f+3ac4acbffc,21.0.0-2-g8dde007+04a6aea1af,21.0.0-2-g8f08a60+9402881886,21.0.0-2-ga326454+3ac4acbffc,21.0.0-2-ga63a54e+81dd751046,21.0.0-2-gc738bc1+5f65c6e7a9,21.0.0-2-gde069b7+26c92b3210,21.0.0-2-gecfae73+0993ddc9bd,21.0.0-2-gfc62afb+e2614ded12,21.0.0-21-gba890a8+5a4f502a26,21.0.0-23-g9966ff26+03098d1af8,21.0.0-3-g357aad2+8ad216c477,21.0.0-3-g4be5c26+e2614ded12,21.0.0-3-g6d51c4a+4d2fe0280d,21.0.0-3-g7d9da8d+75b29a8a7f,21.0.0-3-gaa929c8+522e0f12c2,21.0.0-3-ge02ed75+4d2fe0280d,21.0.0-4-g3300ddd+e70536a077,21.0.0-4-gc004bbf+eac6615e82,21.0.0-4-gccdca77+f94adcd104,21.0.0-4-gd1c1571+18b81799f9,21.0.0-5-g7b47fff+4d2fe0280d,21.0.0-5-gb155db7+d2632f662b,21.0.0-5-gdf36809+637e4641ee,21.0.0-6-g722ad07+28c848f42a,21.0.0-7-g959bb79+522e0f12c2,21.0.0-7-gfd72ab2+cf01990774,21.0.0-9-g87fb7b8d+e2ab11cdd6,w.2021.04
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 /*****************************************************************************/
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) {
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 in afw::geom::ellipses is to order moments (xx, yy, xy),
601  // but the older algorithmic code uses (xx, xy, yy) - the order of
602  // indices here is not a bug.
603  shape->instFluxErr = std::sqrt(cov(0, 0));
604  shape->xxErr = std::sqrt(cov(1, 1));
605  shape->xyErr = std::sqrt(cov(2, 2));
606  shape->yyErr = std::sqrt(cov(3, 3));
607  shape->instFlux_xx_Cov = cov(0, 1);
608  shape->instFlux_xy_Cov = cov(0, 2);
609  shape->instFlux_yy_Cov = cov(0, 3);
610  shape->xx_yy_Cov = cov(1, 3);
611  shape->xx_xy_Cov = cov(1, 2);
612  shape->yy_xy_Cov = cov(2, 3);
613  }
614  }
615  }
616  }
617 
618  return true;
619 }
620 
621 } // namespace
622 
624  : instFlux_xx_Cov(std::numeric_limits<ErrElement>::quiet_NaN()),
625  instFlux_yy_Cov(std::numeric_limits<ErrElement>::quiet_NaN()),
626  instFlux_xy_Cov(std::numeric_limits<ErrElement>::quiet_NaN()) {}
627 
629  bool doMeasurePsf) {
631  r._shapeResult =
632  ShapeResultKey::addFields(schema, name, "elliptical Gaussian adaptive moments", SIGMA_ONLY);
633  r._centroidResult = CentroidResultKey::addFields(schema, name, "elliptical Gaussian adaptive moments",
635  r._instFluxResult = FluxResultKey::addFields(schema, name, "elliptical Gaussian adaptive moments");
636 
637  // Only include PSF shape fields if doMeasurePsf = True
638  if (doMeasurePsf) {
639  r._includePsf = true;
640  r._psfShapeResult = afw::table::QuadrupoleKey::addFields(
641  schema, schema.join(name, "psf"), "adaptive moments of the PSF model at the object position");
642  } else {
643  r._includePsf = false;
644  }
645 
646  r._instFlux_xx_Cov =
647  schema.addField<ErrElement>(schema.join(name, "instFlux", "xx", "Cov"),
648  (boost::format("uncertainty covariance between %s and %s") %
649  schema.join(name, "instFlux") % schema.join(name, "xx"))
650  .str(),
651  "count*pixel^2");
652  r._instFlux_yy_Cov =
653  schema.addField<ErrElement>(schema.join(name, "instFlux", "yy", "Cov"),
654  (boost::format("uncertainty covariance between %s and %s") %
655  schema.join(name, "instFlux") % schema.join(name, "yy"))
656  .str(),
657  "count*pixel^2");
658  r._instFlux_xy_Cov =
659  schema.addField<ErrElement>(schema.join(name, "instFlux", "xy", "Cov"),
660  (boost::format("uncertainty covariance between %s and %s") %
661  schema.join(name, "instFlux") % schema.join(name, "xy"))
662  .str(),
663  "count*pixel^2");
664 
665  // Skip the psf flag if not recording the PSF shape.
666  if (r._includePsf) {
668  } else {
671  }
672  return r;
673 }
674 
676  : _shapeResult(s),
677  _centroidResult(s),
678  _instFluxResult(s),
679  _instFlux_xx_Cov(s["instFlux"]["xx"]["Cov"]),
680  _instFlux_yy_Cov(s["instFlux"]["yy"]["Cov"]),
681  _instFlux_xy_Cov(s["instFlux"]["xy"]["Cov"]) {
682  // The input SubSchema may optionally provide for a PSF.
683  try {
684  _psfShapeResult = afw::table::QuadrupoleKey(s["psf"]);
686  _includePsf = true;
687  } catch (pex::exceptions::NotFoundError &e) {
688  _flagHandler =
690  _includePsf = false;
691  }
692 }
693 
696  static_cast<ShapeResult &>(result) = record.get(_shapeResult);
697  static_cast<CentroidResult &>(result) = record.get(_centroidResult);
698  static_cast<FluxResult &>(result) = record.get(_instFluxResult);
699  result.instFlux_xx_Cov = record.get(_instFlux_xx_Cov);
700  result.instFlux_yy_Cov = record.get(_instFlux_yy_Cov);
701  result.instFlux_xy_Cov = record.get(_instFlux_xy_Cov);
702  for (size_t n = 0; n < SdssShapeAlgorithm::N_FLAGS; ++n) {
703  if (n == SdssShapeAlgorithm::PSF_SHAPE_BAD.number && !_includePsf) continue;
704  result.flags[n] = _flagHandler.getValue(record, n);
705  }
706  return result;
707 }
708 
710  return record.get(_psfShapeResult);
711 }
712 
714  record.set(_shapeResult, value);
715  record.set(_centroidResult, value);
716  record.set(_instFluxResult, value);
717  record.set(_instFlux_xx_Cov, value.instFlux_xx_Cov);
718  record.set(_instFlux_yy_Cov, value.instFlux_yy_Cov);
719  record.set(_instFlux_xy_Cov, value.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  _flagHandler.setValue(record, n, value.flags[n]);
723  }
724 }
725 
727  afw::geom::ellipses::Quadrupole const &value) const {
728  record.set(_psfShapeResult, value);
729 }
730 
732  return _shapeResult == other._shapeResult && _centroidResult == other._centroidResult &&
733  _instFluxResult == other._instFluxResult && _psfShapeResult == other._psfShapeResult &&
734  _instFlux_xx_Cov == other._instFlux_xx_Cov && _instFlux_yy_Cov == other._instFlux_yy_Cov &&
735  _instFlux_xy_Cov == other._instFlux_xy_Cov;
736  // don't bother with flags - if we've gotten this far, it's basically impossible the flags don't match
737 }
738 
740  return _shapeResult.isValid() && _centroidResult.isValid() && _instFluxResult.isValid() &&
741  _psfShapeResult.isValid() && _instFlux_xx_Cov.isValid() && _instFlux_yy_Cov.isValid() &&
742  _instFlux_xy_Cov.isValid();
743  // don't bother with flags - if we've gotten this far, it's basically impossible the flags are invalid
744 }
745 
748  : _ctrl(ctrl),
749  _resultKey(ResultKey::addFields(schema, name, ctrl.doMeasurePsf)),
750  _centroidExtractor(schema, name) {}
751 
752 template <typename ImageT>
754  bool negative, Control const &control) {
755  double xcen = center.getX(); // object's column position
756  double ycen = center.getY(); // object's row position
757 
758  xcen -= image.getX0(); // work in image Pixel coordinates
759  ycen -= image.getY0();
760 
761  float shiftmax = control.maxShift; // Max allowed centroid shift
762  if (shiftmax < 2) {
763  shiftmax = 2;
764  } else if (shiftmax > 10) {
765  shiftmax = 10;
766  }
767 
769  try {
770  result.flags[FAILURE.number] =
771  !getAdaptiveMoments(image, control.background, xcen, ycen, shiftmax, &result, control.maxIter,
772  control.tol1, control.tol2, negative);
773  } catch (pex::exceptions::Exception &err) {
774  result.flags[FAILURE.number] = true;
775  }
776  if (result.flags[UNWEIGHTED.number] || result.flags[SHIFT.number]) {
777  // These are also considered fatal errors in terms of the quality of the results,
778  // even though they do produce some results.
779  result.flags[FAILURE.number] = true;
780  }
781  double IxxIyy = result.getQuadrupole().getIxx() * result.getQuadrupole().getIyy();
782  double Ixy_sq = result.getQuadrupole().getIxy();
783  Ixy_sq *= Ixy_sq;
784  double epsilon = 1.0e-6;
785  if (IxxIyy < (1.0 + epsilon) * Ixy_sq)
786  // We are checking that Ixx*Iyy > (1 + epsilon)*Ixy*Ixy where epsilon is suitably small. The
787  // value of epsilon used here is a magic number subject to future review (DM-5801 was marked won't fix).
788  {
789  if (!result.flags[FAILURE.number]) {
791  (boost::format("computeAdaptiveMoments IxxIxy %d < (1 + eps=%d)*(Ixy^2=%d);"
792  " implying singular moments without any flag set")
793  % IxxIyy % epsilon % Ixy_sq).str());
794  }
795  }
796 
797  // getAdaptiveMoments() just computes the zeroth moment in result.instFlux (and its error in
798  // result.instFluxErr, result.instFlux_xx_Cov, etc.) That's related to the instFlux by some geometric
799  // factors, which we apply here.
800  // This scale factor is just the inverse of the normalization constant in a bivariate normal distribution:
801  // 2*pi*sigma_x*sigma_y*sqrt(1-rho^2), where rho is the correlation.
802  // This happens to be twice the ellipse area pi*sigma_maj*sigma_min, which is identically equal to:
803  // pi*sqrt(I_xx*I_yy - I_xy^2); i.e. pi*sqrt(determinant(I))
804  double instFluxScale = geom::TWOPI * std::sqrt(IxxIyy - Ixy_sq);
805 
806  result.instFlux *= instFluxScale;
807  result.instFluxErr *= instFluxScale;
808  result.x += image.getX0();
809  result.y += image.getY0();
810 
811  if (ImageAdaptor<ImageT>::hasVariance) {
812  result.instFlux_xx_Cov *= instFluxScale;
813  result.instFlux_yy_Cov *= instFluxScale;
814  result.instFlux_xy_Cov *= instFluxScale;
815  }
816 
817  return result;
818 }
819 
820 template <typename ImageT>
822  afw::geom::ellipses::Quadrupole const &shape,
823  geom::Point2D const &center) {
824  // while arguments to computeFixedMomentsFlux are in PARENT coordinates, the implementation is LOCAL.
825  geom::Point2D localCenter = center - geom::Extent2D(image.getXY0());
826 
827  geom::BoxI const bbox = computeAdaptiveMomentsBBox(image.getBBox(afw::image::LOCAL), localCenter,
828  shape.getIxx(), shape.getIxy(), shape.getIyy());
829 
830  std::tuple<std::pair<bool, double>, double, double, double> weights =
831  getWeights(shape.getIxx(), shape.getIxy(), shape.getIyy());
832 
834 
835  if (!std::get<0>(weights).first) {
836  throw pex::exceptions::InvalidParameterError("Input shape is singular");
837  }
838 
839  double const w11 = std::get<1>(weights);
840  double const w12 = std::get<2>(weights);
841  double const w22 = std::get<3>(weights);
842  bool const interp = shouldInterp(shape.getIxx(), shape.getIyy(), std::get<0>(weights).second);
843 
844  double sum0 = 0; // sum of pixel values weighted by a Gaussian
845  if (calcmom<true>(ImageAdaptor<ImageT>().getImage(image), localCenter.getX(), localCenter.getY(), bbox,
846  0.0, interp, w11, w12, w22, NULL, &sum0, NULL, NULL, NULL, NULL, NULL, NULL) < 0) {
847  throw LSST_EXCEPT(pex::exceptions::RuntimeError, "Error from calcmom");
848  }
849 
850  result.instFlux = sum0 * 2.0;
851 
852  if (ImageAdaptor<ImageT>::hasVariance) {
853  int ix = static_cast<int>(center.getX() - image.getX0());
854  int iy = static_cast<int>(center.getY() - image.getY0());
855  if (!image.getBBox(afw::image::LOCAL).contains(geom::Point2I(ix, iy))) {
857  (boost::format("Center (%d,%d) not in image (%dx%d)") % ix % iy %
858  image.getWidth() % image.getHeight())
859  .str());
860  }
861  double var = ImageAdaptor<ImageT>().getVariance(image, ix, iy);
862  // 0th moment (i0) error = sqrt(var / wArea); instFlux (error) = 2 * wArea * i0 (error)
863  double const wArea = geom::PI * std::sqrt(shape.getDeterminant());
864  result.instFluxErr = 2 * std::sqrt(var * wArea);
865  }
866 
867  return result;
868 }
869 
871  afw::image::Exposure<float> const &exposure) const {
872  bool negative = false;
873 
874  try {
875  negative = measRecord.get(measRecord.getSchema().find<afw::table::Flag>("flags_negative").key);
876  } catch (pexExcept::Exception &e) {
877  }
879  exposure.getMaskedImage(), _centroidExtractor(measRecord, _resultKey.getFlagHandler()), negative,
880  _ctrl);
881 
882  if (_ctrl.doMeasurePsf) {
883  // Compute moments of Psf model. In the interest of implementing this quickly, we're just
884  // calling Psf::computeShape(), which delegates to SdssShapeResult::computeAdaptiveMoments
885  // for all nontrivial Psf classes. But this could in theory save the results of a shape
886  // computed some other way as part of base_SdssShape, which might be confusing. We should
887  // fix this eventually either by making Psf shape measurement not part of base_SdssShape, or
888  // by making the measurements stored with shape.sdss always computed via the
889  // SdssShapeAlgorithm instead of delegating to the Psf class.
890  try {
891  PTR(afw::detection::Psf const) psf = exposure.getPsf();
892  if (!psf) {
893  result.flags[PSF_SHAPE_BAD.number] = true;
894  } else {
895  _resultKey.setPsfShape(measRecord, psf->computeShape(geom::Point2D(result.x, result.y)));
896  }
897  } catch (pex::exceptions::Exception &err) {
898  result.flags[PSF_SHAPE_BAD.number] = true;
899  }
900  }
901 
902  measRecord.set(_resultKey, result);
903 }
904 
906  _resultKey.getFlagHandler().handleFailure(measRecord, error);
907 }
908 
909 #define INSTANTIATE_IMAGE(IMAGE) \
910  template SdssShapeResult SdssShapeAlgorithm::computeAdaptiveMoments( \
911  IMAGE const &, geom::Point2D const &, bool, Control const &); \
912  template FluxResult SdssShapeAlgorithm::computeFixedMomentsFlux( \
913  IMAGE const &, afw::geom::ellipses::Quadrupole const &, geom::Point2D const &)
914 
915 #define INSTANTIATE_PIXEL(PIXEL) \
916  INSTANTIATE_IMAGE(afw::image::Image<PIXEL>); \
917  INSTANTIATE_IMAGE(afw::image::MaskedImage<PIXEL>);
918 
919 INSTANTIATE_PIXEL(int);
920 INSTANTIATE_PIXEL(float);
921 INSTANTIATE_PIXEL(double);
922 
925  : BaseTransform{name}, _instFluxTransform{name, mapper}, _centroidTransform{name, mapper} {
926  // If the input schema has a PSF flag -- it's optional -- assume we are also transforming the PSF.
927  _transformPsf = mapper.getInputSchema().getNames().count("sdssShape_flag_psf") ? true : false;
928 
929  // Skip the last flag if not transforming the PSF shape.
930  for (std::size_t i = 0; i < SdssShapeAlgorithm::getFlagDefinitions().size(); i++) {
931  FlagDefinition const &flag = SdssShapeAlgorithm::getFlagDefinitions()[i];
932  if (flag == SdssShapeAlgorithm::FAILURE) continue;
933  if (mapper.getInputSchema().getNames().count(mapper.getInputSchema().join(name, flag.name)) == 0)
934  continue;
936  mapper.getInputSchema().find<afw::table::Flag>(name + "_" + flag.name).key;
937  mapper.addMapping(key);
938  }
939 
940  _outShapeKey = ShapeResultKey::addFields(mapper.editOutputSchema(), name, "Shape in celestial moments",
942  if (_transformPsf) {
943  _outPsfShapeKey = afw::table::QuadrupoleKey::addFields(mapper.editOutputSchema(), name + "_psf",
944  "PSF shape in celestial moments",
946  }
947 }
948 
950  afw::table::BaseCatalog &outputCatalog, afw::geom::SkyWcs const &wcs,
951  afw::image::PhotoCalib const &photoCalib) const {
952  // The instFlux and cetroid transforms will check that the catalog lengths
953  // match and throw if not, so we don't repeat the test here.
954  _instFluxTransform(inputCatalog, outputCatalog, wcs, photoCalib);
955  _centroidTransform(inputCatalog, outputCatalog, wcs, photoCalib);
956 
957  CentroidResultKey centroidKey(inputCatalog.getSchema()[_name]);
958  ShapeResultKey inShapeKey(inputCatalog.getSchema()[_name]);
959  afw::table::QuadrupoleKey inPsfShapeKey;
960  if (_transformPsf) {
961  inPsfShapeKey = afw::table::QuadrupoleKey(
962  inputCatalog.getSchema()[inputCatalog.getSchema().join(_name, "psf")]);
963  }
964 
965  afw::table::SourceCatalog::const_iterator inSrc = inputCatalog.begin();
966  afw::table::BaseCatalog::iterator outSrc = outputCatalog.begin();
967  for (; inSrc != inputCatalog.end(); ++inSrc, ++outSrc) {
968  ShapeResult inShape = inShapeKey.get(*inSrc);
969  ShapeResult outShape;
970 
971  // The transformation from the (x, y) to the (Ra, Dec) basis.
972  geom::AffineTransform crdTr =
973  wcs.linearizePixelToSky(centroidKey.get(*inSrc).getCentroid(), geom::radians);
974  outShape.setShape(inShape.getShape().transform(crdTr.getLinear()));
975 
976  // Transformation matrix from pixel to celestial basis.
978  outShape.setShapeErr((m * inShape.getShapeErr().cast<double>() * m.transpose()).cast<ErrElement>());
979 
980  _outShapeKey.set(*outSrc, outShape);
981 
982  if (_transformPsf) {
983  _outPsfShapeKey.set(*outSrc, inPsfShapeKey.get(*inSrc).transform(crdTr.getLinear()));
984  }
985  }
986 }
987 
988 } // namespace base
989 } // namespace meas
990 } // namespace lsst
lsst::afw::geom::ellipses::Quadrupole::getIxx
double const getIxx() const
Definition: Quadrupole.h:54
y
int y
Definition: SpanSet.cc:49
lsst::meas::base::FULL_COVARIANCE
@ FULL_COVARIANCE
The full covariance matrix is provided.
Definition: constants.h:46
lsst::meas::base::SdssShapeResult
Result object SdssShapeAlgorithm.
Definition: SdssShape.h:224
lsst::afw::image
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Definition: imageAlgorithm.dox:1
lsst::afw::image::Exposure::getPsf
std::shared_ptr< lsst::afw::detection::Psf const > getPsf() const
Return the Exposure's Psf object.
Definition: Exposure.h:322
SdssShape.h
lsst::meas::base::ShapeResult::xy
ShapeElement xy
image or model second moment for xy^2
Definition: ShapeUtilities.h:46
std::make_tuple
T make_tuple(T... args)
lsst::afw::image::LOCAL
@ LOCAL
Definition: ImageBase.h:94
lsst::afw::table::CoordinateType::PIXEL
@ PIXEL
lsst::geom::PI
constexpr double PI
The ratio of a circle's circumference to diameter.
Definition: Angle.h:39
lsst::meas::base::FlagDefinition::number
std::size_t number
Definition: FlagHandler.h:54
lsst::afw::table::QuadrupoleKey::set
void set(BaseRecord &record, geom::ellipses::Quadrupole const &value) const override
Set a Quadrupole in the given record.
Definition: aggregates.cc:114
std::string
STL class.
lsst::meas::base::FluxResult::instFlux
meas::base::Flux instFlux
Measured instFlux in DN.
Definition: FluxUtilities.h:42
Psf.h
lsst::meas::base::SdssShapeControl::tol2
float tol2
"Convergence tolerance for FWHM" ;
Definition: SdssShape.h:58
lsst::meas::base::FluxResult
A reusable result struct for instFlux measurements.
Definition: FluxUtilities.h:41
lsst::log.log.logContinued.error
def error(fmt, *args)
Definition: logContinued.py:213
std::fabs
T fabs(T... args)
lsst::afw::image::positionToIndex
int positionToIndex(double pos)
Convert image position to nearest integer index.
Definition: ImageUtils.h:69
lsst::meas::base::FlagHandler::handleFailure
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
lsst::afw::table::SourceRecord
Record class that contains measurements made on a single exposure.
Definition: Source.h:80
lsst::afw::table::BaseRecord::get
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Definition: BaseRecord.h:151
lsst::meas::base::NO_UNCERTAINTY
@ NO_UNCERTAINTY
Algorithm provides no uncertainy information at all.
Definition: constants.h:44
lsst::afw::image::Exposure< float >
wcs
table::Key< table::Array< std::uint8_t > > wcs
Definition: SkyWcs.cc:71
lsst::meas::base::SIGMA_ONLY
@ SIGMA_ONLY
Only the diagonal elements of the covariance matrix are provided.
Definition: constants.h:45
std::numeric_limits::quiet_NaN
T quiet_NaN(T... args)
lsst::meas::base::SdssShapeResult::flags
std::bitset< SdssShapeAlgorithm::N_FLAGS > flags
Status flags (see SdssShapeAlgorithm).
Definition: SdssShape.h:230
lsst::meas::base::SdssShapeControl::maxIter
int maxIter
"Maximum number of iterations" ;
Definition: SdssShape.h:55
lsst.pex::exceptions::NotFoundError
Reports attempts to access elements using an invalid key.
Definition: Runtime.h:151
psf
Key< int > psf
Definition: Exposure.cc:65
lsst::meas::base::FlagHandler::getValue
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
lsst::afw::table._match.first
first
Definition: _match.py:74
lsst::afw
Definition: imageAlgorithm.dox:1
ellipses.h
lsst::meas::base::SdssShapeResultKey::isValid
bool isValid() const
Return True if the key is valid.
Definition: SdssShape.cc:739
lsst::meas::base::MeasurementError
Exception to be thrown when a measurement algorithm experiences a known failure mode.
Definition: exceptions.h:48
lsst::afw::geom::ellipses::Quadrupole::getIxy
double const getIxy() const
Definition: Quadrupole.h:60
std::tuple
lsst::afw::table::Schema
Defines the fields and offsets for a table.
Definition: Schema.h:50
lsst::meas::base::ShapeResultKey
A FunctorKey for ShapeResult.
Definition: ShapeUtilities.h:113
lsst::afw::table::QuadrupoleKey::get
geom::ellipses::Quadrupole get(BaseRecord const &record) const override
Get a Quadrupole from the given record.
Definition: aggregates.cc:110
lsst.pex.config.history.format
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
lsst::afw::geom::SkyWcs
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Definition: SkyWcs.h:117
lsst::meas::base::CentroidResultKey
A FunctorKey for CentroidResult.
Definition: CentroidUtilities.h:88
lsst::meas::base::SdssShapeResultKey::operator==
bool operator==(SdssShapeResultKey const &other) const
Compare the FunctorKey for equality with another, using the underlying Keys.
Definition: SdssShape.cc:731
lsst.pipe.tasks.processCcdWithFakes.radius
radius
Definition: processCcdWithFakes.py:347
lsst::afw::geom.transform.transformContinued.name
string name
Definition: transformContinued.py:32
lsst::afw::image::X
@ X
Definition: ImageUtils.h:36
lsst::meas::base::SdssShapeResultKey::setPsfShape
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:726
lsst::meas::base::CentroidResult
A reusable struct for centroid measurements.
Definition: CentroidUtilities.h:41
lsst::afw::image::Exposure::getMaskedImage
MaskedImageT getMaskedImage()
Return the MaskedImage.
Definition: Exposure.h:228
std::sqrt
T sqrt(T... args)
lsst::meas::base::SdssShapeResult::instFlux_xx_Cov
ErrElement instFlux_xx_Cov
instFlux, xx term in the uncertainty covariance matrix
Definition: SdssShape.h:226
lsst::afw::table::Schema::find
SchemaItem< T > find(std::string const &name) const
Find a SchemaItem in the Schema by name.
Definition: Schema.cc:656
lsst::meas::base::FluxResultKey::isValid
bool isValid() const
Return True if both the instFlux and instFluxErr Keys are valid.
Definition: FluxUtilities.h:111
lsst::meas::base::CentroidResultKey::addFields
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.
Definition: CentroidUtilities.cc:66
lsst::meas::base::CentroidResultKey::get
virtual CentroidResult get(afw::table::BaseRecord const &record) const
Get a CentroidResult from the given record.
Definition: CentroidUtilities.cc:105
lsst::meas::base::FlagDefinitionList
vector-type utility class to build a collection of FlagDefinitions
Definition: FlagHandler.h:60
lsst::meas::base::ShapeTrMatrix
Eigen::Matrix< ShapeElement, 3, 3, Eigen::DontAlign > ShapeTrMatrix
Definition: constants.h:62
Angle.h
lsst::geom::AffineTransform
An affine coordinate transformation consisting of a linear transformation and an offset.
Definition: AffineTransform.h:75
lsst::meas::base::SdssShapeAlgorithm::computeAdaptiveMoments
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.
lsst::meas::base::SdssShapeAlgorithm::measure
virtual void measure(afw::table::SourceRecord &measRecord, afw::image::Exposure< float > const &exposure) const
Called to measure a single child source in an image.
Definition: SdssShape.cc:870
schema
table::Schema schema
Definition: python.h:134
lsst::meas::base::FlagHandler::setValue
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
std::isnan
T isnan(T... args)
INSTANTIATE_PIXEL
#define INSTANTIATE_PIXEL(PIXEL)
Definition: SdssShape.cc:915
lsst::meas::base::SdssShapeControl
A C++ control class to handle SdssShapeAlgorithm's configuration.
Definition: SdssShape.h:52
lsst::meas::base::SdssShapeAlgorithm::SHIFT
static FlagDefinition const SHIFT
Definition: SdssShape.h:158
lsst::afw::image::MaskedImage
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:73
lsst::meas::base::SdssShapeControl::tol1
float tol1
"Convergence tolerance for e1,e2" ;
Definition: SdssShape.h:57
lsst.pex::exceptions::DomainError
Reports arguments outside the domain of an operation.
Definition: Runtime.h:57
lsst::meas::base::SdssShapeTransform::operator()
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:949
lsst::meas::base::SdssShapeAlgorithm::UNWEIGHTED_BAD
static FlagDefinition const UNWEIGHTED_BAD
Definition: SdssShape.h:156
lsst::meas::base::SdssShapeAlgorithm::N_FLAGS
static unsigned int const N_FLAGS
Definition: SdssShape.h:154
std::printf
T printf(T... args)
image
afw::table::Key< afw::table::Array< ImagePixelT > > image
Definition: HeavyFootprint.cc:216
lsst::afw::table::QuadrupoleKey
A FunctorKey used to get or set a geom::ellipses::Quadrupole from a tuple of constituent Keys.
Definition: aggregates.h:282
lsst::meas::base::BaseTransform::_name
std::string _name
Definition: Transform.h:107
lsst::meas::base::FlagHandler
Utility class for handling flag fields that indicate the failure modes of an algorithm.
Definition: FlagHandler.h:148
x
double x
Definition: ChebyshevBoundedField.cc:277
lsst::afw::table::QuadrupoleKey::addFields
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
lsst::meas::base::ShapeResultKey::addFields
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.
Definition: ShapeUtilities.cc:74
lsst::geom::AffineTransform::getLinear
LinearTransform const & getLinear() const noexcept
Definition: AffineTransform.h:155
lsst::afw::geom::ellipses::Quadrupole::getIyy
double const getIyy() const
Definition: Quadrupole.h:57
lsst::afw::table._source.SourceCatalog
Definition: _source.py:32
other
ItemVariant const * other
Definition: Schema.cc:56
lsst::afw::table::BaseRecord
Base class for all records.
Definition: BaseRecord.h:31
lsst.pex::exceptions::LogicError
Reports errors in the logical structure of the program.
Definition: Runtime.h:46
lsst::afw::table::SchemaMapper
A mapping between the keys of two Schemas, used to copy data between them.
Definition: SchemaMapper.h:21
lsst::meas::base::SdssShapeResultKey::SdssShapeResultKey
SdssShapeResultKey()
Default constructor; instance will not be usuable unless subsequently assigned to.
Definition: SdssShape.h:90
lsst::meas::base::ShapeResult::setShapeErr
void setShapeErr(ShapeCov const &matrix)
Set the struct standard deviation elements from the given matrix, with rows and columns ordered (xx,...
Definition: ShapeUtilities.cc:56
lsst::meas::base::SdssShapeResultKey::set
virtual void set(afw::table::BaseRecord &record, SdssShapeResult const &value) const
Set an SdssShapeResult in the given record.
Definition: SdssShape.cc:713
lsst::afw::table::Key< afw::table::Flag >
LinearTransform.h
lsst::meas::base::ShapeResult::yy
ShapeElement yy
image or model second moment for y^2
Definition: ShapeUtilities.h:45
lsst::afw::table::Key::isValid
bool isValid() const noexcept
Return true if the key was initialized to valid offset.
Definition: Key.h:97
lsst::meas::base::CentroidResultKey::isValid
bool isValid() const
Return True if the centroid key is valid.
Definition: CentroidUtilities.h:138
lsst::afw::image::Y
@ Y
Definition: ImageUtils.h:36
lsst::afw::table::CatalogIterator
Iterator class for CatalogT.
Definition: Catalog.h:39
lsst.pipe.tasks.cli.cmd.commands.str
str
Definition: commands.py:50
ptr
uint64_t * ptr
Definition: RangeSet.cc:88
lsst::meas::base::SdssShapeControl::doMeasurePsf
bool doMeasurePsf
"Whether to also compute the shape of the PSF model" ;
Definition: SdssShape.h:59
image.h
lsst::afw::image::PhotoCalib
The photometric calibration of an exposure.
Definition: PhotoCalib.h:114
PTR
#define PTR(...)
Definition: base.h:41
lsst::geom::TWOPI
constexpr double TWOPI
Definition: Angle.h:40
lsst::meas::base::SdssShapeAlgorithm::UNWEIGHTED
static FlagDefinition const UNWEIGHTED
Definition: SdssShape.h:157
result
py::object result
Definition: _schema.cc:429
weight
afw::table::Key< double > weight
Definition: CoaddBoundedField.cc:166
Source.h
lsst::meas::base::ShapeResultKey::get
virtual ShapeResult get(afw::table::BaseRecord const &record) const
Get a ShapeResult from the given record.
Definition: ShapeUtilities.cc:127
exceptions.h
lsst::meas::base::SdssShapeAlgorithm::getFlagDefinitions
static FlagDefinitionList const & getFlagDefinitions()
Definition: SdssShape.cc:58
lsst::meas::base::SdssShapeResultKey::get
virtual SdssShapeResult get(afw::table::BaseRecord const &record) const
Get an SdssShapeResult from the given record.
Definition: SdssShape.cc:694
lsst
A base class for image defects.
Definition: imageAlgorithm.dox:1
lsst::meas::base::SdssShapeResult::instFlux_xy_Cov
ErrElement instFlux_xy_Cov
instFlux, xy term in the uncertainty covariance matrix
Definition: SdssShape.h:228
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
lsst::afw::detection
Definition: Footprint.h:50
std::min
T min(T... args)
lsst::afw::geom::ellipses::Quadrupole
An ellipse core with quadrupole moments as parameters.
Definition: Quadrupole.h:47
lsst::meas::base::FlagDefinitionList::size
std::size_t size() const
return the current size (number of defined elements) of the collection
Definition: FlagHandler.h:125
photoCalib
Key< int > photoCalib
Definition: Exposure.cc:67
lsst::meas::base::SdssShapeResultKey::getPsfShape
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:709
lsst::meas::base::SdssShapeResultKey::addFields
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:628
lsst::meas::base::SdssShapeAlgorithm::computeFixedMomentsFlux
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:821
lsst.pex::exceptions::InvalidParameterError
Reports invalid arguments.
Definition: Runtime.h:66
std::exp
T exp(T... args)
lsst::meas::base::SdssShapeAlgorithm::PSF_SHAPE_BAD
static FlagDefinition const PSF_SHAPE_BAD
Definition: SdssShape.h:160
lsst::meas::base::makeShapeTransformMatrix
ShapeTrMatrix makeShapeTransformMatrix(geom::LinearTransform const &xform)
Construct a matrix suitable for transforming second moments.
Definition: ShapeUtilities.cc:143
lsst::meas::base::ErrElement
float ErrElement
Definition: constants.h:55
lsst::afw::table::SubSchema
A proxy type for name lookups in a Schema.
Definition: Schema.h:357
std
STL namespace.
key
Key< U > key
Definition: Schema.cc:281
lsst::geom::Point< double, 2 >
lsst::afw::geom::ellipses::BaseCore::transform
Transformer transform(lsst::geom::LinearTransform const &transform)
Return the transform that maps the ellipse to the unit circle.
Definition: Transformer.h:116
lsst::meas::base::SdssShapeAlgorithm::FAILURE
static FlagDefinition const FAILURE
Definition: SdssShape.h:155
lsst::meas::base::ShapeResult::getShape
Shape const getShape() const
Return an afw::geom::ellipses object corresponding to xx, yy, xy.
Definition: ShapeUtilities.cc:41
lsst::meas::base::ShapeResult::getShapeErr
ShapeCov const getShapeErr() const
Return the 3x3 symmetric covariance matrix, with rows and columns ordered (xx, yy,...
Definition: ShapeUtilities.cc:49
lsst::afw::geom::ellipses::Quadrupole::getDeterminant
double getDeterminant() const
Return the determinant of the matrix representation.
Definition: Quadrupole.h:83
lsst::geom::Box2I
An integer coordinate rectangle.
Definition: Box.h:55
mapper
SchemaMapper * mapper
Definition: SchemaMapper.cc:78
lsst::meas::base::ShapeResult::setShape
void setShape(Shape const &shape)
Set struct elements from the given Quadrupole object.
Definition: ShapeUtilities.cc:43
lsst::geom::radians
constexpr AngleUnit radians
constant with units of radians
Definition: Angle.h:108
lsst::geom::Extent2D
Extent< double, 2 > Extent2D
Definition: Extent.h:400
lsst::meas::base::SdssShapeTransform::SdssShapeTransform
SdssShapeTransform(Control const &ctrl, std::string const &name, afw::table::SchemaMapper &mapper)
Definition: SdssShape.cc:923
lsst::meas::base::FluxResultKey::addFields
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...
Definition: FluxUtilities.cc:36
lsst::afw::table::BaseRecord::getSchema
Schema getSchema() const
Return the Schema that holds this record's fields and keys.
Definition: BaseRecord.h:80
std::size_t
lsst::meas::base::SdssShapeResultKey::getFlagHandler
FlagHandler const & getFlagHandler() const
Definition: SdssShape.h:125
std::make_pair
T make_pair(T... args)
lsst::meas::base::SdssShapeResult::SdssShapeResult
SdssShapeResult()
Constructor; initializes everything to NaN.
Definition: SdssShape.cc:623
lsst.pex::exceptions::Exception
Provides consistent interface for LSST exceptions.
Definition: Exception.h:107
lsst::afw::detection::Psf
A polymorphic base class for representing an image's Point Spread Function.
Definition: Psf.h:76
lsst::meas::base::ShapeResult
A reusable struct for moments-based shape measurements.
Definition: ShapeUtilities.h:43
lsst::meas::base::SdssShapeAlgorithm::MAXITER
static FlagDefinition const MAXITER
Definition: SdssShape.h:159
lsst::afw::table::BaseRecord::set
void set(Key< T > const &key, U const &value)
Set value of a field for the given key.
Definition: BaseRecord.h:164
lsst::meas::base::SdssShapeAlgorithm::fail
virtual void fail(afw::table::SourceRecord &measRecord, MeasurementError *error=nullptr) const
Handle an exception thrown by the current algorithm by setting flags in the given record.
Definition: SdssShape.cc:905
std::max
T max(T... args)
lsst::geom::Box2D
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
lsst::meas::base::ShapeResultKey::set
virtual void set(afw::table::BaseRecord &record, ShapeResult const &value) const
Set a ShapeResult in the given record.
Definition: ShapeUtilities.cc:136
lsst::meas::base::SdssShapeControl::maxShift
double maxShift
"Maximum centroid shift, limited to 2-10" ;
Definition: SdssShape.h:56
lsst::meas::base::SdssShapeResult::instFlux_yy_Cov
ErrElement instFlux_yy_Cov
instFlux, yy term in the uncertainty covariance matrix
Definition: SdssShape.h:227
lsst::meas::base::BaseTransform
Abstract base class for all C++ measurement transformations.
Definition: Transform.h:86
lsst::afw::table::CatalogT::begin
iterator begin()
Iterator access.
Definition: Catalog.h:396
lsst::meas::base::CentroidResult::getCentroid
Centroid const getCentroid() const
Return a Point object containing the measured x and y.
Definition: CentroidUtilities.cc:41
lsst::meas::base::SdssShapeResultKey
A FunctorKey that maps SdssShapeResult to afw::table Records.
Definition: SdssShape.h:73
lsst::meas::base::ShapeResult::xx
ShapeElement xx
image or model second moment for x^2
Definition: ShapeUtilities.h:44
lsst.obs.base.makeRawVisitInfo.NaN
NaN
Definition: makeRawVisitInfo.py:43
lsst::meas::base::SdssShapeAlgorithm::SdssShapeAlgorithm
SdssShapeAlgorithm(Control const &ctrl, std::string const &name, afw::table::Schema &schema)
Definition: SdssShape.cc:746
lsst::afw::table::CatalogT< BaseRecord >
lsst::geom::Extent< double, 2 >
std::numeric_limits
astshim.fitsChanContinued.iter
def iter(self)
Definition: fitsChanContinued.py:88
Box.h
lsst::afw::image::MaskedImage::Image
lsst::afw::image::Image< ImagePixelT > Image
Definition: MaskedImage.h:85
lsst::meas::base::ShapeResultKey::isValid
bool isValid() const
Return True if the shape key is valid.
Definition: ShapeUtilities.h:164
lsst::meas::base::FlagHandler::addFields
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
lsst::meas::base::SdssShapeControl::background
double background
"Additional value to add to background" ;
Definition: SdssShape.h:54
bbox
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
lsst.pex::exceptions::RuntimeError
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104
m
int m
Definition: SpanSet.cc:49
lsst::afw::table::QuadrupoleKey::isValid
bool isValid() const noexcept
Return True if all the constituent Keys are valid.
Definition: aggregates.h:342
std::pow
T pow(T... args)