Tools for 4 Dimensional Geometry / 5X5 Matrix Transforms

[Charles Pegge]

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.

Code Select Expand

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


[Theo Gottwald]

How about this?

Code Select Expand

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.

Code Select Expand

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


[Theo Gottwald]

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.

Code Select Expand

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

[Charles Pegge]

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

[Theo Gottwald]

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

Code Select Expand


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


Code Select Expand


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]


[Charles Pegge]

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.

[Theo Gottwald]

Here's the implementation of the interpolate function in O2 Basic, along with variants in x86 and x64 SIMD assembler:

O2 Basic:

Code Select Expand

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):

Code Select Expand

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):

Code Select Expand

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.

Code Select Expand

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" :-)

[Charles Pegge]

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