Haven’t produced any posts in a while (well at least not published any), I guess that’s just an inevitable consequence of personal blogs! However, the explanation, as it usually is, is that I have been extremely busy with my theoretical studies and actually still am. I’m writing my dissertation over the summer. About the topic, well I might as well just post the project formulation. Dissertation Formulation

To be fair, the graphics side of my work has been suspended since last October! Crazy indeed, I guess physics won this interval of time. Anyway, I hope to produce something shiny in the autumn. For now, I can say that I’m pretty excited about what is going on over at Echobit.

And as a final note; I have no clue of what is going on with the use of LaTeX on this site, though I hope to get it fixed.

The derivation of the matrix for a general rotation in three dimensional Euclidean space is often encountered in strange forms. This is sometimes a bit puzzling until you have found your own way around it. In this post, I will not derive the expression, but show how the useful matrix exponential of the SO(3) generators is in fact producing a rotation about an arbitrary axis through some angle. [Note: SO(3) is the group of 3x3 orthogonal matrices with unit determinant].

Let the rotation matrix about an arbitrary rotation axis through an angle be denoted . The matrix satisfies and we can therefore pick and let be arbitrary. This way any double covering is avoided.

The infinitesimal generators of the SO(3) group is can be written as:

; ;

Which for example can be obtained by considering the generators of SO(2) (which is far more easier obtained if you are starting from scratch) and generalize them to three dimension. Note that you can form Hermitian matrices by defining . The generators satisfy the Lie algebra , where is the Levi-Civita symbol.

