The FCC Coordinate System

miller_111.py

"""

a demo of how to convert from cubic to a rhombohedral isometric lattice, which
is equivalent to tetrahedral / octahedral packing, but without the headache
of having to manage two separate primitives or interlaced grid structures
as with the face-centered cubic lattice, which produces an equivalent structure.

an integer conversion from tetrahedral to cartesian coordinates is as easy as

tx = y+z
ty = z+x
tz = x+y

equivalently, the conversion from cartesian to tetrahedral coordinates goes

x = (-tx+ty+tz) / 2
y = ( tx-ty+tz) / 2
z = ( tx+ty-tz) / 2

if unit edge length should be preserved, then the conversions are

tx = (y+z)*D
ty = (z+x)*D
tz = (x+y)*D

and

x = (-tx+ty+tz)*D
y = ( tx-ty+tz)*D
z = ( tx+ty-tz)*D

where D = sqrt(2)/2 = 2^-0.5 = 1/sqrt(2)

the vertex neighborhood of any coordinate consists of 12 points, and
can be reached by the tetrahedral coordinates

(  1  0  0 )
( -1  0  0 )
(  0  1  0 )
(  0 -1  0 )
(  0  0  1 )
(  0  0 -1 )

(  0  1 -1 )
(  0 -1  1 )
( -1  0  1 )
(  1  0 -1 )
(  1 -1  0 )
( -1  1  0 )

which in integer cartesian coordinates are all diagonals of the FCC lattice:

(  0  1  1 )
(  0 -1 -1 )
(  1  0  1 )
( -1  0 -1 )
(  1  1  0 )
( -1 -1  0 )

(  0 -1  1 )
(  0  1 -1 )
(  1  0 -1 )
( -1  0  1 )
( -1  1  0 )
(  1 -1  0 )

about distance metrics: 

the cartesian distance vector between two points A and B in the lattice is

tx = (By-Ay+Bz-Az)*D
ty = (Bz-Az+Bx-Ax)*D
tz = (Bx-Ax+By-Ay)*D

which can also purely be done in tetrahedral metrics

x = Bx-Ax
y = By-Ay
z = Bz-Az

likewise, scaling uniformity is preserved; the lattice is invariant to
translations and scale, but rotations are not preserved; An axis/angle rotation
by Rodrigues' formula nevertheless works if the tetrahedral dot and cross
products are used (see code).

the euclidean norm (length) of a tetrahedral vector is

|v| = sqrt(((y+z)^2 + (z+x)^2 + (x+y)^2)/2)
 
or

|v| = sqrt(x*(x+y+z) + y*(y+z) + z^2)

which saves the final division by 2.

the L1 distance (taxicab, manhattan) is

|v| = |x|+|y|+|z|

the cross product of two tetrahedral vectors in cartesian coordinates is

v = u x v

vx = ( Ux*(Vy - Vz) - Uy*(Vz + Vx) + Uz*(Vx + Vy)) / 2
vy = ( Ux*(Vy + Vz) + Uy*(Vz - Vx) - Uz*(Vx + Vy)) / 2
vz = (-Ux*(Vy + Vz) + Uy*(Vz + Vx) + Uz*(Vx - Vy)) / 2

remaining in tetrahedral coordinates, the cross product is

vx = (- Ux*(  Vy -   Vz) + Uy*(3*Vz +   Vx) - Uz*(  Vx + 3*Vy))*(D/2)
vy = (- Ux*(  Vy + 3*Vz) - Uy*(  Vz -   Vx) + Uz*(3*Vx +   Vy))*(D/2)
vz = (  Ux*(3*Vy +   Vz) - Uy*(  Vz + 3*Vx) - Uz*(  Vx -   Vy))*(D/2)

the dot product of two tetrahedral vectors is

q = u . v

q = ( Ux*(Vz+Vy) + Uy*(Vx+Vz) + Uz*(Vy+Vx) )/2 + Ux*Vx + Uy*Vy + Uz*Vz


"""

from math import *

# these are effectively the same constant, keeping the coordinates at unit 
# length, but due to float bias we need to counteract rounding errors
DIVR2 = 1.0/2**0.5
MULR2 = 2**-0.5

def norm(x,y,z):
    return (x*x+y*y+z*z)**0.5
    
