Geant4  10.01
G4ClippablePolygon.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: G4ClippablePolygon.cc 66356 2012-12-18 09:02:32Z gcosmo $
28 //
29 //
30 // --------------------------------------------------------------------
31 // GEANT 4 class source file
32 //
33 //
34 // G4ClippablePolygon.cc
35 //
36 // Includes code from G4VSolid (P.Kent, V.Grichine, J.Allison)
37 //
38 // --------------------------------------------------------------------
39 
40 #include "G4ClippablePolygon.hh"
41 
42 #include "G4VoxelLimits.hh"
43 #include "G4GeometryTolerance.hh"
44 
45 //
46 // Constructor
47 //
49  : normal(0.,0.,0.)
50 {
52 }
53 
54 
55 //
56 // Destructor
57 //
59 {
60 }
61 
62 
63 //
64 // AddVertexInOrder
65 //
67 {
68  vertices.push_back( vertex );
69 }
70 
71 
72 //
73 // ClearAllVertices
74 //
76 {
77  vertices.clear();
78 }
79 
80 
81 //
82 // Clip
83 //
85 {
86  if (voxelLimit.IsLimited()) {
87  ClipAlongOneAxis( voxelLimit, kXAxis );
88  ClipAlongOneAxis( voxelLimit, kYAxis );
89  ClipAlongOneAxis( voxelLimit, kZAxis );
90  }
91 
92  return (vertices.size() > 0);
93 }
94 
95 
96 //
97 // PartialClip
98 //
99 // Clip, while ignoring the indicated axis
100 //
102  const EAxis IgnoreMe )
103 {
104  if (voxelLimit.IsLimited()) {
105  if (IgnoreMe != kXAxis) ClipAlongOneAxis( voxelLimit, kXAxis );
106  if (IgnoreMe != kYAxis) ClipAlongOneAxis( voxelLimit, kYAxis );
107  if (IgnoreMe != kZAxis) ClipAlongOneAxis( voxelLimit, kZAxis );
108  }
109 
110  return (vertices.size() > 0);
111 }
112 
113 
114 //
115 // GetExtent
116 //
118  G4double &min,
119  G4double &max ) const
120 {
121  //
122  // Okay, how many entries do we have?
123  //
124  G4int noLeft = vertices.size();
125 
126  //
127  // Return false if nothing is left
128  //
129  if (noLeft == 0) return false;
130 
131  //
132  // Initialize min and max to our first vertex
133  //
134  min = max = vertices[0].operator()( axis );
135 
136  //
137  // Compare to the rest
138  //
139  G4int i;
140  for( i=1; i<noLeft; i++ )
141  {
142  G4double component = vertices[i].operator()( axis );
143  if (component < min )
144  min = component;
145  else if (component > max )
146  max = component;
147  }
148 
149  return true;
150 }
151 
152 
153 //
154 // GetMinPoint
155 //
156 // Returns pointer to minimum point along the specified axis.
157 // Take care! Do not use pointer after destroying parent polygon.
158 //
160 {
161  G4int noLeft = vertices.size();
162  if (noLeft==0)
163  G4Exception("G4ClippablePolygon::GetMinPoint()",
164  "GeomSolids0002", FatalException, "Empty polygon.");
165 
166  const G4ThreeVector *answer = &(vertices[0]);
167  G4double min = answer->operator()(axis);
168 
169  G4int i;
170  for( i=1; i<noLeft; i++ )
171  {
172  G4double component = vertices[i].operator()( axis );
173  if (component < min)
174  {
175  answer = &(vertices[i]);
176  min = component;
177  }
178  }
179 
180  return answer;
181 }
182 
183 
184 //
185 // GetMaxPoint
186 //
187 // Returns pointer to maximum point along the specified axis.
188 // Take care! Do not use pointer after destroying parent polygon.
189 //
191 {
192  G4int noLeft = vertices.size();
193  if (noLeft==0)
194  G4Exception("G4ClippablePolygon::GetMaxPoint()",
195  "GeomSolids0002", FatalException, "Empty polygon.");
196 
197  const G4ThreeVector *answer = &(vertices[0]);
198  G4double max = answer->operator()(axis);
199 
200  G4int i;
201  for( i=1; i<noLeft; i++ )
202  {
203  G4double component = vertices[i].operator()( axis );
204  if (component > max)
205  {
206  answer = &(vertices[i]);
207  max = component;
208  }
209  }
210 
211  return answer;
212 }
213 
214 
215 //
216 // InFrontOf
217 //
218 // Decide if this polygon is in "front" of another when
219 // viewed along the specified axis. For our purposes here,
220 // it is sufficient to use the minimum extent of the
221 // polygon along the axis to determine this.
222 //
223 // In case the minima of the two polygons are equal,
224 // we use a more sophisticated test.
225 //
226 // Note that it is possible for the two following
227 // statements to both return true or both return false:
228 // polygon1.InFrontOf(polygon2)
229 // polygon2.BehindOf(polygon1)
230 //
232  EAxis axis ) const
233 {
234  //
235  // If things are empty, do something semi-sensible
236  //
237  G4int noLeft = vertices.size();
238  if (noLeft==0) return false;
239 
240  if (other.Empty()) return true;
241 
242  //
243  // Get minimum of other polygon
244  //
245  const G4ThreeVector *minPointOther = other.GetMinPoint( axis );
246  const G4double minOther = minPointOther->operator()(axis);
247 
248  //
249  // Get minimum of this polygon
250  //
251  const G4ThreeVector *minPoint = GetMinPoint( axis );
252  const G4double min = minPoint->operator()(axis);
253 
254  //
255  // Easy decision
256  //
257  if (min < minOther-kCarTolerance) return true; // Clear winner
258 
259  if (minOther < min-kCarTolerance) return false; // Clear loser
260 
261  //
262  // We have a tie (this will not be all that rare since our
263  // polygons are connected)
264  //
265  // Check to see if there is a vertex in the other polygon
266  // that is behind this one (or vice versa)
267  //
268  G4bool answer;
269  G4ThreeVector normalOther = other.GetNormal();
270 
271  if (std::fabs(normalOther(axis)) > std::fabs(normal(axis)))
272  {
273  G4double minP, maxP;
274  GetPlanerExtent( *minPointOther, normalOther, minP, maxP );
275 
276  answer = (normalOther(axis) > 0) ? (minP < -kCarTolerance)
277  : (maxP > +kCarTolerance);
278  }
279  else
280  {
281  G4double minP, maxP;
282  other.GetPlanerExtent( *minPoint, normal, minP, maxP );
283 
284  answer = (normal(axis) > 0) ? (maxP > +kCarTolerance)
285  : (minP < -kCarTolerance);
286  }
287  return answer;
288 }
289 
290 //
291 // BehindOf
292 //
293 // Decide if this polygon is behind another.
294 // See notes in method "InFrontOf"
295 //
297  EAxis axis ) const
298 {
299  //
300  // If things are empty, do something semi-sensible
301  //
302  G4int noLeft = vertices.size();
303  if (noLeft==0) return false;
304 
305  if (other.Empty()) return true;
306 
307  //
308  // Get minimum of other polygon
309  //
310  const G4ThreeVector *maxPointOther = other.GetMaxPoint( axis );
311  const G4double maxOther = maxPointOther->operator()(axis);
312 
313  //
314  // Get minimum of this polygon
315  //
316  const G4ThreeVector *maxPoint = GetMaxPoint( axis );
317  const G4double max = maxPoint->operator()(axis);
318 
319  //
320  // Easy decision
321  //
322  if (max > maxOther+kCarTolerance) return true; // Clear winner
323 
324  if (maxOther > max+kCarTolerance) return false; // Clear loser
325 
326  //
327  // We have a tie (this will not be all that rare since our
328  // polygons are connected)
329  //
330  // Check to see if there is a vertex in the other polygon
331  // that is in front of this one (or vice versa)
332  //
333  G4bool answer;
334  G4ThreeVector normalOther = other.GetNormal();
335 
336  if (std::fabs(normalOther(axis)) > std::fabs(normal(axis)))
337  {
338  G4double minP, maxP;
339  GetPlanerExtent( *maxPointOther, normalOther, minP, maxP );
340 
341  answer = (normalOther(axis) > 0) ? (maxP > +kCarTolerance)
342  : (minP < -kCarTolerance);
343  }
344  else
345  {
346  G4double minP, maxP;
347  other.GetPlanerExtent( *maxPoint, normal, minP, maxP );
348 
349  answer = (normal(axis) > 0) ? (minP < -kCarTolerance)
350  : (maxP > +kCarTolerance);
351  }
352  return answer;
353 }
354 
355 
356 //
357 // GetPlanerExtent
358 //
359 // Get min/max distance in or out of a plane
360 //
362  const G4ThreeVector &planeNormal,
363  G4double &min,
364  G4double &max ) const
365 {
366  //
367  // Okay, how many entries do we have?
368  //
369  G4int noLeft = vertices.size();
370 
371  //
372  // Return false if nothing is left
373  //
374  if (noLeft == 0) return false;
375 
376  //
377  // Initialize min and max to our first vertex
378  //
379  min = max = planeNormal.dot(vertices[0]-pointOnPlane);
380 
381  //
382  // Compare to the rest
383  //
384  G4int i;
385  for( i=1; i<noLeft; i++ )
386  {
387  G4double component = planeNormal.dot(vertices[i] - pointOnPlane);
388  if (component < min )
389  min = component;
390  else if (component > max )
391  max = component;
392  }
393 
394  return true;
395 }
396 
397 
398 //
399 // Clip along just one axis, as specified in voxelLimit
400 //
402  const EAxis axis )
403 {
404  if (!voxelLimit.IsLimited(axis)) return;
405 
406  G4ThreeVectorList tempPolygon;
407 
408  //
409  // Build a "simple" voxelLimit that includes only the min extent
410  // and apply this to our vertices, producing result in tempPolygon
411  //
412  G4VoxelLimits simpleLimit1;
413  simpleLimit1.AddLimit( axis, voxelLimit.GetMinExtent(axis), kInfinity );
414  ClipToSimpleLimits( vertices, tempPolygon, simpleLimit1 );
415 
416  //
417  // If nothing is left from the above clip, we might as well return now
418  // (but with an empty vertices)
419  //
420  if (tempPolygon.size() == 0)
421  {
422  vertices.clear();
423  return;
424  }
425 
426  //
427  // Now do the same, but using a "simple" limit that includes only the max
428  // extent. Apply this to out tempPolygon, producing result in vertices.
429  //
430  G4VoxelLimits simpleLimit2;
431  simpleLimit2.AddLimit( axis, -kInfinity, voxelLimit.GetMaxExtent(axis) );
432  ClipToSimpleLimits( tempPolygon, vertices, simpleLimit2 );
433 
434  //
435  // If nothing is left, return now
436  //
437  if (vertices.size() == 0) return;
438 }
439 
440 
441 //
442 // pVoxelLimits must be only limited along one axis, and either the maximum
443 // along the axis must be +kInfinity, or the minimum -kInfinity
444 //
446  G4ThreeVectorList& outputPolygon,
447  const G4VoxelLimits& pVoxelLimit )
448 {
449  G4int i;
450  G4int noVertices=pPolygon.size();
451  G4ThreeVector vEnd,vStart;
452 
453  outputPolygon.clear();
454 
455  for (i=0;i<noVertices;i++)
456  {
457  vStart=pPolygon[i];
458  if (i==noVertices-1)
459  {
460  vEnd=pPolygon[0];
461  }
462  else
463  {
464  vEnd=pPolygon[i+1];
465  }
466 
467  if (pVoxelLimit.Inside(vStart))
468  {
469  if (pVoxelLimit.Inside(vEnd))
470  {
471  // vStart and vEnd inside -> output end point
472  //
473  outputPolygon.push_back(vEnd);
474  }
475  else
476  {
477  // vStart inside, vEnd outside -> output crossing point
478  //
479  pVoxelLimit.ClipToLimits(vStart,vEnd);
480  outputPolygon.push_back(vEnd);
481  }
482  }
483  else
484  {
485  if (pVoxelLimit.Inside(vEnd))
486  {
487  // vStart outside, vEnd inside -> output inside section
488  //
489  pVoxelLimit.ClipToLimits(vStart,vEnd);
490  outputPolygon.push_back(vStart);
491  outputPolygon.push_back(vEnd);
492  }
493  else // Both point outside -> no output
494  {
495  }
496  }
497  }
498 }
virtual G4bool InFrontOf(const G4ClippablePolygon &other, EAxis axis) const
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector G4ThreeVector
virtual void ClipAlongOneAxis(const G4VoxelLimits &voxelLimit, const EAxis axis)
virtual G4bool GetPlanerExtent(const G4ThreeVector &pointOnPlane, const G4ThreeVector &planeNormal, G4double &min, G4double &max) const
G4double GetSurfaceTolerance() const
virtual G4bool PartialClip(const G4VoxelLimits &voxelLimit, const EAxis IgnoreMe)
virtual void AddVertexInOrder(const G4ThreeVector vertex)
virtual void ClearAllVertices()
virtual const G4ThreeVector * GetMinPoint(const EAxis axis) const
int G4int
Definition: G4Types.hh:78
virtual G4bool BehindOf(const G4ClippablePolygon &other, EAxis axis) const
const G4ThreeVector GetNormal() const
virtual G4bool GetExtent(const EAxis axis, G4double &min, G4double &max) const
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
virtual const G4ThreeVector * GetMaxPoint(const EAxis axis) const
virtual G4bool Clip(const G4VoxelLimits &voxelLimit)
void AddLimit(const EAxis pAxis, const G4double pMin, const G4double pMax)
std::vector< G4ThreeVector > G4ThreeVectorList
G4bool IsLimited() const
G4ThreeVectorList vertices
bool G4bool
Definition: G4Types.hh:79
G4bool Inside(const G4ThreeVector &pVec) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
T max(const T t1, const T t2)
brief Return the largest of the two arguments
EAxis
Definition: geomdefs.hh:54
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void ClipToSimpleLimits(G4ThreeVectorList &pPolygon, G4ThreeVectorList &outputPolygon, const G4VoxelLimits &pVoxelLimit)
G4bool Empty() const
double G4double
Definition: G4Types.hh:76
G4double GetMaxExtent(const EAxis pAxis) const
G4double GetMinExtent(const EAxis pAxis) const
static G4GeometryTolerance * GetInstance()
G4bool ClipToLimits(G4ThreeVector &pStart, G4ThreeVector &pEnd) const