RayTracing versus PhotonMapping

Post your FreeBASIC tips and tricks here. Please don’t post your code without including an explanation.
D.J.Peters
Posts: 8023
Joined: May 28, 2005 3:28
Contact:

RayTracing versus PhotonMapping

Postby D.J.Peters » Feb 21, 2010 3:39

http://web.cs.wpi.edu/~emmanuel/courses ... pping.html

You can switch between RayTracing and PhotonMapping while rendering

[p] = photon mapping
[r] = ray tracing
[v] = view photons

any other key = quit

Joshy

Code: Select all

' RayTracing versus PhotonMapping

type REAL as single

const as integer Dimension     = 256*2
const as integer Types         = 2
const as integer MaxPhotons    = 3000
const as integer MaxRelections = 3
const as REAL gRadius2         = 0.7
const as REAL gExposure        = 40

dim shared as integer Objects(1) = {1,5}
dim shared as REAL gAmbient    = 0.2
dim shared as REAL gOrigin(2)  = {0 ,0 , 0}
dim shared as REAL gLight(2)   = {0,1.2,3}
dim shared as REAL gSphere(3) =  { -0.5, -0.75, 4, 0.75}
dim shared as REAL gPlanes(4,1) = {{0,  1.5}, _
                                    {1, -1.5}, _
                                    {0, -1.5}, _
                                    {1,  1.5}, _
                                    {2,  5.0}}
dim shared as boolean gPhotonMapping = true
dim shared as integer NumberOfPhotons(1,4)
dim shared as REAL gPhotons(1,4,MaxPhotons*5,2,2)
dim shared as boolean gIntersect
dim shared as integer gType
dim shared as integer gIndex
dim shared as REAL gDist2, gDist
dim shared as REAL gPoint(2)

dim shared as boolean gEmpty = true
dim shared as boolean gView3D
dim shared as integer pRow, pCol, pIteration, pMax

declare sub GatherPhotons(r() as REAL, p() as REAL,typ as integer, index as integer)
declare sub GetColor(r() as REAL, rgbIn() as REAL, typ as integer, index as integer)
declare sub StorePhoton(typ as integer,index as integer, _
                        location() as REAL,direction() as REAL,energy() as REAL)
declare sub ShadowPhoton(ray() as REAL)
declare sub DrawPhoton(frgb() as REAL,p() as REAL)
declare sub ResetRender
declare sub Render

#define min(a,b) iif(a<b,a,b)
#define max(a,b) iif(a>b,a,b)

#macro NormalizeVec(r,v)
  scope
  dim as REAL L2 = v(0)*v(0)+v(1)*v(1)+v(2)*v(2)
  if l2<>0 then L2=1.0/sqr(L2)
  r(0)=v(0)*L2
  r(1)=v(1)*L2
  r(2)=v(2)*L2
  end scope
#endmacro

#macro SubVec(r,a,b)
  r(0)=a(0)-b(0)
  r(1)=a(1)-b(1)
  r(2)=a(2)-b(2)
#endmacro

#macro AddVec(r,a,b)
  r(0)=a(0)+b(0)
  r(1)=a(1)+b(1)
  r(2)=a(2)+b(2)
#endmacro

#macro MulScalarVec(r,a,b)
  r(0)=a(0)*b
  r(1)=a(1)*b
  r(2)=a(2)*b
#endmacro

#define DotProduct(a,b) (a(0)*b(0) + a(1)*b(1) + a(2)*b(2))

#define rnd2 (rnd-rnd)

#macro RandomVec(r)
  r(0)=rnd-rnd
  r(1)=rnd-rnd
  r(2)=rnd-rnd
#endmacro

#define odd(xx) (xx and 1)


function Distance2(a() as REAL, _
                   b() as REAL, _
                   d2 as REAL) as boolean
  dim as REAL c = a(0) - b(0)
  dim as REAL d = c*c
  if (d > d2) then return false
  c = a(1) - b(1)
  d += c*c
  if (d > d2) then return false
  c = a(2) - b(2)
  d += c*c
  if (d > d2) then return false
  gDist2 = d2
  return true
end function

sub RaySphere(idx as integer,r() as REAL, o() as REAL)
  dim as REAL s(2)
  SubVec(s,gSphere,o)
  dim as REAL radius = gSphere(3)
  dim as REAL A = DotProduct(r,r)
  dim as REAL B = -2.0 * DotProduct(s,r)
  dim as REAL C = DotProduct(s,s) - (radius*radius)
  dim as REAL D = B*B - 4*A*C
  if (D < 0.0) then return
  dim as REAL sign  = iif(C < -0.00001,1,-1)
  dim as REAL l = (-B + sign*sqr(D))/(2*A)
  if (l<0) or (l>gDist) then return
  gType  = 0
  gIndex = idx
  gDist  = l
  gIntersect = true
end sub

sub RayPlane(idx as integer, r() as REAL,o() as REAL)
  dim as integer axis = gPlanes(idx,0)
  if r(axis)=0 then return
  dim as REAL l = (gPlanes(idx,1) - o(axis)) / r(axis)
  if (l<0) or (l>gDist) then return
  gType  = 1
  gIndex = idx
  gDist  = l
  gIntersect = true 
end sub

sub RayObject(typ as integer, idx as integer, r()as REAL, o() as REAL)
  if (typ = 0) then
    RaySphere(idx,r(),o())
  else
    RayPlane(idx,r(),o())
  end if
end sub

sub SphereNormal(r() as REAL,P() as REAL)
  dim as REAL v(2)
  SubVec(v,P,gSphere)
  NormalizeVec(r,v)
end sub

sub PlaneNormal(r() as REAL, idx as integer, P() as REAL, O() as REAL)
  dim as integer axis = gPlanes(idx,0)
  dim as REAL N(2)
  N(axis) = O(axis) - gPlanes(idx,1)
  NormalizeVec(r,N)
end sub

sub SurfaceNormal(r() as REAL, _
                  typ as integer, _
                  idx as integer, _
                  P() as REAL, _
                  Inside() as REAL)
  if (typ = 0) then
    SphereNormal(r(),P())
  else
    PlaneNormal(r(),idx,P(),Inside())
  end if
end sub

sub MirrorVec(Ret() as REAL,_
              Ray() as REAL, _
              FromPoint() as REAL)
  dim as REAL N(2)=any,tmp(2)=any
  SurfaceNormal(N(),gType, gIndex, gPoint(), fromPoint())
  MulScalarVec(tmp,N,(2 * DotProduct(ray,N) ))
  SubVec(tmp,ray,tmp)
  NormalizeVec(Ret,tmp)
end sub

'
' Lighting
'
function LightDiffuse(N() as REAL,P() as REAL) as REAL
  dim as REAL L(2)=any
  SubVec(L,gLight,P)
  NormalizeVec(L,L)
  return DotProduct(N,L)
end function

function LightObject(typ as integer, _
                     idx as integer, _
                     P() as REAL, _
                     Ambient as REAL ) as REAL
  dim as REAL N(2)=any
  SurfaceNormal(N(),typ, idx, P(), gLight())
  dim as REAL L = LightDiffuse(N() , P() )
  return min(1.0, max(L, Ambient))
end function

'
' Raytracing
'
sub Raytrace(ray() as REAL,origin() as REAL)
  gIntersect = false
  gDist      = 999999.9
  for typ as integer = 0 to Types-1
    for idx as integer = 0 to Objects(typ)-1
      RayObject(typ,idx,ray(),origin())
    next
  next
end sub

sub AbsorbColor(ret() as REAL, _
                rgbIn() as REAL, _
                r as REAL,g as REAL,b as REAL) ' e.g. White Light Hits Red Wall
  dim as REAL rgbOut(2)={r,g,b}
  for c as integer =0 to 2
    ret(c) = min(rgbOut(c),rgbIn(c))
  next
end sub

sub GetColor(r() as REAL, _
            rgbIn() as REAL, _
            typ as integer, _
            idx as integer)
  if (typ=0) then     ' sphere
    AbsorbColor(r(),rgbIn(), 1, 1, 0.5)
  elseif (typ=1) then ' plane
    if idx=0 then
      AbsorbColor(r(),rgbIn(), 1.0, 0, 0)
    elseif idx=2 then
      AbsorbColor(r(),rgbIn(), 0, 1.0, 0)
    else 
      AbsorbColor(r(),rgbIn(), 1, 1, 1)
    end if
  end if
end sub

sub ComputePixelColor(prgb() as REAL,x as REAL,y as REAL)
  dim as REAL ray(2) = {  x/Dimension - 0.5 , _
                          -(y/Dimension - 0.5), _
                                          1.0}
  Raytrace(ray(), gOrigin())
  if (gIntersect)  then
    MulScalarVec(gPoint,ray,gDist)
    if (gType = 0) then
      MirrorVec(ray(),ray(),gOrigin())
      Raytrace(ray(), gPoint())
      if (gIntersect) then
        dim as REAL tmp(2)=any
        MulScalarVec(tmp,ray,gDist)
        AddVec(gPoint,tmp,gPoint)
      end if
    end if
 
    if (gPhotonMapping) then
      GatherPhotons(prgb(),gPoint(),gType,gIndex)
    else
      dim as integer tType  = gType
      dim as integer tIndex = gIndex
      dim as REAL a = gAmbient
      dim as REAL tmp(2)=any
      SubVec(tmp,gPoint,gLight)
      Raytrace(tmp(),gLight())
      if (tType = gType) and (tIndex = gIndex) then
        a = LightObject(gType, gIndex, gPoint(), gAmbient)
      end if
      prgb(0)=a:prgb(1)=a:prgb(2)=a
      GetColor(prgb(),prgb(),tType,tIndex)
    end if
  end if
end sub

'
' Photon Mapping
'
sub GatherPhotons(energy() as REAL, _
                  p() as REAL, _
                  typ as integer, _
                  idx as integer)
  dim as REAL N(2)=any
  dim as REAL tmp(2)=any
  dim as REAL g(2)=any
  dim as REAL weight=any
  dim as REAL frgb(2)=any
  SurfaceNormal(N(), typ, idx, p(), gOrigin())

  for i as integer = 0  to  NumberOfPhotons(typ,idx)-1
    ' location
    g(0)=gPhotons(typ,idx,i,0,0)
    g(1)=gPhotons(typ,idx,i,0,1)
    g(2)=gPhotons(typ,idx,i,0,2)
    if (Distance2(p(),g(),gRadius2)) then
      ' direction
      g(0)=gPhotons(typ,idx,i,1,0)
      g(1)=gPhotons(typ,idx,i,1,1)
      g(2)=gPhotons(typ,idx,i,1,2)
      weight = max(0.0, -DotProduct(N,g) )
      weight *= (1.0 - sqr(gDist2)) / gExposure
      ' energy
      g(0)=gPhotons(typ,idx,i,2,0)
      g(1)=gPhotons(typ,idx,i,2,1)
      g(2)=gPhotons(typ,idx,i,2,2)
      MulScalarVec(tmp,g,weight)
      AddVec(frgb,frgb,tmp)
    end if
  next
  for j as integer=0 to 2
    energy(j)=max(0,min(1,frgb(j) ) )
  next
