Geant4  10.03
G4TessellatedGeometryAlgorithms.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration and of QinetiQ Ltd, *
20 // * subject to DEFCON 705 IPR conditions. *
21 // * By using, copying, modifying or distributing the software (or *
22 // * any work based on the software) you agree to acknowledge its *
23 // * use in resulting scientific publications, and indicate your *
24 // * acceptance of all terms of the Geant4 Software license. *
25 // ********************************************************************
26 //
27 //
28 // $Id: G4TessellatedGeometryAlgorithms.cc 66356 2012-12-18 09:02:32Z gcosmo $
29 //
30 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
31 //
32 // CHANGE HISTORY
33 // --------------
34 //
35 // 07 August 2007, P R Truscott, QinetiQ Ltd, UK - Created, with member
36 // functions based on the work of Rickard Holmberg.
37 //
38 // 26 September 2007
39 // P R Truscott, qinetiQ Ltd, UK
40 // Updated to assign values of location array, not update
41 // just the pointer.
42 //
43 // 12 October 2012, M Gayer, CERN, - Reviewed optimized implementation.
44 //
46 
48 
49 #include <cfloat>
50 
52 //
53 // IntersectLineAndTriangle2D
54 //
55 // Determines whether there is an intersection between a line defined
56 // by r = p + s.v and a triangle defined by verticies p0, p0+e0 and p0+e1.
57 //
58 // Here:
59 // p = 2D vector
60 // s = scaler on [0,infinity)
61 // v = 2D vector
62 // p0, e0 and e1 are 2D vectors
63 // Information about where the intersection occurs is returned in the
64 // variable location.
65 //
66 // This is based on the work of Rickard Holmberg.
67 //
69  const G4TwoVector &p, const G4TwoVector &v,
70  const G4TwoVector &p0, const G4TwoVector &e0, const G4TwoVector &e1,
71  G4TwoVector location[2])
72 {
73  G4TwoVector loc0[2];
74  G4int e0i = IntersectLineAndLineSegment2D (p,v,p0,e0,loc0);
75  if (e0i == 2)
76  {
77  location[0] = loc0[0];
78  location[1] = loc0[1];
79  return true;
80  }
81 
82  G4TwoVector loc1[2];
83  G4int e1i = IntersectLineAndLineSegment2D (p,v,p0,e1,loc1);
84  if (e1i == 2)
85  {
86  location[0] = loc1[0];
87  location[1] = loc1[1];
88  return true;
89  }
90 
91  if ((e0i == 1) && (e1i == 1))
92  {
93  if ((loc0[0]-p).mag2() < (loc1[0]-p).mag2())
94  {
95  location[0] = loc0[0];
96  location[1] = loc1[0];
97  }
98  else
99  {
100  location[0] = loc1[0];
101  location[1] = loc0[0];
102  }
103  return true;
104  }
105 
106  G4TwoVector p1 = p0 + e0;
107  G4TwoVector DE = e1 - e0;
108  G4TwoVector loc2[2];
109  G4int e2i = IntersectLineAndLineSegment2D (p,v,p1,DE,loc2);
110  if (e2i == 2)
111  {
112  location[0] = loc2[0];
113  location[1] = loc2[1];
114  return true;
115  }
116 
117  if ((e0i == 0) && (e1i == 0) && (e2i == 0)) return false;
118 
119  if ((e0i == 1) && (e2i == 1))
120  {
121  if ((loc0[0]-p).mag2() < (loc2[0]-p).mag2())
122  {
123  location[0] = loc0[0];
124  location[1] = loc2[0];
125  }
126  else
127  {
128  location[0] = loc2[0];
129  location[1] = loc0[0];
130  }
131  return true;
132  }
133 
134  if ((e1i == 1) && (e2i == 1))
135  {
136  if ((loc1[0]-p).mag2() < (loc2[0]-p).mag2())
137  {
138  location[0] = loc1[0];
139  location[1] = loc2[0];
140  }
141  else
142  {
143  location[0] = loc2[0];
144  location[1] = loc1[0];
145  }
146  return true;
147  }
148 
149  return false;
150 }
151 
153 //
154 // IntersectLineAndLineSegment2D
155 //
156 // Determines whether there is an intersection between a line defined
157 // by r = p0 + s.d0 and a line-segment with endpoints p1 and p1+d1.
158 // Here:
159 // p0 = 2D vector
160 // s = scaler on [0,infinity)
161 // d0 = 2D vector
162 // p1 and d1 are 2D vectors
163 //
164 // This function returns:
165 // 0 - if there is no intersection;
166 // 1 - if there is a unique intersection;
167 // 2 - if the line and line-segments overlap, and the intersection is a
168 // segment itself.
169 // Information about where the intersection occurs is returned in the
170 // as ??.
171 //
172 // This is based on the work of Rickard Holmberg as well as material published
173 // by Philip J Schneider and David H Eberly, "Geometric Tools for Computer
174 // Graphics," ISBN 1-55860-694-0, pp 244-245, 2003.
175 //
177  const G4TwoVector &p0, const G4TwoVector &d0,
178  const G4TwoVector &p1, const G4TwoVector &d1,
179  G4TwoVector location[2])
180 {
181  G4TwoVector e = p1 - p0;
182  G4double kross = cross(d0,d1);
183  G4double sqrKross = kross * kross;
184  G4double sqrLen0 = d0.mag2();
185  G4double sqrLen1 = d1.mag2();
186  location[0] = G4TwoVector(0.0,0.0);
187  location[1] = G4TwoVector(0.0,0.0);
188 
189  if (sqrKross > DBL_EPSILON * DBL_EPSILON * sqrLen0 * sqrLen1)
190  {
191  //
192  // The line and line segment are not parallel. Determine if the intersection
193  // is in positive s where r=p0 + s*d0, and for 0<=t<=1 where r=p1 + t*d1.
194  //
195  G4double ss = cross(e,d1)/kross;
196  if (ss < 0) return 0; // Intersection does not occur for positive ss
197  G4double t = cross(e,d0)/kross;
198  if (t < 0 || t > 1) return 0; // Intersection does not occur on line-segment
199  //
200  // Intersection of lines is a single point on the forward-propagating line
201  // defined by r=p0 + ss*d0, and the line segment defined by r=p1 + t*d1.
202  //
203  location[0] = p0 + ss*d0;
204  return 1;
205  }
206  //
207  // Line and line segment are parallel. Determine whether they overlap or not.
208  //
209  G4double sqrLenE = e.mag2();
210  kross = cross(e,d0);
211  sqrKross = kross * kross;
212  if (sqrKross > DBL_EPSILON * DBL_EPSILON * sqrLen0 * sqrLenE)
213  {
214  return 0; //Lines are different.
215  }
216  //
217  // Lines are the same. Test for overlap.
218  //
219  G4double s0 = d0.dot(e)/sqrLen0;
220  G4double s1 = s0 + d0.dot(d1)/sqrLen0;
221  G4double smin = 0.0;
222  G4double smax = 0.0;
223 
224  if (s0 < s1) {smin = s0; smax = s1;}
225  else {smin = s1; smax = s0;}
226 
227  if (smax < 0.0) return 0;
228  else if (smin < 0.0)
229  {
230  location[0] = p0;
231  location[1] = p0 + smax*d0;
232  return 2;
233  }
234  else
235  {
236  location[0] = p0 + smin*d0;
237  location[1] = p0 + smax*d0;
238  return 2;
239  }
240 }
241 
243 //
244 // CrossProduct
245 //
246 // This is just a ficticious "cross-product" function for two 2D vectors...
247 // "ficticious" because such an operation is not relevant to 2D space compared
248 // with 3D space.
249 //
251  const G4TwoVector &v2)
252 {
253  return v1.x()*v2.y() - v1.y()*v2.x();
254 }
static G4double cross(const G4TwoVector &v1, const G4TwoVector &v2)
int G4int
Definition: G4Types.hh:78
const G4int smax
bool G4bool
Definition: G4Types.hh:79
#define DBL_EPSILON
Definition: templates.hh:87
static const G4double d1
static G4bool IntersectLineAndTriangle2D(const G4TwoVector &p, const G4TwoVector &v, const G4TwoVector &p0, const G4TwoVector &e0, const G4TwoVector &e1, G4TwoVector location[2])
CLHEP::Hep2Vector G4TwoVector
Definition: G4TwoVector.hh:42
double G4double
Definition: G4Types.hh:76
static G4int IntersectLineAndLineSegment2D(const G4TwoVector &p0, const G4TwoVector &d0, const G4TwoVector &p1, const G4TwoVector &d1, G4TwoVector location[2])