OpenCascade/occt
1
0
Fork 0
mirror of https://git.dev.opencascade.org/repos/occt.git synced 2025-04-05 18:16:23 +03:00
Code Packages Releases Activity

0023894: Voxel_BooleanOperation (Cut) gives incorrect results

Browse Source 
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
This commit is contained in:
Pawel 2013-06-13 15:19:29 +04:00
parent 3830895866
commit f03671a0c4
14 changed files with 380 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
src/Voxel/Voxel_FastConverter.cdl
View File

@ -85,6 +85,16 @@ is
-- 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;
ithread : Integer from Standard = 1)
@ -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;

342
src/Voxel/Voxel_FastConverter.cxx
View File

@ -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;
}
}
}
}
}
}
}

1
tests/v3d/grids.list
View File

@ -9,3 +9,4 @@
009 vertex_wire
010 wire
011 wire_solid
012 voxel

1
tests/v3d/voxel/A1 Normal file
View File

@ -0,0 +1 @@
voxelboolds

1
tests/v3d/voxel/A2 Normal file
View File

@ -0,0 +1 @@
voxelcolords

1
tests/v3d/voxel/A3 Normal file
View File

@ -0,0 +1 @@
voxelfloatds

1
tests/v3d/voxel/A4 Normal file
View File

@ -0,0 +1 @@
voxeloctboolds

1
tests/v3d/voxel/A5 Normal file
View File

@ -0,0 +1 @@
voxelroctboolds

1
tests/v3d/voxel/A6 Normal file
View File

@ -0,0 +1 @@
voxelfuseboolds

1
tests/v3d/voxel/A7 Normal file
View File

@ -0,0 +1 @@
voxelfusecolords

1
tests/v3d/voxel/A8 Normal file
View File

@ -0,0 +1 @@
voxelfusefloatds

1
tests/v3d/voxel/A9 Normal file
View File

@ -0,0 +1 @@
voxelcutboolds

1
tests/v3d/voxel/B1 Normal file
View File

@ -0,0 +1 @@
voxelcutcolords

1
tests/v3d/voxel/B2 Normal file
View File

@ -0,0 +1 @@
voxelcutfloatds