-
Notifications
You must be signed in to change notification settings - Fork 52
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #44 from Yeicor/import-stl
Convert triangle meshes to SDF3 and voxel-based cache/smoothing
- Loading branch information
Showing
10 changed files
with
423 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,2 @@ | ||
TOP = ../.. | ||
include $(TOP)/mk/example.mk |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1 @@ | ||
d717263395e88df0a1f878b6e51021b02f0656ba monkey-out.stl |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,72 @@ | ||
//----------------------------------------------------------------------------- | ||
/* | ||
Imported monkey model, with modifications | ||
*/ | ||
//----------------------------------------------------------------------------- | ||
|
||
package main | ||
|
||
import ( | ||
"github.com/deadsy/sdfx/obj" | ||
"github.com/deadsy/sdfx/render" | ||
"github.com/deadsy/sdfx/sdf" | ||
"log" | ||
"os" | ||
"time" | ||
) | ||
|
||
func monkeyWithHat() sdf.SDF3 { | ||
// MONKEY | ||
// - Open the STL file | ||
file, err := os.OpenFile("monkey.stl", os.O_RDONLY, 0400) | ||
if err != nil { | ||
file, err = os.OpenFile("examples/monkey_hat/monkey.stl", os.O_RDONLY, 0400) | ||
if err != nil { | ||
panic(err) | ||
} | ||
} | ||
// - Create the SDF from the mesh (a modified Suzanne from Blender with 366 faces) | ||
monkeyImported, err := obj.ImportSTL(file, 20, 3, 5) | ||
if err != nil { | ||
panic(err) | ||
} | ||
|
||
// HAT | ||
hatHeight := 0.5 | ||
hat, err := sdf.Cylinder3D(hatHeight, 0.6, 0) | ||
if err != nil { | ||
panic(err) | ||
} | ||
edge, err := sdf.Cylinder3D(hatHeight*0.4, 1, 0) | ||
if err != nil { | ||
panic(err) | ||
} | ||
edge = sdf.Transform3D(edge, sdf.Translate3d(sdf.V3{Z: -hatHeight / 2})) | ||
fullHat := sdf.Union3D(hat, edge) | ||
|
||
// Union | ||
fullHat = sdf.Transform3D(fullHat, sdf.Translate3d(sdf.V3{Y: 0.15, Z: 1})) | ||
monkeyHat := sdf.Union3D(monkeyImported, fullHat) | ||
|
||
// - Cache the mesh full SDF3 hierarchy for faster evaluation (at the cost of initialization time and memory). | ||
// It also smooths the mesh a little using trilinear interpolation. | ||
// It is actually slower for this mesh (unless meshCells <<< renderer's meshCells), but should be faster for | ||
// more complex meshes (with more triangles) or SDF3 hierarchies that take longer to evaluate. | ||
monkeyHat = sdf.NewVoxelSDF3(monkeyHat, 64, nil) // Use 32 for harder smoothing demo | ||
|
||
return monkeyHat | ||
} | ||
|
||
func main() { | ||
startTime := time.Now() | ||
monkeyHat := monkeyWithHat() | ||
|
||
render.ToSTL(monkeyHat, 128, "monkey-out.stl", &render.MarchingCubesUniform{}) | ||
|
||
// Dual Contouring is very sensitive to noise (produced when close to shared triangle vertices) | ||
//render.ToSTL(monkeyHat, 64, "monkey-out.stl", dc.NewDualContouringDefault()) | ||
|
||
log.Println("Monkey + hat rendered in", time.Since(startTime)) | ||
} |
Binary file not shown.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,210 @@ | ||
//----------------------------------------------------------------------------- | ||
/* | ||
Closed-surface triangle meshes (and STL files) | ||
*/ | ||
//----------------------------------------------------------------------------- | ||
|
||
package obj | ||
|
||
import ( | ||
"github.com/deadsy/sdfx/render" | ||
"github.com/deadsy/sdfx/sdf" | ||
"github.com/dhconnelly/rtreego" | ||
"github.com/hschendel/stl" | ||
"io" | ||
"math" | ||
) | ||
|
||
//----------------------------------------------------------------------------- | ||
|
||
type triMeshSdf struct { | ||
rtree *rtreego.Rtree | ||
numNeighbors int | ||
bb sdf.Box3 | ||
} | ||
|
||
const stlEpsilon = 1e-1 | ||
|
||
func (t *triMeshSdf) Evaluate(p sdf.V3) float64 { | ||
// Check all triangle distances | ||
signedDistanceResult := 1. | ||
closestTriangle := math.MaxFloat64 | ||
// Quickly skip checking most triangles by only checking the N closest neighbours (AABB based) | ||
neighbors := t.rtree.NearestNeighbors(t.numNeighbors, stlToPoint(p)) | ||
for _, neighbor := range neighbors { | ||
triangle := neighbor.(*stlTriangle).Triangle3 | ||
testPointToTriangle := p.Sub(triangle.V[0]) | ||
triNormal := triangle.Normal() | ||
signedDistanceToTriPlane := triNormal.Dot(testPointToTriangle) | ||
// Take this triangle as the source of truth if the projection of the point on the triangle is the closest | ||
distToTri, _ := stlPointToTriangleDistSq(p, triangle) | ||
if distToTri < closestTriangle { | ||
closestTriangle = distToTri | ||
signedDistanceResult = signedDistanceToTriPlane | ||
} | ||
} | ||
return signedDistanceResult | ||
} | ||
|
||
func (t *triMeshSdf) BoundingBox() sdf.Box3 { | ||
return t.bb | ||
} | ||
|
||
// ImportTriMesh converts a triangle-based mesh into a SDF3 surface. minChildren and maxChildren are parameters that can | ||
// affect the performance of the internal data structure (3 and 5 are a good default; maxChildren >= minChildren > 0). | ||
// | ||
// WARNING: Setting a low numNeighbors will consider many fewer triangles for each evaluated point, greatly speeding up | ||
// the algorithm. However, if the count of triangles is too low artifacts will appear on the surface (triangle | ||
// continuations). Setting this value to MaxInt is extremely slow but will provide correct results, so choose a value | ||
// that works for your model. | ||
// | ||
// It is recommended to cache (and/or smooth) its values by using sdf.VoxelSdf3. | ||
// | ||
// WARNING: It will only work on non-intersecting closed-surface(s) meshes. | ||
// NOTE: Fix using blender for intersecting surfaces: Edit mode > P > By loose parts > Add boolean modifier to join them | ||
func ImportTriMesh(tris chan *render.Triangle3, numNeighbors, minChildren, maxChildren int) sdf.SDF3 { | ||
m := &triMeshSdf{ | ||
rtree: nil, | ||
numNeighbors: numNeighbors, | ||
bb: sdf.Box3{ | ||
Min: sdf.V3{X: math.MaxFloat64, Y: math.MaxFloat64, Z: math.MaxFloat64}, | ||
Max: sdf.V3{X: -math.MaxFloat64, Y: -math.MaxFloat64, Z: -math.MaxFloat64}, | ||
}, | ||
} | ||
|
||
// Compute the bounding box | ||
bulkLoad := make([]rtreego.Spatial, 0) | ||
for triangle := range tris { | ||
bulkLoad = append(bulkLoad, &stlTriangle{Triangle3: triangle}) | ||
for _, vertex := range triangle.V { | ||
m.bb = m.bb.Include(vertex) | ||
} | ||
} | ||
if !m.bb.Contains(m.bb.Min) { // Return a valid bounding box if no vertices are found in the mesh | ||
m.bb = sdf.Box3{} // Empty box centered at {0,0,0} | ||
} | ||
//m.bb = m.bb.ScaleAboutCenter(1 + 1e-12) // Avoids missing faces due to inaccurate math operations. | ||
m.rtree = rtreego.NewTree(3, minChildren, maxChildren, bulkLoad...) | ||
|
||
return m | ||
} | ||
|
||
//----------------------------------------------------------------------------- | ||
|
||
func stlPointToTriangleDistSq(p sdf.V3, triangle *render.Triangle3) (float64, bool /* falls outside? */) { | ||
// Compute the closest point | ||
closest, fallsOutside := stlClosestTrianglePointTo(p, triangle) | ||
// Compute distance to the closest point | ||
closestToP := p.Sub(closest) | ||
distance := closestToP.Length2() | ||
// Solve influence (distance) ties, by prioritizing triangles with normals more aligned to `closestToP`. | ||
// This should fix ghost triangle extensions and smooth the field over sharp angles. | ||
if fallsOutside { // <-- This is an optimization, as others have 0 extra influence in this step | ||
distance *= 1 + (1-math.Abs(closestToP.Normalize().Dot(triangle.Normal())))*stlEpsilon | ||
} | ||
//log.Println(distance, closestToP.Normalize().Dot(triangle.Normal())) | ||
return distance, fallsOutside | ||
} | ||
|
||
// https://stackoverflow.com/a/47505833 | ||
func stlClosestTrianglePointTo(p sdf.V3, triangle *render.Triangle3) (sdf.V3, bool /* falls outside? */) { | ||
edgeAbDelta := triangle.V[1].Sub(triangle.V[0]) | ||
edgeCaDelta := triangle.V[0].Sub(triangle.V[2]) | ||
edgeBcDelta := triangle.V[2].Sub(triangle.V[1]) | ||
|
||
// The closest point may be a vertex | ||
uab := stlEdgeProject(triangle.V[0], edgeAbDelta, p) | ||
uca := stlEdgeProject(triangle.V[2], edgeCaDelta, p) | ||
if uca > 1 && uab < 0 { | ||
return triangle.V[0], true | ||
} | ||
ubc := stlEdgeProject(triangle.V[1], edgeBcDelta, p) | ||
if uab > 1 && ubc < 0 { | ||
return triangle.V[1], true | ||
} | ||
if ubc > 1 && uca < 0 { | ||
return triangle.V[2], true | ||
} | ||
|
||
// The closest point may be on an edge | ||
triNormal := triangle.Normal() | ||
planeAbNormal := triNormal.Cross(edgeAbDelta) | ||
planeBcNormal := triNormal.Cross(edgeBcDelta) | ||
planeCaNormal := triNormal.Cross(edgeCaDelta) | ||
if uab >= 0 && uab <= 1 && !stlPlaneIsAbove(triangle.V[0], planeAbNormal, p) { | ||
return stlEdgePointAt(triangle.V[0], edgeAbDelta, uab), true | ||
} | ||
if ubc >= 0 && ubc <= 1 && !stlPlaneIsAbove(triangle.V[1], planeBcNormal, p) { | ||
return stlEdgePointAt(triangle.V[1], edgeBcDelta, ubc), true | ||
} | ||
if uca >= 0 && uca <= 1 && !stlPlaneIsAbove(triangle.V[2], planeCaNormal, p) { | ||
return stlEdgePointAt(triangle.V[2], edgeCaDelta, uca), true | ||
} | ||
|
||
// The closest point is in the triangle so project to the plane to find it | ||
return stlPlaneProject(triangle.V[0], triNormal, p), false | ||
} | ||
|
||
func stlEdgeProject(edge1, edgeDelta, p sdf.V3) float64 { | ||
return p.Sub(edge1).Dot(edgeDelta) / edgeDelta.Length2() | ||
} | ||
|
||
func stlEdgePointAt(edge1, edgeDelta sdf.V3, t float64) sdf.V3 { | ||
return edge1.Add(edgeDelta.MulScalar(t)) | ||
} | ||
|
||
func stlPlaneIsAbove(anyPoint, normal, testPoint sdf.V3) bool { | ||
return normal.Dot(testPoint.Sub(anyPoint)) > 0 | ||
} | ||
|
||
func stlPlaneProject(anyPoint, normal, testPoint sdf.V3) sdf.V3 { | ||
v := testPoint.Sub(anyPoint) | ||
d := normal.Dot(v) | ||
p := testPoint.Sub(normal.MulScalar(d)) | ||
return p | ||
} | ||
|
||
//----------------------------------------------------------------------------- | ||
|
||
type stlTriangle struct { | ||
*render.Triangle3 | ||
} | ||
|
||
func (s *stlTriangle) Bounds() *rtreego.Rect { | ||
bounds := sdf.Box3{Min: s.V[0], Max: s.V[0]} | ||
bounds = bounds.Include(s.V[1]) | ||
bounds = bounds.Include(s.V[2]) | ||
points, err := rtreego.NewRectFromPoints(stlToPoint(bounds.Min), stlToPoint(bounds.Max)) | ||
if err != nil { | ||
panic(err) // Implementation error | ||
} | ||
return points | ||
} | ||
|
||
func stlToPoint(v3 sdf.V3) rtreego.Point { | ||
return rtreego.Point{v3.X, v3.Y, v3.Z} | ||
} | ||
|
||
//----------------------------------------------------------------------------- | ||
|
||
// ImportSTL converts an STL model into a SDF3 surface. See ImportTriMesh. | ||
func ImportSTL(reader io.ReadSeeker, numNeighbors, minChildren, maxChildren int) (sdf.SDF3, error) { | ||
mesh, err := stl.ReadAll(reader) | ||
if err != nil { | ||
return nil, err | ||
} | ||
tris := make(chan *render.Triangle3, 128) // Buffer some triangles and send in batches if scheduler prefers it | ||
go func() { | ||
for _, triangle := range mesh.Triangles { | ||
tri := &render.Triangle3{} | ||
for i, vertex := range triangle.Vertices { | ||
tri.V[i] = sdf.V3{X: float64(vertex[0]), Y: float64(vertex[1]), Z: float64(vertex[2])} | ||
} | ||
tris <- tri | ||
} | ||
close(tris) | ||
}() | ||
return ImportTriMesh(tris, numNeighbors, minChildren, maxChildren), nil | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.