end sub

sub EmitPhotons
  randomize 1
  dim as REAL frgb(2)=any
  dim as REAL ray(2)=any
  dim as REAL p(2)=any
  for typ as integer = 0 to Types-1
    for idx as integer = 0 to Objects(typ)-1
      NumberOfPhotons(typ,idx) = 0
    next
  next
  for i as integer = 0 to MaxPhotons-1
    dim as integer bounces = 1
    ' white photon color
    frgb(0)=1:frgb(1)=1:frgb(2)=1
    RandomVec(ray)
    NormalizeVec(ray,ray)
    p(0) = gLight(0)
    p(1) = gLight(1)
    p(2) = gLight(2)
    while (p(1) >= gLight(1))
      dim as REAL N(2)=any
      RandomVec(N)
      NormalizeVec(N,N)
      MulScalarVec(N,N,0.75)
      AddVec(p,gLight,N)
    wend
   
    Raytrace(ray(), p())
       
    if abs(p(0) > 1.5) then Continue for
    if abs(p(1) > 1.5) then  Continue for
    if Distance2(p(), gSphere(),gSphere(3)*gSphere(3)) then  Continue for
    while (gIntersect<>0) and (bounces <= MaxRelections)
      dim as REAL tmp(2)=any
      MulScalarVec(tmp,ray,gDist)
      AddVec(gPoint, tmp, p)
      GetColor(frgb(),frgb(),gType,gIndex)
      MulScalarVec(frgb,frgb, 1.0/sqr(bounces))
      StorePhoton(gType, gIndex, gPoint(), ray(),frgb())
      DrawPhoton(frgb(), gPoint())
      ShadowPhoton(ray())
      MirrorVec(ray(),ray(),p())
      Raytrace(ray(), gPoint())
      p(0) = gPoint(0)
      p(1) = gPoint(1)
      p(2) = gPoint(2)
      bounces+=1
    wend
  next
end sub

sub StorePhoton(typ as integer, _
                idx as integer, _
                l() as REAL,_
                d() as REAL, _
                e() as REAL)
  dim as integer Photon=NumberOfPhotons(typ,idx)
  for i as integer=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
  NumberOfPhotons(typ,idx)=Photon+1
end sub

sub ShadowPhoton(ray() as REAL)
  dim as REAL shadow(2) = {-0.25,-0.25,-0.25}
  dim as REAL tPoint(2) = {gPoint(0), gPoint(1),gPoint(2)}
  dim as integer tType  = gType ' Save State
  dim as integer tIndex = gIndex
  dim as REAL BumpedPoint(2)
  MulScalarVec(BumpedPoint,ray,0.000001)
  AddVec(BumpedPoint,gPoint,BumpedPoint)
  Raytrace(ray(), BumpedPoint())
  dim as REAL ShadowPoint(2)=any
  MulScalarVec(ShadowPoint,ray,gDist)
  AddVec(ShadowPoint,ShadowPoint, BumpedPoint) ' 3D Point
  StorePhoton(gType, gIndex, ShadowPoint(), ray(), shadow())
  gPoint(0) = tPoint(0)
  gPoint(1) = tPoint(1)
  gPoint(2) = tPoint(2)
  gType  = tType
  gIndex = tIndex
end sub




sub Render
  dim as integer x,y,iterations = 0
  dim as integer nruns
  nruns=256
  while (iterations < nruns)
    nruns=max(pMax, 256)
    if (pCol >= pMax) then
      pRow+=1
      pCol = 0
      if (pRow >= pMax) then
        pIteration+=1
        pRow = 0
        pMax = int(2^pIteration)
      end if
    end if
    dim as integer pNeedsDrawing = (pIteration = 1) or odd(pRow) _
                                 or ( (odd(pRow)=0) and odd(pCol))
    x = pCol * (Dimension/pMax)
    y = pRow * (Dimension/pMax)
   
    pCol+=1
   
    if (pNeedsDrawing) then
      iterations+=1
      dim as REAL b(2)
      ComputePixelColor(b(),x,y)
      dim as uinteger col=rgb(b(0)*255,b(1)*255,b(2)*255)
      line (x,y)-step((Dimension/pMax)-1,(Dimension/pMax)-1),col,BF' Draw the Possibly Enlarged Pixel
      'pset (x,y),col
    end if
  wend
  if (pRow = Dimension-1) then
    gEmpty = false
  end if
end sub

sub ResetRender()
  pRow=0
  pCol=0
  pIteration=1
  pMax=2
  gEmpty=true
  if (gPhotonMapping and not gView3D) then
    EmitPhotons()
  end if
end sub


sub setup()
  screenres Dimension,Dimension,24
  ResetRender
end sub

sub drawit()
  if gPhotonMapping then
    windowtitle "[r]aytracing [v]iew"
  else
    windowtitle "[p]hoton Mapping [v]iew"
  end if
  if (gView3D) then
    if (gEmpty) then
      cls
      emitPhotons()
      gEmpty = false
    end if
  else
    if (gEmpty) then
      render()
    end if
  end if
end sub

sub DrawPhoton(frgb() as REAL,p() as REAL)
  if (gView3D=true) and cbool(p(2) > 0.0) then
    dim as integer x = (Dimension/2) + int(Dimension *  p(0)/p(2))
    dim as integer y = (Dimension/2) + int(Dimension * -p(1)/p(2))
    if (y <= Dimension) then
      pset (x,y),rgb(frgb(0)*255,frgb(1)*255,frgb(2)*255)
    end if
  end if
end sub

sub SwitchToMode(i as string)
  if (i="r") then
    gView3D = false
    gPhotonMapping = false
    resetRender()
  elseif (i="p") then
    gView3D = false
    gPhotonMapping = true 
    resetRender()
  elseif (i="v") then
    gView3D = true
    resetRender()
  end if
end sub

'
' main
'
dim as string Key
Setup()
while key=""
  key=lcase(inkey)
  if key="r" or key="p" or key="v" then
    SwitchToMode(key)
    key=""
  end if 
  DrawIt()
  sleep(100)
wend
Last edited by D.J.Peters on Dec 10, 2015 18:00, edited 2 times in total.
phycowelder
Posts: 74
Joined: Dec 19, 2007 6:55

Postby phycowelder » Feb 21, 2010 4:52

wow! that's awesome! keep up the good work DJ
D.J.Peters
Posts: 8023
Joined: May 28, 2005 3:28
Contact:

Postby D.J.Peters » Feb 21, 2010 10:38

Code: Select all

' PhotonMapping

#define true -1
#define false 0

const as integer Dimension       = 512
const as integer Types           = 2
const as integer MaxPhotons      = 6000
const as integer MaxRelections   = 3
const as single gRadius2         = 0.5
const as single gExposure        = 30

dim shared as integer Objects(1) = {1,5}
dim shared as single gAmbient    = 0.2
dim shared as single gOrigin(2)  = {0,0,0} 'x,y,z
dim shared as single gLight(2)   = {0,0,5} 'x,y,z                                 
dim shared as single gSphere(3) =  {-0.75, -0.75, 5, 0.75} ' x,y,z,radius
dim shared as single gPlanes(4,1) = {{0,  1.5}, _
                                     {1, -1.5}, _
                                     {0, -1.5}, _
                                     {1,  1.5}, _
                                     {2,  7.0}}
dim shared as integer PhotonMapping  = true
dim shared as integer NumberOfPhotons(1,4)
dim shared as single gPhotons(1,4,MaxPhotons*5,2,2)
dim shared as integer gIntersect
dim shared as integer gType
dim shared as integer gIndex
dim shared as single gDist2,gDist,gDistSqr
dim shared as single gPoint(2)


function min(a as single,b as single) as single
  if b<a then return b
  return a
end function

function max(a as single,b as single) as single
  if b>a then return b
  return a
end function

function rnd2 as single
  return rnd*2-rnd*2
end function 

function Distance2(a() as single, _
                   b() as single, _
                   d2 as single) as integer
  dim as single c = a(0) - b(0)         
  dim as single d = c*c
  if (d > d2) then return false
  c = a(1) - b(1)
  d += c*c
  if (d > d2) then return false
  c = a(2) - b(2)
  d += c*c
  if (d > d2) then return false
  gDist2 = d2
  gDistSqr=sqr(d2)
  return true
end function

sub SurfaceNormal(r() as single, _
                  typ as integer, _
                  idx as integer, _
                  P() as single, _
                  O() as single)
  if (typ = 0) then
    r(0)=p(0)-gSphere(0) 
    r(1)=p(1)-gSphere(1) 
    r(2)=p(2)-gSphere(2) 
    dim as single l=r(0)*r(0)+r(1)*r(1)+r(2)*r(2)
    if l<>0 then l=1/sqr(l)
    r(0)*=l:r(1)*=l:r(2)*=l
  else 
    dim as integer axis = gPlanes(idx,0)
    dim as single N(2)
    N(axis) = O(axis) - gPlanes(idx,1)
    dim as single l=n(0)*n(0)+n(1)*n(1)+n(2)*n(2)
    if l<>0.0 then l=1/sqr(l)
    r(0)=n(0)*l:r(1)=n(1)*l:r(2)=n(2)*l
  end if
end sub

sub MirrorVec(Ret() as single,_
              Ray() as single, _
              FromPoint() as single)
  dim as single N(2)
  SurfaceNormal(N(0),gType, gIndex, gPoint(0), fromPoint(0))
  dim as single l=ray(0)*N(0)+ray(1)*N(1)+ray(2)*N(2)
  l+=l
  ret(0)=ray(0)-n(0)*l
  ret(1)=ray(1)-n(1)*l
  ret(2)=ray(2)-n(2)*l
  l=ret(0)*ret(0)+ret(1)*ret(1)+ret(2)*ret(2)
  if l<>0.0 then l=1/sqr(l)
  Ret(0)*=l:Ret(1)*=l:Ret(2)*=l
end sub

'
' Raytracing
'
sub Raytrace(r() as single,o() as single)
  gIntersect = false
  gDist      = 999999.9
       
  for idx as integer =0 to 4
    dim as integer axis = gPlanes(idx,0)
    if r(axis)<>0 then
      dim as single l = (gPlanes(idx,1) - o(axis)) / r(axis)
      if (l>0) and (l<gDist) then
        gType  = 1
        gIndex = idx
        gDist  = l
        gIntersect = true
      end if
    end if     
  next
   
  dim as single s(2)
  s(0)=gSphere(0)-o(0)
  s(1)=gSphere(1)-o(1)
  s(2)=gSphere(2)-o(2)
  dim as single radius = gSphere(3)
  dim as single A = r(0)*r(0)+r(1)*r(1)+r(2)*r(2)
  dim as single B = -2.0 * (s(0)*r(0)+s(1)*r(1)+s(2)*r(2))
  dim as single C = s(0)*s(0)+s(1)*s(1)+s(2)*s(2) -(radius*radius)
  dim as single D = B*B - 4*A*C
  if (D>0.0) then
    dim as single sign  = iif(C < -0.00001,1,-1)
    dim as single l = (-B + sign*sqr(D))/(2*A)
    if (l>0) and (l<gDist) then
      gType  = 0
      gIndex = 0
      gDist  = l
      gIntersect = true
    end if 
  end if   
