#!/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 rectangular workpiece at the specified xy
			# distance around the model (default: none)
z_step = None		# maximum Z step (default: unlimited)
flip = False		# flip around X center (default: don't)
height = None		# height of the workpiece above the Z plane (can be
			# negative). Default: use model dimensions.
align_top = None	# align the Z position of the model to the workpiece
align_bottom = 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, z):
	p = v.Point
	print p[0], " ", p[1], " ", z


# Make a vector from a point. While we're at it, also apply flipping (if
# requested).

def vec(p):
	if flip:
		return Base.Vector(p[0],
		    bb.YMax - p[1] + bb.YMin, bb.ZMax - p[2] + bb.ZMin)
	else:
		return Base.Vector(p[0], p[1], p[2])

#
# Dump the current Z level (plateau or intermediate level).
#


def dump_level(wires, z):
	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 wires:
		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, z)
			last = v
		if first is not None:
			print_vec(first, z)
			print
	print


def usage():
	print >>sys.stderr, "usage:", sys.argv[0], \
	    "[-a (top|bottom)(+|-)offset] [-f] [-h height]" + \
	    "\t[-b piece_distance] [-s max_step] 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:], "a:fh:p:s:")
for opt, arg in opts:
	if opt == "-a":
		if arg[0:3] == "top":
			align_top = float(arg[3:])
		elif arg[0:6] == "bottom":
			align_bottom = float(arg[6:])
		else:
			usage()
	elif opt == "-f":
		flip = True
	elif opt == "-h":
		height = float(arg)
	elif opt == "-p":
		margin = float(arg)
	elif opt == "-s":
		z_step = float(arg)
	else:
		assert False

if len(args) != 1:
	usage()

#
# Read the STL mesh and determine its bounding box
#

mesh = Mesh.Mesh(args[0])
bb = mesh.BoundBox

#
# 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)
		vert.addFacet(vec(facet.Points[0]), vec(facet.Points[1]),
		    vec(facet.Points[2]))

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(vert.Topology, mech_eps)

z_off = 0
if height is not None:
	if align_top is not None:
		z_off = align_top - bb.ZMax
		if height > 0:
			z_off += height
	if align_bottom is not None:
		z_off = align_bottom - bb.ZMin
		if height < 0:
			z_off += height

#
# 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.
#
# We advance by at most z_step and insert intermediate layers if needed. 
#

if height is not None and height > 0:
	last_z = height 
else:
	last_z = None
if height is not None and height < 0 and z_levels[-1] > height:
	z_levels.append(height)
	
for next_z in z_levels:
	wires = shape.slice(Base.Vector(0, 0, 1), next_z + epsilon)
	if z_step is None or last_z is None or last_z + z_step >= next_z:
		dump_level(wires, next_z + z_off)
	else:
		d = next_z - last_z
		n = (d // z_step) + 1
		for i in range(0, n):
			dump_level(wires, last_z + (i + 1) * (d / n) + z_off)
	last_z = next_z
	
#
# That's all, folks !
#