Oliver Klehm
email: oklehm@mpi-inf.mpg.de

Ivo Ihrke
email: ihrke@mpi-inf.mpg.de

Hans-Peter Seidel
email: hpseidel@mpi-inf.mpg.de

Elmar Eisemann
email: e.eisemann@tudelft.nl

Volume Stylizer: Tomography-based Volume Painting

Keywords: CR Categories: I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism-Color, shading, shadowing, and texture artist control, optimization, participating media

Figure 1: Volume stylization of an environmentally lit cloud hovering over a city. On the left, we show the original cloud model, whereas on the right, the volume was stylized to increase the atmospheric tension of the scene. The modifications include the red color cast of the cloud, an increase in contrast, and the addition of a logo. Our technique optimizes for the volume properties emission and albedo from a handful of user defined images. Rendering the new volume reveals details following the user input such as the clearly visible logo.

Introduction

Volume rendering has long been established as a tool to enrich the appearance of virtual scenes and enabling volumetric phenomena [START_REF] Kajiya | Ray Tracing Volume Densities[END_REF][START_REF] Max | Optical Models for Direct Volume Rendering[END_REF][START_REF] Jarosz | Progressive Photon Beams[END_REF]Novák et al. 2012]. Hereby, an atmospheric touch is added to otherwise "sterile" synthetic scenes. Voxels are often used to define volumes by storing physical properties such as extinction, absorption, and scattering behavior. These volumes are typically generated by means of simulation [START_REF] Stam | Stable Fluids[END_REF]] or noise functions [START_REF] Perlin | Hypertexture[END_REF]]. While effective methods for shape control [START_REF] Treuille | Keyframe Control of Smoke Simulations[END_REF][START_REF] Mcnamara | Fluid Control using the Adjoint Method[END_REF]] exist, current solutions do not allow for modifying the appearance of such volumes, especially under complex illumination conditions.

We provide a solution that facilitates appearance control. Starting with a lit volume under environmental lighting, the user modifies (or redefines) its appearance for certain viewpoints using familiar image editing operations. Hereby, we eliminate the need on the user's side to estimate the influence of changes under complex illumination conditions. To achieve plausible results from all possible viewpoints, our approach modifies physically-based volume properties (albedo, emission, or -under special conditions-extinction) to match given appearances under a standard rendering model. While we do restrict ourselves to single-scattering, many important cues are captured and images produced with this model are convincing. Further, our solution is fast enough for interactive editing sessions with update times in the order of seconds.

Precisely, our contributions are:

• We identify properties in the radiative-transport equation and derive conditions that allow us to enable stylization via a fast linear optimization (Sec. 4.1) • We provide an efficient implementation to keep execution time and memory cost practical (Sec. 5) • We show several use cases to illustrate our system (Sec. 6).

Related Work

It is very common for artists to start with a physical simulation, tweaking it to match a desired appearance while maintaining visual plausibility. Light source editing [START_REF] Schoeneman | Painting with Light[END_REF] or modifications to the light transport of spot lights [START_REF] Kerr | Bendy-Lights: Artistic Control of Direct Illumination by Curving Light Rays[END_REF]], indirect light [START_REF] Obert | iCheat: A Representation for Artistic Control of Indirect Cinematic Lighting[END_REF], or shadows [START_REF] Obert | Visibility Editing For All-Frequency Shadow Design[END_REF] have been proposed, but did not focus on participating media. Re-cently, first steps were taken to modify more complex effects like sub-surface scattering [START_REF] Song | SubEdit: a Representation for Editing Measured Heterogeneous Subsurface Scattering[END_REF]]. While the results look impressive, the scattering is relatively strong and objects need to be rather opaque. Light beams [START_REF] Nowrouzezahrai | A Programmable System for Artistic Volumetric Lighting[END_REF]] also include volumetric effects, but the focus is on light modification, not on changing the properties of the medium. Further, both methods require special rendering techniques, whereas we target the application of redefining properties of a volume relevant to the radiative transport equation [START_REF] Chandrasekar | Radiative Transfer[END_REF]] to ensure a plausible look for any view while keeping the rendering method unchanged.

Generating convincing volume data is known to be hard. This fact gave rise to methods that aimed at capturing properties of natural phenomena. These include flames [START_REF] Ihrke | Image-Based Tomographic Reconstruction of Flames[END_REF], smoke [START_REF] Hawkins | Acquisition of Time-Varying Participating Media[END_REF][START_REF] Ihrke | Adaptive Optical Tomography[END_REF]] and refractive elements [START_REF] Ihrke | Reconstructing the Geometry of Flowing Water[END_REF][START_REF] Atcheson | Time-resolved 3d capture of non-stationary gas flows[END_REF]. Many of them rely on tomographic reconstruction [START_REF] Ihrke | Image-Based Tomographic Reconstruction of Flames[END_REF][START_REF] Ihrke | Adaptive Optical Tomography[END_REF][START_REF] Atcheson | Time-resolved 3d capture of non-stationary gas flows[END_REF]] and the ability to convincingly recover volumetric descriptions from a sparse number of views (8 to 16) hints at the possibility of defining volume properties in an image-based fashion. However, the image formation models employed in these techniques are usually simple and restricted to a single phenomenon such as emission [START_REF] Ihrke | Image-Based Tomographic Reconstruction of Flames[END_REF][START_REF] Ihrke | Adaptive Optical Tomography[END_REF], or refraction [START_REF] Atcheson | Time-resolved 3d capture of non-stationary gas flows[END_REF], only Lintu et al. [2007] attempt recovery of both, emission and absorption, for a nebula from astronomical observations. In contrast, we aim at artistic control and modify existing volumes to achieve a certain appearance for a low number of views under complex lighting. This goal shares some characteristics with recent developments in fabrication, where objects are physically built to have a certain interaction with light in the real world. Shadows [START_REF] Pauly | Shadow Art[END_REF][START_REF] Baran | Manufacturing Layered Attenuators for Multiple Prescribed Shadow Images[END_REF], caustics [START_REF] Papas | Goal-based Caustics[END_REF][START_REF] Papas | The magic lens: Refractive steganography[END_REF], and surface reflectance [START_REF] Weyrich | Fabricating Microgeometry for Custom Surface Reflectance[END_REF]] have been investigated. Others [START_REF] Holroyd | Computing and fabricating multilayer models[END_REF][START_REF] Wetzstein | Layered 3D: Tomographic Image Synthesis for Attenuation-based Light Field and High Dynamic Range Displays[END_REF]] made it possible to show specific views when uniformly illuminating from the back and observing under particular viewing angles. The latter techniques lead to light-field display technology [START_REF] Lanman | Polarization Fields: Dynamic Light Field Display using Multi-Layer LCDs[END_REF][START_REF] Wetzstein | Tensor Displays: Compressive Light Field Synthesis using Multilayer Displays with Directional Backlighting[END_REF], or even 6D displays [START_REF] Fuchs | Towards Passive 6D Reflectance Field Displays[END_REF]]. Nonetheless, these methods usually consider only a very small number of three to five layered light modulating planes, which are viewed from a well-defined viewing zone. Also, only one particular effect like absorption [START_REF] Holroyd | Computing and fabricating multilayer models[END_REF][START_REF] Wetzstein | Layered 3D: Tomographic Image Synthesis for Attenuation-based Light Field and High Dynamic Range Displays[END_REF] or change of polarization [START_REF] Lanman | Polarization Fields: Dynamic Light Field Display using Multi-Layer LCDs[END_REF]] is considered. We enable control over appearance and achieve a full surround view under complex illumination with volumes exhibiting full resolution in all three spatial dimensions. Since we work in a rendering context, we can also use more information than for physical object manufacture and display technology. In particular, refracted ray paths are possible, where previous work is restricted to straight rays passing through a volume.

Background

In this section, we discuss the necessary background for our approach. In particular, we review the volumetric rendering equation and show how it, under specific conditions, can be inverted. The inversion leads to a standard tomography problem which is being solved by inverting a large scale linear system.

We start by reviewing volumetric rendering, which simulates how light interacts with participating media. We also introduce the main properties that influence the appearance. This background will be useful to understand how we offer control over appearance.

Rendering complex light transport in participating media is done by solving the radiative transport equation [START_REF] Chandrasekar | Radiative Transfer[END_REF]]. It applies the physical-based cross section volume model where the volume is described by infinitely small particles of random orientation. The particles can absorb, scatter, and emit light. For the moment, we do not consider refraction, which will be discussed in Sec. 5.1. We summarize all symbols in Tab. 1.

