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