Author Topic: PhotoMapping  (Read 1404 times)

Galileo

  • Guest
PhotoMapping
« on: December 21, 2017, 10:44:10 AM »
Hello.

A beautiful program. Merry Christmas for everibody.

Code: [Select]
// 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$()

B+

  • Guest
Re: PhotoMapping
« Reply #1 on: December 22, 2017, 12:06:34 AM »
To All, Merry Christmas!

Galileo, Nice graphic effect. :)