#!/usr/bin/env python3
# -*- coding: utf-8 -*-

#**************************************************************************
#   Copyright (C) 2019, Paul Lutus                                        *
#                                                                         *
#   This program is free software; you can redistribute it and/or modify  *
#   it under the terms of the GNU General Public License as published by  *
#   the Free Software Foundation; either version 2 of the License, or     *
#   (at your option) any later version.                                   *
#                                                                         *
#   This program is distributed in the hope that it will be useful,       *
#   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
#   GNU General Public License for more details.                          *
#                                                                         *
#   You should have received a copy of the GNU General Public License     *
#   along with this program; if not, write to the                         *
#   Free Software Foundation, Inc.,                                       *
#   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
#**************************************************************************

class Types:
  TETRAHEDRON,\
  OCTAHEDRON,\
  ICOSAHEDRON,\
  DODECAHEDRON,\
  CUBE,\
  HYPERCUBE,\
  BUCKYBALL = range(7)
  
# BEGIN user enties

# if DEBUG = True, print state and results

DEBUG = False

# ident must be set to one of the following

#ident = Types.TETRAHEDRON
#ident = Types.OCTAHEDRON
#ident = Types.ICOSAHEDRON
#ident = Types.DODECAHEDRON
#ident = Types.CUBE
#ident = Types.HYPERCUBE
ident = Types.BUCKYBALL

# if 3D printing, rotate objects so they rest on a face
 
printing = True

# rsc = radius in mm of edge elements (cylinders and spheres)
# i.e. 1/2 width of elements in mm

# Recommend rsc >= 3.5 for 3D printing,
# may be smaller in other cases

rsc = 3

# polydim = approximate enclosing volume edge size in mm

polydim = 100

# use savepath to provide a directory path (string) at which
# to auto-save an STL file from Blender containing
# the chosen Platonic solid, otherwise set to False

savepath = False

# END user entries

import os, sys, re
import math
from itertools import permutations

# rendering modes
    
FREECAD,BLENDER,MATPLOTLIB  = range(3)
    
RENDER_MODE = -1
    
try:
  import FreeCAD
  from FreeCAD import Base
  import PartDesignGui
  RENDER_MODE = FREECAD
  if DEBUG:
    print("Mode: FreeCAD.")
except:
  try:
    import bpy
    from mathutils import Vector
    RENDER_MODE = BLENDER
    if DEBUG:
      print("Mode: Blender.")
  except:
    try:
      from mpl_toolkits.mplot3d import Axes3D
      import matplotlib.pyplot as plt
      RENDER_MODE = MATPLOTLIB
      if DEBUG:
        print("Mode: Matplotlib.")
    except:
      print("No required libraries, quitting.")
      quit()

# vertex management class

class Vertex:
  def __init__(self,ix,iy,iz):
    self.x = ix
    self.y = iy
    self.z = iz
  
  # magnitude squared     
  def magsq(self,v):
    return (self.x-v.x)**2 + (self.y-v.y)**2 + (self.z-v.z)**2
  
  # magnitude
  def mag(self,v):
    return math.sqrt(self.magsq(v))

  def equal(self,v):
    return v.x == self.x and v.y == self.y and v.z == self.z
    
  def vec(self):
    return (self.x,self.y,self.z)
       
  # midpoint of two vectors for cylinder placement  
  def mid(self,v):
    return ((self.x+v.x)/2,(self.y+v.y)/2,(self.z+v.z)/2)
        
  def sub(self,v):
    return (v.x-self.x,v.y-self.y,v.z-self.z)
        
  def __repr__(self):
    return "(%f,%f,%f)" % (self.x,self.y,self.z)

# edge class

class Edge:
  def __init__(self,ia,ib):
    self.a = ia
    self.b = ib
    
  def r(self):
    return math.sqrt((self.a.x-self.b.x)**2 + (self.a.y - self.b.y)**2 + (self.a.z - self.b.z)**2)

  def defines(self,ia,ib):
    return (self.a.equal(ia) and self.b.equal(ib)) or (self.a.equal(ib) and self.b.equal(ia))
    
  def __repr__(self):
    return "(%s,%s)" % (self.a,self.b)

# main generator class
    
