Universität Karlsruhe (TH)Forschungsuniversität - gegründet 1825
Fakultät für Informatik
Institut für Betriebs- und Dialogsysteme

Studienarbeit

Accelerating the bidirectional path tracing
algorithm using a dedicated intersection processor
Wenzel Jakob
30.07.2007

Supervised by Prof. em. Alfred Schmitt
Dipl. inform. Stefan PreußAbstract
Even with the processing power of today’s computers, the calculation of global illumination solutions for general scenes remains a lengthy undertaking. This project investigates
the possibility of accelerating such calculations using a dedicated hardware co-processor.
The algorithm of choice, bidirectional path tracing, presents an additional challenge, since
its three-dimensional random walks make the use of coherent ray tracing techniques impracticable. A parallel and deeply pipelined ray-triangle intersection processor including a
kd-tree traversal unit was developed during the course of the project, and numerous benchmarks were performed.

Contents
List of Figures

7

List of Tables

8

Listings

9

1 Overview
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Architecture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 License . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11
11
11
13

2 Global Illumination
2.1 Solid angles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Radiometric quantities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Irradiance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 Radiance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Geometric term . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4 Bidirectional reflectance distribution function . . . . . . . . . . . . . . . . . .
2.4.1 Modified Blinn-Phong BRDF . . . . . . . . . . . . . . . . . . . . . . . .
2.5 Rendering equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5.1 Area formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5.2 Path integral formulation . . . . . . . . . . . . . . . . . . . . . . . . . .

15
15
15
16
16
17
18
19
19
20
20

3 Monte Carlo Integration
3.1 Basic Monte Carlo integration . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.1 Primary estimator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.2 Secondary estimator . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Variance reduction techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Importance Sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 Stratified sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Russian roulette . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 Estimators for Fredholm equations . . . . . . . . . . . . . . . . . . . . . . . . .
3.5 Sampling the hemisphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23
23
23
24
24
24
25
25
26
26

4 Bidirectional path tracing
4.1 Path tracing and light tracing . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29
29

4

C
4.2
4.3
4.4
4.5

C

A bidirectional estimator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Pure software implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The bipartite mutual visibility query . . . . . . . . . . . . . . . . . . . . . . . .

29
31
37
39

5 Target architecture & design flow
5.1 Logic resources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.1 Look-up tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.2 Interconnect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.3 Clocking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.4 DSP48 tiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.5 Block RAM tiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Design flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.1 Logic design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.2 Synthesis and mapping . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.3 Placement and routing . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41
41
43
43
43
44
44
44
45
46
46

6 Building blocks
6.1 Larger multipliers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2 Non-restoring integer divider . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3 Floating-point unit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3.1 Floating-point adder/subtracter . . . . . . . . . . . . . . . . . . . . . .
6.3.2 Floating-point multiplier/divider units . . . . . . . . . . . . . . . . . .
6.4 Memory hierarchy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.4.1 DDR SDRAM Controller . . . . . . . . . . . . . . . . . . . . . . . . . .
6.4.2 Caching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.5 USB Interface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.5.1 Device driver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49
49
50
51
51
52
53
53
54
54
54

7 Intersection unit
7.1 Triangle intersection unit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.1.1 Intersection algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.1.2 Hardware implementation . . . . . . . . . . . . . . . . . . . . . . . . .
7.2 kd-tree intersection unit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2.1 Intersection algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2.2 Generating kd-trees . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2.3 Hardware implementation . . . . . . . . . . . . . . . . . . . . . . . . .
7.3 Bidirectional path-tracing frontend . . . . . . . . . . . . . . . . . . . . . . . . .

57
57
57
58
60
60
61
61
65

8 Software
8.1 Software pipeline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67
67

9 Results
9.1 Projected performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9.2 Conclusion and future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69
71
71

5

C
Bibliography

6C
73

List of Figures
1.1
1.2

Architecture overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
FPGA architecture overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12
13

2.1
2.2
2.3
2.4
2.5

A graphical illustration of radiance . . . . . . . . . . . . . . . . . . . . . . . . .
Illustration of the geometric term . . . . . . . . . . . . . . . . . . . . . . . . . .
The Bidirectional Reflectance Distribution Function . . . . . . . . . . . . . .
Polar plots of the modified Blinn-Phong BRDF . . . . . . . . . . . . . . . . . .
An illustration of the operators K and G [1]. . . . . . . . . . . . . . . . . . . .

17
18
19
19
21

4.1
4.2
4.3

Overview of the bidirectional path tracing algorithm . . . . . . . . . . . . . .
Cornell box rendering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Cornell box rendering grouped by path configuration. Each image contains
the paths generated by one of the partial estimators – “i ” and “ j ” denotethe eye and light path steps, respectively. . . . . . . . . . . . . . . . . . . . . .

30
37

5.1
5.2
5.3
5.4
5.5
5.6
5.7

A generic FPGA array architecture . . . . . . . . . . . . . . . . . . . . . . . . .
Layout of a Configurable Logic Block . . . . . . . . . . . . . . . . . . . . . . . .
Simplified diagram of a Virtex-4 slice . . . . . . . . . . . . . . . . . . . . . . . .
Hierarchical interconnect resources . . . . . . . . . . . . . . . . . . . . . . . . .
DSP48 slice architecture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Diagram of the state machine . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Synthesized parity calculation circuit . . . . . . . . . . . . . . . . . . . . . . . .

41
42
42
43
45
45
46

6.1
6.2
6.3
6.4

24-bit floating-point format . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Block diagram of the floating-point adder/subtracter . . . . . . . . . . . . . .
Block diagram of the floating-point multiplier . . . . . . . . . . . . . . . . . .
Implementation hardware . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51
52
53
55

7.1
7.2
7.3

Ray-triangle intersection data path . . . . . . . . . . . . . . . . . . . . . . . . .
kd-tree traversal data path . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
kd-tree control unit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

59
62
62

9.1

Renderings of the used test scenes . . . . . . . . . . . . . . . . . . . . . . . . . .

70

38

7

List of Tables

8

5.1

Logic resources available in current Virtex-4 devices . . . . . . . . . . . . . . .

45

6.1

Implementation costs of the floating-point cores . . . . . . . . . . . . . . . . .

53

7.1
7.2

Profiler report excerpt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Implementation costs of the intersection units . . . . . . . . . . . . . . . . . .

57
65

9.1
9.2

Average number of accesses per intersection operation . . . . . . . . . . . . .
Intersections per second using different core configurations . . . . . . . . . .

71
71

Listings
3.1

An unbiased primary estimator for Fredholm equations . . . . . . . . . . . .

27

4.1
4.2
4.3
4.4
4.5
4.6

Surface interaction data structure . . . . . . . . . . . . . . . . . . . . . . . . . .
Main rendering loop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Random walk function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Path evaluation function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Path rendering function (part 1) . . . . . . . . . . . . . . . . . . . . . . . . . .
Path rendering function (part 2) . . . . . . . . . . . . . . . . . . . . . . . . . .

32
32
33
34
35
36

5.1
5.2

Entity declaration of the finite state machine . . . . . . . . . . . . . . . . . . .
Architecture declaration of the finite state machine . . . . . . . . . . . . . . .

46
47

6.1
6.2

Restoring division algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Non-restoring division algorithm . . . . . . . . . . . . . . . . . . . . . . . . . .

50
51

7.1
7.2
7.3
7.4

Ray-triangle intersection pseudocode . . . . . . . . . . . . . . . . . . . . . . .
kd-tree traversal pseudocode . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Control unit – stack data type . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Control unit – pipeline stage data type . . . . . . . . . . . . . . . . . . . . . . .

58
60
63
63

8.1

A scene using beam’s XML-based scene description language . . . . . . . . .

68

9L

10

L

1 Overview
1.1 Motivation
The main motivation of this project is the inherent parallelism of global illumination algorithms. Large numbers of rays need to be shot through virtual scenes – an operation,
which seems very appropriate for acceleration using dedicated hardware.
On the other hand, hardware development platforms capable of prototyping such designs have only recently become available. However, successes by other groups, such as thegraphics group of the University of Saarland, encourage further investigation. While mostof the current focus is directed towards real-time Whitted-style ray tracing, little has been
done in the direction of global illumination acceleration using hardware. One reason forthis might be the sheer computational complexity of these algorithms in comparison toplain ray tracing. This project investigates the possibility of accelerating these calculations
using a dedicated hardware co-processor. The used algorithm, bidirectional path tracing,poses an additional problem, since its random walks make the use of any coherent ray
tracing techniques1 impractical. To achieve a significant speedup despite this drawback,the algorithm was analyzed, and a new bipartite mutual visibility query instruction was
implemented in silicon. Each such instruction consists of several intersection tests using
hierarchical data structures. The final hardware implementation is able to execute up to 96
of these intersections tests in parallel.

1.2 Architecture
The project consists of three major parts: A global illumination renderer called beam (Discussed in Chapter 8), a Linux kernel driver, and the FPGA2 processor core. The scope of
this project was to plan and develop these components. Finally, the intersection architecture was to be implemented on actual hardware.
The FPGA core can be divided into five major blocks: The USB host interface, a DDR
memory controller, a multi-ported direct-mapped cache, a ray scheduler, and several combined intersection and kd-tree traversal units. Figure 1.2 contains all FPGA components
and their source code locations.
When beam renders an image, the following operations occur in sequence:
1

Coherent ray tracing techniques exploit the statistical dependence of rays to trace several rays simultaneously with negligible overhead
2
Field Programmable Gate Array - A reconfigurable microchip consisting of programmable logic and programmable interconnect.

111.2 A

1 O

USB 2.0

beam
Global Illumination rendererintersect kernel driverLinux 2.6 kernel

FPGA

Ray scheduler
Parallel intersection units
Cache&DDR Memory Controller

Figure 1.1: Architecture overview(i) beam parses the XML scene description and loads all referenced plugins such as material types, light sources and sampling techniques.
(ii) A kd-tree acceleration data structure is recursively constructed from the scene’s geometry using heuristical methods.
(iii) The renderer converts this kd-tree data structure into a memory image that is compatible with the FPGA circuitry. This involves number format conversion and transforming triangle data into a special redundant storage format.
(iv) This memory image is uploaded to a DDR SDRAM chip connected to the FPGA.
Afterwards, the FPGA cache and intersection units are initialized.
(v) The renderer starts tracing rays using the intersection processors. Multithreading
support is provided by the hardware in order to process packets of rays at the sametime, which helps to reduce the effects of memory fetch and pipeline latencies. Furthermore, the beam renderer uses a software pipeline to keep both the FPGA and thecomputer’s CPU busy at all times.
(vi) The resulting bitmap consisting of spectral samples is developed into a RGB-basedformat. Instead of writing one output image, beam creates several partial images andone combined image to visualize the inner workings of the bidirectional path tracing
algorithm.
The following chapters introduce the relevant foundations in the fields of Global Illumination, Monte Carlo integration, bidirectional path tracing and FPGA architecture. The
processor implementation is discussed in chapters 6 and 7 and the software implementation is explained in chapter 8. Chapter 9 concludes with benchmark results and ideas forfurther work.

12

1 O

1.3 L
Host System

32 MB DDR
SDRAM
100MhzCypress USB
2.0 PHY

3,2 Gb/s

40Mhz

Direct-mapped cachestate machine &round-robin cacherequest scheduler

SSTL I/O buffers,
data recapture &delay control

Communication statemachine & 512-deepcross clock-domain
FIFOs

rcache.vhd

ddr_pads.vhd
Open row manager,command scheduler &initialization statemachineusb_iface.vhd
usb_iface_pkg.vhd

Multi-ported cachememory(2 × 1024 × 398 bits)

ddr_memctrl.vhd

Organic LED display
communicationinterfacercache_mpram.vhd

ddr_wrapper.vhd

cache_pkg.vhd

OSRAM
Organic LED
display (Status
messages)oled_ctrl.vhd
oled_ctrl_pkg.vhd
Command interpreter
finite state machine

Ray scheduler

Fetch stage
kd-treetraversal
pipeline(19 stages)

Ray scheduler &Mutual visibility query
supportTriangleintersection
pipeline(30 stages)

intersect_bd.vhd

Forward
Decision stage

Internal clock signal
generation,
USB+DDR clock
signal generation,
Power-on initialization

Multi-stack(32×32-deep)

intersect_main.vhd

200 Mhzdifferential
clock source

clockmgr_virtex4.vhd

Figure 1.2: FPGA architecture overview

