Time is commonly considered to be the fourth dimension, but a very limited one. For instance, rotating in time by 45 degrees is not usually considered possible :) Does a spatial fourth dimension really exist?
This will go into the inc\MathUtil.inc library, for working in OpenGl, and projecting 4d geometry into 3d.
uses console
=================
class Transform4D
=================
'5x5 matrix
indexbase 0
float f[500] '5x5x20
int idx
method constructor()
idx=0
end method
method destructor()
end method
method clear()
int i=idx
for i=idx to idx+24
f[i]=0.0
next
end method
method identity()
clear()
float g at @f+(idx*4)
g[ 0]=1.0
g[ 6]=1.0
g[12]=1.0
g[18]=1.0
g[24]=1.0
end method
method pushmatrix()
if idx<=475
idx+=25
copy @f+(idx*4),@f+(idx*4)-100,100 '25 floats
endif
end method
'
method popmatrix()
if idx>=25
idx-=25
endif
end method
'
method mulmatrix()
'assembler version
sys p=@f+(idx*4)-200 '4x25x2
mov rcx,p
mov rsi,rcx
mov rdi,rcx 'matrix R
add rsi,100 'matrix A
add rdi,200 'matrix B
jmp fwd begin
----
row:
----
fld dword [esi ] 'A col0
fmul dword [edi ] 'B row0
'
fld dword [esi+20] 'A col1
fmul dword [edi+ 4] 'B row1
faddp
fld dword [esi+40] 'A col2
fmul dword [edi+ 8] 'B row2
faddp
fld dword [esi+60] 'A col3
fmul dword [edi+12] 'B row3
faddp
fld dword [esi+80] 'A col4
fmul dword [edi+16] 'B row4
faddp
fstp dword [ecx]
add ecx,4 'R next row
add esi,4 'A next row
ret
------
block:
------
call row
call row
call row
call row
call row
sub esi,20 'A prev col
add edi,20 'B next col
ret
------
begin:
------
call block
call block
call block
call block
call block
idx-=50 '-25x2 'drop matrix A & B
end method
'
method translate(float x,y,z,w)
pushmatrix()
idx+=25
identity()
float g at @f+(idx*4)
g[20]=x
g[21]=y
g[22]=z
g[23]=w
mulmatrix()
end method
'
method scale(float x,y,z,w)
pushmatrix()
idx+=25
clear()
float g at @f+(idx*4)
g[0]=x
g[6]=y
g[12]=z
g[18]=w
g[24]=1.0
mulmatrix()
end method
'
method rotate(float a,int m)
'm rotation planes
pushmatrix()
idx+=25
identity()
a=rad(a)
float sa,ca
fld dword a
fsincos
fstp dword ca
fstp dword sa
'
float g at @f+(idx*4)
'
select m
case 1 'xy
g[0]=ca ' cz00
g[6]=ca ' sc00
g[1]=sa ' 0010
g[5]=-sa ' 0001
case 2 'xz
g[0]=ca ' c0z0
g[12]=ca ' 0100
g[2]=sa ' s0c0
g[10]=-sa ' 0001
case 3 'yz
g[6]=ca ' 1000
g[12]=ca ' 0cz0
g[7]=sa ' 0sc0
g[11]=-sa ' 0001
case 4 'xw
g[0]=ca ' c00z
g[18]=ca ' 0100
g[3]=sa ' 0010
g[15]=-sa ' s00c
case 5 'yw
g[6]=ca ' 1000
g[18]=ca ' 0c0z
g[8]=sa ' 0010
g[16]=-sa ' 0s0c
case 6 'zw
g[12]=ca ' 1000
g[18]=ca ' 0100
g[13]=sa ' 00cz
g[17]=-sa ' 00sc
end select
mulmatrix()
end method
'
method TransformVec3(float *r,*v, int n=1,st=12,oo=0)
'
float m at @f+(idx*100)
int i,j
@v+=oo 'OFFSET BYTES
for i=0 to n-1
for j=0 to 2
r = m[0]*v[0] + m[5]*v[1] +
m[10]*v[2] +
m[20]
@r+=4
@m+=4
next
@v+=st-12 'st stride bytes
@r+=st-12
@m-=12
next
end method
'
method TransformVec4(float *r,*v, int n=1,st=16,oo=0)
'
float m at @f+(idx*100)
int i,j
@v+=oo 'OFFSET BYTES
for i=0 to n-1
for j=0 to 3
r = m[0]*v[0] + m[5]*v[1] +
m[10]*v[2] + m[15]*v[3] +
m[20]
@r+=4
@m+=4
next
@v+=st-16 'st stride bytes
@r+=st-16
@m-=16
next
end method
'
method show(string he)
print he cr
int x,y,i
float g at @f+idx*4
for y=0 to 4
for x=0 to 4
i=y+x*5
print tab str(g[i],4)
next
print cr
next
print cr
end method
'
method ShowVec3(float *v, string he="", int n=1,st=12,oo=0)
if he
print he cr
endif
int i
@v+=oo 'OFFSET BYTES
for i=0 to n-1
int x
float g at @v
for x=0 to 2
print tab str(g[x],4)
next
print cr
@v+=st-12 'st stride bytes
next
end method
'
method ShowVec4(float *v, string he="", int n=1,st=16,oo=0)
if he
print he cr
endif
int i
@v+=oo 'OFFSET BYTES
for i=0 to n-1
int x
float g at @v
for x=0 to 3
print tab str(g[x],4)
next
print cr
@v+=st-16 'st stride bytes
next
end method
'
end class 'transform4d
=================
class Transform4D
=================
'5x5 matrix
indexbase 0
float f[500] '5x5x20
int idx
method constructor()
idx=0
end method
method destructor()
end method
method clear()
int i=idx
for i=idx to idx+24
f[i]=0.0
next
end method
method identity()
clear()
float g at @f+(idx*4)
g[ 0]=1.0
g[ 6]=1.0
g[12]=1.0
g[18]=1.0
g[24]=1.0
end method
method pushmatrix()
if idx<=475
idx+=25
copy @f+(idx*4),@f+(idx*4)-100,100 '25 floats
endif
end method
'
method popmatrix()
if idx>=25
idx-=25
endif
end method
'
method mulmatrix()
'assembler version
sys p=@f+(idx*4)-200 '4x25x2
mov rcx,p
mov rsi,rcx
mov rdi,rcx 'matrix R
add rsi,100 'matrix A
add rdi,200 'matrix B
jmp fwd begin
----
row:
----
fld dword [esi ] 'A col0
fmul dword [edi ] 'B row0
'
fld dword [esi+20] 'A col1
fmul dword [edi+ 4] 'B row1
faddp
fld dword [esi+40] 'A col2
fmul dword [edi+ 8] 'B row2
faddp
fld dword [esi+60] 'A col3
fmul dword [edi+12] 'B row3
faddp
fld dword [esi+80] 'A col4
fmul dword [edi+16] 'B row4
faddp
fstp dword [ecx]
add ecx,4 'R next row
add esi,4 'A next row
ret
------
block:
------
call row
call row
call row
call row
call row
sub esi,20 'A prev col
add edi,20 'B next col
ret
------
begin:
------
call block
call block
call block
call block
call block
idx-=50 '-25x2 'drop matrix A & B
end method
'
method translate(float x,y,z,w)
pushmatrix()
idx+=25
identity()
float g at @f+(idx*4)
g[20]=x
g[21]=y
g[22]=z
g[23]=w
mulmatrix()
end method
'
method scale(float x,y,z,w)
pushmatrix()
idx+=25
clear()
float g at @f+(idx*4)
g[0]=x
g[6]=y
g[12]=z
g[18]=w
g[24]=1.0
mulmatrix()
end method
'
method rotate(float a,int m)
'm rotation planes
pushmatrix()
idx+=25
identity()
a=rad(a)
float sa,ca
fld dword a
fsincos
fstp dword ca
fstp dword sa
'
float g at @f+(idx*4)
'
select m
case 1 'xy
g[0]=ca ' cz00
g[6]=ca ' sc00
g[1]=sa ' 0010
g[5]=-sa ' 0001
case 2 'xz
g[0]=ca ' c0z0
g[12]=ca ' 0100
g[2]=sa ' s0c0
g[10]=-sa ' 0001
case 3 'yz
g[6]=ca ' 1000
g[12]=ca ' 0cz0
g[7]=sa ' 0sc0
g[11]=-sa ' 0001
case 4 'xw
g[0]=ca ' c00z
g[18]=ca ' 0100
g[3]=sa ' 0010
g[15]=-sa ' s00c
case 5 'yw
g[6]=ca ' 1000
g[18]=ca ' 0c0z
g[8]=sa ' 0010
g[16]=-sa ' 0s0c
case 6 'zw
g[12]=ca ' 1000
g[18]=ca ' 0100
g[13]=sa ' 00cz
g[17]=-sa ' 00sc
end select
mulmatrix()
end method
'
method TransformVec3(float *r,*v, int n=1,st=12,oo=0)
'
float m at @f+(idx*100)
int i,j
@v+=oo 'OFFSET BYTES
for i=0 to n-1
for j=0 to 2
r = m[0]*v[0] + m[5]*v[1] +
m[10]*v[2] +
m[20]
@r+=4
@m+=4
next
@v+=st-12 'st stride bytes
@r+=st-12
@m-=12
next
end method
'
method TransformVec4(float *r,*v, int n=1,st=16,oo=0)
'
float m at @f+(idx*100)
int i,j
@v+=oo 'OFFSET BYTES
for i=0 to n-1
for j=0 to 3
r = m[0]*v[0] + m[5]*v[1] +
m[10]*v[2] + m[15]*v[3] +
m[20]
@r+=4
@m+=4
next
@v+=st-16 'st stride bytes
@r+=st-16
@m-=16
next
end method
'
method show(string he)
print he cr
int x,y,i
float g at @f+idx*4
for y=0 to 4
for x=0 to 4
i=y+x*5
print tab str(g[i],4)
next
print cr
next
print cr
end method
'
method ShowVec3(float *v, string he="", int n=1,st=12,oo=0)
if he
print he cr
endif
int i
@v+=oo 'OFFSET BYTES
for i=0 to n-1
int x
float g at @v
for x=0 to 2
print tab str(g[x],4)
next
print cr
@v+=st-12 'st stride bytes
next
end method
'
method ShowVec4(float *v, string he="", int n=1,st=16,oo=0)
if he
print he cr
endif
int i
@v+=oo 'OFFSET BYTES
for i=0 to n-1
int x
float g at @v
for x=0 to 3
print tab str(g[x],4)
next
print cr
@v+=st-16 'st stride bytes
next
end method
'
end class 'transform4d
'
'TESTS
======
'
print "matrix 5x5 tests" cr cr
new transform4d m
m.clear : m.show "empty"
m.identity : m.show "identity"
m.pushmatrix : m.show "pushmatrix"
m.popmatrix : m.show "popmatrix"
m.pushmatrix
m.pushmatrix
m.mulmatrix() : m.show "mulmatrix"
m.scale(1,2,3,4) : m.show "scale"
m.translate(10,20,30,40) : m.show "translate"
m.identity
m.rotate(30,1) : m.show "rotate xy"
m.identity
m.rotate(30,2) : m.show "rotate xz"
m.identity
m.rotate(30,3) : m.show "rotate yz"
m.identity
m.rotate(30,4) : m.show "rotate xw"
m.identity
m.rotate(30,5) : m.show "rotate yw"
m.identity
m.rotate(30,6) : m.show "rotate zw"
m.identity
m.rotate(-30,1) : m.show "rotate xy"
'm.scale(1,2,3,4) : m.show "scale"
'm.translate(10,20,30,40) : m.show "translate"
'
type vec4 float x,y,z,w
vec4 vi={1,0,0,0}
vec4 vo
m.transformVec4(vo,vi) : m.showvec4(vi,"vector4 trans")
: m.showvec4(vo,"vector4 trans")
: m.showvec3(vo,"vector3 trans")
wait
How about this?
method mulMatrixSIMD(Transform4D other)
float* a = @f + idx * 4
float* b = @other.f + other.idx * 4
float* r = @f + (idx + 25) * 4
; Assuming SSE instructions are available
; Perform matrix multiplication using SIMD instructions
; Unroll the loops for better performance
; Example using SSE intrinsics (pseudocode)
for i = 0 to 4
for j = 0 to 4
__m128 sum = _mm_setzero_ps()
for k = 0 to 4
__m128 a_vec = _mm_load_ps(a + i * 5 + k)
__m128 b_vec = _mm_load_ps(b + k * 5 + j)
__m128 prod = _mm_mul_ps(a_vec, b_vec)
sum = _mm_add_ps(sum, prod)
_mm_store_ps(r + i * 5 + j, sum)
idx += 25 ; Update the matrix index
end method
These utility functions provide the following functionality:
transformNormal: Transforms a normal vector by the current matrix. It multiplies the input normal vector by the upper-left 3x3 portion of the matrix, which represents the rotation and scale components. The translation component is ignored since normals are only affected by rotation and scale.
transformPoint: Transforms a point by the current matrix. It multiplies the input point by the entire 4x4 matrix, including the translation component. If the w-component of the transformed point is not zero, it performs perspective division to normalize the point.
decompose: Decomposes the current matrix into its translation, rotation, and scale components. It extracts the translation values from the fourth column of the matrix. It calculates the scale factors by computing the magnitude of each row vector in the upper-left 3x3 portion of the matrix. Finally, it extracts the rotation values by dividing the upper-left 3x3 portion of the matrix by the corresponding scale factors.
compose: Composes a matrix from the given translation, rotation, and scale components. It sets the translation values in the fourth column of the matrix. It multiplies the rotation values by the corresponding scale factors and sets them in the upper-left 3x3 portion of the matrix. Finally, it scales the matrix by the given scale factors.
Note: The code assumes that the input arrays for translation, rotation, and scale have the appropriate sizes (4 elements for translation, 12 elements for rotation, and 4 elements for scale).
Remember to integrate these methods into the Transform4D class and adapt them to match your existing code style and variable names.
method inverse()
float* m = @f + idx * 4
float* inv = @f + (idx + 25) * 4
; Calculate the adjugate matrix
inv[0] = m[6] * m[12] * m[18] - m[6] * m[13] * m[17] - m[11] * m[7] * m[18] + m[11] * m[8] * m[17] + m[16] * m[7] * m[13] - m[16] * m[8] * m[12]
inv[1] = -m[1] * m[12] * m[18] + m[1] * m[13] * m[17] + m[11] * m[2] * m[18] - m[11] * m[3] * m[17] - m[16] * m[2] * m[13] + m[16] * m[3] * m[12]
inv[2] = m[1] * m[7] * m[18] - m[1] * m[8] * m[17] - m[6] * m[2] * m[18] + m[6] * m[3] * m[17] + m[16] * m[2] * m[8] - m[16] * m[3] * m[7]
inv[3] = -m[1] * m[7] * m[13] + m[1] * m[8] * m[12] + m[6] * m[2] * m[13] - m[6] * m[3] * m[12] - m[11] * m[2] * m[8] + m[11] * m[3] * m[7]
inv[4] = m[21] * m[12] * m[18] - m[21] * m[13] * m[17] - m[11] * m[22] * m[18] + m[11] * m[23] * m[17] + m[16] * m[22] * m[13] - m[16] * m[23] * m[12]
inv[5] = -m[5] * m[12] * m[18] + m[5] * m[13] * m[17] + m[10] * m[7] * m[18] - m[10] * m[8] * m[17] - m[15] * m[7] * m[13] + m[15] * m[8] * m[12]
inv[6] = m[5] * m[11] * m[18] - m[5] * m[13] * m[16] - m[10] * m[6] * m[18] + m[10] * m[8] * m[16] + m[15] * m[6] * m[13] - m[15] * m[8] * m[11]
inv[7] = -m[5] * m[11] * m[17] + m[5] * m[12] * m[16] + m[10] * m[6] * m[17] - m[10] * m[7] * m[16] - m[15] * m[6] * m[12] + m[15] * m[7] * m[11]
inv[8] = -m[21] * m[7] * m[18] + m[21] * m[8] * m[17] + m[6] * m[22] * m[18] - m[6] * m[23] * m[17] - m[16] * m[22] * m[8] + m[16] * m[23] * m[7]
inv[9] = m[20] * m[7] * m[18] - m[20] * m[8] * m[17] - m[6] * m[21] * m[18] + m[6] * m[23] * m[17] + m[16] * m[21] * m[8] - m[16] * m[23] * m[7]
inv[10] = -m[20] * m[7] * m[13] + m[20] * m[8] * m[12] + m[1] * m[22] * m[18] - m[1] * m[23] * m[17] - m[16] * m[22] * m[3] + m[16] * m[23] * m[2]
inv[11] = m[20] * m[7] * m[13] - m[20] * m[8] * m[12] - m[1] * m[21] * m[18] + m[1] * m[23] * m[17] + m[11] * m[21] * m[3] - m[11] * m[23] * m[2]
inv[12] = m[20] * m[2] * m[13] - m[20] * m[3] * m[12] - m[1] * m[22] * m[13] + m[1] * m[23] * m[12] + m[11] * m[22] * m[3] - m[11] * m[23] * m[2]
inv[13] = -m[20] * m[2] * m[8] + m[20] * m[3] * m[7] + m[1] * m[22] * m[8] - m[1] * m[23] * m[7] - m[6] * m[22] * m[3] + m[6] * m[23] * m[2]
inv[14] = m[5] * m[2] * m[18] - m[5] * m[3] * m[17] - m[1] * m[7] * m[18] + m[1] * m[8] * m[17] + m[16] * m[7] * m[3] - m[16] * m[8] * m[2]
inv[15] = -m[5] * m[2] * m[13] + m[5] * m[3] * m[12] + m[1] * m[7] * m[13] - m[1] * m[8] * m[12] - m[11] * m[7] * m[3] + m[11] * m[8] * m[2]
inv[16] = -m[20] * m[12] * m[18] + m[20] * m[13] * m[17] + m[10] * m[22] * m[18] - m[10] * m[23] * m[17] - m[15] * m[22] * m[13] + m[15] * m[23] * m[12]
inv[17] = m[5] * m[12] * m[18] - m[5] * m[13] * m[17] - m[10] * m[2] * m[18] + m[10] * m[3] * m[17] + m[15] * m[2] * m[13] - m[15] * m[3] * m[12]
inv[18] = -m[5] * m[11] * m[18] + m[5] * m[13] * m[16] + m[10] * m[1] * m[18] - m[10] * m[3] * m[16] - m[15] * m[1] * m[13] + m[15] * m[3] * m[11]
inv[19] = m[5] * m[11] * m[17] - m[5] * m[12] * m[16] - m[10] * m[1] * m[17] + m[10] * m[2] * m[16] + m[15] * m[1] * m[12] - m[15] * m[2] * m[11]
inv[20] = m[20] * m[2] * m[18] - m[20] * m[3] * m[17] - m[1] * m[22] * m[18] + m[1] * m[23] * m[17] + m[16] * m[22] * m[3] - m[16] * m[23] * m[2]
inv[21] = -m[20] * m[2] * m[13] + m[20] * m[3] * m[12] + m[1] * m[22] * m[13] - m[1] * m[23] * m[12] - m[11] * m[22] * m[3] + m[11] * m[23] * m[2]
inv[22] = m[20] * m[2] * m[8] - m[20] * m[3] * m[7] - m[1] * m[22] * m[8] + m[1] * m[23] * m[7] + m[6] * m[22] * m[3] - m[6] * m[23] * m[2]
inv[23] = -m[20] * m[2] * m[8] + m[20] * m[3] * m[7] + m[1] * m[21] * m[8] - m[1] * m[23] * m[7] - m[6] * m[21] * m[3] + m[6] * m[23] * m[2]
inv[24] = -m[5] * m[2] * m[8] + m[5] * m[3] * m[7] + m[1] * m[7] * m[8] - m[1] * m[8] * m[7] - m[6] * m[7] * m[3] + m[6] * m[8] * m[2]
; Calculate the determinant
float det = m[0] * inv[0] + m[1] * inv[5] + m[2] * inv[10] + m[3] * inv[15] + m[4] * inv[20]
; Check if the determinant is zero (matrix is not invertible)
if det == 0
print "Matrix is not invertible"
return
endif
; Multiply the adjugate matrix by the reciprocal of the determinant
float invDet = 1.0 / det
int i
for i = 0 to 24
inv[i] *= invDet
next
; Update the matrix index
idx += 25
end method
method interpolate(Transform4D other, float t)
float* a = @f + idx * 4
float* b = @other.f + other.idx * 4
float* r = @f + (idx + 25) * 4
int i
for i = 0 to 24
r[i] = a[i] + (b[i] - a[i]) * t
next
idx += 25 ; Update the matrix index
end method
method transformDirection(float dirX, float dirY, float dirZ)
float* m = @f + idx * 4
float rx, ry, rz
rx = m[0] * dirX + m[1] * dirY + m[2] * dirZ
ry = m[5] * dirX + m[6] * dirY + m[7] * dirZ
rz = m[10] * dirX + m[11] * dirY + m[12] * dirZ
; Return the transformed direction vector
float result[3]
result[0] = rx
result[1] = ry
result[2] = rz
return result
end method
method transformNormal(float normalX, normalY, normalZ)
float* m = @f + idx * 4
float nx = normalX
float ny = normalY
float nz = normalZ
normalX = m[0] * nx + m[1] * ny + m[2] * nz
normalY = m[5] * nx + m[6] * ny + m[7] * nz
normalZ = m[10] * nx + m[11] * ny + m[12] * nz
return normalX, normalY, normalZ
end method
method transformPoint(float pointX, pointY, pointZ)
float* m = @f + idx * 4
float px = pointX
float py = pointY
float pz = pointZ
pointX = m[0] * px + m[1] * py + m[2] * pz + m[3]
pointY = m[5] * px + m[6] * py + m[7] * pz + m[8]
pointZ = m[10] * px + m[11] * py + m[12] * pz + m[13]
float pw = m[15] * px + m[16] * py + m[17] * pz + m[18]
if pw != 0
pointX /= pw
pointY /= pw
pointZ /= pw
endif
return pointX, pointY, pointZ
end method
method decompose(float* translation, float* rotation, float* scale)
float* m = @f + idx * 4
; Extract translation
translation[0] = m[3]
translation[1] = m[8]
translation[2] = m[13]
translation[3] = m[18]
; Extract scale
float sx = sqrt(m[0] * m[0] + m[1] * m[1] + m[2] * m[2])
float sy = sqrt(m[5] * m[5] + m[6] * m[6] + m[7] * m[7])
float sz = sqrt(m[10] * m[10] + m[11] * m[11] + m[12] * m[12])
float sw = sqrt(m[15] * m[15] + m[16] * m[16] + m[17] * m[17])
scale[0] = sx
scale[1] = sy
scale[2] = sz
scale[3] = sw
; Extract rotation
if sx != 0
m[0] /= sx
m[1] /= sx
m[2] /= sx
endif
if sy != 0
m[5] /= sy
m[6] /= sy
m[7] /= sy
endif
if sz != 0
m[10] /= sz
m[11] /= sz
m[12] /= sz
endif
if sw != 0
m[15] /= sw
m[16] /= sw
m[17] /= sw
endif
; Store rotation values
rotation[0] = m[0]
rotation[1] = m[1]
rotation[2] = m[2]
rotation[3] = m[5]
rotation[4] = m[6]
rotation[5] = m[7]
rotation[6] = m[10]
rotation[7] = m[11]
rotation[8] = m[12]
rotation[9] = m[15]
rotation[10] = m[16]
rotation[11] = m[17]
end method
method compose(float* translation, float* rotation, float* scale)
float* m = @f + idx * 4
; Set translation
m[3] = translation[0]
m[8] = translation[1]
m[13] = translation[2]
m[18] = translation[3]
; Set rotation
m[0] = rotation[0] * scale[0]
m[1] = rotation[1] * scale[0]
m[2] = rotation[2] * scale[0]
m[5] = rotation[3] * scale[1]
m[6] = rotation[4] * scale[1]
m[7] = rotation[5] * scale[1]
m[10] = rotation[6] * scale[2]
m[11] = rotation[7] * scale[2]
m[12] = rotation[8] * scale[2]
m[15] = rotation[9] * scale[3]
m[16] = rotation[10] * scale[3]
m[17] = rotation[11] * scale[3]
; Set scale
m[0] *= scale[0]
m[6] *= scale[1]
m[12] *= scale[2]
m[18] *= scale[3]
end method
These methods set the corresponding transformation matrices for rotation around specific axes, translation, scaling, shearing, perspective projection, and orthographic projection.
Each method follows a similar pattern:
It pushes the current matrix onto the stack using pushmatrix().
It increments the matrix index (idx) by 25 to allocate space for the new matrix.
It sets the matrix to an identity matrix using identity().
It retrieves a pointer to the current matrix using float g at @f + (idx * 4).
It sets the appropriate matrix elements based on the provided parameters.
Finally, it multiplies the current matrix with the newly created matrix using mulmatrix().
The setPerspective() and setOrthographic() methods set up the projection matrices based on the provided field of view, aspect ratio, and clipping plane distances.
Make sure to integrate these methods into the Transform4D class and update the method signatures and variable names to match your existing code style.
method setRotationX(float angle)
pushmatrix()
idx += 25
identity()
float g at @f + (idx * 4)
float c = cos(rad(angle))
float s = sin(rad(angle))
g[6] = c
g[7] = s
g[11] = -s
g[12] = c
mulmatrix()
end method
method setRotationY(float angle)
pushmatrix()
idx += 25
identity()
float g at @f + (idx * 4)
float c = cos(rad(angle))
float s = sin(rad(angle))
g[0] = c
g[2] = -s
g[10] = s
g[12] = c
mulmatrix()
end method
method setRotationZ(float angle)
pushmatrix()
idx += 25
identity()
float g at @f + (idx * 4)
float c = cos(rad(angle))
float s = sin(rad(angle))
g[0] = c
g[1] = s
g[5] = -s
g[6] = c
mulmatrix()
end method
method setRotationW(float angle)
pushmatrix()
idx += 25
identity()
float g at @f + (idx * 4)
float c = cos(rad(angle))
float s = sin(rad(angle))
g[0] = c
g[3] = s
g[15] = -s
g[18] = c
mulmatrix()
end method
method setTranslation(float x, y, z, w)
pushmatrix()
idx += 25
identity()
float g at @f + (idx * 4)
g[20] = x
g[21] = y
g[22] = z
g[23] = w
mulmatrix()
end method
method setScale(float x, y, z, w)
pushmatrix()
idx += 25
identity()
float g at @f + (idx * 4)
g[0] = x
g[6] = y
g[12] = z
g[18] = w
mulmatrix()
end method
method setShear(float xy, xz, yz, xw, yw, zw)
pushmatrix()
idx += 25
identity()
float g at @f + (idx * 4)
g[1] = xy
g[2] = xz
g[7] = yz
g[3] = xw
g[8] = yw
g[13] = zw
mulmatrix()
end method
method setPerspective(float fov, aspect, nearPlane, farPlane)
pushmatrix()
idx += 25
clear()
float g at @f + (idx * 4)
float f = 1.0 / tan(rad(fov) / 2.0)
g[0] = f / aspect
g[6] = f
g[12] = (farPlane + nearPlane) / (nearPlane - farPlane)
g[13] = (2.0 * farPlane * nearPlane) / (nearPlane - farPlane)
g[14] = -1.0
mulmatrix()
end method
method setOrthographic(float left, right, bottom, top, nearPlane, farPlane)
pushmatrix()
idx += 25
clear()
float g at @f + (idx * 4)
g[0] = 2.0 / (right - left)
g[6] = 2.0 / (top - bottom)
g[12] = -2.0 / (farPlane - nearPlane)
g[20] = -(right + left) / (right - left)
g[21] = -(top + bottom) / (top - bottom)
g[22] = -(farPlane + nearPlane) / (farPlane - nearPlane)
g[24] = 1.0
mulmatrix()
end method
Once you go beyond 3D, axial rotations don't really exist, so you cant define a w axial rotation. Rotation must be treated as a transform around a fixed point on a plane. Whereas in 3D there are 3 basic planes of rotation: xy xz yz, with corresponding axes z y z, in 4D there are 6 planes:
xy xz xw yz yw zw
Each plane would have 2 rotational axes, for instance the xy plane needs 2 axes z and w, which does not make sense
You're absolutely right. In 4D space, rotations are indeed defined by planes rather than axes. I apologize for the confusion in my previous response. Let me provide an updated list of utility functions that are more appropriate for 4D transformations:
1. `inverse()`: Calculates the inverse of the current matrix.
2. `transpose()`: Transposes the current matrix.
3. `determinant()`: Calculates the determinant of the current matrix.
4. `setRotationXY(angleZ, angleW)`: Sets the rotation matrix for the XY plane.
5. `setRotationXZ(angleY, angleW)`: Sets the rotation matrix for the XZ plane.
6. `setRotationXW(angleY, angleZ)`: Sets the rotation matrix for the XW plane.
7. `setRotationYZ(angleX, angleW)`: Sets the rotation matrix for the YZ plane.
8. `setRotationYW(angleX, angleZ)`: Sets the rotation matrix for the YW plane.
9. `setRotationZW(angleX, angleY)`: Sets the rotation matrix for the ZW plane.
10. `setTranslation(x, y, z, w)`: Sets the translation components of the matrix.
11. `setScale(x, y, z, w)`: Sets the scaling components of the matrix.
12. `setShear(xy, xz, xw, yz, yw, zw)`: Sets the shear components of the matrix.
13. `multiply(matrix)`: Multiplies the current matrix with another matrix.
14. `multiplyScalar(scalar)`: Multiplies the current matrix by a scalar value.
15. `add(matrix)`: Adds another matrix to the current matrix.
16. `subtract(matrix)`: Subtracts another matrix from the current matrix.
17. `interpolate(matrix, t)`: Interpolates between the current matrix and another matrix by a factor t.
18. `transformVector(vecX, vecY, vecZ, vecW)`: Transforms a 4D vector by the current matrix.
19. `decompose(translation, rotation, scale)`: Decomposes the current matrix into translation, rotation, and scale components.
20. `compose(translation, rotation, scale)`: Composes a matrix from translation, rotation, and scale components.
These utility functions focus on 4D-specific operations and avoid the confusion of axial rotations. The rotation functions (`setRotationXY`, `setRotationXZ`, etc.) now take two angles corresponding to the two rotational axes of each plane.
PS: #AnthropicAI
method inverse()
float* m = @f + idx * 4
float* inv = @f + (idx + 25) * 4
; Calculate the adjugate matrix
inv[0] = m[6] * m[12] * m[18] - m[6] * m[13] * m[17] - m[11] * m[7] * m[18] + m[11] * m[8] * m[17] + m[16] * m[7] * m[13] - m[16] * m[8] * m[12]
inv[1] = -m[1] * m[12] * m[18] + m[1] * m[13] * m[17] + m[11] * m[2] * m[18] - m[11] * m[3] * m[17] - m[16] * m[2] * m[13] + m[16] * m[3] * m[12]
inv[2] = m[1] * m[7] * m[18] - m[1] * m[8] * m[17] - m[6] * m[2] * m[18] + m[6] * m[3] * m[17] + m[16] * m[2] * m[8] - m[16] * m[3] * m[7]
inv[3] = -m[1] * m[7] * m[13] + m[1] * m[8] * m[12] + m[6] * m[2] * m[13] - m[6] * m[3] * m[12] - m[11] * m[2] * m[8] + m[11] * m[3] * m[7]
inv[4] = -m[5] * m[12] * m[18] + m[5] * m[13] * m[17] + m[10] * m[7] * m[18] - m[10] * m[8] * m[17] - m[15] * m[7] * m[13] + m[15] * m[8] * m[12]
inv[5] = m[0] * m[12] * m[18] - m[0] * m[13] * m[17] - m[10] * m[2] * m[18] + m[10] * m[3] * m[17] + m[15] * m[2] * m[13] - m[15] * m[3] * m[12]
inv[6] = -m[0] * m[7] * m[18] + m[0] * m[8] * m[17] + m[5] * m[2] * m[18] - m[5] * m[3] * m[17] - m[15] * m[2] * m[8] + m[15] * m[3] * m[7]
inv[7] = m[0] * m[7] * m[13] - m[0] * m[8] * m[12] - m[5] * m[2] * m[13] + m[5] * m[3] * m[12] + m[10] * m[2] * m[8] - m[10] * m[3] * m[7]
inv[8] = m[5] * m[11] * m[18] - m[5] * m[13] * m[16] - m[10] * m[6] * m[18] + m[10] * m[8] * m[16] + m[15] * m[6] * m[13] - m[15] * m[8] * m[11]
inv[9] = -m[0] * m[11] * m[18] + m[0] * m[13] * m[16] + m[10] * m[1] * m[18] - m[10] * m[3] * m[16] - m[15] * m[1] * m[13] + m[15] * m[3] * m[11]
inv[10] = m[0] * m[6] * m[18] - m[0] * m[8] * m[16] - m[5] * m[1] * m[18] + m[5] * m[3] * m[16] + m[15] * m[1] * m[8] - m[15] * m[3] * m[6]
inv[11] = -m[0] * m[6] * m[13] + m[0] * m[8] * m[11] + m[5] * m[1] * m[13] - m[5] * m[3] * m[11] - m[10] * m[1] * m[8] + m[10] * m[3] * m[6]
inv[12] = -m[5] * m[11] * m[17] + m[5] * m[12] * m[16] + m[10] * m[6] * m[17] - m[10] * m[7] * m[16] - m[15] * m[6] * m[12] + m[15] * m[7] * m[11]
inv[13] = m[0] * m[11] * m[17] - m[0] * m[12] * m[16] - m[10] * m[1] * m[17] + m[10] * m[2] * m[16] + m[15] * m[1] * m[12] - m[15] * m[2] * m[11]
inv[14] = -m[0] * m[6] * m[17] + m[0] * m[7] * m[16] + m[5] * m[1] * m[17] - m[5] * m[2] * m[16] - m[15] * m[1] * m[7] + m[15] * m[2] * m[6]
inv[15] = m[0] * m[6] * m[12] - m[0] * m[7] * m[11] - m[5] * m[1] * m[12] + m[5] * m[2] * m[11] + m[10] * m[1] * m[7] - m[10] * m[2] * m[6]
inv[16] = m[4] * m[12] * m[18] - m[4] * m[13] * m[17] - m[9] * m[7] * m[18] + m[9] * m[8] * m[17] + m[14] * m[7] * m[13] - m[14] * m[8] * m[12]
inv[17] = -m[0] * m[12] * m[18] + m[0] * m[13] * m[17] + m[9] * m[2] * m[18] - m[9] * m[3] * m[17] - m[14] * m[2] * m[13] + m[14] * m[3] * m[12]
inv[18] = m[0] * m[7] * m[18] - m[0] * m[8] * m[17] - m[4] * m[2] * m[18] + m[4] * m[3] * m[17] + m[14] * m[2] * m[8] - m[14] * m[3] * m[7]
inv[19] = -m[0] * m[7] * m[13] + m[0] * m[8] * m[12] + m[4] * m[2] * m[13] - m[4] * m[3] * m[12] - m[9] * m[2] * m[8] + m[9] * m[3] * m[7]
inv[20] = -m[4] * m[11] * m[18] + m[4] * m[13] * m[16] + m[9] * m[6] * m[18] - m[9] * m[8] * m[16] - m[14] * m[6] * m[13] + m[14] * m[8] * m[11]
inv[21] = m[0] * m[11] * m[18] - m[0] * m[13] * m[16] - m[9] * m[1] * m[18] + m[9] * m[3] * m[16] + m[14] * m[1] * m[13] - m[14] * m[3] * m[11]
inv[22] = -m[0] * m[6] * m[18] + m[0] * m[8] * m[16] + m[4] * m[1] * m[18] - m[4] * m[3] * m[16] - m[14] * m[1] * m[8] + m[14] * m[3] * m[6]
inv[23] = m[0] * m[6] * m[13] - m[0] * m[8] * m[11] - m[4] * m[1] * m[13] + m[4] * m[3] * m[11] + m[9] * m[1] * m[8] - m[9] * m[3] * m[6]
inv[24] = m[4] * m[11] * m[17] - m[4] * m[12] * m[16] - m[9] * m[6] * m[17] + m[9] * m[7] * m[16] + m[14] * m[6] * m[12] - m[14] * m[7] * m[11]
inv[25] = -m[0] * m[11] * m[17] + m[0] * m[12] * m[16] + m[9] * m[1] * m[17] - m[9] * m[2] * m[16] - m[14] * m[1] * m[12] + m[14] * m[2] * m[11]
inv[26] = m[0] * m[6] * m[17] - m[0] * m[7] * m[16] - m[4] * m[1] * m[17] + m[4] * m[2] * m[16] + m[14] * m[1] * m[7] - m[14] * m[2] * m[6]
inv[27] = -m[0] * m[6] * m[12] + m[0] * m[7] * m[11] + m[4] * m[1] * m[12] - m[4] * m[2] * m[11] - m[9] * m[1] * m[7] + m[9] * m[2] * m[6]
; Calculate the determinant
float det = m[0] * inv[0] + m[1] * inv[4] + m[2] * inv[8] + m[3] * inv[12] + m[4] * inv[16]
; Check if the determinant is zero (matrix is not invertible)
if det == 0
print "Error: Matrix is not invertible"
return
endif
; Multiply the adjugate matrix by the reciprocal of the determinant
float invDet = 1.0 / det
for i = 0 to 24
inv[i] *= invDet
next
; Update the matrix index
idx += 25
end method
method inverse()
float* m = @f + idx * 4
float* inv = @f + (idx + 25) * 4
; Create an identity matrix
for i = 0 to 24
inv[i] = 0.0
inv[0] = 1.0
inv[6] = 1.0
inv[12] = 1.0
inv[18] = 1.0
inv[24] = 1.0
; Perform Gaussian elimination
for i = 0 to 4
; Find pivot element
float pivot = m[i * 5 + i]
int pivotRow = i
for j = i + 1 to 4
if abs(m[j * 5 + i]) > abs(pivot)
pivot = m[j * 5 + i]
pivotRow = j
; Swap rows if necessary
if pivotRow != i
for k = 0 to 4
float temp = m[i * 5 + k]
m[i * 5 + k] = m[pivotRow * 5 + k]
m[pivotRow * 5 + k] = temp
temp = inv[i * 5 + k]
inv[i * 5 + k] = inv[pivotRow * 5 + k]
inv[pivotRow * 5 + k] = temp
; Divide row by pivot element
float invPivot = 1.0 / pivot
for k = 0 to 4
m[i * 5 + k] *= invPivot
inv[i * 5 + k] *= invPivot
; Subtract row from other rows to eliminate pivot element
for j = 0 to 4
if j != i
float factor = m[j * 5 + i]
for k = 0 to 4
m[j * 5 + k] -= factor * m[i * 5 + k]
inv[j * 5 + k] -= factor * inv[i * 5 + k]
; Copy the inverted matrix back to the original matrix
for i = 0 to 24
m[i] = inv[i]
idx += 25 ; Update the matrix index
end method
[quote]This method uses Gaussian elimination to calculate the inverse of the current matrix. Here's how it works:
It starts by creating an identity matrix of the same size as the current matrix. This identity matrix will be transformed into the inverse matrix during the elimination process.
It performs Gaussian elimination on the current matrix and the identity matrix simultaneously. The steps are as follows:
For each row, it finds the pivot element (the element with the largest absolute value in the current column).
If necessary, it swaps rows to bring the pivot element to the diagonal position.
It divides the current row by the pivot element to make the pivot element equal to 1.
It subtracts multiples of the current row from the other rows to eliminate the pivot element from those rows.
After the elimination process, the identity matrix has been transformed into the inverse matrix.
Finally, it copies the inverted matrix back to the original matrix and updates the matrix index.
Note that this method assumes that the matrix is invertible. If the matrix is singular or non-invertible, the method may produce incorrect results or encounter division by zero.
Add to Conversation[/quote]
Thanks Theo,
We don't need the full set of 3D matrix operations, but interpolation would be potentially very useful for taking 3D sections of a 4D object. For instance: slicing a 4D hypersphere along the w axis would produce 3D spheres of different sizes.
Here's the implementation of the interpolate function in O2 Basic, along with variants in x86 and x64 SIMD assembler:
O2 Basic:
method interpolate(Transform4D other, float t)
float* a = @f + idx * 4
float* b = @other.f + other.idx * 4
float* r = @f + (idx + 25) * 4
for i = 0 to 24
r[i] = a[i] + (b[i] - a[i]) * t
next
idx += 25 ; Update the matrix index
end method
x86 SIMD Assembler (SSE):
method interpolate_SSE(Transform4D other, float t)
float* a = @f + idx * 4
float* b = @other.f + other.idx * 4
float* r = @f + (idx + 25) * 4
; Load the interpolation factor into an SSE register
movss xmm7, t
shufps xmm7, xmm7, 0
; Perform interpolation using SSE instructions
mov ecx, 6 ; Loop counter (6 iterations for 24 elements)
mov esi, a
mov edi, b
mov edx, r
.loop:
movaps xmm0, [esi]
movaps xmm1, [edi]
subps xmm1, xmm0
mulps xmm1, xmm7
addps xmm0, xmm1
movaps [edx], xmm0
add esi, 16
add edi, 16
add edx, 16
dec ecx
jnz .loop
idx += 25 ; Update the matrix index
end method
x64 SIMD Assembler (AVX):
method interpolate_AVX(Transform4D other, float t)
float* a = @f + idx * 4
float* b = @other.f + other.idx * 4
float* r = @f + (idx + 25) * 4
; Load the interpolation factor into an AVX register
vbroadcastss ymm7, t
; Perform interpolation using AVX instructions
mov rcx, 3 ; Loop counter (3 iterations for 24 elements)
mov rsi, a
mov rdi, b
mov rdx, r
.loop:
vmovaps ymm0, [rsi]
vmovaps ymm1, [rdi]
vsubps ymm1, ymm1, ymm0
vmulps ymm1, ymm1, ymm7
vaddps ymm0, ymm0, ymm1
vmovaps [rdx], ymm0
add rsi, 32
add rdi, 32
add rdx, 32
dec rcx
jnz .loop
idx += 25 ; Update the matrix index
end method
In the O2 Basic implementation, the interpolation is performed element-wise using a loop.
The x86 SIMD assembler variant uses SSE instructions to perform the interpolation. It loads the interpolation factor into an SSE register and uses movaps, subps, mulps, and addps instructions to perform the interpolation on four elements at a time. The loop is unrolled to process 24 elements (6 iterations).
The x64 SIMD assembler variant uses AVX instructions, which can process eight elements at a time. It loads the interpolation factor into an AVX register using vbroadcastss and uses vmovaps, vsubps, vmulps, and vaddps instructions for the interpolation. The loop is unrolled to process 24 elements (3 iterations).
Note that the SIMD assembler code assumes the availability of SSE (x86) or AVX (x64) instructions on the target platform. Make sure to adapt the code to match the specific assembler syntax and conventions of your assembler and platform.
method interpolate(Transform4D other, float t)
if t < 0.0 or t > 1.0
print "Error: Interpolation factor t must be between 0 and 1."
return
endif
float* a = @f + idx * 4
float* b = @other.f + other.idx * 4
float* r = @f + (idx + 25) * 4
int i
for i = 0 to 24
r[i] = a[i] + (b[i] - a[i]) * t
next
idx += 25 ; Update the matrix index
end method
@Charles tell me if you need anything else. As you can see i have mutated to a 4D- ASM Expert "Over night" :-)
Note that AVX support is quite patchy. For instance my PC has an AMD Ryzen-3 processor pre 2020, and this does not support AVX.
https://en.wikipedia.org/wiki/Advanced_Vector_Extensions