usage : minilight scenefile
Joshy
scene file: cornellbox.txt
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
