Modelling and Rendering of Teeth A Study in Sub Surface Scattering

Christian Thode Larsen

Kongens Lyngby 2011 IMM-M.Sc.-2011-31

Technical University of Denmark Informatics and Mathematical Modelling Building 321, DK-2800 Kongens Lyngby, Denmark Phone +45 45253351, Fax +45 45882673 [email protected] www.imm.dtu.dk

IMM-M.Sc.: ISSN 0909-3192

Summary

This thesis presents an appearance model for teeth. The model is based on an in depth study of subsurface scattering and diffusion theory. Theoretical components that have relevance for the model are presented. The physical components of teeth that have significance for image synthesis are also covered. Based on the appearance model for teeth, a rendering solution has been designed. The solution is able to render images of teeth in real-time (interactive frame rates on a Nvidia GTX 580) and without precomputation. The components of the solution and how they relate to the studied theory are covered. Results of image synthesis achieved by implementing the solution is presented. The advantages and disadvantages of the model and solution, found in the process of making these results, are discussed. Identified potential for further work is also presented.

ii

Acknowledgements

This project has not only required an in depth study of subsurface scattering theory, but also a study of diffusion. Together, they form a very broad area of knowledge. I would like to thank DTU supervisors Allan Peter Engsig-Karup, Jeppe Revall Frisvad and Jakob Andreas Bærentzen for their continued input and guidance. They all proved to be invaluable in helping me advance up a very steep learning curve. Without their help, I would not have been able to cover such a large area of knowledge. I spent a greal deal of time working on the project at 3shape, and I would like to thank them for providing a nice working environment that consists of very friendly people. Feeling welcome there proved to be a great motivation for my continued efforts in the project. I would especially like to thank Peter Dahl Ejby Jensen, the technical supervisor assigned to me from 3shape. Peter showed an amazing drive and interest in facilitating the project. His knowledge and supervision made it much easier for me to focus on the academic content and the assignment at hand. Peter especially proved to be invaluable during the periods where progress was slow and where results seemed far away. Without his support, also during out-of-office hours, and without his continued belief in my abilities, the project would likely have developed quite differently. Working on a master thesis is hard work, and it has required a lot of my attention within the last eight months. Thanks goes to family and friends, who always show a great deal of understanding and support when I am occupied with work. Finally I would like to thank, in no particular order, coffee, TV-series, computer games, books, karate practice and parties. Sometimes you just have to lay back, relax and forget about work.

iv

Contents

Summary

i

Acknowledgements

iii

1 Introduction

1

2 Structure of Teeth

5

3 Propagation of light 3.1 Radiance . . . . . . . . 3.2 Geometrical Optics . . . 3.3 Radiative Transfer . . . 3.4 Scattering Components

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

11 12 15 20 30

4 Diffusion 39 4.1 Diffusion Approximation . . . . . . . . . . . . . . . . . . . . . . . 40 4.2 Diffusion Profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.3 Interlayer Scattering . . . . . . . . . . . . . . . . . . . . . . . . . 58 5 Modelling and Rendering 5.1 Rendering Components . . . . . . 5.2 Parameter Analysis . . . . . . . . . 5.3 Diffusion Profiles for Dentine . . . 5.4 Sampling the Diffusion Profile . . . 5.5 Dentine Geometry . . . . . . . . . 5.6 Single Scattering . . . . . . . . . . 5.7 Transmission of Back Illumination 5.8 The Shader Algorithm . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

63 64 65 68 69 69 71 74 74

vi

CONTENTS

6 Results 6.1 Dentine Diffusion Profiles . . . . . . . . . . 6.2 Shader Components . . . . . . . . . . . . . 6.3 12-tap BSSRDF sampling/BRDF . . . . . . 6.4 Roughness and Wetness . . . . . . . . . . . 6.5 Enamel and Dentine Scattering Parameters 6.6 Enamel Thickness . . . . . . . . . . . . . . 6.7 Translucency . . . . . . . . . . . . . . . . . 6.8 Real Teeth/Crowns . . . . . . . . . . . . . . 6.9 Sets of teeth . . . . . . . . . . . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

77 77 81 84 86 88 88 91 91 97

7 Discussion 7.1 Scattering in Dentine (and Pulp) . . . . 7.2 12-tap sampling . . . . . . . . . . . . . . 7.3 Scattering in Enamel . . . . . . . . . . . 7.4 Enamel Thickness, Dentine Layer Model 7.5 Transmission of Backlight . . . . . . . . 7.6 Absorption and Scattering Parameters . 7.7 Age and Wearing . . . . . . . . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

101 101 103 103 104 105 106 107

8 Further Work 8.1 Predictive Rendering . . . . . . . 8.2 Dentine Layer Modelling . . . . . 8.3 Internal Reflection . . . . . . . . 8.4 Enamel-Dentine Layer Noise . . . 8.5 Parameters . . . . . . . . . . . . 8.6 Scattering in Dentine . . . . . . . 8.7 Single Scattering Approximations 8.8 Texturing . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

109 109 110 110 111 111 112 113 113

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

9 Conclusion

115

Index of Notation

119

Chapter

1 Introduction

This thesis presents an apperaance model for teeth, as well as a rendering solution that is able to produce synthesised images of teeth in real-time. The work is a result of a master project proposed by 3shape, a company that specializes in scanners and software for dental products. The project was proposed due to a need for improving their visualization of teeth, a visualization intended to run in real-time and preferably without a dependency on any major precomputation step or artistic work. The need was further bolstered by the current lack of appearance models for the particulate material. The purpose of an appearance model is to capture the characteristics of a material, or a group of similar materials, as they appear in real-world photographs. In the context of teeth, no existing appearance model captures their very special characteristics well, thereby motivating the need for a better model and consequently the possibility of performing better image synthesis. To date, published work on appearance models for teeth are very limited. Only recently was a layered model for rendering of human teeth presented [Shetty 2011]. However, the presented model depends on artistic work in order to model one of the physical components of teeth, and the presented work in general did not meet the requirements of this project. Teeth are complex and are composed of several layers of different materials that

2

Introduction

exhibits translucency, a phenomenon that occurs as a result of light propagation inside the materials. Specifically, light enters and emerges at different points due to subsurface transport of light, commonly referred to as subsurface scattering. When light propagates between several translucent materials of different properties, or alternatively a single translucent material with variations in its properties, the transport of light also results in color bleeding. One way to create an appearance model is to describe the propagation of light, its interaction with objects and the materials they are composed of, with models that are physically correct. Such models closely relates to the area of predictive rendering, where high accuracy in all aspects of the rendering process (physically based appearance model, rendering technique and its display) ensures that it is possible to predict the result – the appearance of the rendered material [Wilkie et al. 2009]. Another way to synthesise images is to employ models, given one or several mathematical expressions based solely on preference or observation, that produces a convincing visual result. This is referred to as believable rendering. The model for teeth presented in this thesis has its roots firmly placed in physically based models for subsurface scattering of light. This choice, resembling the choices that are necessary for predictive rendering, is motivated by the general idea that physically based models should provide the best tool set for good and accurate image synthesis of teeth. This should be seen in contrast to a choice closer to believable rendering, which is assumed likely to depend on a high degree of artistic modelling. It was previously mentioned how all rendering aspects must be accurate in order to make predictive rendering. The rendering step in predictive rendering requires adequate sampling of the physically based models that simulate light-material interaction, which is typically done by Monte Carlo- or finite element-methods. This type of rendering is generally too costly to perform in real-time except for very simple materials. The rendering solution presented in this thesis is consequently one of believable rendering, because it is necessary to employ simplifications and effect tweaks in order to maintain real-time compability. But the solution can also be thought of as a hybrid between the areas of predictive and believable rendering, because it models the physical light-material interaction much more closely than other believable rendering models. This small introduction to the areas of predictive and believable rendering is intended to make the reader familiar with the train of thoughts and the line of reasoning that lead to, and motivated, the work presented in this thesis. It is also intended to illustrate how image quality and performance (in the form of maintaining interactive frame-rates and absence of precomputation) are two factors that oppose each other. In order to improve image quality it is generally necessary to improve the underlying models that describe propagation of light

3 and/or to increase their sampling. This in turn is likely to increase the cost of each image, and may penalize performance. Likewise, in order to improve frame-rates it is often required that the implemented models are approximated or simplified, which in turn may have a notable effect on the resulting images. The solution proposed in thesis is thought to be a good compromise between image quality and performance, as it captures the primary traits of teeth while maintaining interactive frame rates. This thesis has been clearly partitioned in chapters that aims to familiarise the reader with the concepts that must be studied in order to implement the proposed shader, or alternatively to design ones own. The reader will initially be familiarized with the components of teeth that are believed to be main contributors to their appearance. From this follows the theory and relevant approximations that can be used to model light interaction with said components. This leads to a presentation of the shader solution and the key components that it contains. With the shader presented, the results found during the performed study of subsurface scattering, teeth and implementation of the shader are presented. Finally these results are discussed, and the identified potentials for further work are covered.

4

Introduction

Chapter

2 Structure of Teeth

