topo_operators.coffee

# Polyhédronisme
#===================================================================================================
#
# A toy for constructing and manipulating polyhedra and other meshes
#
# Includes implementation of the conway polyhedral operators derived
# from code by mathematician and mathematical sculptor
# George W. Hart http://www.georgehart.com/
#
# Copyright 2011, Anselm Levskaya
# Released under the MIT License


#===================================================================================================
# Polyhedron Flagset Construct
#
# A Flag is an associative triple of a face index and two adjacent vertex indices,
# listed in geometric clockwise order (staring into the normal)
#
# Face_i -> V_i -> V_j
#
# They are a useful abstraction for defining topological transformations of the polyhedral mesh, as
# one can refer to vertices and faces that don't yet exist or haven't been traversed yet in the
# transformation code.
#
# A flag is similar in concept to a directed halfedge in halfedge data structures.
#
class polyflag
  constructor: ->
    @flags= new Object() # flags[face][vertex] = next vertex of flag; symbolic triples
    @verts= new Object() # XYZ coordinates
    @xyzs = new Object() # [symbolic names] holds vertex index

  # Add a new vertex named "name" with coordinates "xyz".
  newV: (name, xyz) ->
    if @verts[name] is undefined
      @verts[name] = 0
      @xyzs[name] = xyz

  newFlag: (facename, v1, v2) ->
    if @flags[facename] is undefined
      @flags[facename] = {}
    @flags[facename][v1] = v2

  topoly: () ->
    poly = new polyhedron()

    ctr=0 # first number the vertices
    for i,v of @verts
      poly.xyz[ctr]=@xyzs[i] # store in array
      @verts[i] = ctr
      ctr++

    ctr=0
    for i,f of @flags
      poly.face[ctr] = [] # new face
      # grab _any_ vertex as starting point
      for j,v of f
        v0 = v
        break  # need just one
      # build face out of all the edge relations in the flag assoc array
      v = v0 # v moves around face
      poly.face[ctr].push @verts[v] #record index
      v = @flags[i][v] # goto next vertex
      faceCTR=0
      while v isnt v0 # loop until back to start
        poly.face[ctr].push @verts[v]
        v = @flags[i][v]
        faceCTR++
        if faceCTR>1000 # necessary to prevent browser hangs on badly formed flagsets!
          console.log "Bad flag spec, have a neverending face:", i, @flags[i]
          break
      ctr++

    poly.name = "unknown polyhedron"
    poly


#===================================================================================================
# Polyhedron Operators
#===================================================================================================
#          for each vertex of new polyhedron:
#              call newV(Vname, xyz) with a symbolic name and coordinates
#          for each flag of new polyhedron:
#              call newFlag(Fname, Vname1, Vname2) with a symbolic name for the new face
#              and the symbolic name for two vertices forming an oriented edge
#          ORIENTATION -must- be dealt with properly to make a manifold (correct) mesh.
#          Specifically, no edge v1->v2 can ever be crossed in the -same direction- by
#          two different faces
#
#          call topoly() to assemble flags into polyhedron structure by following the orbits
#          of the vertex mapping stored in the flagset for each new face
#
#          set name as appropriate

# Kis(N)
# ------------------------------------------------------------------------------------------
# Kis (abbreviated from triakis) transforms an N-sided face into an N-pyramid rooted at the
# same base vertices.
# only kis n-sided faces, but n==0 means kiss all.
#
kisN = (poly, n, apexdist)->
  n or= 0
  apexdist or= 0.1
  console.log "Taking kis of #{if n==0 then "" else n}-sided faces of #{poly.name}..."

  flag = new polyflag()
  for p,i in poly.xyz
    # each old vertex is a new vertex
    flag.newV "v#{i}", p

  normals = poly.normals()
  centers = poly.centers()
  foundAny = false                 # alert if don't find any
  for f,i in poly.face
    v1 = "v"+f[f.length-1]
    for v in f
      v2 = "v"+v
      if f.length is n or n is 0
        foundAny = true
        apex = "apex#{i}"
        fname = "#{i}#{v1}"
        flag.newV apex, add(centers[i],mult(apexdist,normals[i])) # new vertices in centers of n-sided face
        flag.newFlag fname,   v1,   v2 # the old edge of original face
        flag.newFlag fname,   v2, apex # up to apex of pyramid
        flag.newFlag fname, apex,   v1 # and back down again
      else
        flag.newFlag "#{i}", v1, v2  # same old flag, if non-n
      v1=v2  # current becomes previous

  if not foundAny
    console.log "No #{n}-fold components were found."

  newpoly = flag.topoly()
  newpoly.name = "k" + (if n is 0 then "" else n) + poly.name
  #newpoly.xyz = adjustXYZ(newpoly, 3)
  #newpoly.xyz = canonicalXYZ(newpoly, 3)  # this tends to make results look like shit
  newpoly


