LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
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/format.hpp"
28#include "Eigen/LU"
29#include "lsst/geom/Angle.h"
30#include "lsst/geom/Box.h"
32#include "lsst/afw/image.h"
38
39namespace lsst {
40namespace meas {
41namespace base {
42namespace {
43FlagDefinitionList flagDefinitions;
44} // namespace
45
46FlagDefinition const SdssShapeAlgorithm::FAILURE = flagDefinitions.addFailureFlag();
47FlagDefinition const SdssShapeAlgorithm::UNWEIGHTED_BAD =
48 flagDefinitions.add("flag_unweightedBad", "Both weighted and unweighted moments were invalid");
49FlagDefinition const SdssShapeAlgorithm::UNWEIGHTED = flagDefinitions.add(
50 "flag_unweighted", "Weighted moments converged to an invalid value; using unweighted moments");
51FlagDefinition const SdssShapeAlgorithm::SHIFT =
52 flagDefinitions.add("flag_shift", "centroid shifted by more than the maximum allowed amount");
53FlagDefinition const SdssShapeAlgorithm::MAXITER =
54 flagDefinitions.add("flag_maxIter", "Too many iterations in adaptive moments");
55FlagDefinition const SdssShapeAlgorithm::PSF_SHAPE_BAD =
56 flagDefinitions.add("flag_psf", "Failure in measuring PSF model shape");
57
59
60namespace { // anonymous
61
62typedef 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 */
76Matrix4d 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//
128template <typename ImageT> // general case
129struct 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
139template <typename T> // specialise to a MaskedImage
140struct 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
153std::tuple<std::pair<bool, double>, double, double, double> getWeights(double sigma11, double sigma12,
154 double sigma22) {
155 double const NaN = std::numeric_limits<double>::quiet_NaN();
156 if (std::isnan(sigma11) || std::isnan(sigma12) || std::isnan(sigma22)) {
157 return std::make_tuple(std::make_pair(false, NaN), NaN, NaN, NaN);
158 }
159 double const det = sigma11 * sigma22 - sigma12 * sigma12; // determinant of sigmaXX matrix
160 if (std::isnan(det) || det < std::numeric_limits<float>::epsilon()) { // a suitably small number
161 return std::make_tuple(std::make_pair(false, det), NaN, NaN, NaN);
162 }
163 return std::make_tuple(std::make_pair(true, det), sigma22 / det, -sigma12 / det, sigma11 / det);
164}
165
167bool 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).
175geom::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 */
193template <typename ImageT>
194static void calcmom(ImageT const &image, // the image data
195 float xcen, float ycen, // centre of object
196 const 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 &psum) { // sum w*I
201
202 float tmod, ymod;
203 float X, Y; // sub-pixel interpolated [xy]
204 float weight;
205 float tmp;
206 double sum, sumx, sumy, sumxx, sumyy, sumxy, sums4;
207
208 if (w11 < 0 || w11 > 1e6 || fabs(w12) > 1E6 || w22 < 0 || w22 > 1e6) {
209 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "Invalid weight parameter(s)");
210 }
211
212 sum = sumx = sumy = sumxx = sumxy = sumyy = sums4 = 0;
213
214 int const ix0 = bbox.getMinX(); // corners of the box being analyzed
215 int const ix1 = bbox.getMaxX();
216 int const iy0 = bbox.getMinY(); // corners of the box being analyzed
217 int const iy1 = bbox.getMaxY();
218
219 if (ix0 < 0 || ix1 >= image.getWidth() || iy0 < 0 || iy1 >= image.getHeight()) {
220 throw LSST_EXCEPT(pex::exceptions::LengthError, "Invalid image dimensions");
221 }
222
223 for (int i = iy0; i <= iy1; ++i) {
224 typename ImageT::x_iterator ptr = image.x_at(ix0, i);
225 float const y = i - ycen;
226 float const y2 = y * y;
227 float const yl = y - 0.375;
228 float const yh = y + 0.375;
229 for (int j = ix0; j <= ix1; ++j, ++ptr) {
230 float x = j - xcen;
231 if (interpflag) {
232 float const xl = x - 0.375;
233 float const xh = x + 0.375;
234
235 float expon = xl * xl * w11 + yl * yl * w22 + 2.0 * xl * yl * w12;
236 tmp = xh * xh * w11 + yh * yh * w22 + 2.0 * xh * yh * w12;
237 expon = (expon > tmp) ? expon : tmp;
238 tmp = xl * xl * w11 + yh * yh * w22 + 2.0 * xl * yh * w12;
239 expon = (expon > tmp) ? expon : tmp;
240 tmp = xh * xh * w11 + yl * yl * w22 + 2.0 * xh * yl * w12;
241 expon = (expon > tmp) ? expon : tmp;
242
243 if (expon <= 9.0) {
244 tmod = *ptr - bkgd;
245 for (Y = yl; Y <= yh; Y += 0.25) {
246 double const interpY2 = Y * Y;
247 for (X = xl; X <= xh; X += 0.25) {
248 double const interpX2 = X * X;
249 double const interpXy = X * Y;
250 expon = interpX2 * w11 + 2 * interpXy * w12 + interpY2 * w22;
251 weight = std::exp(-0.5 * expon);
252
253 ymod = tmod * weight;
254 sum += ymod;
255 }
256 }
257 }
258 } else {
259 float x2 = x * x;
260 float xy = x * y;
261 float expon = x2 * w11 + 2 * xy * w12 + y2 * w22;
262
263 if (expon <= 14.0) {
264 weight = std::exp(-0.5 * expon);
265 tmod = *ptr - bkgd;
266 ymod = tmod * weight;
267 sum += ymod;
268 }
269 }
270 }
271 }
272
273 psum = sum;
274
275}
276
277template <typename ImageT>
278static int calcmom(ImageT const &image, // the image data
279 float xcen, float ycen, // centre of object
280 const geom::BoxI& bbox, // bounding box to consider
281 float bkgd, // data's background level
282 bool interpflag, // interpolate within pixels?
283 double w11, double w12, double w22, // weights
284 double &pI0, // amplitude of fit (if !NULL)
285 double &psum, // sum w*I (if !NULL)
286 double &psumx, double &psumy, // sum [xy]*w*I (if !instFluxOnly)
287 double &psumxx, double &psumxy, double &psumyy, // sum [xy]^2*w*I (if !instFluxOnly)
288 double &psums4, // sum w*I*weight^2 (if !instFluxOnly && !NULL)
289 bool negative = false) {
290 float tmod, ymod;
291 float X, Y; // sub-pixel interpolated [xy]
292 float weight;
293 float tmp;
294 double sum, sumx, sumy, sumxx, sumyy, sumxy, sums4;
295
296 if (w11 < 0 || w11 > 1e6 || fabs(w12) > 1E6 || w22 < 0 || w22 > 1e6) {
297 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "Invalid weight parameter(s)");
298 }
299
300 sum = sumx = sumy = sumxx = sumxy = sumyy = sums4 = 0;
301
302 int const ix0 = bbox.getMinX(); // corners of the box being analyzed
303 int const ix1 = bbox.getMaxX();
304 int const iy0 = bbox.getMinY(); // corners of the box being analyzed
305 int const iy1 = bbox.getMaxY();
306
307 if (ix0 < 0 || ix1 >= image.getWidth() || iy0 < 0 || iy1 >= image.getHeight()) {
308 throw LSST_EXCEPT(pex::exceptions::LengthError, "Invalid image dimensions");
309 }
310
311 for (int i = iy0; i <= iy1; ++i) {
312 typename ImageT::x_iterator ptr = image.x_at(ix0, i);
313 float const y = i - ycen;
314 float const y2 = y * y;
315 float const yl = y - 0.375;
316 float const yh = y + 0.375;
317 for (int j = ix0; j <= ix1; ++j, ++ptr) {
318 float x = j - xcen;
319 if (interpflag) {
320 float const xl = x - 0.375;
321 float const xh = x + 0.375;
322
323 float expon = xl * xl * w11 + yl * yl * w22 + 2.0 * xl * yl * w12;
324 tmp = xh * xh * w11 + yh * yh * w22 + 2.0 * xh * yh * w12;
325 expon = (expon > tmp) ? expon : tmp;
326 tmp = xl * xl * w11 + yh * yh * w22 + 2.0 * xl * yh * w12;
327 expon = (expon > tmp) ? expon : tmp;
328 tmp = xh * xh * w11 + yl * yl * w22 + 2.0 * xh * yl * w12;
329 expon = (expon > tmp) ? expon : tmp;
330
331 if (expon <= 9.0) {
332 tmod = *ptr - bkgd;
333 for (Y = yl; Y <= yh; Y += 0.25) {
334 double const interpY2 = Y * Y;
335 for (X = xl; X <= xh; X += 0.25) {
336 double const interpX2 = X * X;
337 double const interpXy = X * Y;
338 expon = interpX2 * w11 + 2 * interpXy * w12 + interpY2 * w22;
339 weight = std::exp(-0.5 * expon);
340
341 ymod = tmod * weight;
342 sum += ymod;
343 sumx += ymod * (X + xcen);
344 sumy += ymod * (Y + ycen);
345 sumxx += interpX2 * ymod;
346 sumxy += interpXy * ymod;
347 sumyy += interpY2 * ymod;
348 sums4 += expon * expon * ymod;
349 }
350 }
351 }
352 } else {
353 float x2 = x * x;
354 float xy = x * y;
355 float expon = x2 * w11 + 2 * xy * w12 + y2 * w22;
356
357 if (expon <= 14.0) {
358 weight = std::exp(-0.5 * expon);
359 tmod = *ptr - bkgd;
360 ymod = tmod * weight;
361 sum += ymod;
362 sumx += ymod * j;
363 sumy += ymod * i;
364 sumxx += x2 * ymod;
365 sumxy += xy * ymod;
366 sumyy += y2 * ymod;
367 sums4 += expon * expon * ymod;
368 }
369 }
370 }
371 }
372
373
374 std::tuple<std::pair<bool, double>, double, double, double> const weights = getWeights(w11, w12, w22);
375 double const detW = std::get<1>(weights) * std::get<3>(weights) - std::pow(std::get<2>(weights), 2);
376 pI0 = sum / (geom::PI * sqrt(detW));
377
378
379 psum = sum;
380
381 psumx = sumx;
382 psumy = sumy;
383 psumxx = sumxx;
384 psumxy = sumxy;
385 psumyy = sumyy;
386
387 psums4 = sums4;
388
389 if (negative) {
390 return (sum < 0 && sumxx < 0 && sumyy < 0) ? 0 : -1;
391 } else {
392 return (sum > 0 && sumxx > 0 && sumyy > 0) ? 0 : -1;
393 }
394}
395
396/*
397 * Workhorse for adaptive moments
398 *
399 * All inputs are expected to be in LOCAL image coordinates
400 */
401template <typename ImageT>
402bool getAdaptiveMoments(ImageT const &mimage, double bkgd, double xcen, double ycen, double shiftmax,
403 SdssShapeResult *shape, int maxIter, float tol1, float tol2, bool negative) {
404 double I0 = 0; // amplitude of best-fit Gaussian
405 double sum; // sum of intensity*weight
406 double sumx, sumy; // sum ((int)[xy])*intensity*weight
407 double sumxx, sumxy, sumyy; // sum {x^2,xy,y^2}*intensity*weight
408 double sums4; // sum intensity*weight*exponent^2
409 float const xcen0 = xcen; // initial centre
410 float const ycen0 = ycen; // of object
411
412 double sigma11W = 1.5; // quadratic moments of the
413 double sigma12W = 0.0; // weighting fcn;
414 double sigma22W = 1.5; // xx, xy, and yy
415
416 double w11 = -1, w12 = -1, w22 = -1; // current weights for moments; always set when iter == 0
417 float e1_old = 1e6, e2_old = 1e6; // old values of shape parameters e1 and e2
418 float sigma11_ow_old = 1e6; // previous version of sigma11_ow
419
420 typename ImageAdaptor<ImageT>::Image const &image = ImageAdaptor<ImageT>().getImage(mimage);
421
422 if (std::isnan(xcen) || std::isnan(ycen)) {
423 // Can't do anything
424 shape->flags[SdssShapeAlgorithm::UNWEIGHTED_BAD.number] = true;
425 return false;
426 }
427
428 bool interpflag = false; // interpolate finer than a pixel?
430 int iter = 0; // iteration number
431 for (; iter < maxIter; iter++) {
432 bbox = computeAdaptiveMomentsBBox(image.getBBox(afw::image::LOCAL), geom::Point2D(xcen, ycen),
433 sigma11W, sigma12W, sigma22W);
434 std::tuple<std::pair<bool, double>, double, double, double> weights =
435 getWeights(sigma11W, sigma12W, sigma22W);
436 if (!std::get<0>(weights).first) {
437 shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
438 break;
439 }
440
441 double const detW = std::get<0>(weights).second;
442
443 if( sigma11W * sigma22W < sigma12W * sigma12W - std::numeric_limits<float>::epsilon()) return false;
444
445 {
446 const double ow11 = w11; // old
447 const double ow12 = w12; // values
448 const double ow22 = w22; // of w11, w12, w22
449
450 w11 = std::get<1>(weights);
451 w12 = std::get<2>(weights);
452 w22 = std::get<3>(weights);
453
454 if (shouldInterp(sigma11W, sigma22W, detW)) {
455 if (!interpflag) {
456 interpflag = true; // N.b.: stays set for this object
457 if (iter > 0) {
458 sigma11_ow_old = 1.e6; // force at least one more iteration
459 w11 = ow11;
460 w12 = ow12;
461 w22 = ow22;
462 iter--; // we didn't update wXX
463 }
464 }
465 }
466 }
467 if (calcmom(image, xcen, ycen, bbox, bkgd, interpflag, w11, w12, w22, I0, sum, sumx, sumy,
468 sumxx, sumxy, sumyy, sums4, negative) < 0) {
469 shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
470 break;
471 }
472
473 shape->x = sumx / sum; // update centroid. N.b. we're not setting errors here
474 shape->y = sumy / sum;
475
476 if (fabs(shape->x - xcen0) > shiftmax || fabs(shape->y - ycen0) > shiftmax) {
477 shape->flags[SdssShapeAlgorithm::SHIFT.number] = true;
478 }
479 /*
480 * OK, we have the centre. Proceed to find the second moments.
481 */
482 float const sigma11_ow = sumxx / sum; // quadratic moments of
483 float const sigma22_ow = sumyy / sum; // weight*object
484 float const sigma12_ow = sumxy / sum; // xx, xy, and yy
485
486 if (sigma11_ow <= 0 || sigma22_ow <= 0) {
487 shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
488 break;
489 }
490
491 float const d = sigma11_ow + sigma22_ow; // current values of shape parameters
492 float const e1 = (sigma11_ow - sigma22_ow) / d;
493 float const e2 = 2.0 * sigma12_ow / d;
494 /*
495 * Did we converge?
496 */
497 if (iter > 0 && fabs(e1 - e1_old) < tol1 && fabs(e2 - e2_old) < tol1 &&
498 fabs(sigma11_ow / sigma11_ow_old - 1.0) < tol2) {
499 break; // yes; we converged
500 }
501
502 e1_old = e1;
503 e2_old = e2;
504 sigma11_ow_old = sigma11_ow;
505 /*
506 * Didn't converge, calculate new values for weighting function
507 *
508 * The product of two Gaussians is a Gaussian:
509 * <x^2 exp(-a x^2 - 2bxy - cy^2) exp(-Ax^2 - 2Bxy - Cy^2)> =
510 * <x^2 exp(-(a + A) x^2 - 2(b + B)xy - (c + C)y^2)>
511 * i.e. the inverses of the covariances matrices add.
512 *
513 * We know sigmaXX_ow and sigmaXXW, the covariances of the weighted object
514 * and of the weights themselves. We can estimate the object's covariance as
515 * sigmaXX_ow^-1 - sigmaXXW^-1
516 * and, as we want to find a set of weights with the _same_ covariance as the
517 * object we take this to be the an estimate of our correct weights.
518 *
519 * N.b. This assumes that the object is roughly Gaussian.
520 * Consider the object:
521 * O == delta(x + p) + delta(x - p)
522 * the covariance of the weighted object is equal to that of the unweighted
523 * object, and this prescription fails badly. If we detect this, we set
524 * the UNWEIGHTED flag, and calculate the UNweighted moments
525 * instead.
526 */
527 {
528 float n11, n12, n22; // elements of inverse of next guess at weighting function
529 float ow11, ow12, ow22; // elements of inverse of sigmaXX_ow
530
531 std::tuple<std::pair<bool, double>, double, double, double> weights =
532 getWeights(sigma11_ow, sigma12_ow, sigma22_ow);
533 if (!std::get<0>(weights).first) {
534 shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
535 break;
536 }
537
538 ow11 = std::get<1>(weights);
539 ow12 = std::get<2>(weights);
540 ow22 = std::get<3>(weights);
541
542 n11 = ow11 - w11;
543 n12 = ow12 - w12;
544 n22 = ow22 - w22;
545
546 weights = getWeights(n11, n12, n22);
547 if (!std::get<0>(weights).first) {
548 // product-of-Gaussians assumption failed
549 shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
550 break;
551 }
552
553 sigma11W = std::get<1>(weights);
554 sigma12W = std::get<2>(weights);
555 sigma22W = std::get<3>(weights);
556 }
557
558 if (sigma11W <= 0 || sigma22W <= 0) {
559 shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
560 break;
561 }
562 }
563
564 if (iter == maxIter) {
565 shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
566 shape->flags[SdssShapeAlgorithm::MAXITER.number] = true;
567 }
568
569 if (sumxx + sumyy == 0.0) {
570 shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = true;
571 }
572 /*
573 * Problems; try calculating the un-weighted moments
574 */
575 if (shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number]) {
576 w11 = w22 = w12 = 0;
577 double ignored;
578 if (calcmom(image, xcen, ycen, bbox, bkgd, interpflag, w11, w12, w22, I0, sum, sumx, sumy,
579 sumxx, sumxy, sumyy, ignored, negative) < 0 ||
580 (!negative && sum <= 0) || (negative && sum >= 0)) {
581 shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number] = false;
582 shape->flags[SdssShapeAlgorithm::UNWEIGHTED_BAD.number] = true;
583
584 if (sum > 0) {
585 shape->xx = 1 / 12.0; // a single pixel
586 shape->xy = 0.0;
587 shape->yy = 1 / 12.0;
588 }
589
590 return false;
591 }
592
593 sigma11W = sumxx / sum; // estimate of object moments
594 sigma12W = sumxy / sum; // usually, object == weight
595 sigma22W = sumyy / sum; // at this point
596 }
597
598 shape->instFlux = I0;
599 shape->xx = sigma11W;
600 shape->xy = sigma12W;
601 shape->yy = sigma22W;
602
603 if (shape->xx + shape->yy != 0.0) {
604 int const ix = afw::image::positionToIndex(xcen);
605 int const iy = afw::image::positionToIndex(ycen);
606
607 if (ix >= 0 && ix < mimage.getWidth() && iy >= 0 && iy < mimage.getHeight()) {
608 float const bkgd_var = ImageAdaptor<ImageT>().getVariance(
609 mimage, ix, iy); // XXX Overestimate as it includes object
610
611 if (bkgd_var > 0.0) { // NaN is not > 0.0
612 if (!(shape->flags[SdssShapeAlgorithm::UNWEIGHTED.number])) {
613 Matrix4d fisher = calc_fisher(*shape, bkgd_var); // Fisher matrix
614 Matrix4d cov = fisher.inverse();
615 // convention is to order moments (xx, yy, xy)
616 shape->instFluxErr = std::sqrt(cov(0, 0));
617 shape->xxErr = std::sqrt(cov(1, 1));
618 shape->yyErr = std::sqrt(cov(2, 2));
619 shape->xyErr = std::sqrt(cov(3, 3));
620 shape->instFlux_xx_Cov = cov(0, 1);
621 shape->instFlux_yy_Cov = cov(0, 2);
622 shape->instFlux_xy_Cov = cov(0, 3);
623 shape->xx_yy_Cov = cov(1, 2);
624 shape->xx_xy_Cov = cov(1, 3);
625 shape->yy_xy_Cov = cov(2, 3);
626 }
627 }
628 }
629 }
630
631 return true;
632}
633
634} // namespace
635
637 : instFlux_xx_Cov(std::numeric_limits<ErrElement>::quiet_NaN()),
638 instFlux_yy_Cov(std::numeric_limits<ErrElement>::quiet_NaN()),
639 instFlux_xy_Cov(std::numeric_limits<ErrElement>::quiet_NaN()) {}
640
642 bool doMeasurePsf) {
644 r._shapeResult =
645 ShapeResultKey::addFields(schema, name, "elliptical Gaussian adaptive moments", SIGMA_ONLY);
646 r._centroidResult = CentroidResultKey::addFields(schema, name, "elliptical Gaussian adaptive moments",
648 r._instFluxResult = FluxResultKey::addFields(schema, name, "elliptical Gaussian adaptive moments");
649
650 // Only include PSF shape fields if doMeasurePsf = True
651 if (doMeasurePsf) {
652 r._includePsf = true;
653 r._psfShapeResult = afw::table::QuadrupoleKey::addFields(
654 schema, schema.join(name, "psf"), "adaptive moments of the PSF model at the object position");
655 } else {
656 r._includePsf = false;
657 }
658
659 r._instFlux_xx_Cov =
660 schema.addField<ErrElement>(schema.join(name, "instFlux", "xx", "Cov"),
661 (boost::format("uncertainty covariance between %s and %s") %
662 schema.join(name, "instFlux") % schema.join(name, "xx"))
663 .str(),
664 "count*pixel^2");
665 r._instFlux_yy_Cov =
666 schema.addField<ErrElement>(schema.join(name, "instFlux", "yy", "Cov"),
667 (boost::format("uncertainty covariance between %s and %s") %
668 schema.join(name, "instFlux") % schema.join(name, "yy"))
669 .str(),
670 "count*pixel^2");
671 r._instFlux_xy_Cov =
672 schema.addField<ErrElement>(schema.join(name, "instFlux", "xy", "Cov"),
673 (boost::format("uncertainty covariance between %s and %s") %
674 schema.join(name, "instFlux") % schema.join(name, "xy"))
675 .str(),
676 "count*pixel^2");
677
678 // Skip the psf flag if not recording the PSF shape.
679 if (r._includePsf) {
680 r._flagHandler = FlagHandler::addFields(schema, name, SdssShapeAlgorithm::getFlagDefinitions());
681 } else {
682 r._flagHandler = FlagHandler::addFields(schema, name, SdssShapeAlgorithm::getFlagDefinitions(),
684 }
685 return r;
686}
687
689 : _shapeResult(s),
690 _centroidResult(s),
691 _instFluxResult(s),
692 _instFlux_xx_Cov(s["instFlux"]["xx"]["Cov"]),
693 _instFlux_yy_Cov(s["instFlux"]["yy"]["Cov"]),
694 _instFlux_xy_Cov(s["instFlux"]["xy"]["Cov"]) {
695 // The input SubSchema may optionally provide for a PSF.
696 try {
697 _psfShapeResult = afw::table::QuadrupoleKey(s["psf"]);
699 _includePsf = true;
700 } catch (pex::exceptions::NotFoundError &e) {
701 _flagHandler =
703 _includePsf = false;
704 }
705}
706
709 static_cast<ShapeResult &>(result) = record.get(_shapeResult);
710 static_cast<CentroidResult &>(result) = record.get(_centroidResult);
711 static_cast<FluxResult &>(result) = record.get(_instFluxResult);
712 result.instFlux_xx_Cov = record.get(_instFlux_xx_Cov);
713 result.instFlux_yy_Cov = record.get(_instFlux_yy_Cov);
714 result.instFlux_xy_Cov = record.get(_instFlux_xy_Cov);
715 for (size_t n = 0; n < SdssShapeAlgorithm::N_FLAGS; ++n) {
716 if (n == SdssShapeAlgorithm::PSF_SHAPE_BAD.number && !_includePsf) continue;
717 result.flags[n] = _flagHandler.getValue(record, n);
718 }
719 return result;
720}
721
723 return record.get(_psfShapeResult);
724}
725
727 record.set(_shapeResult, value);
728 record.set(_centroidResult, value);
729 record.set(_instFluxResult, value);
730 record.set(_instFlux_xx_Cov, value.instFlux_xx_Cov);
731 record.set(_instFlux_yy_Cov, value.instFlux_yy_Cov);
732 record.set(_instFlux_xy_Cov, value.instFlux_xy_Cov);
733 for (size_t n = 0; n < SdssShapeAlgorithm::N_FLAGS; ++n) {
734 if (n == SdssShapeAlgorithm::PSF_SHAPE_BAD.number && !_includePsf) continue;
735 _flagHandler.setValue(record, n, value.flags[n]);
736 }
737}
738
740 afw::geom::ellipses::Quadrupole const &value) const {
741 record.set(_psfShapeResult, value);
742}
743
745 return _shapeResult == other._shapeResult && _centroidResult == other._centroidResult &&
746 _instFluxResult == other._instFluxResult && _psfShapeResult == other._psfShapeResult &&
747 _instFlux_xx_Cov == other._instFlux_xx_Cov && _instFlux_yy_Cov == other._instFlux_yy_Cov &&
748 _instFlux_xy_Cov == other._instFlux_xy_Cov;
749 // don't bother with flags - if we've gotten this far, it's basically impossible the flags don't match
750}
751
753 return _shapeResult.isValid() && _centroidResult.isValid() && _instFluxResult.isValid() &&
754 _psfShapeResult.isValid() && _instFlux_xx_Cov.isValid() && _instFlux_yy_Cov.isValid() &&
755 _instFlux_xy_Cov.isValid();
756 // don't bother with flags - if we've gotten this far, it's basically impossible the flags are invalid
757}
758
760 afw::table::Schema &schema)
761 : _ctrl(ctrl),
762 _resultKey(ResultKey::addFields(schema, name, ctrl.doMeasurePsf)),
763 _centroidExtractor(schema, name) {}
764
765template <typename ImageT>
767 bool negative, Control const &control) {
768 double xcen = center.getX(); // object's column position
769 double ycen = center.getY(); // object's row position
770
771 xcen -= image.getX0(); // work in image Pixel coordinates
772 ycen -= image.getY0();
773
774 float shiftmax = control.maxShift; // Max allowed centroid shift
775 if (shiftmax < 2) {
776 shiftmax = 2;
777 } else if (shiftmax > 10) {
778 shiftmax = 10;
779 }
780
782 try {
784 !getAdaptiveMoments(image, control.background, xcen, ycen, shiftmax, &result, control.maxIter,
785 control.tol1, control.tol2, negative);
786 } catch (pex::exceptions::Exception &err) {
787 result.flags[FAILURE.number] = true;
788 }
789 if (result.flags[UNWEIGHTED.number] || result.flags[SHIFT.number]) {
790 // These are also considered fatal errors in terms of the quality of the results,
791 // even though they do produce some results.
792 result.flags[FAILURE.number] = true;
793 }
794 double IxxIyy = result.getQuadrupole().getIxx() * result.getQuadrupole().getIyy();
795 double Ixy_sq = result.getQuadrupole().getIxy();
796 Ixy_sq *= Ixy_sq;
797 double epsilon = 1.0e-6;
798 if (IxxIyy < (1.0 + epsilon) * Ixy_sq)
799 // We are checking that Ixx*Iyy > (1 + epsilon)*Ixy*Ixy where epsilon is suitably small. The
800 // value of epsilon used here is a magic number subject to future review (DM-5801 was marked won't fix).
801 {
802 if (!result.flags[FAILURE.number]) {
804 (boost::format("computeAdaptiveMoments IxxIxy %d < (1 + eps=%d)*(Ixy^2=%d);"
805 " implying singular moments without any flag set")
806 % IxxIyy % epsilon % Ixy_sq).str());
807 }
808 }
809
810 // getAdaptiveMoments() just computes the zeroth moment in result.instFlux (and its error in
811 // result.instFluxErr, result.instFlux_xx_Cov, etc.) That's related to the instFlux by some geometric
812 // factors, which we apply here.
813 // This scale factor is just the inverse of the normalization constant in a bivariate normal distribution:
814 // 2*pi*sigma_x*sigma_y*sqrt(1-rho^2), where rho is the correlation.
815 // This happens to be twice the ellipse area pi*sigma_maj*sigma_min, which is identically equal to:
816 // pi*sqrt(I_xx*I_yy - I_xy^2); i.e. pi*sqrt(determinant(I))
817 double instFluxScale = geom::TWOPI * std::sqrt(IxxIyy - Ixy_sq);
818
819 result.instFlux *= instFluxScale;
820 result.instFluxErr *= instFluxScale;
821 result.x += image.getX0();
822 result.y += image.getY0();
823
824 if (ImageAdaptor<ImageT>::hasVariance) {
825 result.instFlux_xx_Cov *= instFluxScale;
826 result.instFlux_yy_Cov *= instFluxScale;
827 result.instFlux_xy_Cov *= instFluxScale;
828 }
829
830 return result;
831}
832
833template <typename ImageT>
836 geom::Point2D const &center) {
837 // while arguments to computeFixedMomentsFlux are in PARENT coordinates, the implementation is LOCAL.
838 geom::Point2D localCenter = center - geom::Extent2D(image.getXY0());
839
840 geom::BoxI const bbox = computeAdaptiveMomentsBBox(image.getBBox(afw::image::LOCAL), localCenter,
841 shape.getIxx(), shape.getIxy(), shape.getIyy());
842
843 std::tuple<std::pair<bool, double>, double, double, double> weights =
844 getWeights(shape.getIxx(), shape.getIxy(), shape.getIyy());
845
847
848 if (!std::get<0>(weights).first) {
849 throw pex::exceptions::InvalidParameterError("Input shape is singular");
850 }
851
852 double const w11 = std::get<1>(weights);
853 double const w12 = std::get<2>(weights);
854 double const w22 = std::get<3>(weights);
855 bool const interp = shouldInterp(shape.getIxx(), shape.getIyy(), std::get<0>(weights).second);
856
857 double sum0 = 0; // sum of pixel values weighted by a Gaussian
858 calcmom(ImageAdaptor<ImageT>().getImage(image), localCenter.getX(), localCenter.getY(), bbox,
859 0.0, interp, w11, w12, w22, sum0);
860
861 result.instFlux = sum0 * 2.0;
862
863 if (ImageAdaptor<ImageT>::hasVariance) {
864 int ix = static_cast<int>(center.getX() - image.getX0());
865 int iy = static_cast<int>(center.getY() - image.getY0());
866 if (!image.getBBox(afw::image::LOCAL).contains(geom::Point2I(ix, iy))) {
868 (boost::format("Center (%d,%d) not in image (%dx%d)") % ix % iy %
869 image.getWidth() % image.getHeight())
870 .str());
871 }
872 double var = ImageAdaptor<ImageT>().getVariance(image, ix, iy);
873 // 0th moment (i0) error = sqrt(var / wArea); instFlux (error) = 2 * wArea * i0 (error)
874 double const wArea = geom::PI * std::sqrt(shape.getDeterminant());
875 result.instFluxErr = 2 * std::sqrt(var * wArea);
876 }
877
878 return result;
879}
880
882 afw::image::Exposure<float> const &exposure) const {
883 bool negative = false;
884
885 try {
886 negative = measRecord.get(measRecord.getSchema().find<afw::table::Flag>("flags_negative").key);
887 } catch (pexExcept::Exception &e) {
888 }
890 exposure.getMaskedImage(), _centroidExtractor(measRecord, _resultKey.getFlagHandler()), negative,
891 _ctrl);
892
893 if (_ctrl.doMeasurePsf) {
894 // Compute moments of Psf model. In the interest of implementing this quickly, we're just
895 // calling Psf::computeShape(), which delegates to SdssShapeResult::computeAdaptiveMoments
896 // for all nontrivial Psf classes. But this could in theory save the results of a shape
897 // computed some other way as part of base_SdssShape, which might be confusing. We should
898 // fix this eventually either by making Psf shape measurement not part of base_SdssShape, or
899 // by making the measurements stored with shape.sdss always computed via the
900 // SdssShapeAlgorithm instead of delegating to the Psf class.
901 try {
902 std::shared_ptr<afw::detection::Psf const> psf = exposure.getPsf();
903 if (!psf) {
904 result.flags[PSF_SHAPE_BAD.number] = true;
905 } else {
906 _resultKey.setPsfShape(measRecord, psf->computeShape(geom::Point2D(result.x, result.y)));
907 }
908 } catch (pex::exceptions::Exception &err) {
909 result.flags[PSF_SHAPE_BAD.number] = true;
910 }
911 }
912
913 measRecord.set(_resultKey, result);
914}
915
917 _resultKey.getFlagHandler().handleFailure(measRecord, error);
918}
919
920#define INSTANTIATE_IMAGE(IMAGE) \
921 template SdssShapeResult SdssShapeAlgorithm::computeAdaptiveMoments( \
922 IMAGE const &, geom::Point2D const &, bool, Control const &); \
923 template FluxResult SdssShapeAlgorithm::computeFixedMomentsFlux( \
924 IMAGE const &, afw::geom::ellipses::Quadrupole const &, geom::Point2D const &)
925
926#define INSTANTIATE_PIXEL(PIXEL) \
927 INSTANTIATE_IMAGE(afw::image::Image<PIXEL>); \
928 INSTANTIATE_IMAGE(afw::image::MaskedImage<PIXEL>);
929
931INSTANTIATE_PIXEL(float);
932INSTANTIATE_PIXEL(double);
933
936 : BaseTransform{name}, _instFluxTransform{name, mapper}, _centroidTransform{name, mapper} {
937 // If the input schema has a PSF flag -- it's optional -- assume we are also transforming the PSF.
938 _transformPsf = mapper.getInputSchema().getNames().count("sdssShape_flag_psf") ? true : false;
939
940 // Skip the last flag if not transforming the PSF shape.
943 if (flag == SdssShapeAlgorithm::FAILURE) continue;
944 if (mapper.getInputSchema().getNames().count(mapper.getInputSchema().join(name, flag.name)) == 0)
945 continue;
947 mapper.getInputSchema().find<afw::table::Flag>(name + "_" + flag.name).key;
948 mapper.addMapping(key);
949 }
950
951 _outShapeKey = ShapeResultKey::addFields(mapper.editOutputSchema(), name, "Shape in celestial moments",
953 if (_transformPsf) {
954 _outPsfShapeKey = afw::table::QuadrupoleKey::addFields(mapper.editOutputSchema(), name + "_psf",
955 "PSF shape in celestial moments",
957 }
958}
959
961 afw::table::BaseCatalog &outputCatalog, afw::geom::SkyWcs const &wcs,
962 afw::image::PhotoCalib const &photoCalib) const {
963 // The instFlux and cetroid transforms will check that the catalog lengths
964 // match and throw if not, so we don't repeat the test here.
965 _instFluxTransform(inputCatalog, outputCatalog, wcs, photoCalib);
966 _centroidTransform(inputCatalog, outputCatalog, wcs, photoCalib);
967
968 CentroidResultKey centroidKey(inputCatalog.getSchema()[_name]);
969 ShapeResultKey inShapeKey(inputCatalog.getSchema()[_name]);
970 afw::table::QuadrupoleKey inPsfShapeKey;
971 if (_transformPsf) {
972 inPsfShapeKey = afw::table::QuadrupoleKey(
973 inputCatalog.getSchema()[inputCatalog.getSchema().join(_name, "psf")]);
974 }
975
976 afw::table::SourceCatalog::const_iterator inSrc = inputCatalog.begin();
977 afw::table::BaseCatalog::iterator outSrc = outputCatalog.begin();
978 for (; inSrc != inputCatalog.end(); ++inSrc, ++outSrc) {
979 ShapeResult inShape = inShapeKey.get(*inSrc);
980 ShapeResult outShape;
981
982 // The transformation from the (x, y) to the (Ra, Dec) basis.
984 wcs.linearizePixelToSky(centroidKey.get(*inSrc).getCentroid(), geom::radians);
985 outShape.setShape(inShape.getShape().transform(crdTr.getLinear()));
986
987 // Transformation matrix from pixel to celestial basis.
989 outShape.setShapeErr((m * inShape.getShapeErr().cast<double>() * m.transpose()).cast<ErrElement>());
990
991 _outShapeKey.set(*outSrc, outShape);
992
993 if (_transformPsf) {
994 _outPsfShapeKey.set(*outSrc, inPsfShapeKey.get(*inSrc).transform(crdTr.getLinear()));
995 }
996 }
997}
998
999} // namespace base
1000} // namespace meas
1001} // namespace lsst
AmpInfoBoxKey bbox
Definition Amplifier.cc:117
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
afw::table::Key< afw::table::Array< ImagePixelT > > image
uint64_t * ptr
Definition RangeSet.cc:95
SchemaMapper * mapper
int y
Definition SpanSet.cc:48
int m
Definition SpanSet.cc:48
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Definition SkyWcs.h:117
Transformer transform(lsst::geom::LinearTransform const &transform)
Return the transform that maps the ellipse to the unit circle.
An ellipse core with quadrupole moments as parameters.
Definition Quadrupole.h:47
double getDeterminant() const
Return the determinant of the matrix representation.
Definition Quadrupole.h:83
A class to contain the data, WCS, and other information needed to describe an image of the sky.
Definition Exposure.h:72
A class to manipulate images, masks, and variance as a single object.
Definition MaskedImage.h:74
lsst::afw::image::Image< ImagePixelT > Image
Definition MaskedImage.h:86
The photometric calibration of an exposure.
Definition PhotoCalib.h:114
Tag types used to declare specialized field types.
Definition misc.h:31
Base class for all records.
Definition BaseRecord.h:31
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Definition BaseRecord.h:151
Schema getSchema() const
Return the Schema that holds this record's fields and keys.
Definition BaseRecord.h:80
void set(Key< T > const &key, U const &value)
Set value of a field for the given key.
Definition BaseRecord.h:164
iterator begin()
Iterator access.
Definition Catalog.h:400
bool isValid() const noexcept
Return true if the key was initialized to valid offset.
Definition Key.h:97
A FunctorKey used to get or set a geom::ellipses::Quadrupole from a tuple of constituent Keys.
Definition aggregates.h:369
geom::ellipses::Quadrupole get(BaseRecord const &record) const override
Get a Quadrupole from the given record.
void set(BaseRecord &record, geom::ellipses::Quadrupole const &value) const override
Set a Quadrupole in the given record.
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.
bool isValid() const noexcept
Return True if all the constituent Keys are valid.
Definition aggregates.h:429
Defines the fields and offsets for a table.
Definition Schema.h:51
SchemaItem< T > find(std::string const &name) const
Find a SchemaItem in the Schema by name.
Definition Schema.cc:467
A mapping between the keys of two Schemas, used to copy data between them.
typename Base::const_iterator const_iterator
Record class that contains measurements made on a single exposure.
Definition Source.h:78
A proxy type for name lookups in a Schema.
Definition Schema.h:367
An affine coordinate transformation consisting of a linear transformation and an offset.
LinearTransform const & getLinear() const noexcept
A floating-point coordinate rectangle geometry.
Definition Box.h:413
An integer coordinate rectangle.
Definition Box.h:55
Abstract base class for all C++ measurement transformations.
Definition Transform.h:86
A FunctorKey for CentroidResult.
virtual CentroidResult get(afw::table::BaseRecord const &record) const
Get a CentroidResult from the given record.
bool isValid() const
Return True if the centroid key is valid.
static CentroidResultKey addFields(afw::table::Schema &schema, std::string const &name, std::string const &doc, UncertaintyEnum uncertainty)
Add the appropriate fields to a Schema, and return a CentroidResultKey that manages them.
vector-type utility class to build a collection of FlagDefinitions
Definition FlagHandler.h:60
std::size_t size() const
return the current size (number of defined elements) of the collection
Utility class for handling flag fields that indicate the failure modes of an algorithm.
void handleFailure(afw::table::BaseRecord &record, MeasurementError const *error=nullptr) const
Handle an expected or unexpected Exception thrown by a measurement algorithm.
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.
void setValue(afw::table::BaseRecord &record, std::size_t i, bool value) const
Set the flag field corresponding to the given flag index.
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.
bool isValid() const
Return True if both the instFlux and instFluxErr Keys are valid.
static FluxResultKey addFields(afw::table::Schema &schema, std::string const &name, std::string const &doc)
Add a pair of _instFlux, _instFluxErr fields to a Schema, and return a FluxResultKey that points to t...
Exception to be thrown when a measurement algorithm experiences a known failure mode.
Definition exceptions.h:48
static FlagDefinition const SHIFT
Definition SdssShape.h:158
static FlagDefinitionList const & getFlagDefinitions()
Definition SdssShape.cc:58
virtual void measure(afw::table::SourceRecord &measRecord, afw::image::Exposure< float > const &exposure) const
Called to measure a single child source in an image.
Definition SdssShape.cc:881
static Result computeAdaptiveMoments(ImageT const &image, geom::Point2D const &position, bool negative=false, Control const &ctrl=Control())
Compute the adaptive Gaussian-weighted moments of an image.
static FlagDefinition const FAILURE
Definition SdssShape.h:155
static unsigned int const N_FLAGS
Definition SdssShape.h:154
static FluxResult computeFixedMomentsFlux(ImageT const &image, afw::geom::ellipses::Quadrupole const &shape, geom::Point2D const &position)
Compute the instFlux within a fixed Gaussian aperture.
Definition SdssShape.cc:834
static FlagDefinition const MAXITER
Definition SdssShape.h:159
static FlagDefinition const UNWEIGHTED
Definition SdssShape.h:157
static FlagDefinition const PSF_SHAPE_BAD
Definition SdssShape.h:160
static FlagDefinition const UNWEIGHTED_BAD
Definition SdssShape.h:156
SdssShapeAlgorithm(Control const &ctrl, std::string const &name, afw::table::Schema &schema)
Definition SdssShape.cc:759
virtual void fail(afw::table::SourceRecord &measRecord, MeasurementError *error=nullptr) const
Handle an exception thrown by the current algorithm by setting flags in the given record.
Definition SdssShape.cc:916
A C++ control class to handle SdssShapeAlgorithm's configuration.
Definition SdssShape.h:52
double maxShift
"Maximum centroid shift, limited to 2-10" ;
Definition SdssShape.h:56
int maxIter
"Maximum number of iterations" ;
Definition SdssShape.h:55
float tol2
"Convergence tolerance for FWHM" ;
Definition SdssShape.h:58
bool doMeasurePsf
"Whether to also compute the shape of the PSF model" ;
Definition SdssShape.h:59
float tol1
"Convergence tolerance for e1,e2" ;
Definition SdssShape.h:57
double background
"Additional value to add to background" ;
Definition SdssShape.h:54
Result object SdssShapeAlgorithm.
Definition SdssShape.h:224
SdssShapeResult()
Constructor; initializes everything to NaN.
Definition SdssShape.cc:636
std::bitset< SdssShapeAlgorithm::N_FLAGS > flags
Status flags (see SdssShapeAlgorithm).
Definition SdssShape.h:230
A FunctorKey that maps SdssShapeResult to afw::table Records.
Definition SdssShape.h:73
SdssShapeResultKey()
Default constructor; instance will not be usuable unless subsequently assigned to.
Definition SdssShape.h:90
virtual SdssShapeResult get(afw::table::BaseRecord const &record) const
Get an SdssShapeResult from the given record.
Definition SdssShape.cc:707
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:641
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:722
bool isValid() const
Return True if the key is valid.
Definition SdssShape.cc:752
bool operator==(SdssShapeResultKey const &other) const
Compare the FunctorKey for equality with another, using the underlying Keys.
Definition SdssShape.cc:744
virtual void set(afw::table::BaseRecord &record, SdssShapeResult const &value) const
Set an SdssShapeResult in the given record.
Definition SdssShape.cc:726
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:739
FlagHandler const & getFlagHandler() const
Definition SdssShape.h:125
virtual void operator()(afw::table::SourceCatalog const &inputCatalog, afw::table::BaseCatalog &outputCatalog, afw::geom::SkyWcs const &wcs, afw::image::PhotoCalib const &photoCalib) const
Definition SdssShape.cc:960
SdssShapeTransform(Control const &ctrl, std::string const &name, afw::table::SchemaMapper &mapper)
Definition SdssShape.cc:934
A FunctorKey for ShapeResult.
static ShapeResultKey addFields(afw::table::Schema &schema, std::string const &name, std::string const &doc, UncertaintyEnum uncertainty, afw::table::CoordinateType coordType=afw::table::CoordinateType::PIXEL)
Add the appropriate fields to a Schema, and return a ShapeResultKey that manages them.
virtual ShapeResult get(afw::table::BaseRecord const &record) const
Get a ShapeResult from the given record.
bool isValid() const
Return True if the shape key is valid.
virtual void set(afw::table::BaseRecord &record, ShapeResult const &value) const
Set a ShapeResult in the given record.
Reports arguments outside the domain of an operation.
Definition Runtime.h:57
Provides consistent interface for LSST exceptions.
Definition Exception.h:107
Reports invalid arguments.
Definition Runtime.h:66
Reports errors in the logical structure of the program.
Definition Runtime.h:46
Reports attempts to access elements using an invalid key.
Definition Runtime.h:151
Reports errors that are due to events beyond the control of the program.
Definition Runtime.h:104
T exp(T... args)
T fabs(T... args)
T isnan(T... args)
T make_pair(T... args)
T make_tuple(T... args)
T max(T... args)
T min(T... args)
int positionToIndex(double pos)
Convert image position to nearest integer index.
Definition ImageUtils.h:69
Extent< double, 2 > Extent2D
Definition Extent.h:400
double constexpr TWOPI
Definition Angle.h:41
AngleUnit constexpr radians
constant with units of radians
Definition Angle.h:109
double constexpr PI
The ratio of a circle's circumference to diameter.
Definition Angle.h:40
@ SIGMA_ONLY
Only the diagonal elements of the covariance matrix are provided.
Definition constants.h:45
@ FULL_COVARIANCE
The full covariance matrix is provided.
Definition constants.h:46
@ NO_UNCERTAINTY
Algorithm provides no uncertainy information at all.
Definition constants.h:44
Eigen::Matrix< ShapeElement, 3, 3, Eigen::DontAlign > ShapeTrMatrix
Definition constants.h:62
ShapeTrMatrix makeShapeTransformMatrix(geom::LinearTransform const &xform)
Construct a matrix suitable for transforming second moments.
STL namespace.
T pow(T... args)
T quiet_NaN(T... args)
T sqrt(T... args)
afw::table::Key< double > weight
#define INSTANTIATE_PIXEL(PIXEL)
Definition SdssShape.cc:926
A reusable struct for centroid measurements.
Centroid const getCentroid() const
Return a Point object containing the measured x and y.
Simple class used to define and document flags The name and doc constitute the identity of the FlagDe...
Definition FlagHandler.h:40
A reusable result struct for instFlux measurements.
meas::base::Flux instFlux
Measured instFlux in DN.
A reusable struct for moments-based shape measurements.
Shape const getShape() const
Return an afw::geom::ellipses object corresponding to xx, yy, xy.
ShapeCov const getShapeErr() const
Return the 3x3 symmetric covariance matrix, with rows and columns ordered (xx, yy,...
void setShape(Shape const &shape)
Set struct elements from the given Quadrupole object.
ShapeElement xy
image or model second moment for xy^2
ShapeElement xx
image or model second moment for x^2
void setShapeErr(ShapeCov const &matrix)
Set the struct standard deviation elements from the given matrix, with rows and columns ordered (xx,...
ShapeElement yy
image or model second moment for y^2