gtfs-actions/generate-shapes
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

get rid of boost

Browse source
This commit is contained in:
Patrick Brosi 2018-07-22 17:17:53 +02:00 committed by Patrick Brosi
parent 63dfda9749
commit 13fb21f626
12 changed files with 717 additions and 174 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

34
README.md
View file

@ -1,15 +1,4 @@
```
_-====-__-=-__-===-__-=======-__
_( _)
OO( )
. o '===-______-===-____-==-__-====='
.o
. ______ _______________
_()_||__|| __o^o___ | [] [] [] [] |
( | | | |o
/-OO----OO""="OO--OO"="OO---------OO"
############################################################
```
![Map-Matched path of a single train through Switzerland](geo/schweiz_mmatched.png?raw=true)
# pfaedle
@ -18,7 +7,7 @@ Precise map-matching for public transit schedules (GTFS data).
## Requirements
* `cmake`
* `gcc` >= 4.8
* `gcc` >= 4.8 (may work on lower versions, untested)
## Building and Installation
@ -49,16 +38,29 @@ pfaedle -c <CFG FILE> -x <OSM FILE> <GTFS INPUT FOLDER>
A shape'd version of the input GTFS feed will be written to `./gtfs-out`.
A default configuration file `pfaedle.cfg` can be found in this repo.
By default, shapes are only calculated for trips that don't have a shape in the
input feed. To drop all existing shapes, use the `-D` flag.
For example, you may generate (and replace existing, see -D parameter) shapes for the GTFS dataset for Freiburg like this:
```
$ mkdir freiburg_gtfs && cd freiburg_gtfs
$ wget https://fritz.freiburg.de/csv_Downloads/VAGFR.zip
$ unzip VAGFR.zip
$ wget http://download.geofabrik.de/europe/germany/baden-wuerttemberg/freiburg-regbez-latest.osm.bz2
$ bunzip2 freiburg-regbez-latest.osm.bz2
$ mkdir gtfs-out
$ pfaedle -D -c pfaedle.cfg -x freiburg-regbez-latest.osm .
```
A default configuration file `pfaedle.cfg` can be found in this repo.
## Generating shapes for a specific MOT
To generate shapes only for a specific mot, use the `-m` option. Possible
values are either `tram`, `bus`, `rail`, `subway`, `ferry`, `funicular`,
`gondola`, `all`.
`gondola`, `all` (default).
Multiple values can be specified (comma separated).

BIN
geo/schweiz_mmatched.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 1 MiB

2
src/pfaedle/PfaedleMain.cpp
View file

