https://github.com/lhf/dk
Tip revision: 55f7c5bfa52cea3f4305e3a23b4e2a8e34bd91c3 authored by lhf on 29 March 2024, 13:38:20 UTC
Update README.md
Update README.md
Tip revision: 55f7c5b
dk.py
# A vertex-centric representation for adaptive diamond-kite meshes
#
# usage: python R INPUT -- all arguments optional
from __future__ import print_function
from sys import argv
from random import random
# depth of refinement
R= int(argv[1]) if len(argv)>1 else 6
# file containing initial mesh
INPUT= argv[2] if len(argv)>2 else "-"
# size of base mesh
N=6
# make array
def array(n):
return list(range(n))
# normalize 3-adic lattice coordinates [a,b,m]=(a+b*w^2)/(3^m)
def normalize(a,b,m):
a,b=int(a),int(b)
while m>0 and a%3==0 and b%3==0:
a=a/3
b=b/3
m=m-1
return [a,b,m]
# add 3-adic lattice points with scaling
def add(v1,v2,n):
a1,b1,m1=v1
a2,b2,m2=v2
m2=m2+n//2
d1=3**m1
d2=3**m2
a=d2*a1+d1*a2
b=d2*b1+d1*b2
m=m1+m2
return normalize(a,b,m)
# multiply 3-adic lattice points
# PolynomialRemainder[(a_1 + b_1 z)(a_2 + b_2 z),z^2-z+1,z]
# z (a_2 b_1 + a_1 b_2 + b_2 b_1) + a_1 a_2 - b_1 b_2
def multiply(v1,v2):
a1,b1,m1=v1
a2,b2,m2=v2
a=a1*a2-b1*b2
b=a2*b1+a1*b2+b2*b1
m=m1+m2
return normalize(a,b,m)
# rotate star
def rotate(a,z):
n=len(a)
b=array(n)
for i in range(n):
b[i]=multiply(a[i],z)
return b
# 3-adic lattice coordinates of basic directions; odd directions are shorter
W=array(12)
W[0]=[1,0,0] # 1
W[1]=[1,1,1] # (1+w^2)/3
W[2]=[0,1,0] # w^2
for i in range(3,12):
W[i]=multiply(W[i-2],W[2])
# stars
STAR=array(6+1)
for d in range(3,6+1):
STAR[d]=array(12)
# basic stars from Eppstein (2014) fig 4
STAR[3][0]=[W[0],W[4],W[8]]
STAR[4][0]=[W[0],W[4],W[7],W[9]]
STAR[5][0]=[W[0],W[3],W[5],W[7],W[9]]
STAR[6][0]=[W[0],W[2],W[4],W[6],W[8],W[10]]
# rotate basic stars
for d in range(3,6+1):
STAR[d][1]=rotate(STAR[d][0],W[1])
for k in range(2,12):
STAR[d][k]=rotate(STAR[d][k-2],W[2])
# opposites
OPP=array(6+1)
for d in range(3,6+1):
OPP[d]=array(12)
OPP[d][0]=array(2*d)
# basic opposites for degree 3
OPP[3][0][0]=[0,1,0]
OPP[3][0][1]=[0,2,0]
OPP[3][0][2]=multiply(OPP[3][0][0],W[4])
OPP[3][0][3]=multiply(OPP[3][0][1],W[4])
OPP[3][0][4]=multiply(OPP[3][0][2],W[4])
OPP[3][0][5]=multiply(OPP[3][0][3],W[4])
# basic opposites for degree 4
OPP[4][0][0]=[0,1,0]
OPP[4][0][1]=[0,2,0]
OPP[4][0][2]=multiply([0,1,0],W[4])
OPP[4][0][3]=[0,0,0]
OPP[4][0][4]=multiply([2,2,1],W[7])
OPP[4][0][5]=multiply([1,1,0],W[7])
OPP[4][0][6]=[1,-1,0]
OPP[4][0][7]=[0,0,0]
# basic opposites for degree 5
OPP[5][0][0]=[0,1,0]
OPP[5][0][1]=[0,0,0]
OPP[5][0][2]=multiply([2,2,1],W[3])
OPP[5][0][3]=multiply([1,1,0],W[3])
OPP[5][0][4]=multiply([2,2,1],W[5])
OPP[5][0][5]=multiply([1,1,0],W[5])
OPP[5][0][6]=multiply([2,2,1],W[7])
OPP[5][0][7]=multiply([1,1,0],W[7])
OPP[5][0][8]=[1,-1,0]
OPP[5][0][9]=[0,0,0]
# basic opposites for degree 6
OPP[6][0][0]=[2,2,1]
OPP[6][0][1]=[1,1,0]
for i in range(2,2*6):
OPP[6][0][i]=multiply(OPP[6][0][i-2],W[2])
# rotate basic opposites
for d in range(3,6+1):
OPP[d][1]=rotate(OPP[d][0],W[1])
for k in range(2,12):
OPP[d][k]=rotate(OPP[d][k-2],W[2])
# Cartesian coordinates
w2=complex(0.5,0.866025403784438646763723170752936183471402626905190314027)
# vertex cloud
V={}
# topological data
U={}
E={}
F={}
nU=0
# use dot notation for dict -- https://stackoverflow.com/a/74214556/107090
class DOTTED(dict):
__getattr__ = dict.get
__setattr__ = dict.__setitem__
__delattr__ = dict.__delitem__
# add vertex to cloud
def addvertex(a,b,m,d,k,n):
global nU
a,b,m=normalize(a,b,m)
if not (a,b,m) in V:
z=(a+b*w2)/(3**m)
t=(a,b,m)
U[nU]=V[a,b,m]=DOTTED({'t':t,'d':d,'k':k,'n':n,'z':z,'w':False,'id':nU})
nU=nU+1
return V[a,b,m]
# add vertex in base mesh
def basevertex(w,d,k):
a,b=w.real,w.imag
v=addvertex(a,b,0,d,k,0)
v.w=w
# base mesh: hexagonal grid
def basemesh(N):
# use complex numbers for easy addition
for i in range(12):
a,b,m=W[i]
W[i]=complex(a,b)
for i in range(N):
for j in range(N):
c=0+i*(W[2]+W[4])+j*3*W[0]
basevertex(c,3,0)
for k in range(0,12,4):
basevertex(c+W[k],6,0)
for k in range(2,12,4):
basevertex(c+W[k],3,2)
if i<N-1 and j<N-1:
basevertex(c+W[0]+W[2],3,0)
boundary()
# base mesh: mark boundary vertices
def boundary():
for k in V:
v=V[k]
if isboundary(v):
v.d=0
v.k=0
def isboundary(v):
for k in range(0,12,2):
w=v.w+W[k]
a,b=w.real,w.imag
a,b,m=normalize(a,b,0)
if not (a,b,m) in V:
return True
return False
# load mesh from csv file
def loadmesh(filename):
N=0
for line in open(filename):
if N>0:
a,b,m,d,k,n=list(map(int,line.split(",")))
addvertex(a,b,m,d,k,n)
N=N+1
addredundant()
# add vertices absent in file
# vertices of degree 3 that are adjacent to a vertex of degree 6
def addredundant():
for k in range(len(U)):
v=U[k]
if v.d==6:
for i in range(v.d):
a,b,m=adj(v,i)
if not (a,b,m) in V:
k=(6+2*i+v.k)%12
addvertex(a,b,m,3,k,v.n)
# stars
def adda(v,i,A):
w=A[v.d][v.k][i]
return add(v.t,w,v.n)
def adj(v,i):
return adda(v,i,STAR)
def adjacent(v,i):
a,b,m=adj(v,i)
return V[a,b,m]
def star(v):
s=array(v.d)
for i in range(v.d):
s[i]=adjacent(v,i)
return s
def makestar(v):
s=array(v.d)
for i in range(v.d):
a,b,m=adj(v,i)
assert not (a,b,m) in V,(a,b,m)
s[i]=addvertex(a,b,m,0,0,0)
return s
# local subdivision step
def subdivide(v):
assert v.d==6,v.d
k=v.k
n=v.n
s0=star(v)
v.k=(k+1)%12
v.n=n+1
s1=makestar(v)
for j in range(v.d):
w=s1[j]
w.d=3
w.k=(6+2*j+k+1)%12
w.n=n+1
for j in range(v.d):
w=s0[j]
kk=(6+2*j+k)%12
if w.d==0:
pass
elif w.d==3:
w.d=4
w.k=(kk+4)%12
elif w.d==4:
w.d=5
if w.k==kk:
w.k=(kk+4)%12
else:
w.k=(kk-4)%12
elif w.d==5:
w.d=6
w.k=(kk-1)%12
w.n=n+1
else:
assert False,w.d
return s0
# opposites
def opp(v,i):
return adda(v,i,OPP)
def opposites(v):
s=array(v.d)
for i in range(0,2*v.d,2):
a,b,m=opp(v,i)
if not (a,b,m) in V:
a,b,m=opp(v,i+1)
s[i//2]=V[a,b,m]
return s
# faces and edges
def addedge(v1,v2):
a=[v1,v2]
j=a.index(min(a))
a=(a[j],a[(j+1)%2])
E[a]=True
def addface(v1,v2,v3,v4):
a=[v1.id,v2.id,v3.id,v4.id]
j=a.index(min(a))
a=(a[j],a[(j+1)%4],a[(j+2)%4],a[(j+3)%4])
F[a]=True
addedge(a[0],a[1])
addedge(a[1],a[2])
addedge(a[2],a[3])
addedge(a[3],a[0])
def addfaces(v):
s=star(v)
o=opposites(v)
for i in range(v.d):
addface(v,s[i],o[i],s[(i+1)%v.d])
# diamond or kite?
def facetype(k):
v1,v2,v3,v4=k
a1,b1,m1=U[v1].t
a2,b2,m2=U[v2].t
a3,b3,m3=U[v3].t
a4,b4,m4=U[v4].t
m=max(m1,m2,m3,m4)
d1=3**(m-m1); a1=a1*d1; b1=b1*d1
d2=3**(m-m2); a2=a2*d2; b2=b2*d2
d3=3**(m-m3); a3=a3*d3; b3=b3*d3
d4=3**(m-m4); a4=a4*d4; b4=b4*d4
if (a1+a3)==(a2+a4) and (b1+b3)==(b2+b4):
return "qd"
else:
return "qk"
# color face according to level
def facetype(k):
v1,v2,v3,v4=k
n=max(U[v1].n,U[v2].n,U[v3].n,U[v4].n)
return str(n)+" qt"
# refinement
queue=set()
# Eppstein's prerequisite structure of replacement operations
# must select all adjacent vertices beforehand: refine(w) can change v.d
# v.d!=6 can happen if object is near the boundary when base mesh is too small
def refine(v):
if v.d==4:
w0=adjacent(v,0)
w1=adjacent(v,1)
refine(w0)
refine(w1)
elif v.d==5:
w0=adjacent(v,0)
refine(w0)
#assert v.d==6,v
if v.d!=6:
return
if v.n>=R:
return
s=subdivide(v)
for w in s:
if w.d==4:
queue.add(w.t)
queue.add(v.t)
# refinement criterion: uniform
def needsrefinement(v):
return v.n < R
# refinement criterion: random
def needsrefinement(v):
return v.n < R and random() < 0.45
# refinement criterion: implicit curve
def needsrefinement(v):
z = min([f(w)*f(v) for w in star(v)])
return v.n < R and z<=0
# implicit curve (Taubin, 1994)
def f(v):
z=v.z-complex(8,4)
x,y=z.real,z.imag
z=0.004+0.110*x-0.177*y-0.174*x*x+0.224*x*y-0.303*y*y-0.168*x*x*x+0.327*x*x*y-0.087*x*y*y-0.013*y*y*y+0.235*x*x*x*x-0.667*x*x*x*y+0.745*x*x*y*y-0.029*x*y*y*y+0.072*y*y*y*y;
return z
# reduction: vertices of degree 3 that are adjacent to a vertex of degree 6
# same as vertices of degree 3 that are adjacent to an internal vertex
def redundant(v):
if v.d!=3:
return False
for w in star(v):
#if w.d!=0:
if w.d==6:
return True
return False
# main -----------------------------------------------------------------------
# initial mesh
if INPUT=="-":
basemesh(N)
else:
loadmesh(INPUT)
# refine mesh
queue = {k for k in V if V[k].d==6}
while queue:
k=next(iter(queue))
queue.remove(k)
v=V[k]
if needsrefinement(v):
refine(v)
# reconstruct mesh
for k in V:
if V[k].d>0:
addfaces(V[k])
# output PostScript
print('''\
%!PS-Adobe-2.0 EPSF-2.0
%%BoundingBox: 0 0 500 300
%%BoundingBox: 20 25 455 295
%%BoundingBox: 135 80 390 240
50 50 translate
25 dup scale
0 setlinewidth
1 setlinejoin
1 setlinecap
/c { 0.02 0 360 arc fill } bind def
/p { 0.04 0 360 arc fill } bind def
/l { moveto lineto stroke } bind def
/q { moveto lineto lineto lineto closepath fill } bind def
/p0 { 0.5 0.5 0.5 setrgbcolor p } bind def
/p3 { 0 0 1 setrgbcolor p } bind def
/p4 { 0 1 0 setrgbcolor p } bind def
/p5 { 1 0 1 setrgbcolor p } bind def
/p6 { 1 0 0 setrgbcolor p } bind def
/qd { 1 0.8 0.8 setrgbcolor q } bind def
/qk { 0.8 0.8 1 setrgbcolor q } bind def
/qt { 0.9 6 div mul 0.1 add 1 sub neg dup 1 setrgbcolor q } bind def
''')
print("%=ARG",R,INPUT)
print("%=VEF",R,len(V),len(U),len(E),len(F))
print("")
print("% faces")
print("1 0.8 0.8 setrgbcolor")
for k in F:
z1=U[k[0]].z
z2=U[k[1]].z
z3=U[k[2]].z
z4=U[k[3]].z
print(z1.real,z1.imag, z2.real,z2.imag, z3.real,z3.imag, z4.real,z4.imag,facetype(k))
print("")
print("% edges")
print("0 0 0 setrgbcolor")
print("0.01 setlinewidth")
for k in E:
z1=U[k[0]].z
z2=U[k[1]].z
print(z1.real,z1.imag, z2.real,z2.imag, "l")
print("")
print("% curve")
print("1 0 0 setrgbcolor")
for k in E:
v0=U[k[0]]; z0=v0.z; f0=f(v0)
v1=U[k[1]]; z1=v1.z; f1=f(v1)
if f0*f1<=0:
t=(f0-0)/(f0-f1)
print((1-t)*z0.real+t*z1.real,(1-t)*z0.imag+t*z1.imag, "c")
"""
print ""
print "% vertices"
print "0 0 0 setrgbcolor"
for k in U:
v=U[k]
print v.z.real, v.z.imag, "p" + str(v.d) # + "0" # + str(v.k)
"""
print("")
print("showpage")
print("%%EOF")
#exit() # if extra files not needed
# output CSV
print("")
PREFIX="%=CSV"
print(PREFIX, "a,b,m,d,k,n")
for k in range(len(U)):
v=U[k]
print(PREFIX, ','.join(map(str,[v.t[0],v.t[1],v.t[2],v.d,v.k,v.n])))
# output OBJ
print("")
PREFIX="%=OBJ"
print(PREFIX, "OBJ")
print(PREFIX, len(U),len(F),0)
for k in range(len(U)):
v=U[k]
print(PREFIX, "v", v.z.real, v.z.imag, 0)
for k in F:
print(PREFIX, "f",k[0]+1,k[1]+1,k[2]+1,k[3]+1)
# output OFF
print("")
PREFIX="%=OFF"
print(PREFIX, "OFF")
print(PREFIX, len(U),len(F),0)
for k in range(len(U)):
v=U[k]
print(PREFIX, v.z.real, v.z.imag, 0)
for k in F:
print(PREFIX, "4",k[0],k[1],k[2],k[3])
# done
print("")
print("%=EOF")