def tetnorm(x,y,z):
    return (x*(x+z+y) + y*(y+z) + z*z)**0.5

def tet2cart(x,y,z):
    return (y+z)*DIVR2,(x+z)*DIVR2,(x+y)*DIVR2

def cart2tet(tx,ty,tz):
    return (-tx+ty+tz)*MULR2,(tx-ty+tz)*MULR2,(tx+ty-tz)*MULR2
    
coords = [
(  1, 0, 0 ),
( -1, 0, 0 ),
(  0, 1, 0 ),
(  0,-1, 0 ),
(  0, 0, 1 ),
(  0, 0,-1 ),

(  0, 1,-1 ),
(  0,-1, 1 ),
( -1, 0, 1 ),
(  1, 0,-1 ),
(  1,-1, 0 ),
( -1, 1, 0 ),
]

for x,y,z in coords:
    print (x,y,z),tetnorm(x,y,z),tet2cart(x,y,z),cart2tet(*tet2cart(x,y,z))

U = (1,2,3)
V = (0,4,1)

def cross(U,V):
    Ux,Uy,Uz = U
    Vx,Vy,Vz = V
    return Uy*Vz - Uz*Vy, Uz*Vx - Ux*Vz, Ux*Vy - Uy*Vx

def dot(U,V):
    Ux,Uy,Uz = U
    Vx,Vy,Vz = V
    return Ux*Vx+Uy*Vy+Uz*Vz

def tetdot(U,V):
    Ux,Uy,Uz = U
    Vx,Vy,Vz = V
    return ( Ux*(Vz+Vy) + Uy*(Vx+Vz) + Uz*(Vy+Vx) )/2.0 + Ux*Vx + Uy*Vy + Uz*Vz

def tetcross2cart(U,V):
    Ux,Uy,Uz = U
    Vx,Vy,Vz = V
    vx = ( Ux*(Vy - Vz) - Uy*(Vz + Vx) + Uz*(Vx + Vy)) / 2.0
    vy = ( Ux*(Vy + Vz) + Uy*(Vz - Vx) - Uz*(Vx + Vy)) / 2.0
    vz = (-Ux*(Vy + Vz) + Uy*(Vz + Vx) + Uz*(Vx - Vy)) / 2.0
    return vx,vy,vz

def tetcross(U,V):
    Ux,Uy,Uz = U
    Vx,Vy,Vz = V
    D = MULR2
    Wx = (- Ux*(  Vy -   Vz) + Uy*(3*Vz +   Vx) - Uz*(  Vx + 3*Vy))*(D/2)
    Wy = (- Ux*(  Vy + 3*Vz) - Uy*(  Vz -   Vx) + Uz*(3*Vx +   Vy))*(D/2)
    Wz = (  Ux*(3*Vy +   Vz) - Uy*(  Vz + 3*Vx) - Uz*(  Vx -   Vy))*(D/2)
    return Wx,Wy,Wz
    
def taxicab(x,y,z):
    return abs(x)+abs(y)+abs(z)

def tetaxisangle(a,U,W):
    s,c = sin(a),cos(a)
    Sx,Sy,Sz = tetcross(W,U)
    t = (1-c)*tetdot(W,U)
    Ux,Uy,Uz = U
    Wx,Wy,Wz = W
    vx = c*Ux + s*Sx + t*Wx
    vy = c*Uy + s*Sy + t*Wy
    vz = c*Uz + s*Sz + t*Wz
    return vx,vy,vz

def axisangle(a,U,W):
    s,c = sin(a),cos(a)
    Sx,Sy,Sz = cross(W,U)
    t = (1-c)*dot(W,U)
    Ux,Uy,Uz = U
    Wx,Wy,Wz = W
    vx = c*Ux + s*Sx + t*Wx
    vy = c*Uy + s*Sy + t*Wy
    vz = c*Uz + s*Sz + t*Wz
    return vx,vy,vz

print cross(tet2cart(*U),tet2cart(*V))
print tetcross2cart(U,V)
print cart2tet(*tetcross2cart(U,V))
print tetcross(U,V)

print dot(tet2cart(*U),tet2cart(*V))

print cart2tet(*axisangle(radians(33.0),tet2cart(*U),tet2cart(0,-1,1)))
print tetaxisangle(radians(33.0),U,(0,-1,1))