In this (brief) chapter the structure and composition of teeth will be described, with emphasis on the components that plays a role in their visual appearance. Since the reader has not yet been familiarized with related computer graphics concepts, physical models for light and the parameters they rely on, the chapter will be (fairly) non-technical. These models will be covered in chapter 3 and 4, and details about parameters for teeth will be disclosed in chapter 5 where the proposed solution to an appearance model and real-time based shader will be presented. Human teeth are categorised in two groups, anterior teeth (incisor and canine), the front teeth and the four sharp teeth that are used to rip food apart, and molar teeth, positioned further into the mouth and used to chew food. As can be seen in figure 2.1 on the following page1 , teeth consist of several layers. Of these layers, the enamel and dentine layers are believed to be major contributors to the visual appearance of teeth, while any contribution from the pulp may be neglected. The work in this thesis is based on this assumption for two reasons. The first reason is that the dental references located through the course of this master project generally seem to focus entirely on the propagation of light in enamel 1 From

http://en.wikipedia.org/wiki/Tooth_(human)

6

Structure of Teeth

Figure 2.1: The different layers in teeth illustrated for a human molar tooth. The properties of the enamel and dentine layers are key components in creating an appearance model for teeth, since light is able to propagate in both.

7

Figure 2.2: A sliced human incisor tooth. The thickness of the enamel layer varies from bottom to to top, from 0.5 mm to 1.5 mm respectively.

and dentine, which implies that contributions from other components of the teeth are minor at best. From a general point of view this seems to be a valid assumption, since the enamel and dentine layers are the two outermost layers of teeth. The second reason is that the available (limited) measurements for teeth (enamel and dentine) places a natural limit on the components of teeth that can be taken into account in the solution, if one relies solely on a study of publicised work. The alternative would be to measure the different components of teeth directly, but this is a time consuming process and would require an entirely different study. Returning to the composition of teeth, the enamel layer is mostly composed of hydroxyapatite crystals, together with water and proteins. The hydroxyapatite crystals are oriented in keyhole shaped prisms [Zijp et al. 1995]. This structure is the reason why the enamel layer is hard, and serves to protect the teeth from damage when eating while it makes it easier to break the food into pieces at the same time. Light is able to propagate in the enamel layer [Spitzer and Bosch 1975; Fried et al. 1995; Zijp et al. 1995], which consequently makes it a key component in modelling the appearance of teeth, since light always passes through this layer before any interaction with other layers. The thickness of the enamel layer is known to vary quite a lot. It depends on the type of tooth (incisor or molar), sex, and where on the tooth the thickness is measured [Alvesalo and Tammisalo 1981]. This is also apparent for a molar tooth in figure 2.1 on the preceding page. Specifically, the enamel layer is thinnest near the gums down to 0.5 mm, and becomes thicker as it gets closer

8

Structure of Teeth

Figure 2.3: An anterior tooth (incisor or canine). The enamel and dentine layers have been revealed due to grinding of teeth, which provides an insight into the thickness of the enamel and dentine layers.

to the top. Molar teeth are known to be as thick as up to 2.5 mm near the cusps, while incisor teeth may be as thick 1.5 mm near the top [Alvesalo and Tammisalo 1981; Ten Cate 1985]. An illustration of a human incisor tooth with a ruler aligned can be seen in figure 2.2 on the preceding page2 . The dentine layer works as a supporting layer for enamel. It is not as hard as the enamel layer and also contains hydroxyapatite crystals [Zijp and Ten Bosch 1993]. Furthermore it contains tubules, cylindrical rods that are oriented along the dentine normals outward from the pulp towards the interface towards the enamel layer. Like the enamel layer, the dentine layer also propagate light, which is known to be affected by the orientation of the dentine tubules [Zijp and Ten Bosch 1993; Vaarkamp et al. 1995; Fried et al. 1995; Kienle et al. 2003; Kienle et al. 2006]. The dentine layer is sensitive to physical stimuli, because it is connected to the pulp by cells called odontoblasts that extend into the tubules. Together with the enamel layer, it is believed to be a key contributor to the real-world appearance of teeth. No measurements of dentine thickness has been found in this project. However, 2 From

http://en.wikipedia.org/wiki/Tooth_enamel

9

Figure 2.4: An illustration of the variation in translucent effect near the edges of incisor teeth. Translucency is a primary charecteristica of teeth. Another important trait is color bleeding, which occurs due to light that interacts with the dentine and enamel layers.

inspection of figure 2.3 on the preceding page3 that illustrates the cross section of a human molar tooth aligned with a ruler, suggests that it varies along the normal between approximately 1 mm and 3 mm. Inspection of figure 2.2 on page 7 seems to support this claim for an incisor tooth. The pulp is composed by tissue, blood vessels and nerves, and connects the tooth to the blood vesssels and nerve system of the body. The previously mentioned odontoblasts that extends from the pulp serve to produce new dentine and to feed the biological material with nutrients as well as moisture. The pulp will be neglected in the work presented in this thesis, in the sense that it is assumed either to absorb any entering light, or alternatively that no light reaches the pulp layer. Figure 2.44 illustrates the primary characteristics of teeth, specifically the semi transparent effect (translucency) and color bleeding. These effects are the result of the light interaction between the different layers as well as their individual properties, and capturing these effects (in some way) is fairly essential in order to achieve a believable tooth image. Note that these visual effects are the result of the interaction between light and the pure tooth layers, and that they do not include factors such as smooking or bad healthcare, which affects how teeth 3 From 4 From

http://en.wikipedia.org/wiki/Tooth_enamel http://www.drgregkerbel.com/Dental-Smile-Catalog-Photos.html.

10

Structure of Teeth

appear. Modelling such factors require specific attention, and it is not a focus in this thesis, as the goal is to model teeth that looks nice. With this brief introduction to the composition of teeth, it is time to cover physically based models for propagation of light.

Chapter

3 Propagation of light

Appearance modelling is a subject that aims to produce models that can be used to make convincing image synthesis of objects in the real-world. This is a question of how the object is perceived by the human eye, which is affected by the shape of the object and by its color. While the shape plays an important role in this perception, the focus of this chapter is why different objects have different colors, and how this is accurately modelled. The color of an object appears as a result of interaction between light and a material (or several), which in turn is a result of the intensity of the light that hits the object and the properties and composition of the material. Slight variations in these properties can be very important for the perception of the object, and this illustrates why it is important that models that are used for realistic image synthesis are able to accurately capture these variations. In order to do this, it is necessary with mathematical formulations of the propagation of light as well as abstract parameters that can be used to represent the properties of any material. One large set of materials are those that exhibit translucency. This is a phenomenon that occurs due to something called subsurface scattering. Specifically, light enters at one point on the surface and emerges at another due to interaction with the particles of the medium. This is commonly referred to as radiative transfer and subsurface transport of light, and the theory on this subject presents the models that are needed in order to capture the appearance of many real-world materials.

12

Propagation of light

Figure 3.1: The electromagnetic spectrum. Visible light lies in the interval of wavelengths between approximately 400 and 700 nanometers.

The theory on radiative transfer is a complicated subject that is closely related to quantum mechanics, electromagnetism and fundamental physics. In order to get into the details of a subject that has such a large base of scientific knowledge, some sort of starting point is needed. A base of understanding from where it becomes possible to delve further into the various details and complications of the subject, specifically the mathematical models, their advantages as well as the assumptions about light and and materials that they rely on. In order to form this base, this chapter will begin by defining the most basic quantity that are involved in radiative transfer, namely radiance.

3.1

Radiance

Before the term radiance is defined, it makes sense to begin with the term light, which should be much more familiar to most people. In a somewhat philosophical formulation, light is used intuitively to express that which has to be present in order for us to be able to see things around us, as opposed to darkness, a term that is used to refer to its absence. Light is also often referred to as beams of light, which are created because the sun shines, when a candle is lit or when a lamp is turned on. From a more practical point of view, light is an electro magnetic force. More precisely, light is electromagnetic radiation, hence the term radiance. This radiation consist of waves of energy that are emitted from sources of light, such as the sun or a lamp, and travels at approximately 3 · 108 meters per second in

3.1 Radiance

13

Figure 3.2: An illustration of the additive RGB color model, and how adding red, green and blue together forms the various colors of the visible electromagnetic spectrum.

vacuum. The radiated waves comes at different wavelengths within the electromagnetic spectrum as shown in figure 3.1 on the facing page1 . The human eye is able to perceive the difference between wavelengths in the visible spectrum, and this is the reason why different colors exist. The visible spectrum lies in a range between approximately 400 and 700 nano meters, and as a rule of thumb, a wavelength of 450 can be considered to be be a fully blue color, 525 fully green and 700 fully red. Perhaps also familiar to most people are the non-visible ultraviolet and infrared wavelengths, which lies below and above the wavelengths of the visible spectrum, respectively. In computer graphics, the visible spectrum is normally approximated through blending of the three colors red, green and blue as defined by the additive RGB color model, which can be observed in figure 3.22 . By mixing different intensities of these three colors, it is possible to approximate most of the colors that the human eye is able to perceive. If the intensity of all three colors are zero, it corresponds to black, and if the intensities are one, to white. This is a computationally cheap way to represent the spectrum of colors, and much more simple than simulating light that is emitted at a multitude of different wavelengths. Also, for most purposes blending between three channels is sufficient for visual realism. The behaviour of light light does not only resemble the behaviour of waves, but also that of particles. These particles consist of energy and are referred to as photons. In total the theory behind this particle-wave behaviour (the particle1 From 2 From

http://en.wikipedia.org/wiki/Light http://en.wikipedia.org/wiki/RGB_color_model

14

Propagation of light

Figure 3.3: An illustration of one solid angle (steradian), a two dimensional span given by the square of the radius of sphere in three dimensions. One solid angle corresponds to a fractional area of 1/4π on the sphere.

wave duality) and light-matter interaction is called quantum electro dynamics, and forms the basis for all mechanical, electrical and chemical laws [Feynman et al. 1964]. With a general picture of what light is, focus can be turned to radiance. Radiance is measured (in SI units) as W sr−1 m−2 , watts per steradian per square meter [Ishimaru 1978]. A steradian is an SI unit that describes solid angle, a two-dimensional angular span on a unit sphere in three dimensions, centred around a specific direction w ~ towards a point of interest on a (possibly imaginary) surface a. A steradian has been illustrated by figure 3.33 , and it follows that the size of the angular span depends on the distance to the center point of the sphere, or in other words its radius. In other words, radiance is the intensity of energy W that moves within one solid angle sr−1 through a unit area m−2 on a surface. Radiance is also referred to as the specific intensity, or the average power flux density, both of which signifies that the energy moves in a specific direction. If the radiance is wavelength dependent, it is referred to as spectral radiance, which is measured as W sr−1 m−2 Hz −1 , watts per steradian per square meter per hertz. As previously mentioned, radiance is directional, a flow of energy within one solid angle, e.g. on the unit hemisphere, while the terms irradiance and radiant emittance (or exitance) are used to describe the flux density of light that enters or exits, respectively, a material over the entire hemisphere at a point on the surface. These two quantities are measured in watts per square meter (W m−2 ), or alternatively W m−2 Hz −1 if wavelength dependent. As was previously covered, wavelength dependency is often omitted in computer graphics, for which reason spectral radiance will be omitted in the following. 3 From

