LSSTApplications  19.0.0-14-gb0260a2+72efe9b372,20.0.0+7927753e06,20.0.0+8829bf0056,20.0.0+995114c5d2,20.0.0+b6f4b2abd1,20.0.0+bddc4f4cbe,20.0.0-1-g253301a+8829bf0056,20.0.0-1-g2b7511a+0d71a2d77f,20.0.0-1-g5b95a8c+7461dd0434,20.0.0-12-g321c96ea+23efe4bbff,20.0.0-16-gfab17e72e+fdf35455f6,20.0.0-2-g0070d88+ba3ffc8f0b,20.0.0-2-g4dae9ad+ee58a624b3,20.0.0-2-g61b8584+5d3db074ba,20.0.0-2-gb780d76+d529cf1a41,20.0.0-2-ged6426c+226a441f5f,20.0.0-2-gf072044+8829bf0056,20.0.0-2-gf1f7952+ee58a624b3,20.0.0-20-geae50cf+e37fec0aee,20.0.0-25-g3dcad98+544a109665,20.0.0-25-g5eafb0f+ee58a624b3,20.0.0-27-g64178ef+f1f297b00a,20.0.0-3-g4cc78c6+e0676b0dc8,20.0.0-3-g8f21e14+4fd2c12c9a,20.0.0-3-gbd60e8c+187b78b4b8,20.0.0-3-gbecbe05+48431fa087,20.0.0-38-ge4adf513+a12e1f8e37,20.0.0-4-g97dc21a+544a109665,20.0.0-4-gb4befbc+087873070b,20.0.0-4-gf910f65+5d3db074ba,20.0.0-5-gdfe0fee+199202a608,20.0.0-5-gfbfe500+d529cf1a41,20.0.0-6-g64f541c+d529cf1a41,20.0.0-6-g9a5b7a1+a1cd37312e,20.0.0-68-ga3f3dda+5fca18c6a4,20.0.0-9-g4aef684+e18322736b,w.2020.45
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
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  std::tuple<std::pair<bool, double>, double, double, double> const weights = getWeights(w11, w12, w22);
333  double const detW = std::get<1>(weights) * std::get<3>(weights) - std::pow(std::get<2>(weights), 2);
334  *pI0 = sum / (geom::PI * sqrt(detW));
335  if (psum) {
336  *psum = sum;
337  }
338  if (!instFluxOnly) {
339  *psumx = sumx;
340  *psumy = sumy;
341  *psumxx = sumxx;
342  *psumxy = sumxy;
343  *psumyy = sumyy;
344  if (psums4 != NULL) {
345  *psums4 = sums4;
346  }
347  }
348 
349 #if RECALC_W
350  if (wsum > 0 && !instFluxOnly) {
351  double det = w11 * w22 - w12 * w12;
352  wsumxx /= wsum;
353  wsumxy /= wsum;
354  wsumyy /= wsum;
355  printf("%g %g %g %g %g %g\n", w22 / det, -w12 / det, w11 / det, wsumxx, wsumxy, wsumyy);
356  }
357 #endif
358 
359  if (negative) {
360  return (instFluxOnly || (sum < 0 && sumxx < 0 && sumyy < 0)) ? 0 : -1;
361  } else {
362  return (instFluxOnly || (sum > 0 && sumxx > 0 && sumyy > 0)) ? 0 : -1;
363  }
364 }
365 
366 /*
367  * Workhorse for adaptive moments
368  *
369  * All inputs are expected to be in LOCAL image coordinates
370  */
371 template <typename ImageT>
372 bool getAdaptiveMoments(ImageT const &mimage, double bkgd, double xcen, double ycen, double shiftmax,
373  SdssShapeResult *shape, int maxIter, float tol1, float tol2, bool negative) {
374  double I0 = 0; // amplitude of best-fit Gaussian
375  double sum; // sum of intensity*weight
376  double sumx, sumy; // sum ((int)[xy])*intensity*weight
377  double sumxx, sumxy, sumyy; // sum {x^2,xy,y^2}*intensity*weight
378  double sums4; // sum intensity*weight*exponent^2
379  float const xcen0 = xcen; // initial centre
380  float const ycen0 = ycen; // of object
381 
382  double sigma11W = 1.5; // quadratic moments of the
383  double sigma12W = 0.0; // weighting fcn;
384  double sigma22W = 1.5; // xx, xy, and yy
385 
386  double w11 = -1, w12 = -1, w22 = -1; // current weights for moments; always set when iter == 0
387  float e1_old = 1e6, e2_old = 1e6; // old values of shape parameters e1 and e2
388  float sigma11_ow_old = 1e6; // previous version of sigma11_ow
389 
390  typename ImageAdaptor<ImageT>::Image const &image = ImageAdaptor<ImageT>().getImage(mimage);
391 
392  if (std::isnan(xcen) || std::isnan(ycen)) {
393  // Can't do anything
394  shape->flags[SdssShapeAlgorithm::UNWEIGHTED_BAD.number] = true;
395  return false;
396  }
397 
398  bool interpflag = false; // interpolate finer than a pixel?
400  int iter = 0; // iteration number
401  for (; iter < maxIter; iter++) {
402  bbox = computeAdaptiveMomentsBBox(image.getBBox(afw::image::LOCAL), geom::Point2D(xcen, ycen),
403  sigma11W, sigma12W, sigma22W);
404  std::tuple<std::pair<bool, double>, double, double, double> weights =
405  getWeights(sigma11W, sigma12W, sigma22W);
406  if (!std::get<0>(weights).first) {
407  shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
408  break;
409  }
410 
411  double const detW = std::get<0>(weights).second;
412 
413 #if 0 // this form was numerically unstable on my G4 powerbook
414  assert(detW >= 0.0);
415 #else
416  assert(sigma11W * sigma22W >= sigma12W * sigma12W - std::numeric_limits<float>::epsilon());
417 #endif
418 
419  {
420  const double ow11 = w11; // old
421  const double ow12 = w12; // values
422  const double ow22 = w22; // of w11, w12, w22
423 
424  w11 = std::get<1>(weights);
425  w12 = std::get<2>(weights);
426  w22 = std::get<3>(weights);
427 
428  if (shouldInterp(sigma11W, sigma22W, detW)) {
429  if (!interpflag) {
430  interpflag = true; // N.b.: stays set for this object
431  if (iter > 0) {
432  sigma11_ow_old = 1.e6; // force at least one more iteration
433  w11 = ow11;
434  w12 = ow12;
435  w22 = ow22;
436  iter--; // we didn't update wXX
437  }
438  }
439  }
440  }
441 
442  if (calcmom<false>(image, xcen, ycen, bbox, bkgd, interpflag, w11, w12, w22, &I0, &sum, &sumx, &sumy,
443  &sumxx, &sumxy, &sumyy, &sums4, negative) < 0) {
444  shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
445  break;
446  }
447 
448 #if 0
449 /*
450  * Find new centre
451  *
452  * This is only needed if we update the centre; if we use the input position we've already done the work
453  */
454  xcen = sumx/sum;
455  ycen = sumy/sum;
456 #endif
457  shape->x = sumx / sum; // update centroid. N.b. we're not setting errors here
458  shape->y = sumy / sum;
459 
460  if (fabs(shape->x - xcen0) > shiftmax || fabs(shape->y - ycen0) > shiftmax) {
461  shape->flags[SdssShapeAlgorithm::SHIFT.number] = true;
462  }
463  /*
464  * OK, we have the centre. Proceed to find the second moments.
465  */
466  float const sigma11_ow = sumxx / sum; // quadratic moments of
467  float const sigma22_ow = sumyy / sum; // weight*object
468  float const sigma12_ow = sumxy / sum; // xx, xy, and yy
469 
470  if (sigma11_ow <= 0 || sigma22_ow <= 0) {
471  shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
472  break;
473  }
474 
475  float const d = sigma11_ow + sigma22_ow; // current values of shape parameters
476  float const e1 = (sigma11_ow - sigma22_ow) / d;
477  float const e2 = 2.0 * sigma12_ow / d;
478  /*
479  * Did we converge?
480  */
481  if (iter > 0 && fabs(e1 - e1_old) < tol1 && fabs(e2 - e2_old) < tol1 &&
482  fabs(sigma11_ow / sigma11_ow_old - 1.0) < tol2) {
483  break; // yes; we converged
484  }
485 
486  e1_old = e1;
487  e2_old = e2;
488  sigma11_ow_old = sigma11_ow;
489  /*
490  * Didn't converge, calculate new values for weighting function
491  *
492  * The product of two Gaussians is a Gaussian:
493  * <x^2 exp(-a x^2 - 2bxy - cy^2) exp(-Ax^2 - 2Bxy - Cy^2)> =
494  * <x^2 exp(-(a + A) x^2 - 2(b + B)xy - (c + C)y^2)>
495  * i.e. the inverses of the covariances matrices add.
496  *
497  * We know sigmaXX_ow and sigmaXXW, the covariances of the weighted object
498  * and of the weights themselves. We can estimate the object's covariance as
499  * sigmaXX_ow^-1 - sigmaXXW^-1
500  * and, as we want to find a set of weights with the _same_ covariance as the
501  * object we take this to be the an estimate of our correct weights.
502  *
503  * N.b. This assumes that the object is roughly Gaussian.
504  * Consider the object:
505  * O == delta(x + p) + delta(x - p)
506  * the covariance of the weighted object is equal to that of the unweighted
507  * object, and this prescription fails badly. If we detect this, we set
508  * the UNWEIGHTED flag, and calculate the UNweighted moments
509  * instead.
510  */
511  {
512  float n11, n12, n22; // elements of inverse of next guess at weighting function
513  float ow11, ow12, ow22; // elements of inverse of sigmaXX_ow
514 
515  std::tuple<std::pair<bool, double>, double, double, double> weights =
516  getWeights(sigma11_ow, sigma12_ow, sigma22_ow);
517  if (!std::get<0>(weights).first) {
518  shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
519  break;
520  }
521 
522  ow11 = std::get<1>(weights);
523  ow12 = std::get<2>(weights);
524  ow22 = std::get<3>(weights);
525 
526  n11 = ow11 - w11;
527  n12 = ow12 - w12;
528  n22 = ow22 - w22;
529 
530  weights = getWeights(n11, n12, n22);
531  if (!std::get<0>(weights).first) {
532  // product-of-Gaussians assumption failed
533  shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
534  break;
535  }
536 
537  sigma11W = std::get<1>(weights);
538  sigma12W = std::get<2>(weights);
539  sigma22W = std::get<3>(weights);
540  }
541 
542  if (sigma11W <= 0 || sigma22W <= 0) {
543  shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
544  break;
545  }
546  }
547 
548  if (iter == maxIter) {
549  shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
550  shape->flags[SdssShapeAlgorithm::MAXITER.number] = true;
551  }
552 
553  if (sumxx + sumyy == 0.0) {
554  shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
555  }
556  /*
557  * Problems; try calculating the un-weighted moments
558  */
559  if (shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number]) {
560  w11 = w22 = w12 = 0;
561  if (calcmom<false>(image, xcen, ycen, bbox, bkgd, interpflag, w11, w12, w22, &I0, &sum, &sumx, &sumy,
562  &sumxx, &sumxy, &sumyy, NULL, negative) < 0 ||
563  (!negative && sum <= 0) || (negative && sum >= 0)) {
564  shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = false;
565  shape->flags[SdssShapeAlgorithm::UNWEIGHTED_BAD.number] = true;
566 
567  if (sum > 0) {
568  shape->xx = 1 / 12.0; // a single pixel
569  shape->xy = 0.0;
570  shape->yy = 1 / 12.0;
571  }
572 
573  return false;
574  }
575 
576  sigma11W = sumxx / sum; // estimate of object moments
577  sigma12W = sumxy / sum; // usually, object == weight
578  sigma22W = sumyy / sum; // at this point
579  }
580 
581  shape->instFlux = I0;
582  shape->xx = sigma11W;
583  shape->xy = sigma12W;
584  shape->yy = sigma22W;
585 
586  if (shape->xx + shape->yy != 0.0) {
587  int const ix = afw::image::positionToIndex(xcen);
588  int const iy = afw::image::positionToIndex(ycen);
589 
590  if (ix >= 0 && ix < mimage.getWidth() && iy >= 0 && iy < mimage.getHeight()) {
591  float const bkgd_var = ImageAdaptor<ImageT>().getVariance(
592  mimage, ix, iy); // XXX Overestimate as it includes object
593 
594  if (bkgd_var > 0.0) { // NaN is not > 0.0
595  if (!(shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number])) {
596  Matrix4d fisher = calc_fisher(*shape, bkgd_var); // Fisher matrix
597  Matrix4d cov = fisher.inverse();
598  // convention in afw::geom::ellipses is to order moments (xx, yy, xy),
599  // but the older algorithmic code uses (xx, xy, yy) - the order of
600  // indices here is not a bug.
601  shape->instFluxErr = std::sqrt(cov(0, 0));
602  shape->xxErr = std::sqrt(cov(1, 1));
603  shape->xyErr = std::sqrt(cov(2, 2));
604  shape->yyErr = std::sqrt(cov(3, 3));
605  shape->instFlux_xx_Cov = cov(0, 1);
606  shape->instFlux_xy_Cov = cov(0, 2);
607  shape->instFlux_yy_Cov = cov(0, 3);
608  shape->xx_yy_Cov = cov(1, 3);
609  shape->xx_xy_Cov = cov(1, 2);
610  shape->yy_xy_Cov = cov(2, 3);
611  }
612  }
613  }
614  }
615 
616  return true;
617 }
618 
619 } // namespace
620 
622  : instFlux_xx_Cov(std::numeric_limits<ErrElement>::quiet_NaN()),
623  instFlux_yy_Cov(std::numeric_limits<ErrElement>::quiet_NaN()),
624  instFlux_xy_Cov(std::numeric_limits<ErrElement>::quiet_NaN()) {}
625 
627  bool doMeasurePsf) {
629  r._shapeResult =
630  ShapeResultKey::addFields(schema, name, "elliptical Gaussian adaptive moments", SIGMA_ONLY);
631  r._centroidResult = CentroidResultKey::addFields(schema, name, "elliptical Gaussian adaptive moments",
633  r._instFluxResult = FluxResultKey::addFields(schema, name, "elliptical Gaussian adaptive moments");
634 
635  // Only include PSF shape fields if doMeasurePsf = True
636  if (doMeasurePsf) {
637  r._includePsf = true;
638  r._psfShapeResult = afw::table::QuadrupoleKey::addFields(
639  schema, schema.join(name, "psf"), "adaptive moments of the PSF model at the object position");
640  } else {
641  r._includePsf = false;
642  }
643 
644  r._instFlux_xx_Cov =
645  schema.addField<ErrElement>(schema.join(name, "instFlux", "xx", "Cov"),
646  (boost::format("uncertainty covariance between %s and %s") %
647  schema.join(name, "instFlux") % schema.join(name, "xx"))
648  .str(),
649  "count*pixel^2");
650  r._instFlux_yy_Cov =
651  schema.addField<ErrElement>(schema.join(name, "instFlux", "yy", "Cov"),
652  (boost::format("uncertainty covariance between %s and %s") %
653  schema.join(name, "instFlux") % schema.join(name, "yy"))
654  .str(),
655  "count*pixel^2");
656  r._instFlux_xy_Cov =
657  schema.addField<ErrElement>(schema.join(name, "instFlux", "xy", "Cov"),
658  (boost::format("uncertainty covariance between %s and %s") %
659  schema.join(name, "instFlux") % schema.join(name, "xy"))
660  .str(),
661  "count*pixel^2");
662 
663  // Skip the psf flag if not recording the PSF shape.
664  if (r._includePsf) {
666  } else {
669  }
670  return r;
671 }
672 
674  : _shapeResult(s),
675  _centroidResult(s),
676  _instFluxResult(s),
677  _instFlux_xx_Cov(s["instFlux"]["xx"]["Cov"]),
678  _instFlux_yy_Cov(s["instFlux"]["yy"]["Cov"]),
679  _instFlux_xy_Cov(s["instFlux"]["xy"]["Cov"]) {
680  // The input SubSchema may optionally provide for a PSF.
681  try {
682  _psfShapeResult = afw::table::QuadrupoleKey(s["psf"]);
684  _includePsf = true;
685  } catch (pex::exceptions::NotFoundError &e) {
686  _flagHandler =
688  _includePsf = false;
689  }
690 }
691 
694  static_cast<ShapeResult &>(result) = record.get(_shapeResult);
695  static_cast<CentroidResult &>(result) = record.get(_centroidResult);
696  static_cast<FluxResult &>(result) = record.get(_instFluxResult);
697  result.instFlux_xx_Cov = record.get(_instFlux_xx_Cov);
698  result.instFlux_yy_Cov = record.get(_instFlux_yy_Cov);
699  result.instFlux_xy_Cov = record.get(_instFlux_xy_Cov);
700  for (size_t n = 0; n < SdssShapeAlgorithm::N_FLAGS; ++n) {
701  if (n == SdssShapeAlgorithm::PSF_SHAPE_BAD.number && !_includePsf) continue;
702  result.flags[n] = _flagHandler.getValue(record, n);
703  }
704  return result;
705 }
706 
708  return record.get(_psfShapeResult);
709 }
710 
712  record.set(_shapeResult, value);
713  record.set(_centroidResult, value);
714  record.set(_instFluxResult, value);
715  record.set(_instFlux_xx_Cov, value.instFlux_xx_Cov);
716  record.set(_instFlux_yy_Cov, value.instFlux_yy_Cov);
717  record.set(_instFlux_xy_Cov, value.instFlux_xy_Cov);
718  for (size_t n = 0; n < SdssShapeAlgorithm::N_FLAGS; ++n) {
719  if (n == SdssShapeAlgorithm::PSF_SHAPE_BAD.number && !_includePsf) continue;
720  _flagHandler.setValue(record, n, value.flags[n]);
721  }
722 }
723 
725  afw::geom::ellipses::Quadrupole const &value) const {
726  record.set(_psfShapeResult, value);
727 }
728 
730  return _shapeResult == other._shapeResult && _centroidResult == other._centroidResult &&
731  _instFluxResult == other._instFluxResult && _psfShapeResult == other._psfShapeResult &&
732  _instFlux_xx_Cov == other._instFlux_xx_Cov && _instFlux_yy_Cov == other._instFlux_yy_Cov &&
733  _instFlux_xy_Cov == other._instFlux_xy_Cov;
734  // don't bother with flags - if we've gotten this far, it's basically impossible the flags don't match
735 }
736 
738  return _shapeResult.isValid() && _centroidResult.isValid() && _instFluxResult.isValid() &&
739  _psfShapeResult.isValid() && _instFlux_xx_Cov.isValid() && _instFlux_yy_Cov.isValid() &&
740  _instFlux_xy_Cov.isValid();
741  // don't bother with flags - if we've gotten this far, it's basically impossible the flags are invalid
742 }
743 
746  : _ctrl(ctrl),
747  _resultKey(ResultKey::addFields(schema, name, ctrl.doMeasurePsf)),
748  _centroidExtractor(schema, name) {}
749 
750 template <typename ImageT>
752  bool negative, Control const &control) {
753  double xcen = center.getX(); // object's column position
754  double ycen = center.getY(); // object's row position
755 
756  xcen -= image.getX0(); // work in image Pixel coordinates
757  ycen -= image.getY0();
758 
759  float shiftmax = control.maxShift; // Max allowed centroid shift
760  if (shiftmax < 2) {
761  shiftmax = 2;
762  } else if (shiftmax > 10) {
763  shiftmax = 10;
764  }
765 
767  try {
768  result.flags[FAILURE.number] =
769  !getAdaptiveMoments(image, control.background, xcen, ycen, shiftmax, &result, control.maxIter,
770  control.tol1, control.tol2, negative);
771  } catch (pex::exceptions::Exception &err) {
772  result.flags[FAILURE.number] = true;
773  }
774  if (result.flags[UNWEIGHTED.number] || result.flags[SHIFT.number]) {
775  // These are also considered fatal errors in terms of the quality of the results,
776  // even though they do produce some results.
777  result.flags[FAILURE.number] = true;
778  }
779  double IxxIyy = result.getQuadrupole().getIxx() * result.getQuadrupole().getIyy();
780  double Ixy_sq = result.getQuadrupole().getIxy();
781  Ixy_sq *= Ixy_sq;
782  double epsilon = 1.0e-6;
783  if (IxxIyy < (1.0 + epsilon) * Ixy_sq)
784  // We are checking that Ixx*Iyy > (1 + epsilon)*Ixy*Ixy where epsilon is suitably small. The
785  // value of epsilon used here is a magic number subject to future review (DM-5801 was marked won't fix).
786  {
787  if (!result.flags[FAILURE.number]) {
789  (boost::format("computeAdaptiveMoments IxxIxy %d < (1 + eps=%d)*(Ixy^2=%d);"
790  " implying singular moments without any flag set")
791  % IxxIyy % epsilon % Ixy_sq).str());
792  }
793  }
794 
795  // getAdaptiveMoments() just computes the zeroth moment in result.instFlux (and its error in
796  // result.instFluxErr, result.instFlux_xx_Cov, etc.) That's related to the instFlux by some geometric
797  // factors, which we apply here.
798  // This scale factor is just the inverse of the normalization constant in a bivariate normal distribution:
799  // 2*pi*sigma_x*sigma_y*sqrt(1-rho^2), where rho is the correlation.
800  // This happens to be twice the ellipse area pi*sigma_maj*sigma_min, which is identically equal to:
801  // pi*sqrt(I_xx*I_yy - I_xy^2); i.e. pi*sqrt(determinant(I))
802  double instFluxScale = geom::TWOPI * std::sqrt(IxxIyy - Ixy_sq);
803 
804  result.instFlux *= instFluxScale;
805  result.instFluxErr *= instFluxScale;
806  result.x += image.getX0();
807  result.y += image.getY0();
808 
809  if (ImageAdaptor<ImageT>::hasVariance) {
810  result.instFlux_xx_Cov *= instFluxScale;
811  result.instFlux_yy_Cov *= instFluxScale;
812  result.instFlux_xy_Cov *= instFluxScale;
813  }
814 
815  return result;
816 }
817 
818 template <typename ImageT>
820  afw::geom::ellipses::Quadrupole const &shape,
821  geom::Point2D const &center) {
822  // while arguments to computeFixedMomentsFlux are in PARENT coordinates, the implementation is LOCAL.
823  geom::Point2D localCenter = center - geom::Extent2D(image.getXY0());
824 
825  geom::BoxI const bbox = computeAdaptiveMomentsBBox(image.getBBox(afw::image::LOCAL), localCenter,
826  shape.getIxx(), shape.getIxy(), shape.getIyy());
827 
828  std::tuple<std::pair<bool, double>, double, double, double> weights =
829  getWeights(shape.getIxx(), shape.getIxy(), shape.getIyy());
830 
832 
833  if (!std::get<0>(weights).first) {
834  throw pex::exceptions::InvalidParameterError("Input shape is singular");
835  }
836 
837  double const w11 = std::get<1>(weights);
838  double const w12 = std::get<2>(weights);
839  double const w22 = std::get<3>(weights);
840  bool const interp = shouldInterp(shape.getIxx(), shape.getIyy(), std::get<0>(weights).second);
841 
842  double i0 = 0; // amplitude of Gaussian
843  if (calcmom<true>(ImageAdaptor<ImageT>().getImage(image), localCenter.getX(), localCenter.getY(), bbox,
844  0.0, interp, w11, w12, w22, &i0, NULL, NULL, NULL, NULL, NULL, NULL, NULL) < 0) {
845  throw LSST_EXCEPT(pex::exceptions::RuntimeError, "Error from calcmom");
846  }
847 
848  double const wArea = geom::PI * std::sqrt(shape.getDeterminant());
849 
850  result.instFlux = i0 * 2 * wArea;
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  result.instFluxErr = 2 * std::sqrt(var * wArea);
864  }
865 
866  return result;
867 }
868 
870  afw::image::Exposure<float> const &exposure) const {
871  bool negative = false;
872 
873  try {
874  negative = measRecord.get(measRecord.getSchema().find<afw::table::Flag>("flags_negative").key);
875  } catch (pexExcept::Exception &e) {
876  }
878  exposure.getMaskedImage(), _centroidExtractor(measRecord, _resultKey.getFlagHandler()), negative,
879  _ctrl);
880 
881  if (_ctrl.doMeasurePsf) {
882  // Compute moments of Psf model. In the interest of implementing this quickly, we're just
883  // calling Psf::computeShape(), which delegates to SdssShapeResult::computeAdaptiveMoments
884  // for all nontrivial Psf classes. But this could in theory save the results of a shape
885  // computed some other way as part of base_SdssShape, which might be confusing. We should
886  // fix this eventually either by making Psf shape measurement not part of base_SdssShape, or
887  // by making the measurements stored with shape.sdss always computed via the
888  // SdssShapeAlgorithm instead of delegating to the Psf class.
889  try {
890  PTR(afw::detection::Psf const) psf = exposure.getPsf();
891  if (!psf) {
892  result.flags[PSF_SHAPE_BAD.number] = true;
893  } else {
894  _resultKey.setPsfShape(measRecord, psf->computeShape(geom::Point2D(result.x, result.y)));
895  }
896  } catch (pex::exceptions::Exception &err) {
897  result.flags[PSF_SHAPE_BAD.number] = true;
898  }
899  }
900 
901  measRecord.set(_resultKey, result);
902 }
903 
905  _resultKey.getFlagHandler().handleFailure(measRecord, error);
906 }
907 
908 #define INSTANTIATE_IMAGE(IMAGE) \
909  template SdssShapeResult SdssShapeAlgorithm::computeAdaptiveMoments( \
910  IMAGE const &, geom::Point2D const &, bool, Control const &); \
911  template FluxResult SdssShapeAlgorithm::computeFixedMomentsFlux( \
912  IMAGE const &, afw::geom::ellipses::Quadrupole const &, geom::Point2D const &)
913 
914 #define INSTANTIATE_PIXEL(PIXEL) \
915  INSTANTIATE_IMAGE(afw::image::Image<PIXEL>); \
916  INSTANTIATE_IMAGE(afw::image::MaskedImage<PIXEL>);
917 
918 INSTANTIATE_PIXEL(int);
919 INSTANTIATE_PIXEL(float);
920 INSTANTIATE_PIXEL(double);
921 
924  : BaseTransform{name}, _instFluxTransform{name, mapper}, _centroidTransform{name, mapper} {
925  // If the input schema has a PSF flag -- it's optional -- assume we are also transforming the PSF.
926  _transformPsf = mapper.getInputSchema().getNames().count("sdssShape_flag_psf") ? true : false;
927 
928  // Skip the last flag if not transforming the PSF shape.
929  for (std::size_t i = 0; i < SdssShapeAlgorithm::getFlagDefinitions().size(); i++) {
930  FlagDefinition const &flag = SdssShapeAlgorithm::getFlagDefinitions()[i];
931  if (flag == SdssShapeAlgorithm::FAILURE) continue;
932  if (mapper.getInputSchema().getNames().count(mapper.getInputSchema().join(name, flag.name)) == 0)
933  continue;
935  mapper.getInputSchema().find<afw::table::Flag>(name + "_" + flag.name).key;
936  mapper.addMapping(key);
937  }
938 
939  _outShapeKey = ShapeResultKey::addFields(mapper.editOutputSchema(), name, "Shape in celestial moments",
941  if (_transformPsf) {
942  _outPsfShapeKey = afw::table::QuadrupoleKey::addFields(mapper.editOutputSchema(), name + "_psf",
943  "PSF shape in celestial moments",
945  }
946 }
947 
949  afw::table::BaseCatalog &outputCatalog, afw::geom::SkyWcs const &wcs,
950  afw::image::PhotoCalib const &photoCalib) const {
951  // The instFlux and cetroid transforms will check that the catalog lengths
952  // match and throw if not, so we don't repeat the test here.
953  _instFluxTransform(inputCatalog, outputCatalog, wcs, photoCalib);
954  _centroidTransform(inputCatalog, outputCatalog, wcs, photoCalib);
955 
956  CentroidResultKey centroidKey(inputCatalog.getSchema()[_name]);
957  ShapeResultKey inShapeKey(inputCatalog.getSchema()[_name]);
958  afw::table::QuadrupoleKey inPsfShapeKey;
959  if (_transformPsf) {
960  inPsfShapeKey = afw::table::QuadrupoleKey(
961  inputCatalog.getSchema()[inputCatalog.getSchema().join(_name, "psf")]);
962  }
963 
964  afw::table::SourceCatalog::const_iterator inSrc = inputCatalog.begin();
965  afw::table::BaseCatalog::iterator outSrc = outputCatalog.begin();
966  for (; inSrc != inputCatalog.end(); ++inSrc, ++outSrc) {
967  ShapeResult inShape = inShapeKey.get(*inSrc);
968  ShapeResult outShape;
969 
970  // The transformation from the (x, y) to the (Ra, Dec) basis.
971  geom::AffineTransform crdTr =
972  wcs.linearizePixelToSky(centroidKey.get(*inSrc).getCentroid(), geom::radians);
973  outShape.setShape(inShape.getShape().transform(crdTr.getLinear()));
974 
975  // Transformation matrix from pixel to celestial basis.
977  outShape.setShapeErr((m * inShape.getShapeErr().cast<double>() * m.transpose()).cast<ErrElement>());
978 
979  _outShapeKey.set(*outSrc, outShape);
980 
981  if (_transformPsf) {
982  _outPsfShapeKey.set(*outSrc, inPsfShapeKey.get(*inSrc).transform(crdTr.getLinear()));
983  }
984  }
985 }
986 
987 } // namespace base
988 } // namespace meas
989 } // namespace lsst
lsst::afw::geom::ellipses::Quadrupole::getIxx
double const getIxx() const
Definition: Quadrupole.h:54
y
int y
Definition: SpanSet.cc:49
schema
table::Schema schema
Definition: Amplifier.cc:115
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:307
SdssShape.h
lsst::meas::base::ShapeResult::xy
ShapeElement xy
image or model second moment for xy^2
Definition: ShapeUtilities.h:46
lsst::afw::table::CoordinateType::CELESTIAL
@ CELESTIAL
std::make_tuple
T make_tuple(T... args)
lsst::afw::image::LOCAL
@ LOCAL
Definition: ImageBase.h:94
lsst::geom::PI
constexpr double PI
The ratio of a circle's circumference to diameter.
Definition: Angle.h:39
lsst::meas::base::FlagDefinition::number
std::size_t number
Definition: FlagHandler.h:54
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 >
Box.h
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:737
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:729
lsst.pipe.tasks.processCcdWithFakes.radius
radius
Definition: processCcdWithFakes.py:345
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:724
lsst::meas::base::CentroidResult
A reusable struct for centroid measurements.
Definition: CentroidUtilities.h:41
PTR
#define PTR(...)
Definition: base.h:41
lsst::afw::image::Exposure::getMaskedImage
MaskedImageT getMaskedImage()
Return the MaskedImage.
Definition: Exposure.h:228
Angle.h
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
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:869
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:914
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:948
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
exceptions.h
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:711
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
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
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
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:692
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:707
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:626
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:819
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:922
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:621
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:904
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:744
lsst::afw::table::CatalogT< BaseRecord >
lsst::geom::Extent< double, 2 >
std::numeric_limits
astshim.fitsChanContinued.iter
def iter(self)
Definition: fitsChanContinued.py:88
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)