## raytrace try

General FreeBASIC programming questions.
dodicat
Posts: 6645
Joined: Jan 10, 2006 20:30
Location: Scotland

### Re: raytrace try

The main thing is udt's are passed byref by default.
So you should stipulate byval if you are messing with them inside a function.
You should check the cone hit maths.
Why create constructors when they are not needed, they add a function call each time.
I think operators on vectors slow things, although they are handy.
I am not keen on using not.

Code: Select all

'' bluatigro 3 okt 2018
'' slow ray 1
'' this can be speeded up if i undestand :
''https://www.cl.cam.ac.uk/teaching/1999/AGraphHCI/SMAG/node2.html#SECTION00023100000000000000
type d3d
dim as double x , y , z
'declare constructor ()
'declare constructor ( a as double , b as double , c as double )
declare function tocolor() as ulong
end type
'constructor d3d()
' x = 0
' y = 0
' z = 0
'end constructor
'constructor d3d( a as double , b as double , c as double )
'x = a
'y = b
' z = c
'end constructor
operator + ( a as d3d , b as d3d ) as d3d
return type<d3d>( a.x + b.x , a.y + b.y , a.z + b.z )
end operator
operator - ( a as d3d , b as d3d ) as d3d
return type<d3d>( a.x - b.x , a.y - b.y , a.z - b.z )
end operator
function d3d.tocolor() as ulong
return rgb( x * 255 , y * 255 , z * 255 )
end function

type tsphere
dim as d3d c , kl
dim as double d
declare function hit(byval nu as d3d ) as integer
end type
function tsphere.hit(byval nu as d3d ) as integer
nu = nu - c
'nu.x-=c.x
'nu.y-=c.y
'nu.z-=c.z
return ( nu.x * nu.x + nu.y * nu.y + nu.z * nu.z ) < d * d
end function
type tcylinder
dim as d3d c , kl
dim as double d , h
declare function hit(byval nu as d3d ) as integer
end type
function tcylinder.hit(byval nu as d3d ) as integer
nu = nu - c
'nu.x-=c.x
' nu.y-=c.y
'nu.z-=c.z
return ( nu.x * nu.x + nu.z  * nu.z ) < d * d and abs( nu.y ) < h
end function
type tcone
dim as d3d c , kl
dim as double d , h
declare function hit( byval nu as d3d ) as integer
end type
function tcone.hit(byval nu as d3d ) as integer
nu = nu - c
'nu.x-=c.x
' nu.y-=c.y
'nu.z-=c.z
return ( nu.x * nu.x + nu.z * nu.z ) < ( nu.z * nu.z ) / (( h * h ) * ( d * d )) and abs( nu.y ) < h
end function

dim as tsphere sphere
sphere.c = type<d3d>( -150 , 0 , 1000 )
sphere.d = 50
sphere.kl= type<d3d>( 1 , 0 , 0 ) ''red
dim as tcylinder cylinder
cylinder.c = type<d3d>( 0 , 0 , 1000 )
cylinder.d = 30
cylinder.h = 70
cylinder.kl = type<d3d>( 0 , 1 , 0 ) ''green
dim as tcone cone
cone.c = type<d3d>( 150 , 0 , 1000 )
cone.d = 50
cone.h = 100
cone.kl = type<d3d>( 0 , 0 , 1 ) ''blue

screen 20 , 32
dim as integer winx , winy
screeninfo winx , winy
dim as d3d nu , dnu
dim as double x , y
dim as ulong kl
paint ( 1 , 1 ) , &h7f7f7f
for x = -250 to 250
for y = -120 to 120
nu = type<d3d>( x , y , 930 )
dnu = type<d3d>( 0 , 0 , 2 )

while sphere.hit( nu )=0 _
andalso cylinder.hit( nu ) =0 _
andalso nu.z < 1070
nu += dnu
dnu += type<d3d>( 0 , 0 , 1 )
wend

'while not sphere.hit( nu ) _
'and not cylinder.hit( nu ) _
'and not cone.hit( nu ) _
' and nu.z < 1070
' nu += dnu
'dnu += type<d3d>( 0 , 0 , 1 )
'wend
kl = rgb(0,0,0)
if sphere.hit( nu ) then kl = sphere.kl.tocolor()
if cylinder.hit( nu ) then kl = cylinder.kl.tocolor()
if cone.hit( nu ) then kl = cone.kl.tocolor()
pset ( x + winx / 2 , y + winy / 2 ) , kl
next y
next x
sleep
bluatigro
Posts: 651
Joined: Apr 25, 2012 10:35
Location: netherlands

### Re: raytrace try

i have to use operator's or function's
i m a litle layzy : i do not like to type
i m not used to macro's
dafhi
Posts: 1357
Joined: Jun 04, 2005 9:51

### Re: raytrace try

dodicat that's very nice. i'll have a look

ray tracing is a lot different than ray marching
I only know the sphere

Code: Select all

'' bluatigro 3 okt 2018
'' slow ray 1
'' this can be speeded up if i undestand :
''https://www.cl.cam.ac.uk/teaching/1999/AGraphHCI/SMAG/node2.html#SECTION00023100000000000000

type d3d
dim as double x , y , z
declare constructor ()
declare constructor ( a as double , b as double , c as double )
declare function tocolor() as ulong
declare property lenSq as single
declare sub normalize
end type
constructor d3d()
x = 0
y = 0
z = 0
end constructor
constructor d3d( a as double , b as double , c as double )
x = a
y = b
z = c
end constructor
operator + ( a as d3d , b as d3d ) as d3d
return d3d( a.x + b.x , a.y + b.y , a.z + b.z )
end operator
operator - ( a as d3d , b as d3d ) as d3d
return d3d( a.x - b.x , a.y - b.y , a.z - b.z )
end operator
operator * ( a as d3d , b as d3d ) as single
return a.x * b.x + a.y * b.y + a.z * b.z
end operator
function d3d.tocolor() as ulong
return rgb( x * 255 , y * 255 , z * 255 )
end function

property d3d.lenSq as Single
return x*x+y*y+z*z
end property
sub d3d.normalize
dim as single s=1/sqr(lenSq): x*=s:y*=s:z*=s
end sub

type tray
as d3d o, d
end type

type tsphere
dim as d3d c , kl
dim as double d
'declare function hit( nu as d3d ) as integer
declare function intersect(byref as tray ptr) as single
end type
'function tsphere.hit( nu as d3d ) as integer
'  nu = nu - c
'  return ( nu.x * nu.x + nu.y * nu.y + nu.z * nu.z ) < d * d
'end function
function tsphere.intersect(byref ray as tray ptr) as single
/' sphere: (p-c)^2 = r^2

(o + t*d - c)^2 = r^2:  find t

o(o+td-c) + td(o+td-c) - c(o+td-c) = rr
oo+tdo-co + tdo+ttdd-tdc - tdc-co+cc
oo+2tdo-2co +ttdd -2tdc + cc

ttdd + 2t(do-dc) - 2co + oo + cc - rr = 0:  quadratic
tt*A + 2t*B - C = 0

A = dd
B = d(o-c)
C = oo + cc - rr

'/

var rayocen = ray->o - c
dim as single b = ray->d*rayocen
var _c = rayocen*rayocen - d*d
if _c > b*b then return -1
dim as single sqrdisc = sqr((b*b - _c)*4)/2
if sqrdisc-b < 0 then return -1
return -b-sqrdisc

end function

type tcylinder
dim as d3d c , kl
dim as double d , h
declare function hit( nu as d3d ) as integer
end type
function tcylinder.hit( nu as d3d ) as integer
nu = nu - c
return ( nu.x * nu.x + nu.z  * nu.z ) < d * d and abs( nu.y ) < h
end function

type tcone
dim as d3d c , kl
dim as double d , h
declare function hit( nu as d3d ) as integer
end type
function tcone.hit( nu as d3d ) as integer
nu = nu - c
return ( nu.x * nu.x + nu.z * nu.z ) < ( nu.z * nu.z ) / ( h * h ) * ( d * d ) and abs( nu.y ) < h
end function

