Gyromesh (hyperbolic mesh) with Julia

gyromesh.jl

import Meshes
using LinearAlgebra
using MeshViz
using GLMakie
using Makie
using Printf

# beta factor
function betaF(A, s)
  1.0 / √(1 + LinearAlgebra.norm_sqr(A) / (s*s))
end

# gyroaddition
function gyroadd(A, B, s)
  betaA = betaF(A, s)
  betaB = betaF(B, s)
  (
    1 + (betaA / (1 + betaA)) * dot(A, B) / (s*s) + (1 - betaB) / betaB
  ) * A + B
end

# scalar gyromultiplication
function gyroscalar(r, A, s)
  Anorm = norm(A)
  (s / Anorm) * sinh(r * asinh(Anorm / s)) * A
end

# point of coordinate t on the gyroline (AB)
function gyroABt(A, B, t, s)
  gyroadd(A, gyroscalar(t, gyroadd(-A, B, s), s), s)
end

# gyromidpoint
function gyromidpoint(A, B, s)
  gyroABt(A, B, 0.5, s)
end

# gyrosegment
function gyrosegment(A, B, s; n=50)
  A, B = Meshes.coordinates(A), Meshes.coordinates(B)
  [Meshes.Point(gyroABt(A, B, t, s)) for t in LinRange(0, 1, n)]
end

# edge index (I don't remember what is it...)
function edgeindex(from, to, size)
  row = min(from, to)
  col = max(from, to)
  Int(row * size - (row * (row + 1)) / 2 - (size - col))
end

# gyrosubdivision
function gyrosubdiv(s, vertices, triangles)
  nv = length(vertices)
  nt = length(triangles)
  nvmax = Int(nv + (nv * (nv + 1)) / 2)
  newvertices = fill([0.0, 0.0, 0.0], nvmax)
  for i in 1:nv
    newvertices[i] = vertices[i]
  end
  vcnt = nv
  em = fill(-1, Int(nv * (nv + 1) / 2))
  newtriangles = Vector{Tuple{Int,Int,Int}}(undef, 4 * nt)
  for i in 1:nt
    local ithis, iprev, inext
    visited = Int[]
    for j in 1:3
      iprev = triangles[i][((j + 1) % 3) + 1]
      ithis = triangles[i][j]
      inext = triangles[i][(j % 3) + 1]
      mindex = edgeindex(ithis, inext, nv)
      enext = em[mindex]
      if enext == -1
        vcnt = vcnt + 1
        enext = vcnt
        em[mindex] = enext
      end
      mindex = edgeindex(iprev, ithis, nv)
      eprev = em[mindex]
      if eprev == -1
        vcnt = vcnt + 1
        eprev = vcnt
        em[mindex] = eprev
      end
      newtriangles[(i-1)*4 + j] = (enext, ithis, eprev)
      newvertices[enext] = in(enext, visited) ? 
        gyromidpoint(newvertices[enext], vertices[ithis], s) : vertices[ithis]
      newvertices[eprev] = in(eprev, visited) ? 
        gyromidpoint(newvertices[eprev], vertices[ithis], s) : vertices[ithis]
      push!(visited, enext, eprev)
    end
    newtriangles[(i-1)*4 + 4] = (
      em[edgeindex(inext, iprev, nv)],
      em[edgeindex(ithis, inext, nv)], 
      em[edgeindex(iprev, ithis, nv)]
    )
  end
  return (vertices = newvertices[1:vcnt], triangles = newtriangles)
end

# gyrotriangle
function gyrotriangle(A, B, C, s; depth=3)
  ABC = [Meshes.coordinates(A), Meshes.coordinates(B), Meshes.coordinates(C)]
  triangle = (1, 2, 3)
  if iseven(depth)
    triangle = reverse(triangle)
  end
  subdiv = gyrosubdiv(s, ABC, [triangle])
  for _ in 1:depth
    subdiv = gyrosubdiv(s, subdiv...)
  end
  vertices = [Meshes.Point(v...) for v in subdiv.vertices]
  connectivities = Meshes.connect.(subdiv.triangles, Meshes.Triangle)
  Meshes.SimpleMesh(vertices, connectivities)
end