Consider the matrix , where the repeated index i is to be summed over. Now, since the generators are antisymmetric we can write: for some matrix . This is a general result from linear algebra, that any antisymmetric matrix can be decomposed like this (c.f. http://en.wikipedia.org/wiki/Antisymmetric_matrix). Hence,

The neat thing is that the exponential matrix containing is easy to evaluate. In fact, we probably already know the answer from the derivation of the generators in the first place.

This matrix is the well known matrix which produces rotations about the z-axis. The matrix can be expressed in terms of the generator by realizing that , where is the identity matrix:

Putting this result into the relation:

Which is the expression for a general rotation in three dimensional Euclidean space! Oh, that was easy to show, if you knew all the right algebra! Cheers.

In group theory Lagrange’s theorem is one of the most important in the theory of finite groups.

Theorem 2.6 (Lagrange’s theorem) The order of a subgroup H of a finite group G is a divisor of the order of G, i.e. |H| divides |G|.

One evident, but funny, implication of the theorem is in the answer of the following question: List all of the subgroups of any group whose order is a prime number.

Solution: According to Lagrange’s theorem, the order of a subgroup H of a group G must be a divisor of |G|. Since the divisors of a prime number are only the number itself and unity, the subgroups of a group of prime order must be either the unit element alone, H = {e}, or the group G itself, H = G, both of which are improper subgroups. Therefore, a group of prime order has no proper subgroups.

I decided that it was probably time for a release of a graphical demo! Since I did not code up anything in relation to the participating media research I thought that this time I would do the research the other way around! This time I have done some research related to light scattering from large spheres. When considering large spheres, it is vaguely speaking possible to use geometrical optics to describe the physics of the scattering pattern. Geometrical optics is also know as ray tracing in computer graphics.

Although, the use of geometrical optics, instead of the full solution to Maxwell’s equations (known as Mie Theory in this relation), does not provide the classical exact solution to the scattering pattern, it will provide an almost correct one. The scattering pattern will contain the forward scattering known as caustics, the primary rainbow, secondary rainbow and so on.

Having computed the scattering pattern from a sphere and stored in a simple 1D table (say floating point texture), it is pretty straight forward to imagine the use in a real-time application, since the scattering of light in any direction from the sphere is now known from a single look up. Actually it is only a 1D table for a infinite distance light source. In order to handle a nearer point light source correctly, it is necessary to compute a two dimensional table, because the finite distance to the sphere excludes some possibilites of the incoming angles depending on the radius of the sphere. Also to get the entire range of rainbow colors one would have to compute the scattering for different wavelengths.

Demo Info

It is possible to use either the 1D or the 2D apporach in the implementation. Make sure you try both, because there is a pretty rainbow hiding! In the demo the sphere has the same properties as a large water drop and the scattering was only computed for three different wavelengths red, green, and blue. No indirect lighting was taken care of (hey it should be able to run real-time). This means that when the light source (which is modeled as being a point light source) is close to the sphere there will be a large ugly unrealistic shadow on the other side of the sphere. In addition, only light coming directly from the light source will be scattered through the sphere. In principle one could add secondary bounces.

Please, be patient with the loading time, the scattering tables are stored as simple ASCII data. The sphere itself is rendered by a combination of the reflection of light and the first bounce of light coming from the environment. The reflection is computed using the scattering tables ones again, while the light from the environment is estimated by rendering a cube map and looking up in the direction of the double refracted ray through the sphere.

Note that the reference is really enjoyable reading for folks with that sense of humor and is highly recommendable.

Download

Sphere Caustics Demo (2.4MB, Windows only)

References

[1] H.C.Van De Hulst, Light Scattering by Small Particles

I am proud to finally publish my article concerning analytic approaches for rendering participating medias. Enjoy!

Abstract

An analytical approach toward visualization of participating media is presented. The volume rendering equation is considered under the assumption that the in-scattering integral can be sampled as a Dirac delta function in the direction of the illuminating light source. This assumption leads to a completely general result for isotropic homogeneous single scattering medias for which the airlight integral is shown to be a special case of. In addition, the limit for which a light source positioned far away is considered, which is particularly interesting in the case of ocean rendering.

The properties of the result are investigated and an approximate model for point light sources is proposed. The presented model is shown to approximate well for thin medias with a close light source and is suitable for real-time rendering. It is capable of simulating homogeneous translucent medias with a low albedo like marble, wax, skin, and smoke.

Download
Analytic Approaches to Single Scattering in Participating Medias (0.3MB)

References
[1] GREEN, S. 2004. Real-time approximations to subsurface scattering,
http://http.developer.nvidia.com/GPUGems/gpugems_ch16.html.
GPU Gems, Chap 16.

[2] JENSEN, H. W. 2001. Realistic image synthesis using photon mapping.
A. K. Peters, Ltd., Natick, MA, USA.

[3] SUN, B., RAMAMOORTHI, R., NARASIMHAN, S. G., AND NAYAR, S. K. 2005.
A practical analytic single scattering model for real time rendering.
ACM Trans. Graph. 24, 3, 1040–1049.

[4] WENZEL, C. 2007. Real-time atmospheric effects in games revisited,
http://ati.amd.com/developer/gdc/2007/d3dtutorial_crytek.pdf.
GDC 2007.

In which direction is the development of relief mapping techniques heading?

Techniques for visualizing surfaces of highly detailed geometry without exploiting the bandwidth to the stress limits by passing a huge amount of primitives down the graphics pipe have been a research topic for quite a while. Currently, techniques such as bump mapping, parallax mapping, and parallax occlusion mapping seem to be the standard techniques that are present in every pipeline of a 3D engine on the market today. Among newer detailed surface techniques we find a large group of relief mapping techniques:

- Per-Pixel Displacement Mapping with Distance Function (GPU GEMS 2)
- Interval Mapping
- Cone Step Mapping
- Relaxed cone step mapping For Relief Mapping (GPU GEMS 3).

All these techniques are based on fixed point iteration schemes (sometimes called ray tracing), which are solved numerically in order to find the actually displacement/offset of the texture coordinates for each pixel. Such approaches are of course much more accurate than the approximations techniques bump mapping and parallax mapping. The advantage of viewing the problem as a fixed point problem is that the problem can be mapped directly onto the pixel shader, which is intuitive and easy!

However, I think techniques which consider the problem more globally will occur in the future. A global approach will solve for all pixels at once and then use the solution to render. I am of course working on such an approach, which most likely will be called inverse mapping. The technique has some relation to my previous post “The Inverse Function Theorem”. The meaning with this post was however to discuss some of the advantages/disadvantages of global approaches versus the fixed point schemes (I will denote the fixed point schemes local schemes from now on), and write about the idea of using mapping techniques for skybox rendering. First a list of some significant advantages and disadvantages:

Advantages: coherent solution, no artifacts, no iteration in pixel shaders, better solutions of deep surfaces, more general computing minded, only needs to be recomputed when the viewpoint changes.
Disadvantages: same time complexity, no culling.

Obtaining a global coherent solution will not contain artifacts as seen in for example the local approach cone step mapping, where solutions sort of blend around the edge of ridges on the surface. The solving technique will differ much from the fixed point techniques and does not match the fix graphical pipeline as well as the these, however GPU’s are beginning to support better languages for general purposes (I am especially referring to CUDA), which the first robust implementation most likely will be targeted towards. One disadvantage is the lack of culling. It is hard to do know which pixels that should be culling since the solution is obtained independent of what is going on in the color buffer, where as the local approaches are performed directly into the color buffer and culling is therefore easy.

I have done some simple CPU and shader implementations which executes fine! I have included a sample picture and a movie of a parabolic cylinder surface, which are rendered using my technique. The white dot shows the position of a point light. Animation Clip (17MB)

The time complexity of the global technique seems to be very similar to the local approaches - but it provides much better quality of the rendering. Furthermore, it should be noted that both local and global techniques lack the support of multiple depth maps, which is a research problem for itself. Allowing multiple depth maps enables the possibility of the surface to fold in the depth direction.

Using a global mapping technique for skyboxes

Everyone that has tried to play an ordinary open scenario FPS game have noticed that the levels of the game are cut off at some point in space. Standing at this cut off boundary one can often see what is called the skybox. My favorite example of a skybox can be found in Half Life 2 Episode 2 (see picture), where the background is not a simple static background, but contains dynamical elements which helps making it a part of the whole scenario. Another good example is the skyboxes used in the game Riddick EFBB – which the Swedish developers Starbreeze actually acquired out of house, but never the less it is some real nice artwork.

Now, imagine having a technique that could present as much geometry you would like at the far distances of the scenario. This technique could be the global mapping technique! It will be able to handle such things and will provide backgrounds of same quality as the rest of the scene! This means that city scenes could render houses of the same quality in the far distances and open spaces with a lot of vegetation could be rendered at the same quality in the far distances. Imagine such techniques in games like: Assassins Creed, Alan Wake, Far Cry 2, GTA IV, the city seen in the Half Life screen, Hitman, The Witcher, to name a few. Of course, the technique that I have discussed have some immediate technical problems, such as the resolution of surface details (normal map, depth maps), the requirement of multiple depth maps, and some difficulties with dynamic environments. But seen in the light of the skybox and LOD techniques (seen in far cry) that are the standards today – this idea would certainly be a step forward!

I was working on methods for determining the inverse of a specific function. Namely, a two dimensional bijective vector function of the following kind:

In the particular problem the function is sampled on a regular grid so only a discrete map is available. It is therefore tempting to find the inverse by simply exchanging the arguments:

The problem with this naive approach is that the inverse function will properly not be represented on a regular grid and even if the function is bijective in its continuous representation the discrete sampling does not have to be so. Finally, neighbors in the discrete representation of the function does not have to be neighbors in the inverse representation, because the function is not necessarily conformal (does not preserve angles). This essentially means that if you want the inverse function represented on a regular grid you will have to do some linear interpolation between the non-regular grid points, which you do not know the neighbors of. It is therefore difficult and computational expensive to reach good results on a regular grid! (I have tried this by using a BFS algorithm to locate the neighbors and then a bresenhams rasterization algorithm in order to compute the linear interpolation at the regular grid points.)

Instead I turned to another method which would completely neglect these problems! Namely, exploiting the inverse function theorem. You might have encountered this theorem in its one dimensional version.

Let f: I->R be a monotonic function, which is differentiable in the interval I with f’(X) ≠ 0 for all x in I. The inverse function f^-1: J -> R is then differentiable and for all y in J the following holds:

In my application I need the generalized version to include the nature of higher dimensions

The Inverse Function Theorem
Let f: R^n -> R^n be continuously differentiable on some open set containing a, and suppose det(Jf(a))≠0 (J being the Jacobian). Then there is some open set V containing a and an open W containing f(a) such that f: V -> W has a continuous inverse f^-1: W -> V which is differentiable for all y in W.

In matrix representation this can be expressed as:

where J denotes the Jacobian.

Using this magnificent theorem it is possible to find the inverse function by simply solving the differential equation. In practice this can be done numerically using various methods (for example using Runge-Kutta methods).

For the specific case of two dimensions one will obtain two pairs of equations for which one will need to solve only one of the pair’s equations simultaneously. This means that it is possible to solve independently in either the x direction or the y direction. The two pairs of equations takes the form:

or

where

The inverse function can then be obtained on a regular grid and the only worry is whether the determinant of the Jacobian of the function is zero, because then the approach breaks down!

Useful Links
http://en.wikipedia.org/wiki/Inverse_function
http://en.wikipedia.org/wiki/Runge_kutta
http://en.wikipedia.org/wiki/Bresenham%27s_line_algorithm

I really wanted to post my newest research a few weeks ago entitled “Analytic Single Scattering in Participating Medias”, but then I got stuck in a dilemma concerning the use of LaTeX in WordPress blogs, which I thought I had solved with mimetex. I will reveal the pre-abstract of the post:

Abstract
Computing participating media efficiently is a tough task, however recent research in real-time computer graphics brings some innovative examples of how to do this. For example the CryENGINE 2 shows some innovation images. Although these approaches might not be correct in a physical sense - the way they fake inhomogeneous medias such as ice, jade, vegetation and skin still produces convincing results into the gaming world.

An analytic approach towards single scattering is explored and an exact expression on how to simulate volumetric effects, like ocean rendering where the surface is illuminated by a distant light source is derived. In addition, a novel approach of how to simulate homogeneous translucent objects illuminated by a near point light source is derived. Both approaches is done correctly in the physical sense.

Once it was done I discovered that I had inserted a lot of mathematical expressions in the text using $$, which is not supported by WordPress! So I thought that a search and replace with the public mimetex link would do the job, but doing so yielded a very poor result! Iam currently trying different things to get a good result, because posting the post as pdf file kind of eliminates the idea of blogging!

However, while trying to solve the problem I thought that I might as well write another post on a completely different topic, namely how to exploit the inverse function theorem to find the inverse of a function. Enjoy!

Before reading this very first blog entry I will like to inform the reader that knowledge about the radiosity algorithm may be helpful in order to follow the various discussions and derivations.

Introduction

In this entry I will show that the classical radiosity equation always can be rewritten such that the linear system has a system matrix, which is symmetric and positive definite. A symmetric positive definite matrix can be Cholesky decomposed, which enables the system to be solved with an arbitrary emission configuration using backsubstitution. This way the radiosity algorithm gets extremely stable numerically and faster than traditionally matrix solvers.

The motivation for making research in the rather old radiosity equation [6] is that a rearrangement of the equation of this type makes the radiosity algorithm competitive to the irradiance caching method [2] used in the photon mapping technique. I will leave that discussion for a later blog entry and for now just show how to rewrite the radiosity equation. I will however leave some final words in the conclusion.

From the Rendering Equation to the Classical Radiosity Solution

The derivation from the Rendering equation [5] to the Classical radiosity solution can be found in many text books (see for example [2,3,4]). To keep it short the idea is to express the rendering equation as an integral over surface locations and exploit that the radiance is constant in all directions making it possible to express the equation using radiant existance instead of radiance. Finally, it is discretized into a linear system of equations yielding:

where , , and denotes the radiosity (radiant exitance), the emitted power per area, and the reflectance of element i, respectively, while denotes the form factor from element i to element k. The radiosity equation can be written in matrix notation as:

This equation expresses then the energy transfer between n pure diffuse elements.

Rewriting the Radiosity Equation

I will show that the system matrix can be rearranged so it is symmetric; thereby reducing the storage requirements to about n^2/2 elements. Furthermore the result can be shown not only to be symmetric, but also positive definite, which enables us to perform a Cholesky decomposition [1]. This is decomposition is really neat, since it is extremely stable numerically and enables one of the fastest ways to solve a linear system of equations exact.
First we note that the equation can be written as,

where is a matrix containing the diffuse reflectance’s in its diagonal. Second, we exploit the reciprocity relation,

,

to form a new matrix, which we know is symmetric:

Using this new matrix it is possible to rewrite the linear system of equations such that the system matrix is symmetric:

where is a matrix containing the areas of the elements in the diagonal. Note that the right hand side is known and very easy to compute since and are diagonal matrices. The symmetric system matrix and right hand side vector now have the form:

and

The reader has maybe by now noticed that I have chosen for all i. There is no particular reason for this, since the matrix still would be symmetric if the self-interaction was included (I have chosen this to simplify further derivations).

The new system matrix is positive definite if and only if all leading principal minors of the matrix are shown to be positive [7]. Let us consider the first three minors:

1. , is always true since the reflectance is positive and the area is positive.

2. , is always true if the reflectance’s and the form factors are constrained to lie in the domain .

3.
, unfortunately this determinant is not necessarily positive!

Fortunately, it is possible to exploit the structure of the matrix equation and some basic rules of matrix determinants to ensure that that this determinant and the higher order leading minors always are positive. Lets consider the case where the determinant of the 3×3 matrix is negative (given by case 3). In this case we can multiple row 3 or column 3 of the matrix without making any impact on the previously computed determinants. It is thereby possible to use the theorem saying:

If the matrix B comes from matrix A by multiplication of a row with a constant k, then

In the present case choosing k=-1 will imply a positive determinant of the 3×3 matrix. A similar argument can be used for the n’th leading minor. It is now just a matter of bookkeeping the multiplied rows, which can be done by introducing a diagonal matrix:

Multiplying the new radiosity equation on both sides yields:

This system looks complicated, but is in fact pretty straight forward and has a system matrix, which always is symmetric positive definite!

Solving Using Cholesky Decomposition

It is possible to perform a Cholesky decomposition on the new system matrix
,
since it always will be symmetric positive definite. However, one important note is that the coefficients of the matrix P should be in cooperated into the Cholesky decomposition and will demand that the arbitrary right hand side vector is multiplied by the respective coefficients. This in cooperation is, however, easily done and the cost of multiplying the P matrix with the right hand side is O(N), which can be combined with the inevitable multiplication of the area and inverse reflectance.

The overall operation count for the Cholesky decomposition is N^3/6 executions of the inner loop (consisting of one multiply and one subtract), with also N square roots [1]. This is not a very pleasing count. Fortunately, the system can be solved with an arbitrary patch emission vector once the Cholesky decomposition is known by performing two executions of the backsubstitution algorithm.

Instead of solving the full radiosity matrix for the exact solution for all elements in a given scene using the above method or some other matrix solving technique, the problem has been confronted with different numerical matrix solvers, which I just want to mention the names of so you can Google them: Jacobi iteration, Gauss-Seidel iteration, Southwell iteration, and progressive refinement. Note that solving the complete radiosity system exact has a storage requirement, which goes as O(N^2), which is unacceptable nowadays where the polygon count (and tessellation into elements) is greater than the memory can handle.

Final Words

Even through the execution count and memory consumptions for the Cholesky decomposition does not seem pleasing at first sight it can be shown that the radiosity at an arbitrary location in the scene can be computed in O(N) time if the full system has been solved in before hand. (where N is the size of the system). The storage requirements can be shown to be in the order of O(N). These facts are the main topics in one of my next blog entries and are the reasons why the radiosity algorithm gets competitive to the irradiance caching method. The whole discussion about form factor evaluation has been left out and I will maybe consider bringing them into the discussion also (since they are in fact the must computationally expensive). Much more about that later!

Another topic that could be interesting is how an arbitrary BRDF can be combined with the pure diffuse lighting solution obtained by solving the radiosity equation without violating the energy conservation. This is actually rather easy, but is not physically correct and I do not know whether it will provide nice images. However, introducing a specular component can not provide anything less interesting to an image consisting of pure diffuse surfaces.

It could also be interesting to consider how point light sources can be introduced correctly into the radiosity algorithm. Points in a finite element method does not make much sense, but in relation to real-time rendering, where point light sources are used a lot, this could be an interesting topic.

References

[1] http://www.nrbook.com/, Section 2.9 Cholesky Decomposition
[2] Henrik Wann Jensen, Realistic Image Synthesis Using Photon Mapping, 2001 http://graphics.ucsd.edu/~henrik/
[3] P. Hanrahan, Rendering concepts, In Cohen and Wallace, Radiosity and Realistic Image Synthesis, pp. 13-40
[4] Andrew Glassner, Principles of Digital Image Synthesis¸ January 15, 1995
[5] James T. Kajiya, The Rendering Equation, 1986
[6] Cindy M. Goral, Modeling the Interaction of Light Between Diffuse Surfaces, 1984
[7] http://en.wikipedia.org/wiki/Positive-definite_matrix

Copyright © Jakob Gath 2007

This is the first entry of this blog! In the time to come I will try and post various thoughts in the subjects of physics, mathematics, and computer graphics regularly, which I hope people will gain some motivation and benefits from. Since our server do not have LaTeX installed I will have to stick with mimetex, which makes the time-dependent Schrödinger equation look like this: