Skip to content

Commit f1e2374

Browse files
unixpickleunixpickle
committed
convex hull
1 parent 96ec17e commit f1e2374

2 files changed

Lines changed: 318 additions & 0 deletions

File tree

‎model2d/hull_mesh.go‎

Lines changed: 111 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
1+
package model2d
2+
3+
import (
4+
"math"
5+
"sort"
6+
)
7+
8+
// ConvexHullMesh computes the convex hull of points and returns it as a mesh.
9+
//
10+
// The hull edges are ordered clockwise (assuming the y-axis points up), so
11+
// segment normals face outward. Colinear points on the hull are removed.
12+
func ConvexHullMesh(points []Coord) *Mesh {
13+
unique := make([]Coord, 0, len(points))
14+
seen := NewCoordToNumber[int]()
15+
for _, p := range points {
16+
if _, ok := seen.Load(p); !ok {
17+
seen.Store(p, 1)
18+
unique = append(unique, p)
19+
}
20+
}
21+
22+
if len(unique) < 2 {
23+
return NewMesh()
24+
}
25+
26+
pivotIdx := 0
27+
for i := 1; i < len(unique); i++ {
28+
p, q := unique[i], unique[pivotIdx]
29+
if p.Y < q.Y || (p.Y == q.Y && p.X < q.X) {
30+
pivotIdx = i
31+
}
32+
}
33+
unique[0], unique[pivotIdx] = unique[pivotIdx], unique[0]
34+
pivot := unique[0]
35+
36+
scale := maxCoordAbs(unique)
37+
if scale == 0 {
38+
scale = 1
39+
}
40+
eps := scale * scale * 1e-12
41+
42+
others := append([]Coord{}, unique[1:]...)
43+
sort.Slice(others, func(i, j int) bool {
44+
a, b := others[i], others[j]
45+
cross := cross2(pivot, a, b)
46+
if math.Abs(cross) > eps {
47+
return cross > 0
48+
}
49+
da := pivot.SquaredDist(a)
50+
db := pivot.SquaredDist(b)
51+
if da != db {
52+
return da < db
53+
}
54+
if a.X != b.X {
55+
return a.X < b.X
56+
}
57+
return a.Y < b.Y
58+
})
59+
60+
hull := []Coord{pivot}
61+
for _, p := range others {
62+
for len(hull) >= 2 {
63+
c := cross2(hull[len(hull)-2], hull[len(hull)-1], p)
64+
if c > eps {
65+
break
66+
}
67+
hull = hull[:len(hull)-1]
68+
}
69+
hull = append(hull, p)
70+
}
71+
72+
if len(hull) < 2 {
73+
return NewMesh()
74+
}
75+
if len(hull) == 2 {
76+
m := NewMesh()
77+
m.Add(&Segment{hull[0], hull[1]})
78+
return m
79+
}
80+
81+
for i, j := 0, len(hull)-1; i < j; i, j = i+1, j-1 {
82+
hull[i], hull[j] = hull[j], hull[i]
83+
}
84+
85+
m := NewMesh()
86+
for i := 0; i < len(hull); i++ {
87+
p1 := hull[i]
88+
p2 := hull[(i+1)%len(hull)]
89+
m.Add(&Segment{p1, p2})
90+
}
91+
return m
92+
}
93+
94+
func cross2(a, b, c Coord) float64 {
95+
return (b.X-a.X)*(c.Y-a.Y) - (b.Y-a.Y)*(c.X-a.X)
96+
}
97+
98+
func maxCoordAbs(points []Coord) float64 {
99+
var maxAbs float64
100+
for _, p := range points {
101+
ax := math.Abs(p.X)
102+
ay := math.Abs(p.Y)
103+
if ax > maxAbs {
104+
maxAbs = ax
105+
}
106+
if ay > maxAbs {
107+
maxAbs = ay
108+
}
109+
}
110+
return maxAbs
111+
}

‎model2d/hull_mesh_test.go‎