1.3 License
All source code, including the hardware description code, is available under an open source
license and may be distributed accordingly. This document is placed into the public domain.

13

1.3 L

14

1 O

2 Global Illumination
Global illumination (GI) tries to visualize a three-dimensional scene in a way that the appearance of a single object can be influenced by any other object present in the scene. Thisbehavior makes it possible to account for complex illumination features such as indirect
lighting or caustics. Images rendered using global illumination look amazingly photorealistic – but they also are very costly to create in terms of required computational power.
This chapter discusses some of the theoretical foundations of global illumination.

2.1 Solid angles
A sound theoretical description of global illumination algorithms requires the integrationof functions over the hemisphere above a point. Thus, a measure for the contents of ahemispherical area is needed in order to define a meaningful integral. This is called asolid angle and is defined as the surface area on the hemisphere divided by the squared
hemisphere radius (usually r = 1). Its unit is named steradian (sr).
Using this measure and a transformation to spherical coordinates, a solid angle is represented as Θ = (ϕ, θ ), with ϕ and θ being the azimuth and elevation, respectively. We thendefine a differential solid angledωΘ := sin(θ ) dθ dϕ(2.1)

where dθ and dϕ are the differential elevation and azimuth. The decrease in hemisphericalarea near the horizon is the reason for the additional sine factor. We defineπ(2.2)Ω+ := [0, 2π] × [0, ]2as the set comprising the entire hemisphere. The integral can now be expressed:Z
Z 2π Z π/2Ω+

f (Θ) dωΘ =f (ϕ, θ ) · sin(θ ) dθ dϕ

0

(2.3)

0

2.2 Radiometric quantities
Global illumination algorithms attempt to converge to a stable, physically correct distribution of light – therefore, the development of such an algorithm requires an understandingof the underlying physical quantities. Radiometry is the measurement of optical radiation,
including the visible range between 400 and 700nm. This section will not go into the complex physical details and instead only mention the parts that are important for computer15

2.2 R 

2 G I

graphics. A thorough introduction can be found in [2]. As an additional simplification,
this section will leave out the concept of wavelengths and spectral radiance.
The most fundamental radiometric quantity is called radiant power or flux and is usuallydenoted by the letter Φ. It expresses how much total energy flows through a subset of R3
per unit time and is measured in Watt (Joules per second).

2.2.1 Irradiance
The irradiance with respect to a location x denotes the radiant power per unit surface area
and is measured in Watt/m2 . It describes the density of radiant power arriving at a surface
point.
dΦE(x) :=(2.4)
d Ax
where A x is an infinitesimal surface patch containing x . The word irradiance usually signifies that only incident radiant power is considered. Otherwise, the word radiosity is used.

2.2.2 Radiance
The radiance for a given point x and a direction of interest Θ is defined as the flux per
projected area per unit solid angle. The projected area A⊥x is a differential surface area
projected perpendicularly to the direction of interest and ωΘ is a differential solid angle
around Θ. For an illustration, refer to figure 2.1. Radiance is measured in Watt/(sr · m2 ).
d2 ΦL(x, Θ) :=d A⊥x dωΘ

(2.5)

If we want to describe the radiance that is incident to or exitant1 from a given differential
surface A x with the surface normal Nx at x , the following equality holds:d2 ΦL(x, Θ) =cos(Nx , Θ) d A x dωΘ(2.6)

where cos(Nx , Θ) is the cosine of the angle between Θ and the surface normal. This result makes sense when imagining a circular shaft of light shining on a surface with a nonperpendicular angle (cos(Nx , Θ) < 1). Since the surface area receiving light resembles astretched circle, it is larger than a perpendicularly projected area receiving this light. The
cosine term in (2.6) is required to cancel out this increase in area.
From now on, we will categorize radiance by its orientation to a surface point and write
L o for exitant radiance and L i for incident radiance. In space, these functions satisfy
L i (x, Θ) = L o (x, −Θ).

(2.7)

Thus, there is no distinction between incident or outgoing radiance. However, at a surface
points, these functions measure different quantities.
1

The usage of this word stems from Eric Veach’s Ph. D. thesis [10]16

2 G I

2.3 G 
ωΘ

Θ

A⊥
x
x

α
(a) Definition of radiance(b) Non-perpendicular incident light

Figure 2.1: A graphical illustration of radiance
The irradiance at a surface point x can be obtained by integrating the radiance function
over the hemisphere above x . In (2.8), the cosine factor models the fact that incident lightat a grazing angle will deposit less radiant energy than light at a perpendicular angle.
Z
E(x) =Ω+

L i (x, Θ) cos(Nx , Θ) dωΘ

(2.8)

One important property of radiance is that, in absence of participating media such asfog or flame, its value is invariant along straight paths.

2.3 Geometric term
Let us consider the situation of figure 2.3: there are two differential surfaces A x and A y anda normalized direction vector pointing from x to y.
Θ :=

y−x
k y − xk.

Dutré et al. [2] use a transformation for expressing a solid angle ωΘ as a differential surface:dωΘ =cos(N y , −Θ) d A y
kx − yk2(2.9)

Suppose that we want to express the radiance traveling on a ray from x to y . (2.6) and (2.9)
imply thatd2 ΦL(x, Θ) =cos(Nx , Θ) d A x dωΘ
d2 Φ
1
=
·
d A x d A y G(x, y)
where
G(x, y) :=cos(Nx , Θ) cos(N y , −Θ)kx − yk2(2.10)(2.11)(2.12)

17

2.4 B   

2 G I

y
−Θ

Nx

Ax

x

Ny

Θ

Ay

Figure 2.2: Illustration of the geometric termG(x, y) is called the geometric term and describes how much radiance is passed from A x
to A y . The denominator expresses how the subtended solid angle of the differential surface
A y on the hemisphere of x will shrink quadratically with increasing distance.

2.4 Bidirectional reflectance distribution function
The bidirectional reflectance distribution function (BRDF) describes what fraction of theincident irradiance from a certain direction (Ψ) leaves as reflected radiance in anotherdirection (Φ). It was introduced by Nicodemus in 1970 and models how light interacts
with a material. Using (2.8), we can derive the irradiance function for a solid angle:
d E(x, Θ) = L i (x, Θ) cos(Θ, Nx ) dωΘ

(2.13)

Now we can write down the definition of f r : R3 × Ω+ × Ω+ → R.
f r (x, Ψ, Θ) :=d L o (x, Θ)d L o (x, Θ)=d E(x, Ψ)
L i (x, Ψ) cos(Ψ, Nx ) dωΨ

(2.14)In order to be physically plausible, the BRDF needs to have several properties:(i) Range: The value of a BRDF is always greater than or equal to zero.
(ii) Reciprocity: The BRDF value does not change when Θ and Ψ are exchanged. The
notation f r (x, Ψ ↔ Θ) denotes this property.
(iii) Energy conservation: The law of conservation of energy demands that, given incident irradiance along an angle Ψ, the total reflected radiance is less than or equal tothis amount.
ZΩ+f r (x, Ψ, Θ) cos(Nx , Θ) dωΘ ≤ 1(2.15)The bidirectional scattering distribution function (BSDF) is an extension of the BRDF to
the whole sphere around a surface point and as such, it can model refraction in addition toreflection.

18

2 G I

2.5 R 

Ψ

Θ

Nx

x
A

Figure 2.3: The Bidirectional Reflectance Distribution Function

2.4.1 Modified Blinn-Phong BRDF
Since the commonly used Phong BRDF does not satisfy the physical plausibility propertiesof energy conservation and reciprocity, a modified [5] Blinn-Phong model can be usedinstead:1
n+2
cosn (Nx , H) + kd
(2.16)f r (x, Ψ ↔ Θ) = ks
2ππ
where H is a half-vector in the plane of Ψ and Θ and kd + ks ≤ 1 (energy conservation). kd
and ks respectively control how diffuse and specular the surface is. This model raises the
cosine of the angle between H and the surface normal to a power n, which determines how
sharp the specular reflection is – a low value causes glossy reflection.

0.35

0.18
0.16

0.3

0.14
0.25

0.12

0.2

0.1

0.15

0.08
0.06

0.1
0.3
0.05

0.25
0.2

0
-0.05

0.15
0

0.05

0.1
0.1

0.05
0.15

0.2

0
0.25

0.3

-0.05

(a) Perfect specular reflection (kd = 0, n = 160)

0.04

0.2
0.15

0.02
0
-0.15

0.1
0.05
-0.1

-0.05

0
0

-0.05
0.05

0.1

-0.1
0.15

0.2

-0.15

(b) Diffuse and glossy specular reflection (kd = 0.4, n = 40)

Figure 2.4: Polar plots of the modified Blinn-Phong BRDF (ks = 0.02 in both cases). The
black arrow denotes the (constant) direction of incoming irradiance. The hemisphere radius in a direction denotes the amount of reflectance.

2.5 Rendering equation
When solving a global illumination problem, we are most interested in a radiance function
describing the state after an energy equilibrium has been reached. In 1986, Kajiya presented

19

2.5 R 

2 G I

the rendering equation [3] as a generalization of various global illumination algorithms. It
is based on similar equations known from radiative heat transfer research. We will use aslightly modified [2] version expressed using angles and a BRDF.
The exitant radiance L o (x, Θ) can be split up into an emitted radiance term L e (emitted by the surface itself) and a reflected radiance term L r (created as a result of incident
radiance).
L o (x, Θ) = L e (x, Θ) + L r (x, Θ)(2.17)
Using the definition of radiance (2.5) and the definition of the BRDF (2.14), we can calculate L o using an intuitive integral equation:Z
L o (x, Θ) = L e (x, Θ) +Ω+

L i (x, Ψ) f r (x, Ψ, Θ) cos(Nx , Ψ) dωΨ

(2.18)

This is the rendering equation in its basic hemispherical form – being a Fredholm equationof the second kind, it cannot be analytically solved except for some absolutely trivial cases.

