src/Voxel/Voxel_FastConverter.cxx

   1 // Created on: 2008-05-30
   2 // Created by: Vladislav ROMASHKO
   3 // Copyright (c) 2008-2012 OPEN CASCADE SAS
   4 //
   5 // The content of this file is subject to the Open CASCADE Technology Public
   6 // License Version 6.5 (the "License"). You may not use the content of this file
   7 // except in compliance with the License. Please obtain a copy of the License
   8 // at http://www.opencascade.org and read it completely before using this file.
   9 //
  10 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
  11 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
  12 //
  13 // The Original Code and all software distributed under the License is
  14 // distributed on an "AS IS" basis, without warranty of any kind, and the
  15 // Initial Developer hereby disclaims all such warranties, including without
  16 // limitation, any warranties of merchantability, fitness for a particular
  17 // purpose or non-infringement. Please see the License for the specific terms
  18 // and conditions governing the rights and limitations under the License.
  19
  20
  21 #include <Voxel_FastConverter.ixx>
  22
  23 #include <Bnd_Box.hxx>
  24 #include <BRep_Tool.hxx>
  25 #include <BRepBndLib.hxx>
  26 #include <BRepMesh.hxx>
  27
  28 #include <TopoDS.hxx>
  29 #include <TopoDS_Face.hxx>
  30 #include <TopExp_Explorer.hxx>
  31
  32 #include <gp_Lin2d.hxx>
  33 #include <gce_MakePln.hxx>
  34
  35 #include <ElSLib.hxx>
  36 #include <Poly_Triangulation.hxx>
  37 #include <IntAna2d_AnaIntersection.hxx>
  38
  39 // Printing the progress in stdout.
  40 //#define CONV_DUMP
  41
  42 Voxel_FastConverter::Voxel_FastConverter(const TopoDS_Shape&    shape,
  43                                          Voxel_ROctBoolDS&      voxels,
  44                                          const Standard_Real    deflection,
  45                                          const Standard_Integer nbx,
  46                                          const Standard_Integer nby,
  47                                          const Standard_Integer nbz,
  48                                          const Standard_Integer nbthreads,
  49                                          const Standard_Boolean useExistingTriangulation)
  50 :myShape(shape),myVoxels(&voxels),
  51  myDeflection(deflection),
  52  myNbX(nbx),myNbY(nby),myNbZ(nbz),
  53  myNbThreads(nbthreads),myIsBool(2),
  54  myNbTriangles(0),
  55  myUseExistingTriangulation(useExistingTriangulation)
  56 {
  57   Init();
  58 }
  59
  60 Voxel_FastConverter::Voxel_FastConverter(const TopoDS_Shape&    shape,
  61                                          Voxel_BoolDS&          voxels,
  62                                          const Standard_Real    deflection,
  63                                          const Standard_Integer nbx,
  64                                          const Standard_Integer nby,
  65                                          const Standard_Integer nbz,
  66                                          const Standard_Integer nbthreads,
  67                                          const Standard_Boolean useExistingTriangulation)
  68 :myShape(shape),myVoxels(&voxels),
  69  myDeflection(deflection),
  70  myNbX(nbx),myNbY(nby),myNbZ(nbz),
  71  myNbThreads(nbthreads),myIsBool(1),
  72  myNbTriangles(0),
  73  myUseExistingTriangulation(useExistingTriangulation)
  74 {
  75   Init();
  76 }
  77
  78 Voxel_FastConverter::Voxel_FastConverter(const TopoDS_Shape&    shape,
  79                                          Voxel_ColorDS&         voxels,
  80                                          const Standard_Real    deflection,
  81                                          const Standard_Integer nbx,
  82                                          const Standard_Integer nby,
  83                                          const Standard_Integer nbz,
  84                                          const Standard_Integer nbthreads,
  85                                          const Standard_Boolean useExistingTriangulation)
  86 :myShape(shape),myVoxels(&voxels),
  87  myDeflection(deflection),
  88  myNbX(nbx),myNbY(nby),myNbZ(nbz),
  89  myNbThreads(nbthreads),myIsBool(0),
  90  myNbTriangles(0),
  91  myUseExistingTriangulation(useExistingTriangulation)
  92 {
  93   Init();
  94 }
  95
  96 void Voxel_FastConverter::Init()
  97 {
  98   if (myShape.IsNull())
  99     return;
 100   if (myNbThreads < 1)
 101     return;
 102
 103   // Check number of splits.
 104   Voxel_DS* voxels = (Voxel_DS*) myVoxels;
 105   if (voxels->GetNbX() != myNbX || voxels->GetNbY() != myNbY || voxels->GetNbZ() != myNbZ)
 106   {
 107     // Compute boundary box of the shape
 108     Bnd_Box box;
 109     BRepBndLib::Add(myShape, box);
 110
 111     // Define the voxel model by means of the boundary box of shape
 112     Standard_Real xmin, ymin, zmin, xmax, ymax, zmax;
 113     box.Get(xmin, ymin, zmin, xmax, ymax, zmax);
 114
 115     // Initialize the voxels.
 116     if (myIsBool == 2)
 117       ((Voxel_ROctBoolDS*) voxels)->Init(xmin, ymin, zmin, xmax - xmin, ymax - ymin, zmax - zmin, myNbX, myNbY, myNbZ);
 118     else if (myIsBool == 1)
 119       ((Voxel_BoolDS*) voxels)->Init(xmin, ymin, zmin, xmax - xmin, ymax - ymin, zmax - zmin, myNbX, myNbY, myNbZ);
 120     else if (myIsBool == 0)
 121       ((Voxel_ColorDS*) voxels)->Init(xmin, ymin, zmin, xmax - xmin, ymax - ymin, zmax - zmin, myNbX, myNbY, myNbZ);
 122   }
 123
 124   // Check presence of triangulation.
 125   TopLoc_Location L;
 126   Standard_Boolean triangulate = Standard_False;
 127   TopExp_Explorer expl(myShape, TopAbs_FACE);
 128   if(myUseExistingTriangulation == Standard_False)
 129   {
 130     for (; expl.More(); expl.Next())
 131     {
 132       const TopoDS_Face & F = TopoDS::Face(expl.Current());
 133       Handle(Poly_Triangulation) T = BRep_Tool::Triangulation(F, L);
 134       if (T.IsNull() || (T->Deflection() > myDeflection))
 135       {
 136         triangulate = Standard_True;
 137         break;
 138       }
 139     }
 140   }
 141
 142   // Re-create the triangulation.
 143   if (triangulate)
 144   {
 145     BRepMesh::Mesh(myShape, myDeflection);
 146   }
 147
 148   // Compute the number of triangles.
 149   myNbTriangles = 0;
 150   expl.Init(myShape, TopAbs_FACE);
 151   for (; expl.More(); expl.Next())
 152   {
 153     const TopoDS_Face & F = TopoDS::Face(expl.Current());
 154     Handle(Poly_Triangulation) T = BRep_Tool::Triangulation(F, L);
 155     if (T.IsNull() == Standard_False)
 156       myNbTriangles += T->NbTriangles();
 157   }
 158 }
 159
 160 // Destructor
 161 void Voxel_FastConverter::Destroy()
 162 {
 163
 164 }
 165
 166 Standard_Boolean Voxel_FastConverter::Convert(Standard_Integer&      progress,
 167                                               const Standard_Integer ithread)
 168 {
 169   if (ithread == 1)
 170     progress = 0;
 171 #ifdef CONV_DUMP
 172   if (ithread == 1)
 173     printf("Progress = %d   \r", progress);
 174 #endif
 175
 176   if (myNbX <= 0 || myNbY <= 0 || myNbZ <= 0)
 177     return Standard_False;
 178
 179   if(myNbTriangles == 0)
 180     return Standard_False;
 181
 182   // Half of diagonal of a voxel
 183   Voxel_DS* ds = (Voxel_DS*) myVoxels;
 184   Standard_Real dx = ds->GetXLen() / (Standard_Real) ds->GetNbX(),
 185                 dy = ds->GetYLen() / (Standard_Real) ds->GetNbY(),
 186                 dz = ds->GetZLen() / (Standard_Real) ds->GetNbZ();
 187   Standard_Real hdiagonal = sqrt(dx * dx + dy * dy + dz * dz);
 188   hdiagonal /= 2.0;
 189
 190   // Compute the scope of triangles for current thread
 191   Standard_Integer start_thread_triangle = 1, end_thread_triangle = myNbTriangles, ithread_triangle = 0;
 192   start_thread_triangle = (ithread - 1) * (myNbTriangles / myNbThreads) + 1;
 193   end_thread_triangle   = (ithread - 0) * (myNbTriangles / myNbThreads);
 194
 195   // Convert
 196   TopLoc_Location L;
 197   Standard_Integer iprogress = 0, prev_progress = 0;
 198   Standard_Integer n1, n2, n3;
 199   Standard_Integer ixmin, iymin, izmin, ixmax, iymax, izmax;
 200   Standard_Real xmin, ymin, zmin, xmax, ymax, zmax;
 201   TopExp_Explorer expl(myShape, TopAbs_FACE);
 202   for (; expl.More(); expl.Next())
 203   {
 204     const TopoDS_Face & F = TopoDS::Face(expl.Current());
 205     Handle(Poly_Triangulation) T = BRep_Tool::Triangulation(F, L);
 206     if (T.IsNull())
 207       continue;
 208
 209     gp_Trsf trsf;
 210     Standard_Boolean transform = !L.IsIdentity();
 211     if (transform)
 212       trsf = L.Transformation();
 213
 214     const TColgp_Array1OfPnt& nodes = T->Nodes();
 215     const Poly_Array1OfTriangle& triangles = T->Triangles();
 216     Standard_Integer itriangle = triangles.Lower(), nb_triangles = triangles.Upper();
 217     for (; itriangle <= nb_triangles; itriangle++)
 218     {
 219       ithread_triangle++;
 220       if (ithread_triangle < start_thread_triangle || ithread_triangle > end_thread_triangle)
 221         continue;
 222
 223       const Poly_Triangle& t = triangles.Value(itriangle);
 224       t.Get(n1, n2, n3);
 225       gp_Pnt p1 = nodes.Value(n1);
 226       gp_Pnt p2 = nodes.Value(n2);
 227       gp_Pnt p3 = nodes.Value(n3);
 228       if (transform)
 229       {
 230         p1.Transform(trsf);
 231         p2.Transform(trsf);
 232         p3.Transform(trsf);
 233       }
 234
 235       // Get boundary box of the triangle
 236       GetBndBox(p1, p2, p3, xmin, ymin, zmin, xmax, ymax, zmax);
 237
 238       // Find the range of voxels inside the boudary box of the triangle.
 239       if (!ds->GetVoxel(xmin, ymin, zmin, ixmin, iymin, izmin))
 240         continue;
 241       if (!ds->GetVoxel(xmax, ymax, zmax, ixmax, iymax, izmax))
 242         continue;
 243
 244       // Refuse voxels for whom distance from their center to plane of triangle is greater than half of diagonal.
 245       // Make a line from center of each voxel to the center of triangle and
 246       // compute intersection of the line with sides of triangle.
 247       // Refuse the voxel in case of intersection.
 248       gce_MakePln mkPlane(p1, p2, p3);
 249       if (!mkPlane.IsDone())
 250         continue;
 251       gp_Pln plane = mkPlane.Value();
 252       ComputeVoxelsNearTriangle(plane, p1, p2, p3, hdiagonal, ixmin, iymin, izmin, ixmax, iymax, izmax);
 253
 254       // Progress
 255       if (ithread == 1)
 256       {
 257         iprogress++;
 258         progress = (Standard_Integer) ( (Standard_Real) iprogress / (Standard_Real) myNbTriangles * 100.0 );
 259       }
 260 #ifdef CONV_DUMP
 261       if (ithread == 1 && prev_progress != progress)
 262       {
 263         printf("Progress = %d  \r", progress);
 264         prev_progress = progress;
 265       }
 266 #endif
 267
 268     } // iteration of triangles
 269   } // iteration of faces
 270
 271   if (ithread == 1)
 272     progress = 100;
 273 #ifdef CONV_DUMP
 274   if (ithread == 1)
 275     printf("Progress = %d  \r", progress);
 276 #endif
 277   return Standard_True;
 278 }
 279
 280 Standard_Boolean Voxel_FastConverter::FillInVolume(const Standard_Byte inner,
 281                                                    const Standard_Integer ithread)
 282 {
 283   Voxel_DS* ds = (Voxel_DS*) myVoxels;
 284   Standard_Integer ix, iy, iz, nbx = ds->GetNbX(), nby = ds->GetNbY(), nbz = ds->GetNbZ();
 285   Standard_Boolean prev_surface, surface, volume;
 286
 287   if (inner)
 288   {
 289     // Fill-in internal voxels by the value "inner"
 290     for (ix = 0; ix < nbx; ix++)
 291     {
 292       for (iy = 0; iy < nby; iy++)
 293       {
 294         // Check existence of volume.
 295         volume = Standard_False;
 296         surface = Standard_False;
 297         prev_surface = Standard_False;
 298         for (iz = 0; iz < nbz; iz++)
 299         {
 300           surface = (myIsBool == 1) ?
 301             ((Voxel_BoolDS*)myVoxels)->Get(ix, iy, iz) == Standard_True :
 302               ((Voxel_ColorDS*)myVoxels)->Get(ix, iy, iz) > 0;
 303           if (prev_surface && !surface)
 304           {
 305             volume = !volume;
 306           }
 307           prev_surface = surface;
 308         }
 309         if (volume)
 310           continue;
 311
 312         // Fill-in the volume.
 313         volume = Standard_False;
 314         surface = Standard_False;
 315         prev_surface = Standard_False;
 316         for (iz = 0; iz < nbz; iz++)
 317         {
 318           surface = (myIsBool == 1) ?
 319             ((Voxel_BoolDS*)myVoxels)->Get(ix, iy, iz) == Standard_True :
 320               ((Voxel_ColorDS*)myVoxels)->Get(ix, iy, iz) > 0;
 321           if (prev_surface && !surface)
 322           {
 323             volume = !volume;
 324           }
 325           if (volume && !surface)
 326           {
 327             (myIsBool == 1) ? ((Voxel_BoolDS*)myVoxels)->Set(ix, iy, iz, inner) :
 328               ((Voxel_ColorDS*)myVoxels)->Set(ix, iy, iz, inner);
 329           }
 330           prev_surface = surface;
 331         }
 332       }
 333     }
 334   }
 335   else
 336   {
 337     // Set value of interbal voxels to 0 ("inner" = 0)
 338     Standard_Boolean next_surface;
 339     for (ix = 0; ix < nbx; ix++)
 340     {
 341       for (iy = 0; iy < nby; iy++)
 342       {
 343         volume = Standard_False;
 344         surface = Standard_False;
 345         prev_surface = Standard_False;
 346         next_surface = Standard_False;
 347         for (iz = 0; iz < nbz; iz++)
 348         {
 349           surface = (myIsBool == 1) ?
 350             ((Voxel_BoolDS*)myVoxels)->Get(ix, iy, iz) == Standard_True :
 351               ((Voxel_ColorDS*)myVoxels)->Get(ix, iy, iz) > 0;
 352           if (prev_surface != surface)
 353           {
 354             volume = !volume;
 355           }
 356           if (volume && iz + 1 < nbz)
 357           {
 358             next_surface = (myIsBool == 1) ?
 359               ((Voxel_BoolDS*)myVoxels)->Get(ix, iy, iz + 1) == Standard_True :
 360               ((Voxel_ColorDS*)myVoxels)->Get(ix, iy, iz + 1) > 0;
 361           }
 362           if (volume && prev_surface == surface && next_surface)
 363           {
 364             (myIsBool == 1) ? ((Voxel_BoolDS*)myVoxels)->Set(ix, iy, iz, inner) :
 365               ((Voxel_ColorDS*)myVoxels)->Set(ix, iy, iz, inner);
 366           }
 367           prev_surface = surface;
 368         }
 369       }
 370     }
 371   }
 372
 373   return Standard_True;
 374 }
 375
 376 void Voxel_FastConverter::GetBndBox(const gp_Pnt&  p1,
 377                                     const gp_Pnt&  p2,
 378                                     const gp_Pnt&  p3,
 379                                     Standard_Real& xmin,
 380                                     Standard_Real& ymin,
 381                                     Standard_Real& zmin,
 382                                     Standard_Real& xmax,
 383                                     Standard_Real& ymax,
 384                                     Standard_Real& zmax) const
 385 {
 386   // P1:
 387   xmin = p1.X();
 388   ymin = p1.Y();
 389   zmin = p1.Z();
 390   xmax = p1.X();
 391   ymax = p1.Y();
 392   zmax = p1.Z();
 393   // P2:
 394   if (xmin > p2.X())
 395     xmin = p2.X();
 396   if (ymin > p2.Y())
 397     ymin = p2.Y();
 398   if (zmin > p2.Z())
 399     zmin = p2.Z();
 400   if (xmax < p2.X())
 401     xmax = p2.X();
 402   if (ymax < p2.Y())
 403     ymax = p2.Y();
 404   if (zmax < p2.Z())
 405     zmax = p2.Z();
 406   // P3:
 407   if (xmin > p3.X())
 408     xmin = p3.X();
 409   if (ymin > p3.Y())
 410     ymin = p3.Y();
 411   if (zmin > p3.Z())
 412     zmin = p3.Z();
 413   if (xmax < p3.X())
 414     xmax = p3.X();
 415   if (ymax < p3.Y())
 416     ymax = p3.Y();
 417   if (zmax < p3.Z())
 418     zmax = p3.Z();
 419 }
 420
 421 // This method is copied from Voxel_ShapeIntersector.cxx
 422 static Standard_Boolean mayIntersect(const gp_Pnt2d& p11, const gp_Pnt2d& p12,
 423                                      const gp_Pnt2d& p21, const gp_Pnt2d& p22)
 424 {
 425     if (p11.X() > p21.X() && p11.X() > p22.X() && p12.X() > p21.X() && p12.X() > p22.X())
 426         return Standard_False;
 427     if (p11.X() < p21.X() && p11.X() < p22.X() && p12.X() < p21.X() && p12.X() < p22.X())
 428         return Standard_False;
 429     if (p11.Y() > p21.Y() && p11.Y() > p22.Y() && p12.Y() > p21.Y() && p12.Y() > p22.Y())
 430         return Standard_False;
 431     if (p11.Y() < p21.Y() && p11.Y() < p22.Y() && p12.Y() < p21.Y() && p12.Y() < p22.Y())
 432         return Standard_False;
 433     return Standard_True;
 434 }
 435
 436 void Voxel_FastConverter::ComputeVoxelsNearTriangle(const gp_Pln&          plane,
 437                                                     const gp_Pnt&          p1,
 438                                                     const gp_Pnt&          p2,
 439                                                     const gp_Pnt&          p3,
 440                                                     const Standard_Real    hdiagonal,
 441                                                     const Standard_Integer ixmin,
 442                                                     const Standard_Integer iymin,
 443                                                     const Standard_Integer izmin,
 444                                                     const Standard_Integer ixmax,
 445                                                     const Standard_Integer iymax,
 446                                                     const Standard_Integer izmax) const
 447 {
 448   gp_Pnt pc;
 449   Standard_Real xc, yc, zc, uc, vc, u1, v1, u2, v2, u3, v3;
 450   Standard_Integer ix, iy, iz;
 451   IntAna2d_AnaIntersection intersector2d;
 452
 453   // Project points of triangle onto the plane
 454   ElSLib::Parameters(plane, p1, u1, v1);
 455   ElSLib::Parameters(plane, p2, u2, v2);
 456   ElSLib::Parameters(plane, p3, u3, v3);
 457
 458   // Make lines of triangle
 459   gp_Pnt2d p2d1(u1, v1), p2d2(u2, v2), p2d3(u3, v3), p2dt((u1+u2+u3)/3.0,(v1+v2+v3)/3.0), p2dc;
 460   gp_Vec2d v2d12(p2d1, p2d2), v2d23(p2d2, p2d3), v2d31(p2d3, p2d1);
 461   gp_Lin2d L1(p2d1, v2d12), L2(p2d2, v2d23), L3(p2d3, v2d31), Lv;
 462   Standard_Real d1 = p2d1.Distance(p2d2) - Precision::Confusion(),
 463                 d2 = p2d2.Distance(p2d3) - Precision::Confusion(),
 464                 d3 = p2d3.Distance(p2d1) - Precision::Confusion(), dv;
 465
 466   Voxel_DS* ds = (Voxel_DS*) myVoxels;
 467   for (ix = ixmin; ix <= ixmax; ix++)
 468   {
 469     for (iy = iymin; iy <= iymax; iy++)
 470     {
 471       for (iz = izmin; iz <= izmax; iz++)
 472       {
 473         ds->GetCenter(ix, iy, iz, xc, yc, zc);
 474         pc.SetCoord(xc, yc, zc);
 475         if (plane.Distance(pc) < hdiagonal)
 476         {
 477           ElSLib::Parameters(plane, pc, uc, vc);
 478           p2dc.SetCoord(uc, vc);
 479
 480           gp_Vec2d v2dct(p2dc, p2dt);
 481           dv = v2dct.Magnitude() - Precision::Confusion();
 482           Lv.SetLocation(p2dc);
 483           Lv.SetDirection(v2dct);
 484
 485           // Side 1:
 486           if (mayIntersect(p2d1, p2d2, p2dc, p2dt))
 487           {
 488             intersector2d.Perform(Lv, L1);
 489             if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
 490             {
 491               const IntAna2d_IntPoint& i2d = intersector2d.Point(1);
 492               Standard_Real param1 = i2d.ParamOnFirst();
 493               Standard_Real param2 = i2d.ParamOnSecond();
 494               if (param1 > Precision::Confusion() && param1 < dv &&
 495                   param2 > Precision::Confusion() && param2 < d1)
 496               {
 497                 continue;
 498               }
 499             }
 500           }
 501
 502           // Side 2:
 503           if (mayIntersect(p2d2, p2d3, p2dc, p2dt))
 504           {
 505             intersector2d.Perform(Lv, L2);
 506             if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
 507             {
 508               const IntAna2d_IntPoint& i2d = intersector2d.Point(1);
 509               Standard_Real param1 = i2d.ParamOnFirst();
 510               Standard_Real param2 = i2d.ParamOnSecond();
 511               if (param1 > Precision::Confusion() && param1 < dv &&
 512                   param2 > Precision::Confusion() && param2 < d2)
 513               {
 514                 continue;
 515               }
 516             }
 517           }
 518
 519           // Side 3:
 520           if (mayIntersect(p2d3, p2d1, p2dc, p2dt))
 521           {
 522             intersector2d.Perform(Lv, L3);
 523             if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
 524             {
 525               const IntAna2d_IntPoint& i2d = intersector2d.Point(1);
 526               Standard_Real param1 = i2d.ParamOnFirst();
 527               Standard_Real param2 = i2d.ParamOnSecond();
 528               if (param1 > Precision::Confusion() && param1 < dv &&
 529                   param2 > Precision::Confusion() && param2 < d3)
 530               {
 531                 continue;
 532               }
 533             }
 534           }
 535
 536           // Set positive value to this voxel:
 537           switch (myIsBool)
 538           {
 539             case 0:
 540               ((Voxel_ColorDS*) myVoxels)->Set(ix, iy, iz, 15);
 541               break;
 542             case 1:
 543               ((Voxel_BoolDS*) myVoxels)->Set(ix, iy, iz, Standard_True);
 544               break;
 545             case 2:
 546             {
 547               //((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, Standard_True);
 548
 549               // Check intersection between the triangle & sub-voxels of the voxel.
 550               Standard_Real hdiagonal2 = hdiagonal / 2.0, hdiagonal4 = hdiagonal / 4.0;
 551               for (Standard_Integer i = 0; i < 8; i++)
 552               {
 553                 ((Voxel_ROctBoolDS*) myVoxels)->GetCenter(ix, iy, iz, i, xc, yc, zc);
 554                 pc.SetCoord(xc, yc, zc);
 555                 if (plane.Distance(pc) < hdiagonal2)
 556                 {
 557                   ElSLib::Parameters(plane, pc, uc, vc);
 558                   p2dc.SetCoord(uc, vc);
 559
 560                   gp_Vec2d v2dct(p2dc, p2dt);
 561                   dv = v2dct.Magnitude() - Precision::Confusion();
 562                   Lv.SetLocation(p2dc);
 563                   Lv.SetDirection(v2dct);
 564
 565                   // Side 1:
 566                   if (mayIntersect(p2d1, p2d2, p2dc, p2dt))
 567                   {
 568                     intersector2d.Perform(Lv, L1);
 569                     if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
 570                     {
 571                       const IntAna2d_IntPoint& i2d = intersector2d.Point(1);
 572                       Standard_Real param1 = i2d.ParamOnFirst();
 573                       Standard_Real param2 = i2d.ParamOnSecond();
 574                       if (param1 > Precision::Confusion() && param1 < dv &&
 575                           param2 > Precision::Confusion() && param2 < d1)
 576                       {
 577                         continue;
 578                       }
 579                     }
 580                   }
 581
 582                   // Side 2:
 583                   if (mayIntersect(p2d2, p2d3, p2dc, p2dt))
 584                   {
 585                     intersector2d.Perform(Lv, L2);
 586                     if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
 587                     {
 588                       const IntAna2d_IntPoint& i2d = intersector2d.Point(1);
 589                       Standard_Real param1 = i2d.ParamOnFirst();
 590                       Standard_Real param2 = i2d.ParamOnSecond();
 591                       if (param1 > Precision::Confusion() && param1 < dv &&
 592                           param2 > Precision::Confusion() && param2 < d2)
 593                       {
 594                         continue;
 595                       }
 596                     }
 597                   }
 598
 599                   // Side 3:
 600                   if (mayIntersect(p2d3, p2d1, p2dc, p2dt))
 601                   {
 602                     intersector2d.Perform(Lv, L3);
 603                     if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
 604                     {
 605                       const IntAna2d_IntPoint& i2d = intersector2d.Point(1);
 606                       Standard_Real param1 = i2d.ParamOnFirst();
 607                       Standard_Real param2 = i2d.ParamOnSecond();
 608                       if (param1 > Precision::Confusion() && param1 < dv &&
 609                           param2 > Precision::Confusion() && param2 < d3)
 610                       {
 611                         continue;
 612                       }
 613                     }
 614                   }
 615
 616                   //((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, i, Standard_True);
 617
 618                   // Check intersection between the triangle & sub-voxels of the sub-voxel.
 619                   for (Standard_Integer j = 0; j < 8; j++)
 620                   {
 621                     ((Voxel_ROctBoolDS*) myVoxels)->GetCenter(ix, iy, iz, i, j, xc, yc, zc);
 622                     pc.SetCoord(xc, yc, zc);
 623                     if (plane.Distance(pc) < hdiagonal4)
 624                     {
 625                       ElSLib::Parameters(plane, pc, uc, vc);
 626                       p2dc.SetCoord(uc, vc);
 627
 628                       gp_Vec2d v2dct(p2dc, p2dt);
 629                       dv = v2dct.Magnitude() - Precision::Confusion();
 630                       Lv.SetLocation(p2dc);
 631                       Lv.SetDirection(v2dct);
 632
 633                       // Side 1:
 634                       if (mayIntersect(p2d1, p2d2, p2dc, p2dt))
 635                       {
 636                         intersector2d.Perform(Lv, L1);
 637                         if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
 638                         {
 639                           const IntAna2d_IntPoint& i2d = intersector2d.Point(1);
 640                           Standard_Real param1 = i2d.ParamOnFirst();
 641                           Standard_Real param2 = i2d.ParamOnSecond();
 642                           if (param1 > Precision::Confusion() && param1 < dv &&
 643                               param2 > Precision::Confusion() && param2 < d1)
 644                           {
 645                             continue;
 646                           }
 647                         }
 648                       }
 649
 650                       // Side 2:
 651                       if (mayIntersect(p2d2, p2d3, p2dc, p2dt))
 652                       {
 653                         intersector2d.Perform(Lv, L2);
 654                         if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
 655                         {
 656                           const IntAna2d_IntPoint& i2d = intersector2d.Point(1);
 657                           Standard_Real param1 = i2d.ParamOnFirst();
 658                           Standard_Real param2 = i2d.ParamOnSecond();
 659                           if (param1 > Precision::Confusion() && param1 < dv &&
 660                               param2 > Precision::Confusion() && param2 < d2)
 661                           {
 662                             continue;
 663                           }
 664                         }
 665                       }
 666
 667                       // Side 3:
 668                       if (mayIntersect(p2d3, p2d1, p2dc, p2dt))
 669                       {
 670                         intersector2d.Perform(Lv, L3);
 671                         if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
 672                         {
 673                           const IntAna2d_IntPoint& i2d = intersector2d.Point(1);
 674                           Standard_Real param1 = i2d.ParamOnFirst();
 675                           Standard_Real param2 = i2d.ParamOnSecond();
 676                           if (param1 > Precision::Confusion() && param1 < dv &&
 677                               param2 > Precision::Confusion() && param2 < d3)
 678                           {
 679                             continue;
 680                           }
 681                         }
 682                       }
 683
 684                       ((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, i, j, Standard_True);
 685                     }
 686                   } // End of "Check level 2".
 687
 688                 }
 689               } // End of "Check level 1".
 690
 691               break;
 692             }
 693           }
 694         }
 695       }
 696     }
 697   }
 698 }