Aleph-w 3.0
A C++ Library for Data Structures and Algorithms
Loading...
Searching...
No Matches
geom_algorithms_test_mec.cc
Go to the documentation of this file.
1/*
2 Aleph_w
3
4 Data structures & Algorithms
5 version 2.0.0b
6 https://github.com/lrleon/Aleph-w
7
8 This file is part of Aleph-w library
9
10 Copyright (c) 2002-2026 Leandro Rabindranath Leon
11
12 Permission is hereby granted, free of charge, to any person obtaining a copy
13 of this software and associated documentation files (the "Software"), to deal
14 in the Software without restriction, including without limitation the rights
15 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
16 copies of the Software, and to permit persons to whom the Software is
17 furnished to do so, subject to the following conditions:
18
19 The above copyright notice and this permission notice shall be included in all
20 copies or substantial portions of the Software.
21
22 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
23 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
25 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
27 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
28 SOFTWARE.
29*/
30
32#include <cmath>
33#include <random>
34
35using namespace Aleph;
36
37// ---------- Basic tests ----------
38
40{
42 pts.append(Point(3, 7));
43
45 auto c = mec(pts);
46
47 EXPECT_EQ(c.center, Point(3, 7));
48 EXPECT_EQ(c.radius_squared, 0);
49}
50
52{
54 pts.append(Point(0, 0));
55 pts.append(Point(4, 0));
56
58 auto c = mec(pts);
59
60 EXPECT_EQ(c.center, Point(2, 0));
61 EXPECT_EQ(c.radius_squared, 4);
62 EXPECT_TRUE(c.contains(Point(0, 0)));
63 EXPECT_TRUE(c.contains(Point(4, 0)));
64}
65
67{
68 // Equilateral triangle with vertices at (0,0), (4,0), (2, 2*sqrt(3))
69 // Circumradius = 2*sqrt(3)/sqrt(3)*2/sqrt(3) ... let's just check containment
71 pts.append(Point(0, 0));
72 pts.append(Point(4, 0));
74 pts.append(Point(2, h));
75
77 auto c = mec(pts);
78
79 // All points must be contained
80 EXPECT_TRUE(c.contains(Point(0, 0)));
81 EXPECT_TRUE(c.contains(Point(4, 0)));
82 EXPECT_TRUE(c.contains(Point(2, h)));
83
84 // Radius should equal circumradius = 4/sqrt(3)
86 const Geom_Number tol = Geom_Number(1, 1000000);
87 const Geom_Number actual_r = c.radius();
90 EXPECT_TRUE(diff < tol)
91 << "radius " << actual_r << " != expected " << expected_r;
92}
93
95{
96 // Right triangle: hypotenuse is diameter of circumscribed circle
98 pts.append(Point(0, 0));
99 pts.append(Point(6, 0));
100 pts.append(Point(0, 8));
101
103 auto c = mec(pts);
104
105 // Center should be midpoint of hypotenuse = (3, 4)
106 EXPECT_EQ(c.center, Point(3, 4));
107 // Radius^2 = 9 + 16 = 25
108 EXPECT_EQ(c.radius_squared, 25);
109
110 EXPECT_TRUE(c.contains(Point(0, 0)));
111 EXPECT_TRUE(c.contains(Point(6, 0)));
112 EXPECT_TRUE(c.contains(Point(0, 8)));
113}
114
116{
118 pts.append(Point(0, 0));
119 pts.append(Point(4, 0));
120 pts.append(Point(4, 4));
121 pts.append(Point(0, 4));
122
124 auto c = mec(pts);
125
126 // Center at (2, 2), radius^2 = 8
127 EXPECT_EQ(c.center, Point(2, 2));
128 EXPECT_EQ(c.radius_squared, 8);
129
130 for (DynList<Point>::Iterator it(pts); it.has_curr(); it.next_ne())
131 EXPECT_TRUE(c.contains(it.get_curr()));
132}
133
135{
136 // Regular pentagon centered at origin, circumradius = 5
138 const double r = 5.0;
139 for (int i = 0; i < 5; ++i)
140 {
141 const double angle = 2.0 * M_PI * i / 5.0;
142 pts.append(Point(Geom_Number(r * std::cos(angle)),
143 Geom_Number(r * std::sin(angle))));
144 }
145
147 auto c = mec(pts);
148
149 // All points should be contained
150 for (DynList<Point>::Iterator it(pts); it.has_curr(); it.next_ne())
151 EXPECT_TRUE(c.contains(it.get_curr()));
152
153 // Radius should be close to 5
154 const Geom_Number tol = Geom_Number(1, 100000);
155 const Geom_Number diff = c.radius() > Geom_Number(5) ?
156 c.radius() - Geom_Number(5) : Geom_Number(5) - c.radius();
157 EXPECT_TRUE(diff < tol)
158 << "radius " << c.radius() << " expected ~5";
159}
160
161// ---------- Edge cases ----------
162
169
171{
173 for (int i = 0; i < 10; ++i)
174 pts.append(Point(5, 5));
175
177 auto c = mec(pts);
178
179 EXPECT_EQ(c.center, Point(5, 5));
180 EXPECT_EQ(c.radius_squared, 0);
181}
182
184{
186 pts.append(Point(0, 0));
187 pts.append(Point(2, 0));
188 pts.append(Point(4, 0));
189 pts.append(Point(6, 0));
190 pts.append(Point(10, 0));
191
193 auto c = mec(pts);
194
195 // Diameter from (0,0) to (10,0): center = (5,0), r^2 = 25
196 EXPECT_EQ(c.center, Point(5, 0));
197 EXPECT_EQ(c.radius_squared, 25);
198
199 for (DynList<Point>::Iterator it(pts); it.has_curr(); it.next_ne())
200 EXPECT_TRUE(c.contains(it.get_curr()));
201}
202
204{
206 pts.append(Point(0, 0));
207 pts.append(Point(4, 0));
208 pts.append(Point(4, 0));
209 pts.append(Point(0, 0));
210 pts.append(Point(2, 3));
211
213 auto c = mec(pts);
214
215 for (DynList<Point>::Iterator it(pts); it.has_curr(); it.next_ne())
216 EXPECT_TRUE(c.contains(it.get_curr()));
217}
218
220{
221 // Triangle with an interior point
223 pts.append(Point(0, 0));
224 pts.append(Point(10, 0));
225 pts.append(Point(5, 10));
226 pts.append(Point(5, 3)); // interior
227
229 auto c = mec(pts);
230
231 for (DynList<Point>::Iterator it(pts); it.has_curr(); it.next_ne())
232 EXPECT_TRUE(c.contains(it.get_curr()));
233
234 // Interior point should not affect the circle (boundary should have 3 pts)
237 pts_no_interior.append(Point(10, 0));
238 pts_no_interior.append(Point(5, 10));
239
240 auto c2 = mec(pts_no_interior);
241 EXPECT_EQ(c.center, c2.center);
242 EXPECT_EQ(c.radius_squared, c2.radius_squared);
243}
244
245// ---------- Property tests ----------
246
248{
250 pts.append(Point(1, 1));
251 pts.append(Point(3, 7));
252 pts.append(Point(8, 2));
253 pts.append(Point(5, 9));
254 pts.append(Point(0, 4));
255 pts.append(Point(7, 6));
256 pts.append(Point(2, 8));
257
259 auto c = mec(pts);
260
261 EXPECT_TRUE(c.radius_squared >= 0);
262
263 for (DynList<Point>::Iterator it(pts); it.has_curr(); it.next_ne())
264 EXPECT_TRUE(c.contains(it.get_curr()));
265}
266
268{
270 auto c = mec({Point(0, 0), Point(4, 0), Point(4, 4), Point(0, 4)});
271
272 EXPECT_EQ(c.center, Point(2, 2));
273 EXPECT_EQ(c.radius_squared, 8);
274}
275
276// ---------- Static helpers ----------
277
279{
281 EXPECT_EQ(c.center, Point(3, 4));
282 // r^2 = 9 + 16 = 25
283 EXPECT_EQ(c.radius_squared, 25);
284}
285
287{
289 Point(0, 0), Point(5, 0), Point(10, 0));
290 EXPECT_EQ(c.center, Point(5, 0));
291 EXPECT_EQ(c.radius_squared, 25);
292}
293
294// ---------- Stress tests ----------
295
297{
298 // Place 20 points on a circle of radius 10 centered at (5, 5)
300 const double cx = 5.0, cy = 5.0, r = 10.0;
301 for (int i = 0; i < 20; ++i)
302 {
303 const double angle = 2.0 * M_PI * i / 20.0;
304 pts.append(Point(Geom_Number(cx + r * std::cos(angle)),
305 Geom_Number(cy + r * std::sin(angle))));
306 }
307
309 auto c = mec(pts);
310
311 // All points should be contained
312 for (DynList<Point>::Iterator it(pts); it.has_curr(); it.next_ne())
313 EXPECT_TRUE(c.contains(it.get_curr()));
314
315 // Radius should be close to 10
316 const Geom_Number tol = Geom_Number(1, 100000);
317 const Geom_Number diff = c.radius() > Geom_Number(10) ?
318 c.radius() - Geom_Number(10) : Geom_Number(10) - c.radius();
319 EXPECT_TRUE(diff < tol)
320 << "radius " << c.radius() << " expected ~10";
321}
322
324{
325 // 10x10 grid
327 for (int x = 0; x < 10; ++x)
328 for (int y = 0; y < 10; ++y)
329 pts.append(Point(x, y));
330
332 auto c = mec(pts);
333
334 for (DynList<Point>::Iterator it(pts); it.has_curr(); it.next_ne())
335 EXPECT_TRUE(c.contains(it.get_curr()));
336}
337
339{
340 std::mt19937 gen(42);
341 std::uniform_real_distribution<double> dist(-100.0, 100.0);
342
344 for (int i = 0; i < 200; ++i)
345 pts.append(Point(Geom_Number(dist(gen)), Geom_Number(dist(gen))));
346
348 auto c = mec(pts);
349
350 EXPECT_TRUE(c.radius_squared >= 0);
351
352 for (DynList<Point>::Iterator it(pts); it.has_curr(); it.next_ne())
353 EXPECT_TRUE(c.contains(it.get_curr()));
354}
long double h
Definition btreepic.C:154
Iterator on the items of list.
Definition htlist.H:1420
Dynamic singly linked list with functional programming support.
Definition htlist.H:1155
T & append(const T &item)
Definition htlist.H:1271
bool has_curr() const noexcept
Definition htlist.H:930
Smallest circle enclosing a point set (Welzl's algorithm).
static Circle from_two_points(const Point &a, const Point &b)
Smallest circle with a and b on its boundary (diameter).
static Circle from_three_points(const Point &a, const Point &b, const Point &c)
Circumscribed circle through three points.
Represents a point with rectangular coordinates in a 2D plane.
Definition point.H:229
TEST_F(GeomAlgorithmsTest, MEC_SinglePoint)
static mpfr_t y
Definition mpfr_mul_d.c:3
Main namespace for Aleph-w library functions.
Definition ah-arena.H:89
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.
Geom_Number square_root(const Geom_Number &x)
Square root of x (wrapper over mpfr).
Definition point.H:205
bool diff(const C1 &c1, const C2 &c2, Eq e=Eq())
Check if two containers differ.
mpq_class Geom_Number
Numeric type used by the geometry module.
Definition point.H:115
gsl_rng * r