MiniLight 1.7 a minimal global illumination renderer.

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

MiniLight 1.7 a minimal global illumination renderer.

Postby D.J.Peters » Mar 21, 2015 3:46

MiniLight 1.7 in FreeBASIC C-Version: http://www.hxa.name/minilight/

usage : minilight scenefile

Joshy

scene file: cornellbox.txt
Image
fbc -gen gcc -O 3 minilight.bas
file: minilight.bas

Code: Select all

' MiniLight 1.7 in C from here: http://www.hxa.name/minilight/
' MiniLight 1.7 (optimized version) in FreeBASIC
#include once "crt/time.bi"
#include once "crt/limits.bi"
#include once "crt/math.bi"
#include once "crt/stdio.bi"
#include once "crt/stdlib.bi"
#include once "crt/string.bi"


#define DEBUG
#ifdef DEBUG
 #macro dprint(txt)
  scope
   dim as integer hFile=FreeFile()
   open err for output as #hFile:
   print #hFile,txt
   close #hFile
  end scope
  #endmacro
#else
 #define dprint(txt) :
#endif

#ifndef USE_DOUBLE
  type real as single
  const as real REAL_MAX = FLT_MAX
  #define realSQRT(x)  sqrtf(x)
  #define realPOW(x,y) powf(x,y)
  #define realTAN(x)   tanf(x)
  #define realCOS(x)   cosf(x)
  #define realSIN(x)   sinf(x)
#else
  type real as double
  const as real REAL_MAX = DBL_MAX
  #define realSQRT(x)  sqrt(x)
  #define realPOW(x,y) pow(x,y)
  #define realTAN(x)   tan_(x)
  #define realCOS(x)   cos_(x)
  #define realSIN(x)   sin_(x)
#endif

type bool as integer
const as real EPSILON               = 1.0 / 1048576.0
const as real TOLERANCE             = 1.0 / 1024.0
const as real PI                    = atn(1)*4.0
const as real PI2                   = PI*2.0
const as real VIEW_ANGLE_MIN        =  10.0
const as real VIEW_ANGLE_MAX        = 160.0
const as real DISPLAY_LUMINANCE_MAX = 200.0
const as real GAMMA_ENCODE          = 0.45

const as integer MAX_NODE_LEVELS    = 44
const as integer MAX_NODE_ITEMS     =  8
const as real    MIN_NODE_TOLERANCE = TOLERANCE * 4.0
const as integer MAX_TRIANGLES      = &H1000000
const as integer IMAGE_DIM_MAX      = 4000


type realVector3
  as real xyz(2)
end type
dim shared as realVector3 realVector3ZERO = Type({0,0,0})
dim shared as realVector3 realVector3ONE  = Type({1,1,1})
' ITU-R BT.709 standard RGB luminance weighting
dim shared as realVector3 RGB_LUMINANCE = type({0.2126,0.7152,0.0722})

function realVector3Read(byval pIn as FILE ptr) as realVector3
  dim as realVector3 v=any
  dim as zstring * 2 l,r
  dim as single f(2)=any
  fscanf( pIn, "%1s %g %g %g %1s",@l, @f(0), @f(1), @f(2),@r)
  if l[0]<>asc("(") or r[0]<>asc(")") then
    print "error: realVector3Read() !"
    beep:sleep:end
  end if
  for i as integer=0 to 2:v.xyz(i)=f(i):next
  return v
end function

function realVector3Add(byval pV0 as const realVector3 ptr,byval pV1 as const realVector3 ptr) as realVector3
   dim as realVector3 v=any
   v.xyz(0) = pV0->xyz(0) + pV1->xyz(0)
   v.xyz(1) = pV0->xyz(1) + pV1->xyz(1)
   v.xyz(2) = pV0->xyz(2) + pV1->xyz(2)
   return v
end function

function realVector3Sub(byval pV0 as const realVector3 ptr,byval pV1 as const realVector3 ptr) as realVector3
   dim as realVector3 v=any
   v.xyz(0) = pV0->xyz(0) - pV1->xyz(0)
   v.xyz(1) = pV0->xyz(1) - pV1->xyz(1)
   v.xyz(2) = pV0->xyz(2) - pV1->xyz(2)
   return v
end function

function realVector3MulV(byval pV0 as const realVector3 ptr,byval pV1 as const realVector3 ptr) as realVector3
   dim as realVector3 v=any
   v.xyz(0) = pV0->xyz(0) * pV1->xyz(0)
   v.xyz(1) = pV0->xyz(1) * pV1->xyz(1)
   v.xyz(2) = pV0->xyz(2) * pV1->xyz(2)
   return v
end function

function realVector3MulF(byval pV as const realVector3 ptr,byval f as real) as realVector3
  dim as realVector3 v=any
  v.xyz(0) = pV->xyz(0) * f
  v.xyz(1) = pV->xyz(1) * f
  v.xyz(2) = pV->xyz(2) * f
  return v
end function

function realVector3Negative(byval pV as const realVector3 ptr) as realVector3
  dim as realVector3 v=any
  v.xyz(0) = -pV->xyz(0)
  v.xyz(1) = -pV->xyz(1)
  v.xyz(2) = -pV->xyz(2)
  return v
end function

function realVector3Dot(byval pV0 as const realVector3 ptr,byval pV1 as const realVector3 ptr) as real
  return pV0->xyz(0) * pV1->xyz(0) + pV0->xyz(1) * pV1->xyz(1) + pV0->xyz(2) * pV1->xyz(2)
end function

function realVector3Unitized(byval pV as const realVector3 ptr) as realVector3
  dim as real length = realSQRT( realVector3Dot( pV, pV ) )
  length = iif (length<>0.0, 1.0 / length,0.0)
  return realVector3MulF( pV, length )
end function

function realVector3Cross(byval pV0 as const realVector3 ptr,byval pV1 as const realVector3 ptr) as realVector3
  dim as realVector3 v=any
  v.xyz(0) = pV0->xyz(1)*pV1->xyz(2) - pV0->xyz(2)*pV1->xyz(1)
  v.xyz(1) = pV0->xyz(2)*pV1->xyz(0) - pV0->xyz(0)*pV1->xyz(2)
  v.xyz(2) = pV0->xyz(0)*pV1->xyz(1) - pV0->xyz(1)*pV1->xyz(0)
  return v
end function

function realVector3IsZero(byval pV as const realVector3 ptr) as bool
  if pV->xyz(0)<>0.0 orelse _
     pV->xyz(1)<>0.0 orelse _
     pV->xyz(2)<>0.0 then return 0
  return 1
end function

function realVector3Clamped(byval pV   as const realVector3 ptr, _
                            byval pMin as const realVector3 ptr, _
                            byval pMax as const realVector3 ptr) as realVector3
  dim as realVector3 v=any
  v.xyz(0) = iif(pV->xyz(0) < pMin->xyz(0), pMin->xyz(0),iif(pV->xyz(0) > pMax->xyz(0),pMax->xyz(0),pV->xyz(0)))
  v.xyz(1) = iif(pV->xyz(1) < pMin->xyz(1), pMin->xyz(1),iif(pV->xyz(1) > pMax->xyz(1),pMax->xyz(1),pV->xyz(1)))
  v.xyz(2) = iif(pV->xyz(2) < pMin->xyz(2), pMin->xyz(2),iif(pV->xyz(2) > pMax->xyz(2),pMax->xyz(2),pV->xyz(2)))
  return v