http://en.wikipedia.org/wiki/Solid_angle

3.2 Geometrical Optics

15

Radiance is a basic quantity in physically based rendering. It signifies the intensity of light that travels, and how this energy is distributed when light interacts with different materials. Modelling the travel of light and its material interactions is then the challenge, and appropriate models must be used. Initially it is necessary to describe how light behaves as it interacts with the surface of a material. One way of describing this is given by geometrical optics.

3.2

Geometrical Optics

Geometrical optics is a set of geometrical observations of how light, treated as rays, interacts with the surface (or interface) of objects that consist of different materials. In [Feynman et al. 1964], it is pointed out that geometrical optics is purely ”geometrical approximations intended to describe how light behaves on a large scale”, and that it does not take into consideration the theory of quantum electro dynamics, the fundamental theory that provides the full picture of what light really is (e.g. wave propagation). However, even though geometrical optics do not capture the full physical picture of light, treating light as rays has proven to be quite effective for the purpose of producing convincing image synthesis. A commonly used off-line technique that relies on geometrical optics it is path tracing, where rays initially are traced from the camera into the virtual scene. At some point a ray hits an object in the scene, and the portion of light that will travel towards the origin of the ray (the camera) at that specific point on the object is computed, as well as a new path for that particular ray, given the shape and properties of the object it hit. This process then repeats itself recursively until it is terminated probabilistically (using Russian Roulette). Most real-time based shader solutions also relies on geometrical optics, with the exception that the number of material interactions for a single ray is limited in order to maintain speed. In both cases (path tracing and real-time shaders) the path that the rays follow after interaction with the surface of the object material is given by the definitions presented in this section. The first observation is that of reflection, which can be observed when light is reflected by a completely smooth surface, such as a mirror or a still pool of water. It has been illustrated by figure 3.4 on the next page and is given by [Feynman et al. 1964]

θi = θr .

(3.1)

16

Propagation of light

~n

ω ~i

θi

θr

ω ~r

Figure 3.4: An illustration of the reflection of light. A beam of light with direction ω ~ i hits a mirror and is then reflected along direction ω ~ r . The angles of incidence θi and reflection θr are equal, relative to the normal ~n of the mirror surface. The reflected direction lies in the plane of incidence.

The relation simply tells us that the angle of light θi incident on a smooth surface is equal to the angle of the reflected light θr in the plane of incidence. Note that in this context, terms such as incident angle or simply angle are commonly used as implicit references to angles formed by some ray of light and the surface normal. The next relation is observed when light refracts through the interface of a material. The light penetrates the surface object, and follows a new direction in the material it has entered due to the refraction. If the object is closed, the light may then refract through the interface from the inside. Alternatively it is internally reflected, and continue to travel inside the object. Refraction is modelled by Snells law, it can be observed in figure 3.5 on the facing page and is given by [Feynman et al. 1964]

sin(θi )ηi = sin(θt )ηt .

(3.2)

The relation says that the product of the sine of the incident angle sin(θi ) and the refractive index ηi of the first material is equal to the sine of the refracted angle sin(θt ) and the refractive index of the second material ηt . As with reflection, the refracted direction lies in the plane of incidence. Note that the refracted angle is the angle between the direction of the refracted light, and the inward surface normal, not the outward normal. The refractive index is used as a measure of how light refracts between different materials. Specifically it is a measure of how the light changes speed as it cross the interface between different materials [Born and Wolf 2003]. Examples of refractive indices are that of air (an index of 1.0) or water (1.33).

3.2 Geometrical Optics

17

interface ηi ηt ω ~i θi ~n

θt

-~n ω ~t

Figure 3.5: An illustration of the refraction of light. A beam of light with direction ω ~ i hits an interface and is then refracted along direction ω ~ t in the plane of incidence. The angle of refraction θt is determined according to Snells Law, which takes into account the refractive indices ηi and ηt of the two materials as well as as the angle θi incident on the interface between them(the surface of the material the light enters).

18

Propagation of light

It can be seen from the relation (or by experimentation) that as the incident angle on the interface becomes more and more glancing on the surface, the refracted angle increases. It also follows that as the index of the material light refracts into becomes higher, relative to the material light arrives from, the refracted angle will be much smaller, i.e. the refracted light will follow a direction that is closer to the inward surface normal. It follows from this case (light that moves from a high to a low index of refraction) that there is a set of directions that the light will never follow after refraction. The opposite case is particularly important, i.e. when the incident material has a higher index of refraction than the refracted material. In this case total internal reflection may occur at or near incident glancing angles, in which case no light refracts across the interface. This is important when the path of light is simulated, and is most easily tested by evaluation of the sine or cosine to the refracted angle that has been computed by Snells Law. If the sine or cosine is higher than 1, or below 0, respectively, the angle is actually one of reflection, and no refraction has occurred. Previously it was stated that light either reflects, or refracts when it hits interface of a material. Which of the two depends on the probability for given action to occur, something which can be determined by calculating Fresnel reflectance. Alternatively, the Fresnel reflectance can be viewed as fraction of light that is reflected, and one minus the Fresnel reflectance as fraction of light that is refracted across the interface.

the the the the the

The Fresnel reflectance F (unit-less) is defined by [Born and Wolf 2003]

Fr =

Fk + F⊥ . 2

(3.3)

Here Fk is the reflectance of the parallel component of the incident electromagnetic wave, relative to the interface, and F⊥ is the reflectance of the perpendicular component. Note that this expression assumes unpolarised light; the light is reflected equally parallel and perpendicular to the interface. This assumption is normally made in the context of computer graphics, but if not, then the two components of the incident wave must be handled separately. Further details about polarisation of light can be found in [Born and Wolf 2003]. As seen in equation 3.3, in order to obtain the Fresnel reflectance at an interface, it is necessary to compute Fk and F⊥ , where the two are given by [Ishimaru 1978]

3.2 Geometrical Optics

19

 Fk =

cos(θi )ηt − cos(θt )ηi cos(θi )ηt + cos(θt )ηi

2

cos(θi )ηi − cos(θt )ηt cos(θi )ηi + cos(θt )ηt

2

(3.4)

and

 F⊥ =

.

(3.5)

Note how these equations utilize the cosine to the angle and not the sine. The cosine is easily found through trigonometry: sin(θ)2 + cos(θ)2 = 1. Given the Fresnel reflectance, the transmittance Ft (the fraction of light refracted across the interface) can now be found

Ft = 1 − Fr .

(3.6)

It was previously mentioned how refraction affects the direction of light that transmits across a surface between materials of mismatched refractive indices. This is a consequence of Snells law, and specifically has the consequence that the solid angle light moves within is compressed. This can be understood intuitively by considering the relation between the entire hemisphere of 2π solid angles of incident radiance on a surface, which is compressed upon refraction to a smaller area of solid angles on the opposite hemisphere. Conservation of energy states that energy incident at a point on a surface must equal the sum of reflected and transmitted energy. Accounting for energy conservation at all points on a surface, and adjusting refracted radiance according to the compression, the following relation can be found

 Li = Lr +

ηi ηt

2 Lt .

(3.7)

The relation shows that transmission into a medium of a higher refractive index results in an increase in radiance. Conversely, it is lowered upon transmission into a medium of a lower refractive index. A more complete description of compression of solid angle and how it is derived can be seen in [Ishimaru 1978; Frisvad 2008]. In practice, compression of solid angle tends to cancel out in simulations, as light that refracts into a material normally also refracts out again. Nevertheless, it is

20

Propagation of light

presented for completeness, as there may be cases where it remains necessary to include it during image synthesis. An example of this could be light that travels from a light source through air towards a camera placed below the surface of water. The previous example raises the question of how light behaves when it travels inside materials other than air or vacuum. Such materials that allow transport of light are referred to as participating media, as they participate in the light transport together with air. This phenomenon of light transport is closely tied with subsurface scattering, and finally brings the chapter to the actual topic of radiative transfer in and between participating media.

3.3

Radiative Transfer

