#! /usr/bin/python3
# Last edited on 2024-07-26 08:21:53 by stolfi

# Generic definitions for creating solid polyhedral object models
# for the topological slicing project.

from math import sin, cos, tan, atan2, asin, hypot, sqrt, pi, inf, floor, sqrt
import rn, rmxn
import sys
import re
import random


def reven(C):
  # Rounds {C} to even multiple of {eps}.
  return eps*2*floor(C/eps/2 + 0.5)

def radians(deg):
  # Converts degrees to radians.
  return pi*deg/180

def checkpow(n, msg):
  # Checks that {n} is a power of 2:
  pp = n;
  while pp % 2 == 0: pp = pp//2
  assert pp == 1, msg


def tilt_matrix(xrot,yrot):
  
  if xrot == 0 and yrot == 0:
    rmat = rmxn.ident_matrix(3,3)
  else:
    t, et = rn.dir((tan(xrot*pi/180), tan(yrot*pi/180), 1))
    assert et != 0
    s, es = rn.dir(rn.cross3d((0,0,1), t))
    assert es != 0
    r, er = rn.dir(rn.cross3d(t, s))
    assert abs(er - 1.0) <= 0.001
    rmat = (r, s, t)
    
  return rmat
  
def face_normal_and_barycenter(Fp):
  # Compute normal {nrm} and barycenter {bar}
  # given a list of vertices {Fp[0..deg-1]}
  deg = len(Fp)
  assert deg >= 3
  po = Fp[0][:3]
  pa = Fp[1][:3]
  pb = Fp[2][:3]
  u = rn.sub(pa,po)
  v = rn.sub(pb,po)
  nrm, sz = rn.dir(rn.cross3d(u, v))

  bar = (0, 0, 0)
  for jv in range(deg):
    pj = Fp[jv][:3]
    bar = rn.mix(1.0, bar, 1.0/deg, pj)
    
  # Check planarity:
  for jv in range(deg):
    pj = Fp[jv][:3]
    u = rn.sub(pj, bar)
    s = rn.dot(u, nrm)
    assert abs(s) < 0.5*eps, f"face is not planar s = {s}"
  return nrm, bar


  ??
  ??vref = ref_vertices(S, pock)
  ??vtot, Nv = all_vertices(vref)
  ??perturb_and_round_vertices(vtot, rmat,pert)

  
  for key,v in vtot.items():
    for iv in range(len(v)):
      vi = v[iv]
      ci = rmxn.map_row(vi[:3], rmat)
      v[iv] = tuple([ pertround(ci[j],pert) for j in range(3) ] + [ vi[3], ])
  return None

def pertround(c,pert):
  # Adds a random amount in {[-pert _ pert]*2*eps} to {c}, then 
  # rounds to an even multiple of {eps}.
  
  if pert != 0: 
    d = random.uniform(-pert*2*eps, pert*2*eps)
    c = 2*eps * floor(c/(2*eps) + 0.5)
  return c

def write(wro,obj):
  write_vertices_OBJ(wro, obj)
  write_faces(wro, obj)

    # Exclude real and ghost edges from the {D} list:
    EGs = set((E[1:] + G[1:]))
    Ds = set(D[1:])
    Ds = Ds - EGs
    D = [None] + list(Ds)
    Nd = len(D) - 1 
  return
  
def perturb_and_round_vertices(vtot, rmat,pert):
  # Multiplies eack coord vector of each vertex of {vtot} (as row) by the 3x3 matrix {rmat}.
  # Then adds a random amount in {[-pert _ pert]} to each coordinate.
  # Then rounds each coordinate to an even multiple of {eps}.
  
  for key,v in vtot.items():
    for iv in range(len(v)):
      vi = v[iv]
      ci = rmxn.map_row(vi[:3], rmat)
      v[iv] = tuple([ pertround(ci[j],pert) for j in range(3) ] + [ vi[3], ])
  return None