end function

function realVector3Interpolate(byval pV0 as const realVector3 ptr, _
                                byval pV1 as const realVector3 ptr, _
                                byval pV2 as const realVector3 ptr, _
                                byval u as real, _
                                byval v as real) as  realVector3
  dim as realVector3 r=any
  dim as real uv=1.0 - (u + v)
  r.xyz(0) = -(uv * pV1->xyz(0) + pV0->xyz(0) * u + pV2->xyz(0) * v)
  r.xyz(1) = -(uv * pV1->xyz(1) + pV0->xyz(1) * u + pV2->xyz(1) * v)
  r.xyz(2) = -(uv * pV1->xyz(2) + pV0->xyz(2) * u + pV2->xyz(2) * v)
  return r
end function

function realVector3RGB(byval pV as const realVector3 ptr) as ulong
  dim as integer r=pV->xyz(0)
  dim as integer g=pV->xyz(1)
  dim as integer b=pV->xyz(2)
  return RGB(r,g,b)
end function

type Triangle
  as realVector3 aVertexs(2)
  as realVector3 edge1
  as realVector3 edge2
  as realVector3 edge3
  as realVector3 normal
  as realVector3 tangent
  as real        area
  as real        reflectdot3
  as realVector3 reflectcolor
  as realVector3 reflectivity
  as realVector3 emitivity
end type

function TriangleNormalV(byval pT as const Triangle ptr) as realVector3
   dim as realVector3 edge1 = realVector3Sub( @pT->aVertexs(1), @pT->aVertexs(0) )
   dim as realVector3 edge3 = realVector3Sub( @pT->aVertexs(2), @pT->aVertexs(1) )
   return realVector3Cross( @edge1, @edge3 )
end function

function TriangleNormal(byval pT as const Triangle ptr) as realVector3
  dim as realVector3 normalV = TriangleNormalV( pT )
  return realVector3Unitized( @normalV )
end function

function TriangleTangent(byval pT as const Triangle ptr) as realVector3
  dim as realVector3 edge1 = realVector3Sub( @pT->aVertexs(1), @pT->aVertexs(0) )
  return realVector3Unitized( @edge1 )
end function

function TriangleArea(byval pT as const Triangle ptr) as real
  dim as realVector3 normalV = TriangleNormalV( pT )
  return realSQRT( realVector3Dot( @normalV, @normalV ) ) * 0.5
end function


function TriangleCreate(byval pIn as FILE ptr) as Triangle
  dim as Triangle t=any
  for i as integer = 0 to 2
    t.aVertexs(i) = realVector3Read(pIn)
  next
  t.edge1        = realVector3Sub(@t.aVertexs(1),@t.aVertexs(0))
  t.tangent      = realVector3Unitized( @t.edge1 )
  t.edge2        = realVector3Sub(@t.aVertexs(2),@t.aVertexs(0))
  t.edge3        = realVector3Sub(@t.aVertexs(2),@t.aVertexs(1))
  t.reflectivity = realVector3Read(pIn)
  t.emitivity    = realVector3Read(pIn)
  t.normal       = TriangleNormal(@t)
  t.area         = TriangleArea(@t)
  t.reflectivity = realVector3Clamped(@t.reflectivity,@realVector3ZERO,@realVector3ONE)
  t.reflectdot3  = realVector3Dot(@t.reflectivity,@realVector3ONE ) / 3.0
  t.reflectcolor = realVector3MulF(@t.reflectivity,1.0 / t.reflectdot3)
  t.emitivity    = realVector3Clamped(@t.emitivity,@realVector3ZERO,@t.emitivity)
  return t
end function

sub TriangleBound(byval pT as const Triangle ptr, _
                  aBound_o() as real)
  ' initialise to one vertex
  for i as integer=0 to 5
    aBound_o(i) = pT->aVertexs(2).xyz(i mod 3)
  next
  ' expand to surround all vertexs
  for i as integer= 0 to 2
    dim as integer d,m
    for j as integer = 0 to 5
      d = j \ 3
      m = j mod 3
      dim as real v = pT->aVertexs(i).xyz(m) + iif (d,1.0,-1.0) * TOLERANCE
      dim as integer g=iif(aBound_o(j)>v,1,0)
      if (g xor d) then aBound_o(j)=v
    next
  next
end sub


function TriangleIntersection(byval pTriangle     as const Triangle    ptr, _
                              byval pRayOrigin    as const realVector3 ptr , _
                              byval pRayDirection as const realVector3 ptr, _
                              byval pOutDistance  as real ptr) as bool
  dim as realVector3 pvec = realVector3Cross( pRayDirection, @pTriangle->edge2 )
  dim as real det = pTriangle->edge1.xyz(0)*pvec.xyz(0) _
                  + pTriangle->edge1.xyz(1)*pvec.xyz(1) _
                  + pTriangle->edge1.xyz(2)*pvec.xyz(2)

  if abs(det)<EPSILON then return 0
  det = 1.0 / det
  dim as realVector3 tvec = realVector3Sub( pRayOrigin, @pTriangle->aVertexs(0))
  dim as real u = tvec.xyz(0)*pvec.xyz(0) _
                + tvec.xyz(1)*pvec.xyz(1) _
                + tvec.xyz(2)*pvec.xyz(2)
  u *= det
  if u<0 then return 0
  if u>1 then return 0
  dim as realVector3 qvec = realVector3Cross( @tvec, @pTriangle->edge1 )
  dim as real v = pRayDirection->xyz(0)*qvec.xyz(0) _
                + pRayDirection->xyz(1)*qvec.xyz(1) _
                + pRayDirection->xyz(2)*qvec.xyz(2)
  v *= det
  if v<0 then return 0
  if u+v >1 then return 0
  *pOutDistance = pTriangle->edge2.xyz(0)*qvec.xyz(0) _
                + pTriangle->edge2.xyz(1)*qvec.xyz(1) _
                + pTriangle->edge2.xyz(2)*qvec.xyz(2)
  *pOutDistance *= det
  if *pOutDistance< 0 then return 0
  return 1
end function

function TriangleRandomSamplePoint(byval pTriangle as const Triangle ptr) as realVector3
  dim as realVector3 SamplePoint=any
  dim as real rndSQR = realSQRT(rnd())
  dim as real u =  1.0 - rndSQR
  dim as real v = (1.0 - rnd()) * rndSQR
  SamplePoint.xyz(0)=pTriangle->aVertexs(0).xyz(0) + pTriangle->edge1.xyz(0)*u + pTriangle->edge2.xyz(0)*v
  SamplePoint.xyz(1)=pTriangle->aVertexs(0).xyz(1) + pTriangle->edge1.xyz(1)*u + pTriangle->edge2.xyz(1)*v
  SamplePoint.xyz(2)=pTriangle->aVertexs(0).xyz(2) + pTriangle->edge1.xyz(2)*u + pTriangle->edge2.xyz(2)*v
  return SamplePoint
