src/IntCurveSurface/IntCurveSurface_Polyhedron.gxx

   1 // Created on: 1993-02-03
   2 // Created by: Laurent BUCHARD
   3 // Copyright (c) 1993-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 under
   9 // the terms of the GNU Lesser General Public License 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 <gp.hxx>
  18 #include <gp_Pnt.hxx>
  19 #include <gp_Vec.hxx>
  20 #include <gp_Dir.hxx>
  21 #include <gp_Lin.hxx>
  22 #include <TColgp_Array2OfPnt.hxx>
  23 #include <TColStd_Array2OfReal.hxx>
  24
  25 #include <Bnd_Array1OfBox.hxx>
  26
  27 #include <Standard_ConstructionError.hxx>
  28
  29 #include <stdio.h>
  30
  31 #define CHECKBOUNDS 0
  32
  33
  34 //-----------------------------------------------------
  35 #define LONGUEUR_MINI_EDGE_TRIANGLE 1e-15
  36
  37
  38 //=======================================================================
  39 //function : IntCurveSurface_Polyhedron
  40 //purpose  :
  41 //=======================================================================
  42 IntCurveSurface_Polyhedron::IntCurveSurface_Polyhedron (const ThePSurface&     Surface,
  43                                                         const Standard_Integer nbdU,
  44                                                         const Standard_Integer nbdV,
  45                                                         const Standard_Real    u1,
  46                                                         const Standard_Real    v1,
  47                                                         const Standard_Real    u2,
  48                                                         const Standard_Real    v2)
  49      : nbdeltaU((nbdU<3)? 3 : nbdU),
  50        nbdeltaV((nbdV<3)? 3 : nbdV),
  51        TheDeflection(Epsilon(100.)),
  52        C_MyPnts(NULL),C_MyU(NULL),C_MyV(NULL),C_MyIsOnBounds(NULL)
  53 {
  54   Standard_Integer t = (nbdeltaU+1)*(nbdeltaV+1)+1;
  55   gp_Pnt           *CMyPnts = new gp_Pnt[t];         C_MyPnts = (void *)CMyPnts;
  56   Standard_Real    *CMyU    = new Standard_Real[t];  C_MyU    = (void *)CMyU;
  57   Standard_Real    *CMyV    = new Standard_Real[t];  C_MyV    = (void *)CMyV;
  58
  59 //  Modified by Sergey KHROMOV - Fri Dec  7 12:03:46 2001 Begin
  60   Standard_Boolean *CMyIsOnBounds = new Standard_Boolean[t];
  61
  62   C_MyIsOnBounds = (void *)CMyIsOnBounds;
  63 //  Modified by Sergey KHROMOV - Fri Dec  7 12:03:47 2001 End
  64   Init(Surface,u1,v1,u2,v2);
  65 }
  66
  67 //=======================================================================
  68 //function : IntCurveSurface_Polyhedron
  69 //purpose  :
  70 //=======================================================================
  71 IntCurveSurface_Polyhedron::IntCurveSurface_Polyhedron (const ThePSurface&     Surface,
  72                                                         const TColStd_Array1OfReal&    Upars,
  73                                                         const TColStd_Array1OfReal&    Vpars)
  74     : nbdeltaU(Upars.Length()-1),
  75       nbdeltaV(Vpars.Length()-1),
  76       TheDeflection(Epsilon(100.)),
  77       C_MyPnts(NULL),C_MyU(NULL),C_MyV(NULL),C_MyIsOnBounds(NULL)
  78 {
  79   Standard_Integer t = (nbdeltaU+1)*(nbdeltaV+1)+1;
  80   gp_Pnt           *CMyPnts = new gp_Pnt[t];         C_MyPnts = (void *)CMyPnts;
  81   Standard_Real    *CMyU    = new Standard_Real[t];  C_MyU    = (void *)CMyU;
  82   Standard_Real    *CMyV    = new Standard_Real[t];  C_MyV    = (void *)CMyV;
  83
  84 //  Modified by Sergey KHROMOV - Fri Dec  7 12:03:46 2001 Begin
  85   Standard_Boolean *CMyIsOnBounds = new Standard_Boolean[t];
  86
  87   C_MyIsOnBounds = (void *)CMyIsOnBounds;
  88 //  Modified by Sergey KHROMOV - Fri Dec  7 12:03:47 2001 End
  89   Init(Surface, Upars, Vpars);
  90 }
  91
  92
  93 void IntCurveSurface_Polyhedron::Destroy() {
  94   //-- printf("\n IntCurveSurface_Polyhedron::Destroy\n");
  95   gp_Pnt *CMyPnts     = (gp_Pnt *)C_MyPnts;       if(C_MyPnts) delete [] CMyPnts;
  96   Standard_Real *CMyU = (Standard_Real *)C_MyU;   if(C_MyU)    delete [] CMyU;
  97   Standard_Real *CMyV = (Standard_Real *)C_MyV;   if(C_MyV)    delete [] CMyV;
  98
  99 //  Modified by Sergey KHROMOV - Fri Dec  7 12:03:46 2001 Begin
 100   Standard_Boolean *CMyIsOnBounds = (Standard_Boolean *)C_MyIsOnBounds;
 101
 102   if(C_MyIsOnBounds)    delete [] CMyIsOnBounds;
 103 //  Modified by Sergey KHROMOV - Fri Dec  7 12:03:47 2001 End
 104
 105   C_MyPnts=C_MyU=C_MyV=C_MyIsOnBounds=NULL;
 106 }
 107
 108 //=======================================================================
 109 //function : Init
 110 //purpose  :
 111 //=======================================================================
 112 void IntCurveSurface_Polyhedron::Init(const ThePSurface&     Surface,
 113                                       const Standard_Real U0,
 114                                       const Standard_Real V0,
 115                                       const Standard_Real U1,
 116                                       const Standard_Real V1) {
 117   //#define DEBUGDUMP
 118   Standard_Integer i1,i2;
 119   Standard_Real    U,V;
 120   Standard_Real    U1mU0sNbdeltaU = (U1-U0)/(Standard_Real)nbdeltaU;
 121   Standard_Real    V1mV0sNbdeltaV = (V1-V0)/(Standard_Real)nbdeltaV;
 122   gp_Pnt TP;
 123   Standard_Integer Index=1;
 124   //-- --------------------------------------------------
 125   //-- Index varie de 1 -> (nbdu+1)*(nbdv+1)
 126   //-- V est la colonne
 127   //-- U est la ligne
 128   //-- --------------------------------------------------
 129   gp_Pnt *CMyPnts     = (gp_Pnt *)C_MyPnts;
 130   Standard_Real *CMyU = (Standard_Real *)C_MyU;
 131   Standard_Real *CMyV = (Standard_Real *)C_MyV;
 132   Standard_Boolean *CMyIsOnBounds = (Standard_Boolean *)C_MyIsOnBounds;
 133
 134   for (i1 = 0, U = U0; i1 <= nbdeltaU; i1++, U+= U1mU0sNbdeltaU) {
 135     for (i2 = 0, V = V0; i2 <= nbdeltaV; i2++, V+= V1mV0sNbdeltaV ) {
 136       ThePSurfaceTool::D0(Surface,U,V,TP);
 137       //-- Point(TP,i1, i2,U,V);
 138       CMyPnts[Index] = TP;
 139       CMyU[Index]    = U;
 140       CMyV[Index]    = V;
 141 //  Modified by Sergey KHROMOV - Fri Dec  7 12:07:51 2001
 142       CMyIsOnBounds[Index] = (i1 == 0 || i1 == nbdeltaU ||
 143                               i2 == 0 || i2 == nbdeltaV);
 144 //  Modified by Sergey KHROMOV - Fri Dec  7 12:07:52 2001
 145       TheBnd.Add(TP);
 146       Index++;
 147     }
 148   }
 149   //-- Calcul de la deflection Triangle <-> Point milieu
 150   Standard_Real tol=0.0; Standard_Integer nbtriangles = NbTriangles();
 151   for (i1=1; i1<=nbtriangles; i1++) {
 152     Standard_Real tol1 = DeflectionOnTriangle(Surface,i1);
 153     if(tol1>tol) tol=tol1;
 154   }
 155   //-- Calcul de la deflection Bord <-> Point milieu
 156
 157
 158   DeflectionOverEstimation(tol*1.2);
 159   FillBounding();
 160
 161 //  Modified by Sergey KHROMOV - Fri Dec  7 11:23:33 2001 Begin
 162   Standard_Real aDeflection;
 163
 164   TheBorderDeflection = RealFirst();
 165
 166 // Compute the deflection on the lower bound (U-isoline) of the surface.
 167   aDeflection = ComputeBorderDeflection(Surface, U0, V0, V1, Standard_True);
 168
 169   if (aDeflection > TheBorderDeflection)
 170     TheBorderDeflection = aDeflection;
 171
 172 // Compute the deflection on the upper bound (U-isoline) of the surface.
 173   aDeflection = ComputeBorderDeflection(Surface, U1, V0, V1, Standard_True);
 174
 175   if (aDeflection > TheBorderDeflection)
 176     TheBorderDeflection = aDeflection;
 177
 178 // Compute the deflection on the lower bound (V-isoline) of the surface.
 179   aDeflection = ComputeBorderDeflection(Surface, V0, U0, U1, Standard_False);
 180
 181   if (aDeflection > TheBorderDeflection)
 182     TheBorderDeflection = aDeflection;
 183
 184 // Compute the deflection on the upper bound (V-isoline) of the surface.
 185   aDeflection = ComputeBorderDeflection(Surface, V1, U0, U1, Standard_False);
 186
 187   if (aDeflection > TheBorderDeflection)
 188     TheBorderDeflection = aDeflection;
 189
 190 //  Modified by Sergey KHROMOV - Fri Dec  7 11:23:34 2001 End
 191
 192   #ifdef OCCT_DEBUG_DUMP
 193     Dump();
 194   #endif
 195 }
 196 //=======================================================================
 197 //function : Init
 198 //purpose  :
 199 //=======================================================================
 200 void IntCurveSurface_Polyhedron::Init(const ThePSurface&     Surface,
 201                                       const TColStd_Array1OfReal& Upars,
 202                                       const TColStd_Array1OfReal& Vpars) {
 203   //#define DEBUGDUMP
 204   Standard_Integer i1,i2;
 205   Standard_Real    U,V;
 206   gp_Pnt TP;
 207   Standard_Integer Index=1;
 208   //-- --------------------------------------------------
 209   //-- Index varie de 1 -> (nbdu+1)*(nbdv+1)
 210   //-- V est la colonne
 211   //-- U est la ligne
 212   //-- --------------------------------------------------
 213   gp_Pnt *CMyPnts     = (gp_Pnt *)C_MyPnts;
 214   Standard_Real *CMyU = (Standard_Real *)C_MyU;
 215   Standard_Real *CMyV = (Standard_Real *)C_MyV;
 216   Standard_Boolean *CMyIsOnBounds = (Standard_Boolean *)C_MyIsOnBounds;
 217   Standard_Integer i0 = Upars.Lower(), j0 = Vpars.Lower();
 218
 219   for (i1 = 0; i1 <= nbdeltaU; i1++) {
 220     U = Upars(i1+i0);
 221     for (i2 = 0; i2 <= nbdeltaV; i2++) {
 222       V = Vpars(i2+j0);
 223       ThePSurfaceTool::D0(Surface,U,V,TP);
 224       //-- Point(TP,i1, i2,U,V);
 225       CMyPnts[Index] = TP;
 226       CMyU[Index]    = U;
 227       CMyV[Index]    = V;
 228 //  Modified by Sergey KHROMOV - Fri Dec  7 12:07:51 2001
 229       CMyIsOnBounds[Index] = (i1 == 0 || i1 == nbdeltaU ||
 230                               i2 == 0 || i2 == nbdeltaV);
 231 //  Modified by Sergey KHROMOV - Fri Dec  7 12:07:52 2001
 232       TheBnd.Add(TP);
 233       Index++;
 234     }
 235   }
 236   //-- Calcul de la deflection Triangle <-> Point milieu
 237   Standard_Real tol=0.0; Standard_Integer nbtriangles = NbTriangles();
 238   for (i1=1; i1<=nbtriangles; i1++) {
 239     Standard_Real tol1 = DeflectionOnTriangle(Surface,i1);
 240     if(tol1>tol) tol=tol1;
 241   }
 242   //-- Calcul de la deflection Bord <-> Point milieu
 243
 244
 245   DeflectionOverEstimation(tol*1.2);
 246   FillBounding();
 247
 248 //  Modified by Sergey KHROMOV - Fri Dec  7 11:23:33 2001 Begin
 249   Standard_Real aDeflection;
 250
 251   TheBorderDeflection = RealFirst();
 252   Standard_Real U0 = Upars(i0);
 253   Standard_Real V0 = Vpars(j0);
 254   Standard_Real U1 = Upars(Upars.Upper());
 255   Standard_Real V1 = Vpars(Vpars.Upper());
 256
 257 // Compute the deflection on the lower bound (U-isoline) of the surface.
 258   aDeflection = ComputeBorderDeflection(Surface, U0, V0, V1, Standard_True);
 259
 260   if (aDeflection > TheBorderDeflection)
 261     TheBorderDeflection = aDeflection;
 262
 263 // Compute the deflection on the upper bound (U-isoline) of the surface.
 264   aDeflection = ComputeBorderDeflection(Surface, U1, V0, V1, Standard_True);
 265
 266   if (aDeflection > TheBorderDeflection)
 267     TheBorderDeflection = aDeflection;
 268
 269 // Compute the deflection on the lower bound (V-isoline) of the surface.
 270   aDeflection = ComputeBorderDeflection(Surface, V0, U0, U1, Standard_False);
 271
 272   if (aDeflection > TheBorderDeflection)
 273     TheBorderDeflection = aDeflection;
 274
 275 // Compute the deflection on the upper bound (V-isoline) of the surface.
 276   aDeflection = ComputeBorderDeflection(Surface, V1, U0, U1, Standard_False);
 277
 278   if (aDeflection > TheBorderDeflection)
 279     TheBorderDeflection = aDeflection;
 280
 281 //  Modified by Sergey KHROMOV - Fri Dec  7 11:23:34 2001 End
 282
 283   #ifdef OCCT_DEBUG_DUMP
 284     Dump();
 285   #endif
 286 }
 287 //=======================================================================
 288 //function : DeflectionOnTriangle
 289 //purpose  :
 290 //=======================================================================
 291 Standard_Real IntCurveSurface_Polyhedron::DeflectionOnTriangle (const ThePSurface& Surface,
 292                                                                 const Standard_Integer Triang) const
 293 {
 294   Standard_Integer i1,i2,i3;
 295   Triangle(Triang,i1,i2,i3);
 296   //-- Calcul de l equation du plan
 297   Standard_Real u1,v1,u2,v2,u3,v3;
 298   gp_Pnt P1,P2,P3;
 299   P1 = Point(i1,u1,v1);
 300   P2 = Point(i2,u2,v2);
 301   P3 = Point(i3,u3,v3);
 302   if(P1.SquareDistance(P2)<=LONGUEUR_MINI_EDGE_TRIANGLE) return(0);
 303   if(P1.SquareDistance(P3)<=LONGUEUR_MINI_EDGE_TRIANGLE) return(0);
 304   if(P2.SquareDistance(P3)<=LONGUEUR_MINI_EDGE_TRIANGLE) return(0);
 305   gp_XYZ XYZ1=P2.XYZ()-P1.XYZ();
 306   gp_XYZ XYZ2=P3.XYZ()-P2.XYZ();
 307   gp_XYZ XYZ3=P1.XYZ()-P3.XYZ();
 308   gp_Vec NormalVector((XYZ1^XYZ2)+(XYZ2^XYZ3)+(XYZ3^XYZ1));
 309   Standard_Real aNormLen = NormalVector.Magnitude();
 310   if (aNormLen < gp::Resolution()) {
 311     return 0.;
 312   }
 313   //
 314   NormalVector.Divide(aNormLen);
 315   //-- Standard_Real PolarDistance = NormalVector * P1.XYZ();
 316   //-- Calcul du point u,v  au centre du triangle
 317   Standard_Real u = (u1+u2+u3)/3.0;
 318   Standard_Real v = (v1+v2+v3)/3.0;
 319   gp_Pnt P =  ThePSurfaceTool::Value(Surface,u,v);
 320   gp_Vec P1P(P1,P);
 321   return(Abs(P1P.Dot(NormalVector)));
 322 }
 323 //=======================================================================
 324 //function : Parameters
 325 //purpose  :
 326 //=======================================================================
 327 void IntCurveSurface_Polyhedron::Parameters( const Standard_Integer Index
 328                                             ,Standard_Real &U
 329                                             ,Standard_Real &V) const
 330 {
 331 #if CHECKBOUNDS
 332   if(Index<0 || Index>((nbdeltaU+1)*(nbdeltaV+1))) {
 333     printf("\n Erreur IntCurveSurface_Polyhedron::Parameters\n");
 334   }
 335 #endif
 336   Standard_Real *CMyU = (Standard_Real *)C_MyU;
 337   U = CMyU[Index];
 338   Standard_Real *CMyV = (Standard_Real *)C_MyV;
 339   V = CMyV[Index];
 340 }
 341 //=======================================================================
 342 //function : DeflectionOverEstimation
 343 //purpose  : Set
 344 //=======================================================================
 345 void IntCurveSurface_Polyhedron::DeflectionOverEstimation(const Standard_Real flec)
 346 {
 347   if(flec<0.0001) {
 348     TheDeflection=0.0001;
 349     TheBnd.Enlarge(0.0001);
 350   }
 351   else {
 352     TheDeflection=flec;
 353     TheBnd.Enlarge(flec);
 354   }
 355 }
 356 //=======================================================================
 357 //function : DeflectionOverEstimation
 358 //purpose  : Get
 359 //=======================================================================
 360 Standard_Real IntCurveSurface_Polyhedron::DeflectionOverEstimation() const
 361 {
 362   return TheDeflection;
 363 }
 364 //=======================================================================
 365 //function : Bounding
 366 //purpose  :
 367 //=======================================================================
 368 const Bnd_Box& IntCurveSurface_Polyhedron::Bounding() const
 369 {
 370   return TheBnd;
 371 }
 372 //=======================================================================
 373 //function : FillBounding
 374 //purpose  :
 375 //=======================================================================
 376 void IntCurveSurface_Polyhedron::FillBounding()
 377 {
 378   TheComponentsBnd=new Bnd_HArray1OfBox(1, NbTriangles());
 379   Bnd_Box Boite;
 380   Standard_Integer np1, np2, np3;
 381   Standard_Integer nbtriangles = NbTriangles();
 382   for (Standard_Integer iTri=1; iTri<=nbtriangles; iTri++) {
 383     Triangle(iTri, np1, np2, np3);
 384     gp_Pnt p1(Point(np1));
 385     gp_Pnt p2(Point(np2));
 386     gp_Pnt p3(Point(np3));
 387     Boite.SetVoid();
 388     if(p1.SquareDistance(p2)>LONGUEUR_MINI_EDGE_TRIANGLE) {
 389       if(p1.SquareDistance(p3)>LONGUEUR_MINI_EDGE_TRIANGLE) {
 390         if(p2.SquareDistance(p3)>LONGUEUR_MINI_EDGE_TRIANGLE) {
 391           Boite.Add(p1);
 392           Boite.Add(p2);
 393           Boite.Add(p3);
 394           Boite.Enlarge(TheDeflection);
 395         }
 396       }
 397     }
 398     Boite.Enlarge(TheDeflection);
 399     TheComponentsBnd->SetValue(iTri,Boite);
 400   }
 401 }
 402 //=======================================================================
 403 //function : ComponentsBounding
 404 //purpose  :
 405 //=======================================================================
 406 const Handle(Bnd_HArray1OfBox)&
 407       IntCurveSurface_Polyhedron::ComponentsBounding() const
 408 {
 409  return TheComponentsBnd;
 410 }
 411 //=======================================================================
 412 //function : NbTriangles
 413 //purpose  :
 414 //=======================================================================
 415 Standard_Integer IntCurveSurface_Polyhedron::NbTriangles  () const
 416 {
 417   return nbdeltaU*nbdeltaV*2;
 418 }
 419 //=======================================================================
 420 //function : NbPoints
 421 //purpose  :
 422 //=======================================================================
 423 Standard_Integer IntCurveSurface_Polyhedron::NbPoints () const
 424 {
 425   return (nbdeltaU+1)*(nbdeltaV+1);
 426 }
 427 //=======================================================================
 428 //function : TriConnex
 429 //purpose  :
 430 //=======================================================================
 431 Standard_Integer IntCurveSurface_Polyhedron::TriConnex
 432   (const Standard_Integer Triang,
 433    const Standard_Integer Pivot,
 434    const Standard_Integer Pedge,
 435    Standard_Integer&      TriCon,
 436    Standard_Integer&      OtherP) const
 437 {
 438   Standard_Integer Pivotm1    = Pivot-1;
 439   Standard_Integer nbdeltaVp1 = nbdeltaV+1;
 440   Standard_Integer nbdeltaVm2 = nbdeltaV + nbdeltaV;
 441
 442 // Pivot position in the MaTriangle :
 443   Standard_Integer ligP = Pivotm1/nbdeltaVp1;
 444   Standard_Integer colP = Pivotm1 - ligP * nbdeltaVp1;
 445
 446 // Point sur Edge position in the MaTriangle and edge typ :
 447   Standard_Integer ligE =0, colE =0, typE =0;
 448   if (Pedge!=0) {
 449     ligE= (Pedge-1)/nbdeltaVp1;
 450     colE= (Pedge-1) - (ligE * nbdeltaVp1);
 451   // Horizontal
 452     if      (ligP==ligE) typE=1;
 453   // Vertical
 454     else if (colP==colE) typE=2;
 455   // Oblique
 456     else                 typE=3;
 457   }
 458   else {
 459     typE=0;
 460   }
 461
 462 // Triangle position General case :
 463
 464   Standard_Integer linT =0, colT =0;
 465   Standard_Integer linO =0, colO =0;
 466   Standard_Integer t =0, tt =0;
 467
 468   if (Triang!=0) {
 469     t = (Triang-1)/(nbdeltaVm2);
 470     tt= (Triang-1)-t*nbdeltaVm2;
 471     linT= 1+t;
 472     colT= 1+tt;
 473     if (typE==0) {
 474       if (ligP==linT) {
 475         ligE=ligP-1;
 476         colE=colP-1;
 477         typE=3;
 478       }
 479       else {
 480         if (colT==ligP+ligP) {
 481           ligE=ligP;
 482           colE=colP-1;
 483           typE=1;
 484         }
 485         else {
 486           ligE=ligP+1;
 487           colE=colP+1;
 488           typE=3;
 489         }
 490       }
 491     }
 492     switch (typE) {
 493     case 1:  // Horizontal
 494       if (linT==ligP) {
 495         linT++;
 496         linO=ligP+1;
 497         colO=(colP>colE)? colP : colE; //--colO=Max(colP, colE);
 498       }
 499       else {
 500         linT--;
 501         linO=ligP-1;
 502         colO=(colP<colE)? colP : colE;  //--colO=Min(colP, colE);
 503       }
 504       break;
 505     case 2:  // Vertical
 506       if (colT==(colP+colP)) {
 507         colT++;
 508         linO=(ligP>ligE)? ligP : ligE;  //--linO=Max(ligP, ligE);
 509         colO=colP+1;
 510       }
 511       else {
 512         colT--;
 513         linO=(ligP<ligE)? ligP : ligE;  //--linO=Min(ligP, ligE);
 514         colO=colP-1;
 515       }
 516       break;
 517     case 3:  // Oblique
 518       if ((colT&1)==0) {
 519         colT--;
 520         linO=(ligP>ligE)? ligP : ligE;  //--linO=Max(ligP, ligE);
 521         colO=(colP<colE)? colP : colE;  //--colO=Min(colP, colE);
 522       }
 523       else {
 524         colT++;
 525         linO=(ligP<ligE)? ligP : ligE;  //--linO=Min(ligP, ligE);
 526         colO=(colP>colE)? colP : colE;  //--colO=Max(colP, colE);
 527       }
 528       break;
 529     }
 530   }
 531   else {
 532     // Unknown Triangle position :
 533     if (Pedge==0) {
 534       // Unknown edge :
 535       linT=(1>ligP)? 1 : ligP;      //--linT=Max(1, ligP);
 536       colT=(1>(colP+colP))? 1 : (colP+colP);       //--colT=Max(1, colP+colP);
 537       if (ligP==0) linO=ligP+1;
 538       else         linO=ligP-1;
 539       colO=colP;
 540     }
 541     else {
 542       // Known edge We take the left or down connectivity :
 543       switch (typE) {
 544       case 1:  // Horizontal
 545         linT=ligP+1;
 546         colT=(colP>colE)? colP : colE;  //--colT=Max(colP,colE);
 547         colT+=colT;
 548         linO=ligP+1;
 549         colO=(colP>colE)? colP : colE;  //--colO=Max(colP,colE);
 550         break;
 551       case 2:  // Vertical
 552         linT=(ligP>ligE)? ligP : ligE;  //--linT=Max(ligP, ligE);
 553         colT=colP+colP;
 554         linO=(ligP<ligE)? ligP : ligE;  //--linO=Min(ligP, ligE);
 555         colO=colP-1;
 556         break;
 557       case 3:  // Oblique
 558         linT=(ligP>ligE)? ligP : ligE;  //--linT=Max(ligP, ligE);
 559         colT=colP+colE;
 560         linO=(ligP>ligE)? ligP : ligE;  //--linO=Max(ligP, ligE);
 561         colO=(colP<colE)? colP : colE;  //--colO=Min(colP, colE);
 562         break;
 563       }
 564     }
 565   }
 566
 567   TriCon=(linT-1)*nbdeltaVm2 + colT;
 568
 569   if (linT<1) {
 570     linO=0;
 571     colO=colP+colP-colE;
 572     if (colO<0) {colO=0;linO=1;}
 573     else if (colO>nbdeltaV) {colO=nbdeltaV;linO=1;}
 574     TriCon=0;
 575   }
 576   else if (linT>nbdeltaU) {
 577     linO=nbdeltaU;
 578     colO=colP+colP-colE;
 579     if (colO<0) {colO=0;linO=nbdeltaU-1;}
 580     else if (colO>nbdeltaV) {colO=nbdeltaV;linO=nbdeltaU-1;}
 581     TriCon=0;
 582   }
 583
 584   if (colT<1) {
 585     colO=0;
 586     linO=ligP+ligP-ligE;
 587     if (linO<0) {linO=0;colO=1;}
 588     else if (linO>nbdeltaU) {linO=nbdeltaU;colO=1;}
 589     TriCon=0;
 590   }
 591   else if (colT>nbdeltaV) {
 592     colO=nbdeltaV;
 593     linO=ligP+ligP-ligE;
 594     if (linO<0) {linO=0;colO=nbdeltaV-1;}
 595     else if (linO>nbdeltaU) {linO=nbdeltaU;colO=nbdeltaV-1;}
 596     TriCon=0;
 597   }
 598
 599   OtherP=linO*nbdeltaVp1 + colO+1;
 600
 601   return TriCon;
 602 }
 603
 604
 605 //=======================================================================
 606 //function : PlaneEquation
 607 //purpose  :
 608 //=======================================================================
 609 void IntCurveSurface_Polyhedron::PlaneEquation (const Standard_Integer Triang,
 610                                         gp_XYZ&        NormalVector,
 611                                         Standard_Real& PolarDistance)  const
 612 {
 613   Standard_Integer i1,i2,i3;
 614   Triangle(Triang,i1,i2,i3);
 615
 616   //--  gp_XYZ v1=Point(i2).XYZ()-Point(i1).XYZ();
 617   //--  gp_XYZ v2=Point(i3).XYZ()-Point(i2).XYZ();
 618   //--  gp_XYZ v3=Point(i1).XYZ()-Point(i3).XYZ();
 619
 620   gp_XYZ Pointi1(Point(i1).XYZ());
 621   gp_XYZ Pointi2(Point(i2).XYZ());
 622   gp_XYZ Pointi3(Point(i3).XYZ());
 623
 624
 625   gp_XYZ v1= Pointi2 - Pointi1;
 626   gp_XYZ v2= Pointi3 - Pointi2;
 627   gp_XYZ v3= Pointi1 - Pointi3;
 628
 629   if(v1.SquareModulus()<=LONGUEUR_MINI_EDGE_TRIANGLE) { NormalVector.SetCoord(1.0,0.0,0.0); return; }
 630   if(v2.SquareModulus()<=LONGUEUR_MINI_EDGE_TRIANGLE) { NormalVector.SetCoord(1.0,0.0,0.0); return; }
 631   if(v3.SquareModulus()<=LONGUEUR_MINI_EDGE_TRIANGLE) { NormalVector.SetCoord(1.0,0.0,0.0); return; }
 632
 633   NormalVector= (v1^v2)+(v2^v3)+(v3^v1);
 634   Standard_Real aNormLen = NormalVector.Modulus();
 635   if (aNormLen < gp::Resolution()) {
 636     PolarDistance = 0.;
 637   }
 638   else {
 639     NormalVector.Divide(aNormLen);
 640     PolarDistance = NormalVector * Point(i1).XYZ();
 641   }
 642 }
 643 //=======================================================================
 644 //function : Contain
 645 //purpose  :
 646 //=======================================================================
 647 Standard_Boolean IntCurveSurface_Polyhedron::Contain (const Standard_Integer Triang,
 648                                                       const gp_Pnt& ThePnt) const
 649 {
 650   Standard_Integer i1,i2,i3;
 651   Triangle(Triang,i1,i2,i3);
 652   gp_XYZ Pointi1(Point(i1).XYZ());
 653   gp_XYZ Pointi2(Point(i2).XYZ());
 654   gp_XYZ Pointi3(Point(i3).XYZ());
 655
 656   gp_XYZ v1=(Pointi2-Pointi1)^(ThePnt.XYZ()-Pointi1);
 657   gp_XYZ v2=(Pointi3-Pointi2)^(ThePnt.XYZ()-Pointi2);
 658   gp_XYZ v3=(Pointi1-Pointi3)^(ThePnt.XYZ()-Pointi3);
 659   if (v1*v2 >= 0. && v2*v3 >= 0. && v3*v1>=0.)
 660     return Standard_True;
 661   else
 662     return Standard_False;
 663 }
 664 //=======================================================================
 665 //function : Dump
 666 //purpose  :
 667 //=======================================================================
 668 void IntCurveSurface_Polyhedron::Dump() const
 669 {
 670
 671 }
 672 //=======================================================================
 673 //function : Size
 674 //purpose  :
 675 //=======================================================================
 676 void IntCurveSurface_Polyhedron::Size (Standard_Integer& nbdu,
 677                               Standard_Integer& nbdv) const
 678 {
 679   nbdu=nbdeltaU;
 680   nbdv=nbdeltaV;
 681 }
 682 //=======================================================================
 683 //function : Triangle
 684 //purpose  :
 685 //=======================================================================
 686 void IntCurveSurface_Polyhedron::Triangle (const Standard_Integer Index,
 687                                   Standard_Integer& P1,
 688                                   Standard_Integer& P2,
 689                                   Standard_Integer& P3) const
 690 {
 691   Standard_Integer line=1+((Index-1)/(nbdeltaV*2));
 692   Standard_Integer colon=1+((Index-1)%(nbdeltaV*2));
 693   Standard_Integer colpnt=(colon+1)/2;
 694
 695 // General formula = (line-1)*(nbdeltaV+1)+colpnt
 696
 697 //  Position of P1 = MesXYZ(line,colpnt);
 698   P1= (line-1)*(nbdeltaV+1) + colpnt;
 699
 700 //  Position of P2= MesXYZ(line+1,colpnt+((colon-1)%2));
 701   P2= line*(nbdeltaV+1) + colpnt+((colon-1)%2);
 702
 703 //  Position of P3= MesXYZ(line+(colon%2),colpnt+1);
 704   P3= (line-1+(colon%2))*(nbdeltaV+1) + colpnt + 1;
 705 }
 706 //=======================================================================
 707 //function : Point
 708 //=======================================================================
 709 const gp_Pnt& IntCurveSurface_Polyhedron::Point(const Standard_Integer Index
 710                                                 ,Standard_Real& U
 711                                                 ,Standard_Real& V) const
 712 {
 713 #if CHECKBOUNDS
 714   if(Index<0 || Index>((nbdeltaU+1)*(nbdeltaV+1))) {
 715     printf("\n Erreur IntCurveSurface_Polyhedron::Parameters\n");
 716   }
 717 #endif
 718   gp_Pnt *CMyPnts     = (gp_Pnt *)C_MyPnts;
 719   Standard_Real *CMyU = (Standard_Real *)C_MyU;
 720   Standard_Real *CMyV = (Standard_Real *)C_MyV;
 721   U=CMyU[Index];
 722   V=CMyV[Index];
 723   return CMyPnts[Index];
 724 }
 725 //=======================================================================
 726 //function : Point
 727 //=======================================================================
 728 const gp_Pnt& IntCurveSurface_Polyhedron::Point(const Standard_Integer Index) const {
 729 #if CHECKBOUNDS
 730   if(Index<0 || Index>((nbdeltaU+1)*(nbdeltaV+1))) {
 731     printf("\n Erreur IntCurveSurface_Polyhedron::Parameters\n");
 732   }
 733 #endif
 734
 735   gp_Pnt *CMyPnts     = (gp_Pnt *)C_MyPnts;
 736   return CMyPnts[Index];
 737 }
 738 //=======================================================================
 739 //function : Point
 740 //=======================================================================
 741 //void IntCurveSurface_Polyhedron::Point (const gp_Pnt& p,
 742 //                                      const Standard_Integer lig,
 743 //                                      const Standard_Integer col,
 744 //                                      const Standard_Real u,
 745 //                                      const Standard_Real v)
 746 void IntCurveSurface_Polyhedron::Point (const gp_Pnt& ,
 747                                         const Standard_Integer ,
 748                                         const Standard_Integer ,
 749                                         const Standard_Real ,
 750                                         const Standard_Real )
 751 {
 752   printf("\n IntCurveSurface_Polyhedron::Point : Ne dois pas etre appelle\n");
 753 }
 754 //=======================================================================
 755 //function : Point
 756 //=======================================================================
 757 void IntCurveSurface_Polyhedron::Point (const Standard_Integer Index,gp_Pnt& P) const
 758 {
 759 #if CHECKBOUNDS
 760   if(Index<0 || Index>((nbdeltaU+1)*(nbdeltaV+1))) {
 761     printf("\n Erreur IntCurveSurface_Polyhedron::Parameters\n");
 762   }
 763 #endif
 764
 765   gp_Pnt *CMyPnts     = (gp_Pnt *)C_MyPnts;
 766   P = CMyPnts[Index];
 767 }
 768
 769 //  Modified by Sergey KHROMOV - Fri Dec  7 10:12:47 2001 Begin
 770
 771 //=======================================================================
 772 //function : IsOnBound
 773 //purpose  : This method returns true if the edge based on points with
 774 //           indices Index1 and Index2 represents a boundary edge.
 775 //=======================================================================
 776
 777 Standard_Boolean IntCurveSurface_Polyhedron::IsOnBound
 778                                  (const Standard_Integer Index1,
 779                                   const Standard_Integer Index2) const
 780 {
 781 #if CHECKBOUNDS
 782   if(Index1<0 || Index1>((nbdeltaU+1)*(nbdeltaV+1))) {
 783     printf("\n Erreur IntCurveSurface_Polyhedron::IsOnBound\n");
 784   }
 785   if(Index2<0 || Index2>((nbdeltaU+1)*(nbdeltaV+1))) {
 786     printf("\n Erreur IntCurveSurface_Polyhedron::IsOnBound\n");
 787   }
 788 #endif
 789   Standard_Boolean *CMyIsOnBounds = (Standard_Boolean *)C_MyIsOnBounds;
 790   Standard_Integer  aDiff         = Abs(Index1 - Index2);
 791   Standard_Integer  i;
 792
 793 // Check if points are neighbour ones.
 794   if (aDiff != 1 && aDiff != nbdeltaV + 1)
 795     return Standard_False;
 796
 797   for (i = 0; i <= nbdeltaU; i++) {
 798     if ((Index1 == 1 + i*(nbdeltaV + 1)) && (Index2 == Index1 - 1))
 799       return Standard_False;
 800
 801     if ((Index1 == (1 + i)*(nbdeltaV + 1)) && (Index2 == Index1 + 1))
 802       return Standard_False;
 803   }
 804
 805   return (CMyIsOnBounds[Index1] && CMyIsOnBounds[Index2]);
 806 }
 807
 808 //=======================================================================
 809 //function : ComputeBorderDeflection
 810 //purpose  : This method computes and returns a deflection of isoline
 811 //           of given parameter on Surface.
 812 //=======================================================================
 813
 814 Standard_Real IntCurveSurface_Polyhedron::ComputeBorderDeflection
 815                               (const ThePSurface      &Surface,
 816                                const Standard_Real     Parameter,
 817                                const Standard_Real     PMin,
 818                                const Standard_Real     PMax,
 819                                const Standard_Boolean  isUIso) const
 820 {
 821   Standard_Integer aNbSamples;
 822   Standard_Integer i;
 823
 824   if (isUIso)
 825     aNbSamples = nbdeltaV;
 826   else
 827     aNbSamples = nbdeltaU;
 828
 829   Standard_Real aDelta      = (PMax - PMin)/aNbSamples;
 830   Standard_Real aPar        = PMin;
 831   Standard_Real aDeflection = RealFirst();
 832   gp_XYZ        aP1;
 833   gp_XYZ        aP2;
 834   gp_XYZ        aPMid;
 835   gp_XYZ        aPParMid;
 836
 837   for (i = 0; i <= aNbSamples; i++, aPar+= aDelta) {
 838     if (isUIso) {
 839       aP1      = ThePSurfaceTool::Value(Surface, Parameter, aPar).XYZ();
 840       aP2      = ThePSurfaceTool::Value(Surface, Parameter, aPar + aDelta).XYZ();
 841       aPParMid = ThePSurfaceTool::Value(Surface, Parameter, aPar + aDelta/2.).XYZ();
 842     } else {
 843       aP1      = ThePSurfaceTool::Value(Surface, aPar,             Parameter).XYZ();
 844       aP2      = ThePSurfaceTool::Value(Surface, aPar + aDelta,    Parameter).XYZ();
 845       aPParMid = ThePSurfaceTool::Value(Surface, aPar + aDelta/2., Parameter).XYZ();
 846     }
 847     aPMid = (aP2 + aP1)/2.;
 848
 849     Standard_Real aDist = (aPMid - aPParMid).Modulus();
 850
 851     if (aDist > aDeflection)
 852       aDeflection = aDist;
 853   }
 854
 855   return aDeflection;
 856 }
 857
 858 //  Modified by Sergey KHROMOV - Fri Dec  7 11:21:52 2001 End