Loading [MathJax]/extensions/tex2jax.js
LSST Applications 28.0.0,g1653933729+a8ce1bb630,g1a997c3884+a8ce1bb630,g28da252d5a+5bd70b7e6d,g2bbee38e9b+638fca75ac,g2bc492864f+638fca75ac,g3156d2b45e+07302053f8,g347aa1857d+638fca75ac,g35bb328faa+a8ce1bb630,g3a166c0a6a+638fca75ac,g3e281a1b8c+7bbb0b2507,g4005a62e65+17cd334064,g414038480c+5b5cd4fff3,g41af890bb2+4ffae9de63,g4e1a3235cc+0f1912dca3,g6249c6f860+3c3976f90c,g80478fca09+46aba80bd6,g82479be7b0+77990446f6,g858d7b2824+78ba4d1ce1,g89c8672015+f667a5183b,g9125e01d80+a8ce1bb630,ga5288a1d22+2a6264e9ca,gae0086650b+a8ce1bb630,gb58c049af0+d64f4d3760,gc22bb204ba+78ba4d1ce1,gc28159a63d+638fca75ac,gcf0d15dbbd+32ddb6096f,gd6b7c0dfd1+3e339405e9,gda3e153d99+78ba4d1ce1,gda6a2b7d83+32ddb6096f,gdaeeff99f8+1711a396fd,gdd5a9049c5+b18c39e5e3,ge2409df99d+a5e4577cdc,ge33fd446bb+78ba4d1ce1,ge79ae78c31+638fca75ac,gf0baf85859+64e8883e75,gf5289d68f6+e1b046a8d7,gfa443fc69c+91d9ed1ecf,gfda6b12a05+8419469a56
LSST Data Management Base Package
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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
30namespace lsst {
31namespace afw {
32namespace geom {
33namespace 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//
63public:
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
93private:
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
132Span 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
#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
An ellipse defined by an arbitrary BaseCore and a center point.
Definition Ellipse.h:51
lsst::geom::Box2D computeBBox() const
Return the bounding box of the ellipse.
Definition Ellipse.cc:55
BaseCore const & getCore() const
Return the ellipse core.
Definition Ellipse.h:71
PixelRegion(Ellipse const &ellipse)
Construct a PixelRegion from an Ellipse.
Span const getSpanAt(int y) const
Return the span at the given y coordinate value.
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)
T sqrt(T... args)