Subsurface scattering is a phenomenon that can be observed in translucent materials, which generally produces an appearance that tends to be blurry or muddy. A few examples are milk, skin, blood, marble and whine. Translucency occurs as a result of transport of light inside the material. More specifically, light enters at one point on the surface due to refraction given equation 3.2 on page 16, and (possibly) emerges at another due to interaction with the particles of the material. Translucent materials are also referred to as participating media, a term used as a general description for any material that participates in the propagation of light together with air. As light travels in the medium, some of the radiance will be attenuated due to absorption. This depends on the absorption properties of the participating medium, which is modelled by the absorption coefficient σa . The absorption coefficient is a macroscopic parameter that is used to capture the generally observed absorption behaviour of a specific material. It is defined per meter, with an SI-unit of m−1 [Ishimaru 1978]. The absorption coefficient relates to the absorption cross section, which can be thought of as a two dimensional planar slice of a three dimensional particle that corresponds to a certain amount of observed absorption of radiance [Ishimaru 1978]. Put a bit differently, it couples the size of a particle to its absorption; as the size of the particle increases, so does the area of its cross section and thereby the amount of radiance that is absorbed by the particle. The absorption coefficient can also be thought of as the probability per unit length that light will be absorbed by the participating medium [Hanrahan and Krueger 1993]. The radiance is not only absorbed by the medium, it is also scattered as it

3.3 Radiative Transfer

21

hits the particles that reside in it. Similarly to the absorption coefficient, the scattering of light by a medium is defined by the scattering coefficient σs . The scattering coefficient also relates to the scattering cross section with an interpretation similar to that for absorption; the larger the particle is, the larger is its cross section and thereby the amount of scattered radiance. Similarly, it can be thought of as the probability of a scattering event per unit length. Note that it is important to distinguish between a scattering event in computer graphics and simulation of light that actually hits a particle. The former express the macroscopically observed scattering behaviour of the medium, and is an abstract model for scattering due to the scattering parameters and the phase function (covered in the following). The latter is what happens on microscopic level, which is computationally infeasible to model in computer graphics. The reason is, that a material easily consists of millions of particles. Together, the absorption coefficient and scattering coefficients express the total amount of attenuated radiance due to absorption and radiance that is scattered in other directions by the particles of the medium. This attenuation is represented by the extinction coefficient σt , given by [Ishimaru 1978]

σt = σa + σs ,

(3.8)

and it has definitions similar to those given for absorption and scattering. It is also called to total cross section, and as previously can be viewed as the total probability for a scattering or absorption event to take place per unit length. Related to the total cross section, or extinction coefficient, is the mean free path, which is abbreviated mf p. This is the mean distance light can be assumed to travel before an absorption or scattering event occur. The mean free path is defined by [Ishimaru 1978]

mf p =

1 . σt

(3.9)

This relation should be quite intuitive, since the extinction coefficient has an SI-unit of m−1 . Consequently, the mean free path is given by one over the probability per unit length for an event to happen. Another important definition is the albedo, which signifies the ratio between the scattering cross section and the sum of the cross sections for absorption and scattering. It is given by [Ishimaru 1978]

22

Propagation of light

ω ~0

ω ~0

ω ~0

Figure 3.6: Illustrations of isotropic, forward and backward scattering (respectively) of radiance by particles.

α=

σs . σt

(3.10)

The albedo expresses the relation between absorption and scattering. For highly scattering media where σs  σa , the albedo will have a value close to one, where as more absorbing media have albedos closer to zero. Intuitively it can be thought of as the probability for a scattering event to occur relative to the total probability of an absorption or scattering event taking place. The albedo is of particular importance for the validity of many approximations to the radiative transfer equation (equation 3.14) which will be described shortly. This importance of the albedo will be covered in further detail in chapter 4. The size, shape and orientation of the particles in the material are important, as they determine how the interacting light is scattered. For example, some materials may exhibit even scattering in all directions (isotropic scattering), while other materials scatter light predominantly in the forward or backward direction (anisotropic scattering), relative to the direction incident to the scattering event. More specifically, this is called directionally independent anisotropic scattering, as the orientation of the medium relative to the light incident on its surface does not affect the scattering. The scattering across the particle can be said to be rotationally invariant, the particle can be oriented in any direction, without it having any influence on the scattering. More complex particles may also depend on the orientation the medium and the light incident on its surface. The scattering behaviour of such particles is also anisotropic, but differs from forward and backward scattering due to the significance of orientation of medium. Consequently, the scattering is rotationally variant, and this scattering behaviour is referred to as directionally dependent anisotropic scattering. Isotropic, forward and backward scattering of radiance has been illustrated by figure 3.6. The scattering behaviour of a medium is macroscopically modelled by a phase

3.3 Radiative Transfer

23

function that forms a probability distribution with unit sr−1 , the probability that light will scatter in a certain direction ω ~ at point x given the incident direction ω ~ 0 [Pharr and Humphreys 2004]. This function is generally formulated as p(x, ω ~ 0, ω ~ ) for directionally dependent scattering media. For directionally independent (isotropic, forward and backward) scattering media it can be modelled by a much more simple function p(µs ) that only takes into account the cosine of the scattering angle µs between the incident direction and the outgoing direction. It is important to understand, that the phase function models the observed scattering behaviour of the entire medium, and not that of the individual particles as they may vary in size, shape and orientation. It is also assumed that particles have some distance between them, which makes it possible to neglect inter-particle interactions at point x. Assuming the simplest possible shape of a particle, a sphere, this distance is required to a be at least a few radii of the spherical particle. From this it should be intuitively clear, that choosing an appropriate phase function for a given medium can be of big importance for image synthesis. The actual effects of choosing different phase functions in order to model different media is an area of study in itself, and only a couple of phase functions will be covered here. The simplest possible phase function is simply a constant that models isotropic scattering, and is given by

p=

1 . 4π

(3.11)

As it can be seen, the function does not depend on any parameters, it simply signifies an even distribution, or probability, of scattering in all directions (4π solid angles, i.e. a sphere). A somewhat more complicated phase function is the Henyey-Greenstein Function, which is given by [Henyey and Greenstein 1941] p(θs ) =

1 1 − g2 . 4π (1 + g − 2g cos(θs ))3/2

(3.12)

This function models the forward and back scattering behaviour of a medium, which takes into account the scattering angle θs between the incident and scattered direction, where its cosine is given by ω ~ ·ω ~ 0 , and the asymmetry parameter g, which is defined as the mean cosine of the scattering angle [Ishimaru 1978]

24

Propagation of light

Z g= 4π

(~ ω·ω ~ 0 )p(~ ω·ω ~ 0 )dω 0 .

(3.13)

Note how none of the listed phase functions have a directional dependency. For an example of a more complex phase function that does have a directional dependency, for the purpose of modelling scattering of light by hair and cloth, see [Jakob et al. 2010]. Appropriate choices of phase functions for the purpose of modelling teeth, as well as assumptions made for the specific work presented in this thesis, will be covered in chapter 5. The asymmetry parameter g plays an important role in models for subsurface scattering, especially for approximations that models multiple scattering. These approximations will be covered more in depth in chapter 4 Relevant definitions are now in place for an introduction to the radiative transfer equation, which models propagation of light in participating media. It is also commonly referred to as the volume rendering equation [Kajiya and Von Herzen 1984], and is defined by (for non-emissive materials) [Chandrasekhar 1960]

~ (~ ω · ∇)L(x, ω ~ ) = −σt (x)L(x, ω ~ ) + σs (x)

Z

p(x, ω ~ 0, ω ~ )L(x, ω ~ 0 )dω 0 .

(3.14)

4π

It follows that it is a differential equation, and that the solution is radiance L(x, ω ~ ) with direction ω ~ at point x. σt , σs and p(x, ω ~ 0, ω ~ ) were previously ~ covered. ∇ signifies the gradient of the solution to the equation, i.e. the direction that has the largest change in radiance. This direction is related to the direction of the solution radiance ω ~ through the dot product, which yields the magnitude of change in radiance in direction ω ~ . The integral over all incoming directions (solid angles) dω 0 Z

p(x, ω ~ 0, ω ~ )L(x, ω ~ 0 )dω 0 ,

4π

signifies radiance that travels in direction ω ~ due to inscattering of radiance, according to the probability given by the phase function p(x, ω ~ 0, ω ~ ), and L(x, ω ~ 0) is the quantity of inscattered radiance. In total the equation describes the change in radiance that travels in direction ω ~ at point x as the sum of attenuation of radiance due to outscattering and absorption (hence the negative sign), and irradiance that is scattered in direction

3.3 Radiative Transfer

25

ω ~ due to inscattering. As such it models how the radiative transfer at a given point in a material is affected by, and also affects, the radiative transfer at surrounding points due to scattering of light. Most materials are non-emissive, i.e. they do not contain light sources. Should a material be emissive, a volumetric light source term Q(x, ω ~ ) may be added to the equation. The radiative transfer equation is a fundamental model for radiative transfer, and it is generally difficult, if not impossible to determine a full solution analytically, unless certain assumptions are made about medium and incident light. Full solutions are normally determined by numerical methods, either by sampling (Monte Carlo based techniques such as path tracing [Pattanaik and Mudur 1993]) or by techniques where the problem is formulated as a system of linear equations (examples to this are finite element, finite difference or finite volume methods[LeVeque 2007]). The degree of accuracy of the solution depends on the chosen technique. However, methods like those that use sampling or finite element approaches are costly. For this reason many approximations exist to the radiative transfer equation. These approximations are made under certain assumptions about the radiative transfer, and the participating medium in question, in order to lower the cost of the image synthesis. This is of particular importance for shaders that are intended to work in real-time, and relevant approximations for this purpose are covered in chapter 4. The full picture of radiative transfer is now in place, which makes it possible to summarize the transport of light in a participating medium. Obviously, the light begins its travel at some light source. Under the assumption that no scattering events takes place in the material that contains the light source (e.g. air), the light then hits and refracts into a participating medium. If the medium is sufficiently thin, it is likely that the light is transmitted through the object and only undergoes attenuation due to outscattering and absorption (referred to as 0th order multiple scattering). The light then hits the interface of the object from the inside, and is refracted out into the air again (under the assumption that no internal reflection occurs). Alternatively, the light is scattered due to interaction with the particles of the medium, which is referred to as 1th order multiple scattering. The probability this to happen is related to the thickness of the medium by the mean free path, given by equation 3.9 on page 21, and is referred to as a scattering event (due to the probabilistic nature of scattering). The probability for scattering in a specific direction is given by the phase function. After the scattering event the process repeats itself until the light exits the