dim as tsphere sphere
sphere.c = d3d( 0 , 0 , 1000 )
sphere.d = 2
sphere.kl= d3d( 1 , 0 , 0 ) ''red
dim as tcylinder cylinder
cylinder.c = d3d( 0 , 0 , 1000 )
cylinder.d = 30
cylinder.h = 70
cylinder.kl = d3d( 0 , 1 , 0 ) ''green
dim as tcone cone
cone.c = d3d( 150 , 0 , 1000 )
cone.d = 50
cone.h = 100
cone.kl = d3d( 0 , 0 , 1 ) ''blue

screen 20 , 32
dim as integer winx , winy
screeninfo winx , winy
dim as d3d nu , dnu
dim as double x , y
dim as ulong kl
paint ( 1 , 1 ) , &h7f7f7f

dim as tray ray
ray.o = d3d( 0, 0, -50 )
sphere.d = 40

for z as single = -10 to 40 step .5
sphere.c = d3d(0,0,z)
'for x = -0 to 0
for x = -250 to 250
'for y = -0 to 0
for y = -120 to 120
ray.d = d3d( x, y, 0 ) - ray.o
ray.d.normalize
'nu = d3d( x , y , 930 )
'dnu = d3d( 0 , 0 , 2 )

/'
while not sphere.hit( nu ) _
and not cylinder.hit( nu ) _
and not cone.hit( nu ) _
and nu.z < 1070
nu += dnu
dnu += d3d( 0 , 0 , 1 )
wend
'/
var t = sphere.intersect( @ray )
dim as ubyte b = t
'if t>0 andalso t < INF then kl = rgb(255,0,0)
kl = rgb(b,b,b)
'if sphere.hit( nu ) then kl = sphere.kl.tocolor()
'if cylinder.hit( nu ) then kl = cylinder.kl.tocolor()
'if cone.hit( nu ) then kl = cone.kl.tocolor()
pset ( x + winx / 2 , y + winy / 2 ) , kl
next y
next x
sleep 1
if inkey<>"" then end
next
sleep
bluatigro
Posts: 651
Joined: Apr 25, 2012 10:35
Location: netherlands

### Re: raytrace try

update :
sphere's whit texture's

you shoot see a 'earth' and a 'moon'

cone's , cylinder's , tori

i know triangles but i did not het good results whit texture's on that

Code: Select all

'' bluatigro 14 nov 2018
'' bluaray 1
#include "noise.bas"

dim shared as integer winx , winy
screen 18 , 32
screeninfo winx , winy
type d3d
dim as double x , y , z
declare constructor ()
declare constructor ( _
a as double , b as double , c as double )
declare function tocolor() as ulong
end type
constructor d3d()
x = 0
y = 0
z = 0
end constructor
constructor d3d( _
a as double , b as double , c as double )
x = a
y = b
z = c
end constructor
function d3d.tocolor() as ulong
return rgb( x * 255 , y * 255 , z * 255 )
end function
operator + ( a as d3d , b as d3d ) as d3d
return d3d( a.x + b.x , a.y + b.y , a.z + b.z )
end operator
operator - ( a as d3d , b as d3d ) as d3d
return d3d( a.x - b.x , a.y - b.y , a.z - b.z )
end operator
operator * ( a as d3d , b as double ) as d3d
return d3d( a.x * b , a.y * b , a.z * b )
end operator
operator / ( a as d3d , b as double ) as d3d
return d3d( a.x / b , a.y / b , a.z / b )
end operator
operator \ ( a as d3d , b as d3d ) as d3d
return d3d( a.y * b.x - a.z * b.y _
, a.z * b.x - a.x * b.z _
, a.x * b.z - a.y * b.x )
end operator
function dot( a as d3d , b as d3d ) as double
return a.x * b.x + a.y * b.y + a.z * b.z
end function
function lenght( a as d3d ) as double
return sqr( dot( a , a ) )
end function
function getangle( a as d3d , b as d3d ) as double
return acos( dot( a , b ) _
/ ( lenght( a ) * lenght( b ) ) )
end function

type m33
dim as double e11 , e12 , e13 _
, e21 , e22 , e23 _
, e31 , e32 , e33
declare constructor()
declare constructor( _
a11 as double , a12 as double , a13 as double _
, a21 as double , a22 as double , a23 as double _
, a31 as double , a32 as double , a33 as double )
declare function inverse() as m33
end type
constructor m33()
e11 = 1
e12 = 0
e13 = 0
e21 = 0
e22 = 1
e23 = 0
e31 = 0
e32 = 0
e33 = 1
end constructor
constructor m33( _
a11 as double , a12 as double , a13 as double _
, a21 as double , a22 as double , a23 as double _
, a31 as double , a32 as double , a33 as double )
e11 = a11
e12 = a12
e13 = a13
e21 = a21
e22 = a22
e23 = a23
e31 = a31
e32 = a32
e33 = a33
end constructor
operator * ( a as m33 , b as m33 ) as m33
return m33( _
a.e11 * b.e11 + a.e12 * b.e21 + a.e13 * b.e31 _
, a.e11 * b.e12 + a.e12 * b.e22 + a.e13 * b.e32 _
, a.e11 * b.e13 + a.e12 * b.e23 + a.e13 * b.e33 _
, a.e21 * b.e11 + a.e22 * b.e21 + a.e23 * b.e31 _
, a.e21 * b.e12 + a.e22 * b.e22 + a.e23 * b.e32 _
, a.e21 * b.e13 + a.e22 * b.e23 + a.e23 * b.e33 _
, a.e31 * b.e11 + a.e32 * b.e21 + a.e33 * b.e31 _
, a.e31 * b.e12 + a.e32 * b.e22 + a.e33 * b.e32 _
, a.e31 * b.e13 + a.e32 * b.e23 + a.e33 * b.e33 )
end operator
function m33.inverse() as m33
dim as double d = e11 * e22 * e33 _
- e11 * e32 * e23 _
+ e21 * e32 * e13 _
- e21 * e12 * e33 _
+ e31 * e12 * e23 _
- e31 * e22 * e13
if d = 0 then d = 1
return m33(  ( e22 * e33 - e23 * e32 ) / d _
, -( e12 * e33 - e13 * e32 ) / d _
,  ( e12 * e23 - e13 * e22 ) / d _
, -( e21 * e33 - e23 * e31 ) / d _
,  ( e11 * e33 - e13 * e31 ) / d _
, -( e11 * e23 - e13 * e21 ) / d _
,  ( e21 * e32 - e22 * e31 ) / d _
, -( e11 * e32 - e12 * e31 ) / d _
,  ( e11 * e22 - e12 * e21 ) / d )
end function

operator * ( m as m33 , v as d3d ) as d3d
return d3d( _
m.e11 * v.x + m.e12 * v.y + m.e13 * v.z _
, m.e21 * v.x + m.e22 * v.y + m.e23 * v.z _
, m.e31 * v.x + m.e32 * v.y + m.e33 * v.z )
end operator

type mplek
dim as m33 matrix
dim as d3d plek
end type
operator * ( m as mplek , v as d3d ) as d3d
return m.matrix * v + m.plek
end operator
operator * ( a as mplek , b as mplek ) as mplek
dim as mplek uit
uit.matrix = a.matrix * b.matrix
uit.plek = a * b.plek
return uit
end operator

dim shared as mplek mv( 20 )
dim shared as integer matrix_number

sub spot( byref x as double _
, byref y as double _
, byref z as double )
dim as d3d temp = d3d( x , y , z )
temp = mv( matrix_number ) * temp
x = temp.x
y = temp.y
z = temp.z
end sub

