3d without openGL

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

Re: 3d without openGL

Postby BasicCoder2 » Oct 10, 2017 19:32

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

Re: 3d without openGL

Postby dafhi » Oct 10, 2017 21:02

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

Re: 3d without openGL

Postby dodicat » Oct 10, 2017 21:28

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: 3349
Joined: Jan 01, 2009 7:03

Re: 3d without openGL

Postby BasicCoder2 » Oct 11, 2017 3:40

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

Re: 3d without openGL

Postby paul doe » Oct 11, 2017 4:33

owen wrote:@ paul doe

This problem disappears if you use matrices


I couldn't agree with you more.
And not only that, using matrices avoids gimbal lock

in my example you can see it encounter gimbal lock by first rotating 90 degrees on the x axis, then 90 on y axis and then stop at 89 on the z.
when you get to 89 on the z axis and proceed to rotate to 90 look closely and notice a tiny glitch.

what a wonderful experience it was to try and avoid gimbal lock to no avail. well actually i kinda did it but what length i went to.

im glad that my example can serve the perfect example of what not to do. LOL

Exactly. Thanks for pointing it out. I can imagine that you had a lot of fun trying to avoid the lock with Euler angles X-D. Don't fret. We all used Euler, at one point of another. In physics is the same, you stick with an Euler integrator until the simulation craps out on you ('explodes' to infinity).

Fun fact: the Apollo 11 mission almost failed thanks to the cheap a$$es of the NASA, because they didn't wanted to install another gimbal in the landing module to cut costs. Look at the Wikipedia entry for Gimbal Lock, quite amusing LOL.

@dodicat: pretty nice!. Reminds me of Point Clouds. I had some code on that, will look for it.

@BasicCoder2: aha. I see it now. I think that you are confusing rotation about an arbitrary axis with rotation around the object's own axes. One is extrinsic (outside the object's coordinate system) and the other is intrinsic. I will now post some code that will help you clarify that (hopefully).
paul doe
Posts: 912
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: 3d without openGL

Postby paul doe » Oct 11, 2017 4:35

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

Hi dafhi

Glad to hear that. Which one in particular sparkled your interest?
Last edited by paul doe on Oct 11, 2017 11:33, edited 1 time in total.
paul doe
Posts: 912
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: 3d without openGL

Postby paul doe » Oct 11, 2017 4:45

BasicCoder2 wrote:Demonstration without enough explanation for my programming limitations.
What I tried to do was modify your code to produce the three axes of the object (without the object). That is the axes were the object. I tried to make the z blue the x red and the y green but they appear to be random colors. I wanted to be able to rotate around any of those three axes.

BasicCoder2:

As I promised, here is some code that will help to understand how coordinate systems work. You will need five more files (besides the ones that I already posted). Check the next post.
Last edited by paul doe on Oct 11, 2017 11:35, edited 2 times in total.
paul doe
Posts: 912
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: 3d without openGL

Postby paul doe » Oct 11, 2017 4:49

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

Re: 3d without openGL

Postby paul doe » Oct 11, 2017 4:55

vec4.bi

Code: Select all

