sfc/slicer.py: FreeCAD-based slicer, first commit

sfc/slicer.py

@@ -0,0 +1,161 @@
+#!/usr/bin/python
+#
+# slicer.py - FreeCAD-based STL to Gnuplot slicer
+#
+# Written 2015 by Werner Almesberger
+# Copyright 2015 by Werner Almesberger
+#
+# This program//library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License, or (at your option) any later version.
+#
+
+
+import sys
+
+sys.path.append("/usr/lib/freecad/lib")
+
+
+import FreeCAD, Part, Mesh
+import os, getopt
+from FreeCAD import Base
+from math import hypot
+
+
+epsilon = 0.0001	# acceptable math rounding error and slicing offset
+mech_eps = 0.01		# acceptable mechanical deviation
+margin = None		# draw a workpiece at the specified xy distance around
+			# the model (default: none)
+
+
+def dist(a, b):
+	pa = a.Point
+	pb = b.Point
+	return hypot(pa[0] - pb[0], pa[1] - pb[1])
+
+
+def print_vec(v):
+	p = v.Point
+	print p[0], " ", p[1], " ", p[2] - epsilon
+
+
+def usage():
+	print >>sys.stderr, "usage:", sys.argv[0], "file.stl"
+	sys.exit(1)
+
+
+#
+# FreeCAD prints progress information to stdout instead of stderr.
+# We don't want that ...
+#
+
+stdout = os.dup(1)
+os.dup2(2, 1)
+sys.stdout = os.fdopen(stdout, "w")
+
+opts, args = getopt.getopt(sys.argv[1:], "b:")
+for opt, arg in opts:
+	if opt == "-b":
+		margin = float(arg)
+	else:
+		assert False
+
+if len(args) != 1:
+	usage()
+
+#
+# Read the STL mesh
+#
+
+mesh = Mesh.Mesh(args[0])
+
+#
+# The 2.5D model consists of "plateaus" (facets parallel to the xy plane) and
+# "walls" (facets parallel to the z axis). Anything else is an error and will
+# produce incorrect results.
+#
+# We use plateau facets only for their z position, as indication where to mill
+# a plateau. Wall facets are kept for later use.
+#
+
+vert = Mesh.Mesh()
+z_raw = {}
+max_nz = 0
+inclined = 0
+for facet in mesh.Facets:
+	if abs(facet.Normal.z) >= 1 - epsilon:
+		z_raw[facet.Points[0][2]] = 1
+	else:
+		nz = abs(facet.Normal.z)
+		if nz > epsilon:
+			inclined += 1
+			max_nz = max(max_nz, nz)
+		v1 = Base.Vector(facet.Points[0])
+		v2 = Base.Vector(facet.Points[1])
+		v3 = Base.Vector(facet.Points[1])
+		vert.addFacet(v1, v2, v3)
+
+if inclined:
+	print >>sys.stderr	# FreeCAD progress reporting messes up newlines
+	print >>sys.stderr, inclined, "inclined facets, maximum normal", max_nz
+
+#
+# @@@ This is perhaps a bit too paranoid
+#
+# I wrote the Z noise filtering because I had mis-read perfectly good
+# distinct coordinates as being essentially the same value but with
+# rounding errors.
+#
+
+z_levels = []
+last = None
+for z in sorted(z_raw.keys(), reverse = True):
+	if last is None or last - z > epsilon:
+		z_levels.append(z)
+		last = z
+
+#
+# Convert the walls to a FreeCAD shape
+#
+
+shape = Part.Shape()
+shape.makeShapeFromMesh(mesh.Topology, mech_eps)
+bb = shape.BoundBox
+
+#
+# Iterate over all plateaus and determine how they intersect with the walls.
+# For this, we add a small offset to the z position so that we intersect above
+# the plateau.
+#
+
+for z in z_levels:
+	print "# level z = ", z
+
+	if margin is not None:
+		print bb.XMin - margin, " ", bb.YMin - margin, " ", z
+		print bb.XMax + margin, " ", bb.YMin - margin, " ", z
+		print bb.XMax + margin, " ", bb.YMax + margin, " ", z
+		print bb.XMin - margin, " ", bb.YMax + margin, " ", z
+		print bb.XMin - margin, " ", bb.YMin - margin, " ", z
+		print
+
+	for wire in shape.slice(Base.Vector(0, 0, 1), z + epsilon):
+		print "# wire = ", wire
+		first = None
+		last = None
+		for e in wire.Edges:
+			v = e.Vertexes[0]
+			if first is None:
+				first = v
+			if last is None or dist(v, last) >= mech_eps:
+				print_vec(v)
+			last = v
+		if first is not None:
+			print_vec(first)
+			print
+	print
+
+#
+# That's all, folks !
+#