end sub

sub AbsorbColor(ret() as single, _
                rgbIn() as single, _
                r as single,g as single,b as single)
  dim as single rgbOut(2)={r,g,b}
  for c as integer =0 to 2
    if rgbOut(c)<rgbIn(c) then
      ret(c)=rgbOut(c)
    else
      ret(c)=rgbIn(c)
    end if
  next
end sub

sub GetColor(r() as single, _
            rgbIn() as single, _
            typ as integer, _
            idx as integer)
  if (typ=0) then     ' sphere
    AbsorbColor(r(0),rgbIn(0), 1, 1, 1)
  elseif (typ=1) then ' plane
    if idx=0 then
      AbsorbColor(r(0),rgbIn(0), 1, 0, 0)
    elseif idx=1 then
      AbsorbColor(r(0),rgbIn(0), 1, 1, 1)
    elseif idx=2 then
      AbsorbColor(r(0),rgbIn(0), 0, 1, 0)
    elseif idx=3 then
      AbsorbColor(r(0),rgbIn(0), 0, 0, 1)   
    else 
      AbsorbColor(r(0),rgbIn(0), 0, 1, 1)
    end if
  end if
end sub

'
' Photon Mapping
'
sub GatherPhotons(energy() as single, _
                  p() as single, _
                  typ as integer, _
                  idx as integer)
  dim as single N(2)
  'dim as single tmp(2)
  dim as single g(2)
  dim as single weight
  dim as single frgb(2)
  SurfaceNormal(N(0), typ, idx, p(0), gOrigin(0))
  for i as integer = 0  to  NumberOfPhotons(typ,idx)-1
    ' location
    g(0)=gPhotons(typ,idx,i,0,0)
    g(1)=gPhotons(typ,idx,i,0,1)
    g(2)=gPhotons(typ,idx,i,0,2)
    if (Distance2(p(0),g(0),gRadius2)) then
      ' direction
      g(0)=gPhotons(typ,idx,i,1,0)
      g(1)=gPhotons(typ,idx,i,1,1)
      g(2)=gPhotons(typ,idx,i,1,2)
      weight = max(0.0, -(N(0)*g(0)+N(1)*g(1)+N(2)*g(2)) )
      weight *= (1.0 - gDistSqr) / gExposure
      ' energy
      frgb(0)+=gPhotons(typ,idx,i,2,0)*weight
      frgb(1)+=gPhotons(typ,idx,i,2,1)*weight
      frgb(2)+=gPhotons(typ,idx,i,2,2)*weight
    end if
  next
  for j as integer=0 to 2
    energy(j)=max(0,min(1,frgb(j) ) )
  next
end sub

sub StorePhoton(typ as integer, _
                idx as integer, _
                l() as single,_
                d() as single, _
                e() as single)
  dim as integer Photon=NumberOfPhotons(typ,idx)
  for i as integer=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
  NumberOfPhotons(typ,idx)=Photon+1
end sub

sub ShadowPhoton(ray() as single)
  dim as single shadow(2) = {-0.25,-0.25,-0.25}
  dim as single tPoint(2) = {gPoint(0), gPoint(1),gPoint(2)}
  dim as integer tType  = gType ' Save State
  dim as integer tIndex = gIndex
  dim as single BumpedPoint(2)
  BumpedPoint(0)=gpoint(0)+ray(0)*0.000001
  BumpedPoint(1)=gpoint(1)+ray(1)*0.000001
  BumpedPoint(2)=gpoint(2)+ray(2)*0.000001
  Raytrace(ray(0), BumpedPoint(0))
  dim as single ShadowPoint(2)
  ShadowPoint(0)=ray(0)*gDist+BumpedPoint(0)
  ShadowPoint(1)=ray(1)*gDist+BumpedPoint(1)
  ShadowPoint(2)=ray(2)*gDist+BumpedPoint(2)
  StorePhoton(gType, gIndex, ShadowPoint(0), ray(0), shadow(0))
  gPoint(0) = tPoint(0)
  gPoint(1) = tPoint(1)
  gPoint(2) = tPoint(2)
  gType  = tType
  gIndex = tIndex
end sub

sub EmitPhotons
  randomize 1
  dim as single frgb(2)
  dim as single ray(2)
  dim as single p(2),l
  for typ as integer = 0 to Types-1
    for idx as integer = 0 to Objects(typ)-1
      NumberOfPhotons(typ,idx) = 0
    next
  next
  for i as integer = 0 to MaxPhotons-1
    dim as integer bounces = 1
    ' white photon color
    frgb(0)=1:frgb(1)=1:frgb(2)=1
    ray(0)=rnd2:ray(1)=rnd2:ray(2)=rnd2:
    l=ray(0)*ray(0)+ray(1)*ray(1)+ray(2)*ray(2)
    if l<>0 then l=1/sqr(l)
    ray(0)*=l:ray(1)*=l:ray(2)*=l
    p(0)=gLight(0):p(1)=gLight(1):p(2)=gLight(2)
    while (p(1) >= gLight(1))
      dim as single N(2)
      N(0)=rnd2:N(1)=rnd2:N(2)=rnd2
      l=N(0)*N(0)+N(1)*N(1)+N(2)*N(2)
      if l<>0 then l=1/sqr(l)
      n(0)*=l:n(1)*=l:n(2)*=l
      p(0)=gLight(0)+N(0)*0.75
      p(1)=gLight(1)+N(1)*0.75
      p(2)=gLight(2)+N(2)*0.75
    wend
    Raytrace(ray(0), p(0))
    if abs(p(0) > 1.5) then Continue for
    if abs(p(1) > 1.5) then Continue for
    if abs(p(2) >   5) then Continue for
    if Distance2(p(0), gSphere(0),gSphere(3)*gSphere(3)) then  Continue for
           
    while (gIntersect<>0) and (bounces <= MaxRelections)
      dim as single tmp(2)
      gpoint(0)=ray(0)*gDist+p(0)
      gpoint(1)=ray(1)*gDist+p(1)
      gpoint(2)=ray(2)*gDist+p(2)
      GetColor(frgb(0),frgb(0),gType,gIndex)
      dim as single sq=1.0/sqr(bounces)
      frgb(0)*=sq
      frgb(1)*=sq
      frgb(2)*=sq
      StorePhoton(gType, gIndex, gPoint(0), ray(0),frgb(0))
      ShadowPhoton(ray(0))
      MirrorVec(ray(0),ray(0),p(0))
      Raytrace(ray(0), gPoint(0))
      p(0)=gPoint(0):p(1)=gPoint(1):p(2)=gPoint(2)
      bounces+=1
    wend
  next
end sub

sub ComputePixelColor(prgb() as single,x as single,y as single)
  dim as single ray(2) = {  x/Dimension - 0.5 , _
                          -(y/Dimension - 0.5), _
                                          1.0}
  Raytrace(ray(0), gOrigin(0))
  if (gIntersect)  then
    gPoint(0)=ray(0)*gDist
    gPoint(1)=ray(1)*gDist
    gPoint(2)=ray(2)*gDist
    if gType=0 then
      MirrorVec(ray(0),ray(0),gOrigin(0))
      Raytrace(ray(0), gPoint(0))
      if (gIntersect) then
        gPoint(0)=ray(0)*gDist+gPoint(0)
        gPoint(1)=ray(1)*gDist+gPoint(1)
        gPoint(2)=ray(2)*gDist+gPoint(2)
      end if
    end if
    GatherPhotons(prgb(0),gPoint(0),gType,gIndex)
  end if
end sub

Sub Render
  Dim As Integer x,y,s,s2
  Dim As Uinteger col
  Dim As Single b(2)
  Dim As Single t1=Timer
  s=16
  While s>0
    WindowTitle " PhotonMapper " & s
    s2=s-1
    For y=0 To Dimension-1 Step s
      Line (0,y+s)-Step(Dimension,s2),0,BF
      screenlock
      For x=0 To Dimension-1 Step s
        ComputePixelColor(b(0),x,y)
        Line (x,y)-Step(s2,s2),rgb(b(0)*255,b(1)*255,b(2)*255),bf
      Next
      screenunlock
      Sleep(50)
      If Inkey=Chr(27) Then Return
    Next
    s=s\2
  Wend
  WindowTitle " PhotonMapper ready "   
End Sub

'
' main
'
ScreenRes Dimension,Dimension,24
WindowTitle " PhotonMapper ... "   
EmitPhotons
Render
bsave "PhotonMapper2.bmp",0
sleep


feel the dirty energy :-)

Joshy
Image
Image
Image
one wrong number
but isn't an error it's computer art :lol:
Image
Last edited by D.J.Peters on Sep 25, 2017 16:40, edited 2 times in total.
rdc
Posts: 1713
Joined: May 27, 2005 17:22
Location: Texas, USA
Contact:

Postby rdc » Feb 21, 2010 14:07

Extremely awesome D.J. dirty or not. :)
cha0s
Site Admin
Posts: 5317
Joined: May 27, 2005 6:42
Location: Illinois
Contact:

Postby cha0s » Feb 21, 2010 15:39

Really sexy.
D.J.Peters
Posts: 8023
Joined: May 28, 2005 3:28
Contact:

Postby D.J.Peters » Feb 21, 2010 16:22

the -gen gas build needs 129 seconds and the -gen gcc build 128 :lol:

edit:
sorry my fault -O 3 was missing and i removed sleep and inkey from the render loop

now i get gas 115 and gcc 43 seconds :-)

Joshy
dodicat
Posts: 6357
Joined: Jan 10, 2006 20:30
Location: Scotland

Postby dodicat » Feb 21, 2010 17:11

Similar idea, different method, simpler example.
I don't know whether it's photon mapping or ray tracing or whatever, in fact I don't think it's any of those, so I'll just call it little squares.

Code: Select all



sub drawpolygon(x() as double,y() as double,colour as uinteger,img as any pointer=0)
    dim k as integer=ubound(x)+1
    dim as integer index,nextindex
    dim as double xc,yc
    for n as integer=1 to ubound(x)'+1
        xc=xc+x(n):yc=yc+y(n)
        index=n mod k:nextindex=(n+1) mod k
        if nextindex=0 then nextindex=1
    line img,(x(index),y(index))-(x(nextindex),y(nextindex)),colour
    next
  xc=xc/ubound(x):yc=yc/ubound(y)
  paint img,(xc,yc),colour,colour
end sub

function inpolygon(colour as uinteger,x as double,y as double) as integer
    if point(x,y)=colour then
        return -1
    else
        return 0
        endif
    end function
    ' ******* ROTATOR ***************
    Sub mv(m1() As Double,m2() As Double,ans() As Double) 'MATRIX x VECTOR
   Dim s As Double
    For i As Integer=1 To 3
        s=0
        For k As Integer = 1 To 3
            s=s+m1(i,k)*m2(k)
        Next k
        ans(i)=s
        Next i
    End Sub
    Dim Shared np(1 To 6) As Double 'For a position rotation, no lines or circles
