Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement getOrientedEnvelope using rotating calipers #1

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Implement getOrientedEnvelope using rotating calipers #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
141 changes: 105 additions & 36 deletions geo/Geo.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,17 +119,29 @@ inline Point<T> rotate(const Point<T>& p, double deg) {

// _____________________________________________________________________________
template <typename T>
inline Point<T> rotate(Point<T> p, double deg, const Point<T>& c) {
deg *= -RAD;
double si = sin(deg);
double co = cos(deg);
inline Point<T> rotateRAD(const Point<T>& p, double deg) {
UNUSED(deg);
return p;
}

// _____________________________________________________________________________
template <typename T>
inline Point<T> rotateRAD(Point<T> p, double rad, const Point<T>& c) {
double si = sin(rad);
double co = cos(rad);
p = p - c;

return Point<T>(p.getX() * co - p.getY() * si,
p.getX() * si + p.getY() * co) +
c;
}

// _____________________________________________________________________________
template <typename T>
inline Point<T> rotate(Point<T> p, double deg, const Point<T>& c) {
return rotateRAD(p, deg * -RAD, c);
}

// _____________________________________________________________________________
template <typename T>
inline LineSegment<T> rotate(LineSegment<T> geo, double deg,
Expand All @@ -139,10 +151,25 @@ inline LineSegment<T> rotate(LineSegment<T> geo, double deg,
return geo;
}

// _____________________________________________________________________________
template <typename T>
inline LineSegment<T> rotateRAD(LineSegment<T> geo, double deg,
const Point<T>& c) {
geo.first = rotateRAD(geo.first, deg, c);
geo.second = rotateRAD(geo.second, deg, c);
return geo;
}

// _____________________________________________________________________________
template <typename T>
inline LineSegment<T> rotate(LineSegment<T> geo, double deg) {
return (geo, deg, centroid(geo));
return rotate(geo, deg, centroid(geo));
}

// _____________________________________________________________________________
template <typename T>
inline LineSegment<T> rotateRAD(LineSegment<T> geo, double deg) {
return rotateRAD(geo, deg, centroid(geo));
}

// _____________________________________________________________________________
Expand All @@ -152,11 +179,32 @@ inline Line<T> rotate(Line<T> geo, double deg, const Point<T>& c) {
return geo;
}

// _____________________________________________________________________________
template <typename T>
inline Line<T> rotateRAD(Line<T> geo, double deg, const Point<T>& c) {
for (size_t i = 0; i < geo.size(); i++) geo[i] = rotateRAD(geo[i], deg, c);
return geo;
}

// _____________________________________________________________________________
template <typename T>
inline Polygon<T> rotate(Polygon<T> geo, double deg, const Point<T>& c) {
for (size_t i = 0; i < geo.getOuter().size(); i++)
geo.getOuter()[i] = rotate(geo.getOuter()[i], deg, c);
for (size_t i = 0; i < geo.getInners().size(); i++)
for (size_t j = 0; j < geo.getInners()[i].size(); j++)
geo.getInners()[i][j] = rotate(geo.getInners()[i][j], deg, c);
return geo;
}

// _____________________________________________________________________________
template <typename T>
inline Polygon<T> rotateRAD(Polygon<T> geo, double deg, const Point<T>& c) {
for (size_t i = 0; i < geo.getOuter().size(); i++)
geo.getOuter()[i] = rotateRAD(geo.getOuter()[i], deg, c);
for (size_t i = 0; i < geo.getInners().size(); i++)
for (size_t j = 0; j < geo.getInners()[i].size(); j++)
geo.getInners()[i][j] = rotateRAD(geo.getInners()[i][j], deg, c);
return geo;
}

Expand All @@ -169,6 +217,15 @@ inline std::vector<Geometry<T>> rotate(std::vector<Geometry<T>> multigeo,
return multigeo;
}

// _____________________________________________________________________________
template <template <typename> class Geometry, typename T>
inline std::vector<Geometry<T>> rotateRAD(std::vector<Geometry<T>> multigeo,
double deg, const Point<T>& c) {
for (size_t i = 0; i < multigeo.size(); i++)
multigeo[i] = rotateRAD(multigeo[i], deg, c);
return multigeo;
}

// _____________________________________________________________________________
template <template <typename> class Geometry, typename T>
inline std::vector<Geometry<T>> rotate(std::vector<Geometry<T>> multigeo,
Expand All @@ -179,6 +236,16 @@ inline std::vector<Geometry<T>> rotate(std::vector<Geometry<T>> multigeo,
return multigeo;
}

// _____________________________________________________________________________
template <template <typename> class Geometry, typename T>
inline std::vector<Geometry<T>> rotateRAD(std::vector<Geometry<T>> multigeo,
double deg) {
auto c = centroid(multigeo);
for (size_t i = 0; i < multigeo.size(); i++)
multigeo[i] = rotateRAD(multigeo[i], deg, c);
return multigeo;
}

// _____________________________________________________________________________
template <typename T>
inline Point<T> move(const Point<T>& geo, double x, double y) {
Expand Down Expand Up @@ -1445,50 +1512,44 @@ inline double parallelity(const Box<T>& box, const MultiLine<T>& multiline) {

// _____________________________________________________________________________
template <template <typename> class Geometry, typename T>
inline RotatedBox<T> getOrientedEnvelope(std::vector<Geometry<T>> pol,
double step) {
// TODO: implement this nicer, works for now, but inefficient
// see
// https://geidav.wordpress.com/tag/gift-wrapping/#fn-1057-FreemanShapira1975
// for a nicer algorithm

inline RotatedBox<T> getOrientedEnvelope(std::vector<Geometry<T>> pol) {
Point<T> center = centroid(pol);
Box<T> tmpBox = getBoundingBox(pol);
double rotateDeg = 0;

// rotate in steps
double i = step;
while (i < 360) {
pol = rotate(pol, step, center);
Box<T> e = getBoundingBox(pol);
Line<T> hull = convexHull(pol).getOuter();
double rotateAngle = 0;

// check each segment
for (size_t i = 1; i < hull.size(); i++) {
// rotate segment such that it is parallel to the x axis
const auto s = hull[i] - hull[i-1];
const double angle = -std::atan2(s.getY(), s.getX());
const auto p = rotateRAD(pol, angle, center);
Box<T> e = getBoundingBox(p);
if (area(tmpBox) > area(e)) {
tmpBox = e;
rotateDeg = i;
rotateAngle = angle;
}
}
// Check segment between the ends of the hull line
if (hull[0] != hull[hull.size()-1]) {
const auto s = hull[0] - hull[hull.size() - 1];
const double angle = -std::atan2(s.getY(), s.getX());
const auto p = rotateRAD(pol, angle, center);
Box<T> e = getBoundingBox(p);
if (area(tmpBox) > area(e)) {
tmpBox = e;
rotateAngle = angle;
}
i+= step;
}

return RotatedBox<T>(tmpBox, -rotateDeg, center);
}

// _____________________________________________________________________________
template <template <typename> class Geometry, typename T>
inline RotatedBox<T> getOrientedEnvelope(std::vector<Geometry<T>> pol) {
return getOrientedEnvelope(pol, 1);
}

// _____________________________________________________________________________
template <template <typename> class Geometry, typename T>
inline RotatedBox<T> getOrientedEnvelope(Geometry<T> pol, double step) {
std::vector<Geometry<T>> mult{pol};
return getOrientedEnvelope(mult, step);
return RotatedBox<T>(tmpBox, -rotateAngle * -IRAD, center);
}

// _____________________________________________________________________________
template <template <typename> class Geometry, typename T>
inline RotatedBox<T> getOrientedEnvelope(Geometry<T> pol) {
std::vector<Geometry<T>> mult{pol};
return getOrientedEnvelope(mult, 1);
return getOrientedEnvelope(mult);
}

// _____________________________________________________________________________
Expand Down Expand Up @@ -1715,6 +1776,14 @@ inline Polygon<T> convexHull(const Polygon<T>& p) {
return convexHull(p.getOuter());
}

// _____________________________________________________________________________
template <typename T>
inline Polygon<T> convexHull(const MultiPolygon<T>& ps) {
MultiPoint<T> mp;
for (const auto& p : ps) mp.insert(mp.end(), p.getOuter().begin(), p.getOuter().end());
return convexHull(mp);
}

// _____________________________________________________________________________
template <typename T>
inline Polygon<T> convexHull(const MultiLine<T>& ls) {
Expand Down
13 changes: 13 additions & 0 deletions tests/TestMain.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1791,6 +1791,19 @@ int main(int argc, char** argv) {
geo::Polygon<double>({{0.0, 0.0}, {1.0, 1.0}, {1.5, 0.5}, {0.5, -0.5}}),
geo::convexHull(obox)));

obox =
geo::getOrientedEnvelope(geo::Polygon<double>({{0.0, 1.0}, {0.9, 1.1},
{1.0, 2.0}, {1.1, 1.1},
{2.0, 1.0}, {1.1, 0.9},
{1.0, 0.0}, {0.9, 0.9}}));
std::cerr << geo::getWKT(geo::convexHull(obox)) << std::endl << std::flush;
TEST(geo::contains(
geo::convexHull(obox),
geo::Polygon<double>({{0.0, 1.0}, {1.0, 2.0}, {2.0, 1.0}, {1.0, 0.0}})));
TEST(geo::contains(
geo::Polygon<double>({{0.0, 1.0}, {1.0, 2.0}, {2.0, 1.0}, {1.0, 0.0}}),
geo::convexHull(obox)));

TEST(geo::dist(geo::LineSegment<double>{{1, 1}, {3, 1}},
geo::LineSegment<double>{{2, 2}, {2, 0}}) == approx(0));
TEST(geo::dist(geo::LineSegment<double>{{1, 1}, {3, 1}},
Expand Down