diff options
| -rw-r--r-- | articles/2021-09-26_noise_and_ray_marching.org | 254 | ||||
| -rw-r--r-- | flake.lock | 8 | ||||
| l--------- | tags/development/2021-09-26_noise_and_ray_marching.org | 1 | ||||
| l--------- | tags/english/2021-09-26_noise_and_ray_marching.org | 1 | ||||
| l--------- | tags/math/2021-09-26_noise_and_ray_marching.org | 1 |
5 files changed, 261 insertions, 4 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..c5ae921 --- /dev/null +++ b/articles/2021-09-26_noise_and_ray_marching.org @@ -0,0 +1,254 @@ +* 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 gradients. 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 +Let $C$ be the sample density along a ray $r$ with length $L$ and given an absorption $\mu$. Then +\[ C(r) = \int_0^L c(x) \mu(x) \exp\left(-\int_0^x \mu(t) dt\right) dx \] +is the volume rendering equation yielding per-pixel density values that are easily mapped to +some convenient color palette. + +#+BEGIN_EXPORT html +<video style="width:100%" src="https://literatelb.org/media/nozzle.webm" controls="controls"> +</video> +#+END_EXPORT + +I.e. we integrate over the values $c(x)$ along the ray weighted by the current absorption $\mu$ +and the accumulated absorption up to the current point. This way samples that are closer to +the view origin will be more prominent than samples of the same magnitude that are farther +away which of course is also a rough approximation of how a real volume behaves, hence +the name. + +In practice this integral is approximated by the sum +$$C(r) = \sum_{i=0}^N c(i \Delta x) \alpha (i \Delta x) \prod_{j=0}^{i-1} \left(1 - \alpha(j\Delta x)\right)$$ +for step width \(\Delta x \in \mathbb{R}^+\) in the most basic case. + +This approach may be extended arbitrarily, e.g. it is only 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 samples \(c(x)\) in the volume rendering +equation. 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 style="display:flex;flex-direction:row;flex-wrap:wrap;width:100%"> +<div style="display:flex;flex-direction:column;flex-basis:100%;flex:1"> +<span style="text-align:center">Q Criterion</span> +<img src="https://static.kummerlaender.eu/media/q_criterion_default.png"/> +</div> +<div style="display:flex;flex-direction:column;flex-basis:100%;flex:1"> +<span style="text-align:center">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 +Volumetric renderings that are produced when using a plain ray marching algorithm may contain +undesired artifacts. Specifically the choice of start offsets w.r.t. to the view origin can lead +to slicing artifacts. While this tends to become less noticable for smaller step widths these are +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 produce a distracting static-like pattern. A common choice in practice is to use so +called /blue noise/ instead. 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 (doi:10.1117/12.152707)]] provides a straight forward method for +pre-computing tileable blue noise textures to be reused during the visualization. + +The first ingredient for this algorithm is a =filteredPattern= function that applies a +plain Gaussian filter with given $\sigma$ to a cyclic 2d array. Using cyclic wrapping here 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): + n = np.prod(shape) + n_start = int(n / 10) + + initial_pattern = initialPattern(shape, n_start, sigma) + noise = np.zeros(shape) + + 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 + + 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 + + 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 + + return noise / (n-1) +#+END_SRC + +The actual algorithm utilizes these three helper functions in four steps: +1. Initial pattern generation +2. Eliminiation of =n_start= tightest clusters +3. Elimination of =n/2-n_start= largest voids +4. Elimination of =n-n/2= tightest clusters of inverse pattern +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]$. + +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 style="display:flex;flex-direction:row;flex-wrap:wrap;width:100%"> +<div style="display:flex;flex-direction:column;flex-basis:100%;flex:1"> +<span style="text-align:center">Blue noise texture</span> +<img src="https://static.kummerlaender.eu/media/blue_noise.png"/> +</div> +<div style="display:flex;flex-direction:column;flex-basis:100%;flex:1"> +<span style="text-align:center">Fourier transformation</span> +<img src="https://static.kummerlaender.eu/media/blue_noise_fourier.png"/> +</div> +</div> +#+END_EXPORT + +One can see clearly 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 style="display:flex;flex-direction:row;flex-wrap:wrap;width:100%"> +<div style="display:flex;flex-direction:column;flex-basis:100%;flex:1"> +<span style="text-align:center">White noise texture</span> +<img src="https://static.kummerlaender.eu/media/white_noise.png"/> +</div> +<div style="display:flex;flex-direction:column;flex-basis:100%;flex:1"> +<span style="text-align:center">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 to one produced using blue noise jittering +followed by a soft blurring shader, we can see that the slicing artifacts mostly vanish. +While the jittering is apparent if one looks closely, the result is more visually pleasing +also arguably more faithful to the underlying scalar field. + +#+BEGIN_EXPORT html +<div style="display:flex;flex-direction:row;flex-wrap:wrap;width:100%"> +<div style="display:flex;flex-direction:column;flex-basis:100%;flex:1"> +<span style="text-align:center">Simple ray marching</span> +<img src="https://static.kummerlaender.eu/media/q_criterion_default.png"/> +</div> +<div style="display:flex;flex-direction:column;flex-basis:100%;flex:1"> +<span style="text-align:center">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 the additional lower frequency components +produce a more obvious pattern in the resulting image compared to blue noise. + +#+BEGIN_EXPORT html +<div style="display:flex;flex-direction:row;flex-wrap:wrap;width:100%"> +<div style="display:flex;flex-direction:column;flex-basis:100%;flex:1"> +<span style="text-align:center">Blue noise</span> +<img src="https://static.kummerlaender.eu/media/q_criterion_blue_noise_close.png"/> +</div> +<div style="display:flex;flex-direction:column;flex-basis:100%;flex:1"> +<span style="text-align:center">White noise</span> +<img src="https://static.kummerlaender.eu/media/q_criterion_white_noise_close.png"/> +</div> +</div> +#+END_EXPORT + +#+BEGIN_EXPORT html +<video style="width:100%" src="https://static.kummerlaender.eu/media/nozzle_q_criterion.webm" controls="controls"> +</video> +#+END_EXPORT @@ -6,11 +6,11 @@ "pkgs-personal": "pkgs-personal" }, "locked": { - "lastModified": 1632338037, - "narHash": "sha256-Uj/6702xXVq0UDrY2v/LtmUQBguNpfuVmNyKisc4+Gg=", + "lastModified": 1632431502, + "narHash": "sha256-ZTauOeZwYIzKuQJoM4xmofwdNg9derVFz9vc2Lifh1E=", "ref": "master", - "rev": "24373f1e7d2f536efac22caa71de50223aca12d7", - "revCount": 249, + "rev": "3fe48fff2bb624400a015e734bcecee245df5949", + "revCount": 250, "type": "git", "url": "https://code.kummerlaender.eu/blog.kummerlaender.eu" }, diff --git a/tags/development/2021-09-26_noise_and_ray_marching.org b/tags/development/2021-09-26_noise_and_ray_marching.org new file mode 120000 index 0000000..8b90651 --- /dev/null +++ b/tags/development/2021-09-26_noise_and_ray_marching.org @@ -0,0 +1 @@ +../../articles/2021-09-26_noise_and_ray_marching.org
\ No newline at end of file diff --git a/tags/english/2021-09-26_noise_and_ray_marching.org b/tags/english/2021-09-26_noise_and_ray_marching.org new file mode 120000 index 0000000..8b90651 --- /dev/null +++ b/tags/english/2021-09-26_noise_and_ray_marching.org @@ -0,0 +1 @@ +../../articles/2021-09-26_noise_and_ray_marching.org
\ No newline at end of file diff --git a/tags/math/2021-09-26_noise_and_ray_marching.org b/tags/math/2021-09-26_noise_and_ray_marching.org new file mode 120000 index 0000000..8b90651 --- /dev/null +++ b/tags/math/2021-09-26_noise_and_ray_marching.org @@ -0,0 +1 @@ +../../articles/2021-09-26_noise_and_ray_marching.org
\ No newline at end of file |