LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
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
int y
Definition: SpanSet.cc: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)