Hello.
A beautiful program. Merry Christmas for everibody.
// PhotonMapping
// Adapted from FreeBASIC to Yabasic by Galileo 12/2017
// Original code in https://www.freebasic.net/forum/viewtopic.php?f=7&t=15154&sid=29297a206f6e6ae000da4fed47242502&start=15
clear screen
xdimension = 256
ydimension = 256
ntypes = 2 // two types of objects sphere and plane
maxphotons = 200
maxphotonrefections = 1
photonsize = 0.1
ginversesquared = 1.0 / photonsize * photonsize
milliseconds = 500
gexposure = 1.0 / (1000 / milliseconds)
maxraytracerefections = 12
nspheres = 3 // number of spheres
nplanes = 6 // number of planes (a box)
dim gnobjects(1)
gnobjects(0) = nspheres
gnobjects(1) = nplanes
gambient = 0
dim gorigin(2)
dim glight(2)
glight(0) = 0
glight(1) = 1
glight(2) = 4 //x,y,z position
dim gspheres(nspheres - 1, 3)
data -1.0, -1.0, 4.0, 0.5, 6.4, 0.0, 2.8, 5.0, 1.0, -1.0, 4.0, 0.5 // x,y,z, position and radius
for i = 0 to nspheres - 1
for j = 0 to 3
read gspheres(i, j)
next j
next i
dim gplanes(nplanes - 1, 1)
data 0, 1.5, 0, -1.5, 1, -1.5, 1, 1.5, 2, 5.0, 2, -1.0
for i = 0 to nplanes - 1
for j = 0 to 1
read gplanes(i, j)
next j
next i
dim numberofphotons(1, 5)
dim gphotons(ntypes, nplanes + nspheres, maxphotons + maxphotons * maxphotonrefections, 2, 2) //3 * x,y,z
dim gpoint(2)
sub rnd2()
return ran()-ran()
end sub
sub squareddistance(a(), b(), d2)
local ab, d
ab = a(0) - b(0)
d = ab * ab
if (d > d2) return false
ab = a(1) - b(1)
d = d + ab * ab
if (d > d2) return false
ab = a(2) - b(2)
d = d + ab * ab
if (d > d2) return false
gdist2 = d
return true
end sub
sub surfacenormal(r(), typ, idx, hitpoint(), o())
local normal(2), l, axis
if (typ = 0) then
normal(0) = hitpoint(0) - gspheres(idx, 0)
normal(1) = hitpoint(1) - gspheres(idx, 1)
normal(2) = hitpoint(2) - gspheres(idx, 2)
elseif (typ = 1) then
axis = gplanes(idx, 0)
normal(axis) = o(axis) - gplanes(idx, 1)
else
//beep
end if
l = normal(0) * normal(0) + normal(1) * normal(1) + normal(2) * normal(2)
if l l = 1 / sqrt(l)
r(0) = normal(0) * l
r(1) = normal(1) * l
r(2) = normal(2) * l
end sub
sub mirror2(ret(), direction(), hitpoint(), normale())
local l
hitpoint(0) = hitpoint(0) + direction(0) * -0.1
hitpoint(1) = hitpoint(1) + direction(1) * -0.1
hitpoint(2) = hitpoint(2) + direction(2) * -0.1
l = normale(0) * direction(0) + normale(1) * direction(1) + normale(2) * direction(2)
l = l * -2
ret(0) = direction(0) + l * normale(0)
ret(1) = direction(1) + l * normale(1)
ret(2) = direction(2) + l * normale(2)
end sub
sub mirrorvec(ret(), direction(), origen())
local l, normale(2)
surfacenormal(normale(), gtype, gindex, gpoint(), origen())
gpoint(0) = gpoint(0) + direction(0) * -0.1
gpoint(1) = gpoint(1) + direction(1) * -0.1
gpoint(2) = gpoint(2) + direction(2) * -0.1
l = normale(0) * direction(0) + normale(1) * direction(1) + normale(2) * direction(2)
l = l * -2
ret(0) = direction(0) + l * normale(0)
ret(1) = direction(1) + l * normale(1)
ret(2) = direction(2) + l * normale(2)
end sub
//
// raytracing
//
sub raytrace(raydirection(), rayorigin())
local idx, axis, l, spheredirection(2), a, r, b, c, d, sign
gintersect = false
gdist = 999999.9
if gnobjects(1) > 0 then
for idx = 0 to gnobjects(1) - 1
axis = gplanes(idx, 0)
if raydirection(axis) <> 0 then
l = (gplanes(idx, 1) - rayorigin(axis)) / raydirection(axis)
if (l > 0) and (l < gdist) then
gtype = 1
gindex = idx
gdist = l
gintersect = true
end if
end if
next idx
end if
if gnobjects(0) > 0 then
a = raydirection(0) * raydirection(0) + raydirection(1) * raydirection(1) + raydirection(2) * raydirection(2)
a = a + a
for idx = 0 to gnobjects(0) - 1
spheredirection(0) = rayorigin(0) - gspheres(idx, 0)
spheredirection(1) = rayorigin(1) - gspheres(idx, 1)
spheredirection(2) = rayorigin(2) - gspheres(idx, 2)
r = gspheres(idx, 3) * gspheres(idx, 3)
b = spheredirection(0) * raydirection(0) + spheredirection(1) * raydirection(1) + spheredirection(2) * raydirection(2)
b = b + b
c = spheredirection(0) * spheredirection(0) + spheredirection(1) * spheredirection(1) + spheredirection(2) * spheredirection(2)
c = c - r
d = b * b - 2 * a * c
if d > 0.0 then
if c < -0.0001 then
sign = 1
else
sign = -1
end if
l = (-b + sign * sqrt(d)) / a
if (l > 0.0) and (l < gdist) then
gtype = 0
gindex = idx
gdist = l
gintersect = true
end if
end if
next idx
end if
end sub
sub absorbcolor(ret(), rgbin(), r, g, b)
local rgbout(2), c
rgbout(0) = r
rgbout(1) = g
rgbout(2) = b
for c = 0 to 2
if rgbout(c) < rgbin(c) then
ret(c) = rgbout(c)
else
ret(c) = rgbin(c)
end if
next c
end sub
sub getcolor(r(), rgbin(), typ, idx)
if (typ = 0) then // spheres
if idx = 0 then
absorbcolor(r(), rgbin(), 0.0, 0.0, 0.0)
elseif idx = 1 then
absorbcolor(r(), rgbin(), 0.0, 0.0, 0.0)
elseif idx = 2 then
absorbcolor(r(), rgbin(), 0.2, 0.2, 0.8)
end if
elseif (typ = 1) then // planes
if idx = 0 then
absorbcolor(r(), rgbin(), 0.1, 0.8, 0.1) // right
elseif idx = 1 then
absorbcolor(r(), rgbin(), 0.8, 0.1, 0.1) // left
elseif idx = 2 then
absorbcolor(r(), rgbin(), 0.5, 0.5, 0.0) // floor
elseif idx = 3 then
absorbcolor(r(), rgbin(), 0.2, 0.2, 0.2) // ceil
elseif idx = 4 then
absorbcolor(r(), rgbin(), 0.0, 0.0, 0.0) // front
elseif idx = 5 then
absorbcolor(r(), rgbin(), 0.0, 0.0, 0.0) // behind camera
end if
end if
end sub
//
// photon mapping
//
sub gatherphotons(energy(), hitpoint(), typ, idx)
local n(2), v(2), weight, frgb(2), i, cosin
surfacenormal(n(), typ, idx, hitpoint(), gorigin())
for i = 0 to numberofphotons(typ, idx) - 1
// photon location
v(0) = gphotons(typ, idx, i, 0, 0)
v(1) = gphotons(typ, idx, i, 0, 1)
v(2) = gphotons(typ, idx, i, 0, 2)
// in the near of an active photon ?
if (squareddistance(hitpoint(), v(), ginversesquared)) then
// photon direction
cosin = n(0) * gphotons(typ, idx, i, 1, 0) + n(1) * gphotons(typ, idx, i, 1, 1) + n(2) * gphotons(typ, idx, i, 1, 2)
weight = max(0, -cosin)
if weight then
weight = weight * (1 - sqrt(gdist2)) * gexposure
if weight then
// photon energy
frgb(0) = frgb(0) + gphotons(typ, idx, i, 2, 0) * weight
frgb(1) = frgb(1) + gphotons(typ, idx, i, 2, 1) * weight
frgb(2) = frgb(2) + gphotons(typ, idx, i, 2, 2) * weight
end if
end if
end if
next i
for i = 0 to 2
energy(i) = max(0, min(1, frgb(i)))
next i
end sub
sub storephoton(typ, idx, l(), d(), e())
local photon, i
photon = numberofphotons(typ, idx)
for i = 0 to 2
gphotons(typ, idx, photon, 0, i) = l(i) // location
gphotons(typ, idx, photon, 1, i) = d(i) // direction
gphotons(typ, idx, photon, 2, i) = e(i) // energy
next i
numberofphotons(typ, idx) = photon + 1
end sub
sub shadowphoton(raydirection())
local bumpedpoint(2), shadowpoint(2), shadowenerg(2), oldpoint(2), oldtype, oldindex, olddist
shadowenerg(0) = -0.2
shadowenerg(1) = -0.2
shadowenerg(2) = -0.2
oldpoint(0) = gpoint(0)
oldpoint(1) = gpoint(1)
oldpoint(2) = gpoint(2)
oldtype = gtype : oldindex = gindex : olddist = gdist
bumpedpoint(0) = gpoint(0) + raydirection(0) * 0.000001
bumpedpoint(1) = gpoint(1) + raydirection(1) * 0.000001
bumpedpoint(2) = gpoint(2) + raydirection(2) * 0.000001
raytrace(raydirection(), bumpedpoint())
shadowpoint(0) = bumpedpoint(0) + raydirection(0) * gdist
shadowpoint(1) = bumpedpoint(1) + raydirection(1) * gdist
shadowpoint(2) = bumpedpoint(2) + raydirection(2) * gdist
storephoton(gtype, gindex, shadowpoint(), raydirection(), shadowenerg())
gpoint(0) = oldpoint(0)
gpoint(1) = oldpoint(1)
gpoint(2) = oldpoint(2)
gtype = oldtype
gindex = oldindex
gdist = olddist
end sub
sub emitphotons()
local photonenergy(2), photondirection(2), photonorigin(2), l, typ, idx, i, reflection
for typ = 0 to ntypes - 1
for idx = 0 to gnobjects(typ) - 1
numberofphotons(typ, idx) = 0
next idx
next typ
for i = 0 to maxphotons - 1
reflection = 0
// random photon energy
photonenergy(0) = ran()
photonenergy(1) = ran()
photonenergy(2) = ran()
// normalize energy
l = photonenergy(0) * photonenergy(0) + photonenergy(1) * photonenergy(1) + photonenergy(2) * photonenergy(2)
if l then
l = 1.0 / sqrt(l)
photonenergy(0) = photonenergy(0) * l
photonenergy(1) = photonenergy(1) * l
photonenergy(2) = photonenergy(2) * l
end if
// radom photon direction
photondirection(0) = rnd2()
photondirection(1) = rnd2() * 2
photondirection(2) = rnd2()
// normalize direction
l = photondirection(0) * photondirection(0) + photondirection(1) * photondirection(1) + photondirection(2) * photondirection(2)
if l then
l = 1.0 / sqrt(l)
photondirection(0) = photondirection(0) * l
photondirection(1) = photondirection(1) * l
photondirection(2) = photondirection(2) * l
end if
// photon position origin from light
photonorigin(0) = glight(0)
photonorigin(1) = glight(1)
photonorigin(2) = glight(2)
raytrace(photondirection(), photonorigin())
while((gintersect<>0) and (reflection < maxphotonrefections))
reflection = reflection + 1
gpoint(0) = photonorigin(0) + photondirection(0) * gdist
gpoint(1) = photonorigin(1) + photondirection(1) * gdist
gpoint(2) = photonorigin(2) + photondirection(2) * gdist
getcolor(photonenergy(), photonenergy(), gtype, gindex)
l = 1.0 / reflection
photonenergy(0) = photonenergy(0) * l
photonenergy(1) = photonenergy(1) * l
photonenergy(2) = photonenergy(2) * l
storephoton(gtype, gindex, gpoint(), photondirection(), photonenergy())
shadowphoton(photondirection())
mirrorvec(photondirection(), photondirection(), photonorigin())
raytrace(photondirection(), gpoint())
photonorigin(0) = gpoint(0)
photonorigin(1) = gpoint(1)
photonorigin(2) = gpoint(2)
wend
next i
end sub
sub getpixelcolor(pixelrgb(), x, y, z)
local raydirection(2), ndivs, nreflection, mirrorsrgb(2), mirrgb(2)
raydirection(0) = x
raydirection(1) = y
raydirection(2) = z
raytrace(raydirection(), gorigin())
if (gintersect) then
gpoint(0) = gorigin(0) + raydirection(0) * gdist
gpoint(1) = gorigin(1) + raydirection(1) * gdist
gpoint(2) = gorigin(2) + raydirection(2) * gdist
gatherphotons(pixelrgb(), gpoint(), gtype, gindex)
while(((gtype = 0 and gindex < 2) or gindex > 3) and gintersect and (nreflection < maxraytracerefections))
nreflection = nreflection + 1
mirrorvec(raydirection(), raydirection(), gorigin())
raytrace(raydirection(), gpoint())
if (gintersect) then
ndivs = ndivs + 1
gpoint(0) = gpoint(0) + raydirection(0) * gdist
gpoint(1) = gpoint(1) + raydirection(1) * gdist
gpoint(2) = gpoint(2) + raydirection(2) * gdist
gatherphotons(mirrgb(), gpoint(), gtype, gindex)
mirrorsrgb(0) = mirrorsrgb(0) + mirrgb(0)
mirrorsrgb(1) = mirrorsrgb(1) + mirrgb(1)
mirrorsrgb(2) = mirrorsrgb(2) + mirrgb(2)
end if
wend
if ndivs > 0 then
mirrorsrgb(0) = mirrorsrgb(0) / ndivs
mirrorsrgb(1) = mirrorsrgb(1) / ndivs
mirrorsrgb(2) = mirrorsrgb(2) / ndivs
pixelrgb(0) = pixelrgb(0) * 0.25 + mirrorsrgb(0) * 0.75
pixelrgb(1) = pixelrgb(1) * 0.25 + mirrorsrgb(1) * 0.75
pixelrgb(2) = pixelrgb(2) * 0.25 + mirrorsrgb(2) * 0.75
end if
else
pixelrgb(0) = 1 // !!! debug only !!!
pixelrgb(1) = 0
pixelrgb(2) = 1
end if
end sub
sub render()
local h, m, s, l, t, b(2), x, y
//dim as double t1=timer
//windowtitle " photonmapper rendering ..."
for t = 0 to ydimension - 1
y = -(t / ydimension - 0.5)
for l = 0 to xdimension - 1
x = (l / xdimension - 0.5) * ratio
getpixelcolor(b(), x, y, 1)
color b(0) * 255, b(1) * 255, b(2) * 255
dot l, t
next l
next t
end sub
//
// main
//
if xdimension > ydimension then
ratio = xdimension / ydimension
else
ratio = ydimension / xdimension
end if
open window xdimension, ydimension
emitphotons()
render()
print time$()
