Maths3D

BlitzMax Forums/BlitzMax Programming/Maths3D

JoshK(Posted 2008) [#1]
Maths3D code for games. Originally based on code posted by Mark Sibly:
Strict

Import BRL.Math
Import BRL.Blitz

Function Vec2:TVec2(x#=0.0,y#=0.0/0.0)
	Local t:TVec2=New TVec2
	If IsNan(y) y=x
	t.x=x
	t.y=y
	Return t
EndFunction

Function Vec3:TVec3( x#=0.0,y#=0.0/0.0,z#=0.0/0.0 )
	Local t:TVec3=New TVec3
	If IsNan(y) y=x
	If IsNan(z) z=y
	t.x=x
	t.y=y
	t.z=z
	Return t
EndFunction

Function Vec4:TVec4( x#=0.0,y#=0.0/0.0,z#=0.0/0.0,w#=0.0/0.0 )
	Local t:TVec4=New TVec4
	If IsNan(y) y=x
	If IsNan(z) z=y
	If IsNan(w) w=z
	t.x=x
	t.y=y
	t.z=z
	t.w=w
	Return t
EndFunction

Function Line:TLine( x#,y#,z#,dx#,dy#,dz# )
	Local t:TLine=New TLine
	t.x=x
	t.y=y
	t.z=z
	t.dx=dx
	t.dy=dy
	t.dz=dz
	Return t
EndFunction

Function Quat:TQuat( x#,y#,z#,w# )
	Local t:TQuat=New TQuat
	t.x=x
	t.y=y
	t.z=z
	t.w=w
	Return t
EndFunction

Function Plane:TPlane( nx#,ny#,nz#,d# )
	Local t:TPlane=New TPlane
	t.nx=nx
	t.ny=ny
	t.nz=nz
	t.d=d
	Return t
EndFunction

Function UlpsEqual( x#,y#,maxUlps=8 )
	Return DeltaUlps( x,y )<=maxUlps
EndFunction

Function DeltaUlps( x#,y# )
	Local ix=(Int Ptr Varptr x)[0]
	Local iy=(Int Ptr Varptr y)[0]
	If ix<0 ix=$80000000-ix
	If iy<0 iy=$80000000-iy
	Return Abs( ix-iy )
EndFunction

Function PointDistance#(p1:TVec3,p2:TVec3)
	Local dx#,dy#,dz#
	dx=p1.x-p2.x
	dy=p1.y-p2.y
	dz=p1.z-p2.z
	Return Sqr(dx*dx+dy*dy+dz*dz)
EndFunction

Function PointDistanceSqrd#(p1:TVec3,p2:TVec3)
	Local dx#,dy#,dz#
	dx=p1.x-p2.x
	dy=p1.y-p2.y
	dz=p1.z-p2.z
	Return dx*dx+dy*dy+dz*dz
EndFunction

Type TVec2
	Field x#,y#

	Method ToString$()
		Return "Vec3("+x+","+y+")"
	End Method

EndType

Type TVec3 {value}

	Field x#,y#,z# {attribute}
	
	Method Copy:TVec3()
		Return Vec3( x,y,z )
	End Method
	
	Method ToString$()
		Return "Vec3("+x+","+y+","+z+")"
	End Method
	
	Method ToVec4:TVec4()
		Return Vec4( x,y,z,1 )
	End Method
	
	Method ToField$()
		Return x+","+y+","+z
	End Method
	
	Method FromField:TVec3( t$ )
		Local bits$[]=t.Split( "," )
		If bits.length<>3 Throw "Format error"
		x=bits[0].ToFloat()
		y=bits[1].ToFloat()
		z=bits[2].ToFloat()
		Return Self
	End Method

	Method Pointer:Float Ptr()
		Return Varptr x
	EndMethod

	Method Length#()
		Return Sqr( x*x+y*y+z*z )
	End Method
	
	Method Dot#( v:TVec3 )
		Return x*v.x+y*v.y+z*v.z
	End Method
	
	Method Inverse:TVec3()
		Return Vec3( -x,-y,-z )
	End Method
	
	Method Reciprocal:TVec3()
		Return Vec3( 1/x,1/y,1/z )
	End Method
	
	Method Normalize:TVec3()
		Local t#=Sqr( x*x+y*y+z*z )
		If t=0.0 Return vec3(0.0,0.0,0.0)
		Return Vec3( x/t,y/t,z/t )
	End Method
	
	Method Scale:TVec3( scale# )
		Return Vec3( x*scale,y*scale,z*scale )
	End Method 
	
	Method DistanceTo#( v:TVec3 )
		Local dx#=x-v.x,dy#=y-v.y,dz#=z-v.z
		Return Sqr( dx*dx+dy*dy+dz*dz )
	End Method
	
	Method Plus:TVec3( v:TVec3 )
		Return Vec3( x+v.x,y+v.y,z+v.z )
	End Method
	
	Method Minus:TVec3( v:TVec3 )
		Return Vec3( x-v.x,y-v.y,z-v.z )
	End Method
	
	Method Times:TVec3( v:TVec3 )
		Return Vec3( x*v.x,y*v.y,z*v.z )
	End Method
	
	Method DividedBy:TVec3( v:TVec3 )
		Return Vec3( x/v.x,y/v.y,z/v.z )
	End Method
	
	Method Cross:TVec3( v:TVec3 )
		Return Vec3( y*v.z-z*v.y,z*v.x-x*v.z,x*v.y-y*v.x )
	End Method
	
	Method Blend:TVec3( v:TVec3,Alpha# )
		Local beta#=1-Alpha
		Return Vec3( x*beta+v.x*Alpha,y*beta+v.y*Alpha,z*beta+v.z*Alpha )
	End Method
	
	Method ToQuat:TQuat()
		Local c1#=Cos(-z/2)
		Local s1#=Sin(-z/2)
		Local c2#=Cos(-x/2)
		Local s2#=Sin(-x/2)
		Local c3#=Cos( y/2)
		Local s3#=Sin( y/2)
		Local c1_c2#=c1*c2
		Local s1_s2#=s1*s2
		Local t:TQuat=New TQuat
		t.x=c1*s2*c3 - s1*c2*s3;
		t.y=c1_c2*s3 + s1_s2*c3;
		t.z=s1*c2*c3 + c1*s2*s3;
		t.w=c1_c2*c3 - s1_s2*s3;
		Return t
	EndMethod
	
	Function FromQuat:TVec3(quat:TQuat)
		Return quat.toeuler()
	EndFunction
	
	Method Compare:Int(o:Object)
		Local t:TVec3=TVec3(o)
		If x<t.x Return -1
		If y<t.y Return -1
		If z<t.z Return -1
		If x>t.x Return 1
		If y>t.y Return 1
		If z>t.z Return 1
		Return 0
	EndMethod
	
EndType

Type TVec4' Extends TVec3

	Field x#,y#,z#,w# {attribute}
	
	Method Copy:TVec4()
		Return Vec4( x,y,z,w )
	End Method
	
	Method ToString$()
		Return "Vec4("+x+","+y+","+z+","+w+")"
	End Method
	
	Method ToVec3:TVec3()
		Return vec3(x,y,z)
	End Method
	
	Method XYZ:TVec3()
		Return vec3(x,y,z)
	EndMethod
	
	Method Plus:TVec4( v:TVec4 )
		Return Vec4( x+v.x,y+v.y,z+v.z,w+v.w )
	End Method

	Method Times:TVec4( v:TVec4 )
		Return Vec4( x*v.x,y*v.y,z*v.z,w*v.w )
	End Method	

	Method Pointer:Float Ptr()
		Return Varptr x
	EndMethod

EndType

Type TLine

	Field x#,y#,z#,dx#,dy#,dz# {attribute}
	
	Method Copy:TLine()
		Return Line( x,y,z,dx,dy,dz )
	End Method
	
	Method ToString$()
		Return "Line("+x+","+y+","+z+","+dx+","+dy+","+dz+")"
	End Method
	
	Method Origin:TVec3()
		Return Vec3( x,y,z )
	End Method
	
	Method Delta:TVec3()
		Return Vec3( dx,dy,dz )
	End Method
	
	Method EndPoint:TVec3()
		Return Vec3( x+dx,y+dy,z+dz )
	End Method
	
	Method Evaluate:TVec3( time# )
		Return Vec3( x+dx*time,y+dy*time,z+dz*time )
	End Method
	
	Method DistanceTo#( t:TLine )
		Local wx#=x-t.x
		Local wy#=y-t.y
		Local wz#=z-t.z
		Local a#=dx*dx+dy*dy+dz*dz
		Local b#=dx*t.dx+dy*t.dy+dz*t.dz
		Local c#=t.dx*t.dx+t.dy*t.dy+t.dz*t.dz
		Local d#=dx*wx+dy*wy+dz*wz
		Local e#=t.dx*wx+t.dy*wy+t.dz*wz
		Local q#=a*c-b*b
		Local st#=b*e-c*d
		Local tt#=a*e-b*d
		Local tx#=wx+(st*dx-tt*t.dx)/q
		Local ty#=wy+(st*dy-tt*t.dy)/q
		Local tz#=wz+(st*dz-tt*t.dz)/q
		Return Sqr( tx*tx+ty*ty+tz*tz )
	End Method

	Method Intersects:TVec3( t:TLine,tolerance#=0.001 )
		Local wx#=x-t.x
		Local wy#=y-t.y
		Local wz#=z-t.z
		Local a#=dx*dx+dy*dy+dz*dz
		Local b#=dx*t.dx+dy*t.dy+dz*t.dz
		Local c#=t.dx*t.dx+t.dy*t.dy+t.dz*t.dz
		Local d#=dx*wx+dy*wy+dz*wz
		Local e#=t.dx*wx+t.dy*wy+t.dz*wz
		Local q#=a*c-b*b
		Local st#=b*e-c*d
		Local tt#=a*e-b*d
		Local tx#=wx+(st*dx-tt*t.dx)/q
		Local ty#=wy+(st*dy-tt*t.dy)/q
		Local tz#=wz+(st*dz-tt*t.dz)/q
		If tolerance<>-1
			If tx*tx+ty*ty+tz*tz>=tolerance*tolerance Return
		EndIf	
		tx#=wx+(st*dx)/q
		ty#=wy+(st*dy)/q
		tz#=wz+(st*dz)/q
		Return vec3(tx,ty,tz).plus(vec3(t.x,t.y,t.z))
	EndMethod
	
	Function FromEndPoints:TLine( p0:TVec3,p1:TVec3 )
		Local t:TLine=New TLine
		t.x=p0.x;t.y=p0.y;t.z=p0.z
		t.dx=p1.x-p0.x;t.dy=p1.y-p0.y;t.dz=p1.z-p0.z;
		Return t
	EndFunction
	
EndType

Type TQuat

	Field x#,y#,z#,w#=1.0 {attribute}
	
	Method Copy:TQuat()
		Return Quat( x,y,z,w )
	End Method
	
	Method ToString$()
		Return "Quat("+x+","+y+","+z+","+w+")"
	End Method
	
	Method Inverse:TQuat()
		Return Quat( -x,-y,-z,w )
	End Method

	Method ToAngleAxis( angle# Var,axis:TVec3 Var )
		Local t#=ACos(w)
		angle=t*2
		If angle>0
			t=1/Sin(t)
			axis=Vec3( x*t,y*t,z*t )
		Else
			axis=Vec3(0,0,0)
		EndIf
	End Method
	
	Method Yaw#()
		Return ATan2( 2*y*w-2*z*x,1-2*y*y-2*x*x )
	End Method
	
	Method Pitch#()
		Return -ASin( 2*z*y+2*x*w )
	End Method
	
	Method Roll#()
		Return -ATan2( 2*z*w-2*y*x,1-2*z*z-2*x*x )
	End Method
	
	Method ToEuler:TVec3()
		Local x2#=x*x,y2#=y*y,z2#=z*z
		Local euler:TVec3=New TVec3
		euler.y#=ATan2( 2*y*w-2*z*x,1-2*y2-2*x2 )
		euler.x#=-ASin( 2*z*y+2*x*w )
		euler.z#=-ATan2( 2*z*w-2*y*x,1-2*z2-2*x2 )
		Return euler
	End Method
	
	Method Times:TQuat( q:TQuat )
		Local result:TQuat=New TQuat
		result.x#=w*q.x + x*q.w + y*q.z - z*q.y
		result.y#=w*q.y + y*q.w + z*q.x - x*q.z
		result.z#=w*q.z + z*q.w + x*q.y - y*q.x
		result.w#=w*q.w - x*q.x - y*q.y - z*q.z
		Return result
	EndMethod
	
	Method Slerp:TQuat( q:TQuat,a# )
		Local b#=1-a,f
		Local d#=x*q.x+y*q.y+z*q.z+w*q.w
		If d<0
			d=-d
			f=True
		EndIf
		If d<1
			Local om#=ACos( d )'-
			Local si#=Sin( om )'-
			a=Sin( a*om )/si'-
			b=Sin( b*om )/si'-
		EndIf
		If f a=-a
		Return Quat( x*b + q.x*a , y*b + q.y*a , z*b + q.z*a , w*b + q.w*a )
	EndMethod
	
	Method Normalize:TQuat()
		Local q:TQuat=New TQuat
		Local m#=Sqr(x*x+y*y+z*z+w*w)
		q.x=x/m
		q.y=y/m
		q.z=z/m
		q.w=w/m
		Return q
	EndMethod
	
	Function FromAngleAxis:TQuat( angle#,axis:TVec3 )
		angle:/2
		Local s#=Sin(angle)
		Local t:TQuat=New TQuat
		t.x=s*axis.x
		t.y=s*axis.y
		t.z=s*axis.z
		t.w=Cos(angle)
		Return t
	EndFunction
	
	Function FromEuler:TQuat(euler:TVec3)
		Return euler.toquat()
	EndFunction	

EndType

Type TPlane

	Field nx#,ny#,nz#,d# {attribute}
	
	Method Copy:TPlane()
		Return Plane( nx,ny,nz,d )
	End Method
	
	Method ToString$()
		Return "Plane("+nx+","+ny+","+nz+","+d+")"
	End Method
	
	Method MajorAxis()
		Local ax#=Abs(nx),ay#=Abs(ny),az#=Abs(nz)
		If ax>ay And ax>az Return 0
		If ay>az Return 1
		Return 2
	End Method
	
	Method SolveX:TVec3( p:TVec3 )
		Return Vec3( (p.y*ny+p.z*nz+d)/-nx,p.y,p.z )
	End Method
	
	Method SolveY:TVec3( p:TVec3 )
		Return Vec3( p.x,(p.x*nx+p.z*nz+d)/-ny,p.z )
	End Method
	
	Method SolveZ:TVec3( p:TVec3 )
		Return Vec3( p.x,p.y,(p.x*nx+p.y*ny+d)/-nz )
	End Method
	
	Method Inverse:TPlane()
		Return Plane( -nx,-ny,-nz,-d )
	End Method
	
	Method Normal:TVec3()
		Return Vec3( nx,ny,nz )
	End Method
	
	Method Normalize:TPlane()
		Local i#=1/Sqr(nx*nx+ny*ny+nz*nz)
		Return Plane( nx*i,ny*i,nz*i,d*i )
	End Method
	
	Method t_IntersectLine#( t:TLine )
		Return -DistanceToPoint( t.Origin() ) / Normal().Dot( t.Delta() )
	End Method
	
	Method IntersectsLine:TVec3(l:TLine)
		Local v:TVec3=New TVec3
		Local u#
		u#=(nx*l.x+ny*l.y+nz*l.z+d)/(nx*l.dx+ny*l.dy+nz*l.dz)
		v.x=l.x-l.dx*u
		v.y=l.y-l.dy*u
		v.z=l.z-l.dz*u
		Return v
	EndMethod
	
	Method DistanceToPoint#( p:TVec3 )
		Return nx*p.x+ny*p.y+nz*p.z+d
	End Method
	
	Function FromPointNormal:TPlane( p:TVec3,n:TVec3 )
		Return Plane( n.x,n.y,n.z,-p.Dot(n) )
	EndFunction
	
	Function FromTriangle:TPlane( v0:TVec3,v1:TVec3,v2:TVec3 )
		Local va:TVec3=v2.Minus(v1)
		Local vb:TVec3=v0.Minus(v1)
		Return FromPointNormal( v1,va.Cross(vb).Normalize() )
	EndFunction
	
	Function Ground:TPlane()
		Return Plane( 0,1,0,0 )
	EndFunction

EndType

Type TMat4

	Field ix#,iy#,iz#,iw# {attribute}
	Field jx#,jy#,jz#,jw# {attribute}
	Field kx#,ky#,kz#,kw# {attribute}
	Field tx#,ty#,tz#,tw# {attribute}

	Method Copy:TMat4()
		Local t:TMat4=New TMat4
		MemCopy t,Self,SizeOf Self
		Return t
	End Method
	
	Method ToString$()
		Local t$="Mat4{~n"
		t:+ix+","+iy+","+iz+","+iw+"~n"
		t:+jx+","+jy+","+jz+","+jw+"~n"
		t:+kx+","+ky+","+kz+","+kw+"~n"
		t:+tx+","+ty+","+tz+","+tw+"~n"
		Return t+"}~n"
	End Method
	
	Method Pointer:Float Ptr()
		Return Varptr ix
	EndMethod
	
	Method Translation:TVec3()
		Return Vec3(tx,ty,tz)
	End Method

	Rem
	Local sc:TVec3=Scale()
	Local iv:TVec3=Vec3(ix,iy,iz).Scale(1/sc.x)
	Local jv:TVec3=Vec3(jx,jy,jz).Scale(1/sc.y)
	Local kv:TVec3=Vec3(kx,ky,kz).Scale(1/sc.z)
	Local x#,y#,z#,w#
	Local t#=iv.x+jv.y+kv.z+1
	If t>0
		t=Sqr(t)*2
		x=(kv.y-jv.z)/t
		y=(iv.z-kv.x)/t
		z=(jv.x-iv.y)/t
		w=t/4
	Else If iv.x>jv.y And iv.x>kv.z
		t=Sqr(iv.x-jv.y-kv.z+1)*2.0
		x=t/4
		y=(jv.x+iv.y)/t
		z=(iv.z+kv.x)/t
		w=(kv.y-jv.z)/t
	Else If jv.y>kv.z
		t=Sqr(jv.y-kv.z-iv.x+1)*2.0
		x=(jv.x+iv.y)/t
		y=t/4
		z=(kv.y+jv.z)/t
		w=(iv.z-kv.x)/t
	Else
		t=Sqr(kv.z-jv.y-iv.x+1)*2.0
		x=(iv.z+kv.x)/t
		y=(kv.y+jv.z)/t
		z=t/4
		w=(jv.x-iv.y)/t
	EndIf
	Return Quat( x,y,z,w )
	EndRem
	
	Method Rotation:TQuat()
		Return normalize().Rotation2()
	EndMethod
	
	Method Rotation2:TQuat()	
		Local mult:Float
		Local x#,y#,z#,w#
		Local fourWSquaredMinus1:Float = +ix+jy+kz
		Local fourXSquaredMinus1:Float = +ix-jy-kz
		Local fourYSquaredMinus1:Float = +jy-ix-kz
		Local fourZSquaredMinus1:Float = +kz-ix-jy
		Local biggestIndex=0
		Local fourBiggestSquaredMinus1:Float=fourWSquaredMinus1
		If fourXSquaredMinus1 > fourBiggestSquaredMinus1
			fourBiggestSquaredMinus1 = fourXSquaredMinus1 
			biggestIndex=1
		EndIf
		If fourYSquaredMinus1 > fourBiggestSquaredMinus1
			fourBiggestSquaredMinus1 = fourYSquaredMinus1 
			biggestIndex=2
		EndIf
		If fourZSquaredMinus1 > fourBiggestSquaredMinus1
			fourBiggestSquaredMinus1 = fourZSquaredMinus1 
			biggestIndex=3
		EndIf
		Local biggestValue:Float=Sqr(fourBiggestSquaredMinus1+1.0)*0.5
		mult=0.25/biggestValue
		Select BiggestIndex
			Case 0
				w#=BiggestValue
				x#=(jz-ky)*mult
				y#=-(kx-iz)*mult
				z#=-(iy-jx)*mult
			Case 1
				x#=-BiggestValue
				w#=-( jz - ky )*mult
				y#=( iy + jx )*mult
				z#=( kx + iz )*mult
			Case 2
				y#=BiggestValue
				w#=-(kx-iz)*mult
				x#=-(iy+jx)*mult
				z#=(jz+ky)*mult
			Case 3
				z#=-BiggestValue
				w#=(iy-jx)*mult
				x#=(kx+iz)*mult
				y#=-(jz+ky)*mult
		EndSelect
		Return quat(-x,y,z,w)
	EndMethod
	
	Function KeyHit(blah)
	EndFunction
	
	Rem
	'testing code:
	Global mode=0
	If KeyHit(KEY_N)
		mode:+1
		'Notify mode
	EndIf
	If KeyHit(KEY_L) Notify mode
	Select mode		
		Case 0 Return quat(x,y,z,w)
		Case 1 Return quat(x,y,z,-w)
		Case 2 Return quat(x,-y,z,w)
		Case 3 Return quat(x,y,-z,w)
		Case 4 Return quat(-x,-y,z,w)
		Case 5 Return quat(x,-y,-z,w)
		Case 6 Return quat(-x,y,z,w)
	EndSelect
	EndRem
	
	Method Scale:TVec3()
		Return Vec3( Vec3(ix,iy,iz).Length(),Vec3(jx,jy,jz).Length(),Vec3(kx,ky,kz).Length() )
	End Method
	
	Method Times:TMat4( m:TMat4 )
		Local t:TMat4=New TMat4
		t.ix= ix*m.ix + jx*m.iy + kx*m.iz + tx*m.iw
		t.iy= iy*m.ix + jy*m.iy + ky*m.iz + ty*m.iw
		t.iz= iz*m.ix + jz*m.iy + kz*m.iz + tz*m.iw
		t.iw= iw*m.ix + jw*m.iy + kw*m.iz + tw*m.iw
		t.jx= ix*m.jx + jx*m.jy + kx*m.jz + tx*m.jw
		t.jy= iy*m.jx + jy*m.jy + ky*m.jz + ty*m.jw
		t.jz= iz*m.jx + jz*m.jy + kz*m.jz + tz*m.jw
		t.jw= iw*m.jx + jw*m.jy + kw*m.jz + tw*m.jw
		t.kx= ix*m.kx + jx*m.ky + kx*m.kz + tx*m.kw
		t.ky= iy*m.kx + jy*m.ky + ky*m.kz + ty*m.kw
		t.kz= iz*m.kx + jz*m.ky + kz*m.kz + tz*m.kw
		t.kw= iw*m.kx + jw*m.ky + kw*m.kz + tw*m.kw
		t.tx= ix*m.tx + jx*m.ty + kx*m.tz + tx*m.tw
		t.ty= iy*m.tx + jy*m.ty + ky*m.tz + ty*m.tw
		t.tz= iz*m.tx + jz*m.ty + kz*m.tz + tz*m.tw
		t.tw= iw*m.tx + jw*m.ty + kw*m.tz + tw*m.tw
		Return t
	End Method
	
	Method TimesPoint:TVec3( v:TVec3 )
'	Rem
		Return Vec3(..
			ix*v.x + jx*v.y + kx*v.z + tx,..
			iy*v.x + jy*v.y + ky*v.z + ty,..
			iz*v.x + jz*v.y + kz*v.z + tz )
'	End Rem
'		Local x#=ix*v.x + jx*v.y + kx*v.z + tx
'		Local y#=iy*v.x + jy*v.y + ky*v.z + ty
'		Local z#=iz*v.x + jz*v.y + kz*v.z + tz
'		Local w#=iw*v.x + jw*v.y + kw*v.z + tw
'		Return Vec3( x/w,y/w,z/w )
	End Method
	
	Method TimesVector:TVec3( v:TVec3 )
		Return Vec3(..
			ix*v.x + jx*v.y + kx*v.z,..
			iy*v.x + jy*v.y + ky*v.z,..
			iz*v.x + jz*v.y + kz*v.z )
	End Method

	Method TimesNormal:TVec3( v:TVec3 )
		Local m:TMat4=Inverse()
		Return Vec3(..
			m.ix*v.x + m.iy*v.y + m.iz*v.z,..
			m.jx*v.x + m.jy*v.y + m.jz*v.z,..
			m.kx*v.x + m.ky*v.y + m.kz*v.z )
	End Method
	
	Method TimesPlane:TPlane( p:TPlane )
		Local m:TMat4=Inverse()
		Return Plane(..
			m.ix*p.nx + m.iy*p.ny + m.iz*p.nz + m.iw*p.d,..
			m.jx*p.nx + m.jy*p.ny + m.jz*p.nz + m.jw*p.d,..
			m.kx*p.nx + m.ky*p.ny + m.kz*p.nz + m.kw*p.d,..
			m.tx*p.nx + m.ty*p.ny + m.tz*p.nz + m.tw*p.d ).Normalize()
	End Method
	
	Method TimesLine:TLine( t:TLine )
		Return TLine.FromEndPoints( TimesPoint( t.Origin() ),TimesPoint( t.EndPoint() ) )
	End Method
	
	Method Determinant#()
		Assert Abs(iw)<=.001 And Abs(jw)<=.001 And Abs(kw<=.001) And tw>=1-.001
		Return ix*(jy*kz-jz*ky) - iy*(jx*kz-jz*kx) + iz*(jx*ky-jy*kx)
	End Method
	
	Method I:TVec3()
		Return Vec3(ix,iy,iz)
	End Method
	
	Method J:TVec3()
		Return Vec3(jx,jy,jz)
	End Method
		
	Method K:TVec3()
		Return Vec3(kx,ky,kz)
	End Method
		
	Method Cofactor:TMat4()
		Local t:TMat4=New TMat4
		t.ix= (jy*kz-jz*ky) ; t.iy=-(jx*kz-jz*kx) ; 	t.iz= (jx*ky-jy*kx)
		t.jx=-(iy*kz-iz*ky) ; t.jy= (ix*kz-iz*kx) ; 	t.jz=-(ix*ky-iy*kx)
		t.kx= (iy*jz-iz*jy) ; t.ky=-(ix*jz-iz*jx) ; 	t.kz= (ix*jy-iy*jx)
		Return t
	End Method
	
	Method Inverse:TMat4()
		Local c#=1.0/Determinant()
		Local t:TMat4=New TMat4
		t.ix= c * ( jy*kz - jz*ky )
		t.iy=-c * ( iy*kz - iz*ky )
		t.iz= c * ( iy*jz - iz*jy )
		t.jx=-c * ( jx*kz - jz*kx )
		t.jy= c * ( ix*kz - iz*kx )
		t.jz=-c * ( ix*jz - iz*jx )
		t.kx= c * ( jx*ky - jy*kx )
		t.ky=-c * ( ix*ky - iy*kx )
		t.kz= c * ( ix*jy - iy*jx )
		t.tx=-( tx*t.ix + ty*t.jx + tz*t.kx )
		t.ty=-( tx*t.iy + ty*t.jy + tz*t.ky )
		t.tz=-( tx*t.iz + ty*t.jz + tz*t.kz )
		t.tw=1
		Return t
	End Method
	
	Method MyInverse:TMat4()
		Local t:TMat4=Transpose()
		t.iw=0
		t.jw=0
		t.kw=0
		t.tw=1
		t.tx=-( tx*t.ix + ty*t.jx + tz*t.kx )
		t.ty=-( tx*t.iy + ty*t.jy + tz*t.ky )
		t.tz=-( tx*t.iz + ty*t.jz + tz*t.kz )
		Return t
		'mat.tx = -( (ix * tx) + (grid[0,1] * ty) + (grid[0,2] * tz) )
		'mat.ty = -( (jx * tx) + (grid[1,1] * ty) + (grid[1,2] * tz) )
		'mat.tz = -( (kx * tx) + (grid[2,1] * ty) + (grid[2,2] * tz) )		
	EndMethod
	
	Method Transpose:TMat4()
		Local t:TMat4=New TMat4
		t.ix=ix ; t.iy=jx ; t.iz=kx ; t.iw=tx
		t.jx=iy ; t.jy=jy ; t.jz=ky ; t.jw=ty
		t.kx=iz ; t.ky=jz ; t.kz=kz ; t.kw=tz
		t.tx=iw ; t.ty=jw ; t.tz=kw ; t.tw=tw
		Return t
	End Method
	
	Function FromDirection:TMat4( dir:TVec3,up:TVec3=Null )
		If Not up up=Vec3(0,1,0)
		Local j:TVec3=up.Normalize()
		Local k:TVec3=dir.Normalize()
		Local i:TVec3=j.Cross(k).Normalize()
		j=k.Cross(i).Normalize()
		Local t:TMat4=New TMat4
		t.ix=i.x ; t.iy=i.y ; t.iz=i.z
		t.jx=j.x ; t.jy=j.y ; t.jz=j.z
		t.kx=k.x ; t.ky=k.y ; t.kz=k.z
		t.tw=1
		Return t
	EndFunction
	
	Function FromTranslation:TMat4( v:TVec3 )
		Local t:TMat4=New TMat4
		t.ix=1 ; t.jy=1 ;	t.kz=1
		t.tx=v.x ; t.ty=v.y ; t.tz=v.z
		t.tw=1
		Return t
	EndFunction
	
	Function FromRotation:TMat4( q:TQuat )
		Local t:TMat4=New TMat4
		Local xx#=q.x*q.x,yy#=q.y*q.y,zz#=q.z*q.z
		Local xy#=q.x*q.y,xz#=q.x*q.z,yz#=q.y*q.z
		Local wx#=q.w*q.x,wy#=q.w*q.y,wz#=q.w*q.z
		t.ix=1-2*(yy+zz) ; t.iy=  2*(xy-wz) ; t.iz=  2*(xz+wy)
		t.jx=  2*(xy+wz) ; t.jy=1.0-2*(xx+zz) ; t.jz=  2*(yz-wx)
		t.kx=  2*(xz-wy) ; t.ky=  2*(yz+wx) ; t.kz=1.0-2*(xx+yy)
		t.tw=1
		Return t
	EndFunction
	
	Function FromScale:TMat4( v:TVec3 )
		Local t:TMat4=New TMat4
		t.ix=v.x
		t.jy=v.y
		t.kz=v.z
		t.tw=1
		Return t
	EndFunction
	
	Function FromUniformScale:TMat4( sc# )
		Local t:TMat4=New TMat4
		t.ix=sc
		t.jy=sc
		t.kz=sc
		t.tw=1
		Return t
	EndFunction
	
	Function FromTransRotScale:TMat4( trans:TVec3,rot:TQuat,scale:TVec3 )
		Local t:TMat4
		If trans
			t=FromTranslation(trans)
			If rot t=t.Times(FromRotation(rot))
			If scale t=t.Times(FromScale(scale))
		Else If rot
			t=FromRotation(rot)
			If scale t=t.Times(FromScale(scale))
		Else If scale
			t=FromScale(scale)
		Else
			t=Identity()
		EndIf
		Return t
	EndFunction
	
	Function FromYaw:TMat4( yaw# )	'untested!
		Local t:TMat4=New TMat4
		Local s#=Sin(yaw),c#=Cos(yaw)
		t.ix=c  ; t.iz=s ; t.jy=1
		t.kx=-s ; t.kz=c ; t.tw=1
		Return t
	EndFunction
	
	Function FromPitch:TMat4( pitch# )	'untested!
		Local t:TMat4=New TMat4
		Local s#=Sin(pitch),c#=Cos(pitch)
		t.ix=1  ; t.jy=c ; t.jz=s
		t.ky=-s ; t.kz=c ; t.tw=1
		Return t
	EndFunction
	
	Function FromRoll:TMat4( roll# )	'untested!
		Local t:TMat4=New TMat4
		Local s#=Sin(roll),c#=Cos(roll)
		t.ix=c ; t.iy=s ; t.jx=-s
		t.jy=c ; t.kz=1 ; t.tw=1
		Return t
	EndFunction
	
	Function FromYawPitchRoll:TMat4( yaw#,pitch#,roll# )
		Local t:TMat4
		If yaw
			t=FromYaw(yaw)
			If pitch t=t.Times(FromPitch(pitch))
			If roll t=t.Times(FromRoll(roll))
		Else If pitch
			t=FromPitch(pitch)
			If roll t=t.Times(FromRoll(roll))
		Else If roll
			t=FromRoll(roll)
		Else
			t=Identity()
		EndIf
		Return t
	EndFunction
	
	Function FromOrtho:TMat4( l#,r#,b#,t#,n#,f# )
		Local q:TMat4=New TMat4
		q.ix=2/(r-l)
		q.tx=-(r+l)/(r-l)
		q.jy=2/(t-b)
		q.ty=-(t+b)/(t-b)
		q.kz=2/(f-n)
		q.tz=-(f+n)/(f-n)
		q.tw=1
		Return q
	EndFunction
	
	Function FromFrustum:TMat4( near_left#,near_right#,near_bottom#,near_top#,near#,far# )
		Local t:TMat4=New TMat4
		Local near2#=near*2
		Local w#=near_right-near_left
		Local h#=near_top-near_bottom
		Local d#=far-near
		t.ix=near2/w
		t.jy=near2/h
		t.kx=(near_right+near_left)/w
		t.ky=(near_top+near_bottom)/h
		t.kz=(far+near)/d
		t.kw=1
		t.tz=-(far*near2)/d
		Return t
	EndFunction
	
	Function FromInfiniteFrustum:TMat4( near_left#,near_right#,near_bottom#,near_top#,near# )
		Local t:TMat4=New TMat4
		Local near2#=near*2
		Local w#=near_right-near_left
		Local h#=near_top-near_bottom
		Local e#=.001
		t.ix=near2/w
		t.jy=near2/w
		t.kx=(near_right+near_left)/w
		t.ky=(near_top+near_bottom)/h
		t.kz=1-e
		t.kw=1
		t.tz=near*(e-1)
		Return t
	EndFunction

	Function Identity:TMat4()
		Local t:TMat4=New TMat4
		t.ix=1 ; t.jy=1 ; t.kz=1 ; t.tw=1
		Return t
	EndFunction
	
	Method Normalize:TMat4()
		Local mat:TMat4=New TMat4
		Local vec:TVec3
		vec=i().normalize()
		mat.ix=vec.x
		mat.iy=vec.y
		mat.iz=vec.z
		vec=j().normalize()
		mat.jx=vec.x
		mat.jy=vec.y
		mat.jz=vec.z
		vec=k().normalize()
		mat.kx=vec.x
		mat.ky=vec.y
		mat.kz=vec.z
		mat.tx=tx
		mat.ty=ty
		mat.tz=tz
		mat.tw=1.0
		Return mat
	EndMethod
	
EndType

Type TAABB
	
	Field x0#,y0#,z0#,x1#,y1#,z1#,cx#,cy#,cz#,w#,h#,d#,radius# { attribute }
	
	Method Reset()
		x0=0.0
		y0=0.0
		z0=0.0
		x1=0.0
		y1=0.0
		z1=0.0
		cx#=0.0
		cy#=0.0
		cz#=0.0
		w#=0.0
		h#=0.0
		d#=0.0
	EndMethod
	
	'Rem
	Method TForm:TAABB(src:TMat4,dst:TMat4)
		Local vec:TVec3
		
		Local aabb:TAABB
		aabb=New taabb
		
		vec=TFormPointm(vec3(x0,y0,z0),src,dst)
		aabb.x0=vec.x
		aabb.y0=vec.y
		aabb.z0=vec.z
		aabb.x1=vec.x
		aabb.y1=vec.y
		aabb.z1=vec.z
		
		vec=TFormPointm(vec3(x0,y0,z1),src,dst)
		aabb.x0=Min(aabb.x0,vec.x)
		aabb.y0=Min(aabb.y0,vec.y)
		aabb.z0=Min(aabb.z0,vec.z)
		aabb.x1=Max(aabb.x1,vec.x)
		aabb.y1=Max(aabb.y1,vec.y)
		aabb.z1=Max(aabb.z1,vec.z)

		vec=TFormPointm(vec3(x0,y1,z0),src,dst)
		aabb.x0=Min(aabb.x0,vec.x)
		aabb.y0=Min(aabb.y0,vec.y)
		aabb.z0=Min(aabb.z0,vec.z)
		aabb.x1=Max(aabb.x1,vec.x)
		aabb.y1=Max(aabb.y1,vec.y)
		aabb.z1=Max(aabb.z1,vec.z)

		vec=TFormPointm(vec3(x1,y0,z0),src,dst)
		aabb.x0=Min(aabb.x0,vec.x)
		aabb.y0=Min(aabb.y0,vec.y)
		aabb.z0=Min(aabb.z0,vec.z)
		aabb.x1=Max(aabb.x1,vec.x)
		aabb.y1=Max(aabb.y1,vec.y)
		aabb.z1=Max(aabb.z1,vec.z)

		vec=TFormPointm(vec3(x1,y1,z0),src,dst)
		aabb.x0=Min(aabb.x0,vec.x)
		aabb.y0=Min(aabb.y0,vec.y)
		aabb.z0=Min(aabb.z0,vec.z)
		aabb.x1=Max(aabb.x1,vec.x)
		aabb.y1=Max(aabb.y1,vec.y)
		aabb.z1=Max(aabb.z1,vec.z)

		vec=TFormPointm(vec3(x0,y1,z1),src,dst)
		aabb.x0=Min(aabb.x0,vec.x)
		aabb.y0=Min(aabb.y0,vec.y)
		aabb.z0=Min(aabb.z0,vec.z)
		aabb.x1=Max(aabb.x1,vec.x)
		aabb.y1=Max(aabb.y1,vec.y)
		aabb.z1=Max(aabb.z1,vec.z)
		
		vec=TFormPointm(vec3(x1,y0,z1),src,dst)
		aabb.x0=Min(aabb.x0,vec.x)
		aabb.y0=Min(aabb.y0,vec.y)
		aabb.z0=Min(aabb.z0,vec.z)
		aabb.x1=Max(aabb.x1,vec.x)
		aabb.y1=Max(aabb.y1,vec.y)
		aabb.z1=Max(aabb.z1,vec.z)

		vec=TFormPointm(vec3(x1,y1,z1),src,dst)
		aabb.x0=Min(aabb.x0,vec.x)
		aabb.y0=Min(aabb.y0,vec.y)
		aabb.z0=Min(aabb.z0,vec.z)
		aabb.x1=Max(aabb.x1,vec.x)
		aabb.y1=Max(aabb.y1,vec.y)
		aabb.z1=Max(aabb.z1,vec.z)
		
		aabb.Update()
		
		Return aabb
	EndMethod
	
	Method Update()
		w=x1-x0
		h=y1-y0
		d=z1-z0
		cx=x0+w*0.5
		cy=y0+h*0.5
		cz=z0+d*0.5
		UpdateRadius()
	EndMethod
	

	
	Method IntersectsPoint:Int(p:TVec3,radius#=0.0)
		If p.x<x0-radius Return False
		If p.y<y0-radius Return False
		If p.z<z0-radius Return False
		If p.x>x1+radius Return False
		If p.y>y1+radius Return False
		If p.z>z1+radius Return False
		Return True
	EndMethod

	Method IntersectsAABB:Int(aabb:TAABB)
		If aabb.x0>x1 Return False
		If aabb.y0>y1 Return False
		If aabb.z0>z1 Return False
		If aabb.x1<x0 Return False
		If aabb.y1<y0 Return False
		If aabb.z1<z0 Return False
		Return True
	EndMethod
	
	Function FromCloud:TAABB(vertex_count,vertex_cloud:Float Ptr)
		Local aabb:TAABB=New TAABB
		Local v:Int
		For v=0 To vertex_count-1
			If v=0 Or vertex_cloud[v*3+0]<aabb.x0 aabb.x0=vertex_cloud[v*3+0]
			If v=0 Or vertex_cloud[v*3+1]<aabb.y0 aabb.y0=vertex_cloud[v*3+1]
			If v=0 Or vertex_cloud[v*3+2]<aabb.z0 aabb.z0=vertex_cloud[v*3+2]
			If v=0 Or vertex_cloud[v*3+0]>aabb.x1 aabb.x1=vertex_cloud[v*3+0]
			If v=0 Or vertex_cloud[v*3+1]>aabb.y1 aabb.y1=vertex_cloud[v*3+1]
			If v=0 Or vertex_cloud[v*3+2]>aabb.z1 aabb.z1=vertex_cloud[v*3+2]
		Next
		aabb.Update()
		Return aabb
	EndFunction
	
	Method UpdateRadius()
		Local dx#,dy#,dz#
		dx=w*0.5
		dy=h*0.5
		dz=d*0.5
		radius=Sqr(dx*dx+dy*dy+dz*dz)
	EndMethod

	Method MinExtent:TVec3()
		Return vec3(x0,y0,z0)
	EndMethod
	
	Method MaxExtent:TVec3()
		Return vec3(x1,y1,z1)
	EndMethod
	
	Method Center:TVec3()
		Return vec3(cx,cy,cz)
	EndMethod
	
	Method Copy:TAABB()
		Local aabb:TAABB
		aabb=New TAABB
		aabb.x0=x0
		aabb.y0=y0
		aabb.z0=z0
		aabb.x1=x1
		aabb.y1=y1
		aabb.z1=z1
		aabb.cx=cx
		aabb.cy=cy
		aabb.cz=cz
		aabb.w=w
		aabb.h=h
		aabb.d=d
		aabb.radius=radius
		Return aabb
	EndMethod
	
	Method IntersectsPlane:Int(plane:TPlane)
		Local d#,n
		Local greater
		Local lesser
		Local p:TVec3
		For n=1 To 8
			Select n
				Case 1 p=vec3(x0,y0,z0)
				Case 2 p=vec3(x1,y0,z0)
				Case 3 p=vec3(x0,y1,z0)
				Case 4 p=vec3(x0,y0,z1)
				Case 5 p=vec3(x1,y1,z1)
				Case 6 p=vec3(x0,y1,z1)
				Case 7 p=vec3(x1,y0,z1)
				Case 8 p=vec3(x1,y1,z1)
			EndSelect
			d#=plane.DistanceToPoint(p)
			If d>0.0 greater=True
			If d<0.0 lesser=True
			If greater=True And lesser=True Return 0
		Next
		If greater=True Return 1
		If lesser=True Return -1
	EndMethod
	
	Method ToString$()
		Return x0+", "+y0+", "+z0+", "+x1+", "+y1+", "+z1+Chr(13)+Chr(10)+cx+", "+cy+", "+cz
	EndMethod
	
	Method Intersects:Int(aabb:TAABB,bias#=0.0)
		If x0>aabb.x1-bias Return
		If y0>aabb.y1-bias Return
		If z0>aabb.z1-bias Return
		If aabb.x0>x1-bias Return
		If aabb.y0>y1-bias Return
		If aabb.z0>z1-bias Return
		Return 1
	EndMethod
	
EndType

Function TFormPointm:TVec3(p:TVec3,mat1:TMat4,mat2:TMat4)
	Local mat:TMat4
	If mat1
		mat=mat1.times(TMat4.FromTranslation(p))
		p=mat.Translation()
	EndIf
	If mat2
		mat=mat2.inverse().times(TMat4.FromTranslation(p))
		p=mat.Translation()
	EndIf
	If mat1=Null And mat2=Null Return p.copy() Else Return p
EndFunction

Function TFormVectorm:TVec3(p:TVec3,mat1:TMat4,mat2:TMat4)
	Local mat:TMat4
	If mat1
		mat=mat1.copy()
		mat.tx=0.0 ; mat.ty=0.0 ; mat.tz=0.0
		mat=mat.times(TMat4.FromTranslation(p))
		p=mat.Translation()
	EndIf
	If mat2
		mat=mat2.inverse()
		mat.tx=0.0 ; mat.ty=0.0 ; mat.tz=0.0
		mat=mat.times(TMat4.FromTranslation(p))
		p=mat.Translation()
	EndIf
	If mat1=Null And mat2=Null Return p.copy() Else Return p
EndFunction

Function TFormNormalm:TVec3(p:TVec3,mat1:TMat4,mat2:TMat4)
	Return TFormVectorm(p,mat1,mat2).normalize()
EndFunction

Function TFormQuatm:TQuat(quat:TQuat,mat1:TMat4,mat2:TMat4)
	Local mat:TMat4
	If mat1
		mat=mat1.copy()
		mat.tx=0.0 ; mat.ty=0.0 ; mat.tz=0.0
		mat=mat.times(TMat4.FromRotation(quat))
		quat=mat.rotation()
	EndIf
	If mat2
		mat=mat2.inverse()
		mat.tx=0.0 ; mat.ty=0.0 ; mat.tz=0.0
		mat=mat.times(TMat4.FromRotation(quat))
		quat=mat.rotation().normalize()
	EndIf
	If mat1=Null And mat2=Null Return quat.copy() Else Return quat	
EndFunction



JoshK(Posted 2008) [#2]
Also added my own AABB class.


ImaginaryHuman(Posted 2008) [#3]
Looks like a good reference. Thanks.

I don't personally understand eulers and quarternians or whatever, at this point. For that matter I don't understand matrices. Maybe someone can point us simpler folk to a very easy to understand tutorial or math primer which lays it all out in very simple terms.


Gabriel(Posted 2008) [#4]
Eulers are what everyone knows and uses by default. Axis Angle rotations between 0 and 360, and three axes, X, Y and Z. It's easy to visualize for us humans, but very poor mathematically speaking because the order you do the rotations in defines where you end up. So the same apparent rotation can have several different meanings. Not good from a mathematical standpoint. Not to mention you cannot interpolate between orientations and the Gimbal lock problem.

Quaternions are a much better, much more straightforward approach, but they're nowhere near as easy for the human brain to visualize. Essentially a quaternion is a set of four values describing a single rotation around one axis. The axis is completely arbitrary ( not just x, y and z, but any axis ) and so a single rotation can achieve any rotation possible in the Euler system. X, Y and Z define the axis, and W defines the size of rotation around that axis. Values are between 0 and 1.

Quaternions can be interpolated, and do not suffer from Gimbal lock, so they're very suitable for games use. Also, because they only have four values, they're faster than using matrices.

Let someone else explain matrices, I'm allergic.

There are countless internet resources, and I wouldn't want to pick one out in particular as I found none were exhaustive and I actually needed to pick up bits and pieces from various places. The 3D Math Primer is an excellent book on maths if you're interested.


Dreamora(Posted 2008) [#5]
Very nice, thank you.
But some points on it:

1. There is a GL call in. Not really the "nice" programming style for a pure math lib

2. The AABB looks nice but is bound against TEntity which is not present and basically could be from 5 engines but most likely yours :)

3. The { } in there could easily be removed. You don't assign anything to the attributes so they are worthless reflection wise (which is their only purpose)


JoshK(Posted 2008) [#6]
Mark added the gl command, not me.

I prefer using this to return a float pointer for different commands:

Method Pointer:Float Ptr()
Return Varptr ix
EndMethod


BlitzSupport(Posted 2008) [#7]

For that matter I don't understand matrices. Maybe someone can point us simpler folk to a very easy to understand tutorial or math primer which lays it all out in very simple terms.



This is really good... even I understood it when I read it, though it's now all forgotten:

http://chortle.ccsu.edu/VectorLessons/vectorIndex.html


Azathoth(Posted 2008) [#8]
3. The { } in there could easily be removed. You don't assign anything to the attributes so they are worthless reflection wise (which is their only purpose)
Not so, reflection can be used to tell if the attribute exists. My clone function makes use of this http://www.blitzbasic.com/codearcs/codearcs.php?code=2132


Dreamora(Posted 2008) [#9]
OK, that indeed is right ...
But it isn't used here at all


JoshK(Posted 2008) [#10]
I just used the original code. Don't complain to me.


JoshK(Posted 2008) [#11]
Here is a method for a plane to intersect a line:
	Method IntersectsLine:TVec3(l:TLine)
		Local v:TVec3=New TVec3
		Local u#
		u#=(-nx*l.x-ny*l.y-nz*l.z-d)/(-nx*(l.x-l.dx)-ny*(l.x-l.dy)-nz*(l.x-l.dz))
		v.x=l.x+l.dx*u
		v.y=l.y+l.dy*u
		v.z=l.z+l.dz*u
		Return v
	EndMethod



JoshK(Posted 2008) [#12]
Added transformation commands for points, vectors, normals, and quaternions.