class Platonic:
  
  def main(self,ident = Types.BUCKYBALL,rsc = 3,polydim = 100,printing = True,savepath = False):
    
    # Cartesian coordinates for icosahedron vertices
    # from http://eusebeia.dyndns.org/4d/icosahedron
    # (0, ±1, ±φ)
    # (±1, ±φ, 0)
    # (±φ, 0, ±1)
    
    # Cartesian coordinates for dodecahedron vertices
    # from https://math.wikia.org/wiki/Dodecahedron
    # (±1, ±1, ±1)
    # (0, ±ϕ, ±1/ϕ)
    # (±1/ϕ, 0, ±ϕ)
    # (±ϕ, ±1/ϕ, 0)
    
    # Cartesian coordinates for Buckyball (truncated icosahedron) vertices
    # from https://en.wikipedia.org/wiki/Truncated_icosahedron
    # this object's construction requires a special algorithm
    # (0, ±1, ±3φ)
    # (±1, ±(2 + φ), ±2φ)
    # (±φ, ±2, ±(2φ + 1))
    
    # ϕ = https://en.wikipedia.org/wiki/Golden_ratio
           
    ϕ = (math.sqrt(5)+1)/2
    
    # inverse golden ratio
    # note that ϕ has a unique
    # property: ϕ - 1/ϕ = 1
    
    iϕ = 1 / ϕ
    
    # now we can get started:
    
    # 1. generate vertices
    
    # use a dict to eliminate duplicate vertices
    
    td = {}
    
    polyscaled = polydim/(2*ϕ)
    
    # tuple for most generator loops
    polyv = (-polyscaled,polyscaled)
    # special tuple for Buckyball
    polyb = (-polyscaled/3,polyscaled/3)
    
    # Buckyball requires special algorithm
    # see https://en.wikipedia.org/wiki/Truncated_icosahedron
    
    if ident == Types.BUCKYBALL:  
      cc = (
        (0,1,3*ϕ),
        (1,2+ϕ,2*ϕ),
        (ϕ,2,2*ϕ+1)
      )
      # each row in cc
      for row in cc:
          # each column in row
          for n,cols in enumerate(permutations(row)):
            # even permutations only
            if n in (0,3,4):
              x,y,z = cols
              # generate all forms
              for a in polyb:
                for b in polyb:
                  for c in polyb:
                    td[(a*x,b*y,c*z)] = 1
                    
    elif ident == Types.TETRAHEDRON:
      a = polydim / 2
      td[(a,  a,  a)] = 1
      td[(a, -a, -a)] = 1
      td[(-a, a, -a)] = 1
      td[(-a, -a, a)] = 1
                          
    else:
      for n in (range(0,1),range(0,2))[ident == Types.HYPERCUBE]:
        # used only by cube/hypercube
        w = (ϕ,ϕ/2.0)[n]
        # generate all forms
        for a in polyv:
          for b in polyv:
            for c in polyv:
              if ident == Types.OCTAHEDRON:
                td[(a * ϕ, 0, 0)] = 1
                td[(0, b * ϕ, 0)] = 1
                td[(0, 0, c * ϕ)] = 1
              elif ident == Types.DODECAHEDRON:
                td[(a,b,c)] = 1 # base cube
                td[(0,b * ϕ,c * iϕ)] = 1
                td[(a * iϕ,0,c * ϕ)] = 1
                td[(a * ϕ,b * iϕ,0)] = 1
              elif ident == Types.ICOSAHEDRON:
                td[(0,a,b * ϕ)] = 1
                td[(a,b * ϕ,0)] = 1
                td[(b * ϕ,0,a)] = 1
              elif ident == Types.HYPERCUBE or ident == Types.CUBE:
                td[(a*w,b*w,c*w)] = 1
                td[(b*w,c*w,a*w)] = 1
                td[(c*w,a*w,b*w)] = 1
              else:
                print("Option setting error!")
                quit()
    
    vertices = [Vertex(*v) for v in td.keys()]

    # 2. rotate vertices for printing
          
    # when needed, a matrix rotation to allow the object to rest
    # flat on a face for stable 3D printing
    
    def rot(x,y,ra):
      return (
        x * math.cos(ra) + y * math.sin(ra),
        x * -math.sin(ra) + y * math.cos(ra)
      )
    
    if printing:
    
      # for 3D printing, rotate so one face
      # is flush with printing surface
      # columns are in object ID order  
      
      axy = (  0,  0,    0,    0,0,0, 0)[ident]
      axz = (-45,-45,-20.9,    0,0,0,21)[ident]
      ayz = (-35,-35,    0,-31.7,0,0, 0)[ident]
      
      # only perform rotation if needed
      
      if axy != 0 or axz != 0 or ayz != 0:
      
        rxy = axy * math.pi/180
        rxz = axz * math.pi/180
        ryz = ayz * math.pi/180
      
        tvert = []
      
        for v in vertices:
          x,y,z = v.vec()
          if rxy != 0: x,y = rot(x,y,rxy)
          if rxz != 0: x,z = rot(x,z,rxz)
          if ryz != 0: y,z = rot(y,z,ryz)
          tvert += [[x,y,z]]
      
        vertices = [Vertex(*v) for v in tvert]
      
    if DEBUG:
      print("vertices: %d" % len(vertices))
      
    # 3. set up edge search parameters
    
    if ident == Types.HYPERCUBE:
      # hypercube requires more test edge lengths
      psch = polyscaled * ϕ
      pd2 = psch/2.0
      tests = (
        4.0*pd2**2, # small inner cube
        4.0*psch**2, # large outer cube
        3.0*(psch-pd2)**2 # diagonal between inner & outer
        )
    else:
      # dodecahdron, icosahedron, Buckyball
      # have edges all same length
      v1 = vertices[0]
      if ident == Types.BUCKYBALL:
        v2 = vertices[2]
      else:
        v2 = vertices[1]
      tests = (v1.magsq(v2),)
           
    # 4. define edges using vertex list
    
    delta = 0.01 # acceptance bound
    
    edges = []
    
    for a in vertices:
      for b in vertices:
        r = a.magsq(b)
        for t in tests:
          # test for a valid edge
          if abs(r - t) < delta:
            # found a candidate edge
            # see if this edge is already present
            present = False
            for edge in edges:
              present |= edge.defines(a,b)
              if present: break
            if not present:
              edges.append(Edge(a,b))
            break
    
    if DEBUG:
      print("edges: %d" % len(edges))
              
    titlestr = (
      'tetrahedron',
      'octahedron',
      'icosahedron',
      'dodecahedron',
      'cube',
      'hypercube',
      'buckyball'
    )[ident]
    
    # 5. generate renderable object
    
    if RENDER_MODE == FREECAD:
      
      Gui.activateWorkbench("PartWorkbench")
    
      App.newDocument(titlestr)
      App.setActiveDocument(titlestr)
      doc = App.getDocument(titlestr)
      
      for n,vert in enumerate(vertices):
        sph = Part.makeSphere(rsc,Base.Vector(*vert.vec()))
        Part.show(sph)
         
      for n,edge in enumerate(edges):
        veca = Base.Vector(*edge.a.vec())
        vecb = Base.Vector(*edge.b.vec())
        
        rr = edge.a.mag(edge.b)
        
        cyl = Part.makeCylinder(rsc,rr,veca,vecb.sub(veca))
        Part.show(cyl)
        
      all_objects = FreeCAD.ActiveDocument.Objects
      
      if DEBUG:
        print("list length: %d" % len(all_objects))
     
      # FreeCAD can't complete this operation for lack of memory
      # on larger objects with lots of vertices and edges
      #doc.addObject("Part::MultiFuse","Fusion")  
      #doc.Fusion.Shapes = all_objects
      #for n,obj in enumerate(all_objects):
      #  obj.ViewObject.Visibility = False
      #doc.recompute()
      #doc.Fusion.Placement=App.Placement(App.Vector(0,0,0),
      #App.Rotation(App.Vector(1,0,0),0), App.Vector(0,0,0))
      Gui.SendMsgToActiveView("ViewFit")
    
    elif RENDER_MODE == BLENDER:
      
      # user-settable default values
      segments = 32
      rings = 16
       
      def make_sphere(loc_vec,radius=1,us = segments,vs = rings):
        bpy.ops.mesh.primitive_uv_sphere_add(
          segments=us,
          ring_count=vs,
          radius=radius,
          location=loc_vec
        )
      
      def make_cylinder(loc_vec,rot_vec,radius=1,length=1,us = segments):
        bpy.ops.mesh.primitive_cylinder_add(
          radius=radius,
          depth=length,
          vertices=us,
          location=loc_vec,
          rotation=rot_vec
        ) 
      
      bpy.ops.object.select_by_type(type='MESH')
    
      bpy.ops.object.delete(use_global=False, confirm=False)
      
      for vert in vertices:
        make_sphere(vert.vec(),rsc)
        
      for edge in edges:
        # calculate cylinder length and orientation
        v1 = Vector((0,0,1))
        dd = edge.b.sub(edge.a)
        rr = edge.a.mag(edge.b)
        v2 = Vector(dd)
        rot =  v1.rotation_difference(v2).to_euler()
        make_cylinder(edge.a.mid(edge.b),rot,rsc,rr)
      
      bpy.ops.object.select_by_type(type='MESH')
      
      bpy.ops.object.join()
      
      # adjust position of new combined object
      bpy.ops.object.origin_set(type='ORIGIN_CENTER_OF_MASS', center='MEDIAN')
      bpy.context.object.location = (0,0,0)
      #bpy.context.object.rotation_euler = (0,0,0)  
      bpy.context.object.name = titlestr
      
      if savepath:
        path = savepath + "/" + re.sub(' +','_',titlestr) + ".stl"
        bpy.ops.export_mesh.stl(filepath=path,check_existing=False,use_selection=True)
    
    
    # fall-back option: draw graph in matplotlib
    elif RENDER_MODE == MATPLOTLIB:
      fig = plt.figure()
      ax = fig.add_subplot(111, projection='3d')
      
      for edge in edges:
        ax.plot((edge.a.x,edge.b.x),(edge.a.y,edge.b.y),zs=(edge.a.z,edge.b.z))
       
      for v in vertices:
        ax.scatter(*v.vec(), marker='o')
        
      ax.set_xlabel('X axis')
      ax.set_ylabel('Y axis')
      ax.set_zlabel('Z axis')
    
      plt.show()
    else:
      print("Render option Setting error, quitting.")
      return

if __name__ == "__main__":
  # execute only if run as a script
  Platonic().main(ident,rsc,polydim,printing,savepath)