2.5.1 Area formulation
There are several equivalent formulations of this equation – the second important one does
not integrate over the hemisphere of a point but does so over all geometry contained in thescene. For this, we need a ray casting function:r(x, Θ) := x + Θ · min{t ∈ R+ | x + t · Θ ∈ A }where A ⊆ R3 contains all surface points of the scene. Thus, r finds the closest intersection
point between a surface and the ray starting at x in direction Θ. Using r , we can also define
a mutual visibility function V :(1, if r(x, −x→y) = yV (x, y) :=0, otherwise
Using V and the geometric term (2.12), the L i term in (2.18) can be substituted.
Z
L (x, Θ) = L (x, Θ) +L ( y, −y→x) f (x, −x→y, Θ) V (x, y) G(x, y) d A
o

e

o

r

A

2.5.2 Path integral formulation
The third approach [1] introduces three new integral operators K , G and T :Z(Kh)(x, Θ) :=h(x, Ψ) f r (x, Ψ, Θ) cos(Nx , Ψ) dωΨ
Ω+

20

y

(2.19)

2 G I

2.5 R 

K is called the local scattering operator because it matches the reflectance integral L r of therendering equation.
(h(r(x, Θ), −Θ), if r(x, Θ) is defined(Gh)(x, Θ) :=0,otherwise
G is called the propagation operator, since it converts exitant radiance on distant surfaces toincident radiance on the local surface such that G L o = L i . For a visual illustration of theoperators, refer to figure 2.5.

G

K

r

r

r

Figure 2.5: An illustration of the operators K and G [1].
T := K ◦ G

The third operator T is the concatenation of the previous two operators and is called thelight transport operator. Using T , the rendering equation simplifies to
Lo = Le + K Li⇔

Lo = Le + T Lo⇔

L e = (I − T )L o⇔

L o = (I − T )−1 L e(2.20)

assuming that (I − T ) has an inverse. Similarly to the geometric series, the inverse exists if
kT k < 1, where k · k is the operator norm given by:
kSk := sup kS f k∼
k f k∼ ≤1

for some norm k · k∼ . In this case, the domain of this norm is the set of all radiancefunctions, which becomes a banach space using the L p norm defined by:k f kp =

Z Z
R3

1
p

| f (x, Θ)| p dωΘ d xΩ

In 1995, Arvo [1] has shown for one-sided reflective surfaces with physically plausible
BRDFs that kGk ≤ 1 and kKk ≤ 1 if k · k∼ is a L p -norm with 1 ≤ p ≤ ∞. By futherassuming that no surface in the scene is perfectly reflective, he has furthermore shown thatkKk < 1 and thuskT k = kK Gk ≤ kKk kGk < 1(2.21)

21

2.5 R 

2 G I

Therefore, (I − T ) is invertible and its inverse is the Neumann series of T
(I − T )−1= I + T + T + ··· =
2

∞
X

Tn

(2.22)Lo = Le + T Le + T 2 Le + T 3 Le + · · ·(2.23)

n=0

Combining (2.20) and (2.22) yields

In other words, the rendering equation can also be expressed as the sum of contributions
by paths of lengths 1, 2, etc. coming from an emissive scene object. This is the integral we
will try to solve.

22

3 Monte Carlo Integration
Monte Carlo Integration tries to solve integrals which might otherwise be hard or impossible to solve using analytical methods. Consider the one-dimensional definite integral overthe interval [0, 1].
Z
1

I :=

f (x) d x.

(3.1)

0

Many functions from the domain of computer graphics do not possess elementary antiderivatives. In these cases, numerical quadrature methods are usually applied to approximate I . However, these methods require the evaluation of f at many sample points.
Furthermore, when integrals of higher dimension are to be solved, quadrature methods
become increasingly impractical since they require the evaluation of f over a multidimensional grid with a polynomial amount of sample points. This phenomenon is known as
the curse of dimensionality. In these cases, Monte Carlo Integration may still provide very
accurate numerical solutions while not suffering from any of these problems. Monte Carlo
Integration is covered in depth in [2], [4] and [10].

3.1 Basic Monte Carlo integration
One important quality of Monte Carlo methods is that their convergence rates are completely independent of the integration domain’s dimensionality. This property makes themthe method of choice for solving integrals like (2.23), which require the integration over an
infinite amount of dimensions. For simplicity, the following estimators are restricted to
one dimension. However, they could easily be extended to multiple dimensions withoutchanging their fundamental properties.

3.1.1 Primary estimator
Given n independent and identically-distributed random variables X 1 , . . . , X n ∼ U(0, 1)
with U being the standard uniform distribution, we can define n primary estimators〈I i 〉 := f (X i )In this chapter, pX i denotes the probability distribution function of X i . The expected valueof these estimators is the value of the definite integral (3.1):
E[〈I i 〉] =Z10f (x) pX i (x) d x =Z1

f (x) d x = I

(3.2)

0

23

3.2 V  

3 M C I

The estimator is unbiased because the expected values matches the quantity which is to be
estimated. The variance σ〈I2 i 〉 of 〈I i 〉 can also be calculated:
σ〈I2 i 〉 = E[( f (x) − I)2 ] =Z1

f 2 (x) d x − I 2

(3.3)

0

In practice, one would never use this estimator by itself because the variance is constantand generally too high to yield acceptable results.

3.1.2 Secondary estimator
A secondary estimator can be created by averaging the contributions of the primary estimators. An actual implementation generates a series of pseudo-random numbers andaverages the function evaluations at these points.
〈J〉 :=n
1X〈I i 〉n i=0

The expected value of this estimator also matches the integral (3.1) and it is thus unbiased.
n
n
1X
1X
E[〈J〉] = E[
〈I i 〉] =E[〈I i 〉] = I
n i=0
n i=0

(3.4)

Compared to (3.3), the variance is reduced by a factor of n.
iid

2
=
σ〈J〉

n
1 X
2

n

σ〈I2 i 〉 =

i=0

1
n

σ〈I2 0 〉(3.5)

p
However, the standard deviation σ is only reduced by a factor of n. This is also the main
practical problem with 〈J〉 and Monte-Carlo integration in general. In order to improve
the quality of the result I by a factor of two, the number of samples has to be increased
four-fold.

3.2 Variance reduction techniques
Many techniques have been developed to mitigate this property of Monte Carlo methods.
Prominent examples are importance sampling and stratified sampling.

3.2.1 Importance Sampling
The underlying idea of importance sampling is that some parts on the interval [0, 1] aremore important for the calculation of I than others, specifically the areas where f (x) has arelatively large value. Ideal random variables X 1 , . . . , X i would only sample f within theseareas while still resulting in an unbiased estimator.

24

3 M C I

3.3 R 

Suppose that X 0 . . . X n are independent and distributed according to a probability density
function p and p(x) 6= 0 (x ∈ [0, 1]). We can define new estimators〈I i0 〉 :=f (X i )andp(X i )n
1X
〈J 〉 :=〈I i 〉n i=0
0

Again, these estimators are unbiased, but they have different variance terms:
Z1
f (x)
0
p(x) d x = I
E[〈I i 〉] =
p(x)
0
Z 1
Z1 2
2
f (x)
f (x)
2
2
σ〈I 0 〉 =p(x) d x − I =dx − I2
i
p(x)
p(x)
0
0
!
Z1 2
n
n
X
X
1
f
(x)
1 2
1
iid
2
2
σ〈J
σ
d
x
−
I
0 =
0〉 =
n2 i=0 〈Ii 〉 n2 i=0
p(x)
n
0

(3.6)

(3.7)
(3.8)

2
As can be seen in (3.8), the variance σ〈J0 〉 is zero if p(x) = f (x)/I . Unfortunately, samplingaccording to this probability density function involves both normalizing p and calculating
the inverse of its cumulative distribution function. Both of these steps require the integration of f over [0, 1], which we were trying to avoid in the first place. In practice, any
probability density function f with a “similar shape” will cause a major variance reduction.

3.2.2 Stratified sampling
Sampling uniformly distributed random variables often causes the accumulation of many
sample points in one location (“clumping”). This effect slows down convergence of Monte
Carlo methods and can be avoided by distributing samples more uniformly. The fundamental idea of stratified sampling is to subdivide the integration domain into several disjoint intervals (“strata”) Di such thatZZ
Z1
f (x) d x + · · · +f (x) d x =0D0

f (x) d x

(3.9)

Dm

and then equally distributing the samples amongst these intervals, which causes negligible
computational overhead. According to Lafortune [4], the variance of stratified sampling isalways less than or equal to the variance of the estimator 〈J〉. This technique is effective,
easy to implement, and can be combined with importance sampling.

3.3 Russian rouletteA naïve algorithm numerically approximating each integral in the path formulation of therendering equation (2.23) would never terminate. For a well-defined algorithmic implementation, a stop condition, where evaluation terminates, is needed. Care must be taken
as not to introduce any bias into the estimator.

25

3.4 E  F 

3 M C I

Russian-roulette defines an estimator 〈J〉 r r based on 〈J〉 , which evaluates to zero with aprobability a (as in absorption).
(0,with probability a〈J〉 r r :=−1
(1 − a) · 〈J〉, with probability 1 − aThis estimator is still unbiased, although the variance is obviously increased.
E[〈J〉 r r ] = a · 0 + (1 − a) ·11−a· E[〈J〉] = E[〈J〉](3.10)

3.4 Estimators for Fredholm equations
Integral equations of the formf (x) = g(x) +Z1

K(x, y) f ( y) d y(3.11)

0where K is the kernel of the associated integral operator, g is a known function and f is
unknown are called inhomogeneous Fredholm equations of the second kind. Lafortune [4]uses an unbiased estimator for f using importance sampling:
K(x, X 1 )〈 f (X 1 )〉p1 (X 1 )K(X 1 , X 2 )K(x, X 1 )g(X 1 ) +〈 f (X 3 )〉= g(x) +p1 (X 1 )p2 (X 2 )K(x, X 1 )K(x, X 1 ) K(X 1 , X 2 )= g(x) +g(X 1 ) +g(X 2 ) + · · ·p1 (X 1 )p1 (X 1 ) p2 (X 2 )
!
∞
i
X
Y
K(X j−1 , X j )=g(X i )
p
(X
)
j
j
i=0
j=1

〈 f (x)〉 = g(x) +(3.12)

where X 1 , . . . , X n are random variables distributed according to the probability density
functions p1 , . . . , pn and X 0 = x . This estimator can be combined with russian roulette
for an algorithmic implementation (Listing 3.1).

3.5 Sampling the hemisphere
There are several different approaches to reduce variance when sampling the hemisphereabove a surface point.
(i) Cosine lobe sampling: The rendering equation (2.18) includes a cosine term since
light at grazing angles contributes less radiant energy to the irradiance at the differential surface position. Cosine lobe sampling uses this fact to prefer angles that arecloser to the surface normal.

26

3 M C I
1
2
3
4
5
6
7
8

3.5 S  

MC-I (x)s ← Uniform sample from [0, 1]v←0if (s < a)return // Absorptions ← Sample [0,1] using the PDF pv ← g(s) + (MC-I(s) · K(x, s))/p(s)return v / (1 − a)Listing 3.1: An unbiased primary estimator for Fredholm equations

To create random numbers with the distribution pcl , we can, for example, calculate
the inverse of the cumulative distribution function Fcl . The probability density function pcl is given bycos(Nx , Θ) cos(θ )pcl (Θ) ==ππwhere Θ = (ϕ, θ ) and θ is the elevation on the hemisphere.
Fcl (ϕ, θ ) =Z ϕZ θ
0

=

Zϕ
0

=

ϕ

2π

cos(θ 0 )
π

0

1
π



sin(θ 0 ) dθ 0 dϕ 0

cos2 (θ 0 )

θ

2

dϕ 0
0

· (1 − cos2 (θ )) =: Fϕ (ϕ) · Fθ (θ )The cumulative distribution function Fcl is separable and the inverse of its parts caneasily be calculated by setting ϕ to 2π or θ to π2 .
Fϕ−1 (ξ1 ) = 2πξ1Fθ−1 (ξ2 ) = cos−1 (p

1 − ξ2 )(3.13)(3.14)

Since ξi will be uniform random numbers on the interval [0, 1], the term 1 − ξ2 in(3.14) can be replaced by ξ2 . Applying these transformations to ξ1 and ξ2 will resultin directions that are distributed according to pcl .
(ii) BRDF sampling: When a surface reflection occurs, BRDF sampling prefers directions
which transport much radiant energy. Unfortunately, creating samples with such a
distribution is not always possible using analytic methods. While BRDF sampling is
superior to cosine lobe sampling, it generally has a higher computational cost. This
project uses the modified Blinn-Phong BRDF because proper samples can be generated using a simple closed form.

27

3.5 S  

28

3 M C I

4 Bidirectional path tracing
Bidirectional path tracing was independently developed by Veach and Guibas [10] andLafortune and Willems [4]. It unites and generalizes the path tracing and light tracing
algorithms while substantially reducing their variance.

4.1 Path tracing and light tracing
The path tracing technique proposed by Kajiya [3] iteratively solves the rendering equation
by recursively tracing rays through the scene. Whenever a ray hits a non-emissive surface,
it is randomly either absorbed, reflected or refracted. When a ray hits a light source, thecontribution of the hitting path is accumulated on the camera’s film. The construction ofthese so-called random walks constitutes an estimator of the rendering equation’s solution.
The reciprocity laws make this approach physically correct even though photons actually
travel in the opposite direction. A major improvement of the technique samples the lightsources at every intersection to rectify cases where the probability of this event is almost
zero. A significant drawback of the path tracing technique is its low convergence rate. Even
with the mentioned modifications, more than 1024 samples per pixel are needed for themost basic of scenes.
The light tracing technique is the dual version of path tracing – photons are traced fromthe light source to the viewpoint. Light tracing can handle some scenes well where pathtracing takes a very long time to converge. However, the number of of photons required for
an image without noise is staggeringly high. Instead of tracing radiance, the light tracing
technique propagates a quantity called importance. Importance describes how much influence an infinitesimal light source positioned at some point has on the scene. Importancecan be handled like radiance and all presented equations can also be applied to it. Tracing
radiance and importance are essentially symmetric operations and it can be proven thatthe linear operators expressing them are mutually adjoint.

4.2 A bidirectional estimator
Bidirectional path tracing works by performing two random walks at the same time –one from the light source and another one from the viewpoint. Both radiance and importance are traced and later combined. Importance sampling, stratified sampling andrussian roulette can all be used to create samples. The path’s vertices are denoted byx i (i = 0, . . . , S) on the eye path and yi (i = 0, . . . , T ) on the light path.

29

4.2 A  

x2

4 B  
Eye path
x1
Image plane

y0

x0

y2

y1

Light path
Mutal visibility test
Path rayFigure 4.1: Overview of the bidirectional path tracing algorithm
In a second step, the mutual visibility of hit points on the two paths is tested and stored(Figure 4.1). Finally, the contribution of every possible composite path is calculated. Let
C t,s be the path created by concatenating the vertices y0 , . . . , y t and x s , . . . , x 0 . It consistsof t steps on the light path, s steps on the eye path and one additional deterministic step
connecting the two sub-paths. We would like to create a partial estimator 〈C t,s 〉 to calculate
the flux over this composite path. Lafortune [4] differentiates several cases:(i) s = 0, t = 0: The light source is directly visible from the eye and the affected pixel
coordinates have to be determined a posteriori. G is the geometric term known from
(2.12).
→〈C0,0 〉 = L e ( y0 , −y−
0 x 0 ) G(x 0 , y0 )(ii) s > 0, t = 0: Paths with this configuration correspond to classic path tracing – thelight has to be explicitly sampled for this path configuration. Let the normal vectors
at the surface interaction points be named n x i and n yi :→−−→ −−−→−−→〈C0,s 〉 =L e ( y0 , −y−
0 x s ) f r (x s , x s y0 , x s x S−1 ) cos(n x s , x s y0 )S−1Y
→ −−−→−−−→f r (x i , −x−−
i x i+1 , x i x i−1 ) cos(n x i , x i x i+1 )
i=1This complicated-looking estimator is simply the product of the BRDF and exitant cosine values between surface interactions on the eye path and between the light source
and the last eye path interaction.
(iii) s = 0, t > 0: Paths with this configuration correspond to classic light tracing. The30

4 B  

4.3 P  

affected pixel coordinates have to be determined a posteriori.
T −1
Y

→
〈C t,0 〉 =L e ( y0 , −y−
0 y1 )→ −−−→−−−→f r ( yi , −y−−
i yi+1 , yi y−1 ) cos(n yi , yi yi+1 )!

i=1→ −−−−→−−→f r ( yt , −y−t x 0 , y t y T −1 ) cos(n ys , ys x 0 )(iv) s > 0, t > 0: The composite path consists of t reflections on the light path and s
reflections on the eye path.
!
T −1
Y〈C 〉 =L ( y , −y−→y )f (y ,−y−−
y →, −y−−y→) cos(n , −
y−−
y →)
t,s

e

0

0 1

r

i

i

i+1

i

−1

yi

i

i+1

i=1

→ −−−−→−−→ −−−→f r ( yt , −y−t x s , y t y T −1 ) G( y t , x t ) f r (x s , x s y t , x s x S−1 )!
S−1Yf (x , −x−−x →, −
x−−x →) cos(n , −
x−−
x →)
r

i

i

i+1

i

i−1

xi

i

i+1

i=1

When russian roulette or importance sampling are used, additional factors have to be taken
into account. The flux arriving at the eye location x 0 can be expressed as a weighted sumof the partial estimators multiplied by the visibility function V :
Φ=

T X
S
X

α t,s · V ( y t , x s ) · 〈C t,s 〉t=0 s=0

These weights cannot be chosen arbitrarily – some constraints need to be satisfied or the
resulting estimator might be biased:NX

α t,N −t = 1(N = 0, 1, . . .)t=0

The rationale is that the groups of composite paths with lengths 1, 2, . . . each receive
weights that add up to one. For this project, simple coefficients of the form
α t,s =

1
t +s+1

were chosen. It can be seen that they satisfy the above condition.

4.3 Pure software implementation
The following code consists of simplified C++ excerpts of a pure software implementationof the bidirectional path tracing algorithm.
The main rendering loop (Listing 4.2) uses a loadable sampling technique to generate animage plane sample and transforms it into an eye ray. Afterwards, a light ray is sampled31

4.3 P  
1
2
3
4
5
6
7
8

4 B  

struct Vertex {Point p;// Surface point in world coordinates
Normal n;// Surface normal (world coordinates)
Vector dgdu, dgdv;// Tangent and binormal vectors
Vector wi, wo;// Incident and exitant directions (local)
const BSDF *bsdf;// Pointer to the surface’s BSDF
Spectrum cumulative; // Cumulative attenuation factor};

Listing 4.1: Surface interaction data structure
1
2
3
4
5
6
7
8
9
10
11
12
13

Sample eyeSample;
Ray eyeRay;
while (m_sampler->getNextSample(eyeSample)) {m_camera->generateRay(eyeSample, eyeRay);
Spectrum Le =sampleLight(lightRay, lightPdf);

if (lightPdf != 0.0f) {Le /= lightPdf;
renderPaths(eyeSample, eyeRay, lightRay, Le);}}

Listing 4.2: Main rendering loop

and the method BiDir::renderPaths is invoked on these values. These two rays form
the start of the eye and light path random walks. Before the renderPaths function is
discussed, we will first introduce two utility functions:
The method BiDir::randomWalk in listing 4.3 performs a random walk either on aneye or a light path. This is done by following one initial ray and then sampling the BSDF
at each intersection point. A cumulative attenuation factor is stored in the cumulative
variable. The evaluation is only stopped once an absorption event occurs.
The method BiDir::evalPath in listing 4.4 takes two arrays containing the surface interactions of previously created eye and light paths. The paths are then evaluated according
to case (iv) of the algorithm.
The BiDir::renderPaths function uses these methods to evaluate cases (ii) to (iv).
Case (i) is omitted here for simplicity. It can be seen how the source code implementation
corresponds exactly to the previously introduced estimator types.

32

4 B  
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

4.3 P  

void BiDir::randomWalk(Ray ray,// Initial ray
std::vector<Vertex> &vertices, // Target arraybool eye) const {// Eye or light path?
Intersection its;
Spectrum cumulative(1.0f); // Cumulative attenuation factor
while (true) { // Repeat until an absorption event occurs
if (!rayIntersect(ray, its)) // Intersect the ray withbreak;// the scene
const BSDF *bsdf = its.mesh->getBSDF();float bsdfPdf;
Vertex v;// Store the local surface geometry:
v.p = its.p;// - intersection point
v.n = its.n;// - normal vector
v.dgdu = its.dgdu; // - tangent vector
v.dgdv = its.dgdv; // - binormal vector
v.bsdf = bsdf;// - local BSDF// Ensure wi/wo are stored with regard to the requested typeif (eye) {/* Convert incoming dir. into local coordinate system */v.wo = BSDF::toLocal(-ray.d, its);
cumulative *= bsdf->sample(v.wo, v.wi, bsdfPdf);
ray.d = BSDF::toWorld(v.wi, its); // New ray direction/* Multiply by the absolute value of the outgoing cosineand divide by the PDF value (importance sampling) */cumulative *= std::abs(BSDF::cosTheta(v.wi)) / bsdfPdf;} else {/* Convert incoming dir. into local coordinate system */v.wi = BSDF::toLocal(-ray.d, its);
cumulative *= bsdf->sample(v.wi, v.wo, bsdfPdf);
ray.d = BSDF::toWorld(v.wo, its); // New ray direction/* Multiply by the absolute value of the outgoing cosineand divide by the PDF value (importance sampling) */cumulative *= std::abs(BSDF::cosTheta(v.wo)) / bsdfPdf;}
v.cumulative = cumulative;if (cumulative.isBlack()) // Absorption
break;/* Store the surface interaction */vertices.push_back(v);ray.o = its.p; // New ray origin: current intersection point}}

Listing 4.3: Random walk function33

4.3 P  

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

4 B  

Spectrum BiDir::evalPath(
const std::vector<Vertex> &eyePath,// Eye and light path
const std::vector<Vertex> &lightPath, // surface interactions
int nEye, int nLight) const {// Number of steps// on each path
const Vertex &ev = eyePath[nEye-1];
const Vertex &lv = lightPath[nLight-1];
Spectrum L(1.0f); /* Start with an importance of 1.0 *//* Account for reflections on the eye and light paths */if (nEye > 1)
L *= eyePath[nEye-2].cumulative;
if (nLight > 1)
L *= lightPath[nLight-2].cumulative;/* Calculate the geometric term */Vector etl = lv.p - ev.p;
float lengthSquared = etl.lengthSquared();etl /= std::sqrt(lengthSquared);float geoTerm = absDot(etl, ev.n) * absDot(-etl, lv.n)/ lengthSquared;/* Extremely close points cause numerical problems */if (lengthSquared < 0.05)
return Spectrum(0.0f);/* Evaluate the BSDFs at the last light and eye pathsurface interaction points */L *= ev.bsdf->f(ev.wi, BSDF::toLocal(etl, ev.dgdu, ev.dgdv, ev.n));
L *= lv.bsdf->f(
BSDF::toLocal(-etl, lv.dgdu, lv.dgdv, lv.n), lv.wo);
L *= geoTerm; // Attenuate by the geometric term
return L;}

Listing 4.4: Path evaluation function34

4 B  

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

4.3 P  

void BiDir::renderPaths(const Sample &eyeSample,
const Ray &eyeRay, const Ray &lightRay,
const Spectrum &Le) {Vector wi;
Point onLight;/* Perform two random walks */randomWalk(eyeRay, m_eyePath, true);randomWalk(lightRay, m_lightPath, false);/* Implementation of cases (ii) and (iv) */for (unsigned int i=1; i<m_eyePath.size()+1; i++) {const Vertex &ev = m_eyePath[i-1];// Sample the light source from the current surface location
Spectrum localLe = m_lights[0]->sample(ev.p, onLight, wi);// Convert the incident light dir. to the local coord. sys.
Vector localWi = BSDF::toLocal(wi, ev.dgdu, ev.dgdv, ev.n);
float cosTheta = std::abs(BSDF::cosTheta(localWi));
if (V(ev.p, onLight) && cosTheta > 0) { // Visibility test/* Account for multiple reflections */if (i > 1)
localLe *= m_eyePath[i-2].cumulative;/* Store the contribution on the film */m_film->addSample(eyeSample,
localLe * ev.bsdf->f(ev.wo, localWi) * cosTheta *pathWeight(i, 0));}/* Case (iv) */for (unsigned int j=1; j<m_lightPath.size()+1; j++) {const Vertex &lv = m_lightPath[j-1];
if (V(ev.p, lv.p)) { // Mutual visibility test/* Store the contribution on the film */m_film->addSample(eyeSample, Le *evalPath(m_eyePath, i, lightRay, m_lightPath, j)* pathWeight(i, j));}}}// Continued on the next page ..

Listing 4.5: Path rendering function (part 1)

35

4.3 P  

41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79

4 B  

// BiDir::renderPaths continued ..
/* Implementation of case (iii) */for (unsigned int j=1; j<m_lightPath.size()+1; j++) {const Vertex &lv = m_lightPath[j-1];
if (V(lv.p, m_camera->getPosition())) { // Visibility test
Spectrum localLe = Le;/* Calculate the exitant direction */Vector wo = m_camera->getPosition() - lv.p;
float lengthSquared = wo.lengthSquared();wo /= std::sqrt(lengthSquared);/* Extremely close points cause numerical problems */if (lengthSquared < .05)
continue;/* Convert exitant vector to the local coord. sys. */Vector localWo = BSDF::toLocal(
wo, lv.dgdu, lv.dgdv, lv.n);/* Account for multiple reflections */if (j > 1)
localLe *= m_lightPath[j-2].cumulative;
localLe *= lv.bsdf->f(lv.wi, localWo)* std::abs(BSDF::cosTheta(localWo))/ lengthSquared;/* Determine the position on the image planeand store the contribution if the positionis valid */Sample s;if (m_camera->positionToSample(lv.p, s))m_film->addSample(s, localLe * pathWeight(0, j));}}}Listing 4.6: Path rendering function (part 2)36

4 B  

4.4 R

4.4 Results
Example renderings using the pure software implementation of the bidirectional path tracing algorithm showed very good results at 128 to 256 samples per pixel. The noise wasconsiderably lower compared to path tracing renderings created with an equal amount ofsamples. The rendering time, however, was quite long – about 30 minutes were required torender a Cornell box scene at 512x512 without visible noise.
Caustics caused by specular transmission and refraction are usually notoriously hard to
render – as expected, the algorithm was able to reproduce a caustic without problems. The
result of the Cornell box scene and the image representations of the partial estimators canbe seen in figures 4.2 and 4.3. Note the color bleeding on the walls and the light reflections
on the ceiling and floor caused by the mirror object.

Figure 4.2: Cornell box rendering (256 samples, ray depth: 3). Monkey courtesy of theBlender project.

37

4.4 R

4 B  

(a) i=0, j=0(b) i=0, j=1(c) i=0, j=2(d) i=0, j=3(e) i=1, j=0(f) i=1, j=1(g) i=1, j=2(h) i=1, j=3(i) i=2, j=0(j) i=2, j=1(k) i=2, j=2(l) i=2, j=3(m) i=3, j=0(n) i=3, j=1(o) i=3, j=2(p) i=3, j=3

Figure 4.3: Cornell box rendering grouped by path configuration. Each image contains thepaths generated by one of the partial estimators – “i ” and “ j ” denote the eye
and light path steps, respectively.

38

4 B  

4.5 T    

4.5 The bipartite mutual visibility query
Since (n + 1)2 mutual visibility checks have to be performed for a pair of paths of lengthn, these costs will eventually dominate the total computation time with increasing valuesof n. Therefore, this operation was chosen for acceleration using dedicated hardware. The
implemented operation works as follows:
1
2
3
4
5
6

B-M-V-Q(X, Y, S, T)
for i ← 0 to S dofor j ← 0 to T do
check whether x i and y j are mutually visible
store the result in bit vector
return bit vector

This operation has the highly beneficial property of transporting a computationally complex problem with minimal communication overhead – especially, when using long lightand eye random walks. Only 6n floating-point values need to be transferred to the intersection co-processor which, in turn, sends a small result bit vector. On the other hand, sucha request causes (n + 1)2 intersection tests on the co-processor. This makes it possible tobefefit from such acceleration using only low-bandwith links.

394.5 T    

40

4 B  

5 Target architecture & design flow
A Field-Programmable Gate Array (FPGA) is a semiconductor device consisting of reconfigurable logic resources embedded within a reconfigurable interconnect. The high degreeof programmability enables the implementation of almost any logic circuit, including recent developments such as fully asynchronous circuits. The size and speed improvements
of modern FPGAs make them viable for very complex designs, especially in the domain ofhigh performance computing, where the parallelism of their logic resources pays off compared to current general-purpose CPUs. This chapter discusses the FPGA architecture anddesign flow of Xilinx [13] parts, however, parts by other manufacturers are usually very
similar in structure and some of the information provided here will also apply to them.

5.1 Logic resources
The Configurable Logic Block (CLB), is the main resource for the implementation of combinatorial and sequential circuits. It can be further subdivided into four slices and a switch
matrix, as shown in figure 5.2. The switch matrix provides very fast local connectivity
amongst the four slices in a CLB and somewhat slower connectivity to adjacent CLBs. Eachslice (Figure 5.3) contains two registers and two 4-input look-up tables (LUTs), some ofwhich can also be configured to operate as 16-bit RAM/ROM or as a shift register. Additionally, each slice contains dedicated logic for the creation of extremely fast ripple-carry 1
adders.

IOB
CLB

Memory
Cell

Mult.

CLB

CLB

CLB

IOB

IOB
IOB

IOB

IOB

Figure 5.1: A generic FPGA array architecture – including I/O blocks, on-chip memory
slices and multipliers embedded in a reconfigurable interconnect
1

These adders consist of a cascade of 1-bit full adders, where the carry bits “ripple” from the last significant
digit to the most significant digit.

415.1 L 

5 T  &  Logic/RAM/ Logic-onlyslices
Shift-registerslices
COUT

COUT

Slice
X1Y1
Slice
X1Y0

SwitchMatrix

Slice
X0Y1

Neighboring
CLB

Slice
X0Y0
CIN

CIN

Figure 5.2: Layout of a Configurable Logic Block

COUT (to upper slice)
YB

MUX
G4

A4

G3

A3

4-Input

G2

A2

LUT

G1

A1

D
D

Q

YQ

CE
D-FF/
Latch
Y

From
Interconnect

ToInterconnect

1
0

XB

MUX
F4

A4

F3

A3

F2

A2

F1

A1

D

4-Input
D

Q

XQ

CE

LUT

D-FF/
Latch
X

1
0
CE
CLK

CIN (from lower slice)

Figure 5.3: Simplified diagram of a Virtex-4 slice (wide function circuitry not drawn). The
highlighted elements can be configured by a bitstream.

42

5 T  &  5.1 L 

5.1.1 Look-up tables
A look-up table can implement an arbitrary boolean function with four inputs. The timerequired for the function generator to react to a change in inputs is very short ( TI LO =0.2ns) and independent of the implemented function. By combining look-up tables withdedicated wide function multiplexers available in each slice (not drawn in figure 5.3), fastboolean functions with up to 8 inputs are possible.

5.1.2 Interconnect
Routing is one of the most complex tasks of FPGA design and this reflects in the abundanceof interconnect types available to the developer: There are CLB-internal lines, direct lines,
double lines, hex lines, long lines (Figure 5.4), each possessing individual electrical capacitance and skew characteristics. Fortunately, there is software which will perform automatic
placement and routing of logic on the FPGA. However, it is necessary to be aware of skew,
which might be introduced by inappropriate usage of the interconnect resources. Skew is
the difference in path delays when a signal simultaneously travels to several different destinations. This is a particular problem when occurring on nets distributing clock pulses tothousands of flip-flops because it removes the synchronization that is so essential to any
synchronous circuit. Thus, it has to be reduced at all costs.

CLB

CLB

CLB
CLB

CLB

CLB

CLB

CLB

(b) Double lines
CLB

CLB

CLB

CLB

CLB

CLB

CLB

CLB

6

(a) Local lines

CLB

CLB

(c) Hex lines

CLB

CLB

CLB

6

CLB

6

(d) Long lines

Figure 5.4: Hierarchical interconnect resources

5.1.3 Clocking
Since the most of today’s logic circuits are synchronous designs, they have to be supported
by an external clock source. Virtex-4 devices include an extensive number of systems tofacilitate this task. A total of 32 dedicated low-skew clock distribution networks exist on
the chip, a subset of which leads to the clock input of every logic resource on the FPGA.
These networks can only be driven by special clock buffers, which ensures that a clock
impulse can be delivered within hundreds of picoseconds and almost no skew.

43

5.2 D 5 T  &  Another important clocking component is the Digital Clock Manager (DCM), which isjust a Xilinx-specific term for delay-locked loop (DLL). High-speed systems running at frequencies in excess of 100 Mhz frequently encounter synchronization anomalies as a resultof this speed. Suppose two components C1 and C2 located on the same circuit board exchange data through a communications link. Additionally, C1 forwards a clock signal toC2 , which C2 uses to internally clock its synchronous logic resources. However, the copper
traces connecting C1 and C2 have a capacitive value, which causes a propagation delay inthe order of several nanoseconds when forwarding the clock signal. As a result, C1 andC2 receive the rising edge of the clock signal at different times and thus, are not properly
synchronized from a system perspective.
This is where the DCM comes in. In order to remove this anomaly, the circuit board
needs to provide a feedback loop trace leading from component C1 to component C2 andback again. The DCM will send a clock signal through the feedback loop and measure
the delay caused by I/O drivers and capacitive values. Through observation of the clock
signal, it is able to “anticipate” its future behavior and, using this information, constructs apredicted version of the clock signal. This signal is then forwarded to C2 and again delayedby the trace capacitance. However, this time, the offset into the future and path delays
cancel each other out, causing C1 and C2 to be in perfect synchronization.

5.1.4 DSP48 tiles
Newer Xilinx FPGAs include “DSP48” tiles, which basically consist of a 18×18 two’s complement multiplier followed by an 48-bit sign-extended adder/subtracter/accumulator [14].
On the chip surface, these are arranged in vertical columns and use a dedicated 48-bit internal bus so that a tile can “cascade” data upwards to the next tile. The pipelining of
input operands, intermediate products and the adder/subtracter-outputs is completely programmable. The arithmetic mode of operation of a tile can be changed dynamically using
an opcode input. An optional 17-bit right shift makes it possible to perform any of these
operations with an arbitrary precision.

5.1.5 Block RAM tiles
On the target FPGA, there is also an ample amount of 18Kb block RAM tiles, which canbe configured to use “aspect ratios” from 512×36, to 16K×1. Each tile has two totally
independent ports supporting full read and write access. Some very interesting applications
are possible thanks to the dual-port capability — including fast direct–mapped caches andshared memory communication. Furthermore, the memory block can be “pre-configured”
using the configuration bitstream and can thus be used as ROM or as a microsequencer.

5.2 Design flow
This section documents the creation of a tiny piece of hardware to illustrate the typical design flow. The circuit is a moore state machine, which reads a stream of bits and determines

44

5 T  &  BCOUT

5.2 D Cascade out to next slice

Conc.

18

PCOUT

X

36

48

18

A

B

48

36

18

×

72
36

Outputs two 36-bit partial
products
18

Y

48

±

P

48

0

C

48

17-bit right shift

Shared withneighboring
DSP48 sliceZ48

17-bit right shift

18

48

Cascade from previous slice

BCIN

PCIN

Figure 5.5: DSP48 slice architecture as described in [14].
Device

CLB matrix

SlicesTotal LUTs

Slice RAM (Kb)

Block RAM (Kb)

DSP48sXC4VLX15
XC4VLX25
XC4VLX60 (∗)
XC4VLX100
XC4VLX200

64×24
96×28
128×52
192×64
192×116

6,144
10,752
26,624
49,12589,08812,28821,50453,24898,394
178,176

96
168
416
768
1,392

864
1,296
2,880
4,320
6,048

32
48
64
96
96

Table 5.1: Logic resources available in current Virtex-4 devices [12].
whether the amount read so far contained an even or odd number of ones (parity).

1
0

Even

Odd

0

1

Figure 5.6: Diagram of the state machine

5.2.1 Logic design
Logic design is done using the Very High Speed Integrated Hardware Description Language(VHDL), a language commonly used for the simulation and synthesis of digital circuits.
VHDL supports modularity, re-use and encourages the developer to partition his design
into many interconnected entities. An entity possibly possesses multiple architectures. This45

5.2 D 5 T  &  can be compared to Java’s interfaces and implementation classes. Exemplary code for thestate machine example is provided in listings 5.1 and 5.2.
1
2
3
4
5
6
7
8
9
10
11

library ieee;use ieee.std_logic_1164.all;-- Import std_logic data type-- (corresponds to a single bit)

entity my_fsm is-- Declare the entity ’my_fsm’
port (-- Interface to the outside worldclk: in std_logic; -- Clock input
rst: in std_logic; -- Synchronous reset input
input : in std_logic; -- FSM input
output : out std_logic -- FSM output);
end my_fsm;

Listing 5.1: Entity declaration of the finite state machine

5.2.2 Synthesis and mapping
The synthesis step collects a number of VHDL-files and transforms this behavioral circuitdescription into an equivalent logic description — while the result has very little resemblance with the original source code, it is completely faithful to the specified functionality.
This device-independent synthesized logic description is then translated and mapped
onto the available device resources (e.g. LUTs, flip flops). Figure 5.7 is the result of invoking
the Xilinx tools on the provided state machine code.
IBUF
input

I

LUT3
FD

O

OBUF
I

O

output

IBUF
rst

I

O

IBUFG
clk

I

O

BUFG
I

O

Figure 5.7: Synthesized parity calculation circuit

5.2.3 Placement and routing
The placement and routing steps continue with the results from synthesis and mapping.
Furthermore, they need a so–called constraints file, which forces them to adhere to additional functional specifications such as operating frequency in MHz or more fine-grained
timing requirements. Unfortunately, placement and routing of a complex circuit can take46

5 T  &  
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

5.2 D architecture rtl of my_fsm is-- Declare an architecture ’rtl’
type state_t is-- ’RTL’=Register Transfer Level(STATE_EVEN, STATE_ODD); -- New state data typesignal state_r : state_t;-- Declare the state register
begin
output <= ’0’ when state_r = STATE_EVEN else ’1’;
fsm_update: process (clk)-- Sensitive to clock signal
begin
if rising_edge(clk) then -- When a rising clock pulse occurs
if rst = ’1’ then-- Synchronous reset detection
state_r <= STATE_EVEN;
else-- otherwise invoke the statecase state_r is -- transition function
when STATE_EVEN =>if input = ’1’ thenstate_r <= STATE_ODD;
end if;when STATE_ODD =>if input = ’1’ thenstate_r <= STATE_EVEN;
end if;end case;end if;end if;end process fsm_update;end rtl;

Listing 5.2: Architecture declaration of the finite state machinefrom hours to days — dependent on how “strict” these constraints are. The reason is thatthe underlying problems are NP-hard and have to be solved heuristically. If the whole
design is supposed to run at 100 Mhz, signals need to travel from flip-flop through combinatorial logic to flip-flop within 10 nanoseconds! Otherwise, the previous result is stored
and erroneous circuit operation will ensue. As a consequence, these steps need to ensurethat all components are placed appropriately (e.g. close enough on the chip surface) tosatisfy the constraints.
While there are many third-party synthesis and mapping solutions, the final placement
and routing has to be performed with proprietary tools made by the device manufacturer.
As last part of the design flow, these tools produce a configuration bitstream, which canthen be uploaded onto the target FPGA.

47

5.2 D 48

5 T  &  

6 Building blocks
When developing a core of this complexity, it becomes essential to make use of resources
already embedded in the FPGA fabric. This project was designed to fit into a Xilinx Virtex4 LX60 FPGA and extensively uses its internal dual-ported memory blocks and DSP48
pipelined multipliers. For reasons of cost and simplicity, a USB 2.0 interface was designed
to connect the intersection processor to a CPU. While a PCI or PCI Express interface would
likely involve much less CPU overhead and latency, the USB 2.0 variant still allows fast
data transfers and serves its purpose of demonstrating the performance of such a hardware
acceleration device.

6.1 Larger multipliers
The 18-bit operand width limitation of the DSP48 multiply-accumulate blocks would usually severely limit their use. However, by cascading four of these multiplier units, it ispossible to build a larger (35×35-bit) signed integer multiplier. The exact multiplication
scheme can be derived using divide-and-conquer: Let ai and bi be the i -th digit of the two’s
complement representation of a and b such that:
33
33
X
X
34
i
34
a = −2 · a34 +2 · ai and b = −2 · b34 +2i · bii=0

i=0

By rearranging the summation order, the sum can be split up into the sum of a shifted,
signed 18-bit value and an unsigned 17-bit value.
a = (−217 · a34 +16
X

2i ai+17 ) · 217 +

{z

}

=:ah

b = (−217 · b34 +16
X

2i bi+17 ) · 217 +i=0|2i ai .

i=0

i=0

|

16
X

{z

=:bh

| {z }=:al16
X

2i bi .

i=0}

| {z }=:bl

The multiplication can then be rewritten as:a · b = (ah · 217 + al )(bh · 217 + bl )= ah bh · 234 + 217 · ah bl + 217 · al bh + al bl

The summands are now composed of 17-bit unsigned multiplications and 18-bit signed
multiplications, which can be realized using the built-in DSP48-blocks.

496.2 N-  

6 B 

6.2 Non-restoring integer divider
Of the four basic arithmetic operations, division is the slowest and most complex in termsof hardware. Even on modern processors such as the Pentium-4, division takes in excess
of 40 clock cycles. Signed integer division is defined by the equation z = d · q + s and theconstraint sgn(z) = sgn(s), where z is the dividend, d is divisor, q is the quotient and s is
the remainder.
A fully pipelined, non-restoring divider is used to realize the division operation occurring inside the kd-tree traversal and triangle intersection units. This core is characterized
by a short critical path and moderate size. To clarify its operation, we will first examine thesimpler, restoring divider: Both are digit recurrence schemes, generating one quotient bit
during each cycle. Restoring division is very intuitive in that it is similar to how a humanwould perform division of two large numbers. If the trial difference of the shifted divisor
and the partial remainder from the n − 1th iteration is positive, the result is stored for the
next iteration. Otherwise, the previous value is restored and a 0 quotient bit is generated.
1
2
3
4
5
6
7
8
9
10

R-D(, )r ← z
for i ← n−1 downto 0 do
if r ≤ 2i · d thenr ← r − 2i · dq ← 2·q+1
else
q ← 2·q+0
end
end

Listing 6.1: Restoring division algorithm
However, restoring division has two disadvantages: Because the decision whether or not
to store the subtraction result depends on its sign, the timing-critical paths of the resulting logic are very long. Furthermore, additional sign conversion steps are required if thedivision core needs to handle two’s-complement numbers.
The non-restoring division algorithm avoids these drawbacks using a binary signed-digit(BSD) number representation. Depending on the partial reminder and dividend signs,
either an addition or a subtraction is performed during each cycle. The resulting quotient
bits are, unlike restoring division, not in {0, 1}, but {−1, 1} – as a result, this intermediate
quotient is always an odd number. This method can be compared to adding weights to bothsides of a scale until balance is reached. After the intermediate quotient bits in BSD form
have been calculated, two additional correction cycles are required for two’s-complement
conversion and remainder calculation. A thorough review of various hardware divisionalgorithms can be found in [7].

50

6 B 
1
2
3
4
5
6
7
8
9
10
11

6.3 F- N-D(, )r ← z
for i ← n−1 downto 0 do
if sgn(r) 6= sgn(d) thenr ← r − 2i · dqi ← 1elser ← r + 2i · dq i ← −1
end
end

Listing 6.2: Non-restoring division algorithm

6.3 Floating-point unit
During simulation runs, an implementation of the triangle intersector using 32-bit fixedpoint arithmetic showed signs of serious numerical problems. When a distant ray is transformed into the unit triangle space of a very small triangle, overflows occur. The inability of
fixed-point arithmetic to represent very large and very small numbers makes it unsuitable
for this design.
Instead, a simple 24-bit floating point format was chosen. It is neither compatible with
the IEEE 754 standard, nor does it support denormalized numbers or rounding modes.
However, the simplicity reflects in a simple and efficient hardware implementation, whichstill achieves an amazing precision during triangle intersection. Since this design’s floating
point calculation results are only stored temporarily, no rounding error accumulation will
occur.
23 22
S

E

16 15
E

E

E

E

E

0

E M M M M M M M M M M M M M M M M

7-bit exponent

16-bit mantissa, implicit 1

Figure 6.1: 24-bit floating-point format

6.3.1 Floating-point adder/subtracter
When comparing the integer versus floating-point complexity of the four basic arithmetic
operations, the adder/subtracter incurs the greatest hardware overhead. The reason forthis complexity is that the operands may need to be swapped and costly alignment preand post-shifts have to be performed. Given two floating point values a = (−1)sa · 2ea · ma
and b = (−1)s b · 2e b · m b , the core performs the following operations in sequence:(i) Calculate the difference of the operand’s exponents: δ ← |ea − e b |.

51

6.3 F- 6 B 

Unpack operands
fp_addsub_preShiftcalculation

Magnitude comparison& Operand swap

Alignment pre-shiftfp_addsub_shift

Mantissa magnitude comparison &Operation selection
fp_addsub_main
AddSub

Leading zero counter &Alignment post-shiftfp_addsub_post

Figure 6.2: Block diagram of the floating-point adder/subtracter. Dashed lines denotepipeline stages.
(ii) Possibly swap the operands so that a = 〈Operand with large exponent〉 and b =〈Operand with small Exponent〉.
(iii) Limit the alignment pre-shift: δ ← max{δ, 18}. If the exponent difference is greaterthan 18, b can be considered zero.
(iv) Perform an alignment pre-shift of the operand with the smaller exponent: m b ←m b · 2−δ(v) Operation selection depending on the operand signs and the requested operation(vi) Either c ← a ± b or c ← b ± a, so that the mantissa cm is non-negative.
(vii) Perform a normalization post-shift and correct the exponent
On Virtex4, steps 1-4 can be performed within 10 ns and form the first pipeline stage.
Steps 5-6 and step 7 are implemented in pipeline stages 2 and 3, respectively.

6.3.2 Floating-point multiplier/divider units
The floating point multiplier is much simpler. Given a and b, it performs the following
operations in sequence:
1. Calculate the result sign: sc ← sa ⊕ s b
2. Calculate an intermediate exponent: ec ← ea + e b
3. Calculate an intermediate mantissa: mc ← ma · m b52

6 B 

6.4 M 

Unpack operands
fp_mul_pre
Calculatesign (XOR)

Exponent
addition

Overflow /Underflow test

fp_mul_main

Alignment post-shiftfp_mul_post

Figure 6.3: Block diagram of the floating-point multiplier
4. In case of an exponent underflow: mc ← 0. In case of an overflow, flag the error.
5. The multiplication result of the two mantissas is in [1, 4). If mc > 2, an alignment
shift has to be performed and the exponent needs to be incremented.
The floating point divider works analogous to the multiplier: Instead of adding the exponents, they are now subtracted and the alignment post-shift performs a shift in the otherdirection. The pipelined non-restoring divider previously mentioned is used by the divider
core.
Core

4-Input LUTs

Flip-flops

SRL161 DSP48

FP add/sub
FP multiplier
FP divider
FP comparator

376
69
533
66

139
92
509
49

8
1
29
0

0
1
0
0

Available

53248

53248

26624

64

1

16-bit linear shift register

Table 6.1: Implementation costs of the floating-point cores

6.4 Memory hierarchy
6.4.1 DDR SDRAM Controller
Intersect uses a custom, pipelined DDR SDRAM controller capable of achieving 3,2Gb/s(381 MiB/s) at 100 Mhz. It takes care of refresh cycles, command timeouts, data strobes,
open rows and delay/skew issues. It supports read and write bursts of length 2 (16 bit
width), thus requiring one command every clock cycle to sustain full read/write speed.

53

6.5 USB I

6 B 

6.4.2 Caching
The project uses a 1024-line direct-mapped cache to speed up read accesses. Each cacheline stores sixteen 24-bit words amounting to a total of 48 KB of addressable memory.
Internally, the cache uses a series of dual ported RAM blocks to provide a single write
port and a read port for each of the three intersection units.
After a cache miss occurs, a read request is sent to the cache miss handler state machine.
The cache miss handler then instructs the DDR controller to bring in 16 consecutive words
using a burst transfer and updates the cache memory, which takes about 25 clock cycles
altogether. Since the cache uses no FIFO to store cache misses occurring while the statemachine is busy, special care must be taken to avoid cases where several threads fight toload data into the same cache line, which gets overwritten just before each thread is readyto receive its requested data.

6.5 USB Interface
The USB interface makes use of a Cypress EZ-USB FX2 controller chip on the circuit board.
This chip has been programmed with a firmware putting the device into FIFO mode. A simple state machine on the FPGA moves arriving/leaving data between the FX2’s “slave” FIFO
and internal FPGA FIFOs, which enables a theoretical peak performance of 480MBit/s. The
firmware of the USB controller is stored on an external EEPROM and automatically loadedon power-up. On the FPGA side, there is also some endianness detection/conversion logic.
As a consequence, the different possible hosts (e.g. Linux/i386, Linux/PowerPC) can use
the intersection processor without spending any CPU cycles on endianness conversion.

6.5.1 Device driver
As for now, there is only a Linux device driver. It uses DMA and asynchronous I/O for highspeed data transfers. A special USB scheduler manages a list of ongoing requests to keep
the USB bus saturated at all times. Up to sixteen simultaneous requests can be scheduled;this only works very well with the ehci-hcd host controller and current Linux versions(e.g. ≥ 2.6.17). The driver stores received packets in a 128-deep temporary packet buffer
and combines them to reduce the number of system calls (and context switches) required
to read data. The optimum transfer speed was achieved when scheduling up to eight bulk
output USB requests and eight inbound interrupt transfers at the same time. To evaluate
different modes of operation and chip configurations, a simple benchmark was designed. It
measures the raw read/write performance without any additional processing on the FPGA’s
side. The device driver achieved 27 MiB/s sustained write bursts and 31 MiB/s sustained
read bursts on a Pentium-4 PC running Linux 2.6.17 at 3 Ghz.

54

6 B 

6.5 USB I

Figure 6.4: Implementation hardware55

6.5 USB I

56

6 B 

7 Intersection unit
Table 7.1 shows an excerpt of a profiler report generated while running the pure software
implementation of the bidirectional path tracing algorithm. Since most time is spent doing ray-triangle intersection and kd-tree traversal, an acceleration device could drasticallyreduce the run-time of the algorithm.
% CPU time

Self seconds

Calls

Name

67.69
11.68
4.90
4.80
2.00
1.60
1.35
1.13
1.05
0.45
0.36
0.23

90.71
15.65
6.57
6.43
2.68
2.15
1.81
1.52
1.41
0.60
0.48
0.31

402295180
12519060
12540538
12540538
1048576
8347638
3312701
9723139
20971524192484
14799230
1

Triangle::rayIntersect
KDTree::rayIntersect
AABB::rayIntersect
KDTree::rayIntersect
BiDir::renderPaths
BiDir::V
BiDir::evalPathcoordinateSystem
BiDir::randomWalk
BSDF::toWorld
BSDF::toLocal
KDTree::buildTree

Table 7.1: A profiler report excerpt after rendering a scene using the pure-software implementation. Note that almost all CPU time is spent tracing rays!

7.1 Triangle intersection unit
Various intersection schemes were compared in the context of hardware implementation
feasibility. Established algorithms such as the fast intersection test by Möller and Trumbore [6] are unfortunately too complex for this purpose. Instead, a modified of Arenberg’s
algorithm [9] was chosen, which is also used on the SaarCOR chip by the University of
Saarland’s graphics group.

7.1.1 Intersection algorithm
For high performance operation, the ray-triangle intersection unit uses a redundant triangle storage format, which first has to be generated from the raw triangle data. This requires
two steps for every triangle : First, an affine transformation is created to transform the unit
triangle {(0, 0, 0), (1, 0, 0), (0, 1, 0)} onto the specified triangle. During the second step,57

7.1 T  7 I the affine transform’s inverse is calculated and converted to the processor’s 24-bit floatingpoint format. This pre-computation speeds up the actual intersection calculations andreduces hardware complexity by a great deal.
Let n be the normal vector of the triangle defined by A, B , C ∈ R3 .
n :=(B − A) × (C − A)k(B − A) × (C − A)k.

The affinity T : R3 → R3 maps the unit triangle onto 4ABC and maps (0, 0, 1) onto theorthogonal vector n.

Bx − Ax Cx − Ax nx
T := B y − A y C y − A y n y  + A
Bz − Az Cz − Az nz
If the triangle is well-formed (i.e. the points are not collinear), the normal vector n completes {B − A, C − A} to a base of R3 . Thus, T is an isomorphism and its inverse can becalculated – for example using the Gauss-Jordan algorithm.
All that is left to be done by the hardware in order to calculate an intersection point is totranslate the ray into unit triangle space with the help of the inverse affine transformationand then perform the intersection there, which is now a much easier problem.
1
2
3
4
5
6
7
8
9
10
11

I-T(T −1 , O, D)
O0 ← T −1 · O // Affine transformD0 ← T −1 · D // Linear transformt ← −Oz0 /Dz0
u ← Ox0 + D0x · tv ← O0y + D0y · t
if t ≥ 0 and u ≥ 0 and v ≥ 0 and u + v ≤ 1 thenreturn (t, u, v)
else
return nil

Listing 7.1: Ray-triangle intersection pseudocode

7.1.2 Hardware implementation
A diagram of the ray-triangle intersection data-path can be found in figure 7.1. This design
utilizes a massive amount of parallelism by unwrapping and pipelining the whole calculation. This results in the use of 20 fp24-multipliers, 18 fp24-adders, one fully pipelined
fp24-divider, and many register stages. A consequence of the pipelined design is that atriangle intersection has a latency of thirty clock cycles. However, since up to thirty triangle intersection requests can be simultaneously contained inside the pipeline, one request
is completed during each clock cycle, giving the unit a throughput of 100M Triangles per
second at 100Mhz.

58

7 I 7.1 T  

×

rdz

valid

×

+

a32
rdx

&

×

mint, maxt
check

+

a33
rdy

a31

/

t

+

a34

Legend

×

roz

+

a33
roy

Register stage++
×fp24 Multiplier

×

+

a23
rdy

×

×

a31
rdz

fp24 Adder

×

×
a32
rox

+

×

+

a22
rdx
a21

v

+

a24

×

roz

+

×

+

a22
rox

×

+

a13
rdy

×

×

a21
rdz

≤ 1 check

×

+

a23
roy

+

×

+

a12
rdx
a11

u

+

a14

×

roz

×

+

a13
roy

×

+

a12
rox
a11
2

3

3

13

2

3

3

1

Total: 30 pipeline stages

Figure 7.1: Ray-triangle intersection data path59

7.2 -  7 I 7.2 kd-tree intersection unit
A kd-tree is the geometric equivalent of a binary search tree. It is a data structure frequently
used to organize points or surfaces in a k-dimensional space. In comparison to the wellknown quadtree or octree data structures, kd-trees benefit from their simplicity and goodcache behavior, giving them better overall performance. A kd-tree is created by recursively
partitioning the scene geometry using a generated hyperplane. In contrast to the BSP-tree,
the hyperplanes are always axis-aligned, therefore the only information needed to uniquely
characterize a hyperplane is an axis index and a single coordinate. kd-trees that are used toorganize geometry usually store all primitives inside the tree’s leaf nodes.

7.2.1 Intersection algorithm
kd-trees are very efficient data structures for elementary operations such as nearest-neighborqueries or ray intersections. Listing 7.2 contains some pseudo-code to illustrate how thepartitioning information makes the intersection operation very efficient. Only triangles in
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

T-T(node, ray, mint, maxt)
if node.isLeaf() thenreturn the closest intersection in node.triangles or −1
if ray.d[node.axis] = 0 thennoSplitIntersection = true // Corner case: The ray is parallel to the intersection planeelse // Calculate the distance along the ray to the ray−plane intersection
splitIntersection = (node.splitPos − ray.o[node.axis]) / ray.d[node.axis]
if ray.o[node.axis] ≤ node.splitPos then // Determine which cell contains the ray’s origin
near = node.left
far = node.rightelse
near = node.rightfar = node.left
if noSplitIntersection or splitIntersection < 0 or splitIntersection >= maxt thenreturn T-T(near, ray, mint, maxt) // Whole interval on near cellelse if splitIntersection ≤ mint thenreturn T-T(far, ray, mint, maxt) // Whole interval on far cellelse // Traverse both cells in the proper ordert ← T-T(near, ray, mint, splitIntersection)
if t >= 0 then // An intersection has been found in the near cell, the far
return t// cell does not need to be checked.
return T-T(far, ray, splitIntersection, maxt)

Listing 7.2: kd-tree traversal pseudocode
leaf nodes, which are actually encountered by the ray, are checked for intersections. Also,
once the closest intersection in a node has been found, the evaluation can be stopped since
any intersection in other nodes will be farther away. The mint and maxt parameters specify
the considered search interval along the ray.

60

7 I 7.2 -  7.2.2 Generating kd-trees
The kd-tree construction algorithm is based on the Surface Area Heuristic, which is explained in detail in [11]. At each recursion level, a split plane is chosen using greedy selection. It can be proven that such a “perfect” greedy split plane passes through a triangle
vertex, restricting the set of possible candidates to a finite subset.
Given the set of candidate planes, the algorithm tries to calculate the overall cost ofintersecting arbitrary rays against each of the possible trees resulting from the partitioning
operation. It then selects the split plane with the lowest cost function value. The costfunction is defined asc(A, B) := ckd + ctri · (pA · nA + pB · nB )where ctri and ckd are the costs of intersecting a single triangle and performing one traversal
step, respectively. nA and nB are the number of triangles contained on each side of the partition (A, B) – overlapping triangles are counted twice. Finally, the conditional probabilities
pA and pB that, given an intersection of the parent node, A or B are also intersected, arecalculated using the following approximation:
pA =

SA
S

and

pB =

SB
S

where S is the surface area of the parent node and SA, SB are the surface areas of A and
B , respectively. This heuristic generally produces very good results. It can be proven thatthe algorithmic average case complexity of intersecting a ray against a tree containing nprimitives is in O(log n). Separate trees are built for the host CPU and intersection processor, since the relation between intersection and traversal costs is so fundamentally different
amongst the two.

7.2.3 Hardware implementation
Designing a kd-tree unit for hardware implementation is a challenge because of the complex control flow and stack requirements. This project used a bottom-up approach: first,
flat data-paths were developed for the triangle and kd-tree intersection units. Later, thesewere embedded into a control unit handling stack accesses and control flow.
Data path
The kd-tree data path (Figure 7.2) is a straight-forward implementation of the function
body in listing 7.2. Again, there is one expensive division step, which requires 13 pipeline
stages and makes up a big part of this unit’s hardware consumption. Since the length of
19 stages is much smaller than the triangle unit’s 30 stages, 11 empty register stages areappended so that the results of both units are available at the same time. This may be seen
as a waste of available FPGA resources – however, such register stages can be implemented
very cheaply using the built-in SRL16 linear shift registers.

61

7.2 -  7 I 

maxt

<

maxt

split

1

3

13

Decision
logic

−

MUX

axis
roz
roy
rox

midt

/

MUX

rdz
rdy
rdx

mint

<

mint

decision1

Total: 19 stages

1

Figure 7.2: kd-tree traversal data path
Control unit
The final core has a pipeline length of 32 stages and embeds the triangle and kd-tree intersection data paths. Each pipeline stage has an associated thread ID between 0 and 31, which
is used to ensure that each intersection thread only sees its own, private stack memory area.
The stack was implemented in a way such that the top item is always kept inside registers
and the remaining data is stored inside memory blocks.

Ray scheduler

Fetch stage
kd-treetraversal
pipeline(19 stages)

Triangleintersection
pipeline(30 stages)

Forward
Decision stage
Multi-stack(32×32-deep)

Figure 7.3: kd-tree control unit
When the decision stage pushes a value onto the stack, the registered top item is replacedwith new data and the old value is written out into the private stack memory area. When the627 I 7.2 -  decision stage pops an item from the stack, the top value is overwritten and will be available
during the next iteration. This all works because only one push or pop operation can occur
per thread per iteration. Here, the term iteration specifies the event of a single thread
passing through pipeline stages 1 to 30. The VHDL stack data type declaration is shown inlisting 7.3. The data type in listing 7.4 describes the complete state of an intersection thread
and exists as register storage inside each pipeline stage. This creates a massive amount ofregisters but, again, these registers can be cheaply implemented using linear shift registers.
1
2
3
4
5

type stack_t is record
addr: memptr_t; -- Node address in memorymint: fpval;-- Start of the checked intervalmaxt: fpval;-- End of the checked interval
end record;

Listing 7.3: Control unit – stack data type

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

type intersection_t is record
id: thread_id; -- Thread id
reqid: reqid_t;-- Request id
top: stack_t;-- Top of the stackvalid: boolean;-- Is this stage’s data valid or a bubble?
shadow: boolean;-- Does this stage contain a shadow ray?
fetch_ok : boolean;-- Was the data fetch sucessful?
kd: boolean;-- Does the addr. point to a kd-tree node?
sp: thraddr_t; -- Stack pointerro: fpval3;-- Ray originrd: fpval3;-- Ray directionu: fpval;-- U coord. of the closest intersectionv: fpval;-- V coord. of the closest intersection
hit: boolean;-- Was there an intersection so far?
tri_idx : triidx_t; -- Index of the current triangle
best_idx : triidx_t; -- Index of the best triangle
end record;

Listing 7.4: Control unit – pipeline stage data typeThe kd-tree control unit performs the following operations in each stage:(i) Schedule:1• If there is an active thread in the last pipeline stage, schedule it. This operation
basically copies all registers from the last pipeline stage.
1

This is technically not a pipeline stage because the schedule & fetch operations are performed within thesame clock cycle – however, it is listed separately here for clarity.

637.2 -  7 I • If the last stage contains a bubble and the input queue contains a ray to be intersected, schedule it. This operation initializes the thread’s stack pointer to zeroand pushes the ray and the memory address 0x0000002 on top of the stack. No
stack memory access is actually performed since the top was previously unused.
(ii) Fetch:• Fetch 16×24 bits = 48 bytes from the memory address specified by the top of
the stack using the direct-mapped cache. 48 bytes can contain one triangle or
up to four kd-tree nodes.
• If there is a cache miss, notify the cache and ignore this thread for the current
iteration by setting the decision stage to read-only when the thread arrives there.
This behavior repeats until the stack has finally brought in the requested cacheline.
(iii) Intersect or traverse (30 stages):• If the fetched memory block contains a triangle, check for an intersection withthe thread’s ray along the current search interval.
• If the top item references a kd-tree node, determine the order in which the childnodes should be traversed, and check for corner cases.
(iv) Decision:• If the previously performed memory fetch of this thread was unsuccessful, do
nothing.
• If the top item of the stack referenced a kd-tree node, push the far node onto thestack and overwrite the top of the stack with the near node. The search interval
is also reduced for both nodes. If a corner case occurred and one of the nodes
does not need to be checked at all, only overwrite the top of the stack.
• If the top item of the stack referenced the last triangle of a leaf kd-tree node,
there are several possibilities:– If the remaining stack is empty, the search is done and the core announces
whether an intersection has been found so far.
– If the remaining stack is not empty, the search is only terminated if an intersection has been found so far. If this is not the case, a node is popped
from the stack memory and the search continues with the next iteration.
• If the top item of the stack referenced a triangle of a leaf kd-tree node – but not
the last one, there are also several possibilities:– If no intersection was found, the top of the stack is changed to point to thenext triangle of the leaf node.
2

The root of the kd-tree is located at this position in memory64

7 I 7.3 B - 

– If an intersection was found and the thread does not contain a shadow ray,the core checks whether this intersection is closer than all previous intersection. The t , u and v coordinates are stored as well as the triangle index,
if this is found to be true.
– If an intersection was found and the thread contains a shadow ray, thesearch is aborted because an obstruction has been found.
Table 7.2 lists the resource requirements of the presented VHDL cores in terms of Virtex-4
resources. As expected, the triangle unit uses a relatively large amount of logic resources,
while the kd-tree control unit uses mostly storage elements. Because of the core’s size,
three intersection units are included on the final FPGA version, which allows 96 rays to beprocessed in parallel with a peak performance of 300M intersection operations per second
at 100Mhz.
Core

4-Input LUTs

Flip-flops

SRL161

DSP482

RAMB163

Triangle unit
kd-tree data pathkd-tree control

8676
1392
3625

3221
949
8465

606
206
2788

20
0
0

0
0
4

Combined
Available

13693 (26%)
53248

12635 (24%) 3600 (14%)
53248
26624

20 (31%) 4 (3%)
64
160

1

16-bit linear shift register
18-bit pipelined multiplier3
16kbit dual-port memory block
2

Table 7.2: Implementation costs of the intersection units

7.3 Bidirectional path-tracing frontend
The bidirectional path-tracing frontend is located between the USB interface and the intersection cores. It reads commands sent from the host CPU and executes them. This coreautomatically performs fp32 to fp24 conversion and vice versa to free the host CPU fromthis burden. Commands carry a request identification number, which is included whensending the response back to the CPU. The rationale is that responses may not be received
in the order, in which the request were previously sent. When using a software pipeline on
the host CPU, this information is needed to associate responses with the correct data.
This core also implements the bipartite mutual visibility query operation. When the CPU
sends such an operation to the FPGA, the FPGA first expects the three-dimensional coordinates of the eye path surface interactions followed by the light-path surface interactions.
The core issues shadow ray intersection requests for each possible pair (ei , l j ) and sends
visibility results back as soon as they exit the pipeline.

657.3 B - 

66

7 I 

8 Software
The design of the beam global illumination renderer was strongly influenced by the bookPhysically Based Rendering [8] by Matt Pharr and Greg Humphreys. beam uses the sameplugin-based approach, where almost anything is an external plugin: sampling techniques,light types, camera types and BRDFs are loaded on demand as the scene references them.
The plugin interfaces are specified in a very genereral way so that almost any technique canbe implemented without a change in the rendering core. beam uses a simple XML-based
scene format – an example can be seen in listing 8.1.
One additional feature is the use of the EXR high dynamic-range format to store image
data. When images turn out too bright or too dark after a long rendering time, they can be
re-developed using a different exposure time.

8.1 Software pipeline
One major difference, however, is the main rendering loop. Since the rendering task now
involves communication and external processing latencies, a software pipeline using 4
threads of execution was introduced. These are as follows:(i) Ray generator: This thread uses a custom sampling technique, a light source, and acamera to sample the start rays of a pair of random walks. In an additional step, theserandom walks are continued until a maximal path length is reached or an absorption
event occurs.
(ii) Write thread: This thread receives a packet of rays from the ray generator and sends
a mutual visibility query to the FPGA.
(iii) Read thread: This thread receives query responses from the FPGA, writes the results
to the affected packets, and then passes them to the processor thread.
(iv) Processor thread: The processor thread receives packets containing a number of random walks and visibility information. The only thing left to be done is to write thecontribution of each possible path onto the camera’s film. Afterwards, the packet isrecycled and given back to the ray generator.
The write and read threads have an elevated priority to get data into and out of the FPGA
at the earliest possible time. The pipeline works quite efficiently when many packets areinserted at the start. This approach keeps both the host CPU and the FPGA busy andmitigates the effects of external processing latencies.

678.1 S 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

8 S

<scene><!-- Specifies a triangle mesh using the OBJ file loader --><trimesh type="obj"><string name="filename" value="monkey.obj"/><!-- Sets the surface properties of the monkey mesh --><bsdf type="blinnphong"><float name="kd" value="0.7"/><float name="ks" value="0.0"/></bsdf></trimesh><!-- Instantiates a point light source --><light type="point"><transform name="lightToWorld"> <!-- Affine transform. --><translate x="0" y="1" z="0"/></transform><spectrum name="intensity" value="30, 30, 30"/></light><!-- Instantiates a perspective --><camera type="perspective"><transform name="cameraToWorld"><translate x="0" y="0" z="2.8"/></transform><float name="fov" value="50"/><!-- The camera uses the stratified sampling technique --><sampler type="stratified"><integer name="xSamples" value="2"/><integer name="ySamples" value="2"/><boolean name="jitter" value="true"/></sampler><!-- The camera uses a high-dynamic range film --><film type="exrfilm"><integer name="xRes" value="512"/><integer name="yRes" value="512"/></film></camera></scene>Listing 8.1: A scene using beam’s XML-based scene description language

689 Results
During the course of this project, three components were planned and developed: a global
illumination renderer, a linux kernel driver, and an intersection processor architecture. Various building blocks had to be designed in order to facilitate its development. Finally, thepresented intersection architecture was validated and synthesized for the target hardware.
As of this writing, the implementation of the co-processor works very reliable – although
it has a series of drawbacks, which made it impossible to exploit the full advantage of theacceleration. These disadvantages are mostly related to the evaluation board’s hardware
and are listed below:(i) USB 2.0 Link: The USB connection proved to be a significant bottleneck. Even
though the implemented operation was very efficient in terms of required data transfer, the FPGA usage was only about 5 percent under full USB load. Gigabytes ofvisibility information had to be transferred even for the most simple of scenes. Despite this severe bottleneck, a speedup of about 30% was measurable when offloading
the visibility computation.
(ii) Memory bandwidth: The memory connection of the used FPGA board had a bandwidth of 381 MiB/s and was only 16 bits wide, which is tiny compared to the 6,4
GiB/s of bandwidth and 128 bits of connectivity available on commonly used DDR2
memory modules.
(iii) Cache size: The 1024-line cache proved to be too small to achieve good performance
when rendering large scenes. This was especially true for the Sponza Atrium test
scene, which took up almost all of the 32 MB of on-board memory. Cache hit rates
below 10% made the use of a cache practically useless.
With the benefit of hindsight, a different board configuration would have eliminated these
problems. A PCI Express board with an integrated DDR2 memory module socket would
have eliminated the bandwidth problems. Furthermore, a FPGA with more on-chip SRAM
resources or a series of on-board SRAM chips would have reduced the cache thrashingwhen using very large scenes.
Since the project is written in generic VHDL, it can easily be ported to different target
hardware. The following section analyzes the projected performance of the co-processor
when the bandwidth and cache issues are not the limiting factors. While these conditions
are idealized, they are realistic in conjunction with the proposed target hardware and showwhat could be done with proper funding.

69

9 R

(a) Rendering of a caustic(b) Glossy reflection(c) Sponza Atrium(d) Cornell Box

Figure 9.1: Renderings of the used test scenes70

9 R

9.1 P 

9.1 Projected performance
When using a 128-bit memory connection using double data rate I/O, a 384-bit fetch takes
only two clock cycles. To calculate the expected memory access timeout, the worst case wasassumed: Every memory access references an inactive DDR memory row, which results ina “pre-charge” and “activate” command including two no-operation cycles before actually
being able to access the data. Thus, a total of six clock cycles is needed for any external
memory access. Assuming a cache hit ratio of 90%, the expected memory access latency is
0.9 · 1 + 0.1 · 6 = 1.5 cycles

To calculate the projected performance, the average number of kd-tree node and triangle
data accesses per intersection operation were measured (Table 9.1).
Scene name

Number of triangles

Avg. kd-tree node accesses

Avg. triangle data accesses

Caustic
Cornell Box
Sponza Atrium

130
2,270
76,136

6
10
41

4
9
29

Table 9.1: Average number of accesses per intersection operation
These values were then multiplied by the memory access latency to calculate the average
number of clock cycles required for an intersection operation in these scenes. Based onthese numbers, the maximum number of intersections per second at 100Mhz was calculated for different core configurations (Table 9.2).
Scene name
Caustic
Cornell Box
Sponza Atrium

Clock cycles per intersection

1 Core3 Cores

8 Cores

16.5
28.5
105

6,060,606
3,508,771
952,380

18,181,818
10,526,316
2,857,142

48,484,848
28,070,175
7,619,047

Table 9.2: Intersections per second using different core configurations
Further improvements could be achieved by increasing the core’s clock speed beyond
100 Mhz – for example by performing critical path optimizations or switching to a higher
speed grade FPGA device. Implementing this core using an ASIC1 process would also dramatically increase its speed.

9.2 Conclusion and future work
A highly space-efficient processor architecture using floating point arithmetic was presented. The recursive kd-tree intersection code was mapped into an efficient, pipelined
hardware design. A multi-core version of this architecture was implemented on an FPGA
1Application-specific integrated circuit - A custom, highly efficient, non-reconfigurable integrated circuit.

719.2 C   9 R

and both verified using test-benches and a specially designed global illumination renderer.
While the performance of the prototype did not meet the initial expectations, much improvement will be possible using more adequate evaluation hardware. A custom floating
point core design made it possible to fit the exceptional number of 126 floating point cores
onto a single FPGA. The designed architecture proved to be highly efficient in the numberof required clock cycles per intersection operation.
The single most important future work will be to acquire more adequate evaluation
hardware, and to adapt the project to it. With some effort, it might be possible to design
a hardware implementation of the ray generator, which runs directly on the FPGA – thus
significantly lowering the communication bandwidth requirements. Additional areas of
future work include such optimizations as mailboxing2 or studies on how variations inpipelining and cache design affect the performance of the intersection core.

2

Mailboxing is a technique which avoids unnecessary intersection tests when a triangle strikes multiple kdtree nodes.

72Bibliography[1] A, J. R. Analytic Methods for Simulated Light Transport. PhD thesis, Yale University,
CT, USA, 1995.
[2] D, P., B, P.,  B, K. Advanced Global Illumination. A K Peters,
Natick, MA, USA, 2003.
[3] K, J. T. The rendering equation. In SIGGRAPH (1986), pp. 143–150.
[4] L, E. Mathematical Models and Monte Carlo Algorithms for Physically Based
Rendering. PhD thesis, Katholieke Universiteit Leuven, Belgium, 1996.
[5] L, E.,  W, Y. Using the modified phong reflectance model forphysically based rendering. Tech. Rep. CW 197, Department of Computing Science,
Katholieke Universiteit Leuven, 1994.
[6] M, T.,  T, B. Fast, minimum storage ray-triangle intersection.
Journal of Graphics Tools 2, 1 (1997).
[7] P, B. Computer Arithmetic: Algorithms and Hardware Designs. Oxford University Press, New York, NY, USA, 2000, pp. 211–221, 279–307.
[8] P, M.,  H, G. Physically Based Rendering: From Theory to Implementation. Morgan Kaufmann, San Francisco, CA, USA, 2004.
[9] S, J., W, S., W, D., P, W. J.,  S, P. Realtime ray
tracing of dynamic scenes on an FPGA chip. In HWWS ’04: Proceedings of the ACM
SIGGRAPH/EUROGRAPHICS conference on Graphics hardware (New York, NY, USA,
2004), ACM Press, pp. 95–106.
[10] V, E. Robust Monte Carlo Methods for Light Transport Simulation. PhD thesis,
Stanford University, CA, USA, 1997.
[11] W, I.,  H, V. On building fast kd-trees for ray tracing, and on doingthat in O(N log N). In Proceedings of the 2006 IEEE Symposium on Interactive Ray
Tracing (2006), pp. 61–69.
[12] X, I. Virtex-4 Family Overview, v1.5 ed., Feb. 2006. http://www.xilinx.
com.
[13] X, I. Virtex-4 User Guide, v1.5 ed., Mar. 2006. http://www.xilinx.com.

73

B

B

[14] X, I. XtremeDSP for Virtex-4 FPGAs User Guide, v2.2 ed., July 2006. http://www.xilinx.com.

74