'np(1),np(2),np(3)-------np(4),np(5),np(6)  for line ends
'np(4),np(5)                    for circle centres
' use linepointset,circlepointset to draw the pixels
        Sub rotate(Byval pivot_x As Double,_  'x pivot for rotation
                   Byval pivot_y As Double,_  'y pivot for rotation
                   Byval pivot_z As Double,_  'z pivot for rotation
                   Byval first_x As Double,_  'x for line,or centre for circle
                   Byval first_y As Double,_  'y for line,or centre for circle
                   Byval first_z As Double,_  'z for line or circle
                   Byval second_x As Double, _'x for line,or radius for circle
                   Byval second_y As Double, _'y for line,or aspect for circle
                   Byval second_z As Double,_ 'z for line, first arc position circle
                   Byval second_arc As Double,_ 'second arc position circle,0 line
                   Byval angleX As Double, _   'angle to rotate round x axis
                   Byval angleY As Double,_    'angle to rotate round y axis
                   Byval angleZ As Double,_    'angle to rotate round z axis
                   Byval magnifier As Double,_ '1=no magnifacation
                   Byval dilator As Double,_   'times distance from pivot(1=no dilation)
                   Byval colour As Integer,_   'color for line or circle
                   Byval thickness As Double,_ 'thickness line or circle
                   Byref shape As String,_ 'line/circle/circlefill/box/boxfill/linepoint[set],circlepoint[set]"
                   Byref mode As String,_    '2d or 3d
                   Byval perspective As Double=0,_ 'add some 3d perspective 0 to 1 approx
                   image As Any Pointer=0)        'write to an image if required
  shape=Lcase$(shape)
  mode=Lcase$(mode)
  Dim th As Double
  th=thickness
  Dim As Double zval,pp   'used in get_perspective
  Dim sx As Double=second_x
Dim p As Double = 4*Atn(1)  '(pi)
Dim angleX_degrees As Double
Dim angleY_degrees As Double
Dim angleZ_degrees As Double

#Macro thickline(t)
Dim As Double s,h,c
Dim As Uinteger prime=rgb(255,255,255)
h=Sqr(((np(1))-(np(4)))^2+((np(2))-(np(5)))^2)
s=((np(5))-np(2))/h
c=(np(1)-(np(4)))/h
Line image, (np(4)+s*t/2,np(5)+c*t/2)-(np(1)+s*t/2,np(2)+c*t/2),prime
Line image, (np(4)-s*t/2,np(5)-c*t/2)-(np(1)-s*t/2,np(2)-c*t/2),prime
Line image, (np(4)+s*t/2,np(5)+c*t/2)-(np(4)-s*t/2,np(5)-c*t/2),prime
Line image, (np(1)+s*t/2,np(2)+c*t/2)-(np(1)-s*t/2,np(2)-c*t/2),prime
Paint image,((np(4)+np(1))/2, (np(5)+np(2))/2),prime,prime

Line image, (np(4)+s*t/2,np(5)+c*t/2)-(np(1)+s*t/2,np(2)+c*t/2),colour
Line image, (np(4)-s*t/2,np(5)-c*t/2)-(np(1)-s*t/2,np(2)-c*t/2),colour
Line image, (np(4)+s*t/2,np(5)+c*t/2)-(np(4)-s*t/2,np(5)-c*t/2),colour
Line image, (np(1)+s*t/2,np(2)+c*t/2)-(np(1)-s*t/2,np(2)-c*t/2),colour
Paint image,((np(4)+np(1))/2, (np(5)+np(2))/2), colour, colour
#EndMacro

#macro thickcircle(t)
Dim As Uinteger prime=rgb(255,255,255)
Dim As Double xp1,xp2,yp1,yp2
Dim arc1 As Double=second_z*p/180
Dim arc2 As Double=second_arc*p/180
arc1=2*p+(arc1-(anglez_degrees))
arc2=2*p+(arc2-(anglez_degrees))
sx=sx*magnifier
If arc1=arc2 Then
     Circle image,(np(4),np(5)),sx,prime,,,second_y
    Circle image,(np(4),np(5)),sx-t,prime,,,second_y
    Paint image,(np(4),np(5)+sx-t/2),prime,prime
    Paint image,(np(4)+sx-t/2,np(5)),prime,prime
    Circle image,(np(4),np(5)),sx,colour,,,second_y
    Circle image,(np(4),np(5)),sx-t,colour,,,second_y
    Paint image,(np(4),np(5)+sx-t/2),colour,colour
    Paint image,(np(4)+sx-t/2,np(5)),colour,colour
End If
if arc1<>arc2 Then
    xp1=np(4)+(sx-t/2)*Cos(.5*(arc2+arc1))
yp1=np(5)-(sx-t/2)*Sin(.5*(arc2+arc1))
Circle image,(np(4),np(5)),sx,prime,arc1,arc2,second_y
    Circle image,(np(4),np(5)),sx-t,prime,arc1,arc2,second_y
    Line image,(np(4)+sx*Cos(arc1),np(5)-sx*Sin(arc1))-(np(4)+(sx-t)*Cos(arc1),np(5)-(sx-t)*Sin(arc1)),prime
    Line image,(np(4)+sx*Cos(arc2),np(5)-sx*Sin(arc2))-(np(4)+(sx-t)*Cos(arc2),np(5)-(sx-t)*Sin(arc2)),prime
    'pset(xp1,yp1),rgb(255,255,255)
    Paint image,(xp1,yp1),prime,prime

   Circle image,(np(4),np(5)),sx,colour,arc1,arc2,second_y
    Circle image,(np(4),np(5)),sx-t,colour,arc1,arc2,second_y
    Line image,(np(4)+sx*Cos(arc1),np(5)-sx*Sin(arc1))-(np(4)+(sx-t)*Cos(arc1),np(5)-(sx-t)*Sin(arc1)),colour
    Line image,(np(4)+sx*Cos(arc2),np(5)-sx*Sin(arc2))-(np(4)+(sx-t)*Cos(arc2),np(5)-(sx-t)*Sin(arc2)),colour
    'pset(xp1,yp1),rgb(255,255,255)
    Paint image,(xp1,yp1),colour,colour
End If
#endmacro

#macro get_perspective(np3,np6)
For n As Integer=3 To 6 Step 3
zval =np(n)  'for perspective
pp=perspective*((zval+1000)/1000-1)
pp=(1-pp)
If n=3 Then
np(n-2)=np(n-2)-pivot_x
np(n-1)=np(n-1)-pivot_y
np(n-2)=np(n-2)*pp
np(n-1)=np(n-1)*pp
np(n-2)=np(n-2)+pivot_x
np(n-1)=np(n-1)+pivot_y
Endif
If n=6 Then
    np(n-2)=np(n-2)-pivot_x
    np(n-1)=np(n-1)-pivot_y
    np(n-2)=np(n-2)*pp
    np(n-1)=np(n-1)*pp
    np(n-2)=np(n-2)+pivot_x
    np(n-1)=np(n-1)+pivot_y
Endif
Next n
sx=(pp)*sx
#endmacro

Dim pivot_vector(1 To 3) As Double
Dim line_vector(1 To 3) As Double
magnifier=dilator*magnifier
If shape="circle" Then
angleX=angleX Mod 360:angleY=angleY Mod 360:angleZ=angleZ Mod 360
End If
angleX_degrees=(2*p/360)*angleX     
angleY_degrees=(2*p/360)*angleY
angleZ_degrees=(2*p/360)*angleZ
pivot_vector(1)=first_x-pivot_x
pivot_vector(2)=first_y-pivot_y
pivot_vector(3)=first_z-pivot_z
pivot_vector(1)=dilator*pivot_vector(1)
pivot_vector(2)=dilator*pivot_vector(2)
pivot_vector(3)=dilator*pivot_vector(3)

Dim Rx(1 To 3,1 To 3) As Double
Dim Ry(1 To 3,1 To 3) As Double
Dim Rz(1 To 3,1 To 3) As Double
'rotat1on matrices about the three axix
If mode="3d" Then
Rx(1,1)=1:Rx(1,2)=0:Rx(1,3)=0
Rx(2,1)=0:Rx(2,2)=Cos(angleX_degrees):Rx(2,3)=-Sin(angleX_degrees)
Rx(3,1)=0:Rx(3,2)=Sin(angleX_degrees):Rx(3,3)=Cos(angleX_degrees)

Ry(1,1)=Cos(angleY_degrees):Ry(1,2)=0:Ry(1,3)=Sin(angleY_degrees)
Ry(2,1)=0:Ry(2,2)=1:Ry(2,3)=0
Ry(3,1)=-Sin(angleY_degrees):Ry(3,2)=0:Ry(3,3)=Cos(angleY_degrees)
Endif

Rz(1,1)=Cos(angleZ_degrees):Rz(1,2)=-Sin(angleZ_degrees):Rz(1,3)=0
Rz(2,1)=Sin(angleZ_degrees):Rz(2,2)=Cos(angleZ_degrees):Rz(2,3)=0
Rz(3,1)=0:Rz(3,2)=0:Rz(3,3)=1

line_vector(1)=magnifier*(second_x-first_x)'*pp                   'get the vector
line_vector(2)=magnifier*(second_y-first_y)'*pp                   'get the vector
line_vector(3)=magnifier*(second_z-first_z)'*pp

Dim new_pos(1 To 3) As Double
Dim temp1(1 To 3) As Double
Dim temp2(1 To 3) As Double
If mode="3d" Then
mv Rx(),pivot_vector(),temp1()           
mv Ry(),temp1(),temp2()
mv Rz(),temp2(),new_pos()
Endif
If mode="2d" Then
    mv Rz(),pivot_vector(),new_pos()
    Endif
new_pos(1)=new_pos(1)+pivot_x
new_pos(2)=new_pos(2)+pivot_y
new_pos(3)=new_pos(3)+pivot_z


Dim new_one(1 To 3) As Double            'To hold the turned value
If mode="3d" Then
mv Rx(),line_vector(),temp1()              'rotate
mv Ry(),temp1(),temp2()
mv Rz(),temp2(),new_one()
Endif
If mode="2d" Then
    mv Rz(),line_vector(),new_one()
    Endif
new_one(1)=new_one(1)+first_x              'translate
new_one(2)=new_one(2)+first_y
new_one(3)=new_one(3)+first_z

Dim xx As Double   
Dim yy As Double
Dim zz As Double
xx=first_x-new_pos(1)
yy=first_y-new_pos(2)
zz=first_z-new_pos(3)
 np(1)=new_one(1)-xx 
 np(2)=new_one(2)-yy
 np(3)=new_one(3)-zz
 np(4)=first_x-xx
 np(5)=first_y-yy
 np(6)= first_z-zz
If perspective <> 0 Then
get_perspective(np(3),np(6))
End If
Select Case shape
Case "line"
    If th<2 Then
 Line image,(np(4),np(5))-(np(1),np(2)),colour
Else
 thickline(th)   
 End If
Case "circle"
    Dim arc1 As Double=second_z*p/180
