tfEvenly-Spaced Streamlines in 3D
using Illuminated Stream Lines
Motivation
Overview:
The purpose of this project is to combine two methods of flow visualisation. First, we use the "Evenly-Spaced Streamlines" algorithm(ESS) by Bruno Jobard and Wilfrid Lefer to evenly distribute streamlines over a 3D vectorfield. This algorithm produces excellent results for 2D pictures and can be extended to 3D in a rather straight-forward way.
As it is difficult to distinquish length, orientation and position of streamlines in a 3D picture because of the perspective, we want to improve spatial perception by adding illumination from a light source. Unfortunately, popular lighting models like Phong shading cannot easily be applied to streamlines, because lines have abitrary many normal vectors. M. Zöckler, D. Stalling, and H.-C. Hege. developed a illumination method which doesn't depend on the normal vector, but entirely on the tangent vector of the streamline, and they called it "Illuminated Stream Lines". By using texture mapping and texture matrix, a graphics library like OpenGL can be used to render the illumination correctly.
Mapping of scalar values on streamline attributes:
Besides from combining two techniques of flow visualisation, we concentrated our efforts on the effective mapping of scalar values (velocity, pressure) on streamline attributes like line width, opacity, color, streamline material and streamline density.
For mapping streamline density interactively, we developed a method that we called "Incremential streamline
resolutions". This is an
extension to the "Evenly-Spaced Streamlines"-algorithm,
with the number of resolutions as an additional
input parameter. It works exactly like "Evenly-Spaced
Streamlines", with the exeption that it does not calculate all streamlines
for the highest resolution at once.
Instead it first calculates the streamlines for the lowest resolution. Then it
reuses the calculated streamlines for the next resolution(by putting them back
into the streamline queue => see "Evenly-Spaced
Streamlines"-section). This step is repeated until the streamlines for the
highest resolution are calculated.
In pseude-code, the algorithm would look like following:
rn = number of resolutions.
dsep = separating distance between streamlines of highest resolution.
q = ESS-streamline queue.
l = list of already calculated streamlines.
ESS = Evenly-Spaced Streamlines algorithm
q := {};
l := {};
for(int i=1;i<rn + 1;i++)
{
Put already calculated streamlines into q.
Apply ESS for separating distance dsep * rn / i;
Assign newly calculated streamlines a resolution index of i;
Add new streamlines to l;
}
Output:
l = {streamlines with resolution index};
q = {};
For example, if the resolution number is 4, the algorithm calculates streamlines for 25%, 50%, 75% and finally 100% of the highest streamline resolution and assigns them a resolution index of 1-4.
Because streamlines of lower resolutions are reused for higher resolutions, there is no cutoff of lower resolution streamlines in regions of higher streamline density. Therefore, with enough resolutions calculated, smooth and continous changings of the streamline density according to some scalar value can be achieved.
Problems:
The main problem of our project lies again in the perspective. Whereas the ESS-algorithm works good in 2D, the evenly distribution of the streamlines is destroyed in 3D. There are regions with many streamlines and regions with few streamlines, and this regions change as the view direction changes. The ultimate goal is to achieve "Screen-Space Evenly-Spaced Streamlines".
Detailed description
The idea behind the algorithm is to derive all possible streamlines in a distance dsep from an existing streamline. For that reason reason we use a queue. The algorithm goes like this:
1) We place the first streamline in the queue.
2) We take one streamline out of the queue, calculate all possible seedpoints for a new streamline in distance dsep, integrate the streamlines and place them into the queue.
3) Repeat step 2 until the queue is empty.
For the streamline integration we use another parameter dtest
(test distance). Dtest is in the range of 0.0 to 1.0, and is called test
distance, because the streamline is tested on every integration step if it still
is in a distance >= dsep*dtest from all
existing streamlines. If not, the integration of this streamline is stopped. For
the integration a Euler, Runge-Kutta-2 or Runge-Kutta-4 operator can be used,
for example.
To avoid having to make the distance test with every samplepoint of every
streamline, which would be very time consuming, an additional grid called
dsep-grid is used. The grid has cell size of dsep, and (pointers of) all
calculated samplepoints are inserted into this grid on the right position. Now
for every streamline seedpoint, samplepoints closer than dsep can only be in the
cell where the seedpoint itself belongs, or in the neighbouring dsep-grid cells.
Therefore the test must take place only for the samplepoints in that cells(9
cells in 2D and 27 in 3D).
For 3D, the only additional difficulty compared to the 2D case is that of the new streamline seedpoint calculation. In 2D, you just iterate through every point of the streamline, and consider the 2 points in a distance dsep which are normal to the streamline tangent vector as new seedpoint for a streamline. In 3D, as I said before, there are arbitrary many normal vectors and therefore abitrary many seedpoints. Therefore the following approach is made: One normal vector is chosen and then rotated around the streamline for n-1 times, with streamline tangent vector as the rotation axis. So we get n candidates for a streamline seedpoint.
Illuminated Stream Lines
1) Mathematical background
We use following notations:
I = light intensity
L = light vector
V = view vector
N = normal vector
T = tangent vector
R = reflection vector = the vector in the L-N-plane with the same angle to N as L.
n = coefficent regulating the highlight intensity
k = diffuse, ambient and speculat color components
With the Phong shading model, the light intensity on any surface point is calculated by
Unfortunately, for line segments
there are no unique R and N , but infinitely many, so which normal
and reflection vector should we choose? We take the normal vector coplanar to
L and T as N, and we take the reflection vector coplanar to
L and T as R.
With this particular N, we can now project L into the lines
tangent and normal space, yielding an orthogonal decomposition of L:
Because the decomposition is orthogonal, we can totally eliminate N from the diffuse reflection term with the help of Phythagoras's theorem:
We can use a similar approach to eliminate N from the specular reflection term either:
2) Hardware rendering
Diffuse term:
The illumination can be rendered using texture maps and the texture transformation matrix of OpenGL. Remember, the diffuse term of the Phong equation is given by
We want to express the L-part of the dot product with the texture matrix and use the tangent vector as texture coordinate. The texture matrix is written as
Then the entry point into a one dimensional texture map(note that the term is always between 0 and 1) is given by
Specular term:
If we also want to add the specular term of Phong's model, the view vector V comes into play, and we have to rewrite the texture matrix like this:
Then we get an entry point into a two dimensional texture map:
Now we can use following formulas to create a 2-dimensional texture map:
The program: Flow3D
The program is called Flow3D and is written in C++ with MFC-support. For a detailed description please read the source documentation and the user manual. The documentation, the executeable and the source code are free for download.
Files available for download:
1) Source code: Flow3D.zip
2) Executeable: Flow3D_exe.zip
3) Source documentation: Flow3D_doc_html.zip
Rendering results:
The development of the streamline rendering can be seen here.
Except from visualisation of fractal attractors, there are also some visualisations of real datasets which can be inspected here.
Some tool to render volume datasets as flow-field can be inspected here.
Memory consumption:
Every grid point is stored as a samplepoint. A samplepoint contains the following data:
Position: 3 floats
Vector: 3 float for unit vector, 1 float for magnitude.
(Unit vector can be directly used for illumination calculation)
Quadratic minimal distance(to other samplepoints): 1
float.
Additional scalar information: 1 float each.
Consider there is one additional scalar information like pressure => every samplepoint needs 9 floats and one pointer = 38 bytes. An average grid size of 100x100x100, for instance, would need 36 000 000 bytes just for the sample points = ca. 38 MB.
Streamline calculation adds new samplepoints. The number of the new samplepoints depends on the choosen separating distance and integrator step size, so no exact consideration can be made. For every new samplepoint, one pointer to this samplepoint is stored in the dsep-grid.
The actually displayed samplepoints is usually a much lower number that the samplepoints calculated, especially if dt is low, because a minimal distance to previous samplepoint must be given for a samplepoints to be rendered.
Two adjacent displayed samplepoints are put together in a
streamline piece, and in the rendering step the streamlines are rendered by
drawing lines between the two samplepoints of a streamline piece. Besides the
two pointers to the samplepoints, a streamline piece contains a couple of
rendering information about the line, like color, line width, opacity, scalar
region index, resolution index and so on.
Running the program, I made following observations of the memory consumption:
After starting up, the application needed
8148KB. Calculating a grid of
100x100x100, the application needed
55576 KB.
After the streamline calculation with reasonable values, the application needed
58436 KB => Streamline calculation needed
much less than the grid computing.
Oliver Mattausch, 9506125
Department of Computer Science TU WIEN