3d without openGL
-
- Posts: 3906
- Joined: Jan 01, 2009 7:03
- Location: Australia
Re: 3d without openGL
@dodicat,
Yes your 3d demos look great but I can't visualize or understand how your compact code works.
May I suggest you add keyboard control over the rotation and translation of the object? It is only then that it would be visually clear that the computations allow rotation around an arbitrary axis. In my wire frame rocket example I show the absolute world coordinates and also the rocket coordinates. The problem was to rotate around any one of the rocket's coordinates.
.
Yes your 3d demos look great but I can't visualize or understand how your compact code works.
May I suggest you add keyboard control over the rotation and translation of the object? It is only then that it would be visually clear that the computations allow rotation around an arbitrary axis. In my wire frame rocket example I show the absolute world coordinates and also the rocket coordinates. The problem was to rotate around any one of the rocket's coordinates.
.
Re: 3d without openGL
@paul doe - Very intrigued and inspired by your offerings
Re: 3d without openGL
OK basiccoder2.
I use only the necessary for a Saturn type ring of stones.
Use the arrow keys.
I use only the necessary for a Saturn type ring of stones.
Use the arrow keys.
Code: Select all
Type V3
As Single x,y,z
As Uinteger col
End Type
Type angle3D 'FLOATS for angles
As Single sx,sy,sz
As Single cx,cy,cz
Declare static Function construct(As Single,As Single,As Single) As Angle3D
End Type
Type float As V3
#define incircle(cx,cy,radius,x,y) (cx-x)*(cx-x) +(cy-y)*(cy-y)<= radius*radius
Function Angle3D.construct(x As Single,y As Single,z As Single) As Angle3D
Return Type (Sin(x),Sin(y),Sin(z), _
Cos(x),Cos(y),Cos(z))
End Function
Function Rotate(c As V3,p As V3,a As Angle3D,scale As float=type(1,1,1)) As V3
Dim As Single dx=p.x-c.x,dy=p.y-c.y,dz=p.z-c.z
Return Type<V3>((scale.x)*((a.cy*a.cz)*dx+(-a.cx*a.sz+a.sx*a.sy*a.cz)*dy+(a.sx*a.sz+a.cx*a.sy*a.cz)*dz)+c.x,_
(scale.y)*((a.cy*a.sz)*dx+(a.cx*a.cz+a.sx*a.sy*a.sz)*dy+(-a.sx*a.cz+a.cx*a.sy*a.sz)*dz)+c.y,_
(scale.z)*((-a.sy)*dx+(a.sx*a.cy)*dy+(a.cx*a.cy)*dz)+c.z,p.col)
End Function
Function perspective(p As V3,eyepoint As V3) As V3
Dim As Single w=1+(p.z/eyepoint.z)
Return Type<V3>((p.x-eyepoint.x)/w+eyepoint.x,_
(p.y-eyepoint.y)/w+eyepoint.y,_
(p.z-eyepoint.z)/w+eyepoint.z,p.col)
End Function
Sub QsortZ(array() As V3,begin As Long,Finish As Ulong)
Dim As Long i=begin,j=finish
Dim As V3 x =array(((I+J)\2))
While I <= J
While array(I).z > X .z:I+=1:Wend
While array(J).z < X .z:J-=1:Wend
If I<=J Then Swap array(I),array(J): I+=1:J-=1
Wend
If J >begin Then QsortZ(array(),begin,J)
If I <Finish Then QsortZ(array(),I,Finish)
End Sub
Function Regulate(Byval MyFps As long,Byref fps As long) As long
Static As Double timervalue,_lastsleeptime,t3,frames
frames+=1
If (Timer-t3)>=1 Then t3=Timer:fps=frames:frames=0
Var sleeptime=_lastsleeptime+((1/myfps)-Timer+timervalue)*1000
If sleeptime<1 Then sleeptime=1
_lastsleeptime=sleeptime
timervalue=Timer
Return sleeptime
End Function
screen 19,32
const pi=4*atn(1)
redim as v3 a(0)
'set some points in a saturn type disc
for n as long=0 to 100000
var xp=rnd*800,yp=rnd*600
if incircle(400,300,150,xp,yp)=0 then
if incircle(400,300,250,xp,yp) then
var u=ubound(a)
redim preserve a(1 to u+1)
a(ubound(a))=type<v3>(xp,yp,rnd*10-rnd*10,rnd*rgb(255,255,255))
pset(xp,yp),a(ubound(a)).col
end if
end if
next
print "the points are on the x y plane with a small random z component"
print "flip them all 90 degrees around y axis and expand a bit"
print "Press a key"
sleep
cls
redim as V3 rot(lbound(a) to ubound(a))'to hold all rotated points
dim as Angle3D A3d=angle3D.construct(0,pi/2,0) 'flip all points by pi/2 on y axis
for n as long=lbound(a) to ubound(a)
rot(n)=rotate(type(400,300,0),a(n),A3D,type(1.4,1.4,1.4)) '1.4 is a scale up
a(n)=rot(n)
pset(a(n).x,a(n).y),a(n).col
next
print "Thats 'em all flipped without perspctive"
print "Press a key"
sleep
dim as float ang
dim as long fps
'fix these two components at start
ang.z=pi/2 'tilt angle sideways
ang.y=-pi/7 'tilt angle back and forward
dim as string key
do
key=inkey
If key=Chr(255)+"K" Then ang.z-=.02 'left
If key=Chr(255)+"M" Then ang.z+=.02 'right
If key=Chr(255)+"P" Then ang.y-=.02 'down
If key=Chr(255)+"H" Then ang.y+=.02 'up
ang.x+=.01 'the orbiting speed
A3D=Angle3D.construct(ang.x,ang.y,ang.z)'get the six rotate components (sines, coses ...)
screenlock
cls
draw string(50,50),"Framerate "&fps
draw string(50,150),"Use the arrow keys"
for n as long=lbound(a) to ubound(a)
rot(n)=rotate(type(400,300,0),a(n),A3D,type(1,1,1))
rot(n)=perspective(rot(n),type(400,300,1000))
next
qsortZ(rot(),lbound(rot),ubound(rot))
for n as long=lbound(rot) to ubound(rot)
pset(rot(n).x,rot(n).y),rot(n).col
next
screenunlock
sleep regulate(64,fps),1
loop until key=chr(27)
-
- Posts: 3906
- Joined: Jan 01, 2009 7:03
- Location: Australia
Re: 3d without openGL
Demonstration without enough explanation for my programming limitations.
What I tried to do was modify your code to produce the three axes of the object (without the object). That is the axes were the object. I tried to make the z blue the x red and the y green but they appear to be random colors. I wanted to be able to rotate around any of those three axes.
What I tried to do was modify your code to produce the three axes of the object (without the object). That is the axes were the object. I tried to make the z blue the x red and the y green but they appear to be random colors. I wanted to be able to rotate around any of those three axes.
Code: Select all
Type V3
As Single x,y,z
As Uinteger col
End Type
Type angle3D 'FLOATS for angles
As Single sx,sy,sz
As Single cx,cy,cz
Declare static Function construct(As Single,As Single,As Single) As Angle3D
End Type
Type float As V3
#define incircle(cx,cy,radius,x,y) (cx-x)*(cx-x) +(cy-y)*(cy-y)<= radius*radius
Function Angle3D.construct(x As Single,y As Single,z As Single) As Angle3D
Return Type (Sin(x),Sin(y),Sin(z), _
Cos(x),Cos(y),Cos(z))
End Function
Function Rotate(c As V3,p As V3,a As Angle3D,scale As float=type(1,1,1)) As V3
Dim As Single dx=p.x-c.x,dy=p.y-c.y,dz=p.z-c.z
Return Type<V3>((scale.x)*((a.cy*a.cz)*dx+(-a.cx*a.sz+a.sx*a.sy*a.cz)*dy+(a.sx*a.sz+a.cx*a.sy*a.cz)*dz)+c.x,_
(scale.y)*((a.cy*a.sz)*dx+(a.cx*a.cz+a.sx*a.sy*a.sz)*dy+(-a.sx*a.cz+a.cx*a.sy*a.sz)*dz)+c.y,_
(scale.z)*((-a.sy)*dx+(a.sx*a.cy)*dy+(a.cx*a.cy)*dz)+c.z,p.col)
End Function
Function perspective(p As V3,eyepoint As V3) As V3
Dim As Single w=1+(p.z/eyepoint.z)
Return Type<V3>((p.x-eyepoint.x)/w+eyepoint.x,_
(p.y-eyepoint.y)/w+eyepoint.y,_
(p.z-eyepoint.z)/w+eyepoint.z,p.col)
End Function
Sub QsortZ(array() As V3,begin As Long,Finish As Ulong)
Dim As Long i=begin,j=finish
Dim As V3 x =array(((I+J)\2))
While I <= J
While array(I).z > X .z:I+=1:Wend
While array(J).z < X .z:J-=1:Wend
If I<=J Then Swap array(I),array(J): I+=1:J-=1
Wend
If J >begin Then QsortZ(array(),begin,J)
If I <Finish Then QsortZ(array(),I,Finish)
End Sub
Function Regulate(Byval MyFps As long,Byref fps As long) As long
Static As Double timervalue,_lastsleeptime,t3,frames
frames+=1
If (Timer-t3)>=1 Then t3=Timer:fps=frames:frames=0
Var sleeptime=_lastsleeptime+((1/myfps)-Timer+timervalue)*1000
If sleeptime<1 Then sleeptime=1
_lastsleeptime=sleeptime
timervalue=Timer
Return sleeptime
End Function
screenres 1280,600,32
const pi=4*atn(1)
redim as v3 a(0)
'set some points
for n as single =-250 to 250
var xp=0,yp=0,zp = n
var u=ubound(a)
redim preserve a(1 to u+1)
a(ubound(a))=type<v3>(xp, yp, zp, rnd*rgb(0,0,255))
pset(xp,yp),a(ubound(a)).col
next n
for n as single =-250 to 250
var xp=0,yp=n,zp = 0
var u=ubound(a)
redim preserve a(1 to u+1)
a(ubound(a))=type<v3>(xp, yp, zp, rnd*rgb(255,0,0))
pset(xp,yp),a(ubound(a)).col
next n
for n as single =-250 to 250
var xp=n,yp=0,zp = 0
var u=ubound(a)
redim preserve a(1 to u+1)
a(ubound(a))=type<v3>(xp, yp, zp, rnd*rgb(0,255,0))
pset(xp,yp),a(ubound(a)).col
next n
'print "the points are on the x y plane with a small random z component"
'print "flip them all 90 degrees around y axis and expand a bit"
'print "Press a key"
'sleep
'cls
redim as V3 rot(lbound(a) to ubound(a))'to hold all rotated points
dim as Angle3D A3d=angle3D.construct(0,pi/2,0) 'flip all points by pi/2 on y axis
for n as long=lbound(a) to ubound(a)
rot(n)=rotate(type(400,300,0),a(n),A3D,type(0.5,0.5,0.5)) '1.4 is a scale up
a(n)=rot(n)
pset(a(n).x,a(n).y),a(n).col
next
'print "Thats 'em all flipped without perspctive"
'print "Press a key"
'sleep
dim as float ang
dim as long fps
'fix these two components at start
ang.z=pi/2 'tilt angle sideways
ang.y=-pi/7 'tilt angle back and forward
dim as string key
do
key=inkey
If key=Chr(255)+"K" Then ang.z-=.02 'left
If key=Chr(255)+"M" Then ang.z+=.02 'right
If key=Chr(255)+"P" Then ang.y-=.02 'down
If key=Chr(255)+"H" Then ang.y+=.02 'up
if key="z" then ang.x -=.02
if key="x" then ang.x +=.02
A3D = Angle3D.construct(ang.x,ang.y,ang.z)'get the six rotate components (sines, coses ...)
screenlock
cls
draw string(50,50),"Framerate "&fps
draw string(50,150),"Use the arrow keys"
for n as long=lbound(a) to ubound(a)
rot(n)=rotate(type(400,300,0),a(n),A3D,type(1,1,1))
rot(n)=perspective(rot(n),type(400,300,1000))
next
qsortZ(rot(),lbound(rot),ubound(rot))
for n as long=lbound(rot) to ubound(rot)
pset(rot(n).x,rot(n).y),rot(n).col
next
screenunlock
sleep regulate(64,fps),1
loop until key=chr(27)
Re: 3d without openGL
Exactly. Thanks for pointing it out. I can imagine that you had a lot of fun trying to avoid the lock with Euler angles X-D. Don't fret. We all used Euler, at one point of another. In physics is the same, you stick with an Euler integrator until the simulation craps out on you ('explodes' to infinity).owen wrote:@ paul doe
I couldn't agree with you more.This problem disappears if you use matrices
And not only that, using matrices avoids gimbal lock
in my example you can see it encounter gimbal lock by first rotating 90 degrees on the x axis, then 90 on y axis and then stop at 89 on the z.
when you get to 89 on the z axis and proceed to rotate to 90 look closely and notice a tiny glitch.
what a wonderful experience it was to try and avoid gimbal lock to no avail. well actually i kinda did it but what length i went to.
im glad that my example can serve the perfect example of what not to do. LOL
Fun fact: the Apollo 11 mission almost failed thanks to the cheap a$$es of the NASA, because they didn't wanted to install another gimbal in the landing module to cut costs. Look at the Wikipedia entry for Gimbal Lock, quite amusing LOL.
@dodicat: pretty nice!. Reminds me of Point Clouds. I had some code on that, will look for it.
@BasicCoder2: aha. I see it now. I think that you are confusing rotation about an arbitrary axis with rotation around the object's own axes. One is extrinsic (outside the object's coordinate system) and the other is intrinsic. I will now post some code that will help you clarify that (hopefully).
Re: 3d without openGL
Hi dafhidafhi wrote:@paul doe - Very intrigued and inspired by your offerings
Glad to hear that. Which one in particular sparkled your interest?
Last edited by paul doe on Oct 11, 2017 11:33, edited 1 time in total.
Re: 3d without openGL
BasicCoder2:BasicCoder2 wrote:Demonstration without enough explanation for my programming limitations.
What I tried to do was modify your code to produce the three axes of the object (without the object). That is the axes were the object. I tried to make the z blue the x red and the y green but they appear to be random colors. I wanted to be able to rotate around any of those three axes.
As I promised, here is some code that will help to understand how coordinate systems work. You will need five more files (besides the ones that I already posted). Check the next post.
Last edited by paul doe on Oct 11, 2017 11:35, edited 2 times in total.
Re: 3d without openGL
camera.bas
object.bas
UPDATE 2: 13/11/2017 - camera.bas: corrected some more bugs and added a few methods that I forgot to port. Sorry :'(
UPDATE 1: 10/11/2017 - camera.bas: killed a few bugses and sprinkled some more comments
Code: Select all
#include once "platform.bi" '' platform specific definitions
#include once "core.bi" '' core functions and defines
#include once "math.bi" '' math functions and constants
#include once "vec4.bi" '' vec4 type and operations
#include once "mat4.bi" '' mat4 type and operations
/'
this serves to define a projection plane
in the code, x and y coordinates are not used. It only serves to replace
the IRenderContext interface
'/
type rectangle
as double x
as double y
as double width
as double height
end type
/'
simple camera class
responsibility: represents a pinhole-model camera. Can be used to navigate a scene FPS style
very useful for debugging purposes (especially geometry code)
mutability: mutable. Use with caution or solely for debugging. Just don't try to write Doom 5 with this.
10/10/2017: refactored the code because it's a mess.
'/
class camera
public:
declare constructor( _
byval position as vec4, _
byval x as vec4, _
byval y as vec4, _
byval z as vec4, _
byval near as double = 1.0, _
byval far as double = 10000.0, _
byval projPlane as rectangle )
'' lots of crappy functions to play with
declare function getU() as vec4
declare function getV() as vec4
declare function getDir() as vec4
declare function getPos() as vec4
declare sub move( byval offset as vec4 )
declare sub rotate( byval axis as vec4, byval angle as double )
declare function getRatioUV() as double
declare function getRatioVU() as double
declare sub setRatioUV( byval ratio as double )
declare sub setRatioVU( byval ratio as double)
declare function getZoomU() as double
declare function getZoomV() as double
declare sub setZoomU( byval a as double )
declare sub setZoomV( byval a as double )
declare sub zoom( byval a as double )
declare sub zoomU( byval a as double )
declare sub zoomV( byval a as double )
declare function getFOVV() as double
declare function getFOVU() as double
declare sub setFOVU( byval angle as double )
declare sub setFOVV( byval angle as double )
declare function getYaw() as double
declare function getPitch() as double
declare function getRoll() as double
declare sub setYaw( byval angle as double )
declare sub setPitch( byval angle as double )
declare sub setRoll( byval angle as double )
declare sub yawAbsolute( byval angle as double )
declare sub yaw( byval angle as double )
declare sub pitch( byval angle as double )
declare sub roll( byval angle as double )
declare function getCameraMatrix() as mat4
declare function getInverseCameraMatrix() as mat4
declare function transform( byval v as vec4 ) as vec4
declare function orthogonal( byval pnt as vec4, byval plane as rectangle, byref x as double, byref y as double, byref z as double ) as integer
declare function perspective( byval pnt as vec4, byval plane as rectangle, byref x as double, byref y as double, byref z as double ) as integer
declare function projectOnScreen( byval pnt as vec4, byval plane as rectangle, byref x as double, byref y as double, byref z as double ) as integer
declare function projectOnScreenOrtho( byval pnt as vec4, byval plane as rectangle, byref x as double, byref y as double, byref z as double ) as integer
declare function getScale() as double
declare sub setScale( byval dirLength as double )
declare sub scale( byval factor as double )
declare sub resetSkewU()
declare sub resetSkewV()
declare function getDistance( byval pnt as vec4 ) as double
declare sub setDistance( byval pnt as vec4, byval dist as double )
declare sub setLookDirection( byval newDir as vec4 )
declare sub lookAt( byval target as vec4 )
declare sub lookAt( byval target as vec4, byval up as vec4 )
declare sub setU( byval newU as vec4 )
declare sub setV( byval newV as vec4 )
declare sub setDir( byval newDir as vec4 )
declare sub setPos( byval newPos as vec4 )
declare function nearClip() as double
declare function farClip() as double
declare function projectionPlane() as rectangle
private:
declare sub generateMatrix()
'' members
as vec4 m_pos
as vec4 m_u
as vec4 m_v
as vec4 m_dir
as mat4 m_position
as mat4 m_camMatrix
as mat4 m_invCamMatrix
as double m_nearClip
as double m_farClip
as double m_fov
as rectangle m_projPlane
end class
constructor camera( _
byval position as vec4, _
byval x as vec4, _
byval y as vec4, _
byval z as vec4, _
byval near as double = 1.0, _
byval far as double = 10000.0, _
byval projPlane as rectangle )
m_pos = position
m_u = x
m_v = y
m_dir = z
m_nearClip = near
m_farClip = far
m_position = matrices.translation( position )
m_projPlane = projPlane
/'
adjust the aspect ratio of the camera to match the
projection plane resolution
sorry for this crap, I had to reimplement it this way to shorten
the code
'/
with m_projPlane
if( .width < .height ) then
setRatioUV( .width / .height )
else
setRatioVU( .height / .width )
end if
end with
generateMatrix()
end constructor
sub camera.generateMatrix()
'' generates the camera's rotation matrix
with m_camMatrix
.a = m_u.x : .b = m_v.x : .c = m_dir.x : .d = 0.0
.e = m_u.y : .f = m_v.y : .g = m_dir.y : .h = 0.0
.i = m_u.z : .j = m_v.z : .k = m_dir.z : .l = 0.0
.m = m_u.w : .n = m_v.w : .o = m_dir.w : .p = 1.0
end with
'' and calculates the inverse immediately
m_invCamMatrix = inverse( m_camMatrix )
end sub
sub camera.move( byval offset as vec4 )
m_pos = m_pos + offset
end sub
sub camera.rotate( byval axis as vec4, byval angle as double )
'' rotates the camera
m_u = rotateAroundAxis( m_u, axis, angle )
m_v = rotateAroundAxis( m_v, axis, angle )
m_dir = rotateAroundAxis( m_dir, axis, angle )
generateMatrix()
end sub
function camera.projectionPlane() as rectangle
return( m_projPlane )
end function
function camera.getU() as vec4
return( m_u )
end function
function camera.getV() as vec4
return( m_v )
end function
function camera.getDir() as vec4
return( m_dir )
end function
function camera.getPos() as vec4
return( m_pos )
end function
function camera.getRatioUV() as double
return( m_u.length() / m_v.length() )
end function
function camera.getRatioVU() as double
return( m_v.length() / m_u.length() )
end function
sub camera.setRatioUV( byval ratio as double )
m_v.normalize()
m_v = m_v * m_u.length() / ratio
generateMatrix()
end sub
sub camera.setRatioVU(byval ratio as double)
m_u.normalize()
m_u = m_u * m_v.length() / ratio
generateMatrix()
end sub
function camera.getZoomU() as double
return( m_dir.length() / m_u.length() )
end function
function camera.getZoomV() as double
return( m_dir.length() / m_v.length() )
end function
sub camera.setZoomU( byval a as double )
m_u = m_u / ( a / getZoomU() )
generateMatrix()
end sub
sub camera.setZoomV( byval a as double )
m_v = m_v / ( a / getZoomV() )
generateMatrix()
end sub
sub camera.zoom( byval a as double )
zoomU( a )
zoomV( a )
end sub
sub camera.zoomU( byval a as double )
m_u = m_u / a
generateMatrix()
end sub
sub camera.zoomV( byval a as double )
m_v = m_v / a
generateMatrix()
end sub
function camera.getFOVV() as double
return( 2.0 * atan2( m_v.length(), m_dir.length() ) )
end function
function camera.getFOVU() as double
return( 2.0 * atan2( m_u.Length(), m_dir.length() ) )
end function
sub camera.setFOVU( byval angle as double )
setzoomU( 1.0 / tan( angle / 2.0 ) )
end sub
sub camera.setFOVV( byval angle as double )
setzoomV( 1.0 / tan( angle / 2.0 ) )
end sub
function camera.getYaw() as double
return( atan2( m_dir.x, m_dir.z ) )
end function
function camera.getPitch() as double
return( atan2( m_dir.y, sqr( m_dir.x * m_dir.x + m_dir.z * m_dir.z ) ) )
end function
function camera.getRoll() as double
return( vectorAngle( cross( vec4( 0.0, 1.0, 0.0 ), m_dir ), cross( m_v, m_dir ) ) )
end function
sub camera.setYaw( byval angle as double )
'' to change yaw, you have to rotate around the "up" axis of the WORLD = the y axis
dim currentAngle as double = getyaw()
rotate( vec4( 0.0, 1.0, 0.0 ), -angle + currentAngle )
end sub
sub camera.setPitch( byval angle as double )
dim currentAngle as double = getPitch()
rotate( m_u, angle - currentAngle )
end sub
sub camera.setRoll( byval angle as double )
dim currentAngle as double = getRoll()
rotate( m_dir, angle - currentAngle )
end sub
sub camera.yawAbsolute( byval angle as double )
rotate( vec4( 0.0, 1.0, 0.0 ), angle )
end sub
sub camera.yaw( byval angle as double )
rotate( m_v, angle)
end sub
sub camera.pitch(byval angle as double)
rotate( m_u, angle )
end sub
sub camera.roll( byval angle as double )
rotate( m_dir, angle )
end sub
function camera.getCameraMatrix() as mat4
return( m_camMatrix )
end function
function camera.getInverseCameraMatrix() as mat4
return( m_invCamMatrix )
end function
function camera.transform( byval v as vec4 ) as vec4
return( m_invCamMatrix * v )
end function
function camera.orthogonal( byval pnt as vec4, byval plane as rectangle, byref x as double, byref y as double, byref z as double ) as integer
dim as mat4 PM
dim as double ratio
if( plane.width < plane.height ) then
ratio = plane.height
else
ratio = plane.width
end if
with PM
.a = ( 2 / ratio ) * 7 : .b = 0.0 : .c = 0.0 : .d = 0.0
.e = 0.0 : .f = -( 2 / ratio ) * 7 : .g = 0.0 : .h = 0.0
.i = 0.0 : .j = 0.0 : .k = ( -2.0 / ( m_farClip - m_nearClip ) ) * 7 : .l = ( ( m_farClip + m_nearClip ) / ( m_farClip - m_nearClip ) )
.m = 0.0 : .n = 0.0 : .o = 0.0 : .p = 1.0
end with
Dim as vec4 b
b = PM * pnt
dim as double px, py
px = ( b.x * plane.width ) + plane.width / 2
py = ( b.y * plane.height ) + plane.height / 2
if( px >= 0 and px <= plane.width - 1 and py >= 0 and py <= plane.height - 1 ) then
x = px
y = py
z = b.z
return( -1 )
else
x = px
y = py
z = b.z
return( 0 )
end if
end function
function camera.perspective( byval pnt as vec4, byval plane as rectangle, byref x as double, byref y as double, byref z as double ) as integer
'' the current matrix
dim as double fovX = radians( 90 )
dim as double fovY = radians( 90 )
dim as double fovTanX = 1 / tan( fovX / 2 )
dim as double fovTanY = 1 / tan( fovY / 2 )
'' projection matrix (assumes a right-handed coord system)
dim as mat4 PM
with PM
.a = fovTanX : .b = 0.0 : .c = 0.0 : .d = 0.0
.e = 0.0 : .f = fovTanY : .g = 0.0 : .h = 0.0
.i = 0.0 : .j = 0.0 : .k = -( ( farClip + nearClip ) / ( farClip - nearClip ) ) : .l = -( ( 2 * ( farClip * nearClip ) ) / ( farClip - nearClip ) )
.m = 0.0 : .n = 0.0 : .o = -1.0 : .p = 0.0
end with
'' project the point
dim as vec4 b
b = PM * pnt
dim as double px
dim as double py
px = ( b.x * plane.width ) / ( 1.0 * b.w ) + plane.width / 2
py = ( b.y * plane.height ) / ( 1.0 * b.w ) + plane.height / 2
if( px >= 0 and px <= plane.width - 1 and py >= 0 and py <= plane.height - 1 and b.z > m_nearClip ) then
x = px
y = py
z = b.z
return( -1 )
else
x = px
y = py
z = b.z
return( 0 )
end if
end function
function camera.projectOnScreen( byval pnt as vec4, byval plane as rectangle, byref x as double, byref y as double, byref z as double ) as integer
'' first transformation: translation
dim as mat4 TM = matrices.translation( -m_pos )
dim as vec4 a = vec4( TM * pnt )
'' second transformation: rotation
'' transform() multiplies the vector with the INVERSE camera matrix
'' to bring the point from camera space to world space
dim as vec4 b = vec4( transform( a ) )
'' Third transformation: projection
return( perspective( b, plane, x, y, z ) )
end function
function camera.projectOnScreenOrtho( byval pnt as vec4, byval plane as rectangle, byref x as double, byref y as double, byref z as double ) as integer
dim as mat4 TM
'' first transformation: translation
with TM
.a = 1.0 : .b = 0.0 : .c = 0.0 : .d = -m_pos.x
.e = 0.0 : .f = 1.0 : .g = 0.0 : .h = -m_pos.y
.i = 0.0 : .j = 0.0 : .k = 1.0 : .l = -m_pos.z
.m = 0.0 : .n = 0.0 : .o = 0.0 : .p = 1.0
end with
dim as vec4 a = vec4( TM * pnt )
'' transform() multiplies the vector with the INVERSE camera matrix
'' to bring the point from camera space to world space
dim as vec4 b = vec4( transform( a ) )
'' third transformation: projection
return( orthogonal( b, plane, x, y, z ) )
end function
function camera.getScale() as double
return( m_dir.length )
end function
sub camera.setScale( byval dirLength as double )
scale( m_dir.length() / dirLength )
end sub
sub camera.scale( byval factor as double )
m_dir = m_dir * factor
m_u = m_u * factor
m_v = m_v * factor
generateMatrix()
end sub
sub camera.resetSkewU()
dim as double oldzoomU = getzoomU()
dim as double oldzoomV = getzoomV()
m_u = cross( m_dir, m_v )
m_v = cross( m_dir, -m_u )
setzoomU( oldzoomU )
setzoomV( oldzoomV )
end sub
sub camera.resetSkewV()
dim as double oldzoomU = getzoomU()
dim as double oldzoomV = getzoomV()
m_v = cross( m_dir, m_u )
m_u = cross( m_dir, -m_v )
setzoomU( oldzoomU )
setzoomV( oldzoomV )
end sub
sub camera.setLookDirection( byval newDir as vec4 )
dim as vec4 axis = vec4( ( cross( m_dir, newDir ) ) )
if axis.length() = 0 then
exit sub
end if
dim as double angle = vectorAngle( m_dir, newDir )
if angle <> 0 then
rotate( axis, angle )
end if
end sub
sub camera.lookAt( byval target as vec4, byval up as vec4 )
/'
preserve the lengths of the camera vectors
every time one messes with the camera vectors directly, it has to
remember to preserve their length, if you do not do this, the projection
gets distorted, as these vectors control different parameters of the
camera model
REMEMBER
u controls the horizontal FOV
v controls the vertical FOV
dir controls the focal zoom
'/
dim as double dl = m_dir.length()
dim as double ul = m_u.length()
dim as double vl = m_v.length()
m_dir = ( normalize( target - m_pos ) * dl )
m_U = ( normalize( cross( up, m_dir ) ) * ul )
m_V = ( normalize( cross( m_dir, m_u ) ) * vl )
generateMatrix()
end sub
sub camera.lookAt( byval target as vec4 )
/'
preserve the lengths of the camera vectors
every time one messes with the camera vectors directly, it has to
remember to preserve their length, if you do not do this, the projection
gets distorted, as these vectors control different parameters of the
camera model
REMEMBER
u controls the horizontal FOV
v controls the vertical FOV
dir controls the focal zoom
'/
dim as double dl = m_dir.length()
dim as double ul = m_u.length()
dim as double vl = m_v.length()
'' compute the forward vector
dim as vec4 forward = vec4( normalize( target - m_pos ) )
'' compute temporal up vector based on the forward vector
'' watch out when look up/down at 90 degree
'' for example, forward vector is on the Y axis
dim as vec4 up
if( abs( forward.x ) < epsilon and abs( forward.z ) < epsilon ) then
if( forward.y > 0 ) then
'' forward vector is pointing +Y axis
up = vec4( 0.0, 0.0, -1.0 )
else
'' forward vector is pointing -Y axis
up = vec4( 0.0, 0.0, 1.0 )
end if
else
'' in general, up vector is straight up
up = vec4( 0.0, 1.0, 0.0 )
end if
'' compute the left vector
dim as vec4 leftV = vec4( normalize( cross( up, forward ) ) )
'' re-calculate the orthonormal up vector
up = normalize( cross( forward, leftV ) )
m_dir = forward * dl
m_u = leftV * ul
m_v = up * vl
generateMatrix()
end sub
sub camera.setU( byval newU as vec4 )
m_u = newU
generateMatrix()
end sub
sub camera.setV( byval newV as vec4 )
m_v = newV
generateMatrix()
end sub
sub camera.setDir( byval newDir as vec4 )
m_dir = newDir
generateMatrix()
end sub
sub camera.setPos( byval newPos as vec4 )
m_pos = newPos
end sub
function camera.nearClip() as double
return( m_nearClip )
end function
function camera.farClip() as double
return( m_farClip )
end function
function camera.getDistance( byval pnt as vec4 ) as double
return( distance( m_pos, pnt ) )
end function
public sub camera.setDistance( byval pnt as vec4, byval dist as double )
dim currentDist as double = distance( m_pos, pnt )
move( currentDist * normalize( pnt - m_pos ) )
end sub
Code: Select all
'' needed headers
#include once "platform.bi"
#include once "core.bi"
#include once "math.bi"
#include once "vec4.bi"
#include once "mat4.bi"
#include once "arrayList.bi"
type edge
public:
declare constructor( byval v1 as vec4 ptr, byval v2 as vec4 ptr )
declare function vertex1() as vec4 ptr
declare function vertex2() as vec4 ptr
private:
as vec4 ptr m_v1
as vec4 ptr m_v2
end type
constructor edge( byval v1 as vec4 ptr, byval v2 as vec4 ptr )
m_v1 = v1
m_v2 = v2
end constructor
function edge.vertex1() as vec4 ptr
return( m_v1 )
end function
function edge.vertex2() as vec4 ptr
return( m_v2 )
end function
class object3D
public:
declare constructor()
declare constructor( _
byval position as vec4 = vec4( 0.0, 0.0, 0.0 ), _
byval xAxis as vec4, _
byval yAxis as vec4, _
byval zAxis as vec4 )
declare destructor()
declare function transform( byval v as vec4 ) as vec4
declare function objectSpaceToWorldSpace( byval pnt as vec4 ) as vec4
declare function getObjectMatrix() as mat4
declare function getInverseObjectMatrix() as mat4
declare sub move( byval offset as vec4 )
declare sub rotate( byval axis as vec4, byval angle as double )
declare function getyaw() as double
declare function getPitch() as double
declare function getRoll() as double
declare sub setyaw( byval angle as double )
declare sub setPitch( byval angle as double )
declare sub setRoll( byval angle as double )
declare sub yawAbsolute( byval angle as double )
declare sub yaw( byval angle as double )
declare sub pitch( byval angle as double )
declare sub roll( byval angle as double )
declare function getDistance( byval pnt as vec4 ) as double
declare sub setDistance( byval pnt as vec4, byval dist as double )
declare sub setLookDirection( byval newDir as vec4 )
declare sub lookAt( byval lookAtMe as vec4 )
declare sub setU( byval newU as vec4 )
declare sub setV( byval newV as vec4 )
declare sub setDir( byval newDir as vec4 )
declare sub setPos( byval newPos as vec4 )
declare function getU() as vec4
declare function getV() as vec4
declare function getDir() as vec4
declare function getPos() as vec4
'' these two functions are wrong on so many levels that isn't even funny
'' made it this way to simplify a little
declare function vertices() as arrayList ptr
declare function edges() as arrayList ptr
private:
declare sub generateMatrix()
as vec4 m_pos = any
as vec4 m_u = any
as vec4 m_v = any
as vec4 m_dir = any
as mat4 m_objMatrix = any
as mat4 m_invObjMatrix = any
as arrayList m_vertices
as arrayList m_edges
end class
constructor object3d()
'' construct an object with the default axes
m_pos = vec4( 0.0, 0.0, 0.0 )
m_u = vec4( 1.0, 0.0, 0.0 )
m_v = vec4( 0.0, 1.0, 0.0 )
m_dir = vec4( 0.0, 0.0, 1.0 )
generateMatrix()
end constructor
constructor object3d( _
byval position as vec4 = vec4( 0.0, 0.0, 0.0 ), _
byval xAxis as vec4, _
byval yAxis as vec4, _
byval zAxis as vec4 )
m_pos = position
m_u = xAxis
m_v = yAxis
m_dir = zAxis
generateMatrix()
end constructor
destructor object3d()
'' if there are vertices on the object, release them
dim as vec4 ptr v
if( m_vertices.count > 0 ) then
for i as integer = 0 to m_vertices.count - 1
v = m_vertices.get( i )
delete( v )
next
end if
'' same for the edges
dim as edge ptr e
if( m_edges.count > 0 ) then
for i as integer = 0 to m_edges.count - 1
e = m_edges.get( i )
delete( e )
next
end if
end destructor
function object3d.vertices() as arrayList ptr
return( @m_vertices )
end function
function object3d.edges() as arrayList ptr
return( @m_edges )
end function
function object3d.getObjectMatrix() as mat4
return( m_objMatrix )
end function
function object3d.getInverseObjectMatrix() as mat4
return( m_invObjMatrix )
end function
sub object3d.generateMatrix()
'' generates the object rotation matrix
m_objMatrix = mat4( _
m_u.x , m_v.x , m_dir.x , 0.0, _
m_u.y , m_v.y , m_dir.y , 0.0, _
m_u.z , m_v.z , m_dir.z , 0.0, _
m_u.w , m_v.w , m_dir.w , 1.0 )
'' and compute the inverse immediately
m_invObjMatrix = inverse( m_objMatrix )
end sub
sub object3d.move( byval offset as vec4 )
m_pos += offset
end sub
sub object3d.rotate( byval axis as vec4, byval angle as double )
'' rotates the object around an arbitrary axis
m_u = rotateAroundAxis( m_u, axis, angle)
m_v = rotateAroundAxis( m_v, axis, angle)
m_dir = rotateAroundAxis( m_dir, axis, angle)
generateMatrix()
end sub
function object3d.getU() as vec4
return( m_u )
end function
function object3d.getV() as vec4
return( m_v )
end function
function object3d.getDir() as vec4
return( m_dir )
end function
function object3d.getPos() as vec4
return( m_pos )
end function
function object3d.getyaw() as double
return( atan2( m_dir.x, m_dir.z ) )
end function
function object3d.getPitch() as double
return( atan2( m_dir.y, sqr( m_dir.x * m_dir.x + m_dir.z * m_dir.z ) ) )
end function
function object3d.getRoll() as double
return( vectorAngle( cross( vec4( 0.0, 1.0, 0.0 ), m_dir ), cross( m_v, m_dir ) ) )
end function
sub object3d.setyaw( byval angle as double )
'' to change yaw, you have to rotate around the "up" axis of the WORLD = the y axis
dim as double currentAngle = getyaw()
rotate( vec4( 0.0, 1.0, 0.0 ), -angle + currentAngle )
end sub
sub object3d.setPitch( byval angle as double )
dim as double currentAngle = getPitch()
rotate( m_u, angle - currentAngle )
end sub
sub object3d.setRoll( byval angle as double )
dim as double currentAngle = getRoll()
rotate( m_dir, angle - currentAngle )
end sub
sub object3d.yawAbsolute( byval angle as double )
rotate( vec4( 0.0, 1.0, 0.0 ), angle )
end sub
sub object3d.yaw( byval angle as double )
rotate( m_v, angle )
end sub
sub object3d.pitch( byval angle as double )
rotate( m_u, angle )
end sub
sub object3d.roll( byval angle as double )
rotate( m_dir, angle )
end sub
function object3d.getDistance( byval pnt as vec4 ) as double
return( distance( m_pos, pnt ) )
end function
sub object3d.setDistance( byval pnt as vec4, byval dist as double )
dim as double currentDist = distance( m_pos, pnt )
move( currentDist * normalize( pnt - pos ) )
end sub
sub object3d.setLookDirection( byval newDir as vec4 )
dim axis as vec4 = vec4( ( cross( m_dir, newDir ) ) )
if axis.length() = 0 then
exit sub
end if
dim as double angle = vectorAngle( m_dir, newDir )
if( angle <> 0 ) then
rotate( axis, angle )
end if
end sub
sub object3d.lookAt( byval lookAtMe as vec4 )
dim as double dl = m_dir.length()
dim as double ul = m_u.length()
dim as double vl = m_v.length()
dim as vec4 w = vec4( lookAtMe - m_pos )
setDir( Normalize( w ) * dl )
setU( normalize( cross( vec4( 0.0, 1.0, 0.0 ), m_dir ) ) * ul )
setV( normalize( cross( m_dir, m_u ) ) * vl )
end sub
sub object3d.setU( byval newU as vec4 )
m_u = newU
generateMatrix()
end sub
sub object3d.setV( byval newV as vec4 )
m_v = newV
generateMatrix()
end sub
sub object3d.setDir( byval newDir as vec4 )
m_dir = newDir
generateMatrix()
end sub
sub object3d.setPos( byval newPos as vec4 )
m_pos = newPos
end sub
function object3d.transform( byval v as vec4 ) as vec4
return( m_objMatrix * v )
end function
function object3d.objectSpaceToWorldSpace( byval pnt as vec4 ) as vec4
dim as vec4 b = vec4( transform( pnt ) )
dim as vec4 a = vec4( b + m_pos )
return( a )
end function
UPDATE 1: 10/11/2017 - camera.bas: killed a few bugses and sprinkled some more comments
Last edited by paul doe on Oct 13, 2017 10:27, edited 2 times in total.
Re: 3d without openGL
vec4.bi
mat4.bi
utility.bi
Code: Select all
#ifndef __vec4__
#define __vec4__
#include once "math.bi"
/'
homogeneous 4 tuple vector type
it is meant to be used as an homogeneous 3D vector (like the ones
used by OpenGL or Direct3D)
conceptually, they are to be interpreted like this:
| x |
| y |
| z |
| 1 |
which is known as an homogeneous column vector
10/7/2017: fixed a typo in the + operator
now this was the original code. Take a very close look at it (the comment was removed
after the fix out of pure anger, as you can imagine):
operator +( byref v as vec4, byref w as vec4 ) as vec4
'' addition operator (all rocket science stuff...)
return( vec4( v.x + w.x, v.y + w.y, v.x + w.z ) )
end operator
goes to show you what happens when you are a #%$@ imbecile and write code in a hurry.
STAY TUNED FOR THE NExT ExCITING TUTORIAL: how to write unbelievable $%#@ 3D code!!!
'/
type vec4
public:
as double x
as double y
as double z
as double w
declare constructor( byval nx as double = 0.0, byval ny as double = 0.0, byval nz as double = 0.0, byval nw as double = 1.0 )
declare constructor( byref rhs as vec4 )
declare operator let( byref rhs as vec4 )
declare operator cast() as string
declare function length() as double
declare function squaredLength() as double
declare sub normalize()
declare function normal() as vec4
declare sub homogeneize()
declare function homogeneous() as vec4
declare function cross( byref v as vec4 ) as vec4
declare function dot( byref v as vec4 ) as double
declare function distance( byref v as vec4 ) as double
end type
constructor vec4( byval nx as double = 0.0, byval ny as double = 0.0, byval nz as double = 0.0, byval nw as double = 1.0 )
'' default constructor creates an homogeneous vector
x = nx
y = ny
z = nz
w = nw
end constructor
constructor vec4( byref rhs as vec4 )
'' copy constructor
x = rhs.x
y = rhs.y
z = rhs.z
w = rhs.w
end constructor
operator vec4.let( byref rhs as vec4 )
'' assignment constructor
x = rhs.x
y = rhs.y
z = rhs.z
w = rhs.w
end operator
operator vec4.cast() as string
'' human readable string representation (useful for debugging)
return( _
"| " & trim( str( x ) ) & " |" & chr( 13 ) & chr( 10 ) & _
"| " & trim( str( y ) ) & " |" & chr( 13 ) & chr( 10 ) & _
"| " & trim( str( z ) ) & " |" & chr( 13 ) & chr( 10 ) & _
"| " & trim( str( w ) ) & " |" & chr( 13 ) & chr( 10 ) )
end operator
function vec4.squaredLength() as double
/'
returns the squared length of this vector
useful when you just want to compare which one is bigger, as
this avoids having to compute a square root
'/
return( x * x + y * y + z * z )
end function
function vec4.length() as double
'' returns the length of this vector
return( sqr( x * x + y * y + z * z ) )
end function
sub vec4.normalize()
/'
normalizes the vector
note that the homogeneous coordinate (w) is not touched
'/
dim as double l = length()
if l > 0 then
x /= l : y /= l : z /= l
end if
end sub
function normalize( byref v as vec4 ) as vec4
'' for compatibility
return( v.normal() )
end function
function vec4.normal() as vec4
/'
returns this vector normalized but without altering itself
again the homogeneous coordinate is left alone
'/
dim as double l = length()
dim as vec4 v = vec4( this )
if l > 0 then
v.x /= l : v.y /= l : v.z /= l
end if
return( v )
end function
sub vec4.homogeneize()
/'
homogeneizes the vector
this is done by dividing the components by the homogeneous coordinate (w)
'/
x /= w : y /= w : z /= w : w /= w
end sub
function vec4.homogeneous() as vec4
/'
returns this vector homogeneized but without altering it
'/
return( vec4( x / w, y / w, z / w, w / w ) )
end function
function vec4.cross( byref v as vec4 ) as vec4
/'
returns the cross product (aka vectorial product) of this vector and
another vector v
note that there's no cross product defined for a 4 tuple, so
we simply use a standard 3d cross product, and make the w component 1
'/
return( vec4( _
v.y * z - v.z * y, _
v.z * x - v.x * z, _
v.x * y - v.y * x, _
1.0 ) )
end function
function cross( byref v as vec4, byref w as vec4 ) as vec4
'' returns the cross product between vectors v and w
return( v.cross( w ) )
end function
function vec4.dot( byref v as vec4 ) as double
'' returns the dot product (aka scalar product) of this vector and vector v
return( v.x * x + v.y * y + v.z * z )
end function
function dot( byref v as vec4, byref w as vec4 ) as double
'' returns the dot product between two vectors
return( v.dot( w ) )
end function
function vec4.distance( byref v as vec4 ) as double
/'
gets the distance of this vector with vector v
to calculate the distance, substract them and calculate the length of the resultant vector
'/
return( sqr( ( v.x - x ) ^ 2 + ( v.y - y ) ^ 2 + ( v.z - z ) ^ 2 ) )
end function
function distance( byref v as vec4, byref w as vec4 ) as double
return( sqr( ( v.x - w.x ) ^ 2 + ( v.y - w.y ) ^ 2 + ( v.z - w.z ) ^ 2 ) )
end function
operator -( byref v as vec4, byref w as vec4 ) as vec4
'' substraction operator
return( vec4( v.x - w.x, v.y - w.y, v.z - w.z ) )
end operator
operator -( byref v as vec4 ) as vec4
'' negation operator
return( vec4( -v.x, -v.y, -v.z ) )
end operator
operator +( byref v as vec4, byref w as vec4 ) as vec4
return( vec4( v.x + w.x, v.y + w.y, v.z + w.z ) )
end operator
operator *( byref v as vec4, byref a as double ) as vec4
/'
multiplication with a scalar
note that this is not a 'scalar product', but a multiplication with a number
vectors do not define multiplications per se, they define the dot product
and the cross product. To avoid confusion (and also correctness), the
multiplication operator is overloaded to a scaling of the vector
'/
return( vec4( v.x * a, v.y * a, v.z * a ) )
end operator
operator *( byref a as double, byref v as vec4 ) as vec4
'' same as above but with the parameters inversed
return( vec4( v.x * a, v.y * a, v.z * a ) )
end operator
operator /( byref v as vec4, byref a as double ) as vec4
'' division by a scalar. See note above on multiplying a vector
return( vec4( v.x / a, v.y / a, v.z / a ) )
end operator
function vectorAngle( byref v as vec4, byref w as vec4 ) as double
/'
returns the angle between two vectors using the dot product, in radians
note that the result of the dot product used here
should mathematically speaking, always give a result between -1 and 1. Due to imprecisions of
numerical calculations it might sometimes be a little bit outside this range however (especially
if you define Scalar to be float instead of double). If that happens, the acos function will give
an invalid result. So instead a protection was added that sets the value back to 1 or -1
(because, if the value became 1.0000023 for example, it was probably supposed to be 1 anyway).
dotProduct( v, w ) = length( v ) * length( w ) * cos( angle )
angle = aCos( dotProduct / ( length( v ) * length( w ) ) )
thus:
angle = aCos( dot( normal( v ) * normal( w ) ) )
'/
dim as double angleCos = dot( v.normal(), w.normal() )
'' for acos, the value has to be between -1.0 and 1.0, but due to numerical imprecisions it
'' sometimes comes outside this range
angleCos = clamp( -1.0, 1.0, angleCos )
return( -acos( angleCos ) )
end function
function rotateAroundAxis( byref v as vec4, byref axis as vec4, byval angle as double ) as vec4
/'
rotate vector v around arbitrary axis for angle radians
it can only rotate around an axis through our object, to rotate around another axis:
first translate the object to the axis, then use this function, then translate back
in the new direction.
'/
if( ( v.x = 0 ) and ( v.y = 0 ) and ( v.z = 0 ) ) then
return vec4( 0.0, 0.0, 0.0 )
end if
dim nAxis as vec4 = vec4( axis.x, axis.y, axis.z )
nAxis.normalize()
'' calculate parameters of the rotation matrix
dim as double c = cos( angle )
dim as double s = sin( angle )
dim as double t = 1 - c
'' multiply w with rotation matrix
dim w as vec4
w.x = ( t * nAxis.x * nAxis.x + c ) * v.x _
+ ( t * nAxis.x * nAxis.y + s * nAxis.z ) * v.y _
+ ( t * nAxis.x * nAxis.z - s * nAxis.y ) * v.z
w.y = ( t * nAxis.x * nAxis.y - s * nAxis.z ) * v.x _
+ ( t * nAxis.y * nAxis.y + c ) * v.y _
+ ( t * nAxis.y * nAxis.z + s * nAxis.x ) * v.z
w.z = ( t * nAxis.x * nAxis.z + s * nAxis.y ) * v.x _
+ ( t * nAxis.y * nAxis.z - s * nAxis.x ) * v.y _
+ ( t * nAxis.z * nAxis.z + c ) * v.z
'' the vector has to retain its length, so it's normalized and
'' multiplied with the original length
w.normalize()
w = w * v.length()
return( w )
end function
#endif
Code: Select all
#ifndef __mat4__
#define __mat4__
#include once "core.bi"
#include once "platform.bi"
#include once "vec4.bi"
/'
| a b c d |
4x4 matrix type | e f g h |
| i j k l |
| m n o p |
9/30/2017: improved 4x4 matrix inverse calculation. With -gen gcc -O max
(the settings I always use) it is more than 60% faster than the previous
version. A byproduct of correcting the bug mentioned below.
9/29/2017: fixed determinant calculation (it was erroneously computed)
The funny thing was that this library is used in various applications,
including my 3D engine and various tools. When I was implementing a
feature, it kept doing weird things, and the bug was finally tracked to
an erroneous calculation of the determinant. The determinant was
correctly calculated in the 3x3 matrix code, but when I ported the code
to use OpenGL for rendering, it didn't worked as intended. Goes to show
you that one is to be extra careful with the math code, as it is the
foundation of the entire game engine.
the calculations were cross checked with the help of this online resource:
https://www.mathsisfun.com/algebra/matrix-calculator.html
which has a very neat matrix calculator for various dimensions.
'/
type mat4
public:
as double a, b, c, d
as double e, f, g, h
as double i, j, k, l
as double m, n, o, p
declare constructor( _
byval sa as double = 1.0, byval sb as double = 0.0, byval sc as double = 0.0, byval sd as double = 0.0, _
byval se as double = 0.0, byval sf as double = 1.0, byval sg as double = 0.0, byval sh as double = 0.0, _
byval si as double = 0.0, byval sj as double = 0.0, byval sk as double = 1.0, byval sl as double = 0.0, _
byval sm as double = 0.0, byval sn as double = 0.0, byval so as double = 0.0, byval sp as double = 1.0 _
)
declare constructor( byref NM as mat4 )
declare operator let( byref RHS as mat4 )
declare operator cast() as string
declare function determinant() as double
declare sub transpose()
declare sub inverse() '' slower version
'declare sub inverse()
declare sub identity()
end type
constructor mat4( _
byval sa as double = 1.0, byval sb as double = 0.0, byval sc as double = 0.0, byval sd as double = 0.0, _
byval se as double = 0.0, byval sf as double = 1.0, byval sg as double = 0.0, byval sh as double = 0.0, _
byval si as double = 0.0, byval sj as double = 0.0, byval sk as double = 1.0, byval sl as double = 0.0, _
byval sm as double = 0.0, byval sn as double = 0.0, byval so as double = 0.0, byval sp as double = 1.0 _
)
/'
default constructor initializes the matrix to an identity, if no coefficients are specified.
this is far more useful than initializing it to all zeros
'/
a = sa : b = sb : c = sc : d = sd
e = se : f = sf : g = sg : h = sh
i = si : j = sj : k = sk : l = sl
m = sm : n = sn : o = so : p = sp
end constructor
constructor mat4( byref RHS as mat4 )
'' copy constructor
a = RHS.a : b = RHS.b : c = RHS.c : d = RHS.d
e = RHS.e : f = RHS.f : g = RHS.g : h = RHS.h
i = RHS.i : j = RHS.j : k = RHS.k : l = RHS.l
m = RHS.m : n = RHS.n : o = RHS.o : p = RHS.p
end constructor
operator mat4.let( byref RHS as mat4 )
'' assignment construction
a = RHS.a : b = RHS.b : c = RHS.c : d = RHS.d
e = RHS.e : f = RHS.f : g = RHS.g : h = RHS.h
i = RHS.i : j = RHS.j : k = RHS.k : l = RHS.l
m = RHS.m : n = RHS.n : o = RHS.o : p = RHS.p
end operator
operator mat4.cast() as string
'' the matrix in a human readable form (very useful for debugging purposes)
return( _
"| " & trim( str( a ) ) & " | " & trim( str( b ) ) & " | " & trim( str( c ) ) & " | " & trim( str( d ) ) & " |" & chr( 13 ) & chr( 10 ) & _
"| " & trim( str( e ) ) & " | " & trim( str( f ) ) & " | " & trim( str( g ) ) & " | " & trim( str( h ) ) & " |" & chr( 13 ) & chr( 10 ) & _
"| " & trim( str( i ) ) & " | " & trim( str( j ) ) & " | " & trim( str( k ) ) & " | " & trim( str( l ) ) & " |" & chr( 13 ) & chr( 10 ) & _
"| " & trim( str( m ) ) & " | " & trim( str( n ) ) & " | " & trim( str( o ) ) & " | " & trim( str( p ) ) & " |" & chr( 13 ) & chr( 10 ) _
)
end operator
function mat4.determinant() as double
/'
computes the determinant of the matrix using Laplace cofactor expansion
the determinant of a 3x3 matrix is:
a * ( e * i - f * h ) - b * ( d * i - f * g ) + c * ( d * h - e * g )
and a 4x4 matrix determinant is given by:
a * b * c * d *
| f g h | | e g h | | e f h | | e f g |
| j k l | - | i k l | + | i j l | - | i j k |
| n o p | | m o p | | m n p | | m n o |
where the '|' means the determinant of the inner 3x3 matrices. Note that the
cofactors are already factored in the calculation.
the determinant is thus:
+ ( a * ( f * ( k * p - l * o ) - g * ( j * p - l * n ) + h * ( j * o - k * n ) ) )
- ( b * ( e * ( k * p - l * o ) - g * ( i * p - l * m ) + h * ( i * o - k * m ) ) )
+ ( c * ( e * ( j * p - l * n ) - f * ( i * p - l * m ) + h * ( i * n - j * m ) ) )
- ( d * ( e * ( j * o - k * n ) - f * ( i * o - k * m ) + g * ( i * n - j * m ) ) )
'/
dim as double det = _
( a * ( f * ( k * p - l * o ) - g * ( j * p - l * n ) + h * ( j * o - k * n ) ) ) _
- ( b * ( e * ( k * p - l * o ) - g * ( i * p - l * m ) + h * ( i * o - k * m ) ) ) _
+ ( c * ( e * ( j * p - l * n ) - f * ( i * p - l * m ) + h * ( i * n - j * m ) ) ) _
- ( d * ( e * ( j * o - k * n ) - f * ( i * o - k * m ) + g * ( i * n - j * m ) ) )
/'
this isn't matematically correct, just a programmer's dirty hack.
if the determinant of a matrix is 0, it means it has no inverse. In the code for
calculating the inverse, a division by the determinant is performed; and if it is
zero, a division by zero is performed on *every* element of the matrix, filling it
with positive or negative infinity values and rendering it useless. A matrix
without inverse is the matrix itself, so setting the determinant value to 1
does the trick.
'/
if det = 0 then det = 1.0
return( det )
end function
sub mat4.transpose()
/'
transposes the matrix
[ a b c d ]T [ a e i m ]
[ e f g h ] = [ b f j n ]
[ i j k l ] [ c g k o ]
[ m n o p ] [ d h l p ]
why have it, if it is not used by the matrix code itself? Well, there is a
nice property of matrices, which has to do with rotations. If you can
be sure that the matrix contains only rotations, transposing it is the
same as taking its inverse, thus saving you *a lot* of computation
'/
this = mat4( a, e, i, m, b, f, j, n, c, g, k, o, d, h, l, p )
end sub
/'
sub mat4.inverse()
/'
computes the inverse of a 4x4 matrix
this version is 60%+ faster and 400%+ uglier than the previous version.
it was made so by computing the determinant inside the method and
recycling as much calculation as possible
'/
'' list of 2x2 determinants
dim as double kplo = k * p - l * o
dim as double jpln = j * p - l * n
dim as double jokn = j * o - k * n
dim as double iplm = i * p - l * m
dim as double iokm = i * o - k * m
dim as double injm = i * n - j * m
dim as double gpho = g * p - h * o
dim as double fphn = f * p - h * n
dim as double fogn = f * o - g * n
dim as double ephm = e * p - h * m
dim as double eogm = e * o - g * m
dim as double enfm = e * n - f * m
dim as double glhk = g * l - h * k
dim as double flhj = f * l - h * j
dim as double fkgj = f * k - g * j
dim as double elhi = e * l - h * i
dim as double ekgi = e * k - g * i
dim as double ejfi = e * j - f * i
'' list of 3x3 determinants
dim as double d1kplo = f * kplo
dim as double d1jpln = g * jpln
dim as double d1jokn = h * jokn
dim as double d2kplo = e * kplo
dim as double d2iplm = g * iplm
dim as double d2iokm = h * iokm
dim as double d3jpln = e * jpln
dim as double d3iplm = f * iplm
dim as double d3injm = h * injm
dim as double d4jokn = e * jokn
dim as double d4iokm = f * iokm
dim as double d4injm = g * injm
dim as double d5kplo = b * kplo
dim as double d5jpln = c * jpln
dim as double d5jokn = d * jokn
dim as double d6kplo = a * kplo
dim as double d6iplm = c * iplm
dim as double d6iokm = d * iokm
dim as double d7jpln = a * jpln
dim as double d7iplm = b * iplm
dim as double d7injm = d * injm
dim as double d8jokn = a * jokn
dim as double d8iokm = b * iokm
dim as double d8injm = c * injm
dim as double d9gpho = b * gpho
dim as double d9fphn = c * fphn
dim as double d9fogn = d * fogn
dim as double d10gpho = a * gpho
dim as double d10ephm = c * ephm
dim as double d10eogm = d * eogm
dim as double d11fphn = a * fphn
dim as double d11ephm = b * ephm
dim as double d11enfm = d * enfm
dim as double d12fogn = a * fogn
dim as double d12eogm = b * eogm
dim as double d12enfm = c * enfm
dim as double d13glhk = b * glhk
dim as double d13flhj = c * flhj
dim as double d13fkgj = d * fkgj
dim as double d14glhk = a * glhk
dim as double d14elhi = c * elhi
dim as double d14ekgi = d * ekgi
dim as double d15flhj = a * flhj
dim as double d15elhi = b * elhi
dim as double d15ejfi = d * ejfi
dim as double d16fkgj = a * fkgj
dim as double d16ekgi = b * ekgi
dim as double d16ejfi = c * ejfi
'' 4x4 determinant (inversed)
dim as double det = _
( a * ( d1kplo - d1jpln + d1jokn ) _
- ( b * ( d2kplo - d2iplm + d2iokm ) ) _
+ ( c * ( d3jpln - d3iplm + d3injm ) ) _
- ( d * ( d4jokn - d4iokm + d4injm ) ) )
'' if the determinant is 0, the matrix has no inverse
if det = 0 then exit sub
'' Multiplying with the reciprocal is slightly faster than dividing
dim as double invDet = 1.0 / det
'' minors
dim as double Ma = d1kplo - d1jpln + d1jokn
dim as double Mb = d2kplo - d2iplm + d2iokm
dim as double Mc = d3jpln - d3iplm + d3injm
dim as double Md = d4jokn - d4iokm + d4injm
dim as double Me = d5kplo - d5jpln + d5jokn
dim as double Mf = d6kplo - d6iplm + d6iokm
dim as double Mg = d7jpln - d7iplm + d7injm
dim as double Mh = d8jokn - d8iokm + d8injm
dim as double Mi = d9gpho - d9fphn + d9fogn
dim as double Mj = d10gpho - d10ephm + d10eogm
dim as double Mk = d11fphn - d11ephm + d11enfm
dim as double Ml = d12fogn - d12eogm + d12enfm
dim as double Mm = d13glhk - d13flhj + d13fkgj
dim as double Mn = d14glhk - d14elhi + d14ekgi
dim as double Mo = d15flhj - d15elhi + d15ejfi
dim as double Mp = d16fkgj - d16ekgi + d16ejfi
/'
adjugate (the adjugate is the transpose of the cofactored matrix of minors)
Ma -Me Mi -Mm
-Mb Mf -Mj Mn
Mc -Mg Mk -Mo
-Md Mh -Ml Mp
'/
this = mat4( _
Ma * invDet, -Me * invDet, Mi * invDet, -Mm * invDet, _
-Mb * invDet, Mf * invDet, -Mj * invDet, Mn * invDet, _
Mc * invDet, -Mg * invDet, Mk * invDet, -Mo * invDet, _
-Md * invDet, Mh * invDet, -Ml * invDet, Mp * invDet )
end sub
'/
sub mat4.inverse()
/'
calculates the inverse of a 4x4 matrix
the minor, cofactor and adjugate matrices are all calculated and factored into the
code to make it shorter and more efficient
'/
'' multiplying by the reciprocal of the determinant is faster than dividing by it
dim as double invDet = 1 / determinant
this = mat4( _
( f * k * p + g * l * n + h * j * o - f * l * o - g * j * p - h * k * n ) * invDet, _
( b * l * o + c * j * p + d * k * n - b * k * p - c * l * n - d * j * o ) * invDet, _
( b * g * p + c * h * n + d * f * o - b * h * o - c * f * p - d * g * n ) * invDet, _
( b * h * k + c * f * l + d * g * j - b * g * l - c * h * j - d * f * k ) * invDet, _
( e * l * o + g * i * p + h * k * m - e * k * p - g * l * m - h * i * o ) * invDet, _
( a * k * p + c * l * m + d * i * o - a * l * o - c * i * p - d * k * m ) * invDet, _
( a * h * o + c * e * p + d * g * m - a * g * p - c * h * m - d * e * o ) * invDet, _
( a * g * l + c * h * i + d * e * k - a * h * k - c * e * l - d * g * i ) * invDet, _
( e * j * p + f * l * m + h * i * n - e * l * n - f * i * p - h * j * m ) * invDet, _
( a * l * n + b * i * p + d * j * m - a * j * p - b * l * m - d * i * n ) * invDet, _
( a * f * p + b * h * m + d * e * n - a * h * n - b * e * p - d * f * m ) * invDet, _
( a * h * j + b * e * l + d * f * i - a * f * l - b * h * i - d * e * j ) * invDet, _
( e * k * n + f * i * o + g * j * m - e * j * o - f * k * m - g * i * n ) * invDet, _
( a * j * o + b * k * m + c * i * n - a * k * n - b * i * o - c * j * m ) * invDet, _
( a * g * n + b * e * o + c * f * m - a * f * o - b * g * m - c * e * n ) * invDet, _
( a * f * k + b * g * i + c * e * j - a * g * j - b * e * k - c * f * i ) * invDet _
)
'' lovely, isn't it?
end sub
sub mat4.identity()
'' makes the matrix an identity matrix
a = 1.0 : b = 0.0 : c = 0.0 : d = 0.0
e = 0.0 : f = 1.0 : g = 0.0 : h = 0.0
i = 0.0 : j = 0.0 : k = 1.0 : l = 0.0
m = 0.0 : n = 0.0 : o = 0.0 : p = 1.0
end sub
operator *( byref A as mat4, byref B as mat4 ) as mat4
/'
multiply two 4x4 matrices
remember that matrix multiplication is not commutative!
A * B != B * A
'/
return( mat4( _
A.a * B.a + A.b * B.e + A.c * B.i + A.d * B.m, _
A.a * B.b + A.b * B.f + A.c * B.j + A.d * B.n, _
A.a * B.c + A.b * B.g + A.c * B.k + A.d * B.o, _
A.a * B.d + A.b * B.h + A.c * B.l + A.d * B.p, _
A.e * B.a + A.f * B.e + A.g * B.i + A.h * B.m, _
A.e * B.b + A.f * B.f + A.g * B.j + A.h * B.n, _
A.e * B.c + A.f * B.g + A.g * B.k + A.h * B.o, _
A.e * B.d + A.f * B.h + A.g * B.l + A.h * B.p, _
A.i * B.a + A.j * B.e + A.k * B.i + A.l * B.m, _
A.i * B.b + A.j * B.f + A.k * B.j + A.l * B.n, _
A.i * B.c + A.j * B.g + A.k * B.k + A.l * B.o, _
A.i * B.d + A.j * B.h + A.k * B.l + A.l * B.p, _
A.m * B.a + A.n * B.e + A.o * B.i + A.p * B.m, _
A.m * B.b + A.n * B.f + A.o * B.j + A.p * B.n, _
A.m * B.c + A.n * B.g + A.o * B.k + A.p * B.o, _
A.m * B.d + A.n * B.h + A.o * B.l + A.p * B.p _
) )
'return( mat4( _
' A.a * B.a + A.b * B.e + A.c * B.i + A.d * B.m, _
' A.a * B.b + A.b * B.f + A.c * B.j + A.d * B.n, _
' A.a * B.c + A.b * B.g + A.c * B.k + A.d * B.o, _
' A.a * B.d + A.b * B.h + A.c * B.l + A.d * B.p, _
' A.e * B.a + A.f * B.e + A.g * B.i + A.h * B.m, _
' A.e * B.b + A.f * B.f + A.g * B.j + A.h * B.n, _
' A.e * B.c + A.f * B.g + A.g * B.k + A.h * B.o, _
' A.e * B.d + A.f * B.h + A.g * B.l + A.h * B.p, _
' A.i * B.a + A.j * B.e + A.k * B.i + A.l * B.m, _
' A.i * B.b + A.j * B.f + A.k * B.j + A.l * B.n, _
' A.i * B.c + A.j * B.g + A.k * B.k + A.l * B.o, _
' A.i * B.d + A.j * B.h + A.k * B.l + A.l * B.p, _
' A.m * B.a + A.n * B.e + A.o * B.i + A.p * B.m, _
' A.m * B.b + A.n * B.f + A.o * B.j + A.p * B.n, _
' A.m * B.c + A.n * B.g + A.o * B.k + A.p * B.o, _
' A.m * B.d + A.n * B.h + A.o * B.l + A.p * B.p _
') )
end operator
operator +( byref A as mat4, byref B as mat4 ) as mat4
'' adds two 4x4 matrices
return( mat4( _
A.a + B.a, A.b + B.b, A.c + B.c, A.d + B.d, _
A.e + B.e, A.f + B.f, A.g + B.g, A.h + B.h, _
A.i + B.i, A.j + B.j, A.k + B.k, A.l + B.l, _
A.m + B.m, A.n + B.n, A.o + B.o, A.p + B.p _
) )
end operator
operator -( byref A as mat4, byref B as mat4 ) as mat4
'' substracts two 4x4 matrices
return( mat4( _
A.a - B.a, A.b - B.b, A.c - B.c, A.d - B.d, _
A.e - B.e, A.f - B.f, A.g - B.g, A.h - B.h, _
A.i - B.i, A.j - B.j, A.k - B.k, A.l - B.l, _
A.m - B.m, A.n - B.n, A.o - B.o, A.p - B.p _
) )
end operator
operator -( byref A as mat4 ) as mat4
'' negates the matrix
return( mat4( _
-A.a, -A.b, -A.c, -A.d, _
-A.e, -A.f, -A.g, -A.h, _
-A.i, -A.j, -A.k, -A.l, _
-A.m, -A.n, -A.o, -A.p _
) )
end operator
operator *( byref A as mat4, byref s as double ) as mat4
'' scalar multiplication
return( mat4( _
A.a * s, A.b * s, A.c * s, A.d * s, _
A.e * s, A.f * s, A.g * s, A.h * s, _
A.i * s, A.j * s, A.k * s, A.l * s, _
A.m * s, A.n * s, A.o * s, A.p * s _
) )
end operator
operator *( byref s as double, byref A as mat4 ) as mat4
'' scalar multiplication
return( mat4( _
A.a * s, A.b * s, A.c * s, A.d * s, _
A.e * s, A.f * s, A.g * s, A.h * s, _
A.i * s, A.j * s, A.k * s, A.l * s, _
A.m * s, A.n * s, A.o * s, A.p * s _
) )
end operator
operator /( byref s as double, byref A as mat4 ) as mat4
/'
scalar divide
the 'division' of a matrix by another matrix is not defined, the
equivalent operation is multiplying one matrix by the inverse of the other
on scalars though, it can be defined, mostly for convenience purposes
'/
return( mat4( _
A.a / s, A.b / s, A.c / s, A.d / s, _
A.e / s, A.f / s, A.g / s, A.h / s, _
A.i / s, A.j / s, A.k / s, A.l / s, _
A.m / s, A.n / s, A.o / s, A.p / s _
) )
end operator
operator /( byref A as mat4, byref s as double ) as mat4
return( mat4( _
A.a / s, A.b / s, A.c / s, A.d / s, _
A.e / s, A.f / s, A.g / s, A.h / s, _
A.i / s, A.j / s, A.k / s, A.l / s, _
A.m / s, A.n / s, A.o / s, A.p / s _
) )
end operator
operator *( byref v as vec4, byref A as mat4 ) as vec4
/'
multiply a vector with a row matrix, resulting in a row vector (like Direct3D)
a row vector looks like this:
| x y z w |
and is the format that Direct3D uses. What this means, code-wise, is that you
have to pre-multiply the vectors with the matrices, and some other stuff, like
transposing the matrices if you are using column vectors (as this library does)
'/
return( vec4( _
A.a * v.x + A.e * v.y + A.i * v.z + A.m * v.w, _
A.b * v.x + A.f * v.y + A.j * v.z + A.n * v.w, _
A.c * v.x + A.g * v.y + A.k * v.z + A.o * v.w, _
A.d * v.x + A.h * v.y + A.l * v.z + A.p * v.w _
) )
end operator
operator *( byref A as mat4, byref v as vec4 ) as vec4
/'
multiply a vector with a column matrix, resulting in a column vector (like OpenGL)
a column vector looks like this
| x |
| y |
| z |
| w |
and is the format favored by OpenGL. In this library, column vectors are used, for
compatibility.
'/
'w.X = A.a * v.X + A.b * v.Y + A.c * v.Z + A.d * v.W
'w.Y = A.e * v.X + A.f * v.Y + A.g * v.Z + A.h * v.W
'w.Z = A.i * v.X + A.j * v.Y + A.k * v.Z + A.l * v.W
'w.W = A.m * v.X + A.n * v.Y + A.o * v.Z + A.p * v.W
return( vec4( _
A.a * v.X + A.b * v.Y + A.c * v.Z + A.d * v.W, _
A.e * v.X + A.f * v.Y + A.g * v.Z + A.h * v.W, _
A.i * v.X + A.j * v.Y + A.k * v.Z + A.l * v.W, _
A.m * v.X + A.n * v.Y + A.o * v.Z + A.p * v.W _
) )
end operator
'' utility functions
function inverse( byref M as mat4 ) as mat4
'' returns the inverse of the provided matrix
dim as mat4 I = mat4( M )
I.inverse()
return( I )
end function
namespace matrices
function translation( byval t as vec4 ) as mat4
return mat4( _
1.0, 0.0, 0.0, t.x, _
0.0, 1.0, 0.0, t.y, _
0.0, 0.0, 1.0, t.z, _
0.0, 0.0, 0.0, 1.0 )
end function
function identity() as mat4
'' returns an identity matrix
return( mat4( _
1.0, 0.0, 0.0, 0.0, _
0.0, 1.0, 0.0, 0.0, _
0.0, 0.0, 1.0, 0.0, _
0.0, 0.0, 0.0, 1.0 ) )
end function
end namespace
#endif
Code: Select all
/'
auxiliary (read: rushed) functions for performing some drawings on the
3D plane
'/
function clipSegment( byref cam as camera, byref s1 as vec4, byref s2 as vec4 ) as integer
/'
clip both vertices of a line if necessary
note that this 'clipping' is only done on the near plane, not on the scene
the sensible (and correct) way to do this would be to clip them against
the view frustum; but alas, in order to keep the code short and clear
the calculation of the frustum was removed from the camera class. Hence
this crappy function.
Clipping (especially against the near plane) is important, because if you
don't clip, when one of the line endpoints goes behind the near clipping
plane, it gets a negative coordinate (due to the projection), and the line
segment is no longer valid
'/
if( s1.z >= cam.nearClip() and s2.z >= cam.nearClip() ) then
' Neither one is behind the camera, draw them both
return( -1 )
elseIf( s1.z < cam.nearClip() and s2.z >= cam.nearClip() ) then
' First coordinate behind the camera, clip it
s1.x = s1.x + ((cam.nearClip() - s1.z) / (s2.z - s1.z)) * (s2.x - s1.x)
s1.y = s1.y + ((cam.nearClip() - s1.z) / (s2.z - s1.z)) * (s2.y - s1.y)
s1.z = cam.nearClip()
return( -1 )
elseIf( s1.z >= cam.nearClip() and s2.z < cam.nearClip() ) then
' Second coordinate behind the camera, clip it
s2.x = s1.x + ((cam.nearClip() - s1.z) / (s2.z - s1.z)) * (s2.x - s1.x)
s2.y = s1.y + ((cam.nearClip() - s1.z) / (s2.z - s1.z)) * (s2.y - s1.y)
s2.z = cam.nearClip()
return( -1 )
else
' Both coordinates behind the camera, don't draw
return( 0 )
end if
end function
sub drawLine3d( byref cam as camera, byref p1 As vec4, byref p2 As vec4, byval c as uint32 )
dim as vec4 pos1, pos2
pos1 = cam.transform( p1 - cam.getPos() )
pos2 = cam.transform( p2 - cam.getPos() )
dim as double x1, y1, x2, y2
dim as integer visible = clipSegment( cam, pos1, pos2 )
if( visible ) then
dim as double z1, z2
cam.perspective( pos1, cam.projectionPlane, x1, y1, z1 )
cam.perspective( pos2, cam.projectionPlane, x2, y2, z2 )
line( x1, y1 ) - ( x2, y2 ), c
end if
end sub
Re: 3d without openGL
So, you will need to have, in the same folder:
platform.bi
core.bi
math.bi
vec4.bi
mat4.bi
arrayList.bi
utility.bi
camera.bas
object.bas
Once you've got all that assembled, copy and paste this demo code, in the same folder:
camera test.bas
And there you have it, your own 3D playground. Read the comments in the code to see how it works. Should you have any questions, or need some explanation on what the code does, I'm all ears ;-)
Have fun!
Paul
platform.bi
core.bi
math.bi
vec4.bi
mat4.bi
arrayList.bi
utility.bi
camera.bas
object.bas
Once you've got all that assembled, copy and paste this demo code, in the same folder:
camera test.bas
Code: Select all
#include once "fbgfx.bi"
#include once "platform.bi" '' platform specific definitions
#include once "core.bi" '' core functions and defines
#include once "math.bi" '' math functions and constants
#include once "vec4.bi" '' vec4 type and operations
#include once "mat4.bi" '' mat4 type and operations
#include once "arrayList.bi" '' convenience
#include once "camera.bas" '' the camera class definition
#include once "object.bas" '' the object definition
#include once "utility.bi" '' to draw lines an' stuff directly in 3D space
/'
3D playground
intended as a simple tutorial/framework to do 3D stuff without too much complication
it is a bare-bones implementation, actually. Mostly to test all the math/conventions
behind the representation of a 3D scene (using the term 'scene' veeeery loosely here)
conventions used
RIGHT HANDED coordinate system
positive rotations are COUNTER-CLOCKWISE
facets are defined in COUNTER-CLOCKWISE order
angles are in RADIANS - use radians( angle_in_degrees ) for convenience
'/
'' set a screen mode
const as integer scrW = 1200, scrH = 600
screenRes( scrW, scrH, 32, , fb.gfx_alpha_primitives )
windowTitle( fbVersion & " - 3D Playground" )
/'
define a projection plane (in this case the screen itself)
the projection plane is where the image gets projected (duh)
in the pinhole camera model, the projection plane is actually the
so called 'near' plane. The near parameter of the camera is really
the distance of this plane to the origin of coordinates (in the
case of the camera, its position in the world coordinate system
'/
dim as rectangle projectionPlane = ( 0.0, 0.0, scrW, scrH )
/'
instantiate a camera class
remember the parameters for the constructor:
position, in WORLD space
X axis of the camera's coordinate system
Y axis of the camera's coordinate system
Z axis of the camera's coordinate system. If you flip it, you will end with a
left-handed coordinate system. Not that it matters much, save for the fact that
all other methods treat the matrix as a right-handed one. If you use lookAt()
and suddenly the axes get flipped, this is the most probable reason
near clipping plane Z distance
far clipping plane Z distance
projection plane
a word of advice: the axes of the camera have to be of a certain length, for the
projection matrix depends on it. See the code below, where the U and V vectors
get scaled to match the relation of the projection plane. If you fool around with
them and suddenly the image looks like crap, this is the most probable reason.
to remember:
the U vector (the X one) controls the horizontal field of view
the V vector (the Y one) controls the vertical field of view
and the Z vector controls the zoom
of course, they must be perfectly perpendicular to each other, if not the resulting
image is sheared. Anyway, you don't have to mess with the constructor if you don't
like and/or understand what it does for the camera to do its job.
having said all that, play around! See how each function transforms the camera, and
the effect it has on the image. Look at the implementation to see how it's done, but
most important of all, have fun!
'/
dim as camera cam = camera( _
vec4( 0.0, 1.0, 3.0 ), _
vec4( 1.0, 0.0, 0.0 ), _
vec4( 0.0, 1.0, 0.0 ), _
vec4( 0.0, 0.0, -1.0 ), _
0.1, 10000.0, projectionPlane )
/'
scene manager (yeah right)
here, it's just a list that contains the objects to be rendered, but in real code it
is an instance of the scene manager.
'/
dim as arrayList objects
/'
let's create some objects, shall we? One, actually (yeah I'm cheap)
these are the vertices and edges to define a small paper plane that
you can manipulate (see the code in the main loop)
the objects ideally should be loaded from a file, but that will
complicate the implementation. If there's some interest, I will show
how you can load an obj file (or invent your own format if you wish)
'/
dim as object3d ptr obj = new object3d() '' it will be at the origin, at the start
obj->vertices->add( new vec4( -0.1, 0.0, 0.0 ) )
obj->vertices->add( new vec4( 0.1, 0.0, 0.0 ) )
obj->vertices->add( new vec4( 0.0, 0.0, 0.3 ) )
obj->vertices->add( new vec4( 0.0, 0.0, 0.0 ) )
obj->vertices->add( new vec4( 0.0, -0.1, 0.0 ) )
obj->vertices->add( new vec4( 0.0, 0.0, 0.3 ) )
obj->edges->add( new edge( obj->vertices->get( 0 ), obj->vertices->get( 1 ) ) )
obj->edges->add( new edge( obj->vertices->get( 1 ), obj->vertices->get( 2 ) ) )
obj->edges->add( new edge( obj->vertices->get( 2 ), obj->vertices->get( 0 ) ) )
obj->edges->add( new edge( obj->vertices->get( 3 ), obj->vertices->get( 4 ) ) )
obj->edges->add( new edge( obj->vertices->get( 4 ), obj->vertices->get( 5 ) ) )
obj->edges->add( new edge( obj->vertices->get( 5 ), obj->vertices->get( 3 ) ) )
'' add the defined object to the scene
objects.add( obj )
'' move it a little so it doesn't lie exactly in the center of the world
obj->move( vec4( 1.0, 0.0, -1.0 ) )
'' some variables used for interaction
dim as integer mouseX, mouseY, mouseButton
dim as integer oldMouseX, oldMouseY
'' to make movement (somewhat) constant
dim as double frameTime
'' background color (change it if you don't like it)
dim as uint32 backgroundColor = rgba( 8, 8, 8, 255 )
color( , backgroundColor )
'' (extremely crappy) main loop
do
'' grab an object to have fun
dim as object3d ptr o = objects.get( 0 )
'' update interaction
oldMouseX = mouseX
oldMouseY = mouseY
getMouse( mouseX, mouseY, , mouseButton )
/'
camera controls
click and draw the mouse to free look
w: forward
s: backward
a: strafe left
d: strafe right
q: up
e: down
space: look at the paper plane (hold to keep looking at it)
you see, they are very similar to that of a FPS
'/
if( multikey( fb.sc_e ) ) then
'' down
cam.move( -vec4( 0.0, 1.0, 0.0 ) * 20.0 * frameTime )
end if
if( multikey( fb.sc_q ) ) then
'' up
cam.move( vec4( 0.0, 1.0, 0.0 ) * 20.0 * frameTime )
end if
if( multikey( fb.sc_w ) ) then
'' forward
cam.move( normalize( cam.getDir() ) * 20.0 * frameTime )
end if
if( multikey( fb.sc_s ) ) then
'' backwards
cam.move( normalize( -cam.getDir() ) * 20.0 * frameTime )
end if
if( multikey( fb.sc_a ) ) then
'' left
cam.move( normalize( cam.getU() ) * 20.0 * frameTime )
end if
if( multikey( fb.sc_d ) ) then
'' right
cam.move( normalize( -cam.getU() ) * 20.0 * frameTime )
end if
if( multikey( fb.sc_space ) ) then
'' look at the object (in case we lost it)
cam.lookAt( o->getPos() )
end if
/'
paper plane controls
i: forward
k: backwards
j: turn left
l: turn right
u: turn up
o: turn down
'/
if( multikey( fb.sc_u ) ) then
o->rotate( normalize( o->getU() ), radians( 0.5 ) )
end if
if( multikey( fb.sc_o ) ) then
o->rotate( normalize( o->getU() ), radians( -0.5 ) )
end if
if( multikey( fb.sc_i ) ) then
o->move( normalize( o->getDir() ) * 10.0 * frameTime )
end if
if( multikey( fb.sc_k ) ) then
o->move( normalize( -( o->getDir() ) ) * 10.0 * frameTime )
end if
if( multikey( fb.sc_j ) ) then
o->rotate( normalize( o->getV() ), radians( 0.5 ) )
end if
if( multikey( fb.sc_l ) ) then
o->rotate( normalize( o->getV() ), radians( -0.5 ) )
end if
'' if the left mouse button is pressed, activate free look mode
if( mouseButton = 1 ) then
cam.rotate( vec4( 0.0, 1.0, 0.0 ), 5.0 * ( oldMouseX - mouseX ) / cam.projectionPlane.width )
cam.rotate( cam.getU, -5.0 * ( oldMouseY - mouseY ) / cam.projectionPlane.height )
end if
'' render the screen
frameTime = timer()
screenLock()
cls()
'' draw the floor to have a frame of reference
for zz as integer = -10 to 10 step 1
drawLine3d( cam, vec4( -10, 0, zz ), vec4( 10, 0, zz ), rgba( 32, 32, 32, 255 ) )
next
for xx as integer = -10 to 10 step 1
drawLine3d( cam, vec4( xx, 0, -10 ), vec4( xx, 0, 10 ), rgba( 32, 32, 32, 255 ) )
next
/'
draw the absolute axes of the world
x = red
y = green
z = blue
'/
drawLine3d( cam, vec4( 0.0, 0.0, 0.0 ), vec4( 3.0, 0.0, 0.0 ), rgba( 128, 0, 0, 255 ) )
drawLine3d( cam, vec4( 0.0, 0.0, 0.0 ), vec4( 0.0, 3.0, 0.0 ), rgba( 0, 128, 0, 255 ) )
drawLine3d( cam, vec4( 0.0, 0.0, 0.0 ), vec4( 0.0, 0.0, 3.0 ), rgba( 0, 0, 128, 255 ) )
'' renders the objects
for i as integer = 0 to objects.count - 1
dim as object3d ptr o
o = objects.get( i )
'' render the edges
for j as integer = 0 to o->edges->count - 1
dim as edge ptr e = o->edges->get( j )
drawLine3d( cam, o->objectSpaceToWorldSpace( *e->vertex1 ), o->objectSpaceToWorldSpace( *e->vertex2 ), rgba( 255, 255, 0, 255 ) )
next
next
screenUnlock()
/'
it's important not to let this var to be negative, as it gets multiplied with
the camera vectors and could screw some calculations (movement, for example)
told you, the loop implementation is very crappy, but it seems like this is
esier for most people to grasp
'/
frameTime = abs( timer() - frameTime )
sleep( 1 )
loop until multikey( fb.sc_escape ) '' loop until esc is pressed
'' finally, if there's objects in the list, free them
if( objects.count > 0 ) then
dim as object3d ptr o
for i as integer = 0 to objects.count - 1
o = objects.get( i )
delete( o )
next
end if
Have fun!
Paul
-
- Posts: 3906
- Joined: Jan 01, 2009 7:03
- Location: Australia
Re: 3d without openGL
You get an A+ on the code. It does exactly what I would have liked to have done.paul doe wrote:And there you have it, your own 3D playground. Read the comments in the code to see how it works. Should you have any questions, or need some explanation on what the code does, I'm all ears ;-)
It is unlikely I will ever understand it but I wonder to what extent I could use the code myself?
I imagine it could be extended to have a world of filled polygons to do something like the old F/A-18 Interceptor Amiga game.
https://www.youtube.com/watch?v=oHsUbGTIXyE
.
Re: 3d without openGL
It was obvious from your first demo 'main' that you've got a 'complete solution' - based on my own experience one typically needs a 'camera,' an object list, or 'world' or 'object manager' - something many of us haven't got into yet. Your next demo 'grid test,' combined with what it does, and how 'trivially' it does it is, well, unprecedented.paul doe wrote:Hi dafhidafhi wrote:@paul doe - Very intrigued and inspired by your offerings
Glad to hear that. Which one in particular sparkled your interest?
A lot of us are eager to step into 3d programming and you've dropped a deluxe-model scene composition lab on our doorstep
Last edited by dafhi on Oct 11, 2017 18:00, edited 1 time in total.
Re: 3d without openGL
Basiccoder2
In my rotate the first parameter is the fulcrum.
The second parameter is the point to rotate.
The third parameter is the angle which holds all six trigs necessary :
Sin(x),Sin(y),Sin(z)
Cos(x),Cos(y),Cos(z)
The fourth parameter is the scaler which is (1,1,1) default (no scaling)
You have set up your points around (0,0,0), so this is the fulcrum.
Also the eye point for perspective will be .x=0 and .y=0
I set the eye point .z to 600 to show a little more perspective.
In the final showing of points I add 640 to the x's and 300 to the y's so the shape is in the screen centre.
You don't really need the initial rotate by pi/2.
(This was to get the stones around Saturn into a Saturn disc.)
You have only orthogonal axis to show.
You can apply your .5 scale during the rotate in the loop.
Your code:
Others:
Matrices are too slow.
My rotate does the matrix * vector multiplication straight off in three lines.
For axial rotation Rodrigues method can also be used (I have posted many examples).
But it involves normalising and cross products which slows the whole process down.
Although it is easier in a way, you just stipulate an axis (2,-6,8) for example.
In my rotate the first parameter is the fulcrum.
The second parameter is the point to rotate.
The third parameter is the angle which holds all six trigs necessary :
Sin(x),Sin(y),Sin(z)
Cos(x),Cos(y),Cos(z)
The fourth parameter is the scaler which is (1,1,1) default (no scaling)
You have set up your points around (0,0,0), so this is the fulcrum.
Also the eye point for perspective will be .x=0 and .y=0
I set the eye point .z to 600 to show a little more perspective.
In the final showing of points I add 640 to the x's and 300 to the y's so the shape is in the screen centre.
You don't really need the initial rotate by pi/2.
(This was to get the stones around Saturn into a Saturn disc.)
You have only orthogonal axis to show.
You can apply your .5 scale during the rotate in the loop.
Your code:
Code: Select all
Type V3
As Single x,y,z
As Uinteger col
End Type
Type angle3D 'FLOATS for angles
As Single sx,sy,sz
As Single cx,cy,cz
Declare static Function construct(As Single,As Single,As Single) As Angle3D
End Type
Type float As V3
#define incircle(cx,cy,radius,x,y) (cx-x)*(cx-x) +(cy-y)*(cy-y)<= radius*radius
Function Angle3D.construct(x As Single,y As Single,z As Single) As Angle3D
Return Type (Sin(x),Sin(y),Sin(z), _
Cos(x),Cos(y),Cos(z))
End Function
Function Rotate(c As V3,p As V3,a As Angle3D,scale As float=type(1,1,1)) As V3
Dim As Single dx=p.x-c.x,dy=p.y-c.y,dz=p.z-c.z
Return Type<V3>((scale.x)*((a.cy*a.cz)*dx+(-a.cx*a.sz+a.sx*a.sy*a.cz)*dy+(a.sx*a.sz+a.cx*a.sy*a.cz)*dz)+c.x,_
(scale.y)*((a.cy*a.sz)*dx+(a.cx*a.cz+a.sx*a.sy*a.sz)*dy+(-a.sx*a.cz+a.cx*a.sy*a.sz)*dz)+c.y,_
(scale.z)*((-a.sy)*dx+(a.sx*a.cy)*dy+(a.cx*a.cy)*dz)+c.z,p.col)
End Function
Function perspective(p As V3,eyepoint As V3) As V3
Dim As Single w=1+(p.z/eyepoint.z)
Return Type<V3>((p.x-eyepoint.x)/w+eyepoint.x,_
(p.y-eyepoint.y)/w+eyepoint.y,_
(p.z-eyepoint.z)/w+eyepoint.z,p.col)
End Function
Sub QsortZ(array() As V3,begin As Long,Finish As Ulong)
Dim As Long i=begin,j=finish
Dim As V3 x =array(((I+J)\2))
While I <= J
While array(I).z > X .z:I+=1:Wend
While array(J).z < X .z:J-=1:Wend
If I<=J Then Swap array(I),array(J): I+=1:J-=1
Wend
If J >begin Then QsortZ(array(),begin,J)
If I <Finish Then QsortZ(array(),I,Finish)
End Sub
Function Regulate(Byval MyFps As long,Byref fps As long) As long
Static As Double timervalue,_lastsleeptime,t3,frames
frames+=1
If (Timer-t3)>=1 Then t3=Timer:fps=frames:frames=0
Var sleeptime=_lastsleeptime+((1/myfps)-Timer+timervalue)*1000
If sleeptime<1 Then sleeptime=1
_lastsleeptime=sleeptime
timervalue=Timer
Return sleeptime
End Function
screenres 1280,600,32
const pi=4*atn(1)
redim as v3 a(0)
'set some points
for n as single =-250 to 250
var xp=0,yp=0,zp = n
var u=ubound(a)
redim preserve a(1 to u+1)
a(ubound(a))=type<v3>(xp, yp, zp, rnd*rgb(0,0,255))
pset(xp,yp),a(ubound(a)).col
next n
for n as single =-250 to 250
var xp=0,yp=n,zp = 0
var u=ubound(a)
redim preserve a(1 to u+1)
a(ubound(a))=type<v3>(xp, yp, zp, rnd*rgb(255,0,0))
pset(xp,yp),a(ubound(a)).col
next n
for n as single =-250 to 250
var xp=n,yp=0,zp = 0
var u=ubound(a)
redim preserve a(1 to u+1)
a(ubound(a))=type<v3>(xp, yp, zp, rnd*rgb(0,255,0))
pset(xp,yp),a(ubound(a)).col
next n
'print "the points are on the x y plane with a small random z component"
'print "flip them all 90 degrees around y axis and expand a bit"
'print "Press a key"
'sleep
'cls
redim as V3 rot(lbound(a) to ubound(a))'to hold all rotated points
dim as Angle3D A3d'=angle3D.construct(0,pi/2,0) 'flip all points by pi/2 on y axis
'for n as long=lbound(a) to ubound(a)
'rot(n)=rotate(type(0,0,0),a(n),A3D,type(0.5,0.5,0.5)) '1.4 is a scale up
'a(n)=rot(n)
'pset(a(n).x,a(n).y),a(n).col
'next
'print "Thats 'em all flipped without perspctive"
'print "Press a key"
'sleep
dim as float ang
dim as long fps
'fix these two components at start
ang.z=pi/2 'tilt angle sideways
ang.y=-pi/7 'tilt angle back and forward
dim as string key
do
key=inkey
If key=Chr(255)+"K" Then ang.z-=.02 'left
If key=Chr(255)+"M" Then ang.z+=.02 'right
If key=Chr(255)+"P" Then ang.y-=.02 'down
If key=Chr(255)+"H" Then ang.y+=.02 'up
if key="z" then ang.x -=.02
if key="x" then ang.x +=.02
A3D = Angle3D.construct(ang.x,ang.y,ang.z)'get the six rotate components (sines, coses ...)
screenlock
cls
draw string(50,50),"Framerate "&fps
draw string(50,150),"Use the arrow keys"
for n as long=lbound(a) to ubound(a)
rot(n)=rotate(type(0,0,0),a(n),A3D,type(.5,.5,.5))
rot(n)=perspective(rot(n),type(0,0,600))
next
qsortZ(rot(),lbound(rot),ubound(rot))
for n as long=lbound(rot) to ubound(rot)
pset(rot(n).x+640,rot(n).y+300),rot(n).col
next
screenunlock
sleep regulate(64,fps),1
loop until key=chr(27)
Matrices are too slow.
My rotate does the matrix * vector multiplication straight off in three lines.
For axial rotation Rodrigues method can also be used (I have posted many examples).
But it involves normalising and cross products which slows the whole process down.
Although it is easier in a way, you just stipulate an axis (2,-6,8) for example.
-
- Posts: 3906
- Joined: Jan 01, 2009 7:03
- Location: Australia
Re: 3d without openGL
Maybe axial rotation was the term I was looking for? That the controls would rotate the points around the axes of the object not the axes of the world? I would not have recognized at the time which examples of yours does that. I have collected all your demos but of course none of them have user control or any tutorial explanations of how the equations were derived.dodicat wrote:For axial rotation Rodrigues method can also be used (I have posted many examples).
With paul doe's paper plane the keys U and O rotate it around one unmoving axis and the same with the J and L keys. In both cases one axis remains fixed while the other two rotate around it. He uses the I and K to move along the the third axis which is what an plane would do but I would also want to be able to rotate around that axis as well move along the direction of an axis.
By the way with your example I just noticed the rnd*rgb(255,0,0) should have been plain rgb(255,0,0) to get the three pure colors for each axis instead of random values meant for your Saturn ring demo.
.
Last edited by BasicCoder2 on Oct 11, 2017 18:07, edited 1 time in total.
Re: 3d without openGL
Hi basiccoder2.
Here is an example with a grid:
(I have drawn the axis of rotation)
I am not quite sure what you seek.
I presume a shape and an axis to rotate the shape around.
My first shape was Saturn's rings.
My second shape was your x/y/z axis.
My third shape (here) is a grid, using only the end points.
Rodrigues rotation is more complex.
Here is an example, I did it a long while back .
My equations for the ordinary rotations (not rodrigue) were derived from the 3D rotation matrix.
But Matrix multiplication by looping is far too slow, and the matrix is only 3 X 3 anyway, so it can be written straight down.
Thus the contorted looking return in the rotator.
Here is an example with a grid:
(I have drawn the axis of rotation)
Code: Select all
Screen 19,32
Const pi=4*Atn(1)
Type V3
As Single x,y,z
As Ulong col
End Type
Type angle3D 'FLOATS for angles
As Single sx,sy,sz
As Single cx,cy,cz
Declare Static Function construct(As Single,As Single,As Single) As Angle3D
End Type
Type float As V3
Function Angle3D.construct(x As Single,y As Single,z As Single) As Angle3D
Return Type (Sin(x),Sin(y),Sin(z), _
Cos(x),Cos(y),Cos(z))
End Function
Function Rotate(c As V3,p As V3,a As Angle3D,scale As float=Type(1,1,1)) As V3
Dim As Single dx=p.x-c.x,dy=p.y-c.y,dz=p.z-c.z
Return Type<V3>((scale.x)*((a.cy*a.cz)*dx+(-a.cx*a.sz+a.sx*a.sy*a.cz)*dy+(a.sx*a.sz+a.cx*a.sy*a.cz)*dz)+c.x,_
(scale.y)*((a.cy*a.sz)*dx+(a.cx*a.cz+a.sx*a.sy*a.sz)*dy+(-a.sx*a.cz+a.cx*a.sy*a.sz)*dz)+c.y,_
(scale.z)*((-a.sy)*dx+(a.sx*a.cy)*dy+(a.cx*a.cy)*dz)+c.z,p.col)
End Function
Function perspective(p As V3,eyepoint As V3) As V3
Dim As Single w=1+(p.z/eyepoint.z)
Return Type<V3>((p.x-eyepoint.x)/w+eyepoint.x,_
(p.y-eyepoint.y)/w+eyepoint.y,_
(p.z-eyepoint.z)/w+eyepoint.z,p.col)
End Function
Sub QsortZ(array() As V3,begin As Long,Finish As Ulong)'NOT USED HERE.
Dim As Long i=begin,j=finish
Dim As V3 x =array(((I+J)\2))
While I <= J
While array(I).z > X .z:I+=1:Wend
While array(J).z < X .z:J-=1:Wend
If I<=J Then Swap array(I),array(J): I+=1:J-=1
Wend
If J >begin Then QsortZ(array(),begin,J)
If I <Finish Then QsortZ(array(),I,Finish)
End Sub
Function Regulate(Byval MyFps As Long,Byref fps As Long) As Long
Static As Double timervalue,_lastsleeptime,t3,frames
frames+=1
If (Timer-t3)>=1 Then t3=Timer:fps=frames:frames=0
Var sleeptime=_lastsleeptime+((1/myfps)-Timer+timervalue)*1000
If sleeptime<1 Then sleeptime=1
_lastsleeptime=sleeptime
timervalue=Timer
Return sleeptime
End Function
Sub setgrid(sx As Long,bx As Long,sy As Long,by As Long,st As Long,p() As v3)
#define U Ubound(p)
Redim p(0)
For y As Long=sy To by Step st
Redim Preserve p(1 To U+1)
p(U)=Type(sx,y,0)
Redim Preserve p(1 To U+1)
p(U)=Type(bx,y,0)
Line(sx,y)-(bx,y)
Next
For x As Long=sx To bx Step st
Redim Preserve p(1 To U+1)
p(U)=Type(x,sy,0)
Redim Preserve p(1 To U+1)
p(U)=Type(x,by,0)
Line(x,sy)-(x,by)
Next
End Sub
Sub drawgrid(p() As V3)
Dim As Ulong clr=Rgb(0,100,255)
For n As Long=Lbound(p) To Ubound(p)-1
Line (p(n).x,p(n).y)-(p(n+1).x,p(n+1).y),clr
If n=Ubound(p)-3 Then clr=Rgb(0,200,0)'the axis green
n+=1
Next
End Sub
Redim As v3 a()'the main array of 3D points
Print "The initial grid"
setgrid(200,600,100,500,40,a())
'get dead centre of grid -- fulcrum
Dim As Single cx,cy,cz
For n As Long=Lbound(a) To Ubound(a)
cx+=a(n).x
cy+=a(n).y
cz+=a(n).z
Next
Dim As Long sz=(Ubound(a) - Lbound(a) +1)
cx=cx/sz
cy=cy/sz
cz=cz/sz
Print "Grid fulcrum:"
Print cx,cy,cz
'put in an axis (normal to the grid)
Redim Preserve a(Lbound(a) To (Ubound(a)+2))
a(Ubound(a)-1)=Type(cx,cy,0)
a(Ubound(a)) =Type(cx,cy,-200)
Print "Now rotate pi/2 round y axis"
Print "Press a key"
Sleep
Cls
Redim As V3 rot(Lbound(a) To Ubound(a))'to hold all rotated points
Dim As Angle3D A3d=angle3D.construct(0,pi/2,0) 'flip all points by pi/2 on y axis
For n As Long=Lbound(a) To Ubound(a)
rot(n)=rotate(Type(cx,cy,cz),a(n),A3D,Type(1.4,1.4,1.4)) '1.4 is a scale up
a(n)=rot(n)
Pset(a(n).x,a(n).y),a(n).col
Next
Print "Thats the grid rotated, with axis shown"
Print "Press a key"
drawgrid(a())
Sleep
Dim As float ang
Dim As Long fps
'fix these two components at start
ang.z=pi/2 'tilt angle sideways
ang.y=-pi/7 'tilt angle back and forward
Dim As String key
Do
key=Inkey
If key=Chr(255)+"K" Then ang.z-=.02 'left
If key=Chr(255)+"M" Then ang.z+=.02 'right
If key=Chr(255)+"P" Then ang.y-=.02 'down
If key=Chr(255)+"H" Then ang.y+=.02 'up
ang.x+=.01 'the orbiting speed
A3D=Angle3D.construct(ang.x,ang.y,ang.z)'get the six rotate components (sines, coses ...)
Screenlock
Cls
Draw String(50,50),"Framerate "&fps
Draw String(50,150),"Use the arrow keys"
For n As Long=Lbound(a) To Ubound(a)
rot(n)=rotate(Type(cx,cy,cz),a(n),A3D,Type(1,1,1))
rot(n)=perspective(rot(n),Type(cx,cy,1000))
Next
''qsortZ(rot(),lbound(rot),ubound(rot)) 'definately do not sort
drawgrid(rot())
Screenunlock
Sleep regulate(64,fps),1
Loop Until key=Chr(27)
Sleep
I presume a shape and an axis to rotate the shape around.
My first shape was Saturn's rings.
My second shape was your x/y/z axis.
My third shape (here) is a grid, using only the end points.
Rodrigues rotation is more complex.
Here is an example, I did it a long while back .
Code: Select all
'SetEnviron("fbgfx=GDI")
#include "fbgfx.bi"
Using fb
Type v3
As Single x,y,z
Declare Property length As Single
Declare Property unit As v3
Declare Function AxialRotate(As v3,As Single,As v3) As v3
Declare Function PointRotate(As v3,As v3,As v3=Type<v3>(1,1,1)) As v3
Declare Function perspective(eyepoint As v3) As v3
#define vct Type<v3>
#define dot *
#define cross ^
End Type
Type line3d
As v3 v1,v2
Declare property length As Single
End Type
#define map(a,b,x,c,d) ((d)-(c))*((x)-(a))/((b)-(a))+(c)
#define incircle(cx,cy,radius,x,y) (cx-x)*(cx-x) +(cy-y)*(cy-y)<= radius*radius
'degrees to radians
#define r .01745329
#define onscreen (mx>0) and (mx<xres) and (my>0) and (my<yres)
Operator + (v1 As v3,v2 As v3) As v3
Return Type<v3>(v1.x+v2.x,v1.y+v2.y,v1.z+v2.z)
End Operator
Operator -(v1 As v3,v2 As v3) As v3
Return Type<v3>(v1.x-v2.x,v1.y-v2.y,v1.z-v2.z)
End Operator
Operator * (f As Single,v1 As v3) As v3
Return vct(f*v1.x,f*v1.y,f*v1.z)
End Operator
Operator *(v1 As v3,f As Single) As v3
Return f*v1
End Operator
Operator * (v1 As v3,v2 As v3) As Single 'dot product
Return v1.x*v2.x+v1.y*v2.y+v1.z*v2.z
End Operator
Operator ^ (v1 As v3,v2 As v3) As v3 'cross product
Return Type<v3>(v1.y*v2.z-v2.y*v1.z,-(v1.x*v2.z-v2.x*v1.z),v1.x*v2.y-v2.x*v1.y)
End Operator
Property v3.length As Single
Return Sqr(this.x*this.x+this.y*this.y+this.z*this.z)
End Property
property line3d.length As Single
Dim As v3 p=this.v1-this.v2
Return p.length
End property
Property v3.unit As v3
Dim n As Single=this.length
If n=0 Then n=1e-20
Return vct(this.x/n,this.y/n,this.z/n)
End Property
Function v3.AxialRotate(centre As v3,Angle As Single,axis As v3) As v3
Dim As v3 v=This-centre
Dim As v3 norm=axis.unit
Return (V*Cos(Angle)+(Norm cross V)*Sin(Angle)+Norm*(Norm dot V)*(1-Cos(Angle)))+centre
End Function
Function v3.perspective(eyepoint As v3) As v3
Dim As Single w=1+(this.z/eyepoint.z)
If w=0 Then w=1e-20
Return vct((this.x-eyepoint.x)/w+eyepoint.x,(this.y-eyepoint.y)/w+eyepoint.y,(this.z-eyepoint.z)/w+eyepoint.z)
End Function
Function length(v As v3) As Single
Return Sqr(v.x*v.x+v.y*v.y+v.z*v.z)
End Function
Function v3.PointRotate(_point As v3,Angle As v3,scale As v3=Type<v3>(1,1,1)) As v3
Dim As v3 p=vct(This-_point)
Dim As v3 rot,temp
Dim As Single s=Sin(angle.x),c=Cos(angle.x)
temp=vct((p.y)*C+(-p.z)*S,(p.z)*C+(p.y)*S)
rot.y=temp.x
s=Sin(angle.y):c=Cos(angle.y)
temp=vct((temp.y)*C+(-p.x)*S,(p.x)*C+(temp.y)*S)
rot.z=temp.x
s=Sin(angle.z):c=Cos(angle.z)
temp=vct((temp.y)*C+(-rot.y)*S,(rot.y)*C+(temp.y)*S)
rot.x=temp.x:rot.y=temp.y
Return vct((scale.x*rot.x+_point.x),(scale.y*rot.y+_point.y),(scale.z*rot.z+_point.z))
End Function
Function segment_distance(l As line3d,p As v3,Byref ip As v3=vct(0,0,0)) As Single
Dim As Single linelength=length(l.v1-l.v2)
Dim As Single dist= length( (1/linelength)*((l.v1-l.v2) cross (p-l.v1)))
If length(p-l.v1) >= length(p-l.v2) Then
var leg=(p-l.v1)
var part=Sqr((length(leg))^2-dist^2)
var temp=part/linelength
'If temp>=1 Then temp=1:dist=length(p-l.v2)
ip=l.v1+(temp)*(l.v2-l.v1)
Return dist
Else
var leg=(p-l.v2)
var part=Sqr((length(leg))^2-dist^2)
var temp=part/linelength
'If temp>=1 Then temp=1:dist=length(p-l.v1)
ip=l.v2+(temp)*(l.v1-l.v2)
Return dist
End If
Return dist
End Function
Function Regulate(Byval MyFps As Long,Byref fps As Long) As Long
Static As Double timervalue,_lastsleeptime,t3,frames
frames+=1
If (Timer-t3)>=1 Then t3=Timer:fps=frames:frames=0
Var sleeptime=_lastsleeptime+((1/myfps)-Timer+timervalue)*1000
If sleeptime<1 Then sleeptime=1
_lastsleeptime=sleeptime
timervalue=Timer
Return sleeptime
End Function
'=================================================================================
Function keys(a As v3) As Integer
If Multikey(SC_LEFT) Then a.y=a.y+.005':return 1
If Multikey(SC_RIGHT) Then a.y=a.y-.005':return 1
If Multikey(SC_UP ) Then a.x=a.x-.005':return 1
If Multikey(SC_DOWN ) Then a.x=a.x+.005':return 1
If Multikey(SC_Q) Then a.z=a.z-.005':return 1
If Multikey(SC_W) Then a.z=a.z+.005':return 1
If Multikey(SC_SPACE) Then Return 2
Return 1
End Function
#macro display()
ans=pt.AxialRotate(centre,angle*r,normal_line)
Dim As v3 pp
var dist=segment_distance(L1,ans,pp)
ans=ans.perspective(vct(xres/2,yres/2,600))
angle=angle+2
If angle >=360 Then angle=0
Screenlock
Cls
Circle (pt.x,pt.y),5,5,,,,f
pp=pp.perspective(vct(xres/2,yres/2,600))
Circle(pp.x,pp.y),3,6,,,,f
Draw String(pt.x,pt.y),"Point ("&int(pt.x) &"," &int(pt.y) & "," &pt.z &")"
For a As Single=1 To 360 Step 2
Var trace=pt.AxialRotate(centre,a*r,normal_line)
trace=trace.perspective(vct(xres/2,yres/2,600))
Pset(trace.x,trace.y)
Next a
var p=(-250*normal_line+centre).x,q=(-250*normal_line+centre).y
L1.v1=vct(p,q,-250*(normal_line+centre).z)
Line(p,q)-(centre.x,centre.y)
Draw String(p,q),"Normal (" & Int(p) &"," &int(q) &"," &int(-250*(normal_line+centre).z) &")"
p=(250*normal_line+centre).x:q=(250*normal_line+centre).y
L1.v2=vct(p,q,250*(normal_line+centre).z)
Line(p,q)-(centre.x,centre.y)
Draw String(p,q),"Normal (" & Int(p) &"," &int(q) &"," &int(250*(normal_line+centre).z) &")"
''Circle(centre.x,centre.y),5,3,,,,f
Var rad=map(-300,300,ans.z,20,10)
Circle(ans.x,ans.y),rad,,,,,f
Draw String(10,10),"FPS " &fps
Draw String(10,30),"Use Up,Down,Left,Right,q and w to turn the axis"
Draw String(10,50),"Up/down for x axis, Left/right for y axis, q/w for z axis"
Draw String(10,70),"Space resets"
Draw String(10,90),"Use the mouse to move the point"
Draw String(10,110),"Distance of point from normal " &int(dist)
Draw String(10,140),"NOTE - Positive Z is into the screen"
Screenunlock
Sleep snooze,1
#endmacro
#macro mouse()
Dim As Integer x=mx,y=my,dx,dy
While mb = 1
Display()
Getmouse mx,my,,mb
If onscreen Then
If mx<>x Or my<>y Then
dx = mx - x
dy = my - y
x = mx
y = my
pt.x=x+dx
pt.y=y+dy
End If
End If
Wend
#endmacro
Dim As Integer xres,yres
Screen 20
Screeninfo xres,yres
Dim As v3 centre=Type<v3>(xres/2,yres/2,0),pt
Dim As v3 normal_line
Dim As Single angle
Dim As v3 ans,Rangle,rot=vct(0,1,.5)+centre
Dim As v3 firstpoint=centre-vct(200,0,0)
Dim As line3d L1
pt=firstpoint
Dim As long flag,fps
Dim As Integer mx,my,mb,mw
Do
Getmouse mx,my,mw,mb
var snooze=regulate(60,fps)
flag=keys(Rangle)
If flag=1 Then rot=rot.PointRotate(centre,Rangle):Rangle=vct(0,0,0):flag=0
If flag=2 Then rot=vct(0,1,.5)+centre:Rangle=vct(0,0,0):flag=0:pt=firstpoint
normal_line=vct(rot.x,rot.y,rot.z)-centre
display()
If incircle(pt.x,pt.y,10,mx,my) Then
mouse()
End If
Loop Until Multikey(SC_ESCAPE)
Sleep
But Matrix multiplication by looping is far too slow, and the matrix is only 3 X 3 anyway, so it can be written straight down.
Thus the contorted looking return in the rotator.