Geant4  10.01.p03
UTessellatedGeometryAlgorithms.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * This Software is part of the AIDA Unified Solids Library package *
4 // * See: https://aidasoft.web.cern.ch/USolids *
5 // ********************************************************************
6 //
7 // $Id:$
8 //
9 // --------------------------------------------------------------------
10 //
11 // UTessellatedGeometryAlgorithms
12 //
13 // 11.07.12 Marek Gayer
14 // Created from original implementation in Geant4
15 // --------------------------------------------------------------------
16 
18 
19 #include <cfloat>
20 
22 //
23 // IntersectLineAndTriangle2D
24 //
25 // Determines whether there is an intersection between a line defined
26 // by r = p + s.v and a triangle defined by verticies p0, p0+e0 and p0+e1.
27 //
28 // Here:
29 // p = 2D vector
30 // s = scaler on [0,infinity)
31 // v = 2D vector
32 // p0, e0 and e1 are 2D vectors
33 // Information about where the intersection occurs is returned in the
34 // variable location.
35 //
36 // This is based on the work of Rickard Holmberg.
37 //
39  const UVector2& p, const UVector2& v,
40  const UVector2& p0, const UVector2& e0, const UVector2& e1,
41  UVector2 location[2])
42 {
43  UVector2 loc0[2];
44  int e0i = IntersectLineAndLineSegment2D(p, v, p0, e0, loc0);
45  if (e0i == 2)
46  {
47  location[0] = loc0[0];
48  location[1] = loc0[1];
49  return true;
50  }
51 
52  UVector2 loc1[2];
53  int e1i = IntersectLineAndLineSegment2D(p, v, p0, e1, loc1);
54  if (e1i == 2)
55  {
56  location[0] = loc1[0];
57  location[1] = loc1[1];
58  return true;
59  }
60 
61  if ((e0i == 1) && (e1i == 1))
62  {
63  if ((loc0[0] - p).mag2() < (loc1[0] - p).mag2())
64  {
65  location[0] = loc0[0];
66  location[1] = loc1[0];
67  }
68  else
69  {
70  location[0] = loc1[0];
71  location[1] = loc0[0];
72  }
73  return true;
74  }
75 
76  UVector2 p1 = p0 + e0;
77  UVector2 DE = e1 - e0;
78  UVector2 loc2[2];
79  int e2i = IntersectLineAndLineSegment2D(p, v, p1, DE, loc2);
80  if (e2i == 2)
81  {
82  location[0] = loc2[0];
83  location[1] = loc2[1];
84  return true;
85  }
86 
87  if ((e0i == 0) && (e1i == 0) && (e2i == 0)) return false;
88 
89  if ((e0i == 1) && (e2i == 1))
90  {
91  if ((loc0[0] - p).mag2() < (loc2[0] - p).mag2())
92  {
93  location[0] = loc0[0];
94  location[1] = loc2[0];
95  }
96  else
97  {
98  location[0] = loc2[0];
99  location[1] = loc0[0];
100  }
101  return true;
102  }
103 
104  if ((e1i == 1) && (e2i == 1))
105  {
106  if ((loc1[0] - p).mag2() < (loc2[0] - p).mag2())
107  {
108  location[0] = loc1[0];
109  location[1] = loc2[0];
110  }
111  else
112  {
113  location[0] = loc2[0];
114  location[1] = loc1[0];
115  }
116  return true;
117  }
118 
119  return false;
120 }
121 
123 //
124 // IntersectLineAndLineSegment2D
125 //
126 // Determines whether there is an intersection between a line defined
127 // by r = p0 + s.d0 and a line-segment with endpoints p1 and p1+d1.
128 // Here:
129 // p0 = 2D vector
130 // s = scaler on [0,infinity)
131 // d0 = 2D vector
132 // p1 and d1 are 2D vectors
133 //
134 // This function returns:
135 // 0 - if there is no intersection;
136 // 1 - if there is a unique intersection;
137 // 2 - if the line and line-segments overlap, and the intersection is a
138 // segment itself.
139 // Information about where the intersection occurs is returned in the
140 // as ??.
141 //
142 // This is based on the work of Rickard Holmberg as well as material published
143 // by Philip J Schneider and David H Eberly, "Geometric Tools for Computer
144 // Graphics," ISBN 1-55860-694-0, pp 244-245, 2003.
145 //
147  const UVector2& p0, const UVector2& d0,
148  const UVector2& p1, const UVector2& d1,
149  UVector2 location[2])
150 {
151  UVector2 e = p1 - p0;
152  double kross = Cross(d0, d1);
153  double sqrKross = kross * kross;
154  double sqrLen0 = d0.mag2();
155  double sqrLen1 = d1.mag2();
156  location[0] = UVector2(0.0, 0.0);
157  location[1] = UVector2(0.0, 0.0);
158 
159  if (sqrKross > DBL_EPSILON * DBL_EPSILON * sqrLen0 * sqrLen1)
160  {
161  //
162  //
163  // The line and line segment are not parallel. Determine if the intersection
164  // is in positive s where r = p0 + s*d0, and for 0<=t<=1 where r = p1 + t*d1.
165  //
166  double ss = Cross(e, d1) / kross;
167  if (ss < 0) return 0; // Intersection does not occur for positive ss.
168  double t = Cross(e, d0) / kross;
169  if (t < 0 || t > 1) return 0; // Intersection does not occur on line-segment.
170  //
171  //
172  // Intersection of lines is a single point on the forward-propagating line
173  // defined by r = p0 + ss*d0, and the line segment defined by r = p1 + t*d1.
174  //
175  location[0] = p0 + ss * d0;
176  return 1;
177  }
178  //
179  //
180  // Line and line segment are parallel. Determine whether they overlap or not.
181  //
182  double sqrLenE = e.mag2();
183  kross = Cross(e, d0);
184  sqrKross = kross * kross;
185  if (sqrKross > DBL_EPSILON * DBL_EPSILON * sqrLen0 * sqrLenE)
186  {
187  return 0; //Lines are different.
188  }
189  //
190  //
191  // Lines are the same. Test for overlap.
192  //
193  double s0 = d0.dot(e) / sqrLen0;
194  double s1 = s0 + d0.dot(d1) / sqrLen0;
195  double smin = 0.0;
196  double smax = 0.0;
197 
198  if (s0 < s1)
199  {
200  smin = s0;
201  smax = s1;
202  }
203  else
204  {
205  smin = s1;
206  smax = s0;
207  }
208 
209  if (smax < 0.0) return 0;
210  else if (smin < 0.0)
211  {
212  location[0] = p0;
213  location[1] = p0 + smax * d0;
214  return 2;
215  }
216  else
217  {
218  location[0] = p0 + smin * d0;
219  location[1] = p0 + smax * d0;
220  return 2;
221  }
222 }
223 
225  const UVector2& v2)
226 {
227  return v1.x * v2.y - v1.y * v2.x;
228 }
static const G4double d1
static bool IntersectLineAndTriangle2D(const UVector2 &p, const UVector2 &v, const UVector2 &p0, const UVector2 &e0, const UVector2 &e1, UVector2 location[2])
static double Cross(const UVector2 &v1, const UVector2 &v2)
const G4int smax
const G4double p1
static int IntersectLineAndLineSegment2D(const UVector2 &p0, const UVector2 &d0, const UVector2 &p1, const UVector2 &d1, UVector2 location[2])
#define DBL_EPSILON
Definition: templates.hh:87
static const G4double e1
const G4double p0
double x
Definition: UVector2.hh:195
double mag2() const
Definition: UVector2.hh:305
double dot(const UVector2 &p) const
Definition: UVector2.hh:300
double y
Definition: UVector2.hh:196