1 /**********************************************************************
2 *
3 * GEOS - Geometry Engine Open Source
4 * http://geos.osgeo.org
5 *
6 * Copyright (C) 2012 Excensus LLC.
7 *
8 * This is free software; you can redistribute and/or modify it under
9 * the terms of the GNU Lesser General Licence as published
10 * by the Free Software Foundation.
11 * See the COPYING file for more information.
12 *
13 **********************************************************************
14 *
15 * Last port: triangulate/quadedge/TrianglePredicate.java r524
16 *
17 **********************************************************************/
18
19 #include <geos/triangulate/quadedge/TrianglePredicate.h>
20
21 #include <geos/geom/Coordinate.h>
22
23 namespace geos {
24 namespace geom { // geos.geom
25
26 bool
isInCircleNonRobust(const Coordinate & a,const Coordinate & b,const Coordinate & c,const Coordinate & p)27 TrianglePredicate::isInCircleNonRobust(
28 const Coordinate& a, const Coordinate& b, const Coordinate& c,
29 const Coordinate& p)
30 {
31 bool isInCircle =
32 (a.x * a.x + a.y * a.y) * triArea(b, c, p)
33 - (b.x * b.x + b.y * b.y) * triArea(a, c, p)
34 + (c.x * c.x + c.y * c.y) * triArea(a, b, p)
35 - (p.x * p.x + p.y * p.y) * triArea(a, b, c)
36 > 0;
37 return isInCircle;
38 }
39
40 bool
isInCircleNormalized(const Coordinate & a,const Coordinate & b,const Coordinate & c,const Coordinate & p)41 TrianglePredicate::isInCircleNormalized(
42 const Coordinate& a, const Coordinate& b, const Coordinate& c,
43 const Coordinate& p)
44 {
45 // Unfortunately this implementation is not robust either. For robust one see:
46 // https://www.cs.cmu.edu/~quake/robust.html
47 // https://www.cs.cmu.edu/afs/cs/project/quake/public/code/predicates.c
48
49 long double adx = a.x - p.x;
50 long double ady = a.y - p.y;
51 long double bdx = b.x - p.x;
52 long double bdy = b.y - p.y;
53 long double cdx = c.x - p.x;
54 long double cdy = c.y - p.y;
55
56 long double bdxcdy = bdx * cdy;
57 long double cdxbdy = cdx * bdy;
58 long double alift = adx * adx + ady * ady;
59
60 long double cdxady = cdx * ady;
61 long double adxcdy = adx * cdy;
62 long double blift = bdx * bdx + bdy * bdy;
63
64 long double adxbdy = adx * bdy;
65 long double bdxady = bdx * ady;
66 long double clift = cdx * cdx + cdy * cdy;
67 return (alift * bdxcdy + blift * cdxady + clift * adxbdy) >
68 (alift * cdxbdy + blift * adxcdy + clift * bdxady);
69 }
70
71 double
triArea(const Coordinate & a,const Coordinate & b,const Coordinate & c)72 TrianglePredicate::triArea(const Coordinate& a,
73 const Coordinate& b, const Coordinate& c)
74 {
75 return (b.x - a.x) * (c.y - a.y)
76 - (b.y - a.y) * (c.x - a.x);
77 }
78
79 bool
isInCircleRobust(const Coordinate & a,const Coordinate & b,const Coordinate & c,const Coordinate & p)80 TrianglePredicate::isInCircleRobust(
81 const Coordinate& a, const Coordinate& b, const Coordinate& c,
82 const Coordinate& p)
83 {
84 // This implementation is not robust, name is ported from JTS.
85 return isInCircleNormalized(a, b, c, p);
86 }
87
88 } // namespace geos.geom
89 } // namespace geos
90