Thanks, here's the code:
Vector class:
class Vector(){
public:
func Init(float xx, float yy, float zz)
x = xx
y = yy
z = zz
endfunc
float func Dot(Vector other)
return x*other.x + y*other.y + z*other.z
endfunc
func Divide(float f)
x = x / f
y = y / f
z = z / f
endfunc
func Multiply(Vector a, float f)
x = a.x * f
y = a.y * f
z = a.z * f
endfunc
func MultiplyEquals(float f)
x = x * f
y = y * f
z = z * f
endfunc
func Normalise()
float mag = sqrt(x*x + y * y + z*z)
if mag>0
x = x / mag
y = y / mag
z = z / mag
endif
endfunc
float func Magnitude()
return sqrt(x*x + y*y + z*z)
endfunc
func Add(Vector a, Vector b)
x = a.x + b.x
y = a.y + b.y
z = a.z + b.z
endfunc
func AddEquals(Vector a)
x = x + a.x
y = y + a.y
z = z + a.z
endfunc
func Subtract(Vector a, Vector b)
x = a.x - b.x
y = a.y - b.y
z = a.z - b.z
endfunc
func Cross(Vector a, Vector b)
float _x = (a.y * b.z) - (a.z * b.y)
float _y = (a.z * b.x) - (a.x * b.z)
float _z = (a.x * b.y) - (a.y * b.x)
x = _x
y = _y
z = _z
endfunc
float x
float y
float z
}
Initialisation:
Vector corners[4] ; unit length tetrahedron corners
corners[0] = new Vector()
corners[1] = new Vector()
corners[2] = new Vector()
corners[3] = new Vector()
corners[0].Init(0, 1, 0)
corners[1].Init(0, -1/3, sqrt(8/9))
corners[2].Init(-sqrt(2/3), -1/3, -(1/3) * sqrt(2))
corners[3].Init(sqrt(2/3), -1/3, -(1/3) * sqrt(2))
Vector normals[4]
Vector crosses[4,3] ; crosses[face index][corner index]. cross product of corners for inside-outside measures
Vector C = new Vector()
int triangles[4,3] ; triangles[face index, corner of face]
triangles[0, 0] = 0 ; wish I knew how to do an initialiser list
triangles[0, 1] = 1
triangles[0, 2] = 2
triangles[1, 0] = 0
triangles[1, 1] = 2
triangles[1, 2] = 3
triangles[2, 0] = 0
triangles[2, 1] = 3
triangles[2, 2] = 1
triangles[3, 0] = 2
triangles[3, 1] = 1
triangles[3, 2] = 3
; for use later
Vector point = new Vector()
Vector corn0 = new Vector()
Vector corn1 = new Vector()
Vector corn2 = new Vector()
Vector bary = new Vector()
int it = 0
repeat
normals[it] = new Vector()
crosses[it, 0] = new Vector()
crosses[it, 1] = new Vector()
crosses[it, 2] = new Vector()
Vector cr1 = new Vector()
cr1.Subtract(corners[triangles[it, 2]], corners[triangles[it, 0]])
Vector cr2 = new Vector()
cr2.Subtract(corners[triangles[it, 1]], corners[triangles[it, 0]])
normals[it].Cross(cr1, cr2)
normals[it].Normalise()
; leave this divide out below for an alternative where face of tetrahedron has magnitude 1, rather than corners. Produces slightly different shape
normals[it].Divide(normals[it].Dot(corners[triangles[it, 0]])) ; normal is not unit length so that a corner point gives dot of 1
crosses[it, 0].Cross(corners[triangles[it, 2]], corners[triangles[it, 1]])
crosses[it, 1].Cross(corners[triangles[it, 0]], corners[triangles[it, 2]])
crosses[it, 2].Cross(corners[triangles[it, 1]], corners[triangles[it, 0]])
until (it=it+1)>=4
int squareMap[4,4,3] ; squareMap[quadrant, face index, corner of quadrant]
squareMap[0, 0, 0] = 0
squareMap[0, 0, 1] = 1
squareMap[0, 0, 2] = 2
squareMap[0, 1, 0] = 0
squareMap[0, 1, 1] = 2
squareMap[0, 1, 2] = 3
squareMap[0, 2, 0] = 0
squareMap[0, 2, 1] = 3
squareMap[0, 2, 2] = 1
squareMap[0, 3, 0] = 0
squareMap[0, 3, 1] = 3
squareMap[0, 3, 2] = 1
squareMap[1, 0, 0] = 1
squareMap[1, 0, 1] = 0
squareMap[1, 0, 2] = 3
squareMap[1, 1, 0] = 2
squareMap[1, 1, 1] = 0
squareMap[1, 1, 2] = 1
squareMap[1, 2, 0] = 3
squareMap[1, 2, 1] = 0
squareMap[1, 2, 2] = 2
squareMap[1, 3, 0] = 3
squareMap[1, 3, 1] = 0
squareMap[1, 3, 2] = 2
squareMap[2, 0, 0] = 2
squareMap[2, 0, 1] = 3
squareMap[2, 0, 2] = 0
squareMap[2, 1, 0] = 3
squareMap[2, 1, 1] = 1
squareMap[2, 1, 2] = 0
squareMap[2, 2, 0] = 1
squareMap[2, 2, 1] = 2
squareMap[2, 2, 2] = 0
squareMap[2, 3, 0] = 1
squareMap[2, 3, 1] = 2
squareMap[2, 3, 2] = 0
squareMap[3, 0, 0] = 1
squareMap[3, 0, 1] = 3
squareMap[3, 0, 2] = 2
squareMap[3, 1, 0] = 2
squareMap[3, 1, 1] = 1
squareMap[3, 1, 2] = 3
squareMap[3, 2, 0] = 3
squareMap[3, 2, 1] = 2
squareMap[3, 2, 2] = 1
squareMap[3, 3, 0] = 3
squareMap[3, 3, 1] = 2
squareMap[3, 3, 2] = 1
Inner loop:
; convert c and z to vectors, for convenience
C.x = real(cri)
C.y = imag(cri)
C.z = cj
point.x = real(zri)
point.y = imag(zri)
point.z = zj
float magnitude = 1
int ij = 0 ; ij is just a unique name for the iterator over each tetrahedron face
repeat
magnitude = point.Dot(normals[ij]) ; normals isn't unit length so this gives radial magnitude in shape of tetrahedron
if (magnitude > 0) ; on right side of triangle
; TDL super simple way to get barycentric coordinates of triangle
bary.x = crosses[ij, 0].Dot(point)
bary.y = crosses[ij, 1].Dot(point)
bary.z = crosses[ij, 2].Dot(point)
if (bary.x >= 0 && bary.y >= 0 && bary.z >= 0) ; correct face is found
bary.Divide(bary.x+bary.y+bary.z) ; barycentrics sum to 1
int n = 0
if (bary.x >= 0.5) ; quadrant 0
n = 0
bary.x = bary.x - 0.5
elseif (bary.y >= 0.5) ; quadrant 1
n = 1
bary.y = bary.y - 0.5
elseif (bary.x + bary.y <= 0.5) ; quadrant 2
n = 2
else ; quadrant 3 (middle)
n = 3
float baryz = 1 - (bary.x + bary.y)
bary.y = 0.5 - bary.x
bary.x = 0.5 - baryz
endif
; Z^2 + c
bary.x = bary.x * 2 ; double the size here
bary.y = bary.y * 2
bary.z = 1 - (bary.x + bary.y)
magnitude = magnitude*magnitude;
corn0.Multiply(corners[squareMap[n, ij, 0]], bary.x)
corn1.Multiply(corners[squareMap[n, ij, 1]], bary.y)
corn2.Multiply(corners[squareMap[n, ij, 2]], bary.z)
point.Add(corn0, corn1)
point.AddEquals(corn2)
point.MultiplyEquals(magnitude)
point.AddEquals(C)
ij=10 ; break out of loop since the face was found
endif
endif
until (ij = ij+1)>=4
zri = point.x + flip(point.y) ; convert back to z
zj = point.z