src/Voxel/Voxel_FastConverter.cxx

   1 // Created on: 2008-05-30
   2 // Created by: Vladislav ROMASHKO
   3 // Copyright (c) 2008-2014 OPEN CASCADE SAS
   4 //
   5 // This file is part of Open CASCADE Technology software library.
   6 //
   7 // This library is free software; you can redistribute it and/or modify it under
   8 // the terms of the GNU Lesser General Public License version 2.1 as published
   9 // by the Free Software Foundation, with special exception defined in the file
  10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  11 // distribution for complete text of the license and disclaimer of any warranty.
  12 //
  13 // Alternatively, this file may be used under the terms of Open CASCADE
  14 // commercial license or contractual agreement.
  15
  16
  17 #include <Bnd_Box.hxx>
  18 #include <BRep_Tool.hxx>
  19 #include <BRepBndLib.hxx>
  20 #include <BRepClass3d_SolidClassifier.hxx>
  21 #include <BRepMesh_IncrementalMesh.hxx>
  22 #include <ElSLib.hxx>
  23 #include <gce_MakePln.hxx>
  24 #include <gp_Lin2d.hxx>
  25 #include <gp_Pln.hxx>
  26 #include <gp_Pnt.hxx>
  27 #include <IntAna2d_AnaIntersection.hxx>
  28 #include <Poly_Triangulation.hxx>
  29 #include <TopExp_Explorer.hxx>
  30 #include <TopoDS.hxx>
  31 #include <TopoDS_Face.hxx>
  32 #include <TopoDS_Shape.hxx>
  33 #include <Voxel_BoolDS.hxx>
  34 #include <Voxel_ColorDS.hxx>
  35 #include <Voxel_FastConverter.hxx>
  36 #include <Voxel_ROctBoolDS.hxx>
  37
  38 // Printing the progress in stdout.
  39 //#define CONV_DUMP
  40 Voxel_FastConverter::Voxel_FastConverter(const TopoDS_Shape&    shape,
  41                                          Voxel_ROctBoolDS&      voxels,
  42                                          const Standard_Real    deflection,
  43                                          const Standard_Integer nbx,
  44                                          const Standard_Integer nby,
  45                                          const Standard_Integer nbz,
  46                                          const Standard_Integer nbthreads,
  47                                          const Standard_Boolean useExistingTriangulation)
  48 :myShape(shape),myVoxels(&voxels),
  49  myDeflection(deflection),
  50  myIsBool(2),
  51  myNbX(nbx),myNbY(nby),myNbZ(nbz),
  52  myNbThreads(nbthreads),
  53  myNbTriangles(0),
  54  myUseExistingTriangulation(useExistingTriangulation)
  55 {
  56   Init();
  57 }
  58
  59 Voxel_FastConverter::Voxel_FastConverter(const TopoDS_Shape&    shape,
  60                                          Voxel_BoolDS&          voxels,
  61                                          const Standard_Real    deflection,
  62                                          const Standard_Integer nbx,
  63                                          const Standard_Integer nby,
  64                                          const Standard_Integer nbz,
  65                                          const Standard_Integer nbthreads,
  66                                          const Standard_Boolean useExistingTriangulation)
  67 :myShape(shape),myVoxels(&voxels),
  68  myDeflection(deflection),
  69  myIsBool(1),
  70  myNbX(nbx),myNbY(nby),myNbZ(nbz),
  71  myNbThreads(nbthreads),
  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  myIsBool(0),
  89  myNbX(nbx),myNbY(nby),myNbZ(nbz),
  90  myNbThreads(nbthreads),
  91  myNbTriangles(0),
  92  myUseExistingTriangulation(useExistingTriangulation)
  93 {
  94   Init();
  95 }
  96
  97 void Voxel_FastConverter::Init()
  98 {
  99   if (myShape.IsNull())
 100     return;
 101   if (myNbThreads < 1)
 102     return;
 103
 104   // Check number of splits.
 105   Voxel_DS* voxels = (Voxel_DS*) myVoxels;
 106   if (voxels->GetNbX() != myNbX || voxels->GetNbY() != myNbY || voxels->GetNbZ() != myNbZ)
 107   {
 108     // Compute boundary box of the shape
 109     Bnd_Box box;
 110     BRepBndLib::Add(myShape, box);
 111
 112     // Define the voxel model by means of the boundary box of shape
 113     Standard_Real xmin, ymin, zmin, xmax, ymax, zmax;
 114     box.Get(xmin, ymin, zmin, xmax, ymax, zmax);
 115
 116     // Initialize the voxels.
 117     if (myIsBool == 2)
 118       ((Voxel_ROctBoolDS*) voxels)->Init(xmin, ymin, zmin, xmax - xmin, ymax - ymin, zmax - zmin, myNbX, myNbY, myNbZ);
 119     else if (myIsBool == 1)
 120       ((Voxel_BoolDS*) voxels)->Init(xmin, ymin, zmin, xmax - xmin, ymax - ymin, zmax - zmin, myNbX, myNbY, myNbZ);
 121     else if (myIsBool == 0)
 122       ((Voxel_ColorDS*) voxels)->Init(xmin, ymin, zmin, xmax - xmin, ymax - ymin, zmax - zmin, myNbX, myNbY, myNbZ);
 123   }
 124
 125   // Check presence of triangulation.
 126   TopLoc_Location L;
 127   Standard_Boolean triangulate = Standard_False;
 128   TopExp_Explorer expl(myShape, TopAbs_FACE);
 129   if(myUseExistingTriangulation == Standard_False)
 130   {
 131     for (; expl.More(); expl.Next())
 132     {
 133       const TopoDS_Face & F = TopoDS::Face(expl.Current());
 134       Handle(Poly_Triangulation) T = BRep_Tool::Triangulation(F, L);
 135       if (T.IsNull() || (T->Deflection() > myDeflection))
 136       {
 137         triangulate = Standard_True;
 138         break;
 139       }
 140     }
 141   }
 142
 143   // Re-create the triangulation.
 144   if (triangulate)
 145   {
 146     BRepMesh_IncrementalMesh(myShape, myDeflection);
 147   }
 148
 149   // Compute the number of triangles.
 150   myNbTriangles = 0;
 151   expl.Init(myShape, TopAbs_FACE);
 152   for (; expl.More(); expl.Next())
 153   {
 154     const TopoDS_Face & F = TopoDS::Face(expl.Current());
 155     Handle(Poly_Triangulation) T = BRep_Tool::Triangulation(F, L);
 156     if (T.IsNull() == Standard_False)
 157       myNbTriangles += T->NbTriangles();
 158   }
 159 }
 160
 161 // Destructor
 162 void Voxel_FastConverter::Destroy()
 163 {
 164
 165 }
 166
 167 Standard_Boolean Voxel_FastConverter::Convert(Standard_Integer&      progress,
 168                                               const Standard_Integer ithread)
 169 {
 170   if (ithread == 1)
 171     progress = 0;
 172 #ifdef CONV_DUMP
 173   if (ithread == 1)
 174     printf("Progress = %d   \r", progress);
 175 #endif
 176
 177   if (myNbX <= 0 || myNbY <= 0 || myNbZ <= 0)
 178     return Standard_False;
 179
 180   if(myNbTriangles == 0)
 181     return Standard_False;
 182
 183   // Half of diagonal of a voxel
 184   Voxel_DS* ds = (Voxel_DS*) myVoxels;
 185   Standard_Real dx = ds->GetXLen() / (Standard_Real) ds->GetNbX(),
 186                 dy = ds->GetYLen() / (Standard_Real) ds->GetNbY(),
 187                 dz = ds->GetZLen() / (Standard_Real) ds->GetNbZ();
 188   Standard_Real hdiagonal = sqrt(dx * dx + dy * dy + dz * dz);
 189   hdiagonal /= 2.0;
 190
 191   // Compute the scope of triangles for current thread
 192   Standard_Integer start_thread_triangle = 1, end_thread_triangle = myNbTriangles, ithread_triangle = 0;
 193   if(myNbTriangles < myNbThreads)
 194   {
 195     if(ithread != 1)
 196       return Standard_False;
 197     //in case we're in thread one process all triangles
 198   }
 199   else
 200   {
 201     div_t division = div(myNbTriangles, myNbThreads);
 202     start_thread_triangle = (ithread - 1) * division.quot + 1;
 203     end_thread_triangle   = (ithread - 0) * division.quot;
 204
 205     if(ithread == myNbThreads)
 206       end_thread_triangle += division.rem;
 207   }
 208
 209   // Convert
 210   TopLoc_Location L;
 211   Standard_Integer iprogress = 0;
 212   Standard_Integer n1, n2, n3;
 213   Standard_Integer ixmin, iymin, izmin, ixmax, iymax, izmax;
 214   Standard_Real xmin, ymin, zmin, xmax, ymax, zmax;
 215   TopExp_Explorer expl(myShape, TopAbs_FACE);
 216   for (; expl.More(); expl.Next())
 217   {
 218     const TopoDS_Face & F = TopoDS::Face(expl.Current());
 219     Handle(Poly_Triangulation) T = BRep_Tool::Triangulation(F, L);
 220     if (T.IsNull())
 221       continue;
 222
 223     gp_Trsf trsf;
 224     Standard_Boolean transform = !L.IsIdentity();
 225     if (transform)
 226       trsf = L.Transformation();
 227
 228     const TColgp_Array1OfPnt& nodes = T->Nodes();
 229     const Poly_Array1OfTriangle& triangles = T->Triangles();
 230     Standard_Integer itriangle = triangles.Lower(), nb_triangles = triangles.Upper();
 231     for (; itriangle <= nb_triangles; itriangle++)
 232     {
 233       ithread_triangle++;
 234       if (ithread_triangle < start_thread_triangle )
 235         continue;
 236       if (ithread_triangle > end_thread_triangle)
 237       {
 238         if (ithread == 1)
 239           progress = 100;
 240 #ifdef CONV_DUMP
 241         if (ithread == 1)
 242           printf("Progress = %d  \r", progress);
 243 #endif
 244         return Standard_True;
 245       }
 246
 247       const Poly_Triangle& t = triangles.Value(itriangle);
 248       t.Get(n1, n2, n3);
 249       gp_Pnt p1 = nodes.Value(n1);
 250       gp_Pnt p2 = nodes.Value(n2);
 251       gp_Pnt p3 = nodes.Value(n3);
 252       if (transform)
 253       {
 254         p1.Transform(trsf);
 255         p2.Transform(trsf);
 256         p3.Transform(trsf);
 257       }
 258
 259       // Get boundary box of the triangle
 260       GetBndBox(p1, p2, p3, xmin, ymin, zmin, xmax, ymax, zmax);
 261
 262       // Find the range of voxels inside the boudary box of the triangle.
 263       if (!ds->GetVoxel(xmin, ymin, zmin, ixmin, iymin, izmin))
 264         continue;
 265       if (!ds->GetVoxel(xmax, ymax, zmax, ixmax, iymax, izmax))
 266         continue;
 267
 268       // Refuse voxels for whom distance from their center to plane of triangle is greater than half of diagonal.
 269       // Make a line from center of each voxel to the center of triangle and
 270       // compute intersection of the line with sides of triangle.
 271       // Refuse the voxel in case of intersection.
 272       gce_MakePln mkPlane(p1, p2, p3);
 273       if (!mkPlane.IsDone())
 274         continue;
 275       gp_Pln plane = mkPlane.Value();
 276       ComputeVoxelsNearTriangle(plane, p1, p2, p3, hdiagonal, ixmin, iymin, izmin, ixmax, iymax, izmax);
 277
 278       // Progress
 279       if (ithread == 1)
 280       {
 281         iprogress++;
 282         progress = (Standard_Integer) ( (Standard_Real) iprogress / (Standard_Real) myNbTriangles * 100.0 );
 283       }
 284 #ifdef CONV_DUMP
 285       if (ithread == 1 && prev_progress != progress)
 286       {
 287         printf("Progress = %d  \r", progress);
 288         prev_progress = progress;
 289       }
 290 #endif
 291
 292     } // iteration of triangles
 293   } // iteration of faces
 294
 295   if (ithread == 1)
 296     progress = 100;
 297 #ifdef CONV_DUMP
 298   if (ithread == 1)
 299     printf("Progress = %d  \r", progress);
 300 #endif
 301   return Standard_True;
 302 }
 303
 304 Standard_Boolean Voxel_FastConverter::ConvertUsingSAT(Standard_Integer&      progress,
 305                                                       const Standard_Integer ithread)
 306 {
 307   if (ithread == 1)
 308     progress = 0;
 309 #ifdef CONV_DUMP
 310   if (ithread == 1)
 311     printf("Progress = %d   \r", progress);
 312 #endif
 313
 314   if (myNbX <= 0 || myNbY <= 0 || myNbZ <= 0)
 315     return Standard_False;
 316
 317   if(myNbTriangles == 0)
 318     return Standard_False;
 319
 320   // Half of size of a voxel (also for Voxel_ROctBoolDS)
 321   Voxel_DS* ds = (Voxel_DS*) myVoxels;
 322   Standard_Real dx = ds->GetXLen() / (Standard_Real) ds->GetNbX(),
 323                 dy = ds->GetYLen() / (Standard_Real) ds->GetNbY(),
 324                 dz = ds->GetZLen() / (Standard_Real) ds->GetNbZ();
 325   gp_Pnt extents(dx/2.0, dy/2.0, dz/2.0);
 326   gp_Pnt extents2(dx/4.0, dy/4.0, dz/4.0);
 327   gp_Pnt extents4(dx/8.0, dy/8.0, dz/8.0);
 328
 329   // Compute the scope of triangles for current thread
 330   Standard_Integer start_thread_triangle = 1, end_thread_triangle = myNbTriangles, ithread_triangle = 0;
 331   if(myNbTriangles < myNbThreads)
 332   {
 333     if(ithread != 1)
 334       return Standard_False;
 335     //in case we're in thread one process all triangles
 336   }
 337   else
 338   {
 339     div_t division = div(myNbTriangles, myNbThreads);
 340     start_thread_triangle = (ithread - 1) * division.quot + 1;
 341     end_thread_triangle   = (ithread - 0) * division.quot;
 342
 343     if(ithread == myNbThreads)
 344       end_thread_triangle += division.rem;
 345   }
 346
 347   // Convert
 348   TopLoc_Location L;
 349   Standard_Integer iprogress = 0;
 350   Standard_Integer n1, n2, n3;
 351   Standard_Integer ixmin, iymin, izmin, ixmax, iymax, izmax;
 352   Standard_Real xmin, ymin, zmin, xmax, ymax, zmax;
 353   TopExp_Explorer expl(myShape, TopAbs_FACE);
 354   for (; expl.More(); expl.Next())
 355   {
 356     const TopoDS_Face & F = TopoDS::Face(expl.Current());
 357     Handle(Poly_Triangulation) T = BRep_Tool::Triangulation(F, L);
 358     if (T.IsNull())
 359       continue;
 360
 361     gp_Trsf trsf;
 362     Standard_Boolean transform = !L.IsIdentity();
 363     if (transform)
 364       trsf = L.Transformation();
 365
 366     const TColgp_Array1OfPnt& nodes = T->Nodes();
 367     const Poly_Array1OfTriangle& triangles = T->Triangles();
 368     Standard_Integer itriangle = triangles.Lower(), nb_triangles = triangles.Upper();
 369     for (; itriangle <= nb_triangles; itriangle++)
 370     {
 371       ithread_triangle++;
 372       if (ithread_triangle < start_thread_triangle )
 373         continue;
 374       if (ithread_triangle > end_thread_triangle)
 375       {
 376         if (ithread == 1)
 377           progress = 100;
 378 #ifdef CONV_DUMP
 379         if (ithread == 1)
 380           printf("Progress = %d  \r", progress);
 381 #endif
 382         return Standard_True;
 383       }
 384
 385       const Poly_Triangle& t = triangles.Value(itriangle);
 386
 387       t.Get(n1, n2, n3);
 388       gp_Pnt p1 = nodes.Value(n1);
 389       gp_Pnt p2 = nodes.Value(n2);
 390       gp_Pnt p3 = nodes.Value(n3);
 391       if (transform)
 392       {
 393         p1.Transform(trsf);
 394         p2.Transform(trsf);
 395         p3.Transform(trsf);
 396       }
 397
 398       // Get boundary box of the triangle
 399       GetBndBox(p1, p2, p3, xmin, ymin, zmin, xmax, ymax, zmax);
 400
 401       // Find the range of voxels inside the boudary box of the triangle.
 402       if (!ds->GetVoxel(xmin, ymin, zmin, ixmin, iymin, izmin))
 403         continue;
 404       if (!ds->GetVoxel(xmax, ymax, zmax, ixmax, iymax, izmax))
 405         continue;
 406
 407       // Perform triangle-box intersection to find the voxels resulting from the processed triangle.;
 408       // Using SAT theorem to quickly find the intersection.
 409       ComputeVoxelsNearTriangle(p1, p2, p3,
 410                                 extents, extents2, extents4,
 411                                 ixmin, iymin, izmin, ixmax, iymax, izmax);
 412
 413       // Progress
 414       if (ithread == 1)
 415       {
 416         iprogress++;
 417         progress = (Standard_Integer) ( (Standard_Real) iprogress / (Standard_Real) myNbTriangles * 100.0 );
 418       }
 419 #ifdef CONV_DUMP
 420       if (ithread == 1 && prev_progress != progress)
 421       {
 422         printf("Progress = %d  \r", progress);
 423         prev_progress = progress;
 424       }
 425 #endif
 426
 427     } // iteration of triangles
 428   } // iteration of faces
 429
 430   if (ithread == 1)
 431     progress = 100;
 432 #ifdef CONV_DUMP
 433   if (ithread == 1)
 434     printf("Progress = %d  \r", progress);
 435 #endif
 436   return Standard_True;
 437 }
 438
 439 Standard_Boolean Voxel_FastConverter::FillInVolume(const Standard_Byte inner,
 440                                                    const Standard_Integer /*ithread*/)
 441 {
 442   Voxel_DS* ds = (Voxel_DS*) myVoxels;
 443   Standard_Integer ix, iy, iz, nbx = ds->GetNbX(), nby = ds->GetNbY(), nbz = ds->GetNbZ();
 444   Standard_Boolean prev_surface, surface, volume;
 445
 446   if (inner)
 447   {
 448     // Fill-in internal voxels by the value "inner"
 449     for (ix = 0; ix < nbx; ix++)
 450     {
 451       for (iy = 0; iy < nby; iy++)
 452       {
 453         // Check existence of volume.
 454         volume = Standard_False;
 455         surface = Standard_False;
 456         prev_surface = Standard_False;
 457         for (iz = 0; iz < nbz; iz++)
 458         {
 459           surface = (myIsBool == 1) ?
 460             ((Voxel_BoolDS*)myVoxels)->Get(ix, iy, iz) == Standard_True :
 461               ((Voxel_ColorDS*)myVoxels)->Get(ix, iy, iz) > 0;
 462           if (prev_surface && !surface)
 463           {
 464             volume = !volume;
 465           }
 466           prev_surface = surface;
 467         }
 468         if (volume)
 469           continue;
 470
 471         // Fill-in the volume.
 472         volume = Standard_False;
 473         surface = Standard_False;
 474         prev_surface = Standard_False;
 475         for (iz = 0; iz < nbz; iz++)
 476         {
 477           surface = (myIsBool == 1) ?
 478             ((Voxel_BoolDS*)myVoxels)->Get(ix, iy, iz) == Standard_True :
 479               ((Voxel_ColorDS*)myVoxels)->Get(ix, iy, iz) > 0;
 480           if (prev_surface && !surface)
 481           {
 482             volume = !volume;
 483           }
 484           if (volume && !surface)
 485           {
 486             (myIsBool == 1) ? ((Voxel_BoolDS*)myVoxels)->Set(ix, iy, iz, inner) :
 487               ((Voxel_ColorDS*)myVoxels)->Set(ix, iy, iz, inner);
 488           }
 489           prev_surface = surface;
 490         }
 491       }
 492     }
 493   }
 494   else
 495   {
 496     // Set value of interbal voxels to 0 ("inner" = 0)
 497     Standard_Boolean next_surface;
 498     for (ix = 0; ix < nbx; ix++)
 499     {
 500       for (iy = 0; iy < nby; iy++)
 501       {
 502         volume = Standard_False;
 503         surface = Standard_False;
 504         prev_surface = Standard_False;
 505         next_surface = Standard_False;
 506         for (iz = 0; iz < nbz; iz++)
 507         {
 508           surface = (myIsBool == 1) ?
 509             ((Voxel_BoolDS*)myVoxels)->Get(ix, iy, iz) == Standard_True :
 510               ((Voxel_ColorDS*)myVoxels)->Get(ix, iy, iz) > 0;
 511           if (prev_surface != surface)
 512           {
 513             volume = !volume;
 514           }
 515           if (volume && iz + 1 < nbz)
 516           {
 517             next_surface = (myIsBool == 1) ?
 518               ((Voxel_BoolDS*)myVoxels)->Get(ix, iy, iz + 1) == Standard_True :
 519               ((Voxel_ColorDS*)myVoxels)->Get(ix, iy, iz + 1) > 0;
 520           }
 521           if (volume && prev_surface == surface && next_surface)
 522           {
 523             (myIsBool == 1) ? ((Voxel_BoolDS*)myVoxels)->Set(ix, iy, iz, inner) :
 524               ((Voxel_ColorDS*)myVoxels)->Set(ix, iy, iz, inner);
 525           }
 526           prev_surface = surface;
 527         }
 528       }
 529     }
 530   }
 531
 532   return Standard_True;
 533 }
 534
 535 Standard_Boolean Voxel_FastConverter::FillInVolume(const Standard_Byte inner, const TopoDS_Shape & shape, const Standard_Integer /*ithread*/)
 536 {
 537   Voxel_DS* ds = (Voxel_DS*) myVoxels;
 538   Standard_Integer ix, iy, iz, nbx = ds->GetNbX(), nby = ds->GetNbY(), nbz = ds->GetNbZ();
 539   Standard_Boolean prev_surface, surface, volume, isOnVerticalSurface;
 540
 541   BRepClass3d_SolidClassifier solidClassifier(shape);
 542   Standard_Real xc, yc, zc;
 543
 544   if (inner)
 545   {
 546     // Fill-in internal voxels by the value "inner"
 547     for (ix = 0; ix < nbx; ix++)
 548     {
 549       for (iy = 0; iy < nby; iy++)
 550       {
 551         // Check existence of volume.
 552         volume = Standard_False;
 553         surface = Standard_False;
 554         prev_surface = Standard_False;
 555         isOnVerticalSurface = Standard_False;
 556         for (iz = 0; iz < nbz; iz++)
 557         {
 558           surface = (myIsBool == 1) ?
 559             ((Voxel_BoolDS*)myVoxels)->Get(ix, iy, iz) == Standard_True :
 560               ((Voxel_ColorDS*)myVoxels)->Get(ix, iy, iz) > 0;
 561           if (prev_surface && !surface)
 562           {
 563             if(isOnVerticalSurface)
 564             {
 565               isOnVerticalSurface = Standard_False;
 566               ((Voxel_BoolDS*)myVoxels)->GetCenter(ix, iy, iz, xc, yc, zc);
 567               gp_Pnt P(xc, yc, zc);
 568               solidClassifier.Perform(P, Precision::Confusion());
 569
 570               if(solidClassifier.State() == TopAbs_IN)
 571                 volume = Standard_True;
 572               else
 573                 volume = Standard_False;
 574             }
 575             else
 576               volume = !volume;
 577           }
 578           if(prev_surface && surface)
 579             isOnVerticalSurface = Standard_True;
 580           else
 581             isOnVerticalSurface = Standard_False;
 582           prev_surface = surface;
 583         }
 584         if (volume)
 585           continue;
 586
 587         // Fill-in the volume.
 588         volume = Standard_False;
 589         surface = Standard_False;
 590         prev_surface = Standard_False;
 591         isOnVerticalSurface = Standard_False;
 592         for (iz = 0; iz < nbz; iz++)
 593         {
 594           surface = (myIsBool == 1) ?
 595             ((Voxel_BoolDS*)myVoxels)->Get(ix, iy, iz) == Standard_True :
 596               ((Voxel_ColorDS*)myVoxels)->Get(ix, iy, iz) > 0;
 597           if (prev_surface && !surface)
 598           {
 599             if(isOnVerticalSurface)
 600             {
 601               isOnVerticalSurface = Standard_False;
 602               ((Voxel_BoolDS*)myVoxels)->GetCenter(ix, iy, iz, xc, yc, zc);
 603               gp_Pnt P(xc, yc, zc);
 604               solidClassifier.Perform(P, Precision::Confusion());
 605
 606               if(solidClassifier.State() == TopAbs_IN)
 607                 volume = Standard_True;
 608               else
 609                 volume = Standard_False;
 610             }
 611             else
 612               volume = !volume;
 613           }
 614           if (volume && !surface)
 615           {
 616             (myIsBool == 1) ? ((Voxel_BoolDS*)myVoxels)->Set(ix, iy, iz, inner) :
 617               ((Voxel_ColorDS*)myVoxels)->Set(ix, iy, iz, inner);
 618           }
 619           if(prev_surface && surface)
 620             isOnVerticalSurface = Standard_True;
 621           else
 622             isOnVerticalSurface = Standard_False;
 623           prev_surface = surface;
 624         }
 625       }
 626     }
 627   }
 628
 629   return Standard_True;
 630 }
 631
 632 void Voxel_FastConverter::GetBndBox(const gp_Pnt&  p1,
 633                                     const gp_Pnt&  p2,
 634                                     const gp_Pnt&  p3,
 635                                     Standard_Real& xmin,
 636                                     Standard_Real& ymin,
 637                                     Standard_Real& zmin,
 638                                     Standard_Real& xmax,
 639                                     Standard_Real& ymax,
 640                                     Standard_Real& zmax) const
 641 {
 642   // P1:
 643   xmin = p1.X();
 644   ymin = p1.Y();
 645   zmin = p1.Z();
 646   xmax = p1.X();
 647   ymax = p1.Y();
 648   zmax = p1.Z();
 649   // P2:
 650   if (xmin > p2.X())
 651     xmin = p2.X();
 652   if (ymin > p2.Y())
 653     ymin = p2.Y();
 654   if (zmin > p2.Z())
 655     zmin = p2.Z();
 656   if (xmax < p2.X())
 657     xmax = p2.X();
 658   if (ymax < p2.Y())
 659     ymax = p2.Y();
 660   if (zmax < p2.Z())
 661     zmax = p2.Z();
 662   // P3:
 663   if (xmin > p3.X())
 664     xmin = p3.X();
 665   if (ymin > p3.Y())
 666     ymin = p3.Y();
 667   if (zmin > p3.Z())
 668     zmin = p3.Z();
 669   if (xmax < p3.X())
 670     xmax = p3.X();
 671   if (ymax < p3.Y())
 672     ymax = p3.Y();
 673   if (zmax < p3.Z())
 674     zmax = p3.Z();
 675 }
 676
 677 // This method is copied from Voxel_ShapeIntersector.cxx
 678 static Standard_Boolean mayIntersect(const gp_Pnt2d& p11, const gp_Pnt2d& p12,
 679                                      const gp_Pnt2d& p21, const gp_Pnt2d& p22)
 680 {
 681     if (p11.X() > p21.X() && p11.X() > p22.X() && p12.X() > p21.X() && p12.X() > p22.X())
 682         return Standard_False;
 683     if (p11.X() < p21.X() && p11.X() < p22.X() && p12.X() < p21.X() && p12.X() < p22.X())
 684         return Standard_False;
 685     if (p11.Y() > p21.Y() && p11.Y() > p22.Y() && p12.Y() > p21.Y() && p12.Y() > p22.Y())
 686         return Standard_False;
 687     if (p11.Y() < p21.Y() && p11.Y() < p22.Y() && p12.Y() < p21.Y() && p12.Y() < p22.Y())
 688         return Standard_False;
 689     return Standard_True;
 690 }
 691
 692 void Voxel_FastConverter::ComputeVoxelsNearTriangle(const gp_Pln&          plane,
 693                                                     const gp_Pnt&          p1,
 694                                                     const gp_Pnt&          p2,
 695                                                     const gp_Pnt&          p3,
 696                                                     const Standard_Real    hdiagonal,
 697                                                     const Standard_Integer ixmin,
 698                                                     const Standard_Integer iymin,
 699                                                     const Standard_Integer izmin,
 700                                                     const Standard_Integer ixmax,
 701                                                     const Standard_Integer iymax,
 702                                                     const Standard_Integer izmax) const
 703 {
 704   gp_Pnt pc;
 705   Standard_Real xc, yc, zc, uc, vc, u1, v1, u2, v2, u3, v3;
 706   Standard_Integer ix, iy, iz;
 707   IntAna2d_AnaIntersection intersector2d;
 708
 709   // Project points of triangle onto the plane
 710   ElSLib::Parameters(plane, p1, u1, v1);
 711   ElSLib::Parameters(plane, p2, u2, v2);
 712   ElSLib::Parameters(plane, p3, u3, v3);
 713
 714   // Make lines of triangle
 715   gp_Pnt2d p2d1(u1, v1), p2d2(u2, v2), p2d3(u3, v3), p2dt((u1+u2+u3)/3.0,(v1+v2+v3)/3.0), p2dc;
 716   gp_Vec2d v2d12(p2d1, p2d2), v2d23(p2d2, p2d3), v2d31(p2d3, p2d1);
 717   gp_Lin2d L1(p2d1, v2d12), L2(p2d2, v2d23), L3(p2d3, v2d31), Lv;
 718   Standard_Real d1 = p2d1.Distance(p2d2) - Precision::Confusion(),
 719                 d2 = p2d2.Distance(p2d3) - Precision::Confusion(),
 720                 d3 = p2d3.Distance(p2d1) - Precision::Confusion(), dv;
 721
 722   Voxel_DS* ds = (Voxel_DS*) myVoxels;
 723   for (ix = ixmin; ix <= ixmax; ix++)
 724   {
 725     for (iy = iymin; iy <= iymax; iy++)
 726     {
 727       for (iz = izmin; iz <= izmax; iz++)
 728       {
 729         ds->GetCenter(ix, iy, iz, xc, yc, zc);
 730         pc.SetCoord(xc, yc, zc);
 731         if (plane.Distance(pc) < hdiagonal)
 732         {
 733           ElSLib::Parameters(plane, pc, uc, vc);
 734           p2dc.SetCoord(uc, vc);
 735
 736           gp_Vec2d v2dct(p2dc, p2dt);
 737           dv = v2dct.Magnitude() - Precision::Confusion();
 738           Lv.SetLocation(p2dc);
 739           Lv.SetDirection(v2dct);
 740
 741           // Side 1:
 742           if (mayIntersect(p2d1, p2d2, p2dc, p2dt))
 743           {
 744             intersector2d.Perform(Lv, L1);
 745             if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
 746             {
 747               const IntAna2d_IntPoint& i2d = intersector2d.Point(1);
 748               Standard_Real param1 = i2d.ParamOnFirst();
 749               Standard_Real param2 = i2d.ParamOnSecond();
 750               if (param1 > Precision::Confusion() && param1 < dv &&
 751                   param2 > Precision::Confusion() && param2 < d1)
 752               {
 753                 continue;
 754               }
 755             }
 756           }
 757
 758           // Side 2:
 759           if (mayIntersect(p2d2, p2d3, p2dc, p2dt))
 760           {
 761             intersector2d.Perform(Lv, L2);
 762             if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
 763             {
 764               const IntAna2d_IntPoint& i2d = intersector2d.Point(1);
 765               Standard_Real param1 = i2d.ParamOnFirst();
 766               Standard_Real param2 = i2d.ParamOnSecond();
 767               if (param1 > Precision::Confusion() && param1 < dv &&
 768                   param2 > Precision::Confusion() && param2 < d2)
 769               {
 770                 continue;
 771               }
 772             }
 773           }
 774
 775           // Side 3:
 776           if (mayIntersect(p2d3, p2d1, p2dc, p2dt))
 777           {
 778             intersector2d.Perform(Lv, L3);
 779             if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
 780             {
 781               const IntAna2d_IntPoint& i2d = intersector2d.Point(1);
 782               Standard_Real param1 = i2d.ParamOnFirst();
 783               Standard_Real param2 = i2d.ParamOnSecond();
 784               if (param1 > Precision::Confusion() && param1 < dv &&
 785                   param2 > Precision::Confusion() && param2 < d3)
 786               {
 787                 continue;
 788               }
 789             }
 790           }
 791
 792           // Set positive value to this voxel:
 793           switch (myIsBool)
 794           {
 795             case 0:
 796               ((Voxel_ColorDS*) myVoxels)->Set(ix, iy, iz, 15);
 797               break;
 798             case 1:
 799               ((Voxel_BoolDS*) myVoxels)->Set(ix, iy, iz, Standard_True);
 800               break;
 801             case 2:
 802             {
 803               //((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, Standard_True);
 804
 805               // Check intersection between the triangle & sub-voxels of the voxel.
 806               Standard_Real hdiagonal2 = hdiagonal / 2.0, hdiagonal4 = hdiagonal / 4.0;
 807               for (Standard_Integer i = 0; i < 8; i++)
 808               {
 809                 ((Voxel_ROctBoolDS*) myVoxels)->GetCenter(ix, iy, iz, i, xc, yc, zc);
 810                 pc.SetCoord(xc, yc, zc);
 811                 if (plane.Distance(pc) < hdiagonal2)
 812                 {
 813                   ElSLib::Parameters(plane, pc, uc, vc);
 814                   p2dc.SetCoord(uc, vc);
 815
 816                   gp_Vec2d aVec2dct1(p2dc, p2dt);
 817                   dv = aVec2dct1.Magnitude() - Precision::Confusion();
 818                   Lv.SetLocation(p2dc);
 819                   Lv.SetDirection(aVec2dct1);
 820
 821                   // Side 1:
 822                   if (mayIntersect(p2d1, p2d2, p2dc, p2dt))
 823                   {
 824                     intersector2d.Perform(Lv, L1);
 825                     if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
 826                     {
 827                       const IntAna2d_IntPoint& i2d = intersector2d.Point(1);
 828                       Standard_Real param1 = i2d.ParamOnFirst();
 829                       Standard_Real param2 = i2d.ParamOnSecond();
 830                       if (param1 > Precision::Confusion() && param1 < dv &&
 831                           param2 > Precision::Confusion() && param2 < d1)
 832                       {
 833                         continue;
 834                       }
 835                     }
 836                   }
 837
 838                   // Side 2:
 839                   if (mayIntersect(p2d2, p2d3, p2dc, p2dt))
 840                   {
 841                     intersector2d.Perform(Lv, L2);
 842                     if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
 843                     {
 844                       const IntAna2d_IntPoint& i2d = intersector2d.Point(1);
 845                       Standard_Real param1 = i2d.ParamOnFirst();
 846                       Standard_Real param2 = i2d.ParamOnSecond();
 847                       if (param1 > Precision::Confusion() && param1 < dv &&
 848                           param2 > Precision::Confusion() && param2 < d2)
 849                       {
 850                         continue;
 851                       }
 852                     }
 853                   }
 854
 855                   // Side 3:
 856                   if (mayIntersect(p2d3, p2d1, p2dc, p2dt))
 857                   {
 858                     intersector2d.Perform(Lv, L3);
 859                     if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
 860                     {
 861                       const IntAna2d_IntPoint& i2d = intersector2d.Point(1);
 862                       Standard_Real param1 = i2d.ParamOnFirst();
 863                       Standard_Real param2 = i2d.ParamOnSecond();
 864                       if (param1 > Precision::Confusion() && param1 < dv &&
 865                           param2 > Precision::Confusion() && param2 < d3)
 866                       {
 867                         continue;
 868                       }
 869                     }
 870                   }
 871
 872                   //((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, i, Standard_True);
 873
 874                   // Check intersection between the triangle & sub-voxels of the sub-voxel.
 875                   for (Standard_Integer j = 0; j < 8; j++)
 876                   {
 877                     ((Voxel_ROctBoolDS*) myVoxels)->GetCenter(ix, iy, iz, i, j, xc, yc, zc);
 878                     pc.SetCoord(xc, yc, zc);
 879                     if (plane.Distance(pc) < hdiagonal4)
 880                     {
 881                       ElSLib::Parameters(plane, pc, uc, vc);
 882                       p2dc.SetCoord(uc, vc);
 883
 884                       gp_Vec2d aVec2dct2(p2dc, p2dt);
 885                       dv = aVec2dct2.Magnitude() - Precision::Confusion();
 886                       Lv.SetLocation(p2dc);
 887                       Lv.SetDirection(aVec2dct2);
 888
 889                       // Side 1:
 890                       if (mayIntersect(p2d1, p2d2, p2dc, p2dt))
 891                       {
 892                         intersector2d.Perform(Lv, L1);
 893                         if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
 894                         {
 895                           const IntAna2d_IntPoint& i2d = intersector2d.Point(1);
 896                           Standard_Real param1 = i2d.ParamOnFirst();
 897                           Standard_Real param2 = i2d.ParamOnSecond();
 898                           if (param1 > Precision::Confusion() && param1 < dv &&
 899                               param2 > Precision::Confusion() && param2 < d1)
 900                           {
 901                             continue;
 902                           }
 903                         }
 904                       }
 905
 906                       // Side 2:
 907                       if (mayIntersect(p2d2, p2d3, p2dc, p2dt))
 908                       {
 909                         intersector2d.Perform(Lv, L2);
 910                         if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
 911                         {
 912                           const IntAna2d_IntPoint& i2d = intersector2d.Point(1);
 913                           Standard_Real param1 = i2d.ParamOnFirst();
 914                           Standard_Real param2 = i2d.ParamOnSecond();
 915                           if (param1 > Precision::Confusion() && param1 < dv &&
 916                               param2 > Precision::Confusion() && param2 < d2)
 917                           {
 918                             continue;
 919                           }
 920                         }
 921                       }
 922
 923                       // Side 3:
 924                       if (mayIntersect(p2d3, p2d1, p2dc, p2dt))
 925                       {
 926                         intersector2d.Perform(Lv, L3);
 927                         if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
 928                         {
 929                           const IntAna2d_IntPoint& i2d = intersector2d.Point(1);
 930                           Standard_Real param1 = i2d.ParamOnFirst();
 931                           Standard_Real param2 = i2d.ParamOnSecond();
 932                           if (param1 > Precision::Confusion() && param1 < dv &&
 933                               param2 > Precision::Confusion() && param2 < d3)
 934                           {
 935                             continue;
 936                           }
 937                         }
 938                       }
 939
 940                       ((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, i, j, Standard_True);
 941                     }
 942                   } // End of "Check level 2".
 943
 944                 }
 945               } // End of "Check level 1".
 946
 947               break;
 948             }
 949           }
 950         }
 951       }
 952     }
 953   }
 954 }
 955
 956 //! This macro quickly finds the min & max values among 3 variables
 957 #define FINDMINMAX(x0, x1, x2, min, max)        \
 958   min = max = x0;                               \
 959   if(x1<min) min=x1;                            \
 960   if(x1>max) max=x1;                            \
 961   if(x2<min) min=x2;                            \
 962   if(x2>max) max=x2;
 963
 964 static bool planeBoxOverlap(const gp_Vec & normal, const double d, const gp_Pnt & maxbox)
 965 {
 966   gp_Vec vmin, vmax;
 967   if(normal.X() > 0.0)    { vmin.SetX(-maxbox.X()); vmax.SetX(maxbox.X()); }
 968   else                    { vmin.SetX(maxbox.X());  vmax.SetX(-maxbox.X()); }
 969
 970   if(normal.Y() > 0.0)    { vmin.SetY(-maxbox.Y()); vmax.SetY(maxbox.Y()); }
 971   else                    { vmin.SetY(maxbox.Y());  vmax.SetY(-maxbox.Y()); }
 972
 973   if(normal.Z() > 0.0)    { vmin.SetZ(-maxbox.Z()); vmax.SetZ(maxbox.Z()); }
 974   else                    { vmin.SetZ(maxbox.Z());  vmax.SetZ(-maxbox.Z()); }
 975
 976   if((normal.Dot(vmin)) + d > 0.0) return false;
 977   if((normal.Dot(vmax)) + d>= 0.0) return true;
 978
 979   return false;
 980 }
 981
 982 #define AXISTEST_X01(a, b, fa, fb)                            \
 983     min = a*v0.Y() - b*v0.Z();                                    \
 984     max = a*v2.Y() - b*v2.Z();                                    \
 985     if(min>max) {const double tmp=max; max=min; min=tmp;    }    \
 986     rad = fa * extents.Y() + fb * extents.Z();                    \
 987     if(min>rad || max<-rad) return false;
 988
 989 #define AXISTEST_X2(a, b, fa, fb)                            \
 990     min = a*v0.Y() - b*v0.Z();                                    \
 991     max = a*v1.Y() - b*v1.Z();                                    \
 992     if(min>max) {const double tmp=max; max=min; min=tmp;    }    \
 993     rad = fa * extents.Y() + fb * extents.Z();                    \
 994     if(min>rad || max<-rad) return false;
 995
 996 #define AXISTEST_Y02(a, b, fa, fb)                            \
 997     min = b*v0.Z() - a*v0.X();                                    \
 998     max = b*v2.Z() - a*v2.X();                                    \
 999     if(min>max) {const double tmp=max; max=min; min=tmp;    }    \
1000     rad = fa * extents.X() + fb * extents.Z();                    \
1001     if(min>rad || max<-rad) return false;
1002
1003 #define AXISTEST_Y1(a, b, fa, fb)                            \
1004     min = b*v0.Z() - a*v0.X();                                    \
1005     max = b*v1.Z() - a*v1.X();                                    \
1006     if(min>max) {const double tmp=max; max=min; min=tmp;    }    \
1007     rad = fa * extents.X() + fb * extents.Z();                    \
1008     if(min>rad || max<-rad) return false;
1009
1010 #define AXISTEST_Z12(a, b, fa, fb)                            \
1011     min = a*v1.X() - b*v1.Y();                                    \
1012     max = a*v2.X() - b*v2.Y();                                    \
1013     if(min>max) {const double tmp=max; max=min; min=tmp;    }    \
1014     rad = fa * extents.X() + fb * extents.Y();                    \
1015     if(min>rad || max<-rad) return false;
1016
1017 #define AXISTEST_Z0(a, b, fa, fb)                            \
1018     min = a*v0.X() - b*v0.Y();                                    \
1019     max = a*v1.X() - b*v1.Y();                                    \
1020     if(min>max) {const double tmp=max; max=min; min=tmp;    }    \
1021     rad = fa * extents.X() + fb * extents.Y();                    \
1022     if(min>rad || max<-rad) return false;
1023
1024 // compute triangle edges
1025 // - edges lazy evaluated to take advantage of early exits
1026 // - fabs precomputed (half less work, possible since extents are always >0)
1027 // - customized macros to take advantage of the null component
1028 // - axis vector discarded, possibly saves useless movs
1029 #define IMPLEMENT_CLASS3_TESTS                        \
1030   double rad;                                         \
1031   double min, max;                                    \
1032                                                       \
1033   const double fey0 = fabs(e0.Y());                   \
1034   const double fez0 = fabs(e0.Z());                   \
1035   AXISTEST_X01(e0.Z(), e0.Y(), fez0, fey0);           \
1036   const double fex0 = fabs(e0.X());                   \
1037   AXISTEST_Y02(e0.Z(), e0.X(), fez0, fex0);           \
1038   AXISTEST_Z12(e0.Y(), e0.X(), fey0, fex0);           \
1039                                                       \
1040   const double fey1 = fabs(e1.Y());                   \
1041   const double fez1 = fabs(e1.Z());                   \
1042   AXISTEST_X01(e1.Z(), e1.Y(), fez1, fey1);           \
1043   const double fex1 = fabs(e1.X());                   \
1044   AXISTEST_Y02(e1.Z(), e1.X(), fez1, fex1);           \
1045   AXISTEST_Z0(e1.Y(), e1.X(), fey1, fex1);            \
1046                                                       \
1047   const gp_Vec e2 = v2 - v0;                          \
1048   const double fey2 = fabs(e2.Y());                   \
1049   const double fez2 = fabs(e2.Z());                   \
1050   AXISTEST_X2(e2.Z(), e2.Y(), fez2, fey2);            \
1051   const double fex2 = fabs(e2.X());                   \
1052   AXISTEST_Y1(e2.Z(), e2.X(), fez2, fex2);            \
1053   AXISTEST_Z12(e2.Y(), e2.X(), fey2, fex2);
1054
1055 static bool TriBoxOverlap(const gp_Pnt & p1, const gp_Pnt & p2, const gp_Pnt & p3,
1056                           const gp_Pnt & center, const gp_Pnt & extents)
1057 {
1058   // use separating axis theorem to test overlap between triangle and box
1059   // need to test for overlap in these directions:
1060   // 1) the {x,y,z}-directions (actually, since we use the AABB of the triangle
1061   //    we do not even need to test these)
1062   // 2) normal of the triangle
1063   // 3) crossproduct(edge from tri, {x,y,z}-directin)
1064   //    this gives 3x3=9 more tests
1065
1066   // move everything so that the boxcenter is in (0,0,0)
1067   gp_Vec v0(center, p1);
1068   gp_Vec v1(center, p2);
1069   gp_Vec v2(center, p3);
1070
1071   // First, test overlap in the {x,y,z}-directions
1072   double aMin,aMax;
1073   // Find min, max of the triangle in x-direction, and test for overlap in X
1074   FINDMINMAX(v0.X(), v1.X(), v2.X(), aMin, aMax);
1075   if(aMin>extents.X() || aMax<-extents.X()) return false;
1076
1077   // Same for Y
1078   FINDMINMAX(v0.Y(), v1.Y(), v2.Y(), aMin, aMax);
1079   if(aMin>extents.Y() || aMax<-extents.Y()) return false;
1080
1081   // Same for Z
1082   FINDMINMAX(v0.Z(), v1.Z(), v2.Z(), aMin, aMax);
1083   if(aMin>extents.Z() || aMax<-extents.Z()) return false;
1084
1085   // 2) Test if the box intersects the plane of the triangle
1086   // compute plane equation of triangle: normal*x+d=0
1087   // ### could be precomputed since we use the same leaf triangle several times
1088   const gp_Vec e0 = v1 - v0;
1089   const gp_Vec e1 = v2 - v1;
1090   const gp_Vec normal = e0.Crossed(e1);
1091   const double d = -normal.Dot(v0);
1092   if(!planeBoxOverlap(normal, d, extents)) return false;
1093
1094   // 3) "Class III" tests
1095   //if(mFullPrimBoxTest)
1096   {
1097     IMPLEMENT_CLASS3_TESTS
1098   }
1099
1100   return true;
1101 }
1102
1103 void Voxel_FastConverter::ComputeVoxelsNearTriangle(const gp_Pnt& p1,
1104                                                     const gp_Pnt& p2,
1105                                                     const gp_Pnt& p3,
1106                                                     const gp_Pnt& extents,
1107                                                     const gp_Pnt& extents2,
1108                                                     const gp_Pnt& extents4,
1109                                                     const Standard_Integer ixmin,
1110                                                     const Standard_Integer iymin,
1111                                                     const Standard_Integer izmin,
1112                                                     const Standard_Integer ixmax,
1113                                                     const Standard_Integer iymax,
1114                                                     const Standard_Integer izmax) const
1115 {
1116   gp_Pnt pc;
1117   Standard_Real xc, yc, zc;
1118   Standard_Integer ix, iy, iz;
1119
1120   Voxel_DS* ds = (Voxel_DS*) myVoxels;
1121   for (ix = ixmin; ix <= ixmax; ix++)
1122   {
1123     for (iy = iymin; iy <= iymax; iy++)
1124     {
1125       for (iz = izmin; iz <= izmax; iz++)
1126       {
1127         ds->GetCenter(ix, iy, iz, xc, yc, zc);
1128         pc.SetCoord(xc, yc, zc);
1129
1130         if(TriBoxOverlap(p1, p2, p3, pc, extents))
1131         {
1132           // Set positive value to this voxel:
1133           switch (myIsBool)
1134           {
1135             case 0:
1136               ((Voxel_ColorDS*) myVoxels)->Set(ix, iy, iz, 15);
1137               break;
1138             case 1:
1139               ((Voxel_BoolDS*) myVoxels)->Set(ix, iy, iz, Standard_True);
1140               break;
1141             case 2:
1142             {
1143               //((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, Standard_True);
1144
1145               // Check intersection between the triangle & sub-voxels of the voxel.
1146               for (Standard_Integer i = 0; i < 8; i++)
1147               {
1148                 ((Voxel_ROctBoolDS*) myVoxels)->GetCenter(ix, iy, iz, i, xc, yc, zc);
1149                 pc.SetCoord(xc, yc, zc);
1150                 if(TriBoxOverlap(p1, p2, p3, pc, extents2))
1151                 {
1152                   //((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, i, Standard_True);
1153
1154                   // Check intersection between the triangle & sub-voxels of the sub-voxel.
1155                   for (Standard_Integer j = 0; j < 8; j++)
1156                   {
1157                     ((Voxel_ROctBoolDS*) myVoxels)->GetCenter(ix, iy, iz, i, j, xc, yc, zc);
1158                     pc.SetCoord(xc, yc, zc);
1159                     if(TriBoxOverlap(p1, p2, p3, pc, extents4))
1160                     {
1161                       ((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, i, j, Standard_True);
1162                     }
1163                   } // End of "Check level 2".
1164
1165                 }
1166               } // End of "Check level 1".
1167               break;
1168             }
1169           }
1170         }
1171       }
1172     }
1173   }
1174 }