#! /usr/bin/python3
# Last edited on 2024-09-20 11:44:49 by stolfi

from math import floor, pi, inf
import sys
import rn

def make_null_array(dims):
  # Creates an array of {None}s with {dims[k]} elements along index {k}. 
  
  arr = [ None ]*dims[0]
  if len(dims) > 1:
    for ks in range(dims[0]):
      arr[ks] = make_null_array(dims[1:])
  return arr

def reven(C, eps):
  # 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 make_vertex(vx, vy, vz, kv, vlab):
  # Creates a vertex record {v} from given data.
  # 
  # Each vertex record a tuple {(vx,vy,vz,kv,vlab)} where {vx,vy,vz} are
  # the coordinates of the vertex, {kv} is the index (from 1) of the
  # vertex in the global vertex list, and {vlab} is a string that can be
  # used to identify the vertex in comments or printouts.
  v = ( vx, vy, vz, kv, vlab )
  return v

def prtv(v, prec):
  # Prints the point record {v}.
  vx,vy,vz,kv,vlab = v
  sys.stderr.write("  V[%04d] = (" % kv)
  sys.stderr.write(" %+12.*f" % (prec+1, vx))
  sys.stderr.write(" %+12.*f" % (prec+1, vy))
  sys.stderr.write(" %+12.*f" % (prec+1, vz))
  sys.stderr.write(" ) # %s\n" % vlab) 

def make_edge(p_org, p_dst, ke, etype, elab, prec):
  # Each edge is represented by an /edge record/ (tuple) {(kvorg, kvdst,
  # etype, ke, elab)} where {kvorg,kvdst} are the indices of the origin
  # and destination vertices, sorted so that {kvorg < kvdst}; {ke} is a
  # sequential edge index (from 1); and {elab} is a string that
  # identifies legibly the edge, useful for printouts, comments, etc.
  # 
  # The {etype} is 0 for edges of the original object, 1 for "ghost"
  # edges that connect components of the border in a single loop, and 2
  # for subdivision edges.
  #

  # If the index of vertex {pv_org} is less than that of {pv_dst},
  # returns an edge record (tuple) from {p_org} to {p_dst} with type {etype} and label 
  # {elab}.  Otherwise returns {None}.
  #
  # If {elab} is {None}, makes a label up from the two vertex labels.

  kv_org = p_org[3]
  kv_dst = p_dst[3]
  # Avoid duplicate edges (even ghost ones):
  if kv_org < kv_dst:
    if elab == None:
      elab = p_org[4] + "--" + p_dst[4]
    e = ( kv_org, kv_dst, etype, ke, elab )
    return e
  else:
    return None
 
def prte(e, prec):
  # Prints the edge record {e}.
  kv_org,kv_dst,etype,ke,elab = e
  sys.stderr.write("  E[%04d] v%d --> v%d" % (ke, kv_org, kv_dst))
  sys.stderr.write(" type = %d" % etype)
  sys.stderr.write(" # %s\n" % elab)

def min_edge_length(Vlst, Elst):
  # Find length of shortest edge: 
  dmin = +inf
  for ke in range(1,len(Elst)):
    ek = Elst[ke]
    kv_org,kv_dst,etype,ke_chk,elab = ek
    assert ke_chk == ke
    p = Vlst[kv_org][0:3]
    q = Vlst[kv_dst][0:3]
    d = rn.dist(p, q)
    dmin = min(dmin, d)
  return dmin

def check_repeated_edge(ne, Elst, Eset, prec):
  # Checks whether edges {Elst[1..ne-1]} are all
  # distinct from {Elst[ne]}. Assumes {Eset} is a
  # set of all pairs {(kv_org,kv_dst)} of vertices,
  # each sorted. Adds the edge {Elst[ne]} to that set.
  en = Elst[ne]
  kv_org = en[0]
  kv_dst = en[1]
  assert kv_org < kv_dst
  if (kv_org, kv_dst) in Eset:
    # Look for the culprit and print it:
    for je in range(1,ne):
      ej = Elst[je]
      jv_org = ej[0]
      jv_dst = ej[1]
      if kv_org == jv_org and kv_dst == jv_dst:
        sys.stderr.write( "** repeated edge:\n")
        prte(ej, prec)
        prte(en, prec)
        assert False
  else:
    Eset.add((kv_org, kv_dst))
  return None

