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