purept (smallpt) path tracer

Programmation d'applications complexes
Avatar de l’utilisateur
grabiller
Messages : 103
Inscription : lun. 10/sept./2012 11:55
Localisation : France - 89220 Rogny-Les-Septs-Ecluses
Contact :

purept (smallpt) path tracer

Message par grabiller »

20130824-21h07: Mise à jour
- Récursion retirée.
- 8 scènes ajoutées.
- meilleurs infos à l'image et en console.

Hello,

Explorant la méthode de rendu en path tracing j'ai porté smallpt en PureBasic, à titre d'expérience.

Image

C'est vraiment l'un des plus basiques path tracer que l'on peut trouver mais ça marche et - bien que plus lent qu'en C/C++ - les performances restent respectables et encore le code est loin d'être optimisé.

La version originale utilise OpenMP pour le multithreading, celle-ci est 'manuellement' multithreadée et utilise tous les cores disponibles.

Il y a un problème par contre, vu que smallpt utilise la récursion pour sa fonction 'radiance', avec la génération d'une erreur 'Stack Overflow' au bout d'un moment, alors j'ai limité la récursion à 5 dans la section roulette russe pour la terminaison du rayon. A part cela la version PureBasic a toutes les features de la version originale.

Je n'ai pas essayé de porter la fonction erand48(Xi) et j'utilise un générateur de nombres aléatoires (très) basique.

Certaines améliorations sont disponibles sur le site de Kevin Beason que je compte ajouter ultérieurement (sans recursion, sampling de lumière explicite, autres scenes, ..). Plus tard j'essaierai également de créer une version OpenCL.

A propos des performances, Kevin a publié quelques statistiques avec sa config de l'époque:
~5m pour 200 spp (render(50) dans le code plus bas) avec un 2.4 GHz Intel Core 2 Quad CPU utilisant 4 threads.

Si quelqu'un a encore une config similaire, je serai curieux de connaitre les temps de rendu pour comparaison.

Sur mon i7 x980 3.2 GHz j'obtiens ~1m43s avec 12 threads.

N'oubliez pas de compiler le code avec le switch "Create threadsafe executable" activé!

Code : Tout sélectionner

EnableExplicit

OpenConsole("purept")

; ---[ 3d Vector ]------------------------------------------------------------
Structure v3d
  x.d
  y.d
  z.d
EndStructure
Macro v3Null( v_ )
  v_\x = 0.0
  v_\y = 0.0
  v_\z = 0.0
EndMacro
Macro v3Clamp( v_ )
  v_\x = clamp(v_\x)
  v_\y = clamp(v_\y)
  v_\z = clamp(v_\z)
EndMacro
Macro v3Norm( v_, n_ )
  n_ = Sqr( v_\x*v_\x + v_\y*v_\y + v_\z*v_\z )
  If n_ <> 0.0
    n_ = 1.0/n_
    v_\x*n_
    v_\y*n_
    v_\z*n_
  EndIf
EndMacro
Macro v3Dot( v1_, v2_ )
  ( v1_\x * v2_\x + v1_\y * v2_\y + v1_\z * v2_\z )
EndMacro
Macro v3Copy( v_, v1_ )
  v_\x = v1_\x
  v_\y = v1_\y
  v_\z = v1_\z
EndMacro
Macro v3Add( v_, v1_, v2_ )
  v_\x = v1_\x + v2_\x
  v_\y = v1_\y + v2_\y
  v_\z = v1_\z + v2_\z
EndMacro
Macro v3AddScaledInPlace( v_, v1_, s1_ )
  v_\x + v1_\x * s1_
  v_\y + v1_\y * s1_
  v_\z + v1_\z * s1_
EndMacro
Macro v3Sub( v_, v1_, v2_ )
  v_\x = v1_\x - v2_\x
  v_\y = v1_\y - v2_\y
  v_\z = v1_\z - v2_\z
EndMacro
Macro v3Mul( v_, v1_, v2_ )
  v_\x = v1_\x * v2_\x
  v_\y = v1_\y * v2_\y
  v_\z = v1_\z * v2_\z
EndMacro
Macro v3MulInPlace( v_, v1_ )
  v_\x * v1_\x
  v_\y * v1_\y
  v_\z * v1_\z
EndMacro
Macro v3Scale( v_, v1_, s1_ )
  v_\x = v1_\x * s1_
  v_\y = v1_\y * s1_
  v_\z = v1_\z * s1_
EndMacro
Macro v3ScaleInPlace( v_, s_ )
  v_\x * s_
  v_\y * s_
  v_\z * s_
EndMacro
Macro v3Cross( v_, v1_, v2_ )
  v_\x = v1_\y * v2_\z - v1_\z * v2_\y
  v_\y = v1_\z * v2_\x - v1_\x * v2_\z
  v_\z = v1_\x * v2_\y - v1_\y * v2_\x
EndMacro
Macro v3Neg( v_, v1_ )
  v_\x = -v1_\x
  v_\y = -v1_\y
  v_\z = -v1_\z
EndMacro

; ---[ 3d Ray ]---------------------------------------------------------------
Structure r3d
  o.v3d
  d.v3d
EndStructure

; ---[ Material ]-------------------------------------------------------------
Enumeration
  #DIFF
  #SPEC
  #REFR
EndEnumeration

; ---[ Sphere ]---------------------------------------------------------------
Structure s3d
  rad.d ; radius
  p.v3d ; position
  e.v3d ; emission
  c.v3d ; color
  m.i   ; material: #DIFF/#SPEC/#REFR
EndStructure
Procedure.d s3dIntersect( *sphere.s3d, *ray.r3d )
  ; returns distance, 0 if no hit
  ; Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0
  
  ; Vec op = p - r.o
  Protected op.v3d
  v3Sub( op, *sphere\p, *ray\o )
  
  ; double t, eps=1e-4, b=op.dot(r.d), det=b*b-op.dot(op)+rad*rad;
  Protected t.d
  Protected eps.d = 1.0e-4
  Protected b.d = v3Dot( op, *ray\d )
  Protected det.d = b*b - v3Dot( op, op ) + ( *sphere\rad * *sphere\rad )
  
  ; if (det<0) return 0; else det=sqrt(det);
  If det < 0.0
    ProcedureReturn 0.0
  Else
    det = Sqr(det)
  EndIf
  
  ; return (t=b-det)>eps ? t : ((t=b+det)>eps ? t : 0);
  t = b - det
  If t <= eps
    t = b + det
    If t <= eps
      ProcedureReturn 0.0
    EndIf
  EndIf
  
  ProcedureReturn t

EndProcedure

; ---[ Scene ]----------------------------------------------------------------
Macro Sphere( idx_, rad_, posx_, posy_, posz_, emir_, emig_, emib_, colr_, colg_, colb_, mat_ )
  scene(idx_)\rad = rad_
  scene(idx_)\p\x = posx_
  scene(idx_)\p\y = posy_
  scene(idx_)\p\z = posz_
  scene(idx_)\e\x = emir_
  scene(idx_)\e\y = emig_
  scene(idx_)\e\z = emib_
  scene(idx_)\c\x = colr_
  scene(idx_)\c\y = colg_
  scene(idx_)\c\z = colb_
  scene(idx_)\m = mat_
