1 // Created on: 2008-05-30
2 // Created by: Vladislav ROMASHKO
3 // Copyright (c) 2008-2012 OPEN CASCADE SAS
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.
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.
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.
21 #include <Voxel_FastConverter.ixx>
23 #include <Bnd_Box.hxx>
24 #include <BRep_Tool.hxx>
25 #include <BRepBndLib.hxx>
26 #include <BRepMesh.hxx>
29 #include <TopoDS_Face.hxx>
30 #include <TopExp_Explorer.hxx>
32 #include <gp_Lin2d.hxx>
33 #include <gce_MakePln.hxx>
36 #include <Poly_Triangulation.hxx>
37 #include <IntAna2d_AnaIntersection.hxx>
38 #include <BRepClass3d_SolidClassifier.hxx>
40 // Printing the progress in stdout.
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 myIsBool(2),myNbThreads(nbthreads),
56 myUseExistingTriangulation(useExistingTriangulation)
61 Voxel_FastConverter::Voxel_FastConverter(const TopoDS_Shape& shape,
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 myIsBool(1),myNbThreads(nbthreads),
74 myUseExistingTriangulation(useExistingTriangulation)
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 myIsBool(0),myNbThreads(nbthreads),
92 myUseExistingTriangulation(useExistingTriangulation)
97 void Voxel_FastConverter::Init()
104 // Check number of splits.
105 Voxel_DS* voxels = (Voxel_DS*) myVoxels;
106 if (voxels->GetNbX() != myNbX || voxels->GetNbY() != myNbY || voxels->GetNbZ() != myNbZ)
108 // Compute boundary box of the shape
110 BRepBndLib::Add(myShape, box);
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);
116 // Initialize the voxels.
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);
125 // Check presence of triangulation.
127 Standard_Boolean triangulate = Standard_False;
128 TopExp_Explorer expl(myShape, TopAbs_FACE);
129 if(myUseExistingTriangulation == Standard_False)
131 for (; expl.More(); expl.Next())
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))
137 triangulate = Standard_True;
143 // Re-create the triangulation.
146 BRepMesh::Mesh(myShape, myDeflection);
149 // Compute the number of triangles.
151 expl.Init(myShape, TopAbs_FACE);
152 for (; expl.More(); expl.Next())
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();
162 void Voxel_FastConverter::Destroy()
167 Standard_Boolean Voxel_FastConverter::Convert(Standard_Integer& progress,
168 const Standard_Integer ithread)
174 printf("Progress = %d \r", progress);
177 if (myNbX <= 0 || myNbY <= 0 || myNbZ <= 0)
178 return Standard_False;
180 if(myNbTriangles == 0)
181 return Standard_False;
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);
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)
196 return Standard_False;
197 //in case we're in thread one process all triangles
201 div_t division = div(myNbTriangles, myNbThreads);
202 start_thread_triangle = (ithread - 1) * division.quot + 1;
203 end_thread_triangle = (ithread - 0) * division.quot;
205 if(ithread == myNbThreads)
206 end_thread_triangle += division.rem;
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())
218 const TopoDS_Face & F = TopoDS::Face(expl.Current());
219 Handle(Poly_Triangulation) T = BRep_Tool::Triangulation(F, L);
224 Standard_Boolean transform = !L.IsIdentity();
226 trsf = L.Transformation();
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++)
234 if (ithread_triangle < start_thread_triangle )
236 if (ithread_triangle > end_thread_triangle)
242 printf("Progress = %d \r", progress);
244 return Standard_True;
247 const Poly_Triangle& t = triangles.Value(itriangle);
249 gp_Pnt p1 = nodes.Value(n1);
250 gp_Pnt p2 = nodes.Value(n2);
251 gp_Pnt p3 = nodes.Value(n3);
259 // Get boundary box of the triangle
260 GetBndBox(p1, p2, p3, xmin, ymin, zmin, xmax, ymax, zmax);
262 // Find the range of voxels inside the boudary box of the triangle.
263 if (!ds->GetVoxel(xmin, ymin, zmin, ixmin, iymin, izmin))
265 if (!ds->GetVoxel(xmax, ymax, zmax, ixmax, iymax, izmax))
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())
275 gp_Pln plane = mkPlane.Value();
276 ComputeVoxelsNearTriangle(plane, p1, p2, p3, hdiagonal, ixmin, iymin, izmin, ixmax, iymax, izmax);
282 progress = (Standard_Integer) ( (Standard_Real) iprogress / (Standard_Real) myNbTriangles * 100.0 );
285 if (ithread == 1 && prev_progress != progress)
287 printf("Progress = %d \r", progress);
288 prev_progress = progress;
292 } // iteration of triangles
293 } // iteration of faces
299 printf("Progress = %d \r", progress);
301 return Standard_True;
304 Standard_Boolean Voxel_FastConverter::ConvertUsingSAT(Standard_Integer& progress,
305 const Standard_Integer ithread)
311 printf("Progress = %d \r", progress);
314 if (myNbX <= 0 || myNbY <= 0 || myNbZ <= 0)
315 return Standard_False;
317 if(myNbTriangles == 0)
318 return Standard_False;
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);
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)
334 return Standard_False;
335 //in case we're in thread one process all triangles
339 div_t division = div(myNbTriangles, myNbThreads);
340 start_thread_triangle = (ithread - 1) * division.quot + 1;
341 end_thread_triangle = (ithread - 0) * division.quot;
343 if(ithread == myNbThreads)
344 end_thread_triangle += division.rem;
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())
356 const TopoDS_Face & F = TopoDS::Face(expl.Current());
357 Handle(Poly_Triangulation) T = BRep_Tool::Triangulation(F, L);
362 Standard_Boolean transform = !L.IsIdentity();
364 trsf = L.Transformation();
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++)
372 if (ithread_triangle < start_thread_triangle )
374 if (ithread_triangle > end_thread_triangle)
380 printf("Progress = %d \r", progress);
382 return Standard_True;
385 const Poly_Triangle& t = triangles.Value(itriangle);
388 gp_Pnt p1 = nodes.Value(n1);
389 gp_Pnt p2 = nodes.Value(n2);
390 gp_Pnt p3 = nodes.Value(n3);
398 // Get boundary box of the triangle
399 GetBndBox(p1, p2, p3, xmin, ymin, zmin, xmax, ymax, zmax);
401 // Find the range of voxels inside the boudary box of the triangle.
402 if (!ds->GetVoxel(xmin, ymin, zmin, ixmin, iymin, izmin))
404 if (!ds->GetVoxel(xmax, ymax, zmax, ixmax, iymax, izmax))
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);
417 progress = (Standard_Integer) ( (Standard_Real) iprogress / (Standard_Real) myNbTriangles * 100.0 );
420 if (ithread == 1 && prev_progress != progress)
422 printf("Progress = %d \r", progress);
423 prev_progress = progress;
427 } // iteration of triangles
428 } // iteration of faces
434 printf("Progress = %d \r", progress);
436 return Standard_True;
439 Standard_Boolean Voxel_FastConverter::FillInVolume(const Standard_Byte inner,
440 const Standard_Integer /*ithread*/)
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;
448 // Fill-in internal voxels by the value "inner"
449 for (ix = 0; ix < nbx; ix++)
451 for (iy = 0; iy < nby; iy++)
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++)
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)
466 prev_surface = surface;
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++)
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)
484 if (volume && !surface)
486 (myIsBool == 1) ? ((Voxel_BoolDS*)myVoxels)->Set(ix, iy, iz, inner) :
487 ((Voxel_ColorDS*)myVoxels)->Set(ix, iy, iz, inner);
489 prev_surface = surface;
496 // Set value of interbal voxels to 0 ("inner" = 0)
497 Standard_Boolean next_surface;
498 for (ix = 0; ix < nbx; ix++)
500 for (iy = 0; iy < nby; iy++)
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++)
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)
515 if (volume && iz + 1 < nbz)
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;
521 if (volume && prev_surface == surface && next_surface)
523 (myIsBool == 1) ? ((Voxel_BoolDS*)myVoxels)->Set(ix, iy, iz, inner) :
524 ((Voxel_ColorDS*)myVoxels)->Set(ix, iy, iz, inner);
526 prev_surface = surface;
532 return Standard_True;
535 Standard_Boolean Voxel_FastConverter::FillInVolume(const Standard_Byte inner, const TopoDS_Shape & shape, const Standard_Integer /*ithread*/)
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;
541 BRepClass3d_SolidClassifier solidClassifier(shape);
542 Standard_Real xc, yc, zc;
546 // Fill-in internal voxels by the value "inner"
547 for (ix = 0; ix < nbx; ix++)
549 for (iy = 0; iy < nby; iy++)
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++)
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)
563 if(isOnVerticalSurface)
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());
570 if(solidClassifier.State() == TopAbs_IN)
571 volume = Standard_True;
573 volume = Standard_False;
578 if(prev_surface && surface)
579 isOnVerticalSurface = Standard_True;
581 isOnVerticalSurface = Standard_False;
582 prev_surface = surface;
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++)
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)
599 if(isOnVerticalSurface)
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());
606 if(solidClassifier.State() == TopAbs_IN)
607 volume = Standard_True;
609 volume = Standard_False;
614 if (volume && !surface)
616 (myIsBool == 1) ? ((Voxel_BoolDS*)myVoxels)->Set(ix, iy, iz, inner) :
617 ((Voxel_ColorDS*)myVoxels)->Set(ix, iy, iz, inner);
619 if(prev_surface && surface)
620 isOnVerticalSurface = Standard_True;
622 isOnVerticalSurface = Standard_False;
623 prev_surface = surface;
629 return Standard_True;
632 void Voxel_FastConverter::GetBndBox(const gp_Pnt& p1,
640 Standard_Real& zmax) const
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)
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;
692 void Voxel_FastConverter::ComputeVoxelsNearTriangle(const gp_Pln& plane,
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
705 Standard_Real xc, yc, zc, uc, vc, u1, v1, u2, v2, u3, v3;
706 Standard_Integer ix, iy, iz;
707 IntAna2d_AnaIntersection intersector2d;
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);
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;
722 Voxel_DS* ds = (Voxel_DS*) myVoxels;
723 for (ix = ixmin; ix <= ixmax; ix++)
725 for (iy = iymin; iy <= iymax; iy++)
727 for (iz = izmin; iz <= izmax; iz++)
729 ds->GetCenter(ix, iy, iz, xc, yc, zc);
730 pc.SetCoord(xc, yc, zc);
731 if (plane.Distance(pc) < hdiagonal)
733 ElSLib::Parameters(plane, pc, uc, vc);
734 p2dc.SetCoord(uc, vc);
736 gp_Vec2d v2dct(p2dc, p2dt);
737 dv = v2dct.Magnitude() - Precision::Confusion();
738 Lv.SetLocation(p2dc);
739 Lv.SetDirection(v2dct);
742 if (mayIntersect(p2d1, p2d2, p2dc, p2dt))
744 intersector2d.Perform(Lv, L1);
745 if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
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)
759 if (mayIntersect(p2d2, p2d3, p2dc, p2dt))
761 intersector2d.Perform(Lv, L2);
762 if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
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)
776 if (mayIntersect(p2d3, p2d1, p2dc, p2dt))
778 intersector2d.Perform(Lv, L3);
779 if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
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)
792 // Set positive value to this voxel:
796 ((Voxel_ColorDS*) myVoxels)->Set(ix, iy, iz, 15);
799 ((Voxel_BoolDS*) myVoxels)->Set(ix, iy, iz, Standard_True);
803 //((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, Standard_True);
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++)
809 ((Voxel_ROctBoolDS*) myVoxels)->GetCenter(ix, iy, iz, i, xc, yc, zc);
810 pc.SetCoord(xc, yc, zc);
811 if (plane.Distance(pc) < hdiagonal2)
813 ElSLib::Parameters(plane, pc, uc, vc);
814 p2dc.SetCoord(uc, vc);
816 gp_Vec2d v2dct(p2dc, p2dt);
817 dv = v2dct.Magnitude() - Precision::Confusion();
818 Lv.SetLocation(p2dc);
819 Lv.SetDirection(v2dct);
822 if (mayIntersect(p2d1, p2d2, p2dc, p2dt))
824 intersector2d.Perform(Lv, L1);
825 if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
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)
839 if (mayIntersect(p2d2, p2d3, p2dc, p2dt))
841 intersector2d.Perform(Lv, L2);
842 if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
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)
856 if (mayIntersect(p2d3, p2d1, p2dc, p2dt))
858 intersector2d.Perform(Lv, L3);
859 if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
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)
872 //((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, i, Standard_True);
874 // Check intersection between the triangle & sub-voxels of the sub-voxel.
875 for (Standard_Integer j = 0; j < 8; j++)
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)
881 ElSLib::Parameters(plane, pc, uc, vc);
882 p2dc.SetCoord(uc, vc);
884 gp_Vec2d v2dct(p2dc, p2dt);
885 dv = v2dct.Magnitude() - Precision::Confusion();
886 Lv.SetLocation(p2dc);
887 Lv.SetDirection(v2dct);
890 if (mayIntersect(p2d1, p2d2, p2dc, p2dt))
892 intersector2d.Perform(Lv, L1);
893 if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
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)
907 if (mayIntersect(p2d2, p2d3, p2dc, p2dt))
909 intersector2d.Perform(Lv, L2);
910 if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
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)
924 if (mayIntersect(p2d3, p2d1, p2dc, p2dt))
926 intersector2d.Perform(Lv, L3);
927 if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints())
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)
940 ((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, i, j, Standard_True);
942 } // End of "Check level 2".
945 } // End of "Check level 1".
956 //! This macro quickly finds the min & max values among 3 variables
957 #define FINDMINMAX(x0, x1, x2, min, max) \
964 static bool planeBoxOverlap(const gp_Vec & normal, const double d, const gp_Pnt & maxbox)
967 if(normal.X() > 0.0) { vmin.SetX(-maxbox.X()); vmax.SetX(maxbox.X()); }
968 else { vmin.SetX(maxbox.X()); vmax.SetX(-maxbox.X()); }
970 if(normal.Y() > 0.0) { vmin.SetY(-maxbox.Y()); vmax.SetY(maxbox.Y()); }
971 else { vmin.SetY(maxbox.Y()); vmax.SetY(-maxbox.Y()); }
973 if(normal.Z() > 0.0) { vmin.SetZ(-maxbox.Z()); vmax.SetZ(maxbox.Z()); }
974 else { vmin.SetZ(maxbox.Z()); vmax.SetZ(-maxbox.Z()); }
976 if((normal.Dot(vmin)) + d > 0.0) return false;
977 if((normal.Dot(vmax)) + d>= 0.0) return true;
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;
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;
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;
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;
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;
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;
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 \
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); \
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); \
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);
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)
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
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);
1071 // First, test overlap in the {x,y,z}-directions
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;
1078 FINDMINMAX(v0.Y(), v1.Y(), v2.Y(), min, max);
1079 if(min>extents.Y() || max<-extents.Y()) return false;
1082 FINDMINMAX(v0.Z(), v1.Z(), v2.Z(), min, max);
1083 if(min>extents.Z() || max<-extents.Z()) return false;
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;
1094 // 3) "Class III" tests
1095 //if(mFullPrimBoxTest)
1097 IMPLEMENT_CLASS3_TESTS
1103 void Voxel_FastConverter::ComputeVoxelsNearTriangle(const gp_Pnt& p1,
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
1117 Standard_Real xc, yc, zc;
1118 Standard_Integer ix, iy, iz;
1120 Voxel_DS* ds = (Voxel_DS*) myVoxels;
1121 for (ix = ixmin; ix <= ixmax; ix++)
1123 for (iy = iymin; iy <= iymax; iy++)
1125 for (iz = izmin; iz <= izmax; iz++)
1127 ds->GetCenter(ix, iy, iz, xc, yc, zc);
1128 pc.SetCoord(xc, yc, zc);
1130 if(TriBoxOverlap(p1, p2, p3, pc, extents))
1132 // Set positive value to this voxel:
1136 ((Voxel_ColorDS*) myVoxels)->Set(ix, iy, iz, 15);
1139 ((Voxel_BoolDS*) myVoxels)->Set(ix, iy, iz, Standard_True);
1143 //((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, Standard_True);
1145 // Check intersection between the triangle & sub-voxels of the voxel.
1146 for (Standard_Integer i = 0; i < 8; i++)
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))
1152 //((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, i, Standard_True);
1154 // Check intersection between the triangle & sub-voxels of the sub-voxel.
1155 for (Standard_Integer j = 0; j < 8; j++)
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))
1161 ((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, i, j, Standard_True);
1163 } // End of "Check level 2".
1166 } // End of "Check level 1".