template< int DIM>
PolygonAlgorithms<DIM>::PolygonAlgorithms(double epsilon, double precision)//: (0)
{
- _Inter_size = 0;//To do remove des push_back and update of _Inter_size
_Is_in_intersection = false;
_Epsilon = epsilon;
_Precision = precision;
bool PolygonAlgorithms<DIM>::intersect_segment_segment(const double * A, const double * B, const double * C,
const double * D, const double * E, double * V)
{
- double AB[DIM], DC[DIM], AC[DIM];
+ double AB[DIM], DC[DIM], AC[DIM], det, t1, t2, inv_det;
+
+ /******* Initialisation of the linear system t1*AB+t2*DC=AC ***********/
for(int idim=0;idim<DIM;idim++)
{
AB[idim] = B[idim]-A[idim];//B-A
DC[idim] = C[idim]-D[idim];//C-D
AC[idim] = C[idim]-A[idim];//C-A
}
- double det = determinant(AB,DC);//determinant of the first two coefficients
- double t1, t2;
- int z_plane=0;
- if(abs(det) < _Epsilon)
+
+ /******* Resolution of the linear system t1*AB+t2*DC=AC ***********/
+ det = determinant(AB,DC);//determinant of the first two coordinates
+ if(abs(det) >_Epsilon)
+ {
+ inv_det = 1/det;
+ t1 = determinant(AC,DC)*inv_det;//solves the linear system t1*AB+t2*DC=AC
+ t2 = determinant(AB,AC)*inv_det;//solves the linear system t1*AB+t2*DC=AC
+ }
+ else
{
switch(DIM)
{
case 2:
return false;
- case 3:
- det = determinant(&AB[1],&DC[1]);//determinant of the two last coefficients
- if(abs(det) < _Epsilon) return false;//case of paralell segments
- else z_plane=1;
+ 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)
+ {
+ inv_det = 1/det;
+ t1=(AC[1]*DC[2]-AC[2]*DC[1])*inv_det;
+ t2=(AB[1]*AC[2]-AB[2]*AC[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)
+ {
+ inv_det = 1/det;
+ t1=(AC[0]*DC[2]-AC[2]*DC[0])*inv_det;
+ t2=(AB[0]*AC[2]-AB[2]*AC[0])*inv_det;
+ }
+ else return false;//case of paralell segments
+ }
}
}
- t1 = determinant(&AC[z_plane],&DC[z_plane])/det;//solves the linear system t1*AB+t2*DC=AC
- t2 = determinant(&AB[z_plane],&AC[z_plane])/det;//solves the linear system t1*AB+t2*DC=AC
-
if(t1>_Precision && t1<1-_Precision)
{
if( t2>_Precision && t2<1-_Precision)
else if( same_side > _Epsilon ) _Terminus= !_Is_in_intersection;//reflexion
else //separation of overlaping edges
{
- if(_Inter.size() ==0) _Terminus=true;
+ if(_Inter.empty() ) _Terminus=true;
else if(!_Is_in_intersection)
{
for(int idim=0;idim<DIM;idim++) V[idim]=A[idim];
crossprod<DIM>(A,E,B,_Vdouble);
if(dotprod<DIM>(_Vdouble,Vdoublebis) >=_Epsilon )//crossing
{
- if(_Inter.size()==0) _Terminus=true;
+ if(_Inter.empty()) _Terminus=true;
else if(!_Is_in_intersection)
{
for(int idim=0;idim<DIM;idim++) V[idim]=A[idim];
inline void PolygonAlgorithms<DIM>::add_crossing( double * ABCD, pair< int,int > i_i_next,
pair< int,int > j_j_next)
{
- if( _Inter.size()>0 )
+ if(!_Inter.empty() )
{
if(_End_segments[1] ==i_i_next)
{
i_next_glob = i_next+shift;
i_prev_glob = i_prev+shift;
//warning: sign is either 1 or -1;
- //To do test and remove from Convex_intersecor.cxx
- while(distance2<DIM>(&Poly1[DIM*i_loc],&Poly1[DIM*i_next])< _Epsilon && i_next != i_loc)
- i_next =(i_next+sign+N0)%N0;
- while(distance2<DIM>(&Poly1[DIM*i_loc],&Poly1[DIM*i_prev])< _Epsilon && i_prev != i_loc)
- i_prev =(i_prev+sign+N0)%N0;
+ //To do: test and remove from Convex_intersecor.cxx
+// while(distance2<DIM>(&Poly1[DIM*i_loc],&Poly1[DIM*i_next])< _Epsilon && i_next != i_loc)
+// i_next =(i_next+sign+N0)%N0;
+// while(distance2<DIM>(&Poly1[DIM*i_loc],&Poly1[DIM*i_prev])< _Epsilon && i_prev != i_loc)
+// i_prev =(i_prev+sign+N0)%N0;
}
/*******************************************************/
/* computes the vertices of the intersection of two COPLANAR */
i_prev, i_prev_glob, i_next, i_next_glob, nb_prev, sign, idim;
const double * Poly1, * Poly2;
bool four_neighbours=false;
-
- int minN1N2 = N1 < N2 ? N1 : N2;
-// _Inter.resize(2*minN1N2);
- _Terminus = minN1N2 < 3;
+ _Terminus = N1 < 3 || N2<3;
/* list of future events ordered according to their coordinates (x,y,z) (lexicographical order) */
multimap< const double *, int, VertexLess<DIM> > events;
if( _Is_in_intersection ) add_new_vertex(i_loc, i_glob, i_next_glob, i_prev_glob, Poly1);
add_crossings(&Poly1[DIM*i_loc], &Poly1[DIM*i_next], i_glob, i_next_glob,
&Poly2[DIM*j1] , &Poly2[DIM*j2] , j1_glob,j2_glob,
- &Poly2[DIM*j3] , &Poly2[DIM*j4] , j3_glob,j4_glob,
- &Poly1[DIM*i_prev]);
+ &Poly2[DIM*j3] , &Poly2[DIM*j4] , j3_glob,j4_glob, &Poly1[DIM*i_prev]);
break;
case 2 :
- if(_Inter.size()>0)
+ if(!_Inter.empty())
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 ://MN: à virer d'ici
- if(_Inter.size()==0){
+ case 0 ://To do: 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;
define_indices(i_loc,i_next,i_prev, Poly1,Poly2,
}
add_crossings(&Poly1[DIM*i_loc], &Poly1[DIM*i_next], i_glob, i_next_glob,
&Poly2[DIM*j1] , &Poly2[DIM*j2] , j1_glob,j2_glob,
- &Poly2[DIM*j3] , &Poly2[DIM*j4] , j3_glob,j4_glob,
- &Poly1[DIM*i_prev]);
+ &Poly2[DIM*j3] , &Poly2[DIM*j4] , j3_glob,j4_glob, &Poly1[DIM*i_prev]);
add_crossings(&Poly1[DIM*i_loc], &Poly1[DIM*i_prev], i_glob, i_prev_glob,
&Poly2[DIM*j1] , &Poly2[DIM*j2] , j1_glob,j2_glob,
- &Poly2[DIM*j3] , &Poly2[DIM*j4] , j3_glob,j4_glob,
- &Poly1[DIM*i_next]);
+ &Poly2[DIM*j3] , &Poly2[DIM*j4] , j3_glob,j4_glob, &Poly1[DIM*i_next]);
}
else //vertex on an edge
{
+++ /dev/null
-// Copyright (C) 2003 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
-// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
-//
-// This library is free software; you can redistribute it and/or
-// modify it under the terms of the GNU Lesser General Public
-// License as published by the Free Software Foundation; either
-// version 2.1 of the License.
-//
-// This library is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-// Lesser General Public License for more details.
-//
-// You should have received a copy of the GNU Lesser General Public
-// License along with this library; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-//
-// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-//
-//
-//
-// File : testUUnit.cxx
-// Module : MED
-
-#include "InterpolationUtils.hxx"
-using namespace std;
-
-int main ()
-{
- cerr<< endl;
- cerr<< " ************ Test 1 ************ "<< endl;
-
- const double Losange1[12]=
- {
- 1,0,0,
- 0,1,0,
- -1,0,0,
- 0,-1,0
- };
-
- const double Losange2[12]=
- {
- 2,0,0,
- 1,1,0,
- 0,0,0,
- 1,-1,0
- };
-
- deque< double > resultat_test1 = INTERP_UTILS::intersect_convex_polygons(Losange1,Losange2,4,4,1,1e-6);
-
- cerr<< " Test1: Résultat théorique" << endl;
- cerr<< 0.5 << " ,"<< -0.5 << " ," << 0 << endl;
- cerr<< 0 << " ,"<< 0 << " ,"<< 0 << endl;
- cerr<< 0.5 << " ,"<< 0.5 << " ,"<< 0 << endl;
- cerr<< 1 << " ," << 0 << " ," << 0 << endl;
-
- cerr<< " Test1: Résultat obtenu" << endl;
- for (int i=0; i<resultat_test1.size()/3; i++)
- {
- cerr << resultat_test1[3*i] << " ";
- cerr << resultat_test1[3*i+1] << " ";
- cerr << resultat_test1[3*i+2] << " ";
- cerr << endl ;
- }
-
- cerr<< " ************ Test 2 ************ "<< endl;
-
- const double Losange3[12]=
- {
- 2.5,0.5,0,
- 1.5,1.5,0,
- 0.5,0.5,0,
- 1.5,-0.5,0
- };
- deque< double > resultat_test2 = INTERP_UTILS::intersect_convex_polygons(Losange1,Losange3,4,4,1,1e-6);
-
- cerr<< " Test2: Résultat théorique" << endl;
-// cerr<< 0.5 << " ,"<< 0.5 << " ," << 0 << endl;
-// cerr<< 1 << " ,"<< 0 << " ," << 0 << endl;
-
- cerr<< " Test2: Résultat obtenu" << endl;
- for (int i=0; i<resultat_test2.size()/3; i++)
- {
- cerr << resultat_test2[3*i] << " ";
- cerr << resultat_test2[3*i+1] << " ";
- cerr << resultat_test2[3*i+2] << " ";
- cerr << endl ;
- }
-
- cerr<< " ************ Test 3 ************ "<< endl;
-
- const double Carre1[12]=
- {
- -1,-1,0,
- -1,1,0,
- 1,1,0,
- 1,-1,0
- };
- const double Carre2[12]=
- {
- 1,-0.25,0,
- 0,-0.25,0,
- 0,0.25,0,
- 1,0.25,0
- };
- deque< double > resultat_test3 = INTERP_UTILS::intersect_convex_polygons( Carre1, Carre2,4,4,1,1e-6);
-
- cerr<< " Test3: Résultat théorique" << endl;
- cerr<< 0 << " ,"<< 0.25 << " ," << 0 << endl;
- cerr<< 0 << " ,"<< -0.25 << " ," << 0 << endl;
- cerr<< 1 << " ,"<< -0.25 << " ,"<< 0 << endl;
- cerr<< 1 << " ,"<< 0.25 << " ," << 0 << endl;
-
- cerr<< " Test3: Résultat obtenu" << endl;
- for (int i=0; i<resultat_test3.size()/3; i++)
- {
- cerr << resultat_test3[3*i] << " ";
- cerr << resultat_test3[3*i+1] << " ";
- cerr << resultat_test3[3*i+2] << " ";
- cerr << endl ;
- }
- cerr<< " ***************** Test 4 ***************** " << endl;
-
- const double Losange4[12]=
- {
- 3,0,0,
- 2,1,0,
- 1,0,0,
- 2,-1,0
- };
- deque< double > resultat_test4 = INTERP_UTILS::intersect_convex_polygons( Losange1, Losange4,4,4,1,1e-6);
-
- cerr<< " Test4: Résultat théorique" << endl;
- cerr<< 1 << " ,"<< 0 << " ,"<< 0 << endl;
-
- cerr<< " Test4: Résultat obtenu" << endl;
- for (int i=0; i<resultat_test4.size()/3; i++)
- {
- cerr << resultat_test4[3*i] << " ";
- cerr << resultat_test4[3*i+1] << " ";
- cerr << resultat_test4[3*i+2] << " ";
- cerr << endl ;
- }
-
- cerr<< " ***************** Test 5 ***************** " << endl;
- deque< double > resultat_test5 = INTERP_UTILS::intersect_convex_polygons( Carre1, Carre1,4,4,1,1e-6);
-
- cerr<< " Test5: Résultat théorique" << endl;
- cerr<< -1 << " ,"<< 1 << " ," << 0 << endl;
- cerr<< -1 << " ,"<< -1 << " ," << 0 << endl;
- cerr<< 1 << " ,"<< -1 << " ," << 0 << endl;
- cerr<< 1 << " ,"<< 1 << " ," << 0 << endl;
-
- cerr<< " Test5: Résultat obtenu" << endl;
- for (int i=0; i<resultat_test5.size()/3; i++)
- {
- cerr << resultat_test5[3*i] << " ";
- cerr << resultat_test5[3*i+1] << " ";
- cerr << resultat_test5[3*i+2] << " ";
- cerr << endl ;
- }
-
- cerr<< " ***************** Test 6 ***************** " << endl;
-
- const double Losange5[12]=
- {
- 1.5,0,0,
- 0,1.5,0,
- -1.5,0,0,
- 0,-1.5,0
- };
- deque< double > resultat_test6 = INTERP_UTILS::intersect_convex_polygons( Carre1, Losange5,4,4,1,1e-6);
-
- cerr<< " Test6: Résultat théorique" << endl;
- cerr<< 1 << " ,"<< -0.5 << " ," << 0 << endl;
- cerr<< 0.5 << " ,"<< -1 << " ," << 0 << endl;
- cerr<< -0.5 << " ,"<< -1 << " ," << 0 << endl;
- cerr<< -1 << " ,"<< -0.5 << " ," << 0 << endl;
- cerr<< -1 << " ,"<< 0.5 << " ," << 0 << endl;
- cerr<< -0.5 << " ,"<< 1 << " ," << 0 << endl;
- cerr<< 0.5 << " ,"<< 1 << " ,"<< 0 << endl;
- cerr<< 1 << " ,"<< 0.5 << " ," << 0 << endl;
-
- cerr<< " Test6: Résultat obtenu" << endl;
- for (int i=0; i<resultat_test6.size()/3; i++)
- {
- cerr << resultat_test6[3*i] << " ";
- cerr << resultat_test6[3*i+1] << " ";
- cerr << resultat_test6[3*i+2] << " ";
- cerr << endl ;
- }
- cerr<< " ***************** Test 7 ***************** " << endl;
-
- deque< double > resultat_test7 = INTERP_UTILS::intersect_convex_polygons( Losange1, Carre1,4,4,1,1e-6);
-
- cerr<< " Test7: Résultat théorique" << endl;
- cerr<< 0 << " ,"<< -1 << " ," << 0 << endl;
- cerr<< -1 << " ,"<< 0 << " ," << 0 << endl;
- cerr<< 0 << " ,"<< 1 << " ,"<< 0 << endl;
- cerr<< 1 << " ,"<< 0 << " ," << 0 << endl;
-
- cerr<< " Test7: Résultat obtenu" << endl;
- for (int i=0; i<resultat_test7.size()/3; i++)
- {
- cerr << resultat_test7[3*i] << " ";
- cerr << resultat_test7[3*i+1] << " ";
- cerr << resultat_test7[3*i+2] << " ";
- cerr << endl ;
- }
- cerr<< " ************ Test 8 ************ "<< endl;
-
- const double Losange6[18]=
- {
- 2,0,0,
- 1,1,0,
- 0.5,0.5,0,
- 0,0,0,
- 0.5,-0.5,0,
- 1,-1,0
- };
- const double Losange7[15]=
- {
- 1,0,0,
- 0,1,0,
- -1,0,0,
- 0,-1,0,
- 0.5,-0.5,0
- };
-
- deque< double > resultat_test8 = INTERP_UTILS::intersect_convex_polygons(Losange6,Losange7,6,5,1,1e-6);
-
- cerr<< " Test8: Résultat théorique" << endl;
- cerr<< 0.5 << " ,"<< -0.5 << " ," << 0 << endl;
- cerr<< 0 << " ,"<< 0 << " ,"<< 0 << endl;
- cerr<< 0.5 << " ,"<< 0.5 << " ,"<< 0 << endl;
- cerr<< 1 << " ," << 0 << " ," << 0 << endl;
-
- cerr<< " Test8: Résultat obtenu" << endl;
- for (int i=0; i<resultat_test1.size()/3; i++)
- {
- cerr << resultat_test1[3*i] << " ";
- cerr << resultat_test1[3*i+1] << " ";
- cerr << resultat_test1[3*i+2] << " ";
- cerr << endl ;
- }
- cerr<< " ************ Test 9 ************ "<< endl;
- const double Carre3[15]=
- {
- -1,-1,0,
- -1,1,0,
- 0.5,1,0,
- 1,1,0,
- 1,-1,0
- };
- const double Carre4[12]=
- {
- -0.5,-1,0,
- -0.5,1,0,
- 1.5,1,0,
- 1.5,-1,0
- };
-
- deque< double > resultat_test9 = INTERP_UTILS::intersect_convex_polygons(Carre4,Carre3,4,5,1,1e-6);
-
- cerr<< " Test9: Résultat théorique" << endl;
- cerr<< -0.5 << " ,"<< -1 << " ," << 0 << endl;
- cerr<< -0.5 << " ,"<< 1 << " ," << 0 << endl;
- cerr<< 1 << " ,"<< 1 << " ," << 0 << endl;
- cerr<< 1 << " ,"<< -1 << " ," << 0 << endl;
-
- cerr<< " Test9: Résultat obtenu" << endl;
- for (int i=0; i<resultat_test9.size()/3; i++)
- {
- cerr << resultat_test9[3*i] << " ";
- cerr << resultat_test9[3*i+1] << " ";
- cerr << resultat_test9[3*i+2] << " ";
- cerr << endl ;
- }
-
- cerr<< " ************ Test 10 ************ "<< endl;
- const double Carre5[15]=
- {
- -1,-1,0,
- -1,1,0,
- 0,1,0,
- 1,1,0,
- 1,-1,0
- };
- const double Losange8[12]=
- {
- 0,1,0,
- 1,-1,0,
- 0,-1.5,0,
- -0.5,-1,0
- };
- deque< double > resultat_test10 = INTERP_UTILS::intersect_convex_polygons(Losange8,Carre5,4,5,1,1e-6);
- cerr<< " Test10: Résultat théorique" << endl;
- cerr<< 0 << " ,"<< 1 << " ," << 0 << endl;
- cerr<< -0.5 << " ,"<< -1 << " ," << 0 << endl;
- cerr<< 1 << " ,"<< -1 << " ," << 0 << endl;
-
- cerr<< " Test10: Résultat obtenu" << endl;
- for (int i=0; i<resultat_test10.size()/3; i++)
- {
- cerr << resultat_test10[3*i] << " ";
- cerr << resultat_test10[3*i+1] << " ";
- cerr << resultat_test10[3*i+2] << " ";
- cerr << endl ;
- }
-
- cerr<< " ************ Test 11 ************ "<< endl;
-
- const double Losange9[12]=
- {
- 0.5,0,0,
- 0,1,0,
- -1.5,0,0,
- 0,-1,0
- };
-
- deque< double > resultat_test11 = INTERP_UTILS::intersect_convex_polygons(Losange1,Losange9,4,4,1,1e-6);
-
- cerr<< " Test11: Résultat théorique" << endl;
- cerr<< 0 << " ,"<< -1 << " ," << 0 << endl;
- cerr<< -1 << " ,"<< 0 << " ,"<< 0 << endl;
- cerr<< 0 << " ,"<< 1 << " ,"<< 0 << endl;
- cerr<< 0.5 << " ," << 0 << " ," << 0 << endl;
-
- cerr<< " Test11: Résultat obtenu" << endl;
- for (int i=0; i<resultat_test11.size()/3; i++)
- {
- cerr << resultat_test11[3*i] << " ";
- cerr << resultat_test11[3*i+1] << " ";
- cerr << resultat_test11[3*i+2] << " ";
- cerr << endl ;
- }
-};
cerr<< " Test17: Résultat théorique" << endl;
cerr<< 0 << " ,"<< 0.5 << " ,"<< 0 << endl;
cerr<< 0 << " ,"<< -0.5 << " ," << 0 << endl;
- cerr<< "??" << " ," << 1/3 << " ," << 0 << endl;
+ cerr<< 2./3 << " ," << 1/3 << " ," << 0 << endl;
cerr<< " Test17: Résultat obtenu" << endl;
for (int i=0; i<resultat_test17.size()/3; i++)
cerr << resultat_test17[3*i+2] << " ";
cerr << endl ;
}
+ cerr<< " ************ Test 18 ************ "<< endl;
+ INTERP_UTILS::PolygonAlgorithms<3> P_18(1e-6,1e-6);
+
+
+ const double triangle3[9]=
+ {
+ 3.928371006592e-03, 0.000000000000e+00, -1.276720577142e-02 ,
+ 0.000000000000e+00, 0.000000000000e+00, -1.669557677802e-02 ,
+ 3.928371006592e-03, 0.000000000000e+00, -1.669557677802e-02 ,
+ };
+ const double triangle4[9]=
+ {
+ 0.000000000000e+00, 0.000000000000e+00, -1.767766952966e-02,
+ 4.419417382416e-03, 0.000000000000e+00, -1.325825214725e-02,
+ 0.000000000000e+00, 0.000000000000e+00, -1.325825214725e-02,
+ };
+
+ deque< double > resultat_test18 = P_18.intersect_convex_polygons(triangle3,triangle4,3,3);
+
+ cerr<< " Test18: Résultat théorique" << endl;
+ cerr<< 0.001 << " ,"<< 0. << " ,"<< -0.02 << endl;
+ cerr<< 0 << " ,"<< 0. << " ," << -1.669557677802e-02 << endl;
+ cerr<< 0.0034 << " ," << 0. << " ," << -0.0132 << endl;
+ cerr<< 0.0039 << " ," << 0. << " ," << -0.0132 << endl;
+ cerr<< 0.0039 << " ," << 0. << " ," << -0.014 << endl;
+
+ cerr<< " Test18: Résultat obtenu" << endl;
+ for (int i=0; i<resultat_test18.size()/3; i++)
+ {
+ cerr << resultat_test18[3*i] << " ";
+ cerr << resultat_test18[3*i+1] << " ";
+ cerr << resultat_test18[3*i+2] << " ";
+ cerr << endl ;
+ }
};
