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