Aleph-w 3.0
A C++ Library for Data Structures and Algorithms
Loading...
Searching...
No Matches
geom_algorithms_test_cdt.cc
Go to the documentation of this file.
2#include <cmath>
3
4// ============================================================================
5// Helper: check that an indexed edge (u,v) appears in a CDT result
6// ============================================================================
7static bool cdt_has_edge(
9 const size_t u, const size_t v)
10{
11 for (size_t t = 0; t < r.triangles.size(); ++t)
12 {
13 const auto & tri = r.triangles(t);
14 size_t vs[3] = {tri.i, tri.j, tri.k};
15 for (int e = 0; e < 3; ++e)
16 {
17 const size_t a = vs[e], b = vs[(e + 1) % 3];
18 if ((a == u and b == v) or (a == v and b == u))
19 return true;
20 }
21 }
22 return false;
23}
24
25// Helper: find point index in sorted result sites
26static size_t find_site(const Array<Point> & sites, const Point & p)
27{
28 for (size_t i = 0; i < sites.size(); ++i)
29 if (sites(i) == p)
30 return i;
31 return ~static_cast<size_t>(0);
32}
33
34// Helper: check the constrained edge list contains an edge between two points
37 const Point & p, const Point & q)
38{
39 const size_t u = find_site(r.sites, p);
40 const size_t v = find_site(r.sites, q);
41 if (u == ~static_cast<size_t>(0) or v == ~static_cast<size_t>(0))
42 return false;
43
44 for (size_t i = 0; i < r.constrained_edges.size(); ++i)
45 {
46 const auto & e = r.constrained_edges(i);
47 if ((e.u == u and e.v == v) or (e.u == v and e.v == u))
48 return true;
49 }
50 return false;
51}
52
53// Helper: check Delaunay property for non-constrained edges
54// For each triangle, check that no non-adjacent vertex violates circumcircle
55// (only approximate check — checks that mesh triangles don't violate each other)
58{
59 // Build a set of constrained edge pairs for fast lookup
60 struct EdgeKey { size_t u, v; };
62 for (size_t i = 0; i < r.constrained_edges.size(); ++i)
63 {
64 const auto & e = r.constrained_edges(i);
65 const size_t a = e.u < e.v ? e.u : e.v;
66 const size_t b = e.u < e.v ? e.v : e.u;
67 con_edges.append(EdgeKey{a, b});
68 }
69
70 auto is_constrained = [&](size_t a, size_t b) -> bool
71 {
72 if (a > b) { const size_t t = a; a = b; b = t; }
73 for (size_t i = 0; i < con_edges.size(); ++i)
74 if (con_edges(i).u == a and con_edges(i).v == b)
75 return true;
76 return false;
77 };
78
79 // For each pair of triangles sharing a non-constrained edge,
80 // check that neither opposite vertex is inside the other's circumcircle.
81 for (size_t t1 = 0; t1 < r.triangles.size(); ++t1)
82 {
83 const auto & tri1 = r.triangles(t1);
84 size_t vs1[3] = {tri1.i, tri1.j, tri1.k};
85
86 for (int e = 0; e < 3; ++e)
87 {
88 const size_t ea = vs1[e], eb = vs1[(e + 1) % 3];
89 if (is_constrained(ea, eb))
90 continue;
91
92 // Find the other triangle sharing this edge
93 for (size_t t2 = t1 + 1; t2 < r.triangles.size(); ++t2)
94 {
95 const auto & tri2 = r.triangles(t2);
96 size_t vs2[3] = {tri2.i, tri2.j, tri2.k};
97
98 bool shares = false;
99 size_t opp = ~static_cast<size_t>(0);
100 for (int f = 0; f < 3; ++f)
101 {
102 const size_t fa = vs2[f], fb = vs2[(f + 1) % 3];
103 if ((fa == ea and fb == eb) or (fa == eb and fb == ea))
104 {
105 shares = true;
106 opp = vs2[(f + 2) % 3];
107 break;
108 }
109 }
110
111 if (not shares)
112 continue;
113
114 // Check that opp is not inside circumcircle of tri1
115 const Orientation o = orientation(r.sites(tri1.i),
116 r.sites(tri1.j),
117 r.sites(tri1.k));
118 if (o == Orientation::COLLINEAR)
119 continue;
120
122 r.sites(tri1.i), r.sites(tri1.j),
123 r.sites(tri1.k), r.sites(opp));
124
125 if (o == Orientation::CCW and det > 0)
126 return false; // violation
127 if (o == Orientation::CW and det < 0)
128 return false;
129 }
130 }
131 }
132
133 return true;
134}
135
136// ============================================================================
137// Basic tests
138// ============================================================================
139
141{
143 pts.append(Point(0, 0));
144 pts.append(Point(4, 0));
145 pts.append(Point(2, 3));
146
148
150 auto result = cdt(pts, constraints);
151
152 EXPECT_EQ(result.sites.size(), 3u);
153 EXPECT_EQ(result.triangles.size(), 1u);
154 EXPECT_TRUE(result.constrained_edges.is_empty());
155}
156
158{
159 Point a(0, 0), b(4, 0), c(4, 4), d(0, 4);
160
162 pts.append(a); pts.append(b); pts.append(c); pts.append(d);
163
164 // Force diagonal a-c
167
169 auto result = cdt(pts, constraints);
170
171 EXPECT_EQ(result.sites.size(), 4u);
172 EXPECT_EQ(result.triangles.size(), 2u);
173
174 // The constrained diagonal must be present
175 EXPECT_TRUE(has_constrained_edge(result, a, c));
176
177 // The edge a-c should be in the triangulation
178 const size_t ia = find_site(result.sites, a);
179 const size_t ic = find_site(result.sites, c);
180 EXPECT_TRUE(cdt_has_edge(result, ia, ic));
181}
182
184{
185 // Triangle — any edge is already in DT
186 Point a(0, 0), b(6, 0), c(3, 5);
187
189 pts.append(a); pts.append(b); pts.append(c);
190
192 constraints.append(Segment(a, b)); // already a DT edge
193
195 auto result = cdt(pts, constraints);
196
197 EXPECT_EQ(result.triangles.size(), 1u);
198 EXPECT_TRUE(has_constrained_edge(result, a, b));
199}
200
202{
203 // Diamond shape — force the non-Delaunay diagonal
204 Point a(0, 0), b(2, 1), c(4, 0), d(2, -1);
205
207 pts.append(a); pts.append(b); pts.append(c); pts.append(d);
208
209 // The Delaunay diagonal is likely b-d. Force a-c instead.
212
214 auto result = cdt(pts, constraints);
215
216 EXPECT_EQ(result.sites.size(), 4u);
217 EXPECT_EQ(result.triangles.size(), 2u);
218 EXPECT_TRUE(has_constrained_edge(result, a, c));
219
220 const size_t ia = find_site(result.sites, a);
221 const size_t ic = find_site(result.sites, c);
222 EXPECT_TRUE(cdt_has_edge(result, ia, ic));
223}
224
226{
227 Point a(0, 0), b(6, 0), c(6, 6), d(0, 6), e(3, 3);
228
230 pts.append(a); pts.append(b); pts.append(c);
231 pts.append(d); pts.append(e);
232
234 constraints.append(Segment(a, c)); // diagonal
235 constraints.append(Segment(b, d)); // other diagonal
236
238 auto result = cdt(pts, constraints);
239
240 EXPECT_GE(result.sites.size(), 5u);
241 EXPECT_GE(result.triangles.size(), 4u);
242
243 // Both diagonals cross at e(3,3), so constraints are split.
244 // a-c becomes (a,e) + (e,c); b-d becomes (b,e) + (e,d).
245 EXPECT_TRUE(has_constrained_edge(result, a, e));
246 EXPECT_TRUE(has_constrained_edge(result, e, c));
247 EXPECT_TRUE(has_constrained_edge(result, b, e));
248 EXPECT_TRUE(has_constrained_edge(result, e, d));
249}
250
252{
253 // A zig-zag chain as constraints
254 Point p0(0, 0), p1(2, 3), p2(4, 0), p3(6, 3), p4(8, 0);
255
257 pts.append(p0); pts.append(p1); pts.append(p2);
258 pts.append(p3); pts.append(p4);
259
261 constraints.append(Segment(p0, p1));
262 constraints.append(Segment(p1, p2));
263 constraints.append(Segment(p2, p3));
264 constraints.append(Segment(p3, p4));
265
267 auto result = cdt(pts, constraints);
268
269 EXPECT_EQ(result.sites.size(), 5u);
270 EXPECT_GE(result.triangles.size(), 3u);
271
272 EXPECT_TRUE(has_constrained_edge(result, p0, p1));
273 EXPECT_TRUE(has_constrained_edge(result, p1, p2));
274 EXPECT_TRUE(has_constrained_edge(result, p2, p3));
275 EXPECT_TRUE(has_constrained_edge(result, p3, p4));
276}
277
278// ============================================================================
279// Edge cases
280// ============================================================================
281
283{
286
288 auto result = cdt(pts, constraints);
289
290 EXPECT_TRUE(result.triangles.is_empty());
291}
292
294{
296 pts.append(Point(0, 0));
297 pts.append(Point(1, 1));
298
300
302 auto result = cdt(pts, constraints);
303
304 EXPECT_TRUE(result.triangles.is_empty());
305}
306
308{
310 pts.append(Point(0, 0));
311 pts.append(Point(1, 0));
312 pts.append(Point(2, 0));
313
315
317 auto result = cdt(pts, constraints);
318
319 EXPECT_TRUE(result.triangles.is_empty());
320}
321
323{
324 // Points come only from constraint segments
325 Point a(0, 0), b(4, 0), c(2, 3);
326
327 DynList<Point> pts; // empty point list
328
331 constraints.append(Segment(b, c));
332 constraints.append(Segment(c, a));
333
335 auto result = cdt(pts, constraints);
336
337 EXPECT_EQ(result.sites.size(), 3u);
338 EXPECT_EQ(result.triangles.size(), 1u);
339
340 EXPECT_TRUE(has_constrained_edge(result, a, b));
341 EXPECT_TRUE(has_constrained_edge(result, b, c));
342 EXPECT_TRUE(has_constrained_edge(result, c, a));
343}
344
346{
347 // A constraint from a to c, with b lying on the segment
348 Point a(0, 0), b(2, 0), c(4, 0), d(2, 3);
349
351 pts.append(a); pts.append(b); pts.append(c); pts.append(d);
352
354 constraints.append(Segment(a, c)); // b lies on this
355
357 auto result = cdt(pts, constraints);
358
359 EXPECT_EQ(result.sites.size(), 4u);
360 EXPECT_GE(result.triangles.size(), 2u);
361
362 // The constraint should be split into a-b and b-c
363 const size_t ia = find_site(result.sites, a);
364 const size_t ib = find_site(result.sites, b);
365 const size_t ic = find_site(result.sites, c);
366 EXPECT_TRUE(cdt_has_edge(result, ia, ib));
367 EXPECT_TRUE(cdt_has_edge(result, ib, ic));
368}
369
371{
372 Point a(0, 0), b(4, 0), c(2, 3);
373
375 pts.append(a); pts.append(b); pts.append(c);
376 pts.append(a); pts.append(b); // duplicates
377
379
381 auto result = cdt(pts, constraints);
382
383 EXPECT_EQ(result.sites.size(), 3u);
384 EXPECT_EQ(result.triangles.size(), 1u);
385}
386
388{
389 Point a(0, 0), b(4, 0), c(2, 3);
390
392 pts.append(a); pts.append(b); pts.append(c);
393
396 constraints.append(Segment(a, b)); // duplicate
397 constraints.append(Segment(b, a)); // reverse duplicate
398
400 auto result = cdt(pts, constraints);
401
402 EXPECT_EQ(result.sites.size(), 3u);
403 EXPECT_EQ(result.triangles.size(), 1u);
404 EXPECT_TRUE(has_constrained_edge(result, a, b));
405}
406
407// ============================================================================
408// Property tests
409// ============================================================================
410
412{
413 // Larger set with multiple constraints
414 Point p0(0, 0), p1(5, 0), p2(5, 5), p3(0, 5);
415 Point p4(1, 1), p5(4, 1), p6(4, 4), p7(1, 4);
416
418 pts.append(p0); pts.append(p1); pts.append(p2); pts.append(p3);
419 pts.append(p4); pts.append(p5); pts.append(p6); pts.append(p7);
420
421 // Inner square edges as constraints
424 constraints.append(Segment(p5, p6));
425 constraints.append(Segment(p6, p7));
426 constraints.append(Segment(p7, p4));
427
429 auto result = cdt(pts, constraints);
430
431 EXPECT_GE(result.sites.size(), 8u);
432
433 // All constraint edges must be in the triangulation
434 for (size_t i = 0; i < result.constrained_edges.size(); ++i)
435 {
436 const auto & e = result.constrained_edges(i);
437 EXPECT_TRUE(cdt_has_edge(result, e.u, e.v))
438 << "Constrained edge (" << e.u << "," << e.v << ") not in triangulation";
439 }
440
441 EXPECT_TRUE(has_constrained_edge(result, p4, p5));
444 EXPECT_TRUE(has_constrained_edge(result, p7, p4));
445}
446
448{
449 Point p0(0, 0), p1(6, 0), p2(6, 6), p3(0, 6), p4(3, 3);
450
452 pts.append(p0); pts.append(p1); pts.append(p2);
453 pts.append(p3); pts.append(p4);
454
455 // One constraint
457 constraints.append(Segment(p0, p2));
458
460 auto result = cdt(pts, constraints);
461
463 << "Delaunay property violated for non-constrained edges";
464}
465
467{
468 Point a(0, 0), b(4, 0), c(4, 4), d(0, 4), e(2, 2);
469
471 pts.append(a); pts.append(b); pts.append(c);
472 pts.append(d); pts.append(e);
473
476
478 auto result = cdt(pts, constraints);
479
480 EXPECT_NE(find_site(result.sites, a), ~static_cast<size_t>(0));
481 EXPECT_NE(find_site(result.sites, b), ~static_cast<size_t>(0));
482 EXPECT_NE(find_site(result.sites, c), ~static_cast<size_t>(0));
483 EXPECT_NE(find_site(result.sites, d), ~static_cast<size_t>(0));
484 EXPECT_NE(find_site(result.sites, e), ~static_cast<size_t>(0));
485}
486
488{
489 Point p0(0, 0), p1(6, 0), p2(6, 4), p3(0, 4), p4(3, 2);
490
492 pts.append(p0); pts.append(p1); pts.append(p2);
493 pts.append(p3); pts.append(p4);
494
496 constraints.append(Segment(p0, p2));
497
499 auto result = cdt(pts, constraints);
500
501 // Check all triangles that are non-degenerate are CCW
502 for (size_t t = 0; t < result.triangles.size(); ++t)
503 {
504 const auto & tri = result.triangles(t);
505 const Orientation o = orientation(result.sites(tri.i),
506 result.sites(tri.j),
507 result.sites(tri.k));
508 EXPECT_NE(o, Orientation::COLLINEAR) << "Degenerate triangle " << t;
509 }
510}
511
512// ============================================================================
513// Stress tests
514// ============================================================================
515
517{
518 // 4x4 grid of points with diagonal constraints
520 for (int y = 0; y < 4; ++y)
521 for (int x = 0; x < 4; ++x)
522 pts.append(Point(x, y));
523
524 // Diagonal constraints across each cell
526 for (int y = 0; y < 3; ++y)
527 for (int x = 0; x < 3; ++x)
528 constraints.append(Segment(Point(x, y), Point(x + 1, y + 1)));
529
531 auto result = cdt(pts, constraints);
532
533 EXPECT_EQ(result.sites.size(), 16u);
534 EXPECT_GE(result.triangles.size(), 18u); // at least 2 per cell
535
536 // All diagonal constraints present
537 for (int y = 0; y < 3; ++y)
538 for (int x = 0; x < 3; ++x)
540 Point(x + 1, y + 1)));
541}
542
544{
545 // Pentagon boundary as constraints with interior point
546 Point p0(2, 0), p1(4, 1.5), p2(3, 4), p3(1, 4), p4(0, 1.5);
547 Point center(2, 2);
548
550 pts.append(p0); pts.append(p1); pts.append(p2);
551 pts.append(p3); pts.append(p4); pts.append(center);
552
554 constraints.append(Segment(p0, p1));
555 constraints.append(Segment(p1, p2));
556 constraints.append(Segment(p2, p3));
557 constraints.append(Segment(p3, p4));
558 constraints.append(Segment(p4, p0));
559
561 auto result = cdt(pts, constraints);
562
563 EXPECT_EQ(result.sites.size(), 6u);
564 EXPECT_GE(result.triangles.size(), 4u);
565
566 EXPECT_TRUE(has_constrained_edge(result, p0, p1));
567 EXPECT_TRUE(has_constrained_edge(result, p1, p2));
568 EXPECT_TRUE(has_constrained_edge(result, p2, p3));
569 EXPECT_TRUE(has_constrained_edge(result, p3, p4));
570 EXPECT_TRUE(has_constrained_edge(result, p4, p0));
571}
572
574{
575 Point a(0, 0), b(4, 0), c(2, 3);
576
578 pts.append(a); pts.append(b); pts.append(c);
579
582
584 auto result = cdt(pts, constraints);
585
587 EXPECT_EQ(triangles.size(), result.triangles.size());
588}
589
591{
593 auto result = cdt(
594 {Point(0, 0), Point(4, 0), Point(4, 4), Point(0, 4)},
595 {Segment(Point(0, 0), Point(4, 4))});
596
597 EXPECT_EQ(result.sites.size(), 4u);
598 EXPECT_EQ(result.triangles.size(), 2u);
599 EXPECT_TRUE(has_constrained_edge(result, Point(0, 0), Point(4, 4)));
600}
601
603{
604 // 8 points with 2 constraints — verify Delaunay property
606 pts.append(Point(0, 0));
607 pts.append(Point(5, 0));
608 pts.append(Point(10, 0));
609 pts.append(Point(10, 5));
610 pts.append(Point(10, 10));
611 pts.append(Point(5, 10));
612 pts.append(Point(0, 10));
613 pts.append(Point(0, 5));
614
616 constraints.append(Segment(Point(0, 0), Point(10, 10)));
617
619 auto result = cdt(pts, constraints);
620
621 EXPECT_GE(result.triangles.size(), 6u);
622 EXPECT_TRUE(has_constrained_edge(result, Point(0, 0), Point(10, 10)));
624}
Simple dynamic array with automatic resizing and functional operations.
Definition tpl_array.H:139
constexpr size_t size() const noexcept
Return the number of elements stored in the stack.
Definition tpl_array.H:351
T & append(const T &data)
Append a copy of data
Definition tpl_array.H:245
Constrained Delaunay Triangulation via Sloan's flip-based method.
static DynList< Triangle > as_triangles(const Result &result)
Convert indexed triangulation to geometric triangles.
Dynamic singly linked list with functional programming support.
Definition htlist.H:1155
T & append(const T &item)
Definition htlist.H:1271
size_t size() const noexcept
Count the number of elements of the list.
Definition htlist.H:1065
Represents a point with rectangular coordinates in a 2D plane.
Definition point.H:229
Represents a line segment between two points.
Definition point.H:827
TEST_F(GeomAlgorithmsTest, CDT_TriangleNoConstraints)
static bool cdt_has_edge(const ConstrainedDelaunayTriangulation::Result &r, const size_t u, const size_t v)
static bool has_constrained_edge(const ConstrainedDelaunayTriangulation::Result &r, const Point &p, const Point &q)
static bool check_delaunay_for_non_constrained(const ConstrainedDelaunayTriangulation::Result &r)
static size_t find_site(const Array< Point > &sites, const Point &p)
static mpfr_t y
Definition mpfr_mul_d.c:3
and
Check uniqueness with explicit hash + equality functors.
Orientation
Classification of three-point orientation.
Definition point.H:2812
Geom_Number in_circle_determinant(const Point &a, const Point &b, const Point &c, const Point &p)
Return the exact in-circle determinant for (a,b,c,p).
Definition point.H:2835
Divide_Conquer_DP_Result< Cost > divide_and_conquer_partition_dp(const size_t groups, const size_t n, Transition_Cost_Fn transition_cost, const Cost inf=dp_optimization_detail::default_inf< Cost >())
Optimize partition DP using divide-and-conquer optimization.
Orientation orientation(const Point &a, const Point &b, const Point &c)
Return the orientation of the triple (a, b, c).
Definition point.H:2821
gsl_rng * r