const as integer xyz = 0
const as integer xzy = 1
const as integer yxz = 2
const as integer yzx = 3
const as integer zxy = 4
const as integer zyx = 5
const as double pi = atn( 1 ) * 4
function rad( x as double ) as double
return x * pi / 180
end function
function pend( f as double , a as double ) as double
''pendular motion for animation
''looks verry natural
return sin( rad( f ) ) * a
end function
dim shared as double sk( 100 , 2 )
sub skelet( lim as integer , x as double , y as double , z as double )
''set angles of skeletal lim
''for animated avatars
if lim < 0 or lim > 64 then exit sub
sk( lim , 0 ) = x
sk( lim , 1 ) = y
sk( lim , 2 ) = z
end sub
dim shared as double cam( 6 )
sub camera( x as double , y as double , z as double _
, pan as double , tilt as double , rol as double , zoom as double )
''set look from point & look angles
cam( 0 ) = x
cam( 1 ) = y
cam( 2 ) = z
cam( 3 ) = pan
cam( 4 ) = tilt
cam( 5 ) = rol
cam( 6 ) = zoom
end sub
sub rotate( byref k as double , byref l as double , deg as double )
''help sub for rotating coordinates
dim as double hk , hl , s , c
s = sin( rad( deg ) )
c = cos( rad( deg ) )
hk = k * c - l * s
hl = k * s + l * c
k = hk
l = hl
end sub
sub moveCamera( x as double , y as double , z as double _
, pan as double )
rotate x , z , cam( 3 )
cam( 0 ) += x
cam( 1 ) += y
cam( 2 ) += z
cam( 3 ) = ( cam( 3 ) + pan ) mod 360
end sub
Sub link( no As Integer, x As double, y As double, z As double, pan As double, tilt As double, rol As double, ax As Integer, p As Integer )
''set curent matrix wil afect folowing drawing comands
If no < 1 Or no > 20 Then Exit Sub
If p < 0 Or p > 20 Then Exit Sub
If p = no Then Exit Sub
''create some lokal matrix's and fill them
Dim As mplek mp , rotx , roty , rotz , translate
mp = mv( p )
rotz.matrix.e11 = Cos( rad( rol ))
rotz.matrix.e12 = -Sin( rad( rol ))
rotz.matrix.e21 = Sin( rad( rol ))
rotz.matrix.e22 = Cos( rad( rol ))
roty.matrix.e11 = Cos( rad( pan ))
roty.matrix.e13 = -Sin( rad( pan ))
roty.matrix.e31 = Sin( rad( pan ))
roty.matrix.e33 = Cos( rad( pan ))
rotx.matrix.e22 = Cos( rad( tilt ))
rotx.matrix.e23 = -Sin( rad( tilt ))
rotx.matrix.e32 = Sin( rad( tilt ))
rotx.matrix.e33 = Cos( rad( tilt ))
translate.plek.x = x
translate.plek.y = y
translate.plek.z = z
''angles can permutate 6 ways
Select Case ax
Case xyz
mv( no ) = rotx * roty * rotz * translate * mp
Case xzy
mv( no ) = rotx * rotz * roty * translate * mp
Case yxz
mv( no ) = roty * rotx * rotz * translate * mp
Case yzx
mv( no ) = roty * rotz * rotx * translate * mp
Case zxy
mv( no ) = rotz * rotx * roty * translate * mp
Case zyx
mv( no ) = rotz * roty * rotx * translate * mp
Case Else
End Select
matrix_number = no
End Sub
sub child( no as integer , x as double , y as double , z as double _
, lim as integer , ax as integer , p as integer )
''set curent matrix for lim of animated avatar
''wil efect folowing drawings
if lim < 0 or lim > 64 then exit sub
link no , x , y , z _
, sk( lim , 1 ) , sk( lim , 0 ) , sk( lim , 2 ), ax , p
end sub
sub lokal( byref x as double , byref y as double _
, byref z as double , q as m33 )
dim as d3d temp = d3d( x , y , z )
temp = q * temp
x = temp.x
y = temp.y
z = temp.z
end sub

type tmaterial
dim as d3d diffuse( 10 )
dim as double grens( 11 )
dim as integer kleurmax
declare constructor ()
end type
constructor tmaterial()
dim as integer i
for i = 0 to 11
diffuse( i ) = d3d( 1 , 1 , 1 )
grens( i ) = 1
next i
grens( 0 ) = 0
kleurmax = 1
end constructor

dim shared as tmaterial material

dim shared as integer spheretel
const as integer spheremax = 100
type tSphere
dim as d3d center
dim as tmaterial mat
dim as m33 mp
declare sub fill( c as d3d , r as double )
declare function hit( o as d3d , d as d3d ) as double
end type

sub tSphere.fill( c as d3d , r as double )
spot c.x , c.y , c.z
center = c
mat = material
mp = mv( matrix_number ).matrix.inverse()
end sub

const as double infinity = 1e9

function tSphere.hit( o as d3d , d as d3d ) as double
dim as double t , a , b , c , disc
dim as d3d temp = o - center
a = dot( d , d )
b = 2 * dot( temp , d )
c = dot( temp , temp ) - radius2
disc = b ^ 2 - 4 * a * c
if disc < 0 then
return infinity
else
dim as double e = sqr( disc )
dim as double demon = 2 * a
t = ( -b - e ) / demon
if t > 1e-12 then
return t
end if

t = ( -b + e ) / demon
if t > 1e-12 then
return t
end if
end if
return infinity
end function

dim shared as d3d black   = d3d( 0 , 0 , 0 )
dim shared as d3d red     = d3d( 1 , 0 , 0 )
dim shared as d3d green   = d3d( 0 , 1 , 0 )
dim shared as d3d yellow  = d3d( 1 , 1 , 0 )
dim shared as d3d blue    = d3d( 0 , 0 , 1 )
dim shared as d3d magenta = d3d( 1 , 0 , 1 )
dim shared as d3d cyan    = d3d( 0 , 1 , 1 )
dim shared as d3d white   = d3d( 1 , 1 , 1 )

dim shared as tSphere spheres( spheremax )

function pixel( o as d3d , d as d3d ) as d3d
dim as integer i , klno , isph = -1
dim as double dist , slow = infinity
for i = 0 to spheretel
dist = spheres( i ).hit( o , d )
if dist < slow then
slow = dist
isph = i
end if
next i
if isph = -1 then return black
dim as tsphere bol = spheres( isph )
dim as d3d n , p , plek
p = o + d  * slow
n = p - bol.center
plek = n
lokal plek.x , plek.y , plek.z , bol.mp
dim as double dd = turbulence3d( _
, plek.z + bol.radius , 64 )
klno = 0
for i = 0 to bol.mat.kleurmax
if ( bol.mat.grens( i ) <= dd _
and bol.mat.grens( i + 1 ) > dd ) then
klno = i
end if
next i
dim as d3d kll = bol.mat.diffuse( klno )
dim as d3d klh = bol.mat.diffuse( klno + 1 )
dim as double gl = bol.mat.grens( klno )
dim as double gh = bol.mat.grens( klno + 1 )
dim as double f = ( dd - gl ) / ( gh - gl )
dim as d3d kl = kll + ( klh - kll ) * f
dim as double a = getangle( n , d3d(-1,1,0) )
return kl * ( cos( a ) / 2 + .5 )
end function
sub test
material.grens( 0 ) = 0
material.diffuse( 0 ) = blue
material.grens( 1 ) = .4
material.diffuse( 1 ) = cyan
material.grens( 2 ) = .45
material.diffuse( 2 ) = yellow
material.grens( 3 ) = .5
material.diffuse( 3 ) = green
material.grens( 4 ) = .7
material.diffuse( 4 ) = white / 2
material.grens( 5 ) = .8
material.diffuse( 5 ) = white
material.kleurmax = 6
link 1 , 0,0,0 , -30,0,0 , xyz , 0
spheres( 0 ).fill d3d(100,0,0) , 100
material.diffuse( 0 ) = black
material.grens( 1 ) = 1
material.diffuse( 1 ) = white
material.kleurmax = 1
link 1 , 0,0,0 , 30,0,0 , xyz , 0
spheres( 1 ).fill d3d(-100,0,0) , 50
spheretel = 1
dim as d3d o , d
dim as double x , y
for x = -winx / 2 to winx / 2
for y = -winy / 2 to winy / 2
o = d3d( 0 , 0 , -1000 )
d = d3d( x , y , 1000 )
pset ( winx / 2 + x , winy / 2 - y ) _
, pixel( o , d ).tocolor()
next y
next x
end sub
cls
paint ( 1 , 1 ) , rgb( 127 , 127 , 127 )
test
sleep

Code: Select all

'' bluatigro 29 nov 2017
'' noise

#ifndef __SIMPLENOISE_BI__
#define __SIMPLENOISE_BI__

dim shared as integer NOIZESIZE = 63

function noise( x as integer , y as integer , z as integer ) as double
randomize x + 1000 * y + 1e6 * z
return rnd
end function

