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