from geometry import *
from collections import defaultdict
import sys


def sign(x):
   if x > 0: return 1
   if x < 0: return -1
   return 0


class Solid(object):
   def __init__(self, name):
      self.name = name
      self.vertices = []
      self.facets = []

   def read_from_3d_file(self):
      """Import data from my old raw '3d' file format.
      """
      f = open(self.name + '.3d')
      n_facets = int(f.readline())
      for i in range(n_facets):
         n_vertices = int(f.readline())
         assert n_vertices == 3, 'Only triangles please!'
         vtx_idxs = []
         for j in range(n_vertices):
            v = Vector3D(*map(int, f.readline().split()))
            found = False
            for (k, w) in enumerate(self.vertices):
               if v == w:
                  found = True
                  break
            if not found:
               self.vertices.append(v)
               k = len(self.vertices) - 1
            vtx_idxs.append(k)
         self.facets.append(vtx_idxs)

   def Triangle_from_facet(self, i):
      return Triangle(*map(self.vertices.__getitem__, self.facets[i]))

   def triangles(self):
      return map(self.Triangle_from_facet, range(len(self.facets)))

   def write_to_stl_file(self):
      f = open(self.name + '.stl', 'w')
      f.write('solid %s\n' % self.name)
      for T in self.triangles():
         n = T.normal()
         n.normalize()
         f.write('  facet normal %f %f %f\n' % (n.x, n.y, n.z))
         f.write('    outer loop\n')
         f.write('      vertex %f %f %f\n' % (T.v1.x, T.v1.y, T.v1.z))
         f.write('      vertex %f %f %f\n' % (T.v2.x, T.v2.y, T.v2.z))
         f.write('      vertex %f %f %f\n' % (T.v3.x, T.v3.y, T.v3.z))
         f.write('    endloop\n')
         f.write('  endfacet\n')
      f.write('endsolid\n')

   def fix_normals(self):
      """The raycasting algorithm used below to determine if a normal is outward-pointing
      runs into trouble in the event that the ray passes through the edge of a triangle.
      For the geometry I have to hand this actually happens and causes a problem. An easy
      hack for current purposes is to slightly perturb the source of the ray so that it
      does not hit the edge. A better solution would be required for more general use.
      """
      delta = Vector3D(1e-4, 1e-4, 1e-4)
      for (i, T) in enumerate(self.triangles()):
         r = Ray(T.cog() + delta, T.normal())
         n = 0
         for (j, TT) in enumerate(self.triangles()):
            if j == i:
               continue
            if TT.get_intersection(r) is not None:
               n += 1
         if n % 2:
            l = self.facets[i]
            self.facets[i] = [l[1], l[0], l[2]]

   def fix_T_junctions(self):
      for (k, v) in enumerate(self.vertices):
         for (j, f) in enumerate(list(self.facets)): # Use list to get a copy since we modify as we go.
            n_sides = len(f)
            assert n_sides == 3, 'Only triangular facets allowed.'
            for i0 in range(n_sides):
               i1, i2 = (i0 + 1) % n_sides, (i0 + 2) % n_sides
               if strictly_inside_line_segment(v, self.vertices[f[i0]], self.vertices[f[i1]]):
                  self.facets[j] = [f[i0], k, f[i2]]
                  self.facets.append([k, f[i1], f[i2]])

   def create_clearance(self, size, boundary_constant):
      for n in [Vector3D(1, 0, 0), Vector3D(0, 1, 0), Vector3D(0, 0, 1)]:
         planes = defaultdict(list)
         for (i, T) in enumerate(self.triangles()):
            if n.dot(T.v1) == n.dot(T.v2) == n.dot(T.v3):
               planes[n.dot(T.v1)].append(i)
         for (c, fs) in planes.items():
            if abs(c) == boundary_constant:
               """boundary_constant has to be supplied because it is only possible to deduce its
               value by knowing all pieces simultaneously and we're working one piece at a time.
               """
               continue
            last_sign = sign(n.dot(self.Triangle_from_facet(fs[0]).normal()))
            vs = set(self.facets[fs[0]])
            normals_aligned = True
            for f in fs[1:]:
               if last_sign != sign(n.dot(self.Triangle_from_facet(f).normal())):
                  normals_aligned = False
                  break
               vs |= set(self.facets[f])
            if normals_aligned:
               for v in vs:
                  """It is possible that moving a vertex like this could push it through the side of one of
                  the triangles it belongs to. Ideally we should really check for this and abort if so.
                  """
                  self.vertices[v] -= n * last_sign * size
            else:
               print 'Not moving plane with normal', n.x, n.y, n.z, 'constant', c, 'to preserve facet coplanarity'

   def get_bounding_box(self):
      """Note that, as in many other places, we take advantage of the fact that we know all faces are parallel
      to coordinate planes.
      """
      return ((-min(map(lambda v: v.x, self.vertices)), Vector3D(-1, 0, 0)),
              (max(map(lambda v: v.x, self.vertices)), Vector3D(1, 0, 0)),
              (-min(map(lambda v: v.y, self.vertices)), Vector3D(0, -1, 0)),
              (max(map(lambda v: v.y, self.vertices)), Vector3D(0, 1, 0)),
              (-min(map(lambda v: v.z, self.vertices)), Vector3D(0, 0, -1)),
              (max(map(lambda v: v.z, self.vertices)), Vector3D(0, 0, 1)))

   def place_on_z_plane(self, dx = 0, dy = 0):
      max_area = -1
      for (c, m) in self.get_bounding_box():
         area = 0
         for T in self.triangles():
            if m.dot(T.v1) == m.dot(T.v2) == m.dot(T.v3) == c:
               area += T.area()
         if area > max_area:
	    max_area = area
            n = -m
	    dz = c
      Z = Vector3D(0, 0, 1)
      v = Z.cross(n)
      if v.length2() != 0:
         M = Matrix3x3(v, Z, n)
      else:
         Y = Vector3D(0, 1, 0)
         M = Matrix3x3(Y.cross(n), Y, n)
      for i in range(len(self.vertices)):
         self.vertices[i] = M * self.vertices[i]
	 self.vertices[i].x += dx
	 self.vertices[i].y += dy
	 self.vertices[i].z += dz


piece = Solid(sys.argv[1][:-3] if sys.argv[1].endswith('.3d') else sys.argv[1])
piece.read_from_3d_file()
piece.fix_normals()
piece.fix_T_junctions()
piece.create_clearance(0.5, 25)
piece.place_on_z_plane(int(sys.argv[2]), int(sys.argv[3])) # Sorry, horrible hack!
piece.write_to_stl_file()