Skip to content

Commit

Permalink
adopt convex hull from concaveman (#28)
Browse files Browse the repository at this point in the history
* not ready

* fix

* use stable predicates

* fix

* fix

* fix

* not ready

* fix

* fix

* reset

* init orient 2d

* add macro

* update

* clean

* fix

* fix

* not ready

* fix compile

* fix

* fix

* fix

* fix

* fix

* fix

---------

Co-authored-by: tang zhixiong <[email protected]>
  • Loading branch information
district10 and zhixiong-tang authored Dec 14, 2024
1 parent 078c7da commit 5cc777b
Show file tree
Hide file tree
Showing 12 changed files with 518 additions and 13 deletions.
4 changes: 4 additions & 0 deletions docs/about/release-notes.md
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,10 @@ To upgrade `concave_hull` to the latest version, use pip:
pip install -U concave_hull
```

## Version 0.0.9 (2024-12-15)

* Fix broken package (not well tested)

## Version 0.0.8 (2024-10-03)

* Add typing stubs
Expand Down
2 changes: 1 addition & 1 deletion pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ build-backend = "scikit_build_core.build"

[project]
name = "concave_hull"
version = "0.0.8"
version = "0.0.9"
url = "https://concave-hull.readthedocs.io"
description="A very fast 2D concave hull algorithm"
readme = "README.md"
Expand Down
12 changes: 6 additions & 6 deletions src/concave_hull/__init__.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
from typing import List, Tuple, Union
from __future__ import annotations

import numpy as np

Expand All @@ -10,11 +10,11 @@


def concave_hull_indexes(
points: Union[np.ndarray, List, Tuple],
points: np.ndarray | list | tuple,
*,
concavity: float = 2.0,
length_threshold: float = 0.0,
convex_hull_indexes: np.ndarray = None, # noqa
convex_hull_indexes: np.ndarray | None = None, # noqa
is_wgs84: bool = False,
):
"""
Expand Down Expand Up @@ -45,7 +45,7 @@ def concave_hull_indexes(


def concave_hull(
points: Union[np.ndarray, List, Tuple],
points: np.ndarray | list | tuple,
*,
concavity: float = 2.0,
length_threshold: float = 0.0,
Expand All @@ -67,7 +67,7 @@ def concave_hull(


def convex_hull_indexes(
points: Union[np.ndarray, List, Tuple],
points: np.ndarray | list | tuple,
*,
include_colinear: bool = False,
order_only: bool = False,
Expand All @@ -81,7 +81,7 @@ def convex_hull_indexes(


def convex_hull(
points: Union[np.ndarray, List, Tuple],
points: np.ndarray | list | tuple,
*,
include_colinear: bool = False,
order_only: bool = False,
Expand Down
2 changes: 1 addition & 1 deletion src/concave_hull/_core.pyi
Original file line number Diff line number Diff line change
Expand Up @@ -63,4 +63,4 @@ def wgs84_to_east_north(
documents here: https://github.com/mapbox/cheap-ruler
"""

__version__: str = "0.0.8"
__version__: str = "0.0.9"
4 changes: 2 additions & 2 deletions src/concaveman.h
Original file line number Diff line number Diff line change
Expand Up @@ -454,7 +454,7 @@ std::vector<int> concaveman_indexes(
// index the points with an R-tree
rtree<T, 2, MAX_CHILDREN, point_type> tree;
for (size_t index = 0, N = points.size(); index < N; index++) {
point_type p{points[index][0], points[index][1], (T)index};
point_type p{points[index][0], points[index][1], static_cast<T>(index)};
tree.insert(std::move(p), {p[0], p[1], p[0], p[1]});
}

Expand All @@ -467,7 +467,7 @@ std::vector<int> concaveman_indexes(
// queue with the nodes
for (auto &idx : hull) {
auto &pp = points[idx];
point_type p{pp[0], pp[1], (T)idx};
point_type p{pp[0], pp[1], static_cast<T>(idx)};
tree.erase(p, {p[0], p[1], p[0], p[1]});
last = circList.insert(last, p);
queue.push_back(last);
Expand Down
3 changes: 3 additions & 0 deletions src/convex_hull.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,9 @@ convex_hull_indexes(const Eigen::Ref<const RowVectorsNx2> &points,
bool order_only = false)
{
const int N = points.rows();
if (N == 0) {
return Eigen::VectorXi(0);
}
Eigen::Vector2d p0(points(0, 0), points(0, 1));
for (int i = 1; i < N; ++i) {
if (points(i, 1) < p0[1] ||
Expand Down
142 changes: 142 additions & 0 deletions src/convex_hull2.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,142 @@
// https://github.com/mapbox/concaveman/pull/17/files

// https://github.com/mourner/robust-predicates/blob/main/src/orient2d.js
#pragma once

#include "orient2d.hpp"
#include <vector>
#include <array>

namespace cubao
{
namespace convex_hull2
{
// https://github.com/mourner/robust-predicates/blob/main/src/orient2d.js

inline double cross(const double *p1, const double *p2, const double *p3)
{
return robust_predicates::orient2d(p1[0], p1[1], //
p2[0], p2[1], //
p3[0], p3[1]);
}

// https://github.com/mapbox/concaveman/blob/master/index.js
// check if the edges (p1,q1) and (p2,q2) intersect
inline bool intersects(const double *p1, const double *q1, //
const double *p2, const double *q2)
{
return !(p1[0] == q2[0] && p1[1] == q2[1]) &&
!(q1[0] == p2[0] && q1[1] == p2[1]) &&
(cross(p1, q1, p2) > 0) != (cross(p1, q1, q2) > 0) &&
(cross(p2, q2, p1) > 0) != (cross(p2, q2, q1) > 0);
}

inline std::vector<std::array<double, 3>>
convexHull(std::vector<std::array<double, 3>> points)
{
if (points.empty()) {
return {};
}
std::sort(
points.begin(), points.end(),
[](const std::array<double, 3> &a, const std::array<double, 3> &b) {
return a[0] == b[0] ? a[1] < b[1] : a[0] < b[0];
});
points.erase(std::unique(points.begin(), points.end(),
[](const auto &a, const auto &b) {
return a[0] == b[0] && a[1] == b[1];
}),
points.end());

std::vector<std::array<double, 3>> lower;
for (size_t i = 0; i < points.size(); i++) {
while (lower.size() >= 2 &&
cross(lower[lower.size() - 2].data(),
lower[lower.size() - 1].data(), points[i].data()) <= 0) {
lower.pop_back();
}
lower.push_back(points[i]);
}

std::vector<std::array<double, 3>> upper;
for (int i = points.size() - 1; i >= 0; i--) {
while (upper.size() >= 2 &&
cross(upper[upper.size() - 2].data(),
upper[upper.size() - 1].data(), points[i].data()) <= 0) {
upper.pop_back();
}
upper.push_back(points[i]);
}

upper.pop_back();
lower.pop_back();
lower.insert(lower.end(), upper.begin(), upper.end());
return lower;
}

// https://www.npmjs.com/package/point-in-polygon?activeTab=code
inline bool pointInPolygon(const double *point, //
const double *polygon, size_t length, //
int start = -1, int end = -1)
{
double x = point[0], y = point[1];
bool inside = false;
if (start < 0)
start = 0;
if (end < 0)
end = length;
int len = end - start;
for (int i = 0, j = len - 1; i < len; j = i++) {
double xi = polygon[(i + start) * 2], yi = polygon[(i + start) * 2 + 1];
double xj = polygon[(j + start) * 2], yj = polygon[(j + start) * 2 + 1];
bool intersect = ((yi > y) != (yj > y)) &&
(x < (xj - xi) * (y - yi) / (yj - yi) + xi);
if (intersect)
inside = !inside;
}
return inside;
}

// speed up convex hull by filtering out points inside quadrilateral formed by 4
// extreme points
inline std::vector<std::array<double, 3>>
fastConvexHull(const std::vector<std::array<double, 3>> &points)
{
if (points.empty()) {
return {};
}

auto left = points[0];
auto top = points[0];
auto right = points[0];
auto bottom = points[0];

// find the leftmost, rightmost, topmost and bottommost points
for (const auto &p : points) {
if (p[0] < left[0])
left = p;
if (p[0] > right[0])
right = p;
if (p[1] < top[1])
top = p;
if (p[1] > bottom[1])
bottom = p;
}

// filter out points that are inside the resulting quadrilateral
std::vector<std::array<double, 2>> cull = {{left[0], left[1]},
{top[0], top[1]},
{right[0], right[1]},
{bottom[0], bottom[1]}};
std::vector<std::array<double, 3>> filtered = {left, top, right, bottom};
for (const auto &p : points) {
if (!pointInPolygon(&p[0], &cull[0][0], cull.size())) {
filtered.push_back(p);
}
}
// get convex hull around the filtered points
return convexHull(filtered);
}

} // namespace convex_hull2
} // namespace cubao
135 changes: 135 additions & 0 deletions src/macros.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,135 @@
#pragma once

#include <cmath>

namespace cubao
{
namespace robust_predicates
{
// https://github.com/mourner/robust-predicates/blob/312178014c173cfa5ed0a8cb5496baaf82a63663/compile.js#L6
#define Fast_Two_Sum(a, b, x, y) \
do { \
x = a + b; \
y = b - (x - a); \
} while (0)

#define Two_Sum(a, b, x, y) \
do { \
double bvirt; \
x = a + b; \
bvirt = x - a; \
y = a - (x - bvirt) + (b - bvirt); \
} while (0)

#define Two_Diff_Tail(a, b, x, y) \
do { \
double bvirt = a - x; \
y = a - (x + bvirt) + (bvirt - b); \
} while (0)

#define Two_Diff(a, b, x, y) \
do { \
x = a - b; \
Two_Diff_Tail(a, b, x, y); \
} while (0)

#define Split(a, ahi, alo) \
do { \
double c = splitter * a; \
ahi = c - (c - a); \
alo = a - ahi; \
} while (0)

#define Two_Product(a, b, x, y) \
do { \
double ahi, alo, bhi, blo; \
x = a * b; \
Split(a, ahi, alo); \
Split(b, bhi, blo); \
y = alo * blo - (x - ahi * bhi - alo * bhi - ahi * blo); \
} while (0)

#define Two_Product_Presplit(a, b, bhi, blo, x, y) \
do { \
double ahi, alo; \
x = a * b; \
Split(a, ahi, alo); \
y = alo * blo - (x - ahi * bhi - alo * bhi - ahi * blo); \
} while (0)

#define Square(a, x, y) \
do { \
double ahi, alo; \
x = a * a; \
Split(a, ahi, alo); \
y = alo * alo - (x - ahi * ahi - (ahi + ahi) * alo); \
} while (0)

#define Two_One_Sum(a1, a0, b, x2, x1, x0) \
do { \
double _i; \
Two_Sum(a0, b, _i, x0); \
Two_Sum(a1, _i, x2, x1); \
} while (0)

#define Two_One_Diff(a1, a0, b, x2, x1, x0) \
do { \
double _i; \
Two_Diff(a0, b, _i, x0); \
Two_Sum(a1, _i, x2, x1); \
} while (0)

#define Two_Two_Sum(a1, a0, b1, b0, x3, x2, x1, x0) \
do { \
double _j, _0; \
Two_One_Sum(a1, a0, b0, _j, _0, x0); \
Two_One_Sum(_j, _0, b1, x3, x2, x1); \
} while (0)

#define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) \
do { \
double _j, _0; \
Two_One_Diff(a1, a0, b0, _j, _0, x0); \
Two_One_Diff(_j, _0, b1, x3, x2, x1); \
} while (0)

#define Two_One_Product(a1, a0, b, D) \
do { \
double bhi, blo, _i, _j, _0, _k, u3; \
Split(b, bhi, blo); \
Two_Product_Presplit(a0, b, bhi, blo, _i, D[0]); \
Two_Product_Presplit(a1, b, bhi, blo, _j, _0); \
Two_Sum(_i, _0, _k, D[1]); \
Fast_Two_Sum(_j, _k, u3, D[2]); \
D[3] = u3; \
} while (0)

#define Cross_Product(a, b, c, d, D) \
do { \
double s1, s0, t1, t0, u3; \
Two_Product(a, d, s1, s0); \
Two_Product(c, b, t1, t0); \
Two_Two_Diff(s1, s0, t1, t0, u3, D[2], D[1], D[0]); \
D[3] = u3; \
} while (0)

#define Two_Product_Sum(a, b, c, d, D) \
do { \
double s1, s0, t1, t0, u3; \
Two_Product(a, b, s1, s0); \
Two_Product(c, d, t1, t0); \
Two_Two_Sum(s1, s0, t1, t0, u3, D[2], D[1], D[0]); \
D[3] = u3; \
} while (0)

#define Square_Sum(a, b, D) \
do { \
double s1, s0, t1, t0, u3; \
Square(a, s1, s0); \
Square(b, t1, t0); \
Two_Two_Sum(s1, s0, t1, t0, u3, D[2], D[1], D[0]); \
D[3] = u3; \
} while (0)

} // namespace robust_predicates
} // namespace cubao
Loading

0 comments on commit 5cc777b

Please sign in to comment.