Dim arc2 As Double=second_arc*p/180
    If arc1=arc2 Then
    If th<=2 Then
 Circle image,(np(4),np(5)),magnifier*sx,colour,,,second_y
Else
 thickcircle(th)
End If
Endif
If arc1<>arc2 Then
If th<=2 Then
    Circle image,(np(4),np(5)),magnifier*sx,colour,arc1,arc2,second_y
Else
    thickcircle(th)
End If
End If
Case "circlefill"
    Dim As Double xp1,xp2,yp1,yp2
Dim As Uinteger prime=rgb(255,255,255)
Dim arc1 As Double=second_z*p/180
Dim arc2 As Double=second_arc*p/180
If arc1=arc2 Then Circle image,(np(4),np(5)),magnifier*sx,colour,,,second_y,F
If arc1<>arc2 Then
 xp1=np(4)+magnifier*sx*Cos(.5*(arc2+arc1))*3/4
yp1=np(5)-magnifier*sx*Sin(.5*(arc2+arc1))*3/4   
Circle image,(np(4),np(5)),magnifier*sx,prime,arc1,arc2,second_y
Line image,(np(4),np(5))-(np(4)+magnifier*sx*Cos(arc2),np(5)-magnifier*sx*Sin(arc2)),prime
Line image,(np(4),np(5))-(np(4)+magnifier*sx*Cos(arc1),np(5)-magnifier*sx*Sin(arc1)),prime
Paint image,(xp1,yp1),prime,prime

Circle image,(np(4),np(5)),magnifier*sx,colour,arc1,arc2,second_y
Line image,(np(4),np(5))-(np(4)+magnifier*sx*Cos(arc2),np(5)-magnifier*sx*Sin(arc2)),colour
Line image,(np(4),np(5))-(np(4)+magnifier*sx*Cos(arc1),np(5)-magnifier*sx*Sin(arc1)),colour
Paint image,(xp1,yp1),colour,colour
End If
 Case"box"
 
 Line image,(np(4),np(5))-(np(1),np(2)),colour,b
Case "boxfill"
 
 Line image,(np(4),np(5))-(np(1),np(2)),colour,bf
        Case "linepoint","circlepoint"
  'nothing drawn
Case "linepointset","circlepointset"
 If shape="linepointset" Then
 Pset image,(np(1),np(2)),colour
 Pset image,(np(4),np(5)),colour
 Endif
 If shape="circlepointset" Then
     Pset image,(np(4),np(5)),colour
 End If

        Case Else
 Print "unknown rotation shape"
End Select
End Sub
'END OF ROTATOR
   
 screenres 1000,700,32
 windowtitle "Rotating little squares around part of a sphere"
 sub walls(d as double)
     dim as double xpiv,ypiv,zpiv,scx,scy,scz,ax,ay,az,p,a,b
     dim as uinteger cl
     dim pi as double=4*atn(1)
     dim st as double
     dim as integer f=1
     Dim As Double cnp(1 To 8),cznp(1 To 4),cz(1 To f,1 To f) 'copy line end positions
dim as double cx,cy
'draw the walls
Dim As Double Xquad1(1 To 4),Yquad1(1 To 4)
Dim As Double Xquad2(1 To 4),Yquad2(1 To 4)
dim as double Xquad3(1 to 4),Yquad3(1 to 4)
dim as double Xquad4(1 to 4),Yquad4(1 to 4)
Dim As Uinteger leftside=rgb(20,20,20)
Dim As Uinteger rightside=rgb(25,25,25)
dim as uinteger back=rgb(22,22,22)
dim as uinteger bottom=rgb(21,21,21)
    Xquad1(1)=0:Yquad1(1)=300
    Xquad1(2)=0:Yquad1(2)=700
    Xquad1(3)=300:Yquad1(3)=450
    Xquad1(4)=300:Yquad1(4)=350
    drawpolygon(Xquad1(),Yquad1(),leftside)
    Xquad2(1)=1000:Yquad2(1)=300
    Xquad2(2)=1000:Yquad2(2)=700
    Xquad2(3)=700:Yquad2(3)=450
    Xquad2(4)=700:Yquad2(4)=350 
    drawpolygon(Xquad2(),Yquad2(),rightside)
    Xquad3(1)=300:Yquad3(1)=350
    Xquad3(2)=700:Yquad3(2)=350
    Xquad3(3)=700:Yquad3(3)=450
    Xquad3(4)=300:Yquad3(4)=450 
    drawpolygon(Xquad3(),Yquad3(),back)
    Xquad4(1)=0:Yquad4(1)=700
    Xquad4(2)=300:Yquad4(2)=450
    Xquad4(3)=700:Yquad4(3)=450
    Xquad4(4)=1000:Yquad4(4)=700
    drawpolygon(Xquad4(),Yquad4(),bottom)
     #macro edge(number)
 Select Case number
 'define two edges for cross product
 Case 1
   cnp(1)=np(1)
   cnp(2)=np(2)
   cnp(3)=np(4)
   cnp(4)=np(5)
   cznp(1)=np(3)
   cznp(2)=np(6)
 Case 2
   cnp(5)=np(1)
   cnp(6)=np(2)
   cnp(7)=np(4)
   cnp(8)=np(5)
   cznp(3)=np(3)
   cznp(4)=np(6)
   end select
   
cx=(cnp(1)+cnp(3)+cnp(5)+cnp(7))/4
cy=(cnp(2)+cnp(4)+cnp(6)+cnp(8))/4

 #endmacro
 dim as double u1,u2,u3,v1,v2,v3,wx,wy,wz,nw
 #macro crossproduct(of_two_sides)
 'get vectors to origin
 u1=cnp(1)-cnp(3)
 u2=cnp(2)-cnp(4)
 u3=cznp(1)-cznp(2)
 v1=cnp(5)-cnp(7)
 v2=cnp(6)-cnp(8)
 v3=cznp(3)-cznp(4)
 'get the cross product
 wx=(u2*v3-v2*u3)
 wy=-(u1*v3-v1*u3)
 wz=(u1*v2-v1*u2)
 nw=Sqr(wx^2+wy^2+wz^2)
 'normalized cross product components
 wx=wx/nw
 wy=wy/nw
 wz=wz/nw
 #endmacro
 
   #macro surface(cl1,cl2) 
        Paint(cx,cy),cl1,cl2
        #endmacro
   
     xpiv=350:ypiv=350:zpiv=5000:scx=500:scy=350:scz=0:ax=0:ay=12:az=0:p=0
     st=2*atn((d/zpiv)*180/pi)
    dim as double min=2,max =-2
    dim as integer ck
    dim as string action
    dim as uinteger lcl=rgb(0,1,0)

    for b=-4 to 4 step 1*st
        ax=b
         for a=-7 to 7 step st
     ck=250*(wz+.999999955)/(.999999955-.99012835910 )
     cl=rgb(10+ck,10+ck,10+ck)
     lcl=cl
         ay=a
       action="linepoint" 
    rotate(xpiv,ypiv,zpiv,scx-d,scy-d,scz,scx-d,scy+d,scz,.0,ax,ay,az,1,1,lcl,1,action,"3d",p)
     edge(1)
     rotate(xpiv,ypiv,zpiv,scx-d,scy-d,scz,scx+d,scy-d,scz,.0,ax,ay,az,1,1,lcl,1,action,"3d",p)
     edge(2)
     crossproduct(of_edges)
     If inpolygon(leftside,(cx),(cy)) Then
       cl=rgb(10,100+.6*ck,10+.5*ck)
       lcl=cl
       
    elseif inpolygon(rightside,(cx),(cy)) Then
       cl=rgb(10,100+.6*ck,10+.5*ck)
       lcl=cl
   elseif inpolygon(back,cx,cy) then
     cl=rgb(10,100+.5*ck,10+.4*ck)
       lcl=cl
   elseif inpolygon(bottom,cx,cy) then
     cl=rgb(50+.5*ck,50+.5*ck,0)
     lcl=cl
   else
       cl=rgb(10+ck,10+ck,10+ck)
       lcl=cl
       endif
     action="line"
       
     rotate(xpiv,ypiv,zpiv,scx-d,scy-d,scz,scx-d,scy+d,scz,.0,ax,ay,az,1,1,lcl,1,action,"3d",p)
     
     rotate(xpiv,ypiv,zpiv,scx-d,scy-d,scz,scx+d,scy-d,scz,.0,ax,ay,az,1,1,lcl,1,action,"3d",p)
     
     rotate(xpiv,ypiv,zpiv,scx+d,scy-d,scz,scx+d,scy+d,scz,.0,ax,ay,az,1,1,lcl,1,action,"3d",p)
 rotate(xpiv,ypiv,zpiv,scx+d,scy+d,scz,scx-d,scy+d,scz,.0,ax,ay,az,1,1,lcl,1,action,"3d",p)
         
        surface(cl,lcl)
next a
next b
draw string(450,200),"grade"
draw string(550,200),str(d)
'*************** back bars ********************
dim as double thk=5
dim as double yval=350
for count as integer=1 to 2
    if count=1 then yval=350
    if count=2 then yval=452
   for a as double=.08 to .034 step -.0001
    for b as double=0 to 180 step 1
rotate(500,500000,0,0,yval,0,0,yval,0,.0,0,0,a,1,1,rgb(255,255,255),1,"linepoint","2d",0)
rotate(np(1),np(2),0,np(1),np(2),0,np(1),np(2)-(thk+d),0,.0,0,0,b,1,1,rgb(255-b,255-b,200-b),1,"line","2d",0)
next b
next a
next count
' ********top side bar left **************************
yval=309
dim ma as double=1
 for a as double=.034 to 0 step -.0001
     ma=ma+.01
    for b as double=0 to 180 step 1
rotate(-70000,500000,0,0,yval,0,0,yval,0,.0,0,0,a,ma,1,rgb(255,255,255),1,"linepoint","2d",0)
rotate(np(1),np(2),0,np(1),np(2),0,np(1),np(2)-(thk+d),0,.0,0,0,b,ma,1,rgb(255-b,255-b,200-b),1,"line","2d",0)
next b
next a
' ********top right side bar  **************************
yval=448
 for a as double=.115 to .08 step -.0001
     ma=ma-.01
    for b as double=0 to 180 step 1
rotate(70000,500000,0,0,yval,0,0,yval,0,.0,0,0,a,ma,1,rgb(255,255,255),1,"linepoint","2d",0)
rotate(np(1),np(2),0,np(1),np(2),0,np(1),np(2)-(thk+d),0,.0,0,0,b,ma,1,rgb(255-b,255-b,200-b),1,"line","2d",0)
next b
next a
' ********lower side bar left **************************
yval=690
 for a as double=.034 to 0 step -.0001
     ma=ma+.01
    for b as double=0 to 180 step 1
rotate(400000,500000,0,0,yval,0,0,yval,0,.0,0,0,a,ma,1,rgb(255,255,255),1,"linepoint","2d",0)
rotate(np(1),np(2),0,np(1),np(2),0,np(1),np(2)-(thk+d),0,.0,0,0,b,ma,1,rgb(255-b,255-b,200-b),1,"line","2d",0)
next b
next a
' ********lower right bar left **************************
yval=-108

 for a as double=.113 to .08 step -.0001  '.08 to .12
     ma=ma-.01
    for b as double=0 to 180 step 1
