1 // Copyright (c) 1998-2004
2 // Utrecht University (The Netherlands),
3 // ETH Zurich (Switzerland),
4 // INRIA Sophia-Antipolis (France),
5 // Max-Planck-Institute Saarbruecken (Germany),
6 // and Tel-Aviv University (Israel). All rights reserved.
7 //
8 // This file is part of CGAL (www.cgal.org)
9 //
10 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Distance_2/include/CGAL/squared_distance_2_2.h $
11 // $Id: squared_distance_2_2.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot
12 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
13 //
14 //
15 // Author(s) : Geert-Jan Giezeman
16
17
18 #ifndef CGAL_SQUARED_DISTANCE_2_2_H
19 #define CGAL_SQUARED_DISTANCE_2_2_H
20
21 #include <CGAL/user_classes.h>
22
23 #include <CGAL/kernel_assertions.h>
24 #include <CGAL/Point_2.h>
25 #include <CGAL/Segment_2.h>
26 #include <CGAL/Line_2.h>
27 #include <CGAL/Ray_2.h>
28 #include <CGAL/Triangle_2.h>
29 #include <CGAL/enum.h>
30 #include <CGAL/wmult.h>
31 #include <CGAL/squared_distance_utils.h>
32 #include <CGAL/squared_distance_2_1.h>
33
34 namespace CGAL {
35
36 namespace internal {
37
38 template <class K>
39 void
distance_index(int & ind1,int & ind2,const typename K::Point_2 & pt,const typename K::Triangle_2 & triangle,const K & k)40 distance_index(int &ind1,
41 int &ind2,
42 const typename K::Point_2 &pt,
43 const typename K::Triangle_2 &triangle,
44 const K& k )
45 {
46 typename K::Left_turn_2 leftturn = k.left_turn_2_object();
47 typedef typename K::Point_2 Point_2;
48 const Point_2 &vt0 = triangle.vertex(0);
49 const Point_2 &vt1 = triangle.vertex(1);
50 const Point_2 &vt2 = triangle.vertex(2);
51 if (leftturn(vt0, vt1, vt2)) {
52 if (leftturn(pt, vt1, vt0)) {
53 if (!is_acute_angle(vt0, vt1, pt, k)) {
54 if (leftturn(pt, vt2, vt1)) {
55 if (!is_acute_angle(vt1, vt2, pt, k)) {
56 ind1 = 2; ind2 = -1;
57 return;
58 }
59 if (!is_acute_angle(vt2, vt1, pt, k)) {
60 ind1 = 1; ind2 = -1;
61 return;
62 }
63 ind1 = 1; ind2 = 2;
64 return;
65 }
66 ind1 = 1; ind2 = -1;
67 return;
68 }
69 if (!is_acute_angle(vt1, vt0, pt, k)) {
70 if (leftturn(pt, vt0, vt2)) {
71 if (!is_acute_angle(vt0, vt2, pt, k)) {
72 ind1 = 2; ind2 = -1;
73 return;
74 }
75 if (!is_acute_angle(vt2, vt0, pt, k)) {
76 ind1 = 0; ind2 = -1;
77 return;
78 }
79 ind1 = 2; ind2 = 0;
80 return;
81 }
82 ind1 = 0; ind2 = -1;
83 return;
84 }
85 ind1 = 0; ind2 = 1;
86 return;
87 } else {
88 if (leftturn(pt, vt2, vt1)) {
89 if (!is_acute_angle(vt1, vt2, pt, k)) {
90 if (leftturn(pt, vt0, vt2)) {
91 if (!is_acute_angle(vt0, vt2, pt, k)) {
92 ind1 = 2; ind2 = -1;
93 return;
94 }
95 if (!is_acute_angle(vt2, vt0, pt, k)) {
96 ind1 = 0; ind2 = -1;
97 return;
98 }
99 ind1 = 2; ind2 = 0;
100 return;
101 }
102 ind1 = 0; ind2 = -1;
103 return;
104 }
105 if (!is_acute_angle(vt2, vt1, pt, k)) {
106 ind1 = 1; ind2 = -1;
107 return;
108 }
109 ind1 = 1; ind2 = 2;
110 return;
111 } else {
112 if (leftturn(pt, vt0, vt2)) {
113 if (!is_acute_angle(vt2, vt0, pt, k)) {
114 ind1 = 0; ind2 = -1;
115 return;
116 }
117 if (!is_acute_angle(vt0, vt2, pt, k)) {
118 ind1 = 2; ind2 = -1;
119 return;
120 }
121 ind1 = 2; ind2 = 0;
122 return;
123 } else {
124 ind1 = -1; ind2 = -1; // point inside or on boundary.
125 return;
126 }
127 }
128 }
129 } else {
130 if (leftturn(pt, vt2, vt0)) {
131 if (!is_acute_angle(vt0, vt2, pt, k)) {
132 if (leftturn(pt, vt1, vt2)) {
133 if (!is_acute_angle(vt2, vt1, pt, k)) {
134 ind1 = 1; ind2 = -1;
135 return;
136 }
137 if (!is_acute_angle(vt1, vt2, pt, k)) {
138 ind1 = 2; ind2 = -1;
139 return;
140 }
141 ind1 = 2; ind2 = 1;
142 return;
143 }
144 ind1 = 2; ind2 = -1;
145 return;
146 }
147 if (!is_acute_angle(vt2, vt0, pt, k)) {
148 if (leftturn(pt, vt0, vt1)) {
149 if (!is_acute_angle(vt0, vt1, pt, k)) {
150 ind1 = 1; ind2 = -1;
151 return;
152 }
153 if (!is_acute_angle(vt1, vt0, pt, k)) {
154 ind1 = 0; ind2 = -1;
155 return;
156 }
157 ind1 = 1; ind2 = 0;
158 return;
159 }
160 ind1 = 0; ind2 = -1;
161 return;
162 }
163 ind1 = 0; ind2 = 2;
164 return;
165 } else {
166 if (leftturn(pt, vt1, vt2)) {
167 if (!is_acute_angle(vt2, vt1, pt, k)) {
168 if (leftturn(pt, vt0, vt1)) {
169 if (!is_acute_angle(vt0, vt1, pt, k)) {
170 ind1 = 1; ind2 = -1;
171 return;
172 }
173 if (!is_acute_angle(vt1, vt0, pt, k)) {
174 ind1 = 0; ind2 = -1;
175 return;
176 }
177 ind1 = 1; ind2 = 0;
178 return;
179 }
180 ind1 = 0; ind2 = -1;
181 return;
182 }
183 if (!is_acute_angle(vt1, vt2, pt, k)) {
184 ind1 = 2; ind2 = -1;
185 return;
186 }
187 ind1 = 2; ind2 = 1;
188 return;
189 } else {
190 if (leftturn(pt, vt0, vt1)) {
191 if (!is_acute_angle(vt1, vt0, pt, k)) {
192 ind1 = 0; ind2 = -1;
193 return;
194 }
195 if (!is_acute_angle(vt0, vt1, pt, k)) {
196 ind1 = 1; ind2 = -1;
197 return;
198 }
199 ind1 = 1; ind2 = 0;
200 return;
201 } else {
202 ind1 = -1; ind2 = -1; // point inside or on boundary.
203 return;
204 }
205 }
206 }
207 }
208 }
209
210
211
212
213 template <class K>
214 typename K::FT
squared_distance_indexed(const typename K::Point_2 & pt,const typename K::Triangle_2 & triangle,int ind1,int ind2,const K & k)215 squared_distance_indexed(const typename K::Point_2 &pt,
216 const typename K::Triangle_2 &triangle,
217 int ind1, int ind2,
218 const K& k)
219 {
220 typedef typename K::FT FT;
221 typedef typename K::Line_2 Line_2;
222 if (ind1 == -1)
223 return FT(0);
224 if (ind2 == -1)
225 return internal::squared_distance(pt, triangle.vertex(ind1), k);
226 return internal::squared_distance(pt,
227 Line_2(triangle.vertex(ind1), triangle.vertex(ind2)),
228 k);
229 }
230
231
232
233 template <class K>
234 typename K::FT
squared_distance(const typename K::Point_2 & pt,const typename K::Triangle_2 & triangle,const K & k)235 squared_distance(const typename K::Point_2 &pt,
236 const typename K::Triangle_2 &triangle,
237 const K& k)
238 {
239 int ind1,ind2;
240 distance_index<K>(ind1, ind2, pt, triangle, k);
241 return squared_distance_indexed(pt, triangle, ind1, ind2, k);
242 }
243
244
245 template <class K>
246 inline typename K::FT
squared_distance(const typename K::Triangle_2 & triangle,const typename K::Point_2 & pt,const K & k)247 squared_distance(const typename K::Triangle_2 & triangle,
248 const typename K::Point_2 & pt,
249 const K& k)
250 {
251 return internal::squared_distance(pt, triangle, k);
252 }
253
254
255 template <class K>
256 typename K::FT
squared_distance(const typename K::Line_2 & line,const typename K::Triangle_2 & triangle,const K & k)257 squared_distance(const typename K::Line_2 &line,
258 const typename K::Triangle_2 &triangle,
259 const K& k)
260 {
261 typedef typename K::FT FT;
262 Oriented_side side0;
263 side0 = line.oriented_side(triangle.vertex(0));
264 if (line.oriented_side(triangle.vertex(1)) != side0)
265 return FT(0);
266 if (line.oriented_side(triangle.vertex(2)) != side0)
267 return FT(0);
268 FT mindist, dist;
269 int i;
270 mindist = internal::squared_distance(triangle.vertex(0),line,k);
271 for (i=1; i<3; i++) {
272 dist = internal::squared_distance(triangle.vertex(i),line,k);
273 if (dist < mindist)
274 mindist = dist;
275 }
276 return mindist;
277 }
278
279
280 template <class K>
281 inline typename K::FT
squared_distance(const typename K::Triangle_2 & triangle,const typename K::Line_2 & line,const K & k)282 squared_distance(const typename K::Triangle_2 & triangle,
283 const typename K::Line_2 & line,
284 const K& k)
285 {
286 return internal::squared_distance(line, triangle, k);
287 }
288
289
290 template <class K>
291 typename K::FT
squared_distance(const typename K::Ray_2 & ray,const typename K::Triangle_2 & triangle,const K & k)292 squared_distance(const typename K::Ray_2 &ray,
293 const typename K::Triangle_2 &triangle,
294 const K& k)
295 {
296 typedef typename K::FT FT;
297 typedef typename K::Point_2 Point_2;
298 typedef typename K::Line_2 Line_2;
299 int i, ind_tr1, ind_tr2, ind_ray = 0, ind1;
300 FT mindist, dist;
301 distance_index<K>(ind_tr1, ind_tr2, ray.source(), triangle, k);
302 mindist =
303 squared_distance_indexed(ray.source(), triangle, ind_tr1, ind_tr2, k);
304 for (i=0; i<3; i++) {
305 const Point_2& pt = triangle.vertex(i);
306 distance_index<K>(ind1, pt, ray, k);
307 dist = squared_distance_indexed(pt, ray, ind1, k);
308 if (dist < mindist) {
309 ind_ray = ind1;
310 ind_tr1 = i; ind_tr2 = -1;
311 mindist = dist;
312 }
313 }
314 // now check if all vertices are on the right side of the separating line.
315 // In case of vertex-vertex smallest distance this is the case.
316 if (ind_tr2 == -1 && ind_ray != -1)
317 return mindist;
318 if (ind_tr2 != -1) {
319 // Check if all the segment vertices lie at the same side of
320 // the triangle segment.
321 const Point_2 &vt1 = triangle.vertex(ind_tr1);
322 const Point_2 &vt2 = triangle.vertex(ind_tr2);
323 if (clockwise(ray.direction().vector(), vt2-vt1, k)) {
324 mindist = FT(0);
325 }
326 } else {
327 // Check if all the triangle vertices lie
328 // at the same side of the segment.
329 const Line_2 &sl = ray.supporting_line();
330 Oriented_side or_s = sl.oriented_side(triangle.vertex(0));
331 for (i=1; i<3; i++) {
332 if (sl.oriented_side(triangle.vertex(i)) != or_s) {
333 mindist = FT(0);
334 break;
335 }
336 }
337 }
338 return mindist;
339 }
340
341
342 template <class K>
343 inline typename K::FT
squared_distance(const typename K::Triangle_2 & triangle,const typename K::Ray_2 & ray,const K & k)344 squared_distance(const typename K::Triangle_2 & triangle,
345 const typename K::Ray_2 & ray,
346 const K& k)
347 {
348 return internal::squared_distance(ray, triangle, k);
349 }
350
351
352 template <class K>
353 typename K::FT
squared_distance(const typename K::Segment_2 & seg,const typename K::Triangle_2 & triangle,const K & k)354 squared_distance(const typename K::Segment_2 &seg,
355 const typename K::Triangle_2 &triangle,
356 const K& k)
357 {
358 typedef typename K::FT FT;
359 typedef typename K::Point_2 Point_2;
360 typename K::Orientation_2 orientation;
361 int i, ind_tr1 = 0, ind_tr2 = -1, ind_seg = 0, ind1, ind2;
362 FT mindist, dist;
363 mindist = internal::squared_distance(seg.source(), triangle.vertex(0), k);
364 for (i=0; i<2; i++) {
365 const Point_2 &pt = seg.vertex(i);
366 distance_index<K>(ind1, ind2, pt, triangle, k);
367 dist = internal::squared_distance_indexed(pt, triangle, ind1, ind2, k);
368 if (dist < mindist) {
369 ind_seg = i;
370 ind_tr1 = ind1; ind_tr2 = ind2;
371 mindist = dist;
372 }
373 }
374 for (i=0; i<3; i++) {
375 const Point_2& pt = triangle.vertex(i);
376 distance_index<K>(ind1, pt, seg, k);
377 dist = internal::squared_distance_indexed(pt, seg, ind1, k);
378 if (dist < mindist) {
379 ind_seg = ind1;
380 ind_tr1 = i; ind_tr2 = -1;
381 mindist = dist;
382 }
383 }
384 // now check if all vertices are on the right side of the separating line.
385 // In case of vertex-vertex smallest distance this is the case.
386 if (ind_tr2 == -1 && ind_seg != -1)
387 return mindist;
388
389 if (ind_tr2 != -1) {
390 // Check if all the segment vertices lie at the same side of
391 // the triangle segment.
392 const Point_2 &vt1 = triangle.vertex(ind_tr1);
393 const Point_2 &vt2 = triangle.vertex(ind_tr2);
394 Orientation or_s = orientation(vt1, vt2, seg.source());
395 if (orientation(vt1, vt2, seg.target()) != or_s) {
396 mindist = FT(0);
397 }
398 } else {
399 // Check if all the triangle vertices lie
400 // at the same side of the segment.
401 const Point_2 &vt1 = seg.source();
402 const Point_2 &vt2 = seg.target();
403 Orientation or_s = orientation(vt1, vt2, triangle.vertex(0));
404 for (i=1; i<3; i++) {
405 if (orientation(vt1, vt2, triangle.vertex(i)) != or_s) {
406 mindist = FT(0);
407 break;
408 }
409 }
410 }
411 return mindist;
412 }
413
414
415 template <class K>
416 inline typename K::FT
squared_distance(const typename K::Triangle_2 & triangle,const typename K::Segment_2 & seg,const K & k)417 squared_distance(const typename K::Triangle_2 & triangle,
418 const typename K::Segment_2 & seg,
419 const K& k)
420 {
421 return internal::squared_distance(seg, triangle, k);
422 }
423
424
425
426 template <class K>
427 typename K::FT
squared_distance(const typename K::Triangle_2 & triangle1,const typename K::Triangle_2 & triangle2,const K & k)428 squared_distance(const typename K::Triangle_2 &triangle1,
429 const typename K::Triangle_2 &triangle2,
430 const K& k)
431 {
432 typedef typename K::FT FT;
433 typedef typename K::Point_2 Point_2;
434 typename K::Orientation_2 orientation;
435 int i, ind1_1 = 0,ind1_2 = -1, ind2_1 = 0, ind2_2 = -1, ind1, ind2;
436 FT mindist, dist;
437 mindist =
438 internal::squared_distance(triangle1.vertex(0), triangle2.vertex(0), k);
439 for (i=0; i<3; i++) {
440 const Point_2& pt = triangle1.vertex(i);
441 distance_index<K>(ind1, ind2, pt, triangle2, k);
442 dist = squared_distance_indexed(pt, triangle2, ind1, ind2, k);
443 if (dist < mindist) {
444 ind1_1 = i; ind1_2 = -1;
445 ind2_1 = ind1; ind2_2 = ind2;
446 mindist = dist;
447 }
448 }
449 for (i=0; i<3; i++) {
450 const Point_2& pt = triangle2.vertex(i);
451 distance_index<K>(ind1, ind2, pt, triangle1, k);
452 dist = squared_distance_indexed(pt, triangle1, ind1, ind2, k);
453 if (dist < mindist) {
454 ind1_1 = ind1; ind1_2 = ind2;
455 ind2_1 = i; ind2_2 = -1;
456 mindist = dist;
457 }
458 }
459 // now check if all vertices are on the right side of the
460 // separating line.
461 if (ind1_2 == -1 && ind2_2 == -1)
462 return mindist;
463 // In case of point-segment closest distance, there is still the
464 // possibility of overlapping triangles. Check if all the
465 // vertices lie at the same side of the segment.
466 if (ind1_2 != -1) {
467 const Point_2 &vt1 = triangle1.vertex(ind1_1);
468 const Point_2 &vt2 = triangle1.vertex(ind1_2);
469 Orientation or_s = orientation(vt1, vt2, triangle2.vertex(0));
470 for (i=1; i<3; i++) {
471 if (orientation(vt1, vt2, triangle2.vertex(i)) != or_s) {
472 mindist = FT(0);
473 break;
474 }
475 }
476 } else {
477 const Point_2 &vt1 = triangle2.vertex(ind2_1);
478 const Point_2 &vt2 = triangle2.vertex(ind2_2);
479 Orientation or_s = orientation(vt1, vt2, triangle1.vertex(0));
480 for (i=1; i<3; i++) {
481 if (orientation(vt1, vt2, triangle1.vertex(i)) != or_s) {
482 mindist = FT(0);
483 break;
484 }
485 }
486 }
487 return mindist;
488 }
489
490 } // namespace internal
491
492 template <class K>
493 inline typename K::FT
squared_distance(const Point_2<K> & pt,const Triangle_2<K> & triangle)494 squared_distance(const Point_2<K> &pt,
495 const Triangle_2<K> &triangle)
496 {
497 return internal::squared_distance(pt, triangle, K());
498 }
499
500 template <class K>
501 inline typename K::FT
squared_distance(const Triangle_2<K> & triangle,const Point_2<K> & pt)502 squared_distance(const Triangle_2<K> &triangle,
503 const Point_2<K> &pt)
504 {
505 return internal::squared_distance(pt, triangle, K());
506 }
507
508 template <class K>
509 inline typename K::FT
squared_distance(const Line_2<K> & line,const Triangle_2<K> & triangle)510 squared_distance(const Line_2<K> &line,
511 const Triangle_2<K> &triangle)
512 {
513 return internal::squared_distance(line, triangle, K());
514 }
515
516 template <class K>
517 inline typename K::FT
squared_distance(const Triangle_2<K> & triangle,const Line_2<K> & line)518 squared_distance(const Triangle_2<K> &triangle,
519 const Line_2<K> &line)
520 {
521 return internal::squared_distance(line, triangle, K());
522 }
523
524 template <class K>
525 inline typename K::FT
squared_distance(const Ray_2<K> & ray,const Triangle_2<K> & triangle)526 squared_distance(const Ray_2<K> &ray,
527 const Triangle_2<K> &triangle)
528 {
529 return internal::squared_distance(ray, triangle, K());
530 }
531
532 template <class K>
533 inline typename K::FT
squared_distance(const Triangle_2<K> & triangle,const Ray_2<K> & ray)534 squared_distance(const Triangle_2<K> &triangle,
535 const Ray_2<K> &ray)
536 {
537 return internal::squared_distance(ray, triangle, K());
538 }
539
540 template <class K>
541 inline typename K::FT
squared_distance(const Segment_2<K> & seg,const Triangle_2<K> & triangle)542 squared_distance(const Segment_2<K> &seg,
543 const Triangle_2<K> &triangle)
544 {
545 return internal::squared_distance(seg, triangle, K());
546 }
547
548 template <class K>
549 inline typename K::FT
squared_distance(const Triangle_2<K> & triangle,const Segment_2<K> & seg)550 squared_distance(const Triangle_2<K> &triangle,
551 const Segment_2<K> &seg)
552 {
553 return internal::squared_distance(seg, triangle, K());
554 }
555
556 template <class K>
557 inline typename K::FT
squared_distance(const Triangle_2<K> & triangle1,const Triangle_2<K> & triangle2)558 squared_distance(const Triangle_2<K> &triangle1,
559 const Triangle_2<K> &triangle2)
560 {
561 return internal::squared_distance(triangle1, triangle2, K());
562 }
563
564 } //namespace CGAL
565
566 #endif
567