Optimally Fast Soft Shadows on Curved Terrain with Dynamic Programming and Maximum Mipmaps
OOptimally Fast Soft Shadows on Curved Terrain with DynamicProgramming and Maximum Mipmaps
Dawoon Jung , Fridger Schrempp and Seunghee Son Korea Aerospace Research Institute (KARI), Daejeon, Korea{dwjung|seunghee}@kari.re.kr Deutsches Elektronen-Synchrotron (DESY), Hamburg, [email protected] 15, 2020
Abstract
We present a simple, novel method of efficiently ren-dering ray cast soft shadows on curved terrain by us-ing dynamic programming and maximum mipmapsto rapidly find a global minimum shadow cost inconstant runtime complexity. Additionally, we ap-ply a new method of reducing view ray computationtimes that pre-displaces the terrain mesh to boot-strap ray starting positions. Combining these twomethods, our ray casting engine runs in real-timewith more than 200% speed up over uniform raystepping with comparable image quality and withouthardware ray tracing acceleration. To add supportfor accurate planetary ephemerides and interactivefeatures, we integrated the engine into celestia.Sci, ageneral space simulation software. We demonstratethe ability of our engine to accurately handle a largerange of distance scales by using it to generate videosof lunar landing trajectories. The numerical errorwhen compared with real lunar mission imagery issmall, demonstrating the accuracy and efficiency ofour approach.
Modeling soft shadows accurately and efficiently isimportant for interactive rendering of long shadows cast by terrain on curved, planetary surfaces. In ourcase, we are developing an interactive simulator thatmust be able to generate reference imagery for test-ing autonomous landing at the lunar poles. Poten-tial landing sites in these regions are perpetually intwilight, and low sun elevation angles lead to softpenumbrae that are nearly / th the total length ofshadows that are themselves tens of kilometers long.Soft shadows can be computed in a single renderpass using ray casting. While several accelerationtechniques for computing ray-terrain intersections ex-ist [21, 31], such techniques typically are problematicwhen ray casting the surface of a curved base sur-face such as a planet, require expensive precomputa-tion, or result in inaccuracies. Moreover, tracing raysfor soft shadows cannot be terminated early (Sec-tion 3.1), requiring a large number of iterations ofslow texture lookups that dominate the render cost.In this work, we describe how to overcome the slowspeed of shadow ray casting and also introduce asimple optimization for view rays to efficiently andaccurately render planet-scale, curved, rugged, andheavily-shadowed terrain with soft shadows in a real-time, interactive application running on previous-generation hardware. We demonstrate and validateour results using mainly scientific data of the Moon,but later we briefly show how our technique gener-alizes to bodies with atmospheres such as Mars and1 a r X i v : . [ c s . G R ] M a y luto.In essence, we recognized that the problem ofray casting soft shadows is equivalent to finding theglobal minimum of a cost function in as few itera-tions as possible where the cost is proportional tothe shadow ray height above the terrain surface. Wefound that dynamic programming –a recursive opti-mization technique–works well for this case. Ourtechnique relies on a maximum mipmap [31] contain-ing maximum height values in successive mip levelsthat is generated in real-time only for the visible ter-rain area and when the view frustum changes suffi-ciently, in order to guarantee that height above theterrain is minimized for a given subproblem.We also show that view ray iteration step count canbe reduced by a factor of nearly 100 for an equivalentor better intersection accuracy by pre-displacing theterrain mesh in the vertex shader.To summarize our main contributions again:1. Fast, accurate , real-time method of comput-ing soft shadows on spherical geometry using dy-namic programming and maximum mipmaps2. Efficient hybrid vertex displacement and viewray casting3.
Scalable rendering architecture that can handleplanetary distance scales
Section 2 discusses research related to this work. Sec-tion 3 details our soft shadow algorithm. Section 4describes our hybrid view ray casting method. Sec-tion 5 discusses how our ray casting engine fits intocelestia.Sci. Section 6 presents our results, and Sec-tion 7 gives concluding remarks.
Terrain is often rendered using displacement map-ping, a technique that usually involves displacing thevertices of a polygonal mesh or ray casted surfacebased on elevation values looked up from a heightfield texture. Szirmay-Kalos and Umenhoffer [28] provide a comprehensive review of displacement map-ping methods.Self-shadowing on terrain refers to shadows cast onterrain by the terrain itself, and is typically computedusing shadow maps in the case of polygonal displace-ment mapping, or ray casting (hybrid methods alsoexist, e.g., grid tracing [17] which casts rays againsta polygonal grid). Ray marching is a variant of raycasting where the ray is marched in discrete steps.When computing shadows using ray casting, a viewray from the viewer to the terrain surface is cast,and then a shadow ray or bundle of rays is shot fromthe terrain intersection point to each light source todetermine the shadow intensity. Soft shadows areshadows whose edges are softened due to area lightshaving non-zero area, e.g., in the case of the Sun.Completely shadowed portions are defined as umbra,and soft areas are defined as penumbra. The refer-ence method of tracing soft shadows is distributedray tracing [7] (Figure 1a), but this requires castingmany rays to each light source with predictably slowperformance.Recent techniques aim to reduce ray casting costby requiring less ray samples. These include horizonmapping techniques that compute the horizon silhou-ette [20, 26, 32] and cone step mapping [6] in a conestep map normally used for accelerating ray-surfaceintersections is repurposed to compute average visi-bility in a limited set of directions. However, thesemethods are limited by slow precomputation timesand coarse sampling.Uniform sampling [9, 30] takes constant-size stepsalong a single ray until the ray exits the terrainvolume. No precomputation is required, but high-frequency terrain features might be missed unlesssampled at sufficiently high rates (Figure 1b) andeven then the computational complexity is O ( N ) where time to trace the entire shadow interval in-creases linearly with number of steps N .Why then are existing ray-cast soft shadowingmethods so unsatisfactory? Certainly many opti-mized ray-surface intersection methods exist; thereader is referred to maximum mipmaps, relaxedcone step mapping [22] (which must not be confusedwith ref [6] where cone sectors encode horizon in-formation), and precomputed safety shapes [2] for2 rea light (a) (b) (c)Figure 1: Comparison of shadow ray tracing methods. (a) Distributed ray tracing requires tracing multiple,expensive rays. (b) Uniform stepping with analytic occluded light area calculation requires only marchingalong direct line of sight but many samples are still wasted. (c) Our method. Horizontal shading representsmaximum mipmap samples, filtered to prevent aliasing. Arrows indicate points along the shadow ray wheresamples are required. Early samples occur at high mip levels, allowing the algorithm to skip large swaths ofthe total interval.the state-of-the-art. However soft shadows cannot becomputed using these methods on their own, becausethey only solve for the distance to an intersection.However, soft shadow computation is a fundamen-tally different problem than ray-surface intersection.The goal of soft shadow computation is not to findthe distance to an intersection. Instead, it is to findthe occluded fraction of an area light . Indeed,within the solid angle subtended by the area lightthere may be many shadow rays which do not inter-sect a surface.In this work, we show that calculating a softshadow can be cast as a global minimization prob-lem. Global minimization is a well-studied area [5].Classical methods such as dynamic programming [4]and branch and bound are useful when the problemcan be divided into sub-problems [8], while heuristicmethods such as simulated annealing and tabu searchdo not strongly guarantee optimality but can be effi-cient for large problems.Figure 1c illustrates our method. We start withthe same strategy as uniform sampling by tracing asingle ray trajectory toward the center of an arealight source, but use dynamic programming guidedby a maximum mipmap generated cheaply in real-time to find a global minimum shadow cost along thetrajectory in O (log N ) runtime in the worst case, or O (1) in practice. The maximum mipmap generationand optimal ray casting run on the GPU in standard OpenGL vertex and fragment shaders at more than3 times the frame rate of uniform ray stepping withcomparable image quality and without hardware raytracing acceleration. This section describes the method we use to computesoft shadows accurately and efficiently.
Modeling soft shadows is important to us, becauselow solar elevation angles near the lunar poles cre-ate shadows that are tens of kilometers in lengthand penumbrae that are a significant fraction of thelength (Figure 2).Our method computes the fraction s of an arealight that is shadowed by terrain. It does this byshooting a single shadow ray R obj = p + ˆL t from aview ray intersection point p toward the geometriccenter of the area light, and computing the heightdifference between the shadow ray and terrain ∆ h k at strategic points. We observe in Figure 3a that theshadow fraction s k is larger when ∆ h k is small. How-ever, small ∆ h k also always occur trivially just afterthe view ray intersection point ( ∆ h in Figure 3b).Instead, the shadow fraction should also take into ac-count the distance travelled by the shadow ray t k in3 km Figure 2: Color coded illustration of long penumbraecast by lunar crater rims in our renderer (inset: fi-nal render). Sunlight is shining from the top left.Umbrae (blue) can reach over 30 km in length, andpenumbrae (red) over 5 km or nearly 1/7th the totalshadow length in this scene. By comparison, maxi-mum height field resolution is 30 m/texel.the following manner [29, 30]: s k = S (cid:18) ∆ h k t k (cid:19) (1)where S is a function relating ∆ h/t to the sun discocclusion fraction as will be described in Section 3.4. ∆ h/t can be interpreted in terms of the unobscuredsolid angle of a cone of light with apex at the startof the shadow ray.As minimizing ∆ h/t will minimize s , the optimalcost J * ( s * ) might appear to be: J * ( s * ) = min ∆ ht (2)where s * = max( s ) (for brevity we will omit the ( s * ) from henceforth). However, the cost function that weshould use is: J * = min ∆ h (3)Only the numerator should be minimized, otherwise ∆ h/t can be trivially minimized by arbitrarily in-creasing t → ∞ .Since there may be any number of ‘peaks’ along theshadow ray with small but non-negative J k = ∆ h k ,the entire length of the shadow ray from p to thepoint where the ray exits the terrain volume mustbe considered. Next we show how our method solvesthis global optimization problem using dynamic pro-gramming in a single-pass GPU fragment shader. Table 1: Table of symbols. p View ray intersection point ˆL Shadow (light) ray unit vector R obj Object space position vector of tip ofray R tex Texture space position vector of tip ofray h Terrain height in [0 , increasing frombottom to top max h Maximum h within a maximummipmap texel H Height of tip of ray ∆ max h k Difference between max h and HT Length of ray covering one maximummipmap texel t Distance traveled by shadow ray nor-malized to [0 , T ] → [0 , N Height field size N (cid:48) Total number of mipmap levels m Mip level ∆ M Size of mipmap texel k Current iteration step J * Optimal cost over { , , ... , k } As discussed in Section 2, we choose dynamic pro-gramming to solve the shadow ray global minimiza-tion problem. Dynamic programming is an algorithmfor solving global minimization problems recursivelyin cases where the cost to optimize a step dependson previous steps [3]. The solution will be optimalif each step is optimal. In other words, if J * k is opti-mal, then J * k +1 = min J * k will be optimal. The taskthen, is to find a way to generate an efficient policy set π = { J * , J * , ... , J * k } that will satisfy optimality whileguaranteeing that the policy covers the entire shadowray domain. We show that maximum mipmaps canbe used to achieve this. Proof of optimality
The ray parameter t * resultingin the optimal shadow cost can be trivially defined as: t * = arg min t { min ∆ h ( t ) } (4)We will now show how to derive an efficient vari-4 hadow fraction Δh (a) Penumbra geometry Δh t Δh (b) Penumbra cost Figure 3: Penumbra geometry and costant of this that uses dynamic programming to run in O (log N ) time. t * can be calculated by the following dynamic pro-gramming method: t * k = k = 0arg min t (cid:8) min ∆ max h ( t * k +1 + 2 − k ) (cid:9) k > where each step k represents a subdivision by half ofstep.If k = 0 (the base case), t = 1 because theray covers the entire texture. For subsequent k > , we first write t * k = arg min { min ∆ h ( t k,i ) } =arg min { min { ∆ h ( t k, ) , ∆ h ( t k, ) , ... }} where i = { , , ... } and t k,i represents values of t over the in-terval covering all possible texels intersected by theray. ∆ h ( t k,i ) can be expressed in terms of subsequentsubdivisions at k + 1 : min ∆ h ( t k,i ) = min { H ( t k,i ) − h ( t k,i ) } = min { H ( t k,i ) − max { h ( t k +1 , i ) ,h ( t k +1 , i +1 ) , ... }} = min ∆ max h ( t k +1 ) (5)where we write the height of the ray tip in texturespace as a function of t k,i as H ( t k,i ) ≡ R h tex ,k,i where R tex ,k,i = R xy obj ,k,i / R z obj ,k,i + for spherical ter-rain. To see why h ( t k,i ) must take the maximum of { h ( t k +1 , i ) , h ( t k +1 , i +1 ) , ... } , first note that the rayin the interval covered by the i th texel at mip level k necessarily intersects texels i, i + 1 , ... at mip level k + 1 . If ∆ h ( t k,i ) were to minimize t k , then none of ∆ h ( t k +1 ) in the interval covered by the ray in the i thtexel can be greater than ∆ h ( t k,i ) since otherwise theinterval would not minimize t k . Additionally, H ( t k,i ) is monotonically increasing with respect to the ter-rain (otherwise it intersects the terrain and is triviallyshadowed) so h ( t k,i ) must be the largest possible inorder to minimize ∆ h ( t k,i ) .Additionally, when going to the next subdivisionlevel, texel centers become offset by − k . So t * k =arg min (cid:8) min ∆ max h ( t * k +1 + 2 − k ) (cid:9) as claimed. Runtime complexity
To show that t * can be cal-culated in O (log N ) time, note that each subdivisionequates with a mip level and the total number of miplevels is log N where N is the maximum dimension(width or height in pixels) of the height field. Thetotal number of steps required at each mip level isconstant (up to three in our case). Thus the totalruntime complexity is bounded by log N .Empirically however, we find that it is not neces-sary to trace all log N mip levels but only a constant N (cid:48) < log N every time, if we always use a texture sizethat is on the order of the display size (e.g., N (cid:48) = 5 for N = [1024 , ) and noting that most shadowsdo not cover the entire length of the texture. In thiscase, runtime complexity is constant , i.e., O (1) . Maximum mipmaps
On a 1-dimensional terrain,each max h ( t * k +1 + 2 − k ) can be obtained by sam-pling from a maximum mipmap which stores themaximum values of the corresponding collection ofheight field texels in successively decreasing mip lev-els m k = { N (cid:48) − , ... , , } .To prove this, first we define 1-dimensional terrainas terrain that has minimal variation in the directionperpendicular to the shadow ray. In practice this isnot always true but we will later show that the arti-facts resulting from this approximation can be effec-tively mitigated as described in the section on longshadows in Section 3.3.By Equation 5, each h ( t k,i ) = max h ( t k +1 , i ) , ... where h ( t k +1 , i ) , ... are the heights of texels at miplevel k +1 in the interval entirely covered by the texelat mip level k . But this is just the definition of amaximum mipmap texture. Texel traversal
We initially set a shadow ray oflength T to correspond to ∆ M , the vector projectedby the ray into texture space of a maximum mipmap5exel at mip level m = N (cid:48) − , aiming to cover onetexel (see Deriving T in the Appendix for a completedefinition).Then, assuming that the height field curvatureis nearly flat at scales of ∼ T , we use a variant ofDDA [1] to traverse the maximum mipmap texelsthat could be intersected by the ray. This variantdiffers from the ‘canonical’ DDA in that all texelsunderneath the footprint of the ray are visited. Fig-ure 4 compares the two different versions of DDA,and Figure 5 illustrates the artifacts that can resultfrom using canonical DDA. (a) (b) Figure 4: Texels intersected by ray of length T . (a)Canonical DDA sometimes misses a texel. (b) Up tothree texels of size ∆ M/ could be intersected; wevisit them all and find the texel that minimizes J k .As Figure 4b shows, there are up to three inter-secting texels in the interval [0 , ∆ M/ ; we find theone that results in min ∆ max h ( t k +1 ) .The complete algorithm is given in Algorithm 1. Next we discuss several important additional detailsfor scalable rendering including virtual texture coor-dinates, and interval concatenation for handling longshadows.
Virtual texture mapping
For our renderer to scalefrom planetary radii of thousands of kilometers tosurface features of tens of meters (a scale factorof ∼ ) with the limited VRAM of a commoditygraphics card, we employ a tile-based virtual texturescheme and stitch together a subset of visible tiles atruntime. While allowing us to render scenes across awide range of distance scales, this also implies thatcare must be taken to use scaled and offset texturecoordinates of the stitched texture when sampling. (a) (b) -101 (c) (d) -101 Figure 5: Streaks when using canonical DDA (toprow) vs our method (bottom). (b), (d) display J *cost values. Long shadows
We noted in Section 3.2 that wecould use less mip levels than log N as shadows arenot usually very long in practice. But this also has abeneficial effect of mitigating artifacts when approxi-mating max h ( t k +1 ) with a maximum mipmap. Thisis because the effect of an error is amplified at highermip levels k → due to subdivision interval sizesdoubling at each level: at mip level 0 for instance, t = { , , } and supposing that t = 1 was chosen in-stead of t = 0 subsequent steps will always be offsetby 1, while if the algorithm is started at mip level 5then the initial error in t will be ≤ − .However the chosen initial mip level N (cid:48) may notbe sufficiently large to cover very long shadows. In-stead we run Algorithm 1 again, but setting the sam-pling interval immediately after the first interval, i.e., t N (cid:48) + k = t + T + t k . As typically the total number ofrequired iterations N (cid:48) need not be very large, this canbe repeated once or twice as needed. N (cid:48) = { , , } for all scenes in this work unless specified. We implemented a physically-based penumbra calcu-lation that assumes light sources (e.g., Sun, Earth)are circular discs emitting a cone of constant lightflux with the apex at the view ray intersection point.Moreover, we make the simplifying assumption that6 lgorithm 1:
Fast Soft Shadow Algorithm J * ← t ← m ← N (cid:48) − R ← Texel step size for k ← to N (cid:48) − do t * ← − , ∆ max h * ← , i ← Compute height H of ray tip at R ( t k,i ) Sample max h ( t k +1 ) at R ( t k,i ) , mip level m ∆ max h k ← H − max h ( t k +1 ) if ∆ max h k < ∆ max h * then ∆ max h * ← ∆ max h ( t k +1 ) t * ← t k for i ← to do // DDA R ( t k,i ) ← R ( t k,i − ) − − k − ∆ R t k ← t k − − k − Compute height H of ray tip at R ( t k,i ) Sample max h ( t k +1 ) at R ( t k,i ) , level m ∆ max h k ← H − max h ( t k +1 ) if ∆ max h k < ∆ max h * then ∆ max h * ← ∆ max h k t * ← t k m ← m − , k ← k + 1 if t * > − then t k ← t * + 2 − k if ∆ max h * < then J * ← (∆ max h * ) /t *the cross section of the terrain shadowing the lightis a prism, i.e., the terrain height is constant in thedirection perpendicular to the shadow ray. This al-lows us to calculate the shadow fraction s as a cir-cular segment. (Figure 3a). While real light sourcessuch as the Sun exhibit limb darkening, and cross sec-tions of terrain are usually never prisms, neverthelessour disc occlusion method generates more accurateresults over other [29,30], non-disc based approxima-tions. Figure 6 compares disc occlusion and non-discmethods with a calibrated ground-truth image fromthe Lunar Reconnaissance Orbiter (LRO) Wide An-gle Camera (WAC) [27].Defining r L as the apparent light radius, surfacenormal ˆN , normal vector aligned to be perpendicular (a)(b)(c) N o r m a li z ed P i x e l V a l ue s (a) Reference(b) Disc(c) Non-disc Figure 6: Normalized pixel value profiles (scaled be-tween the minimum and maximum value ranges foreach image to highlight slope differences) across lit-shadow boundary for (a) LRO WAC reference im-age [27], (b) our render using disc occlusion and (c)without. (b) tracks the reference profile more closely,while (c) underestimates values near the left side ofthe profile.to the shadow ray ˆN ˆ L , the distance d from the centerof the light disc to the edge of the shadowed segment,normalized to disc radius=1, is: d = 2 J * ( ˆN · ˆN ˆ L ) − r L r L d ∈ [ − , J * < (6)The shadow fraction s is then: s = max (cid:40) π − cos − d + d √ − d π , − J * (cid:41) (7) ˆN · ˆN ˆ L can be computed in the vertex shader forspeed. View ray casting can result in artifacts due to in-accurate ray-surface intersections. To combat this,we developed a novel hybrid ray casting approach(Figure 7), where we start from a cubesphere geom-etry surface, which is a spherical mesh that entirelycontains the lunar terrain at the maximum terrainheight, and first displace vertices inward in the vertexshader with the same height field used to determineray visibility. Then rays are shot from the displacedvertices p in the fragment shader.7tarting from pre-displaced vertices has a coupleof advantages. The biggest is the decreased numberof steps required, because rays are already startingclose to their final intersection points.As Figure 8 shows, ray casting without pre-displacement requires up to 100 times the numberof steps to reach an equivalent level of accuracy andcan still result in artifacts at the object silhouette(Figure 9). Mesh resolution does not need to be veryhigh (we use a cube sphere with up to subdi-visions per cube face depending on the lod), as thepurpose of displacement is simply to bootstrap raystarting points to positions near the actual surface. Height fieldDisplaced meshView rays
Figure 7: Hybrid ray casting scheme. Increased accu-racy is obtained by shooting rays from vertices pre-displaced with the height field used to determine rayvisibility.Aliasing artifacts occur due to point sampling, andalso because of incorrect automatic mip level deduc-tion caused by discontinuous sampling of textures.We were able to reduce temporal flickering by reduc-ing sample rates (view ray step size) at highly oblique(a) (b) (c)
Figure 8: View ray step counts with vertex pre-displacement (a), and (b) without. Altitudes areshown in (c). Without pre-displacement, nearly 50–100 times more steps are required to trace the planetlimb and low altitudes. (a) (b)Figure 9: With vertex pre-displacement (a) and with-out (b). Ray casting without vertex pre-displacementresults in ugly artifacts and missed terrain, problemswhich are gone in our hybrid method.angles where ˆN · ˆV → , in a form of low-pass filtering.We mitigated the second issue by explicitly calculat-ing texture gradients using the method of tracing raydifferentials [13]. celestia.Sci [25] is a real-time, three-dimensional, in-teractive space simulation software that can renderobjects at a large range of scales, from Solar Systemobjects to deep space and galaxies [25]. The softwareis based on Celestia [15]. Celestia has been used byNASA [18] and ESA [10] due to its visualization ac-curacy and extensive astronomical database built onpeer-reviewed scientific data.celestia.Sci provides a scene management andframebuffer display framework, support for readingsolar system and spacecraft ephemerides, and script-ing and interactive features. This work adds shadersand new height and albedo texture types that canoptionally be used in an addon (a directory defininga new or modified celestial body). These new texturetypes consist of tiles organized in cube face directo-ries and their use in an addon automatically activatesthe new shaders. This design adds new features tocelestia.Sci while modifying existing code as little aspossible. Figure 10 outlines how our work integrateswith celestia.Sci.8 Results
Here we present results comparing our simulations toactual Moon imagery.
We compared south pole LRO WAC imagery withour simulations. As our renderer does not producegeoreferenced images, we instead render the WACimages unshaded in celestia.Sci by draping the im-ages onto the same Moon terrain mesh used by ourray casting engine (we use double the mesh frequencyto compensate for the comparatively low resolution).This guarantees that our view of the WAC imagesexactly register with our simulations. Moreover aHapke BRDF [23] and earthshine (indirect light fromthe Earth) was used to closely simulate lunar surfacelighting.WAC images are radiance while our simulationsproduce irradiance values; fortunately as these onlydiffer by a constant multiplicative factor we perform arelative comparison of normalized images by dividingall pixel intensities by their maximum.An example set of results is shown in Figure 11.For this case, the standard deviation of the error σ ≈ . , and 97% of pixels are within σ . Some ofthe ‘errors’ are in fact due to a diffuse glow surround-ing bright pixels in the WAC images that bleeds intoadjacent shadows, brightening them. This is a knownartifact called ghosting that is caused by stray light in This WorkExisting celestia.Sci Framework
Call
NewModifiedModified
Render LoopResource Manager16-bit Texture LoaderFramebuffer ObjectSimulation EnginePlanetary EphemeridesGUI
Cube SphereCube Virtual TextureData StructureMaximum Mipmap ShadersRaycasting ShadersGeometry Patch GeneratorHeight Texture (Tiles)Normal Texture (Tiles)Albedo Texture (Tiles)
Figure 10: Our work within celestia.Sci the camera [16] and which we do not try to simulate.To compensate, we subtract empirically determineddark pixel values from the normalized WAC imagesbefore computing errors. For Figure 11 we subtract adark value of 7/255; additional examples in Figure 15subtract 5/255. We attribute other small errors nearshadow boundaries mainly to artifacts from using amaximum mipmap (see Section 3.2) that have notbeen completely mitigated by the technique describedin Section 3.3.
We compare HDTV frames from the SELENE mis-sion [33] with our simulation. The frames [14] werecaptured during SELENE’s final collision trajectorynear the lunar south pole and are thus a useful ana-logue for images captured during a polar landing.Following Honda et al. [12], we color-calibrate theraw HDTV frames except we omit dark offset sub-traction as no reliable values could be found that didnot create an unnatural color cast. We also appliedirradiance scale factors C rr from Honda et al. as val-ues normalized to [0 , and finally applied the samesRGB transfer (gamma) function used in celestia.Scifor a consistent comparison. The same Hapke BRDFused in Section 6.1 was used here. Results are shownin Figure 12. The distinctive mottled appearance ofshadows and pale pinkish brown color of the terrainis closely matched by our render. Our hardware is a laptop with Intel Core i7-4900MQCPU, 16GB RAM, and NVIDIA GTX 970M graph-ics with 6GB VRAM (we disabled Scalable LinkInterface (SLI)). We ran our performance tests onthe laptop’s SATA III SSD rated at 6Gb/s and the1920 × a) (b) (c) Figure 11: LRO WAC polar imagery (a), Our render (b), Error (c)mipmaps have N in the range of to pixels.A lambertian BRDF with a single light source is usedsince we only want to show the performance of softshadow calculation and not accurate surface shading.Figure 13a shows that mean frame times are . ms( Hz) for N (cid:48) = { , , } (the default for most casesin this work) and . ms ( Hz) for N (cid:48) = { , , , } which is only a . % time penalty for a % increasein step count. By contrast, frame time is . ms foruniform stepping (fixed 100 step count, step size ∆ t =0 . ), indicating that the speedup of our methodis up to %. Another test on a slower laptop withIntel Core i5 CPU, 8 GB RAM, and integrated IntelIris 540 GPU driving an internal 1920 × > Hz. See Table 2 for a comparison of framerates on each hardware at different display sizes (e.g., for machine learning).Figure 13d shows the proportion of total frametime taken up by shadow ray tracing.Only visible height texture tiles are paged in syn-chronously on demand using a virtual texture schemeand stitched into a single texture to prevent visibleTable 2: SELENE scene frame rates (Hz) for differ-ent hardware and window sizes. 1920 × ismore suitable for training machine learning systems.GPU / Dimensions 1920 × × We presented a simple, novel method of efficientlyrendering ray cast soft shadows on curved terrainby using maximum mipmaps computed in real-timeand dynamic programming to find a global minimumshadow cost in constant runtime complexity. Addi-tionally, we demonstrated a method of reducing viewray computation times using vertex pre-displacementto bootstrap ray starting positions. Combining thesetwo methods, our ray casting engine runs in real-time with more than 200% speed up over uniform raystepping with comparable image quality and withouthardware ray tracing acceleration. The ray castingengine is integrated into celestia.Sci, a general, inter-active space simulation software. We demonstratedthe scalability of our engine by generating lunar im-agery across a wide range of distance scales, and accu-rately reproduced real lunar mission imagery contain-ing heavily shadowed scenes with long penumbra. Ar-tifacts due to approximating maximum delta heightsalong ray trajectories with a maximum mipmap were10a)(b) Figure 12: SELENE imagery (a), © JAXA/NHK, Our render (b).11
Uniform Step N ={5,5,5} N ={5,5,5,3} TotalShadow Rays F r a m e T i m e s ( m s ) (a)(b) (c)(d)Figure 13: Left: Frame times for SELENE scene near lunar surface. Histogram (inset) indicates most frametimes are clustered around the mean. Compared to uniform stepping (a), our method (b) is up to %faster with 185 Hz frame rates. The effect on frame times of varying N (cid:48) is also shown. Right: Frame timesfor a descent scene with >100 km vertical travel (c), (d).largely mitigated by starting at finer mip levels andconcatenating several runs of the algorithm to coverlong shadows. Image accuracy and system perfor-mance remain quite acceptable given that we are run-ning on previous-generation, commodity computinghardware.We thank Andrew Tribick for bug reports and feed-back. This work was supported in part by fundingfrom the Korea Ministry of Science and ICT throughthe “Development of Pathfinder Lunar Orbiter andKey Technologies for the Second Stage Lunar Explo-ration” project. Appendix
Definition (Deriving T ) . Letting R be the ob-ject space position vector of the shadow ray start-ing point, R end the position vector corresponding to R projected (spherical distortion) onto the textureplane at z = 1 and offset by ∆ M along the shadowray, we have: R = p + t ˆL (8) R end = RR z + 2 < ˆL xy , > (cid:107) < ˆL xy , > (cid:107) ∆ M (9)(The factor of is due to object space conversion) T can be found by deriving the intersection point between the shadow ray R + ˆL T and R end [11] (Fig-ure 14): a = ˆL b = R end c = − R T = ( c × b ) · ( a × b ) (cid:107) a × b (cid:107) (10) Δ M t z=1 R R + L Δ T - c ˆ R end ba Sun p textureplane t e rr a i n Figure 14: Deriving T . The two lines (position vec-tors) b and − c are known to intersect at the origin,so we can use Equation 10 for intersecting two linesin 3D to find the length a = T .12 .00.10.20.00.10.2 Figure 15: Additional LRO WAC comparisons. Left column: reference images, center: our renders, right:absolute error. N (cid:48) = 5 , , , was used due to long shadow lengths.(a) (b) (b)(a)Figure 16: Apollo 17 photograph AS17-152-23311, color-matched (a) vs our render (b).13a)(b)(c)Figure 17: Valles Marineris on Mars (a), (b) rendered using Mars Orbiter Laser Altimeter463m/px height map [19]. Sputnik Planitia on Pluto (c) rendered with 300m/px height mapfrom New Horizons data [24]. 14 eferences [1] John Amanatides and Andrew Woo. A FastVoxel Traversal Algorithm for Ray Tracing. In EG1987 Proceedings (Technical Papers) , 1987.[2] Lionel Baboud, Elmar Eisemann, and Hans Pe-ter Seidel. Precomputed safety shapes for effi-cient and accurate height-field rendering.
IEEETransactions on Visualization and ComputerGraphics , 18(11):1811–1823, 2012.[3] Richard Bellman. The Theory of Dynamic Pro-gramming.
Bull. Am. Math. Soc. , 60(6):503–515,1954.[4] Dimitri P. Bertsekas.
Dynamic Programmingand Optimal Control, Volume I . Athena Scien-tific, Nashua, NH, 3rd edition, 2005.[5] Edmund Burke and Graham Kendall.
SearchMethodologies: Introductory Tutorials in Op-timization and Decision Support Techniques .Springer, New York, 2nd edition, 2014.[6] Chun-Fa Chang, Bo-Quan Lin, Ying-ChiehChen, and Yung-Feng Chiu. Real-time softshadow for displacement mapped surfaces. In , pages1254–1257, New York, June 2009.[7] Robert L. Cook, Thomas Porter, Loren Carpen-ter, Robert L. Cook, Thomas Porter, and LorenCarpenter. Distributed ray tracing. In
Proc. 11thAnnu. Conf. Comput. Graph. Interact. Tech.- SIGGRAPH ’84 , volume 18, pages 137–145,New York, New York, USA, 1984. ACM Press.[8] Kathryn A. Dowsland. Classical techniques. InEdmund Burke and Graham Kendall, editors,
Search Methodol. Introd. Tutorials Optim. Decis.Support Tech. , chapter 2, pages 19–65. Springer,New York, 2nd edition, 2014.[9] Michał Drobot. Quadtree Displacement Map-ping with Height Blending. In Wolfgang En-gel, editor,
GPU Pro 360 , chapter 1, pages 1–32.Taylor & Francis, Boca Raton, FL, 2018. [10] ESA. Closing in on the Red Planet: Mars Ex-press Orbit Lowered, 2004.[11] F.S. Hill, Jr. The Pleasures of "Perp Dot" Prod-ucts. In Paul S. Heckbert, editor,
Graph. GemsIV , chapter II.5, pages 138–148. Academic Press,San Diego, CA, 1994.[12] Rie Honda, Junichi Yamazaki, Seiji Mitsuhashii,and Junichi Tachino. Calibration of Imagesby High Definition Television System OnboardKaguya (SELENE). In , 2010.[13] Homan Igehy. Tracing Ray Differentials. In
Proc. 26th Annu. Conf. Comput. Graph. Inter-act. Tech. , SIGGRAPH ’99, pages 179–186, LosAngeles, August 1999. ACM.[14] JAXA. SELenological and ENgineering Ex-plorer; SELENE Data Archive.[15] Chris Laurel, Clint Weisbrod, Fridger Schrempp,Bob Ippolito, Christophe Teyssier, Hank Ram-sey, Grant Hutchison, Pat Suwalski, Toti, Da-woon Jung, Vincent Giangiulio, and AndrewTribick. Celestia, 2001.[16] P. Mahanti, D. C. Humm, M. S. Robinson, A. K.Boyd, R. Stelling, H. Sato, B. W. Denevi, S. E.Braden, E. Bowman-Cisneros, S. M. Brylow, andM. Tschimmel. Inflight Calibration of the Lu-nar Reconnaissance Orbiter Camera Wide An-gle Camera.
Space Sci. Rev. , 200(1-4):393–430,April 2016.[17] F. Kenton Musgrave. Grid Tracing: Fast RayTracing for Height Fields. Technical report, YaleUniversity, 1990.[18] NASA. SCaN Network Demonstration, 2010.[19] Gregory A. Neumann, F. G. Lemoine, D. E.Smith, and M. T. Zuber. The Mars OrbiterLaser Altimeter Archive: Final Precision Exper-iment Data Record Release and Status of Ra-diometry. In ,League City, Texas, 2003.1520] Derek Nowrouzezahrai and John Snyder. Fastglobal illumination on dynamic height fields.
Comput. Graph. Forum , 28(4):1131–1139, 2009.[21] Manuel M. Oliveira and Fabio Policarpo. An Ef-ficient Representation for Surface Details. Tech-nical Report RP-351, Universidade Federal doRio Grande do Sul, 2005.[22] Fabio Policarpo and Manuel M. Oliveira. Re-laxed Cone Stepping for Relief Mapping. In Hu-bert Nguyen, editor,
GPU Gems 3 , pages 409–428. Addison-Wesley, 2007.[23] H. Sato, M. S. Robinson, B. Hapke, B. W.Denevi, and A. K. Boyd. Resolved Hapke Pa-rameter Maps of the Moon.
J. Geophys. Res.Planets , 119:1775–1805, August 2014.[24] Paul Michael Schenk, Ross A Beyer, William B.McKinnon, Jeffrey M. Moore, John R. Spencer,Oliver L. White, Kelsi Singer, Francis Nimmo,Carver Thomason, Tod R. Lauer, Stuart Rob-bins, Orkan M. Umurhan, William M. Grundy,S. Alan Stern, Harold A. Weaver, Leslie A.Young, K. Ennico Smith, Cathy Olkin, and theNew Horizons Geology and Geophysics Inves-tigation Team. Basins, Fractures and Volca-noes: Global Cartography and Topography ofPluto from New Horizons.
Icarus , 314:400–433,November 2018.[25] Fridger Schrempp. Welcome: Aims and Statusof celestia.Sci, 2013.[26] John Snyder and Derek Nowrouzezahrai. FastSoft Self-Shadowing on Dynamic Height Fields.In Steve Marschner and Michael Wimmer, ed-itors,
EGSR ’08 Proc. Ninet. EurographicsConf. Render. , pages 1275–1283, Aire-la-Ville,Switzerland, 2008. Eurographics Association.[27] Emerson J. Speyerer and Mark S. Robinson.Persistently Illuminated Regions at the LunarPoles: Ideal Sites for Future Exploration.
Icarus ,222(1):122–136, January 2013.[28] László Szirmay-Kalos and Tamás Umenhoffer.Displacement Mapping on the GPU — State of the Art.
Comput. Graph. Forum , 27(6):1567–1592, September 2008.[29] Márton Támas and Viktor Heisenberger. Prac-tical Screen Space Soft Shadows. In WolfgangEngel, editor,
GPU Pro 6 , chapter 4, pages 297–312. CRC Press, Boca Raton, FL, 2016.[30] Natalya Tatarchuk. Dynamic Parallax OcclusionMapping with Approximate Soft Shadows. In
Proc. 2006 Symp. Interact. 3D Graph. games ,I3D ’06, pages 63–69, Redwood City, CA, 2006.[31] Art Tevs, Ivo Ihrke, and Hans-Peter Seidel.Maximum Mipmaps for Fast, Accurate, andScalable Dynamic Height Field Rendering. In
Proc. 2008 Symp. Interact. 3D Graph. games ,I3D ’08, pages 183–190, New York, 2008. ACMPress.[32] Ville Timonen and Jan Westerholm. Scalableheight field self-shadowing.
Comput. Graph. Fo-rum , 29(2):723–731, 2010.[33] Junichi Yamazaki, Seiji Mitsuhashi, MasahitoYamauchi, Junichi Tachino, Rie Honda, Mo-tomaro Shirao, Kazuo Tanimoto, HiroyukiTanaka, Nobuaki Harajima, Asako Omori,Satoshi Yahagi, Shigehiro Kanayama, YuichiIijima, and Hisashi Ohtake. High-DefinitionTelevision System Onboard Lunar ExplorerKaguya (SELENE) and Imaging of the Moonand the Earth.