#ifndef __vec4__
   #define __vec4__   
   
   #include once "math.bi"
   
   /'
      homogeneous 4 tuple vector type
      
      it is meant to be used as an homogeneous 3D vector (like the ones
      used by OpenGL or Direct3D)
      
      conceptually, they are to be interpreted like this:
      
      | x |
      | y |
      | z |
      | 1 |
      
      which is known as an homogeneous column vector
      
      10/7/2017: fixed a typo in the + operator
         now this was the original code. Take a very close look at it (the comment was removed
         after the fix out of pure anger, as you can imagine):
            
            operator +( byref v as vec4, byref w as vec4 ) as vec4
               '' addition operator (all rocket science stuff...)
               return( vec4( v.x + w.x, v.y + w.y, v.x + w.z ) )
            end operator
         
         goes to show you what happens when you are a #%$@ imbecile and write code in a hurry.
         STAY TUNED FOR THE NExT ExCITING TUTORIAL: how to write unbelievable $%#@ 3D code!!!
   '/   
   type vec4
      public:
         as double x
         as double y
         as double z
         as double w
      
         declare constructor( byval nx as double = 0.0, byval ny as double = 0.0, byval nz as double = 0.0, byval nw as double = 1.0 )
         declare constructor( byref rhs as vec4 )
         declare operator let( byref rhs as vec4 )
         
         declare operator cast() as string
         
         declare function length() as double
         declare function squaredLength() as double
         declare sub normalize()
         declare function normal() as vec4
         declare sub homogeneize()
         declare function homogeneous() as vec4
         declare function cross( byref v as vec4 ) as vec4
         declare function dot( byref v as vec4 ) as double
         declare function distance( byref v as vec4 ) as double
   end type
   
   constructor vec4( byval nx as double = 0.0, byval ny as double = 0.0, byval nz as double = 0.0, byval nw as double = 1.0 )
      '' default constructor creates an homogeneous vector
      x = nx
      y = ny
      z = nz
      w = nw
   end constructor
   
   constructor vec4( byref rhs as vec4 )
      '' copy constructor
      x = rhs.x
      y = rhs.y
      z = rhs.z
      w = rhs.w
   end constructor
   
   operator vec4.let( byref rhs as vec4 )
      '' assignment constructor
      x = rhs.x
      y = rhs.y
      z = rhs.z
      w = rhs.w
   end operator
   
   operator vec4.cast() as string
      '' human readable string representation (useful for debugging)
      return( _
         "| " & trim( str( x ) ) & " |" & chr( 13 ) & chr( 10 ) & _
         "| " & trim( str( y ) ) & " |" & chr( 13 ) & chr( 10 ) & _
         "| " & trim( str( z ) ) & " |" & chr( 13 ) & chr( 10 ) & _
         "| " & trim( str( w ) ) & " |" & chr( 13 ) & chr( 10 ) )
   end operator
   
   function vec4.squaredLength() as double
      /'
         returns the squared length of this vector
         useful when you just want to compare which one is bigger, as
            this avoids having to compute a square root
      '/
      return( x * x + y * y + z * z )
   end function
   
   function vec4.length() as double
      '' returns the length of this vector
      return( sqr( x * x + y * y + z * z ) )
   end function
   
   sub vec4.normalize()
      /'
         normalizes the vector
            note that the homogeneous coordinate (w) is not touched
      '/
      dim as double l = length()
      
      if l > 0 then
         x /= l : y /= l : z /= l
      end if
   end sub
   
   function normalize( byref v as vec4 ) as vec4
      '' for compatibility
      return( v.normal() )
   end function
   
   function vec4.normal() as vec4
      /'
         returns this vector normalized but without altering itself
         again the homogeneous coordinate is left alone
      '/
      dim as double l = length()
      dim as vec4 v = vec4( this )
      
      if l > 0 then
         v.x /= l : v.y /= l : v.z /= l
      end if
      
      return( v )
   end function
   
   sub vec4.homogeneize()
      /'
         homogeneizes the vector
            this is done by dividing the components by the homogeneous coordinate (w)
      '/
      x /= w : y /= w : z /= w : w /= w
   end sub
   
   function vec4.homogeneous() as vec4
      /'
         returns this vector homogeneized but without altering it
      '/
      return( vec4( x / w, y / w, z / w, w / w ) )
   end function
   
   function vec4.cross( byref v as vec4 ) as vec4
      /'
         returns the cross product (aka vectorial product) of this vector and
         another vector v
            note that there's no cross product defined for a 4 tuple, so
            we simply use a standard 3d cross product, and make the w component 1
      '/
      return( vec4( _
         v.y * z - v.z * y, _
         v.z * x - v.x * z, _
         v.x * y - v.y * x, _
         1.0 ) )
   end function
   
   function cross( byref v as vec4, byref w as vec4 ) as vec4
      '' returns the cross product between vectors v and w
      return( v.cross( w ) )
   end function
   
   function vec4.dot( byref v as vec4 ) as double
     '' returns the dot product (aka scalar product) of this vector and vector v
      return( v.x * x + v.y * y + v.z * z )
   end function
   
   function dot( byref v as vec4, byref w as vec4 ) as double
      '' returns the dot product between two vectors
      return( v.dot( w ) )
   end function
   
   function vec4.distance( byref v as vec4 ) as double
     /'
        gets the distance of this vector with vector v
           to calculate the distance, substract them and calculate the length of the resultant vector
      '/
      return( sqr( ( v.x - x ) ^ 2 + ( v.y - y ) ^ 2 + ( v.z - z ) ^ 2 ) )
   end function
   
   function distance( byref v as vec4, byref w as vec4 ) as double
      return( sqr( ( v.x - w.x ) ^ 2 + ( v.y - w.y ) ^ 2 + ( v.z - w.z ) ^ 2 ) )
   end function
   
   operator -( byref v as vec4, byref w as vec4 ) as vec4
     '' substraction operator
     return( vec4( v.x - w.x, v.y - w.y, v.z - w.z ) )
   end operator
   
   operator -( byref v as vec4 ) as vec4
      '' negation operator
      return( vec4( -v.x, -v.y, -v.z ) )
   end operator
   
   operator +( byref v as vec4, byref w as vec4 ) as vec4
      return( vec4( v.x + w.x, v.y + w.y, v.z + w.z ) )
   end operator
   
   operator *( byref v as vec4, byref a as double ) as vec4
      /'
         multiplication with a scalar
            note that this is not a 'scalar product', but a multiplication with a number
            vectors do not define multiplications per se, they define the dot product
            and the cross product. To avoid confusion (and also correctness), the
            multiplication operator is overloaded to a scaling of the vector
      '/
      return( vec4( v.x * a, v.y * a, v.z * a ) )
   end operator
   
   operator *( byref a as double, byref v as vec4 ) as vec4
      '' same as above but with the parameters inversed
      return( vec4( v.x * a, v.y * a, v.z * a ) )
   end operator
   
   operator /( byref v as vec4, byref a as double ) as vec4
      '' division by a scalar. See note above on multiplying a vector
      return( vec4( v.x / a, v.y / a, v.z / a ) )
   end operator
      
   function vectorAngle( byref v as vec4, byref w as vec4 ) as double
      /'
         returns the angle between two vectors using the dot product, in radians
      note that the result of the dot product used here
      should mathematically speaking, always give a result between -1 and 1. Due to imprecisions of
      numerical calculations it might sometimes be a little bit outside this range however (especially
      if you define Scalar to be float instead of double). If that happens, the acos function will give
      an invalid result. So instead a protection was added that sets the value back to 1 or -1
      (because, if the value became 1.0000023 for example, it was probably supposed to be 1 anyway).
       
         dotProduct( v, w ) = length( v ) * length( w ) * cos( angle )
         angle = aCos( dotProduct / ( length( v ) * length( w ) ) )
               
               thus:
               
               angle = aCos( dot( normal( v ) * normal( w ) ) )
      '/
     dim as double angleCos = dot( v.normal(), w.normal() )
    
     '' for acos, the value has to be between -1.0 and 1.0, but due to numerical imprecisions it
     '' sometimes comes outside this range
      
      angleCos = clamp( -1.0, 1.0, angleCos )
   
     return( -acos( angleCos ) )
   end function
   
   function rotateAroundAxis( byref v as vec4, byref axis as vec4, byval angle as double ) as vec4
    /'
       rotate vector v around arbitrary axis for angle radians
       it can only rotate around an axis through our object, to rotate around another axis:
       first translate the object to the axis, then use this function, then translate back
       in the new direction.
      '/
    if( ( v.x = 0 ) and ( v.y = 0 ) and ( v.z = 0 ) ) then
       return vec4( 0.0, 0.0, 0.0 )
    end if

    dim nAxis as vec4 = vec4( axis.x, axis.y, axis.z )
    nAxis.normalize()

    '' calculate parameters of the rotation matrix
    dim as double c = cos( angle )
    dim as double s = sin( angle )
    dim as double t = 1 - c

    '' multiply w with rotation matrix
    dim w as vec4
   
    w.x = ( t * nAxis.x * nAxis.x + c ) * v.x _
        + ( t * nAxis.x * nAxis.y + s * nAxis.z ) * v.y _
        + ( t * nAxis.x * nAxis.z - s * nAxis.y ) * v.z

    w.y = ( t * nAxis.x * nAxis.y - s * nAxis.z ) * v.x _
        + ( t * nAxis.y * nAxis.y + c ) * v.y _
        + ( t * nAxis.y * nAxis.z + s * nAxis.x ) * v.z

    w.z = ( t * nAxis.x * nAxis.z + s * nAxis.y ) * v.x _
        + ( t * nAxis.y * nAxis.z - s * nAxis.x ) * v.y _
        + ( t * nAxis.z * nAxis.z + c ) * v.z
      
      '' the vector has to retain its length, so it's normalized and
      '' multiplied with the original length
    w.normalize()
    w = w * v.length()

    return( w )
   end function
   
