< H X A 7 2 4 1 > : MiniLight minimal global illumination renderer

description

MiniLight is a minimal global illumination renderer. It is primarily an exercise in simplicity. But that makes it a good base and benchmark for development and experimentation. And it just might be the neatest renderer around (on average, about 570 lines). There are translations into several programming languages.

Someone asked me how to develop a renderer in three months. It seemed not nearly enough time. But then I began to think how small such a program could be, whilst still being a complete physically based simulation with a clear, well-formed design. (A somewhat different aim from the ‘Minimal Ray Tracer Programming Contest’ of 1987. Although maybe still competitive, adjusting for code size inflation...) I couldn't resist an attempt.

It is a command-line application. It reads a simple text model file and writes a PPM image file. Supported platforms include Mac , Linux , and Windows .

The rendering features are:

Monte-carlo path-tracing transport
Emitter sampling
Progressive refinement
RGB light
Diffuse materials
Triangle modelling primitives
Octree spatial index
Pin-hole ‘lens’
Ward linear tone-mapping

Full source code is provided under the (new) BSD license.

images

A Cornell Box (of course) throughout the day:

Cornell Box, at morning (24 KB)

Cornell Box, at noon (23 KB)

Cornell Box, at evening (26 KB)

Cornell Box, at night (23 KB)

The models are physically based too. For example, at noon: The sky has a colour temperature of 10000K and a luminance of 1 x 10⁴ cd/m². The sun has a colour temperature of 5400K and a luminance of 1 x 10⁹ cd/m², a diameter of 1.392 x 10⁹ m, and distance of 149.6 x 10⁹ m. (It is just a square though.)

Some of the Smits test scenes:

Smits tests (32 KB)

A room, at noon-ish, and night (made in HDR and tonemapped with p3tonemapper):

Front room, at noon (97 KB)

Front room, at night (47 KB)

All these renders are 1000 paths per pixel. The standard Cornell Box, at 500 x 500 pixels, should take about 15 minutes on a 2007 machine.

comparison

There are now nine translations:

language (and version)	size (lines of code)	size (relative to C++)	speed (relative to C++)
Scala 2.7.7	431	0.45	1 / 2.7
OCaml 3.11.0	457	0.48	1 / 2.3
Python 2.5.1	490	0.52	1 / 180
Python Shed Skin 0.1.1	496	0.52	1 / 1.8
Ruby 1.8.6	498	0.52	1 / 218
Lua 5.1.4	568	0.60	1 / 52
C# 3.5	630	0.66	1 / (1–5 ?)
Flex 2/AS3 (Flash 10,0,22,87)	644	0.68	1 / 55
C++ ISO-98 (LLVM-G++-4.2)	952	1.00	1
C ISO-90 (LLVM-GCC-4.2)	1197	1.26	2.1

Line counts disclude comments and blanks and lines with only single braces/parentheses. Flex/ActionScript3 is 864 lines including extra UI.

Speed measurements were done on a Core2-Duo Mac, using the cornellbox and roomfront (non-sun) models. They will vary with platform and model.

C was much faster than C++ because LLVM's link-time optimisation worked for C, but not for C++ (otherwise, C would likely be about 10% faster). But against LLVM, MS-VC no longer has the lead.

Scala, Lua, Ruby and Flex are twice as slow without some manual optimization (Python too, probably). This is mostly inlining some vector operations (but also a little loop unrolling, and flattening the image pixel array).

Note: The different language download packages each have models that generate different sized images – calculate speed as (image width * image height * iterations) / time.

Try the Flex version here . . .

algorithm

path tracing

The simplest method of global illumination is pure monte-carlo path tracing. For each pixel, a ray is sent into the scene. When it hits something it notes its color and emission. It then bounces off and continues. Eventually all the noted info is sent back, step by step. Each step scales the light by certain factors until a final value is returned to the pixel. Repeat the whole process many times with some randomness added and average all the results. An image slowly emerges.

This strategy is essentially a blind search for light emitters. Unless a ray hits an emitter it will return only blackness. For outdoor scenes most rays will bounce into the sky, and find brightness. But for indoors the emitters are small and many many rays are needed to find them. It is too slow.

emitter sampling

This can be alleviated with emitter sampling. The location of all emitters is known from the model, and this can be used. Each step sends an extra ray to an emitter. But now there are two rays gathering the incoming light they must not duplicate any. To ensure this, ‘local’ emission at each bounce is no longer included.

