Geant4  9.6.p02
 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 67011 2013-01-29 16:17:41Z 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
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 {
152  while( curr )
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;
172  while( curr )
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 {
189  while( curr )
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 {
205  while( curr )
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 {
222  *prev = 0, *next = 0;
223  while( curr )
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 )
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 //
417 // CrossesItself
418 //
419 // Return "true" if the polygon crosses itself
420 //
421 // Warning: this routine is not very fast (runs as N**2)
422 //
424 {
425  G4double tolerance2 = tolerance*tolerance;
426  G4double one = 1.0-tolerance,
427  zero = tolerance;
428  //
429  // Top loop over line segments. By the time we finish
430  // with the second to last segment, we're done.
431  //
432  ABVertex *curr1 = vertexHead, *next1=0;
433  while (curr1->next)
434  {
435  next1 = curr1->next;
436  G4double da1 = next1->a-curr1->a,
437  db1 = next1->b-curr1->b;
438 
439  //
440  // Inner loop over subsequent line segments
441  //
442  ABVertex *curr2 = next1->next;
443  while( curr2 )
444  {
445  ABVertex *next2 = curr2->next;
446  if (next2==0) next2 = vertexHead;
447  G4double da2 = next2->a-curr2->a,
448  db2 = next2->b-curr2->b;
449  G4double a12 = curr2->a-curr1->a,
450  b12 = curr2->b-curr1->b;
451 
452  //
453  // Calculate intersection of the two lines
454  //
455  G4double deter = da1*db2 - db1*da2;
456  if (std::fabs(deter) > tolerance2)
457  {
458  G4double s1, s2;
459  s1 = (a12*db2-b12*da2)/deter;
460 
461  if (s1 >= zero && s1 < one)
462  {
463  s2 = -(da1*b12-db1*a12)/deter;
464  if (s2 >= zero && s2 < one) return true;
465  }
466  }
467  curr2 = curr2->next;
468  }
469  curr1 = next1;
470  }
471  return false;
472 }
473 
474 
475 
476 //
477 // BisectedBy
478 //
479 // Decide if a line through two points crosses the polygon, within tolerance
480 //
482  G4double a2, G4double b2,
483  G4double tolerance )
484 {
485  G4int nNeg = 0, nPos = 0;
486 
487  G4double a12 = a2-a1, b12 = b2-b1;
488  G4double len12 = std::sqrt( a12*a12 + b12*b12 );
489  a12 /= len12; b12 /= len12;
490 
492  do
493  {
494  G4double av = curr->a - a1,
495  bv = curr->b - b1;
496 
497  G4double cross = av*b12 - bv*a12;
498 
499  if (cross < -tolerance)
500  {
501  if (nPos) return true;
502  nNeg++;
503  }
504  else if (cross > tolerance)
505  {
506  if (nNeg) return true;
507  nPos++;
508  }
509  curr = curr->next;
510  } while( curr );
511 
512  return false;
513 }
514 
515 
516 
517 //
518 // Area
519 //
520 // Calculated signed polygon area, where polygons specified in a
521 // clockwise manner (where x==a, y==b) have negative area
522 //
523 // References: [O' Rourke (C)] pp. 18-27; [Gems II] pp. 5-6:
524 // "The Area of a Simple Polygon", Jon Rokne.
525 //
527 {
528  G4double answer = 0;
529 
530  ABVertex *curr = vertexHead, *next;
531  do
532  {
533  next = curr->next;
534  if (next==0) next = vertexHead;
535 
536  answer += curr->a*next->b - curr->b*next->a;
537  curr = curr->next;
538  } while( curr );
539 
540  return 0.5*answer;
541 }
542 
543 
544 //
545 // Print
546 //
548 {
550  do
551  {
552  G4cerr << curr->a << " " << curr->b << G4endl;
553  curr = curr->next;
554  } while( curr );
555 }
556 
557 
558 //
559 // CalculateMaxMin
560 //
561 // To be called when the vertices are changed, this
562 // routine re-calculates global values
563 //
565 {
567  aMin = aMax = curr->a;
568  bMin = bMax = curr->b;
569  curr = curr->next;
570  while( curr )
571  {
572  if (curr->a < aMin)
573  aMin = curr->a;
574  else if (curr->a > aMax)
575  aMax = curr->a;
576 
577  if (curr->b < bMin)
578  bMin = curr->b;
579  else if (curr->b > bMax)
580  bMax = curr->b;
581 
582  curr = curr->next;
583  }
584 }