function smoothNoise1d(x as double) as double
x=abs(x)
dim as integer iX1=int(x)
dim as integer iX2=ix1+1
dim as double tx=x-iX1
iX1 and= NOIZESIZE
iX2 and= NOIZESIZE
dim as double l=noise(iX1,0,0)
dim as double r=noise(iX2,0,0)
dim as double v=l + (r-l)*tx
return v
end function

function smoothNoise2d(x as double, y as double) as double
x=abs(x)
y=abs(y)
dim as integer iX1=int(x)
dim as integer iY1=int(y)
dim as integer iX2=ix1+1
dim as integer iY2=iy1+1
dim as double tx=x-iX1
dim as double ty=y-iY1

iX1 and= NOIZESIZE
iX2 and= NOIZESIZE
iY1 and= NOIZESIZE
iY2 and= NOIZESIZE

dim as double lt=noise(iX1,iY1,0)
dim as double rt=noise(iX2,iY1,0)
dim as double rb=noise(iX2,iY2,0)
dim as double lb=noise(iX1,iY2,0)

dim as double sxt=lt + (rt-lt)*tx
dim as double sxb=lb + (rb-lb)*tx

dim as double v=sxt+(sxb-sxt)*ty
return v
end function

function smoothNoise3d(x as double, y as double, z as double) as double
x=abs(x)
y=abs(y)
z=abs(z)
dim as integer iX1=int(x)
dim as integer iY1=int(y)
dim as integer iZ1=int(z)
dim as integer iX2=ix1+1
dim as integer iY2=iy1+1
dim as integer iZ2=iz1+1
dim as double tx=x-iX1
dim as double ty=y-iY1
dim as double tz=z-iZ1

iX1 and= NOIZESIZE
iX2 and= NOIZESIZE
iY1 and= NOIZESIZE
iY2 and= NOIZESIZE
iZ1 and= NOIZESIZE
iZ2 and= NOIZESIZE

dim as double ltf=noise(iX1,iY1,iZ1)
dim as double rtf=noise(iX2,iY1,iZ1)
dim as double rbf=noise(iX2,iY2,iZ1)
dim as double lbf=noise(iX1,iY2,iZ1)
dim as double sxtf=ltf + (rtf-ltf)*tx
dim as double sxbf=lbf + (rbf-lbf)*tx

dim as double ltb=noise(iX1,iY1,iZ2)
dim as double rtb=noise(iX2,iY1,iZ2)
dim as double rbb=noise(iX2,iY2,iZ2)
dim as double lbb=noise(iX1,iY2,iZ2)
dim as double sxtb = ltb + (rtb-ltb)*tx
dim as double sxbb = lbb + (rbb-lbb)*tx

dim as double vf = sxtf+(sxbf-sxtf)*ty
dim as double vb = sxtb+(sxbb-sxtb)*ty
dim as double v = vf + (vb-vf)*tz

return v
end function

function turbulence1d(x as double, size as double) as double
dim as double value, initialSize=any
if size<1 then size=1
initialSize = size
while(size >= 1)
value += smoothNoise1d(x / size) * size
size /= 2.0
wend
value*=0.5
return value/initialSize
end function

function turbulence2d(x as double, y as double, size as double) as double
dim as double value, initialSize=any
if size<1 then size=1
initialSize = size
while(size >= 1)
value += smoothNoise2d(x / size, y / size) * size
size /= 2.0
wend
value*=0.5
return value/(initialSize-1)
end function

function turbulence3d( _
x as double , y as double , z as double _
, size as double ) as double
dim as double value , initialSize = any
if size < 1 then size = 1
initialSize = size
while( size >= 1 )
value += smoothNoise3d(x / size, y / size, z/size) * size
size /= 2.0
wend
value *= 0.5
return value / ( initialSize - 1 )
end function

#endif ' __SIMPLENOISE_BI__
bluatigro
Posts: 651
Joined: Apr 25, 2012 10:35
Location: netherlands

### Re: raytrace try

update :

now we got :
earth , moon , mars , Jupiter and neptune

Code: Select all

'' bluatigro 16 nov 2018
'' bluaray 1 material

#include "noise.bas"

dim shared as integer winx , winy
screen 18 , 32
screeninfo winx , winy
type d3d
dim as double x , y , z
declare constructor ()
declare constructor ( _
a as double , b as double , c as double )
declare function tocolor() as ulong
end type
constructor d3d()
x = 0
y = 0
z = 0
end constructor
constructor d3d( _
a as double , b as double , c as double )
x = a
y = b
z = c
end constructor
function d3d.tocolor() as ulong
return rgb( x * 255 , y * 255 , z * 255 )
end function
operator + ( a as d3d , b as d3d ) as d3d
return d3d( a.x + b.x , a.y + b.y , a.z + b.z )
end operator
operator - ( a as d3d , b as d3d ) as d3d
return d3d( a.x - b.x , a.y - b.y , a.z - b.z )
end operator
operator * ( a as d3d , b as double ) as d3d
return d3d( a.x * b , a.y * b , a.z * b )
end operator
operator / ( a as d3d , b as double ) as d3d
return d3d( a.x / b , a.y / b , a.z / b )
end operator
operator \ ( a as d3d , b as d3d ) as d3d
return d3d( a.y * b.x - a.z * b.y _
, a.z * b.x - a.x * b.z _
, a.x * b.z - a.y * b.x )
end operator
function dot( a as d3d , b as d3d ) as double
return a.x * b.x + a.y * b.y + a.z * b.z
end function
function lenght( a as d3d ) as double
return sqr( dot( a , a ) )
end function
function getangle( a as d3d , b as d3d ) as double
return acos( dot( a , b ) _
/ ( lenght( a ) * lenght( b ) ) )
end function

type m33
dim as double e11 , e12 , e13 _
, e21 , e22 , e23 _
, e31 , e32 , e33
declare constructor()
declare constructor( _
a11 as double , a12 as double , a13 as double _
, a21 as double , a22 as double , a23 as double _
, a31 as double , a32 as double , a33 as double )
declare function inverse() as m33
end type
constructor m33()
e11 = 1
e12 = 0
e13 = 0
e21 = 0
e22 = 1
e23 = 0
e31 = 0
e32 = 0
e33 = 1
end constructor
constructor m33( _
a11 as double , a12 as double , a13 as double _
, a21 as double , a22 as double , a23 as double _
, a31 as double , a32 as double , a33 as double )
e11 = a11
e12 = a12
e13 = a13
e21 = a21
e22 = a22
e23 = a23
e31 = a31
e32 = a32
e33 = a33
end constructor
operator * ( a as m33 , b as m33 ) as m33
return m33( _
a.e11 * b.e11 + a.e12 * b.e21 + a.e13 * b.e31 _
, a.e11 * b.e12 + a.e12 * b.e22 + a.e13 * b.e32 _
, a.e11 * b.e13 + a.e12 * b.e23 + a.e13 * b.e33 _
, a.e21 * b.e11 + a.e22 * b.e21 + a.e23 * b.e31 _
, a.e21 * b.e12 + a.e22 * b.e22 + a.e23 * b.e32 _
, a.e21 * b.e13 + a.e22 * b.e23 + a.e23 * b.e33 _
, a.e31 * b.e11 + a.e32 * b.e21 + a.e33 * b.e31 _
, a.e31 * b.e12 + a.e32 * b.e22 + a.e33 * b.e32 _
, a.e31 * b.e13 + a.e32 * b.e23 + a.e33 * b.e33 )
end operator
function m33.inverse() as m33
dim as double d = e11 * e22 * e33 _
- e11 * e32 * e23 _
+ e21 * e32 * e13 _
- e21 * e12 * e33 _
+ e31 * e12 * e23 _
- e31 * e22 * e13
if d = 0 then d = 1
return m33(  ( e22 * e33 - e23 * e32 ) / d _
, -( e12 * e33 - e13 * e32 ) / d _
,  ( e12 * e23 - e13 * e22 ) / d _
, -( e21 * e33 - e23 * e31 ) / d _
,  ( e11 * e33 - e13 * e31 ) / d _
, -( e11 * e23 - e13 * e21 ) / d _
,  ( e21 * e32 - e22 * e31 ) / d _
, -( e11 * e32 - e12 * e31 ) / d _
,  ( e11 * e22 - e12 * e21 ) / d )
end function