end function

type AABBNode
  as bool              isBranch
  as real              Bound(5)
  as const any ptr ptr pItems ' subcells or leaf's
  as integer           nItems
end type

sub constructAABBNode(byval pInItems as const Triangle ptr ptr, _
                      byval nInItems as const integer, _
                      byval level    as const integer, _
                      byval pOutNode as AABBNode ptr)
  dim as bool blnToBigForLeaf = iif(nInItems > MAX_NODE_ITEMS,1,0)
  dim as bool blnLevelDeepOk  = iif(level < (MAX_NODE_LEVELS - 1),1,0)
  ' is branch if items overflow leaf and tree not too deep
  pOutNode->isBranch = iif(blnToBigForLeaf and blnLevelDeepOk,1,0)

  if( pOutNode->isBranch) then
    ' make branch: make sub-cells, and recurse construction
    dim as integer q
    pOutNode->nItems  = 8
    pOutNode->pItems = calloc( pOutNode->nItems, sizeof(any ptr))
    for s as integer = 0 to pOutNode->nItems-1
      dim as real SubBound(5)=any
      dim as integer nSubItems = 0
      dim as const Triangle ptr ptr pSubItems = calloc( nSubItems, sizeof(Triangle ptr))
      dim as integer d, m, sm
      ' make subcell bound
      for j as integer = 0 to 5
        d = j \ 3 : m = j mod 3 : sm=s shr m : sm and=1
        if (sm xor d) then
          SubBound(j) = (pOutNode->Bound(m) + pOutNode->Bound(m + 3)) * 0.5
        else
          SubBound(j) = pOutNode->Bound(j)
        end if
      next
      ' collect items that overlap subcell
      for i as integer = 0 to nInItems-1
        dim as integer d1,dg,isOverlap = 1
        dim as real InBound(5)
        TriangleBound( pInItems[i], InBound() )
        for j as integer = 0 to 5
          d = j  \  3 : m = j mod 3 : d1=d xor 1
          dg=iif(InBound( d1 * 3 + m) >= SubBound(j),1,0)
          dg xor= d
          isOverlap and=dg
        next
        ' maybe append to subitems store
        if( isOverlap ) then
          nSubItems+=1
          pSubItems = realloc(cptr(Triangle ptr ptr,pSubItems),nSubItems * sizeof(Triangle ptr) )
          pSubItems[nSubItems - 1] = pInItems[i]
        end if
      next
      if (nSubItems=nInItems) then q+=1
      ' maybe make subcell, if any overlapping subitems
      if( nSubItems > 0 ) then
        dim as AABBNode ptr pS = calloc( 1, sizeof(AABBNode) )
        ' curtail degenerate subdivision by adjusting next level
        ' degenerate if two or more subcells copy entire contents of parent, or if subdivision reaches below mm size
        ' having a model including the sun requires one subcell copying entire contents of parent to be allowed
        dim as integer q1=iif(q>1,1,0)
        ' if xBoundSize < shortestBoundSize  call rerurse with last level (MAX_NODE_LEVELS)
        dim as integer nl=iif((SubBound(3)-SubBound(0)) < MIN_NODE_TOLERANCE,MAX_NODE_LEVELS,level + 1)
        dim as integer nextLevel = q1 or nl
        pOutNode->pItems[s] = pS
        for i as integer = 0 to 5
          pS->Bound(i) = SubBound(i)
        next
        ' recurse
        constructAABBNode(pSubItems, nSubItems, nextLevel, pS)
      end if
      free(pSubItems)
    next
  else
    ' make leaf: store items and exit recursion
    pOutNode->nItems = nInItems
    pOutNode->pItems = calloc( pOutNode->nItems, sizeof(any ptr) )
    for i as integer = 0 to pOutNode->nItems-1
      pOutNode->pItems[i] = pInItems[i]
    next
  end if
end sub


function AABBTreeConstruct(byval pEyePosition as const realVector3 ptr, _
                           byval pTriangles   as const Triangle ptr, _
                           byval nTriangles   as integer) as const AABBNode ptr
  dim as AABBNode ptr pRootNode = calloc( 1, sizeof(AABBNode) )
  ' set overall bound (and convert to collection of pointers)
  dim as const Triangle ptr ptr tempTriangleList = calloc( nTriangles, sizeof(Triangle ptr) )
  ' accommodate eye position (makes tracing algorithm simpler)
  for i as integer = 0 to 5
    pRootNode->Bound(i) = pEyePosition->xyz(i mod 3)
  next
  ' accommodate all items
  for i as integer = 0 to nTriangles-1
    tempTriangleList[i] = @pTriangles[i]
    dim as real aItemBound(5)
    TriangleBound( @pTriangles[i], aItemBound() )
    ' accommodate item
    for j as integer = 0 to 5
      dim as integer jg2=iif(j > 2,1,0)
      dim as integer bg=iif(pRootNode->Bound(j) > aItemBound(j),1,0)
      if ( bg xor jg2 ) then pRootNode->Bound(j) = aItemBound(j)
    next
  next
  ' make cubical
  dim as real maxSize = 0.0
  ' find max dimension
  for i as integer = 0 to 2
    dim as real diff=pRootNode->Bound(3+i) - pRootNode->Bound(i)
    if( maxSize < diff ) then maxSize = diff
  next
  ' set all dimensions to max
  for i as integer = 0 to 2
    dim as real bmax = pRootNode->Bound(i) + maxSize
    if ( pRootNode->Bound(3+i) < bmax ) then pRootNode->Bound(3+i) = bmax
  next
  ' make subcell tree
  constructAABBNode(tempTriangleList,nTriangles,0,pRootNode)
  free (tempTriangleList)
  return pRootNode
end function

sub AABBTreeDestruct(byval pNode as AABBNode ptr)
  dim as integer i=pNode->nItems-1
  while pNode->isBranch and (i>0)
    if ( pNode->pItems[i] ) then
      AABBTreeDestruct( cptr(AABBNode ptr,pNode->pItems[i]) )
    end if
    i-=1
  wend
  free(pNode->pItems)
  free(pNode)
end sub


