#! /usr/bin/python3
# Last edited on 2024-09-19 08:37:55 by stolfi
 
import sys
import rn
from math import inf

  # The macro names will begin with "slicing_".

def write_object(fname, Vlst, Elst, Flst, Np, Zp_min, Zp_max, emag, vmag, prec):
  # Writes to file {fname} a POV-Ray ".inc" file that declares
  # the global parameters {slicing_num_vertices}, {slicing_num_edges}
  # {slicing_num_faces}, as well as the macros {slicing_vertices()},
  # {slicing_edges()}, and {slicing_faces()} that define the 
  # elements of the object.  See the functions below for
  # details.
  # 
  # The parameters {Vlst, Elst, Flst} are vertex, edge, and face lists
  # as explained elsewhere.
  #
  # The parameters {Np,Zp_min,Zp_max} define the number and {Z}-range
  # of the sicing planes.
  #
  # The {emag} is a magnification factor for the radius of edges.
  #
  # The {vmag} is the radius of vertices relative to the
  # radius of edges; sp {vmag=1} makes edges join smoothly but 
  # vertices are not visible.
  
  wrp = open(fname, "w")
  write_parameter_declares(wrp, Vlst, Elst, Flst, Np, Zp_min, Zp_max, emag, vmag, prec)
  wrp.write("\n")
  write_vertices_macro(wrp, Vlst, prec)
  wrp.write("\n")
  write_edges_macro(wrp, Vlst, Elst, prec)
  wrp.write("\n")
  write_faces_macro(wrp, Vlst, Flst, prec)
  wrp.close()
  
  return None

def write_parameter_declares(wrp, Vlst, Elst, Flst, Np, Zp_min, Zp_max, emag, vmag, prec):
  # Writes to the ".inc" file {wrp} the global '#declares' for the POV-ray.
  
  Nv = len(Vlst) - 1
  Ne = len(Elst) - 1
  Nf = len(Flst) - 1
  
  wrp.write("#declare slicing_num_faces = %d;\n" % (Nf))
  wrp.write("#declare slicing_num_edges = %d;\n" % (Ne))
  wrp.write("#declare slicing_num_vertices = %d;\n" % (Nv))
  
  wrp.write("\n")
 
  # Slicing plane parameters:
  wrp.write("#declare slicing_num_planes = %d;\n" % (Np))
  wrp.write("#declare slicing_planes_Z_min = %.*f;\n" % (prec, Zp_min))
  wrp.write("#declare slicing_planes_Z_max = %.*f;\n" % (prec, Zp_max))
 
  # Vertex/edge radius parameters:
  vmin = [ +inf, +inf, +inf ]
  vmax = [ -inf, -inf, -inf ]
  for kv in range(1,Nv+1):
    vk = Vlst[kv][0:3]
    for j in range(3):
      vmin[j] = min(vmin[j], vk[j])
      vmax[j] = max(vmax[j], vk[j])
  obj_radius = rn.dist(vmin, vmax)/2
  edge_radius = emag*0.004*obj_radius
  wrp.write("#declare slicing_edge_radius = %.*f;\n" % (prec+2, edge_radius))
  vertex_radius = vmag*edge_radius
  wrp.write("#declare slicing_vertex_radius = %.*f;\n" % (prec+2, vertex_radius))
  
  return None

def write_vertices_macro(wrp, Vlst, prec):
  # Writes to {wrp} a POV-Ray macro {slicing_vertices} that defines an
  # array {V[1..Nv]} with the coordinates of the vertices {Vlst[1..Nv]}
  # of the object, as 3D POV-ray points. The labels are written as
  # comments.

  Nv = len(Vlst) - 1
  
  # Write the mactro preamble:
  wrp.write("#macro slicing_vertices()\n")
  wrp.write("  #local Nv = %d;\n" % Nv)
  wrp.write("  // Returns an array {V} of {Nv+1} elements.\n")
  wrp.write("  // Element {V[kv]} is the vertex with OBJ index {kv} in {1..Nv}.\n")
  wrp.write("  // Element {V[0]} is not used.\n")
  wrp.write("  #local V = array[Nv+1]\n")
  wrp.write("  #local V[   0] = < -1, -1, -1>; // Not used.\n")

  for kv in range(1,Nv+1):
    v = Vlst[kv]
    vx,vy,vz,kv_chk,vlab = v
    assert kv_chk == kv
    wrp.write("  #local V[%4d] =" % kv)
    wrp.write(" < %.*f, %.*f, %.*f >;" % (prec, vx, prec, vy, prec, vz))
    wrp.write(" // %s\n" % vlab)

  # Write the macro postamble:
  wrp.write("  V\n")
  wrp.write("#end\n\n")
  
  return None
  