#endif


mat4.bi

Code: Select all

#ifndef __mat4__
   #define __mat4__

   #include once "core.bi"
   #include once "platform.bi"
   #include once "vec4.bi"
   /'
                             | a b c d |
      4x4 matrix type      | e f g h |
                            | i j k l |
                            | m n o p |
      
      9/30/2017: improved 4x4 matrix inverse calculation. With -gen gcc -O max
         (the settings I always use) it is more than 60% faster than the previous
         version. A byproduct of correcting the bug mentioned below.
      9/29/2017: fixed determinant calculation (it was erroneously computed)
         The funny thing was that this library is used in various applications,
         including my 3D engine and various tools. When I was implementing a
         feature, it kept doing weird things, and the bug was finally tracked to
         an erroneous calculation of the determinant. The determinant was
         correctly calculated in the 3x3   matrix code, but when I ported the code
         to use OpenGL for rendering, it   didn't worked as intended. Goes to show
         you that one is to be extra   careful with the math code, as it is the
         foundation of the entire game engine.
         
         the calculations were cross checked with the help of this online resource:
            
            https://www.mathsisfun.com/algebra/matrix-calculator.html
         
         which has a very neat matrix calculator for various dimensions.
   '/   
   type mat4
      public:
         as double a, b, c, d
         as double e, f, g, h
         as double i, j, k, l
         as double m, n, o, p
         
         declare constructor( _
            byval sa as double = 1.0, byval sb as double = 0.0, byval sc as double = 0.0, byval sd as double = 0.0, _
            byval se as double = 0.0, byval sf as double = 1.0, byval sg as double = 0.0, byval sh as double = 0.0, _
            byval si as double = 0.0, byval sj as double = 0.0, byval sk as double = 1.0, byval sl as double = 0.0, _
            byval sm as double = 0.0, byval sn as double = 0.0, byval so as double = 0.0, byval sp as double = 1.0 _
            )
         declare constructor( byref NM as mat4 )
         declare operator let( byref RHS as mat4 )
         declare operator cast() as string
         
         declare function determinant() as double
         declare sub transpose()
         declare sub inverse() '' slower version
         'declare sub inverse()
         declare sub identity()
   end type
   
   constructor mat4( _
            byval sa as double = 1.0, byval sb as double = 0.0, byval sc as double = 0.0, byval sd as double = 0.0, _
            byval se as double = 0.0, byval sf as double = 1.0, byval sg as double = 0.0, byval sh as double = 0.0, _
            byval si as double = 0.0, byval sj as double = 0.0, byval sk as double = 1.0, byval sl as double = 0.0, _
            byval sm as double = 0.0, byval sn as double = 0.0, byval so as double = 0.0, byval sp as double = 1.0 _
            )
      /'
         default constructor initializes the matrix to an identity, if no coefficients are specified.
         this is far more useful than initializing it to all zeros
      '/
      a = sa : b = sb : c = sc : d = sd
      e = se : f = sf : g = sg : h = sh
      i = si : j = sj : k = sk : l = sl
      m = sm : n = sn : o = so : p = sp
   end constructor
   
   constructor mat4( byref RHS as mat4 )
      '' copy constructor
      a = RHS.a : b = RHS.b : c = RHS.c : d = RHS.d
      e = RHS.e : f = RHS.f : g = RHS.g : h = RHS.h
      i = RHS.i : j = RHS.j : k = RHS.k : l = RHS.l
      m = RHS.m : n = RHS.n : o = RHS.o : p = RHS.p
   end constructor
   
   operator mat4.let( byref RHS as mat4 )
      '' assignment construction
      a = RHS.a : b = RHS.b : c = RHS.c : d = RHS.d
      e = RHS.e : f = RHS.f : g = RHS.g : h = RHS.h
      i = RHS.i : j = RHS.j : k = RHS.k : l = RHS.l
      m = RHS.m : n = RHS.n : o = RHS.o : p = RHS.p
   end operator
   
   operator mat4.cast() as string
      '' the matrix in a human readable form (very useful for debugging purposes)
      return( _
         "| " & trim( str( a ) ) & " | " & trim( str( b ) ) & " | " & trim( str( c ) ) & " | " & trim( str( d ) ) & " |" & chr( 13 ) & chr( 10 ) & _
         "| " & trim( str( e ) ) & " | " & trim( str( f ) ) & " | " & trim( str( g ) ) & " | " & trim( str( h ) ) & " |" & chr( 13 ) & chr( 10 ) & _
         "| " & trim( str( i ) ) & " | " & trim( str( j ) ) & " | " & trim( str( k ) ) & " | " & trim( str( l ) ) & " |" & chr( 13 ) & chr( 10 ) & _
         "| " & trim( str( m ) ) & " | " & trim( str( n ) ) & " | " & trim( str( o ) ) & " | " & trim( str( p ) ) & " |" & chr( 13 ) & chr( 10 ) _
      )
   end operator
   
   function mat4.determinant() as double
      /'
         computes the determinant of the matrix using Laplace cofactor expansion
         
         the determinant of a 3x3 matrix is:
            a * ( e * i - f * h ) - b * ( d * i - f * g ) + c * ( d * h - e * g )
         
         and a 4x4 matrix determinant is given by:
            a *         b *         c *         d *
            | f g h |   | e g h |   | e f h |   | e f g |
            | j k l | - | i k l | + | i j l | - | i j k |
            | n o p |   | m o p |   | m n p |   | m n o |
            
            where the '|' means the determinant of the inner 3x3 matrices. Note that the
            cofactors are already factored in the calculation.
            
         the determinant is thus:
            + ( a * (   f * ( k * p - l * o ) - g * ( j * p - l * n ) + h * ( j * o - k * n ) ) )
            - ( b * ( e * ( k * p - l * o ) - g * ( i * p - l * m ) + h * ( i * o - k * m ) ) )
            +   ( c * ( e * ( j * p - l * n ) - f * ( i * p - l * m ) + h * ( i * n - j * m ) ) )
            - ( d * ( e * ( j * o - k * n ) - f * ( i * o - k * m ) + g * ( i * n - j * m ) ) )
         
      '/
      dim as double det = _
            ( a * (   f * ( k * p - l * o ) - g * ( j * p - l * n ) + h * ( j * o - k * n ) ) ) _
         - ( b * ( e * ( k * p - l * o ) - g * ( i * p - l * m ) + h * ( i * o - k * m ) ) ) _
         +   ( c * ( e * ( j * p - l * n ) - f * ( i * p - l * m ) + h * ( i * n - j * m ) ) ) _
         - ( d * ( e * ( j * o - k * n ) - f * ( i * o - k * m ) + g * ( i * n - j * m ) ) )
      /'
         this isn't matematically correct, just a programmer's dirty hack.
         if the determinant of a matrix is 0, it means it has no inverse. In the code for
         calculating the inverse, a division by the determinant is performed; and if it is
         zero, a division by zero is performed on *every* element of the matrix, filling it
         with positive or negative infinity values and rendering it useless. A matrix
         without inverse is the matrix   itself, so setting the determinant value to 1
         does the trick.
      '/
      if det = 0 then det = 1.0
      
      return( det )
   end function
   
   sub mat4.transpose()
      /'
        transposes the matrix
        
           [ a b c d ]T    [ a e i m ]
           [ e f g h ]  =  [ b f j n ]
           [ i j k l ]     [ c g k o ]
           [ m n o p ]     [ d h l p ]
       
        why have it, if it is not used by the matrix code itself? Well, there is a
        nice property of matrices, which has to do with rotations. If you can
        be sure that the matrix contains only rotations, transposing it is the
        same as taking its inverse, thus saving you *a lot* of computation
      '/
      this = mat4( a, e, i, m, b, f, j, n, c, g, k, o, d, h, l, p )
   end sub
   
   /'
   sub mat4.inverse()
      /'
         computes the inverse of a 4x4 matrix
         
         this version is 60%+ faster and 400%+ uglier than the previous version.
         it was made so by computing the determinant inside the method and
         recycling as much   calculation as possible
      '/
      '' list of 2x2 determinants
      dim as double kplo = k * p - l * o
      dim as double jpln = j * p - l * n
      dim as double jokn = j * o - k * n
      dim as double iplm = i * p - l * m
      dim as double iokm = i * o - k * m
      dim as double injm = i * n - j * m
      dim as double gpho = g * p - h * o
      dim as double fphn = f * p - h * n
      dim as double fogn = f * o - g * n
      dim as double ephm = e * p - h * m
      dim as double eogm = e * o - g * m
      dim as double enfm = e * n - f * m
      dim as double glhk = g * l - h * k
      dim as double flhj = f * l - h * j
      dim as double fkgj = f * k - g * j
      dim as double elhi = e * l - h * i
      dim as double ekgi = e * k - g * i
      dim as double ejfi = e * j - f * i
   
      '' list of 3x3 determinants
      dim as double d1kplo = f * kplo
      dim as double d1jpln = g * jpln
      dim as double d1jokn = h * jokn
      dim as double d2kplo = e * kplo
      dim as double d2iplm = g * iplm
      dim as double d2iokm = h * iokm
      dim as double d3jpln = e * jpln
      dim as double d3iplm = f * iplm
      dim as double d3injm = h * injm
      dim as double d4jokn = e * jokn
      dim as double d4iokm = f * iokm
      dim as double d4injm = g * injm
      dim as double d5kplo = b * kplo
      dim as double d5jpln = c * jpln
      dim as double d5jokn = d * jokn
      dim as double d6kplo = a * kplo
      dim as double d6iplm = c * iplm
      dim as double d6iokm = d * iokm
      dim as double d7jpln = a * jpln
      dim as double d7iplm = b * iplm
      dim as double d7injm = d * injm
      dim as double d8jokn = a * jokn
      dim as double d8iokm = b * iokm
      dim as double d8injm = c * injm
      dim as double d9gpho = b * gpho
      dim as double d9fphn = c * fphn
      dim as double d9fogn = d * fogn
      dim as double d10gpho = a * gpho
      dim as double d10ephm = c * ephm
      dim as double d10eogm = d * eogm
      dim as double d11fphn = a * fphn
      dim as double d11ephm = b * ephm
      dim as double d11enfm = d * enfm
      dim as double d12fogn = a * fogn
      dim as double d12eogm = b * eogm
      dim as double d12enfm = c * enfm
      dim as double d13glhk = b * glhk
      dim as double d13flhj = c * flhj
      dim as double d13fkgj = d * fkgj
      dim as double d14glhk = a * glhk
      dim as double d14elhi = c * elhi
      dim as double d14ekgi = d * ekgi
      dim as double d15flhj = a * flhj
      dim as double d15elhi = b * elhi
      dim as double d15ejfi = d * ejfi
      dim as double d16fkgj = a * fkgj
      dim as double d16ekgi = b * ekgi
      dim as double d16ejfi = c * ejfi
      
      '' 4x4 determinant (inversed)
      dim as double det = _
            ( a * ( d1kplo - d1jpln + d1jokn ) _
         - ( b * ( d2kplo - d2iplm + d2iokm ) ) _
         + ( c * ( d3jpln - d3iplm + d3injm ) ) _
         - ( d * ( d4jokn - d4iokm + d4injm ) ) )
      
      '' if the determinant is 0, the matrix has no inverse
      if det = 0 then exit sub
      
      '' Multiplying with the reciprocal is slightly faster than dividing
      dim as double invDet = 1.0 / det
      
      '' minors
      dim as double Ma = d1kplo - d1jpln + d1jokn
      dim as double Mb = d2kplo - d2iplm + d2iokm
      dim as double Mc = d3jpln - d3iplm + d3injm
      dim as double Md = d4jokn - d4iokm + d4injm
      dim as double Me = d5kplo - d5jpln + d5jokn
      dim as double Mf = d6kplo - d6iplm + d6iokm
      dim as double Mg = d7jpln - d7iplm + d7injm
      dim as double Mh = d8jokn - d8iokm + d8injm
      dim as double Mi = d9gpho - d9fphn + d9fogn
      dim as double Mj = d10gpho - d10ephm + d10eogm
      dim as double Mk = d11fphn - d11ephm + d11enfm
      dim as double Ml = d12fogn - d12eogm + d12enfm
      dim as double Mm = d13glhk - d13flhj + d13fkgj
      dim as double Mn = d14glhk - d14elhi + d14ekgi
      dim as double Mo = d15flhj - d15elhi + d15ejfi
      dim as double Mp = d16fkgj - d16ekgi + d16ejfi
      
      /'
         adjugate (the adjugate is the transpose of the cofactored matrix of minors)         
         
          Ma   -Me    Mi   -Mm
         -Mb    Mf   -Mj    Mn
          Mc   -Mg    Mk   -Mo
         -Md    Mh   -Ml    Mp
      '/
      this = mat4( _
          Ma * invDet, -Me * invDet,  Mi * invDet, -Mm * invDet, _
         -Mb * invDet,  Mf * invDet, -Mj * invDet,  Mn * invDet, _
          Mc * invDet, -Mg * invDet,  Mk * invDet, -Mo * invDet, _
         -Md * invDet,  Mh * invDet, -Ml * invDet,  Mp * invDet )
   end sub
   '/
   sub mat4.inverse()
     /'
        calculates the inverse of a 4x4 matrix
        
        the minor, cofactor and adjugate matrices are all calculated and factored into the
        code to make it shorter and more efficient
     '/
    
     '' multiplying by the reciprocal of the determinant is faster than dividing by it 
     dim as double invDet = 1 / determinant
    
     this = mat4( _
       ( f * k * p + g * l * n + h * j * o - f * l * o - g * j * p - h * k * n ) * invDet, _
       ( b * l * o + c * j * p + d * k * n - b * k * p - c * l * n - d * j * o ) * invDet, _
       ( b * g * p + c * h * n + d * f * o - b * h * o - c * f * p - d * g * n ) * invDet, _
       ( b * h * k + c * f * l + d * g * j - b * g * l - c * h * j - d * f * k ) * invDet, _
       ( e * l * o + g * i * p + h * k * m - e * k * p - g * l * m - h * i * o ) * invDet, _
       ( a * k * p + c * l * m + d * i * o - a * l * o - c * i * p - d * k * m ) * invDet, _
       ( a * h * o + c * e * p + d * g * m - a * g * p - c * h * m - d * e * o ) * invDet, _
       ( a * g * l + c * h * i + d * e * k - a * h * k - c * e * l - d * g * i ) * invDet, _
       ( e * j * p + f * l * m + h * i * n - e * l * n - f * i * p - h * j * m ) * invDet, _
       ( a * l * n + b * i * p + d * j * m - a * j * p - b * l * m - d * i * n ) * invDet, _
       ( a * f * p + b * h * m + d * e * n - a * h * n - b * e * p - d * f * m ) * invDet, _
       ( a * h * j + b * e * l + d * f * i - a * f * l - b * h * i - d * e * j ) * invDet, _
       ( e * k * n + f * i * o + g * j * m - e * j * o - f * k * m - g * i * n ) * invDet, _
       ( a * j * o + b * k * m + c * i * n - a * k * n - b * i * o - c * j * m ) * invDet, _
       ( a * g * n + b * e * o + c * f * m - a * f * o - b * g * m - c * e * n ) * invDet, _
       ( a * f * k + b * g * i + c * e * j - a * g * j - b * e * k - c * f * i ) * invDet _
      )
      '' lovely, isn't it?
   end sub
   
   sub mat4.identity()
      '' makes the matrix an identity matrix
      a = 1.0 : b = 0.0 : c = 0.0 : d = 0.0
      e = 0.0 : f = 1.0 : g = 0.0 : h = 0.0
      i = 0.0 : j = 0.0 : k = 1.0 : l = 0.0
      m = 0.0 : n = 0.0 : o = 0.0 : p = 1.0
   end sub
   
   operator *( byref A as mat4, byref B as mat4 ) as mat4
      /'
         multiply two 4x4 matrices
         
         remember that matrix multiplication is not commutative!
            A * B != B * A
      '/
      return( mat4( _
      A.a * B.a + A.b * B.e + A.c * B.i + A.d * B.m, _
      A.a * B.b + A.b * B.f + A.c * B.j + A.d * B.n, _
      A.a * B.c + A.b * B.g + A.c * B.k + A.d * B.o, _
      A.a * B.d + A.b * B.h + A.c * B.l + A.d * B.p, _
      A.e * B.a + A.f * B.e + A.g * B.i + A.h * B.m, _
      A.e * B.b + A.f * B.f + A.g * B.j + A.h * B.n, _
      A.e * B.c + A.f * B.g + A.g * B.k + A.h * B.o, _
      A.e * B.d + A.f * B.h + A.g * B.l + A.h * B.p, _
      A.i * B.a + A.j * B.e + A.k * B.i + A.l * B.m, _
      A.i * B.b + A.j * B.f + A.k * B.j + A.l * B.n, _
      A.i * B.c + A.j * B.g + A.k * B.k + A.l * B.o, _
      A.i * B.d + A.j * B.h + A.k * B.l + A.l * B.p, _
      A.m * B.a + A.n * B.e + A.o * B.i + A.p * B.m, _
      A.m * B.b + A.n * B.f + A.o * B.j + A.p * B.n, _
      A.m * B.c + A.n * B.g + A.o * B.k + A.p * B.o, _
      A.m * B.d + A.n * B.h + A.o * B.l + A.p * B.p _
      ) )

      'return( mat4( _
      '   A.a * B.a + A.b * B.e + A.c * B.i + A.d * B.m, _
      '   A.a * B.b + A.b * B.f + A.c * B.j + A.d * B.n, _
      '   A.a * B.c + A.b * B.g + A.c * B.k + A.d * B.o, _
      '   A.a * B.d + A.b * B.h + A.c * B.l + A.d * B.p, _
      '   A.e * B.a + A.f * B.e + A.g * B.i + A.h * B.m, _
      '   A.e * B.b + A.f * B.f + A.g * B.j + A.h * B.n, _
      '   A.e * B.c + A.f * B.g + A.g * B.k + A.h * B.o, _
      '   A.e * B.d + A.f * B.h + A.g * B.l + A.h * B.p, _
      '   A.i * B.a + A.j * B.e + A.k * B.i + A.l * B.m, _
      '   A.i * B.b + A.j * B.f + A.k * B.j + A.l * B.n, _
      '   A.i * B.c + A.j * B.g + A.k * B.k + A.l * B.o, _
      '   A.i * B.d + A.j * B.h + A.k * B.l + A.l * B.p, _
      '   A.m * B.a + A.n * B.e + A.o * B.i + A.p * B.m, _
      '   A.m * B.b + A.n * B.f + A.o * B.j + A.p * B.n, _
      '   A.m * B.c + A.n * B.g + A.o * B.k + A.p * B.o, _
      '   A.m * B.d + A.n * B.h + A.o * B.l + A.p * B.p _
      ') )
   end operator
   
   operator +( byref A as mat4, byref B as mat4 ) as mat4
      '' adds two 4x4 matrices
      return( mat4( _
         A.a + B.a, A.b + B.b, A.c + B.c, A.d + B.d, _
         A.e + B.e, A.f + B.f, A.g + B.g, A.h + B.h, _
         A.i + B.i, A.j + B.j, A.k + B.k, A.l + B.l, _
         A.m + B.m, A.n + B.n, A.o + B.o, A.p + B.p _
      ) )
   end operator
   
   operator -( byref A as mat4, byref B as mat4 ) as mat4
      '' substracts two 4x4 matrices
      return( mat4( _
         A.a - B.a, A.b - B.b, A.c - B.c, A.d - B.d, _
         A.e - B.e, A.f - B.f, A.g - B.g, A.h - B.h, _
         A.i - B.i, A.j - B.j, A.k - B.k, A.l - B.l, _
         A.m - B.m, A.n - B.n, A.o - B.o, A.p - B.p _
      ) )
   end operator
   
   operator -( byref A as mat4 ) as mat4
      '' negates the matrix
      return( mat4( _
         -A.a, -A.b, -A.c, -A.d, _
         -A.e, -A.f, -A.g, -A.h, _
         -A.i, -A.j, -A.k, -A.l, _
         -A.m, -A.n, -A.o, -A.p _
      ) )
   end operator
   
   operator *( byref A as mat4, byref s as double ) as mat4
      '' scalar multiplication
      return( mat4( _
         A.a * s, A.b * s, A.c * s, A.d * s, _
         A.e * s, A.f * s, A.g * s, A.h * s, _
         A.i * s, A.j * s, A.k * s, A.l * s, _
         A.m * s, A.n * s, A.o * s, A.p * s _
      ) )
   end operator
   
   operator *( byref s as double, byref A as mat4 ) as mat4
      '' scalar multiplication
      return( mat4( _
         A.a * s, A.b * s, A.c * s, A.d * s, _
         A.e * s, A.f * s, A.g * s, A.h * s, _
         A.i * s, A.j * s, A.k * s, A.l * s, _
         A.m * s, A.n * s, A.o * s, A.p * s _
      ) )
   end operator
   
   operator /( byref s as double, byref A as mat4 ) as mat4
      /'
         scalar divide
            the 'division' of a matrix by another matrix is not defined, the
            equivalent operation is multiplying one matrix by the inverse of the other
            on scalars though, it can be defined, mostly for convenience purposes
      '/
      return( mat4( _
         A.a / s, A.b / s, A.c / s, A.d / s, _
         A.e / s, A.f / s, A.g / s, A.h / s, _
         A.i / s, A.j / s, A.k / s, A.l / s, _
         A.m / s, A.n / s, A.o / s, A.p / s _
      ) )
   end operator
   
   operator /( byref A as mat4, byref s as double ) as mat4
      return( mat4( _
         A.a / s, A.b / s, A.c / s, A.d / s, _
         A.e / s, A.f / s, A.g / s, A.h / s, _
         A.i / s, A.j / s, A.k / s, A.l / s, _
         A.m / s, A.n / s, A.o / s, A.p / s _
      ) )
   end operator
   
   operator *( byref v as vec4, byref A as mat4 ) as vec4
      /'
         multiply a vector with a row matrix, resulting in a row vector (like Direct3D)
            a row vector looks like this:
            
            | x y z w |
            
            and is the format that Direct3D uses. What this means, code-wise, is that you
            have to pre-multiply the vectors with the matrices, and some other stuff, like
            transposing the matrices if you are using column vectors (as this library does)
      '/
      return( vec4( _
         A.a * v.x + A.e * v.y + A.i * v.z + A.m * v.w, _
         A.b * v.x + A.f * v.y + A.j * v.z + A.n * v.w, _
         A.c * v.x + A.g * v.y + A.k * v.z + A.o * v.w, _
         A.d * v.x + A.h * v.y + A.l * v.z + A.p * v.w _
        ) )
   end operator
   
   operator *( byref A as mat4, byref v as vec4 ) as vec4
      /'
         multiply a vector with a column matrix, resulting in a column vector (like OpenGL)
            a column vector looks like this
            
            | x |
            | y |
            | z |
            | w |
            
            and is the format favored by OpenGL. In this library, column vectors are used, for
            compatibility.
      '/
    'w.X = A.a * v.X + A.b * v.Y + A.c * v.Z + A.d * v.W
    'w.Y = A.e * v.X + A.f * v.Y + A.g * v.Z + A.h * v.W
    'w.Z = A.i * v.X + A.j * v.Y + A.k * v.Z + A.l * v.W
    'w.W = A.m * v.X + A.n * v.Y + A.o * v.Z + A.p * v.W
      return( vec4( _
        A.a * v.X + A.b * v.Y + A.c * v.Z + A.d * v.W, _
        A.e * v.X + A.f * v.Y + A.g * v.Z + A.h * v.W, _
        A.i * v.X + A.j * v.Y + A.k * v.Z + A.l * v.W, _
        A.m * v.X + A.n * v.Y + A.o * v.Z + A.p * v.W _   
      ) )
   end operator
   
   '' utility functions
   function inverse( byref M as mat4 ) as mat4
      '' returns the inverse of the provided matrix
      dim as mat4 I = mat4( M )
      
      I.inverse()
      
      return( I )
   end function
      
   namespace matrices
      function translation( byval t as vec4 ) as mat4
         return mat4( _
            1.0, 0.0, 0.0, t.x, _
            0.0, 1.0, 0.0, t.y, _
            0.0, 0.0, 1.0, t.z, _
            0.0, 0.0, 0.0, 1.0 )
      end function
      
      function identity() as mat4
         '' returns an identity matrix
         return( mat4( _
            1.0, 0.0, 0.0, 0.0, _
            0.0, 1.0, 0.0, 0.0, _
            0.0, 0.0, 1.0, 0.0, _
            0.0, 0.0, 0.0, 1.0 ) )
      end function

   end namespace
   