sub AABBNodeIntersection(byval pSpatialIndex   as const AABBNode ptr, _
                         byval pRayOrigin      as const realVector3 ptr, _
                         byval pRayDirection   as const realVector3 ptr, _
                         byval lastHit         as const any ptr, _
                         byval pStart          as const realVector3 ptr, _
                         byval ppOutHitObject  as const Triangle ptr ptr, _
                         byval pOutHitPosition as realVector3 ptr)
   
  if(pSpatialIndex->isBranch) then
    ' is branch: step through subcells and recurse
    dim as integer subCell
    dim as realVector3 cellPosition=any
    if pStart=0 then pStart=pRayOrigin
    ' find which subcell holds ray origin (ray origin is inside cell)
    for i as integer = 0 to 2
      ' compare dimension with center
      dim as integer blnBit =iif( pStart->xyz(i) >= ((pSpatialIndex->Bound(i) + pSpatialIndex->Bound(i+3)) * 0.5),1,0)
      dim as integer shiftBit=blnBit shl i
      subCell or= shiftBit
    next
    ' step through intersected subcells
    cellPosition = *pStart
    while 1
      dim as integer axis = 2, i
      dim as real steps(2)
      if(pSpatialIndex->pItems[subCell]) then
        ' intersect subcell (by recursing)
        AABBNodeIntersection(cptr(AABBNode ptr,pSpatialIndex->pItems[subCell]), _
                             pRayOrigin, _
                             pRayDirection, _
                             lastHit, _
                             @cellPosition, _
                             ppOutHitObject, _
                             pOutHitPosition)
        ' exit branch if item hit
        if ( *ppOutHitObject ) then exit while
      end if
      ' find next subcell ray moves to (by finding which face of the corner ahead is crossed first)
      '  for( i = 3;  i-- > 0;  axis = step[i] < step[axis] ? i : axis )
      for i as integer = 2 to 0 step -1
        ' find which face (inter-/outer-) the ray is heading for (in this dimension)
        dim as bool blnHigh = iif((subCell shr i) and 1,1,0)
        dim as bool lessZero = iif(pRayDirection->xyz(i)<0.0 ,1,0)
        dim as real face = iif(lessZero xor blnHigh , pSpatialIndex->Bound(i + (blnHigh * 3)) , (pSpatialIndex->Bound(i) + pSpatialIndex->Bound(i + 3)) * 0.5)
        ' calculate distance to face  (div by zero produces infinity, which is later discarded)
        steps(i) = (face - pRayOrigin->xyz(i)) / pRayDirection->xyz(i)
        ' last clause of for-statement notes nearest so far
        axis = iif(steps(i) < steps(axis),i,axis)
      next
      dim as bool IsNegative = iif(pRayDirection->xyz(axis)<0.0,1,0)
      dim as bool IsHigh     = iif((subCell shr axis) and 1,1,0)
      ' leaving branch
      ' if: direction is negative and subcell is low
      '  or direction is positive and subcell is high
      if ( IsHigh xor IsNegative) then exit while
      ' move to outer face of next subcell
      dim as realVector3 rs = realVector3MulF( pRayDirection, steps(axis))
      cellPosition = realVector3Add( pRayOrigin, @rs )
      dim as bool shiftAxis=1 shl axis
      subCell xor= shiftAxis
    wend
  else
    ' is leaf: exhaustively intersect contained items
    dim as real nearestDistance = REAL_MAX
    *ppOutHitObject = 0
    ' step through items
    for i as integer = 0 to pSpatialIndex->nItems-1
      dim as const Triangle ptr pItem = cptr(Triangle ptr,pSpatialIndex->pItems[i])
      ' avoid spurious intersection with surface just come from
      if( pItem <> lastHit ) then
        dim as real distance = REAL_MAX
        ' intersect ray with item, and inspect if nearest so far
        if( TriangleIntersection( pItem, pRayOrigin, pRayDirection,@distance ) andalso (distance < nearestDistance) ) then
          ' check intersection is inside cell bound (with tolerance)
          dim as realVector3 hit = any
          hit.xyz(0)=pRayOrigin->xyz(0)+distance*pRayDirection->xyz(0)
          hit.xyz(1)=pRayOrigin->xyz(1)+distance*pRayDirection->xyz(1)
          hit.xyz(2)=pRayOrigin->xyz(2)+distance*pRayDirection->xyz(2)
          if( ((pSpatialIndex->Bound(0) - hit.xyz(0)) <= TOLERANCE) andalso ((hit.xyz(0) - pSpatialIndex->Bound(3)) <= TOLERANCE) andalso _
              ((pSpatialIndex->Bound(1) - hit.xyz(1)) <= TOLERANCE) andalso ((hit.xyz(1) - pSpatialIndex->Bound(4)) <= TOLERANCE) andalso _
              ((pSpatialIndex->Bound(2) - hit.xyz(2)) <= TOLERANCE) andalso ((hit.xyz(2) - pSpatialIndex->Bound(5)) <= TOLERANCE) ) then
            ' note nearest so far
            *ppOutHitObject  = pItem
            nearestDistance = distance
            *pOutHitPosition = hit
          end if
        end if
      end if
    next
  end if
end sub

type SurfacePoint
  as const Triangle ptr pTriangle
  as realVector3        position
end type
'#define SurfacePointHitId( pS ) ( cptr(any ptr,(pS)->pTriangle) )

function SurfacePointCreate(byval pTriangle as const Triangle ptr, _
                            byval pPosition as const realVector3 ptr) as SurfacePoint
  dim as SurfacePoint s=any
  s.pTriangle = pTriangle
  s.position  = *pPosition
  return s
end function

function SurfacePointEmission(byval pSurfacePoint as const SurfacePoint ptr, _
                              byval pToPosition   as const realVector3 ptr, _
                              byval pOutDirection as const realVector3 ptr, _
                              byval isSolidAngle  as bool) as realVector3
  dim as real cosOut = pOutDirection->xyz(0)*pSurfacePoint->pTriangle->normal.xyz(0) _
                     + pOutDirection->xyz(1)*pSurfacePoint->pTriangle->normal.xyz(1) _
                     + pOutDirection->xyz(2)*pSurfacePoint->pTriangle->normal.xyz(2)
  ' emit from front face of surface only
  if (cosOut<=0) then return realVector3ZERO
  if (isSolidAngle=0) then return pSurfacePoint->pTriangle->emitivity
  dim as realVector3 r=type({pToPosition->xyz(0)-pSurfacePoint->position.xyz(0), _
                             pToPosition->xyz(1)-pSurfacePoint->position.xyz(1), _
                             pToPosition->xyz(2)-pSurfacePoint->position.xyz(2)})
  dim as real dot = r.xyz(0)*r.xyz(0) + r.xyz(1)*r.xyz(1) + r.xyz(2)*r.xyz(2)
  if dot <1e-6 then dot=1e-6
  dim as real solidAngle=(cosOut * pSurfacePoint->pTriangle->area) / dot
  r.xyz(0)=pSurfacePoint->pTriangle->emitivity.xyz(0)*solidAngle
  r.xyz(1)=pSurfacePoint->pTriangle->emitivity.xyz(1)*solidAngle
  r.xyz(2)=pSurfacePoint->pTriangle->emitivity.xyz(2)*solidAngle
  return r
end function