monte carlo

Each interaction along the main path is handled mostly by monte-carlo. The bounce direction is chosen by importance sampling: proportionally more likely in directions the surface allows more reflection through. For ideal diffuse that is lots near perpendicular, little near grazing. The overall reflectivity is controlled by russian-roulette: proportionally more likely to bounce the more reflective the surface is. This allows a non-arbitrary end for the recursion. The colour is just scaling, because there is no advantage in monte-carlo there.

Both uses of monte-carlo sample more where the surface interaction transfers more light. This is a fairly good approach.

Why? A simple analysis is that apparent noisiness is the product of detail/richness and intensity: greater detail means adjacent pixels more likely see different values, greater intensity probably means more contrast. The latter is easiest to address, because interaction's effect on intensity is simple. The treatment—sampling more—reduces noise by averaging-out the effect of the randomness.

Over many iterations this puts more computation into resolving more visibly noisy regions, rather than wasting it on regions that won't noticeably improve. That computes an image of particular quality more quickly. But there are still shortcomings . . .

rendering equation

Here is the specific rendering equation for this algorithm: $L_o(x, \vec w) = \int_\Omega L_i(x, \vec w') (\vec w' \cdot \vec n) \frac{\rho}{\pi}\, d\vec w' \ \ + \ \ \sum_E\int_A L_e(z, \vec w'') \Gamma_A V(x,z) (\vec w'' \cdot \vec n) \frac{\rho}{\pi}\, d a$ Meaning: The radiance outward at a position and direction equals the sum of two terms (the first representing the main path, the second the emitter sampling). The first term integrates over the surrounding hemisphere of directions the inward radiance scaled by its projection and the reflectivity function. The second term sums over all emitters, for each integrating over the area the emission scaled by the solid angle, the visibility, the projection, and the reflectivity function.

unbiased?

Monte-carlo path-tracing is described as ‘unbiased’. The term is now sometimes used in the manner of a marketing tag. Its real meaning is too abstract for that purpose. It denotes a lack of consistent error—the only error is in the effects of the randomness. Such niceties are appropriate for mathematics but for engineering we must live by practicalities. All implementations are biased and incorrect in various ways. The task is to arrange them acceptably.

ultimately

The value of this method is its simplicity, while having full generality. It allows the implementation to be almost proven correct by inspection. (Correct meaning ‘following the rules for a particular standard approximation’.) The images produced can be authoritative references for other renderers.

code

These five functions from two modules, in C, are the core: the light-transport and light-interaction code. They pretty much summarise the algorithm.

The RayTracer module provides a context and services for transport. The SurfacePoint module treats geometric and quality data for a particular interaction.

The pixel sampler calls the RayTracer function RayTracerRadiance. RayTracerRadiance then recursively traces the main ray path, calling sampleEmitters to get illumination from emitters. Both RayTracer functions call the SurfacePoint functions.

RayTracer

Vector3f RayTracerRadiance
(
   const RayTracer* pR,
   const Vector3f*  pRayOrigin,
   const Vector3f*  pRayDirection,
   RandomMwc*       pRandom,
   const void*      lastHit
)
{
   Vector3f radiance;

   const Vector3f rayBackDirection = Vector3fNegative( pRayDirection );

   /* intersect ray with scene */
   const Triangle* pHitObject = 0;
   Vector3f        hitPosition;
   SceneIntersection( pR->pScene, pRayOrigin, pRayDirection, lastHit,
      &pHitObject, &hitPosition );

   if( pHitObject )
   {
      /* make surface point of intersection */
      const SurfacePoint surfacePoint = SurfacePointCreate( pHitObject,
         &hitPosition );

      /* local emission (only for first-hit) */
      const Vector3f localEmission = lastHit ? Vector3fZERO :
         SurfacePointEmission( &surfacePoint, pRayOrigin, &rayBackDirection,
         false );

      /* emitter sample */
      const Vector3f emitterSample = sampleEmitters( pR, &rayBackDirection,
         &surfacePoint, pRandom );

      /* recursed reflection */
      Vector3f recursedReflection = Vector3fZERO;
      {
         /* single hemisphere sample, ideal diffuse BRDF:
            reflected = (inradiance * pi) * (cos(in) / pi * color) * reflectance
            -- reflectance magnitude is 'scaled' by the russian roulette, cos is
            importance sampled (both done by SurfacePoint), and the pi and 1/pi
            cancel out -- leaving just: inradiance * reflectance color */
         Vector3f nextDirection;
         Vector3f color;
         /* check surface bounces ray */
         if( SurfacePointNextDirection( &surfacePoint, pRandom,
            &rayBackDirection, &nextDirection, &color ) )
         {
            /* recurse */
            const Vector3f recursed = RayTracerRadiance( pR,
               &surfacePoint.position, &nextDirection, pRandom,
               SurfacePointHitId( &surfacePoint ) );
            recursedReflection = Vector3fMulV( &recursed, &color );
         }
      }

      /* sum components */
      radiance = Vector3fAdd( &localEmission, &emitterSample );
      radiance = Vector3fAdd( &recursedReflection, &radiance );
   }
   else
   {
      /* no hit: default/background scene emission */
      radiance = SceneDefaultEmission( pR->pScene, &rayBackDirection );
   }

   return radiance;
}


