Views
ReflectanceFunctions
From Odwiki
Under heavy construction. Enter only if wearing a hard hat! 
If you've read any articles or tutorials on shading, then chances are good that you've already come across the acronym BRDF, since it represents a component that just about all shaders must deal with sooner or later. It is one of those ubiquitous yet inscrutable little acronyms for which knowing what each initial stands for, does nothing toward explaining what it actually means.
The initials stand for Bidirectional Reflectance Distribution Function... There; I'm done.
"Huh?!?"
Contents |
[edit] What is a BRDF?
The basic concept is actually quite simple to grasp, but some of the math involved in representing it can get a little hairy at times. My goal here is to try to make it as accessible as possible -- including the math. Let's start by taking this acronym appart.
[edit] Function
The first thing to note is that we're talking about a function. That is to say, a mathematical expression that relates one thing to another -- a device that maps some input to a corresponding output. That's what the Function part of the acronym stands for. In this case, our function is supposed to tell us how much light bounces off a surface in the direction of our eyes (the output). And in order for it to be able to calculate this, it will need to be fed something -- what exactly that something is, will depend on how the function is put together, but at a bare minimum, it will need to know just how much light is hitting the surface to begin with, the direction said light is coming from, and the normal to the surface (the input).
[edit] Reflectance
The purpose of this function then, is to relate incoming to outgoing light intensities -- the amount of light that gets reflected off the surface (and in the viewer's direction). That's what the Reflectance part of the acronym refers to. Figure 1 shows this relationship (in two dimensions). Note that some of the arriving intensity may get absorbed or scattered inside the medium. But our little function doesn't care about that part; it just concerns itself with the amount that gets reflected in a given direction.
[edit] Bidirectional
At this point, we could extend our little experiment and try to see what would happen if we were to swap the viewing direction and the direction to the light. It can be shown that, if our object doesn't emit any light of its own, and if our function is physically correct, then the result should be exactly the same in both cases. In other words, our function should output the same result in both directions. It is this characteristic that earns it the title "Bidirectional" -- identical results in both directions, as shown in figure 2.
[edit] Distribution
But what if the viewer were to move to a different position? If we think of a glossy sunlit surface, it would seem intuitive to think that the amount of light reflected towards the new position could indeed be different than the amount reflected toward the old position, even assuming that everything else stays the same. In fact, it is probably safe to say that all natural surfaces reflect light in different amounts depending on the relationship between all three principal geometric components: the direction of the incoming light, the direction to the viewer, and the surface normal -- sometimes the difference is subtle and sometimes pronounced, but a difference there usually is. Figure 3 shows what our hypothetical function puts out for a number of different viewer directions, while keeping everything else fixed.
Instead, we could choose to represent intensity using the length of the line -- the longer the line, the more intense the reflection. This would make it cleaner and allow us to fit a lot more samples in our plot. Figure 4 plots our function using this method.
Finally, since we already know that all our lines start at the point where the incoming light hits the surface, we don't really need to draw the lines themselves, just the endpoints of each line. This would make it even cleaner and allow us to fit even more samples. As a final aesthetic touch, we can join all our endpoints together to get a better idea of the overall shape, giving us a curve. What we've just done (and what is shown in figure 5) is create a polar plot of our function.
Given an incoming light direction and a surface normal then, we can see by looking at figure 5 that what our function is actually doing is calculating a Distribution of reflected intensities over the semicircle above our surface -- it takes an incoming intensity, and redistributes it along different directions. And besides noting the obvious fact that this explains the last bit of the acronym, I should also point out that this distribution is only applicable to that particular combination of incoming direction and surface normal. If we were to change the relationship between those two components, and then plot another sweep of viewing directions, we can be pretty much guaranteed to get a different distribution... and now you should start to get a feel for this little monster. It should, at least in broad strokes, give you some idea of what the term Bidirectional Reflectance Distribution Function means.
[edit] The Reflection Vector
Having established that a BRDF deals strictly with calculating reflected light, this seems like a good time to delve a little deeper into what a perfect reflection actually is.
The direction of reflected light from a perfect mirror surface follows the very simple rule that the angle of incidence is equal to angle of reflection. Specifically, that the angle between the direction to the light source and the normal to the surface is equal to the angle between the normal to the surface and the direction of perfect reflection.
[edit] Deriving a Reflection Vector
Using the labels shown in figure 6, and assuming that all the vectors shown in the figure have a length of 1 unit (i.e: they are "normalized"), we can derive a formula to calculate the direction of perfect reflection ωr in terms of the known vectors ωi and N as follows:
- Using vector projection, we first calculate the vector component of ωi that is perpendicular to (and directed toward) N , which is given by perpN(ωi) = ( N · ωi ) N − ωi
- Referring back to figure 6, we note that the distance from ωi to the vector we're after, ωr , is exactly twice the lengh of the vector we just computed. Therefore, we can express ωr as
| ωr | = | ωi + 2 perpN(ωi) |
| = | ωi + 2 [ ( N · ωi ) N − ωi ] | |
| = | 2 ( N · ωi ) N − 2 ωi + ωi | |
| = | 2 ( N · ωi ) N − ωi |
[edit] Using VEX's reflect() Function
The good news is that VEX already has a built-in function to calculate a reflection vector, with the signature vector reflect( vector direction, normal). The only thing to keep in mind when using this function is that the first argument (direction) is expected to point toward N -- i.e: in the opposite direction of ωi . So to duplicate our calculation above using VEX's built-in function, we would call it like this:However, it is always good to understand the inner workings of a function such as reflect(), since working with BRDFs is fraught with calculations such as these, and there won't always be a VEX equivalent to hand. Expect to see more of this type of thing as we move on -- it's just the nature of the beast.
[edit] The Lambertian Effect
So far, we've been concentrating mostly on the outgoing side, but what about the incoming side? Exactly how much "light stuff" is striking the surface at the point we're interested in? And just how do we go about getting this information in VEX?
I'm sure it's become painfully clear by now that my terminology for referring to the "thing" that is hitting the surface has been all over the place: "light stuff", "light energy", "light intensity" , and so on. To remedy this, we'll borrow some terminology from the field of radiometry, which concerns itself with studying the transfer of energy via radiation, and incidentally, is the terminology used in all the literature.
[edit] Light, by any other name...
The amount of radiation emitted by a light source (or radiant power) per unit time , is called flux. Flux, is measured in watts, and given the symbol W . The amount of radiation (power) either emitted by a light source or received by a surface per unit area is called flux density and is measured in watts per square meter ( W / m2 ) . When flux density is emitted by a light source, it is called radiance, and when emitted by a surface, radiosity. Finally, when flux density is received by a surface, it is called irradiance.
Good; now we can stop calling it "light stuff" :-)
[edit] Lambert's Cosine Law
Even though the flux emitted by a light source is the same as the one arriving at the surface -- i.e: the amount of watts per hour put out by a light bulb doesn't change just because we are measuring it from some point on a surface --, the flux density can indeed be different between the two locations. To better understand why this is so, let's consider the cross section of a light beam arriving at the surface. Again we'll continue to think about this in 2D because it's easier to visualize, and because the relationship still holds in 3D, even though the terminology changes.
First, recall that flux density is related to area, so let's look at how the area of the emitter and the receiver relate to each other. If we look at figure 8, we notice that the emitter's cross sectional surface area A decreases as the angle of incidence θ increases. Clearly, this means that there must be some relationship between the emitter's area (from now on it is understood that when I say area, I mean the cross sectional 2D line) and the angle of incidence. In figure 8, this unknown factor is labeled as f(θ) -- meaning "some function of the angle of incidence θ ", but we don't know what that is yet.
To find out, we need to take a closer look at the geometry between the two areas. If we imagine our beam detaching itself from the "floor" as it "tips over", then we can interpret the angle of incidence as an angle of inclination away from the normal, and realize that it is equal to the angle of elevation from the "floor". Forming the righ-angle triangle at the bottom of the beam and using basic trigonometry, we come to the conclusion that the emitter's area is "whatever it was when aligned to the surface normal" ( A ), times the cosine of the angle of incidence ( A cos θ ). In fact this relationship is also valid when the emitter is directly perpendicular to the surface, because in that case θ = 0 , meaning that A cos θ = A cos 0 = A . So the relationship holds everywhere and our mystery function is cosine(). And since we're always assuming that the vectors N and ωi are normalized, then we can express the cosine of the incident angle θ as the dot product N · ωi .
Now that we understand how the two areas relate to each other, we can have a look at how this may affect radiance/irradiance. We will denote flux density with the symbol Φ (upper-case phi), so that the flux density emitted by the light will be ΦE , and the one incident to the surface ΦI . Now, the density emitted by the source doesn't depend on the angle of the source to the surface -- it is what it is regardless of angle. But as we discovered, the area of the emitter exposed to the surface is proportional to the cosine of the angle of incidence; and as flux density depends on area, this means that the emitted flux density is inversely proportional to it. So if we had a light source emitting P watts of radiant power toward our surface, we can express the flux density ΦE emitted by the light source (radiance) as
where A is the area of both the emitter and the surface at normal incidence (when the emitter is aligned with the surface normal) and A ( N · ωi ) is the cross sectional area of the emitter. And since the flux density ΦI incident on the surface (irradiance) is equal to P / A , we have the relation
A couple of things that should be mentioned before wrapping this up:
- The property just described is independent of the reflectance properties (BRDF) of a material -- it would apply equally to a polished metal surface as to a cement block. I relates the intensity of the light incident to a surface with the intensity leaving the emitter, and therefore has no relationship to reflectance whatsoever. It simply answers the question "How much light is arriving at the surface?".
- Something I haven't mentioned is that for the above analysis to hold true, we have been assuming that the beam (and the associated surface areas) is infinitesimally small (or that it is generated by a point source which is infinitely far away, such that all rays are parallel). In fact, this assumption will continue to be made as we move into 3D.
[edit] Lambert's Ideal Diffuse BRDF
We haven't properly defined a BRDF yet. To do that requires that we look at the full 3D geometry, and derive a formal mathematical representation for it in those terms. But even with our somewhat limited 2D model, we have enough to start implementing some of the most common BRDFs, so let's start with what we've grown accustomed to calling simply "diffuse", but which is also known as Lambertian diffuse.
[edit] Derivation
A Lambertian surface is defined as one which reflects light equally in all directions. There is no physical basis for this model -- i.e: there are no surfaces in nature which are perfectly diffuse (although the paint industry spends a great deal of effort trying to generate synthetic ones that approximate it). That's why this model is also referred to as an ideal diffuse model. However, there are some surfaces in nature that do come close to this ideal, e.g: chalk, talcum powder, cotton cloth, rough matte paper, etc.
Since the reflected intensity is defined as being the same in all directions, then it follows that it is independent of viewing angle -- i.e: a lambertian surface should look equally as bright from all viewing angles (within the hemisphere of illumination). To come up with a BRDF that would give us such a distribution, we have to look at how the viewer's area projects onto the surface, and realize in doing so that Lambert's cosine law comes into play again, except this time in reverse. As illustrated in figure 10, the amount of surface area seen by an observer increases in inverse proportion to the cosine of the outgoing angle θo . This implies that the perceived brightness would be greater at grazing angles than when viewed straight on. To cancel this out, and thereby achieve equal intensity from all viewing angles, we need to multiply by cos θo . Putting it all together, we have the following:
- Assuming equal reflectance in all directions, the reflected irradiance (radiosity) ΦR seen by the viewer in the direction ωo with angle θo , is inversely proportional to the cosine of that angle, and relates to incident radiance ΦE from direction θi as follows:
- The Lambertian diffuse BRDF ρ ( θo ) = cos θo then modulates this amount by cancelling the cosine term in the denominator, giving us
As we can see, all dependence on the viewing angle has disappeared, leaving us with raw irradiance. But even though this shows the basic relationship between incoming and outgoing radiance, it doesn't represent a full illumination model yet. One of the things we're missing is a way to describe the amount of diffuse reflectance that a given material may have, since not all materials reflect diffusely by the same amount. This is an extra scalar parameter to the BRDF which is usually referred to as the diffuse-reflection coefficient and given the symbol Kd. All it does is scale the final result. Incorporating this coefficient, we get Rd = ΦE Kd ( N · ωi ) . The Kd term can also be interpreted as the diffuse albedo which we'll see later, and there's also a factor of π lurking in the calculations, but we'll ignore those subtleties for now. In mathematical symbols then, this is a complete local illumination model to represent materials whose reflectance properties approximate those of an ideal diffuse reflector. Now let's see if we can turn this into a VEX function that we can use in our shaders.
[edit] Implementation
So far, we've been using "BRDF terminology" -- the naming conventions and choice of symbols have been the same ones used in most of the literature. This notation however, doesn't translate well to programming languages such as VEX or RSL. As a result, a separate set of symbols have, with time, come into common use among shader writers. There is no clear concensus among shader writers though, so the variable names given here are only meant to be a representative sampling of some of the common choices.
[edit] Omega, by any other name...
The following table shows a mapping between the symbols we've seen so far, and the names typically given to them in VEX and RSL. Due to some differences in available data types between the two languages, the variable names are given in the context of a full definition statement. Wherever an actual definition would be uncommon, or where the choice of name for the variable is completely arbitrary, the actual variable name is dropped and only the expression that evaluates to the appropriate quantity is shown, along with the resulting data type.
| Symbol | VEX | RSL |
|---|---|---|
| N | Unit surface normal. | |
| vector Nn = normalize(N); Additionally, as an option:: | normal Nn = normalize(N); Additionally, as an option:: | |
| ωi | Unit direction to the light source. | |
| vector Ln = normalize(L); ¹ | vector Ln = normalize(L); ¹ | |
| ωo | Unit direction to the viewer. | |
| vector V = normalize(−I); | vector V = normalize(−I); | |
| ωr | Unit reflection vector, i.e: mirror reflection of ωi about N. | |
| vector R = reflect(−Ln,Nn); ¹ | vector R = reflect(−Ln,Nn); ¹ | |
| cos θi | Cosine of the angle between ωi and N. | |
| [float] dot(Nn,V) | [float] Nn.V | |
| cos θo | Cosine of the angle between ωo and N. | |
| [float] dot(Nn,Ln) | [float] Nn.Ln | |
| ΦE | Flux density emitted by the light source. | |
| [vector, global] Cl ¹ | [color , global] Cl ¹ | |
|
¹ The global variables L and Cl are only available within an illuminance() construct. | ||
[edit] Seeing the Light
The incoming radiance that a BRDF is put in charge of distributing, and which could potentially contribute to reflectance toward the viewer, could come from any direction within the hemisphere about N. Therefore, we need some way to account for all sources of irradiance within that set of directions. Once we have the means of acquiring both the direction and intensity of each (infinitesimally thin) incoming ray, we repeat, for each one, a very simple procedure: pass this information along with the viewer's direction to a BRDF, which tells us what percentage of that irradiance gets reflected toward the viewer; add this single sample of weighted irradiance to an accumulator, and repeat. The final accumulated value is then the reflected radiance (radiosity) bound in the viewer's direction. In mathematical terms, we're integrating BRDF-weighted irradiance over the hemisphere; not as a continuous function, but through a sum of discreet samples.
The problem is that, since each ray is "infinitesimally thin", then there are an infinite number of rays that could fit in the hemisphere. And if we had to apply our procedure to each one, we'd be rendering for a very long time indeed. Luckily, there are many methods available that will cut this process down to acceptable time spans. In our case, things couldn't be simpler: the "stepping through each ray" is done for us; only along the directions that "matter the most", and as efficiently as possible (within a programmable shader architecture). Nice.
Both languages make the same mechanism available to us for computing irradiance: the illuminance() construct. We can think of it as a parameterized code block -- it has the syntactic feel of a function, but is really a block of code that gets run for every light source that is of importance to us. Just what is "of importance to us" is what we define using the parameters, the rest is the snippet of code that will run for each iteration. For each light source within our volume of interest, the renderer will give us its direction L and intensity (radiance) Cl. The two vector-valued variables L and Cl are global, and only accessible inside illuminance code blocks.
illuminance( /* parameters that define our hemisphere of interest */ )
{
// Code block that runs for each light source in that hemisphere.
// For the light source currently being visited, we are given its
// direction (L) and radiance (Cl). We use these two quantities
// inside this block to calculate our BRDF and accumulate the
// result into a running sum, which is the final reflectance.
}
[edit] Show Me the Code
Using the reference table to map symbols to expressions, and given the generous amount of scaffolding provided by both VEX and RSL, there isn't really much left to do except to map the mathematical form of the BRDF N · ωi directly, term by term, and then wrap it up into a function. In VEX, we get
float brdf_diffuse(vector l,n) {
return max(dot(l,n),0);
}
float brdf_diffuse_fast(vector l,n) {
return dot(l,n);
}
And for RSL:
float brdf_diffuse(vector l; normal n) {
return max(l.n,0);
}
float brdf_diffuse_fast(vector l; normal n) {
return l.n;
}
There are two flavours of the BRDF, the "_fast" version doesn't perform the check to ensure that ωi and N are in the same hemisphere -- it can only be called safely from any context where you categorically know that this condition is satisfied, and possibly save a few cycles. Other than that, there isn't much going on here beyond the dot product, which is the totality of this particular BRDF. As I said, it's a simple one.
Now we can build the local illumination model for this BRDF. The only important part that's brough into play here, is the ΦE term which, as mentioned, is acquired through the use of the illuminance() construct. For VEX, it can be implemented as follows:
#include <math.h> // for the definition of M_PI_2
#include <shading.h> // for the definition of LIGHT_DIFFUSE
vector illum_diffuse(vector p, n; string mask) {
vector Rd = 0;
illuminance(p,n,M_PI_2,LIGHT_DIFFUSE,"lightmask",mask) {
shadow(Cl);
Rd += Cl * brdf_diffuse_fast(normalize(L),n);
}
return Rd;
}
And for RSL:
color illum_diffuse(point p; normal n; string category) {
color Rd = 0;
illuminance( (category==""?"-_":category), p, n, PI/2 ) {
Rd += Cl * brdf_diffuse_fast(normalize(L),n);
}
return Rd;
}
The naming convention used here is to prefix the BRDF functions with "brdf_", and the local illumination functions with "illum_". In both cases, the accumulator is Rd, and the incoming flux density is given by the global variable Cl.
The only thing left to do is to write a shader to test these functions. As it turns out, both RSL and VEX have built-in functions to calculate the BRDF and local illumination model for Lambertian diffuse. We can take advantage of this and use them to compare their results with our own, so we'll include this in our test shader.
For VEX:
#pragma label brdf "Diffuse BRDF"
#pragma choice brdf "0" "Custom"
#pragma choice brdf "1" "VEX"
#pragma label Cdiff "Diffuse Color"
#pragma hint Cdiff color
#pragma label Kd "Diffuse Amplitude"
#pragma label lmask "Light Mask"
surface VEX_Diffuse(
int brdf = 0;
vector Cdiff = 1;
float Kd = 1;
string lmask = "";
)
{
vector Nn = normalize(frontface(N,I));
Cf = Cdiff * Kd *
(brdf ? diffuse(Nn,lmask) : illum_diffuse(P,Nn,lmask));
}
And for RSL:
surface RSL_Diffuse(
float brdf = 0;
color Cdiff = 1;
float Kd = 1;
string lmask = "";
)
{
normal Nn = normalize(faceforward(normalize(N),I));
Ci = Cdiff * Kd *
(brdf==1 ? diffuse(Nn) : illum_diffuse(P,Nn,lmask));
}
Note: I will make various otls available to reproduce these results when I get around to packaging them up.
[edit] An Analysis of Our First BRDF
Now that we've managed to put together our first BRDF, it's time to talk about some items that have suddenly appeared in the implementation but haven't been covered so far.
[edit] But I thought you said...
What is that N · ωi term doing inside the BRDF definition?!
It certainly has no business being there since it's part of the calculation for irradiance, and as was mentioned earlier, is completely independent of reflectance -- i.e: it applies to all BRDFs and is not somehow magically attached to Lambertian diffuse. The only reason for it being there instead of in its rightful location, which is in the illumination model, has to do with how these models are typically built for shaders. Inside an illuminance loop, it is customary to apply the BRDF as a direct weighting factor for Cl:
illuminance(...) {
accumulator += Cl * BRDF;
}
The main advantage of this is that the shader writer doesn't need to keep the N · ωi factor in mind (or even know or be aware of it), which would frankly seem quite foreign unless s/he knew how BRDFs were put together -- and even then, it would be inconvenient to keep having to write the same quantity over and over. The reality is that it doesn't matter where the factor is applied, as long as it is. Additionally, it is not uncommon for some BRDFs to cancel out this term for other reasons, which would mean wasted computation if it were left outside the BRDF function itself. For practical reasons then, all implementations from now on will include this term in the BRDF function itself. However, just to state it clearly, the mathematical representation of reflectance (at least as we know it thus far) is this:
Φo ( ωo ) = ΦE ( ωi ) ρ ( ωi , ωo ) ( N · ωi )
The BRDF in that expression is ρ ( ωi , ωo ) and the term ( N · ωi ) is clearly outside of it. If we were reading it out loud, we would say "The outgoing reflected radiance Φo in the viewer's direction ωo is equal to the radiance ΦE emitted by the light source from the incoming direction ωi times the reflectance distribution function ρ ( ωi , ωo ) times the irradiance factor ( N · ωi ) ". Incidentally, the symbol ρ is the Greek letter "rho", which is equivalent to our 'r', for "reflectance" -- I didn't pick it, it is just the symbol frequently used to denote a BRDF.
[edit] Then what is the real Lambertian BRDF?
Lambertian surfaces have the property that the amount of light reflected from a unit differential area toward the viewer is directly proportional to the cosine of the angle between the viewer and the surface normal. Therefore, the diffuse BRDF is cos θ0 = ( N · ωo ) (please note that this is not the same as the irradiance factor ( N · ωi ) ). But, as mentioned earlier, the amount of surface area seen by the viewer is inversely proportional to cos θ0 and so the two factors end up canceling each other out, leaving just the constant Kd. This canceling out of the two factors is not an accident; it is designed that way so that the perceived brightness is the same everywhere (constant) -- see figure 13.
To convince ourselves of this fact, we can illuminate a flat plane with a directional light (parallel rays) and measure the reflected intensity as the viewer and/or the plane rotates with respect to the light. This is illustrated in figure 14. Notice that 1) the reflectance in both rows is even across the entire patch, even though the camera is not orthographic. 2) as the plane rotates away from the light (top row) it reflects less and less light as there is less irradiance to begin with (due to the cosine law), and 3) as the viewer rotates away from the surface normal, and the light and surface remain fixed (bottom row), reflectance remains constant since the relationship between the light and the surface does not change.
[edit] And what about color?
So far in our discussions we haven't talked about color. BRDFs are meant to be wavelength dependent, so that a more refined definition of our reflectance formula would look like this:
Φo ( ωo , λ ) = ΦE ( ωi , λ ) ρ ( ωi , ωo , λ ) ( N · ωi )
where λ (lambda) represents one particular wavelength. But it is not usually notated this way because the dependency is implicitly understood to exist, so the extra parameter is dropped for the sake of clarity. However, what this implies is that there could potentially be different distributions for different wavelengths (for the same material); and this is indeed the case for most natural materials. If you're curious, you can look at plots for several wavelengths of some measured BRDFs from real surfaces that have been published, and you'll see slightly different distributions between the sampled wavelengths. For the majority of materials however, this difference is not very pronounced and we can get away with using a single distribution across the spectrum. For some materials, the difference is very pronounced and this spectral variation needs to be taken into account -- the classical example of such a material is the underside of a Compact Disc, whose reflectance has a prismatic effect.
For the moment, we will stick to the way color is commonly treated in shaders. The BRDF is constant across the spectrum, and is simply modulated by a color coefficient that is meant to represent the chromatic response of the material for that component in the illumination model (the parameter "Cdiff" in the code), and finally weighted by a scalar "dimmer" (the parameter "Kd" in the code) to control the overall contribution of the component. A generic representation of this treatment then looks like this:
| Cfinal = K1*C1*ILLUM1 + K2*C2*ILLUM2 + ... + Kn*Cn*ILLUMn; | |
|
where each ILLUMn term expands to: | |
| ILLUMn = 0; illuminance() { | |
[edit] Specular Reflectance
In addition to a uniform diffuse reflectance, opaque materials with smooth surfaces tend to show strong reflectivity along the reflection vector. This specularity, which shows up as a bright spot when the source is a point light, can be thought of as a blurry mirror reflection of the incoming light. Because it is closely related to a perfect mirror reflection, it makes use (either directly or indirectly) of the reflection vector, which, for specular BRDFs, becomes an important axis of the distribution.
[edit] The Phong Model
One of the early models for this type of reflectance, which is still in use today, is the Phong model. It is not physically plausible in that it is not energy conserving, and is also not reciprocal -- i.e: not bidirectional. But it is reasonably convincing, fast to compute and, more important to us at this stage of our exploration, easy to implement.
The model simply dampens reflectance according to the angular distance (cosine of the angle, actually) between the viewing direction ωo and the reflection vector ωr . To add some control over the shape of the falloff, the angular distance is raised to some power, determined by the user through the parameter "shinyness". In the form of a BRDF, and using the symbol m to denote shinyness, we can write it as follows:
And the lighting model as
Again, when implementing it as a BRDF for shading, we will include the last term ( N · ωi ) in the BRDF function itself, even though it doesn't really belong there. However, if you look closely, you'll notice that the BRDF ρph(ωi,N,ωo) also has this term in the denominator, and so they will cancel out leaving just (ωr · ωo)m , where ωr = reflect(−ωi,N) is the reflection vector, and the exponent m is the user-supplied "shinyness" parameter.
Again, we are lucky that both VEX and RSL have a Phong lighting model built in, so we can compare. This time, I won't split the listings for each function as it should hopefully be clear by now which functions play which roles. Each listing includes a tiny test shader at the end.
#include <math.h>
#include <shading.h>
//------------------------------------------------------------
// BRDF and local illumination functions
// All vectors are expected to be normalized!
//------------------------------------------------------------
float brdf_phong(vector l,n,v; float shine) {
return pow(max(dot(v,reflect(-l,n)),0),shine);
}
vector illum_phong(vector p, n, v; float shine; string mask) {
vector Rph = 0;
illuminance(p,n,M_PI_2,LIGHT_SPECULAR,"lightmask",mask) {
shadow(Cl);
Rph += Cl * brdf_phong(normalize(L),n,v,shine);
}
return Rph;
}
//------------------------------------------------------------
// Test Shader
//------------------------------------------------------------
#pragma label brdf "Phong BRDF"
#pragma choice brdf "0" "Custom"
#pragma choice brdf "1" "VEX"
#pragma label Cph "Phong Color"
#pragma hint Cph color
#pragma label Kph "Phong Amplitude"
#pragma label shine "Phong Shinyness"
#pragma range shine 1 100
#pragma label lmask "Light Mask"
surface VEX_Phong (
int brdf = 0;
vector Cph = 1;
float Kph = 1;
float shine = 20;
string lmask = "";
)
{
vector Nn = normalize(frontface(N,I));
vector V = normalize(-I);
Cf = Cph * Kph *
(
brdf ?
phong(Nn,V,shine,lmask) :
illum_phong(P,Nn,V,shine,lmask)
);
}
//------------------------------------------------------------
// BRDF and local illumination functions
// All vectors are expected to be normalized!
//------------------------------------------------------------
float brdf_phong(vector l; normal n; vector v; float shine) {
return pow(max(v.reflect(-l,n),0),shine);
}
color illum_phong(point p; normal n; vector v;
float shine; string category)
{
color Rph = 0;
illuminance( (category==""?"-_":category), p, n, PI/2 ) {
Rph += Cl * brdf_phong(normalize(L),n,v,shine);
}
return Rph;
}
//------------------------------------------------------------
// Test Shader
//------------------------------------------------------------
surface RSL_Phong(
float brdf = 0;
color Cph = 1;
float Kph = 1;
float shine = 20;
string lmask = "";
)
{
normal Nn = normalize(faceforward(normalize(N),I));
vector V = normalize(-I);
Ci = Cph * Kph *
(
brdf==1 ?
phong(Nn,V,shine) :
illum_phong(P,Nn,V,shine,lmask)
);
}
[edit] The Halfway Vector
Originally introduced by James Blinn as an alternative to the reflection vector because it is cheaper to compute, the halfway vector is often preferred to the reflection vector when calculating specular reflections. Usually denoted with the letter H , the halfway vector is the direction lying exactly half way between the direction to the viewer ωo and the direction to the light source ωi , and is given by H = normalize( ωi + ωo ) .
The angle φ between ωr and ωo (the angular distance used in the Phong model) is exactly twice the angle α between H and N. This factor of 2 can be explained by looking at figure 16 and noticing that
| φ + θ | = | α + β |
| φ | = | α + β − θ |
But β = α + θ , so
| φ | = | α + α + θ − θ |
| φ | = | 2 α |
Given this constant relationship and the fact that the Phong model is not based on any physical properties, there's really no reason why one couldn't use H instead of the more expensive ωr when calculating the specular component. Doing so would produce different results in terms of the rate at which the specular highlight decays, but would retain the general feel of the original Phong model.
Interestingly enough, some recent research that attempts to validate some of the popular analytical models (which are the type being discussed in these notes) against measured BRDFs from real materials, seems to indicate that using the halfway vector gives a better match to the sampled data than using the actual reflection vector. The first analytical model to make use of this vector was Blinn's modification of the original Phong model, known as the Blinn-Phong specular model.
[edit] The Blinn-Phong Model
The only difference between the Blinn-Phong model and the original Phong model, is the use of the halfway vector H in place of the mirror reflection vector ωr . Visually, they are almost identical and so are their implementations.
However, we will take this opportunity to also change the original "shinyness" parametter to the more familiar "roughness" rendition (although this is not part of the original model). As the name implies, "roughness" is simply the inverse of "shinyness", so where we had the axponent m before, we will now have 1 / m instead. The advantage with this formulation is that the parameter values are a little more intuitive and remain in the (0,1] interval for the most part -- although the code should enforce the condition m > 0 to safeguard against division by zero and root of a negative number. Additionally, we can somewhat compensate for the fact that α is half the size of φ (see figure 16), by changing 1 / m to 2.83 / m , where 2 √2 = ~ 2.83 is an arbitrary scaling factor to get the size of the specular dot to match Phong a little better (although this is not strictly necessary).
#include <math.h>
#include <shading.h>
//------------------------------------------------------------
// BRDF and local illumination functions
// All vectors are expected to be normalized!
//------------------------------------------------------------
float brdf_blph(vector l,n,v; float rough) {
return pow(max(dot(n,normalize(l+v)),0.),2.83/max(rough,1e-4));
}
vector illum_blph(vector p, n, v; float rough; string mask) {
vector Rph = 0;
illuminance(p,n,M_PI_2,LIGHT_SPECULAR,"lightmask",mask) {
shadow(Cl);
Rph += Cl * brdf_blph(normalize(L),n,v,rough);
}
return Rph;
}
//------------------------------------------------------------
// Test Shader
//------------------------------------------------------------
#pragma label brdf "BRDF"
#pragma choice brdf "0" "Custom"
#pragma choice brdf "1" "VEX"
#pragma label Cph "Color"
#pragma hint Cph color
#pragma label Kph "Amplitude"
#pragma label rough "Roughness"
#pragma range rough 1e-5 1
#pragma label lmask "Light Mask"
surface VEX_BlinnPhong (
int brdf = 0;
vector Cph = 1;
float Kph = 1;
float rough = 0.05;
string lmask = "";
)
{
vector Nn = normalize(frontface(N,I));
vector V = normalize(-I);
Cf = Cph * Kph * (
brdf ?
phong(Nn,V,1./rough,lmask) :
illum_blph(P,Nn,V,rough,lmask)
);
}
//------------------------------------------------------------
// BRDF and local illumination functions
// All vectors are expected to be normalized!
//------------------------------------------------------------
float brdf_blph(vector l; normal n; vector v; float rough) {
return pow(max(normalize(l+v).n,0),2.83/max(rough,1e-4));
}
color illum_blph(point p; normal n; vector v;
float rough; string category)
{
color Rph = 0;
illuminance( (category==""?"-_":category), p, n, PI/2 ) {
Rph += Cl * brdf_blph(normalize(L),n,v,rough);
}
return Rph;
}
//------------------------------------------------------------
// Test Shader
//------------------------------------------------------------
surface RSL_BlinnPhong(
color Cph = 1;
float Kph = 1;
float rough = 0.05;
string lmask = "";
)
{
normal Nn = normalize(faceforward(normalize(N),I));
vector V = normalize(-I);
Ci = Cph * Kph * illum_blph(P,Nn,V,rough,lmask);
}
[edit] Enter The Third Dimension
So far, we've been analyzing the BRDF in two dimensions. Specifically, we've been visualizing the geometry as though it were only occurring in the incidence plane -- the plane described by the incoming vector ωi and the surface normal N. But the observer could be positioned along any arbitrary direction in the hemisphere above the surface (and, by extension, so could any vector derived from it, such as H). What we need are some new tools that will allow us to express our mental model in full three dimensions, and hopefully glean some new insights into the BRDF in the process.
[edit] Solid Angle
We're all fairly comfortable with the concept of angular measure in the plane, but we're likely more used to describing it in terms of degrees than radians. To better understand the 3D equivalent of a 2D angle, we need to remind ourselves what the radian unit of angular measurement represents.
As illustrated in figure 17, the measure of a planar angle θ in radians is given by the arc length l swept out on a circle, divided by the radius r of that circle, giving us the formula θ = l / r . In other words, it is the ratio of the arc length with respect to the radius, which means that a single radian unit represents an arc length that is equal to the radius. Since the longest arc length of a circle is its circumference, there are θ = 2 π r / r = 2 π radians in the central angle representing the entire circle.
As an aside, notice that in the degree system of angular measure, a full circumference is represented by 360 degrees. So we can describe the relationship between radians and degrees as the ratio 2 π / 360 = π / 180 . This means that to convert degrees to radians we can multiply the degree measure by π / 180 , and to go from radians to degrees we'd multiply the radian measure by the inverse of that ratio: 180 / π .
If we add another dimension to the concept of "length", we end up with the concept of "area". Similarly, if we add another dimension to the circle, we end up with a sphere. So if we extend the definition of a planar angle as described above to three dimensions, and call it a "solid angle" to distinguish between the two, we can say that the measure of a solid angle ω corresponding to an area A on the surface of a sphere with radius r is given by ω = A / r2 . The unit of measure for a solid angle is called a steradian, and is abbreviated sr. By direct analogy to the circumference of a circle in the plane, the total surface area of a sphere is equal to 4 π r2 , which means there are A / r2 = 4 π r2 / r2 = 4 π steradians in the solid angle representing the entire sphere.
[edit] Spherical Coordinates
Recall how in the introduction we chose to represent the distribution of a BRDF via a series of vectors of different magnitudes and called the result a "polar plot". Each vector in such a plot can be described by an angle and a radius (magnitude), and we have more or less consistently been referring to the angular component of such a pair with the symbol θ . In fact, given that we could describe any point on a plane using this combination of angle and magnitude, it follows that what we actually have here is a full (two-dimensional) coordinate system in its own right -- one that is particularly useful for describing elements of systems such as BRDFs. The angle θ that we've mentioned so much is actually a polar angle (a measure of rotation from a given pole, or polar direction), and in all our BRDF diagrams so far, the polar direction has always been the direction of the normal to the surface N .
For the purposes of converting between polar coordinates { r, θ } and 2D cartesian coordinates { x, y } , we can use the following relations:
| x | = | r cos θ |
| y | = | r sin θ |
| r | = | sqrt( x2 + y2 ) |
| θ | = | atan( y / x ) |
The three dimensional analog of polar coordinates in the plane is known as spherical coordinates, and the extension into the third dimension is fairly straight forward -- since we're adding another dimension, we need an extra angular measure in that dimension. As illustrated in figure 18, this new angle φ is a measure of rotation along the plane tangential to the surface and away from the unit vector -- or axis -- lying on the "floor" of the dome (labeled "u") which is orthogonal to the polar angle which measures angular distance relative to the direction N . The radius r of the spherical coordinate system is called the "radial" component, but the two angular components θ and φ are unfortunately known by several different names. The angle θ is variously referred to as "polar angle", "colatitude", or "zenith", and the angle φ as "azimuthal angle" (or simply "azimuth"), or "longitude". To confuse things even further, the convension where the symbols θ and φ are reversed is also frequently used. In these notes however, we will consistently use the symbol θ to refer to the polar angle, and φ to refer to the azimuthal angle.
For the purposes of converting between spherical coordinates { r, θ, φ } and 3D cartesian coordinates { x, y , z } with the familiar orientation used in Houdini, we can use the following relations:
| x | = | r sin θ sin φ |
| y | = | r cos θ |
| z | = | r sin θ cos φ |
| r | = | sqrt( x2 + y2 + z2 ) |
| θ | = | acos( y / r ) |
| φ | = | atan( z / x ) |
[edit] Differential Solid Angle and Area
Here, the term "differential" is used to denote the difference between two samples of the same function where the samples are an infinitesimal distance apart in the function's domain -- or simply, the change in one variable with respect to a very small change in another, dependent variable.
A differential solid angle dω can be thought of as a small solid angular "spread" around a given direction, and can be written using spherical coordinates in terms of the differential azimuthal angle dφ and the differential polar angle dθ . Recall that a solid angle depends on the area it subtends on a sphere (in the same way that a planar angle depends on the arc length it subtends on a circle) and is given by A / r2 which, in the case of a unit sphere (with r=1) reduces to A . With this in mind, we're going to define a differential solid angle dω by first defining the differential area dA it subtends.
Looking at figure 19, we see that the circle at the polar angle θ that lies on a plane parallel to the x-z plane and passes through the point p with spherical coordinates { r, θ, φ } has a radius equal to r sin θ . Therefore we can say that the differential arc length dlφ in the azimuthal direction is equal to r sin θ dφ . We arrive at this result by using the definition of a planar angle θ = l / r , which in this case is dφ = dlφ / r sin θ , and which we then solve for dlφ to get dlφ = r sin θ dφ . We then use the same method to arrive at the differential arc length dlθ in the polar direction, which results in dlθ = r dθ . Multiplying these together gives us the differential surface area
dA = r2 sin θ dφ dθ
Referring back to the definition of a solid angle ω = A / r2 and substituting the differential area dA that we just calculated, gives us dω = r2 sin θ dφ dθ / r2 , so that the differential solid angle dω corresponding to the differential area dA is
dω = sin θ dφ dθ
[edit] Redefining the BRDF
Earlier, we discussed irradiance (incident flux density) in terms of the energy arriving at a point on the surface from a single direction in the hemisphere about N. To calculate the total amount of irradiance incident from every direction in the hemisphere, we then looked at the special case of sampling several point-light sources using the illuminance() construct. Because all these point sources (including point, cone, and directional lights) can only ever illuminate the surface from one direction in the hemisphere, it allowed us to describe the total amount of irradiance striking the surface simply as a weighted sum of these sources (weighted by the irradiance factor).
From the point of view of writing shaders, we could simply stop there and call it a day since, even if we were dealing with complex area sources, both VEX and RSL provide us with functions that will calculate irradiance from them (as well as indirect sources, such as light reflected from other objects). But for the purposes of gaining a better understanding of BRDFs and in an effort to acquire the tools needed to convert published material to working code, we need to look at calculating irradiance as an integral over the hemisphere.
[edit] Irradiance as an Integral
Here we're going to modify our symbols slightly, so that we will use LI and LE to indicate irradiance and radiance in a single direction ωi , and ΦI and ΦE to refer to those quantities over the whole hemisphere.
The total irradiance (incident flux density) ΦI received by a differential area dA is equal to the integral of the irradiance LI(ωi) received from every direction ωi over the unit hemisphere Ω above the surface:
So we have a double integral. First we integrate Li(θ, φ) sin θ over 2 π radians (360°) in the azimuthal direction (φ), and then we integrate this result over π/2 radians (90°) in the polar direction (θ) completing the hemisphere. In order to do this, we had to change our vector representation of the incoming direction ωi to its polar version { θ, φ } so as to match the variables of integration. But as we found before when looking at this in 2D, the radiance received by a surface relates to the radiance emitted by the light source in direct proportion to the cosine of the angle of incidence, the same applies here. So the radiance LI received by the differential surface area dA and the radiance LE emitted by a light source are related by LI = LE cos θ , which means that we can rewrite the integral for irradiance as
[edit] Bidirectional Reflectance
Good; now we have a formal definition of irradiance. But what about reflectance? Reflectance is just the portion of irradiance (as we just defined it) that bounces off the surface. In other words, it is the ratio of the differential reflected radiance dLR to the differential incident irradiance dLI . And the function ρ(ωi , ωo) that calculates this ratio is the BRDF itself, which is a function of the direction to the light ωi and the direction to the viewer ωo . In symbols, we have
Which means that the differential reflected radiance dLR (the portion that bouces off in a single direction) is
And finally, we can write the total bidirectional reflectance ΦR for the entire hemisphere as
where the two directions ωi and ωo should really be written in their spherical forms as { θi , φi } and { θo , φo } , but which were left as vectors to make it easier to parse.
Notice that the only difference between ΦI and ΦR is the precense of the BRDF. Intuitively, this makes sense since reflectance is just a re-distribution of irradiance; and this is distribution is completely defined by the BRDF -- reflections don't generate light, they just distribute it.
[edit] Physical Plausibility
For a BRDF to be physically plausible, it has to adhere to the following physical rules (all of which should make sense at an intuitive level):
- Physical BRDFs must be non-negative.
In other words. reflected light should not be negative. - Physical BRDFs must be reciprocal (bidirectional).
As described earlier, this means that the BRDF should return the same result even if we exchanged the direction to the viewer ωi and the direction to the light source ωo . - Physical BRDFs must conserve energy.
The amount of reflected energy should not be greater than the amount of energy received by the surface.
This is the topic of our next section.
Having said that, there's nothing to stop us from using a non physically plausible BRDF if we think it is aesthetically correct. Nevertheless, it is useful to know what properties render a BRDF as "physically plausible".
[edit] Normalization and The Enigmatic π
Let's say we have some amount of incident radiant power LE coming in from direction ωi and we wanted to distribute it evenly in all directions of the hemisphere (as a Lambertian surface would). This means our distribution function should be a constant (same amount in all directions) KR , but what should that constant be? Given that the radiant power we're trying to distribute is coming from a single direction, we have an incident radiance that is also constant, and so we can pull both terms out of the integral.
This tells us that if we received 1 unit of emitted radiant power ( LE = 1 ) , and then reflected (or radiated) 1 unit of radiant power for each direction in the hemisphere ( KR = 1 ) , we would end up radiating a total (over the whole hemisphere) of π units -- i.e: we would be radiating π units more than we received, and as a result, our BRDF would not conserve energy. For our BRDF to conserve energy (and be "normalized") we should instead reflect a maximum of KR / π { 0≤ KR ≥1 } units in each direction. And this is the reason behind that mysterious 1 / π factor that keeps showing up in the published material -- it is a normalizing factor.
So if this is so important, then why haven't we included it in our BRDF implementations so far?
The answer to that question has more to do with "usage" than with "correctness". As users, we instinctively expect that if we illuminate a diffuse surface whose surface colour is set to 1 with a point light that has an intensity of 1, we would see a reflected intensity of 1 (at its brightest spot), and not 0.31830... (1/π). To achieve this intuitive result, we can pretend that our illuminance() summation is silently being multiplied by π and that our surface colour parameter (or "albedo") is implicitely divided by π, canceling out the π term. Or alternatively, that our units for light intensity are really multiples of π.
However; even though these assumptions have the desirable effect of simplifying usage, they are never taken into account in published papers on BRDFs, so that we shouldn't be surprised to see those BRDFs include a 1 / π term.
Note that a similar process can be applied to derive a normalizing factor for other analytical BRDFs. If we took the Phong BRDF for example, which is given by (cos θ)m and assumed that m > 0 , we would get
So that the normalizing factor for the original Phong model would be 2+m / 2π -- the inverse of the above result -- and the revised model would then be (cos θ)m (2+m) / 2π
. In fact, this particular formulation is known as the "Modified Phong" specular BRDF which, unlike the original, is energy conserving.
[edit] The Fresnel Effect
So far, we've been assuming that all the incident radiance is available for us to distribute using some BRDF. However, the reality of what goes on when light waves hit a surface is conciderably more complex than that. Part of the incoming light is reflected yes, but another portion may get refracted into the surface. Which portions get reflected and which refracted depends on the angle of incidence and the surface's refractive properties. The Fresnel equations are mathematical models that can be used to account for these two principal effects. In their original form, the equations deal with light polarization and complex indices of refraction, but since both VEX and RSL provide functions to calculate the two Fresnel terms, we won't go into any detail here except to present, without derivation, the formula for the Fresnel reflection factor for unpolarized light.
This equation is somewhat expensive to compute, so it is customary to use Schlick's approximation to to it, which is given by
We can make a couple of useful observations about the behaviour of this function. First, as the angle of incidence approaches 90° (light is at grazing angle), the value of ω · H approaches 0, and thus the value of Fr approaches 1. This means that at grazing angles all light is reflected (and none refracted). Second, for normal incidence (where the direction of incident light is aligned with the surace normal), the value of ω · H is 1, and Fr becomes ((η-1) / (η+1))2 -- this is shown as the symbol Rs in the equation for Schlick's approximation above.
This is useful in that we can use it to derive an approximate value for η even if all we know about a surface is the specular color it reflects at normal incidence. All we have to do is solve for η (in the equation for Rs , which yields
We can use this to calculate the relative indices of refraction ηλ for each channel { R, G, B } by only requiring that the user provide the specular color at normal incidence (i.e: at its darkest value), and avoid a request for an index of refraction, which is a fairly unintuitive quantity to have to deal with from the user's point of view. Nevertheless, requesting "the darkest value of the specular color" is also unintuitive, so we'll need to massge things somewhat -- more on this when we get to the implementation side of things; for now, it is good enough to know that this useful relationship exists.
[edit] Complex Fresnel
The problem with implementing the full unpolarized Fresnel reflection function directly, as given above, is that it can't be used with complex indices of refraction (it would produce non real-valued results in its natural angular domain [−π/2,+π/2]). In the form given above then, η should always be greater than 1 since any number between 0 and 1 would cause the expression η2−1+c2 inside the radical for the definition of g to dip into negative values for some portion of the range, depending on the index of refraction. We can still use it if we impose the constraint max(0, η2−1+c2) , but this would restrict the function to being mostly suited to representing non-conductive materials ("dielectric" materials) such as glass or water. If we also want to properly represent conductive materials such as metals, we need to use a complex index of refraction.
Complex numbers and their properties are beyond the scope of these notes, but for now, we can simply think of them as two dimentional coordinates (either polar or cartesian) -- i.e: each complex number is made up of two elements → {x,y} or {θ,r}, whichever way you prefer to think about it. So now our index of refraction will be made up of two values instead of just one, and for the purposes of interpreting their effect on the reflectivity of a material, we can think of them this way: the first component (x) represents the material's index of refraction in exactly the sense that is normally given to that term, and we can think of the second component (y) as the "spread" of reflectivity -- as this value increases, the reflectivity becomes more and more homogenuous (i.e: more spread out; less difference between parallel and perpendicular incidence angles). If the spread (or y-component) is left at zero, then our complex index is identical to the more familiar real-valued one.
In more technical terms, the second component of the complex-valued refractive index (its "imaginary" part) represents the extinction coefficient of the material (usually denoted with the letter k ), and it determines how fast the amplitude of the wave decreases. The extiction coefficient is directly related to the absorption of the material, which is how fast the material absorbs light, and is given by the absorption coefficient a . These two coefficients, as well as the index of refraction η (the real part of the complex index) are all optical constants that define how a material reacts to incident light. The extiction k and absorption a coefficients relate to each other as a = 4πk / λ , where λ is the wavelength.
In order to be able to use this new complex index of refraction, we need to rewrite the equation in such a way that it can use it, which in turn means that we need to avail ourselves of the necessary tools to manipulate this new type of number. In fact, there's a whole algebra defined for complex numbers, so that we can add them, subtract them, multiply them, and so on, in the same way that we are used to applying these operations to the real numbers. There is now a package of VEX functions available in the wiki library that can be used to do complex number arithmetic. We will make use of these functions in our implementation later on. But first, here is a re-formulation of the Fresnel function for unpolarized light that will accept a complex index of refraction:
All materials have an index of refraction larger than 1.0 (which represents absolute vacuum), but most of them are also no higher than 2.5, with the mineral "eglestonite" at η = 2.49 and diamond not far behind. With this in mind, and looking at the plots in figure 21, it should be quite apparent the difference between increasing reflectivity across the angular range (approaching properties of conductive materials like metals) by incrementing η to some arbitrarily large amount as in figure 20 (which shows plots for values of η of up to 11.2), and via the exctinction coefficient k as shown in figure 21. Looking at these last plots, it should also become apparent that a linear increment in k does not translate to a linear increment in "spread", which is something we may choose to linearize for ease of use when working on the implemenation.
[edit] Implementation
There really isn't much more to the implementation of these formulae outside of dealing with the syntax of the target shading language -- i.e: they are direct translations of the symbolic forms given above. The only thing we should ensure is that the angle α is in the range [−π/2,+π/2]. In our test shader we should also include VEX's and RSL's built-in versions of these functions for comparison purposes. But keep in mind that when comparing the results of our custom formula to the built-in ones, we should keep the second component of our complex index of refraction (the "spread" parameter) at zero so that it matches the real-valued index expected by those functions.
#include <wikiComplex.h>
//------------------------------------------------------
// Fresnel factors (kr and kt) using full formula for
// unpolarized light.
//------------------------------------------------------
float fresnel_kr(float cosalpha; complex ior) {
float kr = 1.0;
if(cosalpha>0) {
float c2 = cosalpha*cosalpha;
complex g2 = cx_addR(cx_subR(cx_mul(ior,ior), 1.0), c2);
float ag2 = cx_abs(g2),
a = sqrt(0.5 * (ag2 + g2.x)),
b = sqrt(0.5 * (ag2 - g2.x)),
b2 = b*b,
s = (1.0 - c2) / cosalpha,
amc = a-cosalpha,
apc = a+cosalpha,
ams = a-s,
aps = a+s;
kr = 0.5 * ((amc*amc+b2)/(apc*apc+b2)) *
(1.0 + ((ams*ams+b2)/(aps*aps+b2)));
}
return kr;
}
float fresnel_kt(float cosalpha, ior) {
return max(0.0,1.0-fresnel_kr(cosalpha,ior));
}
//------------------------------------------------------
// Fresnel factors using Schlick's approximation
//------------------------------------------------------
float schlick_kr(float cosalpha, ior) {
float kr = (ior-1.0)/(ior+1.0); kr*=kr;
return kr + (1.0-kr) * pow(1.0-cosalpha,5);
}
float schlick_kt(float cosalpha, ior) {
return max(0.0,1.0-schlick_kr(cosalpha,ior));
}
//------------------------------------------------------
// Test Shader
//------------------------------------------------------
#pragma label model "Fresnel Model"
#pragma choice model "0" "Built-in Fresnel"
#pragma choice model "1" "Custom Fresnel"
#pragma choice m