function SurfacePointReflection(byval pSurfacePoint as const SurfacePoint ptr, _
                                byval pInDirection  as const realVector3 ptr, _
                                byval pInRadiance   as const realVector3 ptr, _
                                byval pOutDirection as const realVector3 ptr) as realVector3
  dim as real inDot  = pInDirection->xyz(0)*pSurfacePoint->pTriangle->normal.xyz(0) _
                     + pInDirection->xyz(1)*pSurfacePoint->pTriangle->normal.xyz(1) _
                     + pInDirection->xyz(2)*pSurfacePoint->pTriangle->normal.xyz(2)

  dim as real outDot = pOutDirection->xyz(0)*pSurfacePoint->pTriangle->normal.xyz(0) _
                     + pOutDirection->xyz(1)*pSurfacePoint->pTriangle->normal.xyz(1) _
                     + pOutDirection->xyz(2)*pSurfacePoint->pTriangle->normal.xyz(2)
  ' directions must be on same side of surface (no transmission)
  dim as bool isSameSide = iif(iif(inDot<0, 1,0) xor iif(outDot<0,1,0),0,1)
  if isSameSide=0 then return realVector3ZERO
  if inDot<0 then inDot=-inDot
  inDot/=PI
  dim as realVector3 r=any
  r.xyz(0)=pInRadiance->xyz(0)*pSurfacePoint->pTriangle->reflectivity.xyz(0)*inDot
  r.xyz(1)=pInRadiance->xyz(1)*pSurfacePoint->pTriangle->reflectivity.xyz(1)*inDot
  r.xyz(2)=pInRadiance->xyz(2)*pSurfacePoint->pTriangle->reflectivity.xyz(2)*inDot
  return r
end function

function SurfacePointNextDirection(byval pSurfacePoint as const SurfacePoint ptr, _
                                   byval pInDirection  as const realVector3 ptr, _
                                   byval pOutDirection as realVector3 ptr, _
                                   byval pOutColor     as realVector3 ptr) as bool
  dim as bool isAlive = iif(rnd()< pSurfacePoint->pTriangle->reflectdot3,1,0)
  if isAlive=0 then return 0
  'if rnd()<.5 then return 0
  dim as real rndPI2 = rnd()* PI2
  dim as real rndSQR = realSQRT(rnd())
  dim as real x = realCOS(rndPI2) * rndSQR
  dim as real y = realSIN(rndPI2) * rndSQR
  dim as real z = realSQRT( 1.0 - (rndSQR * rndSQR) )
  dim as realVector3 t = pSurfacePoint->pTriangle->tangent
  dim as realVector3 n = pSurfacePoint->pTriangle->normal
  dim as real dot =n.xyz(0)*pInDirection->xyz(0) _
                  +n.xyz(1)*pInDirection->xyz(1) _
                  +n.xyz(2)*pInDirection->xyz(2)
  if dot<0 then
    n.xyz(0) = -n.xyz(0)
    n.xyz(1) = -n.xyz(1)
    n.xyz(2) = -n.xyz(2)
  end if
  dim as realVector3 c = realVector3Cross(@n,@t)
  pOutDirection->xyz(0)=t.xyz(0)*x+c.xyz(0)*y+n.xyz(0)*z
  if pOutDirection->xyz(0) then
    pOutDirection->xyz(1)=t.xyz(1)*x+c.xyz(1)*y+n.xyz(1)*z
    if pOutDirection->xyz(1) then
      pOutDirection->xyz(2)=t.xyz(2)*x+c.xyz(2)*y+n.xyz(2)*z
      if pOutDirection->xyz(2) then
        *pOutColor = pSurfacePoint->pTriangle->reflectcolor
        return 1
      end if
    end if
  end if
  return 0
end function

type Scene
  as Triangle ptr     pTriangles
  as integer          nTriangles
  as Triangle ptr ptr pEmitters
  as integer          nEmitters
  as AABBNode ptr     pRoot
  as realVector3      skyEmission
  as realVector3      groundReflection
  as realVector3      groundEmission
end type

function SceneConstruct(byval pIn as FILE ptr, _
                        byval pEyePosition as const realVector3 ptr) as const Scene ptr
  dim as Scene ptr pScene = calloc( 1, sizeof(Scene))
  ' read and condition background sky and ground values
  scope
    pScene->skyEmission      = realVector3Read( pIn)
    pScene->groundReflection = realVector3Read( pIn)
    pScene->skyEmission      = realVector3Clamped(@pScene->skyEmission     ,@realVector3ZERO, @pScene->skyEmission )
    pScene->groundReflection = realVector3Clamped(@pScene->groundReflection,@realVector3ZERO, @realVector3ONE )
    pScene->groundEmission   = realVector3MulV   (@pScene->skyEmission     ,@pScene->groundReflection)
  end scope
  ' read objects, until end of file or until maximum reached
  scope
    pScene->pTriangles = calloc( 0, sizeof(Triangle) )
    pScene->nTriangles = 0
    for i as integer = 0 to MAX_TRIANGLES-1
      ' stop reading if no more objects
      scope
        ' read next non blank char
        dim as zstring * 2 s
        dim as integer r = fscanf( pIn, "%1s", s )
        ' throw non-EOF failure
        dim as integer code = iif(ferror( pIn ),1,0)
        clearerr( pIn )
        if (code) then
          print "error: SceneConstruct()"
          beep:sleep:end
        end if
        ' if char was found, put back, else end reading
        if( 1 = r ) then
          if ungetc( s[0], pIn )=EOF_ then
            print "error: SceneConstruct() ungetc()"
            beep:sleep:end
          end if
        else
          exit for
        end if
      end scope
      ' read an object
      scope
        dim as Triangle t = TriangleCreate(pIn)
        ' append to objects storage
        pScene->nTriangles+=1
        pScene->pTriangles = realloc( pScene->pTriangles, pScene->nTriangles * sizeof(Triangle))
        pScene->pTriangles[pScene->nTriangles - 1] = t
      end scope
    next
  end scope
  ' find emitting objects
  scope
    pScene->pEmitters = calloc( 0, sizeof(Triangle ptr))
    pScene->nEmitters = 0
    if pScene->nTriangles>0 then
      for i as integer = 0 to pScene->nTriangles-1
        ' has non-zero emission and area
        if( (realVector3IsZero( @pScene->pTriangles[i].emitivity )=0) andalso  (pScene->pTriangles[i].area > 0.0) ) then
          ' append to emitters storage
          pScene->nEmitters+=1
          pScene->pEmitters = realloc(pScene->pEmitters, pScene->nEmitters * sizeof(Triangle ptr))
          pScene->pEmitters[pScene->nEmitters - 1] = @pScene->pTriangles[i]
        end if
      next
    else
      print "warning: SceneConstruct() pScene->trianglesLength=0"
      beep:sleep
    end if
  end scope
  ' make index of objects
  pScene->pRoot = cptr(AABBNode ptr,AABBTreeConstruct( pEyePosition, pScene->pTriangles, pScene->nTriangles))
  return pScene
end function

sub SceneDestruct(byval pScene as Scene ptr)
  if pScene then
    if pScene->pRoot      then AABBTreeDestruct( pScene->pRoot )
    if pScene->pEmitters  then free(pScene->pEmitters)
    if pScene->pTriangles then free(pScene->pTriangles)
    free(pScene)
  end if
end sub