operator * ( m as m33 , v as d3d ) as d3d
return d3d( _
m.e11 * v.x + m.e12 * v.y + m.e13 * v.z _
, m.e21 * v.x + m.e22 * v.y + m.e23 * v.z _
, m.e31 * v.x + m.e32 * v.y + m.e33 * v.z )
end operator

type mplek
dim as m33 matrix
dim as d3d plek
end type
operator * ( m as mplek , v as d3d ) as d3d
return m.matrix * v + m.plek
end operator
operator * ( a as mplek , b as mplek ) as mplek
dim as mplek uit
uit.matrix = a.matrix * b.matrix
uit.plek = a * b.plek
return uit
end operator

dim shared as mplek mv( 20 )
dim shared as integer matrix_number

sub spot( byref x as double _
, byref y as double _
, byref z as double )
dim as d3d temp = d3d( x , y , z )
temp = mv( matrix_number ) * temp
x = temp.x
y = temp.y
z = temp.z
end sub

const as integer xyz = 0
const as integer xzy = 1
const as integer yxz = 2
const as integer yzx = 3
const as integer zxy = 4
const as integer zyx = 5
const as double pi = atn( 1 ) * 4
function rad( x as double ) as double
return x * pi / 180
end function
function pend( f as double , a as double ) as double
''pendular motion for animation
''looks verry natural
return sin( rad( f ) ) * a
end function
dim shared as double sk( 100 , 2 )
sub skelet( lim as integer , x as double , y as double , z as double )
''set angles of skeletal lim
''for animated avatars
if lim < 0 or lim > 64 then exit sub
sk( lim , 0 ) = x
sk( lim , 1 ) = y
sk( lim , 2 ) = z
end sub
dim shared as double cam( 6 )
sub camera( x as double , y as double , z as double _
, pan as double , tilt as double , rol as double , zoom as double )
''set look from point & look angles
cam( 0 ) = x
cam( 1 ) = y
cam( 2 ) = z
cam( 3 ) = pan
cam( 4 ) = tilt
cam( 5 ) = rol
cam( 6 ) = zoom
end sub
sub rotate( byref k as double , byref l as double , deg as double )
''help sub for rotating coordinates
dim as double hk , hl , s , c
s = sin( rad( deg ) )
c = cos( rad( deg ) )
hk = k * c - l * s
hl = k * s + l * c
k = hk
l = hl
end sub
sub moveCamera( x as double , y as double , z as double _
, pan as double )
rotate x , z , cam( 3 )
cam( 0 ) += x
cam( 1 ) += y
cam( 2 ) += z
cam( 3 ) = ( cam( 3 ) + pan ) mod 360
end sub
Sub link( no As Integer, x As double, y As double, z As double, pan As double, tilt As double, rol As double, ax As Integer, p As Integer )
''set curent matrix wil afect folowing drawing comands
If no < 1 Or no > 20 Then Exit Sub
If p < 0 Or p > 20 Then Exit Sub
If p = no Then Exit Sub
''create some lokal matrix's and fill them
Dim As mplek mp , rotx , roty , rotz , translate
mp = mv( p )
rotz.matrix.e11 = Cos( rad( rol ))
rotz.matrix.e12 = -Sin( rad( rol ))
rotz.matrix.e21 = Sin( rad( rol ))
rotz.matrix.e22 = Cos( rad( rol ))
roty.matrix.e11 = Cos( rad( pan ))
roty.matrix.e13 = -Sin( rad( pan ))
roty.matrix.e31 = Sin( rad( pan ))
roty.matrix.e33 = Cos( rad( pan ))
rotx.matrix.e22 = Cos( rad( tilt ))
rotx.matrix.e23 = -Sin( rad( tilt ))
rotx.matrix.e32 = Sin( rad( tilt ))
rotx.matrix.e33 = Cos( rad( tilt ))
translate.plek.x = x
translate.plek.y = y
translate.plek.z = z
''angles can permutate 6 ways
Select Case ax
Case xyz
mv( no ) = rotx * roty * rotz * translate * mp
Case xzy
mv( no ) = rotx * rotz * roty * translate * mp
Case yxz
mv( no ) = roty * rotx * rotz * translate * mp
Case yzx
mv( no ) = roty * rotz * rotx * translate * mp
Case zxy
mv( no ) = rotz * rotx * roty * translate * mp
Case zyx
mv( no ) = rotz * roty * rotx * translate * mp
Case Else
End Select
matrix_number = no
End Sub
sub child( no as integer , x as double , y as double , z as double _
, lim as integer , ax as integer , p as integer )
''set curent matrix for lim of animated avatar
''wil efect folowing drawings
if lim < 0 or lim > 64 then exit sub
link no , x , y , z _
, sk( lim , 1 ) , sk( lim , 0 ) , sk( lim , 2 ), ax , p
end sub
sub lokal( byref x as double , byref y as double _
, byref z as double , q as m33 )
dim as d3d temp = d3d( x , y , z )
temp = q * temp
x = temp.x
y = temp.y
z = temp.z
end sub
const as integer turbulence = 1
const as integer wave = 2

dim shared as d3d black   = d3d( 0 , 0 , 0 )
dim shared as d3d red     = d3d( 1 , 0 , 0 )
dim shared as d3d green   = d3d( 0 , 1 , 0 )
dim shared as d3d yellow  = d3d( 1 , 1 , 0 )
dim shared as d3d blue    = d3d( 0 , 0 , 1 )
dim shared as d3d magenta = d3d( 1 , 0 , 1 )
dim shared as d3d cyan    = d3d( 0 , 1 , 1 )
dim shared as d3d white   = d3d( 1 , 1 , 1 )

type tmaterial
dim as d3d diffuse( 10 ) , scale
dim as double grens( 11 ) , turbsize , wavesize
dim as integer kleurmax , mattype
declare constructor ()
declare sub filldiffuse( a as d3d , b as d3d , c as d3d , d as d3d _
, e as d3d , f as d3d , g as d3d , h as d3d )
declare sub fillgrens( a as double , b as double , c as double , d as double _
, e as double , f as double , g as double , h as double )
end type
constructor tmaterial()
dim as integer i
for i = 0 to 10
diffuse( i ) = d3d( 1 , 1 , 1 )
grens( i + 1 ) = 1
next i
kleurmax = 1
scale = d3d( 1 , 1 , 1 )
mattype = turbulence
turbsize = 0
wavesize = 0
end constructor
sub tmaterial.filldiffuse( a as d3d , b as d3d , c as d3d , d as d3d _
, e as d3d , f as d3d , g as d3d , h as d3d )
diffuse( 0 ) = a
diffuse( 1 ) = b
diffuse( 2 ) = c
diffuse( 3 ) = d
diffuse( 4 ) = e
diffuse( 5 ) = f
diffuse( 6 ) = g
diffuse( 7 ) = h
end sub
sub tmaterial.fillgrens( a as double , b as double , c as double , d as double _
, e as double , f as double , g as double , h as double )
grens( 0 ) = 0
grens( 1 ) = a
grens( 2 ) = b
grens( 3 ) = c
grens( 4 ) = d
grens( 5 ) = e
grens( 6 ) = f
grens( 7 ) = g
grens( 8 ) = h
end sub

dim shared as tmaterial material

dim shared as integer spheretel
const as integer spheremax = 100
type tSphere
dim as d3d center
dim as tmaterial mat
dim as m33 mp
declare sub fill( c as d3d , r as double )
declare function hit( o as d3d , d as d3d ) as double
end type

sub tSphere.fill( c as d3d , r as double )
spot c.x , c.y , c.z
center = c
mat = material
mp = mv( matrix_number ).matrix.inverse()
end sub

const as double infinity = 1e9

function tSphere.hit( o as d3d , d as d3d ) as double
dim as double t , a , b , c , disc
dim as d3d temp = o - center
a = dot( d , d )
b = 2 * dot( temp , d )
c = dot( temp , temp ) - radius2
disc = b ^ 2 - 4 * a * c
if disc < 0 then
return infinity
else
dim as double e = sqr( disc )
dim as double demon = 2 * a
t = ( -b - e ) / demon
if t > 1e-12 then
return t
end if

t = ( -b + e ) / demon
if t > 1e-12 then
return t
end if
end if
return infinity
end function

dim shared as tSphere spheres( spheremax )
sub addsphere( x as double , y as double , z as double , d as double , m as tmaterial )
material = m
spheres( spheretel ).fill d3d( x , y , z  ) , d
spheretel += 1
end sub

