3d without openGL

User projects written in or related to FreeBASIC.
BasicCoder2
Posts: 3400
Joined: Jan 01, 2009 7:03

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.
.
dafhi
Posts: 1251
Joined: Jun 04, 2005 9:51

Re: 3d without openGL

@paul doe - Very intrigued and inspired by your offerings
dodicat
Posts: 5951
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: 3d without openGL

OK basiccoder2.
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

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)

BasicCoder2
Posts: 3400
Joined: Jan 01, 2009 7:03

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.

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

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)
paul doe
Posts: 922
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: 3d without openGL

owen wrote:@ paul doe

This problem disappears if you use matrices

I couldn't agree with you more.
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

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

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).
paul doe
Posts: 922
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: 3d without openGL

dafhi wrote:@paul doe - Very intrigued and inspired by your offerings

Hi dafhi

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.
paul doe
Posts: 922
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: 3d without openGL

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.

BasicCoder2:

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.
paul doe
Posts: 922
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: 3d without openGL

camera.bas

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

object.bas

Code: Select all

#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 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
Last edited by paul doe on Oct 13, 2017 10:27, edited 2 times in total.
paul doe
Posts: 922
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: 3d without openGL

vec4.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

mat4.bi

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

utility.bi

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
paul doe
Posts: 922
Joined: Jul 25, 2017 17:22
Location: Argentina

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

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

'' 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

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
BasicCoder2
Posts: 3400
Joined: Jan 01, 2009 7:03

Re: 3d without openGL

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

You get an A+ on the code. It does exactly what I would have liked to have done.
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.
.
dafhi
Posts: 1251
Joined: Jun 04, 2005 9:51

Re: 3d without openGL

paul doe wrote:
dafhi wrote:@paul doe - Very intrigued and inspired by your offerings

Hi dafhi

Glad to hear that. Which one in particular sparkled your interest?

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.

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.
dodicat
Posts: 5951
Joined: Jan 10, 2006 20:30
Location: Scotland

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.

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

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)

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.
BasicCoder2
Posts: 3400
Joined: Jan 01, 2009 7:03

Re: 3d without openGL

dodicat wrote:For axial rotation Rodrigues method can also be used (I have posted many examples).

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.

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.
dodicat
Posts: 5951
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: 3d without openGL

Hi basiccoder2.
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 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 .

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

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.