# Ambo
# ------------------------------------------------------------------------------------------
# The best way to think of the ambo operator is as a topological "tween" between a polyhedron
# and its dual polyhedron.  Thus the ambo of a dual polyhedron is the same as the ambo of the
# original. Also called "Rectify".
#
ambo = (poly)->
  console.log "Taking ambo of #{poly.name}..."

  # helper func to insure unique names of midpoints
  midName = (v1, v2) -> if v1<v2 then v1+"_"+v2 else v2+"_"+v1

  flag = new polyflag()

  # For each face f in the original poly
  for f,i in poly.face
    [v1, v2] = f[-2..-1]
    for v3 in f
      if v1 < v2 # vertices are the midpoints of all edges of original poly
        flag.newV(midName(v1,v2), midpoint(poly.xyz[v1], poly.xyz[v2]))
      # two new flags:
      # One whose face corresponds to the original f:
      flag.newFlag("orig"+i,  midName(v1,v2), midName(v2,v3))
      # Another flag whose face  corresponds to (the truncated) v2:
      flag.newFlag("dual"+v2, midName(v2,v3), midName(v1,v2))
      # shift over one
      [v1, v2] = [v2, v3]

  newpoly = flag.topoly()
  newpoly.name = "a" + poly.name
  #newpoly.xyz = adjustXYZ(newpoly, 2)
  newpoly


# Gyro
# ----------------------------------------------------------------------------------------------
# This is the dual operator to "snub", i.e dual*Gyro = Snub.  It is a bit easier to implement
# this way.
#
# Snub creates at each vertex a new face, expands and twists it, and adds two new triangles to
# replace each edge.
#
gyro = (poly)->
  console.log "Taking gyro of #{poly.name}..."

  flag = new polyflag()

  for v,i in poly.xyz
    flag.newV "v"+i, unit(v)  # each old vertex is a new vertex

  centers = poly.centers() # new vertices in center of each face
  for f,i in poly.face
    flag.newV "center"+i, unit(centers[i])

  for f,i in poly.face
    [v1, v2] = f[-2..-1]
    for v,j in f
      v3 = v
      flag.newV(v1+"~"+v2, oneThird(poly.xyz[v1],poly.xyz[v2]))  # new v in face
      fname = i+"f"+v1
      flag.newFlag(fname, "center"+i,      v1+"~"+v2) # five new flags
      flag.newFlag(fname, v1+"~"+v2,  v2+"~"+v1)
      flag.newFlag(fname, v2+"~"+v1,  "v"+v2)
      flag.newFlag(fname, "v"+v2,     v2+"~"+v3)
      flag.newFlag(fname, v2+"~"+v3,  "center"+i)
      [v1, v2] = [v2, v3]                       # shift over one

  newpoly = flag.topoly()
  newpoly.name = "g" + poly.name
  #newpoly.xyz = adjustXYZ(newpoly, 3)
  newpoly