function pixel( o as d3d , d as d3d ) as d3d
dim as integer i , klno , isph = -1
dim as double dist , slow = infinity
for i = 0 to spheretel
dist = spheres( i ).hit( o , d )
if dist < slow then
slow = dist
isph = i
end if
next i
if isph = -1 then return black
dim as tsphere bol = spheres( isph )
dim as d3d n , p , plek
p = o + d  * slow
n = p - bol.center
plek = n
lokal plek.x , plek.y , plek.z , bol.mp
dim as double dd
if bol.mat.mattype and turbulence then
dd = turbulence3d( _
( plek.x + bol.radius ) / bol.mat.scale.x _
, ( plek.y + bol.radius ) / bol.mat.scale.y _
, ( plek.z + bol.radius ) / bol.mat.scale.z , 64 )
dd = dd * bol.mat.turbsize
end if
if bol.mat.mattype and wave then
dd = sin( plek.y * bol.mat.wavesize + dd ) / 2 + .5
end if
klno = 0
for i = 0 to bol.mat.kleurmax
if ( bol.mat.grens( i ) <= dd _
and bol.mat.grens( i + 1 ) > dd ) then
klno = i
end if
next i
dim as d3d kll = bol.mat.diffuse( klno )
dim as d3d klh = bol.mat.diffuse( klno + 1 )
dim as double gl = bol.mat.grens( klno )
dim as double gh = bol.mat.grens( klno + 1 )
dim as double f = ( dd - gl ) / ( gh - gl )
dim as d3d kl = kll + ( klh - kll ) * f
dim as double a = getangle( n , d3d(-1,1,0) )
return kl * ( cos( a ) / 2 + .5 )
end function
dim shared as tmaterial earth , moon , mars , jupiter , neptune
earth.filldiffuse blue , cyan , yellow , green , white / 2 , white , red , red
earth.fillgrens     .5 , .6   , .7     , .8    , .9        , 1     ,  1  , 1
earth.kleurmax = 6
earth.scale = d3d( 1/3 , 1/3 , 1/3 )
earth.mattype = turbulence
earth.turbsize = 1
earth.wavesize = 0

moon.grens( 0 ) = 0
moon.filldiffuse black , white , red , red ,red , red , red , red
moon.fillgrens 1,1,1,1,1,1,1,1
moon.kleurmax = 1
moon.scale = d3d( .1 , .1 , .1 )
moon.mattype = turbulence
moon.turbsize = 1
moon.wavesize = 0

mars.filldiffuse red , yellow , red , red , red , red , red , red
mars.fillgrens 1,1,1,1,1,1,1,1
mars.kleurmax = 1
mars.scale = d3d( .1 , .1 , .1 )
mars.mattype = turbulence
mars.turbsize = 1
mars.wavesize = 0

jupiter.filldiffuse red / 2 , red , yellow , white / 2 , white , white , red , red
jupiter.fillgrens .5,.6,.7,.9,1,1,1,1
jupiter.kleurmax = 5
jupiter.scale = d3d( 1 , .5 , 1 )
jupiter.mattype = turbulence or wave
jupiter.wavesize = 1 / 10
jupiter.turbsize = 10

neptune.grens( 0 ) = 0
neptune.filldiffuse blue , white , white , white , red , red , red , red
neptune.fillgrens .9 , 1 ,1,1,1,1,1,1
neptune.kleurmax = 1
neptune.scale = d3d( 1 , 1 , 1 )
neptune.mattype = turbulence ''or wave
neptune.wavesize = 0 ''1 / 5
neptune.turbsize = 1 ''0
sub test
spheretel = 0
link 1 , -200,100,0 , 0,0,0 , xyz , 0
addsphere 0,0,0 , 100 , jupiter
link 1 , 200,100,0 , 0,0,0 , xyz , 0
addsphere 0,0,0 , 20 , mars
link 1 , 200,-100,0 , 0,0,0 , xyz , 0
addsphere 0,0,0 , 15 , moon
link 1 , -200,-100,0 , 0,0,0 , xyz , 0
addsphere 0,0,0 , 30 , earth
link 1 , 0,100,0 , 0,0,0 , xyz , 0
addsphere 0,0,0 , 70 , neptune

dim as d3d o , d
dim as double x , y
for x = -winx / 2 to winx / 2
for y = -winy / 2 to winy / 2
o = d3d( 0 , 0 , -1000 )
d = d3d( x , y , 1000 )
pset ( winx / 2 + x , winy / 2 - y ) _
, pixel( o , d ).tocolor()
next y
next x
end sub
cls
paint ( 1 , 1 ) , rgb( 127 , 127 , 127 )
test
bsave "5-planets.bmp" , 0
sleep

dafhi
Posts: 1357
Joined: Jun 04, 2005 9:51

### Re: raytrace try

nice!

with the help of a video by inigo quilez "iq", and my own skills, i made an optimized sphere hit.

but to use it, you'll need normalize

Code: Select all

declare property norm as d3d '' or change to function
..
property d3d.norm() as d3d
var s = 1 / sqr( x * x + y * y + z * z )
return d3d( s * x , s * y , s * z )
end property
...
o = d3d( 0 , 0 , -1000 )
d = d3d( x , y , 1000 ).norm

you may eliminate sphere hit "A"

D.J.Peters
Posts: 8132
Joined: May 28, 2005 3:28
Contact:

### Re: raytrace try

Nice simplex noise textures.

Joshy
paul doe
Posts: 1255
Joined: Jul 25, 2017 17:22
Location: Argentina

### Re: raytrace try

Nice. Another simple optimization that helps shaving a few cycles: instead of dividing by the determinant, you can multiply by its reciprocal, like this:

Code: Select all

'' ...
function m33.inverse() as m33
dim as double d = e11 * e22 * e33 _
- e11 * e32 * e23 _
+ e21 * e32 * e13 _
- e21 * e12 * e33 _
+ e31 * e12 * e23 _
- e31 * e22 * e13
if d = 0 then d = 1

'' Use the reciprocal of the determinant
dim as double invD = 1.0 / d

'return m33(  ( e22 * e33 - e23 * e32 ) / d _
'          , -( e12 * e33 - e13 * e32 ) / d _
'          ,  ( e12 * e23 - e13 * e22 ) / d _
'          , -( e21 * e33 - e23 * e31 ) / d _
'          ,  ( e11 * e33 - e13 * e31 ) / d _
'          , -( e11 * e23 - e13 * e21 ) / d _
'          ,  ( e21 * e32 - e22 * e31 ) / d _
'          , -( e11 * e32 - e12 * e31 ) / d _
'          ,  ( e11 * e22 - e12 * e21 ) / d )

return m33(  ( e22 * e33 - e23 * e32 ) * invD _
, -( e12 * e33 - e13 * e32 ) * invD _
,  ( e12 * e23 - e13 * e22 ) * invD _
, -( e21 * e33 - e23 * e31 ) * invD _
,  ( e11 * e33 - e13 * e31 ) * invD _
, -( e11 * e23 - e13 * e21 ) * invD _
,  ( e21 * e32 - e22 * e31 ) * invD _
, -( e11 * e32 - e12 * e31 ) * invD _
,  ( e11 * e22 - e12 * e21 ) * invD )
end function
'' ...

9 multiplications vs 9 divisions always help ;)
dodicat
Posts: 6645
Joined: Jan 10, 2006 20:30
Location: Scotland

### Re: raytrace try

Really nice bluatigro.
bluatigro
Posts: 651
Joined: Apr 25, 2012 10:35
Location: netherlands

### Re: raytrace try

update :
the 3Dengine is now good

you can safely link lim's to avatar's again

Code: Select all

'' bluatigro 7 dec 2018
'' bluaray 2 material

#include "noise.bas"