26

Propagation of light

medium. In the case of more scattering events taking place, the total number of scattering events is referred to as nth order multiple scattering. In total the radiance L that hits the eye from a specific point on a surface is the sum of radiance that is reflected by the surface towards the eye Lr , transmission of radiance through the medium towards the eye Lt , and radiance that exits the medium towards the eye after a single and multiple scattering events, Ls and Lm respectively. Determining radiance that enters and exits at different points on the surface due to scattering in the medium is generally modelled through the use of a Bidirectional Surface Scattering Reflection Distribution Function (BSSRDF). Determining this function involves finding partial or full solutions to the radiative transfer equation, either analytically or numerically. The BSSRDF will be covered in the following section. In some cases the BSSRDF can be simplified to a Bidirectional Reflectance Distribution Function (BRDF), under the assumption that the incident light does not vary over the surface. This simplification is described in section 3.3.2 on the facing page.

3.3.1

BSSRDF

The Bidirectional Surface Scattering Reflection Distribution Function (referred to as BSSRDF in the following) was proposed by [Nicodemus et al. 1977]. It relates the incident differential radiative flux dΦi at a point on a surface to the outgoing differential radiance dLo that exits at a different point on the surface due to subsurface scattering with the function S

S(xi , ω ~ i , xo , ω ~ o) =

dLo (xo , ω ~ o) , dΦi (xi , ω ~ i)

(3.15)

where S describes the radiative transfer in the medium between the two points xo and xi . The differential radiative flux is further given by [Nicodemus et al. 1977]

dΦi = L(xi , ω~i ) cos θi dωi da(xi ),

(3.16)

where the term cos θi is attenuation of light due to incident angle (projection of the light onto the surface plane) at the point xi with differential surface area da(xi ).

3.3 Radiative Transfer

27

Since S only describes the transport between a single incident and exiting direction, the total radiance that exits in direction ω ~ o at point xo on a surface depends on the light incident from all directions over all incident points on the surface of the medium [Nicodemus et al. 1977] Z Z L(xo , ω~o ) =

S(xi , ω~i , xo , ω~o )L(xi , ω~i ) cos θi dωi da(xi ). A

(3.17)

2π

The equation is fairly straight forward; the differential radiance that exits in direction ω ~ o at point xo is given by a double integral over light incident from all angles over all points on the surface a. The function S then describes the propagation of light inside the participating medium between the two points xi and xo . The orientation of the surface area a is given by the surface normal ~n. Note that the formulation of outgoing radiance assumes a planar surface. In the case of a non planar surface the normal ~n must be defined locally for each point xi . Determining a solution function S is not trivial, since radiance may pass directly through the medium without scattering, or alternatively scatter one or several times. As such, there is an infinite amount of paths that light may follow as it propagates through the medium. In practice the function S is formed by splitting it into multiple BSSRDF terms, where each of these approximates a specific chain of light interactions. The sum of these terms then form the total BSSRDF. When the BSSRDF has been found, it is possible to determine the radiance that travels towards the eye due to surface and subsurface interactions with different objects, by using techniques such as path tracing. In the case of path tracing, each direction from the eye into the scene is sampled several times based on Monte Carlo techniques, where probability is used to determine the path of light transport for each sample due to reflection, refraction and scattering.

3.3.2

BRDF

The Bidirectional Reflectance and Transmittance Distribution Functions (referred to as BRDF and BTDF in the following) describe the specific cases where outgoing radiance in a given direction ω ~ at point x on a surface can be evaluated at x without considering the entire surface. This is a valid assumption if the incident light does not vary over the surface, or alternatively if there is no subsurface scattering of light. The BRDF fr is defined (in a way similar to the BSSRDF) by [Nicodemus 1965]

28

Propagation of light

fr (x, ω~i , ω~o ) =

dL(x,~ωo ) . L(x, ω~i ) cos θi dωi

(3.18)

dL(x,~ωo ) . L(x, ω~i ) cos θi dωi

(3.19)

The BTDF has a similar definition

ft (x, ω~i , ω~o ) =

The BRDF and BTDF relates the radiance at point x moving in direction ω~i , attenuated due to the angle of incidence, to the differential radiance reflected or transmitted in a single outgoing direction ω ~ o . The total amount of radiance in direction ω ~ o is found by insertion and integration of the BRDF or BTDF in the rendering equation, which is given by [Kajiya 1986] Z L(x, ω~o ) =

fr (x, ω~i , ω~o )L(x, ω~i ) cos θi dωi .

(3.20)

2π

The rendering equation can be seen as a simplification to the volume rendering equation 3.14 on page 24, as it resembles the inscattering term, but is limited to the surface and with the phase function replaced by the BRDF or BTDF. The rendering equation states that the radiative transfer in direction ω~o at point x on the surface is equal to an integral over all radiance reflected according to fr from all incoming directions (2π solid angles) of radiance over the hemisphere. As such, the equation can be said to express the transport of light in air or vacuum and its interaction with the surface of non-translucent materials (such as metals) or highly transparent materials (such as glass). As was the case for the radiative transport equation 3.14 on page 24, methods for finding a solution to the rendering equation involves sampling methods such as path tracing that is based on Monte Carlo techniques [Pharr and Humphreys 2004]. Another method that produces nice results is radiosity, which utilises finite element methods to simulate exchange of radiance between surfaces of objects that are split into patches [Watt 1993]. For an entirely smooth surface that does not exhibit any translucency, the BRDF fr and BTDF ft are simply given by the product of the Fresnel term, as presented in section 3.2, and a delta function. For a material that is not entirely smooth, for example due to the orientation of microfacets in the surface [Pharr and Humphreys 2004], it is necessary to use more complex functions that account for the microfacet distribution. One example of a microfacet model is the

3.3 Radiative Transfer

29

~n ω ~0

Figure 3.7: Lamberts cosine Law. Radiance is reflected evenly in all directions by the surface. The amount of reflected radiance is given by the incident radiance attenuated by the cosine to the incident angle, relative to the surface normal.

Kelemen and Szirmay-Kalos BRDF [Kelemen and Szirmay-Kalos 2001], another is the Cook-Torrance BRDF [Cook and Torrance 1981]. As covered in section 3.3.1, The light that leaves the surface of a material after interaction is generally modelled as the sum of a reflectance term and terms that constitute the subsurface transport of light. The reflectance term is generally represented by a BRDF, such as the previously mentioned Fresnel reflectance or a more complicated micro-facet distribution model. This is normally sufficient, as it is not necessary to consider subsurface scattering for light that does not refract across the surface. For materials that exhibit very little or no scattering, a BTDF sufficiently models light that refracts through the interface and continues its travel. A popular example of this, and often used for the purpose of illustrating refraction of light, is refraction of light through a glass panel. In some cases the subsurface scattering term is also appropriately modelled by a BRDF, which for a Lambertian surface is given by ρd /π (the observed brightness is the constant, regardless on the position of the viewer). Here ρd is the diffuse reflectance, which depends on the absorption properties of the material. For a non-absorbing material the diffuse reflectance is simply one. This can be a fair approximation in the cases where distribution of incident illumination over the entire surface is constant, and if the subsurface transport of light (if any) can be assumed to distribute uniformly in all directions. The diffuse distribution of light has generally been modelled by Lamberts cosine law, which states that the radiance observed at a point on a surface is proportional to the area projected onto the surface by the light source. In other

30

Propagation of light

words, the observed radiance only depends on the incident angle of illumination. Lamberts cosine law is given by [Born and Wolf 2003]

Lo =

ρd cos θi L(ω~i ). π

(3.21)

Here θi denotes the angle of incidence between the incident direction of light and the surface normal, and the cosine term originates from the rendering equation. A couple of examples of models that utilise a simple diffuse term are the Modified Blinn-Phong BRDF [Lafortune and Willems 1994] and the Interfaced Lambertian BRDF [Simonot 2009]. Both models use a Lambertian BRDF to model any light contribution due to subsurface scattering, while the Interfaced Lambertian models an additional interface between air and the diffuse light with Fresnel reflectance and transmittance.

3.4

Scattering Components

With the definitions of the BSSRDF and BRDF covered, the theoretical groundwork is in place for a description of expressions for light-material interactions in a participating medium, which generally are separated into three categories; 0th , 1th and nth order multiple scattering. In the following, BRDFs will be presented for 0th and 1th order multiple scattering, based on derivations by [Ishimaru 1978] and [Hanrahan and Krueger 1993]. nth order multiple scattering will also be mentioned, but only briefly as it has no analytical solution. The absence of a simple expression for nth order multiple scattering then forms the basis for the following chapter about the diffusion approximation.