Lines changed: 207 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,207 @@
1+
package model2d
2+
3+
import (
4+
"math"
5+
"math/rand"
6+
"testing"
7+
)
8+
9+
func TestConvexHullMeshBasic(t *testing.T) {
10+
cases := []struct {
11+
name string
12+
points []Coord
13+
want []Segment
14+
}{
15+
{
16+
name: "single",
17+
points: []Coord{XY(1, 2)},
18+
want: nil,
19+
},
20+
{
21+
name: "two_points",
22+
points: []Coord{XY(0, 0), XY(2, 0)},
23+
want: []Segment{{XY(0, 0), XY(2, 0)}},
24+
},
25+
{
26+
name: "triangle",
27+
points: []Coord{XY(0, 0), XY(1, 0), XY(0, 1)},
28+
want: []Segment{
29+
{XY(0, 0), XY(1, 0)},
30+
{XY(1, 0), XY(0, 1)},
31+
{XY(0, 1), XY(0, 0)},
32+
},
33+
},
34+
{
35+
name: "square_with_interior",
36+
points: []Coord{XY(0, 0), XY(1, 0), XY(1, 1), XY(0, 1), XY(0.5, 0.5)},
37+
want: []Segment{
38+
{XY(0, 0), XY(1, 0)},
39+
{XY(1, 0), XY(1, 1)},
40+
{XY(1, 1), XY(0, 1)},
41+
{XY(0, 1), XY(0, 0)},
42+
},
43+
},
44+
{
45+
name: "colinear",
46+
points: []Coord{XY(0, 0), XY(1, 0), XY(2, 0), XY(0.5, 0), XY(2, 0)},
47+
want: []Segment{{XY(0, 0), XY(2, 0)}},
48+
},
49+
}
50+
51+
for _, tc := range cases {
52+
t.Run(tc.name, func(t *testing.T) {
53+
got := ConvexHullMesh(tc.points)
54+
want := NewMesh()
55+
for i := range tc.want {
56+
seg := tc.want[i]
57+
want.Add(&seg)
58+
}
59+
if !meshEdgeSetEqual(got, want) {
60+
t.Fatalf("edge sets differ: got=%v want=%v",
61+
meshEdgeSet(got), meshEdgeSet(want))
62+
}
63+
})
64+
}
65+
}
66+
67+
func TestConvexHullMeshRandomMatchesNaive(t *testing.T) {
68+
rng := rand.New(rand.NewSource(1337))
69+
for i := 0; i < 5; i++ {
70+
points := make([]Coord, 10)
71+
for j := range points {
72+
points[j] = XY(rng.Float64()*2-1, rng.Float64()*2-1)
73+
}
74+
got := ConvexHullMesh(points)
75+
want := naiveConvexHullMesh(points)
76+
MustValidateMesh(t, got, true)
77+
MustValidateMesh(t, want, true)
78+
if !meshEdgeSetEqual(got, want) {
79+
t.Fatalf("random draw %d mismatch: got=%v want=%v",
80+
i, meshEdgeSet(got), meshEdgeSet(want))
81+
}
82+
}
83+
}
84+
85+
func naiveConvexHullMesh(points []Coord) *Mesh {
86+
unique := make([]Coord, 0, len(points))
87+
seen := NewCoordToNumber[int]()
88+
for _, p := range points {
89+
if _, ok := seen.Load(p); !ok {
90+
seen.Store(p, 1)
91+
unique = append(unique, p)
92+
}
93+
}
94+
if len(unique) < 2 {
95+
return NewMesh()
96+
}
97+
98+
start := unique[0]
99+
for _, p := range unique[1:] {
100+
if p.Y < start.Y || (p.Y == start.Y && p.X < start.X) {
101+
start = p
102+
}
103+
}
104+
105+
scale := maxCoordAbs(unique)
106+
if scale == 0 {
107+
scale = 1
108+
}
109+
eps := scale * scale * 1e-12
110+
111+
hull := []Coord{start}
112+
curr := start
113+
for {
114+
var next Coord
115+
found := false
116+
for _, p := range unique {
117+
if p == curr {
118+
continue
119+
}
120+
if !found {
121+
next = p
122+
found = true
123+
continue
124+
}
125+
cross := cross2(curr, next, p)
126+
if cross < -eps {
127+
next = p
128+
} else if math.Abs(cross) <= eps {
129+
if curr.SquaredDist(p) > curr.SquaredDist(next) {
130+
next = p
131+
}
132+
}
133+
}
134+
if !found || next == start {
135+
break
136+
}
137+
hull = append(hull, next)
138+
curr = next
139+
if len(hull) > len(unique)+1 {
140+
break
141+
}
142+
}
143+
144+
if len(hull) < 2 {
145+
return NewMesh()
146+
}
147+
if len(hull) == 2 {
148+
m := NewMesh()
149+
m.Add(&Segment{hull[0], hull[1]})
150+
return m
151+
}
152+
if !isClockwise(hull) {
153+
for i, j := 0, len(hull)-1; i < j; i, j = i+1, j-1 {
154+
hull[i], hull[j] = hull[j], hull[i]
155+
}
156+
}
157+
158+
m := NewMesh()
159+
for i := 0; i < len(hull); i++ {
160+
p1 := hull[i]
161+
p2 := hull[(i+1)%len(hull)]
162+
m.Add(&Segment{p1, p2})
163+
}
164+
return m
165+
}
166+
167+
func meshEdgeSet(m *Mesh) map[Segment]struct{} {
168+
res := map[Segment]struct{}{}
169+
for _, seg := range m.SegmentSlice() {
170+
a, b := seg[0], seg[1]
171+
if lessCoord(b, a) {
172+
a, b = b, a
173+
}
174+
res[Segment{a, b}] = struct{}{}
175+
}
176+
return res
177+
}
178+
179+
func meshEdgeSetEqual(m1, m2 *Mesh) bool {
180+
s1 := meshEdgeSet(m1)
181+
s2 := meshEdgeSet(m2)
182+
if len(s1) != len(s2) {
183+
return false
184+
}
185+
for k := range s1 {
186+
if _, ok := s2[k]; !ok {
187+
return false
188+
}
189+
}
190+
return true
191+
}
192+
193+
func lessCoord(a, b Coord) bool {
194+
if a.X != b.X {
195+
return a.X < b.X
196+
}
197+
return a.Y < b.Y
198+
}
199+
200+
func isClockwise(poly []Coord) bool {
201+
var area2 float64
202+
for i, p := range poly {
203+
q := poly[(i+1)%len(poly)]
204+
area2 += p.X*q.Y - q.X*p.Y
205+
}
206+
return area2 < 0
207+
}

0 commit comments

Comments
 (0)