# Propellor
# ------------------------------------------------------------------------------------------
# builds a new 'skew face' by making new points along edges, 1/3rd the distance from v1->v2,
# then connecting these into a new inset face.  This breaks rotational symmetry about the
# faces, whirling them into gyres
#
propellor = (poly) ->
  console.log "Taking propellor of #{poly.name}..."

  flag = new polyflag()

  for v,i in poly.xyz
    flag.newV("v"+i, unit(v))  # each old vertex is a new vertex

  for f,i in poly.face
    [v1, v2] = f[-2..-1]
    for v in f
      v3 = "#{v}"
      flag.newV(v1+"~"+v2, oneThird(poly.xyz[v1], poly.xyz[v2]))  # new v in face, 1/3rd along edge
      fname = "#{i}f#{v2}"
      flag.newFlag("v#{i}", v1+"~"+v2,  v2+"~"+v3) # five new flags
      flag.newFlag(fname,   v1+"~"+v2,  v2+"~"+v1)
      flag.newFlag(fname,   v2+"~"+v1,     "v"+v2)
      flag.newFlag(fname,      "v"+v2,  v2+"~"+v3)
      flag.newFlag(fname,   v2+"~"+v3,  v1+"~"+v2)
      [v1, v2] = [v2, v3]                       # shift over one

  newpoly = flag.topoly()
  newpoly.name = "p" + poly.name
  #newpoly.xyz  = adjustXYZ(newpoly, 3)
  newpoly


# Reflection
# ------------------------------------------------------------------------------------------
# geometric reflection through origin
reflect = (poly) ->
  console.log "Taking reflection of #{poly.name}..."
  for i in [0..poly.xyz.length-1]
       poly.xyz[i] = mult(-1, poly.xyz[i])         # reflect each point through origin
  for i in [0..poly.face.length-1]
       poly.face[i] = poly.face[i].reverse()       # repair clockwise-ness of faces!
  poly.name = "r" + poly.name
  poly


# Dual
# ------------------------------------------------------------------------------------------------
# The dual of a polyhedron is another mesh wherein:
# - every face in the original becomes a vertex in the dual
# - every vertex in the original becomes a face in the dual
#
# So N_faces, N_vertices = N_dualfaces, N_dualvertices
#
# The new vertex coordinates are convenient to set to the original face centroids.
#
dual = (poly) ->
  console.log "Taking dual of #{poly.name}..."

  flag = new polyflag()

  face = [] # make table of face as fn of edge
  for i in [0..poly.xyz.length-1]
    face[i] = {} # create empty associative table

  for f,i in poly.face
    v1 = f[f.length-1] #previous vertex
    for v2 in f
      # THIS ASSUMES that no 2 faces that share an edge share it in the same orientation!
      # which of course never happens for proper manifold meshes, so get your meshes right.
      face[v1]["v#{v2}"] = "#{i}"
      v1=v2 # current becomes previous

  centers = poly.centers()
  for i in [0..poly.face.length-1]
    flag.newV("#{i}",centers[i])

  for f,i in poly.face
    v1 = f[f.length-1] #previous vertex
    for v2 in f
      flag.newFlag(v1, face[v2]["v#{v1}"], "#{i}")
      v1=v2 # current becomes previous

  dpoly = flag.topoly() # build topological dual from flags

  # match F index ordering to V index ordering on dual
  sortF = []
  for f in dpoly.face
    k = intersect(poly.face[f[0]],poly.face[f[1]],poly.face[f[2]])
    sortF[k] = f
  dpoly.face = sortF

  if poly.name[0] isnt "d"
    dpoly.name = "d"+poly.name
  else
    dpoly.name = poly.name[1..]

  dpoly


# Chamfer
# ----------------------------------------------------------------------------------------
# A truncation along a polyhedron's edges.
# Chamfering or edge-truncation is similar to expansion, moving faces apart and outward,
# but also maintains the original vertices. Adds a new hexagonal face in place of each
# original edge.
# A polyhedron with e edges will have a chamfered form containing 2e new vertices,
# 3e new edges, and e new hexagonal faces. -- Wikipedia
# See also http://dmccooey.com/polyhedra/Chamfer.html
#
# The dist parameter could control how deeply to chamfer.
# But I'm not sure about implementing that yet.
#
# Q: what is the dual operation of chamfering? I.e.
# if cX = dxdX, and xX = dcdX, what operation is x?

