return Standard_True;
}
+Standard_Boolean Voxel_FastConverter::ConvertUsingSAT(Standard_Integer& progress,
+ const Standard_Integer ithread)
+{
+ if (ithread == 1)
+ progress = 0;
+#ifdef CONV_DUMP
+ if (ithread == 1)
+ printf("Progress = %d \r", progress);
+#endif
+
+ if (myNbX <= 0 || myNbY <= 0 || myNbZ <= 0)
+ return Standard_False;
+
+ if(myNbTriangles == 0)
+ return Standard_False;
+
+ // Half of size of a voxel (also for Voxel_ROctBoolDS)
+ Voxel_DS* ds = (Voxel_DS*) myVoxels;
+ Standard_Real dx = ds->GetXLen() / (Standard_Real) ds->GetNbX(),
+ dy = ds->GetYLen() / (Standard_Real) ds->GetNbY(),
+ dz = ds->GetZLen() / (Standard_Real) ds->GetNbZ();
+ gp_Pnt extents(dx/2.0, dy/2.0, dz/2.0);
+ gp_Pnt extents2(dx/4.0, dy/4.0, dz/4.0);
+ gp_Pnt extents4(dx/8.0, dy/8.0, dz/8.0);
+
+ // Compute the scope of triangles for current thread
+ Standard_Integer start_thread_triangle = 1, end_thread_triangle = myNbTriangles, ithread_triangle = 0;
+ start_thread_triangle = (ithread - 1) * (myNbTriangles / myNbThreads) + 1;
+ end_thread_triangle = (ithread - 0) * (myNbTriangles / myNbThreads);
+
+ // Convert
+ TopLoc_Location L;
+ Standard_Integer iprogress = 0, prev_progress = 0;
+ Standard_Integer n1, n2, n3;
+ Standard_Integer ixmin, iymin, izmin, ixmax, iymax, izmax;
+ Standard_Real xmin, ymin, zmin, xmax, ymax, zmax;
+ TopExp_Explorer expl(myShape, TopAbs_FACE);
+ for (; expl.More(); expl.Next())
+ {
+ const TopoDS_Face & F = TopoDS::Face(expl.Current());
+ Handle(Poly_Triangulation) T = BRep_Tool::Triangulation(F, L);
+ if (T.IsNull())
+ continue;
+
+ gp_Trsf trsf;
+ Standard_Boolean transform = !L.IsIdentity();
+ if (transform)
+ trsf = L.Transformation();
+
+ const TColgp_Array1OfPnt& nodes = T->Nodes();
+ const Poly_Array1OfTriangle& triangles = T->Triangles();
+ Standard_Integer itriangle = triangles.Lower(), nb_triangles = triangles.Upper();
+ for (; itriangle <= nb_triangles; itriangle++)
+ {
+ ithread_triangle++;
+ if (ithread_triangle < start_thread_triangle )
+ continue;
+ if (ithread_triangle > end_thread_triangle)
+ {
+ if (ithread == 1)
+ progress = 100;
+#ifdef CONV_DUMP
+ if (ithread == 1)
+ printf("Progress = %d \r", progress);
+#endif
+ return Standard_True;
+ }
+
+ const Poly_Triangle& t = triangles.Value(itriangle);
+
+ t.Get(n1, n2, n3);
+ gp_Pnt p1 = nodes.Value(n1);
+ gp_Pnt p2 = nodes.Value(n2);
+ gp_Pnt p3 = nodes.Value(n3);
+ if (transform)
+ {
+ p1.Transform(trsf);
+ p2.Transform(trsf);
+ p3.Transform(trsf);
+ }
+
+ // Get boundary box of the triangle
+ GetBndBox(p1, p2, p3, xmin, ymin, zmin, xmax, ymax, zmax);
+
+ // Find the range of voxels inside the boudary box of the triangle.
+ if (!ds->GetVoxel(xmin, ymin, zmin, ixmin, iymin, izmin))
+ continue;
+ if (!ds->GetVoxel(xmax, ymax, zmax, ixmax, iymax, izmax))
+ continue;
+
+ // Perform triangle-box intersection to find the voxels resulting from the processed triangle.;
+ // Using SAT theorem to quickly find the intersection.
+ ComputeVoxelsNearTriangle(p1, p2, p3,
+ extents, extents2, extents4,
+ ixmin, iymin, izmin, ixmax, iymax, izmax);
+
+ // Progress
+ if (ithread == 1)
+ {
+ iprogress++;
+ progress = (Standard_Integer) ( (Standard_Real) iprogress / (Standard_Real) myNbTriangles * 100.0 );
+ }
+#ifdef CONV_DUMP
+ if (ithread == 1 && prev_progress != progress)
+ {
+ printf("Progress = %d \r", progress);
+ prev_progress = progress;
+ }
+#endif
+
+ } // iteration of triangles
+ } // iteration of faces
+
+ if (ithread == 1)
+ progress = 100;
+#ifdef CONV_DUMP
+ if (ithread == 1)
+ printf("Progress = %d \r", progress);
+#endif
+ return Standard_True;
+}
+
Standard_Boolean Voxel_FastConverter::FillInVolume(const Standard_Byte inner,
const Standard_Integer ithread)
{
}
}
}
+
+//! This macro quickly finds the min & max values among 3 variables
+#define FINDMINMAX(x0, x1, x2, min, max) \
+ min = max = x0; \
+ if(x1<min) min=x1; \
+ if(x1>max) max=x1; \
+ if(x2<min) min=x2; \
+ if(x2>max) max=x2;
+
+static bool planeBoxOverlap(const gp_Vec & normal, const double d, const gp_Pnt & maxbox)
+{
+ gp_Vec vmin, vmax;
+ if(normal.X() > 0.0) { vmin.SetX(-maxbox.X()); vmax.SetX(maxbox.X()); }
+ else { vmin.SetX(maxbox.X()); vmax.SetX(-maxbox.X()); }
+
+ if(normal.Y() > 0.0) { vmin.SetY(-maxbox.Y()); vmax.SetY(maxbox.Y()); }
+ else { vmin.SetY(maxbox.Y()); vmax.SetY(-maxbox.Y()); }
+
+ if(normal.Z() > 0.0) { vmin.SetZ(-maxbox.Z()); vmax.SetZ(maxbox.Z()); }
+ else { vmin.SetZ(maxbox.Z()); vmax.SetZ(-maxbox.Z()); }
+
+ if((normal.Dot(vmin)) + d > 0.0) return false;
+ if((normal.Dot(vmax)) + d>= 0.0) return true;
+
+ return false;
+}
+
+#define AXISTEST_X01(a, b, fa, fb) \
+ min = a*v0.Y() - b*v0.Z(); \
+ max = a*v2.Y() - b*v2.Z(); \
+ if(min>max) {const double tmp=max; max=min; min=tmp; } \
+ rad = fa * extents.Y() + fb * extents.Z(); \
+ if(min>rad || max<-rad) return false;
+
+#define AXISTEST_X2(a, b, fa, fb) \
+ min = a*v0.Y() - b*v0.Z(); \
+ max = a*v1.Y() - b*v1.Z(); \
+ if(min>max) {const double tmp=max; max=min; min=tmp; } \
+ rad = fa * extents.Y() + fb * extents.Z(); \
+ if(min>rad || max<-rad) return false;
+
+#define AXISTEST_Y02(a, b, fa, fb) \
+ min = b*v0.Z() - a*v0.X(); \
+ max = b*v2.Z() - a*v2.X(); \
+ if(min>max) {const double tmp=max; max=min; min=tmp; } \
+ rad = fa * extents.X() + fb * extents.Z(); \
+ if(min>rad || max<-rad) return false;
+
+#define AXISTEST_Y1(a, b, fa, fb) \
+ min = b*v0.Z() - a*v0.X(); \
+ max = b*v1.Z() - a*v1.X(); \
+ if(min>max) {const double tmp=max; max=min; min=tmp; } \
+ rad = fa * extents.X() + fb * extents.Z(); \
+ if(min>rad || max<-rad) return false;
+
+#define AXISTEST_Z12(a, b, fa, fb) \
+ min = a*v1.X() - b*v1.Y(); \
+ max = a*v2.X() - b*v2.Y(); \
+ if(min>max) {const double tmp=max; max=min; min=tmp; } \
+ rad = fa * extents.X() + fb * extents.Y(); \
+ if(min>rad || max<-rad) return false;
+
+#define AXISTEST_Z0(a, b, fa, fb) \
+ min = a*v0.X() - b*v0.Y(); \
+ max = a*v1.X() - b*v1.Y(); \
+ if(min>max) {const double tmp=max; max=min; min=tmp; } \
+ rad = fa * extents.X() + fb * extents.Y(); \
+ if(min>rad || max<-rad) return false;
+
+// compute triangle edges
+// - edges lazy evaluated to take advantage of early exits
+// - fabs precomputed (half less work, possible since extents are always >0)
+// - customized macros to take advantage of the null component
+// - axis vector discarded, possibly saves useless movs
+#define IMPLEMENT_CLASS3_TESTS \
+ double rad; \
+ double min, max; \
+ \
+ const double fey0 = fabs(e0.Y()); \
+ const double fez0 = fabs(e0.Z()); \
+ AXISTEST_X01(e0.Z(), e0.Y(), fez0, fey0); \
+ const double fex0 = fabs(e0.X()); \
+ AXISTEST_Y02(e0.Z(), e0.X(), fez0, fex0); \
+ AXISTEST_Z12(e0.Y(), e0.X(), fey0, fex0); \
+ \
+ const double fey1 = fabs(e1.Y()); \
+ const double fez1 = fabs(e1.Z()); \
+ AXISTEST_X01(e1.Z(), e1.Y(), fez1, fey1); \
+ const double fex1 = fabs(e1.X()); \
+ AXISTEST_Y02(e1.Z(), e1.X(), fez1, fex1); \
+ AXISTEST_Z0(e1.Y(), e1.X(), fey1, fex1); \
+ \
+ const gp_Vec e2 = v2 - v0; \
+ const double fey2 = fabs(e2.Y()); \
+ const double fez2 = fabs(e2.Z()); \
+ AXISTEST_X2(e2.Z(), e2.Y(), fez2, fey2); \
+ const double fex2 = fabs(e2.X()); \
+ AXISTEST_Y1(e2.Z(), e2.X(), fez2, fex2); \
+ AXISTEST_Z12(e2.Y(), e2.X(), fey2, fex2);
+
+static bool TriBoxOverlap(const gp_Pnt & p1, const gp_Pnt & p2, const gp_Pnt & p3,
+ const gp_Pnt & center, const gp_Pnt & extents)
+{
+ // use separating axis theorem to test overlap between triangle and box
+ // need to test for overlap in these directions:
+ // 1) the {x,y,z}-directions (actually, since we use the AABB of the triangle
+ // we do not even need to test these)
+ // 2) normal of the triangle
+ // 3) crossproduct(edge from tri, {x,y,z}-directin)
+ // this gives 3x3=9 more tests
+
+ // move everything so that the boxcenter is in (0,0,0)
+ gp_Vec v0(center, p1);
+ gp_Vec v1(center, p2);
+ gp_Vec v2(center, p3);
+
+ // First, test overlap in the {x,y,z}-directions
+ double min,max;
+ // Find min, max of the triangle in x-direction, and test for overlap in X
+ FINDMINMAX(v0.X(), v1.X(), v2.X(), min, max);
+ if(min>extents.X() || max<-extents.X()) return false;
+
+ // Same for Y
+ FINDMINMAX(v0.Y(), v1.Y(), v2.Y(), min, max);
+ if(min>extents.Y() || max<-extents.Y()) return false;
+
+ // Same for Z
+ FINDMINMAX(v0.Z(), v1.Z(), v2.Z(), min, max);
+ if(min>extents.Z() || max<-extents.Z()) return false;
+
+ // 2) Test if the box intersects the plane of the triangle
+ // compute plane equation of triangle: normal*x+d=0
+ // ### could be precomputed since we use the same leaf triangle several times
+ const gp_Vec e0 = v1 - v0;
+ const gp_Vec e1 = v2 - v1;
+ const gp_Vec normal = e0.Crossed(e1);
+ const double d = -normal.Dot(v0);
+ if(!planeBoxOverlap(normal, d, extents)) return false;
+
+ // 3) "Class III" tests
+ //if(mFullPrimBoxTest)
+ {
+ IMPLEMENT_CLASS3_TESTS
+ }
+
+ return true;
+}
+
+void Voxel_FastConverter::ComputeVoxelsNearTriangle(const gp_Pnt& p1,
+ const gp_Pnt& p2,
+ const gp_Pnt& p3,
+ const gp_Pnt& extents,
+ const gp_Pnt& extents2,
+ const gp_Pnt& extents4,
+ const Standard_Integer ixmin,
+ const Standard_Integer iymin,
+ const Standard_Integer izmin,
+ const Standard_Integer ixmax,
+ const Standard_Integer iymax,
+ const Standard_Integer izmax) const
+{
+ gp_Pnt pc;
+ Standard_Real xc, yc, zc;
+ Standard_Integer ix, iy, iz;
+
+ Voxel_DS* ds = (Voxel_DS*) myVoxels;
+ for (ix = ixmin; ix <= ixmax; ix++)
+ {
+ for (iy = iymin; iy <= iymax; iy++)
+ {
+ for (iz = izmin; iz <= izmax; iz++)
+ {
+ ds->GetCenter(ix, iy, iz, xc, yc, zc);
+ pc.SetCoord(xc, yc, zc);
+
+ if(TriBoxOverlap(p1, p2, p3, pc, extents))
+ {
+ // Set positive value to this voxel:
+ switch (myIsBool)
+ {
+ case 0:
+ ((Voxel_ColorDS*) myVoxels)->Set(ix, iy, iz, 15);
+ break;
+ case 1:
+ ((Voxel_BoolDS*) myVoxels)->Set(ix, iy, iz, Standard_True);
+ break;
+ case 2:
+ {
+ //((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, Standard_True);
+
+ // Check intersection between the triangle & sub-voxels of the voxel.
+ for (Standard_Integer i = 0; i < 8; i++)
+ {
+ ((Voxel_ROctBoolDS*) myVoxels)->GetCenter(ix, iy, iz, i, xc, yc, zc);
+ pc.SetCoord(xc, yc, zc);
+ if(TriBoxOverlap(p1, p2, p3, pc, extents2))
+ {
+ //((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, i, Standard_True);
+
+ // Check intersection between the triangle & sub-voxels of the sub-voxel.
+ for (Standard_Integer j = 0; j < 8; j++)
+ {
+ ((Voxel_ROctBoolDS*) myVoxels)->GetCenter(ix, iy, iz, i, j, xc, yc, zc);
+ pc.SetCoord(xc, yc, zc);
+ if(TriBoxOverlap(p1, p2, p3, pc, extents4))
+ {
+ ((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, i, j, Standard_True);
+ }
+ } // End of "Check level 2".
+
+ }
+ } // End of "Check level 1".
+ break;
+ }
+ }
+ }
+ }
+ }
+ }
+}
projects / occt.git / commitdiff