91 return x < 0. ? -1 : x > 0. ? 1 : 0;
113 Tuple2(){SetValue(0, 0);}
114 Tuple2(
const Double_t *vv){SetValue(vv);}
125 const Double_t &U()
const{
return fCo[0];}
127 const Double_t &V()
const{
return fCo[1];}
130 const Double_t *GetValue()
const{
return fCo;}
131 void GetValue(
Double_t *vv)
const{vv[0] = fCo[0]; vv[1] = fCo[1];}
133 void SetValue(
const Double_t *vv){fCo[0] = vv[0]; fCo[1] = vv[1];}
139 return t1[0] == t2[0] && t1[1] == t2[1];
172 fCo[0] += vv[0]; fCo[1] += vv[1];
181 fCo[0] -= vv[0]; fCo[1] -= vv[1];
190 fCo[0] *= s; fCo[1] *= s;
return *
this;
198 return *
this *= 1. / s;
206 return TVector2(v1[0] + v2[0], v1[1] + v2[1]);
214 return TVector2(v1[0] - v2[0], v1[1] - v2[1]);
222 return TVector2(-v[0], -v[1]);
230 return TVector2(v[0] * s, v[1] * s);
246 return v * (1.0 / s);
254 return fCo[0] * vv[0] + fCo[1] * vv[1];
284 Bool_t TVector2::FuzzyZero()
const
300 TVector2 TVector2::Normalized()
const
302 return *
this / Length();
310 fCo[0] *= xx; fCo[1] *= yy;
318 return TVector2(fCo[0] * xx, fCo[1] * yy);
356 return v.FuzzyZero();
379 TPoint2 &operator += (
const TVector2 &v);
380 TPoint2 &operator -= (
const TVector2 &v);
384 Double_t Distance2(
const TPoint2 &p)
const;
385 TPoint2 Lerp(
const TPoint2 &p,
Double_t t)
const;
391 TPoint2 &TPoint2::operator += (
const TVector2 &v)
393 fCo[0] += v[0]; fCo[1] += v[1];
400 TPoint2 &TPoint2::operator-=(
const TVector2& v)
402 fCo[0] -= v[0]; fCo[1] -= v[1];
411 fCo[0] = v[0]; fCo[1] = v[1];
420 return (p - *
this).Length();
426 Double_t TPoint2::Distance2(
const TPoint2& p)
const
428 return (p - *
this).Length2();
434 TPoint2 TPoint2::Lerp(
const TPoint2 &p,
Double_t t)
const
436 return TPoint2(fCo[0] + (p[0] - fCo[0]) * t,
437 fCo[1] + (p[1] - fCo[1]) * t);
445 return TPoint2(p[0] + v[0], p[1] + v[1]);
453 return TPoint2(p[0] - v[0], p[1] - v[1]);
461 return TVector2(p1[0] - p2[0], p1[1] - p2[1]);
469 return p1.Distance(p2);
477 return p1.Distance2(p2);
485 return p1.Lerp(p2, t);
492 Tuple3(){SetValue(0, 0, 0);}
493 Tuple3(
const Double_t *v){SetValue(v);}
507 const Double_t *GetValue()
const{
return fCo; }
520 fCo[0] = xx; fCo[1] = yy; fCo[2] = zz;
528 return t1[0] == t2[0] && t1[1] == t2[1] && t1[2] == t2[2];
556 Int_t ClosestAxis()
const;
566 fCo[0] += v[0]; fCo[1] += v[1]; fCo[2] += v[2];
575 fCo[0] -= v[0]; fCo[1] -= v[1]; fCo[2] -= v[2];
584 fCo[0] *= s; fCo[1] *= s; fCo[2] *= s;
601 return fCo[0] * v[0] + fCo[1] * v[1] + fCo[2] * v[2];
631 Bool_t TVector3::FuzzyZero()
const
639 void TVector3::NoiseGate(
Double_t threshold)
641 if (Length2() < threshold) SetValue(0., 0., 0.);
656 return TVector3(v[0] * s, v[1] * s, v[2] * s);
669 TVector3 TVector3::Normalized()
const
671 return *
this / Length();
677 TVector3 TVector3::SafeNormalized()
const
688 fCo[0] *= xx; fCo[1] *= yy; fCo[2] *= zz;
696 return TVector3(fCo[0] * xx, fCo[1] * yy, fCo[2] * zz);
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]);
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]);
731 Int_t TVector3::ClosestAxis()
const
734 return a[0] < a[1] ? (a[1] < a[2] ? 2 : 1) : (a[0] < a[2] ? 2 : 0);
741 return TVector3(v1[0] + v2[0], v1[1] + v2[1], v1[2] + v2[2]);
748 return TVector3(v1[0] - v2[0], v1[1] - v2[1], v1[2] - v2[2]);
755 return TVector3(-v[0], -v[1], -v[2]);
769 return TVector3(v1[0] * v2[0], v1[1] * v2[1], v1[2] * v2[2]);
797 return v.FuzzyZero();
816 TVector3
cross(
const TVector3 &v1,
const TVector3 &v2)
825 return v1.Triple(v2, v3);
834 TPoint3 &operator += (
const TVector3 &v);
835 TPoint3 &operator -= (
const TVector3 &v);
839 Double_t Distance2(
const TPoint3 &p)
const;
840 TPoint3 Lerp(
const TPoint3 &p,
Double_t t)
const;
846 TPoint3 &TPoint3::operator += (
const TVector3 &v)
848 fCo[0] += v[0]; fCo[1] += v[1]; fCo[2] += v[2];
855 TPoint3 &TPoint3::operator-=(
const TVector3 &v)
857 fCo[0] -= v[0]; fCo[1] -= v[1]; fCo[2] -= v[2];
866 fCo[0] = v[0]; fCo[1] = v[1]; fCo[2] = v[2];
875 return (p - *
this).Length();
881 Double_t TPoint3::Distance2(
const TPoint3 &p)
const
883 return (p - *
this).Length2();
889 TPoint3 TPoint3::Lerp(
const TPoint3 &p,
Double_t t)
const
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);
900 return TPoint3(p[0] + v[0], p[1] + v[1], p[2] + v[2]);
907 return TPoint3(p[0] - v[0], p[1] - v[1], p[2] - v[2]);
914 return TVector3(p1[0] - p2[0], p1[1] - p2[1], p1[2] - p2[2]);
921 return p1.Distance(p2);
928 return p1.Distance2(p2);
935 return p1.Lerp(p2, t);
943 Tuple4(){SetValue(0, 0, 0, 0);}
944 Tuple4(
const Double_t *v){SetValue(v);}
947 SetValue(xx, yy, zz, ww);
960 const Double_t &W()
const{
return fCo[3];}
963 const Double_t *GetValue()
const{
return fCo;}
967 v[0] = fCo[0]; v[1] = fCo[1]; v[2] = fCo[2]; v[3] = fCo[3];
972 fCo[0] = v[0]; fCo[1] = v[1]; fCo[2] = v[2]; fCo[3] = v[3];
976 fCo[0] = xx; fCo[1] = yy; fCo[2] = zz; fCo[3] = ww;
984 return t1[0] == t2[0] && t1[1] == t2[1] && t1[2] == t2[2] && t1[3] == t2[3];
993 TMatrix3x3(
const Double_t *m){SetValue(m);}
994 TMatrix3x3(
const TVector3 &euler){SetEuler(euler);}
997 SetEuler(euler); Scale(s[0], s[1], s[2]);
1003 SetValue(xx, xy, xz, yx, yy, yz, zx, zy, zz);
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;
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;
1023 void SetEuler(
const TVector3 &euler)
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);
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;
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);
1055 SetValue(1., 0., 0., 0., 1., 0., 0., 0., 1.);
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;
1064 TMatrix3x3 &operator *= (
const TMatrix3x3 &m);
1067 return fEl[0][
c] * v[0] + fEl[1][
c] * v[1] + fEl[2][
c] * v[2];
1075 TMatrix3x3 Adjoint()
const;
1077 TMatrix3x3 Transposed()
const;
1086 TMatrix3x3 &TMatrix3x3::operator *= (
const TMatrix3x3 &m)
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]));
1097 Double_t TMatrix3x3::Determinant()
const
1099 return triple((*
this)[0], (*
this)[1], (*
this)[2]);
1115 TMatrix3x3 TMatrix3x3::Transposed()
const
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]);
1127 *
this = Transposed();
1133 TMatrix3x3 TMatrix3x3::Adjoint()
const
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));
1145 TVector3 co(Cofac(1, 1, 2, 2), Cofac(1, 2, 2, 0), Cofac(1, 0, 2, 1));
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);
1165 return TVector3(
dot(m[0], v),
dot(m[1], v),
dot(m[2], v));
1172 return TVector3(m.Tdot(0, v), m.Tdot(1, v), m.Tdot(2, v));
1177 TMatrix3x3
operator * (
const TMatrix3x3 &m1,
const TMatrix3x3 &m2)
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]));
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]);
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]));
1217 TLine3(
const TPoint3 &
p1,
const TPoint3 &
p2);
1218 TLine3(
const TPoint3 &
p1,
const TVector3 &v);
1221 static TLine3 InfiniteRay(
const TPoint3 &
p1,
const TVector3 &v);
1222 const TVector3 &Direction()
const {
return fDir;}
1223 const TPoint3 &Origin()
const{
return fOrigin;}
1227 return (i == 0 ? fBounds[0] : fBounds[1]);
1231 return (i == 0 ? fBounds[0] : fBounds[1]);
1235 return (i == 0 ? fParams[0] : fParams[1]);
1239 return (i == 0 ? fParams[0] : fParams[1]);
1241 TVector3 UnboundSmallestVector(
const TPoint3 &point)
const
1244 return diff - fDir * diff.
Dot(fDir);
1246 Double_t UnboundClosestParameter(
const TPoint3 &point)
const
1249 return diff.Dot(fDir);
1251 Double_t UnboundDistance(
const TPoint3& point)
const
1253 return UnboundSmallestVector(point).Length();
1257 return ((fParams[0] - epsilon < t) || (!fBounds[0])) && ((fParams[1] > t +
epsilon) || (!fBounds[1]));
1264 TLine3::TLine3() : fOrigin(0,0,0), fDir(1,0,0)
1267 fParams[0] = 0; fParams[1] = 1;
1273 TLine3::TLine3(
const TPoint3 &
p1,
const TPoint3 &
p2) : fOrigin(p1), fDir(p2-p1)
1276 fParams[0] = 0; fParams[1] = 1;
1282 TLine3::TLine3(
const TPoint3 &
p1,
const TVector3 &v): fOrigin(p1), fDir(v)
1285 fParams[0] = 0; fParams[1] = 1;
1292 : fOrigin(p1), fDir(v)
1294 fBounds[0] = bound1; fBounds[1] = bound2;
1295 fParams[0] = 0; fParams[1] = 1;
1298 class TPlane3 :
public Tuple4 {
1303 TPlane3(
const TPlane3 & p):Tuple4(p){}
1321 n = n.SafeNormalized();
1323 fCo[0] = n.
X(); fCo[1] = n.
Y(); fCo[2] = n.
Z(); fCo[3] = -d;
1333 fCo[0] = mn.
X(); fCo[1] = mn.
Y(); fCo[2] = mn.
Z(); fCo[3] = -md;
1339 TPlane3::TPlane3() : Tuple4()
1341 fCo[0] = 1.; fCo[1] = 0.;
1342 fCo[2] = 0.; fCo[3] = 0.;
1350 return TVector3(fCo[0],fCo[1],fCo[2]);
1366 fCo[0] = -fCo[0]; fCo[1] = -fCo[1]; fCo[2] = -fCo[2]; fCo[3] = -fCo[3];
1374 fCo[0] = rhs.fCo[0]; fCo[1] = rhs.fCo[1]; fCo[2] = rhs.fCo[2]; fCo[3] = rhs.fCo[3];
1383 return Normal().Dot(v) + fCo[3];
1395 TBBox(
const TPoint3 &mini,
const TPoint3 &maxi){SetValue(mini,maxi);}
1397 const TPoint3 &Center()
const{
return fCenter;}
1398 const TVector3 &Extent()
const{
return fExtent;}
1399 TPoint3 &Center(){
return fCenter;}
1400 TVector3 &Extent(){
return fExtent;}
1402 void SetValue(
const TPoint3 &mini,
const TPoint3 &maxi)
1404 fExtent = (maxi - mini) / 2;
1405 fCenter = mini + fExtent;
1407 void Enclose(
const TBBox &a,
const TBBox &b)
1424 fCenter.SetValue(0., 0., 0.);
1425 fExtent.SetValue(-infinity, -infinity, -infinity);
1427 void Include(
const TPoint3 &p)
1442 void Include(
const TBBox &b)
1448 return fCenter[i] - fExtent[i];
1452 return fCenter[i] + fExtent[i];
1454 TPoint3 Lower()
const
1456 return fCenter - fExtent;
1458 TPoint3 Upper()
const
1460 return fCenter + fExtent;
1466 Int_t LongestAxis()
const
1468 return fExtent.ClosestAxis();
1470 Bool_t IntersectXRay(
const TPoint3 &xBase)
const
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;}
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];
1492 enum ETagType {kLeaf, kInternal};
1497 class TBBoxLeaf :
public TBBoxNode {
1501 TBBoxLeaf() : fPolyIndex(0) {}
1502 TBBoxLeaf(
Int_t polyIndex,
const TBBox &bbox) : fPolyIndex(polyIndex)
1512 class TBBoxInternal :
public TBBoxNode {
1515 NodePtr_t fRightSon;
1516 TBBoxInternal() : fLeftSon(0) ,fRightSon(0) {}
1517 TBBoxInternal(
Int_t n, LeafPtr_t leafIt);
1526 InternalPtr_t fInternals;
1530 TBBoxTree() : fBranch(0), fLeaves(0), fInternals(0), fNumLeaves(0) {}
1531 NodePtr_t RootNode()
const{
return fInternals;}
1535 delete[] fInternals;
1540 void RecursiveTreeBuild(
Int_t n, LeafPtr_t leafIt);
1546 TBBoxInternal::TBBoxInternal(
Int_t n, LeafPtr_t leafIt) :
1547 fLeftSon(0) ,fRightSon(0)
1552 fBBox.Include(leafIt[i].fBBox);
1562 fNumLeaves = numLeaves;
1563 fInternals =
new TBBoxInternal[numLeaves];
1564 RecursiveTreeBuild(fNumLeaves,fLeaves);
1570 void TBBoxTree::RecursiveTreeBuild(
Int_t n, LeafPtr_t leafIt)
1572 fInternals[fBranch] = TBBoxInternal(n,leafIt);
1573 TBBoxInternal &aBBox = fInternals[fBranch];
1576 Int_t axis = aBBox.fBBox.LongestAxis();
1580 if (leafIt[i].fBBox.Center()[axis] < aBBox.fBBox.Center()[axis]) {
1588 if (mid == 0 || mid == n) {
1592 aBBox.fRightSon = fInternals + fBranch;
1593 RecursiveTreeBuild(mid,leafIt);
1595 aBBox.fRightSon = leafIt;
1598 aBBox.fLeftSon = fInternals + fBranch;
1599 RecursiveTreeBuild(n - mid, leafIt + mid);
1601 aBBox.fLeftSon = leafIt + mid;
1605 class TBlenderVProp {
1610 TBlenderVProp(
Int_t vIndex) : fVertexIndex(vIndex){}
1611 TBlenderVProp(
Int_t vIndex,
const TBlenderVProp &,
1612 const TBlenderVProp &,
const Double_t &)
1614 fVertexIndex = vIndex;
1616 TBlenderVProp() : fVertexIndex(-1){}
1617 operator Int_t()
const
1619 return fVertexIndex;
1623 fVertexIndex = i;
return *
this;
1627 template <
class TMesh>
1628 class TPolygonGeometry {
1630 typedef typename TMesh::Polygon TPolygon;
1634 const TPolygon &fPoly;
1636 TPolygonGeometry(
const TMesh &mesh,
Int_t pIndex)
1637 : fMesh(mesh), fPoly(mesh.Polys()[pIndex])
1639 TPolygonGeometry(
const TMesh &mesh,
const TPolygon &poly)
1640 : fMesh(mesh), fPoly(poly)
1642 const TPoint3 &operator [] (
Int_t i)
const
1644 return fMesh.Verts()[fPoly[i]].Pos();
1648 return fPoly.Size();
1652 TPolygonGeometry(
const TPolygonGeometry &);
1653 TPolygonGeometry&
operator = (TPolygonGeometry &);
1656 template <
class TPolygon,
class TVertex>
1657 class TMesh :
public TBaseMesh {
1659 typedef std::vector<TVertex> VLIST;
1660 typedef std::vector<TPolygon> PLIST;
1661 typedef TPolygon Polygon;
1663 typedef TPolygonGeometry<TMesh> TGBinder;
1670 VLIST &Verts(){
return fVerts;}
1671 const VLIST &Verts()
const{
return fVerts;}
1672 PLIST &Polys(){
return fPolys;}
1673 const PLIST &Polys()
const{
return fPolys;}
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();}
1683 return fPolys[polyNum][vertexNum];
1687 const Int_t cofacTable[3][2] = {{1,2}, {0,2}, {0,1}};
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]);
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];
1715 l1Param = (l2.Direction()[ind2] * zX - l2.Direction()[ind1] * zY)/det;
1716 l2Param = -(-l1.Direction()[ind2] * zX + l1.Direction()[ind1] * zY)/det;
1726 if (!isect)
return kFALSE;
1727 return l1.IsParameterOnLine(l1Param) && l2.IsParameterOnLine(l2Param);
1735 else return distance < 0 ? 1 : 2;
1740 template<
typename TGBinder>
1744 Int_t majAxis = plane.Normal().ClosestAxis();
1745 Int_t lastInd = p1.Size()-1;
1747 Double_t isectParam(0.), isectParam2(0.);
1750 Int_t isectsFound(0);
1751 for (i=0;i<=lastInd; j=i,i++ ) {
1752 TLine3 testLine(p1[j],p1[i]);
1760 return (isectsFound > 0);
1765 template<
typename TGBinder>
1769 Double_t determinant = l.Direction().Dot(plane.Normal());
1771 a = -plane.Scalar() - plane.Normal().Dot(l.Origin());
1773 if (a <= 0 )
return kFALSE;
1774 if (!l.IsParameterOnLine(a))
return kFALSE;
1775 TPoint3 pointOnPlane = l.Origin() + l.Direction() *
a;
1781 template<
typename TGBinder>
1783 const TPoint3 &pointOnPlane)
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;
1800 template <
typename TGBinder>
1803 TPoint3 midPoint(0., 0., 0.);
1805 for (i=0; i < p1.Size(); i++)
1807 return TPoint3(midPoint[0] / i, midPoint[1] / i, midPoint[2] / i);
1812 template <
typename TGBinder>
1817 for (i=0; i<p1.Size(); i++) {
1818 Double_t signedDistance = plane1.SignedDistance(p1[i]);
1820 signedDistance < 0 ? (output |= 1) : (output |=2);
1828 template <
typename TGBinder>
1836 template <
typename TGBinder>
1839 TPoint3 plast(poly[poly.Size()-1]);
1843 for (j = 0; j < poly.Size(); j++) {
1845 edge = pivot - plast;
1846 if (!edge.FuzzyZero())
break;
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);
1860 template <
typename TGBinder>
1863 TBBox bbox; bbox.SetEmpty();
1864 for (
Int_t i = 0; i < p1.Size(); ++i)
1865 bbox.Include(p1[i]);
1871 template<
typename TGBinderA,
typename TGBinderB>
1873 const TPlane3 &plane1,
const TPlane3 &plane2)
1875 TLine3 intersectLine;
1876 if (!
intersect(plane1, plane2, intersectLine))
1888 return (maxOMin <= minOMax);
1891 template <
class TMesh,
class TSplitFunctionBinder>
1892 class TSplitFunction {
1895 TSplitFunctionBinder &fFunctionBinder;
1898 TSplitFunction(TMesh &mesh, TSplitFunctionBinder &functionBindor)
1899 : fMesh(mesh), fFunctionBinder(functionBindor)
1901 void SplitPolygon(
const Int_t p1Index,
const TPlane3 &plane,
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();
1913 Int_t totalClassification(lastClassification);
1916 for (i = 0; i < p.Size(); j = i, ++i)
1918 Int_t newIndex = p[i];
1919 TPoint3 aVertex = fMesh.Verts()[newIndex].Pos();
1921 if ((newClassification != lastClassification) && newClassification && lastClassification)
1923 Int_t newVertexIndex = fMesh.Verts().size();
1925 fMesh.Verts().push_back(VERTEX_t());
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);
1935 Classify(inP.Verts(),outP.Verts(),newClassification, p.VertexProps(i));
1936 lastClassification = newClassification;
1937 totalClassification |= newClassification;
1938 lastVertex = aVertex;
1939 lastIndex = newIndex;
1941 if (totalClassification == 3) {
1943 outPiece = fMesh.Polys().size();
1944 fMesh.Polys()[p1Index] = inP;
1945 fMesh.Polys().push_back(outP);
1946 fFunctionBinder.ConnectPolygon(inPiece);
1947 fFunctionBinder.ConnectPolygon(outPiece);
1949 fFunctionBinder.ConnectPolygon(p1Index);
1950 if (totalClassification == 1) {
1960 void Classify(
typename TMesh::Polygon::TVPropList &inGroup,
1961 typename TMesh::Polygon::TVPropList &outGroup,
1962 Int_t classification,
1963 typename TMesh::Polygon::TVProp prop)
1965 switch (classification) {
1967 inGroup.push_back(prop);
1968 outGroup.push_back(prop);
1971 inGroup.push_back(prop);
1974 outGroup.push_back(prop);
1983 template <
typename PROP>
1984 class TDefaultSplitFunctionBinder {
1986 void DisconnectPolygon(
Int_t){}
1987 void ConnectPolygon(
Int_t){}
1988 void InsertVertexAlongEdge(
Int_t,
Int_t,
const PROP &){}
1991 template <
typename TMesh>
1992 class TMeshWrapper {
1997 typedef typename TMesh::Polygon Polygon;
1999 typedef typename TMesh::VLIST VLIST;
2000 typedef typename TMesh::PLIST PLIST;
2001 typedef TPolygonGeometry<TMeshWrapper> TGBinder;
2002 typedef TMeshWrapper<TMesh> MyType;
2005 TMeshWrapper(TMesh &mesh):fMesh(mesh){}
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();}
2012 void ComputePlanes();
2013 TBBox ComputeBBox()
const;
2014 void SplitPolygon(
Int_t p1Index,
const TPlane3 &plane,
2021 template <
typename TMesh>
2022 void TMeshWrapper<TMesh>::ComputePlanes()
2024 PLIST& polyList = Polys();
2026 for (i=0;i < polyList.size(); i++) {
2027 TGBinder binder(*
this, i);
2035 template <
typename TMesh>
2036 TBBox TMeshWrapper<TMesh>::ComputeBBox()
const
2038 const VLIST &vertexList = Verts();
2042 for (i=0;i<vertexList.size(); i++)
2043 bbox.Include(vertexList[i].Pos());
2049 template<
typename TMesh>
2050 void TMeshWrapper<TMesh>::SplitPolygon(
Int_t p1Index,
const TPlane3 &plane,
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);
2061 template <
typename AVProp,
typename AFProp>
2062 class TPolygonBase {
2064 typedef AVProp TVProp;
2065 typedef AFProp TFProp;
2066 typedef std::vector<TVProp> TVPropList;
2067 typedef typename TVPropList::iterator TVPropIt;
2073 Int_t fClassification;
2076 const TVPropList &Verts()
const{
return fVerts;}
2077 TVPropList &Verts(){
return fVerts;}
2080 Int_t operator[](
Int_t i)
const {
return fVerts[i];}
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();}
2092 std::reverse(fVerts.begin(),fVerts.end());
2096 TFProp &FProp(){
return fFaceProp;}
2097 const TFProp &FProp()
const{
return fFaceProp;}
2098 void AddProp(
const TVProp &prop){fVerts.push_back(prop);}
2105 template <
typename TMesh>
2106 class TreeIntersector {
2108 OverlapTable_t *fAoverlapsB;
2109 OverlapTable_t *fBoverlapsA;
2110 const TMesh *fMeshA;
2111 const TMesh *fMeshB;
2114 TreeIntersector(
const TBBoxTree &a,
const TBBoxTree &b,
2115 OverlapTable_t *aOverlapsB, OverlapTable_t *bOverlapsA,
2116 const TMesh *meshA,
const TMesh *meshB)
2118 fAoverlapsB = aOverlapsB;
2119 fBoverlapsA = bOverlapsA;
2122 MarkIntersectingPolygons(a.RootNode(),b.RootNode());
2126 void MarkIntersectingPolygons(
const TBBoxNode *a,
const TBBoxNode *b)
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;
2133 TPolygonGeometry<TMesh> pg1(*fMeshA,la->fPolyIndex);
2134 TPolygonGeometry<TMesh> pg2(*fMeshB,lb->fPolyIndex);
2137 fMeshB->Polys()[lb->fPolyIndex].Plane())) {
2138 (*fAoverlapsB)[lb->fPolyIndex].push_back(la->fPolyIndex);
2139 (*fBoverlapsA)[la->fPolyIndex].push_back(lb->fPolyIndex);
2141 }
else if ( a->fTag == TBBoxNode::kLeaf || (b->fTag != TBBoxNode::kLeaf && a->fBBox.Size() < b->fBBox.Size()))
2143 MarkIntersectingPolygons(a, ((
const TBBoxInternal *)b)->fLeftSon);
2144 MarkIntersectingPolygons(a, ((
const TBBoxInternal *)b)->fRightSon);
2146 MarkIntersectingPolygons(((
const TBBoxInternal *)a)->fLeftSon, b);
2147 MarkIntersectingPolygons(((
const TBBoxInternal *)a)->fRightSon, b);
2152 template<
typename TMesh>
2153 class TRayTreeIntersector {
2155 const TMesh *fMeshA;
2160 TRayTreeIntersector(
const TBBoxTree &a,
const TMesh *meshA,
const TLine3 &xRay,
Int_t &polyIndex)
2161 : fMeshA(meshA), fLastIntersectValue(infinity), fPolyIndex(-1)
2163 FindIntersectingPolygons(a.RootNode(),xRay);
2164 polyIndex = fPolyIndex;
2168 void FindIntersectingPolygons(
const TBBoxNode*a,
const TLine3& xRay)
2170 if ((xRay.Origin().X() + fLastIntersectValue < a->fBBox.Lower(0)) ||!a->fBBox.IntersectXRay(xRay.Origin()))
2172 if (a->fTag == TBBoxNode::kLeaf) {
2173 const TBBoxLeaf *la = (
const TBBoxLeaf *)a;
2175 TPolygonGeometry<TMesh> pg(*fMeshA, la->fPolyIndex);
2178 if (testParameter < fLastIntersectValue) {
2179 fLastIntersectValue = testParameter;
2180 fPolyIndex = la->fPolyIndex;
2184 FindIntersectingPolygons(((
const TBBoxInternal*)a)->fLeftSon, xRay);
2185 FindIntersectingPolygons(((
const TBBoxInternal*)a)->fRightSon, xRay);
2197 TVertexBase():fVertexMap(-1) {}
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();}
2208 class TCVertex :
public TVertexBase {
2210 PIndexList_t fPolygons;
2212 TCVertex() : TVertexBase(){}
2213 TCVertex(
const TVertexBase&
vertex) : TVertexBase(vertex){}
2215 TCVertex &
operator = (
const TVertexBase &other)
2220 const PIndexList_t &Polys()
const{
return fPolygons;}
2221 PIndexList_t &Polys(){
return fPolygons;}
2223 Int_t &operator [] (
Int_t i) {
return fPolygons[i];}
2224 const Int_t &operator [] (
Int_t i)
const{
return fPolygons[i];}
2226 void AddPoly(
Int_t polyIndex){fPolygons.push_back(polyIndex);}
2227 void RemovePolygon(
Int_t polyIndex)
2229 PIndexIt_t foundIt = std::find(fPolygons.begin(), fPolygons.end(), polyIndex);
2230 if (foundIt != fPolygons.end()) {
2232 fPolygons.pop_back();
2237 template <
typename TMesh>
2238 class TConnectedMeshWrapper {
2241 UInt_t fUniqueEdgeTestId;
2243 typedef typename TMesh::Polygon Polygon;
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;
2251 TConnectedMeshWrapper(TMesh &mesh) : fMesh(mesh), fUniqueEdgeTestId(0){}
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);
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,
2267 template <
class CMesh>
class TSplitFunctionBinder {
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)
2277 fMesh.InsertVertexAlongEdge(lastIndex, newIndex,prop);
2284 template <
typename TMesh>
2285 void TConnectedMeshWrapper<TMesh>::BuildVertexPolyLists()
2288 for (i=0; i < Polys().size(); i++)
2295 template <
typename TMesh>
2296 void TConnectedMeshWrapper<TMesh>::DisconnectPolygon(
Int_t polyIndex)
2298 const Polygon &poly = Polys()[polyIndex];
2300 for (j=0;j<poly.Verts().size(); j++) {
2301 Verts()[poly[j]].RemovePolygon(polyIndex);
2308 template <
typename TMesh>
2309 void TConnectedMeshWrapper<TMesh>::ConnectPolygon(
Int_t polyIndex)
2311 const Polygon &poly = Polys()[polyIndex];
2313 for (j=0;j<poly.Verts().size(); j++) {
2314 Verts()[poly[j]].AddPoly(polyIndex);
2321 template <
typename TMesh>
2322 void TConnectedMeshWrapper<TMesh>::EdgePolygons(
Int_t v1,
Int_t v2, PIndexList_t &polys)
2324 ++fUniqueEdgeTestId;
2325 Vertex &vb1 = Verts()[
v1];
2327 for (i=0;i < vb1.Polys().size(); ++i){Polys()[vb1[i]].Classification() = fUniqueEdgeTestId;}
2328 Vertex &vb2 = Verts()[v2];
2330 for (j=0;j < vb2.Polys().size(); ++j) {
2331 if ((
UInt_t)Polys()[vb2[j]].Classification() == fUniqueEdgeTestId) {
2332 polys.push_back(vb2[j]);
2340 template <
typename TMesh>
2341 void TConnectedMeshWrapper<TMesh>::InsertVertexAlongEdge(
Int_t v1,
Int_t v2,
const VProp &prop)
2343 PIndexList_t npolys;
2344 EdgePolygons(v1,v2,npolys);
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);
2360 Verts()[newVertex].AddPoly(npolys[i]);
2369 template <
typename TMesh>
2370 void TConnectedMeshWrapper<TMesh>::SplitPolygon(
Int_t p1Index,
const TPlane3 &plane,
2374 TSplitFunctionBinder<MyType> functionBindor(*
this);
2375 TSplitFunction<MyType,TSplitFunctionBinder<MyType> > splitFunction(*
this,functionBindor);
2376 splitFunction.SplitPolygon(p1Index, plane, inPiece, outPiece, onEpsilon);
2379 struct NullType_t{};
2390 template <
class TMesh>
2392 const TBBoxTree &treeA,
const TBBoxTree &treeB,
2393 OverlapTable_t &aOverlapsB, OverlapTable_t &bOverlapsA)
2397 TreeIntersector<TMesh>(treeA, treeB, &aOverlapsB, &bOverlapsA, &meshA, &meshB);
2402 template <
class CMesh,
class TMesh>
2407 for (i = 0; i < table.size(); i++) {
2408 if (table[i].size()) {
2409 PIndexList_t fragments;
2410 fragments.push_back(i);
2412 for (j =0 ; j <table[i].size(); ++j) {
2413 PIndexList_t newFragments;
2414 TPlane3 splitPlane = mesh2.Polys()[table[i][j]].Plane();
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();
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);
2428 newFragments.push_back(fragments[k]);
2431 fragments = newFragments;
2439 template <
typename CMesh,
typename TMesh>
2443 for (i = 0; i < meshB.Polys().size(); i++) {
2444 typename CMesh::TGBinder pg(meshB,i);
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;
2453 meshB.Polys()[i].Classification()= 2;
2456 meshB.Polys()[i].Classification()= 2;
2463 template <
typename CMesh,
typename TMesh>
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();
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;
2479 newPolygon.VertexProps(j) = meshA.Verts()[newPolygon[j]].VertexMap();
2487 template <
typename MeshA,
typename MeshB>
2490 Int_t vertexNum = source.Verts().size();
2491 Int_t polyNum = source.Polys().size();
2493 typedef typename MeshB::VLIST VLIST_t;
2494 typedef typename MeshB::PLIST PLIST_t;
2496 output.Verts() = VLIST_t(vertexNum);
2497 output.Polys() = PLIST_t(polyNum);
2499 std::copy(source.Verts().begin(), source.Verts().end(), output.Verts().begin());
2500 std::copy(source.Polys().begin(), source.Polys().end(), output.Polys().begin());
2507 Int_t numLeaves = mesh.Polys().size();
2508 TBBoxLeaf *aLeaves =
new TBBoxLeaf[numLeaves];
2510 for (i=0;i< mesh.Polys().size(); i++) {
2511 TPolygonGeometry<AMesh_t> pg(mesh,i);
2512 aLeaves[i] = TBBoxLeaf(i,
fit_bbox(pg));
2514 tree.BuildTree(aLeaves,numLeaves);
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,
2531 AConnectedMesh_t meshAPartitioned;
2532 AConnectedMesh_t meshBPartitioned;
2535 AConnectedMeshWrapper_t meshAWrapper(meshAPartitioned);
2536 AConnectedMeshWrapper_t meshBWrapper(meshBPartitioned);
2537 meshAWrapper.BuildVertexPolyLists();
2538 meshBWrapper.BuildVertexPolyLists();
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,
2561 AMesh_t meshAPartitioned(meshA);
2562 AMesh_t meshBPartitioned(meshB);
2563 AMeshWrapper_t meshAWrapper(meshAPartitioned);
2564 AMeshWrapper_t meshBWrapper(meshBPartitioned);
2577 TBBoxTree aTree, bTree;
2580 OverlapTable_t bOverlapsA(meshA.Polys().size());
2581 OverlapTable_t aOverlapsB(meshB.Polys().size());
2586 meshA, meshB, aTree, bTree,
2587 aOverlapsB, bOverlapsA,
2592 meshA, meshB, aTree, bTree,
2593 aOverlapsB, bOverlapsA,
2604 TBBoxTree aTree, bTree;
2607 OverlapTable_t bOverlapsA(meshA.Polys().size());
2608 OverlapTable_t aOverlapsB(meshB.Polys().size());
2613 meshA, meshB, aTree, bTree,
2614 aOverlapsB, bOverlapsA,
2619 meshA, meshB, aTree, bTree,
2620 aOverlapsB, bOverlapsA,
2631 TBBoxTree aTree, bTree;
2634 OverlapTable_t bOverlapsA(meshA.Polys().size());
2635 OverlapTable_t aOverlapsB(meshB.Polys().size());
2640 meshA, meshB, aTree, bTree,
2641 aOverlapsB, bOverlapsA,
2646 meshA, meshB, aTree, bTree,
2647 aOverlapsB, bOverlapsA,
2658 AMesh_t *newMesh =
new AMesh_t;
2661 newMesh->Verts().resize(buff.
NbPnts());
2664 newMesh->Verts()[i] = TVertexBase(v[i * 3], v[i * 3 + 1], v[i * 3 + 2]);
2669 newMesh->Polys().resize(buff.
NbPols());
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];
2677 Int_t s2 = pols[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};
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];
2690 numPnts[0] = segEnds[0], numPnts[1] = segEnds[1], numPnts[2] = segEnds[2];
2693 currPoly.AddProp(TBlenderVProp(numPnts[0]));
2694 currPoly.AddProp(TBlenderVProp(numPnts[1]));
2695 currPoly.AddProp(TBlenderVProp(numPnts[2]));
2697 Int_t lastAdded = numPnts[2];
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];
2707 currPoly.AddProp(TBlenderVProp(segEnds[0]));
2708 lastAdded = segEnds[0];
2711 j += segmentCol + 2;
2714 AMeshWrapper_t
wrap(*newMesh);
2716 wrap.ComputePlanes();
2725 return build_union(*static_cast<const AMesh_t *>(l), *static_cast<const AMesh_t *>(r),
kFALSE);
TMatrix3x3 mmult_transpose_right(const TMatrix3x3 &m1, const TMatrix3x3 &m2)
void build_split_group(const TMesh &meshA, const TMesh &meshB, const TBBoxTree &treeA, const TBBoxTree &treeB, OverlapTable_t &aOverlapsB, OverlapTable_t &bOverlapsA)
TMeshWrapper< AMesh_t > AMeshWrapper_t
int Invert(LASymMatrix &)
AMesh_t * build_difference(const AMesh_t &meshA, const AMesh_t &meshB, Bool_t preserve)
Bool_t fuzzy_zero(const TVector3 &v)
Float_t Normalize(Float_t v[3])
Normalize a vector v in place.
void build_tree(const AMesh_t &mesh, TBBoxTree &tree)
TVector2 operator/(const TVector2 &v, Double_t s)
TVector3 cross(const TVector3 &v1, const TVector3 &v2)
AMesh_t * build_intersection(const AMesh_t &meshA, const AMesh_t &meshB, Bool_t preserve)
Bool_t intersect_2d_no_bounds_check(const TLine3 &l1, const TLine3 &l2, Int_t majAxis, Double_t &l1Param, Double_t &l2Param)
Double_t distance(const TPoint2 &p1, const TPoint2 &p2)
void swap(ROOT::THist< DIMENSIONS, PRECISION > &a, ROOT::THist< DIMENSIONS, PRECISION > &b) noexcept
Swap two histograms.
Bool_t instersect_poly_with_line_3d(const TLine3 &l, const TGBinder &p1, const TPlane3 &plane, Double_t &a)
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)
TPoint3 polygon_mid_point(const TGBinder &p1)
Double_t triple(const TVector3 &v1, const TVector3 &v2, const TVector3 &v3)
Double_t Dot(const TVector3 &) const
TVector2 operator+(const TVector2 &v1, const TVector2 &v2)
Bool_t point_in_polygon_test_3d(const TGBinder &p1, const TPlane3 &plane, const TPoint3 &origin, const TPoint3 &pointOnPlane)
Double_t Angle(const TVector3 &v1, const TVector3 &v2)
Double_t Angle(const TVector3 &) const
return the angle w.r.t. another 3-vector
Bool_t intersect_polygons(const TGBinderA &p1, const TGBinderB &p2, const TPlane3 &plane1, const TPlane3 &plane2)
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.
std::vector< Int_t > PIndexList_t
Double_t dot(const TVector2 &v1, const TVector2 &v2)
Bool_t operator==(const Tuple2 &t1, const Tuple2 &t2)
TVector2 operator-(const TVector2 &v1, const TVector2 &v2)
TBaseMesh * BuildDifference(const TBaseMesh *leftOperand, const TBaseMesh *rightOperand)
TVector3 & operator*=(Double_t)
AMesh_t * build_union(const AMesh_t &meshA, const AMesh_t &meshB, Bool_t preserve)
static double p2(double t, double a, double b, double c)
void partition_mesh(CMesh &mesh, const TMesh &mesh2, const OverlapTable_t &table)
TBaseMesh * BuildIntersection(const TBaseMesh *leftOperand, const TBaseMesh *rightOperand)
TVector2 & operator/=(Double_t s)
void copy_mesh(const MeshA &source, MeshB &output)
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)
Bool_t intersect_2d_bounds_check(const TLine3 &l1, const TLine3 &l2, Int_t majAxis, Double_t &l1Param, Double_t &l2Param)
SEXP wrap(const TString &s)
std::vector< PIndexList_t > OverlapTable_t
Double_t operator*=(TVector2 const &v)
TPoint2 lerp(const TPoint2 &p1, const TPoint2 &p2, Double_t t)
TPlane3 compute_plane(const TGBinder &poly)
Double_t length(const TVector2 &v)
Double_t distance(const TPoint3 &p1, const TPoint3 &p2)
unsigned int r1[N_CITIES]
TVector3 & operator+=(const TVector3 &)
TConnectedMeshWrapper< AConnectedMesh_t > AConnectedMeshWrapper_t
Generic 3D primitive description class.
TVector3 & operator-=(const TVector3 &)
void classify_mesh(const TMesh &meshA, const TBBoxTree &aTree, CMesh &meshB)
static double p1(double t, double a, double b)
Bool_t intersect(const TPlane3 &p1, const TPlane3 &p2, TLine3 &output)
Double_t Dot(const TGLVector3 &v1, const TGLVector3 &v2)
Bool_t intersect_poly_with_line_2d(const TLine3 &l, const TGBinder &p1, const TPlane3 &plane, Double_t &a, Double_t &b)
I've modified some very nice bounding box tree code from Gino van der Bergen's Free Solid Library bel...
TVector2 & operator-=(TVector2 const &v)
TMesh< TestPolygon_t, TVertexBase > AMesh_t
PIndexList_t::iterator PIndexIt_t
TBBoxInternal * InternalPtr_t
Double_t distance2(const TPoint2 &p1, const TPoint2 &p2)
TMesh< TestPolygon_t, TCVertex > AConnectedMesh_t
Binding & operator=(OUT(*fun)(void))
TBaseMesh * ConvertToMesh(const TBuffer3D &buff)
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)
TLine3 polygon_mid_point_ray(const TGBinder &p1, const TPlane3 &plane)
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)
TVector3 Cross(const TVector3 &) const
Short_t Max(Short_t a, Short_t b)
Bool_t fuzzy_zero2(Double_t x)
AxisAngle::Scalar Distance(const AxisAngle &r1, const R &r2)
Distance between two rotations.
TVector2 & operator+=(TVector2 const &v)
ClassImp(TSlaveInfo) Int_t TSlaveInfo const TSlaveInfo * si
Used to sort slaveinfos by ordinal.
Bool_t fuzzy_zero(Double_t x)
TGLVector3 Cross(const TGLVector3 &v1, const TGLVector3 &v2)
Double_t Sqrt(Double_t x)
Bool_t fuzzy_equal(const TVector2 &v1, const TVector2 &v2)
std::complex< float_v > Z
TPolygonBase< TBlenderVProp, NullType_t > TestPolygon_t
static void output(int code)
Int_t which_side(const TGBinder &p1, const TPlane3 &plane1)
TVector2 operator*(const TVector2 &v, Double_t s)
Double_t Angle(const TVector2 &v1, const TVector2 &v2)
TBBox fit_bbox(const TGBinder &p1)
Double_t length2(const TVector2 &v)
Int_t compute_classification(const Double_t &distance, const Double_t &epsil)
unsigned int r2[N_CITIES]
TBaseMesh * BuildUnion(const TBaseMesh *leftOperand, const TBaseMesh *rightOperand)
void Vertex(const Double_t *v)