sub SceneRandomEmitter(byval pScene       as const Scene ptr, _
                       byval pOutPosition as realVector3 ptr, _
                       byval ppOutID      as const Triangle ptr ptr)
  if( pScene->nEmitters > 0) then
    ' select emitter
    dim as integer index = floor(Rnd() * pScene->nEmitters)
    index = iif(index < pScene->nEmitters ,index,pScene->nEmitters - 1)
    ' choose random position on emitter
    *pOutPosition = TriangleRandomSamplePoint(pScene->pEmitters[index])
    *ppOutID       = pScene->pEmitters[index]
   else
    *pOutPosition = realVector3ZERO
    *ppOutID       = 0
  end if
end sub

sub Scene_Intersection(byval pScene         as const Scene ptr, _
                      byval pRayOrigin     as const realVector3 ptr, _
                      byval pRayDirection  as const realVector3 ptr, _
                      byval lastHit        as const any ptr, _
                      byval ppHitObject_o  as const Triangle ptr ptr, _
                      byval pHitPosition_o as realVector3 ptr)
   AABBNodeIntersection( pScene->pRoot,pRayOrigin,pRayDirection,lastHit,0,ppHitObject_o,pHitPosition_o)
end sub

type RayTracer
  as const Scene ptr pScene
end type

function RayTracerCreate(byval pScene as const Scene ptr) as RayTracer
  dim as RayTracer r
  r.pScene = pScene
  return r
end function

function sampleEmitters(byval pRayTracer        as const RayTracer ptr, _
                        byval pRayBackDirection as const realVector3 ptr, _
                        byval pSurfacePoint     as const SurfacePoint ptr) as realVector3
  dim as realVector3 radiance '= realVector3ZERO
  ' get position on an emitter
  dim as realVector3 emitterPosition=any
  dim as const Triangle ptr emitterId = 0
  SceneRandomEmitter(pRayTracer->pScene,@emitterPosition,@emitterId )
  ' check an emitter was found
  if (emitterId) then
    ' make direction to emit point
    dim as realVector3 emitVector    = realVector3Sub(@emitterPosition,@pSurfacePoint->position)
    dim as realVector3 emitDirection = realVector3Unitized(@emitVector)
    ' send shadow ray
    dim as const Triangle ptr pHitObject = 0
    dim as realVector3 hitPosition
    AABBNodeIntersection(pRayTracer->pScene->pRoot,@pSurfacePoint->position,@emitDirection,pSurfacePoint->pTriangle,0,@pHitObject,@hitPosition)
    '  SceneIntersection(pRayTracer->pScene       ,@pSurfacePoint->position,@emitDirection,pSurfacePoint->pTriangle  ,@pHitObject,@hitPosition)
    ' check if unshadowed
    if( (pHitObject=0) orelse (emitterId = pHitObject) ) then
      ' get inward emission value
      dim as SurfacePoint sp = SurfacePointCreate(emitterId,@emitterPosition)
      dim as realVector3 backEmitDirection = realVector3Negative(@emitDirection)
      dim as realVector3 emissionIn = SurfacePointEmission(@sp,@pSurfacePoint->position,@backEmitDirection,TRUE)
      dim as realVector3 emissionAll = realVector3MulF(@emissionIn,pRayTracer->pScene->nEmitters)
      ' get amount reflected by surface
      radiance = SurfacePointReflection(pSurfacePoint,@emitDirection,@emissionAll,pRayBackDirection)
    end if
  end if
  return radiance
end function

function RayTracerRadiance(byval pRayTracer    as const RayTracer ptr, _
                           byval pRayOrigin    as const realVector3 ptr, _
                           byval pRayDirection as const realVector3 ptr, _
                           byval lastHit       as const any ptr) as realVector3
  dim as realVector3 radiance=any
  dim as realVector3 rayBackDirection =any ' realVector3Negative( pRayDirection )
  rayBackDirection.xyz(0)=-pRayDirection->xyz(0)
  rayBackDirection.xyz(1)=-pRayDirection->xyz(1)
  rayBackDirection.xyz(2)=-pRayDirection->xyz(2)
  ' intersect ray with scene
  dim as const Triangle ptr pHitObject = 0
  dim as realVector3 hitPosition=any
  AABBNodeIntersection(pRayTracer->pScene->pRoot,pRayOrigin,pRayDirection,lastHit,0,@pHitObject,@hitPosition)
     'SceneIntersection(pRayTracer->pScene       ,pRayOrigin,pRayDirection,lastHit,  @pHitObject,@hitPosition)
  if( pHitObject ) then
    ' make surface point of intersection
    dim as SurfacePoint surfaceHitPoint = SurfacePointCreate(pHitObject,@hitPosition)
    ' local emission (only for first-hit)
    dim as realVector3 localEmission =any
    if (lastHit) then
      localEmission = realVector3ZERO
    else
      localEmission = SurfacePointEmission(@surfaceHitPoint,pRayOrigin,@rayBackDirection,FALSE)
    end if
    ' emitter sample
    dim as realVector3 emitterSample = sampleEmitters(pRayTracer,@rayBackDirection,@surfaceHitPoint)
    ' recursed reflection
    dim as realVector3 recursedReflection = realVector3ZERO
    scope
      dim as realVector3 nextDirection=any
      dim as realVector3 colour=any
      ' check surface reflects ray
      if SurfacePointNextDirection(@surfaceHitPoint,@rayBackDirection,@nextDirection,@colour) then
        ' !!! recurse !!!
        dim as realVector3 recursed = RayTracerRadiance(pRayTracer,@surfaceHitPoint.position,@nextDirection,surfaceHitPoint.pTriangle)
        recursedReflection.xyz(0)=recursed.xyz(0)*colour.xyz(0)
        recursedReflection.xyz(1)=recursed.xyz(1)*colour.xyz(1)
        recursedReflection.xyz(2)=recursed.xyz(2)*colour.xyz(2)
      end if
    end scope
    ' sum components
    radiance.xyz(0)=recursedReflection.xyz(0)+localEmission.xyz(0)+emitterSample.xyz(0)
    radiance.xyz(1)=recursedReflection.xyz(1)+localEmission.xyz(1)+emitterSample.xyz(1)
    radiance.xyz(2)=recursedReflection.xyz(2)+localEmission.xyz(2)+emitterSample.xyz(2)
  else
    ' no hit: default/background scene emission sky for downward ray, ground for upward ray
    if rayBackDirection.xyz(1) < 0.0 then
      radiance = pRayTracer->pScene->skyEmission
    else
      radiance = pRayTracer->pScene->groundEmission
    end if
  end if
  return radiance
end function

type Image
  as integer w
  as integer h
  as realVector3 ptr aPixels
end type

function ImageConstruct(byval pIn as FILE ptr) as Image ptr
  dim as Image ptr pImg = calloc(1,sizeof(Image))

  ' read width and height (int32)
  dim as long w,h
  if fscanf( pIn, "%i %i",@w,@h )<>2 then
    print "error: ImageConstruct() read width,height"
    beep:sleep:end
  end if

  ' condition width and height
  pImg->w = iif(w<1,1,w)
  pImg->h = iif(h<1,1,h)
  ' allocate pixels
  pImg->aPixels = calloc(pImg->w * pImg->h,sizeof(realVector3))
  return pImg
