LSSTApplications
20.0.0
LSSTDataManagementBasePackage
stack
1a1d771
Linux64
afw
20.0.0
src
geom
ellipses
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"
27
#include "
lsst/afw/geom/ellipses/PixelRegion.h
"
28
#include "
lsst/afw/geom/ellipses/Quadrupole.h
"
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
//
62
class
EllipseHorizontalLineIntersection
{
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
100
PixelRegion::PixelRegion
(
Ellipse
const
& ellipse)
101
: _bbox(ellipse.computeBBox(),
lsst
::
geom
::Box2I::EXPAND)
102
{
103
// Initial temporary bounding box that may be larger than the final one.
104
lsst::geom::Box2I
envelope(ellipse.
computeBBox
(),
lsst::geom::Box2I::EXPAND
);
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
(
135
pex::exceptions::OutOfRangeError
,
136
"PixelRegion is empty."
137
);
138
}
139
if
(
y
< _bbox.
getMinY
() ||
y
> _bbox.
getMaxY
()) {
140
throw
LSST_EXCEPT
(
141
pex::exceptions::OutOfRangeError
,
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
y
int y
Definition:
SpanSet.cc:49
std::floor
T floor(T... args)
lsst::afw::geom::Span
A range of pixels within one row of an Image.
Definition:
Span.h:47
ellipses
lsst::afw::geom::ellipses::Ellipse::computeBBox
lsst::geom::Box2D computeBBox() const
Return the bounding box of the ellipse.
Definition:
Ellipse.cc:55
pex.config.history.format
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition:
history.py:174
lsst::afw
Definition:
imageAlgorithm.dox:1
Quadrupole.h
std::sqrt
T sqrt(T... args)
lsst::geom::Box2I::isEmpty
bool isEmpty() const noexcept
Return true if the box contains no points.
Definition:
Box.h:213
lsst::geom::Box2I::getEndY
int getEndY() const noexcept
Definition:
Box.h:177
lsst::geom::Box2I::EXPAND
@ EXPAND
Definition:
Box.h:63
lsst::afw::geom::ellipses::EllipseHorizontalLineIntersection::xAt
boost::optional< std::pair< double, double > > xAt(double y) const
Definition:
PixelRegion.cc:81
lsst::afw::geom::ellipses::PixelRegion::PixelRegion
PixelRegion(Ellipse const &ellipse)
Construct a PixelRegion from an Ellipse.
Definition:
PixelRegion.cc:100
lsst::geom::Box2I::getBeginY
int getBeginY() const noexcept
Definition:
Box.h:173
lsst::afw::geom::ellipses::EllipseHorizontalLineIntersection::EllipseHorizontalLineIntersection
EllipseHorizontalLineIntersection(Ellipse const &ellipse)
Definition:
PixelRegion.cc:65
x
double x
Definition:
ChebyshevBoundedField.cc:277
lsst::geom::Box2I::include
void include(Point2I const &point)
Expand this to ensure that this->contains(point).
Definition:
Box.cc:152
lsst::afw::geom::ellipses::EllipseHorizontalLineIntersection
Definition:
PixelRegion.cc:62
std::ceil
T ceil(T... args)
lsst
A base class for image defects.
Definition:
imageAlgorithm.dox:1
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition:
Exception.h:48
lsst::afw::geom::ellipses::PixelRegion::getSpanAt
Span const getSpanAt(int y) const
Return the span at the given y coordinate value.
Definition:
PixelRegion.cc:132
lsst::afw::geom::ellipses::Ellipse
An ellipse defined by an arbitrary BaseCore and a center point.
Definition:
Ellipse.h:51
PixelRegion.h
lsst::afw::geom::ellipses::Quadrupole
An ellipse core with quadrupole moments as parameters.
Definition:
Quadrupole.h:47
lsst::geom
Definition:
geomOperators.dox:4
lsst::geom::Box2I::getMaxY
int getMaxY() const noexcept
Definition:
Box.h:162
lsst.pipe.tasks.mergeDetections.converted
converted
Definition:
mergeDetections.py:377
lsst::geom::Point< double, 2 >
lsst::pex::exceptions::OutOfRangeError
Reports attempts to access elements outside a valid range of indices.
Definition:
Runtime.h:89
lsst::geom::Box2I
An integer coordinate rectangle.
Definition:
Box.h:55
lsst::afw::geom::ellipses::Ellipse::getCore
BaseCore const & getCore() const
Return the ellipse core.
Definition:
Ellipse.h:71
std::make_pair
T make_pair(T... args)
std::numeric_limits
lsst::geom::Box2I::getMinY
int getMinY() const noexcept
Definition:
Box.h:158
Generated on Wed Jun 24 2020 18:10:02 for LSSTApplications by
1.8.18