LSST Applications  21.0.0+c4f5df5339,21.0.0+e70536a077,21.0.0-1-ga51b5d4+7c60f8a6ea,21.0.0-10-g560fb7b+411cd868f8,21.0.0-10-gcf60f90+8c49d86aa0,21.0.0-13-gc485e61d+38156233bf,21.0.0-16-g7a993c7b9+1041c3824f,21.0.0-2-g103fe59+d9ceee3e5a,21.0.0-2-g1367e85+0b2f7db15a,21.0.0-2-g45278ab+e70536a077,21.0.0-2-g5242d73+0b2f7db15a,21.0.0-2-g7f82c8f+feb9862f5e,21.0.0-2-g8f08a60+9c9a9cfcc8,21.0.0-2-ga326454+feb9862f5e,21.0.0-2-gde069b7+bedfc5e1fb,21.0.0-2-gecfae73+417509110f,21.0.0-2-gfc62afb+0b2f7db15a,21.0.0-3-g21c7a62+a91f7c0b59,21.0.0-3-g357aad2+062581ff1a,21.0.0-3-g4be5c26+0b2f7db15a,21.0.0-3-g65f322c+85aa0ead76,21.0.0-3-g7d9da8d+c4f5df5339,21.0.0-3-gaa929c8+411cd868f8,21.0.0-3-gc44e71e+fd4029fd48,21.0.0-3-ge02ed75+5d9b90b8aa,21.0.0-38-g070523fc+44fda2b515,21.0.0-4-g591bb35+5d9b90b8aa,21.0.0-4-g88306b8+3cdc83ea97,21.0.0-4-gc004bbf+d52368b591,21.0.0-4-gccdca77+a5c54364a0,21.0.0-5-g7ebb681+81e2098694,21.0.0-5-gdf36809+87b8d260e6,21.0.0-6-g2d4f3f3+e70536a077,21.0.0-6-g4e60332+5d9b90b8aa,21.0.0-6-g5ef7dad+3f4e29eeae,21.0.0-7-gc8ca178+0f5e56d48f,21.0.0-9-g9eb8d17+cc2c7a81aa,master-gac4afde19b+5d9b90b8aa,w.2021.07
LSST Data Management Base Package
PixelRegion.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
2 /*
3  * LSST Data Management System
4  * Copyright 2008, 2009, 2010 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 <cmath>
25 #include <limits>
26 #include "boost/optional.hpp"
29 
30 namespace lsst {
31 namespace afw {
32 namespace geom {
33 namespace ellipses {
34 
35 //
36 // Functor for computing the x boundary points of an ellipse given y.
37 //
38 // The equation for the boundary of the ellipse is:
39 //
40 // F_{xx} x^2 + 2 F_{xy} x y + F_{yy} y^2 = 1
41 //
42 // where F is the matrix inverse of the quadrupole matrix Q and x and y are
43 // relative to the center.
44 //
45 // Given some fixed y, we want to solve this pretty vanilla quadratic equation
46 // for x. Algebra is not shown, but the solution is:
47 //
48 // x = p y ± \sqrt{r - t y^2}
49 //
50 // where p = -F_{xy}/F_{xx} = Q_{xy}/Q_{yy}
51 // r = 1/F_{xx} = det(Q)/Q_{yy} = Q_{xx} - Q_{xy}p
52 // t = det(F)/F_{xx}^2 = det(Q)/Q_{yy}^2 = r/Q_{yy}
53 //
54 // This runs into divide-by-zero problems as Q_{yy} approaches zero. That
55 // corresponds to the ellipse approaching a horizontal line segment. The
56 // equation for that is just
57 //
58 // x = ± \sqrt{r}
59 //
60 // with r = Q_{xx}. That's equivalent to just setting p and t to zero.
61 //
63 public:
64 
65  explicit EllipseHorizontalLineIntersection(Ellipse const& ellipse) :
66  _center(ellipse.getCenter()),
67  _p(0.0), _t(0.0), _r(0.0)
68  {
69  Quadrupole converted(ellipse.getCore());
70  converted.normalize();
71  auto q = converted.getMatrix();
72  if (q(1, 1) < q(0, 0)*std::numeric_limits<double>::epsilon()) {
73  _r = q(0, 0);
74  } else {
75  _p = q(0, 1)/q(1, 1);
76  _r = q(0, 0) - q(0, 1)*_p;
77  _t = _r/q(1, 1);
78  }
79  }
80 
81  boost::optional<std::pair<double, double>> xAt(double y) const {
82  double yc = y - _center.getY();
83  double xc = _p*yc + _center.getX();
84  double s = _r - _t*yc*yc;
85  if (s < 0) {
86  return boost::none;
87  } else {
88  double d = std::sqrt(s);
89  return std::make_pair(xc - d, xc + d);
90  }
91  }
92 
93 private:
94  lsst::geom::Point2D _center;
95  double _p;
96  double _t;
97  double _r;
98 };
99 
101  : _bbox(ellipse.computeBBox(), lsst::geom::Box2I::EXPAND)
102 {
103  // Initial temporary bounding box that may be larger than the final one.
105 
106  if (envelope.isEmpty()) {
107  // If the outer bbox is empty, we know there can't be any spans that
108  // contain ellipse pixels.
109  return;
110  }
111 
112  // Helper class that does the hard work: compute the boundary points of the
113  // ellipse in x that intersect a horizontal line at some given y.
114  EllipseHorizontalLineIntersection intersection(ellipse);
115 
116  // Iterate over pixel rows in the bounding box, computing the intersection
117  // of the ellipse with that y coordinate.
118  int const yEnd = envelope.getEndY();
119  for (int y = envelope.getBeginY(); y != yEnd; ++y) {
120  auto x = intersection.xAt(y);
121  if (x) {
122  int xMin = std::ceil(x->first);
123  int xMax = std::floor(x->second);
124  if (xMax < xMin) continue;
125  _spans.emplace_back(y, xMin, xMax);
126  _bbox.include(lsst::geom::Point2I(xMin, y));
127  _bbox.include(lsst::geom::Point2I(xMax, y));
128  }
129  }
130 }
131 
132 Span const PixelRegion::getSpanAt(int y) const {
133  if (_bbox.isEmpty()) {
134  throw LSST_EXCEPT(
136  "PixelRegion is empty."
137  );
138  }
139  if (y < _bbox.getMinY() || y > _bbox.getMaxY()) {
140  throw LSST_EXCEPT(
142  (boost::format("No span at y=%s in pixel region with rows between y=%s and y=%s")
143  % y % _bbox.getMinY() % _bbox.getMaxY()).str()
144  );
145  }
146  return _spans[y - _bbox.getBeginY()];
147 }
148 
149 } // namespace ellipses
150 } // namespace geom
151 } // namespace afw
152 } // namespace lsst
double x
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
T ceil(T... args)
A range of pixels within one row of an Image.
Definition: Span.h:47
boost::optional< std::pair< double, double > > xAt(double y) const
Definition: PixelRegion.cc:81
An ellipse defined by an arbitrary BaseCore and a center point.
Definition: Ellipse.h:51
BaseCore const & getCore() const
Return the ellipse core.
Definition: Ellipse.h:71
lsst::geom::Box2D computeBBox() const
Return the bounding box of the ellipse.
Definition: Ellipse.cc:55
PixelRegion(Ellipse const &ellipse)
Construct a PixelRegion from an Ellipse.
Definition: PixelRegion.cc:100
Span const getSpanAt(int y) const
Return the span at the given y coordinate value.
Definition: PixelRegion.cc:132
An ellipse core with quadrupole moments as parameters.
Definition: Quadrupole.h:47
An integer coordinate rectangle.
Definition: Box.h:55
int getMinY() const noexcept
Definition: Box.h:158
void include(Point2I const &point)
Expand this to ensure that this->contains(point).
Definition: Box.cc:152
bool isEmpty() const noexcept
Return true if the box contains no points.
Definition: Box.h:213
int getEndY() const noexcept
Definition: Box.h:177
int getBeginY() const noexcept
Definition: Box.h:173
int getMaxY() const noexcept
Definition: Box.h:162
Reports attempts to access elements outside a valid range of indices.
Definition: Runtime.h:89
T floor(T... args)
T make_pair(T... args)
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
T sqrt(T... args)
int y
Definition: SpanSet.cc:49