end function

sub ImageDestruct(byval pImg as Image ptr)
  ' free pixels
  free( pImg->aPixels )
  free( pImg )
end sub

sub ImageAddToPixel(byval pImg as const Image ptr,byval x as long,byval y as long,byval pRadiance as const realVector3 ptr)
  ' only inside image bounds
  if x<0        then exit sub
  if y<0        then exit sub
  if x>=pImg->w then exit sub
  if y>=pImg->h then exit sub
  dim as integer index = x + ((pImg->h - 1 - y) * pImg->w)
  pImg->aPixels[index] = realVector3Add(@pImg->aPixels[index],pRadiance)
end sub

type Camera
  ' eye definition
  as realVector3 viewPosition
  as real   viewAngle
  ' view frame
  as realVector3 viewDirection
  as realVector3 right
  as realVector3 up
end type
#define CameraEyePoint( pC ) ((pC)->viewPosition)

function CameraConstruct(byval pIn as FILE ptr) as Camera ptr
  dim as Camera ptr pCamera=calloc(1,sizeof(Camera))
  dim as realVector3 Y = type({0.0,1.0,0.0})
  dim as realVector3 Z = type({0.0,0.0,1.0})
  ' read and condition view definition
  dim as single viewAngleF = 0.0
  pCamera->viewPosition  = realVector3Read(pIn)
  pCamera->viewDirection = realVector3Read(pIn)
  if fscanf( pIn, "%g",@viewAngleF )<>1 then
    print "error: CameraCreate() read fov"
    beep:sleep:end 1
  end if
  ' clamp fov
  if viewAngleF<VIEW_ANGLE_MIN then
    viewAngleF=VIEW_ANGLE_MIN
  elseif viewAngleF>VIEW_ANGLE_MAX then
    viewAngleF=VIEW_ANGLE_MAX
  end if
  pCamera->viewDirection = realVector3Unitized(@pCamera->viewDirection)
  ' if degenerate, default to Z
  if realVector3IsZero(@pCamera->viewDirection) then pCamera->viewDirection = Z
  ' convert to radians
  pCamera->viewAngle = viewAngleF * (PI / 180.0)
  ' make other directions of view coord frame
  'make trial 'right', using viewDirection and assuming 'up' is Y
  dim as realVector3 vRight = realVector3Cross(@Y,@pCamera->viewDirection)
  pCamera->up    = Y
  pCamera->right = realVector3Unitized(@vRight)
  'check 'right' is valid i.e. viewDirection was not co-linear with 'up'
  if ( realVector3IsZero(@pCamera->right)=0) then
    ' use 'right', and make 'up' properly orthogonal
    dim as realVector3 vUp = realVector3Cross(@pCamera->viewDirection,@pCamera->right)
    pCamera->up = realVector3Unitized(@vUp)
  else
    ' else, assume a different 'up' and redo
    ' 'up' is Z if viewDirection is down, otherwise -Z
    dim as realVector3 vUp = iif(pCamera->viewDirection.xyz(1) < 0.0,Z,realVector3Negative(@Z))
    ' remake 'right'
    dim as realVector3 vRight = realVector3Cross(@vUp,@pCamera->viewDirection)
    pCamera->up    = vUp
    pCamera->right = realVector3Unitized(@vRight)
  end if
  return pCamera
end function

sub CameraDestruct(byval pCamera as Camera ptr)
  free(pCamera)
end sub

sub CameraFrame(byval pImage_o as const Image  ptr, _
                byval pCamera  as const Camera ptr, _
                byval pScene   as const Scene  ptr)
  dim as RayTracer rayTracer = RayTracerCreate(pScene)
  dim as real w = pImage_o->w
  dim as real h = pImage_o->h
  ' step through image pixels, sampling them
  dim as integer iW=pImage_o->w-1,iH=pImage_o->h-1
  dim as real tanViewX = realTAN( pCamera->viewAngle * 0.5 )
  dim as real tanViewY = tanViewX
  dim as real cx,cy,ratio = w / h
  dim as real inv2W=2*1.0/w
  dim as real inv2H=2*1.0/h
  dim as realVector3 sdv = any
  tanViewX*=ratio
  for y as integer = 0 to iH
    for x as integer = 0 to iW
      ' make sample ray direction, stratified by pixels
      ' make image plane XY displacement vector [-1,+1) coefficients, with sub-pixel jitter
      cx = (( (x + Rnd()) * inv2W) - 1.0) * tanViewX
      cy = (( (y + Rnd()) * inv2H) - 1.0) * tanViewY
      sdv.xyz(0)=pCamera->viewDirection.xyz(0)+pCamera->right.xyz(0)*cx+pCamera->up.xyz(0)*cy
      sdv.xyz(1)=pCamera->viewDirection.xyz(1)+pCamera->right.xyz(1)*cx+pCamera->up.xyz(1)*cy
      sdv.xyz(2)=pCamera->viewDirection.xyz(2)+pCamera->right.xyz(2)*cx+pCamera->up.xyz(2)*cy
      sdv=realVector3Unitized(@sdv)
      ' get radiance from RayTracer
      dim as realVector3 radiance = RayTracerRadiance(@rayTracer,@pCamera->viewPosition,@sdv,0)
      ' add radiance to image
      ImageAddToPixel(pImage_o,x,y,@radiance)
    next
  next
end sub

sub ReadMiniLightFile(byval sModelFile  as string, _
                      byval pIterations as integer ptr, _
                      byval ppImage     as Image ptr ptr, _
                      byval ppCamera    as Camera ptr ptr, _
                      byval ppScene     as const Scene ptr ptr)
  dim as zstring ptr MODEL_FORMAT_ID = @"#MiniLight"
  ' open model file
  dim as FILE ptr pModelFile = fopen( sModelFile, "r" )
  if pModelFile=0 then
    print "error: open '" & sModelFile & "'"
    beep:sleep:end
  end if
  ' check model file format identifier at start of first line
  scope
    ' read chars until a mismatch
    if fscanf( pModelFile, MODEL_FORMAT_ID )<>0 then
      print "error: read model file format identifier"
      beep:sleep:end
    end if
    dim as long p = ftell( pModelFile )
    if fscanf( pModelFile, MODEL_FORMAT_ID )=-1 then
      print "error: read model file format identifier"
      beep:sleep:end
    end if
    ' check if all chars were read */
    if strlen( MODEL_FORMAT_ID )<>p then
      print "error: read model file format identifier"
      beep:sleep:end
    end if
  end scope
  '  read and condition frame iterations (int32)
  dim as long itera
  if fscanf( pModelFile, "%i", @itera )<>1 then
    print "error: read number of iterations !"
    beep:sleep:end
  end if
  *pIterations = iif(itera<1,1,itera)
  ' create main rendering objects, from model file
  *ppImage  = ImageConstruct (pModelFile)
  *ppCamera = CameraConstruct(pModelFile)
  *ppScene  = SceneConstruct (pModelFile,@CameraEyePoint(*ppCamera))
  ' close model file
  if EOF_= fclose( pModelFile ) then
    print "error: close scene file !"
    beep:sleep:end
  end if