#endif


utility.bi

Code: Select all

/'
   auxiliary (read: rushed) functions for performing some drawings on the
   3D plane
'/
function clipSegment( byref cam as camera, byref s1 as vec4, byref s2 as vec4 ) as integer
   /'
      clip both vertices of a line if necessary
      
      note that this 'clipping' is only done on the near plane, not on the scene
      
      the sensible (and correct) way to do this would be to clip them against
      the view frustum; but alas, in order to keep the code short and clear
      the calculation of the frustum was removed from the camera class. Hence
      this crappy function.
      
      Clipping (especially against the near plane) is important, because if you
      don't clip, when one of the line endpoints goes behind the near clipping
      plane, it gets a negative coordinate (due to the projection), and the line
      segment is no longer valid
   '/
   if( s1.z >= cam.nearClip() and s2.z >= cam.nearClip() ) then
       ' Neither one is behind the camera, draw them both
       return( -1 )
   elseIf( s1.z < cam.nearClip() and s2.z >= cam.nearClip() ) then
       ' First coordinate behind the camera, clip it
       s1.x = s1.x + ((cam.nearClip() - s1.z) / (s2.z - s1.z)) * (s2.x - s1.x)
       s1.y = s1.y + ((cam.nearClip() - s1.z) / (s2.z - s1.z)) * (s2.y - s1.y)
       s1.z = cam.nearClip()
   
       return( -1 )
   elseIf( s1.z >= cam.nearClip() and s2.z < cam.nearClip() ) then
       ' Second coordinate behind the camera, clip it
       s2.x = s1.x + ((cam.nearClip() - s1.z) / (s2.z - s1.z)) * (s2.x - s1.x)
       s2.y = s1.y + ((cam.nearClip() - s1.z) / (s2.z - s1.z)) * (s2.y - s1.y)
       s2.z = cam.nearClip()
   
       return( -1 )
   else
       ' Both coordinates behind the camera, don't draw
       return( 0 )
   end if
