|
ATLAS Offline Software
|
Go to the documentation of this file.
110 #define sqr(t) (t)*(t)
119 typedef std::map<unsigned int, Linebase*>
LineMap;
120 typedef std::priority_queue<Pointbase>
PQueue;
126 typedef std::map<unsigned int, std::set<unsigned int> >
AdjEdgeMap;
270 #define Absolute(a) ((a) >= 0.0 ? (a) : -(a))
286 #define Fast_Two_Sum_Tail(a, b, x, y) \
290 #define Fast_Two_Sum(a, b, x, y) \
291 x = (REAL) (a + b); \
292 Fast_Two_Sum_Tail(a, b, x, y)
294 #define Two_Sum_Tail(a, b, x, y) \
295 bvirt = (REAL) (x - a); \
297 bround = b - bvirt; \
298 around = a - avirt; \
301 #define Two_Sum(a, b, x, y) \
302 x = (REAL) (a + b); \
303 Two_Sum_Tail(a, b, x, y)
305 #define Two_Diff_Tail(a, b, x, y) \
306 bvirt = (REAL) (a - x); \
308 bround = bvirt - b; \
309 around = a - avirt; \
312 #define Two_Diff(a, b, x, y) \
313 x = (REAL) (a - b); \
314 Two_Diff_Tail(a, b, x, y)
319 #define Split(a, ahi, alo) \
320 c = (REAL) (1.0 * a); \
321 abig = (REAL) (c - a); \
325 #define Two_Product_Tail(a, b, x, y) \
326 Split(a, ahi, alo); \
327 Split(b, bhi, blo); \
328 err1 = x - (ahi * bhi); \
329 err2 = err1 - (alo * bhi); \
330 err3 = err2 - (ahi * blo); \
331 y = (alo * blo) - err3
333 #define Two_Product(a, b, x, y) \
334 x = (REAL) (a * b); \
335 Two_Product_Tail(a, b, x, y)
340 #define Two_One_Diff(a1, a0, b, x2, x1, x0) \
341 Two_Diff(a0, b , x_i, x0); \
342 Two_Sum( a1, x_i, x2, x1)
344 #define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) \
345 Two_One_Diff(a1, a0, b0, x_j, x_0, x0); \
346 Two_One_Diff(x_j, x_0, b1, x3, x2, x1)
379 REAL avirt, bround, around;
380 int eindex, findex, hindex;
386 if ((fnow > enow) == (fnow > -enow)) {
394 if ((eindex < elen) && (findex < flen)) {
395 if ((fnow > enow) == (fnow > -enow)) {
406 while ((eindex < elen) && (findex < flen)) {
407 if ((fnow > enow) == (fnow > -enow)) {
420 while (eindex < elen) {
428 while (findex < flen) {
436 if ((Q != 0.0) || (hindex == 0)) {
457 for (eindex = 1; eindex < elen; ++eindex) {
492 REAL acxtail, acytail, bcxtail, bcytail;
494 REAL detlefttail, detrighttail;
496 REAL B[4], C1[8], C2[12], D[16];
498 int C1length, C2length, Dlength;
505 REAL avirt, bround, around;
508 REAL ahi, alo, bhi, blo;
509 REAL err1, err2, err3;
513 acx = (
REAL) (pa[0] -
pc[0]);
515 acy = (
REAL) (pa[1] -
pc[1]);
521 Two_Two_Diff(detleft, detlefttail, detright, detrighttail,
522 B3,
B[2],
B[1],
B[0]);
528 if ((
det >= errbound) || (-
det >= errbound)) {
537 if ((acxtail == 0.0) && (acytail == 0.0)
538 && (bcxtail == 0.0) && (bcytail == 0.0)) {
544 det += (acx * bcytail + bcy * acxtail)
545 - (acy * bcxtail + bcx * acytail);
546 if ((
det >= errbound) || (-
det >= errbound)) {
568 return(D[Dlength - 1]);
574 REAL detsum, errbound;
576 detleft = (pa[0] -
pc[0]) * (
pb[1] -
pc[1]);
577 detright = (pa[1] -
pc[1]) * (
pb[0] -
pc[0]);
578 det = detleft - detright;
581 if (detright <= 0.0) {
584 detsum = detleft + detright;
586 }
else if (detleft < 0.0) {
587 if (detright >= 0.0) {
590 detsum = -detleft - detright;
598 if ((
det >= errbound) || (-
det >= errbound)) {
631 template <
class T,
class KeyType>
634 template <
class T,
class KeyType>
659 template <
class T,
class KeyType>
728 template <
class T,
class KeyType>
737 template <
class T,
class KeyType>
746 template <
class T,
class KeyType>
756 m_root = newNode; ++m_size;
760 KeyType
keys=
x->keyValue();
761 splay(
keys, m_root );
762 KeyType rootk=m_root->keyValue();
771 else if(
keys > rootk )
785 x->increaseKeyValue(1.0
e-10);
795 template <
class T,
class KeyType>
800 splay(
keys, m_root );
801 if( m_root->keyValue() !=
keys ) {
res=NULL;
return; }
805 if( m_root->m_left == NULL )
813 newTree->m_right = m_root->m_right;
823 template <
class T,
class KeyType>
828 splay(
keys, m_root );
829 KeyType rootk=m_root->keyValue();
830 if( rootk !=
keys ) {
return; }
832 if( m_root->m_left == NULL )
newTree = m_root->m_right;
839 newTree->m_right = m_root->m_right;
852 template <
class T,
class KeyType>
855 if( IsEmpty( ) ) {
min=NULL;
return; }
858 splay(
keys, m_root );
863 if( m_root->m_left == NULL )
newTree = m_root->m_right;
868 newTree->m_right = m_root->m_right;
879 template <
class T,
class KeyType>
882 if( IsEmpty( ) ) {
max=NULL;
return; }
885 splay(
keys, m_root );
890 if( m_root->m_left == NULL )
newTree = m_root->m_right;
895 newTree->m_right = m_root->m_right;
905 template <
class T,
class KeyType>
908 if( IsEmpty( ) ) {
min=NULL;
return; }
911 while(
ptr->m_left != NULL )
ptr =
ptr->m_left;
912 splay(
ptr->keyValue(), m_root );
919 template <
class T,
class KeyType>
922 if( IsEmpty( ) ) {
max=NULL;
return; }
925 while(
ptr->m_right != NULL )
ptr =
ptr->m_right;
926 splay(
ptr->keyValue(), m_root );
934 template <
class T,
class KeyType>
937 if( IsEmpty( ) ) {
res=NULL;
return; }
938 splay(
keys, m_root );
939 if( m_root->keyValue() !=
keys ) {
res=NULL;
return; }
948 template <
class T,
class KeyType>
951 if( IsEmpty( ) ) {
res=NULL;
return; }
952 splay(
keys, m_root );
954 if( m_root->data()->keyValue() <
keys)
res=m_root;
955 else if(m_root->m_left)
969 template <
class T,
class KeyType>
983 template <
class T,
class KeyType>
986 return m_root == NULL;
992 template <
class T,
class KeyType>
1010 template <
class T,
class KeyType>
1018 leftTreeMax = rightTreeMin = &
header;
1022 KeyType rkey=
t->keyValue();
1025 if(
t->m_left == NULL)
break;
1026 if( keys < t->m_left->keyValue() ) rotateWithLeftChild(
t );
1027 if(
t->m_left == NULL )
break;
1034 else if(
keys > rkey )
1036 if(
t->m_right == NULL )
break;
1037 if(
keys >
t->m_right->keyValue() ) rotateWithRightChild(
t );
1038 if(
t->m_right == NULL )
break;
1049 rightTreeMin->
m_left =
t->m_right;
1058 template <
class T,
class KeyType>
1070 template <
class T,
class KeyType>
1083 template <
class T,
class KeyType>
1086 if(
t !=
t->m_left )
1088 reclaimMemory(
t->m_left );
1089 reclaimMemory(
t->m_right );
1098 template <
class T,
class KeyType>
1101 if(
t ==
t->m_left )
1110 template<
class T,
class KeyType>
1116 PreOrder(Visit,
t->m_left);
1117 PreOrder(Visit,
t->m_right);
1125 template<
class T,
class KeyType>
1130 InOrder(Visit,
t->m_left);
1132 InOrder(Visit,
t->m_right);
1140 template<
class T,
class KeyType>
1146 InOrder(Visit,
t->m_left,
y);
1148 InOrder(Visit,
t->m_right,
y);
1157 template<
class T,
class KeyType>
1162 PostOrder(Visit,
t->m_left);
1163 PostOrder(Visit,
t->m_right);
1171 template<
class T,
class KeyType>
1174 if(subtree==NULL)
return 0;
1176 int rh=Height(subtree->
m_right);
1178 return (
lh>rh)?(++
lh):(++rh);
1324 node->data()->setKeyValue(
y);
1344 return (pa.
x==
pb.x && pa.
y==
pb.y);
1350 return( (pa.
y >
pb.y) || ( (pa.
y==
pb.y) && (pa.
x <
pb.x)) );
1356 return( (pa.
y <
pb.y) || ( (pa.
y==
pb.y) && (pa.
x >
pb.x)) );
1362 return !(pa.
x==
pb.x && pa.
y==
pb.y);
1422 Polygon(
const std::vector<double>&
x,
const std::vector<double>&
y);
1440 void set_contour (
const std::vector<double>&
x,
const std::vector<double>&
y);
1444 unsigned int prev(
const unsigned int&
i);
1445 unsigned int next(
const unsigned int&
i);
1456 void addDiagonal(
const unsigned int&
i,
const unsigned int& j);
1523 assert(
x.size()==
y.size());
1525 m_nVertices.push_back(
x.size() );
1526 unsigned int i = m_points.size()+1;
1529 for (
unsigned int j = 0; j < m_nVertices[m_ncontours]; ++j, ++
i)
1535 if(xx > m_xmax ) m_xmax=xx;
1536 if(xx < m_xmin ) m_xmin=xx;
1537 if(
yy > m_ymax ) m_ymax=
yy;
1538 if(
yy < m_ymin ) m_ymin=
yy;
1557 m_xmin = m_ymin = std::numeric_limits<double>::infinity();
1558 m_xmax = m_ymax = - std::numeric_limits<double>::infinity();
1562 init_vertices_and_lines();
1573 for(;
itp!=m_points.end(); ++
itp)
1579 for(; itl!=m_edges.end(); ++itl)
1591 unsigned int j(0),prevLoop(0),currentLoop(0);
1593 while (
i > m_nVertices[currentLoop] )
1595 prevLoop=currentLoop;
1599 if(
i==1 || (
i==m_nVertices[prevLoop]+1) ) j=m_nVertices[currentLoop];
1600 else if(
i <= m_nVertices[currentLoop] ) j=
i-1;
1610 unsigned int j(0),prevLoop(0),currentLoop(0);
1612 while (
i > m_nVertices[currentLoop] )
1614 prevLoop=currentLoop;
1618 if(
i < m_nVertices[currentLoop] ) j=
i+1;
1619 else if (
i==m_nVertices[currentLoop] )
1621 if( currentLoop==0) j=1;
1622 else j=m_nVertices[prevLoop]+1;
1638 for(;
it!=m_points.end(); ++
it)
1645 if(
p > pnext && pprev >
p )
1647 else if (
p > pprev && pnext >
p)
1651 double pa[2],
pb[2],
pc[2];
1653 pa[0]=m_points[idp]->x;
1654 pa[1]=m_points[idp]->y;
1656 pb[0]=m_points[
id]->x;
1657 pb[1]=m_points[
id]->y;
1659 pc[0]=m_points[idn]->x;
1660 pc[1]=m_points[idn]->y;
1669 m_qpoints.push(*(
it->second));
1671 m_startAdjEdgeMap[
id].insert(
id);
1683 m_edges[diag->
id()]=diag;
1685 m_startAdjEdgeMap[
i].insert(diag->
id());
1686 m_startAdjEdgeMap[j].insert(diag->
id());
1688 m_diagonals[diag->
id()]=diag;
1697 double y=m_points[
i]->y;
1700 m_edges[
i]->setHelper(
i);
1701 m_edges[
i]->setKeyValue(
y);
1702 m_edgebst.Insert(m_edges[
i]);
1711 double y=m_points[
i]->y;
1714 unsigned int previ=prev(
i);
1716 unsigned int helper=m_edges[previ]->helper();
1720 m_edgebst.Delete(edge->
keyValue());
1729 double x=m_points[
i]->x,
y=m_points[
i]->y;
1733 m_edgebst.FindMaxSmallerThan(
x, leftnode);
1740 m_edges[
i]->setHelper(
i);
1741 m_edges[
i]->setKeyValue(
y);
1742 m_edgebst.Insert(m_edges[
i]);
1750 double x=m_points[
i]->x,
y=m_points[
i]->y;
1753 unsigned int previ=prev(
i);
1754 unsigned int helper=m_edges[previ]->helper();
1756 m_edgebst.Delete(m_edges[previ]->
keyValue());
1759 m_edgebst.FindMaxSmallerThan(
x, leftnode);
1773 double y=m_points[
i]->y;
1776 unsigned int previ=prev(
i);
1777 unsigned int helper=m_edges[previ]->helper();
1780 m_edgebst.Delete(m_edges[previ]->
keyValue());
1781 m_edges[
i]->setHelper(
i);
1782 m_edges[
i]->setKeyValue(
y);
1783 m_edgebst.Insert(m_edges[
i]);
1792 double x=m_points[
i]->x,
y=m_points[
i]->y;
1796 m_edgebst.FindMaxSmallerThan(
x, leftnode);
1814 std::cout<<
"Please check your input polygon:\n1)orientations?\n2)duplicated points?\n";
1818 while(!m_qpoints.empty())
1822 unsigned int id=
vertex.id;
1833 std::cout<<
"No duplicated points please!\n";
1849 double dxab = pa[0] -
pb[0];
1850 double dyab = pa[1] -
pb[1];
1852 double dxcb =
pc[0] -
pb[0];
1853 double dycb =
pc[1] -
pb[1];
1855 double dxab2 = dxab * dxab;
1856 double dyab2 = dyab * dyab;
1857 double dxcb2 = dxcb * dxcb;
1858 double dycb2 = dycb * dycb;
1859 double ab = dxab2 + dyab2;
1860 double cb = dxcb2 + dycb2;
1862 double cosb = dxab * dxcb + dyab * dycb;
1863 double denom = sqrt( ab * cb);
1878 std::set<unsigned int> edges=m_startAdjEdgeMap[eid];
1879 assert(!edges.empty());
1881 unsigned int nexte=0;
1882 if( edges.size() == 1 ) nexte=*(edges.begin());
1883 else if( edges.size() > 1 )
1885 unsigned int nexte_ccw(0), nexte_cw(0);
1890 for(;
it!=edges.end(); ++
it)
1892 if(*
it==edge->
id())
continue;
1893 double A[2],
B[2],
C[2];
1897 if(edge->
endPoint(1)!=m_edges[*
it]->endPoint(0)) m_edges[*
it]->reverse();
1898 C[0]=m_edges[*
it]->endPoint(1)->x;
C[1]=m_edges[*
it]->endPoint(1)->y;
1901 double cosb=angleCosb(
A,
B,
C);
1903 if(
area > 0 &&
max < cosb ) { nexte_ccw=*
it;
max=cosb; }
1904 else if(
min > cosb ) { nexte_cw=*
it;
min=cosb; }
1907 nexte = (nexte_ccw!=0) ? nexte_ccw : nexte_cw;
1920 while( edges.size() > m_diagonals.size() )
1928 poly.push_back(startp->
id);
1932 endp=
next->endPoint(1);
1935 edges.erase(
next->id());
1936 m_startAdjEdgeMap[
next->endPoint(0)->id].erase(
next->id());
1938 if(endp==startp)
break;
1939 poly.push_back(endp->
id);
1941 unsigned int nexte=selectNextEdge(
next);
1945 std::cout<<
"Please check your input polygon:\n";
1946 std::cout<<
"1)orientations?\n2)with duplicated points?\n3)is a simple one?\n";
1951 if(
next->endPoint(0) !=endp )
next->reverse();
1954 m_mpolys.push_back(
poly);
1966 for(; itnext=
it,
it!=mpoly.end(); ++
it)
1969 if(itnext==mpoly.end()) itnext=mpoly.begin();
1971 point.
left=(point > pointnext)?
true:
false;
1972 qvertex.push(point);
1975 std::stack<internal_poltrig::Pointbase> spoint;
1976 for(
int i=0;
i<2; ++
i) { spoint.push(qvertex.top()); qvertex.pop(); }
1978 while ( qvertex.size() > 1 )
1983 if(topQueuePoint.
left!=topStackPoint.
left)
1985 while ( spoint.size() > 1 )
1991 v[0]=topQueuePoint.
id;
1994 sort(
v.begin(),
v.end());
1995 m_triangles.push_back(
v);
1999 spoint.push(topStackPoint);
2000 spoint.push(topQueuePoint);
2004 while( spoint.size() > 1 )
2009 spoint.push(stack1Point);
2010 double pa[2],
pb[2],
pc[2];
2011 pa[0]=topQueuePoint.
x; pa[1]=topQueuePoint.
y;
2012 pb[0]=stack2Point.
x;
pb[1]=stack2Point.
y;
2013 pc[0]=stack1Point.
x;
pc[1]=stack1Point.
y;
2016 bool left=stack1Point.
left;
2017 if( (
area > 0 && left) || (
area < 0 && !left ) )
2020 v[0]=topQueuePoint.
id;
2021 v[1]=stack2Point.
id;
2022 v[2]=stack1Point.
id;
2023 sort(
v.begin(),
v.end());
2024 m_triangles.push_back(
v);
2029 spoint.push(topQueuePoint);
2038 while( spoint.size() !=1 )
2045 v[0]=lastQueuePoint.
id;
2048 sort(
v.begin(),
v.end());
2049 m_triangles.push_back(
v);
2058 partition2Monotone();
2061 for(;
it!=m_mpolys.end(); ++
it)
2062 triangulateMonotone(*
it);
2073 const std::vector<double>& polygon_ycoords)
JetConstituentVector::iterator iterator
void Find(const KeyType &keys, BTreeNode< T, KeyType > *&res)
void InOrder(void(*Visit)(BTreeNode< T, KeyType > *u, double y), double y)
const Triangles * triangles()
std::map< unsigned int, Linebase * > LineMap
#define Two_Product(a, b, x, y)
void PostOrder(void(*Visit)(BTreeNode< T, KeyType > *u))
unsigned int next(const unsigned int &i)
std::list< Monopoly > Monopolys
std::map< unsigned int, Pointbase * > PointbaseMap
BTreeNode< T, KeyType > * clone(BTreeNode< T, KeyType > *t) const
void DeleteMin(BTreeNode< T, KeyType > *&min)
std::list< Triangle > Triangles
std::vector< unsigned int > Triangle
void Delete(const KeyType &keys)
Pointbase(const double &xx, const double &yy)
internal_poltrig::PQueue m_qpoints
const Triangles * triangles() const
std::vector< ALFA_RawData_p1 > t0
internal_poltrig::EdgeBST m_edgebst
std::vector< ALFA_RawDataCollection_p1 > t1
void handleMergeVertex(const unsigned int &)
void set_contour(const std::vector< double > &x, const std::vector< double > &y)
BTreeNode< T, KeyType > * Right(BTreeNode< T, KeyType > *node)
std::list< Triangle > Triangles
internal_poltrig::LineMap m_edges
std::map< unsigned int, std::set< unsigned int > > AdjEdgeMap
std::priority_queue< Pointbase > PQueue
void handleStartVertex(const unsigned int &)
double dist_sqr(const Pointbase &sp, const Pointbase &ep)
BTreeNode(const T &data, BTreeNode *lt, BTreeNode *rt)
REAL orient2d(REAL *pa, REAL *pb, REAL *pc)
void triangulateMonotone(internal_poltrig::Monopoly &mpoly)
bool operator!=(const Pointbase &pa, const Pointbase &pb)
internal_poltrig::AdjEdgeMap m_startAdjEdgeMap
Pointbase * endPoint(const int &i) const
void handleEndVertex(const unsigned int &)
#define Two_Sum(a, b, x, y)
BTreeNode< T, KeyType > * Left(BTreeNode< T, KeyType > *node)
@ u
Enums for curvilinear frames.
void PreOrder(void(*Visit)(BTreeNode< T, KeyType > *u))
void InOrder(void(*Visit)(BTreeNode< T, KeyType > *u))
friend bool operator==(const Pointbase &, const Pointbase &)
void init_vertices_and_lines()
PolygonTriangulator(const std::vector< double > &polygon_xcoords, const std::vector< double > &polygon_ycoords)
constexpr double poly(const double x)
void FindMaxSmallerThan(const KeyType &keys, BTreeNode< T, KeyType > *&res)
internal_poltrig::Monopolys m_mpolys
void partition2Monotone()
double angleCosb(double *A, double *B, double *C)
void setKeyValue(const double &y)
internal_poltrig::PointbaseMap m_points
bool operator>(const Pointbase &pa, const Pointbase &pb)
Pointbase(const int &idd, const double &xx, const double &yy)
void UpdateKey(BTreeNode< Linebase *, double > *node, double y)
BTreeNode< T, KeyType > * Root()
unsigned int selectNextEdge(internal_poltrig::Linebase *edge)
std::pair< std::vector< unsigned int >, bool > res
void handleSplitVertex(const unsigned int &)
void increaseKeyValue(const double &diff)
void SetVisited(const bool &visited)
REAL estimate(const int &elen, REAL *e)
BTreeNode< T, KeyType > * m_root
internal_poltrig::PointbaseMap & points()
void rotateWithRightChild(BTreeNode< T, KeyType > *&k1) const
unsigned int prev(const unsigned int &i)
void FindMax(BTreeNode< T, KeyType > *&max)
bool operator==(const Pointbase &pa, const Pointbase &pb)
const SplayTree & operator=(const SplayTree &rhs)
#define Fast_Two_Sum(a, b, x, y)
internal_poltrig::LineMap & edges()
void addDiagonal(const unsigned int &i, const unsigned int &j)
REAL orient2dadapt(REAL *pa, REAL *pb, REAL *pc, const REAL &detsum)
Pointbase(const double &xx, const double &yy, const Type &ttype)
#define Two_Diff_Tail(a, b, x, y)
Pointbase(const int &idd, const double &xx, const double &yy, const Type &ttype)
std::list< unsigned int > Monopoly
void InOrder(void(*Visit)(BTreeNode< T, KeyType > *, double y), BTreeNode< T, KeyType > *t, double y)
int fast_expansion_sum_zeroelim(const int &elen, REAL *e, const int &flen, REAL *f, REAL *h)
#define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0)
void newTree(MVAUtils::ForestTMVA *forest, const TMVA::DecisionTreeNode *node, float weight, bool isRegression, bool useYesNoLeaf)
void DeleteMax(BTreeNode< T, KeyType > *&max)
friend bool operator!=(const Pointbase &, const Pointbase &)
Linebase & operator=(const Linebase &)=delete
SplayTree< Linebase *, double > EdgeBST
Polygon(const std::vector< double > &x, const std::vector< double > &y)
Pointbase & operator=(const Pointbase &)=default
friend bool operator>(const Pointbase &, const Pointbase &)
bool operator<(const Pointbase &pa, const Pointbase &pb)
void rotateWithLeftChild(BTreeNode< T, KeyType > *&k2) const
void reclaimMemory(BTreeNode< T, KeyType > *t) const
std::vector< unsigned int > m_nVertices
internal_poltrig::LineMap m_diagonals
void setHelper(const unsigned int &i)
std::vector< unsigned > Triangle
void splay(const KeyType &keys, BTreeNode< T, KeyType > *&t) const
void handleRegularVertexDown(const unsigned int &)
void FindMin(BTreeNode< T, KeyType > *&min)
void handleRegularVertexUp(const unsigned int &)
friend bool operator<(const Pointbase &, const Pointbase &)