RayTracing versus PhotonMapping

Post your FreeBASIC source, examples, tips and tricks here. Please don’t post code without including an explanation.
Post Reply
rolliebollocks
Posts: 2655
Joined: Aug 28, 2008 10:54
Location: new york

Post by rolliebollocks »

5 seconds, AMD Phenom II
Top
APIStudios
Posts: 15
Joined: Oct 17, 2010 10:13
Contact:
Contact APIStudios

Post by APIStudios »

Just for fun, we optimized it a bit to use multiple threads (To take advantage of multi-core CPUs) then we ran it on a 256 core server cluster.

Ran too fast to measure, increased the resolution to 4096x4096 and it still rendered in under 2 sec.

361KB Medium quality JPG (Click to open):
Image
2'337 KB Pixel perfect PNG (Click to open):
Image

Neat program, photon mapping is awesome!
Top
joseywales72
Posts: 206
Joined: Aug 27, 2005 2:02
Location: Istanbul, Turkey

Post by joseywales72 »

2 seconds with fbc -gen gas, 1 second with fbc -gen gcc -O max
Archlinux 32, 2.6.35 kernel, AMD Athlon X2 5200

Very nice looking shot...
Top
D.J.Peters
Posts: 8631
Joined: May 28, 2005 3:28
Contact:
Contact D.J.Peters

Post by D.J.Peters »

some more short videos from same source code

http://www.youtube.com/playlist?list=PL ... re=viewall
Top
HACK3R ADI
Posts: 27
Joined: Jun 17, 2012 10:16

Re: RayTracing versus PhotonMapping

Post by HACK3R ADI »

ample will help me most in OS development....Bookmarked
Top
dafhi
Posts: 1741
Joined: Jun 04, 2005 9:51

Re: RayTracing versus PhotonMapping

Post by dafhi »

unmangled

Code: Select all

' PhotonMapping

#ifndef true
#define true -1
#define false 0
  #endif

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(),gType, gIndex, gPoint(), Origin())
  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(),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() 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(), typ, idx, HitPoint(), gOrigin())

  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(),v(),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(), 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
  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(), PhotonOrigin())
    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(),PhotonEnergy(),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(), 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
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(), 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)
   ' 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(),RayDirection(),gOrigin())
      Raytrace(RayDirection(), gPoint())
      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(),gPoint(),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(),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
Top
D.J.Peters
Posts: 8631
Joined: May 28, 2005 3:28
Contact:
Contact D.J.Peters

Re: RayTracing versus PhotonMapping

Post by D.J.Peters »

After some changes the old code woks with current fbc version now.

Joshy

Code: Select all

' PhotonMapping

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(),gType, gIndex, gPoint(), Origin())
  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(),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() 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(), typ, idx, HitPoint(), gOrigin())

  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(),v(),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(), 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
  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(), PhotonOrigin())
    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(),PhotonEnergy(),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(), 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
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(), 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)
   ' 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(),RayDirection(),gOrigin())
      Raytrace(RayDirection(), gPoint())
      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(),gPoint(),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(),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
Top
UEZ
Posts: 1079
Joined: May 05, 2017 19:59
Location: Germany

Re: RayTracing versus PhotonMapping

Post by UEZ »

Looks wonderful - awesome!

Thanks for the update.
Top
jj2007
Posts: 2326
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:
Contact jj2007

Re: RayTracing versus PhotonMapping

Post by jj2007 »

Awesome indeed, thumbs up!
Top
TheRaven
Posts: 10
Joined: Mar 09, 2019 18:42

Re: RayTracing versus PhotonMapping

Post by TheRaven »

2019 --much later, and...still AWESOME!
Top
Luxan
Posts: 273
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: RayTracing versus PhotonMapping

Post by Luxan »

For ages now I've been using BASIC of one type or another, always I've been
attempting to incorporate photo realistic scenes of data and objects in general.

This code of yours is a significant step towards that.

With xDimension = 1000 and yDimension = 800, this renders in less than a second,
rather incredible .
Top
Post Reply

Return to “Sources, Examples, Tips and Tricks”