Extinction volume σ e L s Sampled Environment Map Camera x c T r (x c ,x t ) x t Q o (x t ) L s'' L s' L s''' L s'''' L i /L d (x c ,ω i ) k

Symbol

Volume properties σt Extinction coefficient (particle per , in range

Table 1: Symbols used for volume rendering

The light incident at a point in the scene x (e.g., the camera position) from direction ωi consists of two terms:

Li(x, ωi) = Tr(x, xs) Lo(xs, ωi) + c Tr(x, xt) Qo(xt, ωi) dt. (1)
The first summand is the reflected radiance at the first visible surface or emitted radiance from the background, attenuated by the volume Tr(x, xs). The second is the medium radiance Qo(x, ωo), which is emitted directly or out-scattered at each point xt along the light path c and again attenuated by the volume. Consequently, we split Qo(x, ωo) into emission and scattering:

Qo(x, ωo) = Qemit(x, ωo) + Qscat(x, ωo). (2)
Emission is the simpler part; each particle emits a fixed radiance independent of direction and regardless of other properties of the

v 1 image 1 o k o k-1 o k+1 image 2 o n-1 o n o n+1
(o k-1 , o k , o k+1 ) and in image2 (o n-1 , o n , o n+1 ) volume: Qemit(x, ωo) = σt(x) Lemit(x). (3)
Simulating scattering requires integrating over the entire unit sphere surrounding the point to gather incoming light, which is then modulated by the phase function (the volumetric equivalent to BRDFs for surfaces; while not required by our method, we use an isotropic phase function for convenience).

Qscat(x, ωo) = σs(x) Lscat(x, ωo) = σt(x) ρ(x) Ω f (x, (ωi • ωo)) Li(x, ωi) dωi.
It is important to note that only particles that scatter light σs are considered. The relation between these and all particles is defined by the albedo ρ. To compute the incoming light in Lscat(x, ωo), one has to apply Eq. 1 recursively, which makes volume rendering hard.

For simplification, multiple volume scattering is often ignored, i. e., Li(x, ωi) := Tr(x, xs) Lo(xs, ωi); only light originating outside the volume (environmental, direct light sources. . . ) is considered, see Fig. 2. In the following, we also apply this approximation.

Volume Reconstruction

To control appearance, we want to invert the volume-rendering process, upon receiving a single or multiple user-defined images defining the desired target views. We treat images as a collection of N constraint pixels. Given the corresponding camera view, a pixel with index k ∈ [1; N ] corresponds to a ray (origin x k and direction ω k ). We denote its value L k i (x k , ω k ) and, if applicable, the position of the first hit surface x k s , which otherwise is the background at ∞. We want to modify properties of a volume V , such that the pixel constraint is matched, namely Li(x k , ω k ) = L k i (x k , ω k ). For example, one could modify the emission field of V to match a given appearance.

Volume Reconstruction

A single pixel constraint can influence many voxels and, inversely, two pixel constraints might imply changes on one and the same voxel. We illustrate this situation in Fig. 3. Consequently, a perfect solution might not always be possible. Instead, we seek to find a coefficient vector a := (. . . , a i , . . .) such that a linear combination i a i υ i (x) of known basis functions υ i defines a property of the volume such that the constraint pixels are matched best.

From a mathematical point of view, the basis functions could have global support. However, in practice, using a basis with spatially local support speeds up the reconstruction. One convenient choice for a basis are box functions associated to the volume's voxels, which corresponds to nearest-neighbor sampling of a 3D texture. Using triangle functions allows us to consider linearly interpolated solutions as well.

Next, we will show that the best match for emission, scattering, and extinction is the solution of a linear system.

o = Wa, (4)
where a are the basis coefficients to be computed. The observations o := (. . . , o k , . . .) T involve the constraint pixel values, and the matrix W := (. . . , w k , . . .) T can be derived from the volume rendering equation.

Property Reconstruction

The volume parameter reconstructions imply that Eq. 1 needs to be linearized. In the following, we derive the entries of matrix W necessary to estimate specific volume properties. The derivation is carried out for a single constraint pixel, i. e., for one row of matrix W.

Emission reconstruction starts with Eq. 1 in conjunction with Eqs. 2 and 3. For a constraint pixel k, we obtain:

L k i (x k , ω k ) = Tr(x k , x k s ) Lo(x k s , ω k ) + s 0 Tr(x k , xt) (σt(xt) Lemit(xt) + Qscat(xt, ω k )) dt.
Assuming single scattering, we can split the integral into emitted and scattered light and unify all known values in o k :

o k := L k i (x k , ω k ) -Tr(x k , x k s ) Lo(x k s , ω k )- s 0 Tr(x k , xt) Qscat(xt, ω k ) dt = s 0 Tr(x k , xt) σt(xt) Lemit(xt) dt. (5)
We represent Lemit(xt) by a linear combination of basis functions, whose integral can be computed:

o k = s 0 Tr(x k , xt) σt(xt) i L i emit υ i (xt) dt = i L i emit s 0 Tr(x k , xt) σt(xt) υ i (xt) dt =: ( Lemit • w k ). (6)
The coefficients of the above equation are defined by an integral corresponding to one ray passing through the volume. Combining all equations defined by the constraint pixels, we obtain a linear system.

Albedo is reconstructed by establishing a similar linearization as for emission. Combining volume rendering (Eq. 1) and the constraints from a pixel k, we obtain -similar to Eq. 5:

o k := L k i (x k , ω k ) -Tr(x k , x k s ) Lo(x k s , ω k )- s 0 Tr(x k , xt) Qemit(xt) dt = s 0 Tr(x k , xt) σt(xt) Lscat(xt, ω k ) ρ(xt) dt.
With single scattering, Lscat(xt, ω k ) is independent of ρ, and we can solve for it. More precisely, representing the field of ρ with a linear combination of basis functions, we obtain:

o k = i ρ i s 0 Tr(x k , xt) σt(xt) Lscat(xt, ω k ) υ i (xt) dt =: ( ρ • w k ).
Again, the coefficients are defined via an integral that translates to a ray marching process involving the known properties of the volume.

Emission & Albedo can be jointly optimized because Lemit and ρ are linearly independent.

Extinction can only be reconstructed if we assume that the volume's outgoing radiance is constant, i. e., Qo(x, ω) = Qo(x ′ , ω ′ ) ∀x ′ , ω ′ . As extinction is mostly used to define the overall shape of the volume and usually the first property to be derived, this restriction is usually not too problematic. Starting with Eq. 1:

L k i (x k , ω k ) = e -s 0 σ t (x t )dt Lo(x k s , ω k )+ s 0 e -t 0 σ t (x t ′ )dt ′ Qo(xt, ω k ) dt.
We employ the linear combination of basis functions and can rewrite the first summand:

Lo(x k s , ω k ) e -s 0 σ i t υ i (x t )dt = Lo(x k s , ω k ) i e -σ i t s 0 υ i (x t )dt .
Concerning the second summand, we exploit the constant outgoing radiance and remove Qo from the integral. The remainder can be integrated, yielding 1e -s 0 σ t (x t )dt . Mathematically, the result can be shown by decomposing σt into a piecewise-constant approximation and splitting the outer integral accordingly [START_REF] Max | Optical Models for Direct Volume Rendering[END_REF]]. Then each integral can be solved and the result recombined. This proof is valid for any Riemann-integrable extinction function. The proof also follows logically; the above remainder is the probability that a ray from the camera passes through the volume without hitting a particle, so it is one minus the probability that a ray is stopped. Now, we can use a transformation similar to the first summand to obtain:

L k i (x k , ω k ) = Lo(x k s , ω k ) i e -σ i t s 0 υ i (x t )dt + Qo 1 - i e -σ i t s 0 υ i (x t )dt = Qo + (Lo(x k s , ω k ) -Qo) i e -σ i t s 0 υ i (x t )dt .
Finally, we apply a logarithm to obtain a linear system in σ i t :

o k := ln( L k i (x k , ω k ) -Qo Lo(x k s , ω k ) -Qo ) = ln i e -σ i t s 0 υ i (x t )dt = i ln(e -σ i t s 0 υ i (x t )dt ) = i -σ i t s 0 υ i (xt)dt .

Discussion of Reconstruction Schemes

In all three cases that we discussed previously, we start with the volume rendering integral that contains the volume property of interest. The target field is represented as a linear combination of basis functions, which allows us to move the coefficients out of the integral, and, hereby, to isolate the unknowns. The same process is applied in computed tomography [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF], which corresponds to our reconstruction of the extinction coefficient. However, different from medical computed tomography, we have to deal with inconsistent user input, i. e., for which there may be no solution that would satisfy o = Wa. Therefore, we opt for a solution in the least squares sense W T o = W T Wa, which means that the optimal solution minimizes the quadratic function ||Wa -o|| 2 .

Extinction describes the overall shape of a volume, but is hard to optimize for. Emission and albedo are sufficient to change the appearance of an existing volume. Extinction cannot simultaneously be estimated in combination with emission or albedo, since it involves the solution of a complex non-linear problem. However, it is possible to begin with the reconstruction of extinction if Qo is constant. Hence, one can solve this simplified problem in a first step.

Based on this result, one can then add scattering light and refine the appearance by estimating albedo and emission while lifting the constraint that Qo is constant.

Implementation

In order to ensure a quick feedback to the user, we map our optimization to the GPU via compute shaders in OpenGL. However, this mapping is not direct, as special care is needed regarding memory and multi-threading management.

Further, the matrix W is large (total count of constraints pixels × number of voxels), the linear system is ill-conditioned, and, finally, we seek a physically plausible, i. e., a non-negative solution.

Fortunately, W is sparse because each row is derived by a ray passing through the volume, intersecting only a low number of basis functions υ. Still it would be too large to keep in memory. Instead, we implicitly solve the system by performing a conjugate gradient minimization of the quadratic function ||Wa -o|| 2 . All necessary steps of the conjugate gradient method are carried out on the GPU and involve 3D textures to represent the vectors. Operations on these vectors are implemented as shaders.

We employ the conjugate gradient method [START_REF] Shewchuk | An introduction to the conjugate gradient method without the agonizing pain[END_REF]] due to its fast convergence, but, for clarity, we illustrate the required operations by describing a standard gradient descent which yields the same minimum. The main ingredients are the computation of the gradient W T (Wao), an update of our current solution a by adding a scaled version of it and iteration. Performing the update is straightforward, but the computation of the gradient is not.

To understand W T (Wa-o), we examine its elements step by step.

The matrix W encodes volume rendering (e.g., for Eq. ( 6)):

Wa = s 0 Tr(x k , xt) σt(xt) i a i υ i (xt)
. Hence, each entry in Wa can be determined via rendering using a volume defined by a. Next, computing c := (. . . , c k , . . .) T := Wao is straightforward, but applying W T to c is a data scattering operation. Each value c k is associated to a constraint pixel k. It needs to be scattered to those voxels that are traversed by a ray marching procedure for the ray associated with pixel k, weighted according to the voxel's influence on k. In other words, we have to perform a ray marching that is very similar to the case for computing the product with W.

In fact, both operations are so similar, that almost the same code can be reused. During rendering, which implements the multiplication by W, we use ray-marching to sum the voxel contributions along the ray following outRadiance += weight • texRead(a, rayPos). The variable weight includes the transmittance between the current position rayPos and the ray origin due to the volumes' extinction, as well as the value of the basis function at the current position. To implement the multiplication by W T , this single line is exchanged by texWriteAdd(a, rayPos, weight • c k ). Reusing the code is also beneficial as the weight values match up perfectly, and, consequently, the implicitly constructed matrices W and W T agree with each other. The full pseudo code can be found in the supplementary material.

Optimizations are possible to speed up the rendering shaders: while multiplying with W is efficient and directly parallelized as usual when stepping through the volume per pixel, the situation is different for multiplication by W T . As we scatter data (texture writes are realized via shader image load store), synchronization issues may occur. Rays of neighboring pixels will likely write to the same voxel. To avoid the resulting stall, we use an interleaved pattern of 6 × 6 pixels. In each round, only one ray of these sub-windows is shot, decreasing the number of conflicts and speeding up the computation by ≈ 10%. Further, the use of 16bit instead of 32bit textures leads to a speedup of ≈ 25%.

Optimization Extensions

Higher precision is obtained when using accurate ray traversal [START_REF] Amanatides | A Fast Voxel Traversal Algorithm for Ray Tracing[END_REF] instead of marching. For several basis-function choices, e.g., nearest-neighbor or linear (which we use), an accurate integral can be computed. This applies for the accumulated transmittance value during volume traversal as well as the weight of the basis function itself.

Regularization is a standard way of stabilizing and controlling the optimization. Usually, additional quadratic terms, modeling prior knowledge about the solution space, are added to the quadratic error function to address numerical ill-conditioning. A Laplacian term La 2 can smooth the overall volume. a 2 minimizes the solution, while a -1 2 biases it towards one. Each of these regularizers is controlled by user-defined weights. Thus, when optimizing for albedo and emission, we can specify which element to favor and, e.g., minimize emission. In practice, we use weak weights (i.e., 2 -9 ), which proved sufficient for a stable solution.

Weights are often desirable to give different parts of an image more importance than others. They are also handy when creating transitions, e.g., when editing only a part of the volume. These perconstraint pixel weights can be easily integrated into the conjugate gradient method and need to be multiplied with the vector Wao.

Refraction occurs for non-constant refractive indices and rays are bent during the traversal. Our reconstruction scheme can handle arbitrary integration curves c of known geometry (cf. Eq. 1), which enables us to use more complex paths than a straight line. To compute the refractive ray paths in our volumetric setting, we resort to the Euler forward scheme of [START_REF] Ihrke | Eikonal Rendering: Efficient Light Transport in Refractive Objects[END_REF].

Visual Hulls can be used to limit the domain of solution [START_REF] Ihrke | Image-Based Tomographic Reconstruction of Flames[END_REF], which increases quality and performance. In case of emission/albedo optimization, we only consider voxels with non-zero extinction, as the other areas have no influence on the rendering. For estimating extinction itself, the input views of the user can be transformed into masks. In all cases, the user can specify masks (projections of the desired visual hull) as additional input to the optimization. Those masks need not to align with any of the input images, the parts of the volume that are not inside the visual hull are simply ignored during the optimization.

Results and Discussion

We used an Intel Xeon x5650 with NVIDIA 560Ti. We obtain fast execution times, even for medium sized voxel grids and high precision. The execution times for different numbers of ray marching steps, voxel resolution, and input images are given in Fig. 4. The process converges after three to five conjugate gradient steps. For 128 3 voxels, 256 ray marching steps, 524K constraint pixels (2 views), the result (Fig. 5) is computed in less than two seconds. With 2.1M constraint pixels (9 views) the computation takes seven seconds (Fig. 6). In practice, a good choice is twice the number of ray marching steps for the length of the diagonal as the number of voxels along an axis. In the following, we illustrate our method on three examples. It took a user roughly 2 min to produce the smoke bunny, 15 minutes for the cloud scene and 4 min for the refraction.

Smoke Bunny shows our approach with imprinted logos. These were blended with the original appearance in the user-provided views -a purely artistic choice. When optimizing albedo and emission, the result is a close match to the input (Fig. 5, top), while intermediate views look appealing and plausible (middle). One can also modify emission or albedo only (bottom). The result can no longer match the input. Emission can only lighten, albedo only darken the appearance (following physical constraints). Depending on the required realism and desired material, these could be adequate choices. Cloud Stylization illustrates the expressiveness of our system. The scene shows a cloud over a city and the goal is to achieve a "frightening" look when on one side of the cloud, a calm appearance on the other side. Modifying appearance in 10 views (640 × 360 pixel), lead to the result in Fig. 6.

The number of views is slightly increased for two reasons; to define a good appearance all around the object and to avoid visual inconsistencies. The system only optimizes for parts of the volume that are specified in at least a single view. Other parts are unchanged and may result in a soft edge after the optimization. Also, if insufficient constraints are used, stripes can become visible, as the color of a pixel constraint only defines the voxels along its corresponding ray. One could increase the weight of the Laplacian regularizer to avoid these artifacts, but we found that this unnecessarily destroys much of the intricate detail.

The view modifications were done by contrast enhancement in the according views in conjunction with a blended red mask. The resulting input images were inconsistent in parts, as each was independently designed. Yet, the least-square solution maintains the most faithful reconstruction.

The effect of rendering the stylized volume under multiple scattering is shown in Fig. 7. While multiple scattering acts like a blur, that reduces details, the reconstructed logo remains clearly visible.

Extinction and Refraction is shown in Fig. 8. Here, we illustrate the use of our extinction optimization to show the construction of complex shapes and also to illustrate the compatibility with bent rays due to continuous refraction. The user provided four in- put views (growth of a tree) and our reconstruction computed extinction coefficients (per wavelength) that attenuate the background light such as to match the provided images. In the example, the volume does not scatter or emit any light, i. e., Qo = 0. In this reconstruction process, the shape of the volume is implicitly defined by the volume of varying refractive index. The refraction indices were generated using a random process restricted to the inside of a glass sphere. We also used per-pixel weights to concentrate the importance on the main parts of the user images, which leads to a less constrained boundary and a more variable shape. Due to the refraction, the intermediate views appear as random colors. We imagine, such a combination of refraction and extinction could e.g. be used in a game, where the user has to find the correct viewpoint to see a certain image, in order to obtain hints for the solution of a puzzle.

Conclusions

We presented a novel approach to stylize volumes. We showed that solving for a specific volume appearance is possible by adjusting the desired appearance in a number of reference views. Our formulation as an optimization problem results in a large scale linear system. We make this solution practical by employing an efficient GPU-friendly implementation that avoids the explicit construction of the equation system. Several cases illustrated the variability of the method.

In the future, we plan to consider multiple scattering, the current problem being the cost of derivative computations and the demands on memory. Fully GPU-based out-of-core rendering approaches hint at a possible direction [START_REF] Crassin | GPU Pro. AK Peters, ch. X.3 Efficient Rendering of Highly Detailed Volumetric Scenes with GigaVoxels[END_REF]].

Figure 2 :

 2 Figure 2: Single scattering: Radiance is accumulated along the view ray. At each sample point xt emitted and incoming attenuated light is scattered towards the observer (Qo = Qemit + Qscat).

  [0; ∞], it holds: σt = σs + σa) σs Scattering coefficient (ratio of particles scattering incoming light) σa Absorption coefficient (ratio of particles absorbing incoming light) ρ Albedo, i. e., σs(x) = σt(x) ρ(x). Lemit Radiance emitted by the volume f Volumetric phase function Description x, xs Position, first surface hit by view ray or infinity ω, ω i/o , Ω Direction, incoming/outgoing, unit sphere of directions L i/o Radiance incoming/outgoing Qo Medium radiance outgoing from volume point Tr(xa, x b ) := e -c σ t (x t )dt -Transmittance or probability, that a photon travels along an arc-length parameterized curve c from xa to x b

Figure 3 :

 3 Figure 3: Volume from two viewpoints. Problem: Changing one voxel v1 influences several pixels in image1 (o k-1 , o k , o k+1 ) and in image2 (o n-1 , o n , o n+1 )

Figure 4 :

 4 Figure 4: Timings in ms for Bunny per iteration (three for convergence). Left: volume resolution vs. ray marching steps (relative to volume diagonal), more steps estimate the integrals more accurately. Right: volume resolution vs. # of views. In addition to the number of inputs, the performance of the optimization depends on volume resolution as well as the number of ray marching steps.

Figure 5 :

 5 Figure 5: Logos on a smoke bunny. Top: optimized emission and albedo and 8x difference (false color: blue to red); Middle: Intermediate views appear plausible. Bottom: left two images use only emission, the other two only albedo.

Figure 6 :

 6 Figure 6: Cloud stylization; original volume (left), user input (middle), result after optimization (right), except for the last row, which shows random intermediate views.

Figure 7 :

 7 Figure 7: Single (left) vs. multiple scattering (right). To account for the additional energy, we scaled the environmental light such that both images have similar brightness.

Figure 8 :

 8 Figure 8: Reconstruction of extinction in a refractive volume to absorb light of the background. Left, middle: results for user-defined views, rightmost column: intermediate views.

Acknowledgements

We would like to thank Bernhard Reinert, and colleagues at MPI for helpful discussions and the reviewers for their valuable feedback. This work was supported by the German Research Foundation (DFG) through the Emmy-Noether fellowship IH 114/1-1 and by the Intel Visual Computing Institute at Saarland University. The city model is from 3DRT.