def write_edges_macro(wrp, Vlst, Elst, prec):
  # Writes to {wrp} a POV-Ray macro {slicing_edges} that defines an array
  # {E[1..Ne]} describing the edges {Elst[1..Ne]} of the object. Each
  # element {E[ke]} is a POV-ray 3D point with the indices of the origin
  # and destination vertices as {.x} and {.y} components, and the edge
  # type as the {.z} component.
  
  eps = 0.1 ** prec
  
  Nv = len(Vlst) - 1
  Ne = len(Elst) - 1

  # Write the macro preamble:
  wrp.write("#macro slicing_edges()\n")
  wrp.write("  #local Ne = %d;\n" % Ne)
  wrp.write("  // Returns an array {E} of {Ne+1} elements.\n")
  wrp.write("  // Element {E[ke]} is the triple {<ko[ke],kd[ke],ty[ke]>}, for {ke} in {1..Ne}.\n")
  wrp.write("  // Here {ko[ke],kd[ke]} are the indices of the vertices\n")
  wrp.write("  // which are the endpoints of edge {E[ke]},\n")
  wrp.write("  // and {ty[ke]} is an edge type code:\n")
  wrp.write("  // \n")
  wrp.write("  //   0 original object edge.\n")
  wrp.write("  //   1 refinement (triangulation, monotonization, etc) edge.\n")
  wrp.write("  //   2 ghost edge.\n")
  wrp.write("  // \n")
  wrp.write("  #local E = array[Ne+1]\n")
  wrp.write("  #local E[   0] = < -1, -1, -1 >; // Not used.\n")
  
  for ke in range(1, Ne+1):
    e = Elst[ke]
    kv_org, kv_dst, etype, ke_chk, elab = e;
    assert ke_chk == ke

    wrp.write("  #local E[%4d] = < %d, %d, %d >; // %s;\n" % (ke, kv_org, kv_dst, etype, elab))
    elen = rn.norm(rn.sub(Vlst[kv_dst][0:3], Vlst[kv_org][0:3]))
    assert elen >= 5*eps, f"edge {elab}--{lab_dst} too short"
    
  # Write the macro postamble:
  wrp.write("  E\n")
  wrp.write("#end\n\n")
  
  return None
 
def write_faces_macro(wrp, Vlst, Flst, prec):
  # Writes to {wrp} a POV-Ray macro {slicing_faces} that defines an array
  # {F[1..Nf]} with the planes of the faces {Flst[1..Nf]}. Each element
  # {F[kf]} is an instance of the 'prism' POV-Ray
  # primitive. 
  
  Nv = len(Vlst) - 1
  Nf = len(Flst) - 1

  # Write the macro preamble:
  wrp.write("#macro slicing_faces(explode)\n")
  wrp.write("  #local Nf = %d;\n" % Nf)
  wrp.write("  // Returns an array {F} of {Nf+1} elements.\n")
  wrp.write("  // Element {F[kf]} is an instance of the\n")
  wrp.write("  // 'prism' POV-ray primitive, for {kf} in {1..Nf}.\n")
  wrp.write("  // The {explode} parameter is a displacement to be\n")
  wrp.write("  // applied to each face in the normal direction of its normal.\n")
  wrp.write("  #local F = array[Nf+1]\n")
  wrp.write("  #local F[   0] = sphere{ <0,0,0>, 1000} // Not used.\n")
  
  Hth = (0.1 ** prec)  # Half-thickness of face.
  
  for kf in range(1, Nf+1):
    f = Flst[kf]
    Nx,Ny,Nz,Fiv,kf_chk,flab = f
    N = (Nx,Ny,Nz)
    deg = len(Fiv)
    assert kf_chk == kf
    v = Vlst[Fiv[0]]  # A vertex of the face.
    
    # Choose origin {o} two vectors {U,V} orthogonal to the face normal:
    o = Vlst[Fiv[0]][0:3]
    R, er = rn.dir(rn.sub(Vlst[Fiv[1]][0:3], o))
    assert er > 1.0e-12
    S, es = rn.dir(rn.cross3d(R,N))
    assert es > 1.0e-12
    wrp.write("  #local F[%4d] = prism{ // %s\n" % (kf, flab))
    wrp.write("    %+.*f, %+.*f, %d,\n" % (prec+1, +Hth, prec+1, -Hth, deg+1))
    for jv in range(deg):
      vj = Vlst[Fiv[jv]]
      pj = vj[0:3]
      plab = vj[4] 
      opj = rn.sub(pj, o)
      rj = rn.dot(opj, R)
      sj = rn.dot(opj, S)
      if jv == 0: assert rj == 0 and sj == 0
      wrp.write("    < %+.*f, %+.*f >,\n" % (prec+1, rj, prec+1, sj)) 
    wrp.write("    <0,0>\n") 
    wrp.write("    texture{ slicing_tx_face }\n")
    wrp.write("    translate explode*y\n")
    wrp.write("    matrix<\n")
    wrp.write("      %+.8f, %+.8f, %+.8f,\n" % (R[0], R[1], R[2]))
    wrp.write("      %+.8f, %+.8f, %+.8f,\n" % (N[0], N[1], N[2]))
    wrp.write("      %+.8f, %+.8f, %+.8f,\n" % (S[0], S[1], S[2]))
    wrp.write("      %+.8f, %+.8f, %+.8f\n"  % (o[0], o[1], o[2]))
    wrp.write("    >\n")
    wrp.write("  }\n")
    
  # Write the macro postamble:
  wrp.write("  F\n")
  wrp.write("#end\n\n")
  
  return None