PolyIntersect.cpp
1 #include "PolyIntersect.hpp"
2 
3 using std::vector;
4 using std::max;
5 
6 IntersectFlags SegmentIntersection(Vector2D const& p0,Vector2D const& p1,
7  Vector2D const& q0,Vector2D const& q1,Vector2D &Intersection)
8 {
9  /*int n=static_cast<int>(p.size());
10  int m=static_cast<int>(q.size());
11  Vector2D p0(p[ploc]);
12  Vector2D p1(p[(ploc+1)%n]);
13  Vector2D q0(q[qloc]);
14  Vector2D q1(q[(qloc+1)%m]);*/
15  if(min(p0.x,p1.x)>max(q0.x,q1.x)||min(q1.x,q0.x)>max(p0.x,p1.x)||
16  min(p0.y,p1.y)>max(q1.y,q0.y)||min(q1.y,q0.y)>max(p0.y,p1.y))
17  return False;
18  boost::array<Vector2D,3> points;
19  points[0]=p0;
20  points[1]=p1;
21  points[2]=q0;
22  double temp=orient2d(TripleConstRef<Vector2D>(p0,p1,q0));
23  points[2]=q1;
24  if(temp*orient2d(TripleConstRef<Vector2D>(p0,p1,q1))>0)
25  return False;
26  points[0]=q0;
27  points[1]=q1;
28  points[2]=p0;
29  temp=orient2d(TripleConstRef<Vector2D>(q0,q1,p0));
30  points[2]=p1;
31  if(temp*orient2d(TripleConstRef<Vector2D>(q0,q1,p1))>0)
32  return False;
33  points[0]=Vector2D(0,0);
34  points[1]=p1-p0;
35  points[2]=q1-q0;
36  double d=orient2d(TripleConstRef<Vector2D>(Vector2D(0,0),
37  p1-p0,
38  q1-q0));
39  double eps = 1e-7*abs(p1 - p0)*abs(q1 - q0);
40  if(fabs(d)<0)
41  return Par;
42  double xi=((q0.x-q1.x)*(p0.x*p1.y-p1.x*p0.y)-(p0.x-p1.x)*
43  (q0.x*q1.y-q1.x*q0.y))/d;
44  double yi=((q0.y-q1.y)*(p0.x*p1.y-p1.x*p0.y)-(p0.y-p1.y)*
45  (q0.x*q1.y-q1.x*q0.y))/d;
46  Intersection.Set(xi,yi);
47  eps=1e-7*sqrt(abs(p1-p0)*abs(q1-q0));
48  if((xi+eps)<min(p0.x,p1.x)||(xi-eps)>max(p0.x,p1.x))
49  return Par;
50  if((xi+eps)<min(q0.x,q1.x)||(xi-eps)>max(q0.x,q1.x))
51  return Par;
52  if((yi+eps)<min(p0.y,p1.y)||(yi-eps)>max(p0.y,p1.y))
53  return Par;
54  if((yi+eps)<min(q0.y,q1.y)||(yi-eps)>max(q0.y,q1.y))
55  return Par;
56  return True;
57 }
58 
59 /*
60 vector<Vector2D> GetParEdge(Vector2D const& p0,Vector2D const& p1,
61  Vector2D const& q0,Vector2D const& q1)
62 {
63  vector<Vector2D> res(2);
64  if((p1.x-p0.x)*(p1.y-p0.y)>0)
65  {
66  res[0].Set(min(max(p0.x,p1.x),max(q0.x,q1.x)),
67  min(max(p0.y,p1.y),max(q0.y,q1.y)));
68  res[1].Set(max(min(p0.x,p1.x),min(q0.x,q1.x)),
69  max(min(p0.y,p1.y),min(q0.y,q1.y)));
70  }
71  else
72  {
73  res[0].Set(max(min(p0.x,p1.x),min(q0.x,q1.x)),
74  min(max(p0.y,p1.y),max(q0.y,q1.y)));
75  res[1].Set(min(max(p0.x,p1.x),max(q0.x,q1.x)),
76  max(min(p0.y,p1.y),min(q0.y,q1.y)));
77  }
78  return res;
79 }
80 */
81 
82 vector<Vector2D> ConvexIntersect(vector<Vector2D> const& poly0,vector<Vector2D>
83  const& poly1)
84 {
85  // poly0 is a poly1 is b
86  vector<Vector2D> res;
87  InFlags flag=UnKnown;
88  bool FirsPoint=true;
89  const int n=static_cast<int>(poly0.size());
90  const int m=static_cast<int>(poly1.size());
91  int p0index=0,p1index=0;
92  boost::array<Vector2D,3> AreaSignCheck;
93  int p0counter=0,p1counter=0;
94  do
95  {
96  AreaSignCheck[0]=poly0[static_cast<size_t>((p0index+n-1)%n)];
97  AreaSignCheck[1]=poly0[static_cast<size_t>(p0index)];
98  AreaSignCheck[2]=poly1[static_cast<size_t>(p1index)];
99  double bLeftOfa=orient2d
100  (TripleConstRef<Vector2D>
101  (poly0[static_cast<size_t>((p0index+n-1)%n)],
102  poly0[static_cast<size_t>(p0index)],
103  poly1[static_cast<size_t>(p1index)]));
104  AreaSignCheck[0]=poly1[static_cast<size_t>((p1index+m-1)%m)];
105  AreaSignCheck[1]=poly1[static_cast<size_t>(p1index)];
106  AreaSignCheck[2]=poly0[static_cast<size_t>(p0index)];
107  double aLeftOfb=orient2d
108  (TripleConstRef<Vector2D>
109  (poly1[static_cast<size_t>((p1index+m-1)%m)],
110  poly1[static_cast<size_t>(p1index)],
111  poly0[static_cast<size_t>(p0index)]));
112 
113  Vector2D intersect;
114  IntersectFlags inter=SegmentIntersection(poly0[static_cast<size_t>((p0index+n-1)%n)],
115  poly0[static_cast<size_t>(p0index)],
116  poly1[static_cast<size_t>(p1index)],
117  poly1[static_cast<size_t>((p1index+m-1)%m)],intersect);
118  if(inter==True)
119  {
120  res.push_back(intersect);
121  if(aLeftOfb>0)
122  flag=Pi;
123  else
124  flag=Qi;
125  if(FirsPoint)
126  {
127  FirsPoint=false;
128  p0counter=0;
129  p1counter=0;
130  }
131  }
132 
133  AreaSignCheck[0]=Vector2D(0,0);
134  AreaSignCheck[1]=poly0[static_cast<size_t>(p0index)]-poly0[static_cast<size_t>((p0index+n-1)%n)];
135  AreaSignCheck[2]=poly1[static_cast<size_t>(p1index)]-poly1[static_cast<size_t>((p1index+m-1)%m)];
136  double areasign=orient2d
137  (TripleConstRef<Vector2D>
138  (Vector2D(0,0),
139  poly0[static_cast<size_t>(p0index)]-poly0[static_cast<size_t>((p0index+n-1)%n)],
140  poly1[static_cast<size_t>(p1index)]-poly1[static_cast<size_t>((p1index+m-1)%m)]));
141  if(inter==Par)
142  {
143  if(ScalarProd(poly0[static_cast<size_t>(p0index)]-poly0[static_cast<size_t>((p0index+n-1)%n)],
144  poly1[static_cast<size_t>(p1index)]-poly1[static_cast<size_t>((p1index+m-1)%m)])<0)
145  {
146  // Do we need this? The area will be zero
147  /*res=GetParEdge(poly0[p0index],poly0[(p0index+n-1)%n],
148  poly1[p1index],poly1[(p1index+m-1)%m]);*/
149  return vector<Vector2D> ();
150  }
151  }
152  if((fabs(areasign)<1e-9)&&(aLeftOfb<0)&&(bLeftOfa<0))
153  return vector<Vector2D> ();
154  if((fabs(areasign)<1e-9)&&(fabs(aLeftOfb)<1e-9)&&(fabs(bLeftOfa)<1e-9))
155  {
156  if(flag==Pi)
157  {
158  ++p1counter;
159  p1index=(p1index+1)%m;
160  }
161  else
162  {
163  ++p0counter;
164  p0index=(p0index+1)%n;
165  }
166  }
167  if(areasign>=0)
168  {
169  if(bLeftOfa>0)
170  {
171  if(flag==Pi)
172  res.push_back(poly0[static_cast<size_t>(p0index)]);
173  ++p0counter;
174  p0index=(p0index+1)%n;
175  }
176  else
177  {
178  if(flag==Qi)
179  res.push_back(poly1[static_cast<size_t>(p1index)]);
180  ++p1counter;
181  p1index=(p1index+1)%m;
182  }
183  }
184  else
185  {
186  if(aLeftOfb>0)
187  {
188  if(flag==Qi)
189  res.push_back(poly1[static_cast<size_t>(p1index)]);
190  ++p1counter;
191  p1index=(p1index+1)%m;
192  }
193  else
194  {
195  if(flag==Pi)
196  res.push_back(poly0[static_cast<size_t>(p0index)]);
197  ++p0counter;
198  p0index=(p0index+1)%n;
199  }
200  }
201  }
202  while(((p0counter<n)||(p1counter<m))&&(p0counter<2*n)&&(p1counter<2*m));
203  return res;
204 }
SegmentIntersection
bool SegmentIntersection(Edge const &edge1, Edge const &edge2, Vector2D &Intersection, double eps=1e-8)
Calculates the intersection of two edges.
Definition: Edge.cpp:45
InFlags
InFlags
Flags for if the segment is inner or outer.
Definition: PolyIntersect.hpp:10
max
double max(vector< double > const &v)
returns the maximal term in a vector
Definition: utils.cpp:52
Vector2D::y
double y
Component in the y direction.
Definition: geometry.hpp:48
TripleConstRef
A collection of three identical references.
Definition: triplet.hpp:12
ConvexIntersect
vector< Vector2D > ConvexIntersect(vector< Vector2D > const &poly0, vector< Vector2D > const &poly1)
Calculates the intersection between two convex polygons.
Definition: PolyIntersect.cpp:82
PolyIntersect.hpp
Finds the intersection points between two polygons.
ScalarProd
double ScalarProd(Vector3D const &v1, Vector3D const &v2)
Scalar product of two vectors.
Definition: Vector3D.cpp:185
min
double min(vector< double > const &v)
Returns the minimal term in a vector.
Definition: utils.cpp:44
Vector2D::Set
void Set(double ix, double iy)
Set vector components.
Definition: geometry.cpp:39
abs
double abs(Vector3D const &v)
Norm of a vector.
Definition: Vector3D.cpp:44
IntersectFlags
IntersectFlags
Flags for there is an intersction or not or if the segments are parallel.
Definition: PolyIntersect.hpp:12
Vector2D
2D Mathematical vector
Definition: geometry.hpp:15
orient2d
double orient2d(const TripleConstRef< Vector2D > &points)
Checks whether 3 given points are on a counterclockwise circle, clockwise circle or colinear...
Definition: geotests.cpp:76
Vector2D::x
double x
Component in the x direction.
Definition: geometry.hpp:45
Generated on Thu Jun 4 2020 18:12:17 by   doxygen 1.8.13