# We could "almost" do this in terms of already-implemented operations:
# cC = t4daC = t4jC, cO = t3daO, cD = t5daD, cI = t3daI
# But it doesn't work for cases like T.

chamfer = (poly, dist) ->
  console.log "Taking chamfer of #{poly.name}..."

  dist or= 0.5

  flag = new polyflag()

  normals = poly.normals()

  # For each face f in the original poly
  for f,i in poly.face
    v1 = f[f.length-1]
    v1new = i + "_" + v1

    for v2 in f
      # TODO: figure out what distances will give us a planar hex face.
      # Move each old vertex further from the origin.
      flag.newV(v2, mult(1.0 + dist, poly.xyz[v2]))
      # Add a new vertex, moved parallel to normal.
      v2new = i + "_" + v2
      flag.newV(v2new, add(poly.xyz[v2], mult(dist*1.5, normals[i])))
      # Four new flags:
      # One whose face corresponds to the original face:
      flag.newFlag("orig"+i, v1new, v2new)
      # And three for the edges of the new hexagon:
      facename = (if v1<v2 then "hex"+v1+"_"+v2 else "hex"+v2+"_"+v1)
      flag.newFlag(facename, v2, v2new)
      flag.newFlag(facename, v2new, v1new)
      flag.newFlag(facename, v1new, v1)
      v1 = v2
      v1new = v2new

  newpoly = flag.topoly()
  newpoly.name = "c" + poly.name
  newpoly


# Whirl
# ----------------------------------------------------------------------------------------------
# Gyro followed by truncation of vertices centered on original faces.
# This create 2 new hexagons for every original edge.
# (https://en.wikipedia.org/wiki/Conway_polyhedron_notation#Operations_on_polyhedra)
#
# Possible extension: take a parameter n that means only whirl n-sided faces.
# If we do that, the flags marked #* below will need to have their other sides
# filled in one way or another, depending on whether the adjacent face is
# whirled or not.

whirl = (poly, n) ->
  console.log "Taking whirl of #{poly.name}..."
  n or= 0
  
  flag = new polyflag()

  for v,i in poly.xyz
    flag.newV "v"+i, unit(v)  # each old vertex is a new vertex

  centers = poly.centers() # new vertices around center of each face
  #for f,i in poly.face
  #  # Whirl: use "center"+i+"~"+v1
  #  flag.newV "center"+i+"~"+v1, unit(centers[i])

  for f,i in poly.face
    [v1, v2] = f[-2..-1]
    for v,j in f
      v3 = v
      # New vertex along edge
      v1_2 = oneThird(poly.xyz[v1],poly.xyz[v2])
      flag.newV(v1+"~"+v2, v1_2)
      # New vertices near center of face
      cv1name = "center"+i+"~"+v1
      cv2name = "center"+i+"~"+v2
      flag.newV(cv1name, unit(oneThird(centers[i], v1_2))) 
      fname = i+"f"+v1
      # New hexagon for each original edge
      flag.newFlag(fname, cv1name,      v1+"~"+v2)
      flag.newFlag(fname, v1+"~"+v2,  v2+"~"+v1) #*
      flag.newFlag(fname, v2+"~"+v1,  "v"+v2)    #*
      flag.newFlag(fname, "v"+v2,     v2+"~"+v3) #*
      flag.newFlag(fname, v2+"~"+v3,  cv2name)
      flag.newFlag(fname, cv2name, cv1name)
      # New face in center of each old face      
      flag.newFlag("c"+i, cv1name, cv2name)
      
      [v1, v2] = [v2, v3]                       # shift over one

  newpoly = flag.topoly()
  newpoly.name = "w" + poly.name
  #newpoly.xyz = adjustXYZ(newpoly, 3)
  newpoly
  