static Vector3f sampleEmitters
(
   const RayTracer*    pR,
   const Vector3f*     pRayBackDirection,
   const SurfacePoint* pSurfacePoint,
   RandomMwc*          pRandom
)
{
   Vector3f radiance = Vector3fZERO;

   /* single emitter sample, ideal diffuse BRDF:
      reflected = (emitivity * solidangle) * (emitterscount) *
      (cos(emitdirection) / pi * reflectivity)
      -- SurfacePoint does the first and last parts (in separate methods) */

   /* get position on an emitter */
   Vector3f        emitterPosition;
   const Triangle* emitterId = 0;
   SceneEmitter( pR->pScene, pRandom, &emitterPosition, &emitterId );

   /* check an emitter was found */
   if( emitterId )
   {
      /* make direction to emit point */
      const Vector3f emitVector    = Vector3fSub( &emitterPosition,
         &pSurfacePoint->position );
      const Vector3f emitDirection = Vector3fUnitized( &emitVector );

      /* send shadow ray */
      const Triangle* pHitObject = 0;
      Vector3f        hitPosition;
      SceneIntersection( pR->pScene, &pSurfacePoint->position, &emitDirection,
         SurfacePointHitId( pSurfacePoint ), &pHitObject, &hitPosition );

      /* check if unshadowed */
      if( !pHitObject | (emitterId == pHitObject) )
      {
         /* get inward emission value */
         const SurfacePoint sp            = SurfacePointCreate( emitterId,
            &emitterPosition );
         const Vector3f backEmitDirection = Vector3fNegative( &emitDirection );
         const Vector3f emissionIn        = SurfacePointEmission( &sp,
            &pSurfacePoint->position, &backEmitDirection, true );
         const Vector3f emissionAll       = Vector3fMulF( &emissionIn,
            (float)SceneEmittersCount( pR->pScene ) );

         /* get amount reflected by surface */
         radiance = SurfacePointReflection( pSurfacePoint, &emitDirection,
            &emissionAll, pRayBackDirection );
      }
   }

   return radiance;
}

SurfacePoint

Vector3f SurfacePointEmission
(
   const SurfacePoint* pS,
   const Vector3f*     pToPosition,
   const Vector3f*     pOutDirection,
   bool                isSolidAngle
)
{
   const Vector3f ray       = Vector3fSub( pToPosition, &pS->position );
   const float    distance2 = Vector3fDot( &ray, &ray );
   const Vector3f normal    = TriangleNormal( pS->pTriangle );
   const float    cosOut    = Vector3fDot( pOutDirection, &normal );
   const float    area      = TriangleArea( pS->pTriangle );

   /* emit from front face of surface only */
   const float solidAngle = (float)(cosOut > 0.0f) * (isSolidAngle ?
      /* clamp-out infinity */
      (cosOut * area) / (distance2 >= 1e-6f ? distance2 : 1e-6f) : 1.0f);

   return Vector3fMulF( &pS->pTriangle->emitivity, solidAngle );
}


Vector3f SurfacePointReflection
(
   const SurfacePoint* pS,
   const Vector3f*     pInDirection,
   const Vector3f*     pInRadiance,
   const Vector3f*     pOutDirection
)
{
   const Vector3f normal = TriangleNormal( pS->pTriangle );
   const float    inDot  = Vector3fDot( pInDirection,  &normal );
   const float    outDot = Vector3fDot( pOutDirection, &normal );

   /* directions must be on same side of surface (no transmission) */
   const bool isSameSide = !( (inDot < 0.0f) ^ (outDot < 0.0f) );

   /* ideal diffuse BRDF:
      radiance scaled by reflectivity, cosine, and 1/pi  */
   const Vector3f r = Vector3fMulV( pInRadiance, &pS->pTriangle->reflectivity );
   return Vector3fMulF( &r, ((float)fabs( inDot ) / PI) * (float)isSameSide );
}