# edges of a mesh
function edgeslist(mesh)
  topo = convert(Meshes.HalfEdgeTopology, Meshes.topology(mesh))
  edges_map = Meshes.Boundary{1,0}(topo)
  edges = [edges_map(i) for i in 1:Meshes.nfacets(topo)]
  vertices = mesh.vertices
  [(vertices[edge[1]], vertices[edge[2]]) for edge in edges]
end

# gyrotube, made of spheres
function gyrotube(A, B, s, radius; n = 200)
  centers = gyrosegment(A, B, s; n = n)
  spheres = [
    Meshes.discretize(
      Meshes.Sphere(ctr, radius), Meshes.RegularDiscretization(20)
    )
    for ctr in centers
  ]
  reduce(Meshes.merge, spheres)
end

# merge duplicated vertices of a mesh
function gather(mesh)
  # get mesh vertices
  points = Meshes.vertices(mesh)
  npoints = length(points)

  # Boolean vector indicating the duplicated vertices
  duplicated = fill(false, npoints)

  # dictionary to map the old indices to the new indices
  newindices = Dict{Int64,Int64}()

  # initialize new indices
  newindex = 1

  # iterate over vertices
  for index in 1:(npoints-1)
    if !duplicated[index]
      tail = points[(index+1):npoints]
      duplicates = index .+ findall(==(points[index]), tail)
      duplicated[duplicates] .= true
      push!(duplicates, index)
      for i in duplicates
        newindices[i] = newindex
      end
      newindex = newindex + 1
    end
  end

  # if the last vertex is not a duplicate, add it to the dictionary
  if npoints ∉ keys(newindices) #!haskey(newindices, npoints)
    newindices[npoints] = newindex
  end

  # get the faces
  topo = convert(Meshes.HalfEdgeTopology, Meshes.topology(mesh))
  ∂₂₀ = Meshes.Boundary{2,0}(topo)
  nfaces = Meshes.nelements(topo)

  # vector to store the connectivities
  connec = Vector{Meshes.Connectivity}(undef, 0)

  # iterate over the faces
  for f in 1:nfaces
    face = ∂₂₀(f)
    toconnect = [newindices[i] for i in face]
    c = Meshes.connect(tuple(toconnect...))
    push!(connec, c)
  end

  Meshes.SimpleMesh(points[.!duplicated], connec)
end


"""
  gyromesh(triangles, curvature, edgeradius; depth=2)

Gyromesh from a vector of triangles.

# Arguments
- `triangles`: vector made of vectors of three points
- `curvature`: hyperbolic curvature
- `edgeradius`: radius of the gyroedges (gyrotubes)
- `depth`: number of iterations
"""
function gyromesh(triangles::Vector, curvature, edgeradius; depth = 2)
  meshes = [
    Meshes.SimpleMesh(tr, Meshes.connect.([(1, 2, 3)])) 
    for tr in triangles
  ]
  mesh = gather(reduce(Meshes.merge, meshes)) 
  edges = edgeslist(mesh)
  gyroedges = [
      gyrotube(edge[1], edge[2], curvature, edgeradius)
      for edge in edges
  ]
  gyrotriangles = [
    gyrotriangle(tr..., curvature; depth=depth) 
    for tr in triangles
  ]
  spheres = [
    Meshes.discretize(Meshes.Sphere(pt, edgeradius*1.5)) 
    for pt in Meshes.vertices(mesh)
  ]
  (gyroedges = gyroedges, gyrotriangles = gyrotriangles, spheres = spheres)
end

"""
  gyromesh(m, curvature, edgeradius; depth=2)

Gyromesh from a triangle mesh.

# Arguments
- `m`: a triangle mesh
- `curvature`: hyperbolic curvature
- `edgeradius`: radius of the gyroedges
- `depth`: number of iterations
"""
function gyromesh(m::Meshes.Mesh, curvature, edgeradius; depth = 2)
  triangles = Meshes.vertices.(collect(m))
  gyromesh(triangles, curvature, edgeradius; depth = depth)
end


# test ####
#= 
a = 1
p1 = Meshes.Point(a, -a, -a)
p2 = Meshes.Point(a, a, a)
p3 = Meshes.Point(-a, -a, a)
p4 = Meshes.Point(-a, a, -a)
mesh = Meshes.SimpleMesh(
  [p1, p2, p3, p4], 
  Meshes.connect.([(1,2,3), (1,2,4), (1,3,4), (2,3,4)], Meshes.Triangle)
)