end sub


function calculateToneMapping(byval pImg as const Image ptr,byval divider as const real) as real
  ' calculate estimate of world-adaptation luminanceas log mean luminance of scene
  dim as integer nPixels=pImg->w * pImg->h
  dim as real sumOfLogs,invnPixels=1.0/nPixels
  for  i as integer = 0 to nPixels-1
    dim as real Y = realVector3Dot(@pImg->aPixels[i],@RGB_LUMINANCE ) *  divider
    ' clamp luminance to a perceptual minimum
    if Y<1e-4 then Y=1e-4
    sumOfLogs += log10(Y)
  next
  dim as real adaptLuminance = realPOW( 10.0, sumOfLogs *invnPixels)
  ' make scale-factor from:
  ' ratio of minimum visible differences in luminance, in display-adapted
  ' and world-adapted perception (discluding the constant that cancelled),
  ' divided by display max to yield a [0,1] range
  dim as real a = 1.219 + realPOW(DISPLAY_LUMINANCE_MAX * 0.25, 0.4 )
  dim as real b = 1.219 + realPOW(adaptLuminance, 0.4 )
  dim as real ab = a/b
  return realPOW(ab,2.5) / DISPLAY_LUMINANCE_MAX
end function

function isPowerOfTwo(byval x as integer) as integer 
  return (x<>0) andalso ((x and -x) = x)
end function

sub DrawFrame(byval pImage as const Image ptr, _
              byval itera  as integer)
  dim as real divider=iif(itera<1,1.0,1.0/itera)
  dim as real tonemapScaling = calculateToneMapping(pImage, divider)
  dim as real mapped=any
  dim as real gamma=any
  dim as integer iW=pImage->W-1
  dim as integer iH=pImage->h-1
  dim as integer index
  for y as integer=0 to iH
    for x as integer=0 to iW
      dim as realVector3 c=pImage->aPixels[index]:index+=1
      for i as integer=0 to 2
       ' tonemap
        mapped = c.xyz(i) * divider * tonemapScaling
        if mapped<0.0 then mapped=0.0
        ' gamma encode
        gamma = realPOW(mapped,GAMMA_ENCODE)*255.0
        ' quantize
        c.xyz(i) = floor(gamma+0.5)
        if c.xyz(i)>255.0 then c.xyz(i)=255.0
      next
      pset(x,y),realVector3RGB(@c)
    next
  next
end sub

function renderProgressively(byval nIterations as long, _
                             byval pImage  as const Image  ptr, _
                             byval pCamera as const Camera ptr, _
                             byval pScene  as const Scene  ptr) as integer
  ' do progressive refinement render loop
  dim as integer itera
  for itera = 0 to nIterations-1
    ' display current iteration number
    windowtitle "iteration:" & itera
    ' render a frame
    CameraFrame(pImage,pCamera,pScene)
    ' draw frame
    if itera mod 8=0 then DrawFrame(pImage,itera)
    if len(inkey()) then return itera
  next
  return nIterations
end function

'
' main
'
dim as Camera ptr pCamera     ' position, direction and field of view from scene file
dim as Image  ptr pImage      ' size from scene file
dim as Scene  ptr pScene      ' triangle list from scene file
dim as integer    nIterations ' number of iterations from scene file
dim as string scenefilename=command(1)
if len(scenefilename)=0 then
  print "usage: minilight scenefile"
  beep:sleep:end
end if

' create scene,image and camera from MiniLightFile
ReadMiniLightFile(scenefilename,@nIterations,@pImage,@pCamera,@pScene)
screenres pImage->W,pImage->H,32

' render nIterations
dim as real tStart=Timer()
dim as integer lastIteration = renderProgressively(nIterations,pImage,pCamera,pScene)
dim as real tTime=Timer()-tStart
WindowTitle "finished " & int(tTime) & " sec ..."
DrawFrame(pImage,lastIteration)
bSave "tmp.bmp",0

CameraDestruct(pCamera)
SceneDestruct(pScene)
ImageDestruct(pImage)
sleep
Last edited by D.J.Peters on Jul 17, 2016 22:43, edited 7 times in total.
D.J.Peters
Posts: 7852
Joined: May 28, 2005 3:28

Re: MiniLight 1.7 a minimal global illumination renderer.

Postby D.J.Peters » Mar 21, 2015 3:57

scene file: roomlightoff.txt
Image

scene file: roomlighton.txt
Image

MiniLight scene format:

Code: Select all

#MiniLight

128

512 512

(0.278 0.275 -0.789) (0 0 1) 40


(0.0906 0.0943 0.1151) (0.1 0.09 0.07)

(0.343 0.545 0.332) (0.213 0.545 0.227) (0.343 0.545 0.227)  (0.7 0.7 0.7) (1000 1000 1000)
(0.213 0.545 0.227) (0.343 0.545 0.332) (0.213 0.545 0.332)  (0.7 0.7 0.7) (1000 1000 1000)

(0.556 0.000 0.000) (0.006 0.000 0.559) (0.556 0.000 0.559)  (0.7 0.7 0.7) (0 0 0)
(0.006 0.000 0.559) (0.556 0.000 0.000) (0.003 0.000 0.000)  (0.7 0.7 0.7) (0 0 0)
...

required header: #MiniLight
required number of iterations integer >0 : 128
required image width height integer integer: 512 512
required camera position direction field of view (vector) (vector) float: (0.278 0.275 -0.789) (0 0 1) 40
required skyemission groundreflection (vector) (vector): (0.0906 0.0943 0.1151) (0.1 0.09 0.07)
triangles (vertex1) (vertex2) (vertex3) (facereflectioncolor) (lightemission):
(0.343 0.545 0.332) (0.213 0.545 0.227) (0.343 0.545 0.227) (0.7 0.7 0.7) (1000 1000 1000)
...
D.J.Peters
Posts: 7852
Joined: May 28, 2005 3:28

Re: MiniLight 1.7 a minimal global illumination renderer.

Postby D.J.Peters » Mar 22, 2015 22:11

I optimized the whole code see first post
one frame on old Pentium 2.6 seconds (before ~10.0).

Joshy

Pinhole camera with 1024 interations per scene.
scene: roomlightoff.txt
Image
scene: roomlighton.txt
Image
dafhi
Posts: 1275
Joined: Jun 04, 2005 9:51

Re: MiniLight 1.7 a minimal global illumination renderer.

Postby dafhi » Mar 23, 2015 1:09

"usage: .. .ext"
oyster
Posts: 183
Joined: Oct 11, 2005 10:46

Re: MiniLight 1.7 a minimal global illumination renderer.

Postby oyster » May 12, 2015 13:17

it would be nice to see it on the official page, and be compared with other lagnuages :)

Return to “Tips and Tricks”

Who is online

Users browsing this forum: No registered users and 1 guest