bool SurfacePointNextDirection
(
   const SurfacePoint* pS,
   RandomMwc*          pRandom,
   const Vector3f*     pInDirection,
   Vector3f*           pOutDirection_o,
   Vector3f*           pColor_o
)
{
   const float reflectivityMean =
      Vector3fDot( &pS->pTriangle->reflectivity, &Vector3fONE ) / 3.0f;

   /* russian-roulette for reflectance 'magnitude' */
   const bool isAlive = RandomMwcFloat( pRandom ) < reflectivityMean;

   if( isAlive )
   {
      /* cosine-weighted importance sample hemisphere */

      const float _2pr1 = PI * 2.0f * RandomMwcFloat( pRandom );
      const float sr2   = (float)sqrt( RandomMwcFloat( pRandom ) );

      /* make coord frame coefficients (z in normal direction) */
      const float x = (float)cos( _2pr1 ) * sr2;
      const float y = (float)sin( _2pr1 ) * sr2;
      const float z = (float)sqrt( 1.0f - (sr2 * sr2) );

      /* make coord frame */
      const Vector3f t = TriangleTangent( pS->pTriangle );
      Vector3f       n = TriangleNormal( pS->pTriangle );
      Vector3f       c;
      /* put normal on inward ray side of surface (no transmission) */
      if( Vector3fDot( &n, pInDirection ) < 0.0f )
      {
         n = Vector3fNegative( &n );
      }
      c = Vector3fCross( &n, &t );

      {
         /* scale frame by coefficients */
         const Vector3f tx = Vector3fMulF( &t, x );
         const Vector3f cy = Vector3fMulF( &c, y );
         const Vector3f nz = Vector3fMulF( &n, z );

         /* make direction from sum of scaled components */
         const Vector3f sum = Vector3fAdd( &tx, &cy );
         *pOutDirection_o = Vector3fAdd( &sum, &nz );
      }

      /* make color by dividing-out mean from reflectivity */
      *pColor_o = Vector3fMulF( &pS->pTriangle->reflectivity,
         1.0f / reflectivityMean );
   }

   return isAlive && !Vector3fIsZero( pOutDirection_o );
}

downloads

Version 1.5.2 updated 2008-02-17:
(C++ and OCaml 1.5.2.1 updated 2009-04-01)
(C 1.5.2 added 2009-06-20)
(Scala 1.5.2.1 updated 2010-01-09)

minilight15c.zip	C ISO-90 source code and Mac, Linux, and Windows executables (154 KB)
minilight15cpp.zip	C++ ISO-98 source code (205 KB)
minilight15ocaml.zip	OCaml 3.11 source code (30 KB)
minilight15scala.zip	Scala 2.7 source code (32 KB)
minilight15lua.zip	Lua 5.1 source code (28 KB)
minilight15python.zip	Python 2.5 source code (18 KB) (Translation by Juraj Sukop)
minilight15ruby.zip	Ruby 1.8 source code (25 KB)
minilight15flex.zip	Flex 2 source code (44 KB)

external

MiniLight OCaml	adjusted and varied by Mauricio Fernández
MiniLight C# 3.5	translated (with small extras) by Chris Lomont
MiniLight Clojure 2009	in progress, by Mark Reid

contribution

If you can't resist doing a translation, I can probably add/link to it. It is a good way to learn 80% of a new language. The current ‘rules’ are:

Follow the main module separation design.
Make it compact and simple (these are the priority), by good use of the language.
Only optimize where really needed, and not excessively (you want to find out how well the compiler/language can do for itself with no manual help).
Write the best code possible!

Sacrifice some clarity for compactness in the SpatialIndex module. It is a bit large (it deserves careful testing, too). And, only to make the sizes comparable:

lines must be <= 80 chars
indent ought to be 3-ish chars (but mixing tabs and spaces is instant disqualification!)

The reference version is the C, for functionality and logic. Otherwise, the OCaml, for a more functional form.

hxa7241+www (ατ) googlem‌ail (dοτ) com

MINILIGHT