'================
'MATRIX MATHS ETC
'================
string tab=chr(9)
string cr=chr(13,10)
macro ShowBuffer(n)
===================
string bu=space n
end macro
macro ShowTitle(v)
==================
mid bu,c,v cr : c+=2+len v
end macro
macro ShowElement(v)
====================
mid bu,c," "+str(v,4) : c+=14
end macro
macro ShowNextLine()
====================
mid bu,c,cr : c+=2
end macro
macro ShowReturn()
==================
return left bu,c-1
end macro
function ShowVec(float *v, sys n=3) as string
=============================================
indexbase 0
ShowBuffer 256
sys c=1
n--
ShowTitle "Vector:"
for i=0 to n
ShowElement v[i]
next
ShowNextLine
ShowReturn
end function
function ShowMat(float*m,sys r=4) as string
===========================================
indexbase 0
ShowBuffer 512
string s
sys c,e,i,j,k,n
static double min=1.0e-7
c=1
n=r-1
e=r*r-1
ShowTitle "Matrix:"
for i=0 to n 'column
for j=0 to e step r 'row
k=i+j
ShowElement m[k]
next
ShowNextLine
next
ShowReturn
end function
function Transpose(float *r, *m, int n=4)
=========================================
indexbase 0
sys i,j,k1,k2
sys e=n-1
for i=0 to e
k2=i
for j=0 to e
r[k1]=m[k2] : k1++ : k2+=n
next
next
end function
function Identity(float *m, int n=4)
====================================
indexbase 0
sys h=n+1, e=n*n-1
for i=0 to e
m[i]=0.0
next
for i=0 to e step h
m[i]=1.0
next
end function
function det2x2sub(float *m, int i0,i1,i2,i3) as float
======================================================
indexbase 0
return m[i0] * m[i3] - m[i2] * m[i1]
end function
function InverseMat3x3(float *r, *m)
====================================
indexbase 0
float det = 0.0f
'
det += m[0] * det2x2sub(m, 4, 5, 7, 8)
det -= m[3] * det2x2sub(m, 1, 2, 7, 8)
det += m[6] * det2x2sub(m, 1, 2, 4, 5)
'
r[0] = det2x2sub(m, 4, 5, 7, 8) / det
r[1] = -det2x2sub(m, 1, 2, 7, 8) / det
r[2] = det2x2sub(m, 1, 2, 4, 5) / det
r[3] = -det2x2sub(m, 3, 5, 6, 8) / det
r[4] = det2x2sub(m, 0, 2, 6, 8) / det
r[5] = -det2x2sub(m, 0, 2, 3, 5) / det
r[6] = det2x2sub(m, 3, 4, 6, 7) / det
r[7] = -det2x2sub(m, 0, 1, 6, 7) / det
r[8] = det2x2sub(m, 0, 1, 3, 4) / det
end function
function Mat4x4ToMat3x3(float *r,*m)
====================================
indexbase 0
'copy top left section
r[0] = m[0]
r[3] = m[4]
r[6] = m[8]
r[1] = m[1]
r[4] = m[5]
r[7] = m[9]
r[2] = m[2]
r[5] = m[6]
r[8] = m[10]
end function
function mat4x4Mul(float *r,*m1,*m2)
====================================
indexbase 0
r[0] = m1[0] * m2[0] + m1[4] * m2[1] + m1[8] * m2[2] + m1[12] * m2[3]
r[1] = m1[1] * m2[0] + m1[5] * m2[1] + m1[9] * m2[2] + m1[13] * m2[3]
r[2] = m1[2] * m2[0] + m1[6] * m2[1] + m1[10] * m2[2] + m1[14] * m2[3]
r[3] = m1[3] * m2[0] + m1[7] * m2[1] + m1[11] * m2[2] + m1[15] * m2[3]
r[4] = m1[0] * m2[4] + m1[4] * m2[5] + m1[8] * m2[6] + m1[12] * m2[7]
r[5] = m1[1] * m2[4] + m1[5] * m2[5] + m1[9] * m2[6] + m1[13] * m2[7]
r[6] = m1[2] * m2[4] + m1[6] * m2[5] + m1[10] * m2[6] + m1[14] * m2[7]
r[7] = m1[3] * m2[4] + m1[7] * m2[5] + m1[11] * m2[6] + m1[15] * m2[7]
r[8] = m1[0] * m2[8] + m1[4] * m2[9] + m1[8] * m2[10] + m1[12] * m2[11]
r[9] = m1[1] * m2[8] + m1[5] * m2[9] + m1[9] * m2[10] + m1[13] * m2[11]
r[10] = m1[2] * m2[8] + m1[6] * m2[9] + m1[10] * m2[10] + m1[14] * m2[11]
r[11] = m1[3] * m2[8] + m1[7] * m2[9] + m1[11] * m2[10] + m1[15] * m2[11]
r[12] = m1[0] * m2[12] + m1[4] * m2[13] + m1[8] * m2[14] + m1[12] * m2[15]
r[13] = m1[1] * m2[12] + m1[5] * m2[13] + m1[9] * m2[14] + m1[13] * m2[15]
r[14] = m1[2] * m2[12] + m1[6] * m2[13] + m1[10] * m2[14] + m1[14] * m2[15]
r[15] = m1[3] * m2[12] + m1[7] * m2[13] + m1[11] * m2[14] + m1[15] * m2[15]
end function
function transformVec4(float *r,*m,*v, sys n=0)
===============================================
'using mat4x4
indexbase 0
sys i
for i=0 to n
r[0] = m[0] * v[0] + m[4] * v[1] + m[8] * v[2] + m[12] * v[3]
r[1] = m[1] * v[0] + m[5] * v[1] + m[9] * v[2] + m[13] * v[3]
r[2] = m[2] * v[0] + m[6] * v[1] + m[10] * v[2] + m[14] * v[3]
r[3] = m[3] * v[0] + m[7] * v[1] + m[11] * v[2] + m[15] * v[3]
@v+=16
@r+=16
next
end function
function transformVec3(float *r,*m,*v, sys n=1)
===============================================
'using mat4x4
indexbase 0
sys i
n--
for i=0 to n
r[0] = m[0] * v[0] + m[4] * v[1] + m[8] * v[2] + m[12]
r[1] = m[1] * v[0] + m[5] * v[1] + m[9] * v[2] + m[13]
r[2] = m[2] * v[0] + m[6] * v[1] + m[10] * v[2] + m[14]
@v+=12
@r+=12
next
end function
macro rotsetup
==============
indexbase 0
float m[16]
a = rad(a)
float ca = cos(a), sa = sin(a)
end macro
function mat4x4rotateVecX(float *r, *v, a)
==========================================
'
' 1 0 0 0
' 0 c -s 0
' 0 s c 0
' 0 0 0 0
'
rotsetup
m[0] = 1
'm[1] = 0
'm[2] = 0
'm[3] = 0
----------
'm[4] = 0
m[5] = ca
m[6] = sa
'm[7] = 0
----------
'm[8] = 0
m[9] = -sa
m[10] = ca
'm[11] = 0
transformVec3 r,m,v
end function
function mat4x4rotateVecY(float *r, *v, a)
==========================================
'
' c s 0 0
' 0 1 0 0
'-s 0 c 0
' 0 0 0 0
'
rotsetup
m[0] = ca
'm[1] = 0
m[2] = -sa
'm[3] = 0
----------
'm[4] = 0
m[5] = 1
'm[6] = 0
'm[7] = 0
----------
m[8] = sa
'm[9] = 0
m[10] = ca
'm[11] = 0
transformVec3 r,m,v
end function
function mat4x4rotateVecZ(float *r, *v, a)
==========================================
'
' c -s 0 0
' s c 0 0
' 0 0 1 0
' 0 0 0 0
'
rotsetup
m[0] = ca
m[1] = sa
'm[2] = 0
'm[3] = 0
----------
m[4] = -sa
m[5] = ca
'm[6] = 0
'm[7] = 0
----------
'm[8] = 0
'm[9] = 0
m[10] = 1
'm[11] = 0
transformVec3 r,m,v
end function
function det3x3sub(float *m, int i0,i1,i2,i3,i4,i5,i6,i7,i8) as float
=====================================================================
indexbase 0
float det
det += m[i0] * det2x2sub(m, i4, i5, i7, i8)
det -= m[i3] * det2x2sub(m, i1, i2, i7, i8)
det += m[i6] * det2x2sub(m, i1, i2, i4, i5)
return det;
end function
function mat4x4Inverse(float *r,*m)
===================================
indexbase 0
float det
det += m[0] * det3x3sub(m, 5, 6, 7, 9, 10, 11, 13, 14, 15)
det -= m[4] * det3x3sub(m, 1, 2, 3, 9, 10, 11, 13, 14, 15)
det += m[8] * det3x3sub(m, 1, 2, 3, 5, 6, 7, 13, 14, 15)
det -= m[12] * det3x3sub(m, 1, 2, 3, 5, 6, 7, 9, 10, 11)
r[0] = det3x3sub(m, 5, 6, 7, 9, 10, 11, 13, 14, 15) / det;
r[1] = -det3x3sub(m, 1, 2, 3, 9, 10, 11, 13, 14, 15) / det;
r[2] = det3x3sub(m, 1, 2, 3, 5, 6, 7, 13, 14, 15) / det;
r[3] = -det3x3sub(m, 1, 2, 3, 5, 6, 7, 9, 10, 11) / det;
r[4] = -det3x3sub(m, 4, 6, 7, 8, 10, 11, 12, 14, 15) / det;
r[5] = det3x3sub(m, 0, 2, 3, 8, 10, 11, 12, 14, 15) / det;
r[6] = -det3x3sub(m, 0, 2, 3, 4, 6, 7, 12, 14, 15) / det;
r[7] = det3x3sub(m, 0, 2, 3, 4, 6, 7, 8, 10, 11) / det;
r[8] = det3x3sub(m, 4, 5, 7, 8, 9, 11, 12, 13, 15) / det;
r[9] = -det3x3sub(m, 0, 1, 3, 8, 9, 11, 12, 13, 15) / det;
r[10] = det3x3sub(m, 0, 1, 3, 4, 5, 7, 12, 13, 15) / det;
r[11] = -det3x3sub(m, 0, 1, 3, 4, 5, 7, 8, 9, 11) / det;
r[12] = -det3x3sub(m, 4, 5, 6, 8, 9, 10, 12, 13, 14) / det;
r[13] = det3x3sub(m, 0, 1, 2, 8, 9, 10, 12, 13, 14) / det;
r[14] = -det3x3sub(m, 0, 1, 2, 4, 5, 6, 12, 13, 14) / det;
r[15] = det3x3sub(m, 0, 1, 2, 4, 5, 6, 8, 9, 10) / det;
end function
function mat4x4scale(float *r, x, y, z)
=======================================
indexbase 0
r[0] = x
r[5] = y
r[10] = z
r[15] = 1
end function
function mat4x4Translate(float *r, x, y, z)
===========================================
indexbase 0
r[12] = x
r[13] = y
r[14] = z
r[0] = 1
r[5] = 1
r[10] = 1
r[15] = 1
end function
function mat4x4rotate(float *r, a, x, y, z)
===========================================
indexbase 0
a = rad(a)
double f=1/sqr(x*x+y*y+z*z)
x*=f : y*=f : z*=f 'normalise
float ca = 1.0 - cos(a), sa = sin(a)
r[0] = 1.0 + ca * (x * x - 1.0)
r[1] = ca * x * y + z * sa
r[2] = ca * x * z - y * sa
r[4] = ca * x * y - z * sa
r[5] = 1.0 + ca * (y * y - 1.0)
r[6] = ca * y * z + x * sa
r[8] = ca * x * z + y * sa
r[9] = ca * y * z - x * sa
r[10] = 1.0 + ca * (z * z - 1.0)
r[15] = 1.0
end function
type vec3 float x,y,z
=====================
function Normalize(vec3 *r,*v)
==============================
double f=1/sqr(v.x*v.x+v.y*v.y+v.z*v.z)
r.x=v.x*f
r.y=v.y*f
r.z=v.z*f
end function
function dot(vec3 *u,*v) as float
=================================
return u.x * v.x + u.y * v.y + u.z * v.z
end function
function cross(vec3 *r,*u,*v)
=============================
r.x = u.y * v.z - u.z * v.y
r.y = u.z * v.x - u.x * v.z
r.z = u.x * v.y - u.y * v.x
end function
function diff(vec3 *r,*u,*v)
============================
r.x=u.x - v.x
r.y=u.y - v.y
r.z=u.z - v.z
end function
function mat4x4Look(float *r, vec3 *eye,*center,*up)
====================================================
indexbase 0
vec3 vx,vy,vz,vt
diff(eye,center,vt)
normalize(vt,vz)
'
cross(up,vz,vt)
normalize(vt,vx)
'
cross(vz,vx,vy)
'
r[0] = vx.x
r[1] = vy.x
r[2] = vz.x
r[4] = vx.y
r[5] = vy.y
r[6] = vz.y
r[8] = vx.z
r[9] = vy.z
r[10] = vz.z
r[12] = -dot(vx, eye)
r[13] = -dot(vy, eye)
r[14] = -dot(vz, eye)
end function
function CopyFloat(float *r,*s, int n=1)
========================================
indexbase 0
sys i
n--
for i=0 to n
r[i] = s[i]
next
end function
function CopyVec(float *r,*s, int n=1)
======================================
indexbase 0
sys i
n*=3
n--
for i=0 to n
r[i] = s[i]
next
end function
function CopyMat(float *r,*s, int n=1)
======================================
indexbase 0
sys i
n*=16
n--
for i=0 to n
r[i] = s[i]
next
end function
'
'CAMERA SPACE
'http://ogldev.atspace.co.uk/www/tutorial13/tutorial13.html
'
function SetViewMatrix(float *r, vec3 *vp,*vu,*vv,*vn)
======================================================
indexbase 0
vec3 u,v,n
float p[16], m[16]
'
identity p,4
p[12]={-vp.x, -vp.y, -vp.z}
'
normalize u,vu
normalize v,vv
normalize n,vn
'itr cross product u from n and v
'
' U V N
' UP RIGHT LOOKAT
--------------------------------------
m[0] = {u.x, v.x, n.x, 0.0} 'column 0
m[4] = {u.y, v.y, n.y, 0.0} 'column 1
m[8] = {u.z, v.z, n.z, 0.0} 'column 2
m[12]= {0.0, 0.0, 0.0, 1.0} 'column 3
--------------------------------------
'
Mat4x4Mul r,m,p
'
end function
'http://stackoverflow.com/questions/18404890/how-to-build-perspective-projection-matrix-no-api
'http://gamedev.stackexchange.com/questions/56609/how-to-create-a-projection-matrix-in-opengl-es-2-0
'**
'http://www.songho.ca/opengl/gl_projectionmatrix.html
'http://www.terathon.com/gdc07_lengyel.pdf
'
'http://www.swiftless.com/tutorials/opengl4/3-opengl-4-matrices.html
'
function SetProjectionMatrix(float *r, fov, aspect, near, far, int orth=0, leftHanded=1)
========================================================================================
'FOV RECEIVED IN DEGREES
'ASPECT == VIEWPORT WIDTH/HEIGHT
'
indexbase 0
identity r,4
float h = 1 : if not lefthanded then h=-1
float i = 1 / (far-near)
float f = 1 / tan(rad(fov/2))
'
r[5] = aspect * f 'uh
r[0] = h * f 'uw
'
if orth then
'
'ORTHOGRAPHIC PROJECTION
-------------------------------------------------
' uw 0 0 0
' 0 uh 0 0
' 0 0 -2/(f-n) -(f+n)/(f-n)
' 0 0 0 1
-------------------------------------------------
'
r[10] = -2.0 * i
r[11] = 0
r[14] = -(far+near) * i
r[15] = 1
exit function
end if
'
'PERSPECTIVE PROJECTION
-------------------------------------------------
' uw 0 0 0
' 0 uh 0 0
' 0 0 -(f+n)/(f-n) -2fn/(f-n)
' 0 0 -1 0
-------------------------------------------------
r[10] = -(far+near) * i
r[11] = -1
r[14] = -2.0 * far * near * i
r[15] = 0
end function
'===========
'QUATERNIONS
'===========
'
'http://ogldev.atspace.co.uk/www/tutorial15/tutorial15.html
'http://www.cprogramming.com/tutorial/3d/quaternions.html
'
' Hamilton
'
'SLERP, or Spherical Linear intERPolation
'
' Q = xi + yj + zk + w
'
' i*i == j*j == k*k = ijk == -1
'
' 1/Q == -xi - yi - zk +w
'
'ROTATION QUATERNION
'
' W = QVQ^-1
'
' Q = {
' v.x * sin(a/2),
' v.y * sin(a/2),
' v.z * sin(a/2),
' cos(a/2)
' }
function NormQuaternion(float *r)
=================================
indexbase 0
double f = 1 / sqr(r[0]*r[0] + r[1]*r[1] + r[2]*r[2] + r[3]*r[3])
r[0]*=f
r[1]*=f
r[2]*=f
r[3]*=f
end function
function SetQuaternion(float *r,a,x,y,z)
========================================
indexbase 0
a=rad(a/2) 'degrees to radians
float sa=sin(a)
r[0]=x*sa
r[1]=y*sa
r[2]=z*sa
r[3]=cos(a)
NormQuaternion r
end function
function ResetQuaternion(float *r)
==================================
r={0,0,0,1}
end function
function MulQuaternion(float *r,*a,*b)
======================================
'0 x
'1 y
'2 z
'3 w
indexbase 0
'itr test first for deviation from unit
NormQuaternion a
NormQuaternion b
'
r[3] = a[3]*b[3] - a[0]*b[0] - a[1]*b[1] - a[2]*[b[2]
r[0] = a[3]*b[0] + a[0]*b[3] + a[1]*b[2] - a[2]*[b[1]
r[1] = a[3]*b[1] - a[0]*b[2] + a[1]*b[3] + a[2]*[b[0]
r[2] = a[3]*b[2] + a[0]*b[1] - a[1]*b[0] + a[2]*[b[3]
end function
'generate local_rotation (details later)
'Qtotal = Qlocal * Qtotal //multiplication order matters on this line
'
'
function MulQuaternion(float *qTotal,*qLocal)
=============================================
float q[4]
MulQuaternion q,qLocal,qTotal
copyfloat qTotal,q,4
end function
function MatFromQuat(float *r,*q)
=================================
indexbase 0
r[0] = 1-2*q[1]*q[1] - 2*q[2]*q[2]
r[1] = 2*q[0]*q[1] + 2*q[3]*q[2]
r[2] = 2*q[0]*q[2] - 2*q[3]*q[1]
r[3] = 0
r[4] = 2*q[0]*q[1] - 2 * q[3]*q[2]
r[5] = 1 - 2*q[0]*q[0] - 2*q[2]*q[2]
r[6] = 2*q[1]*q[2] - 2*q[3]*q[0]
r[7] = 0
r[8] = 2*q[0]*q[2] + 2*q[3]*q[1]
r[9] = 2*q[1]*q[2] + 2*q[3]*q[0]
r[10] = 1-2*q[0]*q[0] - 2*q[1]*q[1]
r[11] = 0
r[12] = 0
r[13] = 0
r[14] = 0
r[15] = 1
end function
==============================================
'FUNCTIONS EMULATING CLASSIC OPENGL TRANSFORMS
==============================================
float MatStack[800], sys MatIdx 'stack depth: 50
identity MatStack,4 : MatIdx=0
function scalef(float x,y,z)
============================
indexbase 0
float m[32]
mat4x4scale m[0],x,y,z
mat4x4mul m[16],MatStack[MatIdx],m[0]
copymat MatStack[MatIdx],m[16]
end function
function translatef(float x,y,z)
================================
indexbase 0
float m[32]
mat4x4translate m[0],x,y,z
mat4x4mul m[16], MatStack[MatIdx],m[0]
copymat MatStack[MatIdx],m[16]
end function
function rotatef(float a,x,y,z)
===============================
indexbase 0
float m[32]
mat4x4rotate m[0],a,x,y,z
mat4x4mul m[16],MatStack[MatIdx],m[0]
copymat MatStack[MatIdx],m[16]
end function
function PushMatrix()
=====================
indexbase 0
sys msd=MatIdx+16
if msd<800 then
copymat matstack[msd], matstack[MatIdx]
MatIdx=msd
end if
end function
function PopMatrix()
====================
if MatIdx>0 then MatIdx-=16
end function
function LoadIdentity()
=======================
indexbase 0
identity matstack[MatIdx],4
end function
'SUPPORTING FUNCTIONS
function ResetMatrix()
======================
MatIdx=0
LoadIdentity
end function
function MulMatrix(float *m)
============================
indexbase 0
float r[16]
mat4x4mul r,MatStack[MatIdx],m
copymat MatStack[MatIdx],r
end function
function inversef()
===================
indexbase 0
float m[16]
mat4x4inverse m[0],MatStack[MatIdx]
copymat MatStack[MatIdx],m[0]
end function
function RotateAxis(float a,x,y,z)
==================================
float q[4], m[16]
SetQuaternion q,a, x,y,z
MatFromQuat m,q
MulMatrix m
end function
function transformf(float *v, int n=1,st=12,oo=0)
=================================================
'
indexbase 0
float m at (@MatStack[MatIdx])
float r[3]
sys i
@v+=oo 'OFFSET
n-- 'FINAL VERTEX INDEX
for i=0 to n
r[0] = m[0] * v[0] + m[4] * v[1] + m[8] * v[2] + m[12]
r[1] = m[1] * v[0] + m[5] * v[1] + m[9] * v[2] + m[13]
r[2] = m[2] * v[0] + m[6] * v[1] + m[10] * v[2] + m[14]
v={r[0],r[1],r[2]}
@v+=st 'STRIDE
next
end function
