0023894: Voxel_BooleanOperation (Cut) gives incorrect results
authorPawel <pawel-kowalski@wp.pl>
Thu, 13 Jun 2013 11:19:29 +0000 (15:19 +0400)
committerPawel <pawel-kowalski@wp.pl>
Thu, 13 Jun 2013 11:19:29 +0000 (15:19 +0400)
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

14 files changed:
src/Voxel/Voxel_FastConverter.cdl
src/Voxel/Voxel_FastConverter.cxx
tests/v3d/grids.list
tests/v3d/voxel/A1 [new file with mode: 0644]
tests/v3d/voxel/A2 [new file with mode: 0644]
tests/v3d/voxel/A3 [new file with mode: 0644]
tests/v3d/voxel/A4 [new file with mode: 0644]
tests/v3d/voxel/A5 [new file with mode: 0644]
tests/v3d/voxel/A6 [new file with mode: 0644]
tests/v3d/voxel/A7 [new file with mode: 0644]
tests/v3d/voxel/A8 [new file with mode: 0644]
tests/v3d/voxel/A9 [new file with mode: 0644]
tests/v3d/voxel/B1 [new file with mode: 0644]
tests/v3d/voxel/B2 [new file with mode: 0644]

index 14b4f02..ae95d59 100755 (executable)
@@ -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;
index 6d5b35a..9ed4411 100755 (executable)
@@ -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(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;
+            }
+          }
+        }
+      }
+    }
+  }
+}
index 5d5c12c..f80c845 100644 (file)
@@ -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 (file)
index 0000000..9ad45e1
--- /dev/null
@@ -0,0 +1 @@
+voxelboolds
diff --git a/tests/v3d/voxel/A2 b/tests/v3d/voxel/A2
new file mode 100644 (file)
index 0000000..75607d0
--- /dev/null
@@ -0,0 +1 @@
+voxelcolords
diff --git a/tests/v3d/voxel/A3 b/tests/v3d/voxel/A3
new file mode 100644 (file)
index 0000000..04941bd
--- /dev/null
@@ -0,0 +1 @@
+voxelfloatds
diff --git a/tests/v3d/voxel/A4 b/tests/v3d/voxel/A4
new file mode 100644 (file)
index 0000000..f032943
--- /dev/null
@@ -0,0 +1 @@
+voxeloctboolds
diff --git a/tests/v3d/voxel/A5 b/tests/v3d/voxel/A5
new file mode 100644 (file)
index 0000000..e9b030c
--- /dev/null
@@ -0,0 +1 @@
+voxelroctboolds
diff --git a/tests/v3d/voxel/A6 b/tests/v3d/voxel/A6
new file mode 100644 (file)
index 0000000..16b89e1
--- /dev/null
@@ -0,0 +1 @@
+voxelfuseboolds
diff --git a/tests/v3d/voxel/A7 b/tests/v3d/voxel/A7
new file mode 100644 (file)
index 0000000..3ecc2fb
--- /dev/null
@@ -0,0 +1 @@
+voxelfusecolords
diff --git a/tests/v3d/voxel/A8 b/tests/v3d/voxel/A8
new file mode 100644 (file)
index 0000000..03b7baf
--- /dev/null
@@ -0,0 +1 @@
+voxelfusefloatds
diff --git a/tests/v3d/voxel/A9 b/tests/v3d/voxel/A9
new file mode 100644 (file)
index 0000000..94112ae
--- /dev/null
@@ -0,0 +1 @@
+voxelcutboolds
diff --git a/tests/v3d/voxel/B1 b/tests/v3d/voxel/B1
new file mode 100644 (file)
index 0000000..fa9616e
--- /dev/null
@@ -0,0 +1 @@
+voxelcutcolords
diff --git a/tests/v3d/voxel/B2 b/tests/v3d/voxel/B2
new file mode 100644 (file)
index 0000000..9280088
--- /dev/null
@@ -0,0 +1 @@
+voxelcutfloatds