0023934: Compiler warnings in MS VC++ 10
[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 #include <BRepClass3d_SolidClassifier.hxx>
39
40 // Printing the progress in stdout.
41 //#define CONV_DUMP
42
43 Voxel_FastConverter::Voxel_FastConverter(const TopoDS_Shape&    shape,
44                                          Voxel_ROctBoolDS&      voxels, 
45                                          const Standard_Real    deflection,
46                                          const Standard_Integer nbx,
47                                          const Standard_Integer nby,
48                                          const Standard_Integer nbz,
49                                          const Standard_Integer nbthreads,
50                                          const Standard_Boolean useExistingTriangulation)
51 :myShape(shape),myVoxels(&voxels),
52  myDeflection(deflection),
53  myNbX(nbx),myNbY(nby),myNbZ(nbz),
54  myNbThreads(nbthreads),myIsBool(2),
55  myNbTriangles(0),
56  myUseExistingTriangulation(useExistingTriangulation)
57 {
58   Init();
59 }
60
61 Voxel_FastConverter::Voxel_FastConverter(const TopoDS_Shape&    shape,
62                                          Voxel_BoolDS&          voxels, 
63                                          const Standard_Real    deflection,
64                                          const Standard_Integer nbx,
65                                          const Standard_Integer nby,
66                                          const Standard_Integer nbz,
67                                          const Standard_Integer nbthreads,
68                                          const Standard_Boolean useExistingTriangulation)
69 :myShape(shape),myVoxels(&voxels),
70  myDeflection(deflection),
71  myNbX(nbx),myNbY(nby),myNbZ(nbz),
72  myNbThreads(nbthreads),myIsBool(1),
73  myNbTriangles(0),
74  myUseExistingTriangulation(useExistingTriangulation)
75 {
76   Init();
77 }
78
79 Voxel_FastConverter::Voxel_FastConverter(const TopoDS_Shape&    shape,
80                                          Voxel_ColorDS&         voxels, 
81                                          const Standard_Real    deflection,
82                                          const Standard_Integer nbx,
83                                          const Standard_Integer nby,
84                                          const Standard_Integer nbz,
85                                          const Standard_Integer nbthreads,
86                                          const Standard_Boolean useExistingTriangulation)
87 :myShape(shape),myVoxels(&voxels),
88  myDeflection(deflection),
89  myNbX(nbx),myNbY(nby),myNbZ(nbz),
90  myNbThreads(nbthreads),myIsBool(0),
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::Mesh(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 v2dct(p2dc, p2dt);
817                   dv = v2dct.Magnitude() - Precision::Confusion();
818                   Lv.SetLocation(p2dc);
819                   Lv.SetDirection(v2dct);
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 v2dct(p2dc, p2dt);
885                       dv = v2dct.Magnitude() - Precision::Confusion();
886                       Lv.SetLocation(p2dc);
887                       Lv.SetDirection(v2dct);
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 min,max;
1073   // Find min, max of the triangle in x-direction, and test for overlap in X
1074   FINDMINMAX(v0.X(), v1.X(), v2.X(), min, max);
1075   if(min>extents.X() || max<-extents.X()) return false;
1076   
1077   // Same for Y
1078   FINDMINMAX(v0.Y(), v1.Y(), v2.Y(), min, max);
1079   if(min>extents.Y() || max<-extents.Y()) return false;
1080   
1081   // Same for Z
1082   FINDMINMAX(v0.Z(), v1.Z(), v2.Z(), min, max);
1083   if(min>extents.Z() || max<-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 }