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