end function

sub drawLine3d( byref cam as camera, byref p1 As vec4, byref p2 As vec4, byval c as uint32 )
    dim as vec4 pos1, pos2

    pos1 = cam.transform( p1 - cam.getPos() )
    pos2 = cam.transform( p2 - cam.getPos() )

      dim as double x1, y1, x2, y2
      
    dim as integer visible = clipSegment( cam, pos1, pos2 )

    if( visible ) then
      dim as double z1, z2

      cam.perspective( pos1, cam.projectionPlane, x1, y1, z1 )
      cam.perspective( pos2, cam.projectionPlane, x2, y2, z2 )

      line( x1, y1 ) - ( x2, y2 ), c

    end if
end sub
paul doe
Posts: 912
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: 3d without openGL

Postby paul doe » Oct 11, 2017 5:02

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: 3349
Joined: Jan 01, 2009 7:03

Re: 3d without openGL

Postby BasicCoder2 » Oct 11, 2017 5:36

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

Re: 3d without openGL

Postby dafhi » Oct 11, 2017 8:01

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

Re: 3d without openGL

Postby dodicat » Oct 11, 2017 9:46

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: 3349
Joined: Jan 01, 2009 7:03

Re: 3d without openGL

Postby BasicCoder2 » Oct 11, 2017 11:40

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

Re: 3d without openGL

Postby dodicat » Oct 11, 2017 11:55

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.

Return to “Projects”

Who is online

Users browsing this forum: Haubitze and 2 guests