| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4ReduciblePolygon implementation; a utilit
27 // test, reduce, and/or otherwise manipulate a
28 // See G4ReduciblePolygon.hh for more info.
29 //
30 // Author: David C. Williams (davidw@scipp.ucs
31 // -------------------------------------------
32
33 #include "G4ReduciblePolygon.hh"
34 #include "globals.hh"
35
36 // Constructor: with simple arrays
37 //
38 G4ReduciblePolygon::G4ReduciblePolygon( const
39 const
40
41 : aMin(0.), aMax(0.), bMin(0.), bMax(0.)
42 {
43 //
44 // Do all of the real work in Create
45 //
46 Create( a, b, n );
47 }
48
49 // Constructor: special PGON/PCON case
50 //
51 G4ReduciblePolygon::G4ReduciblePolygon( const
52 const
53 const
54 : aMin(0.), aMax(0.), bMin(0.), bMax(0.)
55 {
56 //
57 // Translate
58 //
59 auto a = new G4double[n*2];
60 auto b = new G4double[n*2];
61
62 G4double *rOut = a + n,
63 *zOut = b + n,
64 *rIn = rOut-1,
65 *zIn = zOut-1;
66
67 for( G4int i=0; i < n; ++i, ++rOut, ++zOut,
68 {
69 *rOut = rmax[i];
70 *rIn = rmin[i];
71 *zOut = *zIn = z[i];
72 }
73
74 Create( a, b, n*2 );
75
76 delete [] a;
77 delete [] b;
78 }
79
80 // Create
81 //
82 // To be called by constructors, fill in the l
83 // polygon
84 //
85 void G4ReduciblePolygon::Create( const G4doubl
86 const G4doubl
87 {
88 if (n<3)
89 G4Exception("G4ReduciblePolygon::Create()",
90 FatalErrorInArgument, "Less tha
91
92 const G4double *anext = a, *bnext = b;
93 ABVertex* prev = nullptr;
94 do // Loop checking, 13.08.2015, G.Cosmo
95 {
96 auto newVertex = new ABVertex;
97 newVertex->a = *anext;
98 newVertex->b = *bnext;
99 newVertex->next = nullptr;
100 if (prev==nullptr)
101 {
102 vertexHead = newVertex;
103 }
104 else
105 {
106 prev->next = newVertex;
107 }
108
109 prev = newVertex;
110 } while( ++anext, ++bnext < b+n );
111
112 numVertices = n;
113
114 CalculateMaxMin();
115 }
116
117 // Fake default constructor - sets only member
118 // for usage restri
119 //
120 G4ReduciblePolygon::G4ReduciblePolygon( __void
121 : aMin(0.), aMax(0.), bMin(0.), bMax(0.)
122 {
123 }
124
125
126 //
127 // Destructor
128 //
129 G4ReduciblePolygon::~G4ReduciblePolygon()
130 {
131 ABVertex* curr = vertexHead;
132 while( curr != nullptr ) // Loop checking
133 {
134 ABVertex* toDelete = curr;
135 curr = curr->next;
136 delete toDelete;
137 }
138 }
139
140 // CopyVertices
141 //
142 // Copy contents into simple linear arrays.
143 // ***** CAUTION ***** Be care to declare the
144 // enough size!
145 //
146 void G4ReduciblePolygon::CopyVertices( G4doubl
147 {
148 G4double *anext = a, *bnext = b;
149 ABVertex *curr = vertexHead;
150 while( curr != nullptr ) // Loop checking
151 {
152 *anext++ = curr->a;
153 *bnext++ = curr->b;
154 curr = curr->next;
155 }
156 }
157
158 // ScaleA
159 //
160 // Multiply all a values by a common scale
161 //
162 void G4ReduciblePolygon::ScaleA( G4double scal
163 {
164 ABVertex* curr = vertexHead;
165 while( curr != nullptr ) // Loop checking
166 {
167 curr->a *= scale;
168 curr = curr->next;
169 }
170 }
171
172 // ScaleB
173 //
174 // Multiply all b values by a common scale
175 //
176 void G4ReduciblePolygon::ScaleB( G4double scal
177 {
178 ABVertex* curr = vertexHead;
179 while( curr != nullptr ) // Loop checking
180 {
181 curr->b *= scale;
182 curr = curr->next;
183 }
184 }
185
186 // RemoveDuplicateVertices
187 //
188 // Remove adjacent vertices that are equal. Re
189 // is a problem (too few vertices remaining).
190 //
191 G4bool G4ReduciblePolygon::RemoveDuplicateVert
192 {
193 ABVertex *curr = vertexHead,
194 *prev = nullptr, *next = nullptr;
195 while( curr != nullptr ) // Loop checking
196 {
197 next = curr->next;
198 if (next == nullptr) next = vertexHead;
199
200 if (std::fabs(curr->a-next->a) < tolerance
201 std::fabs(curr->b-next->b) < tolerance
202 {
203 //
204 // Duplicate found: do we have > 3 verti
205 //
206 if (numVertices <= 3)
207 {
208 CalculateMaxMin();
209 return false;
210 }
211
212 //
213 // Delete
214 //
215 ABVertex* toDelete = curr;
216 curr = curr->next;
217 delete toDelete;
218
219 numVertices--;
220
221 if (prev != nullptr)
222 prev->next = curr;
223 else
224 vertexHead = curr;
225 }
226 else
227 {
228 prev = curr;
229 curr = curr->next;
230 }
231 }
232
233 //
234 // In principle, this is not needed, but why
235 //
236 CalculateMaxMin();
237
238 return true;
239 }
240
241 // RemoveRedundantVertices
242 //
243 // Remove any unneeded vertices, i.e. those ve
244 // are on the line connecting the previous and
245 //
246 G4bool G4ReduciblePolygon::RemoveRedundantVert
247 {
248 //
249 // Under these circumstances, we can quit no
250 //
251 if (numVertices <= 2) return false;
252
253 G4double tolerance2 = tolerance*tolerance;
254
255 //
256 // Loop over all vertices
257 //
258 ABVertex *curr = vertexHead, *next = nullptr
259 while( curr != nullptr ) // Loop checking
260 {
261 next = curr->next;
262 if (next == nullptr) next = vertexHead;
263
264 G4double da = next->a - curr->a,
265 db = next->b - curr->b;
266
267 //
268 // Loop over all subsequent vertices, up t
269 //
270 for(;;)
271 {
272 //
273 // Get vertex after next
274 //
275 ABVertex* test = next->next;
276 if (test == nullptr) test = vertexHead;
277
278 //
279 // If we are back to the original vertex
280 //
281 if (test==curr) break;
282
283 //
284 // Test for parallel line segments
285 //
286 G4double dat = test->a - curr->a,
287 dbt = test->b - curr->b;
288
289 if (std::fabs(dat*db-dbt*da)>tolerance2)
290
291 //
292 // Redundant vertex found: do we have >
293 //
294 if (numVertices <= 3)
295 {
296 CalculateMaxMin();
297 return false;
298 }
299
300 //
301 // Delete vertex pointed to by next. Car
302 //
303 if (curr->next != nullptr)
304 { // next is not head
305 if (next->next != nullptr)
306 curr->next = test; // next is not t
307 else
308 curr->next = nullptr; // New tail
309 }
310 else
311 vertexHead = test; // New head
312
313 if ((curr != next) && (next != test)) de
314
315 --numVertices;
316
317 //
318 // Replace next by the vertex we just te
319 // and keep on going...
320 //
321 next = test;
322 da = dat; db = dbt;
323 }
324 curr = curr->next;
325 }
326
327 //
328 // In principle, this is not needed, but why
329 //
330 CalculateMaxMin();
331
332 return true;
333 }
334
335 // ReverseOrder
336 //
337 // Reverse the order of the vertices
338 //
339 void G4ReduciblePolygon::ReverseOrder()
340 {
341 //
342 // Loop over all vertices
343 //
344 ABVertex* prev = vertexHead;
345 if (prev==nullptr) return; // No vertices
346
347 ABVertex* curr = prev->next;
348 if (curr==nullptr) return; // Just one ve
349
350 //
351 // Our new tail
352 //
353 vertexHead->next = nullptr;
354
355 for(;;)
356 {
357 //
358 // Save pointer to next vertex (in origina
359 //
360 ABVertex *save = curr->next;
361
362 //
363 // Replace it with a pointer to the previo
364 // (in original order)
365 //
366 curr->next = prev;
367
368 //
369 // Last vertex?
370 //
371 if (save == nullptr) break;
372
373 //
374 // Next vertex
375 //
376 prev = curr;
377 curr = save;
378 }
379
380 //
381 // Our new head
382 //
383 vertexHead = curr;
384 }
385
386
387 // StartWithZMin
388 //
389 // Starting alway with Zmin=bMin
390 // This method is used for GenericPolycone
391 //
392 void G4ReduciblePolygon::StartWithZMin()
393 {
394 ABVertex* curr = vertexHead;
395 G4double bcurr = curr->b;
396 ABVertex* prev = curr;
397 while( curr != nullptr) // Loop checking,
398 {
399 if(curr->b < bcurr)
400 {
401 bcurr = curr->b;
402 ABVertex* curr1 = curr;
403 while( curr1 != nullptr ) // Loop che
404 {
405 if(curr1->next == nullptr) { curr1->ne
406 curr1 = curr1->next;
407 }
408 vertexHead = curr;
409 prev->next = nullptr;
410 }
411 prev = curr;
412 curr = curr->next;
413 }
414 }
415
416 // CrossesItself
417 //
418 // Return "true" if the polygon crosses itself
419 //
420 // Warning: this routine is not very fast (run
421 //
422 G4bool G4ReduciblePolygon::CrossesItself( G4do
423 {
424 G4double tolerance2 = tolerance*tolerance;
425 G4double one = 1.0-tolerance,
426 zero = tolerance;
427 //
428 // Top loop over line segments. By the time
429 // with the second to last segment, we're do
430 //
431 ABVertex *curr1 = vertexHead, *next1 = nullp
432 while (curr1->next != nullptr) // Loop ch
433 {
434 next1 = curr1->next;
435 G4double da1 = next1->a-curr1->a,
436 db1 = next1->b-curr1->b;
437
438 //
439 // Inner loop over subsequent line segment
440 //
441 ABVertex* curr2 = next1->next;
442 while( curr2 != nullptr ) // Loop check
443 {
444 ABVertex* next2 = curr2->next;
445 if (next2==nullptr) next2 = vertexHead;
446 G4double da2 = next2->a-curr2->a,
447 db2 = next2->b-curr2->b;
448 G4double a12 = curr2->a-curr1->a,
449 b12 = curr2->b-curr1->b;
450
451 //
452 // Calculate intersection of the two lin
453 //
454 G4double deter = da1*db2 - db1*da2;
455 if (std::fabs(deter) > tolerance2)
456 {
457 G4double s1, s2;
458 s1 = (a12*db2-b12*da2)/deter;
459
460 if (s1 >= zero && s1 < one)
461 {
462 s2 = -(da1*b12-db1*a12)/deter;
463 if (s2 >= zero && s2 < one) return t
464 }
465 }
466 curr2 = curr2->next;
467 }
468 curr1 = next1;
469 }
470 return false;
471 }
472
473 // BisectedBy
474 //
475 // Decide if a line through two points crosses
476 //
477 G4bool G4ReduciblePolygon::BisectedBy( G4doubl
478 G4doubl
479 G4doubl
480 {
481 G4int nNeg = 0, nPos = 0;
482
483 G4double a12 = a2-a1, b12 = b2-b1;
484 G4double len12 = std::sqrt( a12*a12 + b12*b1
485 a12 /= len12; b12 /= len12;
486
487 ABVertex* curr = vertexHead;
488 do // Loop checking, 13.08.2015, G.Cosmo
489 {
490 G4double av = curr->a - a1,
491 bv = curr->b - b1;
492
493 G4double cross = av*b12 - bv*a12;
494
495 if (cross < -tolerance)
496 {
497 if (nPos != 0) return true;
498 ++nNeg;
499 }
500 else if (cross > tolerance)
501 {
502 if (nNeg != 0) return true;
503 ++nPos;
504 }
505 curr = curr->next;
506 } while( curr != nullptr );
507
508 return false;
509 }
510
511 // Area
512 //
513 // Calculated signed polygon area, where polyg
514 // clockwise manner (where x==a, y==b) have ne
515 //
516 // References: [O' Rourke (C)] pp. 18-27; [
517 // "The Area of a Simple Polygon", Jon Rokn
518 //
519 G4double G4ReduciblePolygon::Area()
520 {
521 G4double answer = 0;
522
523 ABVertex *curr = vertexHead, *next = nullptr
524 do // Loop checking, 13.08.2015, G.Cosmo
525 {
526 next = curr->next;
527 if (next==nullptr) next = vertexHead;
528
529 answer += curr->a*next->b - curr->b*next->
530 curr = curr->next;
531 } while( curr != nullptr );
532
533 return 0.5*answer;
534 }
535
536 // Print
537 //
538 void G4ReduciblePolygon::Print()
539 {
540 ABVertex* curr = vertexHead;
541 do // Loop checking, 13.08.2015, G.Cosmo
542 {
543 G4cerr << curr->a << " " << curr->b << G4e
544 curr = curr->next;
545 } while( curr != nullptr );
546 }
547
548 // CalculateMaxMin
549 //
550 // To be called when the vertices are changed,
551 // routine re-calculates global values
552 //
553 void G4ReduciblePolygon::CalculateMaxMin()
554 {
555 ABVertex* curr = vertexHead;
556 aMin = aMax = curr->a;
557 bMin = bMax = curr->b;
558 curr = curr->next;
559 while( curr != nullptr ) // Loop checking
560 {
561 if (curr->a < aMin)
562 aMin = curr->a;
563 else if (curr->a > aMax)
564 aMax = curr->a;
565
566 if (curr->b < bMin)
567 bMin = curr->b;
568 else if (curr->b > bMax)
569 bMax = curr->b;
570
571 curr = curr->next;
572 }
573 }
574