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