gmesh = gyromesh(mesh, 0.7, 0.01)

function draw(gmesh)
  [MeshViz.viz!(gtriangle; color=:green) for gtriangle in gmesh.gyrotriangles]
  [MeshViz.viz!(gedge; color=:yellow) for gedge in gmesh.gyroedges]
  [MeshViz.viz!(sphere; color=:yellow) for sphere in gmesh.spheres]
  return 0
end
MeshViz.viz(gmesh.gyrotriangles[1]; color = :green)
 =#

C0 = (1 + sqrt(5)) / 4

vertices = Meshes.Point.([
  (0.5, 0.0, C0),
  (0.5, 0.0, -C0),
  (-0.5, 0.0, C0),
  (-0.5, 0.0, -C0),
  (C0, 0.5, 0.0),
  (C0, -0.5, 0.0),
  (-C0, 0.5, 0.0),
  (-C0, -0.5, 0.0),
  (0.0, C0, 0.5),
  (0.0, C0, -0.5),
  (0.0, -C0, 0.5),
  (0.0, -C0, -0.5)
])

faces = map(f -> f .+ 1, [
  [0, 2, 10],
  [0, 10, 5],
  [0, 5, 4],
  [0, 4, 8],
  [0, 8, 2],
  [3, 1, 11],
  [3, 11, 7],
  [3, 7, 6],
  [3, 6, 9],
  [3, 9, 1],
  [2, 6, 7],
  [2, 7, 10],
  [10, 7, 11],
  [10, 11, 5],
  [5, 11, 1],
  [5, 1, 4],
  [4, 1, 9],
  [4, 9, 8],
  [8, 9, 6],
  [8, 6, 2]
])
 
triangles = [vertices[f] for f in faces]
 
icosahedron = Meshes.SimpleMesh(
  vertices, Meshes.connect.([tuple(f...) for f in faces], Meshes.Triangle)
)
 
function drawGI(s, r, alpha1, alpha2)
  gmesh = gyromesh(icosahedron, s, r)
  unitsphere = Meshes.discretize(Meshes.Sphere((0, 0, 0), 1))
  fig, ax, plt = MeshViz.viz( # invisible bounding sphere
    unitsphere; alpha = 0
  )
  ax.show_axis = false
  [MeshViz.viz!(gtriangle; color=:navy) for gtriangle in gmesh.gyrotriangles]
  [MeshViz.viz!(gedge; color=:gold) for gedge in gmesh.gyroedges]
  [MeshViz.viz!(sphere; color=:gold) for sphere in gmesh.spheres]
  Makie.scale!(ax.scene, 1.75, 1.75, 1.75)
  Makie.rotate!(fig.scene, Meshes.Vec3f(1, 0, 0), alpha1)
  Makie.rotate!(ax.scene, Meshes.Vec3f(0, 0, 1), alpha2)
  Makie.save("gyroIcosahedron.png", fig)
end

function animGI(nplots, r, alpha1, alpha2)
  s_ = LinRange(0.6, 0.8, nplots)
  for i in 1:nplots
    gmesh = gyromesh(icosahedron, s_[i], r)
    unitsphere = Meshes.discretize(Meshes.Sphere((0, 0, 0), 1))
    fig, ax, plt = MeshViz.viz( # invisible bounding sphere
      unitsphere; alpha = 0
    )
    ax.show_axis = false
    [MeshViz.viz!(gtriangle; color=:navy) for gtriangle in gmesh.gyrotriangles]
    [MeshViz.viz!(gedge; color=:gold) for gedge in gmesh.gyroedges]
    [MeshViz.viz!(sphere; color=:gold) for sphere in gmesh.spheres]
    Makie.scale!(ax.scene, 1.75, 1.75, 1.75)
    Makie.rotate!(fig.scene, Meshes.Vec3f(1, 0, 0), alpha1)
    Makie.rotate!(ax.scene, Meshes.Vec3f(0, 0, 1), alpha2)
    png = @sprintf "zzpic%03d.png" i
    Makie.save(png, fig)
  end
end