3.4.1

0th order Multiple Scattering

0th order multiple scattering of light, also referred to as transmission or the reduced intensity, is the specific case where radiance is attenuated due to absorption and outscattering by the medium, but where the direction of transfer is unaffected by the scattering. The amount of radiance that will be transmitted without a change to its direction directly relates to the thickness of the medium it is being transferred in. Intuitively, as the thickness of the medium increases, more radiance is attenuated due to scattering and absorption until the point where the transmitted amount can be neglected. This relation between the thickness of the medium and its scattering and absorption properties is called

3.4 Scattering Components

31

~n Layer 1

ω ~i

η1 η2 ω ~

Layer 2

~n

-~n η2 η3 ω ~t

Layer 3 -~n

Figure 3.8: Transmission of radiance across a participating medium. The light refracts from one medium into another, travels some distance and is then refracted into a different medium. The direction of transfer does not change due to scattering.

32

Propagation of light

the optical depth (also referred to as the optical thickness, or optical distance). It is given by [Ishimaru 1978] Z

z

σt dz.

τ (z) =

(3.22)

0

This formulation of the optical depth assumes that the extinction coefficient may vary along the direction z (the depth in the medium). If the extinction coefficient is constant along z, or alternatively constant along a number n of intervals i with length d in the direction z, it can be formulated as a sum of products between the coefficients σt (i) and the corresponding distances d(i)

τd =

n X

σt (i).

(3.23)

i=1

Transmission of radiance in a medium has been illustrated by figure 3.8 on the previous page, and is given analytically by [Ishimaru 1978] L(0) (~ ω , z) = e−τ /µ0 Li (z = 0).

(3.24)

here µ0 signifies the cosine to the angle of transmission relative to the inward surface normal of the participating medium, an optical thickness of τ in the direction of the normal, and Li (z = 0) the light incident on the surface at depth z = 0. Consequently, division of τ by µ0 yields the optical distance in the transmitted direction by geometric calculus under the assumption of a planeparallel medium. In practice this is often not the case, for which reason τd in the following signifies the optical depth along the transmitted direction, and not along the direction of the inward surface normal, thus leaving out the cosine term. Accounting for attenuation of radiance due to refraction across the border of the medium, the reduced intensity through the medium (volume) is given by [Ishimaru 1978; Hanrahan and Krueger 1993] L(0) (xo , ω~o ) = Ft12 Ft23 e−τd L(xi , ω~i ),

(3.25)

where τd = σt k xo − xi k, and Ft is Fresnel transmittance terms due to refraction. Under the assumption of constant illumination over the surface, and again assuming a plane-parallel medium, it is possible to simplify this equation with a BRDF, which is then given by

3.4 Scattering Components

33

~n Layer 1

ω ~i

η1

θ0

η2

ω ~0

Layer 2 θ θs

η2

ω ~

η3 Layer 3

-~n

ω ~t

Figure 3.9: An example illustration of single forward scattering of light. The cosine to the scattering angle θs is given by the dot product of the directions ω ~ 0 and ω ~ . The cosine to the scattering angle is always positive in the case of forward scattering.

L(0) (ω~o ) = Ft12 Ft23 e−τ /µ0 L(ω~i ).

(3.26)

Both expressions for reduced intensity can be implemented directly for the purpose of rendering with off-line techniques like path tracing, and the BRDF is particularly useful in a real-time shader based solution, as it is not necessary to calculate the intersection point xo between the direction of transmission and the opposite surface. This covers the case of direction transmission of radiance across a participating medium, and forms the basis for the case of single scattering.

3.4.2

1st Order Multiple Scattering

In the case of 1st order multiple scattering, also referred to as single scattering, the transport of light changes direction due to a scattering event in the participating medium. In a more intuitive formulation, consider incident light along

34

Propagation of light

~n Layer 1

ω ~t

ω ~i

η1 η2

θ0

Layer 2

ω ~0

θ

ω ~ θs

η2 Layer 3

η3

Figure 3.10: An illustration of single backward scattering of light. The case is similar to forward scattering, but light is scattered back towards the interface that it already crossed once. In the case of backward scattering, the cosine to the scattering angle is always negative.

direction ω ~ i that refracts across some interface into a participating medium. It travels some distance in direction ω ~ 0 , at which point is scatters into a different direction ω ~ due to interaction with the particles of the medium. This distance is given by the mean free path, the average distance light can be assumed to travel before an event occurs. The probability for the light to scatter in a given direction is given by the choice of phase function, e.g. the Henyey-Greenstein function listed by equation 3.12 on page 23. After the scattering event light transport continues in the new direction ω ~ , until it hits either the initial interface or the interface on the interface on the other side of the medium. Depending on the indices of refraction at the layers, the radiance may refract across the boundary or reflect internally. (1)

Generally, single scattering is divided into two cases; forward scattering L+ and (1) backward scattering L− of light. As with 0th -order multiple scattering there exist analytical expressions that describe both cases. Forward scattering, illustrated by figure 3.9 on the previous page, is the case where light continues in a direction towards the interface opposite the one that it initially crossed into the participating medium. This implies that the cosine to the scattering angle is always positive. [Hanrahan and Krueger 1993] presents analytical BRDF solutions for single

3.4 Scattering Components

35

forward scattering (and single backward scattering that follows), which is based on a few assumptions. The first assumption is, that the incident illumination is constant over the surface. This makes it possible to formulate single scattering in terms of a BRDF, since incident illumination from the same direction is the same at all points on the surface. This also means that the scattered light is affected only by the incident and exiting directions at the surface, and the thickness of the medium (which relates to the reduced intensity, as mentioned in section 3.4.1). The expression for forward scattering of light is given by [Ishimaru 1978; Hanrahan and Krueger 1993] (1)

L+ (ω~o ) = Ft12 Ft23 αp(µs )

 µ0  −τd /µ0 e − e−τd /µ L(ω~i ). µ0 − µ

(3.27)

Here µ0 and µ denote the cosines to the angles θ0 and θ before and after the scattering event correspondingly, as illustrated by figure 3.9 on page 33. The phase function is given by equation 3.12 on page 23, as the medium is assumed to exhibit directionally independent scattering of light. The cosine µs to the scattering angle θs is given by the dot product of the direction of transfer in the participating medium before and after the scattering event. In total the equation expresses the attenuation of light due the probability for a scattering event to occur (the albedo α), refraction across both interfaces of the medium (the Fresnel transmittance terms Ft ), scattering and absorption in the directions relative to the scattering event (the cosine and exponential terms) and the probability for light to scatter in the given direction (the phase function). Note that this formulation assumes the medium to be plane-parallel, which can be seen implicitly from the exponential terms where the optical depth is divided by µ0 and µ. In the case where this assumption is invalid, the optical depth must be computed directly for the incident and outgoing directions to the scattering event. Equation 3.27 is not robust in the case where the cosine to the incident angle corresponds to the cosine to refracted angle, in which case it produces a division by zero. It is necessary to test for this special case in practical implementations. In these cases a different expression should be used, derived by using L’Hospitals rule, and given by [Ishimaru 1978; Hanrahan and Krueger 1993] (1)

L+ (ω~o ) = Ft12 Ft23 αp(µs )

τd −τd /µ e L(ω~i ). µ

(3.28)

The formulation for single backward scattering is similar to that of forward scattering, in the sense that it relies on the same parameters. However, as

36

Propagation of light

the light is scattered towards the interface it already crossed once, the cosine to the scattering angle is always negative. Single backward scattering has been illustrated by figure 3.10 and is given by [Ishimaru 1978; Hanrahan and Krueger 1993] (1)

L− (ω~o ) = Ft12 Ft21 αp(µs )

 µ0  1 − e−τd (1/µ0 +1/µ) L(ω~i ). µ0 + µ

(3.29)

The different terms here have similar interpretations as those given for single forward scattering. As with the reduced intensity BRDF, the BRDF expressions for single forward and backward scattering is quite useful for real-time based shader solutions. If accuracy is very important, specifically if there is (even slight) variations in incident illumination and when the medium is not plane-parallel, the expressions may not be appropriate. In this case, the radiative transfer equation 3.14 on page 24 may be sampled directly with off-line techniques such as path-tracing. It it also important to note that the outgoing radiance for both forward and backward scattering is based on a single incident direction. In order to find the total contribution of radiance in direction ω~o , it is necessary to insert and integrate the BRDF over all incoming directions with the rendering equation 3.20. For a real-time rendering purpose, this integral requires either precomputation, or a reasonable approximation in terms of multiplication of the BRDF by a constant (with the added risk of breaking energy conservation).

3.4.3

N th Order Multiple Scattering

BSSRDF formulations for nth order multiple scattering, or simply multiple scattering, are generally impossible to make analytically, because they require that a full solution to the radiative transfer equation 3.14 on page 24 is found. Finding a solution through brute force Monte Carlo methods such as path tracing is a very costly simulation, as it requires a large amount of sampling and simulated scattering events in order for the solution to converge. For this reason methods that estimates light contributions due to multiple scattering have received a lot of attention, and many approximations have been made for this specific purpose. Some rely on assumptions such as a relatively constant incident illumination over the surface, and others that the medium is plane-parallel. Some rely on diffusion theory, where it is assumed that scattering of light becomes totally random after a few scattering events. For this reason

3.4 Scattering Components

Layer 1

ω ~i

37

ω ~t

η1 η2

Layer 2 η2 Layer 3

η3

