Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ReduciblePolygon.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: G4ReduciblePolygon.cc 92024 2015-08-13 14:16:00Z gcosmo $
28 //
29 //
30 // --------------------------------------------------------------------
31 // GEANT 4 class source file
32 //
33 //
34 // G4ReduciblePolygon.cc
35 //
36 // Implementation of a utility class used to specify, test, reduce,
37 // and/or otherwise manipulate a 2D polygon.
38 //
39 // See G4ReduciblePolygon.hh for more info.
40 //
41 // --------------------------------------------------------------------
42 
43 #include "G4ReduciblePolygon.hh"
44 #include "globals.hh"
45 
46 //
47 // Constructor: with simple arrays
48 //
50  const G4double b[],
51  G4int n )
52  : aMin(0.), aMax(0.), bMin(0.), bMax(0.),
53  vertexHead(0)
54 {
55  //
56  // Do all of the real work in Create
57  //
58  Create( a, b, n );
59 }
60 
61 
62 //
63 // Constructor: special PGON/PCON case
64 //
66  const G4double rmax[],
67  const G4double z[], G4int n )
68  : aMin(0.), aMax(0.), bMin(0.), bMax(0.),
69  vertexHead(0)
70 {
71  //
72  // Translate
73  //
74  G4double *a = new G4double[n*2];
75  G4double *b = new G4double[n*2];
76 
77  G4double *rOut = a + n,
78  *zOut = b + n,
79  *rIn = rOut-1,
80  *zIn = zOut-1;
81 
82  G4int i;
83  for( i=0; i < n; i++, rOut++, zOut++, rIn--, zIn-- )
84  {
85  *rOut = rmax[i];
86  *rIn = rmin[i];
87  *zOut = *zIn = z[i];
88  }
89 
90  Create( a, b, n*2 );
91 
92  delete [] a;
93  delete [] b;
94 }
95 
96 
97 //
98 // Create
99 //
100 // To be called by constructors, fill in the list and statistics for a new
101 // polygon
102 //
104  const G4double b[], G4int n )
105 {
106  if (n<3)
107  G4Exception("G4ReduciblePolygon::Create()", "GeomSolids0002",
108  FatalErrorInArgument, "Less than 3 vertices specified.");
109 
110  const G4double *anext = a, *bnext = b;
111  ABVertex *prev = 0;
112  do // Loop checking, 13.08.2015, G.Cosmo
113  {
114  ABVertex *newVertex = new ABVertex;
115  newVertex->a = *anext;
116  newVertex->b = *bnext;
117  newVertex->next = 0;
118  if (prev==0)
119  {
120  vertexHead = newVertex;
121  }
122  else
123  {
124  prev->next = newVertex;
125  }
126 
127  prev = newVertex;
128  } while( ++anext, ++bnext < b+n );
129 
130  numVertices = n;
131 
132  CalculateMaxMin();
133 }
134 
135 
136 //
137 // Fake default constructor - sets only member data and allocates memory
138 // for usage restricted to object persistency.
139 //
141  : aMin(0.), aMax(0.), bMin(0.), bMax(0.), numVertices(0), vertexHead(0)
142 {
143 }
144 
145 
146 //
147 // Destructor
148 //
150 {
151  ABVertex *curr = vertexHead;
152  while( curr ) // Loop checking, 13.08.2015, G.Cosmo
153  {
154  ABVertex *toDelete = curr;
155  curr = curr->next;
156  delete toDelete;
157  }
158 }
159 
160 
161 //
162 // CopyVertices
163 //
164 // Copy contents into simple linear arrays.
165 // ***** CAUTION ***** Be care to declare the arrays to a large
166 // enough size!
167 //
169 {
170  G4double *anext = a, *bnext = b;
171  ABVertex *curr = vertexHead;
172  while( curr ) // Loop checking, 13.08.2015, G.Cosmo
173  {
174  *anext++ = curr->a;
175  *bnext++ = curr->b;
176  curr = curr->next;
177  }
178 }
179 
180 
181 //
182 // ScaleA
183 //
184 // Multiply all a values by a common scale
185 //
187 {
188  ABVertex *curr = vertexHead;
189  while( curr ) // Loop checking, 13.08.2015, G.Cosmo
190  {
191  curr->a *= scale;
192  curr = curr->next;
193  }
194 }
195 
196 
197 //
198 // ScaleB
199 //
200 // Multiply all b values by a common scale
201 //
203 {
204  ABVertex *curr = vertexHead;
205  while( curr ) // Loop checking, 13.08.2015, G.Cosmo
206  {
207  curr->b *= scale;
208  curr = curr->next;
209  }
210 }
211 
212 
213 //
214 // RemoveDuplicateVertices
215 //
216 // Remove adjacent vertices that are equal. Returns "false" if there
217 // is a problem (too few vertices remaining).
218 //
220 {
221  ABVertex *curr = vertexHead,
222  *prev = 0, *next = 0;
223  while( curr ) // Loop checking, 13.08.2015, G.Cosmo
224  {
225  next = curr->next;
226  if (next == 0) next = vertexHead;
227 
228  if (std::fabs(curr->a-next->a) < tolerance &&
229  std::fabs(curr->b-next->b) < tolerance )
230  {
231  //
232  // Duplicate found: do we have > 3 vertices?
233  //
234  if (numVertices <= 3)
235  {
236  CalculateMaxMin();
237  return false;
238  }
239 
240  //
241  // Delete
242  //
243  ABVertex *toDelete = curr;
244  curr = curr->next;
245  delete toDelete;
246 
247  numVertices--;
248 
249  if (prev) prev->next = curr; else vertexHead = curr;
250  }
251  else
252  {
253  prev = curr;
254  curr = curr->next;
255  }
256  }
257 
258  //
259  // In principle, this is not needed, but why not just play it safe?
260  //
261  CalculateMaxMin();
262 
263  return true;
264 }
265 
266 
267 //
268 // RemoveRedundantVertices
269 //
270 // Remove any unneeded vertices, i.e. those vertices which
271 // are on the line connecting the previous and next vertices.
272 //
274 {
275  //
276  // Under these circumstances, we can quit now!
277  //
278  if (numVertices <= 2) return false;
279 
280  G4double tolerance2 = tolerance*tolerance;
281 
282  //
283  // Loop over all vertices
284  //
285  ABVertex *curr = vertexHead, *next = 0;
286  while( curr ) // Loop checking, 13.08.2015, G.Cosmo
287  {
288  next = curr->next;
289  if (next == 0) next = vertexHead;
290 
291  G4double da = next->a - curr->a,
292  db = next->b - curr->b;
293 
294  //
295  // Loop over all subsequent vertices, up to curr
296  //
297  for(;;)
298  {
299  //
300  // Get vertex after next
301  //
302  ABVertex *test = next->next;
303  if (test == 0) test = vertexHead;
304 
305  //
306  // If we are back to the original vertex, stop
307  //
308  if (test==curr) break;
309 
310  //
311  // Test for parallel line segments
312  //
313  G4double dat = test->a - curr->a,
314  dbt = test->b - curr->b;
315 
316  if (std::fabs(dat*db-dbt*da)>tolerance2) break;
317 
318  //
319  // Redundant vertex found: do we have > 3 vertices?
320  //
321  if (numVertices <= 3)
322  {
323  CalculateMaxMin();
324  return false;
325  }
326 
327  //
328  // Delete vertex pointed to by next. Carefully!
329  //
330  if (curr->next)
331  { // next is not head
332  if (next->next)
333  curr->next = test; // next is not tail
334  else
335  curr->next = 0; // New tail
336  }
337  else
338  vertexHead = test; // New head
339 
340  if ((curr != next) && (next != test)) delete next;
341 
342  numVertices--;
343 
344  //
345  // Replace next by the vertex we just tested,
346  // and keep on going...
347  //
348  next = test;
349  da = dat; db = dbt;
350  }
351  curr = curr->next;
352  }
353 
354  //
355  // In principle, this is not needed, but why not just play it safe?
356  //
357  CalculateMaxMin();
358 
359  return true;
360 }
361 
362 
363 //
364 // ReverseOrder
365 //
366 // Reverse the order of the vertices
367 //
369 {
370  //
371  // Loop over all vertices
372  //
373  ABVertex *prev = vertexHead;
374  if (prev==0) return; // No vertices
375 
376  ABVertex *curr = prev->next;
377  if (curr==0) return; // Just one vertex
378 
379  //
380  // Our new tail
381  //
382  vertexHead->next = 0;
383 
384  for(;;)
385  {
386  //
387  // Save pointer to next vertex (in original order)
388  //
389  ABVertex *save = curr->next;
390 
391  //
392  // Replace it with a pointer to the previous one
393  // (in original order)
394  //
395  curr->next = prev;
396 
397  //
398  // Last vertex?
399  //
400  if (save == 0) break;
401 
402  //
403  // Next vertex
404  //
405  prev = curr;
406  curr = save;
407  }
408 
409  //
410  // Our new head
411  //
412  vertexHead = curr;
413 }
414 
415 
416 // StartWithZMin
417 //
418 // Starting alway with Zmin=bMin
419 // This method is used for GenericPolycone
420 //
422 {
423  ABVertex *curr = vertexHead;
424  G4double bcurr = curr->b;
425  ABVertex *prev = curr;
426  while( curr ) // Loop checking, 13.08.2015, G.Cosmo
427  {
428  if(curr->b < bcurr)
429  {
430  bcurr = curr->b;
431  ABVertex *curr1 = curr;
432  while( curr1 ) // Loop checking, 13.08.2015, G.Cosmo
433  {
434  if(curr1->next == 0) { curr1->next = vertexHead; break; }
435  curr1 = curr1->next;
436  }
437  vertexHead = curr;
438  prev->next = 0;
439  }
440  prev = curr;
441  curr = curr->next;
442  }
443 }
444 
445 
446 //
447 // CrossesItself
448 //
449 // Return "true" if the polygon crosses itself
450 //
451 // Warning: this routine is not very fast (runs as N**2)
452 //
454 {
455  G4double tolerance2 = tolerance*tolerance;
456  G4double one = 1.0-tolerance,
457  zero = tolerance;
458  //
459  // Top loop over line segments. By the time we finish
460  // with the second to last segment, we're done.
461  //
462  ABVertex *curr1 = vertexHead, *next1=0;
463  while (curr1->next) // Loop checking, 13.08.2015, G.Cosmo
464  {
465  next1 = curr1->next;
466  G4double da1 = next1->a-curr1->a,
467  db1 = next1->b-curr1->b;
468 
469  //
470  // Inner loop over subsequent line segments
471  //
472  ABVertex *curr2 = next1->next;
473  while( curr2 ) // Loop checking, 13.08.2015, G.Cosmo
474  {
475  ABVertex *next2 = curr2->next;
476  if (next2==0) next2 = vertexHead;
477  G4double da2 = next2->a-curr2->a,
478  db2 = next2->b-curr2->b;
479  G4double a12 = curr2->a-curr1->a,
480  b12 = curr2->b-curr1->b;
481 
482  //
483  // Calculate intersection of the two lines
484  //
485  G4double deter = da1*db2 - db1*da2;
486  if (std::fabs(deter) > tolerance2)
487  {
488  G4double s1, s2;
489  s1 = (a12*db2-b12*da2)/deter;
490 
491  if (s1 >= zero && s1 < one)
492  {
493  s2 = -(da1*b12-db1*a12)/deter;
494  if (s2 >= zero && s2 < one) return true;
495  }
496  }
497  curr2 = curr2->next;
498  }
499  curr1 = next1;
500  }
501  return false;
502 }
503 
504 
505 
506 //
507 // BisectedBy
508 //
509 // Decide if a line through two points crosses the polygon, within tolerance
510 //
512  G4double a2, G4double b2,
513  G4double tolerance )
514 {
515  G4int nNeg = 0, nPos = 0;
516 
517  G4double a12 = a2-a1, b12 = b2-b1;
518  G4double len12 = std::sqrt( a12*a12 + b12*b12 );
519  a12 /= len12; b12 /= len12;
520 
521  ABVertex *curr = vertexHead;
522  do // Loop checking, 13.08.2015, G.Cosmo
523  {
524  G4double av = curr->a - a1,
525  bv = curr->b - b1;
526 
527  G4double cross = av*b12 - bv*a12;
528 
529  if (cross < -tolerance)
530  {
531  if (nPos) return true;
532  nNeg++;
533  }
534  else if (cross > tolerance)
535  {
536  if (nNeg) return true;
537  nPos++;
538  }
539  curr = curr->next;
540  } while( curr );
541 
542  return false;
543 }
544 
545 
546 
547 //
548 // Area
549 //
550 // Calculated signed polygon area, where polygons specified in a
551 // clockwise manner (where x==a, y==b) have negative area
552 //
553 // References: [O' Rourke (C)] pp. 18-27; [Gems II] pp. 5-6:
554 // "The Area of a Simple Polygon", Jon Rokne.
555 //
557 {
558  G4double answer = 0;
559 
560  ABVertex *curr = vertexHead, *next;
561  do // Loop checking, 13.08.2015, G.Cosmo
562  {
563  next = curr->next;
564  if (next==0) next = vertexHead;
565 
566  answer += curr->a*next->b - curr->b*next->a;
567  curr = curr->next;
568  } while( curr );
569 
570  return 0.5*answer;
571 }
572 
573 
574 //
575 // Print
576 //
578 {
579  ABVertex *curr = vertexHead;
580  do // Loop checking, 13.08.2015, G.Cosmo
581  {
582  G4cerr << curr->a << " " << curr->b << G4endl;
583  curr = curr->next;
584  } while( curr );
585 }
586 
587 
588 //
589 // CalculateMaxMin
590 //
591 // To be called when the vertices are changed, this
592 // routine re-calculates global values
593 //
595 {
596  ABVertex *curr = vertexHead;
597  aMin = aMax = curr->a;
598  bMin = bMax = curr->b;
599  curr = curr->next;
600  while( curr ) // Loop checking, 13.08.2015, G.Cosmo
601  {
602  if (curr->a < aMin)
603  aMin = curr->a;
604  else if (curr->a > aMax)
605  aMax = curr->a;
606 
607  if (curr->b < bMin)
608  bMin = curr->b;
609  else if (curr->b > bMax)
610  bMax = curr->b;
611 
612  curr = curr->next;
613  }
614 }
G4bool CrossesItself(G4double tolerance)
void ScaleB(G4double scale)
void Create(const G4double a[], const G4double b[], G4int n)
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
int G4int
Definition: G4Types.hh:78
G4bool RemoveDuplicateVertices(G4double tolerance)
tuple b
Definition: test.py:12
G4bool RemoveRedundantVertices(G4double tolerance)
void ScaleA(G4double scale)
bool G4bool
Definition: G4Types.hh:79
const G4int n
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
void CopyVertices(G4double a[], G4double b[]) const
tuple z
Definition: test.py:28
G4bool BisectedBy(G4double a1, G4double b1, G4double a2, G4double b2, G4double tolerance)
G4ReduciblePolygon(const G4double a[], const G4double b[], G4int n)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
program test
Definition: Main_HIJING.f:1
G4GLOB_DLL std::ostream G4cerr