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