rotate(-400000,500000,0,0,yval,0,0,yval,0,.0,0,0,a,ma,1,rgb(255,255,255),1,"linepointset","2d",0)
rotate(np(1),np(2),0,np(1),np(2),0,np(1),np(2)-(thk+d),0,.0,0,0,b,ma,1,rgb(255-b,255-b,200-b),1,"line","2d",0)

next b
next a

 
end sub
'************************** RUN ***************************
dim n as integer =27
 do
     n=n-3
     screenlock
     cls
 walls(n)
 screenunlock
 sleep 1,1
 loop until n=3
 draw string(450,250),"END, PRESS ANY KEY"

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

Postby D.J.Peters » Feb 23, 2010 14:42

~20 seconds with 6 vector and 3 photon refelctions

Joshy

Code: Select all

' PhotonMapping

#define true -1
#define false 0

type float as single ' double

Const As Integer xDimension          = 640
Const As Integer yDimension          = 480
Const As Integer nTypes              = 2
Const As Integer MaxPhotons          = 500
Const As Integer MaxPhotonRefections = 3
Const as float  PhotonSize           = 0.1
Const As float gInverseSquared       = 1.0 / PhotonSize*PhotonSize
Const As float  MilliSeconds         = 200
Const As float gExposure             = 1000 / MilliSeconds
Const As Integer MaxRaytraceRefections = 6

#define nSpheres 2
#define nPlanes  6
dim shared As float Ratio
Dim Shared As Integer gNObjects(1) = {nSpheres,nPlanes}
Dim Shared As float  gAmbient      = 0
Dim Shared As float  gOrigin(2)
Dim Shared As float  gLight(2)     = {0,0,4} 'x,y,z

Dim Shared As float  gSpheres(nSpheres-1,3) =  {_
                            {-1.0,-1.00,4.0, 0.5} , _
                            { 1.5,-1.8,3.9, 1} }
                             ' x  ,y    ,z, radius

Dim Shared As float  gPlanes(nPlanes-1,1)  = {_
                                      {0,  1.5}, _
                                      {0, -1.5}, _
                                      {1, -1.5}, _
                                      {1,  1.5}, _
                                      {2,  5.0}, _
                                      {2,  -1.0}}
Dim Shared As Integer PhotonMapping  = true
Dim Shared As Integer NumberOfPhotons(1,5)
Dim Shared As float   gPhotons(nTypes,nPlanes+nSpheres, _
                      MaxPhotons + MaxPhotons*MaxPhotonRefections, _
                      2,2) '3 * x,y,z
Dim Shared As Integer gIntersect
Dim Shared As Integer gType
Dim Shared As Integer gIndex
Dim Shared As float   gDist2,gDist,gDistSqr
Dim Shared As float   gPoint(2)

Function min(a As float ,b As float ) As float
  If b<a Then Return b
  Return a
End Function

Function max(a As float ,b As float ) As float
  If b>a Then Return b
  Return a
End Function

Function rnd2 As float
  Return (Rnd-Rnd)
End Function

Function SquaredDistance(a() As float , _
                   b() As float , _
                   d2 As float ) As Integer
  Dim As float  c = a(0) - b(0)
  Dim As float  d = c*c
  If (d > d2) Then Return false
  c = a(1) - b(1)
  d += c*c
  If (d > d2) Then Return false
  c = a(2) - b(2)
  d += c*c
  If (d > d2) Then Return false
  gDist2 = d'!!!
  gDistSqr=sqr(d)
  Return true
End Function

Sub SurfaceNormal(r() As float , _
                  typ As Integer, _
                  idx As Integer, _
                  HitPoint() As float , _
                  O() As float )
  dim as float  Normal(2),L
  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
    Dim As Integer 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 Then l=1/Sqr(l)
  r(0)=Normal(0)*l
  r(1)=Normal(1)*l
  r(2)=Normal(2)*l
End Sub

sub Mirror2(Ret()    as float, _
            Direction() as float, _
            HitPoint()  as float, _
            Normale()   as float)
  dim as float L
  HitPoint(0)+=Direction(0)*-0.1
  HitPoint(1)+=Direction(1)*-0.1
  HitPoint(2)+=Direction(2)*-0.1
  L=-2*(Normale(0)*Direction(0) + Normale(1)*Direction(1) + Normale(2)*Direction(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() As float ,_
              Direction() As float , _
              Origin() As float )
  Dim As float  L,Normal(2)
  SurfaceNormal(Normal(0),gType, gIndex, gPoint(0), Origin(0))
  gPoint(0)+=Direction(0)*-0.1
  gPoint(1)+=Direction(1)*-0.1
  gPoint(2)+=Direction(2)*-0.1
  L=-2*( Normal(0)*Direction(0)  _
        +Normal(1)*Direction(1) _
        +Normal(2)*Direction(2))
  Ret(0)=Direction(0)+L*Normal(0)
  Ret(1)=Direction(1)+L*Normal(1)
  Ret(2)=Direction(2)+L*Normal(2)

End Sub

'
' Raytracing
'
Sub Raytrace(RayDirection() As float ,RayOrigin() As float )
  gIntersect = false
  gDist      = 999999.9
 
  if gNObjects(1)>0 then
    For idx As Integer =0 To gNObjects(1)-1
      Dim As Integer axis = gPlanes(idx,0)
      If RayDirection(axis)<>0 Then
        Dim As float  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
  end if
 
  if gNObjects(0)>0 then
    For idx As Integer =0 To gNObjects(0)-1
      Dim As float  s(2)
      s(0)=RayOrigin(0)-gSpheres(idx,0)
      s(1)=RayOrigin(1)-gSpheres(idx,1)
      s(2)=RayOrigin(2)-gSpheres(idx,2)
      Dim As float  radius = gSpheres(idx,3)
      Dim As float  A = RayDirection(0)*RayDirection(0) _
                      + RayDirection(1)*RayDirection(1) _
                      + RayDirection(2)*RayDirection(2)
      Dim As float  B = 2.0 * (s(0)*RayDirection(0)+ _
                               s(1)*RayDirection(1)+ _
                               s(2)*RayDirection(2))
      Dim As float  C =       s(0)*s(0)+s(1)*s(1)+s(2)*s(2) -(radius*radius)
      Dim As float  D = B*B - 4*A*C
      If (D>0.0) Then
        Dim As float  sign  = iif(C < -0.00001,1,-1)
        Dim As float  l = (-B + sign*Sqr(D))/(2*A)
        If (l>0) And (l<gDist) Then
          gType  = 0
          gIndex = idx
          gDist  = l
          gIntersect = true
        End If
      End If
    next
  End If
End Sub

Sub AbsorbColor(ret() As float , _
                rgbIn() As float , _
                r As float ,g As float ,b As float )
  Dim As float  rgbOut(2)={r,g,b}
  For c As Integer =0 To 2
    If rgbOut(c)<rgbIn(c) Then
      ret(c)=rgbOut(c)
    Else
      ret(c)=rgbIn(c)
    End If
  Next
End Sub

Sub GetColor(r() As float , _
            rgbIn() As float , _
            typ As Integer, _
            idx As Integer)
  'AbsorbColor(r(0),rgbIn(0), 0.2, 0.2, 0.2)
  'return

  If (typ=0) Then     ' spheres
    if idx=0 then
      AbsorbColor(r(0),rgbIn(0), 0,0,0)
    elseif idx=1 then
      AbsorbColor(r(0),rgbIn(0), 0.2,0.2,1)
    end if
  Elseif (typ=1) Then ' planes
    If idx=0 Then
      AbsorbColor(r(0),rgbIn(0), 0.1, 0.8, 0.1) ' right
    Elseif idx=1 Then
      AbsorbColor(r(0),rgbIn(0), 0.8, 0.1, 0.1) ' left
    Elseif idx=2 Then
      AbsorbColor(r(0),rgbIn(0), 0.5, 0.5, 0.0) ' floor
    Elseif idx=3 Then
      AbsorbColor(r(0),rgbIn(0), 0.1, 0.1, 0.8) ' ceil
    Elseif idx=4 then
      AbsorbColor(r(0),rgbIn(0), 0.0, 0.0, 0.0) ' back side
    Elseif idx=5 then
      AbsorbColor(r(0),rgbIn(0), 0.0, 0.0, 1.0) ' behind camera
    End If
  End If
End Sub

'
' Photon Mapping
'
Sub GatherPhotons(energy() As float , _
                  HitPoint() As float , _
                  typ As Integer, _
                  idx As Integer)
  Dim As float  N(2)
  Dim As float  v(2)
  Dim As float  weight
  Dim As float  frgb(2)
  SurfaceNormal(N(0), typ, idx, HitPoint(0), gOrigin(0))

  For i As Integer = 0  To  NumberOfPhotons(typ,idx)-1
    ' 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)
   
    If (SquaredDistance(HitPoint(0),v(0),gInverseSquared)) Then
     
      dim as float  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.0, -cosin)
      weight *= (1.0 - gDistSqr) / gExposure
      ' energy
      frgb(0)+=gPhotons(typ,idx,i,2,0)*weight
      frgb(1)+=gPhotons(typ,idx,i,2,1)*weight
      frgb(2)+=gPhotons(typ,idx,i,2,2)*weight
    End If
  Next
  For j As Integer=0 To 2
    energy(j)=max(0,min(1,frgb(j) ) )
  Next
End Sub

Sub StorePhoton(typ As Integer, _
                idx As Integer, _
                l() As float ,_
                d() As float , _
                e() As float )
  Dim As Integer Photon=NumberOfPhotons(typ,idx)
 
  For i As Integer=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
  NumberOfPhotons(typ,idx)=Photon+1
End Sub

Sub ShadowPhoton(RayDirection() As float )
  static As float ShadowEnerg(2) = {-0.2,-0.2,-0.2}
  Dim As float   OldPoint(2) = {gPoint(0), gPoint(1),gPoint(2)}
  Dim As Integer OldType     = gType
  Dim As Integer OldIndex    = gIndex
  Dim as float   OldDist     = gDist
  Dim As float  BumpedPoint(2)
 
  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(0), BumpedPoint(0))

  Dim As float  ShadowPoint(2)

  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(0), RayDirection(0), ShadowEnerg(0))

  gPoint(0) = OldPoint(0)
  gPoint(1) = OldPoint(1)
  gPoint(2) = OldPoint(2)
  gType     = OldType
  gIndex    = OldIndex
  gDist     = OldDist ' !!!

End Sub