dim shared as integer winx , winy
screen 20 , 32
screeninfo winx , winy
type d3d
dim as double x , y , z
declare constructor ()
declare constructor ( _
a as double , b as double , c as double )
declare function tocolor() as ulong
end type
constructor d3d()
x = 0
y = 0
z = 0
end constructor
constructor d3d( _
a as double , b as double , c as double )
x = a
y = b
z = c
end constructor
function d3d.tocolor() as ulong
return rgb( x * 255 , y * 255 , z * 255 )
end function
operator + ( a as d3d , b as d3d ) as d3d
return d3d( a.x + b.x , a.y + b.y , a.z + b.z )
end operator
operator - ( a as d3d , b as d3d ) as d3d
return d3d( a.x - b.x , a.y - b.y , a.z - b.z )
end operator
operator * ( a as d3d , b as double ) as d3d
return d3d( a.x * b , a.y * b , a.z * b )
end operator
operator / ( a as d3d , b as double ) as d3d
return d3d( a.x / b , a.y / b , a.z / b )
end operator
operator \ ( a as d3d , b as d3d ) as d3d
return d3d( a.y * b.x - a.z * b.y _
, a.z * b.x - a.x * b.z _
, a.x * b.z - a.y * b.x )
end operator
function dot( a as d3d , b as d3d ) as double
return a.x * b.x + a.y * b.y + a.z * b.z
end function
function lenght( a as d3d ) as double
return sqr( dot( a , a ) )
end function
function getangle( a as d3d , b as d3d ) as double
return acos( dot( a , b ) _
/ ( lenght( a ) * lenght( b ) ) )
end function

type m33
dim as double e( 3 , 3 )
declare constructor()
declare constructor( _
a11 as double , a12 as double , a13 as double _
, a21 as double , a22 as double , a23 as double _
, a31 as double , a32 as double , a33 as double )
declare function inverse() as m33
end type
constructor m33()
e(0,0) = 1
e(0,1) = 0
e(0,2) = 0
e(0,3) = 0
e(1,0) = 0
e(1,1) = 1
e(1,2) = 0
e(1,3) = 0
e(2,0) = 0
e(2,1) = 0
e(2,2) = 1
e(2,3) = 0
e(3,0) = 0
e(3,1) = 0
e(3,2) = 0
e(3,3) = 1
end constructor
constructor m33( _
a11 as double , a12 as double , a13 as double _
, a21 as double , a22 as double , a23 as double _
, a31 as double , a32 as double , a33 as double )
e(0,0) = a11
e(0,1) = a12
e(0,2) = a13
e(0,3) = 0
e(1,0) = a21
e(1,1) = a22
e(1,2) = a23
e(1,3) = 0
e(2,0) = a31
e(2,1) = a32
e(2,2) = a33
e(2,3) = 0
e(3,0) = 0
e(3,1) = 0
e(3,2) = 0
e(3,3) = 1
end constructor
operator * ( a as m33 , b as m33 ) as m33
dim as m33 uit
dim as integer i , j , k
for i = 0 to 3
for j = 0 to 3
uit.e( i , j ) = 0
for k = 0 to 3
uit.e( i , j ) += a.e( i , k ) * b.e( k , j )
next k
next j
next i
return uit
end operator
function m33.inverse() as m33
dim as double d = e(0,0) * e(1,1) * e(2,2) _
- e(0,0) * e(2,1) * e(1,2) _
+ e(1,0) * e(2,1) * e(0,2) _
- e(1,0) * e(0,1) * e(2,2) _
+ e(2,0) * e(0,1) * e(1,2) _
- e(2,0) * e(1,1) * e(0,2)
if d = 0 then d = 1
d = 1 / d
return m33(  ( e(1,1) * e(2,2) - e(1,2) * e(2,1) ) * d _
, -( e(0,1) * e(2,2) - e(0,2) * e(2,1) ) * d _
,  ( e(0,1) * e(1,2) - e(0,2) * e(1,1) ) * d _
, -( e(1,0) * e(2,2) - e(1,2) * e(2,0) ) * d _
,  ( e(0,0) * e(2,2) - e(0,2) * e(2,0) ) * d _
, -( e(0,0) * e(1,2) - e(0,2) * e(1,0) ) * d _
,  ( e(1,0) * e(2,1) - e(1,1) * e(2,0) ) * d _
, -( e(0,0) * e(2,1) - e(0,1) * e(2,0) ) * d _
,  ( e(0,0) * e(1,1) - e(0,1) * e(1,0) ) * d )
end function

operator * ( m as m33 , v as d3d ) as d3d
return d3d( _
m.e(0,0) * v.x + m.e(0,1) * v.y + m.e(0,2) * v.z _
, m.e(1,0) * v.x + m.e(1,1) * v.y + m.e(1,2) * v.z _
, m.e(2,0) * v.x + m.e(2,1) * v.y + m.e(2,2) * v.z )
end operator

dim shared as m33 mv( 20 )
dim shared as integer matrix_number

sub spot( byref x as double _
, byref y as double _
, byref z as double )
dim as d3d temp = d3d( x , y , z )
temp = mv( matrix_number ) * temp
x = temp.x
y = temp.y
z = temp.z
end sub

const as integer xyz = 0
const as integer xzy = 1
const as integer yxz = 2
const as integer yzx = 3
const as integer zxy = 4
const as integer zyx = 5
const as double pi = atn( 1 ) * 4
function rad( x as double ) as double
return x * pi / 180
end function
function pend( f as double , a as double ) as double
''pendular motion for animation
''looks verry natural
return sin( rad( f ) ) * a
end function
dim shared as double sk( 100 , 2 )
sub skelet( lim as integer , x as double , y as double , z as double )
''set angles of skeletal lim
''for animated avatars
if lim < 0 or lim > 64 then exit sub
sk( lim , 0 ) = x
sk( lim , 1 ) = y
sk( lim , 2 ) = z
end sub
dim shared as double cam( 6 )
sub camera( x as double , y as double , z as double _
, pan as double , tilt as double , rol as double , zoom as double )
''set look from point & look angles
cam( 0 ) = x
cam( 1 ) = y
cam( 2 ) = z
cam( 3 ) = pan
cam( 4 ) = tilt
cam( 5 ) = rol
cam( 6 ) = zoom
end sub
sub rotate( byref k as double , byref l as double , deg as double )
''help sub for rotating coordinates
dim as double hk , hl , s , c
s = sin( rad( deg ) )
c = cos( rad( deg ) )
hk = k * c - l * s
hl = k * s + l * c
k = hk
l = hl
end sub
sub moveCamera( x as double , y as double , z as double _
, pan as double )
rotate x , z , cam( 3 )
cam( 0 ) += x
cam( 1 ) += y
cam( 2 ) += z
cam( 3 ) = ( cam( 3 ) + pan ) mod 360
end sub
Sub link( no As Integer, x As double, y As double, z As double, pan As double, tilt As double, rol As double, ax As Integer, p As Integer )
''set curent matrix wil afect folowing drawing comands
If no < 1 Or no > 20 Then Exit Sub
If p < 0 Or p > 20 Then Exit Sub
If p = no Then Exit Sub
''create some lokal matrix's and fill them
Dim As m33 mp , rotx , roty , rotz , translate
mp = mv( p )
rotz.e(0,0) = Cos( rad( rol ))
rotz.e(0,1) = -Sin( rad( rol ))
rotz.e(1,0) = Sin( rad( rol ))
rotz.e(1,1) = Cos( rad( rol ))
roty.e(0,0) = Cos( rad( pan ))
roty.e(0,2) = -Sin( rad( pan ))
roty.e(2,0) = Sin( rad( pan ))
roty.e(2,2) = Cos( rad( pan ))
rotx.e(1,1) = Cos( rad( tilt ))
rotx.e(1,2) = -Sin( rad( tilt ))
rotx.e(2,1) = Sin( rad( tilt ))
rotx.e(2,2) = Cos( rad( tilt ))
translate.e(3,0) = x
translate.e(3,1) = y
translate.e(3,2) = z
''angles can permutate 6 ways
Select Case ax
Case xyz
mv( no ) = rotx * roty * rotz * translate * mp
Case xzy
mv( no ) = rotx * rotz * roty * translate * mp
Case yxz
mv( no ) = roty * rotx * rotz * translate * mp
Case yzx
mv( no ) = roty * rotz * rotx * translate * mp
Case zxy
mv( no ) = rotz * rotx * roty * translate * mp
Case zyx
mv( no ) = rotz * roty * rotx * translate * mp
Case Else
End Select
matrix_number = no
End Sub
sub child( no as integer , x as double , y as double , z as double _
, lim as integer , ax as integer , p as integer )
''set curent matrix for lim of animated avatar
''wil efect folowing drawings
if lim < 0 or lim > 64 then exit sub
link no , x , y , z _
, sk( lim , 1 ) , sk( lim , 0 ) , sk( lim , 2 ), ax , p
end sub
sub lokal( byref x as double , byref y as double _
, byref z as double , q as m33 )
dim as d3d temp = d3d( x , y , z )
temp = q * temp
x = temp.x
y = temp.y
z = temp.z
end sub
const as integer turbulence = 1
const as integer wave = 2