def make_face(Fiv, kf, flab, Vlst, prec):
  # Each face is represented by a tuple {(Nx, Ny, Nz, Fv, kf, flab)} where
  # {(Nx,Ny,Nz)} is the face's outwards-pointing normal direction
  # vector, {Fiv} is a list of the indices (from 1) of its corners in
  # {Vlst}, {kf} is a sequential face index (from 1), and {flab} is a
  # string that identifies legibly the face, useful for printouts,
  # comments, etc.
  #
  # This procedure returns a face record (tuple) as above from the 
  # given parameters {Fiv,kf,flab}.  Parameter {Vlst} should be a
  # list of vertex records (tuples).
  #
  # The list {Fiv} must have at least three vertices. Computes the face
  # normal {N} from the first three vertices.

  eps = 0.1**prec   # Coordinate quantum.

  deg = len(Fiv);
  assert deg >= 3, ("face has only %d vertices" % deg)
  
  # Compute barycenter {bar}:
  Vsum = ( 0, 0, 0)
  for jv in range(0,deg):
    pj = Vlst[Fiv[jv]][0:3]
    Vsum = rn.add(Vsum, pj)
  bar = rn.scale(1.0/deg, Vsum)

  # Compute normal {Fn} and area {Area}:
  Asum = ( 0, 0, 0 )
  for jv in range(0,deg-1):
    pa = Vlst[Fiv[jv]][0:3]
    pb = Vlst[Fiv[jv+1]][0:3]
    u = rn.sub(pa,bar)
    v = rn.sub(pb,bar)
    A = rn.cross3d(u, v)
    Asum = rn.add(Asum, A)
  N, Nsz = rn.dir(Asum)
  assert Nsz > 1.0e-15
  Area = Nsz/2

  # Check planarity:
  for jv in range(3,deg):
    vj = Vlst[Fiv[jv]][0:3]
    vjlab = Vlst[Fiv[jv]][4]
    u = rn.sub(vj, bar)
    s = rn.dot(u, N)
    if abs(s) > 1.0e-7 + 3*eps:
      sys.stderr.write(f"** face is not planar s = {s:.12f} area = {Area:.12f}\n")
      rn.write(sys.stderr, "  bar = ( ", bar, "%.12f", " ", " )\n")
      rn.write(sys.stderr, "  vj  = ( ", vj,  "%.12f", " ", " )\n")
      rn.write(sys.stderr, "  N   = ( ", N,   "%.12f", " ", " )\n")
      assert False

  f = ( N[0], N[1], N[2], Fiv, kf, flab)
  return f
  
def prtf(f, Vlst, prec):
  # Prints the face record {f}.
  Nx,Ny,Nz,Fiv,kf,flab = f
  sys.stderr.write("  F[%04d] N = (" % kf)
  sys.stderr.write(" %+12.*f" % (prec+1, Nx))
  sys.stderr.write(" %+12.*f" % (prec+1, Ny))
  sys.stderr.write(" %+12.*f" % (prec+1, Nz))
  sys.stderr.write(" ) corners = (")
  for kv in Fiv: 
    pk = Vlst[kv]
    pklab = pk[4]
    sys.stderr.write(" %s" % pklab)
  sys.stderr.write(" )# %s\n" % flab)
  
def min_triangle_width(Vlst, Flst):
  # Find width of narrowest triangle: 
  hmin = +inf
  for kf in range(1,len(Flst)):
    fk = Flst[kf]
    Nx,Ny,Nz,Fiv,kf_chk,flab = fk
    assert kf_chk == kf
    assert len(Fiv) == 3, "non-triangular face"
    p = Vlst[Fiv[0]][0:3]
    q = Vlst[Fiv[1]][0:3]
    r = Vlst[Fiv[2]][0:3]
    dmax = max(rn.dist(p,q), max(rn.dist(q,r), rn.dist(r,p))) # Longest side.
    area = 0.5*rn.norm(rn.cross3d(rn.sub(q,p),rn.sub(r,p)))
    h = area/dmax # Height over longest side.
    hmin = min(hmin, h)
  return hmin