Figure 3.11: Multiple scattering. Light changes its direction several times due to scattering events in the participating medium.

these approximations are covered separately in chapter 4. An illustration of multiple scattering can be seen in figure 3.11.

38

Propagation of light

Chapter

4 Diffusion

It was mentioned in chapter 3 how finding a full solution to the radiative transfer equation 3.14 on page 24 is costly, because it requires adequate sampling. Specifically, this is a question of simulating the scattering events in the participating medium, which becomes costly for nth order multiple scattering. This means that sampling enough light paths and scattering events is both time consuming and necessary to produce sufficiently converged images. This motivates finding cheaper (physically based) approximations to nth order multiple scattering that produces acceptable results. Approximations normally rely on one or several assumptions about the phenomenon or physical process that is being modelled. This imposes extra constraints on the usability of the approximation as compared to the actual model. This leaves two choices when the level of approximation increases. One is to respect the constraints of the approximation and use it appropriately, the other is to neglect the error that is introduced when the constraints are not met. Depending on the constraint, the second case may be acceptable, but it may also introduce significant error in the result. This chapter presents approximations based on the theory of diffusion, which forms the basis for most, if not all, current approximations to nth order multiple scattering.

40

Diffusion

Figure 4.1: Diffusion of particles in a glass of water. The particles are initially concentrated in a small area, but will distribute evenly into areas of lower concentration given enough time. The system is said to be in a steady state when the flow becomes constant. The concept is the same for the diffusion approximation to nth order multiple scattering.

4.1

Diffusion Approximation

To reiterate for the case of nth order multiple scattering, the problem lies in the computational cost due to the many scattering events. The approximation that reduces this cost is the diffusion approximation. Diffusion is best described by considering a volume of particles. Initially, the particles enter the volume, resulting in an area with a high concentration of particles, while the remaining volume has low concentration. Diffusion then describes the flow of these particles from the area of high concentration to those of lower concentration, until the volume at some point reaches a steady state. The volume is in a steady state when the flow of particles through the volume does not change over time. When, or if, a steady state is reached depends on the radiative flux across the boundary (surface) of the volume. In other words, if the flux constantly changes, i.e. the number of particles entering the volume varies, the flow of particles in the volume will never enter a steady state. This implies that a steady state is reached only when the flux is constant across the border, and after the particles in the volumes have had a sufficient amount of time to distribute in the volume. The process of diffusion has been illustrated for particles in a glass of water by figure 4.11 . Recall the differential nature of the radiative transfer equation 3.14 and its 1 From

http://en.wikipedia.org/wiki/Diffusion

4.1 Diffusion Approximation

41

description; a change in radiance is due to in- and out-scattering as well as absorption. The intuition of this description is is equally valid if made for particles (photons) instead of radiance, due to the particle-wave duality. Also, the description bears quite some resemblance to the description of diffusion. Under the assumption that photons (radiance) are able to participate in a diffusion process (travel from higher to lower areas of concentration), it becomes possible to relate diffusion theory and the radiative transfer equation (3.14), thereby allowing nth order multiple scattering in a participating medium to be modelled by diffusion. The diffusion approximation relies on a general observation that scattering of light in a highly scattering media becomes almost entirely random after sufficiently many scattering events [Ishimaru 1978]. From this statement follows two important observations of its accuracy. First, the medium must indeed be highly scattering, which is the case when σs  σa or more precisely when the albedo α is close to one. The meaning of these parameters and the relation to the radiative transfer equation was covered in section 3.3. Second, the medium must be thick enough to allow for a sufficient number of scattering events. The implication is that the diffusion approximation is inaccurate in a region close to the light source, as too few scattering events have taken place. In this context, the light source is intuitively interpreted as the point on the volume surface where light enters, but it can be any point in the volume that is treated as a source. These are the general constraints and assumptions that must be met in order to use diffusion theory to model nth order multiple scattering, the details will be covered more thoroughly later. It is now necessary to introduce another important parameter, the reduced scattering coefficient, given by [Ishimaru 1978]

σs0 = σs (1 − g).

(4.1)

Here g is the asymmetry parameter described by equation 3.13 on page 24. The intuition here is, that when nth order multiple scattering is modelled by diffusion, i.e. an isotropic propagation of photons, the fraction of light that is asymmetrically scattered, given by the asymmetry parameter, loses its significance due to the random motion of photons. In other words, it relates the scattering cross section of a particle that exhibits asymmetric scattering of light to the cross section of particle that scatters light isotropically with approximately the same result for nth order multiple scattering. Roughly put, it is a hack that forms a relation between a light model where the direction of scattered light has a significance, and the diffusion model where this significance can be ignored. Consequently, the reduced scattering coefficient is directly related

42

Diffusion

to the previously listed constraints as to when the diffusion approximation is valid. A rough rule of thumb given by [Jacques and Pogue 2008] is, that the diffusion approximation is accurate when σs0 /σa > 10 and preferably > 20. Another parameter that follows directly from the reduced scattering coefficient is the reduced extinction coefficient, given by σt0 = σa + σs0 ,

(4.2)

with similar definitions to those of the extinction coefficient given by equation 3.8 on page 21, only now transformed into the isotropic approximation. From this follows the reduced mean free path mf p0 for the isotropic approximation

mf p0 =

1 , σt0

which is of particular importance. Consider the effect of reducing the scattering coefficient by the asymmetry parameter. This effectively models a particle that exhibit less scattering of light, since (approximately) only the fraction of isotropically scattered light is considered. Consequently, the mean free path changes and becomes a direct measure of how thick the material must be, in order for the isotropic approximation to be valid. It follows, by intuition, that the material should be at least one mf p0 thick, preferably two or more. This makes a lot of sense. If the material is not this thick, there is not sufficient volume to allow for the scattering of photons to become random. Finally follows the reduced albedo, which has similar definitions to the albedo given by equation 3.10 on page 22

α0 =

σs0 . σt0

(4.3)

The diffusion approximation in itself describes propagation of light in an infinite medium, and does not take the boundary (volume surface) into account for semiinfinite or finite volumes. This requires specific boundary conditions, which will be covered in section 4.1.1. One derivation of the diffusion approximation given by [Ishimaru 1978] considers the 1st order scattering of radiance that has entered a semi-infinite volume as a reduced intensity source Qri . In the following, the major steps in the derivation will be presented. This is done to provide some intuitive notion to the relation between the diffusion approximation and the radiative transfer equation, and because it forms a base for the dipole and

4.1 Diffusion Approximation

43

multipole approximations which will be covered in section 4.2. [Keijzer et al. 1988; Jensen et al. 2001; Donner 2006] have been used as additional references for the derivation. It proves to be useful to reformulate the solution radiance L(x, ω ~ ). This is done by splitting L(x, ω ~ ) into two components; radiance that has been scattered at least once, the diffuse radiance Ld (x, ω ~ ), and reduced intensity radiance Lri (x, ω ~) that arrives at point x from the surface of the semi-infinite volume. L(x, ω ~ ) is given by [Ishimaru 1978]

L(x, ω ~ ) = Ld (x, ω ~ ) + Lri (x, ω ~ ).

(4.4)

Of these two components, Ld (x, ω ~ ) is inserted into the radiative transfer equation (3.14). The resulting expression is given by [Ishimaru 1978]

~ d (x, ω (~ ω · ∇)L ~ ) = −σt (x)Ld (x, ω ~ ) + σs (x)

Z

p(x, ω ~ 0, ω ~ )Ld (x, ω ~ 0 )dω 0 + Qri (x, ω ~ ),

4π

(4.5) where Qri is the previously mentioned reduced intensity source, given by [Ishimaru 1978] Z Qri = σs (x)

p(x, ω ~,ω ~ 0 )Lri (x, ω ~ 0 )dw0 .

(4.6)

4π

Note that Lri (x, ω ~ ) can be computed analytically according to equation 3.24 on page 32, and signifies radiance that has travelled some distance into a volume, attenuated due to absorption and outscattering, but not yet been scattered in a different direction. This forms the basis for the derivation. The derivation involves reformulating the solution to the radiative transfer equation (radiance) in terms of fluence, the total quantity of radiance that passes through a point within a unit time interval. The expression for fluence is found by integrating equation 4.5 over all directions, and is given by [Ishimaru 1978] ~ · E = −σa (x)φ(x) + Q0 (x). ∇

(4.7)

~ signifies the direction where the change in fluence is greatest, i.e. the direc∇ tion that has the largest radiative transfer. E describes vector irradiance, or

44

Diffusion

the diffuse flux vector, at point x, i.e. a net radiative flux due to inscattered radiance. It is given by [Ishimaru 1978] Z Ld (x, ω ~ )~ ω dw.

E=

(4.8)

4π

The divergence of the vector irradiance signifies the magnitude of change in fluence, which is equal to the sum of absorption −σa φ(x) and emission from source Q0 , given by [Ishimaru 1978] Z Q0 (x) =

Qri (x, ω ~ )dw.

(4.9)

4π

It follows that Q0 (x) is the net magnitude (i.e. fluence) of the outscattered reduced intensity light source. The next step is to formulate the diffuse radiance Ld (x, ω ~ ) in terms of the diffusion approximation, namely that the propagation of radiance in a highly scattering medium tends to become isotropic. The diffuse radiance is represented by a two term spherical harmonics expansion, given by [Ishimaru 1978]

Ld (x, ω ~) '

3 1 φ(x) + ω ~ · E. 4π 4π

(4.10)