dim shared as d3d black   = d3d( 0 , 0 , 0 )
dim shared as d3d red     = d3d( 1 , 0 , 0 )
dim shared as d3d green   = d3d( 0 , 1 , 0 )
dim shared as d3d yellow  = d3d( 1 , 1 , 0 )
dim shared as d3d blue    = d3d( 0 , 0 , 1 )
dim shared as d3d magenta = d3d( 1 , 0 , 1 )
dim shared as d3d cyan    = d3d( 0 , 1 , 1 )
dim shared as d3d white   = d3d( 1 , 1 , 1 )

type tmaterial
dim as d3d diffuse( 10 ) , scale
dim as double grens( 11 ) , turbsize , wavesize
dim as integer kleurmax , mattype
declare constructor ()
declare sub filldiffuse( a as d3d , b as d3d , c as d3d , d as d3d _
, e as d3d , f as d3d , g as d3d , h as d3d )
declare sub fillgrens( a as double , b as double , c as double , d as double _
, e as double , f as double , g as double , h as double )
end type
constructor tmaterial()
dim as integer i
for i = 0 to 10
diffuse( i ) = d3d( 1 , 1 , 1 )
grens( i + 1 ) = 1
next i
kleurmax = 1
scale = d3d( 1 , 1 , 1 )
mattype = turbulence
turbsize = 0
wavesize = 0
end constructor
sub tmaterial.filldiffuse( a as d3d , b as d3d , c as d3d , d as d3d _
, e as d3d , f as d3d , g as d3d , h as d3d )
diffuse( 0 ) = a
diffuse( 1 ) = b
diffuse( 2 ) = c
diffuse( 3 ) = d
diffuse( 4 ) = e
diffuse( 5 ) = f
diffuse( 6 ) = g
diffuse( 7 ) = h
end sub
sub tmaterial.fillgrens( a as double , b as double , c as double , d as double _
, e as double , f as double , g as double , h as double )
grens( 0 ) = 0
grens( 1 ) = a
grens( 2 ) = b
grens( 3 ) = c
grens( 4 ) = d
grens( 5 ) = e
grens( 6 ) = f
grens( 7 ) = g
grens( 8 ) = h
end sub

dim shared as tmaterial material

dim shared as integer spheretel
const as integer spheremax = 100
type tSphere
dim as d3d center
dim as tmaterial mat
dim as m33 mp
declare sub fill( c as d3d , r as double )
declare function hit( o as d3d , d as d3d ) as double
end type

sub tSphere.fill( c as d3d , r as double )
spot c.x , c.y , c.z
center = c
mat = material
mp = mv( matrix_number ).inverse()
end sub

const as double infinity = 1e9

function tSphere.hit( o as d3d , d as d3d ) as double
dim as double t , a , b , c , disc
dim as d3d temp = o - center
a = dot( d , d )
b = 2 * dot( temp , d )
c = dot( temp , temp ) - radius2
disc = b ^ 2 - 4 * a * c
if disc < 0 then
return infinity
else
dim as double e = sqr( disc )
dim as double demon = 2 * a
t = ( -b - e ) / demon
if t > 1e-12 then
return t
end if

t = ( -b + e ) / demon
if t > 1e-12 then
return t
end if
end if
return infinity
end function

dim shared as tSphere spheres( spheremax )
sub addsphere( x as double , y as double , z as double , d as double , m as tmaterial )
material = m
spheres( spheretel ).fill d3d( x , y , z  ) , d
spheretel += 1
end sub

function pixel( o as d3d , d as d3d ) as d3d
dim as integer i , klno , isph = -1
dim as double dist , slow = infinity
for i = 0 to spheretel
dist = spheres( i ).hit( o , d )
if dist < slow then
slow = dist
isph = i
end if
next i
if isph = -1 then return black
dim as tsphere bol = spheres( isph )
dim as d3d n , p , plek
p = o + d  * slow
n = p - bol.center
plek = n
lokal plek.x , plek.y , plek.z , bol.mp
dim as double dd
if bol.mat.mattype and turbulence then
dd = turbulence3d( _
( plek.x + bol.radius ) / bol.mat.scale.x _
, ( plek.y + bol.radius ) / bol.mat.scale.y _
, ( plek.z + bol.radius ) / bol.mat.scale.z , 64 )
dd = dd * bol.mat.turbsize
end if
if bol.mat.mattype and wave then
dd = sin( plek.y * bol.mat.wavesize + dd ) / 2 + .5
end if
klno = 0
for i = 0 to bol.mat.kleurmax
if ( bol.mat.grens( i ) <= dd _
and bol.mat.grens( i + 1 ) > dd ) then
klno = i
end if
next i
dim as d3d kll = bol.mat.diffuse( klno )
dim as d3d klh = bol.mat.diffuse( klno + 1 )
dim as double gl = bol.mat.grens( klno )
dim as double gh = bol.mat.grens( klno + 1 )
dim as double f = ( dd - gl ) / ( gh - gl )
dim as d3d kl = kll + ( klh - kll ) * f
dim as double a = getangle( n , d3d(-1,1,0) )
return kl * ( cos( a ) / 2 + .5 )
end function
dim shared as tmaterial earth , moon , mars , jupiter , neptune
earth.filldiffuse blue , cyan , yellow , green , white / 2 , white , red , red
earth.fillgrens     .5 , .6   , .7     , .8    , .9        , 1     ,  1  , 1
earth.kleurmax = 6
earth.scale = d3d( 1/3 , 1/3 , 1/3 )
earth.mattype = turbulence
earth.turbsize = 1
earth.wavesize = 0

moon.grens( 0 ) = 0
moon.filldiffuse black , white , red , red ,red , red , red , red
moon.fillgrens 1,1,1,1,1,1,1,1
moon.kleurmax = 1
moon.scale = d3d( .1 , .1 , .1 )
moon.mattype = turbulence
moon.turbsize = 1
moon.wavesize = 0

mars.filldiffuse red , yellow , red , red , red , red , red , red
mars.fillgrens 1,1,1,1,1,1,1,1
mars.kleurmax = 1
mars.scale = d3d( .1 , .1 , .1 )
mars.mattype = turbulence
mars.turbsize = 1
mars.wavesize = 0

jupiter.filldiffuse red / 2 , red , yellow , white / 2 , white , white , red , red
jupiter.fillgrens .5,.6,.7,.9,1,1,1,1
jupiter.kleurmax = 5
jupiter.scale = d3d( 1 , .5 , 1 )
jupiter.mattype = turbulence or wave
jupiter.wavesize = 1 / 10
jupiter.turbsize = 10

neptune.grens( 0 ) = 0
neptune.filldiffuse blue , white , white , white , red , red , red , red
neptune.fillgrens .9 , 1 ,1,1,1,1,1,1
neptune.kleurmax = 1
neptune.scale = d3d( 1 , 1 , 1 )
neptune.mattype = turbulence ''or wave
neptune.wavesize = 0 ''1 / 5
neptune.turbsize = 1 ''0
sub test
spheretel = 0
link 1 , 0,0,0 , 0,0,0 , xyz , 0
addsphere 0,0,0 , 200 , jupiter
dim as integer i
for i = 0 to 10
link 1 , 0,0,0 , i*36,0,0 , xyz , 0
addsphere 350,0,0 , 70 , earth
next i

dim as d3d o , d
dim as double x , y
for x = -winx / 2 to winx / 2
for y = -winy / 2 to winy / 2
o = d3d( 0 , 0 , -1000 )
d = d3d( x , y , 1000 )
pset ( winx / 2 + x , winy / 2 - y ) _
, pixel( o , d ).tocolor()
next y
next x
end sub
cls
paint ( 1 , 1 ) , rgb( 127 , 127 , 127 )
test
''bsave "earth.bmp" , 0
sleep