# insetN
# ------------------------------------------------------------------------------------------
insetN = (poly, n, inset_dist, popout_dist)->
  n or= 0
  inset_dist  or= 0.5
  popout_dist or= -0.2

  console.log "Taking inset of #{if n==0 then "" else n}-sided faces of #{poly.name}..."

  flag = new polyflag()
  for p,i in poly.xyz
    # each old vertex is a new vertex
    flag.newV "v#{i}", p

  normals = poly.normals()
  centers = poly.centers()
  for f,i in poly.face #new inset vertex for every vert in face
    if f.length is n or n is 0
      for v in f
        flag.newV "f"+i+"v"+v, add(tween(poly.xyz[v],centers[i],inset_dist),mult(popout_dist,normals[i]))

  foundAny = false                 # alert if don't find any
  for f,i in poly.face
    v1 = "v"+f[f.length-1]
    for v in f
      v2 = "v"+v
      if f.length is n or n is 0
        foundAny = true
        fname = i + v1
        flag.newFlag fname,      v1,       v2
        flag.newFlag fname,      v2,       "f"+i+v2
        flag.newFlag fname, "f"+i+v2,  "f"+i+v1
        flag.newFlag fname, "f"+i+v1,  v1
        #new inset, extruded face
        flag.newFlag "ex"+i, "f"+i+v1,  "f"+i+v2
      else
        flag.newFlag i, v1, v2  # same old flag, if non-n
      v1=v2  # current becomes previous

  if not foundAny
    console.log "No #{n}-fold components were found."

  newpoly = flag.topoly()
  newpoly.name = "n" + (if n is 0 then "" else n) + poly.name
  #newpoly.xyz = adjustXYZ(newpoly, 3)
  #newpoly.xyz = canonicalXYZ(newpoly, 3)  # this tends to make results look like shit
  newpoly

# ExtrudeN
# ------------------------------------------------------------------------------------------
extrudeN = (poly, n)->
  n or= 0
  console.log "Taking extrusion of #{if n==0 then "" else n}-sided faces of #{poly.name}..."

  flag = new polyflag()
  for p,i in poly.xyz
    # each old vertex is a new vertex
    flag.newV "v#{i}", p

  normals = poly.normals()
  centers = poly.centers()
  for f,i in poly.face #new inset vertex for every vert in face
    if f.length is n or n is 0
      for v in f
        #flag.newV "f"+i+"v"+v, add(midpoint(poly.xyz[v],centers[i]),mult(-0.2,normals[i]))
        flag.newV "f"+i+"v"+v, add(poly.xyz[v], mult(0.3,normals[i]))

  foundAny = false                 # alert if don't find any
  for f,i in poly.face
    v1 = "v"+f[f.length-1]
    for v in f
      v2 = "v"+v
      if f.length is n or n is 0
        foundAny = true
        #fname = i+v1
        flag.newFlag i+v1,       v1,       v2
        flag.newFlag i+v1,       v2, "f"+i+v2
        flag.newFlag i+v1, "f"+i+v2, "f"+i+v1
        flag.newFlag i+v1, "f"+i+v1,       v1
        #new inset, extruded face
        flag.newFlag "ex"+i, "f"+i+v1,  "f"+i+v2
      else
        flag.newFlag i, v1, v2  # same old flag, if non-n
      v1=v2  # current becomes previous

  if not foundAny
    console.log "No #{n}-fold components were found."

  newpoly = flag.topoly()
  newpoly.name = "x" + (if n is 0 then "" else n) + poly.name
  #console.log newpoly
  #newpoly.xyz = adjustXYZ(newpoly, 3)
  #newpoly.xyz = canonicalXYZ(newpoly, 3)  # this tends to make results look like shit
  newpoly