@ -154,7 +154,7 @@ int main(int argc, char** argv) {
if (cfg.evaluate) ecoll.printStats(&std::cout);
if (cfg.feedPaths.size()) {
LOG(INFO) << "Writing output GTFS...";
LOG(INFO) << "Writing output GTFS to " << cfg.outputPath << " ...";
ad::cppgtfs::Writer w;
w.write(&gtfs, cfg.outputPath);
}

1
src/pfaedle/config/PfaedleConfig.h
View file

@ -21,6 +21,7 @@ struct Config {
: dbgOutputPath("geo"),
solveMethod("global"),
evalPath("."),
outputPath("gtfs-out"),
dropShapes(false),
useHMM(false),
writeGraph(false),

37
src/util/geo/Box.h
View file

@ -41,10 +41,41 @@ class Box {
template <typename T>
class RotatedBox {
public:
RotatedBox() : _box() {}
RotatedBox(const Box<T>& box) : _box(box), _deg(0), _center() {}
RotatedBox() : _box(), _deg(0), _center() {}
RotatedBox(const Box<T>& box)
: _box(box),
_deg(0),
_center(Point<T>(
(box.getUpperRight().getX() - box.getLowerLeft().getX()) / T(2),
(box.getUpperRight().getY() - box.getLowerLeft().getY()) / T(2))) {}
RotatedBox(const Point<T>& ll, const Point<T>& ur)
: _box(ll, ur), _deg(0), _center() {}
: _box(ll, ur),
_deg(0),
_center(Point<T>((ur.getX() - ll.getX()) / T(2),
(ur.getY() - ll.getY()) / T(2))) {}
RotatedBox(const Box<T>& box, double deg)
: _box(box),
_deg(deg),
_center(Point<T>(
(box.getUpperRight().getX() - box.getLowerLeft().getX()) / T(2),
(box.getUpperRight().getY() - box.getLowerLeft().getY()) / T(2))) {}
RotatedBox(const Point<T>& ll, const Point<T>& ur, double deg)
: _box(ll, ur),
_deg(deg),
_center(Point<T>((ur.getX() - ll.getX()) / T(2),
(ur.getY() - ll.getY()) / T(2))) {}
RotatedBox(const Box<T>& box, double deg, const Point<T>& center)
: _box(box), _deg(deg), _center(center) {}
RotatedBox(const Point<T>& ll, const Point<T>& ur, double deg,
const Point<T>& center)
: _box(ll, ur), _deg(deg), _center(center) {}
const Box<T>& getBox() const { return _box; }
Box<T>& getBox() { return _box; }
double getDegree() const { return _deg; }
const Point<T>& getCenter() const { return _center; }
Point<T>& getCenter() { return _center; }
private:
Box<T> _box;

517
src/util/geo/Geo.h
View file

@ -8,6 +8,8 @@
#define _USE_MATH_DEFINES
#include <math.h>
#include <algorithm>
#include <cassert>
#include <iostream>
#include <sstream>
#include "util/Misc.h"
@ -45,7 +47,7 @@ typedef Polygon<double> DPolygon;
typedef Polygon<float> FPolygon;
typedef Polygon<int> IPolygon;
const static double EPSILON = 0.0000001;
const static double EPSILON = 0.00000000001;
const static double RAD = 0.017453292519943295; // PI/180
// _____________________________________________________________________________
@ -93,6 +95,14 @@ inline Point<T> centroid(const Box<T> box) {
return centroid(LineSegment<T>(box.getLowerLeft(), box.getUpperRight()));
}
// _____________________________________________________________________________
template <typename T, template <typename> typename Geometry>
inline Point<T> centroid(std::vector<Geometry<T>> multigeo) {
Line<T> a;
for (const auto& g : multigeo) a.push_back(centroid(g));
return centroid(a);
}
// _____________________________________________________________________________
template <typename T>
inline Point<T> rotate(const Point<T>& p, double deg) {
@ -142,12 +152,63 @@ inline Polygon<T> rotate(Polygon<T> geo, double deg, const Point<T>& c) {
return geo;
}
// _____________________________________________________________________________
template <typename T, template <typename> typename Geometry>
inline std::vector<Geometry<T>> rotate(std::vector<Geometry<T>> multigeo,
double deg, const Point<T>& c) {
for (size_t i = 0; i < multigeo.size(); i++)
multigeo[i] = rotate(multigeo[i], deg, c);
return multigeo;
}
// _____________________________________________________________________________
template <typename T, template <typename> typename Geometry>
inline std::vector<Geometry<T>> rotate(std::vector<Geometry<T>> multigeo,
double deg) {
auto c = centroid(multigeo);
for (size_t i = 0; i < multigeo.size(); i++)
multigeo[i] = rotate(multigeo[i], deg, c);
return multigeo;
}
// _____________________________________________________________________________
template <typename T>
inline Point<T> move(const Point<T>& geo, T x, T y) {
return Point<T>(geo.getX() + x, geo.getY() + y);
}
// _____________________________________________________________________________
template <typename T>
inline Line<T> move(Line<T> geo, T x, T y) {
for (size_t i = 0; i < geo.size(); i++) geo[i] = move(geo[i], x, y);
return geo;
}
// _____________________________________________________________________________
template <typename T>
inline LineSegment<T> move(LineSegment<T> geo, T x, T y) {
geo.first = move(geo.first, x, y);
geo.second = move(geo.second, x, y);
return geo;
}
// _____________________________________________________________________________
template <typename T>
inline Polygon<T> move(Polygon<T> geo, T x, T y) {
for (size_t i = 0; i < geo.getOuter().size(); i++)
geo.getOuter()[i] = move(geo.getOuter()[i], x, y);
return geo;
}
// _____________________________________________________________________________
template <typename T, template <typename> typename Geometry>
inline std::vector<Geometry<T>> move(std::vector<Geometry<T>> multigeo, T x,
T y) {
for (size_t i = 0; i < multigeo.size(); i++)
multigeo[i] = move(multigeo[i], x, y);
return multigeo;
}
// _____________________________________________________________________________
template <typename T>
inline Box<T> minbox() {
@ -155,23 +216,20 @@ inline Box<T> minbox() {
}
// _____________________________________________________________________________
// template <typename T>
// inline RotatedBox<T> shrink(const RotatedBox<T>& b, double d) {
// double xd =
// b.b.getUpperRight().getX() - b.b.getLowerLeft().getX();
// double yd =
// b.b.getUpperRight().getY() - b.b.getLowerLeft().getY();
template <typename T>
inline RotatedBox<T> shrink(const RotatedBox<T>& b, double d) {
double xd = b.b.getUpperRight().getX() - b.b.getLowerLeft().getX();
double yd = b.b.getUpperRight().getY() - b.b.getLowerLeft().getY();
// if (xd <= 2 * d) d = xd / 2 - 1;
// if (yd <= 2 * d) d = yd / 2 - 1;
if (xd <= 2 * d) d = xd / 2 - 1;
if (yd <= 2 * d) d = yd / 2 - 1;
// Box<T> r(Point<T>(b.b.getLowerLeft().getX() + d,
// b.b.getLowerLeft().getY() + d),
// Point<T>(b.b.getUpperRight().getX() - d,
// b.b.getUpperRight().getY() - d));
Box<T> r(
Point<T>(b.b.getLowerLeft().getX() + d, b.b.getLowerLeft().getY() + d),
Point<T>(b.b.getUpperRight().getX() - d, b.b.getUpperRight().getY() - d));
// return RotatedBox<T>(r, b.rotateDeg, b.center);
// }
return RotatedBox<T>(r, b.getDegree(), b.getCenter());
}
// _____________________________________________________________________________
inline bool doubleEq(double a, double b) { return fabs(a - b) < 0.000001; }
@ -199,6 +257,12 @@ inline bool contains(const LineSegment<T>& l, const Box<T>& box) {
return contains(l.first, box) && contains(l.second, box);
}
// _____________________________________________________________________________
template <typename T>
inline bool contains(const Box<T>& b, const Box<T>& box) {
return contains(b.getLowerLeft(), box) && contains(b.getUpperRight(), box);
}
// _____________________________________________________________________________
template <typename T>
inline bool contains(const Point<T>& p, const LineSegment<T>& ls) {
@ -236,13 +300,17 @@ inline int8_t polyContCheck(const Point<T>& a, Point<T> b, Point<T> c) {
if (a.getY() == b.getY() && a.getY() == c.getY())
return (!((b.getX() <= a.getX() && a.getX() <= c.getX()) ||
(c.getX() <= a.getX() && a.getX() <= b.getX())));
if (a.getY() == b.getY() && a.getX() == b.getX()) return 0;
if (fabs(a.getY() - b.getY()) < EPSILON &&
fabs(a.getX() - b.getX()) < EPSILON)
return 0;
if (b.getY() > c.getY()) {
Point<T> tmp = b;
b = c;
c = tmp;
}
if (a.getY() <= b.getY() || a.getY() > c.getY()) return 1;
if (a.getY() <= b.getY() || a.getY() > c.getY()) {
return 1;
}
double d = (b.getX() - a.getX()) * (c.getY() - a.getY()) -
(b.getY() - a.getY()) * (c.getX() - a.getX());
@ -255,7 +323,9 @@ inline int8_t polyContCheck(const Point<T>& a, Point<T> b, Point<T> c) {
template <typename T>
inline bool contains(const Polygon<T>& polyC, const Polygon<T>& poly) {
for (const auto& p : polyC.getOuter()) {
if (!contains(p, poly)) return false;
if (!contains(p, poly)) {
return false;
}
}
return true;
}
@ -264,7 +334,20 @@ inline bool contains(const Polygon<T>& polyC, const Polygon<T>& poly) {
template <typename T>
inline bool contains(const Line<T>& l, const Polygon<T>& poly) {
for (const auto& p : l) {
if (!contains(p, poly)) return false;
if (!contains(p, poly)) {
return false;
}
}
return true;
}
// _____________________________________________________________________________
template <typename T>
inline bool contains(const Line<T>& l, const Line<T>& other) {
for (const auto& p : l) {
if (!contains(p, other)) {
return false;
}
}
return true;
}
@ -298,6 +381,16 @@ inline bool contains(const Polygon<T>& poly, const Line<T>& l) {
return true;
}
// _____________________________________________________________________________
template <typename T, template <typename> typename GeometryA,
template <typename> typename GeometryB>
inline bool contains(const std::vector<GeometryA<T>>& multigeo,
const GeometryB<T>& geo) {
for (const auto& g : multigeo)
if (!contains(g, geo)) return false;
return true;
}
// _____________________________________________________________________________
template <typename T>
inline bool intersects(const LineSegment<T>& ls1, const LineSegment<T>& ls2) {
@ -463,6 +556,13 @@ inline Point<T> intersection(const Point<T>& p1, const Point<T>& q1,
p2.getY(), q2.getX(), q2.getY());
}
// _____________________________________________________________________________
template <typename T>
inline Point<T> intersection(const LineSegment<T>& s1,
const LineSegment<T>& s2) {
return intersection(s1.first, s1.second, s2.first, s2.second);
}
// _____________________________________________________________________________
template <typename T>
inline bool lineIntersects(T p1x, T p1y, T q1x, T q1y, T p2x, T p2y, T q2x,
@ -534,11 +634,10 @@ inline double crossProd(const Point<T>& a, const Point<T>& b) {
// _____________________________________________________________________________
template <typename T>
inline double crossProd(const Point<T>& p, const LineSegment<T>& ls) {
LineSegment<T> lss(Point<T>(0, 0),
Point<T>(ls.second.getX() - ls.first.getX(),
ls.second.getY() - ls.first.getY()));
return crossProd(lss.second, Point<T>(p.getX() - ls.first.getX(),
p.getY() - ls.first.getY()));
return crossProd(
Point<T>(ls.second.getX() - ls.first.getX(),
ls.second.getY() - ls.first.getY()),
Point<T>(p.getX() - ls.first.getX(), p.getY() - ls.first.getY()));
}
// _____________________________________________________________________________
@ -555,6 +654,19 @@ inline std::string getWKT(const Point<T>& p) {
return ss.str();
}
// _____________________________________________________________________________
template <typename T>
inline std::string getWKT(const std::vector<Point<T>>& p) {
std::stringstream ss;
ss << "MULTIPOINT (";
for (size_t i = 0; i < p.size(); i++) {
if (i) ss << ", ";
ss << "(" << p[i].getX() << " " << p[i].getY() << ")";
}
ss << ")";
return ss.str();
}
// _____________________________________________________________________________
template <typename T>
inline std::string getWKT(const Line<T>& l) {
@ -568,6 +680,26 @@ inline std::string getWKT(const Line<T>& l) {
return ss.str();
}
// _____________________________________________________________________________
template <typename T>
inline std::string getWKT(const std::vector<Line<T>>& ls) {
std::stringstream ss;
ss << "MULTILINESTRING (";
for (size_t j = 0; j < ls.size(); j++) {
if (j) ss << ", ";
ss << "(";
for (size_t i = 0; i < ls[j].size(); i++) {
if (i) ss << ", ";
ss << ls[j][i].getX() << " " << ls[j][i].getY();
}
ss << ")";
}
ss << ")";
return ss.str();
}
// _____________________________________________________________________________
template <typename T>
inline std::string getWKT(const LineSegment<T>& l) {
@ -594,13 +726,35 @@ inline std::string getWKT(const Polygon<T>& p) {
std::stringstream ss;
ss << "POLYGON ((";
for (size_t i = 0; i < p.getOuter().size(); i++) {
if (i) ss << ", ";
ss << p.getOuter()[i].getX() << " " << p.getOuter()[i].getY();
ss << p.getOuter()[i].getX() << " " << p.getOuter()[i].getY() << ", ";
}
ss << p.getOuter().front().getX() << " " << p.getOuter().front().getY();
ss << "))";
return ss.str();
}
// _____________________________________________________________________________
template <typename T>
inline std::string getWKT(const std::vector<Polygon<T>>& ls) {
std::stringstream ss;
ss << "MULTIPOLYGON (";
for (size_t j = 0; j < ls.size(); j++) {
if (j) ss << ", ";
ss << "((";
for (size_t i = 0; i < ls[j].getOuter().size(); i++) {
ss << ls[j].getOuter()[i].getX() << " " << ls[j].getOuter()[i].getY()
<< ", ";
}
ss << ls[j].getOuter().front().getX() << " "
<< ls[j].getOuter().front().getY();
ss << "))";
}
ss << ")";
return ss.str();
}
// _____________________________________________________________________________
template <typename T>
inline double len(const Point<T>& g) {
@ -633,6 +787,12 @@ inline Box<T> simplify(const Box<T>& g, double d) {
return g;
}
// _____________________________________________________________________________
template <typename T>
inline RotatedBox<T> simplify(const RotatedBox<T>& g, double d) {
return g;
}
// _____________________________________________________________________________
template <typename T>
inline Line<T> simplify(const Line<T>& g, double d) {
@ -658,6 +818,15 @@ inline Line<T> simplify(const Line<T>& g, double d) {
return Line<T>{g.front(), g.back()};
}
// _____________________________________________________________________________
template <typename T>
inline Polygon<T> simplify(const Polygon<T>& g, double d) {
auto simple = simplify(g, d);
std::rotate(simple.begin(), simple.begin() + simple.size() / 2, simple.end());
simple = simplify(simple, d);
return Polygon<T>(simple);
}
// _____________________________________________________________________________
inline double distToSegment(double lax, double lay, double lbx, double lby,
double px, double py) {
@ -755,43 +924,40 @@ inline double parallelity(const Box<T>& box, const Line<T>& line) {
}
// _____________________________________________________________________________
// template <typename T>
// inline double parallelity(const Box<T>& box, const MultiLine<T>& multiline) {
// double ret = 0;
// for (const Line<T>& l : multiline) {
// ret += parallelity(box, l);
// }
template <typename T>
inline double parallelity(const Box<T>& box, const MultiLine<T>& multiline) {
double ret = 0;
for (const Line<T>& l : multiline) {
ret += parallelity(box, l);
}
// return ret / multiline.size();
// }
return ret / multiline.size();
}
// _____________________________________________________________________________
// template <typename T, template <typename> typename Geometry>
// inline RotatedBox<T> getOrientedEnvelope(Geometry<T> pol) {
// // TODO: implement this nicer, works for now, but inefficient
// // see
// // https://geidav.wordpress.com/tag/gift-wrapping/#fn-1057-FreemanShapira1975
// // for a nicer algorithm
template <typename T, template <typename> typename Geometry>
inline RotatedBox<T> getOrientedEnvelope(Geometry<T> pol) {
// TODO: implement this nicer, works for now, but inefficient
// see
// https://geidav.wordpress.com/tag/gift-wrapping/#fn-1057-FreemanShapira1975
// for a nicer algorithm
// Point<T> center;
// bgeo::centroid(pol, center);
Point<T> center = centroid(pol);
Box<T> tmpBox = getBoundingBox(pol);
double rotateDeg = 0;
// Box<T> tmpBox = getBoundingBox(pol);
// double rotateDeg = 0;
// rotate in 1 deg steps
for (int i = 1; i < 360; i += 1) {
pol = rotate(pol, 1, center);
Box<T> e = getBoundingBox(pol);
if (area(tmpBox) > area(e)) {
tmpBox = e;
rotateDeg = i;
}
}
// // rotate in 5 deg steps
// for (int i = 1; i < 360; i += 1) {
// pol = rotate(pol, 1, center);
// Box<T> e;
// bgeo::envelope(pol, e);
// if (bgeo::area(tmpBox) > bgeo::area(e)) {
// tmpBox = e;
// rotateDeg = i;
// }
// }
// return RotatedBox<T>(tmpBox, -rotateDeg, center);
// }
return RotatedBox<T>(tmpBox, -rotateDeg, center);
}
// _____________________________________________________________________________
template <typename T>
@ -849,6 +1015,90 @@ inline Box<T> getBoundingBox(const Box<T>& b) {
return b;
}
// _____________________________________________________________________________
template <typename T, template <typename> typename GeometryA>
inline Box<T> getBoundingBox(const std::vector<GeometryA<T>>& multigeo) {
Box<T> b;
b = extendBox(multigeo, b);
return b;
}
// _____________________________________________________________________________
template <typename T>
inline Polygon<T> convexHull(const Point<T>& p) {
return Polygon<T>({p});
}
// _____________________________________________________________________________
template <typename T>
inline Polygon<T> convexHull(const Box<T>& b) {
return Polygon<T>(b);
}
// _____________________________________________________________________________
template <typename T>
inline Polygon<T> convexHull(const LineSegment<T>& b) {
return Polygon<T>(Line<T>{b.first, b.second});
}
// _____________________________________________________________________________
template <typename T>
inline Polygon<T> convexHull(const RotatedBox<T>& b) {
auto p = convexHull(b.getBox());
p = rotate(p, b.getDegree(), b.getCenter());
return p;
}
// _____________________________________________________________________________
template <typename T>
inline size_t convexHullImpl(const Line<T>& a, size_t p1, size_t p2,
Line<double>* h, uint8_t d) {
// quickhull by Barber, Dobkin & Huhdanpaa
Point<T> pa;
bool found = false;
double maxDist = 0;
for (const auto& p : a) {
double tmpDist = distToSegment((*h)[p1], (*h)[p2], p);
double cp = crossProd(p, LineSegment<T>((*h)[p1], (*h)[p2]));
if (((cp > 0 && !d) || (cp < 0 && d)) && tmpDist > maxDist) {
pa = p;
found = true;
maxDist = tmpDist;
}
}
if (!found) return 0;
h->insert(h->begin() + p2 + !d, pa);
size_t in = 1 + convexHullImpl(a, p1, p2 + !d, h, d);
return in + convexHullImpl(a, p2 + in * d + 1 - 2 * d, p2 + in * d, h, d);
}
// _____________________________________________________________________________
template <typename T>
inline Polygon<T> convexHull(const Line<T>& l) {
if (l.size() == 2) return convexHull(LineSegment<T>(l[0], l[1]));
if (l.size() == 1) return convexHull(l[0]);
Point<T> left(std::numeric_limits<T>::max(), 0);
Point<T> right(std::numeric_limits<T>::min(), 0);
for (const auto& p : l) {
if (p.getX() <= left.getX()) left = p;
if (p.getX() >= right.getX()) right = p;
}
Line<T> hull{left, right};
convexHullImpl(l, 0, 1, &hull, 1);
convexHullImpl(l, 0, hull.size() - 1, &hull, 0);
return Polygon<T>(hull);
}
// _____________________________________________________________________________
template <typename T>
inline Polygon<T> convexHull(const Polygon<T>& p) {
return convexHull(p.getOuter());
}
// _____________________________________________________________________________
template <typename T>
inline Box<T> extendBox(const Line<T>& l, Box<T> b) {
@ -864,6 +1114,55 @@ inline Box<T> extendBox(const LineSegment<T>& ls, Box<T> b) {
return b;
}
// _____________________________________________________________________________
template <typename T, template <typename> typename Geometry>
inline Box<T> extendBox(const std::vector<Geometry<T>>& multigeom, Box<T> b) {
for (const auto& g : multigeom) b = extendBox(g, b);
return b;
}
// _____________________________________________________________________________
template <typename T>
inline double area(const Point<T>& b) {
UNUSED(b);
return 0;
}
// _____________________________________________________________________________
template <typename T>
inline double area(const LineSegment<T>& b) {
UNUSED(b);
return 0;
}
// _____________________________________________________________________________
template <typename T>
inline double area(const Line<T>& b) {
UNUSED(b);
return 0;
}
// _____________________________________________________________________________
template <typename T>
inline double area(const Box<T>& b) {
return (b.getUpperRight().getX() - b.getLowerLeft().getX()) *
(b.getUpperRight().getY() - b.getLowerLeft().getY());
}
// _____________________________________________________________________________
template <typename T>
inline double area(const Polygon<T>& b) {
double ret = 0;
size_t j = b.getOuter().size() - 1;
for (size_t i = 0; i < b.getOuter().size(); i++) {
ret += (b.getOuter()[j].getX() + b.getOuter()[i].getX()) *
(b.getOuter()[j].getY() - b.getOuter()[i].getY());
j = i;
}
return fabs(ret / 2.0);
}
// _____________________________________________________________________________
template <typename T>
inline double commonArea(const Box<T>& ba, const Box<T>& bb) {
@ -877,76 +1176,68 @@ inline double commonArea(const Box<T>& ba, const Box<T>& bb) {
}
// _____________________________________________________________________________
// template <typename T, template <typename> typename Geometry>
// inline RotatedBox<T> getFullEnvelope(Geometry<T> pol) {
// Point<T> center;
// bgeo::centroid(pol, center);
template <typename T, template <typename> typename Geometry>
inline RotatedBox<T> getFullEnvelope(Geometry<T> pol) {
Point<T> center = centroid(pol);
Box<T> tmpBox = getBoundingBox(pol);
double rotateDeg = 0;
// Box<T> tmpBox;
// bgeo::envelope(pol, tmpBox);
// double rotateDeg = 0;
std::vector<Polygon<T>> ml;
// MultiPolygon<T> ml;
// rotate in 5 deg steps
for (int i = 1; i < 360; i += 1) {
pol = rotate(pol, 1, center);
Polygon<T> hull = convexHull(pol);
ml.push_back(hull);
Box<T> e = getBoundingBox(pol);
if (area(tmpBox) > area(e)) {
tmpBox = e;
rotateDeg = i;
}
}
// // rotate in 5 deg steps
// for (int i = 1; i < 360; i += 1) {
// pol = rotate(pol, 1, center);
// Polygon<T> hull;
// bgeo::convex_hull(pol, hull);
// ml.push_back(hull);
// Box<T> e;
// bgeo::envelope(pol, e);
// if (bgeo::area(tmpBox) > bgeo::area(e)) {
// tmpBox = e;
// rotateDeg = i;
// }
// }
tmpBox = getBoundingBox(ml);
// bgeo::envelope(ml, tmpBox);
// return RotatedBox<T>(tmpBox, rotateDeg, center);
// }
return RotatedBox<T>(tmpBox, rotateDeg, center);
}
// _____________________________________________________________________________
// template <typename T>
// inline RotatedBox<T> getOrientedEnvelopeAvg(MultiLine<T> ml) {
// MultiLine<T> orig = ml;
// // get oriented envelope for hull
// RotatedBox<T> rbox = getFullEnvelope<T>(ml);
// Point<T> center;
// bgeo::centroid(rbox.b, center);
template <typename T>
inline RotatedBox<T> getOrientedEnvelopeAvg(MultiLine<T> ml) {
MultiLine<T> orig = ml;
// get oriented envelope for hull
RotatedBox<T> rbox = getFullEnvelope(ml);
Point<T> center = centroid(rbox.getBox());
// ml = rotate(ml, -rbox.rotateDeg - 45, center);
ml = rotate(ml, -rbox.getDegree() - 45, center);
// double bestDeg = -45;
// double score = parallelity(rbox.b, ml);
double bestDeg = -45;
double score = parallelity(rbox.b, ml);
// for (double i = -45; i <= 45; i += .5) {
// ml = rotate(ml, -.5, center);
// double p = parallelity(rbox.b, ml);
// if (parallelity(rbox.b, ml) > score) {
// bestDeg = i;
// score = p;
// }
// }
for (double i = -45; i <= 45; i += .5) {
ml = rotate(ml, -.5, center);
double p = parallelity(rbox.b, ml);
if (parallelity(rbox.b, ml) > score) {
bestDeg = i;
score = p;
}
}
// rbox.rotateDeg += bestDeg;
rbox.rotateDeg += bestDeg;
// // move the box along 45deg angles from its origin until it fits the ml
// // = until the intersection of its hull and the box is largest
// Polygon<T> p = rbox.getPolygon();
// p = rotate(p, -rbox.rotateDeg, rbox.center);
// move the box along 45deg angles from its origin until it fits the ml
// = until the intersection of its hull and the box is largest
Polygon<T> p = convexHull(rbox);
p = rotate(p, -rbox.getDegree(), rbox.getCenter());
// Polygon<T> hull;
// bgeo::convex_hull(orig, hull);
// hull = rotate(hull, -rbox.rotateDeg, rbox.center);
Polygon<T> hull = convexHull(orig);
hull = rotate(hull, -rbox.getDegree(), rbox.getCenter());
// Box<T> box;
// bgeo::envelope(hull, box);
// rbox = RotatedBox<T>(box, rbox.rotateDeg, rbox.center);
Box<T> box = getBoundingBox(hull);
rbox = RotatedBox<T>(box, rbox.getDegree(), rbox.getCenter());
// return rbox;
// }
return rbox;
}
// _____________________________________________________________________________
template <typename T>

10
src/util/geo/Line.h
View file

@ -11,12 +11,16 @@
namespace util {
namespace geo {
template<typename T>
using Line = std::vector<Point<T>>;
template <typename T>
class Line : public std::vector<Point<T>> {
using std::vector<Point<T>>::vector;
};
template<typename T>
template <typename T>
using LineSegment = std::pair<Point<T>, Point<T>>;
template <typename T>
using MultiLine = std::vector<Line<T>>;
} // namespace geo
} // namespace util

3
src/util/geo/Point.h
View file

@ -39,6 +39,9 @@ class Point {
T _x, _y;
};
template <typename T>
using MultiPoint = std::vector<Point<T>>;
} // namespace geo
} // namespace util

3
src/util/geo/PolyLine.h
View file

@ -17,6 +17,8 @@ namespace geo {
static const double MAX_EQ_DISTANCE = 15;
static const double AVERAGING_STEP = 20;
// legacy code, will be removed in the future
template <typename T>
struct LinePoint {
LinePoint() : lastIndex(0), totalPos(-1), p() {}
@ -46,7 +48,6 @@ struct SharedSegments {
std::vector<SharedSegment<T>> segments;
};
// TODO: maybe let this class inherit from a more generic geometry class
template <typename T>
class PolyLine {
public:

12
src/util/geo/PolyLine.tpp
View file

@ -93,11 +93,11 @@ void PolyLine<T>::offsetPerp(double units) {
n1 = n1 / n;
n2 = n2 / n;
lastP.template set<0>(lastP.getX() + (n1 * units));
lastP.template set<1>(lastP.getY() + (n2 * units));
lastP.setX(lastP.getX() + (n1 * units));
lastP.setY(lastP.getY() + (n2 * units));
curP.template set<0>(curP.getX() + (n1 * units));
curP.template set<1>(curP.getY() + (n2 * units));
curP.setX(curP.getX() + (n1 * units));
curP.setY(curP.getY() + (n2 * units));
if (lastIns && befLastIns &&
lineIntersects(*lastIns, *befLastIns, lastP, curP)) {
@ -460,8 +460,8 @@ bool PolyLine<T>::contains(const PolyLine<T>& rhs, double dmax) const {
template <typename T>
void PolyLine<T>::move(double vx, double vy) {
for (size_t i = 0; i < _line.size(); i++) {
_line[i].set<0>(_line[i].getX() + vx);
_line[i].set<1>(_line[i].getY() + vy);
_line[i].setX(_line[i].getX() + vx);
_line[i].setY(_line[i].getY() + vy);
}
}

9
src/util/geo/Polygon.h
View file

@ -25,13 +25,16 @@ class Polygon {
b.getUpperRight(),
Point<T>(b.getLowerLeft().getX(), b.getUpperRight().getY())}) {}
const std::vector<Point<T>>& getOuter() const { return _outer; }
std::vector<Point<T>>& getOuter() { return _outer; }
const Line<T>& getOuter() const { return _outer; }
Line<T>& getOuter() { return _outer; }
private:
std::vector<Point<T>> _outer;
Line<T> _outer;
};
template <typename T>
using MultiPolyon = std::vector<Polygon<T>>;
} // namespace geo
} // namespace util

263
src/util/tests/TestMain.cpp
View file

@ -272,38 +272,38 @@ CASE("frechet distance") {
// ___________________________________________________________________________
CASE("geo box alignment") {
// Line<double> a;
// a.push_back(Point<double>(1, 1));
// a.push_back(Point<double>(1, 2));
Line<double> a;
a.push_back(Point<double>(1, 1));
a.push_back(Point<double>(1, 2));
// Line<double> b;
// b.push_back(Point<double>(1, 2));
// b.push_back(Point<double>(2, 2));
Line<double> b;
b.push_back(Point<double>(1, 2));
b.push_back(Point<double>(2, 2));
// Line<double> c;
// c.push_back(Point<double>(2, 2));
// c.push_back(Point<double>(2, 1));
Line<double> c;
c.push_back(Point<double>(2, 2));
c.push_back(Point<double>(2, 1));
// Line<double> d;
// d.push_back(Point<double>(2, 1));
// d.push_back(Point<double>(1, 1));
Line<double> d;
d.push_back(Point<double>(2, 1));
d.push_back(Point<double>(1, 1));
// Box<double> box(Point<double>(2, 3), Point<double>(5, 4));
// MultiLine<double> ml;
// ml.push_back(a);
// ml.push_back(b);
// ml.push_back(c);
// ml.push_back(d);
Box<double> box(Point<double>(2, 3), Point<double>(5, 4));
MultiLine<double> ml;
ml.push_back(a);
ml.push_back(b);
ml.push_back(c);
ml.push_back(d);
// EXPECT(parallelity(box, ml) == approx(1));
// ml = rotate(ml, 45);
// EXPECT(parallelity(box, ml) == approx(0));
// ml = rotate(ml, 45);
// EXPECT(parallelity(box, ml) == approx(1));
// ml = rotate(ml, 45);
// EXPECT(parallelity(box, ml) == approx(0));
// ml = rotate(ml, 45);
// EXPECT(parallelity(box, ml) == approx(1));
EXPECT(parallelity(box, ml) == approx(1));
ml = rotate(ml, 45);
EXPECT(parallelity(box, ml) == approx(0));
ml = rotate(ml, 45);
EXPECT(parallelity(box, ml) == approx(1));
ml = rotate(ml, 45);
EXPECT(parallelity(box, ml) == approx(0));
ml = rotate(ml, 45);
EXPECT(parallelity(box, ml) == approx(1));
},
// ___________________________________________________________________________
@ -944,7 +944,7 @@ CASE("geometry") {
EXPECT(geo::intersects(Box<double>(Point<double>(0, 0), Point<double>(1.5, 1.5)), boxa));
Polygon<double> poly({Point<double>(1, 1), Point<double>(3, 2), Point<double>(4, 3), Point<double>(6, 3), Point<double>(5, 1)});
EXPECT(geo::getWKT(poly) == "POLYGON ((1 1, 3 2, 4 3, 6 3, 5 1))");
EXPECT(geo::getWKT(poly) == "POLYGON ((1 1, 3 2, 4 3, 6 3, 5 1, 1 1))");
EXPECT(geo::contains(Point<double>(4, 2), poly));
EXPECT(!geo::contains(Point<double>(3, 3), poly));
EXPECT(geo::contains(Point<double>(1, 1), poly));
@ -977,8 +977,215 @@ CASE("geometry") {
EXPECT(geo::rotate(rotLine, 90, Point<double>(2, 2))[1].getX() == approx(3));
EXPECT(geo::rotate(rotLine, 90, Point<double>(2, 2))[1].getY() == approx(1));
MultiLine<double> multiRotLine({{{1, 1}, {3, 3}}, {{1, 3}, {3, 1}}});
EXPECT(geo::rotate(multiRotLine, 90, Point<double>(2, 2))[0][0].getX() == approx(1));
EXPECT(geo::rotate(multiRotLine, 90, Point<double>(2, 2))[0][0].getY() == approx(3));
EXPECT(geo::rotate(multiRotLine, 90, Point<double>(2, 2))[0][1].getX() == approx(3));
EXPECT(geo::rotate(multiRotLine, 90, Point<double>(2, 2))[0][1].getY() == approx(1));
EXPECT(geo::rotate(multiRotLine, 90, Point<double>(2, 2))[1][0].getX() == approx(3));
EXPECT(geo::rotate(multiRotLine, 90, Point<double>(2, 2))[1][0].getY() == approx(3));
EXPECT(geo::rotate(multiRotLine, 90, Point<double>(2, 2))[1][1].getX() == approx(1));
EXPECT(geo::rotate(multiRotLine, 90, Point<double>(2, 2))[1][1].getY() == approx(1));
EXPECT(geo::getWKT(multiRotLine) == "MULTILINESTRING ((1 1, 3 3), (1 3, 3 1))");
EXPECT(geo::contains(multiRotLine[0], geo::move(geo::move(multiRotLine, 1.0, 2.0), -1.0, -2.0)[0]));
EXPECT(geo::contains(multiRotLine, geo::getBoundingBox(Line<double>{{1, 1}, {3, 3}, {1, 3}, {3, 1}})));
EXPECT(geo::contains(getBoundingBox(multiRotLine), geo::getBoundingBox(Line<double>{{1, 1}, {3, 3}, {1, 3}, {3, 1}})));
EXPECT(geo::contains(geo::getBoundingBox(Line<double>{{1, 1}, {3, 3}, {1, 3}, {3, 1}}), getBoundingBox(multiRotLine)));
EXPECT(geo::dist(geo::centroid(rotP), rotP) == approx(0));
EXPECT(geo::dist(geo::centroid(rotLine), rotP) == approx(0));
EXPECT(geo::dist(geo::centroid(polybox), Point<double>(3.5, 2.5)) == approx(0));
EXPECT(geo::dist(geo::centroid(Polygon<double>({{0, 0}, {3, 4}, {4,3}})), Point<double>(7.0/3.0,7.0/3.0)) == approx(0));
auto polyy = Polygon<double>({{0, 0}, {3, 4}, {4,3}});
MultiPolyon<double> mpoly{polyy, polyy};
EXPECT(geo::getWKT(polyy) == "POLYGON ((0 0, 3 4, 4 3, 0 0))");
EXPECT(geo::getWKT(mpoly) == "MULTIPOLYGON (((0 0, 3 4, 4 3, 0 0)), ((0 0, 3 4, 4 3, 0 0)))");
auto hull = geo::convexHull(Line<double>{{0.1, 3}, {1, 1}, {2, 2}, {4, 4}, {0, 0}, {1, 2}, {3, 1}, {3, 3}});
EXPECT(hull.getOuter().size() == size_t(4));
EXPECT(hull.getOuter()[0].getX() == approx(0));
EXPECT(hull.getOuter()[0].getY() == approx(0));
EXPECT(hull.getOuter()[1].getX() == approx(3));
EXPECT(hull.getOuter()[1].getY() == approx(1));
EXPECT(hull.getOuter()[2].getX() == approx(4));
EXPECT(hull.getOuter()[2].getY() == approx(4));
EXPECT(hull.getOuter()[3].getX() == approx(0.1));
EXPECT(hull.getOuter()[3].getY() == approx(3));
EXPECT(geo::contains(geo::convexHull(geo::getBoundingBox(poly)), geo::getBoundingBox(poly)));
EXPECT(geo::contains(geo::getBoundingBox(poly), geo::convexHull(geo::getBoundingBox(poly))));
auto hull2 = geo::convexHull(Line<double>{{0.1, 3}, {1, 1}, {2, 2}, {4, 4}, {0, 0}, {1, 2}, {3, 1}, {3, 3}, {-0.1, 1}});
EXPECT(hull2.getOuter().size() == size_t(5));
EXPECT(hull2.getOuter()[0].getX() == approx(-.1));
EXPECT(hull2.getOuter()[0].getY() == approx(1));
EXPECT(hull2.getOuter()[1].getX() == approx(0));
EXPECT(hull2.getOuter()[1].getY() == approx(0));
EXPECT(hull2.getOuter()[2].getX() == approx(3));
EXPECT(hull2.getOuter()[2].getY() == approx(1));
EXPECT(hull2.getOuter()[3].getX() == approx(4));
EXPECT(hull2.getOuter()[3].getY() == approx(4));
EXPECT(hull2.getOuter()[4].getX() == approx(0.1));
EXPECT(hull2.getOuter()[4].getY() == approx(3));
auto hull3 = geo::convexHull(Line<double>{{0.1, 3}, {4, 4}, {0, 0}, {1, 2}, {3, 1}});
EXPECT(hull3.getOuter().size() == size_t(4));
EXPECT(hull3.getOuter()[0].getX() == approx(0));
EXPECT(hull3.getOuter()[0].getY() == approx(0));
EXPECT(hull3.getOuter()[1].getX() == approx(3));
EXPECT(hull3.getOuter()[1].getY() == approx(1));
EXPECT(hull3.getOuter()[2].getX() == approx(4));
EXPECT(hull3.getOuter()[2].getY() == approx(4));
EXPECT(hull3.getOuter()[3].getX() == approx(0.1));
EXPECT(hull3.getOuter()[3].getY() == approx(3));
hull3 = geo::convexHull(Line<double>{{0.1, 3}, {4, 4}, {2, 1}, {3, 2}, {0, 0}, {1, 2}, {3, 1}});
EXPECT(hull3.getOuter().size() == size_t(4));
EXPECT(hull3.getOuter()[0].getX() == approx(0));
EXPECT(hull3.getOuter()[0].getY() == approx(0));
EXPECT(hull3.getOuter()[1].getX() == approx(3));
EXPECT(hull3.getOuter()[1].getY() == approx(1));
EXPECT(hull3.getOuter()[2].getX() == approx(4));
EXPECT(hull3.getOuter()[2].getY() == approx(4));
EXPECT(hull3.getOuter()[3].getX() == approx(0.1));
EXPECT(hull3.getOuter()[3].getY() == approx(3));
hull3 = geo::convexHull(Line<double>{{4, 4}, {1, 2}, {2, 1}, {3, 2}, {0.1, 3}, {0, 0}, {1, 2}, {3, 1}});
EXPECT(hull3.getOuter().size() == size_t(4));
EXPECT(hull3.getOuter()[0].getX() == approx(0));
EXPECT(hull3.getOuter()[0].getY() == approx(0));
EXPECT(hull3.getOuter()[1].getX() == approx(3));
EXPECT(hull3.getOuter()[1].getY() == approx(1));
EXPECT(hull3.getOuter()[2].getX() == approx(4));
EXPECT(hull3.getOuter()[2].getY() == approx(4));
EXPECT(hull3.getOuter()[3].getX() == approx(0.1));
EXPECT(hull3.getOuter()[3].getY() == approx(3));
hull3 = geo::convexHull(Line<double>{{4, 4}, {1, 2}, {3, 1}});
EXPECT(hull3.getOuter().size() == size_t(3));
EXPECT(hull3.getOuter()[0].getX() == approx(1));
EXPECT(hull3.getOuter()[0].getY() == approx(2));
EXPECT(hull3.getOuter()[1].getX() == approx(3));
EXPECT(hull3.getOuter()[1].getY() == approx(1));
EXPECT(hull3.getOuter()[2].getX() == approx(4));
EXPECT(hull3.getOuter()[2].getY() == approx(4));
hull3 = geo::convexHull(Line<double>{{4, 4}, {1, 2}, {3, 10}});
EXPECT(hull3.getOuter().size() == size_t(3));
EXPECT(hull3.getOuter()[0].getX() == approx(1));
EXPECT(hull3.getOuter()[0].getY() == approx(2));
EXPECT(hull3.getOuter()[1].getX() == approx(4));
EXPECT(hull3.getOuter()[1].getY() == approx(4));
EXPECT(hull3.getOuter()[2].getX() == approx(3));
EXPECT(hull3.getOuter()[2].getY() == approx(10));
Line<double> test{{0.3215348546593775, 0.03629583077160248},
{0.02402358131857918, -0.2356728797179394},
{0.04590851212470659, -0.4156409924995536},
{0.3218384001607433, 0.1379850698988746},
{0.11506479756447, -0.1059521474930943},
{0.2622539999543261, -0.29702873322836},
{-0.161920957418085, -0.4055339716426413},
{0.1905378631228002, 0.3698601009043493},
{0.2387090918968516, -0.01629827079949742},
{0.07495888748668034, -0.1659825110491202},
{0.3319341836794598, -0.1821814101954749},
{0.07703635755650362, -0.2499430638271785},
{0.2069242999022122, -0.2232970760420869},
{0.04604079532068295, -0.1923573186549892},
{0.05054295812784038, 0.4754929463150845},
{-0.3900589168910486, 0.2797829520700341},
{0.3120693385713448, -0.0506329867529059},
{0.01138812723698857, 0.4002504701728471},
{0.009645149586391732, 0.1060251100976254},
{-0.03597933197019559, 0.2953639456959105},
{0.1818290866742182, 0.001454397571696298},
{0.444056063372694, 0.2502497166863175},
{-0.05301752458607545, -0.06553921621808712},
{0.4823896228171788, -0.4776170002088109},
{-0.3089226845734964, -0.06356112199235814},
{-0.271780741188471, 0.1810810595574612},
{0.4293626522918815, 0.2980897964891882},
{-0.004796652127799228, 0.382663812844701},
{0.430695573269106, -0.2995073500084759},
{0.1799668387323309, -0.2973467472915973},
{0.4932166845474547, 0.4928094162538735},
{-0.3521487911717489, 0.4352656197131292},
{-0.4907368011686362, 0.1865826865533206},
{-0.1047924716070224, -0.247073392148198},
{0.4374961861758457, -0.001606279519951237},
{0.003256207800708899, -0.2729194320486108},
{0.04310378203457577, 0.4452604050238248},
{0.4916198379282093, -0.345391701297268},
{0.001675087028811806, 0.1531837672490476},
{-0.4404289572876217, -0.2894855991839297}
};
hull3 = geo::convexHull(test);
EXPECT(geo::contains(test, hull3));
EXPECT(hull3.getOuter().size() == size_t(8));
EXPECT(geo::contains(Polygon<double>({{-0.161920957418085, -0.4055339716426413},
{0.05054295812784038, 0.4754929463150845},
{0.4823896228171788, -0.4776170002088109},
{0.4932166845474547, 0.4928094162538735},
{-0.3521487911717489, 0.4352656197131292},
{-0.4907368011686362, 0.1865826865533206},
{0.4916198379282093, -0.345391701297268},
{-0.4404289572876217,
-0.2894855991839297}}), hull3));
EXPECT(geo::contains(hull3, Polygon<double>({{-0.161920957418085, -0.4055339716426413},
{0.05054295812784038, 0.4754929463150845},
{0.4823896228171788, -0.4776170002088109},
{0.4932166845474547, 0.4928094162538735},
{-0.3521487911717489, 0.4352656197131292},
{-0.4907368011686362, 0.1865826865533206},
{0.4916198379282093, -0.345391701297268},
{-0.4404289572876217,
-0.2894855991839297}})));
hull3 = geo::convexHull(Line<double>{{3, 6}, {8, 10}, {3, 5}, {20, -10}, {-4, 5}, {10, 2}, {5, 1}, {45, 1}, {30, -9}, {3, 14}, {25, -5.5}});
EXPECT(hull3.getOuter().size() == size_t(5));
EXPECT(hull3.getOuter()[0].getX() == approx(-4));
EXPECT(hull3.getOuter()[0].getY() == approx(5));
EXPECT(hull3.getOuter()[1].getX() == approx(20));
EXPECT(hull3.getOuter()[1].getY() == approx(-10));
EXPECT(hull3.getOuter()[2].getX() == approx(30));
EXPECT(hull3.getOuter()[2].getY() == approx(-9));
EXPECT(hull3.getOuter()[3].getX() == approx(45));
EXPECT(hull3.getOuter()[3].getY() == approx(1));
EXPECT(hull3.getOuter()[4].getX() == approx(3));
EXPECT(hull3.getOuter()[4].getY() == approx(14));
hull3 = geo::convexHull(Line<double>{{7, 7}, {7, -7}, {-7, -7}, {-7, 7}, {9, 0}, {-9, 0}, {0, 9}, {0, -9}});
EXPECT(hull3.getOuter().size() == size_t(8));
EXPECT(geo::contains(geo::Polygon<double>({{-9, 0}, {-7, -7}, {0, -9}, {7, -7}, {9, 0}, {7, 7}, {0, 9}, {-7, 7}}), hull3));
EXPECT(geo::contains(hull3, geo::Polygon<double>({{-9, 0}, {-7, -7}, {0, -9}, {7, -7}, {9, 0}, {7, 7}, {0, 9}, {-7, 7}})));
hull3 = geo::convexHull(Line<double>{{7, 7}, {7, -7}, {-7, -7}, {-7, 7}, {9, 0}, {-9, 0}, {0, 9}, {0, -9}, {0, 0}, {1, 2}, {-2, 1}, {-1, -1}, {3, 4}, {4, 3}, {-5, 4}, {6, 5}});
EXPECT(hull3.getOuter().size() == size_t(8));
EXPECT(geo::contains(geo::Polygon<double>({{-9, 0}, {-7, -7}, {0, -9}, {7, -7}, {9, 0}, {7, 7}, {0, 9}, {-7, 7}}), hull3));
EXPECT(geo::contains(hull3, geo::Polygon<double>({{-9, 0}, {-7, -7}, {0, -9}, {7, -7}, {9, 0}, {7, 7}, {0, 9}, {-7, 7}})));
hull3 = geo::convexHull(Line<double>{{0, 0}, {1, 2}, {-2, 1}, {-1, -1}, {3, 4}, {4, 3}, {-5, 4}, {6, 5}, {7, 7}, {7, -7}, {-7, -7}, {-7, 7}, {9, 0}, {-9, 0}, {0, 9}, {0, -9}, {-8, 0}, {8, 0}, {-7, 0}, {7, 0}, {-6, 0}, {6, 0}, {-5, 0}, {5, 0}, {-4, 0}, {4, 0}, {-3, 0}, {3, 0}, {-2, 0}, {2, 0}, {-1, 0}, {1, 0}, {0, -8}, {0, 8}, {0, -7}, {0, 7}, {0, -6}, {0, 6}, {0, -5}, {0, 5}, {0, -4}, {0, 4}, {0, -3}, {0, 3}, {0, -2}, {0, 2}, {0, -1}, {0, 1}, {1, 1}, {2, 2}, {3, 3}, {4, 4}, {5, 5}, {6, 6}, {1, -1}, {2, -2}, {3, -3}, {4, -4}, {5, -5}, {6, -6}, {-1, 1}, {-2, 2}, {-3, 3}, {-4, 4}, {-5, 5}, {-6, 6}, {-1, -1}, {-2, -2}, {-3, -3}, {-4, -4}, {-5, -5}, {-6, -6}});
EXPECT(hull3.getOuter().size() == size_t(8));
EXPECT(geo::contains(geo::Polygon<double>({{-9, 0}, {-7, -7}, {0, -9}, {7, -7}, {9, 0}, {7, 7}, {0, 9}, {-7, 7}}), hull3));
EXPECT(geo::contains(hull3, geo::Polygon<double>({{-9, 0}, {-7, -7}, {0, -9}, {7, -7}, {9, 0}, {7, 7}, {0, 9}, {-7, 7}})));
EXPECT(geo::area(geo::Point<double>(1, 2)) == approx(0));
EXPECT(geo::area(geo::Line<double>{{1, 2}, {2, 5}}) == approx(0));
EXPECT(geo::area(geo::Box<double>({0, 0}, {1, 1})) == approx(1));
EXPECT(geo::area(geo::Box<double>({1, 1}, {1, 1})) == approx(0));
EXPECT(geo::area(geo::Box<double>({0, 0}, {2, 2})) == approx(4));
EXPECT(geo::area(geo::Polygon<double>({{0, 0}, {1, 0}, {1, 1}, {0, 1}})) == approx(1));
EXPECT(geo::area(geo::Polygon<double>({{0, 0}, {1, 0}, {1, 1}})) == approx(0.5));
auto obox = geo::getOrientedEnvelope(geo::Line<double>{{0, 0}, {1, 1}, {1.5, 0.5}});
EXPECT(geo::contains(geo::convexHull(obox), geo::Polygon<double>({{0.0, 0.0}, {1.0, 1.0}, {1.5, 0.5}, {0.5, -0.5}})));
EXPECT(geo::contains(geo::Polygon<double>({{0.0, 0.0}, {1.0, 1.0}, {1.5, 0.5}, {0.5, -0.5}}), geo::convexHull(obox)));
}
};