Here φ signifies the solution fluence φ(x), and the diffuse vector flux E at point x. The dot product between the diffuse vector flux and the direction of interest ω ~ yields the fraction of the irradiance at point x that will transfer in direction ω ~. This term is essential for the diffusion approximation, as it signifies radiance that moves from areas of higher concentration to lower concentration. In total this formulation describes a radiative transfer due to diffusion; radiance that travels in direction ω ~ is due to the sum of an isotropic distribution of radiance in all directions (random motion) plus a fraction of irradiance that joins direction ω ~ due to different concentration levels in the volume. The details about the constants can be found in [Ishimaru 1978]. Equation 4.10 is inserted into the radiative transfer equation ( 4.5 on the preceding page). The equation is multiplied by ω ~ , and is again integrated over all directions. The result is given by [Ishimaru 1978] ~ ~ 1 (x). ∇φ(x) = −3σt0 (x)E + Q

(4.11)

4.1 Diffusion Approximation

45

Similar to equation 4.7, equation 4.11 describes how a change in fluence due to radiative transfer is given by the sum of attenuated radiance due to outscattering and absorption, and radiance that is added to that direction due to a volumetric light source, only now expressed in terms of diffuse flux. As before, the constants ~ 1 has a definition similar to Q0 follow from the derivation. Q

~ 1 (x) = Q

Z Qri (x, ω ~ )~ ω dw, 4π

only now the volumetric source due to incident reduced intensity is defined by a net radiative flux, since each direction ω is weighted by the corresponding radiance. Insertion of equation 4.11 on the preceding page into equation 4.7 on page 43 results in the steady state diffusion approximation (or diffusion equation) for a heterogeneous medium, which is given by ~ · (D(x)∇φ(x)) ~ ~ ·Q ~ 1 (x). ∇ = σa (x)φ(x) − Q0 (x) + 3D(x)∇

(4.12)

Here D is called the diffusion coefficient, which is given by

D=

1 , 3σt0

and follows from the derivation of the diffusion approximation. The diffusion equation states that the density of the gradient flow (a change in the flow of radiance) at point x is given by the source or sink at point x, i.e. radiance that leaves the system due to absorption and radiance that enters due to volumetric light sources. ~ 1 represents light sources due to As previously mentioned, the terms Q0 and Q the incident reduced intensity from the boundary of a semi-infinite medium. They can also be used to express respectively an isotropic and a anisotropic ~1 light source in an infinite emissive medium. If terminology is strict, Q0 and Q can be maintained as only reduced intensity sources, and a third source term can be added to express self emission from the volume (the terminology used by [Ishimaru 1978]). ~ 1 will be placed approximately one mf p0 away from the Intuitively, Q0 and Q surface of the semi-infinite volume, since this is the depth where the first scatter-

46

Diffusion

ing event takes place. The definitions of these reduced intensity sources proves to be particularly useful for the dipole and multipole approximations presented in section 4.2. Note that [Keijzer et al. 1988] have a slight error in their description of the derivation. In their description they incorrectly describe equation 4.7 as the result after insertion of equation 4.10 into, and the following integration of, equation 4.5. This may lead to some confusion about the derivation if references [Ishimaru 1978], [Keijzer et al. 1988] and [Jensen et al. 2001] are compared. However, it is only a question of an error in the description of when equation 4.10 should be inserted, it does not affect the final result, making the error easy to miss. For a homogeneous medium, equation 4.12 on the previous page is slightly simplified, due to constants that no longer depends on the position in the volume. It is given by ~ ·Q ~ 1 (x). D∆φ(x) = σa φ(x) − Q0 (x) + 3D∇

(4.13)

The previously mentioned diffusion coefficient depends on an appropriate choice of phase function that describes the scattering of light in a participating medium, as was mentioned in chapter 3 on page 11. If the derivation of the diffusion approximation is made based on the Henyey-Greenstein phase function given by equation 3.12 on page 23 (which is normally the case), the diffusion coefficient is formed through an expansion of the phase function in Legendre polynomials. For practical purposes, the full expansion is truncated in order to obtain the diffusion coefficient with a lower computational expense. The following expression for D is closer to the original expansion, and may better account for media that exhibit anisotropic scattering as well as higher absorption [Ripoll et al. 2005; Frisvad 2008]

1 D= 3σt0



4 σa 1− 5 σs (1 − g 2 ) + σa

−1

.

(4.14)

It was previously mentioned how it takes time for a volume to enter a steady state. The steady state is normally also what is needed for image synthesis, since this process in many ways resembles the one of taking a real photograph with a proper light configuration. The scene and its illumination is assumed to be perfectly still, which also means that the propagation of light in a material should have converged, i.e. there is no change in the concentration of radiance

4.1 Diffusion Approximation

47

(fluence) at different points in the volume. In the case where a visualization of the propagation of light is needed, it would be necessary to formulate the problem in time as well, in other words to let light enter and flow through the volume in small iterative steps while capturing an image of each step. This is not within the scope of this thesis, and estimating nth order multiple scattering therefore is a question of finding a solution to the diffusion equation. In some cases, such as infinite volumes where no volumetric light sources are ~ 1 in equation 4.12 on page 45 vanish. In this case, present, the terms Q0 and Q the equation simplifies to ~ · (D(x)∇φ(x)) ~ ∇ = σa (x)φ(x),

(4.15)

A further simplification is when the medium is non-emissive and also homogeneous, which means the coefficients become constant. The homogeneous equation is given by

D∆φ(x) = σa φ(x),

(4.16)

where ∆ is the divergence to the gradient (also referred to as the Laplacian). One thing that was not yet explicitly covered is how fluence is formulated with diffusion theory at the boundary (surface) of a participating medium. The previous definitions and equations only describe the steady state equations for fluence inside a (possibly infinite) medium. Since computer graphics generally is concerned with images of objects of finite size and shape, the boundary must also be taken into account, in order to find a solution for an entire volume.

4.1.1

Boundary Conditions

The diffusion equation 4.12 on page 45 only describes propagation of light modelled by diffusion due to nth order scattering, and not the radiative flux across the surface of a finite medium. This makes it necessary to model the propagation of light across the surface specifically, which is called a boundary condition. However, the boundary conditions can be formulated a bit differently with respect to the radiance that enters the volume. One approach is to formulate the conditions in terms of light sources that are ~ 1 sources from equation 4.12 placed at the surface of the volume. The Q0 and Q

48

Diffusion

are then neglected, unless the volume exhibit self emission of radiance. This is the approach employed by [Arbree et al. 2009], and is a normal approach for solving the diffusion approximation by finite difference or similar methods [Jacques and Pogue 2008]. A different approach is the one described in the previous section, where Q0 ~ 1 are used to describe radiance that has crossed the boundary into a and Q semi-infinite medium, until the point where it is scattered and thus becomes a reduced intensity light source. In this case, the radiance sources only reside some distance into the volume, while the boundary conditions are formulated entirely in terms of net inward radiative flux due to internal reflection. Only the second approach will be covered here, as it forms the basis for the methods employed by the solution presented in this thesis. The boundary condition is first presented for the case of a radiative flux across the boundary between two materials of equal refractive indices. This case has a simple intuition, the net inward flux must be zero at each point on the surface of the volume, since radiance that leaves the volume due to nth order multiple scattering is assumed not to return. Also, no internal reflection takes place due to the equal indices of refraction. This boundary condition is given by [Duderstadt and Hamilton 1976]

Z 2π(−)

Ld (xs , ω ~ )(~ ω · −~n(xs ))dω = 0,

(4.17)

Here the negative sign represents integration over 2π solid angles on the inward hemisphere, and the directions ω ~ are inward at surface point xs . Equation 4.17 can be written in terms of fluence given the spherical harmonics expansion of radiance, equation 4.10 on page 44. The boundary condition then becomes [Duderstadt and Hamilton 1976]

~ φ(xs ) − 2D(~n · ∇)φ(x s ) = 0.

(4.18)

Here φ(xs ) signifies the solution fluence at the boundary, which should be equal to the magnitude of outward radiative flux. This is a differential boundary condition, since a value is given for the solution to both a function and its normal derivative. This type of boundary condition is referred to as a Robin boundary condition. More details about boundary conditions, as well as solving differential equations numerically in general, can be found in [LeVeque 2007].

4.1 Diffusion Approximation

49

In the case of mismatching indices of refraction, internal refraction may occur. This means that the previous condition should be adjusted such that the net reflected inward radiative flux equals the radiative flux due to internal reflection. In order to write this condition, it is necessary to formulate the internally reflected inward net flux in terms of the diffuse reflectance, which is given by Z Fdr = 2π(+)

Fr (η, ~n · ω ~ 0 )(~n · ω ~ 0 )dω 0 .

(4.19)

Here η is the relative index refraction, the index of the materials light comes from over the material that light refracts into η = ηi /ηt . The definition of the Fresnel reflectance Fr was previously covered in section 3.2, and specifically given by equation 3.3 on page 18. The intuition of this equation should be clear, the integral over the Fresnel reflectance given all incident directions ω 0 , weighted by the cosine to the incident angle, yields the total fraction of energy that is reflected by the surface, hence diffuse reflectance, or average Fresnel reflectance. The diffuse Fresnel reflectance has a simple approximation, found by polynomial fit over the diffuse reflectance of all relative indices of refraction, given by [Egan and Hilgeman 1979]

Fdr ' −1.4399η −2 + 0.7099η −1 + 0.6681 + 0.0636η, ' −0.4399 + 0.7099η

−1

− 0.3319η

−2

+ 0.0636η

−3

,

η>1

(4.20)

η