From f03671a0c4b49105d210d344e0c01f6268585fb5 Mon Sep 17 00:00:00 2001 From: Pawel Date: Thu, 13 Jun 2013 15:19:29 +0400 Subject: [PATCH] 0023894: Voxel_BooleanOperation (Cut) gives incorrect results Implemented a new method using separating axis theorem to compute triangle-box intersection. Based on the intersection result the decision whether to set the voxel is made. Adjustment of lines (removal of extra-spaces). Adding test cases for voxel testing --- src/Voxel/Voxel_FastConverter.cdl | 26 +++ src/Voxel/Voxel_FastConverter.cxx | 342 ++++++++++++++++++++++++++++++ tests/v3d/grids.list | 1 + tests/v3d/voxel/A1 | 1 + tests/v3d/voxel/A2 | 1 + tests/v3d/voxel/A3 | 1 + tests/v3d/voxel/A4 | 1 + tests/v3d/voxel/A5 | 1 + tests/v3d/voxel/A6 | 1 + tests/v3d/voxel/A7 | 1 + tests/v3d/voxel/A8 | 1 + tests/v3d/voxel/A9 | 1 + tests/v3d/voxel/B1 | 1 + tests/v3d/voxel/B2 | 1 + 14 files changed, 380 insertions(+) create mode 100644 tests/v3d/voxel/A1 create mode 100644 tests/v3d/voxel/A2 create mode 100644 tests/v3d/voxel/A3 create mode 100644 tests/v3d/voxel/A4 create mode 100644 tests/v3d/voxel/A5 create mode 100644 tests/v3d/voxel/A6 create mode 100644 tests/v3d/voxel/A7 create mode 100644 tests/v3d/voxel/A8 create mode 100644 tests/v3d/voxel/A9 create mode 100644 tests/v3d/voxel/B1 create mode 100644 tests/v3d/voxel/B2 diff --git a/src/Voxel/Voxel_FastConverter.cdl b/src/Voxel/Voxel_FastConverter.cdl index 14b4f020ac..ae95d59c4d 100755 --- a/src/Voxel/Voxel_FastConverter.cdl +++ b/src/Voxel/Voxel_FastConverter.cdl @@ -84,6 +84,16 @@ is -- "ithread" is the index of the thread for current call of ::Convert(). -- Start numeration of "ithread" with 1, please. returns Boolean from Standard; + + ConvertUsingSAT(me : in out; + progress : out Integer from Standard; + ithread : Integer from Standard = 1) + ---Purpose: Converts a shape into a voxel representation using separating axis theorem. + -- It sets to 0 the outside volume of the shape and + -- 1 for surfacic part of the shape. + -- "ithread" is the index of the thread for current call of ::Convert(). + -- Start numeration of "ithread" with 1, please. + returns Boolean from Standard; FillInVolume(me : in out; inner : Byte from Standard; @@ -128,6 +138,22 @@ is izmax : Integer from Standard) is private; + ComputeVoxelsNearTriangle(me; + + p1 : Pnt from gp; + p2 : Pnt from gp; + p3 : Pnt from gp; + extents : Pnt from gp; + extents2 : Pnt from gp; + extents4 : Pnt from gp; + ixmin : Integer from Standard; + iymin : Integer from Standard; + izmin : Integer from Standard; + ixmax : Integer from Standard; + iymax : Integer from Standard; + izmax : Integer from Standard) + is private; + fields myShape : Shape from TopoDS; diff --git a/src/Voxel/Voxel_FastConverter.cxx b/src/Voxel/Voxel_FastConverter.cxx index 6d5b35a728..9ed441114f 100755 --- a/src/Voxel/Voxel_FastConverter.cxx +++ b/src/Voxel/Voxel_FastConverter.cxx @@ -287,6 +287,128 @@ Standard_Boolean Voxel_FastConverter::Convert(Standard_Integer& progress, 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) { @@ -706,3 +828,223 @@ void Voxel_FastConverter::ComputeVoxelsNearTriangle(const gp_Pln& plane } } } + +//! This macro quickly finds the min & max values among 3 variables +#define FINDMINMAX(x0, x1, x2, min, max) \ + min = max = x0; \ + if(x1max) max=x1; \ + if(x2max) 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; + } + } + } + } + } + } +} diff --git a/tests/v3d/grids.list b/tests/v3d/grids.list index 5d5c12ca97..f80c8451c9 100644 --- a/tests/v3d/grids.list +++ b/tests/v3d/grids.list @@ -9,3 +9,4 @@ 009 vertex_wire 010 wire 011 wire_solid +012 voxel diff --git a/tests/v3d/voxel/A1 b/tests/v3d/voxel/A1 new file mode 100644 index 0000000000..9ad45e1cee --- /dev/null +++ b/tests/v3d/voxel/A1 @@ -0,0 +1 @@ +voxelboolds diff --git a/tests/v3d/voxel/A2 b/tests/v3d/voxel/A2 new file mode 100644 index 0000000000..75607d0bcd --- /dev/null +++ b/tests/v3d/voxel/A2 @@ -0,0 +1 @@ +voxelcolords diff --git a/tests/v3d/voxel/A3 b/tests/v3d/voxel/A3 new file mode 100644 index 0000000000..04941bd170 --- /dev/null +++ b/tests/v3d/voxel/A3 @@ -0,0 +1 @@ +voxelfloatds diff --git a/tests/v3d/voxel/A4 b/tests/v3d/voxel/A4 new file mode 100644 index 0000000000..f0329433a3 --- /dev/null +++ b/tests/v3d/voxel/A4 @@ -0,0 +1 @@ +voxeloctboolds diff --git a/tests/v3d/voxel/A5 b/tests/v3d/voxel/A5 new file mode 100644 index 0000000000..e9b030c546 --- /dev/null +++ b/tests/v3d/voxel/A5 @@ -0,0 +1 @@ +voxelroctboolds diff --git a/tests/v3d/voxel/A6 b/tests/v3d/voxel/A6 new file mode 100644 index 0000000000..16b89e1b74 --- /dev/null +++ b/tests/v3d/voxel/A6 @@ -0,0 +1 @@ +voxelfuseboolds diff --git a/tests/v3d/voxel/A7 b/tests/v3d/voxel/A7 new file mode 100644 index 0000000000..3ecc2fb063 --- /dev/null +++ b/tests/v3d/voxel/A7 @@ -0,0 +1 @@ +voxelfusecolords diff --git a/tests/v3d/voxel/A8 b/tests/v3d/voxel/A8 new file mode 100644 index 0000000000..03b7baf98e --- /dev/null +++ b/tests/v3d/voxel/A8 @@ -0,0 +1 @@ +voxelfusefloatds diff --git a/tests/v3d/voxel/A9 b/tests/v3d/voxel/A9 new file mode 100644 index 0000000000..94112ae079 --- /dev/null +++ b/tests/v3d/voxel/A9 @@ -0,0 +1 @@ +voxelcutboolds diff --git a/tests/v3d/voxel/B1 b/tests/v3d/voxel/B1 new file mode 100644 index 0000000000..fa9616e5ac --- /dev/null +++ b/tests/v3d/voxel/B1 @@ -0,0 +1 @@ +voxelcutcolords diff --git a/tests/v3d/voxel/B2 b/tests/v3d/voxel/B2 new file mode 100644 index 0000000000..9280088cca --- /dev/null +++ b/tests/v3d/voxel/B2 @@ -0,0 +1 @@ +voxelcutfloatds