EndMacro
Global Dim scene.s3d(8)
Sphere( 0, 1.0e5,   1.0e5+1.0,  40.8,       81.6,         0.0, 0.0, 0.0,   0.75,  0.25,  0.25,  #DIFF ) ; Left
Sphere( 1, 1.0e5,  -1.0e5+99.0, 40.8,       81.6,         0.0, 0.0, 0.0,   0.25,  0.25,  0.75,  #DIFF ) ; Rght
Sphere( 2, 1.0e5,   50.0,       40.8,       1.0e5,        0.0, 0.0, 0.0,   0.75,  0.75,  0.75,  #DIFF ) ; Back
Sphere( 3, 1.0e5,   50.0,       40.8,      -1.0e5+170.0,  0.0, 0.0, 0.0,   0.0,   0.0,   0.0,   #DIFF ) ; Frnt
Sphere( 4, 1.0e5,   50.0,       1.0e5,      81.6,         0.0, 0.0, 0.0,   0.75,  0.75,  0.75,  #DIFF ) ; Botm
Sphere( 5, 1.0e5,   50.0,      -1.0e5+81.6, 81.6,         0.0, 0.0, 0.0,   0.75,  0.75,  0.75,  #DIFF ) ; Top
Sphere( 6, 16.5,    27.0,       16.5,       47.0,         0.0, 0.0, 0.0,   0.999, 0.999, 0.999, #SPEC ) ; Mirr
Sphere( 7, 16.5,    73.0,       16.5,       78.0,         0.0, 0.0, 0.0,   0.999, 0.999, 0.999, #REFR ) ; Glas
Sphere( 8, 600.0,   50.0,       681.6-0.27, 81.6,        12.0,12.0,12.0,   0.0,   0.0,   0.0,   #DIFF ) ; Lite

; ---[ Helpers ]--------------------------------------------------------------
Procedure.d clamp( x.d )
  If x < 0.0
    ProcedureReturn 0.0
  ElseIf x > 1.0
    ProcedureReturn 1.0
  EndIf
  ProcedureReturn x
EndProcedure
Procedure.l toInt( x.d )
  ; linear to ~sRGB (not accurate), yet the spheres colors are not
  ; linearized prior rendering.
  ProcedureReturn Int( Pow(clamp(x),1.0/2.2)*255+0.5 )
EndProcedure
Procedure.b intersect( *ray.r3d, *t.Double, *id.Long )
  Protected d.d, inf.d = 1e20
  Protected i.i = ArraySize(scene())
  *t\d = inf
  While i >= 0
    d = s3dIntersect( scene(i), *ray )
    If d And d<*t\d
      *t\d  = d
      *id\l = i
    EndIf
    i - 1
  Wend
  ProcedureReturn Bool( *t\d < inf )
EndProcedure
Procedure.d rand()
  ProcedureReturn Random(16384,0)/16384.0
EndProcedure

; ---[ Radiance ]-------------------------------------------------------------
Procedure radiance( *ret.v3d, *r.r3d, depth.l, *Xi.Unicode )
  
  Protected ntemp.d
  Protected tt_.d
  Protected r2_.r3d
  
  Protected t .d = 0.0 ; distance to intersection
  Protected id.l = 0   ; id of intersected object

  ; if miss, return black
  If Not intersect( *r, @t, @id )
    v3Null(*ret)
    ProcedureReturn
  EndIf

  Protected *obj.s3d = scene(id) ; the hit object
  
  ; x = r.o + r.d*t
  Protected x.v3d
  x\x = *r\o\x + *r\d\x*t
  x\y = *r\o\y + *r\d\y*t
  x\z = *r\o\z + *r\d\z*t
  
  ; n=(x-obj.p).norm()
  Protected n.v3d
  v3Sub( n, x, *obj\p )
  v3Norm( n, ntemp )
  
  ; nl=n.dot(r.d)<0?n:n*-1
  Protected nl.v3d
  If v3Dot( n, *r\d ) < 0.0
    v3Copy( nl, n )
  Else
    v3Neg( nl, n )
  EndIf
  
  ; f=obj.c
  Protected f.v3d
  v3Copy( f, *obj\c )
  
  ; max refl
  Protected p.d
  If ( f\x > f\y ) And ( f\x > f\z )
    p = f\x
  ElseIf f\y > f\z
    p = f\y
  Else
    p = f\z
  EndIf

  ; Russian Roulette
  depth + 1
  If depth > 5
    ;If rand() < p ; erand48(Xi) < p
      ;f\x/p
      ;f\y/p
      ;f\z/p
    ;Else
      v3Copy( *ret, *obj\e )
      ProcedureReturn
    ;EndIf
  EndIf
  
  Protected rr.v3d
  Protected rt.r3d
  Protected rd.v3d
  
  ; Ideal DIFFUSE reflection
  If *obj\m = #DIFF

    Protected r1  .d = 2.0*#PI*rand()
    Protected r2  .d = rand()
    Protected r2s .d = Sqr(r2)
    Protected w.v3d
    v3Copy( w, nl )
    Protected ut.v3d,u.v3d
    If Abs( w\x ) > 0.1
      ut\x = 0.0
      ut\y = 1.0
      ut\z = 0.0
    Else
      ut\x = 1.0
      ut\y = 0.0
      ut\z = 0.0
    EndIf
    v3Cross(u,ut,w)
    v3Norm(u,ntemp)
    Protected v.v3d
    v3Cross(v,w,u)
    
    Protected d.v3d
    d\x = u\x*Cos(r1)*r2s + v\x*Sin(r1)*r2s + w\x*Sqr(1.0-r2)
    d\y = u\y*Cos(r1)*r2s + v\y*Sin(r1)*r2s + w\y*Sqr(1.0-r2)
    d\z = u\z*Cos(r1)*r2s + v\z*Sin(r1)*r2s + w\z*Sqr(1.0-r2)
    v3Norm(d,ntemp)
    
    ; return obj.e + f.mult( radiance( Ray(x,d),depth,Xi) );
    v3Copy( rt\o, x )
    v3Copy( rt\d, d )
    radiance( @rr, @rt, depth, 0 )
    v3MulInPlace( f, rr )    
    v3Add( *ret, *obj\e, f )
    ProcedureReturn
    
  ; Ideal SPECULAR reflection
  ElseIf *obj\m = #SPEC

    ; return obj.e + f.mult(radiance( Ray( x, r.d-n*2*n.dot(r.d) ), depth, Xi ) );
    v3Copy( rt\o, x )
    v3Copy(rd,n)
    tt_ = 2.0*v3Dot(n,*r\d)
    v3ScaleInPlace( rd, tt_ )
    v3Sub( rt\d, *r\d, rd )
    radiance( @rr, @rt, depth, 0 )
    v3MulInPlace( f, rr )
    v3Add( *ret, *obj\e, f )
    ProcedureReturn
    
  EndIf
  
  ; Ideal dielectric REFRACTION
  ; Ray reflRay( x, r.d - n*2*n.dot(r.d) );

  Protected reflRay.r3d
  v3Copy( reflRay\o, x )
  v3Copy(rd,n)
  tt_ = 2.0*v3Dot(n,*r\d)
  v3ScaleInPlace( rd, tt_ )
  v3Sub( reflRay\d, *r\d, rd )
  ; Ray from outside going in? 
  ; bool into = n.dot(nl)>0;
  Protected into.b = Bool( v3Dot(n,nl) > 0.0 )
  ; double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t;
  Protected nc.d=1.0,nt.d=1.5,nnt.d,ddn.d,cos2t.d
  If into
    nnt = nc/nt
  Else
    nnt = nt/nc
  EndIf
  ddn = v3Dot(*r\d,nl)
  cos2t = 1.0 - nnt*nnt*( 1.0 - ddn*ddn )
  If cos2t < 0 ; Total internal reflection
    ; return obj.e + f.mult(radiance(reflRay,depth,Xi));
    radiance( @rr, @reflRay, depth, 0 )
    v3MulInPlace( f, rr )
    v3Add( *ret, *obj\e, f )
    ProcedureReturn
  EndIf
  ; Vec tdir = ( r.d*nnt - n*(  (into?1:-1)*(ddn*nnt+sqrt(cos2t))  ) ).norm();
  Protected intod.d
  If into
    intod = 1.0
  Else
    intod = -1.0
  EndIf
  Protected t_.d,n_.v3d,r_.v3d
  t_ = intod*(ddn*nnt+Sqr(cos2t))
  v3Copy(n_,n)
  v3ScaleInPlace( n_, t_ )
  v3Copy(r_,*r\d)
  v3ScaleInPlace( r_, nnt )
  Protected tdir.v3d
  v3Sub(tdir,r_,n_)
  v3Norm(tdir,ntemp)
  ; double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c = 1-(into?-ddn:tdir.dot(n));
  Protected a.d=nt-nc,b.d=nt+nc,R0.d=a*a/(b*b),c.d
  If into
    c = 1.0 + ddn
  Else
    c = 1.0 - v3Dot(tdir,n)
  EndIf
  ; double Re=R0+(1-R0)*c*c*c*c*c,Tr=1-Re,P=.25+.5*Re,RP=Re/P,TP=Tr/(1-P);
  Protected Re.d=R0+(1.0-R0)*c*c*c*c*c
  Protected Tr.d=1.0-Re
  Protected P_.d=0.25+0.5*Re
  Protected RP.d=Re/P_
  Protected TP.d=Tr/(1.0-P_)
  ; return obj.e + f.mult(depth>2 ? ( erand48(Xi)<P ?   // Russian roulette
  ; radiance(reflRay,depth,Xi)*RP : radiance(Ray(x,tdir),depth,Xi)*TP ) :
  ; radiance(reflRay,depth,Xi)*Re+radiance(Ray(x,tdir),depth,Xi)*Tr);
  If depth > 2
    If rand() < P_
      ; radiance(reflRay,depth,Xi)*RP
      radiance( @rr, @reflRay, depth, 0 )
      v3ScaleInPlace( rr, RP )
    Else
      ; radiance(Ray(x,tdir),depth,Xi)*TP )
      v3Copy(r2_\o,x)
      v3Copy(r2_\d,tdir)
      radiance( @rr, @r2_, depth, 0 )
      v3ScaleInPlace( rr, TP )
    EndIf
  Else
    Protected rr1.v3d
    Protected rr2.v3d
    ; radiance(reflRay,depth,Xi)*Re + radiance(Ray(x,tdir),depth,Xi)*Tr)
    radiance( @rr1, @reflRay, depth, 0 )
    v3ScaleInPlace( rr1, Re )
    v3Copy(r2_\o,x)
    v3Copy(r2_\d,tdir)
    radiance( @rr2, @r2_, depth, 0 )
    v3ScaleInPlace( rr1, Tr )
    v3Add(rr,rr1,rr2)
  EndIf
  v3Mul( f, f, rr )
  v3Add( *ret, *obj\e, f )
  ProcedureReturn
  
EndProcedure

; ---[ Render ]---------------------------------------------------------------
; ...[ Rendering Context ]....................................................
Structure rContext
  h      .l
  w      .l
  row_s  .l
  row_e  .l
  col_s  .l
  col_e  .l
  samps  .l
  tid    .l ; thread id
  isamps .d
  cx     .v3d
  cy     .v3d
  cam    .r3d
  Array c.v3d(0) ; pixels
EndStructure
; ...[ Rendering Synchronization ]............................................
Global rSem = CreateSemaphore()
; ...[ Rendering Thread ].....................................................
Procedure rThread( *ctx.rContext )
  
  ; Local Variables
  Protected ntemp.d,r.v3d
  
  Protected row_s.l = *ctx\row_s
  Protected row_e.l = *ctx\row_e
  Protected tid  .l = *ctx\tid
  SignalSemaphore(rSem)
  
  PrintN( "Starting Thread (tid:"+Str(tid)+") ["+Str(row_s)+":"+Str(row_e)+"]" )
  
  ; Loop over image rows
  Protected y.l
  For y = row_s To row_e
    PrintN( "Rendering ("+Str(*ctx\samps*4)+" spp) "+Str( 100.0*(y-row_s)/(row_e-row_s) )+"% (tid:"+Str(tid)+") ["+Str(row_s)+":"+Str(row_e)+"]        " )
    ; Loop over image columns
    Protected x.l
    For x = *ctx\col_s To *ctx\col_e
      
      ; Pixel Index
      Protected i.l = ( *ctx\h - y - 1 )* *ctx\w + x
      
      ; 2x2 subpixel rows
      Protected sy.l
      For sy = 0 To 1
        ; 2x2 subpixel cols
        Protected sx.l
        For sx = 0 To 1
          
          ; samples
          v3Null(r)
          Protected s.l
          For s=0 To *ctx\samps-1
            
            ; double r1 = 2*erand48(Xi), dx=r1<1 ? sqrt(r1)-1: 1-sqrt(2-r1);
            Protected r1.d,dx.d
            r1 = 2.0*rand()
            If r1 < 1.0
              dx = Sqr(r1)-1.0
            Else
              dx = 1.0-Sqr(2.0-r1)
            EndIf
            
            ; double r2 = 2*erand48(Xi), dy=r2<1 ? sqrt(r2)-1: 1-sqrt(2-r2);
            Protected r2.d,dy.d
            r2 = 2.0*rand()
            If r2 < 1.0
              dy = Sqr(r2)-1.0
            Else
              dy = 1.0-Sqr(2.0-r2)
            EndIf
            
            ; Vec d = cx*( ( (sx+.5 + dx)/2 + x)/w - .5) + cy*( ( (sy+.5 + dy)/2 + y)/h - .5) + cam.d;
            Protected d.v3d
            Protected ccx.d = ( ( sx + 0.5 + dx )/2.0 + x )/*ctx\w - 0.5
            Protected ccy.d = ( ( sy + 0.5 + dy )/2.0 + y )/*ctx\h - 0.5
            d\x = *ctx\cx\x*ccx + *ctx\cy\x*ccy + *ctx\cam\d\x
            d\y = *ctx\cx\y*ccx + *ctx\cy\y*ccy + *ctx\cam\d\y
            d\z = *ctx\cx\z*ccx + *ctx\cy\z*ccy + *ctx\cam\d\z
            
            ; r = r + radiance(Ray(cam.o+d*140,d.norm()),0,Xi)*(1./samps);
            ; Camera rays are pushed ^^^^^ forward to start in interior
            Protected rr.v3d,ray.r3d
            v3Null(rr)
            v3Copy(ray\d,d)
            v3Norm(ray\d,ntemp)
            v3ScaleInPlace(d,138.0)
            v3Add(ray\o,*ctx\cam\o,d)
            radiance(@rr,@ray,0,0)
            v3AddScaledInPlace(r,rr,*ctx\isamps)
            
          Next
          
          v3Clamp(r)
          v3AddScaledInPlace(*ctx\c(i),r,0.25)
          
        Next ; For sx = 0 To 1
      Next ; For sy = 0 To 1
      
    Next ; For x = *ctx\col_s To *ctx\col_e
  Next ; For y = row_s To row_e

EndProcedure
; ...[ Rendering Main ].......................................................
Procedure render( samps.l )
  
  Protected ntemp.d
  
  ; Context
  Protected ctx.rContext
  
  ; Topology
  ctx\w = 1024
  ctx\h =  768
  
  ; Samples
  ctx\samps  = samps
  ctx\isamps = 1.0/samps
  
  ; Cam pos
  ctx\cam\o\x = 50.0
  ctx\cam\o\y = 52.0
  ctx\cam\o\z = 295.6
  ; Cam dir
  ctx\cam\d\x = 0.0
  ctx\cam\d\y = -0.042612
  ctx\cam\d\z = -1.0
  v3Norm(ctx\cam\d,ntemp)
  
  ; Pixels
  ReDim ctx\c(ctx\w*ctx\h-1)
  
  ; Deltas
  Protected cx.v3d,cy.v3d,r.v3d
  ctx\cx\x = ctx\w*0.5135/ctx\h
  ctx\cx\y = 0.0
  ctx\cx\z = 0.0
  v3Cross(ctx\cy,ctx\cx,ctx\cam\d)
  v3Norm(ctx\cy,ntemp)
  v3Scale(ctx\cy,ctx\cy,0.5135)
  
  ; Columns
  ctx\col_s = 0
  ctx\col_e = ctx\w-1
  
  ; Prepare For n Thread(s)
  Dim aThread.l(0)
  Protected ri.l,rThreadCount = CountCPUs(#PB_System_ProcessCPUs)
  ReDim aThread(rThreadCount-1)
  
  ; Generate Threads
  Protected row_step      .l = ctx\h / rThreadCount
  Protected row_last_step .l = ctx\h - row_step*(rThreadCount-1)
  Protected row_s         .l = 0
  Protected row_e         .l = row_step
  Protected riStop        .l = ArraySize(aThread())
  For ri=0 To riStop
    ctx\row_s = row_s
    ctx\row_e = row_e - 1
    ctx\tid   = ri
    aThread(ri) = CreateThread(@rThread(),@ctx)
    WaitSemaphore(rSem)
    row_s + row_step
    If ri = riStop - 1
      row_e + row_last_step
    Else
      row_e + row_step
    EndIf
  Next
  
  ; Wait For Rendering End
  Protected Ta.d = ElapsedMilliseconds()
  For ri=0 To riStop
    WaitThread(aThread(ri))
  Next
  Protected Tb.d = ElapsedMilliseconds()
  
  ; Save Image
  CreateImage(0,ctx\w,ctx\h)
  StartDrawing(ImageOutput(0))
  Box(0,0,ctx\w,ctx\h,RGB($32,$32,$32))
  Protected x.l,y.l,i.l = 0
  For y=0 To ctx\h-1
    For x=0 To ctx\w-1
      Plot(x,y,RGB(toInt(ctx\c(i)\x),toInt(ctx\c(i)\y),toInt(ctx\c(i)\z)))
      i + 1
    Next
  Next
  DrawText( 0, 0, "Render Time : "+ FormatDate( "%hhh%iim%sss", AddDate(0,#PB_Date_Second,(Tb-Ta)/1000.0) ) )
  StopDrawing()
  PrintN( "Saving purept.bmp.." )
  SaveImage(0,"purept.bmp")
  PrintN( "Done." )
  Input()
  
EndProcedure

;render( 10 ) ; 40 spp (~21s)
render( 25 ) ; 100 spp (~59s)
;render( 50 ) ; 200 spp (~1m43s)
;render( 250 ) ; 1000 spp (~9m43s)
;render( 1250 ) ; 5k spp (~50m04s)
;render( 6250 ) ; 25k spp (not tested yet but several hours)

CloseConsole()
Cheers,
Guy.
Dernière modification par grabiller le sam. 24/août/2013 20:09, modifié 1 fois.
guy rabiller | radfac founder / ceo | raa.tel | raafal.org
Backup
Messages : 14526
Inscription : lun. 26/avr./2004 0:40

Re: purept (smallpt) path tracer

Message par Backup »

perso 1 minute 34 sur mon Portable i7
mais avec un resultat plein de "trous"

Image

quoiqu'il en soit c'est interressant :)
Avatar de l’utilisateur
grabiller
Messages : 103
Inscription : lun. 10/sept./2012 11:55
Localisation : France - 89220 Rogny-Les-Septs-Ecluses
Contact :

Re: purept (smallpt) path tracer

Message par grabiller »

Dobro a écrit :perso 1 minute 34 sur mon Portable i7
../..
Je pense que tu as du lancer le source par défaut en render(25). Essayes en render(50).
(j'aurai du le mettre par défaut en render(50), en fait :wink:)

Les 'trous' - ou le 'bruit' - c'est normal, c'est du au path tracing lui-même qui est une méthode qui 'converge' en fonction du nombre de samples.
Avec un nombre de samples infini, on a la solution idéale, mais avec un temps.. infini aussi :wink:
guy rabiller | radfac founder / ceo | raa.tel | raafal.org
Avatar de l’utilisateur
SPH
Messages : 4722
Inscription : mer. 09/nov./2005 9:53

Re: purept (smallpt) path tracer

Message par SPH »

Bravo, beau boulot ! :P
http://HexaScrabble.com/
!i!i!i!i!i!i!i!i!i!
!i!i!i!i!i!i!
!i!i!i!
//// Informations ////
Intel Core i7 4770 64 bits - GTX 650 Ti
Version de PB : 6.00 - 64 bits
Avatar de l’utilisateur
Ar-S
Messages : 9472
Inscription : dim. 09/oct./2005 16:51
Contact :

Re: purept (smallpt) path tracer

Message par Ar-S »

Test en cours sur mon vieux Core2duo e8200 2.6ghz en render(25) aussi.
=> 8m53 !
=> 17.27 en render(50)
Dernière modification par Ar-S le sam. 24/août/2013 23:30, modifié 1 fois.
~~~~Règles du forum ~~~~
⋅.˳˳.⋅ॱ˙˙ॱ⋅.˳Ar-S ˳.⋅ॱ˙˙ॱ⋅.˳˳.⋅
W11x64 PB 6.x
Section HORS SUJET : ICI
LDV MULTIMEDIA : Dépannage informatique & mes Logiciels PB
UPLOAD D'IMAGES : Uploader des images de vos logiciels
Avatar de l’utilisateur
blendman
Messages : 2017
Inscription : sam. 19/févr./2011 12:46

Re: purept (smallpt) path tracer

Message par blendman »

Sur mon ordi pourri, j'ai lancé en render(10) et ça a mis environ 6min40.

Question subsidiaire : serait-il possible d'ouvrir une scène type .obj (ou mieux .blend ^^), pour en faire un raytracer externe pour soft 3D (blender, max, lightwave..) ?
Ce serait topissime :D
Avatar de l’utilisateur
grabiller
Messages : 103
Inscription : lun. 10/sept./2012 11:55
Localisation : France - 89220 Rogny-Les-Septs-Ecluses
Contact :

Re: purept (smallpt) path tracer

Message par grabiller »

Voici une version améliorée.

La récursion a été retirée et le path tracing dépend maintenant réellement de la roulette russe pour la terminaison du rayon, ce qui permet d'obtenir des images quasi identiques avec celles de smallpt (la seule différence provenant du générateur de nombres aléatoires, donc du 'bruit').

Les 8 scènes de la page de Kevin Beason ont été ajoutées: Sky, NightSky, Island, Vista, Overlap, Wada, Wada2 and Forest. Vous pouvez choisir laquelle avec l'id de l'énumération: render( #SCENE_FOREST, 25 ) par exemple. La scène Overlap semble corrompue, j'ai du faire une erreur quelque part je vais regarder cela.

Les infos à l'image sont un peu meilleures ainsi que le retour de la console.

La scène Wada2 (par défaut si vous compilez directement le code) est particulièrement intéressante car elle ne contient que des surfaces spéculaires (réfléchissantes) ce qui permet d'avoir une image quasiment sans bruit avec seulement 40 samples (et même moins) par pixel. (40 spp en ~15s de mon côté).

N'oubliez pas d'activer l'option "Create threadsafe executable" avant compilation!

ImageImageImageImage
ImageImageImageImage

Code : Tout sélectionner

EnableExplicit

; ---[ 3d Vector ]------------------------------------------------------------
Structure v3d
  x.d
  y.d
  z.d
EndStructure
Macro v3Set( v_, x_, y_, z_ )
  v_\x = x_
  v_\y = y_
  v_\z = z_
EndMacro
Macro v3Null( v_ )
  v_\x = 0.0
  v_\y = 0.0
  v_\z = 0.0
EndMacro
Macro v3Clamp( v_ )
  v_\x = clamp(v_\x)
  v_\y = clamp(v_\y)
  v_\z = clamp(v_\z)
EndMacro
Macro v3Norm( v_, n_ )
  n_ = Sqr( v_\x*v_\x + v_\y*v_\y + v_\z*v_\z )
  If n_ <> 0.0
    n_ = 1.0/n_
    v_\x*n_
    v_\y*n_
    v_\z*n_
  EndIf
EndMacro
Macro v3Dot( v1_, v2_ )
  ( v1_\x * v2_\x + v1_\y * v2_\y + v1_\z * v2_\z )
EndMacro
Macro v3Copy( v_, v1_ )
  v_\x = v1_\x
  v_\y = v1_\y
  v_\z = v1_\z
EndMacro
Macro v3Add( v_, v1_, v2_ )
  v_\x = v1_\x + v2_\x
  v_\y = v1_\y + v2_\y
  v_\z = v1_\z + v2_\z
EndMacro
Macro v3AddScaledInPlace( v_, v1_, s1_ )
  v_\x + v1_\x * s1_
  v_\y + v1_\y * s1_
  v_\z + v1_\z * s1_
EndMacro
Macro v3Sub( v_, v1_, v2_ )
  v_\x = v1_\x - v2_\x
  v_\y = v1_\y - v2_\y
  v_\z = v1_\z - v2_\z
EndMacro
Macro v3SubInPlace( v_, v1_ )
  v_\x - v1_\x
  v_\y - v1_\y
  v_\z - v1_\z
EndMacro
Macro v3Mul( v_, v1_, v2_ )
  v_\x = v1_\x * v2_\x
  v_\y = v1_\y * v2_\y
  v_\z = v1_\z * v2_\z
EndMacro
Macro v3MulInPlace( v_, v1_ )
  v_\x * v1_\x
  v_\y * v1_\y
  v_\z * v1_\z
EndMacro
Macro v3Scale( v_, v1_, s1_ )
  v_\x = v1_\x * s1_
  v_\y = v1_\y * s1_
  v_\z = v1_\z * s1_
EndMacro
Macro v3ScaleInPlace( v_, s_ )
  v_\x * s_
  v_\y * s_
  v_\z * s_
EndMacro
Macro v3Cross( v_, v1_, v2_ )
  v_\x = v1_\y * v2_\z - v1_\z * v2_\y
  v_\y = v1_\z * v2_\x - v1_\x * v2_\z
  v_\z = v1_\x * v2_\y - v1_\y * v2_\x
EndMacro
Macro v3Neg( v_, v1_ )
  v_\x = -v1_\x
  v_\y = -v1_\y
  v_\z = -v1_\z
EndMacro

; ---[ 3d Ray ]---------------------------------------------------------------
Structure r3d
  o.v3d
  d.v3d
EndStructure

; ---[ Material ]-------------------------------------------------------------
Enumeration
  #DIFF
  #SPEC
  #REFR
EndEnumeration

; ---[ Sphere ]---------------------------------------------------------------
Structure s3d
  rad.d ; radius
  p.v3d ; position
  e.v3d ; emission
  c.v3d ; color
  m.i   ; material: #DIFF/#SPEC/#REFR
EndStructure
Procedure.d s3dIntersect( *sphere.s3d, *ray.r3d )
  ; returns distance, 0 if no hit
  ; Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0
  
  ; Vec op = p - r.o
  Protected op.v3d
  v3Sub( op, *sphere\p, *ray\o )
  
  ; double t, eps=1e-4, b=op.dot(r.d), det=b*b-op.dot(op)+rad*rad;
  Protected t.d
  Protected eps.d = 1.0e-4
  Protected b.d = v3Dot( op, *ray\d )
  Protected det.d = b*b - v3Dot( op, op ) + ( *sphere\rad * *sphere\rad )
  
  ; if (det<0) return 0; else det=sqrt(det);
  If det < 0.0
    ProcedureReturn 0.0
  Else
    det = Sqr(det)
  EndIf
  
  ; return (t=b-det)>eps ? t : ((t=b+det)>eps ? t : 0);
  t = b - det
  If t <= eps
    t = b + det
    If t <= eps
      ProcedureReturn 0.0
    EndIf
  EndIf
  
  ProcedureReturn t

EndProcedure

; ---[ Scenes ]---------------------------------------------------------------
Macro Sphere( idx_, rad_, posx_, posy_, posz_, emir_, emig_, emib_, colr_, colg_, colb_, mat_ )
  scene(idx_)\rad = rad_
  scene(idx_)\p\x = posx_
  scene(idx_)\p\y = posy_
  scene(idx_)\p\z = posz_
  scene(idx_)\e\x = emir_
  scene(idx_)\e\y = emig_
  scene(idx_)\e\z = emib_
  scene(idx_)\c\x = colr_
  scene(idx_)\c\y = colg_
  scene(idx_)\c\z = colb_
  scene(idx_)\m = mat_
EndMacro
Global Dim scene.s3d(0)
Enumeration
  #SCENE_CORNELL
  #SCENE_SKY
  #SCENE_NIGHTSKY
  #SCENE_ISLAND
  #SCENE_VISTA
  #SCENE_OVERLAP
  #SCENE_WADA
  #SCENE_WADA2
  #SCENE_FOREST
EndEnumeration
Procedure SetupScene( scene_id.l )
  
  Select scene_id
    Case #SCENE_CORNELL
      ReDim scene(8)
      Sphere( 0, 1.0e5,   1.0e5+1.0,  40.8,       81.6,         0.0, 0.0, 0.0,   0.75,  0.25,  0.25,  #DIFF ) ; Left
      Sphere( 1, 1.0e5,  -1.0e5+99.0, 40.8,       81.6,         0.0, 0.0, 0.0,   0.25,  0.25,  0.75,  #DIFF ) ; Rght
      Sphere( 2, 1.0e5,   50.0,       40.8,       1.0e5,        0.0, 0.0, 0.0,   0.75,  0.75,  0.75,  #DIFF ) ; Back
      Sphere( 3, 1.0e5,   50.0,       40.8,      -1.0e5+170.0,  0.0, 0.0, 0.0,   0.0,   0.0,   0.0,   #DIFF ) ; Frnt
      Sphere( 4, 1.0e5,   50.0,       1.0e5,      81.6,         0.0, 0.0, 0.0,   0.75,  0.75,  0.75,  #DIFF ) ; Botm
      Sphere( 5, 1.0e5,   50.0,      -1.0e5+81.6, 81.6,         0.0, 0.0, 0.0,   0.75,  0.75,  0.75,  #DIFF ) ; Top
      Sphere( 6, 16.5,    27.0,       16.5,       47.0,         0.0, 0.0, 0.0,   0.999, 0.999, 0.999, #SPEC ) ; Mirr
      Sphere( 7, 16.5,    73.0,       16.5,       78.0,         0.0, 0.0, 0.0,   0.999, 0.999, 0.999, #REFR ) ; Glas
      Sphere( 8, 600.0,   50.0,       681.6-0.27, 81.6,        12.0,12.0,12.0,   0.0,   0.0,   0.0,   #DIFF ) ; Lite
    Case #SCENE_SKY
      ReDim scene(8)
      Sphere( 0, 1600,    3000, 0,              6000,  37.4399,       33.6959,      29.9519,        0.0,    0.0,    0.0,    #DIFF ) ; sun
      Sphere( 1, 1560,    3500, 0,              7000,  149.7599,      74.8799,      7.4879,         0.0,    0.0,    0.0,    #DIFF ) ; horizon sun2
      Sphere( 2, 10000,   50,   40.8,          -1060,  0.00030644159, 0.0096070944, 0.13883156639,  0.175,  0.175,  0.25,   #DIFF ) ; sky
      Sphere( 3, 100000,  50,  -100000,         0,     0.0,           0.0,          0.0,            0.3,    0.3,    0.3,    #DIFF ) ; grnd
      Sphere( 4, 110000,  50,  -110048.5,       0,     3.6,           2.0,          0.2,            0.0,    0.0,    0.0,    #DIFF ) ; horizon brightener
      Sphere( 5, 4e4,     50,  -4e4-30,        -3000,  0.0,           0.0,          0.0,            0.2,    0.2,    0.2,    #DIFF ) ; mountains
      Sphere( 6, 26.5,    22,   26.5,           42,    0.0,           0.0,          0.0,            0.596,  0.596,  0.596,  #SPEC ) ; white Mirr
      Sphere( 7, 13,      75,   13,             82,    0.0,           0.0,          0.0,            0.9216, 0.9216, 0.9216, #REFR ) ; Glas
      Sphere( 8, 22,      87,   22,             24,    0.0,           0.0,          0.0,            0.4176, 0.4176, 0.4176, #REFR ) ; Glas2
    Case #SCENE_NIGHTSKY
      ReDim scene(11)
      Sphere(  0, 2.5e3,  8200.0,   9200.0,  -20000.0,  0.8e2,   0.8e2,   0.8e2,    0.0,       0.0,       0.0,       #DIFF ) ; moon
      Sphere(  1, 2.5e4,  50.0,     0.0,      0.0,      0.00114, 0.00133, 0.00212,  0.000648,  0.001152,  0.003,     #DIFF ) ; sky
      Sphere(  2, 5e0,    -2000.0,  1600.0,  -10000.0,  100.0,   84.3,    69.8,     0.0,       0.0,       0.0,       #DIFF ) ; star
      Sphere(  3, 5e0,    0.0,      1800.0,  -10000.0,  100.0,   85.1,    71.0,     0.0,       0.0,       0.0,       #DIFF ) ; star
      Sphere(  4, 5e0,    3000.0,   1500.0,  -10000.0,  67.1,    78.0,    100.0,    0.0,       0.0,       0.0,       #DIFF ) ; star
      Sphere(  5, 3.5e4,  600.0,   -3.5e4+1,  300.0,    0.0,     0.0,     0.0,      0.006,     0.008,     0.01,      #REFR ) ; pool
      Sphere(  6, 5e4,    -500.0,  -5e4+0,    0.0,      0.0,     0.0,     0.0,      0.35,      0.35,      0.35,      #DIFF ) ; hill
      Sphere(  7, 16.5,   27.0,     0.0,      47.0,     0.0,     0.0,     0.0,      0.33,      0.33,      0.33,      #DIFF ) ; hut
      Sphere(  8, 7,      38.3137,  0.0,      58.3137,  0.0,     0.0,     0.0,      0.33,      0.33,      0.33,      #DIFF ) ; door
      Sphere(  9, 500,    -1.0e3, -300.0,    -3.0e3,    0.0,     0.0,     0.0,      0.351,     0.351,     0.351,     #DIFF ) ; mnt
      Sphere( 10, 830,    0.0,    -500.0,    -3.0e3,    0.0,     0.0,     0.0,      0.354,     0.354,     0.354,     #DIFF ) ; mnt
      Sphere( 11, 490,    1.0e3,  -300.0,    -3.0e3,    0.0,     0.0,     0.0,      0.352,     0.352,     0.352,     #DIFF ) ; mnt
    Case #SCENE_ISLAND
      ReDim scene(6)
      Sphere( 0, 160.0,    50.0,  580.0, -1360.0,  2e2,    2e2,   2e2,    0.0,     0.0,     0.0,   #DIFF ) ; sun
      Sphere( 1, 800.0,    50.0, -900.0, -9980.0,  2e1,    2e1,   2e1,    0.0,     0.0,     0.0,   #DIFF ) ; horizon
      Sphere( 2, 10000.0,  50.0, -20.0,  -1060.0,  0.0627, 0.188, 0.569,  0.4,     0.4,     0.4,   #DIFF ) ; sky
      Sphere( 3, 800.0,    50.0, -740.0, -1060.0,  0.0,    0.0,   0.0,    0.10956, 0.89441, 0.996, #REFR ) ; water
      Sphere( 4, 790.0,    50.0, -740.0, -1060.0,  0.0,    0.0,   0.0,    0.24,    0.18,    0.024, #DIFF ) ; earth
      Sphere( 5, 325.0,    50.0, -275.0, -910.0,   0.0,    0.0,   0.0,    0.32,    0.24,    0.032, #DIFF ) ; island
      Sphere( 6, 275.0,    50.0, -225.0, -893.0,   0.0,    0.0,   0.0,    0.015,   0.225,   0.015, #DIFF ) ; grass
    Case #SCENE_VISTA
      ReDim scene(23)
      Sphere(  0, 8000,   50.0,  -8020.0, -1760.0,   0.5,    0.2,    0.05,   0.0,   0.0,   0.0,   #DIFF ) ; sun
      Sphere(  1, 1e4,    50.0,  -20.0,   -860.0,    0.1893, 0.2259, 0.3,    0.5,   0.5,   0.5,   #DIFF ) ; sky
      Sphere(  2, 150,   -300.0, -20.0,   -960.0,    0.0,    0.0,    0.0,    0.3,   0.3,   0.3,   #DIFF ) ; mnt
      Sphere(  3, 200,   -160.0, -20.0,   -960.0,    0.0,    0.0,    0.0,    0.3,   0.3,   0.3,   #DIFF ) ; mnt
      Sphere(  4, 145,   -160.0,  65.0,   -960.0,    0.0,    0.0,    0.0,    0.8,   0.8,   0.8,   #DIFF ) ; snow
      Sphere(  5, 150,    0.0,   -20.0,   -960.0,    0.0,    0.0,    0.0,    0.3,   0.3,   0.3,   #DIFF ) ; mnt
      Sphere(  6, 150,    150.0, -20.0,   -960.0,    0.0,    0.0,    0.0,    0.3,   0.3,   0.3,   #DIFF ) ; mnt
      Sphere(  7, 125,    300.0, -20.0,   -960.0,    0.0,    0.0,    0.0,    0.3,   0.3,   0.3,   #DIFF ) ; mnt
      Sphere(  8, 150,    425.0, -20.0,   -960.0,    0.0,    0.0,    0.0,    0.3,   0.3,   0.3,   #DIFF ) ; mnt
      Sphere(  9, 2500,   50.0,  -2420.0, -1360.0,   0.0,    0.0,    0.0,    0.1,   0.1,   0.1,   #DIFF ) ; mnt base
      Sphere( 10, 8000,   50.0,  -8020.0, -660.0,    0.0,    0.0,    0.0,    0.2,   0.2,   1.0,   #REFR ) ; water
      Sphere( 11, 8000,   50.0,  -8020.0,  240.0,    0.0,    0.0,    0.0,    0.0,   0.3,   0.0,   #DIFF ) ; grass
      Sphere( 12, 8,     -25.0,  -25.0,   -10.0,     0.0,    0.0,    0.0,    0.0,   0.3,   0.0,   #DIFF ) ; bush
      Sphere( 13, 30,     50.0,   3.0,    -35.0,     0.0,    0.0,    0.0,    0.996, 0.996, 0.996, #REFR ) ; ball
      Sphere( 14, 30,     250.0,  260.0,  -1260.0,   0.0,    0.0,    0.0,    0.8,   0.8,   0.8,   #DIFF ) ; clouds
      Sphere( 15, 37,     287.0,  260.0,  -1260.0,   0.0,    0.0,    0.0,    0.8,   0.8,   0.8,   #DIFF ) ; clouds
      Sphere( 16, 28,     317.0,  260.0,  -1260.0,   0.0,    0.0,    0.0,    0.8,   0.8,   0.8,   #DIFF ) ; clouds
      Sphere( 17, 40,     200.0,  260.0,  -1860.0,   0.0,    0.0,    0.0,    0.8,   0.8,   0.8,   #DIFF ) ; clouds
      Sphere( 18, 37,     237.0,  260.0,  -1860.0,   0.0,    0.0,    0.0,    0.8,   0.8,   0.8,   #DIFF ) ; clouds
      Sphere( 19, 40,     650.0,  260.0,  -1960.0,   0.0,    0.0,    0.0,    0.8,   0.8,   0.8,   #DIFF ) ; clouds
      Sphere( 20, 37,     687.0,  260.0,  -1960.0,   0.0,    0.0,    0.0,    0.8,   0.8,   0.8,   #DIFF ) ; clouds
      Sphere( 21, 37,    -750.0,  260.0,  -2260.0,   0.0,    0.0,    0.0,    0.8,   0.8,   0.8,   #DIFF ) ; clouds
      Sphere( 22, 37,     50.0,   260.0,  -2460.0,   0.0,    0.0,    0.0,    0.8,   0.8,   0.8,   #DIFF ) ; clouds
      Sphere( 23, 37,     587.0,  260.0,  -2660.0,   0.0,    0.0,    0.0,    0.8,   0.8,   0.8,   #DIFF ) ; clouds
    Case #SCENE_OVERLAP
      ReDim scene(2)
      Sphere( 0, 150.0,  125.0, 28.0, 62.0,  0.0, 0.0, 0.0,  0.93, 0.837, 0.744, #REFR )
      Sphere( 1, 28.0,   55.0, -28.0, 62.0,  1e1, 1e1, 1e1,  0.0,  0.0,   0.0,   #DIFF )
      Sphere( 2, 300.0,  50.0,  28.0, 62.0,  0.0, 0.0, 0.0,  0.93, 0.93,  0.93,  #SPEC )
    Case #SCENE_WADA
      ReDim scene(6)
      Sphere( 0, 1.0e5,  50.0,  100.0,       0.0,       3.0, 3.0, 3.0,  0.0,     0.0,     0.0,     #DIFF ) ; sky
      Sphere( 1, 1.0e5,  50.0, -100129.2820, 0.0,       0.0, 0.0, 0.0,  0.1,     0.1,     0.1,     #DIFF ) ; grnd
      Sphere( 2, 60.0,   110.0, 75.441,      62.0,      0.0, 0.0, 0.0,  0.999,   0.2997,  0.2997,  #SPEC ) ; red
      Sphere( 3, 60.0,  -10.0,  75.441,      62.0,      0.0, 0.0, 0.0,  0.2997,  0.999,   0.2997,  #SPEC ) ; grn
      Sphere( 4, 60.0,   50.0, -28.4820,     62.0,      0.0, 0.0, 0.0,  0.2997,  0.2997,  0.999,   #SPEC ) ; blue
      Sphere( 5, 60.0,   50.0,  40.8,       -7.2820,    0.0, 0.0, 0.0,  0.52947, 0.52947, 0.52947, #SPEC ) ; back
      Sphere( 6, 60.0,   50.0,  40.8,        131.2820,  0.0, 0.0, 0.0,  0.999,   0.999,   0.999,   #REFR ) ; front
    Case #SCENE_WADA2
      ReDim scene(4)
      Sphere( 0, 120.0,     170.0, 97.2822,  62.0,       0.0165, 0.0367, 0.0569,   0.996,  0.996,  0.996,  #SPEC ) ; red
      Sphere( 1, 120.0,    -70.0,  97.2822,  62.0,       0.0165, 0.0367, 0.0569,   0.996,  0.996,  0.996,  #SPEC ) ; grn
      Sphere( 2, 120.0,     50.0, -110.5641, 62.0,       0.0165, 0.0367, 0.0569,   0.996,  0.996,  0.996,  #SPEC ) ; blue
      Sphere( 3, 120.0,     50.0,  28.0,    -133.9592,   0.0,    0.0,    0.0,      0.996,  0.996,  0.996,  #SPEC ) ; back
      Sphere( 4, 718.517,   50.0,  28.0,    -3.3197,     0.0,    0.0,    0.0,      0.5,    0.5,    0.5,    #SPEC ) ; front
    Case #SCENE_FOREST
      ReDim scene(12)
      Sphere(  0, 1e5,   50.0,    1e5+130, 0.0,        1.3,    1.3,    1.3,   0.0,    0.0,   0.0,    #DIFF ) ; lite
      Sphere(  1, 1e2,   50.0,   -1e2+2,   47.0,       0.0,    0.0,    0.0,   0.7,    0.7,   0.7,    #DIFF ) ; grnd
      Sphere(  2, 1e4,  -7610.4444, -30.0, 6727.8761,  0.0,    0.0,    0.0,   0.99,   0.99,  0.99,   #SPEC ) ; mirr L
      Sphere(  3, 1e4,   7710.4444, -30.0, 6727.8761,  0.0,    0.0,    0.0,   0.99,   0.99,  0.99,   #SPEC ) ; mirr R
      Sphere(  4, 1e4,  -4950.0, -30.0,   -8710.254,   0.0,    0.0,    0.0,   0.99,   0.99,  0.99,   #SPEC ) ; mirr FL
      Sphere(  5, 1e4,   5050.0, -30.0,   -8710.254,   0.0,    0.0,    0.0,   0.99,   0.99,  0.99,   #SPEC ) ; mirr
      Sphere(  6, 4,     50.0,    3.6,     47.0,       0.0,    0.0,    0.0,   0.13,   0.066, 0.033,  #DIFF ) ; "tree"
      Sphere(  7, 16,    50.0,    21.6,    47.0,       0.0,    0.0,    0.0,   0.0588, 0.361, 0.0941, #DIFF ) ; "tree"
      Sphere(  8, 11,    50.0,    37.8,    47.0,       0.0,    0.0,    0.0,   0.0588, 0.361, 0.0941, #DIFF ) ; "tree"
      Sphere(  9, 7,     50.0,    48.6,    47.0,       0.0,    0.0,    0.0,   0.0588, 0.361, 0.0941, #DIFF ) ; "tree"
      Sphere( 10, 15.5,  50.0,    23.4,    47.0,       0.0,    0.0,    0.0,   0.7,    0.7,   0.7,    #DIFF ) ; "tree"
      Sphere( 11, 10.5,  50.0,    39.6,    47.0,       0.0,    0.0,    0.0,   0.7,    0.7,   0.7,    #DIFF ) ; "tree"
      Sphere( 12, 6.5,   50.0,    50.4,    47.0,       0.0,    0.0,    0.0,   0.7,    0.7,   0.7,    #DIFF ) ; "tree"
  EndSelect
        
EndProcedure

; ---[ Helpers ]--------------------------------------------------------------
Procedure.d clamp( x.d )
  If x < 0.0
    ProcedureReturn 0.0
  ElseIf x > 1.0
    ProcedureReturn 1.0
  EndIf
  ProcedureReturn x
EndProcedure
Procedure.l toInt( x.d )
  ; linear to ~sRGB (not accurate), yet the spheres colors are not
  ; linearized prior rendering.
  ProcedureReturn Int( Pow(clamp(x),1.0/2.2)*255+0.5 )
EndProcedure
Procedure.b intersect( *ray.r3d, *t.Double, *id.Long )
  Protected d.d, inf.d = 1e20
  Protected i.i = ArraySize(scene())
  *t\d = inf
  While i >= 0
    d = s3dIntersect( scene(i), *ray )
    If d And d<*t\d
      *t\d  = d
      *id\l = i
    EndIf
    i - 1
  Wend
  ProcedureReturn Bool( *t\d < inf )
EndProcedure
Procedure.d rand()
  ProcedureReturn Random(16384,0)/16384.0
EndProcedure

; ---[ Radiance ]-------------------------------------------------------------
Procedure radiance( *ret.v3d, *r.r3d, depth.l )
  
  Protected ntemp.d
  Protected tt_.d
  Protected r2_.r3d
  
  Protected t .d = 0.0 ; distance to intersection
  Protected id.l = 0   ; id of intersected object

  ; if miss, return black
  If Not intersect( *r, @t, @id )
    v3Null(*ret)
    ProcedureReturn
  EndIf
  
  ; the hit object
  Protected *obj.s3d = scene(id)
  
  ; x = r.o + r.d*t
  Protected x.v3d
  x\x = *r\o\x + *r\d\x*t
  x\y = *r\o\y + *r\d\y*t
  x\z = *r\o\z + *r\d\z*t
  
  ; n=(x-obj.p).norm()
  Protected n.v3d
  v3Sub( n, x, *obj\p )
  v3Norm( n, ntemp )
  
  ; nl=n.dot(r.d)<0?n:n*-1
  Protected nl.v3d
  If v3Dot( n, *r\d ) < 0.0
    v3Copy( nl, n )
  Else
    v3Neg( nl, n )
  EndIf
  
  ; f=obj.c
  Protected f.v3d
  v3Copy( f, *obj\c )
  
  ; max refl
  Protected p.d
  If ( f\x > f\y ) And ( f\x > f\z )
    p = f\x
  ElseIf f\y > f\z
    p = f\y
  Else
    p = f\z
  EndIf

  ; Russian Roulette
  depth + 1
  If depth > 5
    ;If rand() < p ; erand48(Xi) < p
      ;f\x/p
      ;f\y/p
      ;f\z/p
    ;Else
      v3Copy( *ret, *obj\e )
      ProcedureReturn
    ;EndIf
  EndIf
  
  Protected rr.v3d
  Protected rt.r3d
  Protected rd.v3d
  
  ; Ideal DIFFUSE reflection
  If *obj\m = #DIFF

    Protected r1  .d = 2.0*#PI*rand()
    Protected r2  .d = rand()
    Protected r2s .d = Sqr(r2)
    Protected w.v3d
    v3Copy( w, nl )
    Protected ut.v3d,u.v3d
    If Abs( w\x ) > 0.1
      v3Set(ut,0.0,1.0,0.0)
    Else
      v3Set(ut,1.0,0.0,0.0)
    EndIf
    v3Cross(u,ut,w)
    v3Norm(u,ntemp)
    Protected v.v3d
    v3Cross(v,w,u)
    
    Protected d.v3d
    d\x = u\x*Cos(r1)*r2s + v\x*Sin(r1)*r2s + w\x*Sqr(1.0-r2)
    d\y = u\y*Cos(r1)*r2s + v\y*Sin(r1)*r2s + w\y*Sqr(1.0-r2)
    d\z = u\z*Cos(r1)*r2s + v\z*Sin(r1)*r2s + w\z*Sqr(1.0-r2)
    v3Norm(d,ntemp)
    
    ; return obj.e + f.mult( radiance( Ray(x,d),depth,Xi) );
    v3Copy( rt\o, x )
    v3Copy( rt\d, d )
    radiance( @rr, @rt, depth )
    v3MulInPlace( f, rr )    
    v3Add( *ret, *obj\e, f )
    ProcedureReturn
    
  ; Ideal SPECULAR reflection
  ElseIf *obj\m = #SPEC

    ; return obj.e + f.mult(radiance( Ray( x, r.d-n*2*n.dot(r.d) ), depth, Xi ) );
    v3Copy( rt\o, x )
    v3Copy(rd,n)
    tt_ = 2.0*v3Dot(n,*r\d)
    v3ScaleInPlace( rd, tt_ )
    v3Sub( rt\d, *r\d, rd )
    radiance( @rr, @rt, depth )
    v3MulInPlace( f, rr )
    v3Add( *ret, *obj\e, f )
    ProcedureReturn
    
  EndIf
  
  ; Ideal dielectric REFRACTION
  ; Ray reflRay( x, r.d - n*2*n.dot(r.d) );
  Protected reflRay.r3d
  v3Copy( reflRay\o, x )
  v3Copy(rd,n)
  tt_ = 2.0*v3Dot(n,*r\d)
  v3ScaleInPlace( rd, tt_ )
  v3Sub( reflRay\d, *r\d, rd )
  ; Ray from outside going in? 
  ; bool into = n.dot(nl)>0;
  Protected into.b = Bool( v3Dot(n,nl) > 0.0 )
  ; double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t;
  Protected nc.d=1.0,nt.d=1.5,nnt.d,ddn.d,cos2t.d
  If into
    nnt = nc/nt
  Else
    nnt = nt/nc
  EndIf
  ddn = v3Dot(*r\d,nl)
  cos2t = 1.0 - nnt*nnt*( 1.0 - ddn*ddn )
  If cos2t < 0 ; Total internal reflection
    ; return obj.e + f.mult(radiance(reflRay,depth,Xi));
    radiance( @rr, @reflRay, depth )
    v3MulInPlace( f, rr )
    v3Add( *ret, *obj\e, f )
    ProcedureReturn
  EndIf
  ; Vec tdir = ( r.d*nnt - n*(  (into?1:-1)*(ddn*nnt+sqrt(cos2t))  ) ).norm();
  Protected intod.d
  If into
    intod = 1.0
  Else
    intod = -1.0
  EndIf
  Protected t_.d,n_.v3d,r_.v3d
  t_ = intod*(ddn*nnt+Sqr(cos2t))
  v3Copy(n_,n)
  v3ScaleInPlace( n_, t_ )
  v3Copy(r_,*r\d)
  v3ScaleInPlace( r_, nnt )
  Protected tdir.v3d
  v3Sub(tdir,r_,n_)
  v3Norm(tdir,ntemp)
  ; double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c = 1-(into?-ddn:tdir.dot(n));
  Protected a.d=nt-nc,b.d=nt+nc,R0.d=a*a/(b*b),c.d
  If into
    c = 1.0 + ddn
  Else
    c = 1.0 - v3Dot(tdir,n)
  EndIf
  ; double Re=R0+(1-R0)*c*c*c*c*c,Tr=1-Re,P=.25+.5*Re,RP=Re/P,TP=Tr/(1-P);
  Protected Re.d=R0+(1.0-R0)*c*c*c*c*c
  Protected Tr.d=1.0-Re
  Protected P_.d=0.25+0.5*Re
  Protected RP.d=Re/P_
  Protected TP.d=Tr/(1.0-P_)
  ; return obj.e + f.mult(depth>2 ? ( erand48(Xi)<P ?   // Russian roulette
  ; radiance(reflRay,depth,Xi)*RP : radiance(Ray(x,tdir),depth,Xi)*TP ) :
  ; radiance(reflRay,depth,Xi)*Re+radiance(Ray(x,tdir),depth,Xi)*Tr);
  If depth > 2
    If rand() < P_
      ; radiance(reflRay,depth,Xi)*RP
      radiance( @rr, @reflRay, depth )
      v3ScaleInPlace( rr, RP )
    Else
      ; radiance(Ray(x,tdir),depth,Xi)*TP )
      v3Copy(r2_\o,x)
      v3Copy(r2_\d,tdir)
      radiance( @rr, @r2_, depth )
      v3ScaleInPlace( rr, TP )
    EndIf
  Else
    Protected rr1.v3d
    Protected rr2.v3d
    ; radiance(reflRay,depth,Xi)*Re + radiance(Ray(x,tdir),depth,Xi)*Tr)
    radiance( @rr1, @reflRay, depth )
    v3ScaleInPlace( rr1, Re )
    v3Copy(r2_\o,x)
    v3Copy(r2_\d,tdir)
    radiance( @rr2, @r2_, depth )
    v3ScaleInPlace( rr1, Tr )
    v3Add(rr,rr1,rr2)
  EndIf
  v3Mul( f, f, rr )
  v3Add( *ret, *obj\e, f )
  ProcedureReturn
  
EndProcedure
Procedure radiancef( *ret.v3d, *r_.r3d, depth_.l )
  
  Protected ntemp.d
  Protected tt_.d
  Protected r2_.r3d
  
  Protected t .d = 0.0 ; distance to intersection
  Protected id.l = 0   ; id of intersected object
  
  Protected r.r3d
  v3Copy(r\o,*r_\o)
  v3Copy(r\d,*r_\d)
  
  Protected depth.l = depth_
  
  ; L0 = Le0 + f0*(L1)
  ;    = Le0 + f0*(Le1 + f1*L2)
  ;    = Le0 + f0*(Le1 + f1*(Le2 + f2*(L3))
  ;    = Le0 + f0*(Le1 + f1*(Le2 + f2*(Le3 + f3*(L4)))
  ;    = ...
  ;    = Le0 + f0*Le1 + f0*f1*Le2 + f0*f1*f2*Le3 + f0*f1*f2*f3*Le4 + ...
  ; 
  ; So:
  ; F = 1
  ; While (1){
  ;   L += F*Lei
  ;   F *= fi
  ; }
  
  Protected cl.v3d ; accumulated color
  v3Null(cl)
  Protected cf.v3d ; accumulated reflectance
  v3Set(cf,1.0,1.0,1.0)
  
  While 1
    
    ; if miss, return black
    If Not intersect( @r, @t, @id )
      v3Copy(*ret,cl)
      ProcedureReturn
    EndIf
    
    ; the hit object
    Protected *obj.s3d = scene(id)
    
    ; x = r.o + r.d*t
    Protected x.v3d
    x\x = r\o\x + r\d\x*t
    x\y = r\o\y + r\d\y*t
    x\z = r\o\z + r\d\z*t
    
    ; n=(x-obj.p).norm()
    Protected n.v3d
    v3Sub( n, x, *obj\p )
    v3Norm( n, ntemp )
    
    ; nl=n.dot(r.d)<0?n:n*-1
    Protected nl.v3d
    If v3Dot( n, r\d ) < 0.0
      v3Copy( nl, n )
    Else
      v3Neg( nl, n )
    EndIf
    
    ; f=obj.c
    Protected f.v3d
    v3Copy( f, *obj\c )
    
    ; max refl
    Protected p.d
    If ( f\x > f\y ) And ( f\x > f\z )
      p = f\x
    ElseIf f\y > f\z
      p = f\y
    Else
      p = f\z
    EndIf
    
    ; cl = cl + cf.mult(obj.e);
    cl\x + cf\x * *obj\e\x
    cl\y + cf\y * *obj\e\y
    cl\z + cf\z * *obj\e\z
    
    ; Russian Roulette
    depth + 1
    If depth > 5
      If rand() < p ; erand48(Xi) < p
        f\x/p
        f\y/p
        f\z/p
      Else
        v3Copy(*ret,cl)
        ProcedureReturn
      EndIf
    EndIf
    
    ; cf = cf.mult(f);
    v3MulInPlace(cf,f)
    
    ;Protected rr.v3d
    ;Protected rt.r3d
    Protected rd.v3d
    
    ; Ideal DIFFUSE reflection
    If *obj\m = #DIFF
  
      Protected r1  .d = 2.0*#PI*rand()
      Protected r2  .d = rand()
      Protected r2s .d = Sqr(r2)
      Protected w.v3d
      v3Copy( w, nl )
      Protected ut.v3d,u.v3d
      If Abs( w\x ) > 0.1
        v3Set(ut,0.0,1.0,0.0)
      Else
        v3Set(ut,1.0,0.0,0.0)
      EndIf
      v3Cross(u,ut,w)
      v3Norm(u,ntemp)
      Protected v.v3d
      v3Cross(v,w,u)
      
      Protected d.v3d
      d\x = u\x*Cos(r1)*r2s + v\x*Sin(r1)*r2s + w\x*Sqr(1.0-r2)
      d\y = u\y*Cos(r1)*r2s + v\y*Sin(r1)*r2s + w\y*Sqr(1.0-r2)
      d\z = u\z*Cos(r1)*r2s + v\z*Sin(r1)*r2s + w\z*Sqr(1.0-r2)
      v3Norm(d,ntemp)
      
      ; return obj.e + f.mult( radiance( Ray(x,d),depth,Xi) );
      ; -> r = Ray(x,d);
      v3Copy(r\o,x)
      v3Copy(r\d,d)
      Continue
      
    ; Ideal SPECULAR reflection
    ElseIf *obj\m = #SPEC
  
      ; return obj.e + f.mult(radiance( Ray( x, r.d-n*2*n.dot(r.d) ), depth, Xi ) );
      ; -> r = Ray(x, r.d - n* 2*n.dot(r.d) );
      v3Copy(r\o,x)
      v3Copy(rd,n)
      tt_ = 2.0*v3Dot(n,r\d)
      v3ScaleInPlace( rd, tt_ )
      v3SubInPlace( r\d, rd )
      Continue
      
    EndIf
    
    ; Ideal dielectric REFRACTION
    ; Ray reflRay( x, r.d - n*2*n.dot(r.d) );
    Protected reflRay.r3d
    v3Copy( reflRay\o, x )
    v3Copy(rd,n)
    tt_ = 2.0*v3Dot(n,r\d)
    v3ScaleInPlace( rd, tt_ )
    v3Sub( reflRay\d, r\d, rd )
    ; Ray from outside going in? 
    ; bool into = n.dot(nl)>0;
    Protected into.b = Bool( v3Dot(n,nl) > 0.0 )
    ; double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t;
    Protected nc.d=1.0,nt.d=1.5,nnt.d,ddn.d,cos2t.d
    If into
      nnt = nc/nt
    Else
      nnt = nt/nc
    EndIf
    ddn = v3Dot(r\d,nl)
    cos2t = 1.0 - nnt*nnt*( 1.0 - ddn*ddn )
    If cos2t < 0 ; Total internal reflection
      ; return obj.e + f.mult(radiance(reflRay,depth,Xi));
      ; -> r = reflRay;
      v3Copy(r\o,reflRay\o)
      v3Copy(r\d,reflRay\d)
      Continue
    EndIf
    ; Vec tdir = ( r.d*nnt - n*(  (into?1:-1)*(ddn*nnt+sqrt(cos2t))  ) ).norm();
    Protected intod.d
    If into
      intod = 1.0
    Else
      intod = -1.0
    EndIf
    Protected t_.d,n_.v3d,r_.v3d
    t_ = intod*(ddn*nnt+Sqr(cos2t))
    v3Copy(n_,n)
    v3ScaleInPlace( n_, t_ )
    v3Copy(r_,r\d)
    v3ScaleInPlace( r_, nnt )
    Protected tdir.v3d
    v3Sub(tdir,r_,n_)
    v3Norm(tdir,ntemp)
    ; double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c = 1-(into?-ddn:tdir.dot(n));
    Protected a.d=nt-nc,b.d=nt+nc,R0.d=a*a/(b*b),c.d
    If into
      c = 1.0 + ddn
    Else
      c = 1.0 - v3Dot(tdir,n)
    EndIf
    ; double Re=R0+(1-R0)*c*c*c*c*c,Tr=1-Re,P=.25+.5*Re,RP=Re/P,TP=Tr/(1-P);
    Protected Re.d=R0+(1.0-R0)*c*c*c*c*c
    Protected Tr.d=1.0-Re
    Protected P_.d=0.25+0.5*Re
    Protected RP.d=Re/P_
    Protected TP.d=Tr/(1.0-P_)
    
    ; return obj.e + f.mult(depth>2 ? ( erand48(Xi)<P ?   // Russian roulette
    ; radiance(reflRay,depth,Xi)*RP : radiance(Ray(x,tdir),depth,Xi)*TP ) :
    ; radiance(reflRay,depth,Xi)*Re+radiance(Ray(x,tdir),depth,Xi)*Tr);
    ; ->
    ; If (erand48(Xi)<P){
    ;   cf = cf*RP;
    ;   r = reflRay;
    ; } Else {
    ;   cf = cf*TP;
    ;   r = Ray(x,tdir);
    ; }
    ; Continue;
    If rand() < P_
      v3ScaleInPlace(cf,RP)
      v3Copy(r\o,reflRay\o)
      v3Copy(r\d,reflRay\d)
    Else
      v3ScaleInPlace(cf,TP)
      v3Copy(r\o,x)
      v3Copy(r\d,tdir)
    EndIf
    Continue

  Wend
  
EndProcedure
; ---[ Render ]---------------------------------------------------------------
; ...[ Rendering Context ]....................................................
Structure rContext
  h      .l
  w      .l
  row_s  .l
  row_e  .l
  col_s  .l
  col_e  .l
  samps  .l
  tid    .l ; thread id
  isamps .d
  cx     .v3d
  cy     .v3d
  cam    .r3d
  Array c.v3d(0) ; pixels
EndStructure
; ...[ Rendering Synchronization ]............................................
Global rSem   = CreateSemaphore()
Global rMutex = CreateMutex()
; ...[ Helpers ]..............................................................
Procedure.s PadLeft( n.l, count.l, p$ )
  Protected s$ = Str(n)
  While Len(s$) < count
    s$ = p$ + s$
  Wend
  ProcedureReturn s$
EndProcedure
; ...[ Rendering Thread ].....................................................
Procedure rThread( *ctx.rContext )
  
  ; Local Variables
  Protected ntemp.d,r.v3d
  
  Protected row_s.l = *ctx\row_s
  Protected row_e.l = *ctx\row_e
  Protected tid  .l = *ctx\tid
  SignalSemaphore(rSem)
  
  Protected log$ = "> (tid:"+PadLeft(tid,2,"0")+") rows ["+PadLeft(row_s,4,"0")+":"+PadLeft(row_e,4,"0")+"] : "
  
  ; Loop over image rows
  Protected y.l
  For y = row_s To row_e
    
    ; Log
    LockMutex(rMutex)
    ConsoleLocate(1,tid+3)
    Print( log$ + PadLeft( 100.0*(y-row_s)/(row_e-row_s),3," ")+"%" )
    UnlockMutex(rMutex)
    
    ; Loop over image columns
    Protected x.l
    For x = *ctx\col_s To *ctx\col_e
      
      ; Pixel Index
      Protected i.l = ( *ctx\h - y - 1 )* *ctx\w + x
      
      ; 2x2 subpixel rows
      Protected sy.l
      For sy = 0 To 1
        ; 2x2 subpixel cols
        Protected sx.l
        For sx = 0 To 1
          
          ; samples
          v3Null(r)
          Protected s.l
          For s=0 To *ctx\samps-1
            
            ; double r1 = 2*erand48(Xi), dx=r1<1 ? sqrt(r1)-1: 1-sqrt(2-r1);
            Protected r1.d,dx.d
            r1 = 2.0*rand()
            If r1 < 1.0
              dx = Sqr(r1)-1.0
            Else
              dx = 1.0-Sqr(2.0-r1)
            EndIf
            
            ; double r2 = 2*erand48(Xi), dy=r2<1 ? sqrt(r2)-1: 1-sqrt(2-r2);
            Protected r2.d,dy.d
            r2 = 2.0*rand()
            If r2 < 1.0
              dy = Sqr(r2)-1.0
            Else
              dy = 1.0-Sqr(2.0-r2)
            EndIf
            
            ; Vec d = cx*( ( (sx+.5 + dx)/2 + x)/w - .5) + cy*( ( (sy+.5 + dy)/2 + y)/h - .5) + cam.d;
            Protected d.v3d
            Protected ccx.d = ( ( sx + 0.5 + dx )/2.0 + x )/*ctx\w - 0.5
            Protected ccy.d = ( ( sy + 0.5 + dy )/2.0 + y )/*ctx\h - 0.5
            d\x = *ctx\cx\x*ccx + *ctx\cy\x*ccy + *ctx\cam\d\x
            d\y = *ctx\cx\y*ccx + *ctx\cy\y*ccy + *ctx\cam\d\y
            d\z = *ctx\cx\z*ccx + *ctx\cy\z*ccy + *ctx\cam\d\z
            
            ; r = r + radiance(Ray(cam.o+d*140,d.norm()),0,Xi)*(1./samps);
            ; Camera rays are pushed ^^^^^ forward to start in interior
            Protected rr.v3d,ray.r3d
            v3Null(rr)
            v3Copy(ray\d,d)
            v3Norm(ray\d,ntemp)
            v3ScaleInPlace(d,138.0)
            v3Add(ray\o,*ctx\cam\o,d)
            ;radiance(@rr,@ray,0)
            radiancef(@rr,@ray,0)
            v3AddScaledInPlace(r,rr,*ctx\isamps)
            
          Next
          
          v3Clamp(r)
          v3AddScaledInPlace(*ctx\c(i),r,0.25)
          
        Next ; For sx = 0 To 1
      Next ; For sy = 0 To 1
      
    Next ; For x = *ctx\col_s To *ctx\col_e
  Next ; For y = row_s To row_e

EndProcedure
; ...[ Rendering Main ].......................................................
Procedure render( scene_id.l, samps.l )
  
  Protected ntemp.d
  
  ; Setup Scene
  SetupScene( scene_id )
  
  ; Context
  Protected ctx.rContext
  
  ; Topology
  ctx\w = 1024
  ctx\h =  768
  
  ; Samples
  ctx\samps  = samps
  ctx\isamps = 1.0/samps
  
  ; Cam pos
  ctx\cam\o\x = 50.0
  ctx\cam\o\y = 52.0
  ctx\cam\o\z = 295.6
  ; Cam dir
  ctx\cam\d\x = 0.0
  ctx\cam\d\y = -0.042612
  ctx\cam\d\z = -1.0
  v3Norm(ctx\cam\d,ntemp)
  
  ; Pixels
  ReDim ctx\c(ctx\w*ctx\h-1)
  
  ; Deltas
  Protected cx.v3d,cy.v3d,r.v3d
  ctx\cx\x = ctx\w*0.5135/ctx\h
  ctx\cx\y = 0.0
  ctx\cx\z = 0.0
  v3Cross(ctx\cy,ctx\cx,ctx\cam\d)
  v3Norm(ctx\cy,ntemp)
  v3Scale(ctx\cy,ctx\cy,0.5135)
  
  ; Columns
  ctx\col_s = 0
  ctx\col_e = ctx\w-1
  
  ; Prepare For n Thread(s)
  Dim aThread.l(0)
  Protected ri.l,rThreadCount = CountCPUs(#PB_System_ProcessCPUs)
  ReDim aThread(rThreadCount-1)
  
  ; Image Name
  Protected img_name$
  Select scene_id
    Case #SCENE_CORNELL  : img_name$ = "cornell"
    Case #SCENE_SKY      : img_name$ = "sky"
    Case #SCENE_NIGHTSKY : img_name$ = "nightsky"
    Case #SCENE_ISLAND   : img_name$ = "island"
    Case #SCENE_VISTA    : img_name$ = "vista"
    Case #SCENE_OVERLAP  : img_name$ = "overlap"
    Case #SCENE_WADA     : img_name$ = "wada"
    Case #SCENE_WADA2    : img_name$ = "wada2"
    Case #SCENE_FOREST   : img_name$ = "forest"
  EndSelect

  ; Log
  OpenConsole("purept")
  EnableGraphicalConsole(1)
  ConsoleLocate(1,1)
  PrintN( "Rendering "+ img_name$ +" with "+Str(samps*4)+" spp" )
  
  ; Generate Threads
  Protected row_step      .l = ctx\h / rThreadCount
  Protected row_last_step .l = ctx\h - row_step*(rThreadCount-1)
  Protected row_s         .l = 0
  Protected row_e         .l = row_step
  Protected riStop        .l = ArraySize(aThread())
  For ri=0 To riStop
    ctx\row_s = row_s
    ctx\row_e = row_e - 1
    ctx\tid   = ri
    aThread(ri) = CreateThread(@rThread(),@ctx)
    WaitSemaphore(rSem)
    row_s + row_step
    If ri = riStop - 1
      row_e + row_last_step
    Else
      row_e + row_step
    EndIf
  Next
  
  ; Wait For Rendering End
  Protected Ta.d = ElapsedMilliseconds()
  For ri=0 To riStop
    WaitThread(aThread(ri))
  Next
  Protected Tb.d = ElapsedMilliseconds()
  
  ; Save Image
  LoadFont( 0, "Tahoma", 8, #PB_Font_Bold|#PB_Font_HighQuality )
  CreateImage(0,ctx\w,ctx\h)
  StartDrawing(ImageOutput(0))
  Protected x.l,y.l,i.l = 0
  For y=0 To ctx\h-1
    For x=0 To ctx\w-1
      Plot(x,y,RGB(toInt(ctx\c(i)\x),toInt(ctx\c(i)\y),toInt(ctx\c(i)\z)))
      i + 1
    Next
  Next
  DrawingMode( #PB_2DDrawing_Transparent|#PB_2DDrawing_AlphaBlend )
  Box(0,0,ctx\w,16,RGBA($32,$32,$32,202))
  DrawingFont( FontID(0) )
  FrontColor( RGBA($9A,$98,$98,$FF) )
  Protected Tf$ = FormatDate( "%hhh%iim%sss", AddDate(0,#PB_Date_Second,(Tb-Ta)/1000.0) )
  DrawText( 0, 0, " purept : "+img_name$+" : "+ Tf$ +" : "+Str(samps*4)+" spp" )
  StopDrawing()
  ConsoleLocate(1,rThreadCount+4)
  PrintN( "Saving purept_"+ img_name$ +".bmp.." )
  SaveImage(0,"purept_"+ img_name$ +".bmp")
  PrintN( " Done in "+ Tf$ +"." )
  Input()
  
  CloseConsole()
  
EndProcedure

render( #SCENE_WADA2,    10 ) ; 40 spp (~15s)

;render( #SCENE_CORNELL,  25 ) ; 100 spp (~1m17s)
;render( #SCENE_SKY,      25 ) ; 100 spp (~42s)
;render( #SCENE_NIGHTSKY, 25 ) ; 100 spp (~48s)
;render( #SCENE_ISLAND,   25 ) ; 100 spp (~41s)
;render( #SCENE_VISTA,    25 ) ; 100 spp (~1m12s) ?? Broken Scene, I have to check..
;render( #SCENE_OVERLAP,  25 ) ; 100 spp (~1m01s) ?? Broken Scene, I have to check..
;render( #SCENE_WADA,     25 ) ; 100 spp (~52s)
;render( #SCENE_WADA2,    25 ) ; 100 spp (~38s)
;render( #SCENE_FOREST,   25 ) ; 100 spp (~1m09s)

;render( #SCENE_CORNELL,  50 ) ; 200 spp (~2m11s)

;render( #SCENE_CORNELL,  250 ) ; 1000 spp (~11m31s)

;render( #SCENE_CORNELL, 1250 ) ; 5k spp (~50m04s)

;render( #SCENE_CORNELL, 6250 ) ; 25k spp (not tested yet but several hours)
Have fun!

Cheers,
Guy.
guy rabiller | radfac founder / ceo | raa.tel | raafal.org
Avatar de l’utilisateur
grabiller
Messages : 103
Inscription : lun. 10/sept./2012 11:55
Localisation : France - 89220 Rogny-Les-Septs-Ecluses
Contact :

Re: purept (smallpt) path tracer

Message par grabiller »

blendman a écrit :Sur mon ordi pourri, j'ai lancé en render(10) et ça a mis environ 6min40.

Question subsidiaire : serait-il possible d'ouvrir une scène type .obj (ou mieux .blend ^^), pour en faire un raytracer externe pour soft 3D (blender, max, lightwave..) ?
Ce serait topissime :D
C'est pas trop à l'ordre du jour.

Ce code est vraiment expérimental et me sert surtout à étudier le path tracing. De plus il ne rend que des sphères :wink:

Un raytracer externe en PureBasic ? C'est pas encore pour tout de suite! 8O
guy rabiller | radfac founder / ceo | raa.tel | raafal.org
Avatar de l’utilisateur
graph100
Messages : 1318
Inscription : sam. 21/mai/2005 17:50

Re: purept (smallpt) path tracer

Message par graph100 »

Pour ton "render( #SCENE_CORNELL, 1250 ) ; 5k spp" ça a pris 1h08m23s sur mon portable.
Il y a 8 thread et le proc est un i7-3630QM @ 2.40GHz, la machine tourne sur W8

le wada2 prend 18s avec 40spp
_________________________________________________
Mon site : CeriseCode (Attention Chantier perpétuel ;))
G-Rom
Messages : 3626
Inscription : dim. 10/janv./2010 5:29

Re: purept (smallpt) path tracer

Message par G-Rom »

Très intéressant ton code, même si ce n'est qu'une adaptation c'est du bon travail de recherche.
il fût un temps ou je m'étais lancer la dedans aussi ( http://www.dailymotion.com/video/x9cx7u ... basic_tech , http://www.purebasic.fr/french/viewtopic.php?t=9424 ) , si la loi de moore continue, dans quelques années ou aura du raytracing en temps réel.
Tu pourrais peu être te penché sur OpenCL pour optimiser les calculs.
Backup
Messages : 14526
Inscription : lun. 26/avr./2004 0:40

Re: purept (smallpt) path tracer

Message par Backup »

G-Rom a écrit :dans quelques années ou aura du raytracing en temps réel..
bienvenu dans le futur :)
http://www.pcworld.fr/carte-graphique/a ... 2137,1.htm
le RayCore :)

Image

la news date du : Publié le 4 octobre 2012
G-Rom
Messages : 3626
Inscription : dim. 10/janv./2010 5:29

Re: purept (smallpt) path tracer

Message par G-Rom »

on est jamais dans le futur , des innovations , y en a plein , y avait aussi cela http://www.3dvf.com/actualite-1423-unli ... etour.html , mais il faut qu'après le marché suive, à ce jour , aucun jeu ne tourne en raytracing ^^ , donc c'est pas demain la veille ;)
Avatar de l’utilisateur
grabiller
Messages : 103
Inscription : lun. 10/sept./2012 11:55
Localisation : France - 89220 Rogny-Les-Septs-Ecluses
Contact :

Re: purept (smallpt) path tracer

Message par grabiller »

G-Rom a écrit :Très intéressant ton code, même si ce n'est qu'une adaptation c'est du bon travail de recherche.
il fût un temps ou je m'étais lancer la dedans aussi ( http://www.dailymotion.com/video/x9cx7u ... basic_tech , http://www.purebasic.fr/french/viewtopic.php?t=9424 ) , si la loi de moore continue, dans quelques années ou aura du raytracing en temps réel.
Tu pourrais peu être te penché sur OpenCL pour optimiser les calculs.
Oui pour l'OpenCL c'est prévu, mais seulement une fois que j'aurai implémenté les features que je désire tester/ajouter, avant ça ne sert pas trop à grand chose.

Pour le raytracing temps réel cela existe déjà mais la vraie question est: qu'entend-on par 'raytracing temps réel', avec quelle(s) features exactement ?

De mon côté ce qui m'intéresse c'est de créer des images pour des films (même avec une approche Machinima qui m'intéresse particulièrement) et pour cela il faut un minimum de features:

- antialiasing évolué/efficace
- sampling de texture évolué/efficace (référence: OpenImageIO)
- displacement
- motion blur (camera,transformation,déformation, les 3)
- depth of field
- shading/lighting physiquement 'plausible' + GI
- en option 'forte' volumetric scattering

Le problème est qu'on ne trouve jamais toutes ces features dans un raytracer temps réel, et même pas dans bon nombre de moteur de rendu connus et open source y compris 'accélérés GPU' (notamment le motion blur de déformation).

Mon soucis est donc d'abord d'implémenter les features dont j'ai besoin + des améliorations que je ne vois pas encore dans les moteurs de rendu (y compris commerciaux) comme par exemple du noise reduction efficace niveau samples (et non niveau image après le rendu).

Ce port de smallpt me sert à cela: expérimenter. Qui sait, petit à petit il pourra peut-être être utilisable réellement.
guy rabiller | radfac founder / ceo | raa.tel | raafal.org
G-Rom
Messages : 3626
Inscription : dim. 10/janv./2010 5:29

Re: purept (smallpt) path tracer

Message par G-Rom »

Code : Tout sélectionner

 qu'entend-on par 'raytracing temps réel', avec quelle(s) features exactement ?
60 fps, comme un moteur qui "rasterize" ( ogre & cie )
- antialiasing évolué/efficace
- sampling de texture évolué/efficace (référence: OpenImageIO)
- displacement
- motion blur (camera,transformation,déformation, les 3)
- depth of field
- shading/lighting physiquement 'plausible' + GI
- en option 'forte' volumetric scattering
Dans ton cas, le temps réel n'est pas primordial , du moins , ce n'est pas un objectif.
le rendu dois être rapide & de très bonne qualité , donc exit le temps réel.
La tu parles de post-effect , il n'ont rien à voir avec le raytracer en lui même , il te faut une autre approche.
tu pourrais rendre la scène dans plusieurs tampons différent ( multipasses ) puis les combiné pour le rendu finale.
par exemple , pour le volumetric scattering , tu rends dans un tampon blanc tous les occluders en noir , ensuite tu appliques
un radial blur à la position de la lumière dans l'espace 2D , tu rends dans un autre tampon ta scène & tu applique le 1° tampon au second ( mode multiplication ) & tu as ton effet godray ;)

Pour le noise reduction , n'étant pas un spécialiste , une recherche sur google te donnera des résultats , surtout anglophone.
Bon courage ^^

@++
Avatar de l’utilisateur
grabiller
Messages : 103
Inscription : lun. 10/sept./2012 11:55
Localisation : France - 89220 Rogny-Les-Septs-Ecluses
Contact :

Re: purept (smallpt) path tracer

Message par grabiller »

Effectivement à 60 fps on a pas grand chose en raytracing actuellement, d'autant que dans le temps d'une image il faut traiter de beaucoup d'autres données que le rendu lui-même dans le cadre d'un jeu.

Pour ce dont je parle, non il ne s'agit pas de 'post-effect' dans le cas d'un raytracer mais bien de features inhérentes au moteur de rendu. Tu parles de 'tricks' utilisés avec les moteurs OpenGL/DirectX pour simuler ce genre de features. Dans le cas d'un raytracer, le DOF, le motion blur, le volumic se font au niveau du sampling, pas une fois que l'image est rendu ou avec des effets de superpositions. Avec les moteurs de jeu 'ça le fait' mais pour une images films la qualité n'est pas suffisante et ces 'tricks' sont facilement visibles.

Pour le noise réduction j'ai déjà trouvé, c'est pour cela que j'ai porté smallpt sous PureBasic, pour pouvoir facilement expérimenter. Je compte implémenter ces deux approches complémentaires, l'une pour les fireflies (conséquence de l'importance sampling des lumière) et l'autre pour le noise lui-même inhérent au Monte Carlo path tracing:
Density-based Outlier Rejection in Monte Carlo Rendering
On Filtering the Noise from the Random Parameters in Monte Carlo Rendering
guy rabiller | radfac founder / ceo | raa.tel | raafal.org
Répondre