804 lines
30 KiB
Plaintext
804 lines
30 KiB
Plaintext
type Face = number[]
|
|
type XYZ = [number, number, number]
|
|
type Polyhedron = {face: Face[], xyz: XYZ[], name?: string}
|
|
|
|
type Notation = string
|
|
|
|
// Useful constants
|
|
rt2 := Math.sqrt 2
|
|
rth := rt2/2
|
|
phi := (1 + Math.sqrt 5)/2
|
|
ihp := 1/phi
|
|
tau := 2*Math.PI
|
|
|
|
// Sadly needs to be early because we initialize the Dodecahedron as the dual of
|
|
// the icosahedron:
|
|
// Only add one direction of each, will auto reverse as well
|
|
specialDuals: Record<string, string> :=
|
|
Tetrahedron: 'Tetrahedron',
|
|
Cube: 'Octahdron',
|
|
Dodecahedron: 'Icosahedron',
|
|
Cuboctahedron: 'Rhombic dodecahedron',
|
|
'truncated Tetrahedron': 'triakis Tetrahedron',
|
|
'truncated Cube': 'triakis Octahedron',
|
|
'truncated Octahedron': 'tetrakis Cube',
|
|
Rhombicuboctahedron: 'Deltoidal icositetrahedron',
|
|
'truncated Cuboctahedron': 'Disdyakis dodecahedron',
|
|
'snub Cube': 'pentagonal icositetrahedron',
|
|
Icosidodecahedron: 'Rhombic triacontahedron',
|
|
'truncated Dodecahedron': 'triakis Icosahedron',
|
|
'truncated Icosahedron': 'pentakis Dodecahedron',
|
|
Rhombicosidodecahedron: 'Deltoidal hexecontahedron',
|
|
'truncated Icosidodecahedron': 'Disdyakis triacontahedron',
|
|
'snub Dodecahedron': 'pentagonal hexecontahedron'
|
|
|
|
// All seeds are canonical, i.e., edges tangent to unit sphere
|
|
icosahedron: Polyhedron :=
|
|
face: [[0,1,9], [0,8,1], [0,4,8], [0,5,4], [0,9,5],
|
|
[4,5,2], [4,2,10], [4,10,8], [8,10,6], [8,6,1],
|
|
[1,6,7], [1,7,9], [9,7,11], [9,11,5], [5,11,2],
|
|
[3,2,11], [3,10,2], [3,6,10], [3,7,6], [3,11,7]]
|
|
xyz: [[0,ihp,1], [0,-ihp,1], [0,ihp,-1], [0,-ihp,-1],
|
|
[ihp,1,0], [-ihp,1,0], [ihp,-1,0], [-ihp,-1,0],
|
|
[1,0,ihp], [-1,0,ihp], [1,0,-ihp], [-1,0,-ihp]]
|
|
name: 'Icosahedron'
|
|
|
|
polyCache: Record<Notation, Polyhedron> :=
|
|
'': face: [], xyz: []
|
|
T:
|
|
face: [[0,2,1], [0,3,2], [0,1,3], [1,2,3]]
|
|
xyz: [[1,1,1], [1,-1,-1], [-1,1,-1], [-1,-1,1]]
|
|
name: 'Tetrahedron'
|
|
O:
|
|
face: [[0,2,1], [0,3,2], [0,4,3], [0,1,4],
|
|
[1,5,4], [1,2,5], [2,3,5], [3,4,5]]
|
|
xyz: [[0,0,rt2], [rt2,0,0], [0,rt2,0], [-rt2,0,0], [0,-rt2,0], [0,0,-rt2]]
|
|
name: 'Octahedron'
|
|
C:
|
|
face: [[0,3,2,1], [0,5,4,3], [0,1,6,5], [1,2,7,6], [2,3,4,7], [4,5,6,7]]
|
|
xyz: [[rth,rth,rth], [-rth,rth,rth], [-rth,-rth,rth], [rth,-rth,rth],
|
|
[rth,-rth,-rth], [rth,rth,-rth], [-rth,rth,-rth], [-rth,-rth,-rth]]
|
|
name: 'Cube'
|
|
I: icosahedron
|
|
D: geomDual icosahedron
|
|
|
|
export function generateVRML(notation: Notation): string
|
|
outputVRML notation, generatePoly notation
|
|
|
|
function generatePoly(notation: Notation): Polyhedron
|
|
getStandardPoly inform standardize notation
|
|
|
|
function getStandardPoly(notation: Notation): Polyhedron
|
|
if notation in polyCache then return polyCache[notation]
|
|
[ , op, arg, rest] := notation.match(/^(.)(\d*)(.*)$/) or []
|
|
parent := getStandardPoly rest
|
|
// may have created what we want by side effect
|
|
if notation in polyCache then return polyCache[notation]
|
|
dispatch op, arg, parent, notation // will do the caching
|
|
|
|
// Convenience tuple maker
|
|
function ð<Tup extends unknown[]>(...arg: Tup): Tup arg
|
|
|
|
// Note we now allow numeric arguments on all of the basic operations,
|
|
// kis/truncate, join/ambo, and gyro/snub. Likely some of the operations
|
|
// we are taking as composite could have reasonable numeric versions, but
|
|
// there didn't seem to be any sensible way to propagate such an argument
|
|
// to the operations in their rewrites. In other words, the numeric-limited
|
|
// operations may not be composite, or at least not in the same way. So
|
|
// we have just left them as applying throughout the polyhedron.
|
|
rawStandardizations :=
|
|
P4$: 'C', A3$: 'O', Y3$: 'T', // Seed synonyms
|
|
e: 'aa', b: 'ta', o: 'jj', m: 'kj', // abbreviations
|
|
[String.raw`t(\d*)`]: 'dk$1d',
|
|
[String.raw`a(\d*)`]: 'dj$1d', // dual operations
|
|
[String.raw`s(\d*)`]: 'dg$1d',
|
|
dd: '', rr: '', jd: 'j', gd: 'rgr', // absorption rules
|
|
rd: 'dr', // these commute; others? If so, move 'r' in to cancel w/ seed
|
|
// Remainder are all simplifications/unique selections for seeds:
|
|
jY: 'dA', dT: 'T', gT: 'D', jT: 'C', dC: 'O', dO: 'C',
|
|
dI: 'D', dD: 'I', rO: 'O', rC: 'C', rI: 'I', rD: 'D',
|
|
jO: 'jC', jI: 'jD', gO: 'gC', gI: 'gD'
|
|
standardizations :=
|
|
(ð RegExp(pat, 'g'), rep for pat, rep in rawStandardizations)
|
|
|
|
function standardize(notation: Notation): Notation
|
|
lastNotation .= ''
|
|
while lastNotation != notation // iterate in case of rdrd, e.g.
|
|
lastNotation = notation
|
|
for [pat, rep] of standardizations
|
|
notation = notation.replace(pat, rep)
|
|
notation
|
|
|
|
function orb(r: number, n: number,
|
|
height: number, t = 0): XYZ[]
|
|
// A regular n-gon inscribed in a horizontal circle of radius r
|
|
// at given height, rotated t*tau/n from standard position.
|
|
theta := tau/n
|
|
rho := t*theta
|
|
for i of [0...n]
|
|
[r*Math.cos(rho + i*theta), r*Math.sin(rho + i*theta), height] as XYZ
|
|
|
|
operator add(a: XYZ, b: XYZ)
|
|
accumulate copy(a), b
|
|
|
|
ngonalNames := 3: 'triangular', 4: 'square', 5: 'pentagonal', 6: 'hexagonal',
|
|
7: 'heptagonal', 8: 'octagonal', 9: 'enneagonal', 10: 'decagonal',
|
|
12: 'dodecagonal'
|
|
type SpecialNgon = keyof typeof ngonalNames
|
|
function ngonal(n: number)
|
|
ngonalNames[n as SpecialNgon] ?? n+`-gonal`
|
|
|
|
seeds :=
|
|
P: (n: number) => // Prism
|
|
unless n then n = 3
|
|
theta := tau/n
|
|
halfEdge := Math.sin theta/2
|
|
xyz := orb(1, n, halfEdge) ++ orb(1, n, -halfEdge)
|
|
face := [[n-1..0], [n...2*n]] // top and bottom
|
|
for i of [0...n]
|
|
ip1 := (i+1)%n // next vertex around
|
|
face.push [i, ip1, ip1+n, i+n]
|
|
{face, xyz, name: ngonal(n)+` prism`}
|
|
A: (n: number) => // Antiprism
|
|
unless n then n = 4
|
|
theta := tau/n
|
|
halfHeight .= Math.sqrt
|
|
1 - 4/(4 + 2*Math.cos(theta/2) - 2*Math.cos(theta))
|
|
faceRadius .= Math.sqrt 1-halfHeight*halfHeight
|
|
// Scale to put edge midpoints on unit sphere
|
|
f := mag [halfHeight, faceRadius*Math.cos(theta/2), 0]
|
|
halfHeight /= f
|
|
faceRadius /= f
|
|
xyz := orb(faceRadius, n, halfHeight)
|
|
++ orb(faceRadius, n, -halfHeight, 0.5)
|
|
face := [[n-1..0], [n...2*n]] // top and bottom
|
|
for i of [0...n]
|
|
face.push [i, (i+1)%n, i+n]
|
|
face.push [i, i+n, n + (n+i-1)%n]
|
|
{face, xyz, name: ngonal(n)+` antiprism`}
|
|
Y: (n: number) => // pYramid
|
|
unless n then n = 4
|
|
// Canonical solution by Intelligenti Pauca and Ed Pegg, see
|
|
// https://math.stackexchange.com/questions/2286628/canonical-pyramid-polynomials
|
|
theta := tau/n
|
|
c := Math.cos theta/2
|
|
baseRadius := Math.sqrt 2/(c*(1+c))
|
|
depth := Math.sqrt (1-c)/(1+c)
|
|
height := 2*Math.sqrt 1/(1 - c*c)
|
|
xyz := orb baseRadius, n, depth
|
|
edgeMid2 := xyz[0] add xyz[1]
|
|
xyz.push [0, 0, depth-height]
|
|
face := ([i, (i+1)%n, n] for i of [0...n])
|
|
face.unshift [n-1..0]
|
|
{face, xyz, name: ngonal(n)+` pyramid`}
|
|
type SeedOp = keyof typeof seeds
|
|
|
|
// Syntactic sugar to deal with weird TypeScript typing:
|
|
operator æ<T>(A: T[], i: number) A.at(i) as T
|
|
function wordsof(s: string) 1 + (s.match(/\s+/g)?.length ?? 0)
|
|
specialJoins :=
|
|
Tetrahedron: 'Cube',
|
|
Cube: 'Rhombic dodecahedron', Octahedron: 'Rhombic dodecahedron',
|
|
Dodecahedron: 'Rhombic triacontahedron', Icosahedron: 'Rhombic triacontahedron',
|
|
'Rhombic dodecahedron': 'Deltoidal icositetrahedron',
|
|
'Rhombic triacontahedron': 'Deltoidal hexecontahedron',
|
|
'pentakis Dodecahedron': 'Rhombic enneacontahedron'
|
|
type SpecialJoin = keyof typeof specialJoins
|
|
specialKis :=
|
|
'Rhombic dodecahedron': 'Disdyakis dodecahedron'
|
|
'Rhombic triacontahedron': 'Disdyakis triacontahedron'
|
|
|
|
type SpecialKis = keyof typeof specialKis
|
|
kisWords:= 3: 'triakis', 4: 'tetrakis', 5: 'pentakis', 6: 'hexakis', 7: 'heptakis'
|
|
type KisNumber = keyof typeof kisWords
|
|
|
|
function kisjoin(P: Polyhedron, notation: string,
|
|
digits: string, join: boolean): Polyhedron
|
|
// kis and join are closely related operations. Both of them add a
|
|
// pyramid on a selection of faces; join then further deletes any
|
|
// _original_ edge bordered by two _new_ triangles, producing a quad.
|
|
// Faces are selected by their numbers of sides, using the given digits.
|
|
// If there are none, all faces are used. Otherwise, the digits are turned
|
|
// into a list of numbers by breaking after every digit except as needed
|
|
// to prevent leading 0s or isolated 1s or 2s (since no face has one or
|
|
// two sides); this way you can list any subset of the numbers 3 - 32,
|
|
// which is plenty.
|
|
// The operation is then applied just to faces with the numbers of edges on
|
|
// the list. e.g. k3412 will add pyramids to the triangles, quads, and
|
|
// dodecagon faces.
|
|
allowed := parseSides digits
|
|
// first collect a directory from face indices to new vertex numbers
|
|
nextVertex .= P.xyz.length
|
|
newVixes :=
|
|
for each f of P.face
|
|
!digits or f.length is in allowed ? nextVertex++ : 0
|
|
if nextVertex is P.xyz.length then return P // nothing to do
|
|
xyz := P.xyz ++ faceCenters(P).filter (f,ix) => newVixes[ix]
|
|
face: Face[] := []
|
|
for each f, ix of P.face
|
|
v := newVixes[ix]
|
|
if v is 0
|
|
face.push f.slice()
|
|
continue
|
|
// Add the pyramid, possibly eliding edges:
|
|
for each w, jx of f
|
|
pw := f æ jx-1
|
|
neighbor .= 0
|
|
if join
|
|
neighbor = P.face.findIndex (g, gx) =>
|
|
gx !== ix and w is in g and pw is in g
|
|
if join and newVixes[neighbor] // elide this edge
|
|
if pw < w // avoid adding same face twice
|
|
face.push [v, pw, newVixes[neighbor], w]
|
|
else face.push [v, pw, w]
|
|
// Nomenclature
|
|
let name: string|undefined
|
|
if P.name and wordsof(P.name) < 3 // don't go hog-wild with naming
|
|
if join
|
|
unless digits
|
|
name = specialJoins[P.name as SpecialJoin] ?? 'joined ' + P.name
|
|
else
|
|
size .= P.face[0].length
|
|
unless P.face.every (f) => f.length is size
|
|
size = 0
|
|
if !size and digits then size = Number(digits) or 0
|
|
// very special case
|
|
if size is 5 and P.name is 'pentagonal antiprism'
|
|
name = 'Icosahedron'
|
|
else if (!digits or Number(digits) is size) and size in kisWords
|
|
name = (specialKis[P.name as SpecialKis]
|
|
?? kisWords[size as KisNumber] + ' ' + P.name)
|
|
// Cheaty super special case
|
|
if notation is 'jk5djP5'
|
|
name = 'dual elongated pentagonal orthobirotunda'
|
|
// Done, fix up the vertices a bit
|
|
adjustXYZ {face, xyz, name}, notation, 3
|
|
|
|
|
|
specialGyro :=
|
|
Cube: 'pentagonal icositetrahedron', Octahedron: 'pentagonal icositetrahedron',
|
|
Icosahedron: 'pentagonal hexecontahedron',
|
|
Dodecahedron: 'pentagonal hexecontahedron'
|
|
type SpecialGyro = keyof typeof specialGyro
|
|
|
|
// how enums ought to work?
|
|
FromCenter := Symbol()
|
|
AlongEdge := Symbol()
|
|
type Gyway = typeof FromCenter | typeof AlongEdge
|
|
|
|
function gyropel(P: Polyhedron, notation: string,
|
|
digits: string, ...ways: Gyway[]): Polyhedron
|
|
// gyro and propellor are closely related operations. Both of them add new
|
|
// vertices one third of the way along each edge of each face selected
|
|
// by the digits argument (see kisjoin). They then differ in what edges
|
|
// are drawn to the new vertices. In gyro, another new vertex is added
|
|
// at the center of each face and connected to each of them; in propellor,
|
|
// they are just connected in sequence. For completeness, we also allow
|
|
// both at the same time, which is equivalent to propellor followed by kis
|
|
// just on the new rotated faces.
|
|
fromCenter := FromCenter is in ways
|
|
alongEdge := AlongEdge is in ways
|
|
unless fromCenter or alongEdge then return P // nothing to do
|
|
allowed := parseSides digits
|
|
// first collect a directory from directed edges to new vertex numbers
|
|
xyz := P.xyz.slice()
|
|
startV := xyz.length
|
|
edgeV: number[][] := []
|
|
for each f of P.face
|
|
if digits and f.length is not in allowed then continue
|
|
for each v, ix of f
|
|
pv := f æ ix-1
|
|
(edgeV[pv] ??= [])[v] = xyz.length
|
|
xyz.push lerp xyz[pv], xyz[v], 1/3
|
|
if xyz.length is startV then return P // nothing to do
|
|
|
|
// Now revisit each face, accumulating the new faces.
|
|
face: Face[] := []
|
|
centers: XYZ[] := fromCenter ? faceCenters P : []
|
|
for each f, fx of P.face
|
|
if digits and f.length is not in allowed
|
|
// Just collect all of the vertices around
|
|
newFace: Face := []
|
|
for each v, ix of f
|
|
pv := f æ ix-1
|
|
reverseV := edgeV[v]?[pv] ?? -1
|
|
if reverseV >= 0 then newFace.push reverseV
|
|
newFace.push v
|
|
face.push newFace
|
|
continue
|
|
centerV := xyz.length
|
|
if fromCenter then xyz.push centers[fx]
|
|
aroundOutside: Face .= []
|
|
for each v, ix of f
|
|
pv := f æ ix-1
|
|
ppv := f æ ix-2
|
|
firstNew := edgeV[ppv][pv]
|
|
newSection := [firstNew]
|
|
reverseV := edgeV[pv][ppv] ?? -1
|
|
if reverseV >= 0 then newSection.push reverseV
|
|
newSection.push pv
|
|
secondNew := edgeV[pv][v]
|
|
if alongEdge
|
|
newSection.push secondNew
|
|
face.push newSection
|
|
if fromCenter then face.push ð centerV, firstNew, secondNew
|
|
else aroundOutside.push firstNew
|
|
else if fromCenter
|
|
newSection.push secondNew, centerV
|
|
face.push newSection
|
|
else aroundOutside ++= newSection
|
|
if aroundOutside.length then face.push aroundOutside
|
|
let name: string|undefined
|
|
nw .= 0
|
|
if ((fromCenter xor alongEdge)
|
|
and P.name and !digits and (nw = wordsof(P.name)) < 3)
|
|
if alongEdge
|
|
if nw is 1 then name = 'propello' + P.name
|
|
else name = 'propellorized '+P.name
|
|
else
|
|
if P.name in specialGyro then name = specialGyro[P.name as SpecialGyro]
|
|
else if nw is 1 then name = 'gyro' + P.name
|
|
else name = 'gyrated '+P.name
|
|
// Done, fix up the vertices a bit
|
|
adjustXYZ {face, xyz, name}, notation, 3
|
|
|
|
function parseSides(digits: string): number[]
|
|
unless digits return []
|
|
tooSmall := ['1', '2']
|
|
last := digits.length - 1
|
|
return := []
|
|
current .= ''
|
|
pos .= 0
|
|
while pos <= last
|
|
current += digits[pos++]
|
|
nextDigit := digits[pos]
|
|
if (current is in tooSmall
|
|
or nextDigit is '0'
|
|
or pos == last and nextDigit is in tooSmall)
|
|
continue
|
|
return.value.push parseInt current
|
|
current = ''
|
|
|
|
transforms :=
|
|
k: (P: Polyhedron, notation: string, digits: string) => // kis[n]
|
|
kisjoin P, notation, digits, false
|
|
j: (P: Polyhedron, notation: string, digits: string) => // join
|
|
kisjoin P, notation, digits, true
|
|
g: (P: Polyhedron, notation: string, digits: string) => // gyro
|
|
gyropel P, notation, digits, FromCenter
|
|
p: (P: Polyhedron, notation: string, digits: string) => // propellor
|
|
gyropel P, notation, digits, AlongEdge
|
|
f: (P: Polyhedron, notation: string, digits: string) => // fan [new? name?]
|
|
gyropel P, notation, digits, AlongEdge, FromCenter
|
|
r: (P: Polyhedron) => // reverse (mirror)
|
|
face: (f.toReversed() for each f of P.face)
|
|
xyz: (scale copy(v), -1 for each v of P.xyz)
|
|
name: P.name // should we record that it's mirrored?
|
|
d: (P: Polyhedron, notation: string, digits: string) => // dual
|
|
if digits
|
|
console.error `Ignoring ${digits} arg of d in ${notation}`
|
|
// Create a "fake" of P and adjust it (so we don't disturb the
|
|
// cached P), which will create the dual
|
|
parentNotation := notation[1+digits.length..]
|
|
adjustXYZ {face: P.face.slice(), P.xyz, P.name}, parentNotation, 1
|
|
polyCache.`d${parentNotation}`
|
|
c: (P: Polyhedron, notation: string, digits: string) => // canonicalize
|
|
face: P.face.slice()
|
|
xyz: canonicalXYZ P, notation, Number(digits) or 10
|
|
name: P.name
|
|
x: approxCanonicalize // iterative direct adjustment algorithm
|
|
type TransformOp = keyof typeof transforms
|
|
|
|
function dispatch(op: string, digits: string,
|
|
P: Polyhedron, mynotation: string): Polyhedron
|
|
// Note mynotation starts with op!
|
|
return .= P
|
|
if op in seeds
|
|
return = seeds[op as SeedOp] Number(digits) or 0
|
|
else if op in transforms
|
|
return = transforms[op as TransformOp] P, mynotation, digits
|
|
else
|
|
console.error `Unknown operation ${op}${digits} in ${mynotation}.`
|
|
return = polyCache.T
|
|
polyCache[mynotation] = return.value
|
|
|
|
function topoDual(P: Polyhedron): Polyhedron
|
|
// Note this maintains correspondence between V and F indices, but
|
|
// places every vertex at the origin! Don't use without geometrizing
|
|
// in some way.
|
|
face:
|
|
for v of [0...P.xyz.length]
|
|
infaces := // gather labeled list of faces contining v
|
|
for f, index of P.face
|
|
unless f.includes v continue
|
|
ð f, index
|
|
start := infaces[0][1];
|
|
current .= start
|
|
newface := []
|
|
do
|
|
verts := P.face[current]
|
|
preV := verts æ verts.indexOf(v)-1
|
|
nextIx := infaces.findIndex ([face, label]) =>
|
|
label !== current and face.includes preV
|
|
current = infaces[nextIx][1]
|
|
newface.push current
|
|
if newface.length > infaces.length
|
|
console.error 'In topoDual: Malformed polyhedron', P
|
|
break
|
|
until current is start
|
|
newface
|
|
xyz:
|
|
Array(P.face.length).fill([0,0,0]) // warning, every vertex is ===
|
|
name: dualName P.name
|
|
|
|
|
|
function geomDual(P: Polyhedron): Polyhedron
|
|
return := topoDual P
|
|
return.value.xyz = approxDualVertices P
|
|
|
|
function approxDualVertices(P: Polyhedron): XYZ[]
|
|
P.face.map (f) => approxDualVertex f, P.xyz
|
|
|
|
operator dot same (/) (v: number[], w: number[])
|
|
v.reduce (l,r,i) => l + r*w[i], 0
|
|
|
|
function approxDualVertex(f: Face, v: XYZ[]): XYZ
|
|
// For each edge of f, there is a plane containing it perpendicular
|
|
// to the line joining the origin to its nearest approach to the origin.
|
|
// This function returns the point closest to being on all of those planes
|
|
// (in the least-squares sense).
|
|
// This method seems to work OK when the neighborhood of f is convex,
|
|
// and very poorly otherwise. So it probably would not provide any better
|
|
// canonicalization than other methods of approximating the dual.
|
|
// It might be better to replace this with the approximate intersection of the
|
|
// reciprocal planes to each of the vertices of the face...
|
|
normals := (tangentPoint(v[f æ i-1], v[f[i]]) for i of [0...f.length])
|
|
sqlens := normals.map mag2
|
|
columns := (normals.map(&[i]) for i of [0..2])
|
|
target := (columns[i] dot sqlens for i of [0..2]) as XYZ
|
|
CMsource := (for c of [0..2]
|
|
(columns[r] dot columns[c] for r of [0..2])) as [XYZ, XYZ, XYZ]
|
|
cramerD := det ...CMsource
|
|
if Math.abs(cramerD) < 1e-6
|
|
console.error `Face ${f} of ${v.map (p) => '['+p+']'} ill conditioned`
|
|
return [0, 0, 0]
|
|
[ det(target,CMsource[1],CMsource[2])/cramerD,
|
|
det(CMsource[0],target,CMsource[2])/cramerD,
|
|
det(CMsource[0],CMsource[1],target)/cramerD ]
|
|
|
|
function det(a: XYZ, b: XYZ, c:XYZ)
|
|
a[0]*b[1]*c[2] + a[1]*b[2]*c[0] + a[2]*b[0]*c[1]
|
|
- a[2]*b[1]*c[0] - a[1]*b[0]*c[2] - a[0]*b[2]*c[1]
|
|
|
|
operator sub(a: XYZ, b: XYZ)
|
|
diminish copy(a), b
|
|
|
|
function tangentPoint(v: XYZ, w: XYZ) // closest point on vw to origin
|
|
d := w sub v
|
|
v sub scale d, d dot v / mag2 d
|
|
|
|
function faceCenters(P: Polyhedron): XYZ[]
|
|
for each face of P.face
|
|
scale face.reduce((ctr,v) => accumulate(ctr, P.xyz[v]), [0,0,0]),
|
|
1/face.length
|
|
|
|
function dualName(p: string|undefined)
|
|
unless 'Octahedron' in specialDuals
|
|
// one-time reversal of all special duals
|
|
specialDuals[dual] = poly for poly, dual in specialDuals
|
|
|
|
unless p return undefined
|
|
if p in specialDuals
|
|
return specialDuals[p]
|
|
words := p.split(' ')
|
|
if words[0] is 'dual' return words[1..].join ' '
|
|
if words.length is 2
|
|
switch words[1]
|
|
'prism'
|
|
return words[0] + ' bipyramid'
|
|
'bipyramid'
|
|
return words[0] + ' prism'
|
|
'antiprism'
|
|
return words[0] + ' trapezohedron'
|
|
'trapezohedron'
|
|
return words[0] + ' antiprism'
|
|
'pyramid'
|
|
return p // self-dual
|
|
return 'dual ' + p
|
|
|
|
function adjustXYZ(P: Polyhedron, notation: string, iterations = 1): Polyhedron
|
|
dualNotation := 'd' + notation
|
|
D .= topoDual P
|
|
if dualNotation in polyCache
|
|
console.error 'Creating', notation, '_after_ its dual'
|
|
D = polyCache[dualNotation]
|
|
for iter of [1..iterations]
|
|
D.xyz = reciprocalC P
|
|
P.xyz = reciprocalC D
|
|
polyCache[dualNotation] = D
|
|
P
|
|
|
|
function reciprocalC(P: Polyhedron): XYZ[]
|
|
return := faceCenters P
|
|
for each v of return.value
|
|
scale v, 1/mag2 v
|
|
|
|
function canonicalXYZ(P: Polyhedron, notation: string, iterations: number): XYZ[]
|
|
dualNotation := 'd' + notation
|
|
D .= topoDual P
|
|
if dualNotation in polyCache
|
|
console.error 'Creating', notation, '_after_ its dual'
|
|
D = polyCache[dualNotation]
|
|
tempP := Object.assign({}, P) // algorithm is read-only on original data
|
|
if iterations < 1 then iterations = 1
|
|
for iter of [1..iterations]
|
|
D.xyz = reciprocalN tempP
|
|
tempP.xyz = reciprocalN D
|
|
polyCache[dualNotation] = D
|
|
tempP.xyz
|
|
|
|
operator cross tighter (*) (v: XYZ, w: XYZ): XYZ
|
|
[ v[1]*w[2] - v[2]*w[1], v[2]*w[0] - v[0]*w[2], v[0]*w[1] - v[1]*w[0] ]
|
|
|
|
function reciprocalN(P: Polyhedron): XYZ[]
|
|
for each f of P.face
|
|
centroid := ð 0, 0, 0
|
|
normal := ð 0, 0, 0
|
|
meanEdgeRadius .= 0
|
|
for each vlabel, ix of f
|
|
v := P.xyz[vlabel]
|
|
accumulate centroid, v
|
|
pv := P.xyz[f æ ix-1]
|
|
ppv := P.xyz[f æ ix-2]
|
|
// original doc says unit normal but below isn't. Didn't help to try, tho
|
|
nextNormal := (pv sub ppv) cross (v sub pv)
|
|
// or instead, just chop down big ones: (didn't work either)
|
|
// magNext := mag nextNormal
|
|
// if magNext > 1 then scale nextNormal, 1/magNext
|
|
accumulate normal, nextNormal
|
|
meanEdgeRadius += mag tangentPoint pv, v
|
|
scale centroid, 1/f.length
|
|
scale normal, 1/mag normal
|
|
meanEdgeRadius /= f.length
|
|
scale normal, centroid dot normal
|
|
scale normal, 1/mag2 normal // invert in unit sphere
|
|
scale normal, (1 + meanEdgeRadius)/2
|
|
|
|
operator dist(a: XYZ, b: XYZ) mag a sub b
|
|
|
|
// adapted from Polyhedronisme https://levskaya.github.io/polyhedronisme/
|
|
function approxCanonicalize(P: Polyhedron, notation: string,
|
|
digits: String): Polyhedron
|
|
THRESHOLD := 1e-6
|
|
// A difficulty is that the planarizing sometimes has the effect of
|
|
// "folding over" neighboring faces, in which case all is lost.
|
|
// Keeping the weight of the edge smoothing high compared to planarizing
|
|
// seems to help with that.
|
|
EDGE_SMOOTH_FACTOR := 0.5
|
|
PLANARIZE_FACTOR := 0.1
|
|
edge := edges P
|
|
xyz := P.xyz.map copy
|
|
V := xyz.length
|
|
normalizeEdges xyz, edge
|
|
for iter of [1..Number(digits) or 10]
|
|
start := xyz.map copy
|
|
smoothEdgeDists xyz, edge, EDGE_SMOOTH_FACTOR
|
|
normalizeEdges xyz, edge
|
|
planarize xyz, P.face, PLANARIZE_FACTOR
|
|
normalizeEdges xyz, edge
|
|
if Math.max(...(xyz[i] dist start[i] for i of [0...V])) < THRESHOLD
|
|
break
|
|
{face: P.face.slice(), xyz, P.name}
|
|
|
|
type Edge = [number, number]
|
|
function edges(P:Polyhedron): Edge[]
|
|
return: Edge[] := []
|
|
for each f of P.face
|
|
for each v, ix of f
|
|
pv := f æ ix-1
|
|
if pv < v then return.value.push ð pv, v
|
|
|
|
function normalizeEdges(xyz: XYZ[], edge: Edge[]): void
|
|
// Adjusts xyz so that edge tangentpoints have centroid at origin and
|
|
// mean radius 1
|
|
edgeP .= edge.map ([a,b]) => tangentPoint xyz[a], xyz[b]
|
|
edgeCentroid := centroid edgeP
|
|
xyz.forEach (pt) => diminish pt, edgeCentroid
|
|
edgeScale := 1/(mean edge.map ([a,b]) => mag tangentPoint xyz[a], xyz[b])
|
|
xyz.forEach (pt) => scale pt, edgeScale
|
|
|
|
function centroid(xyz: XYZ[]): XYZ
|
|
scale xyz.reduce(accumulate, ð 0,0,0), 1/xyz.length
|
|
|
|
function smoothEdgeDists(xyz: XYZ[], edge: Edge[], fudge: number): void
|
|
// Attempts in the most straightforward way possible to reduce the
|
|
// variance of the radii of the edgepoints
|
|
V := xyz.length
|
|
adj := (ð 0,0,0 for i of [1..V])
|
|
edgeDistsStart := edge.map ([a,b]) => mag tangentPoint xyz[a], xyz[b]
|
|
for each [a,b] of edge
|
|
t := tangentPoint xyz[a], xyz[b]
|
|
scale t, (1 - mag t)/2
|
|
accumulate adj[a], t
|
|
accumulate adj[b], t
|
|
for i of [0...V]
|
|
accumulate xyz[i], scale adj[i], fudge
|
|
edgeDistsEnd := edge.map ([a,b]) => mag tangentPoint xyz[a], xyz[b]
|
|
|
|
function summary(data: number[])
|
|
[Math.min(...data), Math.max(...data), mean(data),
|
|
mean(data.map((x) => Math.abs(1-x)))]
|
|
|
|
function planarize(xyz: XYZ[], face: Face[], fudge: number): void
|
|
V := xyz.length
|
|
adj := (ð 0,0,0 for i of [1..V])
|
|
for each f of face
|
|
if f.length is 3 then continue // triangles always planar
|
|
fxyz := (xyz[v] for each v of f)
|
|
c := centroid fxyz
|
|
n := meanNormal fxyz
|
|
if c dot n < 0 then scale n, -1
|
|
for each v of f
|
|
accumulate adj[v], scale copy(n), n dot (c sub xyz[v])
|
|
for i of [0...V]
|
|
accumulate xyz[i], scale adj[i], fudge
|
|
|
|
function meanNormal(xyz: XYZ[]): XYZ
|
|
mNormal := ð 0,0,0
|
|
[v1, v2] .= xyz.slice(-2);
|
|
for each v3 of xyz
|
|
nextNormal := (v2 sub v1) cross (v3 sub v2)
|
|
magNext := mag nextNormal
|
|
// reduce influence of long edges? (didn't seem to help in brief testing)
|
|
// if magNext > 1 then scale nextNormal, 1/Math.sqrt magNext
|
|
accumulate mNormal, nextNormal
|
|
[v1, v2] = [v2, v3] // shift over one
|
|
scale mNormal, 1/mag mNormal
|
|
|
|
function mean(a: number[])
|
|
m .= 0
|
|
m += e for each e of a
|
|
m / a.length
|
|
|
|
// arithmetic on 3-vectors
|
|
function accumulate(basket: XYZ, egg: XYZ)
|
|
basket[0] += egg[0]
|
|
basket[1] += egg[1]
|
|
basket[2] += egg[2]
|
|
basket
|
|
|
|
function diminish(basket: XYZ, egg: XYZ)
|
|
basket[0] -= egg[0]
|
|
basket[1] -= egg[1]
|
|
basket[2] -= egg[2]
|
|
basket
|
|
|
|
function copy(a: XYZ)
|
|
ð a[0], a[1], a[2]
|
|
|
|
function mag2(a: XYZ)
|
|
a[0]*a[0] + a[1]*a[1] + a[2]*a[2]
|
|
|
|
function mag(a: XYZ)
|
|
Math.sqrt mag2 a
|
|
|
|
function normalize(v: XYZ)
|
|
scale v, 1/mag v
|
|
|
|
function unit
|
|
|
|
function scale(subject: XYZ, by: number)
|
|
subject[0] *= by
|
|
subject[1] *= by
|
|
subject[2] *= by
|
|
subject
|
|
|
|
function lerp(start: XYZ, end: XYZ, howFar: number)
|
|
basket := scale copy(start), 1-howFar
|
|
accumulate basket, scale copy(end), howFar
|
|
|
|
// Feedback
|
|
|
|
function inform(x: string)
|
|
$('input[name="inform"]').val(x)
|
|
x
|
|
|
|
|
|
// VRML97 generation
|
|
|
|
function outputVRML(notation: Notation, P: Polyhedron): string
|
|
shortDescrip := P.name or `A ${P.face.length}-hedron`
|
|
```
|
|
#VRML V2.0 utf8
|
|
Group { children [
|
|
WorldInfo { # Generated by GTW's reimplementation of GWH's Conway script.
|
|
title "${notation} ${stats P}"
|
|
info "Generated by GTW's Conway-notation script inspired by GWH's.
|
|
By using this script, you agree that this image is released
|
|
into the public domain, although you are requested to cite
|
|
George Hart's Encyclopedia of Polyhedra as the source." }
|
|
Background {
|
|
groundColor [ .2 .5 1 ] # light blue
|
|
skyColor [ .2 .5 1 ] }
|
|
NavigationInfo { type [ "EXAMINE" ] }
|
|
DirectionalLight {direction -.5 -1 1 intensity 0.75}
|
|
DirectionalLight {direction .5 1 -1 intensity 0.75}
|
|
Viewpoint { position 0 0 4.5 description "${shortDescrip}" }
|
|
${polyVRML P, colorScheme()}
|
|
Shape {
|
|
appearance Appearance {
|
|
material Material {
|
|
diffuseColor 0 0 0 } }
|
|
geometry IndexedLineSet {
|
|
coord ${useVerts}
|
|
coordIndex [
|
|
${edgeIndices P} ] } }
|
|
${showDual() ? polyVRML geomDual(P), '0.5 0.5 0.5' : ''} ] }
|
|
```
|
|
|
|
function stats(P: Polyhedron): string
|
|
edges := P.face.reduce((e,f) => e + f.length, 0) / 2
|
|
`[${P.face.length} faces, ${edges} edges, ${P.xyz.length} vertices]`
|
|
|
|
useVerts := 'USE verts'
|
|
function defVerts(v: XYZ[])
|
|
`DEF verts Coordinate {
|
|
point [
|
|
${v.map(.join ' ').join(",\n ")} ] }`
|
|
|
|
function polyVRML(P: Polyhedron, color: string): string
|
|
facePartition := color ? [P.face]
|
|
: P.face.reduce ((parts, f) =>
|
|
(parts[f.length] ??= []).push f
|
|
parts),
|
|
[] as Face[][]
|
|
emittedCoords .= ''
|
|
shapes :=
|
|
for part of facePartition
|
|
unless part continue
|
|
`
|
|
Shape {
|
|
appearance Appearance {
|
|
material Material {
|
|
diffuseColor ${color or colorBySides part[0].length} } }
|
|
geometry IndexedFaceSet {
|
|
ccw FALSE
|
|
coord ${emittedCoords ? useVerts : (emittedCoords = defVerts P.xyz)}
|
|
coordIndex [
|
|
${part.map(.join ', ').join(", -1,\n ")}, -1 ] }}`
|
|
shapes.join "\n"
|
|
|
|
function colorScheme
|
|
button := document.getElementsByName('color')[0] as HTMLInputElement
|
|
button.checked ? '1 1 1' : ''
|
|
|
|
function showDual
|
|
false
|
|
|
|
faceColors: Record<number, string> :=
|
|
3: '0.9 0.3 0.3' // red
|
|
4: '0.4 0.4 1.0' // blue
|
|
5: '0.2 0.9 0.3' // green
|
|
6: '0.9 0.9 0.2' // yellow
|
|
7: '0.5 0.25 0.25' // brown
|
|
8: '0.8 0.2 0.8' // magenta
|
|
9: '0.5 0.2 0.8' // purple
|
|
10: '0.1 0.9 0.9' // cyan
|
|
12: '1.0 0.6 0.1' // orange
|
|
|
|
function colorBySides(n: number)
|
|
if n in faceColors
|
|
return faceColors[n]
|
|
'0.5 0.5 0.5' // gray
|
|
|
|
function filtmap<T,U>(A: T[], m: (e:T, i: number, arr: T[]) => U)
|
|
A.map(m).filter (e) => !!e
|
|
|
|
function edgeIndices(P: Polyhedron)
|
|
sep := ",\n "
|
|
filtmap(P.face, (thisf) =>
|
|
filtmap(thisf, (v, ix, f) =>
|
|
preV := f æ ix-1
|
|
preV < v ? `${preV}, ${v}, -1` : '').join sep)
|
|
.join sep
|