1 #include <iostream>
2 #include <ctime>
3 #include "testlib/testlib_test.h"
4 //:
5 // \file
6 // \verbatim
7 // Modifications
8 // 2009-03-08 Peter Vanroose - Increased test coverage: added "inside" & "plane intersection" tests
9 // \endverbatim
10
11 #ifdef _MSC_VER
12 # include "vcl_msvc_warnings.h"
13 #endif
14 #include <cassert>
15 #include "vgl/vgl_triangle_3d.h"
16 #include "vgl/vgl_distance.h"
17 #include "vgl/vgl_plane_3d.h"
18 #include "vgl/vgl_line_segment_3d.h"
19
20 //==============================================================================
21 inline void
test_non_intersecting()22 test_non_intersecting()
23 {
24 vgl_point_3d<double> a_p1(0, 0, 0);
25 vgl_point_3d<double> a_p2(0, 4, 0);
26 vgl_point_3d<double> a_p3(4, 0, 0);
27
28 vgl_point_3d<double> b_p1(0, 0, 2);
29 vgl_point_3d<double> b_p2(0, 4, 2);
30 vgl_point_3d<double> b_p3(4, 0, 2);
31
32 vgl_line_segment_3d<double> iline;
33 unsigned edge_p1, edge_p2;
34 vgl_triangle_3d_intersection_t ret =
35 vgl_triangle_3d_triangle_intersection(a_p1, a_p2, a_p3, b_p1, b_p2, b_p3, iline, edge_p1, edge_p2);
36 TEST("Non-intersecting 1", ret, None);
37 ret = vgl_triangle_3d_triangle_intersection(a_p1, a_p2, a_p3, b_p1, b_p2, b_p3);
38 TEST("Non-intersecting 1 faster version", ret, None);
39
40 b_p1.set(-1, 0, 0);
41 b_p2.set(-1, 4, 0);
42 b_p3.set(-4, 0, 0);
43
44 ret = vgl_triangle_3d_triangle_intersection(a_p1, a_p2, a_p3, b_p1, b_p2, b_p3, iline, edge_p1, edge_p2);
45 TEST("Non-intersecting 2", ret, None);
46 ret = vgl_triangle_3d_triangle_intersection(a_p1, a_p2, a_p3, b_p1, b_p2, b_p3);
47 TEST("Non-intersecting 2 faster version", ret, None);
48 }
49
50
51 //==============================================================================
52 inline void
test_intersecting1()53 test_intersecting1()
54 {
55 vgl_point_3d<double> a_p1(0, 0, 0);
56 vgl_point_3d<double> a_p2(0, 4, 0);
57 vgl_point_3d<double> a_p3(4, 0, 0);
58
59 vgl_point_3d<double> b_p1(0, 0, 0);
60 vgl_point_3d<double> b_p2(0, 4, 0);
61 vgl_point_3d<double> b_p3(4, 4, 0);
62
63 vgl_triangle_3d_intersection_t ret = vgl_triangle_3d_triangle_intersection(a_p1, a_p2, a_p3, b_p1, b_p2, b_p3);
64 TEST("Intersecting 1 coplanar", ret, Coplanar);
65 vgl_line_segment_3d<double> iline;
66 unsigned edge_p1, edge_p2;
67 ret = vgl_triangle_3d_triangle_intersection(a_p1, a_p2, a_p3, b_p1, b_p2, b_p3, iline, edge_p1, edge_p2);
68 TEST("Intersecting coplanar with iline", ret, Coplanar);
69 }
70
71
72 //==============================================================================
73 inline void
test_intersecting2()74 test_intersecting2()
75 {
76 vgl_point_3d<double> a_p1(0, 0, 0);
77 vgl_point_3d<double> a_p2(0, 4, 0);
78 vgl_point_3d<double> a_p3(4, 0, 0);
79
80 vgl_point_3d<double> b_p1(1, 1, 4);
81 vgl_point_3d<double> b_p2(1, 1, -4);
82 vgl_point_3d<double> b_p3(0, 0, 0);
83
84 vgl_triangle_3d_intersection_t ret = vgl_triangle_3d_triangle_intersection(a_p1, a_p2, a_p3, b_p1, b_p2, b_p3);
85 TEST("Intersecting 2 Skew", ret, Skew);
86
87 vgl_line_segment_3d<double> i_line;
88 unsigned edge_p1, edge_p2;
89 ret = vgl_triangle_3d_triangle_intersection(a_p1, a_p2, a_p3, b_p1, b_p2, b_p3, i_line, edge_p1, edge_p2);
90 TEST("Intersecting Skew with iline", ret, Skew);
91 double p1_err = vgl_distance(i_line.point1(), vgl_point_3d<double>(0, 0, 0));
92 double p2_err = vgl_distance(i_line.point2(), vgl_point_3d<double>(1, 1, 0));
93 // handle ambiguity of intersecting line segment ordering.
94 if (p1_err > p2_err)
95 std::swap(edge_p1, edge_p2);
96 TEST_NEAR("Intersecting Skew, iline correct", p1_err + p2_err, 0, 1e-8);
97 // handle ambiguity of intersecting line point 1 , which lies on a vertex of each input triangle.
98 TEST("Intersecting Skew, edge label 1 correct", edge_p1 == 0 || edge_p1 == 2 || edge_p1 == 4 || edge_p1 == 5, true);
99 TEST("Intersecting Skew, edge label 2 correct", edge_p2 == 3, true);
100 }
101
102
103 //==============================================================================
104 inline void
test_intersecting3()105 test_intersecting3()
106 {
107 vgl_point_3d<double> a_p1(0, 0, 0);
108 vgl_point_3d<double> a_p2(0, 4, 0);
109 vgl_point_3d<double> a_p3(4, 0, 0);
110
111 vgl_point_3d<double> b_p1(0, 0, 0);
112 vgl_point_3d<double> b_p2(0, 4, 0);
113 vgl_point_3d<double> b_p3(0, 0, 4);
114
115 vgl_triangle_3d_intersection_t ret = vgl_triangle_3d_triangle_intersection(a_p1, a_p2, a_p3, b_p1, b_p2, b_p3);
116 TEST("Intersecting 3 Skew", ret, Skew);
117
118 vgl_line_segment_3d<double> i_line;
119 unsigned edge_p1, edge_p2;
120 ret = vgl_triangle_3d_triangle_intersection(a_p1, a_p2, a_p3, b_p1, b_p2, b_p3, i_line, edge_p1, edge_p2);
121
122 double p1_err = vgl_distance(i_line.point1(), vgl_point_3d<double>(0, 0, 0));
123 double p2_err = vgl_distance(i_line.point2(), vgl_point_3d<double>(0, 4, 0));
124 // handle ambiguity of intersecting line segment ordering.
125 if (p1_err > 1e-4)
126 {
127 std::swap(edge_p1, edge_p2);
128 p1_err = vgl_distance(i_line.point2(), vgl_point_3d<double>(0, 0, 0));
129 p2_err = vgl_distance(i_line.point1(), vgl_point_3d<double>(0, 4, 0));
130 }
131
132 TEST("Intersecting Skew - with isect line", ret, Skew);
133 TEST_NEAR("Intersecting Skew, iline correct", p1_err + p2_err, 0, 1e-8);
134 TEST("Intersecting Skew, edge label 1 correct", edge_p1 == 2 || edge_p1 == 5, true);
135 TEST("Intersecting Skew, edge label 2 correct", edge_p2 == 1 || edge_p2 == 4, true);
136 }
137
138
139 //==============================================================================
140 inline void
test_intersecting4()141 test_intersecting4()
142 {
143 vgl_point_3d<double> a_p1(0, 0, 0);
144 vgl_point_3d<double> a_p2(0, 4, 0);
145 vgl_point_3d<double> a_p3(4, 0, 0);
146
147 vgl_point_3d<double> b_p1(0, 0, 0);
148 vgl_point_3d<double> b_p2(0, 4, 0);
149 vgl_point_3d<double> b_p3(-4, 0, 0);
150
151 vgl_triangle_3d_intersection_t ret = vgl_triangle_3d_triangle_intersection(a_p1, a_p2, a_p3, b_p1, b_p2, b_p3);
152 TEST("Intersecting 4 coplanar", ret, Coplanar);
153
154 vgl_line_segment_3d<double> i_line;
155 unsigned edge_p1, edge_p2;
156 ret = vgl_triangle_3d_triangle_intersection(a_p1, a_p2, a_p3, b_p1, b_p2, b_p3, i_line, edge_p1, edge_p2);
157
158
159 double p1_err = vgl_distance(i_line.point1(), vgl_point_3d<double>(0, 0, 0));
160 double p2_err = vgl_distance(i_line.point2(), vgl_point_3d<double>(0, 4, 0));
161 // handle ambiguity of intersecting line segment ordering.
162 if (p1_err > 1e-4)
163 {
164 std::swap(edge_p1, edge_p2);
165 p1_err = vgl_distance(i_line.point2(), vgl_point_3d<double>(0, 0, 0));
166 p2_err = vgl_distance(i_line.point1(), vgl_point_3d<double>(0, 4, 0));
167 }
168
169 TEST("Intersecting Coplanar - with isect line", ret, Coplanar);
170 TEST_NEAR("Intersecting Coplanar, iline correct", p1_err + p2_err, 0, 1e-8);
171 // Due to coplanarity, edge labels are very ambiguous
172 TEST(
173 "Intersecting Coplanar, edge label 1 correct", edge_p1 == 2 || edge_p1 == 5 || edge_p1 == 0 || edge_p1 == 3, true);
174 TEST(
175 "Intersecting Coplanar, edge label 2 correct", edge_p2 == 1 || edge_p2 == 4 || edge_p2 == 0 || edge_p2 == 3, true);
176 }
177
178
179 inline void
test_intersecting5_arrangement(const vgl_point_3d<double> & a_A,const vgl_point_3d<double> & a_B,const vgl_point_3d<double> & a_C,const vgl_point_3d<double> & b_A,const vgl_point_3d<double> & b_B,const vgl_point_3d<double> & b_C,unsigned expected_edge_1,unsigned expected_edge_2)180 test_intersecting5_arrangement(const vgl_point_3d<double> & a_A,
181 const vgl_point_3d<double> & a_B,
182 const vgl_point_3d<double> & a_C,
183 const vgl_point_3d<double> & b_A,
184 const vgl_point_3d<double> & b_B,
185 const vgl_point_3d<double> & b_C,
186 unsigned expected_edge_1,
187 unsigned expected_edge_2)
188 {
189 vgl_triangle_3d_intersection_t ret = vgl_triangle_3d_triangle_intersection(a_A, a_B, a_C, b_A, b_B, b_C);
190 TEST("Intersecting Skew", ret, Skew);
191
192 vgl_line_segment_3d<double> i_line;
193 unsigned edge_p1, edge_p2;
194 ret = vgl_triangle_3d_triangle_intersection(a_A, a_B, a_C, b_A, b_B, b_C, i_line, edge_p1, edge_p2);
195
196 double p1_err = vgl_distance(i_line.point1(), vgl_point_3d<double>(1, 1, 0));
197 double p2_err = vgl_distance(i_line.point2(), vgl_point_3d<double>(2, 2, 0));
198 // handle ambiguity of intersecting line segment ordering.
199 if (p1_err > 1e-2)
200 {
201 std::swap(edge_p1, edge_p2);
202 p1_err = vgl_distance(i_line.point2(), vgl_point_3d<double>(1, 1, 0));
203 p2_err = vgl_distance(i_line.point1(), vgl_point_3d<double>(2, 2, 0));
204 }
205
206 TEST("Intersecting Skew - with isect line", ret, Skew);
207 TEST_NEAR("Intersecting Skew, iline correct", p1_err + p2_err, 0, 1e-8);
208 TEST("Intersecting Skew, edge label 1 correct", edge_p1 == expected_edge_1, true);
209 TEST("Intersecting Skew, edge label 2 correct", edge_p2 == expected_edge_2, true);
210 }
211
212
213 //==============================================================================
214 inline void
test_intersecting5()215 test_intersecting5()
216 {
217 vgl_point_3d<double> a_p1(0, 0, 0);
218 vgl_point_3d<double> a_p2(0, 4, 0);
219 vgl_point_3d<double> a_p3(4, 0, 0);
220
221 vgl_point_3d<double> b_p1(4, 4, 0);
222 vgl_point_3d<double> b_p2(1, 1, -4);
223 vgl_point_3d<double> b_p3(1, 1, 4);
224
225 std::cout << "Test skew intersection with multiple arrangements.\n";
226
227 test_intersecting5_arrangement(a_p1, a_p2, a_p3, b_p1, b_p2, b_p3, 4, 1);
228 test_intersecting5_arrangement(a_p2, a_p1, a_p3, b_p1, b_p2, b_p3, 4, 2);
229 test_intersecting5_arrangement(a_p2, a_p3, a_p1, b_p1, b_p2, b_p3, 4, 0);
230 test_intersecting5_arrangement(a_p3, a_p1, a_p2, b_p1, b_p2, b_p3, 4, 2);
231 test_intersecting5_arrangement(a_p1, a_p2, a_p3, b_p2, b_p3, b_p1, 3, 1);
232 test_intersecting5_arrangement(a_p2, a_p1, a_p3, b_p2, b_p3, b_p1, 3, 2);
233 test_intersecting5_arrangement(a_p2, a_p3, a_p1, b_p2, b_p3, b_p1, 3, 0);
234 test_intersecting5_arrangement(a_p3, a_p1, a_p2, b_p2, b_p3, b_p1, 3, 2);
235 test_intersecting5_arrangement(a_p1, a_p2, a_p3, b_p3, b_p1, b_p2, 5, 1);
236 test_intersecting5_arrangement(a_p2, a_p1, a_p3, b_p3, b_p1, b_p2, 5, 2);
237 test_intersecting5_arrangement(a_p2, a_p3, a_p1, b_p3, b_p1, b_p2, 5, 0);
238 test_intersecting5_arrangement(a_p3, a_p1, a_p2, b_p3, b_p1, b_p2, 5, 2);
239 test_intersecting5_arrangement(a_p1, a_p2, a_p3, b_p2, b_p1, b_p3, 5, 1);
240 test_intersecting5_arrangement(a_p2, a_p1, a_p3, b_p2, b_p1, b_p3, 5, 2);
241 test_intersecting5_arrangement(a_p2, a_p3, a_p1, b_p2, b_p1, b_p3, 5, 0);
242 test_intersecting5_arrangement(a_p3, a_p1, a_p2, b_p2, b_p1, b_p3, 5, 2);
243
244 test_intersecting5_arrangement(b_p1, b_p2, b_p3, a_p1, a_p2, a_p3, 1, 4);
245 test_intersecting5_arrangement(b_p1, b_p2, b_p3, a_p2, a_p1, a_p3, 1, 5);
246 test_intersecting5_arrangement(b_p1, b_p2, b_p3, a_p2, a_p3, a_p1, 1, 3);
247 test_intersecting5_arrangement(b_p1, b_p2, b_p3, a_p3, a_p1, a_p2, 1, 5);
248 test_intersecting5_arrangement(b_p2, b_p3, b_p1, a_p1, a_p2, a_p3, 0, 4);
249 test_intersecting5_arrangement(b_p2, b_p3, b_p1, a_p2, a_p1, a_p3, 0, 5);
250 test_intersecting5_arrangement(b_p2, b_p3, b_p1, a_p2, a_p3, a_p1, 0, 3);
251 test_intersecting5_arrangement(b_p2, b_p3, b_p1, a_p3, a_p1, a_p2, 0, 5);
252 test_intersecting5_arrangement(b_p3, b_p1, b_p2, a_p1, a_p2, a_p3, 2, 4);
253 test_intersecting5_arrangement(b_p3, b_p1, b_p2, a_p2, a_p1, a_p3, 2, 5);
254 test_intersecting5_arrangement(b_p3, b_p1, b_p2, a_p2, a_p3, a_p1, 2, 3);
255 test_intersecting5_arrangement(b_p3, b_p1, b_p2, a_p3, a_p1, a_p2, 2, 5);
256 test_intersecting5_arrangement(b_p2, b_p1, b_p3, a_p1, a_p2, a_p3, 2, 4);
257 test_intersecting5_arrangement(b_p2, b_p1, b_p3, a_p2, a_p1, a_p3, 2, 5);
258 test_intersecting5_arrangement(b_p2, b_p1, b_p3, a_p2, a_p3, a_p1, 2, 3);
259 test_intersecting5_arrangement(b_p2, b_p1, b_p3, a_p3, a_p1, a_p2, 2, 5);
260 }
261
262
263 //==============================================================================
264 inline void
test_intersecting_degenerate_triangles1()265 test_intersecting_degenerate_triangles1()
266 {
267 std::cout << '\n'
268 << "----------------------------------------------\n"
269 << " test vgl_triangle_3d_triangle_intersection()\n"
270 << "----------------------------------------------\n";
271
272 // the valid triangle
273 vgl_point_3d<double> valid_p1(0, 0, 0);
274 vgl_point_3d<double> valid_p2(0, 4, 0);
275 vgl_point_3d<double> valid_p3(4, 0, 0);
276
277 // Test on a number of different types of degenerate triangle
278
279 constexpr unsigned NUM_TESTS = 7;
280 constexpr unsigned NUM_VARIANTS = 3;
281
282 vgl_point_3d<double> degen[NUM_TESTS][3] = {
283 { vgl_point_3d<double>(-4, 5, 0),
284 vgl_point_3d<double>(-3, 4, 0),
285 vgl_point_3d<double>(-2, 0, 0) }, // Coplanar non-intersecting
286 { vgl_point_3d<double>(4, 5, 0),
287 vgl_point_3d<double>(-3, 4, 0),
288 vgl_point_3d<double>(2, -1, 0) }, // Coplanar non-intersecting
289 { vgl_point_3d<double>(2, 5, 0),
290 vgl_point_3d<double>(2, 1, 0),
291 vgl_point_3d<double>(2, -3, 0) }, // Coplanar intersecting
292 { vgl_point_3d<double>(0, 8, 0),
293 vgl_point_3d<double>(0, 6, 0),
294 vgl_point_3d<double>(0, 5, 0) }, // Coplanar non-intersecting edge collinear
295 { vgl_point_3d<double>(0, 8, 0),
296 vgl_point_3d<double>(0, 3, 0),
297 vgl_point_3d<double>(0, 5, 0) }, // Coplanar intersecting edge collinear
298 { vgl_point_3d<double>(2, 1, 5),
299 vgl_point_3d<double>(2, 1, 3),
300 vgl_point_3d<double>(2, 1, 1) }, // Non-coplanar non-intersecting
301 { vgl_point_3d<double>(2, 1, 3),
302 vgl_point_3d<double>(2, 1, 1),
303 vgl_point_3d<double>(2, 1, -2) } // Non-coplanar intersecting
304 };
305
306 const char * degen_desc[NUM_TESTS] = { "Coplanar non-intersecting possibly non-degenerate",
307 "Coplanar intersecting possibly non-degenerate",
308 "Coplanar intersecting degenerate",
309 "Coplanar non-intersecting edge collinear degenerate",
310 "Coplanar intersecting edge collinear degenerate",
311 "Non-coplanar non-intersecting degenerate",
312 "Non-coplanar intersecting degenerate" };
313
314 vgl_triangle_3d_intersection_t exp_result[NUM_TESTS][NUM_VARIANTS] = {
315 { None, None, None }, { None, Coplanar, Coplanar }, { Coplanar, Coplanar, Coplanar },
316 { None, None, None }, { Coplanar, Coplanar, Coplanar }, { None, None, None },
317 { None, Skew, Skew }
318 };
319
320 // Expected edges [0-5]=>correct values. 7=>Any repreated values would be correct
321 // 8=>Any repeated values from the degenerate triange would be correct
322 // 9=>No expected answer( since no intersection)
323 unsigned exp_edges[NUM_TESTS][NUM_VARIANTS][2] = {
324 { { 9, 9 }, { 9, 9 }, { 9, 9 } }, // unspecified if non-intersecting
325 { { 9, 9 }, { 8, 8 }, { 7, 7 } }, // repeated edges if coplanar
326 { { 3, 3 }, { 8, 8 }, { 8, 8 } }, // repeated edges if coplanar
327 { { 9, 9 }, { 9, 9 }, { 9, 9 } }, // unspecified if non-intersecting
328 { { 3, 3 }, { 8, 8 }, { 8, 8 } }, { { 9, 9 }, { 9, 9 }, { 9, 9 } }, // unspecified if non-intersecting
329 { { 9, 9 }, { 4, 5 }, { 4, 5 } }
330 };
331
332
333 unsigned degen_tests[NUM_VARIANTS][3] = {
334 { 1, 1, 1 }, // all vertices the same point
335 { 1, 1, 2 }, // one vertex distinct
336 { 0, 1, 2 } // all vertices distinct
337 };
338
339 const char * degen_tests_desc[3] = { "point", "one vertex distinct", "all vertices distinct" };
340
341 for (unsigned i = 0; i < NUM_TESTS; ++i)
342 {
343 for (unsigned j = 0; j < NUM_VARIANTS; ++j)
344 {
345 if (exp_result[i][j] != None) // Check for any unwarranted looseness in test.
346 assert(exp_edges[i][j][0] != 9);
347
348 std::string test_desc = degen_desc[i];
349 test_desc += " (";
350 test_desc += degen_tests_desc[j];
351 test_desc += ")";
352
353 vgl_line_segment_3d<double> i_line;
354 auto point1_edge = (unsigned)-1, point2_edge = (unsigned)-2;
355 vgl_triangle_3d_intersection_t ret = vgl_triangle_3d_triangle_intersection(valid_p1,
356 valid_p2,
357 valid_p3,
358 degen[i][degen_tests[j][0]],
359 degen[i][degen_tests[j][1]],
360 degen[i][degen_tests[j][2]],
361 i_line,
362 point1_edge,
363 point2_edge);
364
365 TEST(test_desc.c_str(), ret, exp_result[i][j]);
366 if (exp_edges[i][j][0] != 9)
367 {
368 if (point1_edge > point2_edge) // Canonicalise edge results.
369 std::swap(point1_edge, point2_edge);
370 if (exp_edges[i][j][0] == 8)
371 TEST("Edge 1 == Edge 2 for no definite edge result", point1_edge == point2_edge && point1_edge >= 3, true);
372 else if (exp_edges[i][j][0] == 7)
373 TEST("Edge 1 == Edge 2 for no definite edge result", point1_edge == point2_edge, true);
374 else
375 {
376 TEST_NEAR("Edge 1 correct", point1_edge, exp_edges[i][j][0], 0);
377 TEST_NEAR("Edge 2 correct", point2_edge, exp_edges[i][j][1], 0);
378 }
379 }
380 test_desc += " reversed";
381 ret = vgl_triangle_3d_triangle_intersection(degen[i][degen_tests[j][0]],
382 degen[i][degen_tests[j][1]],
383 degen[i][degen_tests[j][2]],
384 valid_p1,
385 valid_p2,
386 valid_p3,
387 i_line,
388 point1_edge,
389 point2_edge);
390
391
392 TEST(test_desc.c_str(), ret, exp_result[i][j]);
393 if (exp_edges[i][j][0] != 9)
394 {
395 if (point1_edge > point2_edge) // Canonicalise edge results.
396 std::swap(point1_edge, point2_edge);
397 if (exp_edges[i][j][0] == 8)
398 TEST("Edge 1 == Edge 2 for no definite edge result", point1_edge == point2_edge && point1_edge < 3, true);
399 else if (exp_edges[i][j][0] == 7)
400 TEST("Edge 1 == Edge 2 for no definite edge result", point1_edge == point2_edge, true);
401 else
402 {
403 // modify exp_edge values to point to other triangle.
404 TEST_NEAR("Edge 1 correct", point1_edge, (exp_edges[i][j][0] + 3) % 6, 0);
405 TEST_NEAR("Edge 2 correct", point2_edge, (exp_edges[i][j][1] + 3) % 6, 0);
406 }
407 }
408 }
409 }
410 }
411
412
413 //==============================================================================
414 inline void
test_intersecting_degenerate_triangles2()415 test_intersecting_degenerate_triangles2()
416 {
417 std::cout << '\n'
418 << "----------------------------------------------\n"
419 << " test vgl_triangle_3d_triangle_intersection()\n"
420 << "----------------------------------------------\n";
421
422 // test intersection of 2 degenerate triangles
423 vgl_point_3d<double> d1_p1(0, 0, 0);
424 vgl_point_3d<double> d1_p2(2, 0, 0);
425 vgl_point_3d<double> d1_p3(4, 0, 0);
426
427 vgl_point_3d<double> d2_p1(-1, 0, 0);
428 vgl_point_3d<double> d2_p2(-3, 0, 0);
429 vgl_point_3d<double> d2_p3(-2, 0, 0);
430
431 vgl_triangle_3d_intersection_t ret = vgl_triangle_3d_triangle_intersection(d1_p1, d1_p2, d1_p3, d2_p1, d2_p2, d2_p3);
432
433 TEST("Non-Intersecting two degenerate triangles", ret, None);
434
435 ret = vgl_triangle_3d_triangle_intersection(d2_p1, d2_p2, d2_p3, d1_p1, d1_p2, d1_p3);
436
437 TEST("Non-Intersecting two degenerate triangles reversed", ret, None);
438
439 d2_p1.set(0, 0, 0);
440 d2_p2.set(1, 1, 1);
441 d2_p3.set(2, 2, 2);
442
443 ret = vgl_triangle_3d_triangle_intersection(d1_p1, d1_p2, d1_p3, d2_p1, d2_p2, d2_p3);
444
445 TEST("Intersecting two degenerate triangles", ret, Coplanar);
446
447 ret = vgl_triangle_3d_triangle_intersection(d2_p1, d2_p2, d2_p3, d1_p1, d1_p2, d1_p3);
448
449 TEST("Intersecting two degenerate triangles reversed", ret, Coplanar);
450 }
451
452
453 //==============================================================================
454 inline void
test_coincident_edges1()455 test_coincident_edges1()
456 {
457 std::cout << '\n'
458 << "-----------------------------------------\n"
459 << " test vgl_triangle_3d_coincident_edges()\n"
460 << "-----------------------------------------\n";
461
462 vgl_point_3d<double> a_p1(0, 0, 0);
463 vgl_point_3d<double> a_p2(0, 4, 0);
464 vgl_point_3d<double> a_p3(4, 0, 0);
465
466 vgl_point_3d<double> b_p1(0, 0, 0);
467 vgl_point_3d<double> b_p2(-4, 0, 0);
468 vgl_point_3d<double> b_p3(0, -4, 0);
469
470 std::vector<std::pair<unsigned, unsigned>> coinc =
471 vgl_triangle_3d_coincident_edges(a_p1, a_p2, a_p3, b_p1, b_p2, b_p3);
472
473 TEST("Coincident_edges non-coincident", coinc.empty(), true);
474 }
475
476
477 //==============================================================================
478 inline void
test_coincident_edges2()479 test_coincident_edges2()
480 {
481 vgl_point_3d<double> a_p1(0, 0, 0);
482 vgl_point_3d<double> a_p2(0, 4, 0);
483 vgl_point_3d<double> a_p3(4, 0, 0);
484
485 vgl_point_3d<double> b_p1(1, 0, 0);
486 vgl_point_3d<double> b_p2(0, 0, 4);
487 vgl_point_3d<double> b_p3(2, 0, 0);
488
489 std::vector<std::pair<unsigned, unsigned>> coinc =
490 vgl_triangle_3d_coincident_edges(a_p1, a_p2, a_p3, b_p1, b_p2, b_p3);
491
492 std::pair<unsigned, unsigned> exp_edge(2, 2);
493
494 TEST("Coincident_edges 1 coincident edge", coinc.size() == 1 && coinc[0] == exp_edge, true);
495 }
496
497 //==============================================================================
498 inline void
test_intersect_plane()499 test_intersect_plane()
500 {
501 std::cout << '\n'
502 << "-------------------------------------------\n"
503 << " test vgl_triangle_3d_plane_intersection()\n"
504 << "-------------------------------------------\n";
505
506 vgl_point_3d<double> p1(0, 0, 0);
507 vgl_point_3d<double> p2(0, 15, 25);
508 vgl_point_3d<double> p3(0, 45, 10);
509
510 vgl_plane_3d<double> pl(0, 0, 1, -20); // the plane z=20 ==> normal intersection
511 vgl_line_segment_3d<double> l; // the returned intersection
512 vgl_point_3d<double> p(0, 12, 20), q(0, 25, 20);
513 vgl_line_segment_3d<double> l0(p, q); // the correct intersection
514 vgl_triangle_3d_intersection_t ret;
515
516 ret = vgl_triangle_3d_plane_intersection(p1, p2, p3, pl, l);
517 TEST("normal intersection", ret == Skew && l == l0, true);
518
519 pl.set(0, 0, 1, 1); // z=-1 ==> no intersection
520 ret = vgl_triangle_3d_plane_intersection(p1, p2, p3, pl, l);
521 TEST("no intersection", ret, None);
522
523 pl.set(0, 0, 1, -25); // z=25 ==> corner point intersection
524 l0.set(p2, p2);
525 ret = vgl_triangle_3d_plane_intersection(p1, p2, p3, pl, l);
526 TEST("corner point intersection", ret == Skew && l == l0, true);
527
528 pl.set(1, 0, 0, 0); // x=0 ==> coplanar with triangle
529 ret = vgl_triangle_3d_plane_intersection(p1, p2, p3, pl, l);
530 TEST("coplanar intersection", ret, Coplanar);
531 }
532
533 //==============================================================================
534 inline void
test_point_containment()535 test_point_containment()
536 {
537 std::cout << '\n'
538 << "------------------------------------\n"
539 << " test vgl_triangle_3d_test_inside()\n"
540 << "------------------------------------\n";
541
542
543 // three arbitrary coplanar points (note that coplanarity is assumed by the two methods!)
544 {
545 vgl_point_3d<double> p1(0, 0, 0);
546 vgl_point_3d<double> p2(5, 0, 3);
547 vgl_point_3d<double> p3(2, 0, 10);
548 vgl_point_3d<double> in1 = centre(p1, p2, p3); // barycentre
549 vgl_point_3d<double> in2(2.5, 0, 1.51); // point almost on the edge p1--p2
550 vgl_point_3d<double> in3(0.01, 0, 0.01); // almost a corner point
551 vgl_point_3d<double> out1(2, 0, 1); // same plane, outside triangle
552 vgl_point_3d<double> out2(2, 1, 1); // not coplanar, still outside triangle
553 vgl_point_3d<double> in4(0.01, 0.1, 0.01); // not quite coplanar, but otherwise inside triangle
554 TEST("inside 1 - barycentric method", vgl_triangle_3d_test_inside(in1, p1, p2, p3), true);
555 TEST("inside 1 - barycentric method with coplanar tol", vgl_triangle_3d_test_inside(in1, p1, p2, p3, 1e-10), true);
556 TEST("inside 1 - cosine method", vgl_triangle_3d_test_inside_simple(in1, p1, p2, p3), true);
557 TEST("inside 2 - barycentric method", vgl_triangle_3d_test_inside(in2, p1, p2, p3), true);
558 TEST("inside 2 - barycentric method with coplanar tol", vgl_triangle_3d_test_inside(in2, p1, p2, p3, 1e-10), true);
559 TEST("inside 2 - cosine method", vgl_triangle_3d_test_inside_simple(in2, p1, p2, p3), true);
560 TEST("inside 3 - barycentric method", vgl_triangle_3d_test_inside(in3, p1, p2, p3), true);
561 TEST("inside 3 - cosine method", vgl_triangle_3d_test_inside_simple(in3, p1, p2, p3), true);
562 TEST("not inside 1 - barycentric method", vgl_triangle_3d_test_inside(out1, p1, p2, p3), false);
563 TEST("not inside 1 - cosine method", vgl_triangle_3d_test_inside_simple(out1, p1, p2, p3), false);
564 TEST("not inside 2 - barycentric method", vgl_triangle_3d_test_inside(out2, p1, p2, p3), false);
565 TEST("not inside 2 - cosine method", vgl_triangle_3d_test_inside_simple(out2, p1, p2, p3), false);
566 TEST("sort of inside 4 - barycentric method with coplanar tol",
567 vgl_triangle_3d_test_inside(in4, p1, p2, p3, 0.2),
568 true);
569 TEST("sort of inside 4 - barycentric method", vgl_triangle_3d_test_inside(in4, p1, p2, p3), false);
570 }
571 // three collinear points -- the cosine method will always fail in this case
572 {
573 vgl_point_3d<double> p1(0, 0, 0);
574 vgl_point_3d<double> p2(5, 0, 3);
575 vgl_point_3d<double> p3(15, 0, 9);
576 vgl_point_3d<double> in1(2.5, 0, 1.5); // barycentre
577 vgl_point_3d<double> in2(2, 0, 1.2);
578 vgl_point_3d<double> in3(0.05, 0, 0.03); // almost a corner point
579 vgl_point_3d<double> out1(15.05, 0, 9.03); // collinear, almost a corner point
580 vgl_point_3d<double> out2(20, 0, 12); // collinear
581 vgl_point_3d<double> out3(2, 0, 1); // not collinear
582 TEST("inside 1 (collinear points) - barycentric method", vgl_triangle_3d_test_inside(in1, p1, p2, p3), true);
583 TEST("inside 2 (collinear points) - barycentric method", vgl_triangle_3d_test_inside(in2, p1, p2, p3), true);
584 TEST("inside 3 (collinear points) - barycentric method", vgl_triangle_3d_test_inside(in3, p1, p2, p3), true);
585 TEST("not inside 1 (collinear points) - barycentric method", vgl_triangle_3d_test_inside(out1, p1, p2, p3), false);
586 TEST("not inside 2 (collinear points) - barycentric method", vgl_triangle_3d_test_inside(out2, p1, p2, p3), false);
587 TEST("not inside 3 (collinear points) - barycentric method", vgl_triangle_3d_test_inside(out3, p1, p2, p3), false);
588 }
589 // three coincident points
590 {
591 vgl_point_3d<double> p1(5, 0, 3);
592 vgl_point_3d<double> p2(5, 0, 3);
593 vgl_point_3d<double> p3(5, 0, 3);
594 vgl_point_3d<double> in1(5, 0, 3);
595 vgl_point_3d<double> out1(5, 0, 2);
596 TEST("inside (coincident points) - barycentric method", vgl_triangle_3d_test_inside(in1, p1, p2, p3), true);
597 TEST("not inside (coincident points) - barycentric method", vgl_triangle_3d_test_inside(out1, p1, p2, p3), false);
598 }
599 // almost degenerate triangle
600 {
601 vgl_point_3d<double> p1(-250.00000027628425, 307.98258024874724, 359.57738152050911);
602 vgl_point_3d<double> p2(-239.99999999497970, 314.64924710295026, 359.57738152050911);
603 vgl_point_3d<double> p3(-240.00000000000000, 314.64924709960337, 359.57738152050911);
604
605 vgl_point_3d<double> in1(-245.00000013828425, 311.31591367584872, 359.57738152050911);
606 vgl_point_3d<double> out1(-245.00000013828425, 311.31591367584872, 359.575);
607 // TEST("inside (High AR triangle) - barycentric method", vgl_triangle_3d_test_inside(in1,p1,p2,p3), true);
608 TEST("not inside (High AR triangle) - barycentric method", vgl_triangle_3d_test_inside(out1, p1, p2, p3), false);
609 }
610 }
611
612 //==============================================================================
613 //: A simple performance test on the different triangle point containment algorithms
614 inline void
test_point_containment_algo_perf()615 test_point_containment_algo_perf()
616 {
617 vgl_point_3d<double> p1(0, 0, 0);
618 vgl_point_3d<double> p2(5, 0, 3);
619 vgl_point_3d<double> p3(2, 0, 10);
620
621 vgl_point_3d<double> test_pt = centre(p1, p2, p3);
622
623 unsigned tries = 10000000;
624
625 std::clock_t bary_st = std::clock();
626 for (unsigned i = 0; i < tries; ++i)
627 vgl_triangle_3d_test_inside(test_pt, p1, p2, p3);
628 std::clock_t bary_en = std::clock();
629
630 std::clock_t ang_st = std::clock();
631 for (unsigned i = 0; i < tries; ++i)
632 vgl_triangle_3d_test_inside_simple(test_pt, p1, p2, p3);
633 std::clock_t ang_en = std::clock();
634
635 unsigned long bary_time = bary_en - bary_st;
636 unsigned long ang_time = ang_en - ang_st;
637
638 unsigned long bary_ps = (tries * CLOCKS_PER_SEC) / bary_time;
639 unsigned long ang_ps = (tries * CLOCKS_PER_SEC) / ang_time;
640
641 std::cout << "Barycentric method: " << bary_ps << "/sec (" << bary_time << " ticks)\n"
642 << "Angles method: " << ang_ps << "/sec (" << ang_time << " ticks)\n";
643 std::cout.precision(2);
644 std::cout << "Barycentric " << std::fixed << double(bary_ps) / ang_ps << " times faster than angles\n";
645
646 // usually get Barycentric about 2.5 to 3 times faster than angles KOM (Aug 07)
647 }
648
649
650 //==============================================================================
651 // Test function vgl_triangle_3d_closest_point()
652 //==============================================================================
653 inline void
test_closest_point()654 test_closest_point()
655 {
656 std::cout << '\n'
657 << "----------------------\n"
658 << " test_closest_point()\n"
659 << "----------------------\n";
660
661 {
662 // Consider a triangle in the XY plane.
663 vgl_point_3d<double> p1(0, 0, 0);
664 vgl_point_3d<double> p2(1, 0, 0);
665 vgl_point_3d<double> p3(0, 1, 0);
666
667 // Test one of the actual vertices
668 {
669 vgl_point_3d<double> q(0, 0, 0);
670 vgl_point_3d<double> c = vgl_triangle_3d_closest_point(q, p1, p2, p3);
671 TEST_NEAR("Vertex point", (c - p1).length(), 0.0, 1e-6);
672 }
673
674 // Test a point along a triangle edge
675 {
676 vgl_point_3d<double> q(0.5, 0.0, 0.0);
677 vgl_point_3d<double> c = vgl_triangle_3d_closest_point(q, p1, p2, p3);
678 TEST_NEAR("Edge point", (c - q).length(), 0.0, 1e-6);
679 }
680
681 // Test a point within the triangle
682 {
683 vgl_point_3d<double> q(0.1, 0.1, 0.0);
684 vgl_point_3d<double> c = vgl_triangle_3d_closest_point(q, p1, p2, p3);
685 TEST_NEAR("Point inside", (c - q).length(), 0.0, 1e-6);
686 }
687
688 // Test a point outside the triangle (but in same plane)
689 {
690 vgl_point_3d<double> q(-1.0, 0.0, 0.0);
691 vgl_point_3d<double> c = vgl_triangle_3d_closest_point(q, p1, p2, p3);
692 TEST_NEAR("Point outside but same plane", (c - p1).length(), 0.0, 1e-6);
693 }
694
695 // Test a point outside the triangle (in a different plane, but would project inside the triangle)
696 {
697 vgl_point_3d<double> q(0.1, 0.1, 1.0);
698 vgl_point_3d<double> c = vgl_triangle_3d_closest_point(q, p1, p2, p3);
699 vgl_point_3d<double> r(0.1, 0.1, 0.0); // The expected closest point
700 TEST_NEAR("Point inside", (c - r).length(), 0.0, 1e-6);
701 }
702
703 // Test a point outside the triangle (in a different plane, and would project outside the triangle)
704 {
705 vgl_point_3d<double> q(-1.0, 0.0, 1.0);
706 vgl_point_3d<double> c = vgl_triangle_3d_closest_point(q, p1, p2, p3);
707 TEST_NEAR("Point outside and different plane", (c - p1).length(), 0.0, 1e-6);
708 }
709 }
710 TEST("To be sure that we get here (Mac)", 1, 1);
711
712 // Nearly degenerate triangle
713 {
714 vgl_point_3d<double> p1(-250.00000027628425, 307.98258024874724, 359.57738152050911);
715 vgl_point_3d<double> p2(-239.99999999497970, 314.64924710295026, 359.57738152050911);
716 vgl_point_3d<double> p3(-240.00000000000000, 314.64924709960337, 359.57738152050911);
717 std::cout << "Nearly degenerate triangle " << p1 << p2 << p3 << std::endl;
718
719 vgl_point_3d<double> in1(-245.00000013828425, 311.31591367584872, 359.57738152050911);
720 vgl_point_3d<double> out1(-245.00000013828425, 311.31591367584872, 380);
721 vgl_point_3d<double> c = vgl_triangle_3d_closest_point(in1, p1, p2, p3);
722 TEST_NEAR("Point inside High AR triangle", vgl_distance(c, in1), 0, 1e-8);
723 c = vgl_triangle_3d_closest_point(out1, p1, p2, p3);
724 TEST_NEAR("Point outside High AR triangle", vgl_distance(c, in1), 0, 1e-8);
725 }
726 }
727
728
729 //==============================================================================
730 // Test function vgl_triangle_3d_distance()
731 //==============================================================================
732 inline void
test_distance()733 test_distance()
734 {
735 std::cout << '\n'
736 << "---------------------------------\n"
737 << " test vgl_triangle_3d_distance()\n"
738 << "---------------------------------\n";
739
740 // Consider a triangle in the XY plane.
741 vgl_point_3d<double> p1(0, 0, 0);
742 vgl_point_3d<double> p2(1, 0, 0);
743 vgl_point_3d<double> p3(0, 1, 0);
744
745 // Test one of the actual vertices
746 {
747 vgl_point_3d<double> q(0, 0, 0);
748 double d = vgl_triangle_3d_distance(q, p1, p2, p3);
749 double true_d = 0.0;
750 TEST_NEAR("Vertex point", d, true_d, 1e-6);
751 }
752
753 // Test a point along a triangle edge
754 {
755 vgl_point_3d<double> q(0.5, 0.0, 0.0);
756 double d = vgl_triangle_3d_distance(q, p1, p2, p3);
757 double true_d = 0.0;
758 TEST_NEAR("Edge point", d, true_d, 1e-6);
759 }
760
761 // Test a point within the triangle
762 {
763 vgl_point_3d<double> q(0.1, 0.1, 0.0);
764 double d = vgl_triangle_3d_distance(q, p1, p2, p3);
765 double true_d = 0.0;
766 TEST_NEAR("Point inside", d, true_d, 1e-6);
767 }
768
769 // Test a point outside the triangle (but in same plane)
770 {
771 vgl_point_3d<double> q(-1.0, 0.0, 0.0);
772 double d = vgl_triangle_3d_distance(q, p1, p2, p3);
773 double true_d = 1.0;
774 TEST_NEAR("Point outside but same plane", d, true_d, 1e-6);
775 }
776
777 // Test a point outside the triangle (in a different plane, but would project inside the triangle)
778 {
779 vgl_point_3d<double> q(0.1, 0.1, 1.0);
780 double d = vgl_triangle_3d_distance(q, p1, p2, p3);
781 double true_d = 1.0;
782 TEST_NEAR("Point inside", d, true_d, 1e-6);
783 }
784
785 // Test a point outside the triangle (in a different plane, and would project outside the triangle)
786 {
787 vgl_point_3d<double> q(-1.0, 0.0, 1.0);
788 double d = vgl_triangle_3d_distance(q, p1, p2, p3);
789 double true_d = 1.4142135623730950488; // sqrt2
790 TEST_NEAR("Point outside and different plane", d, true_d, 1e-6);
791 }
792 }
793
794
795 //==============================================================================
796 // Test function vgl_triangle_3d_area()
797 //==============================================================================
798 inline void
test_area()799 test_area()
800 {
801 std::cout << '\n'
802 << "-----------------------------\n"
803 << " test vgl_triangle_3d_area()\n"
804 << "-----------------------------\n";
805
806 vgl_point_3d<double> p1(0, 0, 0);
807 vgl_point_3d<double> p2(1, 0, 0);
808 vgl_point_3d<double> p3(0, 1, 0);
809
810 double area = vgl_triangle_3d_area(p1, p2, p3);
811 TEST_NEAR("Triangle area", area, 0.5, 1e-6);
812 }
813
814
815 //==============================================================================
816 // Test function vgl_triangle_3d_aspect_ratio()
817 //==============================================================================
818 inline void
test_aspect_ratio()819 test_aspect_ratio()
820 {
821 std::cout << '\n'
822 << "-------------------------------------\n"
823 << " test vgl_triangle_3d_aspect_ratio()\n"
824 << "-------------------------------------\n";
825
826 vgl_point_3d<double> p1(0, 0, 0);
827 vgl_point_3d<double> p2(1, 0, 0);
828 vgl_point_3d<double> p3(0, 1, 0);
829
830 double ratio = vgl_triangle_3d_aspect_ratio(p1, p2, p3);
831 TEST_NEAR("Triangle area", ratio, std::sqrt(2.0), 1e-6);
832 }
833
834
835 //==============================================================================
836 // Main test function
837 //==============================================================================
838 static void
test_triangle_3d()839 test_triangle_3d()
840 {
841 test_closest_point();
842
843 test_non_intersecting();
844
845 test_intersecting1();
846 test_intersecting2();
847 test_intersecting3();
848 test_intersecting4();
849 test_intersecting5();
850
851 test_intersect_plane();
852
853 test_intersecting_degenerate_triangles1();
854 test_intersecting_degenerate_triangles2();
855
856 test_coincident_edges1();
857 test_coincident_edges2();
858
859 test_distance();
860 test_area();
861 test_aspect_ratio();
862
863 test_point_containment();
864 #ifdef PERFORMANCE_TEST // Performance Test disabled by default - KOM
865 test_point_containment_algo_perf();
866 #endif
867 }
868
869 TESTMAIN(test_triangle_3d);
870