#include <iostream>
+using namespace std;
+
namespace INTERP_KERNEL
{
template<int DIM>
/******* Resolution of the linear system t1*AB+t2*DC=AC ***********/
det = determinant(AB,DC);//determinant of the first two coordinates
- if(abs(det) >_Epsilon)
+ if(fabs(det) >_Epsilon)
{
inv_det = 1/det;
t1 = determinant(AC,DC)*inv_det;//solves the linear system t1*AB+t2*DC=AC
return false;
case 3://beware AB and CD may belong to a vertical plane
det = determinant(&AB[1],&DC[1]);//determinant of the last two coefficients
- if(abs(det) > _Epsilon)
+ if(fabs(det) > _Epsilon)
{
inv_det = 1/det;
t1=(AC[1]*DC[2]-AC[2]*DC[1])*inv_det;
else //beware AB and CD may belong to a plane y = constant
{
det = AB[0]*DC[2]-AB[2]*DC[0];
- if(abs(det) > _Epsilon)
+ if(fabs(det) > _Epsilon)
{
inv_det = 1/det;
t1=(AC[0]*DC[2]-AC[2]*DC[0])*inv_det;
return true;
}
}
- else if(abs(t1) <= _Precision)
+ else if(fabs(t1) <= _Precision)
{
if( t2>_Precision && t2<1-_Precision)//vertex on an edge
{
}
}
}
- else if(abs(t2-1) <= _Precision)//vertex on a vertex (A=D), first run
+ else if(fabs(t2-1) <= _Precision)//vertex on a vertex (A=D), first run
crossprod<DIM>(A,C,E,_Vdouble);//angle between vectors AC and AE (E=vertex preceding A)
- else if(abs(t2) <= _Precision)//vertex on a vertex (A=C), second run
+ else if(fabs(t2) <= _Precision)//vertex on a vertex (A=C), second run
{
double Vdoublebis[DIM];
crossprod<DIM>(A,B,D,Vdoublebis);
for(int idim=0;idim<DIM;idim++) V[idim]=A[idim];
return true;
}
- else if(abs(in_between)<=_Epsilon && dotprod<DIM>(Vdoublebis,Vdoublebis) > _Epsilon)
+ else if(fabs(in_between)<=_Epsilon && dotprod<DIM>(Vdoublebis,Vdoublebis) > _Epsilon)
//ie _Vdouble=0, separation of overlaping edges at a double point
{
crossprod<DIM>(A,E,B,_Vdouble);
_Terminus = N1 < 3 || N2<3;
/* list of future events ordered according to their coordinates (x,y,z) (lexicographical order) */
- std::multimap< const double *, int, VertexLess<DIM> > events;
- typename std::multimap< const double *, int, VertexLess<DIM> >::iterator mi1,mi2;
+ std::multimap< const double *, int, VertexLess<DIM> > mmap_events;
+ typename std::list< pair< const double *, int > >::iterator mi1,mi2;
std::multimap< int, std::pair< int,bool> >::iterator mi;
/********** Initalisation of events with P1 and P2 vertices ************/
for(i_loc=0;i_loc<N1;i_loc++)
- events.insert(std::make_pair(&P_1[DIM*i_loc],i_loc));
+ mmap_events.insert(std::make_pair(&P_1[DIM*i_loc],i_loc));
for(i_loc=0;i_loc<N2;i_loc++)
- events.insert(std::make_pair(&P_2[DIM*i_loc],i_loc+N1));
-
- if(!_Terminus)
+ mmap_events.insert(std::make_pair(&P_2[DIM*i_loc],i_loc+N1));
+
+ std::list< pair< const double *, int > > events(mmap_events.begin(),mmap_events.end());
+
+ if(!_Terminus)
{
/******** Treatment of the first vertex ********/
mi1=events.begin();
_Status.insert(std::make_pair(i_next_glob,std::make_pair(i_glob, false)));
_Status.insert(std::make_pair(i_prev_glob,std::make_pair(i_glob, false)));
mi1++;
-
+ //std::cout<< "nb_prev= "<< 0 << " i_glob= " << i_glob << std::endl;
+
/******* Loop until the second polygon is reached *******/
while( !four_neighbours)
{
i_glob=(* mi1).second;//global index of vertex i
nb_prev = _Status.count(i_glob);//counts the number of segments ending at i
+
+ //std::cout<< "nb_prev= "<< nb_prev << " i_glob= " << i_glob << std::endl;
switch (nb_prev)
{
case 1 :
mi=_Status.find(i_glob);// pointer to the segment ending at i
i_prev_glob = ((*mi).second).first;//starting point of the segment ending at i
- i_next= (i_prev_glob - i_glob > 0) == (abs(i_prev_glob - i_glob) == 1) ? i_glob - 1 : i_glob + 1;
+ i_next= (i_prev_glob - i_glob > 0) == (fabs(i_prev_glob - i_glob) == 1) ? i_glob - 1 : i_glob + 1;
if(i_glob < N1) i_next_glob = (i_next +N1)%N1;
else i_next_glob = (i_next-N1+N2)%N2 + N1;
_Status.erase(mi);
case 0 :
if( (i_glob < N1) != which_start)
{
- four_neighbours = true;
- /* detection of double points */
- const double * Pi1=(* mi1).first;
- mi2=mi1++;
- const double * Pi2=(* mi2).first;
- int i_double= (*mi2).second;
- if(distance2<DIM>(Pi1, Pi2) < _Epsilon && _Status.count(i_double)==1)
+ mi2=mi1;
+ mi2++;
+ /* detection of double points */
+ if(distance2<DIM>((* mi1).first, (*mi2).first) > _Epsilon)
+ four_neighbours = true;
+ else /* Rare pothological case: */
{
- i_glob=i_double;
- nb_prev=1;
+ //std::cout<<"coucou" << std::endl;
+ const pair< const double *, int > next_pt= *mi2;
events.erase(mi2);
- }
- mi1--;
+ mi1=events.insert(mi1,next_pt);
+ }
}
break;
default:
/******** Loop until a terminal point or crossing is reached ************/
while( !_Terminus)
{
+ //std::cout<< "nb_prev= "<< nb_prev<< " nb_inter= " << _Inter.size()/DIM << std::endl;
switch (nb_prev)
{
case 1 :
mi=_Status.find(i_glob);// pointer to the segment ending at i
i_prev_glob = ((*mi).second).first;//starting point of the segment ending at i
- sign = (i_prev_glob - i_glob > 0) == (abs(i_prev_glob - i_glob) == 1) ? - 1 : + 1;
+ sign = (i_prev_glob - i_glob > 0) == (fabs(i_prev_glob - i_glob) == 1) ? - 1 : + 1;
i_next_glob = i_glob+sign;
_Is_in_intersection = ((*mi).second).second;//boolean that tells if i is in the intersection
_Status.erase(mi);
if(i_glob < N1) for(idim=0;idim<DIM;idim++) _Inter.push_back(P_1[DIM*i_glob+idim]);
else for(idim=0;idim<DIM;idim++) _Inter.push_back(P_2[DIM*(i_glob-N1)+idim]);
return _Inter;
- case 0 ://To do: remove this case from here
+ case 0 ://To do if possible : remove this case from here
if(_Inter.empty() && (i_glob < N1) != which_start){
i_next_glob = i_glob+1;
i_prev_glob = i_glob-1;
&Poly2[DIM*j3], &Poly2[DIM*j4],V12, V34);
_Is_in_intersection=( inside < _Epsilon ); // <= epsilon or 0 ?
- if(abs(inside) > _Epsilon)//vertex clearly inside or outside
- {
+ if(fabs(inside) > _Epsilon)//vertex clearly inside or outside
+ {
+ //std::cout<<"coucou1" << std::endl;
if( _Is_in_intersection)
{
for(int idim=0;idim<DIM;idim++) _Inter.push_back(Poly1[DIM*i_loc+idim]);
}
else //vertex on an edge
{
- bool inside_next, inside_prev;
+ //std::cout<<"coucou2" << std::endl;
+ bool is_inside_next, is_inside_prev;
double Vnext[DIM], Vprev[DIM];
for(idim=0;idim<DIM;idim++) _Inter.push_back(Poly1[DIM*i_loc+idim]);
{
crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j2], &Poly1[DIM*i_next],Vnext);
crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j2], &Poly1[DIM*i_prev],Vprev);
- inside_next= (dotprod<DIM>(Vnext,V34)<0);
- inside_prev= (dotprod<DIM>(Vprev,V34)<0);
+ is_inside_next= (dotprod<DIM>(Vnext,V34)<0);
+ is_inside_prev= (dotprod<DIM>(Vprev,V34)<0);
- if(!(inside_next || inside_prev)) return std::deque< double >();
+ if(!(is_inside_next || is_inside_prev)) return std::deque< double >();
- if(inside_next)
+ if(is_inside_next)
{
_End_segments.push_back(std::make_pair(i_glob,i_next_glob));
add_crossing0(&Poly1[DIM*i_loc], &Poly1[DIM*i_next], i_glob, i_next_glob,
mi=_Status.find(j2_glob);
((* mi).second).second= !((* mi).second).second;
}
- if(inside_prev)
+ if(is_inside_prev)
{
_End_segments.push_back(std::make_pair(i_glob,i_prev_glob));
add_crossing0(&Poly1[DIM*i_loc], &Poly1[DIM*i_prev], i_glob, i_prev_glob,
{
crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j4], &Poly1[DIM*i_next],Vnext);
crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j4], &Poly1[DIM*i_prev],Vprev);
- inside_next= dotprod<DIM>(Vnext,V12)<0;
- inside_prev= dotprod<DIM>(Vprev,V12)<0;
+ is_inside_next= dotprod<DIM>(Vnext,V12)<0;
+ is_inside_prev= dotprod<DIM>(Vprev,V12)<0;
- if(!(inside_next || inside_prev)) return std::deque< double >();
+ if(!(is_inside_next || is_inside_prev)) return std::deque< double >();
- if(inside_next)
+ if(is_inside_next)
{
_End_segments.push_back(std::make_pair(i_glob,i_next_glob));
add_crossing0(&Poly1[DIM*i_loc], &Poly1[DIM*i_next], i_glob, i_next_glob,
mi=_Status.find(j4_glob);
((* mi).second).second= ! ((* mi).second).second;
}
- if(inside_prev)
+ if(is_inside_prev)
{
_End_segments.push_back(std::make_pair(i_glob,i_prev_glob));
add_crossing0(&Poly1[DIM*i_loc], &Poly1[DIM*i_prev], i_glob, i_prev_glob,
}
else //vertices i, j1 and j3 share the same coordinates
{
- crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j2], &Poly1[DIM*i_next],Vnext);
+ crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j2], &Poly1[DIM*i_next],Vnext);
crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j2], &Poly1[DIM*i_prev],Vprev);
crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j4], &Poly1[DIM*i_next],V12);
crossprod<DIM>(&Poly1[DIM*i_loc], &Poly2[DIM*j4], &Poly1[DIM*i_prev],V34);
-
- inside_next= dotprod<DIM>(Vnext,V12)<=_Epsilon;
- inside_prev= dotprod<DIM>(Vprev,V34)<=_Epsilon;
- bool inside_j2= dotprod<DIM>(Vnext,Vprev)<-_Epsilon;
-
- if(!( inside_next || inside_prev || inside_j2 )) return std::deque< double >();
- if(inside_next) _End_segments.push_back(std::make_pair(i_glob,i_next_glob));
- if(inside_prev) _End_segments.push_back(std::make_pair(i_glob,i_prev_glob));
- if(inside_j2)
- {
- _End_segments.push_back(std::make_pair(j1_glob,j2_glob));
- mi=_Status.find(j2_glob);
- ((* mi).second).second= !((* mi).second).second;
- }
- if(!( (inside_next && (inside_prev || inside_j2)) || (inside_prev && inside_j2)))
- {
- _End_segments.push_back(std::make_pair(j3_glob,j4_glob));
- mi=_Status.find(j4_glob);
- ((* mi).second).second= !((* mi).second).second;
- }
- _Status.insert(std::make_pair(i_prev_glob,std::make_pair(i_glob,inside_prev)));
- _Status.insert(std::make_pair(i_next_glob,std::make_pair(i_glob,inside_next)));
+
+ double inside_next= dotprod<DIM>(Vnext,V12);
+ double inside_prev= dotprod<DIM>(Vprev,V34);
+ double inside_j2 = dotprod<DIM>(Vnext,Vprev);
+ double inside_j4 = dotprod<DIM>(V12,V34);
+
+ std::map<double, pair<int,int> > which_is_inside;
+ which_is_inside[inside_next] = std:: make_pair(i_glob,i_next_glob);
+ which_is_inside[inside_prev] = std:: make_pair(i_glob,i_prev_glob);
+ which_is_inside[inside_j2] = std::make_pair(j1_glob,j2_glob);
+ which_is_inside[inside_j4] = std::make_pair(j3_glob,j4_glob);
+
+ std::map<double, pair<int,int> >::iterator min = which_is_inside.begin();
+ std::map<double, pair<int,int> >::iterator minext = min;
+ minext++;
+ std::map<double, pair<int,int> >::reverse_iterator max = which_is_inside.rbegin();
+ std::multimap< int, std::pair< int,bool> >::iterator j2_in_status = _Status.find(((*min).second).second);
+ std::multimap< int, std::pair< int,bool> >::iterator j4_in_status = _Status.find(((*minext).second).second);
+
+ if((*min).first < -_Epsilon) //there is someone clearly inside
+ {
+ _End_segments.push_back( (*min).second );
+ _End_segments.push_back((* minext).second);
+ if(j2_in_status != _Status.end())
+ ((*j2_in_status).second).second = ! ((*j2_in_status).second).second;
+ if(j4_in_status != _Status.end())
+ ((*j4_in_status).second).second = ! ((*j4_in_status).second).second;
+ is_inside_next = ((*min).second).second == i_next_glob || ((*minext).second).second == i_next_glob;
+ is_inside_prev = ((*min).second).second == i_prev_glob || ((*minext).second).second == i_prev_glob;
+ }
+ else
+ if(fabs((*min).first) <= _Epsilon) //nobody is clearly inside but two segments are superposed
+ {
+ if(fabs((*max).first) > _Epsilon)
+ return std::deque< double >();
+ else //all four segments are superposed
+ {
+ _End_segments.push_back(std::make_pair(i_glob,i_next_glob));
+ _End_segments.push_back(std::make_pair(i_glob,i_prev_glob));
+ is_inside_next= true;
+ is_inside_prev= true;
+ }
+ }
+ else //there is nobody inside
+ return std::deque< double >();
+
+ _Status.insert(std::make_pair(i_prev_glob,std::make_pair(i_glob,is_inside_prev)));
+ _Status.insert(std::make_pair(i_next_glob,std::make_pair(i_glob,is_inside_next)));
}
}
}
if(NsubP>3)
{
std::map< int,int >::iterator mi_prev = subP.begin();
- std::map< int,int >::iterator mi = mi_prev++;
- std::map< int,int >::iterator mi_next = mi++;
+ std::map< int,int >::iterator mi = mi_prev;
+ mi++;
+ std::map< int,int >::iterator mi_next = mi;
+ mi_next++;
double directframe=0.;
/* Check if the polygon subP is positively oriented */
while(mi1 != subP.end() && distance2<DIM>(&P[DIM*(*subP.begin()).second],&P[DIM*(*mi1).second])< epsilon)
mi1++;
std::map< int,int >::iterator mi2=mi1;
- while(mi2 != subP.end() && abs(directframe)<epsilon)
+ while(mi2 != subP.end() && fabs(directframe)<epsilon)
{
directframe =direct_frame<DIM>(&P[DIM* (*mi1).second],
&P[DIM* (*subP.begin()).second],
&P[DIM* (*mi_prev).second],
&P[DIM* (*mi_next).second], normal);
if(directframe > -epsilon){
- mi = mi++;
- mi_prev=mi_prev++;
- mi_next=mi_next++;
+ mi ++;
+ mi_prev++;
+ mi_next++;
}
else
{
int i1=1;
while(i1<N && distance2<DIM>(&P[0],&P[i1])< epsilon) i1++;
int i2=i1+1;
- while(i2<N && abs(dotprod<DIM>(normal,normal))<epsilon)
+ while(i2<N && fabs(dotprod<DIM>(normal,normal))<epsilon)
{
crossprod<DIM>(&P[i1], &P[0], &P[i2],normal);
i2++;