Sub EmitPhotons
  'Randomize 1
  Dim As float  PhotonEnergy(2)
  Dim As float  PhotonDirection(2)
  Dim As float  PhotonOrigin(2),l
  For typ As Integer = 0 To nTypes-1
    For idx As Integer = 0 To gNObjects(typ)-1
      NumberOfPhotons(typ,idx) = 0
    Next
  Next
 
  For i As Integer = 0 To MaxPhotons-1
    Dim As Integer Reflection = 0
    ' random photon Energy
    PhotonEnergy(0)=rnd
    PhotonEnergy(1)=rnd
    PhotonEnergy(2)=rnd
    l=PhotonEnergy(0)*PhotonEnergy(0) _
     +PhotonEnergy(1)*PhotonEnergy(1) _
     +PhotonEnergy(2)*PhotonEnergy(2)

    ' radom photon direction
    PhotonDirection(0)=rnd2
    PhotonDirection(1)=rnd2
    PhotonDirection(2)=rnd2
    l=PhotonDirection(0)*PhotonDirection(0) _
     +PhotonDirection(1)*PhotonDirection(1) _
     +PhotonDirection(2)*PhotonDirection(2)

    ' photon position origin from light
    PhotonOrigin(0)=gLight(0)
    PhotonOrigin(1)=gLight(1)
    PhotonOrigin(2)=gLight(2)
    ' random photon direction

    Raytrace(PhotonDirection(0), PhotonOrigin(0))
    While (gIntersect<>0) And (Reflection < MaxPhotonRefections)
      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(0),PhotonEnergy(0),gType,gIndex)
      #if 0
      Dim As float  l=1.0/Sqr(Reflection)
      PhotonEnergy(0)*=l
      PhotonEnergy(1)*=l
      PhotonEnergy(2)*=l
      #endif
      StorePhoton(gType, gIndex, gPoint(0), PhotonDirection(0),PhotonEnergy(0))
      ShadowPhoton(PhotonDirection(0))
      MirrorVec(PhotonDirection(0),PhotonDirection(0),PhotonOrigin(0))
      Raytrace(PhotonDirection(0), gPoint(0))
      PhotonOrigin(0)=gPoint(0)
      PhotonOrigin(1)=gPoint(1)
      PhotonOrigin(2)=gPoint(2)
    Wend
  Next
End Sub

Sub GetPixelColor(PixelRGB() As float , _
                  x As float ,y As float , z as float=1)
  Dim As float  RayDirection(2) = {x,y,z}
  Raytrace(RayDirection(0), gOrigin(0))
  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(0),gPoint(0),gType,gIndex)
   ' If gType=0 or gIndex>3 Then
    dim as integer nDivs=0,nReflection=0
    dim as float  MirrorsRGB(2)
    while ((gType=0 and gIndex=0) or gIndex>3) and gIntersect and _
          (nReflection<MaxRaytraceRefections)
      nReflection+=1
      MirrorVec(RayDirection(0),RayDirection(0),gOrigin(0))
      Raytrace(RayDirection(0), gPoint(0))
      If (gIntersect) Then
        dim as float  MirRGB(2)
        nDivs+=1
        gPoint(0)+=RayDirection(0)*gDist
        gPoint(1)+=RayDirection(1)*gDist
        gPoint(2)+=RayDirection(2)*gDist
        GatherPhotons(MirRGB(0),gPoint(0),gType,gIndex)
        MirrorsRGB(0)+=MirRGB(0)
        MirrorsRGB(1)+=MirRGB(1)
        MirrorsRGB(2)+=MirRGB(2)
      End If
    wend
    if nDivs>0 then
      MirrorsRGB(0)/=nDivs
      MirrorsRGB(1)/=nDivs
      MirrorsRGB(2)/=nDivs
      PixelRGB(0)=PixelRGB(0)*0.5+MirrorsRGB(0)*0.5
      PixelRGB(1)=PixelRGB(1)*0.5+MirrorsRGB(1)*0.5
      PixelRGB(2)=PixelRGB(2)*0.5+MirrorsRGB(2)*0.5
    end if
   
    'End If
  else
    PixelRGB(0)=1 ' !!! debug only !!!
    PixelRGB(1)=0
    PixelRGB(2)=1
  End If
End Sub

dim shared as integer interrupt
sub RenderRectangle(l as integer, t as integer,w as integer,h as integer)
  if interrupt=true then return
  if l>=xDimension then return
  if t>=yDimension then return
  if w=1 and h=1 then return
 
  Dim As Uinteger col
  Dim As float  b(2)
  dim as float  X =  (l/xDimension - 0.5)*Ratio
  dim as float  Y = -(t/yDimension - 0.5)
  GetPixelColor(b(0),x,y)
  line (l,t)-step(w-1,h-1),rgb(b(0)*255,b(1)*255,b(2)*255),BF   

  if w>1 then w=(w+0.5)/2
  if h>1 then h=(h+0.5)/2
  if w<1 then w=1
  if h<1 then h=1
  RenderRectangle(l  ,t  ,w,h)
  RenderRectangle(l+w,t  ,w,h)
  RenderRectangle(l  ,t+h,w,h)
  RenderRectangle(l+w,t+h,w,h)

 
  if inkey()<>"" then interrupt = true
end sub

Sub Render
  dim as integer h,m,s
  Dim As double t1=Timer
  WindowTitle " PhotonMapper rendering ..."
  RenderRectangle(0,0,xDimension,yDimension)
  s=timer-t1:h=s\3600:s-=h*3600:m=s\60:s-=m*60
  WindowTitle "PhotonMapper done " & h & ":" & m & ":" & s
End Sub

'
' main
'
Ratio= iif(xDimension>yDimension, _
           xDimension/yDimension, _
           yDimension/xDimension)

Randomize Timer
ScreenRes xDimension,yDimension,24
EmitPhotons
Render
'Bsave "PhotonMapper.bmp",0
Sleep
 
Last edited by D.J.Peters on Feb 23, 2010 19:09, edited 2 times in total.
duke4e
Posts: 717
Joined: Dec 04, 2005 0:16
Location: Varazdin, Croatia, Europe
Contact:

Postby duke4e » Feb 23, 2010 17:49

Cool, now speed it up with OpenCL ;)
D.J.Peters
Posts: 8023
Joined: May 28, 2005 3:28
Contact:

Postby D.J.Peters » Feb 23, 2010 18:59

duke4e wrote:Cool, now speed it up with OpenCL ;)
may be next month
it's cool to use a high level language first and than translate part by part to OpenCL.

Joshy
tinram
Posts: 88
Joined: Nov 30, 2006 13:35
Location: UK

Postby tinram » Feb 25, 2010 12:42

Superb Joshy.

There's an ethereal glow (interference) in your images above.
Double hit.
D.J.Peters
Posts: 8023
Joined: May 28, 2005 3:28
Contact:

Postby D.J.Peters » Oct 25, 2010 5:41

how fast is it on your boxes ?

Joshy
Image

Code: Select all

' PhotonMapping

#define true -1
#define false 0

Type float As single 'double

Const As Integer xDimension          = 512
Const As Integer yDimension          = 512
Const As Integer nTypes              = 2 ' two types of objects sphere and plane
Const As Integer MaxPhotons          = 200
Const As Integer MaxPhotonRefections = 1
Const As float  PhotonSize           = 0.1
Const As float gInverseSquared       = 1.0 / PhotonSize * PhotonSize
Const As float  MilliSeconds         = 500
Const As float gExposure             = 1.0 / (1000 / MilliSeconds)
Const As Integer MaxRaytraceRefections = 12

#define nSpheres 3 ' number of spheres
#define nPlanes  6 ' number of planes (a box)
Dim Shared As float Ratio
Dim Shared As Integer gNObjects(1) = {nSpheres,nPlanes}
Dim Shared As float  gAmbient      = 0
Dim Shared As float  gOrigin(2)
Dim Shared As float  gLight(2)     = {0,1.0,4} 'x,y,z position

Dim Shared As float  gSpheres(nSpheres-1,3) =  {_
                            {-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

Dim Shared As float  gPlanes(nPlanes-1,1)  = {_
                                      {0,  1.5}, _
                                      {0, -1.5}, _
                                      {1, -1.5}, _
                                      {1,  1.5}, _
                                      {2,  5.0}, _
                                      {2, -1.0}}

Dim Shared As Integer NumberOfPhotons(1,5)
Dim Shared As float   gPhotons(nTypes,nPlanes+nSpheres, _
                               MaxPhotons + MaxPhotons*MaxPhotonRefections, _
                               2,2) '3 * x,y,z
Dim Shared As Integer gIntersect
Dim Shared As Integer gType
Dim Shared As Integer gIndex
Dim Shared As float   gDist2,gDist
Dim Shared As float   gPoint(2)

Function min(a As float ,b As float ) As float
  If b<a Then Return b
  Return a
End Function

Function max(a As float ,b As float ) As float
  If b>a Then Return b
  Return a
End Function

Function rnd2 As float
  Return (Rnd-Rnd)
End Function

Function SquaredDistance(a() As float , _
                         b() As float , _
                         d2 As float ) As Integer
  Dim As float  ab = a(0) - b(0)
  Dim As float  d = ab*ab
  If (d > d2) Then Return false
  ab = a(1) - b(1)
  d += ab*ab
  If (d > d2) Then Return false
  ab = a(2) - b(2)
  d += ab*ab
  If (d > d2) Then Return false
  gDist2 = d
  Return true
End Function

Sub SurfaceNormal(r() As float , _
                  typ As Integer, _
                  idx As Integer, _
                  HitPoint() As float , _
                  O() As float )
  Dim As float  Normal(2),L
  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
    Dim As Integer 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 Then l=1/Sqr(l)
  r(0)=Normal(0)*l
  r(1)=Normal(1)*l
  r(2)=Normal(2)*l
End Sub

Sub Mirror2(Ret()       As float, _
            Direction() As float, _
            HitPoint()  As float, _
            Normale()   As float)
  Dim As float L
  HitPoint(0)+=Direction(0)*-0.1
  HitPoint(1)+=Direction(1)*-0.1
  HitPoint(2)+=Direction(2)*-0.1
  L = Normale(0)*Direction(0) _
    + Normale(1)*Direction(1) _
    + Normale(2)*Direction(2)
  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()       As float ,_
              Direction() As float , _
              Origin()    As float )
  Dim As float  L,Normale(2)
  SurfaceNormal(Normale(0),gType, gIndex, gPoint(0), Origin(0))
  gPoint(0)+=Direction(0)*-0.1
  gPoint(1)+=Direction(1)*-0.1
  gPoint(2)+=Direction(2)*-0.1
  L = Normale(0)*Direction(0) _
    + Normale(1)*Direction(1) _
    + Normale(2)*Direction(2)
  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() As float ,RayOrigin() As float )
  gIntersect = false
  gDist      = 999999.9
 
  If gNObjects(1)>0 Then
    For idx As Integer =0 To gNObjects(1)-1
      Dim As Integer axis = gPlanes(idx,0)
      If RayDirection(axis)<>0 Then
        Dim As float  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
  End If
 
  If gNObjects(0)>0 Then
    Dim As float  SphereDirection(2)
    Dim As float  A = RayDirection(0)*RayDirection(0) _
                    + RayDirection(1)*RayDirection(1) _
                    + RayDirection(2)*RayDirection(2)
                  A+=A
    For idx As Integer =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)
      Dim As float  R = gSpheres(idx,3)*gSpheres(idx,3)
      Dim As float  B = SphereDirection(0)*RayDirection(0) _
                      + SphereDirection(1)*RayDirection(1) _
                      + SphereDirection(2)*RayDirection(2)
                    B+=B
                   
      Dim As float  C = SphereDirection(0)*SphereDirection(0) _
                      + SphereDirection(1)*SphereDirection(1) _
                      + SphereDirection(2)*SphereDirection(2)
                    C-=R
      Dim As float  D = B*B - 2*A*C
      If (D>0.0) Then
        Dim As float  sign  = iif(C < -0.0001,1.0,-1.0)
        Dim As float  l = (-B + sign*Sqr(D))/A
        If (l>0.0) And (l<gDist) Then
          gType      = 0
          gIndex     = idx
          gDist      = l
          gIntersect = true
        End If
      End If
    Next
  End If
