GteDistPointTriangle.h
Go to the documentation of this file.
1 // David Eberly, Geometric Tools, Redmond WA 98052
2 // Copyright (c) 1998-2017
3 // Distributed under the Boost Software License, Version 1.0.
4 // http://www.boost.org/LICENSE_1_0.txt
5 // http://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
6 // File Version: 3.0.0 (2016/06/19)
7 
8 #pragma once
9 
10 #include <Mathematics/GteVector.h>
11 #include <Mathematics/GteDCPQuery.h>
12 #include <Mathematics/GteTriangle.h>
13 
14 namespace gte
15 {
16 
17 template <int N, typename Real>
18 class DCPQuery<Real, Vector<N, Real>, Triangle<N, Real>>
19 {
20 public:
21  struct Result
22  {
23  Real distance, sqrDistance;
24  Real parameter[3]; // barycentric coordinates for triangle.v[3]
25  Vector<N, Real> closest;
26  };
27 
28  Result operator()(Vector<N, Real> const& point,
29  Triangle<N, Real> const& triangle);
30 
31 private:
32  inline void GetMinEdge02(Real const& a11, Real const& b1,
33  Vector<2, Real>& p) const;
34 
35  inline void GetMinEdge12(Real const& a01, Real const& a11, Real const& b1,
36  Real const& f10, Real const& f01, Vector<2, Real>& p) const;
37 
38  inline void GetMinInterior(Vector<2, Real> const& p0, Real const& h0,
39  Vector<2, Real> const& p1, Real const& h1, Vector<2, Real>& p) const;
40 };
41 
42 // Template aliases for convenience.
43 template <int N, typename Real>
44 using DCPPointTriangle =
45 DCPQuery<Real, Vector<N, Real>, Triangle<N, Real>>;
46 
47 template <typename Real>
48 using DCPPoint2Triangle2 = DCPPointTriangle<2, Real>;
49 
50 template <typename Real>
51 using DCPPoint3Triangle3 = DCPPointTriangle<3, Real>;
52 
53 
54 template <int N, typename Real> inline
55 void DCPQuery<Real, Vector<N, Real>, Triangle<N, Real>>::GetMinEdge02(
56  Real const& a11, Real const& b1, Vector<2, Real>& p) const
57 {
58  p[0] = (Real)0;
59  if (b1 >= (Real)0)
60  {
61  p[1] = (Real)0;
62  }
63  else if (a11 + b1 <= (Real)0)
64  {
65  p[1] = (Real)1;
66  }
67  else
68  {
69  p[1] = -b1 / a11;
70  }
71 }
72 
73 template <int N, typename Real> inline
74 void DCPQuery<Real, Vector<N, Real>, Triangle<N, Real>>::GetMinEdge12(
75  Real const& a01, Real const& a11, Real const& b1, Real const& f10,
76  Real const& f01, Vector<2, Real>& p) const
77 {
78  Real h0 = a01 + b1 - f10;
79  if (h0 >= (Real)0)
80  {
81  p[1] = (Real)0;
82  }
83  else
84  {
85  Real h1 = a11 + b1 - f01;
86  if (h1 <= (Real)0)
87  {
88  p[1] = (Real)1;
89  }
90  else
91  {
92  p[1] = h0 / (h0 - h1);
93  }
94  }
95  p[0] = (Real)1 - p[1];
96 }
97 
98 template <int N, typename Real> inline
99 void DCPQuery<Real, Vector<N, Real>, Triangle<N, Real>>::GetMinInterior(
100  Vector<2, Real> const& p0, Real const& h0, Vector<2, Real> const& p1,
101  Real const& h1, Vector<2, Real>& p) const
102 {
103  Real z = h0 / (h0 - h1);
104  p = ((Real)1 - z) * p0 + z * p1;
105 }
106 
107 template <int N,typename Real>
108 typename DCPQuery<Real, Vector<N, Real>, Triangle<N, Real>>::Result
109 DCPQuery<Real, Vector<N, Real>, Triangle<N, Real>>::operator()(
110  Vector<N, Real> const& point, Triangle<N, Real> const& triangle)
111 {
112  Vector<N, Real> diff = point - triangle.v[0];
113  Vector<N, Real> edge0 = triangle.v[1] - triangle.v[0];
114  Vector<N, Real> edge1 = triangle.v[2] - triangle.v[0];
115  Real a00 = Dot(edge0, edge0);
116  Real a01 = Dot(edge0, edge1);
117  Real a11 = Dot(edge1, edge1);
118  Real b0 = -Dot(diff, edge0);
119  Real b1 = -Dot(diff, edge1);
120 
121  Real f00 = b0;
122  Real f10 = b0 + a00;
123  Real f01 = b0 + a01;
124 
125  Vector<2, Real> p0, p1, p;
126  Real dt1, h0, h1;
127 
128  // Compute the endpoints p0 and p1 of the segment. The segment is
129  // parameterized by L(z) = (1-z)*p0 + z*p1 for z in [0,1] and the
130  // directional derivative of half the quadratic on the segment is
131  // H(z) = Dot(p1-p0,gradient[Q](L(z))/2), where gradient[Q]/2 = (F,G).
132  // By design, F(L(z)) = 0 for cases (2), (4), (5), and (6). Cases (1) and
133  // (3) can correspond to no-intersection or intersection of F = 0 with the
134  // triangle.
135  if (f00 >= (Real)0)
136  {
137  if (f01 >= (Real)0)
138  {
139  // (1) p0 = (0,0), p1 = (0,1), H(z) = G(L(z))
140  GetMinEdge02(a11, b1, p);
141  }
142  else
143  {
144  // (2) p0 = (0,t10), p1 = (t01,1-t01), H(z) = (t11 - t10)*G(L(z))
145  p0[0] = (Real)0;
146  p0[1] = f00 / (f00 - f01);
147  p1[0] = f01 / (f01 - f10);
148  p1[1] = (Real)1 - p1[0];
149  dt1 = p1[1] - p0[1];
150  h0 = dt1 * (a11 * p0[1] + b1);
151  if (h0 >= (Real)0)
152  {
153  GetMinEdge02(a11, b1, p);
154  }
155  else
156  {
157  h1 = dt1 * (a01 * p1[0] + a11 * p1[1] + b1);
158  if (h1 <= (Real)0)
159  {
160  GetMinEdge12(a01, a11, b1, f10, f01, p);
161  }
162  else
163  {
164  GetMinInterior(p0, h0, p1, h1, p);
165  }
166  }
167  }
168  }
169  else if (f01 <= (Real)0)
170  {
171  if (f10 <= (Real)0)
172  {
173  // (3) p0 = (1,0), p1 = (0,1), H(z) = G(L(z)) - F(L(z))
174  GetMinEdge12(a01, a11, b1, f10, f01, p);
175  }
176  else
177  {
178  // (4) p0 = (t00,0), p1 = (t01,1-t01), H(z) = t11*G(L(z))
179  p0[0] = f00 / (f00 - f10);
180  p0[1] = (Real)0;
181  p1[0] = f01 / (f01 - f10);
182  p1[1] = (Real)1 - p1[0];
183  h0 = p1[1] * (a01 * p0[0] + b1);
184  if (h0 >= (Real)0)
185  {
186  p = p0; // GetMinEdge01
187  }
188  else
189  {
190  h1 = p1[1] * (a01 * p1[0] + a11 * p1[1] + b1);
191  if (h1 <= (Real)0)
192  {
193  GetMinEdge12(a01, a11, b1, f10, f01, p);
194  }
195  else
196  {
197  GetMinInterior(p0, h0, p1, h1, p);
198  }
199  }
200  }
201  }
202  else if (f10 <= (Real)0)
203  {
204  // (5) p0 = (0,t10), p1 = (t01,1-t01), H(z) = (t11 - t10)*G(L(z))
205  p0[0] = (Real)0;
206  p0[1] = f00 / (f00 - f01);
207  p1[0] = f01 / (f01 - f10);
208  p1[1] = (Real)1 - p1[0];
209  dt1 = p1[1] - p0[1];
210  h0 = dt1 * (a11 * p0[1] + b1);
211  if (h0 >= (Real)0)
212  {
213  GetMinEdge02(a11, b1, p);
214  }
215  else
216  {
217  h1 = dt1 * (a01 * p1[0] + a11 * p1[1] + b1);
218  if (h1 <= (Real)0)
219  {
220  GetMinEdge12(a01, a11, b1, f10, f01, p);
221  }
222  else
223  {
224  GetMinInterior(p0, h0, p1, h1, p);
225  }
226  }
227  }
228  else
229  {
230  // (6) p0 = (t00,0), p1 = (0,t11), H(z) = t11*G(L(z))
231  p0[0] = f00 / (f00 - f10);
232  p0[1] = (Real)0;
233  p1[0] = (Real)0;
234  p1[1] = f00 / (f00 - f01);
235  h0 = p1[1] * (a01 * p0[0] + b1);
236  if (h0 >= (Real)0)
237  {
238  p = p0; // GetMinEdge01
239  }
240  else
241  {
242  h1 = p1[1] * (a11 * p1[1] + b1);
243  if (h1 <= (Real)0)
244  {
245  GetMinEdge02(a11, b1, p);
246  }
247  else
248  {
249  GetMinInterior(p0, h0, p1, h1, p);
250  }
251  }
252  }
253 
254  Result result;
255  result.parameter[0] = (Real)1 - p[0] - p[1];
256  result.parameter[1] = p[0];
257  result.parameter[2] = p[1];
258  result.closest = triangle.v[0] + p[0] * edge0 + p[1] * edge1;
259  diff = point - result.closest;
260  result.sqrDistance = Dot(diff, diff);
261  result.distance = sqrt(result.sqrDistance);
262  return result;
263 }
264 
265 
266 }
distance
GLsizei GLsizei GLfloat distance
Definition: glext.h:9704
gte::DCPQuery::operator()
Result operator()(Type0 const &primitive0, Type1 const &primitive1)
gte::Dot
DualQuaternion< Real > Dot(DualQuaternion< Real > const &d0, DualQuaternion< Real > const &d1)
z
GLdouble GLdouble GLdouble z
Definition: glcorearb.h:843
p
GLfloat GLfloat p
Definition: glext.h:11668
result
GLuint64EXT * result
Definition: glext.h:10003

geometric_tools_engine
Author(s): Yijiang Huang
autogenerated on Thu Jul 18 2019 03:59:59