Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4GeomTools.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. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: $
28 //
29 //
30 // class G4GeomTools Implementation
31 //
32 // Author: evgueni.tcherniaev@cern.ch
33 //
34 // 10.10.2016 E.Tcherniaev: initial version.
35 // --------------------------------------------------------------------
36 
37 #include "G4GeomTools.hh"
38 
39 #include "geomdefs.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4GeometryTolerance.hh"
42 
44 //
45 // Calculate area of a triangle in 2D
46 
48  G4double Bx, G4double By,
49  G4double Cx, G4double Cy)
50 {
51  return ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax))*0.5;
52 }
53 
55 //
56 // Calculate area of a triangle in 2D
57 
59  const G4TwoVector& B,
60  const G4TwoVector& C)
61 {
62  G4double Ax = A.x(), Ay = A.y();
63  return ((B.x()-Ax)*(C.y()-Ay) - (B.y()-Ay)*(C.x()-Ax))*0.5;
64 }
65 
67 //
68 // Calculate area of a quadrilateral in 2D
69 
71  const G4TwoVector& B,
72  const G4TwoVector& C,
73  const G4TwoVector& D)
74 {
75  return ((C.x()-A.x())*(D.y()-B.y()) - (C.y()-A.y())*(D.x()-B.x()))*0.5;
76 }
77 
79 //
80 // Calculate area of a polygon in 2D
81 
83 {
84  G4double area = 0.0;
85  G4int n = p.size();
86  for(G4int i=0,k=n-1; i<n; k=i,++i)
87  {
88  area += p[k].x()*p[i].y() - p[i].x()*p[k].y();
89  }
90  return area*0.5;
91 }
92 
94 //
95 // Point inside 2D triangle
96 
98  G4double Bx, G4double By,
99  G4double Cx, G4double Cy,
100  G4double Px, G4double Py)
101 
102 {
103  if ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) > 0.)
104  {
105  if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) < 0.) return false;
106  if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) < 0.) return false;
107  if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) < 0.) return false;
108  }
109  else
110  {
111  if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) > 0.) return false;
112  if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) > 0.) return false;
113  if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) > 0.) return false;
114  }
115  return true;
116 }
117 
119 //
120 // Point inside 2D triangle
121 
123  const G4TwoVector& B,
124  const G4TwoVector& C,
125  const G4TwoVector& P)
126 {
127  G4double Ax = A.x(), Ay = A.y();
128  G4double Bx = B.x(), By = B.y();
129  G4double Cx = C.x(), Cy = C.y();
130  G4double Px = P.x(), Py = P.y();
131  if ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) > 0.)
132  {
133  if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) < 0.) return false;
134  if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) < 0.) return false;
135  if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) < 0.) return false;
136  }
137  else
138  {
139  if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) > 0.) return false;
140  if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) > 0.) return false;
141  if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) > 0.) return false;
142  }
143  return true;
144 }
145 
147 //
148 // Detemine whether 2D polygon is convex or not
149 
151 {
152  static const G4double kCarTolerance =
154 
155  G4bool gotNegative = false;
156  G4bool gotPositive = false;
157  G4int n = polygon.size();
158  if (n <= 0) return false;
159  for (G4int icur=0; icur<n; ++icur)
160  {
161  G4int iprev = (icur == 0) ? n-1 : icur-1;
162  G4int inext = (icur == n-1) ? 0 : icur+1;
163  G4TwoVector e1 = polygon[icur] - polygon[iprev];
164  G4TwoVector e2 = polygon[inext] - polygon[icur];
165  G4double cross = e1.x()*e2.y() - e1.y()*e2.x();
166  if (std::abs(cross) < kCarTolerance) return false;
167  if (cross < 0) gotNegative = true;
168  if (cross > 0) gotPositive = true;
169  if (gotNegative && gotPositive) return false;
170  }
171  return true;
172 }
173 
175 //
176 // Triangulate simple polygon
177 
180 {
181  result.resize(0);
182  std::vector<G4int> triangles;
183  G4bool reply = TriangulatePolygon(polygon,triangles);
184 
185  G4int n = triangles.size();
186  for (G4int i=0; i<n; ++i) result.push_back(polygon[triangles[i]]);
187  return reply;
188 }
189 
191 //
192 // Triangulation of a simple polygon by "ear clipping"
193 
195  std::vector<G4int>& result)
196 {
197  result.resize(0);
198 
199  // allocate and initialize list of Vertices in polygon
200  //
201  G4int n = polygon.size();
202  if (n < 3) return false;
203 
204  // we want a counter-clockwise polygon in V
205  //
206  G4double area = G4GeomTools::PolygonArea(polygon);
207  G4int* V = new G4int[n];
208  if (area > 0.)
209  for (G4int i=0; i<n; ++i) V[i] = i;
210  else
211  for (G4int i=0; i<n; ++i) V[i] = (n-1)-i;
212 
213  // Triangulation: remove nv-2 Vertices, creating 1 triangle every time
214  //
215  G4int nv = n;
216  G4int count = 2*nv; // error detection counter
217  for(G4int b=nv-1; nv>2; )
218  {
219  // ERROR: if we loop, it is probably a non-simple polygon
220  if ((count--) <= 0)
221  {
222  delete[] V;
223  if (area < 0.) std::reverse(result.begin(),result.end());
224  return false;
225  }
226 
227  // three consecutive vertices in current polygon, <a,b,c>
228  G4int a = (b < nv) ? b : 0; // previous
229  b = (a+1 < nv) ? a+1 : 0; // current
230  G4int c = (b+1 < nv) ? b+1 : 0; // next
231 
232  if (CheckSnip(polygon, a,b,c, nv,V))
233  {
234  // output Triangle
235  result.push_back(V[a]);
236  result.push_back(V[b]);
237  result.push_back(V[c]);
238 
239  // remove vertex b from remaining polygon
240  nv--;
241  for(G4int i=b; i<nv; ++i) V[i] = V[i+1];
242 
243  count = 2*nv; // resest error detection counter
244  }
245  }
246  delete[] V;
247  if (area < 0.) std::reverse(result.begin(),result.end());
248  return true;
249 }
250 
252 //
253 // Helper function for "ear clipping" polygon triangulation.
254 // Check for a valid snip
255 
256 G4bool G4GeomTools::CheckSnip(const G4TwoVectorList& contour,
257  G4int a, G4int b, G4int c,
258  G4int n, const G4int* V)
259 {
260  static const G4double kCarTolerance =
262 
263  // check orientation of Triangle
264  G4double Ax = contour[V[a]].x(), Ay = contour[V[a]].y();
265  G4double Bx = contour[V[b]].x(), By = contour[V[b]].y();
266  G4double Cx = contour[V[c]].x(), Cy = contour[V[c]].y();
267  if ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) < kCarTolerance) return false;
268 
269  // check that there is no point inside Triangle
270  G4double xmin = std::min(std::min(Ax,Bx),Cx);
271  G4double xmax = std::max(std::max(Ax,Bx),Cx);
272  G4double ymin = std::min(std::min(Ay,By),Cy);
273  G4double ymax = std::max(std::max(Ay,By),Cy);
274  for (G4int i=0; i<n; ++i)
275  {
276  if((i == a) || (i == b) || (i == c)) continue;
277  G4double Px = contour[V[i]].x();
278  if (Px < xmin || Px > xmax) continue;
279  G4double Py = contour[V[i]].y();
280  if (Py < ymin || Py > ymax) continue;
281  if (PointInTriangle(Ax,Ay,Bx,By,Cx,Cy,Px,Py)) return false;
282  }
283  return true;
284 }
285 
286 
287 
288 
290 //
291 // Remove collinear and coincident points from 2D polygon
292 
294  std::vector<G4int>& iout,
295  G4double tolerance)
296 {
297  iout.resize(0);
298  // set tolerance squared
299  G4double delta = tolerance*tolerance;
300  // set special value to mark vertices for removal
301  G4double removeIt = kInfinity;
302 
303  G4int nv = polygon.size();
304 
305  // Main loop: check every three consecutive points, if the points
306  // are collinear then mark middle point for removal
307  //
308  G4int icur = 0, iprev = 0, inext = 0, nout = 0;
309  for (G4int i=0; i<nv; ++i)
310  {
311  icur = i; // index of current point
312 
313  for (G4int k=1; k<nv+1; ++k) // set index of previous point
314  {
315  iprev = icur - k;
316  if (iprev < 0) iprev += nv;
317  if (polygon[iprev].x() != removeIt) break;
318  }
319 
320  for (G4int k=1; k<nv+1; ++k) // set index of next point
321  {
322  inext = icur + k;
323  if (inext >= nv) inext -= nv;
324  if (polygon[inext].x() != removeIt) break;
325  }
326 
327  if (iprev == inext) break; // degenerate polygon, stop
328 
329  // Calculate parameters of triangle (iprev->icur->inext),
330  // if triangle is too small or too narrow then mark current
331  // point for removal
332  G4TwoVector e1 = polygon[iprev] - polygon[icur];
333  G4TwoVector e2 = polygon[inext] - polygon[icur];
334 
335  // Check length of edges, then check height of the triangle
336  G4double leng1 = e1.mag2();
337  G4double leng2 = e2.mag2();
338  G4double leng3 = (e2-e1).mag2();
339  if (leng1 <= delta || leng2 <= delta || leng3 <= delta)
340  {
341  polygon[icur].setX(removeIt); nout++;
342  }
343  else
344  {
345  G4double lmax = std::max(std::max(leng1,leng2),leng3);
346  G4double area = std::abs(e1.x()*e2.y()-e1.y()*e2.x())*0.5;
347  if (area/std::sqrt(lmax) <= std::abs(tolerance))
348  {
349  polygon[icur].setX(removeIt); nout++;
350  }
351  }
352  }
353 
354  // Remove marked points
355  //
356  icur = 0;
357  if (nv - nout < 3) // degenerate polygon, remove all points
358  {
359  for (G4int i=0; i<nv; ++i) iout.push_back(i);
360  polygon.resize(0);
361  nv = 0;
362  }
363  for (G4int i=0; i<nv; ++i) // move points, if required
364  {
365  if (polygon[i].x() != removeIt)
366  polygon[icur++] = polygon[i];
367  else
368  iout.push_back(i);
369  }
370  if (icur < nv) polygon.resize(icur);
371  return;
372 }
373 
375 //
376 // Find bounding box of a disk sector
377 
379  G4double startPhi, G4double delPhi,
380  G4TwoVector& pmin, G4TwoVector& pmax)
381 {
382  static const G4double kCarTolerance =
384 
385  // check parameters
386  //
387  pmin.set(0,0);
388  pmax.set(0,0);
389  if (rmin < 0) return false;
390  if (rmax <= rmin + kCarTolerance) return false;
391  if (delPhi <= 0 + kCarTolerance) return false;
392 
393  // calculate extent
394  //
395  pmin.set(-rmax,-rmax);
396  pmax.set( rmax, rmax);
397  if (delPhi >= CLHEP::twopi) return true;
398 
399  DiskExtent(rmin,rmax,
400  std::sin(startPhi),std::cos(startPhi),
401  std::sin(startPhi+delPhi),std::cos(startPhi+delPhi),
402  pmin,pmax);
403  return true;
404 }
405 
407 //
408 // Find bounding box of a disk sector, fast version.
409 // No check of parameters !!!
410 
412  G4double sinStart, G4double cosStart,
413  G4double sinEnd, G4double cosEnd,
414  G4TwoVector& pmin, G4TwoVector& pmax)
415 {
416  static const G4double kCarTolerance =
418 
419  // check if 360 degrees
420  //
421  pmin.set(-rmax,-rmax);
422  pmax.set( rmax, rmax);
423 
424  if (std::abs(sinEnd-sinStart) < kCarTolerance &&
425  std::abs(cosEnd-cosStart) < kCarTolerance) return;
426 
427  // get start and end quadrants
428  //
429  // 1 | 0
430  // ---+---
431  // 3 | 2
432  //
433  G4int icase = (cosEnd < 0) ? 1 : 0;
434  if (sinEnd < 0) icase += 2;
435  if (cosStart < 0) icase += 4;
436  if (sinStart < 0) icase += 8;
437 
438  switch (icase)
439  {
440  // start quadrant 0
441  case 0: // start->end : 0->0
442  if (sinEnd < sinStart) break;
443  pmin.set(rmin*cosEnd,rmin*sinStart);
444  pmax.set(rmax*cosStart,rmax*sinEnd );
445  break;
446  case 1: // start->end : 0->1
447  pmin.set(rmax*cosEnd,std::min(rmin*sinStart,rmin*sinEnd));
448  pmax.set(rmax*cosStart,rmax );
449  break;
450  case 2: // start->end : 0->2
451  pmin.set(-rmax,-rmax);
452  pmax.set(std::max(rmax*cosStart,rmax*cosEnd),rmax);
453  break;
454  case 3: // start->end : 0->3
455  pmin.set(-rmax,rmax*sinEnd);
456  pmax.set(rmax*cosStart,rmax);
457  break;
458  // start quadrant 1
459  case 4: // start->end : 1->0
460  pmin.set(-rmax,-rmax);
461  pmax.set(rmax,std::max(rmax*sinStart,rmax*sinEnd));
462  break;
463  case 5: // start->end : 1->1
464  if (sinEnd > sinStart) break;
465  pmin.set(rmax*cosEnd,rmin*sinEnd );
466  pmax.set(rmin*cosStart,rmax*sinStart);
467  break;
468  case 6: // start->end : 1->2
469  pmin.set(-rmax,-rmax);
470  pmax.set(rmax*cosEnd,rmax*sinStart);
471  break;
472  case 7: // start->end : 1->3
473  pmin.set(-rmax,rmax*sinEnd);
474  pmax.set(std::max(rmin*cosStart,rmin*cosEnd),rmax*sinStart);
475  break;
476  // start quadrant 2
477  case 8: // start->end : 2->0
478  pmin.set(std::min(rmin*cosStart,rmin*cosEnd),rmax*sinStart);
479  pmax.set(rmax,rmax*sinEnd);
480  break;
481  case 9: // start->end : 2->1
482  pmin.set(rmax*cosEnd,rmax*sinStart);
483  pmax.set(rmax,rmax);
484  break;
485  case 10: // start->end : 2->2
486  if (sinEnd < sinStart) break;
487  pmin.set(rmin*cosStart,rmax*sinStart);
488  pmax.set(rmax*cosEnd,rmin*sinEnd );
489  break;
490  case 11: // start->end : 2->3
491  pmin.set(-rmax,std::min(rmax*sinStart,rmax*sinEnd));
492  pmax.set(rmax,rmax);
493  break;
494  // start quadrant 3
495  case 12: // start->end : 3->0
496  pmin.set(rmax*cosStart,-rmax);
497  pmax.set(rmax,rmax*sinEnd);
498  break;
499  case 13: // start->end : 3->1
500  pmin.set(std::min(rmax*cosStart,rmax*cosEnd),-rmax);
501  pmax.set(rmax,rmax);
502  break;
503  case 14: // start->end : 3->2
504  pmin.set(rmax*cosStart,-rmax);
505  pmax.set(rmax*cosEnd,std::max(rmin*sinStart,rmin*sinEnd));
506  break;
507  case 15: // start->end : 3->3
508  if (sinEnd > sinStart) break;
509  pmin.set(rmax*cosStart,rmax*sinEnd);
510  pmax.set(rmin*cosEnd,rmin*sinStart);
511  break;
512  }
513  return;
514 }
515 
517 //
518 // Calculate distance between point P and line segment AB in 3D
519 
523 {
524  G4ThreeVector AP = P - A;
525  G4ThreeVector AB = B - A;
526 
527  G4double u = AP.dot(AB);
528  if (u <= 0) return AP.mag(); // closest point is A
529 
530  G4double len2 = AB.mag2();
531  if (u >= len2) return (B-P).mag(); // closest point is B
532 
533  return ((u/len2)*AB - AP).mag(); // distance to line
534 }
535 
536 
538 //
539 // Calculate bounding box of a spherical sector
540 
541 G4bool
543  G4double startTheta, G4double delTheta,
544  G4double startPhi, G4double delPhi,
545  G4ThreeVector& pmin, G4ThreeVector& pmax)
546 {
547  static const G4double kCarTolerance =
549 
550  // check parameters
551  //
552  pmin.set(0,0,0);
553  pmax.set(0,0,0);
554  if (rmin < 0) return false;
555  if (rmax <= rmin + kCarTolerance) return false;
556  if (delTheta <= 0 + kCarTolerance) return false;
557  if (delPhi <= 0 + kCarTolerance) return false;
558 
559  G4double stheta = startTheta;
560  G4double dtheta = delTheta;
561  if (stheta < 0 && stheta > CLHEP::pi) return false;
562  if (stheta + dtheta > CLHEP::pi) dtheta = CLHEP::pi - stheta;
563  if (dtheta <= 0 + kCarTolerance) return false;
564 
565  // calculate extent
566  //
567  pmin.set(-rmax,-rmax,-rmax);
568  pmax.set( rmax, rmax, rmax);
569  if (dtheta >= CLHEP::pi && delPhi >= CLHEP::twopi) return true;
570 
571  G4double etheta = stheta + dtheta;
572  G4double sinStart = std::sin(stheta);
573  G4double cosStart = std::cos(stheta);
574  G4double sinEnd = std::sin(etheta);
575  G4double cosEnd = std::cos(etheta);
576 
577  G4double rhomin = rmin*std::min(sinStart,sinEnd);
578  G4double rhomax = rmax;
579  if (stheta > CLHEP::halfpi) rhomax = rmax*sinStart;
580  if (etheta < CLHEP::halfpi) rhomax = rmax*sinEnd;
581 
582  G4TwoVector xymin,xymax;
583  DiskExtent(rhomin,rhomax,
584  std::sin(startPhi),std::cos(startPhi),
585  std::sin(startPhi+delPhi),std::cos(startPhi+delPhi),
586  xymin,xymax);
587 
588  G4double zmin = std::min(rmin*cosEnd,rmax*cosEnd);
589  G4double zmax = std::max(rmin*cosStart,rmax*cosStart);
590  pmin.set(xymin.x(),xymin.y(),zmin);
591  pmax.set(xymax.x(),xymax.y(),zmax);
592  return true;
593 }
G4double G4ParticleHPJENDLHEData::G4double result
void set(double x, double y, double z)
double y() const
double x() const
static const G4double kInfinity
Definition: geomdefs.hh:42
double dot(const Hep3Vector &) const
static G4double QuadArea(const G4TwoVector &A, const G4TwoVector &B, const G4TwoVector &C, const G4TwoVector &D)
Definition: G4GeomTools.cc:70
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
std::vector< G4TwoVector > G4TwoVectorList
Definition: G4GeomTools.hh:50
const char * p
Definition: xmltok.h:285
G4double GetSurfaceTolerance() const
static G4bool PointInTriangle(G4double Px, G4double Py, G4double Ax, G4double Ay, G4double Bx, G4double By, G4double Cx, G4double Cy)
Definition: G4GeomTools.cc:97
static G4double PolygonArea(const G4TwoVectorList &polygon)
Definition: G4GeomTools.cc:82
static G4double DistancePointSegment(G4ThreeVector P, G4ThreeVector A, G4ThreeVector B)
Definition: G4GeomTools.cc:520
tuple x
Definition: test.py:50
double C(double temp)
double B(double temperature)
static G4bool SphereExtent(G4double rmin, G4double rmax, G4double startTheta, G4double delTheta, G4double startPhi, G4double delPhi, G4ThreeVector &pmin, G4ThreeVector &pmax)
Definition: G4GeomTools.cc:542
int G4int
Definition: G4Types.hh:78
static double P[]
tuple b
Definition: test.py:12
double mag2() const
double A(double temperature)
void set(double x, double y)
bool G4bool
Definition: G4Types.hh:79
const G4double kCarTolerance
const G4int n
static G4bool TriangulatePolygon(const G4TwoVectorList &polygon, G4TwoVectorList &result)
Definition: G4GeomTools.cc:178
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
double mag2() const
double D(double temp)
static G4bool IsConvex(const G4TwoVectorList &polygon)
Definition: G4GeomTools.cc:150
static constexpr double halfpi
Definition: SystemOfUnits.h:56
static G4double TriangleArea(G4double Ax, G4double Ay, G4double Bx, G4double By, G4double Cx, G4double Cy)
Definition: G4GeomTools.cc:47
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13
double mag() const
static constexpr double twopi
Definition: SystemOfUnits.h:55
static constexpr double pi
Definition: SystemOfUnits.h:54
static G4bool DiskExtent(G4double rmin, G4double rmax, G4double startPhi, G4double delPhi, G4TwoVector &pmin, G4TwoVector &pmax)
Definition: G4GeomTools.cc:378
static G4GeometryTolerance * GetInstance()
static void RemoveRedundantVertices(G4TwoVectorList &polygon, std::vector< G4int > &iout, G4double tolerance=0)
Definition: G4GeomTools.cc:293