End Sub

Sub AbsorbColor(ret()   As float , _
                rgbIn() As float , _
                r As float ,g As float ,b As float )
  Dim As float  rgbOut(2)={r,g,b}
  For c As Integer =0 To 2
    If rgbOut(c)<rgbIn(c) Then
      ret(c)=rgbOut(c)
    Else
      ret(c)=rgbIn(c)
    End If
  Next
End Sub

Sub GetColor(r() As float , _
            rgbIn() As float , _
            typ As Integer, _
            idx As Integer)

  If (typ=0) Then     ' spheres
    If idx=0 Then
      AbsorbColor(r(0),rgbIn(0), 0.0,0.0,0.0)
    Elseif idx=1 Then
      AbsorbColor(r(0),rgbIn(0), 0.0,0.0,0.0)
    elseif idx=2 then
      AbsorbColor(r(0),rgbIn(0), 0.2,0.2,0.8)
    End If
  Elseif (typ=1) Then ' planes
    If idx=0 Then
      AbsorbColor(r(0),rgbIn(0), 0.1, 0.8, 0.1) ' right
    Elseif idx=1 Then
      AbsorbColor(r(0),rgbIn(0), 0.8, 0.1, 0.1) ' left
    Elseif idx=2 Then
      AbsorbColor(r(0),rgbIn(0), 0.5, 0.5, 0.0) ' floor
    Elseif idx=3 Then
      AbsorbColor(r(0),rgbIn(0), 0.2, 0.2, 0.2) ' ceil
    Elseif idx=4 Then
      AbsorbColor(r(0),rgbIn(0), 0.0, 0.0, 0.0) ' front
    Elseif idx=5 Then
      AbsorbColor(r(0),rgbIn(0), 0.0, 0.0, 0.0) ' behind camera
    End If
  End If
End Sub

'
' Photon Mapping
'
Sub GatherPhotons(energy() As float , _
                  HitPoint() As float , _
                  typ As Integer, _
                  idx As Integer)
  Dim As float  N(2)
  Dim As float  v(2)
  Dim As float  weight
  Dim As float  frgb(2)
  SurfaceNormal(N(0), typ, idx, HitPoint(0), gOrigin(0))

  For i As Integer = 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(0),v(0),gInverseSquared)) Then
      ' photon direction
      Dim As float  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.0, -cosin)
      if weight then
        weight *= (1.0 - sqr(gDist2)) * gExposure
        if weight then
          ' photon energy
          frgb(0)+=gPhotons(typ,idx,i,2,0)*weight
          frgb(1)+=gPhotons(typ,idx,i,2,1)*weight
          frgb(2)+=gPhotons(typ,idx,i,2,2)*weight
        end if
      end if
    End If
  Next
  For j As Integer=0 To 2
    energy(j)=max(0,min(1,frgb(j) ) )
  Next
End Sub

Sub StorePhoton(typ As Integer, _
                idx As Integer, _
                l() As float ,_
                d() As float , _
                e() As float )
  Dim As Integer Photon=NumberOfPhotons(typ,idx)
 
  For i As Integer=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
  NumberOfPhotons(typ,idx)=Photon+1
End Sub

Sub ShadowPhoton(RayDirection() As float )
  static As float  BumpedPoint(2)
  static As float  ShadowPoint(2)
  Static As float ShadowEnerg(2) = {-0.2,-0.2,-0.2}
  Dim As float   OldPoint(2) = {gPoint(0), gPoint(1),gPoint(2)}
  Dim As Integer OldType     = gType
  Dim As Integer OldIndex    = gIndex
  Dim As float   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(0), BumpedPoint(0))

  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(0), RayDirection(0), ShadowEnerg(0))

  gPoint(0) = OldPoint(0)
  gPoint(1) = OldPoint(1)
  gPoint(2) = OldPoint(2)
  gType     = OldType
  gIndex    = OldIndex
  gDist     = OldDist

End Sub

Sub EmitPhotons
  randomize 1 'timer
  dim As float PhotonEnergy(2)
  dim As float PhotonDirection(2)
  dim As float PhotonOrigin(2),l
  For typ As Integer = 0 To nTypes-1
    For idx As Integer = 0 To gNObjects(typ)-1
      NumberOfPhotons(typ,idx) = 0
    Next
  Next
  For i As Integer = 0 To MaxPhotons-1
    Dim As Integer Reflection = 0
    ' random photon Energy
    PhotonEnergy(0)=rnd
    PhotonEnergy(1)=rnd
    PhotonEnergy(2)=rnd
#if 0
    ' normalize energy
    l = PhotonEnergy(0)*PhotonEnergy(0) _
      + PhotonEnergy(1)*PhotonEnergy(1) _
      + PhotonEnergy(2)*PhotonEnergy(2)
    if l then
      l=1.0/sqr(l)
      PhotonEnergy(0)*=l
      PhotonEnergy(1)*=l
      PhotonEnergy(2)*=l
    end if
#endif

    ' radom photon Direction
    PhotonDirection(0)= rnd2
    PhotonDirection(1)= rnd2*2
    PhotonDirection(2)= rnd2
#if 0
    ' normalize direction
    l = PhotonDirection(0)*PhotonDirection(0) _
      + PhotonDirection(1)*PhotonDirection(1) _
      + PhotonDirection(2)*PhotonDirection(2)
    if l then
      l=1.0/sqr(l)
      PhotonDirection(0)*=l
      PhotonDirection(1)*=l
      PhotonDirection(2)*=l
    end if
#endif

    ' photon position origin from light
    PhotonOrigin(0)=gLight(0)
    PhotonOrigin(1)=gLight(1)
    PhotonOrigin(2)=gLight(2)

    Raytrace(PhotonDirection(0), PhotonOrigin(0))
    While (gIntersect<>0) And (Reflection < MaxPhotonRefections)
      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(0),PhotonEnergy(0),gType,gIndex)
#if 0
      Dim As float l=1.0/Reflection
      PhotonEnergy(0)*=l
      PhotonEnergy(1)*=l
      PhotonEnergy(2)*=l
#endif
      StorePhoton(gType, gIndex, gPoint(0), PhotonDirection(0),PhotonEnergy(0))
      ShadowPhoton(PhotonDirection(0))
      MirrorVec(PhotonDirection(0),PhotonDirection(0),PhotonOrigin(0))
      Raytrace(PhotonDirection(0), gPoint(0))
      PhotonOrigin(0)=gPoint(0)
      PhotonOrigin(1)=gPoint(1)
      PhotonOrigin(2)=gPoint(2)
    Wend
  Next
End Sub

Sub GetPixelColor(PixelRGB() As float , _
                  x As float ,y As float , z As float=1)
  Dim As float  RayDirection(2) = {x,y,z}
  Raytrace(RayDirection(0), gOrigin(0))
  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(0),gPoint(0),gType,gIndex)
   ' If gType=0 or gIndex>3 Then
    Dim As Integer nDivs=0,nReflection=0
    Dim As float  MirrorsRGB(2)
    While ((gType=0 And gIndex<2) Or gIndex>3) And gIntersect And _
          (nReflection<MaxRaytraceRefections)
      nReflection+=1
      MirrorVec(RayDirection(0),RayDirection(0),gOrigin(0))
      Raytrace(RayDirection(0), gPoint(0))
      If (gIntersect) Then
        Dim As float  MirRGB(2)
        nDivs+=1
        gPoint(0)+=RayDirection(0)*gDist
        gPoint(1)+=RayDirection(1)*gDist
        gPoint(2)+=RayDirection(2)*gDist
        GatherPhotons(MirRGB(0),gPoint(0),gType,gIndex)
        MirrorsRGB(0)+=MirRGB(0)
        MirrorsRGB(1)+=MirRGB(1)
        MirrorsRGB(2)+=MirRGB(2)
      End If
    Wend
    If nDivs>0 Then
      MirrorsRGB(0)/=nDivs
      MirrorsRGB(1)/=nDivs
      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
   
    'End If
  Else
    PixelRGB(0)=1 ' !!! debug only !!!
    PixelRGB(1)=0
    PixelRGB(2)=1
  End If
End Sub


Sub Render
  Dim As Integer h,m,s,l,t
  dim as float b(2),x,y
  Dim As Double t1=Timer
  WindowTitle " PhotonMapper rendering ..."
  for t =0 to yDimension-1
    Y = -(t/yDimension - 0.5)
    screenlock
    for l =0 to xDimension-1
      X =  (l/xDimension - 0.5)*Ratio
      GetPixelColor(b(0),x,y)
      pset (l,t),rgb(b(0)*255,b(1)*255,b(2)*255)
    next
    screenunlock
  next
  s=timer-t1:h=s\3600:s-=h*3600:m=s\60:s-=m*60
  WindowTitle "PhotonMapper done " & h & ":" & m & ":" & s
End Sub

'
' main
'
Ratio= iif(xDimension>yDimension, _
           xDimension/yDimension, _
           yDimension/xDimension)

Randomize Timer
ScreenRes xDimension,yDimension,24
EmitPhotons
Render
'Bsave "PhotonMapper.bmp",0
Sleep
 
Last edited by D.J.Peters on Sep 25, 2017 16:41, edited 1 time in total.
h4tt3n
Posts: 693
Joined: Oct 22, 2005 21:12
Location: Denmark

Postby h4tt3n » Oct 25, 2010 6:53

D.J.Peters wrote:how fast is it on your boxes ?


2 secs on my dual core 2 ghz HP laptop

Cheers,
Mike
TJF
Posts: 3546
Joined: Dec 06, 2009 22:27
Location: N47°, E15°
Contact:

Postby TJF » Oct 25, 2010 7:17

AMD Sempron 3800+:
    2 sec (fbc -gen gcc -O max)
    5 sec (fbc -gen gas)
D.J.Peters
Posts: 8023
Joined: May 28, 2005 3:28
Contact:

Postby D.J.Peters » Oct 25, 2010 9:52

Thanks for testing i found the problem
there was a hidden running thread on my box and i got very slow results
but now all works fine.

test of moving photons

http://www.youtube.com/playlist?list=PL ... re=viewall

Joshy
Last edited by D.J.Peters on Sep 06, 2011 1:39, edited 1 time in total.

Return to “Tips and Tricks”

Who is online

Users browsing this forum: No registered users and 4 guests