src/IntPoly/IntPoly_ShapeSection.cxx

   1 // Created on: 1995-08-01
   2 // Created by: Stagiaire Alain JOURDAIN
   3 // Copyright (c) 1995-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and / or modify it
   9 // under the terms of the GNU Lesser General Public version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 #include <IntPoly_ShapeSection.ixx>
  18
  19 #include <IntPoly_SequenceOfSequenceOfPnt.hxx>
  20 #include <IntPoly_IndexedMapOfPnt.hxx>
  21 #include <IntPoly_PntHasher.hxx>
  22 #include <Precision.hxx>
  23 #include <TopExp_Explorer.hxx>
  24 #include <TopoDS_Shape.hxx>
  25 #include <TopoDS_Edge.hxx>
  26 #include <TColStd_Array1OfReal.hxx>
  27 #include <TColStd_Array1OfInteger.hxx>
  28 #include <BRep_Tool.hxx>
  29 #include <BRep_Builder.hxx>
  30 #include <TopoDS.hxx>
  31 #include <gp_Pnt.hxx>
  32 #include <gp_Dir.hxx>
  33 #include <gp_XYZ.hxx>
  34 #include <Poly_Triangulation.hxx>
  35 #include <Poly_Polygon3D.hxx>
  36
  37
  38 //=======================================================================
  39 //function : IntPoly_ShapeSection
  40 //purpose  :
  41 //=======================================================================
  42
  43 IntPoly_ShapeSection::IntPoly_ShapeSection()
  44 {}
  45
  46
  47 //=======================================================================
  48 //function : IntPoly_ShapeSection
  49 //purpose  :
  50 //=======================================================================
  51
  52 IntPoly_ShapeSection::IntPoly_ShapeSection(const TopoDS_Shape& S1,
  53                                            const TopoDS_Shape& S2)
  54 {
  55   myShape1 = S1;
  56   myShape2 = S2;
  57   Section();
  58 }
  59
  60
  61 //=======================================================================
  62 //function : Section
  63 //purpose  :
  64 //=======================================================================
  65
  66 void IntPoly_ShapeSection::Section()
  67 {
  68   Explore();
  69   Standard_Integer NbLinks = myMapBegPoints.Extent();
  70   Standard_Integer Result,i = 1;
  71   Standard_Real Big = Precision::Infinite();
  72   gp_Pnt BegPoint,EndPoint,OutPoint;
  73
  74   while (i <= NbLinks) {
  75     if (!(Precision::IsInfinite((myMapBegPoints.FindKey(i)).X()))) {
  76       BegPoint = myMapBegPoints.FindKey(i);
  77       EndPoint = myMapEndPoints.FindKey(i);
  78       myCpt++;
  79       myMapBegPoints.Substitute(i,gp_Pnt(Big,myCpt,myCpt));
  80       myMapEndPoints.Substitute(i,gp_Pnt(Big,myCpt,myCpt));
  81       Result = Concat(BegPoint,EndPoint,OutPoint);
  82       if (Result == 2)
  83         ForwConstruction(OutPoint);
  84       else
  85         if (Result == 1)
  86           PrevConstruction(OutPoint);
  87         else {
  88           ForwConstruction(EndPoint);
  89           PrevConstruction(BegPoint);
  90         }
  91     }
  92     i++;
  93   }
  94
  95   i = 1;
  96   while (i < mySection.Length()) {
  97     ConcatSection(mySection.ChangeValue(i),mySection.Length(),i+1);
  98     i++;
  99   }
 100   myNbEdges = mySection.Length();
 101 }
 102
 103
 104
 105 //=======================================================================
 106 //function : Explore
 107 //purpose  :
 108 //=======================================================================
 109
 110 inline Standard_Boolean SAMESIGN(Standard_Real a,
 111                                  Standard_Real b,
 112                                  Standard_Real c)
 113 {
 114   return ((a > 0.0 && b > 0.0 && c > 0.0) ||
 115           (a < 0.0 && b < 0.0 && c < 0.0) ||
 116           (a == 0.0 && b == 0.0 && c == 0.0));
 117 }
 118
 119 inline void MINMAX(Standard_Real a,Standard_Real b,Standard_Real c,
 120                    Standard_Real &m,Standard_Real &M)
 121 {
 122   m = M = a;
 123   if (b < c) {
 124     if (b < m) m = b;
 125     if (c > M) M = c;
 126   }
 127   else {
 128     if (c < m) m = c;
 129     if (b > M) M = b;
 130   }
 131 }
 132
 133 void IntPoly_ShapeSection::Explore()
 134 {
 135   TopExp_Explorer ex1,ex2;
 136   Standard_Integer i,i1,i2,i3,j,j1,j2,j3;
 137   Standard_Integer NbTrian1,NbFaces1,NbTriATotal;
 138   Standard_Integer NbTrian2,NbFaces2,NbTriBTotal;
 139   Standard_Real Big = Precision::Infinite();
 140   Standard_Real min,minfx,minfy,minfz,max,maxfx,maxfy,maxfz;
 141   TopLoc_Location Loc1,Loc2;
 142   Handle(Poly_Triangulation) Tr1,Tr2;
 143   myCpt = 0;
 144
 145   //-----------------------------------------------------------------
 146   //--  Traitement pour recuperer le nombre de faces et le nombre
 147   //--  de triangles des shapes 1 et 2
 148   //-----------------------------------------------------------------
 149   NbFaces1 = NbTriATotal = 0;
 150   for (ex1.Init(myShape1,TopAbs_FACE);ex1.More();ex1.Next()) {
 151     Tr1 = (BRep_Tool::Triangulation(TopoDS::Face(ex1.Current()),Loc1));
 152     if (!Tr1.IsNull()) {
 153       NbTriATotal += Tr1->NbTriangles();
 154       NbFaces1++;
 155     }
 156   }
 157   NbFaces2 = NbTriBTotal = 0;
 158   for (ex2.Init(myShape2,TopAbs_FACE);ex2.More();ex2.Next()) {
 159     Tr2 = (BRep_Tool::Triangulation(TopoDS::Face(ex2.Current()),Loc2));
 160     if (!Tr2.IsNull()) {
 161       NbTriBTotal += Tr2->NbTriangles();
 162       NbFaces2++;
 163     }
 164   }
 165
 166   if (NbTriATotal < NbTriBTotal) {
 167     TopoDS_Shape S = myShape2;
 168     myShape2 = myShape1;
 169     myShape1 = S;
 170     Standard_Integer Nb = NbTriBTotal;
 171     NbTriBTotal = NbTriATotal;
 172     NbTriATotal = Nb;
 173     Nb       = NbFaces2;
 174     NbFaces2 = NbFaces1;
 175     NbFaces1 = Nb;
 176   }
 177
 178   //-----------------------------------------------------------------
 179   //--  Tableaux des sommets des NbTriATotal triangles du shape 1
 180   //-----------------------------------------------------------------
 181   TColgp_Array1OfPnt TA1(1,NbTriATotal);
 182   TColgp_Array1OfPnt TA2(1,NbTriATotal);
 183   TColgp_Array1OfPnt TA3(1,NbTriATotal);
 184   //-----------------------------------------------------------------
 185   //-- Tableaux des MinMax des NbTriATotal triangles du shape 1
 186   //-----------------------------------------------------------------
 187   TColStd_Array1OfReal MinAx(1,NbTriATotal);
 188   TColStd_Array1OfReal MinAy(1,NbTriATotal);
 189   TColStd_Array1OfReal MinAz(1,NbTriATotal);
 190   TColStd_Array1OfReal MaxAx(1,NbTriATotal);
 191   TColStd_Array1OfReal MaxAy(1,NbTriATotal);
 192   TColStd_Array1OfReal MaxAz(1,NbTriATotal);
 193   //-----------------------------------------------------------------
 194   //-- Tableaux des MinMax des NbFaces1 faces du shape 1
 195   //-----------------------------------------------------------------
 196   TColStd_Array1OfReal MinFAx(1,NbFaces1);
 197   TColStd_Array1OfReal MinFAy(1,NbFaces1);
 198   TColStd_Array1OfReal MinFAz(1,NbFaces1);
 199   TColStd_Array1OfReal MaxFAx(1,NbFaces1);
 200   TColStd_Array1OfReal MaxFAy(1,NbFaces1);
 201   TColStd_Array1OfReal MaxFAz(1,NbFaces1);
 202   //-----------------------------------------------------------------
 203   //-- Tableau contenant le nombre de triangles des faces du shape 1
 204   //-----------------------------------------------------------------
 205   TColStd_Array1OfInteger TNbTr1(1,NbFaces1);
 206
 207   //-----------------------------------------------------------------
 208   //--  Preparation : Mise a jour des tableaux et calcul des MinMax
 209   //-----------------------------------------------------------------
 210   NbTriATotal = NbFaces1 = 0;
 211   for (ex1.Init(myShape1,TopAbs_FACE);ex1.More();ex1.Next()) {
 212     Tr1 = (BRep_Tool::Triangulation(TopoDS::Face(ex1.Current()),Loc1));
 213     if (!Tr1.IsNull()) {
 214       NbTrian1 = Tr1->NbTriangles();
 215       NbFaces1++;
 216       TNbTr1.SetValue(NbFaces1,NbTrian1);
 217       const Poly_Array1OfTriangle& TabTrian1 = Tr1->Triangles();
 218       const TColgp_Array1OfPnt&    TabNodes1 = Tr1->Nodes();
 219       minfx = minfy = minfz = -Big;
 220       maxfx = maxfy = maxfz =  Big;
 221       for (i = 1;i <= NbTrian1;i++) {
 222         TabTrian1(i).Get(i1,i2,i3);
 223         NbTriATotal++;
 224         TA1.SetValue(NbTriATotal,TabNodes1(i1).Transformed(Loc1.Transformation()));
 225         TA2.SetValue(NbTriATotal,TabNodes1(i2).Transformed(Loc1.Transformation()));
 226         TA3.SetValue(NbTriATotal,TabNodes1(i3).Transformed(Loc1.Transformation()));
 227         MINMAX(TA1.Value(NbTriATotal).X(),
 228                TA2.Value(NbTriATotal).X(),
 229                TA3.Value(NbTriATotal).X(),min,max);
 230         MinAx.SetValue(NbTriATotal,min);
 231         MaxAx.SetValue(NbTriATotal,max);
 232         if (min < minfx) minfx = min;
 233         if (max > maxfx) maxfx = max;
 234         MINMAX(TA1.Value(NbTriATotal).Y(),
 235                TA2.Value(NbTriATotal).Y(),
 236                TA3.Value(NbTriATotal).Y(),min,max);
 237         MinAy.SetValue(NbTriATotal,min);
 238         MaxAy.SetValue(NbTriATotal,max);
 239         if (min < minfy) minfy = min;
 240         if (max > maxfy) maxfy = max;
 241         MINMAX(TA1.Value(NbTriATotal).Z(),
 242                TA2.Value(NbTriATotal).Z(),
 243                TA3.Value(NbTriATotal).Z(),min,max);
 244         MinAz.SetValue(NbTriATotal,min);
 245         MaxAz.SetValue(NbTriATotal,max);
 246         if (min < minfz) minfz = min;
 247         if (max > maxfz) maxfz = max;
 248       }
 249       MinFAx.SetValue(NbFaces1,minfx);
 250       MaxFAx.SetValue(NbFaces1,maxfx);
 251       MinFAy.SetValue(NbFaces1,minfy);
 252       MaxFAy.SetValue(NbFaces1,maxfy);
 253       MinFAz.SetValue(NbFaces1,minfz);
 254       MaxFAz.SetValue(NbFaces1,maxfz);
 255     }
 256   }
 257
 258   //-----------------------------------------------------------------
 259   //--  Tableaux des sommets des NbTriBTotal triangles du shape 2
 260   //-----------------------------------------------------------------
 261   TColgp_Array1OfPnt TB1(1,NbTriBTotal);
 262   TColgp_Array1OfPnt TB2(1,NbTriBTotal);
 263   TColgp_Array1OfPnt TB3(1,NbTriBTotal);
 264   //-----------------------------------------------------------------
 265   //-- Tableaux des MinMax des NbTriBTotal triangles du shape 2
 266   //-----------------------------------------------------------------
 267   TColStd_Array1OfReal MinBx(1,NbTriBTotal);
 268   TColStd_Array1OfReal MinBy(1,NbTriBTotal);
 269   TColStd_Array1OfReal MinBz(1,NbTriBTotal);
 270   TColStd_Array1OfReal MaxBx(1,NbTriBTotal);
 271   TColStd_Array1OfReal MaxBy(1,NbTriBTotal);
 272   TColStd_Array1OfReal MaxBz(1,NbTriBTotal);
 273   //-----------------------------------------------------------------
 274   //-- Tableaux des MinMax des NbFaces2 faces du shape 2
 275   //-----------------------------------------------------------------
 276   TColStd_Array1OfReal MinFBx(1,NbFaces2);
 277   TColStd_Array1OfReal MinFBy(1,NbFaces2);
 278   TColStd_Array1OfReal MinFBz(1,NbFaces2);
 279   TColStd_Array1OfReal MaxFBx(1,NbFaces2);
 280   TColStd_Array1OfReal MaxFBy(1,NbFaces2);
 281   TColStd_Array1OfReal MaxFBz(1,NbFaces2);
 282   //-----------------------------------------------------------------
 283   //-- Tableau contenant le nombre de triangles des faces du shape 2
 284   //-----------------------------------------------------------------
 285   TColStd_Array1OfInteger TNbTr2(1,NbFaces2);
 286
 287   //-----------------------------------------------------------------
 288   //--  Preparation : Mise a jour des tableaux et calcul des MinMax
 289   //-----------------------------------------------------------------
 290   NbTriBTotal = NbFaces2 = 0;
 291   for (ex2.Init(myShape2,TopAbs_FACE);ex2.More();ex2.Next()) {
 292     Tr2 = (BRep_Tool::Triangulation(TopoDS::Face(ex2.Current()),Loc2));
 293     if (!Tr2.IsNull()) {
 294       NbTrian2 = Tr2->NbTriangles();
 295       NbFaces2++;
 296       TNbTr2.SetValue(NbFaces2,NbTrian2);
 297       const Poly_Array1OfTriangle& TabTrian2 = Tr2->Triangles();
 298       const TColgp_Array1OfPnt&    TabNodes2 = Tr2->Nodes();
 299       minfx = minfy = minfz = -Big;
 300       maxfx = maxfy = maxfz =  Big;
 301       for (j = 1;j <= NbTrian2;j++) {
 302         TabTrian2(j).Get(j1,j2,j3);
 303         NbTriBTotal++;
 304         TB1.SetValue(NbTriBTotal,TabNodes2(j1).Transformed(Loc2.Transformation()));
 305         TB2.SetValue(NbTriBTotal,TabNodes2(j2).Transformed(Loc2.Transformation()));
 306         TB3.SetValue(NbTriBTotal,TabNodes2(j3).Transformed(Loc2.Transformation()));
 307         MINMAX(TB1.Value(NbTriBTotal).X(),
 308                TB2.Value(NbTriBTotal).X(),
 309                TB3.Value(NbTriBTotal).X(),min,max);
 310         MinBx.SetValue(NbTriBTotal,min);
 311         MaxBx.SetValue(NbTriBTotal,max);
 312         if (min < minfx) minfx = min;
 313         if (max > maxfx) maxfx = max;
 314         MINMAX(TB1.Value(NbTriBTotal).Y(),
 315                TB2.Value(NbTriBTotal).Y(),
 316                TB3.Value(NbTriBTotal).Y(),min,max);
 317         MinBy.SetValue(NbTriBTotal,min);
 318         MaxBy.SetValue(NbTriBTotal,max);
 319         if (min < minfy) minfy = min;
 320         if (max > maxfy) maxfy = max;
 321         MINMAX(TB1.Value(NbTriBTotal).Z(),
 322                TB2.Value(NbTriBTotal).Z(),
 323                TB3.Value(NbTriBTotal).Z(),min,max);
 324         MinBz.SetValue(NbTriBTotal,min);
 325         MaxBz.SetValue(NbTriBTotal,max);
 326         if (min < minfz) minfz = min;
 327         if (max > maxfz) maxfz = max;
 328       }
 329       MinFBx.SetValue(NbFaces2,minfx);
 330       MaxFBx.SetValue(NbFaces2,maxfx);
 331       MinFBy.SetValue(NbFaces2,minfy);
 332       MaxFBy.SetValue(NbFaces2,maxfy);
 333       MinFBz.SetValue(NbFaces2,minfz);
 334       MaxFBz.SetValue(NbFaces2,maxfz);
 335     }
 336   }
 337
 338   Standard_Integer ii,jj,fi,fj;
 339   i = j = 0;
 340   for (fi = 1;fi <= NbFaces1;fi++) {
 341     NbTrian1 = TNbTr1.Value(fi);
 342     i += NbTrian1;
 343     j = 0;
 344     for (fj = 1;fj <= NbFaces2;fj++) {
 345       NbTrian2 = TNbTr2.Value(fj);
 346       if (MinFAx.Value(fi) > MaxFBx.Value(fj) ||
 347           MinFBx.Value(fj) > MaxFAx.Value(fi)) {
 348         j += NbTrian2;
 349         continue;
 350       }
 351       else {
 352         if (MinFAy.Value(fi) > MaxFBy.Value(fj) ||
 353             MinFBy.Value(fj) > MaxFAy.Value(fi)) {
 354           j += NbTrian2;
 355           continue;
 356         }
 357         else if (MinFAz.Value(fi) > MaxFBz.Value(fj) ||
 358                  MinFBz.Value(fj) > MaxFAz.Value(fi)) {
 359           j += NbTrian2;
 360           continue;
 361         }
 362       }
 363       i -= NbTrian1;
 364       j += NbTrian2;
 365
 366       for (ii = 1;ii <= NbTrian1;ii++) {
 367         i++;
 368         j -= NbTrian2;
 369         const gp_Pnt& A1 = TA1.Value(i);
 370         const gp_Pnt& A2 = TA2.Value(i);
 371         const gp_Pnt& A3 = TA3.Value(i);
 372         gp_Vec OA1(A1.X(), A1.Y(), A1.Z());
 373         gp_Vec VA0 = gp_Vec(A1,A2);
 374         gp_Vec VA  = gp_Vec(A1,A3);
 375         VA0.Cross(VA);
 376         VA0.Normalize();
 377
 378         for (jj = 1;jj <= NbTrian2;jj++) {
 379           j++;
 380           if(MinAx.Value(i) > MaxBx.Value(j) ||
 381              MinBx.Value(j) > MaxAx.Value(i)) {
 382             continue;
 383           }
 384           else {
 385             if (MinAy.Value(i) > MaxBy.Value(j) ||
 386                 MinBy.Value(j) > MaxAy.Value(i)) {
 387               continue;
 388             }
 389             else if (MinAz.Value(i) > MaxBz.Value(j) ||
 390                      MinBz.Value(j) > MaxAz.Value(i)) {
 391               continue;
 392             }
 393           }
 394           const gp_Pnt& B1 = TB1.Value(j);
 395           const gp_Pnt& B2 = TB2.Value(j);
 396           const gp_Pnt& B3 = TB3.Value(j);
 397           gp_Vec VB0(B1,B2);
 398           gp_Vec VB(B1,B3);
 399           VB0.Cross(VB);
 400           VB0.Normalize();
 401           gp_Vec OB1(B1.XYZ());
 402           gp_Vec V1(B1,A1);
 403           gp_Vec V2(B1,A2);
 404           gp_Vec V3(B1,A3);
 405           Standard_Real h1 = VB0.Dot(V1);
 406           Standard_Real h2 = VB0.Dot(V2);
 407           Standard_Real h3 = VB0.Dot(V3);
 408           Standard_Real ah1 = Abs(h1);
 409           Standard_Real ah2 = Abs(h2);
 410           Standard_Real ah3 = Abs(h3);
 411
 412           if(!SAMESIGN(h1,h2,h3)) {
 413             myFirstTime = Standard_True;
 414             if (h1 == 0.0) {
 415               if (IsInside(A1,B1,B2,B3,VB0)) {
 416                 myFirstTime = Standard_False;
 417                 myBegPoint = A1;
 418                 if (h2 == 0.0) {
 419                   if (IsInside(A2,B1,B2,B3,VB0)) {
 420                     if (!IsEqual(A2,myBegPoint)) {
 421                       myEndPoint = A2;
 422                       InsertInMap();
 423                       continue;
 424                     }
 425                   }
 426                 }
 427                 if (h3 == 0.0) {
 428                   if (IsInside(A3,B1,B2,B3,VB0)) {
 429                     if (!IsEqual(A3,myBegPoint)) {
 430                       myEndPoint = A3;
 431                       InsertInMap();
 432                       continue;
 433                     }
 434                   }
 435                 }
 436               }
 437             }
 438             if (h2 == 0.0) {   //  h1 <> 0
 439               if (IsInside(A2,B1,B2,B3,VB0)) {
 440                 myFirstTime = Standard_False;
 441                 myBegPoint = A2;
 442                 if (h3 == 0.0) {
 443                   if (IsInside(A3,B1,B2,B3,VB0)) {
 444                     if (!IsEqual(A3,myBegPoint)) {
 445                       myEndPoint = A3;
 446                       InsertInMap();
 447                       continue;
 448                     }
 449                   }
 450                 }
 451               }
 452             }
 453             if (h3 == 0.0) {  //  h1 <> 0  and  h2 <> 0
 454               if (IsInside(A3,B1,B2,B3,VB0)) {
 455                 myFirstTime = Standard_False;
 456                 myBegPoint = A3;
 457               }
 458             }
 459             if (Intersect(B1,B2,B3,OB1,VB0,V1,V2,h1,h2,ah1,ah2))
 460               continue;
 461             if (Intersect(B1,B2,B3,OB1,VB0,V2,V3,h2,h3,ah2,ah3))
 462               continue;
 463             if (Intersect(B1,B2,B3,OB1,VB0,V3,V1,h3,h1,ah3,ah1))
 464               continue;
 465
 466             V1 = gp_Vec(A1,B1);
 467             V2 = gp_Vec(A1,B2);
 468             V3 = gp_Vec(A1,B3);
 469             h1 = VA0.Dot(V1);
 470             h2 = VA0.Dot(V2);
 471             h3 = VA0.Dot(V3);
 472             ah1 = Abs(h1);
 473             ah2 = Abs(h2);
 474             ah3 = Abs(h3);
 475
 476             if(!SAMESIGN(h1,h2,h3)) {
 477               if (h1 == 0.0) {
 478                 if (IsInside(B1,A1,A2,A3,VA0)) {
 479                   if (myFirstTime) {
 480                     myFirstTime = Standard_False;
 481                     myBegPoint = B1;
 482                   }
 483                   else {
 484                     if (!IsEqual(B1,myBegPoint)) {
 485                       myEndPoint = B1;
 486                       InsertInMap();
 487                       continue;
 488                     }
 489                   }
 490                   if (h2 == 0.0) {
 491                     if (IsInside(B2,A1,A2,A3,VA0)) {
 492                       if (!IsEqual(B2,myBegPoint)) {
 493                         myEndPoint = B2;
 494                         InsertInMap();
 495                         continue;
 496                       }
 497                     }
 498                   }
 499                   if (h3 == 0.0) {
 500                     if (IsInside(B3,A1,A2,A3,VA0)) {
 501                       if (!IsEqual(B3,myBegPoint)) {
 502                         myEndPoint = B3;
 503                         InsertInMap();
 504                         continue;
 505                       }
 506                     }
 507                   }
 508                 }
 509               }
 510               if (h2 == 0.0) {   //  h1 <> 0
 511                 if (IsInside(B2,A1,A2,A3,VA0)) {
 512                   if (myFirstTime) {
 513                     myFirstTime = Standard_False;
 514                     myBegPoint = B2;
 515                   }
 516                   else {
 517                     if (!IsEqual(B2,myBegPoint)) {
 518                       myEndPoint = B2;
 519                       InsertInMap();
 520                       continue;
 521                     }
 522                   }
 523                   if (h3 == 0.0) {
 524                     if (IsInside(B3,A1,A2,A3,VA0)) {
 525                       if (!IsEqual(B3,myBegPoint)) {
 526                         myEndPoint = B3;
 527                         InsertInMap();
 528                         continue;
 529                       }
 530                     }
 531                   }
 532                 }
 533               }
 534               if (h3 == 0.0) {  //  h1 <> 0  and  h2 <> 0
 535                 if (IsInside(B3,A1,A2,A3,VA0)) {
 536                   if (myFirstTime) {
 537                     myBegPoint = B3;
 538                     myFirstTime = Standard_False;
 539                   }
 540                   else {
 541                     if (!IsEqual(B3,myBegPoint)) {
 542                       myEndPoint = B3;
 543                       InsertInMap();
 544                       continue;
 545                     }
 546                   }
 547                 }
 548               }
 549               if (Intersect(A1,A2,A3,OA1,VA0,V1,V2,h1,h2,ah1,ah2))
 550                 continue;
 551               if (Intersect(A1,A2,A3,OA1,VA0,V2,V3,h2,h3,ah2,ah3))
 552                 continue;
 553               if (Intersect(A1,A2,A3,OA1,VA0,V3,V1,h3,h1,ah3,ah1))
 554                 continue;
 555             }
 556           }
 557         }
 558       }
 559     }
 560   }
 561 }
 562
 563 //=======================================================================
 564 //function : Intersect
 565 //purpose  :
 566 //=======================================================================
 567 #define DIFFSIGN(a,b) (((a)>0.0 && (b)<0.0) || ((a)<0.0 && (b)>0.0))
 568
 569 Standard_Boolean IntPoly_ShapeSection::Intersect(const gp_Pnt& S1,
 570                                                  const gp_Pnt& S2,
 571                                                  const gp_Pnt& S3,
 572                                                  const gp_Vec& OS1,
 573                                                  const gp_Vec& VS0,
 574                                                  const gp_Vec& V1,
 575                                                  const gp_Vec& V2,
 576                                                  Standard_Real& h1,
 577                                                  Standard_Real& h2,
 578                                                  Standard_Real& ah1,
 579                                                  Standard_Real& ah2)
 580 {
 581   if (DIFFSIGN(h1,h2)) {
 582     gp_XYZ OP;
 583     OP.SetLinearForm(ah2/(ah1+ah2),V1.XYZ(),
 584                      ah1/(ah1+ah2),V2.XYZ(),
 585                      OS1.XYZ());
 586     gp_Pnt P(OP.X(),OP.Y(),OP.Z());
 587     if (IsInside(P,S1,S2,S3,VS0)) {
 588       if (myFirstTime) {
 589         myFirstTime = Standard_False;
 590         myBegPoint = P;
 591       }
 592       else {
 593         if (!IsEqual(P,myBegPoint)) {
 594           myEndPoint = P;
 595           InsertInMap();
 596           return (Standard_True);
 597         }
 598       }
 599     }
 600   }
 601   return (Standard_False);
 602 }
 603
 604
 605 //=======================================================================
 606 //function : IsEqual
 607 //purpose  :
 608 //=======================================================================
 609 #define epsilon 0.000000000001
 610
 611 Standard_Boolean IntPoly_ShapeSection::IsEqual(const gp_Pnt& Pt1, const gp_Pnt& Pt2)
 612 {
 613   return (((Pt1.X() <= Pt2.X() && Pt2.X() < Pt1.X()+epsilon) ||
 614            (Pt2.X() <= Pt1.X() && Pt1.X() < Pt2.X()+epsilon)) &&
 615           ((Pt1.Y() <= Pt2.Y() && Pt2.Y() < Pt1.Y()+epsilon) ||
 616            (Pt2.Y() <= Pt1.Y() && Pt1.Y() < Pt2.Y()+epsilon)) &&
 617           ((Pt1.Z() <= Pt2.Z() && Pt2.Z() < Pt1.Z()+epsilon) ||
 618            (Pt2.Z() <= Pt1.Z() && Pt1.Z() < Pt2.Z()+epsilon)));
 619 }
 620
 621
 622 //=======================================================================
 623 //function : IsInside
 624 //purpose  :
 625 //=======================================================================
 626
 627 static Standard_Boolean IsInside1(Standard_Real PU,
 628                                   Standard_Real PV,
 629                                   Standard_Real *Pu,
 630                                   Standard_Real *Pv)
 631 {
 632   Standard_Real u = Pu[0] - PU;
 633   Standard_Real v = Pv[0] - PV;
 634   Standard_Real nu,nv;
 635   Standard_Integer SH = (v < 0.0)? -1 : 1;
 636   Standard_Integer NH,i,ip1 = 1;
 637   Standard_Integer NbCrossing = 0;
 638
 639   for (i = 0;i < 3;i++,ip1++) {
 640     nu = Pu[ip1] - PU;
 641     nv = Pv[ip1] - PV;
 642     NH = (nv < 0.0)? -1 : 1;
 643     if (NH != SH) {
 644       if (u > 0.0 && nu > 0.0) NbCrossing++;
 645       else if (u > 0.0 || nu > 0.0)
 646         if ((u - v * (nu - u) / (nv - v)) > 0.0) NbCrossing++;
 647       SH = NH;
 648     }
 649     u = nu;
 650     v = nv;
 651   }
 652   return (NbCrossing & 1);
 653 }
 654
 655 #define TOLBRUIT 0.00000000000001
 656
 657 Standard_Boolean IntPoly_ShapeSection::IsInside(const gp_Pnt& P,
 658                                                 const gp_Pnt& P1,
 659                                                 const gp_Pnt& P2,
 660                                                 const gp_Pnt& P3,
 661                                                 const gp_Vec& N0)
 662 {
 663   if (IsEqual(P,P1) || IsEqual(P,P2) || IsEqual(P,P3))
 664     return(Standard_True);
 665   Standard_Real Nx = Abs(N0.X());
 666   Standard_Real Ny = Abs(N0.Y());
 667   Standard_Real Nz = Abs(N0.Z());
 668   Standard_Real PU,PV,Pu[4],Pv[4];
 669   Standard_Integer Ind = 1;    //-- 1:x  2:y  3:z
 670   if (Nx > Ny) {               //-- x or z
 671     if (Nx > Nz) Ind = 1;
 672     else         Ind = 3;
 673   }
 674   else {                       //-- y or z
 675     if (Ny > Nz) Ind = 2;
 676     else         Ind = 3;
 677   }
 678   if (Ind == 1) {
 679     PU = P.Y();
 680     PV = P.Z();
 681     Pu[0] = Pu[3] = P1.Y();
 682     Pv[0] = Pv[3] = P1.Z();
 683     Pu[1] = P2.Y();
 684     Pv[1] = P2.Z();
 685     Pu[2] = P3.Y();
 686     Pv[2] = P3.Z();
 687   }
 688   else if (Ind == 2) {
 689     PU = P.Z();
 690     PV = P.X();
 691     Pu[0] = Pu[3] = P1.Z();
 692     Pv[0] = Pv[3] = P1.X();
 693     Pu[1] = P2.Z();
 694     Pv[1] = P2.X();
 695     Pu[2] = P3.Z();
 696     Pv[2] = P3.X();
 697   }
 698   else {
 699     PU = P.X();
 700     PV = P.Y();
 701     Pu[0] = Pu[3] = P1.X();
 702     Pv[0] = Pv[3] = P1.Y();
 703     Pu[1] = P2.X();
 704     Pv[1] = P2.Y();
 705     Pu[2] = P3.X();
 706     Pv[2] = P3.Y();
 707   }
 708   if (IsInside1(PU,PV,Pu,Pv))
 709     return(Standard_True);
 710
 711   if (IsInside1(PU+TOLBRUIT,PV,Pu,Pv))
 712     return(Standard_True);
 713
 714   if (IsInside1(PU-TOLBRUIT,PV,Pu,Pv))
 715     return(Standard_True);
 716
 717   if (IsInside1(PU,PV+TOLBRUIT,Pu,Pv))
 718     return(Standard_True);
 719
 720   if (IsInside1(PU,PV-TOLBRUIT,Pu,Pv))
 721     return(Standard_True);
 722
 723   return(Standard_False);
 724 }
 725
 726
 727 //=======================================================================
 728 //function : InsertInMap
 729 //purpose  :
 730 //=======================================================================
 731
 732 void IntPoly_ShapeSection::InsertInMap()
 733 {
 734   Standard_Integer Index;
 735   Standard_Real Big = Precision::Infinite();
 736
 737   if (myMapBegPoints.Contains(myBegPoint)) {
 738     Index = myMapBegPoints.FindIndex(myBegPoint);
 739     Insert(myMapEndPoints.FindKey(Index),myBegPoint,myEndPoint);
 740     myCpt++;
 741     myMapBegPoints.Substitute(Index,gp_Pnt(Big,myCpt,myCpt));
 742     myMapEndPoints.Substitute(Index,gp_Pnt(Big,myCpt,myCpt));
 743   }
 744   else if (myMapEndPoints.Contains(myEndPoint)) {
 745     Index = myMapEndPoints.FindIndex(myEndPoint);
 746     Insert(myMapBegPoints.FindKey(Index),myEndPoint,myBegPoint);
 747     myCpt++;
 748     myMapBegPoints.Substitute(Index,gp_Pnt(Big,myCpt,myCpt));
 749     myMapEndPoints.Substitute(Index,gp_Pnt(Big,myCpt,myCpt));
 750   }
 751   else {
 752     myMapBegPoints.Add(myBegPoint);
 753     myMapEndPoints.Add(myEndPoint);
 754   }
 755 }
 756
 757
 758 //=======================================================================
 759 //function : Insert
 760 //purpose  :
 761 //=======================================================================
 762
 763 void IntPoly_ShapeSection::Insert(const gp_Pnt& OldPnt,
 764                                   const gp_Pnt& ComPnt,
 765                                   const gp_Pnt& NewPnt)
 766 {
 767   Standard_Integer i = 0;
 768   Standard_Integer NbSection = mySection.Length();
 769   Standard_Boolean IsInSection = Standard_False;
 770
 771   while (i < NbSection) {
 772     i++;
 773     TColgp_SequenceOfPnt& CurSection = mySection.ChangeValue(i);
 774     if (IsEqual(OldPnt,CurSection.First())) {
 775       IsInSection = Standard_True;
 776       CurSection.Prepend(ComPnt);
 777       CurSection.Prepend(NewPnt);
 778       break;
 779     }
 780     if (IsEqual(OldPnt,CurSection.Last())) {
 781       IsInSection = Standard_True;
 782       CurSection.Append(ComPnt);
 783       CurSection.Append(NewPnt);
 784       break;
 785     }
 786     if (IsEqual(NewPnt,CurSection.First())) {
 787       IsInSection = Standard_True;
 788       CurSection.Prepend(ComPnt);
 789       CurSection.Prepend(OldPnt);
 790       break;
 791     }
 792     if (IsEqual(NewPnt,CurSection.Last())) {
 793       IsInSection = Standard_True;
 794       CurSection.Append(ComPnt);
 795       CurSection.Append(OldPnt);
 796       break;
 797     }
 798   }
 799   if (!(IsInSection)) {
 800     TColgp_SequenceOfPnt EmptySec;
 801     EmptySec.Append(OldPnt);
 802     EmptySec.Append(ComPnt);
 803     EmptySec.Append(NewPnt);
 804     mySection.Append(EmptySec);
 805   }
 806 }
 807
 808
 809 //=======================================================================
 810 //function : Concat
 811 //purpose  :
 812 //=======================================================================
 813
 814 Standard_Integer IntPoly_ShapeSection::Concat(const gp_Pnt& BegPnt,
 815                                               const gp_Pnt& EndPnt,
 816                                               gp_Pnt& OutPnt)
 817 {
 818   Standard_Integer i = 0;
 819   Standard_Integer NbSection = mySection.Length();
 820   Standard_Integer ConcatIdx = 0;
 821
 822   while (i < NbSection) {
 823     i++;
 824     TColgp_SequenceOfPnt& CurSection = mySection.ChangeValue(i);
 825     if (IsEqual(BegPnt,CurSection.First())) {
 826       ConcatIdx = 1;
 827       myIndex = i;
 828       CurSection.Prepend(EndPnt);
 829       OutPnt = EndPnt;
 830       break;
 831     }
 832     if (IsEqual(BegPnt,CurSection.Last())) {
 833       ConcatIdx = 2;
 834       myIndex = i;
 835       CurSection.Append(EndPnt);
 836       OutPnt = EndPnt;
 837       break;
 838     }
 839     if (IsEqual(EndPnt,CurSection.First())) {
 840       ConcatIdx = 1;
 841       myIndex = i;
 842       CurSection.Prepend(BegPnt);
 843       OutPnt = BegPnt;
 844       break;
 845     }
 846     if (IsEqual(EndPnt,CurSection.Last())) {
 847       ConcatIdx = 2;
 848       myIndex = i;
 849       CurSection.Append(BegPnt);
 850       OutPnt = BegPnt;
 851       break;
 852     }
 853   }
 854   if (ConcatIdx == 0) {
 855     TColgp_SequenceOfPnt EmptySec;
 856     EmptySec.Append(BegPnt);
 857     EmptySec.Append(EndPnt);
 858     mySection.Append(EmptySec);
 859     NbSection++;
 860     myIndex = NbSection;
 861   }
 862   return ConcatIdx;
 863 }
 864
 865
 866 //=======================================================================
 867 //function : ConcatSection
 868 //purpose  :
 869 //=======================================================================
 870
 871 void IntPoly_ShapeSection::ConcatSection(TColgp_SequenceOfPnt& Section,
 872                                          const Standard_Integer NbSection,
 873                                          const Standard_Integer Index)
 874 {
 875   Standard_Integer j;
 876   Standard_Integer i = Index;
 877   gp_Pnt BegPnt = Section.First();
 878   gp_Pnt EndPnt = Section.Last();
 879
 880   while (i <= NbSection) {
 881     TColgp_SequenceOfPnt& CurSection = mySection.ChangeValue(i);
 882     Standard_Integer CurSection_Length = CurSection.Length();
 883     if (IsEqual(BegPnt,CurSection.First())) {
 884       for (j = 2;j <= CurSection_Length;j++)
 885         Section.Prepend(CurSection.Value(j));
 886       mySection.Remove(i);
 887       ConcatSection(Section,NbSection-1,Index);
 888       break;
 889     }
 890     else
 891       if (IsEqual(BegPnt,CurSection.Last())) {
 892         for (j = CurSection_Length-1;j >= 1;j--)
 893           Section.Prepend(CurSection.Value(j));
 894         mySection.Remove(i);
 895         ConcatSection(Section,NbSection-1,Index);
 896         break;
 897       }
 898       else
 899         if (IsEqual(EndPnt,CurSection.First())) {
 900           for (j = 2;j <= CurSection_Length;j++)
 901             Section.Append(CurSection.Value(j));
 902           mySection.Remove(i);
 903           ConcatSection(Section,NbSection-1,Index);
 904           break;
 905         }
 906         else
 907           if (IsEqual(EndPnt,CurSection.Last())) {
 908             for (j = CurSection_Length-1;j >= 1;j--)
 909               Section.Append(CurSection.Value(j));
 910             mySection.Remove(i);
 911             ConcatSection(Section,NbSection-1,Index);
 912             break;
 913           }
 914     i++;
 915   }
 916 }
 917
 918
 919 //=======================================================================
 920 //function : ForwContruction
 921 //purpose  :
 922 //=======================================================================
 923
 924 void IntPoly_ShapeSection::ForwConstruction(const gp_Pnt& Point)
 925 {
 926   if (myMapBegPoints.Contains(Point)) {
 927     Standard_Integer Index = myMapBegPoints.FindIndex(Point);
 928     gp_Pnt Pnt = myMapEndPoints.FindKey(Index);
 929     (mySection.ChangeValue(myIndex)).Append(Pnt);
 930     Standard_Real Big = Precision::Infinite();
 931     myCpt++;
 932     myMapBegPoints.Substitute(Index,gp_Pnt(Big,myCpt,myCpt));
 933     myMapEndPoints.Substitute(Index,gp_Pnt(Big,myCpt,myCpt));
 934     ForwConstruction(Pnt);
 935   }
 936 }
 937
 938
 939 //=======================================================================
 940 //function : PrevContruction
 941 //purpose  :
 942 //=======================================================================
 943
 944 void IntPoly_ShapeSection::PrevConstruction(const gp_Pnt& Point)
 945 {
 946   if (myMapEndPoints.Contains(Point)) {
 947     Standard_Integer Index = myMapEndPoints.FindIndex(Point);
 948     gp_Pnt Pnt = myMapBegPoints.FindKey(Index);
 949     (mySection.ChangeValue(myIndex)).Prepend(Pnt);
 950     Standard_Real Big = Precision::Infinite();
 951     myCpt++;
 952     myMapBegPoints.Substitute(Index,gp_Pnt(Big,myCpt,myCpt));
 953     myMapEndPoints.Substitute(Index,gp_Pnt(Big,myCpt,myCpt));
 954     PrevConstruction(Pnt);
 955   }
 956 }
 957
 958
 959 //=======================================================================
 960 //function : NbEdges
 961 //purpose  :
 962 //=======================================================================
 963
 964 Standard_Integer IntPoly_ShapeSection::NbEdges()
 965 { return myNbEdges; }
 966
 967
 968 //=======================================================================
 969 //function : Edge
 970 //purpose  :
 971 //=======================================================================
 972
 973 TopoDS_Edge IntPoly_ShapeSection::Edge(const Standard_Integer Index)
 974 {
 975   const TColgp_SequenceOfPnt& CurSection = mySection.ChangeValue(Index);
 976   Standard_Integer NbPoints = CurSection.Length();
 977   TColgp_Array1OfPnt TabPnt(1,NbPoints);
 978   //gp_Pnt CurPoint;
 979   for (Standard_Integer i = 1 ; i <= NbPoints ; i++) {
 980     TabPnt.SetValue(i,CurSection.Value(i));
 981   }
 982   Handle(Poly_Polygon3D) Pol = new Poly_Polygon3D(TabPnt);
 983   TopoDS_Edge Edge;
 984   BRep_Builder B;
 985   B.MakeEdge(Edge,Pol);
 986   return Edge;
 987 }
 988
 989
 990
 991
 992
 993
 994