3d without openGL

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

Re: 3d without openGL

Post by BasicCoder2 »

@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: 1641
Joined: Jun 04, 2005 9:51

Re: 3d without openGL

Post by dafhi »

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

Re: 3d without openGL

Post by dodicat »

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

#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)
           
            
         
BasicCoder2
Posts: 3906
Joined: Jan 01, 2009 7:03
Location: Australia

Re: 3d without openGL

Post by BasicCoder2 »

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

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

Re: 3d without openGL

Post by paul doe »

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
Moderator
Posts: 1730
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: 3d without openGL

Post by paul doe »

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
Moderator
Posts: 1730
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: 3d without openGL

Post by paul doe »

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
Moderator
Posts: 1730
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: 3d without openGL

Post by paul doe »

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

'' 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 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
Moderator
Posts: 1730
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: 3d without openGL

Post by paul doe »

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
Moderator
Posts: 1730
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: 3d without openGL

Post by paul doe »

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
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
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: 3906
Joined: Jan 01, 2009 7:03
Location: Australia

Re: 3d without openGL

Post by BasicCoder2 »

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.
https://www.youtube.com/watch?v=oHsUbGTIXyE
.
dafhi
Posts: 1641
Joined: Jun 04, 2005 9:51

Re: 3d without openGL

Post by dafhi »

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: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: 3d without openGL

Post by dodicat »

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:

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) 
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: 3906
Joined: Jan 01, 2009 7:03
Location: Australia

Re: 3d without openGL

Post by BasicCoder2 »

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: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: 3d without openGL

Post by dodicat »

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

 

  
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.
Post Reply