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