ROOT  6.06/09
Reference Guide
CsgOps.cxx
Go to the documentation of this file.
1 // @(#)root/gl:$Id$
2 // Author: Timur Pocheptsov 01/04/2005
3 /*
4  CSGLib - Software Library for Constructive Solid Geometry
5  Copyright (C) 2003-2004 Laurence Bourn
6 
7  This library is free software; you can redistribute it and/or
8  modify it under the terms of the GNU Library General Public
9  License as published by the Free Software Foundation; either
10  version 2 of the License, or (at your option) any later version.
11 
12  This library is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15  Library General Public License for more details.
16 
17  You should have received a copy of the GNU Library General Public
18  License along with this library; if not, write to the Free
19  Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
20 
21  Please send remarks, questions and bug reports to laurencebourn@hotmail.com
22 */
23 
24 /**
25  * I've modified some very nice bounding box tree code from
26  * Gino van der Bergen's Free Solid Library below. It's basically
27  * the same code - but I've hacked out the transformation stuff as
28  * I didn't understand it. I've also made it far less elegant!
29  * Laurence Bourn.
30  */
31 
32 /*
33  SOLID - Software Library for Interference Detection
34  Copyright (C) 1997-1998 Gino van den Bergen
35 
36  This library is free software; you can redistribute it and/or
37  modify it under the terms of the GNU Library General Public
38  License as published by the Free Software Foundation; either
39  version 2 of the License, or (at your option) any later version.
40 
41  This library is distributed in the hope that it will be useful,
42  but WITHOUT ANY WARRANTY; without even the implied warranty of
43  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
44  Library General Public License for more details.
45 
46  You should have received a copy of the GNU Library General Public
47  License along with this library; if not, write to the Free
48  Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
49 
50  Please send remarks, questions and bug reports to gino@win.tue.nl,
51  or write to:
52  Gino van den Bergen
53  Department of Mathematics and Computing Science
54  Eindhoven University of Technology
55  P.O. Box 513, 5600 MB Eindhoven, The Netherlands
56 */
57 
58 /*
59  This file contains compressed version of CSGSolid and SOLID.
60  All stuff is in RootCsg namespace,
61  I used ROOT's own typedefs and math functions, removed resulting triangulation
62  (I use OpenGL's inner implementation), changed names from MT_xxx CSG_xxx
63  to more natural etc.
64  Interface to CSGSolid :
65  ConvertToMesh(to convert TBuffer3D to mesh)
66  BuildUnion
67  BuildIntersection
68  BuildDifference
69 
70  31.03.05 Timur Pocheptsov.
71 */
72 
73 #include <algorithm>
74 #include <vector>
75 
76 #include "TBuffer3D.h"
77 #include "Rtypes.h"
78 #include "TMath.h"
79 #include "CsgOps.h"
80 
81 namespace RootCsg {
82 
83  const Double_t epsilon = 1e-10;
84  const Double_t epsilon2 = 1e-20;
85  const Double_t infinity = 1e50;
86 
87  /////////////////////////////////////////////////////////////////////////////
88 
90  {
91  return x < 0. ? -1 : x > 0. ? 1 : 0;
92  }
93 
94  /////////////////////////////////////////////////////////////////////////////
95 
97  {
98  return TMath::Abs(x) < epsilon;
99  }
100 
101  /////////////////////////////////////////////////////////////////////////////
102 
104  {
105  return TMath::Abs(x) < epsilon2;
106  }
107 
108  class Tuple2 {
109  protected:
110  Double_t fCo[2];
111 
112  public:
113  Tuple2(){SetValue(0, 0);}
114  Tuple2(const Double_t *vv){SetValue(vv);}
115  Tuple2(Double_t xx, Double_t yy){SetValue(xx, yy);}
116 
117  Double_t &operator [] (Int_t i){return fCo[i];}
118  const Double_t &operator [] (Int_t i)const{return fCo[i];}
119 
120  Double_t &X(){return fCo[0];}
121  const Double_t &X()const{return fCo[0];}
122  Double_t &Y(){return fCo[1];}
123  const Double_t &Y()const{return fCo[1];}
124  Double_t &U(){return fCo[0];}
125  const Double_t &U()const{return fCo[0];}
126  Double_t &V(){return fCo[1];}
127  const Double_t &V()const{return fCo[1];}
128 
129  Double_t *GetValue(){return fCo;}
130  const Double_t *GetValue()const{return fCo;}
131  void GetValue(Double_t *vv)const{vv[0] = fCo[0]; vv[1] = fCo[1];}
132 
133  void SetValue(const Double_t *vv){fCo[0] = vv[0]; fCo[1] = vv[1];}
134  void SetValue(Double_t xx, Double_t yy){fCo[0] = xx; fCo[1] = yy;}
135  };
136 
137  Bool_t operator == (const Tuple2 &t1, const Tuple2 &t2)
138  {
139  return t1[0] == t2[0] && t1[1] == t2[1];
140  }
141 
142  class TVector2 : public Tuple2 {
143  public:
144  TVector2(){}
145  TVector2(const Double_t *v) : Tuple2(v) {}
146  TVector2(Double_t xx, Double_t yy) : Tuple2(xx, yy) {}
147 
148  TVector2 &operator += (const TVector2 &v);
149  TVector2 &operator -= (const TVector2 &v);
150  TVector2 &operator *= (Double_t s);
151  TVector2 &operator /= (Double_t s);
152 
153  Double_t Dot(const TVector2 &v)const;
154  Double_t Length2()const;
155  Double_t Length()const;
156  TVector2 Absolute()const;
157  void Normalize();
158  TVector2 Normalized()const;
159  void Scale(Double_t x, Double_t y);
160  TVector2 Scaled(Double_t x, Double_t y)const;
161  Bool_t FuzzyZero()const;
162  Double_t Angle(const TVector2 &v)const;
163  TVector2 Cross(const TVector2 &v)const;
164  Double_t Triple(const TVector2 &v1, const TVector2 &v2)const;
165  };
166 
167  /////////////////////////////////////////////////////////////////////////////
168  ///
169 
171  {
172  fCo[0] += vv[0]; fCo[1] += vv[1];
173  return *this;
174  }
175 
176  /////////////////////////////////////////////////////////////////////////////
177  ///
178 
180  {
181  fCo[0] -= vv[0]; fCo[1] -= vv[1];
182  return *this;
183  }
184 
185  /////////////////////////////////////////////////////////////////////////////
186  ///
187 
189  {
190  fCo[0] *= s; fCo[1] *= s; return *this;
191  }
192 
193  /////////////////////////////////////////////////////////////////////////////
194  ///
195 
197  {
198  return *this *= 1. / s;
199  }
200 
201  /////////////////////////////////////////////////////////////////////////////
202  ///
203 
204  TVector2 operator + (const TVector2 &v1, const TVector2 &v2)
205  {
206  return TVector2(v1[0] + v2[0], v1[1] + v2[1]);
207  }
208 
209  /////////////////////////////////////////////////////////////////////////////
210  ///
211 
212  TVector2 operator - (const TVector2 &v1, const TVector2 &v2)
213  {
214  return TVector2(v1[0] - v2[0], v1[1] - v2[1]);
215  }
216 
217  /////////////////////////////////////////////////////////////////////////////
218  ///
219 
220  TVector2 operator - (const TVector2 &v)
221  {
222  return TVector2(-v[0], -v[1]);
223  }
224 
225  /////////////////////////////////////////////////////////////////////////////
226  ///
227 
228  TVector2 operator * (const TVector2 &v, Double_t s)
229  {
230  return TVector2(v[0] * s, v[1] * s);
231  }
232 
233  /////////////////////////////////////////////////////////////////////////////
234  ///
235 
236  TVector2 operator * (Double_t s, const TVector2 &v)
237  {
238  return v * s;
239  }
240 
241  /////////////////////////////////////////////////////////////////////////////
242  ///
243 
244  TVector2 operator / (const TVector2 & v, Double_t s)
245  {
246  return v * (1.0 / s);
247  }
248 
249  /////////////////////////////////////////////////////////////////////////////
250  ///
251 
252  Double_t TVector2::Dot(const TVector2 &vv)const
253  {
254  return fCo[0] * vv[0] + fCo[1] * vv[1];
255  }
256 
257  /////////////////////////////////////////////////////////////////////////////
258  ///
259 
260  Double_t TVector2::Length2()const
261  {
262  return Dot(*this);
263  }
264 
265  /////////////////////////////////////////////////////////////////////////////
266  ///
267 
268  Double_t TVector2::Length()const
269  {
270  return TMath::Sqrt(Length2());
271  }
272 
273  /////////////////////////////////////////////////////////////////////////////
274  ///
275 
277  {
278  return TVector2(TMath::Abs(fCo[0]), TMath::Abs(fCo[1]));
279  }
280 
281  /////////////////////////////////////////////////////////////////////////////
282  ///
283 
284  Bool_t TVector2::FuzzyZero()const
285  {
286  return fuzzy_zero2(Length2());
287  }
288 
289  /////////////////////////////////////////////////////////////////////////////
290  ///
291 
292  void TVector2::Normalize()
293  {
294  *this /= Length();
295  }
296 
297  /////////////////////////////////////////////////////////////////////////////
298  ///
299 
300  TVector2 TVector2::Normalized()const
301  {
302  return *this / Length();
303  }
304 
305  /////////////////////////////////////////////////////////////////////////////
306  ///
307 
308  void TVector2::Scale(Double_t xx, Double_t yy)
309  {
310  fCo[0] *= xx; fCo[1] *= yy;
311  }
312 
313  /////////////////////////////////////////////////////////////////////////////
314  ///
315 
316  TVector2 TVector2::Scaled(Double_t xx, Double_t yy)const
317  {
318  return TVector2(fCo[0] * xx, fCo[1] * yy);
319  }
320 
321  /////////////////////////////////////////////////////////////////////////////
322  ///
323 
324  Double_t TVector2::Angle(const TVector2 &vv)const
325  {
326  Double_t s = TMath::Sqrt(Length2() * vv.Length2());
327  return TMath::ACos(Dot(vv) / s);
328  }
329 
330  /////////////////////////////////////////////////////////////////////////////
331  ///
332 
333  Double_t dot(const TVector2 &v1, const TVector2 &v2)
334  {
335  return v1.Dot(v2);
336  }
337 
338  /////////////////////////////////////////////////////////////////////////////
339 
340  Double_t length2(const TVector2 &v)
341  {
342  return v.Length2();
343  }
344 
345  /////////////////////////////////////////////////////////////////////////////
346 
347  Double_t length(const TVector2 &v)
348  {
349  return v.Length();
350  }
351 
352  /////////////////////////////////////////////////////////////////////////////
353 
354  Bool_t fuzzy_zero(const TVector2 &v)
355  {
356  return v.FuzzyZero();
357  }
358 
359  /////////////////////////////////////////////////////////////////////////////
360 
361  Bool_t fuzzy_equal(const TVector2 &v1, const TVector2 &v2)
362  {
363  return fuzzy_zero(v1 - v2);
364  }
365 
366  /////////////////////////////////////////////////////////////////////////////
367 
368  Double_t Angle(const TVector2 &v1, const TVector2 &v2)
369  {
370  return v1.Angle(v2);
371  }
372 
373  class TPoint2 : public TVector2 {
374  public:
375  TPoint2(){}
376  TPoint2(const Double_t *v) : TVector2(v){}
377  TPoint2(Double_t x, Double_t y) : TVector2(x, y) {}
378 
379  TPoint2 &operator += (const TVector2 &v);
380  TPoint2 &operator -= (const TVector2 &v);
381  TPoint2 &operator = (const TVector2 &v);
382 
383  Double_t Distance(const TPoint2 &p)const;
384  Double_t Distance2(const TPoint2 &p)const;
385  TPoint2 Lerp(const TPoint2 &p, Double_t t)const;
386  };
387 
388  /////////////////////////////////////////////////////////////////////////////
389  ///
390 
391  TPoint2 &TPoint2::operator += (const TVector2 &v)
392  {
393  fCo[0] += v[0]; fCo[1] += v[1];
394  return *this;
395  }
396 
397  /////////////////////////////////////////////////////////////////////////////
398  ///
399 
400  TPoint2 &TPoint2::operator-=(const TVector2& v)
401  {
402  fCo[0] -= v[0]; fCo[1] -= v[1];
403  return *this;
404  }
405 
406  /////////////////////////////////////////////////////////////////////////////
407  ///
408 
409  TPoint2 &TPoint2::operator=(const TVector2& v)
410  {
411  fCo[0] = v[0]; fCo[1] = v[1];
412  return *this;
413  }
414 
415  /////////////////////////////////////////////////////////////////////////////
416  ///
417 
418  Double_t TPoint2::Distance(const TPoint2& p)const
419  {
420  return (p - *this).Length();
421  }
422 
423  /////////////////////////////////////////////////////////////////////////////
424  ///
425 
426  Double_t TPoint2::Distance2(const TPoint2& p)const
427  {
428  return (p - *this).Length2();
429  }
430 
431  /////////////////////////////////////////////////////////////////////////////
432  ///
433 
434  TPoint2 TPoint2::Lerp(const TPoint2 &p, Double_t t)const
435  {
436  return TPoint2(fCo[0] + (p[0] - fCo[0]) * t,
437  fCo[1] + (p[1] - fCo[1]) * t);
438  }
439 
440  /////////////////////////////////////////////////////////////////////////////
441  ///
442 
443  TPoint2 operator + (const TPoint2 &p, const TVector2 &v)
444  {
445  return TPoint2(p[0] + v[0], p[1] + v[1]);
446  }
447 
448  /////////////////////////////////////////////////////////////////////////////
449  ///
450 
451  TPoint2 operator - (const TPoint2 &p, const TVector2 &v)
452  {
453  return TPoint2(p[0] - v[0], p[1] - v[1]);
454  }
455 
456  /////////////////////////////////////////////////////////////////////////////
457  ///
458 
459  TVector2 operator - (const TPoint2 &p1, const TPoint2 &p2)
460  {
461  return TVector2(p1[0] - p2[0], p1[1] - p2[1]);
462  }
463 
464  /////////////////////////////////////////////////////////////////////////////
465  ///
466 
467  Double_t distance(const TPoint2 &p1, const TPoint2 &p2)
468  {
469  return p1.Distance(p2);
470  }
471 
472  /////////////////////////////////////////////////////////////////////////////
473  ///
474 
475  Double_t distance2(const TPoint2 &p1, const TPoint2 &p2)
476  {
477  return p1.Distance2(p2);
478  }
479 
480  /////////////////////////////////////////////////////////////////////////////
481  ///
482 
483  TPoint2 lerp(const TPoint2 &p1, const TPoint2 &p2, Double_t t)
484  {
485  return p1.Lerp(p2, t);
486  }
487 
488  class Tuple3 {
489  protected:
490  Double_t fCo[3];
491  public:
492  Tuple3(){SetValue(0, 0, 0);}
493  Tuple3(const Double_t *v){SetValue(v);}
494  Tuple3(Double_t xx, Double_t yy, Double_t zz){SetValue(xx, yy, zz);}
495 
496  Double_t &operator [] (Int_t i){return fCo[i];}
497  const Double_t &operator [] (Int_t i)const{return fCo[i];}
498 
499  Double_t &X(){return fCo[0];}
500  const Double_t &X()const{return fCo[0];}
501  Double_t &Y(){return fCo[1];}
502  const Double_t &Y()const{return fCo[1];}
503  Double_t &Z(){return fCo[2];}
504  const Double_t &Z()const{return fCo[2];}
505 
506  Double_t *GetValue(){return fCo;}
507  const Double_t *GetValue()const{ return fCo; }
508 
509  void GetValue(Double_t *v)const
510  {
511  v[0] = Double_t(fCo[0]), v[1] = Double_t(fCo[1]), v[2] = Double_t(fCo[2]);
512  }
513 
514  void SetValue(const Double_t *v)
515  {
516  fCo[0] = Double_t(v[0]), fCo[1] = Double_t(v[1]), fCo[2] = Double_t(v[2]);
517  }
518  void SetValue(Double_t xx, Double_t yy, Double_t zz)
519  {
520  fCo[0] = xx; fCo[1] = yy; fCo[2] = zz;
521  }
522  };
523 
524  /////////////////////////////////////////////////////////////////////////////
525 
526  Bool_t operator==(const Tuple3& t1, const Tuple3& t2)
527  {
528  return t1[0] == t2[0] && t1[1] == t2[1] && t1[2] == t2[2];
529  }
530 
531  class TVector3 : public Tuple3 {
532  public:
533  TVector3(){}
534  TVector3(const Double_t *v) : Tuple3(v){}
535  TVector3(Double_t xx, Double_t yy, Double_t zz) : Tuple3(xx, yy, zz){}
536 
537  TVector3 &operator += (const TVector3& v);
538  TVector3 &operator -= (const TVector3& v);
539  TVector3 &operator *= (Double_t s);
540  TVector3 &operator /= (Double_t s);
541 
542  Double_t Dot(const TVector3& v)const;
543  Double_t Length2()const;
544  Double_t Length()const;
545  TVector3 Absolute()const;
546  void NoiseGate(Double_t threshold);
547  void Normalize();
548  TVector3 Normalized()const;
549  TVector3 SafeNormalized()const;
550  void Scale(Double_t x, Double_t y, Double_t z);
551  TVector3 Scaled(Double_t x, Double_t y, Double_t z)const;
552  Bool_t FuzzyZero()const;
553  Double_t Angle(const TVector3 &v)const;
554  TVector3 Cross(const TVector3 &v)const;
555  Double_t Triple(const TVector3 &v1, const TVector3 &v2)const;
556  Int_t ClosestAxis()const;
557 
558  static TVector3 Random();
559  };
560 
561  /////////////////////////////////////////////////////////////////////////////
562  ///
563 
565  {
566  fCo[0] += v[0]; fCo[1] += v[1]; fCo[2] += v[2];
567  return *this;
568  }
569 
570  /////////////////////////////////////////////////////////////////////////////
571  ///
572 
574  {
575  fCo[0] -= v[0]; fCo[1] -= v[1]; fCo[2] -= v[2];
576  return *this;
577  }
578 
579  /////////////////////////////////////////////////////////////////////////////
580  ///
581 
583  {
584  fCo[0] *= s; fCo[1] *= s; fCo[2] *= s;
585  return *this;
586  }
587 
588  /////////////////////////////////////////////////////////////////////////////
589  ///
590 
591  TVector3 &TVector3::operator /= (Double_t s)
592  {
593  return *this *= Double_t(1.0) / s;
594  }
595 
596  /////////////////////////////////////////////////////////////////////////////
597  ///
598 
599  Double_t TVector3::Dot(const TVector3 &v)const
600  {
601  return fCo[0] * v[0] + fCo[1] * v[1] + fCo[2] * v[2];
602  }
603 
604  /////////////////////////////////////////////////////////////////////////////
605  ///
606 
607  Double_t TVector3::Length2()const
608  {
609  return Dot(*this);
610  }
611 
612  /////////////////////////////////////////////////////////////////////////////
613  ///
614 
615  Double_t TVector3::Length()const
616  {
617  return TMath::Sqrt(Length2());
618  }
619 
620  /////////////////////////////////////////////////////////////////////////////
621  ///
622 
624  {
625  return TVector3(TMath::Abs(fCo[0]), TMath::Abs(fCo[1]), TMath::Abs(fCo[2]));
626  }
627 
628  /////////////////////////////////////////////////////////////////////////////
629  ///
630 
631  Bool_t TVector3::FuzzyZero()const
632  {
633  return fuzzy_zero(Length2());
634  }
635 
636  /////////////////////////////////////////////////////////////////////////////
637  ///
638 
639  void TVector3::NoiseGate(Double_t threshold)
640  {
641  if (Length2() < threshold) SetValue(0., 0., 0.);
642  }
643 
644  /////////////////////////////////////////////////////////////////////////////
645  ///
646 
647  void TVector3::Normalize()
648  {
649  *this /= Length();
650  }
651 
652  /////////////////////////////////////////////////////////////////////////////
653 
654  TVector3 operator * (const TVector3 &v, Double_t s)
655  {
656  return TVector3(v[0] * s, v[1] * s, v[2] * s);
657  }
658 
659  /////////////////////////////////////////////////////////////////////////////
660 
661  TVector3 operator / (const TVector3 &v, Double_t s)
662  {
663  return v * (1. / s);
664  }
665 
666  /////////////////////////////////////////////////////////////////////////////
667  ///
668 
669  TVector3 TVector3::Normalized()const
670  {
671  return *this / Length();
672  }
673 
674  /////////////////////////////////////////////////////////////////////////////
675  ///
676 
677  TVector3 TVector3::SafeNormalized()const
678  {
679  Double_t len = Length();
680  return fuzzy_zero(len) ? TVector3(1., 0., 0.):*this / len;
681  }
682 
683  /////////////////////////////////////////////////////////////////////////////
684  ///
685 
686  void TVector3::Scale(Double_t xx, Double_t yy, Double_t zz)
687  {
688  fCo[0] *= xx; fCo[1] *= yy; fCo[2] *= zz;
689  }
690 
691  /////////////////////////////////////////////////////////////////////////////
692  ///
693 
694  TVector3 TVector3::Scaled(Double_t xx, Double_t yy, Double_t zz)const
695  {
696  return TVector3(fCo[0] * xx, fCo[1] * yy, fCo[2] * zz);
697  }
698 
699  /////////////////////////////////////////////////////////////////////////////
700  ///
701 
702  Double_t TVector3::Angle(const TVector3 &v)const
703  {
704  Double_t s = TMath::Sqrt(Length2() * v.Length2());
705  return TMath::ACos(Dot(v) / s);
706  }
707 
708  /////////////////////////////////////////////////////////////////////////////
709  ///
710 
711  TVector3 TVector3::Cross(const TVector3 &v)const
712  {
713  return TVector3(fCo[1] * v[2] - fCo[2] * v[1],
714  fCo[2] * v[0] - fCo[0] * v[2],
715  fCo[0] * v[1] - fCo[1] * v[0]);
716  }
717 
718  /////////////////////////////////////////////////////////////////////////////
719  ///
720 
721  Double_t TVector3::Triple(const TVector3 &v1, const TVector3 &v2)const
722  {
723  return fCo[0] * (v1[1] * v2[2] - v1[2] * v2[1]) +
724  fCo[1] * (v1[2] * v2[0] - v1[0] * v2[2]) +
725  fCo[2] * (v1[0] * v2[1] - v1[1] * v2[0]);
726  }
727 
728  /////////////////////////////////////////////////////////////////////////////
729  ///
730 
731  Int_t TVector3::ClosestAxis()const
732  {
733  TVector3 a = Absolute();
734  return a[0] < a[1] ? (a[1] < a[2] ? 2 : 1) : (a[0] < a[2] ? 2 : 0);
735  }
736 
737  /////////////////////////////////////////////////////////////////////////////
738 
739  TVector3 operator + (const TVector3 &v1, const TVector3 &v2)
740  {
741  return TVector3(v1[0] + v2[0], v1[1] + v2[1], v1[2] + v2[2]);
742  }
743 
744  /////////////////////////////////////////////////////////////////////////////
745 
746  TVector3 operator - (const TVector3 &v1, const TVector3 &v2)
747  {
748  return TVector3(v1[0] - v2[0], v1[1] - v2[1], v1[2] - v2[2]);
749  }
750 
751  /////////////////////////////////////////////////////////////////////////////
752 
753  TVector3 operator - (const TVector3 &v)
754  {
755  return TVector3(-v[0], -v[1], -v[2]);
756  }
757 
758  /////////////////////////////////////////////////////////////////////////////
759 
760  TVector3 operator * (Double_t s, const TVector3 &v)
761  {
762  return v * s;
763  }
764 
765  /////////////////////////////////////////////////////////////////////////////
766 
767  TVector3 operator * (const TVector3 &v1, const TVector3 &v2)
768  {
769  return TVector3(v1[0] * v2[0], v1[1] * v2[1], v1[2] * v2[2]);
770  }
771 
772  /////////////////////////////////////////////////////////////////////////////
773 
774  Double_t dot(const TVector3 &v1, const TVector3 &v2)
775  {
776  return v1.Dot(v2);
777  }
778 
779  /////////////////////////////////////////////////////////////////////////////
780 
781  Double_t length2(const TVector3 &v)
782  {
783  return v.Length2();
784  }
785 
786  /////////////////////////////////////////////////////////////////////////////
787 
788  Double_t length(const TVector3 &v)
789  {
790  return v.Length();
791  }
792 
793  /////////////////////////////////////////////////////////////////////////////
794 
795  Bool_t fuzzy_zero(const TVector3 &v)
796  {
797  return v.FuzzyZero();
798  }
799 
800  /////////////////////////////////////////////////////////////////////////////
801 
802  Bool_t fuzzy_equal(const TVector3 &v1, const TVector3 &v2)
803  {
804  return fuzzy_zero(v1 - v2);
805  }
806 
807  /////////////////////////////////////////////////////////////////////////////
808 
809  Double_t Angle(const TVector3 &v1, const TVector3 &v2)
810  {
811  return v1.Angle(v2);
812  }
813 
814  /////////////////////////////////////////////////////////////////////////////
815 
816  TVector3 cross(const TVector3 &v1, const TVector3 &v2)
817  {
818  return v1.Cross(v2);
819  }
820 
821  /////////////////////////////////////////////////////////////////////////////
822 
823  Double_t triple(const TVector3 &v1, const TVector3 &v2, const TVector3 &v3)
824  {
825  return v1.Triple(v2, v3);
826  }
827 
828  class TPoint3 : public TVector3 {
829  public:
830  TPoint3(){}
831  TPoint3(const Double_t *v) : TVector3(v) {}
832  TPoint3(Double_t xx, Double_t yy, Double_t zz) : TVector3(xx, yy, zz) {}
833 
834  TPoint3 &operator += (const TVector3 &v);
835  TPoint3 &operator -= (const TVector3 &v);
836  TPoint3 &operator = (const TVector3 &v);
837 
838  Double_t Distance(const TPoint3 &p)const;
839  Double_t Distance2(const TPoint3 &p)const;
840  TPoint3 Lerp(const TPoint3 &p, Double_t t)const;
841  };
842 
843  /////////////////////////////////////////////////////////////////////////////
844  ///
845 
846  TPoint3 &TPoint3::operator += (const TVector3 &v)
847  {
848  fCo[0] += v[0]; fCo[1] += v[1]; fCo[2] += v[2];
849  return *this;
850  }
851 
852  /////////////////////////////////////////////////////////////////////////////
853  ///
854 
855  TPoint3 &TPoint3::operator-=(const TVector3 &v)
856  {
857  fCo[0] -= v[0]; fCo[1] -= v[1]; fCo[2] -= v[2];
858  return *this;
859  }
860 
861  /////////////////////////////////////////////////////////////////////////////
862  ///
863 
864  TPoint3 &TPoint3::operator=(const TVector3 &v)
865  {
866  fCo[0] = v[0]; fCo[1] = v[1]; fCo[2] = v[2];
867  return *this;
868  }
869 
870  /////////////////////////////////////////////////////////////////////////////
871  ///
872 
873  Double_t TPoint3::Distance(const TPoint3 &p)const
874  {
875  return (p - *this).Length();
876  }
877 
878  /////////////////////////////////////////////////////////////////////////////
879  ///
880 
881  Double_t TPoint3::Distance2(const TPoint3 &p)const
882  {
883  return (p - *this).Length2();
884  }
885 
886  /////////////////////////////////////////////////////////////////////////////
887  ///
888 
889  TPoint3 TPoint3::Lerp(const TPoint3 &p, Double_t t)const
890  {
891  return TPoint3(fCo[0] + (p[0] - fCo[0]) * t,
892  fCo[1] + (p[1] - fCo[1]) * t,
893  fCo[2] + (p[2] - fCo[2]) * t);
894  }
895 
896  /////////////////////////////////////////////////////////////////////////////
897 
898  TPoint3 operator + (const TPoint3 &p, const TVector3 &v)
899  {
900  return TPoint3(p[0] + v[0], p[1] + v[1], p[2] + v[2]);
901  }
902 
903  /////////////////////////////////////////////////////////////////////////////
904 
905  TPoint3 operator - (const TPoint3 &p, const TVector3 &v)
906  {
907  return TPoint3(p[0] - v[0], p[1] - v[1], p[2] - v[2]);
908  }
909 
910  /////////////////////////////////////////////////////////////////////////////
911 
912  TVector3 operator - (const TPoint3 &p1, const TPoint3 &p2)
913  {
914  return TVector3(p1[0] - p2[0], p1[1] - p2[1], p1[2] - p2[2]);
915  }
916 
917  /////////////////////////////////////////////////////////////////////////////
918 
919  Double_t distance(const TPoint3 &p1, const TPoint3 &p2)
920  {
921  return p1.Distance(p2);
922  }
923 
924  /////////////////////////////////////////////////////////////////////////////
925 
926  Double_t distance2(const TPoint3 &p1, const TPoint3 &p2)
927  {
928  return p1.Distance2(p2);
929  }
930 
931  /////////////////////////////////////////////////////////////////////////////
932 
933  TPoint3 lerp(const TPoint3 &p1, const TPoint3 &p2, Double_t t)
934  {
935  return p1.Lerp(p2, t);
936  }
937 
938  class Tuple4 {
939  protected:
940  Double_t fCo[4];
941 
942  public:
943  Tuple4(){SetValue(0, 0, 0, 0);}
944  Tuple4(const Double_t *v){SetValue(v);}
945  Tuple4(Double_t xx, Double_t yy, Double_t zz, Double_t ww)
946  {
947  SetValue(xx, yy, zz, ww);
948  }
949 
950  Double_t &operator [] (Int_t i){return fCo[i];}
951  const Double_t &operator [] (Int_t i)const{return fCo[i];}
952 
953  Double_t &X(){return fCo[0];}
954  const Double_t &X()const{return fCo[0];}
955  Double_t &Y(){return fCo[1];}
956  const Double_t &Y()const{return fCo[1];}
957  Double_t &Z(){return fCo[2];}
958  const Double_t &Z()const{return fCo[2];}
959  Double_t &W(){return fCo[3];}
960  const Double_t &W()const{return fCo[3];}
961 
962  Double_t *GetValue(){return fCo;}
963  const Double_t *GetValue()const{return fCo;}
964 
965  void GetValue(Double_t *v)const
966  {
967  v[0] = fCo[0]; v[1] = fCo[1]; v[2] = fCo[2]; v[3] = fCo[3];
968  }
969 
970  void SetValue(const Double_t *v)
971  {
972  fCo[0] = v[0]; fCo[1] = v[1]; fCo[2] = v[2]; fCo[3] = v[3];
973  }
974  void SetValue(Double_t xx, Double_t yy, Double_t zz, Double_t ww)
975  {
976  fCo[0] = xx; fCo[1] = yy; fCo[2] = zz; fCo[3] = ww;
977  }
978  };
979 
980  /////////////////////////////////////////////////////////////////////////////
981 
982  Bool_t operator == (const Tuple4 &t1, const Tuple4 &t2)
983  {
984  return t1[0] == t2[0] && t1[1] == t2[1] && t1[2] == t2[2] && t1[3] == t2[3];
985  }
986 
987  class TMatrix3x3 {
988  protected:
989  TVector3 fEl[3];
990 
991  public:
992  TMatrix3x3(){}
993  TMatrix3x3(const Double_t *m){SetValue(m);}
994  TMatrix3x3(const TVector3 &euler){SetEuler(euler);}
995  TMatrix3x3(const TVector3 &euler, const TVector3 &s)
996  {
997  SetEuler(euler); Scale(s[0], s[1], s[2]);
998  }
999  TMatrix3x3(Double_t xx, Double_t xy, Double_t xz,
1000  Double_t yx, Double_t yy, Double_t yz,
1001  Double_t zx, Double_t zy, Double_t zz)
1002  {
1003  SetValue(xx, xy, xz, yx, yy, yz, zx, zy, zz);
1004  }
1005 
1006  TVector3 &operator [] (Int_t i){return fEl[i];}
1007  const TVector3 &operator [] (Int_t i)const{return fEl[i];}
1008 
1009  void SetValue(const Double_t *m)
1010  {
1011  fEl[0][0] = *m++; fEl[1][0] = *m++; fEl[2][0] = *m++; m++;
1012  fEl[0][1] = *m++; fEl[1][1] = *m++; fEl[2][1] = *m++; m++;
1013  fEl[0][2] = *m++; fEl[1][2] = *m++; fEl[2][2] = *m;
1014  }
1015  void SetValue(Double_t xx, Double_t xy, Double_t xz,
1016  Double_t yx, Double_t yy, Double_t yz,
1017  Double_t zx, Double_t zy, Double_t zz)
1018  {
1019  fEl[0][0] = xx; fEl[0][1] = xy; fEl[0][2] = xz;
1020  fEl[1][0] = yx; fEl[1][1] = yy; fEl[1][2] = yz;
1021  fEl[2][0] = zx; fEl[2][1] = zy; fEl[2][2] = zz;
1022  }
1023  void SetEuler(const TVector3 &euler)
1024  {
1025  Double_t ci = TMath::Cos(euler[0]);
1026  Double_t cj = TMath::Cos(euler[1]);
1027  Double_t ch = TMath::Cos(euler[2]);
1028  Double_t si = TMath::Sin(euler[0]);
1029  Double_t sj = TMath::Sin(euler[1]);
1030  Double_t sh = TMath::Sin(euler[2]);
1031  Double_t cc = ci * ch;
1032  Double_t cs = ci * sh;
1033  Double_t sc = si * ch;
1034  Double_t ss = si * sh;
1035  SetValue(cj * ch, sj * sc - cs, sj * cc + ss,
1036  cj * sh, sj * ss + cc, sj * cs - sc,
1037  -sj, cj * si, cj * ci);
1038  }
1039 
1040  void Scale(Double_t x, Double_t y, Double_t z)
1041  {
1042  fEl[0][0] *= x; fEl[0][1] *= y; fEl[0][2] *= z;
1043  fEl[1][0] *= x; fEl[1][1] *= y; fEl[1][2] *= z;
1044  fEl[2][0] *= x; fEl[2][1] *= y; fEl[2][2] *= z;
1045  }
1046  TMatrix3x3 Scaled(Double_t x, Double_t y, Double_t z)const
1047  {
1048  return TMatrix3x3(fEl[0][0] * x, fEl[0][1] * y, fEl[0][2] * z,
1049  fEl[1][0] * x, fEl[1][1] * y, fEl[1][2] * z,
1050  fEl[2][0] * x, fEl[2][1] * y, fEl[2][2] * z);
1051  }
1052 
1053  void SetIdentity()
1054  {
1055  SetValue(1., 0., 0., 0., 1., 0., 0., 0., 1.);
1056  }
1057  void GetValue(Double_t *m)const
1058  {
1059  *m++ = fEl[0][0]; *m++ = fEl[1][0]; *m++ = fEl[2][0]; *m++ = 0.0;
1060  *m++ = fEl[0][1]; *m++ = fEl[1][1]; *m++ = fEl[2][1]; *m++ = 0.0;
1061  *m++ = fEl[0][2]; *m++ = fEl[1][2]; *m++ = fEl[2][2]; *m = 0.0;
1062  }
1063 
1064  TMatrix3x3 &operator *= (const TMatrix3x3 &m);
1065  Double_t Tdot(Int_t c, const TVector3 &v)const
1066  {
1067  return fEl[0][c] * v[0] + fEl[1][c] * v[1] + fEl[2][c] * v[2];
1068  }
1069  Double_t Cofac(Int_t r1, Int_t c1, Int_t r2, Int_t c2)const
1070  {
1071  return fEl[r1][c1] * fEl[r2][c2] - fEl[r1][c2] * fEl[r2][c1];
1072  }
1073 
1074  Double_t Determinant()const;
1075  TMatrix3x3 Adjoint()const;
1076  TMatrix3x3 Absolute()const;
1077  TMatrix3x3 Transposed()const;
1078  void Transpose();
1079  TMatrix3x3 Inverse()const;
1080  void Invert();
1081  };
1082 
1083  /////////////////////////////////////////////////////////////////////////////
1084  ///
1085 
1086  TMatrix3x3 &TMatrix3x3::operator *= (const TMatrix3x3 &m)
1087  {
1088  SetValue(m.Tdot(0, fEl[0]), m.Tdot(1, fEl[0]), m.Tdot(2, fEl[0]),
1089  m.Tdot(0, fEl[1]), m.Tdot(1, fEl[1]), m.Tdot(2, fEl[1]),
1090  m.Tdot(0, fEl[2]), m.Tdot(1, fEl[2]), m.Tdot(2, fEl[2]));
1091  return *this;
1092  }
1093 
1094  /////////////////////////////////////////////////////////////////////////////
1095  ///
1096 
1097  Double_t TMatrix3x3::Determinant()const
1098  {
1099  return triple((*this)[0], (*this)[1], (*this)[2]);
1100  }
1101 
1102  /////////////////////////////////////////////////////////////////////////////
1103  ///
1104 
1105  TMatrix3x3 TMatrix3x3::Absolute()const
1106  {
1107  return TMatrix3x3(TMath::Abs(fEl[0][0]), TMath::Abs(fEl[0][1]), TMath::Abs(fEl[0][2]),
1108  TMath::Abs(fEl[1][0]), TMath::Abs(fEl[1][1]), TMath::Abs(fEl[1][2]),
1109  TMath::Abs(fEl[2][0]), TMath::Abs(fEl[2][1]), TMath::Abs(fEl[2][2]));
1110  }
1111 
1112  /////////////////////////////////////////////////////////////////////////////
1113  ///
1114 
1115  TMatrix3x3 TMatrix3x3::Transposed()const
1116  {
1117  return TMatrix3x3(fEl[0][0], fEl[1][0], fEl[2][0],
1118  fEl[0][1], fEl[1][1], fEl[2][1],
1119  fEl[0][2], fEl[1][2], fEl[2][2]);
1120  }
1121 
1122  /////////////////////////////////////////////////////////////////////////////
1123  ///
1124 
1125  void TMatrix3x3::Transpose()
1126  {
1127  *this = Transposed();
1128  }
1129 
1130  /////////////////////////////////////////////////////////////////////////////
1131  ///
1132 
1133  TMatrix3x3 TMatrix3x3::Adjoint()const
1134  {
1135  return TMatrix3x3(Cofac(1, 1, 2, 2), Cofac(0, 2, 2, 1), Cofac(0, 1, 1, 2),
1136  Cofac(1, 2, 2, 0), Cofac(0, 0, 2, 2), Cofac(0, 2, 1, 0),
1137  Cofac(1, 0, 2, 1), Cofac(0, 1, 2, 0), Cofac(0, 0, 1, 1));
1138  }
1139 
1140  /////////////////////////////////////////////////////////////////////////////
1141  ///
1142 
1143  TMatrix3x3 TMatrix3x3::Inverse()const
1144  {
1145  TVector3 co(Cofac(1, 1, 2, 2), Cofac(1, 2, 2, 0), Cofac(1, 0, 2, 1));
1146  Double_t det = dot((*this)[0], co);
1147  Double_t s = 1. / det;
1148  return TMatrix3x3(co[0] * s, Cofac(0, 2, 2, 1) * s, Cofac(0, 1, 1, 2) * s,
1149  co[1] * s, Cofac(0, 0, 2, 2) * s, Cofac(0, 2, 1, 0) * s,
1150  co[2] * s, Cofac(0, 1, 2, 0) * s, Cofac(0, 0, 1, 1) * s);
1151  }
1152 
1153  /////////////////////////////////////////////////////////////////////////////
1154  ///
1155 
1156  void TMatrix3x3::Invert()
1157  {
1158  *this = Inverse();
1159  }
1160 
1161  /////////////////////////////////////////////////////////////////////////////
1162 
1163  TVector3 operator * (const TMatrix3x3& m, const TVector3& v)
1164  {
1165  return TVector3(dot(m[0], v), dot(m[1], v), dot(m[2], v));
1166  }
1167 
1168  /////////////////////////////////////////////////////////////////////////////
1169 
1170  TVector3 operator * (const TVector3& v, const TMatrix3x3& m)
1171  {
1172  return TVector3(m.Tdot(0, v), m.Tdot(1, v), m.Tdot(2, v));
1173  }
1174 
1175  /////////////////////////////////////////////////////////////////////////////
1176 
1177  TMatrix3x3 operator * (const TMatrix3x3 &m1, const TMatrix3x3 &m2)
1178  {
1179  return TMatrix3x3(m2.Tdot(0, m1[0]), m2.Tdot(1, m1[0]), m2.Tdot(2, m1[0]),
1180  m2.Tdot(0, m1[1]), m2.Tdot(1, m1[1]), m2.Tdot(2, m1[1]),
1181  m2.Tdot(0, m1[2]), m2.Tdot(1, m1[2]), m2.Tdot(2, m1[2]));
1182  }
1183 
1184  /////////////////////////////////////////////////////////////////////////////
1185 
1186  TMatrix3x3 mmult_transpose_left(const TMatrix3x3 &m1, const TMatrix3x3 &m2)
1187  {
1188  return TMatrix3x3(m1[0][0] * m2[0][0] + m1[1][0] * m2[1][0] + m1[2][0] * m2[2][0],
1189  m1[0][0] * m2[0][1] + m1[1][0] * m2[1][1] + m1[2][0] * m2[2][1],
1190  m1[0][0] * m2[0][2] + m1[1][0] * m2[1][2] + m1[2][0] * m2[2][2],
1191  m1[0][1] * m2[0][0] + m1[1][1] * m2[1][0] + m1[2][1] * m2[2][0],
1192  m1[0][1] * m2[0][1] + m1[1][1] * m2[1][1] + m1[2][1] * m2[2][1],
1193  m1[0][1] * m2[0][2] + m1[1][1] * m2[1][2] + m1[2][1] * m2[2][2],
1194  m1[0][2] * m2[0][0] + m1[1][2] * m2[1][0] + m1[2][2] * m2[2][0],
1195  m1[0][2] * m2[0][1] + m1[1][2] * m2[1][1] + m1[2][2] * m2[2][1],
1196  m1[0][2] * m2[0][2] + m1[1][2] * m2[1][2] + m1[2][2] * m2[2][2]);
1197  }
1198 
1199  /////////////////////////////////////////////////////////////////////////////
1200 
1201  TMatrix3x3 mmult_transpose_right(const TMatrix3x3& m1, const TMatrix3x3& m2)
1202  {
1203  return TMatrix3x3(m1[0].Dot(m2[0]), m1[0].Dot(m2[1]), m1[0].Dot(m2[2]),
1204  m1[1].Dot(m2[0]), m1[1].Dot(m2[1]), m1[1].Dot(m2[2]),
1205  m1[2].Dot(m2[0]), m1[2].Dot(m2[1]), m1[2].Dot(m2[2]));
1206  }
1207 
1208  class TLine3 {
1209  private :
1210  Bool_t fBounds[2];
1211  Double_t fParams[2];
1212  TPoint3 fOrigin;
1213  TVector3 fDir;
1214 
1215  public :
1216  TLine3();
1217  TLine3(const TPoint3 &p1, const TPoint3 &p2);
1218  TLine3(const TPoint3 &p1, const TVector3 &v);
1219  TLine3(const TPoint3 &p1, const TVector3 &v, Bool_t bound1, Bool_t bound2);
1220 
1221  static TLine3 InfiniteRay(const TPoint3 &p1, const TVector3 &v);
1222  const TVector3 &Direction()const {return fDir;}
1223  const TPoint3 &Origin()const{ return fOrigin;}
1224 
1225  Bool_t Bounds(Int_t i)const
1226  {
1227  return (i == 0 ? fBounds[0] : fBounds[1]);
1228  }
1229  Bool_t &Bounds(Int_t i)
1230  {
1231  return (i == 0 ? fBounds[0] : fBounds[1]);
1232  }
1233  const Double_t &Param(Int_t i)const
1234  {
1235  return (i == 0 ? fParams[0] : fParams[1]);
1236  }
1237  Double_t &Param(Int_t i)
1238  {
1239  return (i == 0 ? fParams[0] : fParams[1]);
1240  }
1241  TVector3 UnboundSmallestVector(const TPoint3 &point)const
1242  {
1243  TVector3 diff(fOrigin - point);
1244  return diff - fDir * diff.Dot(fDir);
1245  }
1246  Double_t UnboundClosestParameter(const TPoint3 &point)const
1247  {
1248  TVector3 diff(fOrigin-point);
1249  return diff.Dot(fDir);
1250  }
1251  Double_t UnboundDistance(const TPoint3& point)const
1252  {
1253  return UnboundSmallestVector(point).Length();
1254  }
1255  Bool_t IsParameterOnLine(const Double_t &t) const
1256  {
1257  return ((fParams[0] - epsilon < t) || (!fBounds[0])) && ((fParams[1] > t + epsilon) || (!fBounds[1]));
1258  }
1259  };
1260 
1261  /////////////////////////////////////////////////////////////////////////////
1262  ///
1263 
1264  TLine3::TLine3() : fOrigin(0,0,0), fDir(1,0,0)
1265  {
1266  fBounds[0] = kFALSE; fBounds[1] = kFALSE;
1267  fParams[0] = 0; fParams[1] = 1;
1268  }
1269 
1270  /////////////////////////////////////////////////////////////////////////////
1271  ///
1272 
1273  TLine3::TLine3(const TPoint3 &p1, const TPoint3 &p2) : fOrigin(p1), fDir(p2-p1)
1274  {
1275  fBounds[0] = kTRUE; fBounds[1] = kTRUE;
1276  fParams[0] = 0; fParams[1] = 1;
1277  }
1278 
1279  /////////////////////////////////////////////////////////////////////////////
1280  ///
1281 
1282  TLine3::TLine3(const TPoint3 &p1, const TVector3 &v): fOrigin(p1), fDir(v)
1283  {
1284  fBounds[0] = kFALSE; fBounds[1] = kFALSE;
1285  fParams[0] = 0; fParams[1] = 1;
1286  }
1287 
1288  /////////////////////////////////////////////////////////////////////////////
1289  ///
1290 
1291  TLine3::TLine3(const TPoint3 &p1, const TVector3 &v, Bool_t bound1, Bool_t bound2)
1292  : fOrigin(p1), fDir(v)
1293  {
1294  fBounds[0] = bound1; fBounds[1] = bound2;
1295  fParams[0] = 0; fParams[1] = 1;
1296  }
1297 
1298  class TPlane3 : public Tuple4 {
1299  public :
1300  TPlane3(const TVector3 &a, const TVector3 &b, const TVector3 &c);
1301  TPlane3(const TVector3 &n, const TVector3 &p);
1302  TPlane3();
1303  TPlane3(const TPlane3 & p):Tuple4(p){}
1304 
1305  TVector3 Normal()const;
1306  Double_t Scalar()const;
1307  void Invert();
1308  Double_t SignedDistance(const TVector3 &)const;
1309 
1310  TPlane3 &operator = (const TPlane3 & rhs);
1311  };
1312 
1313  /////////////////////////////////////////////////////////////////////////////
1314  ///
1315 
1316  TPlane3::TPlane3(const TVector3 &a, const TVector3 &b, const TVector3 &c)
1317  {
1318  TVector3 l1 = b-a;
1319  TVector3 l2 = c-b;
1320  TVector3 n = l1.Cross(l2);
1321  n = n.SafeNormalized();
1322  Double_t d = n.Dot(a);
1323  fCo[0] = n.X(); fCo[1] = n.Y(); fCo[2] = n.Z(); fCo[3] = -d;
1324  }
1325 
1326  /////////////////////////////////////////////////////////////////////////////
1327  ///
1328 
1329  TPlane3::TPlane3(const TVector3 &n, const TVector3 &p)
1330  {
1331  TVector3 mn = n.SafeNormalized();
1332  Double_t md = mn.Dot(p);
1333  fCo[0] = mn.X(); fCo[1] = mn.Y(); fCo[2] = mn.Z(); fCo[3] = -md;
1334  }
1335 
1336  /////////////////////////////////////////////////////////////////////////////
1337  ///
1338 
1339  TPlane3::TPlane3() : Tuple4()
1340  {
1341  fCo[0] = 1.; fCo[1] = 0.;
1342  fCo[2] = 0.; fCo[3] = 0.;
1343  }
1344 
1345  /////////////////////////////////////////////////////////////////////////////
1346  ///
1347 
1348  TVector3 TPlane3::Normal()const
1349  {
1350  return TVector3(fCo[0],fCo[1],fCo[2]);
1351  }
1352 
1353  /////////////////////////////////////////////////////////////////////////////
1354  ///
1355 
1356  Double_t TPlane3::Scalar()const
1357  {
1358  return fCo[3];
1359  }
1360 
1361  /////////////////////////////////////////////////////////////////////////////
1362  ///
1363 
1364  void TPlane3::Invert()
1365  {
1366  fCo[0] = -fCo[0]; fCo[1] = -fCo[1]; fCo[2] = -fCo[2]; fCo[3] = -fCo[3];
1367  }
1368 
1369  /////////////////////////////////////////////////////////////////////////////
1370  ///
1371 
1372  TPlane3 &TPlane3::operator = (const TPlane3 &rhs)
1373  {
1374  fCo[0] = rhs.fCo[0]; fCo[1] = rhs.fCo[1]; fCo[2] = rhs.fCo[2]; fCo[3] = rhs.fCo[3];
1375  return *this;
1376  }
1377 
1378  /////////////////////////////////////////////////////////////////////////////
1379  ///
1380 
1381  Double_t TPlane3::SignedDistance(const TVector3 &v)const
1382  {
1383  return Normal().Dot(v) + fCo[3];
1384  }
1385 
1386  class TBBox {
1387  friend Bool_t intersect(const TBBox &a, const TBBox &b);
1388 
1389  private:
1390  TPoint3 fCenter;
1391  TVector3 fExtent;
1392 
1393  public:
1394  TBBox(){}
1395  TBBox(const TPoint3 &mini, const TPoint3 &maxi){SetValue(mini,maxi);}
1396 
1397  const TPoint3 &Center()const{return fCenter;}
1398  const TVector3 &Extent()const{return fExtent;}
1399  TPoint3 &Center(){return fCenter;}
1400  TVector3 &Extent(){return fExtent;}
1401 
1402  void SetValue(const TPoint3 &mini,const TPoint3 &maxi)
1403  {
1404  fExtent = (maxi - mini) / 2;
1405  fCenter = mini + fExtent;
1406  }
1407  void Enclose(const TBBox &a, const TBBox &b)
1408  {
1409  TPoint3 lower(
1410  TMath::Min(a.Lower(0), b.Lower(0)),
1411  TMath::Min(a.Lower(1), b.Lower(1)),
1412  TMath::Min(a.Lower(2), b.Lower(2))
1413  );
1414  TPoint3 upper(
1415  TMath::Max(a.Upper(0), b.Upper(0)),
1416  TMath::Max(a.Upper(1), b.Upper(1)),
1417  TMath::Max(a.Upper(2), b.Upper(2))
1418  );
1419  SetValue(lower, upper);
1420  }
1421 
1422  void SetEmpty()
1423  {
1424  fCenter.SetValue(0., 0., 0.);
1425  fExtent.SetValue(-infinity, -infinity, -infinity);
1426  }
1427  void Include(const TPoint3 &p)
1428  {
1429  TPoint3 lower(
1430  TMath::Min(Lower(0), p[0]),
1431  TMath::Min(Lower(1), p[1]),
1432  TMath::Min(Lower(2), p[2])
1433  );
1434  TPoint3 upper(
1435  TMath::Max(Upper(0), p[0]),
1436  TMath::Max(Upper(1), p[1]),
1437  TMath::Max(Upper(2), p[2])
1438  );
1439  SetValue(lower, upper);
1440  }
1441 
1442  void Include(const TBBox &b)
1443  {
1444  Enclose(*this, b);
1445  }
1446  Double_t Lower(Int_t i)const
1447  {
1448  return fCenter[i] - fExtent[i];
1449  }
1450  Double_t Upper(Int_t i)const
1451  {
1452  return fCenter[i] + fExtent[i];
1453  }
1454  TPoint3 Lower()const
1455  {
1456  return fCenter - fExtent;
1457  }
1458  TPoint3 Upper()const
1459  {
1460  return fCenter + fExtent;
1461  }
1462  Double_t Size()const
1463  {
1464  return TMath::Max(TMath::Max(fExtent[0], fExtent[1]), fExtent[2]);
1465  }
1466  Int_t LongestAxis()const
1467  {
1468  return fExtent.ClosestAxis();
1469  }
1470  Bool_t IntersectXRay(const TPoint3 &xBase)const
1471  {
1472  if (xBase[0] <= Upper(0)) {
1473  if (xBase[1] <= Upper(1) && xBase[1] >= Lower(1)) {
1474  if (xBase[2] <= Upper(2) && xBase[2] >= Lower(2)) {return kTRUE;}
1475  }
1476  }
1477  return kFALSE;
1478  }
1479  };
1480 
1481  /////////////////////////////////////////////////////////////////////////////
1482 
1483  Bool_t intersect(const TBBox &a, const TBBox &b)
1484  {
1485  return TMath::Abs(a.fCenter[0] - b.fCenter[0]) <= a.fExtent[0] + b.fExtent[0] &&
1486  TMath::Abs(a.fCenter[1] - b.fCenter[1]) <= a.fExtent[1] + b.fExtent[1] &&
1487  TMath::Abs(a.fCenter[2] - b.fCenter[2]) <= a.fExtent[2] + b.fExtent[2];
1488  }
1489 
1490  class TBBoxNode {
1491  public:
1492  enum ETagType {kLeaf, kInternal};
1493  TBBox fBBox;
1494  ETagType fTag;
1495  };
1496 
1497  class TBBoxLeaf : public TBBoxNode {
1498  public:
1499  Int_t fPolyIndex;
1500 
1501  TBBoxLeaf() : fPolyIndex(0) {}
1502  TBBoxLeaf(Int_t polyIndex, const TBBox &bbox) : fPolyIndex(polyIndex)
1503  {
1504  fBBox = bbox;
1505  fTag = kLeaf;
1506  }
1507  };
1508 
1509  typedef TBBoxLeaf *LeafPtr_t;
1510  typedef TBBoxNode *NodePtr_t;
1511 
1512  class TBBoxInternal : public TBBoxNode {
1513  public:
1514  NodePtr_t fLeftSon;
1515  NodePtr_t fRightSon;
1516  TBBoxInternal() : fLeftSon(0) ,fRightSon(0) {}
1517  TBBoxInternal(Int_t n, LeafPtr_t leafIt);
1518  };
1519 
1520  typedef TBBoxInternal* InternalPtr_t;
1521 
1522  class TBBoxTree {
1523  private:
1524  Int_t fBranch;
1525  LeafPtr_t fLeaves;
1526  InternalPtr_t fInternals;
1527  Int_t fNumLeaves;
1528 
1529  public :
1530  TBBoxTree() : fBranch(0), fLeaves(0), fInternals(0), fNumLeaves(0) {}
1531  NodePtr_t RootNode()const{return fInternals;}
1532  ~TBBoxTree()
1533  {
1534  delete[] fLeaves;
1535  delete[] fInternals;
1536  }
1537  void BuildTree(LeafPtr_t leaves, Int_t numLeaves);
1538 
1539  private :
1540  void RecursiveTreeBuild(Int_t n, LeafPtr_t leafIt);
1541  };
1542 
1543  /////////////////////////////////////////////////////////////////////////////
1544  ///
1545 
1546  TBBoxInternal::TBBoxInternal(Int_t n, LeafPtr_t leafIt) :
1547  fLeftSon(0) ,fRightSon(0)
1548  {
1549  fTag = kInternal;
1550  fBBox.SetEmpty();
1551  for (Int_t i=0;i<n;i++)
1552  fBBox.Include(leafIt[i].fBBox);
1553  }
1554 
1555  /////////////////////////////////////////////////////////////////////////////
1556  ///
1557 
1558  void TBBoxTree::BuildTree(LeafPtr_t leaves, Int_t numLeaves)
1559  {
1560  fBranch = 0;
1561  fLeaves = leaves;
1562  fNumLeaves = numLeaves;
1563  fInternals = new TBBoxInternal[numLeaves];
1564  RecursiveTreeBuild(fNumLeaves,fLeaves);
1565  }
1566 
1567  /////////////////////////////////////////////////////////////////////////////
1568  ///
1569 
1570  void TBBoxTree::RecursiveTreeBuild(Int_t n, LeafPtr_t leafIt)
1571  {
1572  fInternals[fBranch] = TBBoxInternal(n,leafIt);
1573  TBBoxInternal &aBBox = fInternals[fBranch];
1574  fBranch++;
1575 
1576  Int_t axis = aBBox.fBBox.LongestAxis();
1577  Int_t i = 0, mid = n;
1578 
1579  while (i < mid) {
1580  if (leafIt[i].fBBox.Center()[axis] < aBBox.fBBox.Center()[axis]) {
1581  ++i;
1582  } else {
1583  --mid;
1584  std::swap(leafIt[i], leafIt[mid]);
1585  }
1586  }
1587 
1588  if (mid == 0 || mid == n) {
1589  mid = n / 2;
1590  }
1591  if (mid >= 2) {
1592  aBBox.fRightSon = fInternals + fBranch;
1593  RecursiveTreeBuild(mid,leafIt);
1594  } else {
1595  aBBox.fRightSon = leafIt;
1596  }
1597  if (n - mid >= 2) {
1598  aBBox.fLeftSon = fInternals + fBranch;
1599  RecursiveTreeBuild(n - mid, leafIt + mid);
1600  } else {
1601  aBBox.fLeftSon = leafIt + mid;
1602  }
1603  }
1604 
1605  class TBlenderVProp {
1606  private:
1607  Int_t fVertexIndex;
1608 
1609  public:
1610  TBlenderVProp(Int_t vIndex) : fVertexIndex(vIndex){}
1611  TBlenderVProp(Int_t vIndex, const TBlenderVProp &,
1612  const TBlenderVProp &, const Double_t &)
1613  {
1614  fVertexIndex = vIndex;
1615  }
1616  TBlenderVProp() : fVertexIndex(-1){}
1617  operator Int_t()const
1618  {
1619  return fVertexIndex;
1620  }
1621  TBlenderVProp &operator = (Int_t i)
1622  {
1623  fVertexIndex = i; return *this;
1624  }
1625  };
1626 
1627  template <class TMesh>
1628  class TPolygonGeometry {
1629  public:
1630  typedef typename TMesh::Polygon TPolygon;
1631 
1632  private:
1633  const TMesh &fMesh;
1634  const TPolygon &fPoly;
1635  public:
1636  TPolygonGeometry(const TMesh &mesh, Int_t pIndex)
1637  : fMesh(mesh), fPoly(mesh.Polys()[pIndex])
1638  {}
1639  TPolygonGeometry(const TMesh &mesh, const TPolygon &poly)
1640  : fMesh(mesh), fPoly(poly)
1641  {}
1642  const TPoint3 &operator [] (Int_t i)const
1643  {
1644  return fMesh.Verts()[fPoly[i]].Pos();
1645  }
1646  Int_t Size()const
1647  {
1648  return fPoly.Size();
1649  }
1650 
1651  private:
1652  TPolygonGeometry(const TPolygonGeometry &);
1653  TPolygonGeometry& operator = (TPolygonGeometry &);
1654  };
1655 
1656  template <class TPolygon, class TVertex>
1657  class TMesh : public TBaseMesh {
1658  public:
1659  typedef std::vector<TVertex> VLIST;
1660  typedef std::vector<TPolygon> PLIST;
1661  typedef TPolygon Polygon;
1662  typedef TVertex Vertex;
1663  typedef TPolygonGeometry<TMesh> TGBinder;
1664 
1665  private:
1666  VLIST fVerts;
1667  PLIST fPolys;
1668 
1669  public:
1670  VLIST &Verts(){return fVerts;}
1671  const VLIST &Verts()const{return fVerts;}
1672  PLIST &Polys(){return fPolys;}
1673  const PLIST &Polys()const{return fPolys;}
1674 
1675  //TBaseMesh's final-overriders
1676  UInt_t NumberOfPolys()const{return fPolys.size();}
1677  UInt_t NumberOfVertices()const{return fVerts.size();}
1678  UInt_t SizeOfPoly(UInt_t polyIndex)const{return fPolys[polyIndex].Size();}
1679  const Double_t *GetVertex(UInt_t vertexNum)const{return fVerts[vertexNum].GetValue();}
1680 
1681  Int_t GetVertexIndex(UInt_t polyNum, UInt_t vertexNum)const
1682  {
1683  return fPolys[polyNum][vertexNum];
1684  }
1685  };
1686 
1687  const Int_t cofacTable[3][2] = {{1,2}, {0,2}, {0,1}};
1688 
1689  /////////////////////////////////////////////////////////////////////////////
1690 
1691  Bool_t intersect(const TPlane3 &p1, const TPlane3 &p2, TLine3 &output)
1692  {
1693  TMatrix3x3 mat;
1694  mat[0] = p1.Normal();
1695  mat[1] = p2.Normal();
1696  mat[2] = mat[0].Cross(mat[1]);
1697  if (mat[2].FuzzyZero()) return kFALSE;
1698  TVector3 aPoint(-p1.Scalar(),-p2.Scalar(),0);
1699  output = TLine3(TPoint3(0., 0., 0.) + mat.Inverse() * aPoint ,mat[2]);
1700  return kTRUE;
1701  }
1702 
1703  /////////////////////////////////////////////////////////////////////////////
1704 
1705  Bool_t intersect_2d_no_bounds_check(const TLine3 &l1, const TLine3 &l2, Int_t majAxis,
1706  Double_t &l1Param, Double_t &l2Param)
1707  {
1708  Int_t ind1 = cofacTable[majAxis][0];
1709  Int_t ind2 = cofacTable[majAxis][1];
1710  Double_t zX = l2.Origin()[ind1] - l1.Origin()[ind1];
1711  Double_t zY = l2.Origin()[ind2] - l1.Origin()[ind2];
1712  Double_t det = l1.Direction()[ind1]*l2.Direction()[ind2] -
1713  l2.Direction()[ind1]*l1.Direction()[ind2];
1714  if (fuzzy_zero(det)) return kFALSE;
1715  l1Param = (l2.Direction()[ind2] * zX - l2.Direction()[ind1] * zY)/det;
1716  l2Param = -(-l1.Direction()[ind2] * zX + l1.Direction()[ind1] * zY)/det;
1717  return kTRUE;
1718  }
1719 
1720  /////////////////////////////////////////////////////////////////////////////
1721 
1722  Bool_t intersect_2d_bounds_check(const TLine3 &l1, const TLine3 &l2, Int_t majAxis,
1723  Double_t &l1Param, Double_t &l2Param)
1724  {
1725  Bool_t isect = intersect_2d_no_bounds_check(l1, l2, majAxis, l1Param, l2Param);
1726  if (!isect) return kFALSE;
1727  return l1.IsParameterOnLine(l1Param) && l2.IsParameterOnLine(l2Param);
1728  }
1729 
1730  /////////////////////////////////////////////////////////////////////////////
1731 
1733  {
1734  if (TMath::Abs(distance) < epsil) return 0;
1735  else return distance < 0 ? 1 : 2;
1736  }
1737 
1738  /////////////////////////////////////////////////////////////////////////////
1739 
1740  template<typename TGBinder>
1741  Bool_t intersect_poly_with_line_2d(const TLine3 &l, const TGBinder &p1, const TPlane3 &plane,
1742  Double_t &a, Double_t &b)
1743  {
1744  Int_t majAxis = plane.Normal().ClosestAxis();
1745  Int_t lastInd = p1.Size()-1;
1746  b = (-infinity); a = (infinity);
1747  Double_t isectParam(0.), isectParam2(0.);
1748  Int_t i;
1749  Int_t j = lastInd;
1750  Int_t isectsFound(0);
1751  for (i=0;i<=lastInd; j=i,i++ ) {
1752  TLine3 testLine(p1[j],p1[i]);
1753  if (intersect_2d_bounds_check(l, testLine, majAxis, isectParam, isectParam2)) {
1754  ++isectsFound;
1755  b = TMath::Max(isectParam, b);
1756  a = TMath::Min(isectParam, a);
1757  }
1758  }
1759 
1760  return (isectsFound > 0);
1761  }
1762 
1763  /////////////////////////////////////////////////////////////////////////////
1764 
1765  template<typename TGBinder>
1766  Bool_t instersect_poly_with_line_3d(const TLine3 &l, const TGBinder &p1,
1767  const TPlane3 &plane, Double_t &a)
1768  {
1769  Double_t determinant = l.Direction().Dot(plane.Normal());
1770  if (fuzzy_zero(determinant)) return kFALSE;
1771  a = -plane.Scalar() - plane.Normal().Dot(l.Origin());
1772  a /= determinant;
1773  if (a <= 0 ) return kFALSE;
1774  if (!l.IsParameterOnLine(a)) return kFALSE;
1775  TPoint3 pointOnPlane = l.Origin() + l.Direction() * a;
1776  return point_in_polygon_test_3d(p1, plane, l.Origin(), pointOnPlane);
1777  }
1778 
1779  /////////////////////////////////////////////////////////////////////////////
1780 
1781  template<typename TGBinder>
1782  Bool_t point_in_polygon_test_3d(const TGBinder& p1, const TPlane3& plane, const TPoint3& origin,
1783  const TPoint3 &pointOnPlane)
1784  {
1785  Bool_t discardSign = plane.SignedDistance(origin) < 0 ? kTRUE : kFALSE;
1786  const Int_t polySize = p1.Size();
1787  TPoint3 lastPoint = p1[polySize-1];
1788  for (Int_t i=0;i<polySize; ++i) {
1789  const TPoint3& aPoint = p1[i];
1790  TPlane3 testPlane(origin, lastPoint, aPoint);
1791  if ((testPlane.SignedDistance(pointOnPlane) <= 0) == discardSign) return kFALSE;
1792  lastPoint = aPoint;
1793  }
1794 
1795  return kTRUE;
1796  }
1797 
1798  /////////////////////////////////////////////////////////////////////////////
1799 
1800  template <typename TGBinder>
1801  TPoint3 polygon_mid_point(const TGBinder &p1)
1802  {
1803  TPoint3 midPoint(0., 0., 0.);
1804  Int_t i;
1805  for (i=0; i < p1.Size(); i++)
1806  midPoint += p1[i];
1807  return TPoint3(midPoint[0] / i, midPoint[1] / i, midPoint[2] / i);
1808  }
1809 
1810  /////////////////////////////////////////////////////////////////////////////
1811 
1812  template <typename TGBinder>
1813  Int_t which_side(const TGBinder &p1, const TPlane3 &plane1)
1814  {
1815  Int_t output = 0;
1816  Int_t i;
1817  for (i=0; i<p1.Size(); i++) {
1818  Double_t signedDistance = plane1.SignedDistance(p1[i]);
1819  if (!fuzzy_zero(signedDistance))
1820  signedDistance < 0 ? (output |= 1) : (output |=2);
1821  }
1822 
1823  return output;
1824  }
1825 
1826  /////////////////////////////////////////////////////////////////////////////
1827 
1828  template <typename TGBinder>
1829  TLine3 polygon_mid_point_ray(const TGBinder &p1, const TPlane3 &plane)
1830  {
1831  return TLine3(polygon_mid_point(p1),plane.Normal(),kTRUE,kFALSE);
1832  }
1833 
1834  /////////////////////////////////////////////////////////////////////////////
1835 
1836  template <typename TGBinder>
1837  TPlane3 compute_plane(const TGBinder &poly)
1838  {
1839  TPoint3 plast(poly[poly.Size()-1]);
1840  TPoint3 pivot;
1841  TVector3 edge;
1842  Int_t j;
1843  for (j = 0; j < poly.Size(); j++) {
1844  pivot = poly[j];
1845  edge = pivot - plast;
1846  if (!edge.FuzzyZero()) break;
1847  }
1848  for (; j < poly.Size(); j++) {
1849  TVector3 v2 = poly[j] - pivot;
1850  TVector3 v3 = edge.Cross(v2);
1851  if (!v3.FuzzyZero())
1852  return TPlane3(v3,pivot);
1853  }
1854 
1855  return TPlane3();
1856  }
1857 
1858  /////////////////////////////////////////////////////////////////////////////
1859 
1860  template <typename TGBinder>
1861  TBBox fit_bbox(const TGBinder &p1)
1862  {
1863  TBBox bbox; bbox.SetEmpty();
1864  for (Int_t i = 0; i < p1.Size(); ++i)
1865  bbox.Include(p1[i]);
1866  return bbox;
1867  }
1868 
1869  /////////////////////////////////////////////////////////////////////////////
1870 
1871  template<typename TGBinderA, typename TGBinderB>
1872  Bool_t intersect_polygons (const TGBinderA &p1, const TGBinderB &p2,
1873  const TPlane3 &plane1, const TPlane3 &plane2)
1874  {
1875  TLine3 intersectLine;
1876  if (!intersect(plane1, plane2, intersectLine))
1877  return kFALSE;
1878  Double_t p1A, p1B;
1879  Double_t p2A, p2B;
1880  if (
1881  !intersect_poly_with_line_2d(intersectLine,p1,plane1,p1A,p1B) ||
1882  !intersect_poly_with_line_2d(intersectLine,p2,plane2,p2A,p2B))
1883  {
1884  return kFALSE;
1885  }
1886  Double_t maxOMin = TMath::Max(p1A,p2A);
1887  Double_t minOMax = TMath::Min(p1B,p2B);
1888  return (maxOMin <= minOMax);
1889  }
1890 
1891  template <class TMesh, class TSplitFunctionBinder>
1892  class TSplitFunction {
1893  private:
1894  TMesh &fMesh;
1895  TSplitFunctionBinder &fFunctionBinder;
1896 
1897  public:
1898  TSplitFunction(TMesh &mesh, TSplitFunctionBinder &functionBindor)
1899  : fMesh(mesh), fFunctionBinder(functionBindor)
1900  {}
1901  void SplitPolygon(const Int_t p1Index, const TPlane3 &plane,
1902  Int_t &inPiece, Int_t &outPiece,
1903  const Double_t onEpsilon)
1904  {
1905  const typename TMesh::Polygon &p = fMesh.Polys()[p1Index];
1906  typename TMesh::Polygon inP(p),outP(p);
1907  inP.Verts().clear();
1908  outP.Verts().clear();
1909  fFunctionBinder.DisconnectPolygon(p1Index);
1910  Int_t lastIndex = p.Verts().back();
1911  TPoint3 lastVertex = fMesh.Verts()[lastIndex].Pos();
1912  Int_t lastClassification = compute_classification(plane.SignedDistance(lastVertex),onEpsilon);
1913  Int_t totalClassification(lastClassification);
1914  Int_t i;
1915  Int_t j=p.Size()-1;
1916  for (i = 0; i < p.Size(); j = i, ++i)
1917  {
1918  Int_t newIndex = p[i];
1919  TPoint3 aVertex = fMesh.Verts()[newIndex].Pos();
1920  Int_t newClassification = compute_classification(plane.SignedDistance(aVertex),onEpsilon);
1921  if ((newClassification != lastClassification) && newClassification && lastClassification)
1922  {
1923  Int_t newVertexIndex = fMesh.Verts().size();
1924  typedef typename TMesh::Vertex VERTEX_t;
1925  fMesh.Verts().push_back(VERTEX_t());
1926  TVector3 v = aVertex - lastVertex;
1927  Double_t sideA = plane.SignedDistance(lastVertex);
1928  Double_t epsil = -sideA / plane.Normal().Dot(v);
1929  fMesh.Verts().back().Pos() = lastVertex + (v * epsil);
1930  typename TMesh::Polygon::TVProp splitProp(newVertexIndex,p.VertexProps(j),p.VertexProps(i),epsil);
1931  inP.Verts().push_back( splitProp );
1932  outP.Verts().push_back( splitProp );
1933  fFunctionBinder.InsertVertexAlongEdge(lastIndex,newIndex,splitProp);
1934  }
1935  Classify(inP.Verts(),outP.Verts(),newClassification, p.VertexProps(i));
1936  lastClassification = newClassification;
1937  totalClassification |= newClassification;
1938  lastVertex = aVertex;
1939  lastIndex = newIndex;
1940  }
1941  if (totalClassification == 3) {
1942  inPiece = p1Index;
1943  outPiece = fMesh.Polys().size();
1944  fMesh.Polys()[p1Index] = inP;
1945  fMesh.Polys().push_back(outP);
1946  fFunctionBinder.ConnectPolygon(inPiece);
1947  fFunctionBinder.ConnectPolygon(outPiece);
1948  } else {
1949  fFunctionBinder.ConnectPolygon(p1Index);
1950  if (totalClassification == 1) {
1951  inPiece = p1Index;
1952  outPiece = -1;
1953  } else {
1954  outPiece = p1Index;
1955  inPiece = -1;
1956  }
1957  }
1958  }
1959 
1960  void Classify(typename TMesh::Polygon::TVPropList &inGroup,
1961  typename TMesh::Polygon::TVPropList &outGroup,
1962  Int_t classification,
1963  typename TMesh::Polygon::TVProp prop)
1964  {
1965  switch (classification) {
1966  case 0 :
1967  inGroup.push_back(prop);
1968  outGroup.push_back(prop);
1969  break;
1970  case 1 :
1971  inGroup.push_back(prop);
1972  break;
1973  case 2 :
1974  outGroup.push_back(prop);
1975  break;
1976  default :
1977  break;
1978  }
1979  }
1980 
1981  };
1982 
1983  template <typename PROP>
1984  class TDefaultSplitFunctionBinder {
1985  public :
1986  void DisconnectPolygon(Int_t){}
1987  void ConnectPolygon(Int_t){}
1988  void InsertVertexAlongEdge(Int_t, Int_t, const PROP &){}
1989  };
1990 
1991  template <typename TMesh>
1992  class TMeshWrapper {
1993  private :
1994  TMesh &fMesh;
1995 
1996  public:
1997  typedef typename TMesh::Polygon Polygon;
1998  typedef typename TMesh::Vertex Vertex;
1999  typedef typename TMesh::VLIST VLIST;
2000  typedef typename TMesh::PLIST PLIST;
2001  typedef TPolygonGeometry<TMeshWrapper> TGBinder;
2002  typedef TMeshWrapper<TMesh> MyType;
2003 
2004  public:
2005  TMeshWrapper(TMesh &mesh):fMesh(mesh){}
2006 
2007  VLIST &Verts(){return fMesh.Verts();}
2008  const VLIST &Verts()const{return fMesh.Verts();}
2009  PLIST &Polys(){return fMesh.Polys();}
2010  const PLIST &Polys() const {return fMesh.Polys();}
2011 
2012  void ComputePlanes();
2013  TBBox ComputeBBox()const;
2014  void SplitPolygon(Int_t p1Index, const TPlane3 &plane,
2015  Int_t &inPiece, Int_t &outPiece, Double_t onEpsilon);
2016  };
2017 
2018  /////////////////////////////////////////////////////////////////////////////
2019  ///
2020 
2021  template <typename TMesh>
2022  void TMeshWrapper<TMesh>::ComputePlanes()
2023  {
2024  PLIST& polyList = Polys();
2025  UInt_t i;
2026  for (i=0;i < polyList.size(); i++) {
2027  TGBinder binder(*this, i);
2028  polyList[i].SetPlane(compute_plane(binder));
2029  }
2030  }
2031 
2032  /////////////////////////////////////////////////////////////////////////////
2033  ///
2034 
2035  template <typename TMesh>
2036  TBBox TMeshWrapper<TMesh>::ComputeBBox()const
2037  {
2038  const VLIST &vertexList = Verts();
2039  TBBox bbox;
2040  bbox.SetEmpty();
2041  Int_t i;
2042  for (i=0;i<vertexList.size(); i++)
2043  bbox.Include(vertexList[i].Pos());
2044  return bbox;
2045  }
2046 
2047  /////////////////////////////////////////////////////////////////////////////
2048 
2049  template<typename TMesh>
2050  void TMeshWrapper<TMesh>::SplitPolygon(Int_t p1Index, const TPlane3 &plane,
2051  Int_t &inPiece, Int_t &outPiece,
2052  Double_t onEpsilon)
2053  {
2054  typedef typename TMesh::Polygon::TVProp mesh;
2055  TDefaultSplitFunctionBinder<mesh> defaultSplitFunction;
2056  TSplitFunction<MyType,TDefaultSplitFunctionBinder<mesh> >
2057  splitFunction(*this,defaultSplitFunction);
2058  splitFunction.SplitPolygon(p1Index,plane,inPiece,outPiece,onEpsilon);
2059  }
2060 
2061  template <typename AVProp, typename AFProp>
2062  class TPolygonBase {
2063  public:
2064  typedef AVProp TVProp;
2065  typedef AFProp TFProp;
2066  typedef std::vector<TVProp> TVPropList;
2067  typedef typename TVPropList::iterator TVPropIt;
2068 
2069  private :
2070  TVPropList fVerts;
2071  TPlane3 fPlane;
2072  TFProp fFaceProp;
2073  Int_t fClassification;
2074 
2075  public:
2076  const TVPropList &Verts()const{return fVerts;}
2077  TVPropList &Verts(){return fVerts;}
2078  Int_t Size()const{return Int_t(fVerts.size());}
2079 
2080  Int_t operator[](Int_t i) const {return fVerts[i];}
2081 
2082  const TVProp &VertexProps(Int_t i)const{return fVerts[i];}
2083  TVProp &VertexProps(Int_t i){return fVerts[i];}
2084  void SetPlane(const TPlane3 &plane){fPlane = plane;}
2085  const TPlane3 &Plane()const{return fPlane;}
2086  TVector3 Normal()const{return fPlane.Normal();}
2087  Int_t &Classification(){ return fClassification;}
2088  const Int_t &Classification()const{return fClassification;}
2089 
2090  void Reverse()
2091  {
2092  std::reverse(fVerts.begin(),fVerts.end());
2093  fPlane.Invert();
2094  }
2095 
2096  TFProp &FProp(){return fFaceProp;}
2097  const TFProp &FProp()const{return fFaceProp;}
2098  void AddProp(const TVProp &prop){fVerts.push_back(prop);}
2099  };
2100 
2101  typedef std::vector<Int_t> PIndexList_t;
2102  typedef PIndexList_t::iterator PIndexIt_t;
2103  typedef std::vector< PIndexList_t > OverlapTable_t;
2104 
2105  template <typename TMesh>
2106  class TreeIntersector {
2107  private:
2108  OverlapTable_t *fAoverlapsB;
2109  OverlapTable_t *fBoverlapsA;
2110  const TMesh *fMeshA;
2111  const TMesh *fMeshB;
2112 
2113  public :
2114  TreeIntersector(const TBBoxTree &a, const TBBoxTree &b,
2115  OverlapTable_t *aOverlapsB, OverlapTable_t *bOverlapsA,
2116  const TMesh *meshA, const TMesh *meshB)
2117  {
2118  fAoverlapsB = aOverlapsB;
2119  fBoverlapsA = bOverlapsA;
2120  fMeshA = meshA;
2121  fMeshB = meshB;
2122  MarkIntersectingPolygons(a.RootNode(),b.RootNode());
2123  }
2124 
2125  private :
2126  void MarkIntersectingPolygons(const TBBoxNode *a, const TBBoxNode *b)
2127  {
2128  if (!intersect(a->fBBox, b->fBBox)) return;
2129  if (a->fTag == TBBoxNode::kLeaf && b->fTag == TBBoxNode::kLeaf) {
2130  const TBBoxLeaf *la = (const TBBoxLeaf *)a;
2131  const TBBoxLeaf *lb = (const TBBoxLeaf *)b;
2132 
2133  TPolygonGeometry<TMesh> pg1(*fMeshA,la->fPolyIndex);
2134  TPolygonGeometry<TMesh> pg2(*fMeshB,lb->fPolyIndex);
2135 
2136  if (intersect_polygons(pg1, pg2, fMeshA->Polys()[la->fPolyIndex].Plane(),
2137  fMeshB->Polys()[lb->fPolyIndex].Plane())) {
2138  (*fAoverlapsB)[lb->fPolyIndex].push_back(la->fPolyIndex);
2139  (*fBoverlapsA)[la->fPolyIndex].push_back(lb->fPolyIndex);
2140  }
2141  } else if ( a->fTag == TBBoxNode::kLeaf || (b->fTag != TBBoxNode::kLeaf && a->fBBox.Size() < b->fBBox.Size()))
2142  {
2143  MarkIntersectingPolygons(a, ((const TBBoxInternal *)b)->fLeftSon);
2144  MarkIntersectingPolygons(a, ((const TBBoxInternal *)b)->fRightSon);
2145  } else {
2146  MarkIntersectingPolygons(((const TBBoxInternal *)a)->fLeftSon, b);
2147  MarkIntersectingPolygons(((const TBBoxInternal *)a)->fRightSon, b);
2148  }
2149  }
2150  };
2151 
2152  template<typename TMesh>
2153  class TRayTreeIntersector {
2154  private:
2155  const TMesh *fMeshA;
2156  Double_t fLastIntersectValue;
2157  Int_t fPolyIndex;
2158 
2159  public:
2160  TRayTreeIntersector(const TBBoxTree &a, const TMesh *meshA, const TLine3 &xRay, Int_t &polyIndex)
2161  : fMeshA(meshA), fLastIntersectValue(infinity), fPolyIndex(-1)
2162  {
2163  FindIntersectingPolygons(a.RootNode(),xRay);
2164  polyIndex = fPolyIndex;
2165  }
2166 
2167  private :
2168  void FindIntersectingPolygons(const TBBoxNode*a,const TLine3& xRay)
2169  {
2170  if ((xRay.Origin().X() + fLastIntersectValue < a->fBBox.Lower(0)) ||!a->fBBox.IntersectXRay(xRay.Origin()))
2171  return;
2172  if (a->fTag == TBBoxNode::kLeaf) {
2173  const TBBoxLeaf *la = (const TBBoxLeaf *)a;
2174  Double_t testParameter(0.);
2175  TPolygonGeometry<TMesh> pg(*fMeshA, la->fPolyIndex);
2176  if (instersect_poly_with_line_3d(xRay,pg,fMeshA->Polys()[la->fPolyIndex].Plane(),testParameter))
2177  {
2178  if (testParameter < fLastIntersectValue) {
2179  fLastIntersectValue = testParameter;
2180  fPolyIndex = la->fPolyIndex;
2181  }
2182  }
2183  } else {
2184  FindIntersectingPolygons(((const TBBoxInternal*)a)->fLeftSon, xRay);
2185  FindIntersectingPolygons(((const TBBoxInternal*)a)->fRightSon, xRay);
2186  }
2187  }
2188  };
2189 
2190  class TVertexBase {
2191  protected:
2192  Int_t fVertexMap;
2193  TPoint3 fPos;
2194 
2195  public:
2196  TVertexBase(Double_t x, Double_t y, Double_t z) : fVertexMap(-1), fPos(x, y, z){}
2197  TVertexBase():fVertexMap(-1) {}
2198 
2199  const TPoint3 &Pos()const{return fPos;}
2200  TPoint3 &Pos(){return fPos;}
2201  Int_t &VertexMap(){return fVertexMap;}
2202  const Int_t &VertexMap()const{return fVertexMap;}
2203  const Double_t * GetValue()const{return fPos.GetValue();}
2204 
2205  Double_t operator [] (Int_t ind)const{return fPos[ind];}
2206  };
2207 
2208  class TCVertex : public TVertexBase {
2209  private:
2210  PIndexList_t fPolygons;
2211  public:
2212  TCVertex() : TVertexBase(){}
2213  TCVertex(const TVertexBase& vertex) : TVertexBase(vertex){}
2214 
2215  TCVertex &operator = (const TVertexBase &other)
2216  {
2217  fPos= other.Pos();
2218  return *this;
2219  }
2220  const PIndexList_t &Polys()const{return fPolygons;}
2221  PIndexList_t &Polys(){return fPolygons;}
2222 
2223  Int_t &operator [] (Int_t i) { return fPolygons[i];}
2224  const Int_t &operator [] (Int_t i)const{return fPolygons[i];}
2225 
2226  void AddPoly(Int_t polyIndex){fPolygons.push_back(polyIndex);}
2227  void RemovePolygon(Int_t polyIndex)
2228  {
2229  PIndexIt_t foundIt = std::find(fPolygons.begin(), fPolygons.end(), polyIndex);
2230  if (foundIt != fPolygons.end()) {
2231  std::swap(fPolygons.back(), *foundIt);
2232  fPolygons.pop_back();
2233  }
2234  }
2235  };
2236 
2237  template <typename TMesh>
2238  class TConnectedMeshWrapper {
2239  private:
2240  TMesh &fMesh;
2241  UInt_t fUniqueEdgeTestId;
2242  public:
2243  typedef typename TMesh::Polygon Polygon;
2244  typedef typename TMesh::Vertex Vertex;
2245  typedef typename TMesh::Polygon::TVProp VProp;
2246  typedef typename TMesh::VLIST VLIST;
2247  typedef typename TMesh::PLIST PLIST;
2248  typedef TPolygonGeometry<TConnectedMeshWrapper> TGBinder;
2249  typedef TConnectedMeshWrapper<TMesh> MyType;
2250 
2251  TConnectedMeshWrapper(TMesh &mesh) : fMesh(mesh), fUniqueEdgeTestId(0){}
2252 
2253  VLIST &Verts(){return fMesh.Verts();}
2254  const VLIST &Verts()const{return fMesh.Verts();}
2255  PLIST &Polys() {return fMesh.Polys();}
2256  const PLIST &Polys() const {return fMesh.Polys();}
2257  void BuildVertexPolyLists();
2258  void DisconnectPolygon(Int_t polyIndex);
2259  void ConnectPolygon(Int_t polyIndex);
2260  //return the polygons neibouring the given edge.
2261  void EdgePolygons(Int_t v1, Int_t v2, PIndexList_t &polys);
2262  void InsertVertexAlongEdge(Int_t v1,Int_t v2, const VProp &prop);
2263  void SplitPolygon(Int_t p1Index, const TPlane3 &plane, Int_t &inPiece,
2264  Int_t &outPiece, Double_t onEpsilon);
2265  };
2266 
2267  template <class CMesh> class TSplitFunctionBinder {
2268  private:
2269  CMesh &fMesh;
2270 
2271  public:
2272  TSplitFunctionBinder(CMesh &mesh):fMesh(mesh){}
2273  void DisconnectPolygon(Int_t polyIndex){fMesh.DisconnectPolygon(polyIndex);}
2274  void ConnectPolygon(Int_t polygonIndex){fMesh.ConnectPolygon(polygonIndex);}
2275  void InsertVertexAlongEdge(Int_t lastIndex, Int_t newIndex, const typename CMesh::VProp &prop)
2276  {
2277  fMesh.InsertVertexAlongEdge(lastIndex, newIndex,prop);
2278  }
2279  };
2280 
2281  /////////////////////////////////////////////////////////////////////////////
2282  ///
2283 
2284  template <typename TMesh>
2285  void TConnectedMeshWrapper<TMesh>::BuildVertexPolyLists()
2286  {
2287  UInt_t i;
2288  for (i=0; i < Polys().size(); i++)
2289  ConnectPolygon(i);
2290  }
2291 
2292  /////////////////////////////////////////////////////////////////////////////
2293  ///
2294 
2295  template <typename TMesh>
2296  void TConnectedMeshWrapper<TMesh>::DisconnectPolygon(Int_t polyIndex)
2297  {
2298  const Polygon &poly = Polys()[polyIndex];
2299  UInt_t j;
2300  for (j=0;j<poly.Verts().size(); j++) {
2301  Verts()[poly[j]].RemovePolygon(polyIndex);
2302  }
2303  }
2304 
2305  /////////////////////////////////////////////////////////////////////////////
2306  ///
2307 
2308  template <typename TMesh>
2309  void TConnectedMeshWrapper<TMesh>::ConnectPolygon(Int_t polyIndex)
2310  {
2311  const Polygon &poly = Polys()[polyIndex];
2312  UInt_t j;
2313  for (j=0;j<poly.Verts().size(); j++) {
2314  Verts()[poly[j]].AddPoly(polyIndex);
2315  }
2316  }
2317 
2318  /////////////////////////////////////////////////////////////////////////////
2319  ///
2320 
2321  template <typename TMesh>
2322  void TConnectedMeshWrapper<TMesh>::EdgePolygons(Int_t v1, Int_t v2, PIndexList_t &polys)
2323  {
2324  ++fUniqueEdgeTestId;
2325  Vertex &vb1 = Verts()[v1];
2326  UInt_t i;
2327  for (i=0;i < vb1.Polys().size(); ++i){Polys()[vb1[i]].Classification() = fUniqueEdgeTestId;}
2328  Vertex &vb2 = Verts()[v2];
2329  UInt_t j;
2330  for (j=0;j < vb2.Polys().size(); ++j) {
2331  if ((UInt_t)Polys()[vb2[j]].Classification() == fUniqueEdgeTestId) {
2332  polys.push_back(vb2[j]);
2333  }
2334  }
2335  }
2336 
2337  /////////////////////////////////////////////////////////////////////////////
2338  ///
2339 
2340  template <typename TMesh>
2341  void TConnectedMeshWrapper<TMesh>::InsertVertexAlongEdge(Int_t v1, Int_t v2, const VProp &prop)
2342  {
2343  PIndexList_t npolys;
2344  EdgePolygons(v1,v2,npolys);
2345  Int_t newVertex = Int_t(prop);
2346  UInt_t i;
2347  for (i=0;i < npolys.size(); i++) {
2348  typename Polygon::TVPropList& polyVerts = Polys()[npolys[i]].Verts();
2349  typename Polygon::TVPropIt v1pos = std::find(polyVerts.begin(),polyVerts.end(),v1);
2350  if (v1pos != polyVerts.end()) {
2351  typename Polygon::TVPropIt prevPos = (v1pos == polyVerts.begin()) ? polyVerts.end()-1 : v1pos-1;
2352  typename Polygon::TVPropIt nextPos = (v1pos == polyVerts.end()-1) ? polyVerts.begin() : v1pos+1;
2353  if (*prevPos == v2) {
2354  polyVerts.insert(v1pos,prop);
2355  } else if (*nextPos == v2) {
2356  polyVerts.insert(nextPos, prop);
2357  } else {
2358  //assert(kFALSE);
2359  }
2360  Verts()[newVertex].AddPoly(npolys[i]);
2361  } else {
2362  //assert(kFALSE);
2363  }
2364  }
2365  }
2366 
2367  /////////////////////////////////////////////////////////////////////////////
2368 
2369  template <typename TMesh>
2370  void TConnectedMeshWrapper<TMesh>::SplitPolygon(Int_t p1Index, const TPlane3 &plane,
2371  Int_t &inPiece, Int_t &outPiece,
2372  Double_t onEpsilon)
2373  {
2374  TSplitFunctionBinder<MyType> functionBindor(*this);
2375  TSplitFunction<MyType,TSplitFunctionBinder<MyType> > splitFunction(*this,functionBindor);
2376  splitFunction.SplitPolygon(p1Index, plane, inPiece, outPiece, onEpsilon);
2377  }
2378 
2379  struct NullType_t{};
2380  //Original TestPolygon_t has two parameters, the second is face property
2381 
2382  typedef TPolygonBase<TBlenderVProp, NullType_t> TestPolygon_t;
2383  typedef TMesh<TestPolygon_t,TVertexBase> AMesh_t;
2384  typedef TMesh<TestPolygon_t,TCVertex > AConnectedMesh_t;
2385  typedef TMeshWrapper<AMesh_t> AMeshWrapper_t;
2386  typedef TConnectedMeshWrapper<AConnectedMesh_t> AConnectedMeshWrapper_t;
2387 
2388  /////////////////////////////////////////////////////////////////////////////
2389 
2390  template <class TMesh>
2391  void build_split_group(const TMesh &meshA, const TMesh &meshB,
2392  const TBBoxTree &treeA, const TBBoxTree &treeB,
2393  OverlapTable_t &aOverlapsB, OverlapTable_t &bOverlapsA)
2394  {
2395  aOverlapsB = OverlapTable_t(meshB.Polys().size());
2396  bOverlapsA = OverlapTable_t(meshA.Polys().size());
2397  TreeIntersector<TMesh>(treeA, treeB, &aOverlapsB, &bOverlapsA, &meshA, &meshB);
2398  }
2399 
2400  /////////////////////////////////////////////////////////////////////////////
2401 
2402  template <class CMesh, class TMesh>
2403  void partition_mesh(CMesh &mesh, const TMesh &mesh2, const OverlapTable_t &table)
2404  {
2405  UInt_t i;
2406  Double_t onEpsilon(1e-4);
2407  for (i = 0; i < table.size(); i++) {
2408  if (table[i].size()) {
2409  PIndexList_t fragments;
2410  fragments.push_back(i);
2411  UInt_t j;
2412  for (j =0 ; j <table[i].size(); ++j) {
2413  PIndexList_t newFragments;
2414  TPlane3 splitPlane = mesh2.Polys()[table[i][j]].Plane();
2415  UInt_t k;
2416  for (k = 0; k < fragments.size(); ++k) {
2417  Int_t newInFragment;
2418  Int_t newOutFragment;
2419  typename CMesh::TGBinder pg1(mesh,fragments[k]);
2420  typename TMesh::TGBinder pg2(mesh2,table[i][j]);
2421  const TPlane3 &fragPlane = mesh.Polys()[fragments[k]].Plane();
2422 
2423  if (intersect_polygons(pg1,pg2,fragPlane,splitPlane)) {
2424  mesh.SplitPolygon(fragments[k], splitPlane, newInFragment, newOutFragment, onEpsilon);
2425  if (-1 != newInFragment) newFragments.push_back(newInFragment);
2426  if (-1 != newOutFragment) newFragments.push_back(newOutFragment);
2427  } else {
2428  newFragments.push_back(fragments[k]);
2429  }
2430  }
2431  fragments = newFragments;
2432  }
2433  }
2434  }
2435  }
2436 
2437  /////////////////////////////////////////////////////////////////////////////
2438 
2439  template <typename CMesh, typename TMesh>
2440  void classify_mesh(const TMesh &meshA, const TBBoxTree &aTree, CMesh &meshB)
2441  {
2442  UInt_t i;
2443  for (i = 0; i < meshB.Polys().size(); i++) {
2444  typename CMesh::TGBinder pg(meshB,i);
2445  TLine3 midPointRay = polygon_mid_point_ray(pg,meshB.Polys()[i].Plane());
2446  TLine3 midPointXRay(midPointRay.Origin(),TVector3(1,0,0));
2447  Int_t aPolyIndex(-1);
2448  TRayTreeIntersector<TMesh>(aTree,&meshA,midPointXRay,aPolyIndex);
2449  if (-1 != aPolyIndex) {
2450  if (meshA.Polys()[aPolyIndex].Plane().SignedDistance(midPointXRay.Origin()) < 0) {
2451  meshB.Polys()[i].Classification()= 1;
2452  } else {
2453  meshB.Polys()[i].Classification()= 2;
2454  }
2455  } else {
2456  meshB.Polys()[i].Classification()= 2;
2457  }
2458  }
2459  }
2460 
2461  /////////////////////////////////////////////////////////////////////////////
2462 
2463  template <typename CMesh, typename TMesh>
2464  void extract_classification(CMesh &meshA, TMesh &newMesh, Int_t classification, Bool_t reverse)
2465  {
2466  UInt_t i;
2467  for (i = 0; i < meshA.Polys().size(); ++i) {
2468  typename CMesh::Polygon &meshAPolygon = meshA.Polys()[i];
2469  if (meshAPolygon.Classification() == classification) {
2470  newMesh.Polys().push_back(meshAPolygon);
2471  typename TMesh::Polygon &newPolygon = newMesh.Polys().back();
2472  if (reverse) newPolygon.Reverse();
2473  Int_t j;
2474  for (j=0; j< newPolygon.Size(); j++) {
2475  if (meshA.Verts()[newPolygon[j]].VertexMap() == -1) {
2476  newMesh.Verts().push_back(meshA.Verts()[newPolygon[j]]);
2477  meshA.Verts()[newPolygon[j]].VertexMap() = newMesh.Verts().size() -1;
2478  }
2479  newPolygon.VertexProps(j) = meshA.Verts()[newPolygon[j]].VertexMap();
2480  }
2481  }
2482  }
2483  }
2484 
2485  /////////////////////////////////////////////////////////////////////////////
2486 
2487  template <typename MeshA, typename MeshB>
2488  void copy_mesh(const MeshA &source, MeshB &output)
2489  {
2490  Int_t vertexNum = source.Verts().size();
2491  Int_t polyNum = source.Polys().size();
2492 
2493  typedef typename MeshB::VLIST VLIST_t;
2494  typedef typename MeshB::PLIST PLIST_t;
2495 
2496  output.Verts() = VLIST_t(vertexNum);
2497  output.Polys() = PLIST_t(polyNum);
2498 
2499  std::copy(source.Verts().begin(), source.Verts().end(), output.Verts().begin());
2500  std::copy(source.Polys().begin(), source.Polys().end(), output.Polys().begin());
2501  }
2502 
2503  /////////////////////////////////////////////////////////////////////////////
2504 
2505  void build_tree(const AMesh_t &mesh, TBBoxTree &tree)
2506  {
2507  Int_t numLeaves = mesh.Polys().size();
2508  TBBoxLeaf *aLeaves = new TBBoxLeaf[numLeaves];
2509  UInt_t i;
2510  for (i=0;i< mesh.Polys().size(); i++) {
2511  TPolygonGeometry<AMesh_t> pg(mesh,i);
2512  aLeaves[i] = TBBoxLeaf(i, fit_bbox(pg));
2513  }
2514  tree.BuildTree(aLeaves,numLeaves);
2515  }
2516 
2517  /////////////////////////////////////////////////////////////////////////////
2518 
2519  void extract_classification_preserve(const AMesh_t &meshA,
2520  const AMesh_t &meshB,
2521  const TBBoxTree &aTree,
2522  const TBBoxTree &bTree,
2523  const OverlapTable_t &aOverlapsB,
2524  const OverlapTable_t &bOverlapsA,
2525  Int_t aClassification,
2526  Int_t bClassification,
2527  Bool_t reverseA,
2528  Bool_t reverseB,
2529  AMesh_t &output)
2530  {
2531  AConnectedMesh_t meshAPartitioned;
2532  AConnectedMesh_t meshBPartitioned;
2533  copy_mesh(meshA,meshAPartitioned);
2534  copy_mesh(meshB,meshBPartitioned);
2535  AConnectedMeshWrapper_t meshAWrapper(meshAPartitioned);
2536  AConnectedMeshWrapper_t meshBWrapper(meshBPartitioned);
2537  meshAWrapper.BuildVertexPolyLists();
2538  meshBWrapper.BuildVertexPolyLists();
2539  partition_mesh(meshAWrapper, meshB, bOverlapsA);
2540  partition_mesh(meshBWrapper, meshA, aOverlapsB);
2541  classify_mesh(meshB, bTree, meshAPartitioned);
2542  classify_mesh(meshA, aTree, meshBPartitioned);
2543  extract_classification(meshAPartitioned, output, aClassification, reverseA);
2544  extract_classification(meshBPartitioned, output, bClassification, reverseB);
2545  }
2546 
2547  /////////////////////////////////////////////////////////////////////////////
2548 
2549  void extract_classification(const AMesh_t &meshA,
2550  const AMesh_t &meshB,
2551  const TBBoxTree &aTree,
2552  const TBBoxTree &bTree,
2553  const OverlapTable_t &aOverlapsB,
2554  const OverlapTable_t &bOverlapsA,
2555  Int_t aClassification,
2556  Int_t bClassification,
2557  Bool_t reverseA,
2558  Bool_t reverseB,
2559  AMesh_t &output)
2560  {
2561  AMesh_t meshAPartitioned(meshA);
2562  AMesh_t meshBPartitioned(meshB);
2563  AMeshWrapper_t meshAWrapper(meshAPartitioned);
2564  AMeshWrapper_t meshBWrapper(meshBPartitioned);
2565  partition_mesh(meshAWrapper, meshB, bOverlapsA);
2566  partition_mesh(meshBWrapper, meshA, aOverlapsB);
2567  classify_mesh(meshB, bTree, meshAPartitioned);
2568  classify_mesh(meshA, aTree, meshBPartitioned);
2569  extract_classification(meshAPartitioned, output, aClassification, reverseA);
2570  extract_classification(meshBPartitioned, output, bClassification, reverseB);
2571  }
2572 
2573  /////////////////////////////////////////////////////////////////////////////
2574 
2575  AMesh_t *build_intersection(const AMesh_t &meshA, const AMesh_t &meshB, Bool_t preserve)
2576  {
2577  TBBoxTree aTree, bTree;
2578  build_tree(meshA, aTree);
2579  build_tree(meshB, bTree);
2580  OverlapTable_t bOverlapsA(meshA.Polys().size());
2581  OverlapTable_t aOverlapsB(meshB.Polys().size());
2582  build_split_group(meshA, meshB, aTree, bTree, aOverlapsB, bOverlapsA);
2583  AMesh_t *output = new AMesh_t;
2584  if (preserve) {
2586  meshA, meshB, aTree, bTree,
2587  aOverlapsB, bOverlapsA,
2588  1, 1, kFALSE, kFALSE, *output
2589  );
2590  } else {
2592  meshA, meshB, aTree, bTree,
2593  aOverlapsB, bOverlapsA,
2594  1, 1, kFALSE, kFALSE, *output
2595  );
2596  }
2597  return output;
2598  }
2599 
2600  /////////////////////////////////////////////////////////////////////////////
2601 
2602  AMesh_t *build_union(const AMesh_t &meshA, const AMesh_t &meshB, Bool_t preserve)
2603  {
2604  TBBoxTree aTree, bTree;
2605  build_tree(meshA, aTree);
2606  build_tree(meshB, bTree);
2607  OverlapTable_t bOverlapsA(meshA.Polys().size());
2608  OverlapTable_t aOverlapsB(meshB.Polys().size());
2609  build_split_group(meshA, meshB, aTree, bTree, aOverlapsB, bOverlapsA);
2610  AMesh_t *output = new AMesh_t;
2611  if (preserve) {
2613  meshA, meshB, aTree, bTree,
2614  aOverlapsB, bOverlapsA,
2615  2, 2, kFALSE, kFALSE, *output
2616  );
2617  } else {
2619  meshA, meshB, aTree, bTree,
2620  aOverlapsB, bOverlapsA,
2621  2, 2, kFALSE, kFALSE, *output
2622  );
2623  }
2624  return output;
2625  }
2626 
2627  /////////////////////////////////////////////////////////////////////////////
2628 
2629  AMesh_t *build_difference(const AMesh_t &meshA, const AMesh_t &meshB, Bool_t preserve)
2630  {
2631  TBBoxTree aTree, bTree;
2632  build_tree(meshA, aTree);
2633  build_tree(meshB, bTree);
2634  OverlapTable_t bOverlapsA(meshA.Polys().size());
2635  OverlapTable_t aOverlapsB(meshB.Polys().size());
2636  build_split_group(meshA, meshB, aTree, bTree, aOverlapsB, bOverlapsA);
2637  AMesh_t *output = new AMesh_t;
2638  if (preserve) {
2640  meshA, meshB, aTree, bTree,
2641  aOverlapsB, bOverlapsA,
2642  2, 1, kFALSE, kTRUE, *output
2643  );
2644  } else {
2646  meshA, meshB, aTree, bTree,
2647  aOverlapsB, bOverlapsA,
2648  2, 1, kFALSE, kTRUE, *output
2649  );
2650  }
2651  return output;
2652  }
2653 
2654  /////////////////////////////////////////////////////////////////////////////
2655 
2657  {
2658  AMesh_t *newMesh = new AMesh_t;
2659  const Double_t *v = buff.fPnts;
2660 
2661  newMesh->Verts().resize(buff.NbPnts());
2662 
2663  for (UInt_t i = 0; i < buff.NbPnts(); ++i)
2664  newMesh->Verts()[i] = TVertexBase(v[i * 3], v[i * 3 + 1], v[i * 3 + 2]);
2665 
2666  const Int_t *segs = buff.fSegs;
2667  const Int_t *pols = buff.fPols;
2668 
2669  newMesh->Polys().resize(buff.NbPols());
2670 
2671  for (UInt_t numPol = 0, j = 1; numPol < buff.NbPols(); ++numPol) {
2672  TestPolygon_t &currPoly = newMesh->Polys()[numPol];
2673  Int_t segmentInd = pols[j] + j;
2674  Int_t segmentCol = pols[j];
2675  Int_t s1 = pols[segmentInd];
2676  segmentInd--;
2677  Int_t s2 = pols[segmentInd];
2678  segmentInd--;
2679  Int_t segEnds[] = {segs[s1 * 3 + 1], segs[s1 * 3 + 2],
2680  segs[s2 * 3 + 1], segs[s2 * 3 + 2]};
2681  Int_t numPnts[3] = {0};
2682 
2683  if (segEnds[0] == segEnds[2]) {
2684  numPnts[0] = segEnds[1], numPnts[1] = segEnds[0], numPnts[2] = segEnds[3];
2685  } else if (segEnds[0] == segEnds[3]) {
2686  numPnts[0] = segEnds[1], numPnts[1] = segEnds[0], numPnts[2] = segEnds[2];
2687  } else if (segEnds[1] == segEnds[2]) {
2688  numPnts[0] = segEnds[0], numPnts[1] = segEnds[1], numPnts[2] = segEnds[3];
2689  } else {
2690  numPnts[0] = segEnds[0], numPnts[1] = segEnds[1], numPnts[2] = segEnds[2];
2691  }
2692 
2693  currPoly.AddProp(TBlenderVProp(numPnts[0]));
2694  currPoly.AddProp(TBlenderVProp(numPnts[1]));
2695  currPoly.AddProp(TBlenderVProp(numPnts[2]));
2696 
2697  Int_t lastAdded = numPnts[2];
2698 
2699  Int_t end = j + 1;
2700  for (; segmentInd != end; segmentInd--) {
2701  segEnds[0] = segs[pols[segmentInd] * 3 + 1];
2702  segEnds[1] = segs[pols[segmentInd] * 3 + 2];
2703  if (segEnds[0] == lastAdded) {
2704  currPoly.AddProp(TBlenderVProp(segEnds[1]));
2705  lastAdded = segEnds[1];
2706  } else {
2707  currPoly.AddProp(TBlenderVProp(segEnds[0]));
2708  lastAdded = segEnds[0];
2709  }
2710  }
2711  j += segmentCol + 2;
2712  }
2713 
2714  AMeshWrapper_t wrap(*newMesh);
2715 
2716  wrap.ComputePlanes();
2717 
2718  return newMesh;
2719  }
2720 
2721  /////////////////////////////////////////////////////////////////////////////
2722 
2724  {
2725  return build_union(*static_cast<const AMesh_t *>(l), *static_cast<const AMesh_t *>(r), kFALSE);
2726  }
2727 
2728  /////////////////////////////////////////////////////////////////////////////
2729 
2731  {
2732  return build_intersection(*static_cast<const AMesh_t *>(l), *static_cast<const AMesh_t *>(r), kFALSE);
2733  }
2734 
2735  /////////////////////////////////////////////////////////////////////////////
2736 
2738  {
2739  return build_difference(*static_cast<const AMesh_t *>(l), *static_cast<const AMesh_t *>(r), kFALSE);
2740  }
2741 
2742 }
TMatrix3x3 mmult_transpose_right(const TMatrix3x3 &m1, const TMatrix3x3 &m2)
Definition: CsgOps.cxx:1201
void build_split_group(const TMesh &meshA, const TMesh &meshB, const TBBoxTree &treeA, const TBBoxTree &treeB, OverlapTable_t &aOverlapsB, OverlapTable_t &bOverlapsA)
Definition: CsgOps.cxx:2391
TMeshWrapper< AMesh_t > AMeshWrapper_t
Definition: CsgOps.cxx:2385
int Invert(LASymMatrix &)
Definition: LaInverse.cxx:22
AMesh_t * build_difference(const AMesh_t &meshA, const AMesh_t &meshB, Bool_t preserve)
Definition: CsgOps.cxx:2629
Bool_t fuzzy_zero(const TVector3 &v)
Definition: CsgOps.cxx:795
Float_t Normalize(Float_t v[3])
Normalize a vector v in place.
Definition: TMath.cxx:499
void build_tree(const AMesh_t &mesh, TBBoxTree &tree)
Definition: CsgOps.cxx:2505
Double_t Z() const
Definition: TVector3.h:222
const Double_t * v1
Definition: TArcBall.cxx:33
void testPlane()
Definition: testTMath.cxx:129
TBBoxLeaf * LeafPtr_t
Definition: CsgOps.cxx:1509
TVector2 operator/(const TVector2 &v, Double_t s)
Definition: CsgOps.cxx:244
TVector3 cross(const TVector3 &v1, const TVector3 &v2)
Definition: CsgOps.cxx:816
AMesh_t * build_intersection(const AMesh_t &meshA, const AMesh_t &meshB, Bool_t preserve)
Definition: CsgOps.cxx:2575
TCanvas * c1
Definition: legend1.C:2
const char * Size
Definition: TXMLSetup.cxx:56
TBBoxNode * NodePtr_t
Definition: CsgOps.cxx:1510
Bool_t intersect_2d_no_bounds_check(const TLine3 &l1, const TLine3 &l2, Int_t majAxis, Double_t &l1Param, Double_t &l2Param)
Definition: CsgOps.cxx:1705
Double_t distance(const TPoint2 &p1, const TPoint2 &p2)
Definition: CsgOps.cxx:467
void swap(ROOT::THist< DIMENSIONS, PRECISION > &a, ROOT::THist< DIMENSIONS, PRECISION > &b) noexcept
Swap two histograms.
Definition: THist.h:196
Bool_t instersect_poly_with_line_3d(const TLine3 &l, const TGBinder &p1, const TPlane3 &plane, Double_t &a)
Definition: CsgOps.cxx:1766
Expr< TransposeOp< SMatrix< T, D1, D2, R >, T, D1, D2 >, T, D2, D1, typename TranspPolicy< T, D1, D2, R >::RepType > Transpose(const SMatrix< T, D1, D2, R > &rhs)
Matrix Transpose B(i,j) = A(j,i) returning a matrix expression.
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:170
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TArc * a
Definition: textangle.C:12
const Bool_t kFALSE
Definition: Rtypes.h:92
TPoint3 polygon_mid_point(const TGBinder &p1)
Definition: CsgOps.cxx:1801
Double_t triple(const TVector3 &v1, const TVector3 &v2, const TVector3 &v3)
Definition: CsgOps.cxx:823
Double_t Dot(const TVector3 &) const
Definition: TVector3.h:290
static const float upper
Definition: main.cpp:49
TVector2 operator+(const TVector2 &v1, const TVector2 &v2)
Definition: CsgOps.cxx:204
TLatex * t1
Definition: textangle.C:20
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
Bool_t point_in_polygon_test_3d(const TGBinder &p1, const TPlane3 &plane, const TPoint3 &origin, const TPoint3 &pointOnPlane)
Definition: CsgOps.cxx:1782
Double_t Angle(const TVector3 &v1, const TVector3 &v2)
Definition: CsgOps.cxx:809
Double_t Angle(const TVector3 &) const
return the angle w.r.t. another 3-vector
Definition: TVector3.cxx:239
Bool_t intersect_polygons(const TGBinderA &p1, const TGBinderB &p2, const TPlane3 &plane1, const TPlane3 &plane2)
Definition: CsgOps.cxx:1872
ABObj< sym, MatrixInverse< sym, ABObj< sym, LASymMatrix, double >, double >, double > Inverse(const ABObj< sym, LASymMatrix, double > &obj)
LAPACK Algebra functions specialize the Invert function for LASymMatrix.
Definition: LaInverse.h:27
std::vector< Int_t > PIndexList_t
Definition: CsgOps.cxx:2101
Double_t dot(const TVector2 &v1, const TVector2 &v2)
Definition: CsgOps.cxx:333
Double_t x[n]
Definition: legend1.C:17
Bool_t operator==(const Tuple2 &t1, const Tuple2 &t2)
Definition: CsgOps.cxx:137
TVector2 operator-(const TVector2 &v1, const TVector2 &v2)
Definition: CsgOps.cxx:212
TBaseMesh * BuildDifference(const TBaseMesh *leftOperand, const TBaseMesh *rightOperand)
Definition: CsgOps.cxx:2737
TVector3 & operator*=(Double_t)
Definition: TVector3.h:283
AMesh_t * build_union(const AMesh_t &meshA, const AMesh_t &meshB, Bool_t preserve)
Definition: CsgOps.cxx:2602
static double p2(double t, double a, double b, double c)
Double_t * fPnts
Definition: TBuffer3D.h:114
void partition_mesh(CMesh &mesh, const TMesh &mesh2, const OverlapTable_t &table)
Definition: CsgOps.cxx:2403
TBaseMesh * BuildIntersection(const TBaseMesh *leftOperand, const TBaseMesh *rightOperand)
Definition: CsgOps.cxx:2730
TVector2 & operator/=(Double_t s)
Definition: TVector2.h:121
void copy_mesh(const MeshA &source, MeshB &output)
Definition: CsgOps.cxx:2488
ROOT::Math::KDTree< _DataPoint > * BuildTree(const std::vector< const _DataPoint * > &vDataPoints, const unsigned int iBucketSize)
TMatrix3x3 mmult_transpose_left(const TMatrix3x3 &m1, const TMatrix3x3 &m2)
Definition: CsgOps.cxx:1186
Bool_t intersect_2d_bounds_check(const TLine3 &l1, const TLine3 &l2, Int_t majAxis, Double_t &l1Param, Double_t &l2Param)
Definition: CsgOps.cxx:1722
SEXP wrap(const TString &s)
Definition: RExports.h:81
UInt_t NbPols() const
Definition: TBuffer3D.h:84
REAL * vertex
Definition: triangle.c:512
Int_t * fPols
Definition: TBuffer3D.h:116
std::vector< PIndexList_t > OverlapTable_t
Definition: CsgOps.cxx:2103
Double_t operator*=(TVector2 const &v)
Definition: TVector2.h:118
TPoint2 lerp(const TPoint2 &p1, const TPoint2 &p2, Double_t t)
Definition: CsgOps.cxx:483
TPlane3 compute_plane(const TGBinder &poly)
Definition: CsgOps.cxx:1837
void Classification()
ROOT::R::TRInterface & r
Definition: Object.C:4
Double_t length(const TVector2 &v)
Definition: CsgOps.cxx:347
XPoint xy[kMAXMK]
Definition: TGX11.cxx:122
Double_t distance(const TPoint3 &p1, const TPoint3 &p2)
Definition: CsgOps.cxx:919
unsigned int r1[N_CITIES]
Definition: simanTSP.cxx:321
TVector3 & operator+=(const TVector3 &)
Definition: TVector3.h:265
TConnectedMeshWrapper< AConnectedMesh_t > AConnectedMeshWrapper_t
Definition: CsgOps.cxx:2386
#define Scalar
Definition: global.h:83
unsigned int UInt_t
Definition: RtypesCore.h:42
TMarker * m
Definition: textangle.C:8
Double_t X() const
Definition: TVector3.h:220
Generic 3D primitive description class.
Definition: TBuffer3D.h:19
TLine * l
Definition: textangle.C:4
TVector3 & operator-=(const TVector3 &)
Definition: TVector3.h:272
void classify_mesh(const TMesh &meshA, const TBBoxTree &aTree, CMesh &meshB)
Definition: CsgOps.cxx:2440
Double_t ACos(Double_t)
Definition: TMath.h:445
static double p1(double t, double a, double b)
const Double_t infinity
Definition: CsgOps.cxx:85
Double_t Cos(Double_t)
Definition: TMath.h:424
Bool_t intersect(const TPlane3 &p1, const TPlane3 &p2, TLine3 &output)
Definition: CsgOps.cxx:1691
#define Absolute(a)
Definition: triangle.c:4729
Double_t Dot(const TGLVector3 &v1, const TGLVector3 &v2)
Definition: TGLUtil.h:322
Bool_t intersect_poly_with_line_2d(const TLine3 &l, const TGBinder &p1, const TPlane3 &plane, Double_t &a, Double_t &b)
Definition: CsgOps.cxx:1741
return c2
Definition: legend2.C:14
I've modified some very nice bounding box tree code from Gino van der Bergen's Free Solid Library bel...
Definition: CsgOps.h:13
void Random()
Definition: utils.cpp:154
double Double_t
Definition: RtypesCore.h:55
TVector2 & operator-=(TVector2 const &v)
Definition: TVector2.h:114
TMesh< TestPolygon_t, TVertexBase > AMesh_t
Definition: CsgOps.cxx:2383
Double_t y[n]
Definition: legend1.C:17
PIndexList_t::iterator PIndexIt_t
Definition: CsgOps.cxx:2102
TBBoxInternal * InternalPtr_t
Definition: CsgOps.cxx:1520
Double_t distance2(const TPoint2 &p1, const TPoint2 &p2)
Definition: CsgOps.cxx:475
TMesh< TestPolygon_t, TCVertex > AConnectedMesh_t
Definition: CsgOps.cxx:2384
int currPoly
Definition: X3DBuffer.c:17
Binding & operator=(OUT(*fun)(void))
TBaseMesh * ConvertToMesh(const TBuffer3D &buff)
Definition: CsgOps.cxx:2656
void extract_classification(const AMesh_t &meshA, const AMesh_t &meshB, const TBBoxTree &aTree, const TBBoxTree &bTree, const OverlapTable_t &aOverlapsB, const OverlapTable_t &bOverlapsA, Int_t aClassification, Int_t bClassification, Bool_t reverseA, Bool_t reverseB, AMesh_t &output)
Definition: CsgOps.cxx:2549
Int_t * fSegs
Definition: TBuffer3D.h:115
TLine3 polygon_mid_point_ray(const TGBinder &p1, const TPlane3 &plane)
Definition: CsgOps.cxx:1829
UInt_t NbPnts() const
Definition: TBuffer3D.h:82
void extract_classification_preserve(const AMesh_t &meshA, const AMesh_t &meshB, const TBBoxTree &aTree, const TBBoxTree &bTree, const OverlapTable_t &aOverlapsB, const OverlapTable_t &bOverlapsA, Int_t aClassification, Int_t bClassification, Bool_t reverseA, Bool_t reverseB, AMesh_t &output)
Definition: CsgOps.cxx:2519
TVector3 Cross(const TVector3 &) const
Definition: TVector3.h:294
const Double_t epsilon2
Definition: CsgOps.cxx:84
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:202
Double_t Sin(Double_t)
Definition: TMath.h:421
Bool_t fuzzy_zero2(Double_t x)
Definition: CsgOps.cxx:103
AxisAngle::Scalar Distance(const AxisAngle &r1, const R &r2)
Distance between two rotations.
Definition: AxisAngle.h:326
TVector2 & operator+=(TVector2 const &v)
Definition: TVector2.h:113
Double_t Y() const
Definition: TVector3.h:221
ClassImp(TSlaveInfo) Int_t TSlaveInfo const TSlaveInfo * si
Used to sort slaveinfos by ordinal.
Definition: TProof.cxx:167
Bool_t fuzzy_zero(Double_t x)
Definition: CsgOps.cxx:96
TGLVector3 Cross(const TGLVector3 &v1, const TGLVector3 &v2)
Definition: TGLUtil.h:328
Int_t sign(Double_t x)
Definition: CsgOps.cxx:89
const Double_t epsilon
Definition: CsgOps.cxx:83
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
Bool_t fuzzy_equal(const TVector2 &v1, const TVector2 &v2)
Definition: CsgOps.cxx:361
std::complex< float_v > Z
Definition: main.cpp:120
TPolygonBase< TBlenderVProp, NullType_t > TestPolygon_t
Definition: CsgOps.cxx:2382
static void output(int code)
Definition: gifencode.c:226
Int_t which_side(const TGBinder &p1, const TPlane3 &plane1)
Definition: CsgOps.cxx:1813
const Bool_t kTRUE
Definition: Rtypes.h:91
TVector2 operator*(const TVector2 &v, Double_t s)
Definition: CsgOps.cxx:228
Double_t Angle(const TVector2 &v1, const TVector2 &v2)
Definition: CsgOps.cxx:368
const Int_t n
Definition: legend1.C:16
TBBox fit_bbox(const TGBinder &p1)
Definition: CsgOps.cxx:1861
Double_t length2(const TVector2 &v)
Definition: CsgOps.cxx:340
Int_t compute_classification(const Double_t &distance, const Double_t &epsil)
Definition: CsgOps.cxx:1732
unsigned int r2[N_CITIES]
Definition: simanTSP.cxx:322
TBaseMesh * BuildUnion(const TBaseMesh *leftOperand, const TBaseMesh *rightOperand)
Definition: CsgOps.cxx:2723
static const float lower
Definition: main.cpp:48
void Vertex(const Double_t *v)
segment * segs
Definition: X3DBuffer.c:21