The paper considers the task of real-time rendering of dynamic relief shadows based on ray casting using origin relief data - detailed height map. The solution proposed is based on looking for shadow rays - sun rays, whose tracks on height map are passed through heights occluding the light. GPU-based methods and algorithms for extracting such rays using an accelerating data structure - a mipmap of maximum and minimum relief heights are developed. This structure provides an effective acceleration of shadow rays extraction by skipping long sections of sun ray tracks that are not involved in relief shadowing. In algorithms developed precise traversing such a data structure is implemented, as well as texture filtering is taking into account, which allows the formation of "torn" shadow edges to be prevented. The solution created was implemented in software complex and a number of comparative shadow visualization tests was conducted. The results of the research can be used in virtual environment systems, video simulators, scientific visualization, educational applications, etc. Relief shadows, height map, ray casting, GPU, virtual environment system
One of the most important tasks of modern virtual environment systems is modelling of realistic
light-shadow appearance [
For such systems, an important feature of developing shadow algorithms is their reliability, i.e. the ability to construct correct shadows for any permissible directions of light sources, as well as the absence of freezes of visualization loop in case of unintentional errors in input data.
From the point of ensuring the correctness of obtained shadows, the best results are derived by algorithms based on ray casting. However, this approach requires the implementation of a per pixel ray processing loop, that, firstly, is associated with heavy time costs, and, secondly, carries potential risks of hanging the visualization system (obtaining an infinite loop). In this paper, we propose algorithms of original implementation of GPU-accelerated ray casting, which ensure reliable real-time relief shadow visualization.
In the field of real-time shadow rendering, a significant part of researches is devoted to the
development of shadow mapping approach. One of the most common is the method of cascaded shadow
mapping [
2021 Copyright for this paper by its authors.
applications. At the same time, the method requires careful adjustment of the shadowing algorithm for
a target virtual scene in order to avoid the effect of "blocky" boundaries of shadows (especially
elongated ones), as well as other undesirable artifacts (truncation of shadows, erroneous self-shading,
etc.). To overcome these disadvantages, solutions based on percentage-closer filtering [
Second major approach is based on constructing and rendering of accurate polygonal model of
shadow volume for each object in the scene [
It is important to note that both at first and second approaches, the reliability of relief shadow rendering will depend not only on shadow algorithm, but also on the reliability of the methods of constructing and rendering of a 3D relief model. These methods may contain hidden defects or optimizations that are invisible during the visualization of the 3D model, but revealing when shadows are rendered. In "Results" section we demonstrate such situation using a simple example.
A completely different principle of shadow rendering is used in the third approach, which involves
working directly with the primary data of 3D model, – a height map of its surface [
The authors of the research [17] parallelize on the GPU execution of ray processing loops using the CUDA hardware-software architecture, however, implemented uniform sampling of height map along the ray can lead to missing some texels and loss of shadow quality. In [18], the task of relief mapping based on ray casting is considered, where looking up ray - height map intersections is accelerated by means of pre-calculated cone map. A significant limitation is great increase of cone map construction time depending on height map size (about 15 minutes for 5122 and about 8 hours for 10242 [19]). The paper [19] describes the task of real-time visualization of dynamic height maps based on ray casting accelerated using a maximum mipmaps calculated "on the fly" on the GPU. The authors of the work [20] extended these ideas of ray casting acceleration to shadow rendering, proposing own implementation algorithm. As the authors mention, their algorithm suffers from insufficient accuracy leading to appearing fake shadows and missed terrain parts, which has to be compensated by adding some pre-displacement to terrain vertices.
Our work takes its origins from our earlier researches on ray casting implementation [21, 22]. Based
on the experience gained, in this work, a core of our shadow ray casting system was developed – a
reliable and accurate algorithm of traversing height map texels along the ray, adapted for the GPU. The
key features of the algorithm are original "transit points" method for accurate height map sampling,
which prevents texel missing (and potential loss of shadow information), as well as performs an
incremental integer transition from texel to texel, which ensures guaranteed progress along the ray and
prevents freezes of ray processing loop due to machine error in the representation of real numbers. An
important advantage of our algorithm over [
During the development of shadow ray casting acceleration method, we, regardless studies [19, 20], logically came up to a data structure similar to maximum mipmap [19], where not only maximum heights are stored, but minimum heights too, as well as improved construction method is used, providing better stitching the mipmap with height map. Checking maximum heights allow large parts of height map, not involving in shadow formation, to be effectively skipped, while checking minimum heights allow ray processing to be earlier completed with shadow formation. If ray casting algorithm contains any inaccuracies or errors, the second check will give fake shadows that are immediately noticeable on the image. This unique property of our accelerating data structure allowed many inaccuracies of ray-casting algorithm to be identified and fixed, and, ultimately, a reliable shadow ray casting system performing entirely on the GPU, to be created. The proposed solution is implemented in C++ using the GLSL shading language and the OpenGL graphics library.
First, we will describe definitions of detailed height map, shadow ray and shadow coefficient, as well as some necessary mathematical background, used in our solution.
Detailed height map. Let there be relief section (for simplicity - a unit square), where the height
above sea level for each point is known. Height map is a single-channel 2m × 2m texture, where height
values (in grayscale), sampled by the step of 1 / (2m - 1), are written. Figure 1 shows an example of
constructing 4 × 4 height map. In our work height values are extracted from map by means of bilinear
interpolation [23]. So, for an arbitrary height map point with texture coordinates s, t ∈ [
ℎ( , ) = (ℎ , (1 − ) + ℎ , +1 )(1 − ) + (ℎ +1, (1 − ) + ℎ +1, +1 ) , (1) where ℎ , , ℎ , +1, ℎ +1, , ℎ +1, +1 are height values of four nearest texels forming a square around the point under consideration (see Figure 1). As it can be seen, interpolated height values will be the closer to origin heights, the higher is m detail degree of height map. Height maps with m close to screen resolution (Full HD, 4K) will be referred to as detailed ones (m = 10,..., 12).
j
In this paper, detailed height maps of relatively small terrain areas are considered, so that for shadow construction task the curvature of sea level surface can be neglected, i.e. considered as flat. An example is the Puget Sound height map [24], being mentioned in many papers related to terrain modeling.
Shadow ray. Let there be relief model given by height map and the normal n to sea level surface, as well as infinite light source (the sun) given by direction s to the sun (unit vector). The coordinates of vectors s and n are set in relief model Object Coordinate System (OCS). Let's choose some two points A and B of relief model, lying in the same plane with vectors s and n, and determine the shadow at point A caused by the relief from point B (see Figure 2a). Denote by bshadow boolean flag of shadow presence at point A (0 - no, 1 - yes). Cast the ray v from point A to point B. Figure 2a shows that the point A will be in the shadow, if the vector s is directed no higher than the ray v. Define this using mixed product p [25] of the following vectors
= ∙ ( × ), (2) where r = s × n is unit vector, complementing vectors s and n to the right triple. Then the expression of bshadow flag can be written as ℎ = { 1, ≤ 0 0, ℎ ( ∙ ) ≠ 1, .
(3) If bshadow is 1, then v is shadow ray.
i
ℎ3,0
ℎ0,0
ℎ3,3
ℎ0,3
brightest one, and kdark, kbright ∈ [
→ ⁄2), and, conversely, shadows are brighter, when the sun is at zenith (
→ 0), what is explained by the contribution of ambient light. Note, that thresholds kdark and kbright are specified based on brightness range used for virtual scene rendering. This is especially actual for modern visualization systems, where rendering is performed in high dynamic range [26]. n A
kshadow
Consider the task of constructing shadows on the image of relief 3D model, synthesized at visualization stage. Let’s choose an arbitrary pixel of this image. Chosen pixel corresponds to some point P of relief model and point P’ on height map. Denote by s’ the track of sun ray drawn from the point P’ (in texture space), and by T - the set of texels intersected by the track s’ (see Figure 3). Each Ti,jth texel has corresponded interpolated point Pi,j of relief model (honestly, each (i, j)th section of the track s’ has corresponded interpolated curve on relief model, however, this can be neglected on detailed height map). Denote by vi,j the ray casted from point P to point Pi,j. As it can be seen, the point P will be in the shadow, if at least one of all vi,jth rays is shadow ray. Then to solve the task, every pixel of relief image, having vi,jth shadow ray, should be extracted, and its brightness decreased by kshadow times.
T (0,1)
i (0,0) ′ , j (1,1)
′ (1,0)
S t0
T (0,1) i+1
i (0,0) s0 (1,0)
S
In this work, pixels of relief image that to be shadowed, are extracted completely on the GPU, in parallel, independently of each other using developed fragment shader. This shader checks the texels of the "footprint" of the track s' either until the first Ti,jth pixel having shadow ray vi,jth will be detected, or until the "footprint" will end. This checking includes the following steps: (a) extracting height value from the Ti,jth texel; (b) constructing vi,jth ray and checking its belonging to shadow rays, and (c) moving to the next Ti,jth texel. Step (b) can be found in the previous Section, so we focus on steps (a) and (c).
To avoid the formation of "torn" shadow edges, we extract height value from the Ti,jth texel using bilinear interpolation (see Eq. (1)) and take each sampling point Ps as the average of two transit points: input point Pin and output point Pout of the track s' in the Ti,jth texel (see Figure 4). As it can be seen from the figure, the point Pout of the Ti,jth texel will be the point Pin for the Ti+1,jth texel, thus, when checking each Ti,jth texel, it is enough to calculate only one transit point Pout. Denote by u unit vector of track s’ direction, then coordinates of the point Pout can be written in parametric form Pout = Pin + ut. The calculation of t parameter, as well as shifts (ioffset, joffset) of next texel along the track (step (c)), are implemented in developed algorithm A1.
1. Write texture coordinates array K for 4 corners of Ti,jth texel:
K[
t = -1; (ioffset, joffset) = (0, 0); exit the algorithm. 3. Calculate the number ncorn of the corner of Ti,jth texel, which is intersected by track s’: ncorn = b0 + 2b1, where b0, b1 are boolean flags, b0 = (|ux| ≥ 0), b1 = (|uy| ≥ 0). 4. Write the difference dP = K[ncorn] - Pin. 5. If (|ux| ≤ ), then: t = |dPy|/|uy|; (ioffset, joffset) = (sign(uy), 0); exit the algorithm. 6. If (|uy| ≤ ), then: t = |dPx|/|ux|; (ioffset, joffset) = (0, sign(ux)); exit the algorithm. 7. Write tVec = dP/u and boolean flags b3 = (tVecx ≤ tVecy), b4 = (tVecy ≤ tVecx). 8. t = min(tVecx, tVecy); (ioffset, joffset) = (sign(uy) ∙ b4, sign(ux) ∙ b3).
End of the algorithm.
In practice, visiting all texels is not always necessary and effective, since height maps almost always contain large areas (plains, valleys, etc.) not involved in shadow constructing (provided not very low sun altitude). We localize such areas and skip them when our checking algorithm reaches such area border. To localize the areas, the mipmap {hmax, hmin} of maximum and minimum heights is created. It is an ordered set of two-channel textures (R and G components) with the following properties: 1) For 0th map (the first one): - dimensions are two times smaller than height map - 2m-1 × 2m-1 (similar to classic mipmaps [27]); - each (i0, j0)th texel stores hmax, hmin values sampled from height map using 3 × 3 matrix (see Fig. 5a); 2) For kth map (except 0th): - dimensions are two times smaller than (k-1)th map, but 1 × 1 map is not created; - each (ik, jk)th texel stores hmax, hmin values sampled from (k-1)th map using 2 × 2 matrix (see Fig. 5b). i0
T i 8 7 6 1 0 5 j 2 3 4 (a) j0
S ik
T ik-1 (b) 1 0 jk-1 jk 2 3
S
Note, that in 3 × 3 matrix sampling points are upper right corners of 0,…,8th texels of height map, and sampling is performed using bilinear interpolation (GL_LINEAR, see Eq. (1)). In 2 × 2 matrix sampling points are centers of 0,…,3th texels of (k-1)th map, and sampling is performed using the "nearest neighbor" principle (GL_NEAREST). The mipmap {hmax, hmin} is constructed "on the fly" on the GPU by means of developed fragment shader once before visualization.
Let the mipmap {hmax, hmin} contain just 2 maps (for simplicity). Introduce the following definitions (see Figure 6):
- l is level of detail (LOD): 0 is height map, 1 and 2 are 0th and 1st maps of mipmap {hmax, hmin}; - (i, j) are (row, column) numbers of start texel at 0th LOD; - (icur,l, jcur.l), (inext,l, jnext,l) are (row, column) numbers of current and next texel along the ray at lth LOD; - Pout is output point of the ray in start texel; - Pout,l is output point of the ray in (icur,l, jcur.l)th texel; - P is the point on relief model (in OCS), corresponding to P’; - Pmin,l, Pmax,l are the minimum and the maximum points of relief model (in OCS), corresponding to (icur,l, jcur.l)th texel (at 0th LOD Pmin,l equals to Pmax,l).
The progress along the ray, accelerated by means of the mipmap {hmax, hmin}, is implemented in our developed algorithm A2.
1. Initialization.
Calculate P and (i, j), as well as kshadow by equation (4).
Set Pin = P’ in the algorithm A1.
Calculate Pout and (inext,0, jnext,0) by means of algorithm A1.
Set Pin = Pout in the algorithm A1. Set l = 2, k = 1. 2. Ray-casting loop. height map dimensions: Do loop while at least one of (icur,0, jcur.0) will go beyond
Calculate (icur,l, jcur,l) from (inext,0, jnext,0).
Calculate Pout,l and (inext,l, jnext,l) by means of algorithm A1.
If l > 0, then do minimum height test:
Calculate Pmin,l based on Pout, l and hmin sampled from (icur,l, jcur,l)th texel.
Calculate for the ray v = (Pmin,l - P) the flag bshadow by equations (2) and (3);
If bshadow equals 1, then k = kshadow, exit the loop.
Do maximum height test:
If l equals to 0, then:
texel.
Otherwise:
Calculate Pmax,0 based on Ps (see Fig. 4) and h sampled from (icur,0, jcur,0)th = 2 , 2 = 1 = 0 ,2 ,1 ,0
′ ′ ′
Calculate for the ray v = (Pmax,l - P) the flag bshadow by equations (2) and (3); If bshadow equals to 1, then: Otherwise:
If l equals to 0, then k = kshadow, exit the loop.
Otherwise: l = l-1.
Set (icur,l, jcur,l) = (inext,l, jnext,l). Calculate (icur,0, jcur.0) from Pout,l.
Set Pin = Pout,l in the algorithm A1. Set l = 2. 3. Set pixel color: C = kC.
End of the algorithm.
Based on methods and algorithms proposed, a software complex for real-time rendering dynamic relief shadows was developed. The complex is a part of the VirSim virtual environment system [28], developed at the Federal State Institution "Scientific Research Institute for System Analysis of the Russian Academy of Sciences".
We carried out a number of comparative shadow visualization tests using common cascaded shadow
mapping method (CSM) varied in configurations, and our shadow ray casting implementation. For this,
we integrated our solution into well-known NVidia CSM implementation, available at [
Figure 7 shows examples of cases where "spike"-test reveals hidden CSM method drawbacks: (a-c) some shadows are torn off from spikes or are passed through them; (e-g) sawtooth shadow edges are appeared (when increasing CSM resolution, edge appearance looks better, but sawtoothness remains still noticeable); (i) conical shape of spike shadow is lost (even at 4K CSM). Our ray casting implementation produces accurate spike shadows with smooth edges in all the above cases (d, h, j).
Figure 8 illustrates results of reliability test where shadows of steps model, having front (visible) faces only from one side, are visualized. Due to enabled back face culling optimization, CSM method failed to produce some of the shadows (left image). In our implementation shadows are produced completely (right image), because object's height map is directly used, and there is no dependency on 3D model construction and visualization methods.
As tests on the Grand Canyon and the Puget Sound height maps showed, dynamic shadow loss is possible in CSM method at long (see Figure 9a) or middle (see Figure 10a, along the lake) distances due to insufficiently close/distant location of intermediate clipping planes. For the same reason, a sawtooth silhouette of elongated mountain shadows may appear (see Figure 10c). In our ray casting implementation all relief shadows are accurately drawn with high detailing, regardless of distances they are located at and sun altitude (see Figures 9b, 10b, 10d).
All tests described were performed at 1920×1080 screen resolution on personal computer (Intel Core i7-6800K 3.40GHz, 16Gb RAM, Windows 10 Pro) equipped with graphics card NVidia GeForce RTX 2080 (2944 cores, 8 GB VRAM, 16x anisotropic filtering, 8x anti-aliasing). Table 1 contains results of performance measurement of the Puget Sound shadow rendering for both methods.
The publication is made within the state task of Federal State Institution “Scientific Research Institute for System Analysis of the Russian Academy of Sciences” on “Carrying out basic scientific researches (47 GP)” on topic No. FNEF-2021-0012 “Virtual environment systems: technologies, methods and algorithms of mathematical modeling and visualization. 0580-2021-0012” (Reg. No. 121031300061-2).
[17] T. Aslandere, M. Flatken, A. Gerndt, A Real-Time Physically Based Algorithm for Hard Shadows on Dynamic Height-Fields, in: Proceedings of 12. Workshop der GI-Fachgruppe on Virtuelle und Erweiterte Realität, Aachen Verlag, Bonn, 2015, pp. 101–112. URL: https://elib.dlr.de/101497/. [18] F. Policarpo, M.M. Oliveira, Relaxed Cone Stepping for Relief Mapping, in: H. Nguyen (Ed.), GPU Gems 3, Addison-Wesley Professional, 2007, pp. 409–428. URL: https://developer.nvidia.com/ gpugems/gpugems3/part-iii-rendering/chapter-18-relaxed-conestepping-relief-mapping. [19] A. Tevs, I. Ihrke, H.-P. Seidel, Maximum Mipmaps for Fast, Accurate, and Scalable Dynamic Height Field Rendering, in: Proceedings of the 2008 Symposium on Interactive 3D Graphics and Games (I3D ’08), New York, NY, 2008, pp. 183–190. doi:10.1145/1342250.1342279. [20] D. Jung, F. Schrempp, S. Son, Optimally Fast Soft Shadows on Curved Terrain with Dynamic
Programming and Maximum Mipmaps, 2020. URL: https://arxiv.org/pdf/2005.06671.pdf. [21] M.V. Mikhaylyuk, D.A. Kononov, D.M. Loginov, The method for calculating the set of cells in 2D square grid intersected by a given ray, Trudy NIISI RAN/Proc. SRISA RAS 9.1 (2019): 44-48 (in Russian). doi:10.25682/NIISI.2019.1.0006. [22] P.Yu. Timokhin, K.D. Panteley, E.M. Vozhegov, A.M. Trushin, Constructing on the GPU maximum intensity projection of 3D scalar fields of digital core material model, Trudy NIISI RAN/Proc. SRISA RAS 9.3 (2019): 58-65 (in Russian). doi:10.25682/NIISI.2019.3.0008. [23] M. Segal, K. Akeley, The OpenGL Graphics System: A Specification, Version 4.6, Core Profile, The Khronos Group Inc., 2006-2019. URL: https://www.khronos.org/registry/OpenGL/specs/gl/ glspec46.core.pdf. [24] Large Geometric Models Archive (Puget Sound, Grand Canyon), Georgia Institute of Technology,
URL: https://www.cc.gatech.edu/projects/large_models/. [25] H. Lass, Vector and Tensor Analysis, International series in pure and applied mathematics, 1st.
ed., McGraw-Hill Book Company, Inc., New York, NY, 1950. [26] M.V. Mikhaylyuk, P.Y. Timokhin, M.A. Torgashev, The Method of Real-Time Implementation of Tone Mapping and Bloom Effect, Programming and Computer Software 41.5 (2015): 289-294. doi:10.1134/S0361768815050084. [27] E. Lengyel, Mathematics for 3D Game Programming and Computer Graphics, 3rd. ed., Course
Technology, Boston, MA, 2012. [28] M.V. Mikhaylyuk, A.V. Maltsev, P.Y. Timokhin, E.V. Strashnov, B.I. Kryuchkov, V.M. Usov, The VirSim Virtual Environment System for the Simulation Complexes of Cosmonaut Training, Scientific Journal Manned Spaceflight 4 (2020): 72-95 (in Russian). doi:10.34131/MSF.20.4.72-95. URL: http://www.gctc.ru/media/files/Periodicheskie_izdaniya/ppk_2020_4_total_37/5_stat.a_mihailuk_pq.pdf