aboutsummaryrefslogtreecommitdiff
path: root/articles/2021-09-26_noise_and_ray_marching.org
diff options
context:
space:
mode:
Diffstat (limited to 'articles/2021-09-26_noise_and_ray_marching.org')
-rw-r--r--articles/2021-09-26_noise_and_ray_marching.org276
1 files changed, 276 insertions, 0 deletions
diff --git a/articles/2021-09-26_noise_and_ray_marching.org b/articles/2021-09-26_noise_and_ray_marching.org
new file mode 100644
index 0000000..d3b231d
--- /dev/null
+++ b/articles/2021-09-26_noise_and_ray_marching.org
@@ -0,0 +1,276 @@
+* Noise and Ray Marching
+[[https://literatelb.org][LiterateLB's]] volumetric visualization functionality relies on a simple ray marching implementation
+to sample both the 3D textures produced by the simulation side of things and the signed distance
+functions that describe the obstacle geometry. While this produces surprisingly [[https://www.youtube.com/watch?v=n86GfhhL7sA][nice looking]]
+results in many cases, some artifacts of the visualization algorithm are visible depending on the
+viewport and sample values. Extending the ray marching code to utilize a noise function is
+one possibility of mitigating such issues that I want to explore in this article.
+
+While my [[https://www.youtube.com/watch?v=J2al5tV14M8][original foray]] into just in time visualization of Lattice Boltzmann based simulations
+was only an aftertought to [[https://tree.kummerlaender.eu/projects/symlbm_playground/][playing around]] with [[https://sympy.org][SymPy]] based code generation approaches I have
+since put some work into a more fully fledged code. The resulting [[https://literatelb.org][LiterateLB]] code combines
+symbolic generation of optimized CUDA kernels and functionality for just in time fluid flow
+visualization into a single /literate/ [[http://code.kummerlaender.eu/LiterateLB/tree/lbm.org][document]].
+
+For all fortunate users of the [[https://nixos.org][Nix]] package manager, tangling and building this from the [[https://orgmode.org][Org]]
+document is as easy as executing the following commands on a CUDA-enabled NixOS host.
+
+#+BEGIN_SRC sh
+git clone https://code.kummerlaender.eu/LiterateLB
+nix-build
+./result/bin/nozzle
+#+END_SRC
+
+** Image Synthesis
+The basic ingredient for producing volumetric images from CFD simulation data is to compute
+some scalar field of samples \(s : \mathbb{R}^3 \to \mathbb{R}_0^+\). Each sample \(s(x)\) can be assigned a color
+\(c(x)\) by some convenient color palette mapping scalar values to a tuple of red, green and blue
+components.
+
+[[https://literatelb.org/tangle/asset/palette/4wave_ROTB.png]]
+
+The task of producing an image then consists to sampling the color field along a ray assigned
+to a pixel by e.g. a simple pinhole camera projection. For this purpose a simple discrete
+approximation of the volume rendering equation with constant step size \(\Delta x \in \mathbb{R}^+\) already
+produces suprisingly good pictures. Specifically
+$$C(r) = \sum_{i=0}^N c(i \Delta x) \mu (i \Delta x) \prod_{j=0}^{i-1} \left(1 - \mu(j\Delta x)\right)$$
+is the color along ray \(r\) of length \(N\Delta x\) with local absorption values \(\mu(x)\). This
+local absorption value may be chosen seperately of the sampling function adding an
+additional tweaking point.
+
+#+BEGIN_EXPORT html
+<video style="width:100%" src="https://literatelb.org/media/nozzle.webm" controls="controls">
+</video>
+#+END_EXPORT
+
+The basic approach may also be extended arbitrarily, e.g. it is only the inclusion of a couple
+of phase functions away from being able [[https://tree.kummerlaender.eu/projects/firmament/][recover the color produced by light travelling through the participating media that is our atmosphere]].
+
+** The Problem
+There are many different possibilities for the choice of sampling function \(s(x)\) given the results of a
+fluid flow simulation. E.g. velocity and curl norms, the scalar product of ray direction and shear layer
+normals or vortex identifiers such as the Q criterion
+\[ Q = \|\Omega\|^2 - \|S\|^2 > 0 \text{ commonly thresholded to recover isosurfaces} \]
+that contrasts the local vorticity and strain rate norms. The strain rate tensor \(S\) is easily
+recovered from the non-equilibrium populations \(f^\text{neq}\) of the simulation lattice — and is in
+fact already used for the turbulence model. Similarly, the vorticity \(\Omega = \nabla \times u\) can be
+computed from the velocity field using a finite difference stencil.
+
+The problem w.r.t. rendering when thresholding sampling values to highlight structures in the flow
+becomes apparent in the following picture:
+
+#+BEGIN_EXPORT html
+<div class="flexcolumns">
+<div>
+<span>Q Criterion</span>
+<img src="https://static.kummerlaender.eu/media/q_criterion_default.png"/>
+</div>
+<div>
+<span>Curl Norm</span>
+<img src="https://static.kummerlaender.eu/media/curl_default.png"/>
+</div>
+</div>
+#+END_EXPORT
+
+While the exact same volume discretization was used for both visualizations, the slices are much
+less apparent for the curl norm samples due to the more gradual changes. In general the issue is
+most prominent for scalar fields with large gradients (specifically the sudden jumps that occur
+when restricting sampling to certain value ranges as is the case for the Q criterion).
+
+** Colors of Noise
+The reason for these artifacts is primarily choice of start offsets w.r.t. the traversed volume
+in addition the the step width. While this tends to become less noticable when decreasing said
+steps, this is not desirable from a performance perspective.
+
+What I settled on for LiterateLB's renderer are view-aligned slicing and random jittering to remove
+most visible artifacts. The choice of /randomness/ for jittering the ray origin is critical here as plain
+random numbers tend to produce a distracting static-like pattern. A common choice in practice is
+to use so called /blue noise/ instead. While both kinds of noise eliminate most slicing artifacts, the
+remaining patterns tend to be less noticeable for blue noise. Noise is called /blue/ if it contains only
+higher frequency components which makes it harder for the pattern recognizer that we call brain to
+find patterns where there should be none.
+
+The [[https://www.spiedigitallibrary.org/conference-proceedings-of-spie/1913/0000/Void-and-cluster-method-for-dither-array-generation/10.1117/12.152707.short?SSO=1][void-and-cluster algorithm]][fn:vac] provides a straight forward method for
+pre-computing tileable blue noise textures that can be reused during the actual visualization.
+Tileability is a desirable property for this as we otherwise would either need a noise texture
+large enough to cover the entire image or instead observe jumps at the boundary between
+the tiled texture.
+
+The first ingredient for /void-and-cluster/ is a =filteredPattern= function that applies a
+plain Gaussian filter with given $\sigma$ to a cyclic 2d array. Using cyclic wrapping during the
+application of this filter is what renders the generated texture tileable.
+
+#+BEGIN_SRC python
+def filteredPattern(pattern, sigma):
+ return gaussian_filter(pattern.astype(float), sigma=sigma, mode='wrap', truncate=np.max(pattern.shape))
+#+END_SRC
+
+This function will be used to compute the locations of the largest void and tightest
+cluster in a binary pattern (i.e. a 2D array of 0s and 1s). In this context a /void/ describes
+an area with only zeros and a /cluster/ describes an area with only ones.
+
+#+BEGIN_SRC python
+def largestVoidIndex(pattern, sigma):
+ return np.argmin(masked_array(filteredPattern(pattern, sigma), mask=pattern))
+#+END_SRC
+
+These two functions work by considering the given binary pattern as a float array that is blurred by
+the Gaussian filter. The blurred pattern gives an implicit ordering of the /voidness/ of each pixel, the
+minimum of which we can determine by a simple search. It is important to exclude the initial binary
+pattern here as void-and-cluster depends on finding the largest areas where no pixel is set.
+
+#+BEGIN_SRC python
+def tightestClusterIndex(pattern, sigma):
+ return np.argmax(masked_array(filteredPattern(pattern, sigma), mask=np.logical_not(pattern)))
+#+END_SRC
+
+Computing the tightest cluster works in the same way with the exception of searching the largest array
+element and masking by the inverted pattern.
+
+#+BEGIN_SRC python
+def initialPattern(shape, n_start, sigma):
+ initial_pattern = np.zeros(shape, dtype=np.bool)
+ initial_pattern.flat[0:n_start] = True
+ initial_pattern.flat = np.random.permutation(initial_pattern.flat)
+ cluster_idx, void_idx = -2, -1
+ while cluster_idx != void_idx:
+ cluster_idx = tightestClusterIndex(initial_pattern, sigma)
+ initial_pattern.flat[cluster_idx] = False
+ void_idx = largestVoidIndex(initial_pattern, sigma)
+ initial_pattern.flat[void_idx] = True
+ return initial_pattern
+#+END_SRC
+
+For the initial binary pattern we set =n_start= random locations to one and then repeatedly
+break up the largest void by setting its center to one. This is also done for the tightest cluster
+by setting its center to zero. We do this until the locations of the tightest cluster and largest
+void overlap.
+
+#+BEGIN_SRC python
+def blueNoise(shape, sigma):
+#+END_SRC
+
+The actual algorithm utilizes these three helper functions in four steps:
+1. Initial pattern generation
+ #+BEGIN_SRC python
+ n = np.prod(shape)
+ n_start = int(n / 10)
+
+ initial_pattern = initialPattern(shape, n_start, sigma)
+ noise = np.zeros(shape)
+ #+END_SRC
+3. Eliminiation of =n_start= tightest clusters
+ #+BEGIN_SRC python
+ pattern = np.copy(initial_pattern)
+ for rank in range(n_start,-1,-1):
+ cluster_idx = tightestClusterIndex(pattern, sigma)
+ pattern.flat[cluster_idx] = False
+ noise.flat[cluster_idx] = rank
+ #+END_SRC
+4. Elimination of =n/2-n_start= largest voids
+ #+BEGIN_SRC python
+ pattern = np.copy(initial_pattern)
+ for rank in range(n_start,int((n+1)/2)):
+ void_idx = largestVoidIndex(pattern, sigma)
+ pattern.flat[void_idx] = True
+ noise.flat[void_idx] = rank
+ #+END_SRC
+5. Elimination of =n-n/2= tightest clusters of the inverted pattern
+ #+BEGIN_SRC python
+ for rank in range(int((n+1)/2),n):
+ cluster_idx = tightestClusterIndex(np.logical_not(pattern), sigma)
+ pattern.flat[cluster_idx] = True
+ noise.flat[cluster_idx] = rank
+ #+END_SRC
+
+For each elimination the current =rank= is stored in the noise texture
+producing a 2D arrangement of the integers from 0 to =n=. As the last
+step the array is divided by =n-1= to yield a grayscale texture with values
+in $[0,1]$.
+
+#+BEGIN_SRC python
+return noise / (n-1)
+#+END_SRC
+
+In order to check whether this actually generated blue noise, we can take a
+look at the Fourier transformation for an exemplary \(100 \times 100\) texture:
+
+#+BEGIN_EXPORT html
+<div class="flexcolumns">
+<div>
+<span>Blue noise texture</span>
+<img src="https://static.kummerlaender.eu/media/blue_noise.png"/>
+</div>
+<div>
+<span>Fourier transformation</span>
+<img src="https://static.kummerlaender.eu/media/blue_noise_fourier.png"/>
+</div>
+</div>
+#+END_EXPORT
+
+One can see qualitatively that higher frequency components are significantly more
+prominent than lower ones. Contrasting this to white noise generated using uniformly
+distributed random numbers, no preference for any range of frequencies can be
+observed:
+
+#+BEGIN_EXPORT html
+<div class="flexcolumns">
+<div>
+<span>White noise texture</span>
+<img src="https://static.kummerlaender.eu/media/white_noise.png"/>
+</div>
+<div>
+<span>Fourier transformation</span>
+<img src="https://static.kummerlaender.eu/media/white_noise_fourier.png"/>
+</div>
+</div>
+#+END_EXPORT
+
+** Comparison
+Contasting the original Q criterion visualization with one produced using blue noise jittering
+followed by a soft blurring shader, we can see that the slicing artifacts largely vanish.
+While the jittering is still visible to closer inspection, the result is significantly more pleasing
+to the eye and arguably more faithful to the underlying scalar field.
+
+#+BEGIN_EXPORT html
+<div class="flexcolumns">
+<div>
+<span>Simple ray marching</span>
+<img src="https://static.kummerlaender.eu/media/q_criterion_default.png"/>
+</div>
+<div>
+<span>Ray marching with blue noise jittering</span>
+<img src="https://static.kummerlaender.eu/media/q_criterion_blue_noise.png"/>
+</div>
+</div>
+#+END_EXPORT
+
+While white noise also obcures the slices, its lower frequency components
+produce more obvious static in the resulting image compared to blue noise.
+As both kinds of noise are precomputed we can freely choose the kind of
+noise that will produce the best results for our sampling data.
+
+#+BEGIN_EXPORT html
+<div class="flexcolumns">
+<div>
+<span>Blue noise</span>
+<img src="https://static.kummerlaender.eu/media/q_criterion_blue_noise_close.png"/>
+</div>
+<div>
+<span>White noise</span>
+<img src="https://static.kummerlaender.eu/media/q_criterion_white_noise_close.png"/>
+</div>
+</div>
+#+END_EXPORT
+
+In practice where the noise is applied just-in-time during the visualization of
+a CFD simulation, all remaining artifacts tend to become invisible. This can
+be seen in the following video of the Q criterion evaluated for a simulated
+nozzle flow in LiterateLB:
+
+#+BEGIN_EXPORT html
+<video style="width:100%" src="https://static.kummerlaender.eu/media/nozzle_q_criterion.webm" controls="controls">
+</video>
+#+END_EXPORT
+
+[fn:vac] Ulichney, R. Void-and-cluster method for dither array generation. In Electronic Imaging (1993). DOI: [[https://www.spiedigitallibrary.org/conference-proceedings-of-spie/1913/0000/Void-and-cluster-method-for-dither-array-generation/10.1117/12.152707.short?SSO=1][10.1117/12.152707]].