def pertround(c,pert):
  # Adds a random amount in {[-pert _ pert]*2*eps} to {c}, then 
  # rounds to an even multiple of {eps}.
  
  if pert != 0: 
    d = random.uniform(-pert*2*eps, pert*2*eps)
    c = 2*eps * floor(c/(2*eps) + 0.5)
  return c
  
def write_vertices(wro, vtot):
  # Writes the vertices of the object to the file {wro} in OBJ format.
  # 
  # Expect {vtot} to be a dict with keys 'mmm', 'mmo', 'mmp', 'mom', ..., 'ppm',
  # 'ppo,' 'ppp' with the vertices in each 'octant' of {\RR^3},
  # Writes the vertices to {wro} in OBJ format.
  # 
  # Assumes that each vertex in {vtot} has three {float}
  # coordinates and a fourth integer element that is the OBJ index.
  # Assumes that each coordinate is rounded to an even multiple of
  # {eps}.
  
  if debug_vertices:
    sys.stderr.write("writing vertices...\n")
    
  debug_vertices = True
  
  Nv = 0
  
  if debug_vertices:
    sys.stderr.write("writing vertices...\n")

  def Prtv(p, lab):
    # Writes {p[0..2]} to {wro} as the coordinates of a new vertex.
    # Expects that the fourth element of {p} is the OBJ index of the vertex
    # and vertices are printed in the order of that index.
    # The {lab} is the vertex label "{key}:{index}" in {vtot}.
    nonlocal Nv
    Nv += 1;
    assert Nv == p[3]
    # Writes {p} to {wro}: 
    wro.write("v");
    for i in range(3): wro.write(" %.*f" % (prec, p[i]))
    wro.write(" # %s\n" % lab);
    if debug_vertices:
      sys.stderr.write("  V[%d] = %s = (" % (p[3], lab))
      for i in range(3): sys.stderr.write(" %.*f" % (prec, p[i]))
      sys.stderr.write(" )\n")
 
  for key,v in vtot.items():
    for iv in range(len(v)):
      lab = f"{key}:{iv}"
      Prtv(v[iv], lab)

  wro.write("\n")
 
  return None
  
  ??Tbar, Tnrm, D = write_triangles_OBJ(wrt, vtot, pock)

def write_faces(wro, vtot, ftot, pock):
  # Writes the faces of the object to OBJ file {wro}, assuming that
  # the vertices are in the dict {vtot} and the faces are in the array {ftot}.
  
  # Also returns tables {Fnrm} of normals to the faces and {Fbar} of their
  # barycenters.
  #
  # Also returns tables {E,G} of edges:
  #
  #   {E} Non-ghost edges of the mesh.
  #   {G} Ghost edges of the mesh.
  #
  # Each element of these arrays is a pair {(iorg,idst)} of OBJ indices
  # of the endpoint vertices. Each array {Fnrm,Fbar,E,G} is indexed
  # {[1..N]}; element 0 is not used ({None}).
  
  # Expected face, and edge counts perpendicular or parallel to each axis, and
  # general oblique:
  
  def Sved(p_org, p_dst, ghost):
    # Saves the edge from {p_org} to {p_dst} in {E} or {G} depending on
    # {ghost}, if traversed in increasing vertex index sense. 

    # Write the face:
    wro.write("f")
    deg = len(Fp)
    for rv in range(deg):
      jv = rv if order == +1 else deg-1-rv
      p = Fp[jv]
      assert len(p) == 5
      kv = p[3]
      lab = p[4]
      if debug_face:
        sys.stderr.write("    v%04d = %s = ( %9.*f %9.*f %9.*f )\n" % (kv, lab, prec, p[0], prec, p[1], prec, p[2]));
      wro.write(" %d" % kv);
    wro.write("\n");
    
    # Compute normal, barycenter:
    nrm, bar = face_normal_and_barycenter(Fp)
    Fnrm.append(nrm)
    Fbar.append(bar)
    Fp = []; 
    return None