# hollow / skeletonize
# ------------------------------------------------------------------------------------------
hollow = (poly, n, inset_dist, thickness)->
  n or= 0
  inset_dist  or= 0.5
  thickness or= 0.2

  console.log "Skeletonizing #{if n==0 then "" else n}-sided faces of #{poly.name}..."

  dualnormals = dual(poly).normals()
  normals = poly.normals()
  centers = poly.centers()

  flag = new polyflag()
  for p,i in poly.xyz
    # each old vertex is a new vertex
    flag.newV "v#{i}", p
    flag.newV "downv#{i}",  add(p,mult(-1*thickness,dualnormals[i]))

  for f,i in poly.face #new inset vertex for every vert in face
    #if f.length is n or n is 0
    for v in f
      flag.newV "fin"+i+"v"+v, tween(poly.xyz[v],centers[i],inset_dist)
      flag.newV "findown"+i+"v"+v, add(tween(poly.xyz[v],centers[i],inset_dist),mult(-1*thickness,normals[i]))

  #foundAny = false                 # alert if don't find any
  for f,i in poly.face
    v1 = "v"+f[f.length-1]
    for v in f
      v2 = "v"+v
      #if f.length is n or n is 0
      foundAny = true
      fname = i + v1
      flag.newFlag fname,      v1,       v2
      flag.newFlag fname,      v2,       "fin"+i+v2
      flag.newFlag fname, "fin"+i+v2,  "fin"+i+v1
      flag.newFlag fname, "fin"+i+v1,  v1

      fname = "sides"+i + v1
      flag.newFlag fname, "fin"+i+v1,     "fin"+i+v2
      flag.newFlag fname, "fin"+i+v2,     "findown"+i+v2
      flag.newFlag fname, "findown"+i+v2, "findown"+i+v1
      flag.newFlag fname, "findown"+i+v1, "fin"+i+v1

      fname = "bottom"+i + v1
      flag.newFlag fname,  "down"+v2,      "down"+v1
      flag.newFlag fname,  "down"+v1,      "findown"+i+v1
      flag.newFlag fname,  "findown"+i+v1, "findown"+i+v2
      flag.newFlag fname,  "findown"+i+v2, "down"+v2

      #new inset, extruded face
      #flag.newFlag "ex"+i, "fin"+i+v1,  "fin"+i+v2
      #else
      #  flag.newFlag i, v1, v2  # same old flag, if non-n
      v1=v2  # current becomes previous

  #if not foundAny
  #  console.log "No #{n}-fold components were found."

  newpoly = flag.topoly()
  #newpoly.name = "h" + (if n is 0 then "" else n) + poly.name
  newpoly.name = "h" + poly.name
  #newpoly.xyz = adjustXYZ(newpoly, 3)
  #newpoly.xyz = canonicalXYZ(newpoly, 3)  # this tends to make results look like shit
  newpoly


# StellaN
# ------------------------------------------------------------------------------------------
stellaN = (poly)->
  console.log "Taking stella of #{poly.name}..."

  centers = poly.centers()  # calculate face centers

  flag = new polyflag()
  for p,i in poly.xyz
    flag.newV "v#{i}", p      # each old vertex is a new vertex

  # iterate over triplets of faces v1,v2,v3
  for f,i in poly.face
    v1 = "v"+f[f.length-2]
    v2 = "v"+f[f.length-1]
    vert1 = poly.xyz[f[f.length-2]]
    vert2 = poly.xyz[f[f.length-1]]
    for v in f
      v3 = "v"+v
      vert3 = poly.xyz[v]
      v12=v1+"~"+v2 # names for "oriented" midpoints
      v21=v2+"~"+v1
      v23=v2+"~"+v3

      # on each Nface, N new points inset from edge midpoints towards center = "stellated" points
      flag.newV v12, midpoint( midpoint(vert1,vert2), centers[i] )

      # inset Nface made of new, stellated points
      flag.newFlag "in#{i}",      v12,       v23

      # new tri face constituting the remainder of the stellated Nface
      flag.newFlag "f#{i}#{v2}",      v23,      v12
      flag.newFlag "f#{i}#{v2}",       v12,      v2
      flag.newFlag "f#{i}#{v2}",      v2,      v23

      # one of the two new triangles replacing old edge between v1->v2
      flag.newFlag "f"+v12,     v1,        v21
      flag.newFlag "f"+v12,     v21,       v12
      flag.newFlag "f"+v12,      v12,       v1

      [v1,v2]=[v2,v3]  # current becomes previous
      [vert1,vert2]=[vert2,vert3]

  newpoly = flag.topoly()
  newpoly.name = "l" + poly.name
  #newpoly.xyz = adjustXYZ(newpoly, 3)
  #newpoly.xyz = canonicalXYZ(newpoly, 3)  # this tends to make results look like shit
  newpoly