1 // Copyright (C) 2007-2016 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 #ifndef __POLYGONALGORITHMS_TXX__
20 #define __POLYGONALGORITHMS_TXX__
22 #include "PolygonAlgorithms.hxx"
23 #include "InterpolationUtils.hxx"
28 namespace INTERP_KERNEL
31 PolygonAlgorithms<DIM>::PolygonAlgorithms(double epsilon, double precision)//: (0)
33 _is_in_intersection = false;
35 _precision = precision;
37 /*************************************************************/
38 /* Computes the 3D intersection between two COPLANAR */
39 /* Segments [A,B] and [C,D], stores the result in V. */
40 /* If A belongs to [CD] then the vertex E (preceeding A) */
41 /* is used to decide if the crossing is real. If A coincides */
42 /* with C or D, a special treatment is performed */
43 /*************************************************************/
45 bool PolygonAlgorithms<DIM>::intersectSegmentSegment(const double * A, const double * B, const double * C,
46 const double * D, const double * E, double * V)
48 double AB[DIM], DC[DIM], AC[DIM], det, t1, t2, inv_det;
50 /******* Initialisation of the linear system t1*AB+t2*DC=AC ***********/
51 for(int idim=0;idim<DIM;idim++)
53 AB[idim] = B[idim]-A[idim];//B-A
54 DC[idim] = C[idim]-D[idim];//C-D
55 AC[idim] = C[idim]-A[idim];//C-A
58 /******* Resolution of the linear system t1*AB+t2*DC=AC ***********/
59 det = determinant(AB,DC);//determinant of the first two coordinates
60 if(fabs(det) >_epsilon)
63 t1 = determinant(AC,DC)*inv_det;//solves the linear system t1*AB+t2*DC=AC
64 t2 = determinant(AB,AC)*inv_det;//solves the linear system t1*AB+t2*DC=AC
72 if(distance2<DIM>(A,D)<_epsilon)
73 crossprod<DIM>(A,C,E,_vdouble);//store the crossprod between vectors AC and AE (E=vertex preceding A)
74 return false;//case of paralell segments
76 case 3://beware AB and CD may belong to a vertical plane
77 det = determinant(&AB[1],&DC[1]);//determinant of the last two coefficients
78 if(fabs(det) > _epsilon)
81 t1=(AC[1]*DC[DIM-1]-AC[DIM-1]*DC[1])*inv_det;
82 t2=(AB[1]*AC[DIM-1]-AB[DIM-1]*AC[1])*inv_det;
84 else //beware AB and CD may belong to a plane y = constant
86 det = AB[0]*DC[DIM-1]-AB[DIM-1]*DC[0];
87 if(fabs(det) > _epsilon)
90 t1=(AC[0]*DC[DIM-1]-AC[DIM-1]*DC[0])*inv_det;
91 t2=(AB[0]*AC[DIM-1]-AB[DIM-1]*AC[0])*inv_det;
95 if(distance2<DIM>(A,D)<_epsilon)
96 crossprod<DIM>(A,C,E,_vdouble);//store the crossprod between vectors AC and AE (E=vertex preceding A)
97 return false;//case of paralell segments
103 if(t1>_precision && t1<1-_precision)
105 if( t2>_precision && t2<1-_precision)
107 for(int idim=0;idim<DIM;idim++) V[idim]=t1*AB[idim] + A[idim];
111 else if(fabs(t1) <= _precision)
113 if( t2>_precision && t2<1-_precision)//vertex on an edge
117 crossprod<DIM>(A,D,B,V12);
118 crossprod<DIM>(A,D,E,V34);
119 double same_side =dotprod<DIM>(V12, V34);
120 if( same_side < -_epsilon ) // <= epsilon or 0 ?//crossing
122 for(int idim=0;idim<DIM;idim++) V[idim]=A[idim];
125 else if( same_side > _epsilon ) _terminus= !_is_in_intersection;//reflexion
126 else //separation of overlaping edges
128 if(_Inter.empty() ) _terminus=true;
129 else if(!_is_in_intersection)
131 for(int idim=0;idim<DIM;idim++) V[idim]=A[idim];
136 else if(fabs(t2-1) <= _precision)//vertex on a vertex (A=D), first run
137 crossprod<DIM>(A,C,E,_vdouble);//store the angle between vectors AC and AE (E=vertex preceding A)
138 else if(fabs(t2) <= _precision)//vertex on a vertex (A=C), second run
140 double Vdoublebis[DIM];
141 //crossprod<DIM>(A,C,E,_vdouble);
142 crossprod<DIM>(A,B,D,Vdoublebis);
143 double in_between =dotprod<DIM>(Vdoublebis,_vdouble);
144 if(in_between>_epsilon)//crossing
146 for(int idim=0;idim<DIM;idim++) V[idim]=A[idim];
149 else if(fabs(in_between)<=_epsilon && dotprod<DIM>(Vdoublebis,Vdoublebis) > _epsilon)
150 //ie _vdouble=0, separation of overlaping edges at a double point
152 //crossprod<DIM>(A,E,B,_vdouble);
153 if(dotprod<DIM>(_vdouble,Vdoublebis) >=_epsilon )//crossing
155 if(_Inter.empty()) _terminus=true;
156 else if(!_is_in_intersection)
158 for(int idim=0;idim<DIM;idim++) V[idim]=A[idim];
168 /*************************************************************/
169 /* adds vertex i to the list inter and updates _End_segments */
170 /* i is the local index of the current vertex */
171 /*************************************************************/
173 inline void PolygonAlgorithms<DIM>::addNewVertex( int i, int i_glob, int i_next_glob, int i_prev_glob,
176 /* Question:Should we add vertex i to the front or back ? */
177 if( _End_segments[1].second == i_glob)
179 for(int idim=0;idim<DIM;idim++) _Inter.push_back(P[DIM*i+idim]);
180 _End_segments[1] = std::make_pair(i_glob, i_next_glob);
184 for(int idim=DIM-1;idim>-1;idim--) _Inter.push_front(P[DIM*i+idim]);
185 _End_segments[0] = std::make_pair(i_glob, i_next_glob);
189 /************************************************************/
190 /* adds a crossing between two segments starting at i and j */
191 /* to the double ended list inter in the correct order */
192 /* according to endsegments, updates _End_segments */
193 /************************************************************/
195 inline void PolygonAlgorithms<DIM>::addCrossing( double * ABCD, std::pair< int,int > i_i_next,
196 std::pair< int,int > j_j_next)
200 if(_End_segments[0] ==i_i_next)
202 for(int idim=DIM-1;idim>=0;idim--) _Inter.push_front(ABCD[idim]);
203 _terminus= (_End_segments[1]== j_j_next);
204 _End_segments[0] = j_j_next;
208 if( _End_segments[0]== j_j_next)
210 for(int idim=DIM-1;idim>=0;idim--) _Inter.push_front(ABCD[idim]);
211 _terminus= (_End_segments[1]== i_i_next);
212 _End_segments[0] = i_i_next;
216 for(int idim=0;idim<DIM;idim++) _Inter.push_back(ABCD[idim]);
217 _End_segments[1] = (_End_segments[1]== i_i_next) ? j_j_next : i_i_next;
223 for(int i=0;i<DIM;i++) _Inter.push_back(ABCD[i]);
224 _End_segments.push_back(i_i_next);
225 _End_segments.push_back(j_j_next);
229 /*******************************************************/
230 /* checks the possible crossing between segments [A,B] */
231 /* (with end-point global indices i and i_next) */
232 /* and [C,D] (end-point global indices j and j_next). */
233 /* If no intersection is detected, checks whether B is */
234 /* inside the quadrangle AEDC. */
235 /* Updates the lists inter and _End_segments */
236 /*******************************************************/
239 void PolygonAlgorithms<DIM>::addCrossing0(const double * A, const double * B, int i, int i_next,
240 const double * C, const double * D, int j, int j_next)
243 if(intersectSegmentSegment(A,B,C,D,ABCD, ABCD))
244 //fifth and sixth arguments are useless here
246 /* Updating _End_segments */
247 std::pair< int,int > i_i_next = std::make_pair(i, i_next);
248 std::pair< int,int > j_j_next = std::make_pair(j, j_next);
249 if( _End_segments[0] == i_i_next)
251 for(int idim=DIM-1;idim>-1;idim--) _Inter.push_front(ABCD[idim]);
252 _End_segments[0] = j_j_next;
256 for(int idim=0;idim<DIM;idim++) _Inter.push_back(ABCD[idim]);
257 _End_segments[1] = j_j_next;
258 _terminus = _End_segments[0]== j_j_next;
261 /* Updating _Status */
262 _Status.insert(make_pair(i_next,std::make_pair(i, false)));
263 std::multimap< int, std::pair< int,bool> >::iterator mi =_Status.find(j_next);
264 ((* mi).second).second= !((* mi).second).second;
266 else _Status.insert(std::make_pair(i_next,std::make_pair(i,true)));
269 /*******************************************************/
270 /* adds the possible crossings between segments [A,B] (with end-point global indices i and i_next) */
271 /*and segments [C,D] and [E,F] to the list inter and updates _End_segments */
272 /* In cases of ambiguity, the vertex G is used to decide wether the crossing should be accepted */
273 /*******************************************************/
275 inline void PolygonAlgorithms<DIM>::addCrossings( const double * A, const double * B, int i , int i_next,
276 const double * C, const double * D, int j1, int j2,
277 const double * E, const double * F, int j3, int j4,
282 std::multimap< int, std::pair< int,bool> >::iterator mi;
284 if(intersectSegmentSegment(A,B,C,D,G,ABCD))
286 if(intersectSegmentSegment(A,B,E,F,G,ABEF))
291 addCrossing(ABCD, std::make_pair(i, i_next), std::make_pair(j1, j2));
292 addCrossing(ABEF, std::make_pair(i, i_next), std::make_pair(j3, j4));
296 addCrossing(ABEF, std::make_pair(i, i_next), std::make_pair(j3, j4));
297 addCrossing(ABCD, std::make_pair(i, i_next), std::make_pair(j1, j2));
299 _Status.insert(std::make_pair(i_next,std::make_pair(i, _is_in_intersection)));
301 ((* mi).second).second= !((* mi).second).second;
303 ((* mi).second).second= !((* mi).second).second;
307 addCrossing(ABCD, std::make_pair( i, i_next), std::make_pair(j1,j2));
308 _Status.insert(std::make_pair(i_next,std::make_pair(i, !_is_in_intersection)));
310 ((* mi).second).second= !((* mi).second).second;
315 if(intersectSegmentSegment(A,B,E,F,G, ABEF))
317 addCrossing(ABEF, std::make_pair( i, i_next), std::make_pair( j3, j4));
318 _Status.insert(std::make_pair(i_next,std::make_pair(i, !_is_in_intersection)));
320 ((* mi).second).second= !((* mi).second).second;
322 else _Status.insert(std::make_pair(i_next,std::make_pair(i, _is_in_intersection)));
327 /* define various indices required in the function intersect_conv_polygon */
328 /* vertices from the both polygons are supposed to be present in the status */
330 inline void PolygonAlgorithms<DIM>::defineIndices(int& i_loc, int& i_next, int& i_prev,
331 const double *& Poly1, const double *& Poly2,
332 int& j1, int& j1_glob, int& j2, int& j2_glob,
333 int& j3, int& j3_glob, int& j4, int& j4_glob,
334 int& i_glob, int& i_next_glob, int& i_prev_glob,
335 const double * P_1, const double * P_2,
336 int N1, int N2, int sign)
346 std::multimap< int, std::pair< int,bool> >::reverse_iterator mi1=_Status.rbegin();
347 j1_glob=((*mi1).second).first;
349 j2_glob=(*mi1).first;
352 j3_glob=((*mi1).second).first;
354 j4_glob=(*mi1).first;
364 std::multimap< int, std::pair< int,bool> >::iterator mi2= _Status.begin();
365 j1_glob=((*mi2).second).first;
367 j2_glob=(*mi2).first;
370 j3_glob=((*mi2).second).first;
372 j4_glob=(*mi2).first;
375 i_loc = i_glob-shift;
376 i_next = (i_next_glob-shift+N0)%N0;//end-point of segment starting at i
377 i_prev = (i_prev_glob-shift+N0)%N0;
378 i_next_glob = i_next+shift;
379 i_prev_glob = i_prev+shift;
380 //warning: sign is either 1 or -1;
381 //To do: test and remove from Convex_intersecor.cxx
382 // while(distance2<DIM>(&Poly1[DIM*i_loc],&Poly1[DIM*i_next])< _epsilon && i_next != i_loc)
383 // i_next =(i_next+sign+N0)%N0;
384 // while(distance2<DIM>(&Poly1[DIM*i_loc],&Poly1[DIM*i_prev])< _epsilon && i_prev != i_loc)
385 // i_prev =(i_prev+sign+N0)%N0;
387 /*******************************************************/
388 /* computes the vertices of the intersection of two COPLANAR */
389 /* simple (no dble points)convex polygons using line sweep algorithm */
390 /* P1 and P2 contain the 3D coordinates of the successive vertices */
391 /*******************************************************/
393 std::deque< double > PolygonAlgorithms<DIM>::intersectConvexPolygons(const double* P_1,const double* P_2,
396 int i_loc, i_glob, j1, j1_glob, j2,j2_glob, j3, j3_glob, j4,j4_glob,
397 i_prev, i_prev_glob, i_next, i_next_glob, nb_prev, sign, idim;
398 const double * Poly1, * Poly2;
399 bool four_neighbours=false;
400 _terminus = N1 < 3 || N2<3;
402 /* list of future events ordered according to their coordinates (x,y,z) (lexicographical order) */
403 std::multimap< const double *, int, VertexLess<DIM> > mmap_events;
404 typename std::list< std::pair< const double *, int > >::iterator mi1,mi2;
406 std::multimap< int, std::pair< int,bool> >::iterator mi;
408 /********** Initalisation of events with P1 and P2 vertices ************/
409 for(i_loc=0;i_loc<N1;i_loc++)
410 mmap_events.insert(std::make_pair(&P_1[DIM*i_loc],i_loc));
411 for(i_loc=0;i_loc<N2;i_loc++)
412 mmap_events.insert(std::make_pair(&P_2[DIM*i_loc],i_loc+N1));
414 std::list< std::pair< const double *, int > > events(mmap_events.begin(),mmap_events.end());
418 /******** Treatment of the first vertex ********/
420 i_glob = (* mi1).second;
421 bool which_start = i_glob < N1;
422 if(i_glob < N1){ i_next_glob = (i_glob +1)%N1; i_prev_glob = (i_glob -1+N1)%N1;}
423 else{ i_next_glob = (i_glob-N1+1)%N2 + N1;i_prev_glob = (i_glob-N1-1+N2)%N2 + N1;}
424 _Status.insert(std::make_pair(i_next_glob,std::make_pair(i_glob, false)));
425 _Status.insert(std::make_pair(i_prev_glob,std::make_pair(i_glob, false)));
427 //std::cout<< "nb_prev= "<< 0 << " i_glob= " << i_glob << std::endl;
429 /******* Loop until the second polygon is reached *******/
430 while( !four_neighbours)
432 i_glob=(* mi1).second;//global index of vertex i
433 nb_prev = _Status.count(i_glob);//counts the number of segments ending at i
435 //std::cout<< "nb_prev= "<< nb_prev << " i_glob= " << i_glob << std::endl;
439 mi=_Status.find(i_glob);// pointer to the segment ending at i
440 i_prev_glob = ((*mi).second).first;//starting point of the segment ending at i
441 i_next= (i_prev_glob - i_glob > 0) == (abs(i_prev_glob - i_glob) == 1) ? i_glob - 1 : i_glob + 1;
442 if(i_glob < N1) i_next_glob = (i_next +N1)%N1;
443 else i_next_glob = (i_next-N1+N2)%N2 + N1;
445 _Status.insert(std::make_pair(i_next_glob,std::make_pair(i_glob, false)));
451 if( (i_glob < N1) != which_start)
455 /* detection of double points */
456 if(distance2<DIM>((* mi1).first, (*mi2).first) > _epsilon)
457 four_neighbours = true;
458 else /* Rare pothological case: */
460 const std::pair< const double *, int > next_pt= *mi2;
462 mi1=events.insert(mi1,next_pt);
467 throw Exception("intersectConvexPolygon: sequence of nodes does not describe a simple polygon (1)");
470 /******** Loop until a terminal point or crossing is reached ************/
473 //std::cout<< "nb_prev= "<< nb_prev<< " nb_inter= " << _Inter.size()/DIM << std::endl;
477 mi=_Status.find(i_glob);// pointer to the segment ending at i
478 i_prev_glob = ((*mi).second).first;//starting point of the segment ending at i
479 sign = (i_prev_glob - i_glob > 0) == (abs(i_prev_glob - i_glob) == 1) ? - 1 : + 1;
480 i_next_glob = i_glob+sign;
481 _is_in_intersection = ((*mi).second).second;//boolean that tells if i is in the intersection
483 defineIndices(i_loc,i_next,i_prev, Poly1,Poly2,
484 j1,j1_glob,j2,j2_glob,j3,j3_glob,j4,j4_glob,
485 i_glob,i_next_glob,i_prev_glob, P_1,P_2, N1, N2, sign);
486 if( _is_in_intersection ) addNewVertex(i_loc, i_glob, i_next_glob, i_prev_glob, Poly1);
487 addCrossings(&Poly1[DIM*i_loc], &Poly1[DIM*i_next], i_glob, i_next_glob,
488 &Poly2[DIM*j1] , &Poly2[DIM*j2] , j1_glob,j2_glob,
489 &Poly2[DIM*j3] , &Poly2[DIM*j4] , j3_glob,j4_glob, &Poly1[DIM*i_prev]);
494 if(i_glob < N1) for(idim=0;idim<DIM;idim++) _Inter.push_back(P_1[DIM*i_glob+idim]);
495 else for(idim=0;idim<DIM;idim++) _Inter.push_back(P_2[DIM*(i_glob-N1)+idim]);
498 case 0 ://To do if possible : remove this case from here
499 if(_Inter.empty() && (i_glob < N1) != which_start){
500 i_next_glob = i_glob+1;
501 i_prev_glob = i_glob-1;
502 defineIndices(i_loc,i_next,i_prev, Poly1,Poly2,
503 j1,j1_glob,j2,j2_glob,j3,j3_glob,j4,j4_glob,
504 i_glob,i_next_glob,i_prev_glob, P_1,P_2, N1, N2, 1);
505 double V12[DIM], V34[DIM];
506 double inside = check_inside<DIM>(&Poly1[DIM*i_loc],&Poly2[DIM*j1],&Poly2[DIM*j2],
507 &Poly2[DIM*j3], &Poly2[DIM*j4],V12, V34);
508 _is_in_intersection=( inside < _epsilon ); // <= epsilon or 0 ?
510 if(fabs(inside) > _epsilon)//vertex clearly inside or outside
512 //std::cout<<"coucou1" << std::endl;
513 if( _is_in_intersection)
515 for(int iidim=0;iidim<DIM;iidim++)
516 _Inter.push_back(Poly1[DIM*i_loc+iidim]);
517 _End_segments.push_back(std::make_pair(i_glob,i_next_glob));
518 _End_segments.push_back(std::make_pair(i_glob,i_prev_glob));
520 addCrossings(&Poly1[DIM*i_loc], &Poly1[DIM*i_next], i_glob, i_next_glob,
521 &Poly2[DIM*j1] , &Poly2[DIM*j2] , j1_glob,j2_glob,
522 &Poly2[DIM*j3] , &Poly2[DIM*j4] , j3_glob,j4_glob, &Poly1[DIM*i_prev]);
523 addCrossings(&Poly1[DIM*i_loc], &Poly1[DIM*i_prev], i_glob, i_prev_glob,
524 &Poly2[DIM*j1] , &Poly2[DIM*j2] , j1_glob,j2_glob,
525 &Poly2[DIM*j3] , &Poly2[DIM*j4] , j3_glob,j4_glob, &Poly1[DIM*i_next]);
527 else //vertex on an edge
529 //std::cout<<"coucou2" << std::endl;
530 bool is_inside_next, is_inside_prev;
531 double Vnext[DIM], Vprev[DIM];
532 for(idim=0;idim<DIM;idim++) _Inter.push_back(Poly1[DIM*i_loc+idim]);
534 if(dotprod<DIM>(V34,V34) > _epsilon)//vertex i on edge (j1,j2), not on (j3,j4)
536 crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j2], &Poly1[DIM*i_next],Vnext);
537 crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j2], &Poly1[DIM*i_prev],Vprev);
538 is_inside_next= (dotprod<DIM>(Vnext,V34)<0);
539 is_inside_prev= (dotprod<DIM>(Vprev,V34)<0);
541 if(!(is_inside_next || is_inside_prev)) return std::deque< double >();
545 _End_segments.push_back(std::make_pair(i_glob,i_next_glob));
546 addCrossing0(&Poly1[DIM*i_loc], &Poly1[DIM*i_next], i_glob, i_next_glob,
547 &Poly2[DIM*j3] , &Poly2[DIM*j4] , j3_glob,j4_glob);
551 _End_segments.push_back(std::make_pair(j1_glob,j2_glob));
552 _Status.insert(std::make_pair(i_next_glob,std::make_pair(i_glob, false)));
553 mi=_Status.find(j2_glob);
554 ((* mi).second).second= !((* mi).second).second;
558 _End_segments.push_back(std::make_pair(i_glob,i_prev_glob));
559 addCrossing0(&Poly1[DIM*i_loc], &Poly1[DIM*i_prev], i_glob, i_prev_glob,
560 &Poly2[DIM*j3] , &Poly2[DIM*j4] , j3_glob,j4_glob);
564 _End_segments.push_back(std::make_pair(j1_glob,j2_glob));
565 _Status.insert(std::make_pair(i_prev_glob,std::make_pair(i_glob, false)));
566 mi=_Status.find(j2_glob);
567 ((* mi).second).second= !((* mi).second).second;
570 else if(dotprod<DIM>(V12,V12) > _epsilon)//vertex i on a edge (j3,j4), not on (j1,j2)
572 crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j4], &Poly1[DIM*i_next],Vnext);
573 crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j4], &Poly1[DIM*i_prev],Vprev);
574 is_inside_next= dotprod<DIM>(Vnext,V12)<0;
575 is_inside_prev= dotprod<DIM>(Vprev,V12)<0;
577 if(!(is_inside_next || is_inside_prev)) return std::deque< double >();
581 _End_segments.push_back(std::make_pair(i_glob,i_next_glob));
582 addCrossing0(&Poly1[DIM*i_loc], &Poly1[DIM*i_next], i_glob, i_next_glob,
583 &Poly2[DIM*j1] , &Poly2[DIM*j2] , j1_glob,j2_glob);
587 _End_segments.push_back(std::make_pair(j3_glob,j4_glob));
588 _Status.insert(std::make_pair(i_next_glob,std::make_pair(i_glob, false)));
589 mi=_Status.find(j4_glob);
590 ((* mi).second).second= ! ((* mi).second).second;
594 _End_segments.push_back(std::make_pair(i_glob,i_prev_glob));
595 addCrossing0(&Poly1[DIM*i_loc], &Poly1[DIM*i_prev], i_glob, i_prev_glob,
596 &Poly2[DIM*j1] , &Poly2[DIM*j2] , j1_glob,j2_glob);
600 _End_segments.push_back(std::make_pair(j3_glob,j4_glob));
601 _Status.insert(std::make_pair(i_prev_glob,std::make_pair(i_glob, false)));
602 mi=_Status.find(j4_glob);
603 ((* mi).second).second= !((* mi).second).second;
606 else //vertices i, j1 and j3 share the same coordinates
608 crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j2], &Poly1[DIM*i_next],Vnext);
609 crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j2], &Poly1[DIM*i_prev],Vprev);
610 crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j4], &Poly1[DIM*i_next],V12);
611 crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j4], &Poly1[DIM*i_prev],V34);
613 double inside_next= dotprod<DIM>(Vnext,V12);
614 double inside_prev= dotprod<DIM>(Vprev,V34);
615 double inside_j2 = dotprod<DIM>(Vnext,Vprev);
616 double inside_j4 = dotprod<DIM>(V12,V34);
618 std::map<double, std::pair<int,int> > which_is_inside;
619 which_is_inside[inside_next] = std::make_pair(i_glob,i_next_glob);
620 which_is_inside[inside_prev] = std::make_pair(i_glob,i_prev_glob);
621 which_is_inside[inside_j2] = std::make_pair(j1_glob,j2_glob);
622 which_is_inside[inside_j4] = std::make_pair(j3_glob,j4_glob);
624 std::map<double, std::pair<int,int> >::iterator min = which_is_inside.begin();
625 std::map<double, std::pair<int,int> >::iterator minext = min;
627 std::map<double, std::pair<int,int> >::reverse_iterator max = which_is_inside.rbegin();
628 std::multimap< int, std::pair< int,bool> >::iterator j2_in_status = _Status.find(((*min).second).second);
629 std::multimap< int, std::pair< int,bool> >::iterator j4_in_status = _Status.find(((*minext).second).second);
631 if((*min).first < -_epsilon) //there is someone clearly inside
633 _End_segments.push_back( (*min).second );
634 _End_segments.push_back((* minext).second);
635 if(j2_in_status != _Status.end())
636 ((*j2_in_status).second).second = ! ((*j2_in_status).second).second;
637 if(j4_in_status != _Status.end())
638 ((*j4_in_status).second).second = ! ((*j4_in_status).second).second;
639 is_inside_next = ((*min).second).second == i_next_glob || ((*minext).second).second == i_next_glob;
640 is_inside_prev = ((*min).second).second == i_prev_glob || ((*minext).second).second == i_prev_glob;
643 if(fabs((*min).first) <= _epsilon) //nobody is clearly inside but two segments are superposed
645 if(fabs((*max).first) > _epsilon)
646 return std::deque< double >();
647 else //all four segments are superposed
649 _End_segments.push_back(std::make_pair(i_glob,i_next_glob));
650 _End_segments.push_back(std::make_pair(i_glob,i_prev_glob));
651 is_inside_next= true;
652 is_inside_prev= true;
655 else //there is nobody inside
656 return std::deque< double >();
658 _Status.insert(std::make_pair(i_prev_glob,std::make_pair(i_glob,is_inside_prev)));
659 _Status.insert(std::make_pair(i_next_glob,std::make_pair(i_glob,is_inside_next)));
665 std::cout << "Problem: nbprev= " << nb_prev << " ; i_glob = " << i_glob << std::endl;
666 throw Exception("intersectConvexPolygon: sequence of nodes does not describe a simple polygon (2)");
669 i_glob=(* mi1).second;//global index of vertex i
670 nb_prev = _Status.count(i_glob);
676 /**************************************************************************/
677 /* computes the convex hull of a polygon subP which is a sub polygon of P */
678 /* P is the array of coordinates, subP is a map containing initially the indices of a subpolygon of P */
679 /* in the end, subP contains only the elements belonging to the convex hull, and not_in_hull the others */
680 /**************************************************************************/
682 inline void PolygonAlgorithms<DIM>::convHull(const double *P, int N, double * normal,
683 std::map< int,int >& subP, std::map< int,int >& not_in_hull,
684 int& NsubP, const double epsilon)
688 std::map< int,int >::iterator mi_prev = subP.begin();
689 std::map< int,int >::iterator mi = mi_prev;
691 std::map< int,int >::iterator mi_next = mi;
693 double directframe=0.;
695 /* Check if the polygon subP is positively oriented */
696 std::map< int,int >::iterator mi1=mi;
697 while(mi1 != subP.end() && distance2<DIM>(&P[DIM*(*subP.begin()).second],&P[DIM*(*mi1).second])< epsilon)
699 std::map< int,int >::iterator mi2=mi1;
700 while(mi2 != subP.end() && fabs(directframe)<epsilon)
702 directframe =direct_frame<DIM>(&P[DIM* (*mi1).second],
703 &P[DIM* (*subP.begin()).second],
704 &P[DIM* (*mi2).second], normal);
707 if(directframe < 0) for(int idim=0; idim< DIM; idim++) normal[idim] *= -1;
709 /* Core of the algorithm */
710 while(mi_next != subP.end())
712 directframe = direct_frame<DIM>(&P[DIM* (*mi).second],
713 &P[DIM* (*mi_prev).second],
714 &P[DIM* (*mi_next).second], normal);
715 if(directframe > -epsilon){
722 not_in_hull.insert(*mi);
728 directframe = direct_frame<DIM>(&P[DIM*(*mi).second],
729 &P[DIM*(*mi_prev).second],
730 &P[DIM*(*subP.begin()).second], normal);
731 if(directframe < -epsilon)
733 not_in_hull.insert(*mi);
741 void PolygonAlgorithms<DIM>::convexDecomposition(const double * P, int N, double *normal, std::vector< int > subP, int NsubP,
742 std::vector< std::map< int,int > >& components, std::vector< int >& components_index,
743 int& Ncomp, int sign, const double epsilon)
746 std::map< int, int > hull;
747 std::map< int, int > not_in_hull;
748 std::map< int, int >::iterator mi, mj;
749 std::vector< int > reflex_region;
752 const double * xmax=&P[DIM*subP[0]];
753 /* checking an extremal point of subP */
754 for(i=0; i<NsubP; i++)
756 if(&P[DIM*subP[i]]> xmax)
759 xmax=&P[DIM*subP[i]];
762 /* renumbering of SubP elements for the convex hull*/
763 for(i=0; i<NsubP; i++) hull.insert(hull.end(),std::make_pair(i,subP[(i+i_xmax)%NsubP]));
764 /* compute the convex hull */
765 convHull(P, N, normal, hull, not_in_hull, NsubP,epsilon);
766 /* convex hull is the next component */
767 components.push_back(hull);
768 components_index.push_back(sign*NsubP);
770 /* searching for reflex regions */
771 for(mi=not_in_hull.begin(); mi!=not_in_hull.end(); mi++)
773 reflex_region.clear();
774 reflex_region.push_back(hull[(*mi).first-1]);
775 reflex_region.push_back( (*mi).second );
779 while((mj != not_in_hull.end()) && ((*mj).first == (*mi).first+1))
781 reflex_region.push_back((*mj).second);
786 reflex_region.push_back(hull[(*mi).first+1]);
788 convexDecomposition( P, N,normal, reflex_region, Nreflex, components, components_index, Ncomp, -sign, epsilon);
792 /**************************************************************************/
793 /* decomposes a non convex polygon P with N vertices contained in a plane */
794 /* into a sequence of convex polygons */
795 /* the input vectors 'components' and 'components_index' should be empty */
796 /* returns the number of convex components */
797 /* if P is composed of a single point, then an empty polygon is returned */
798 /**************************************************************************/
800 int PolygonAlgorithms<DIM>::convexDecomposition(const double * P, int N, std::vector< std::map< int,int > >& components,
801 std::vector< int >& components_index, const double epsilon)
804 std::vector< int > subP(N);
805 double normal[3]={0,0,0};
807 for(int i = 0; i<N; i++) subP[i]=i;
809 //Build the normal of polygon P
811 while(i1<N && distance2<DIM>(&P[0],&P[i1])< epsilon) i1++;
813 while(i2<N && fabs(dotprod<DIM>(normal,normal))<epsilon)
815 crossprod<DIM>(&P[i1], &P[0], &P[i2],normal);
819 convexDecomposition(P, N, normal, subP, N, components, components_index, Ncomp, 1, epsilon);