aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--README.md7
-rw-r--r--articles/2020-05-26_lambda_tuple_swallowing.md2
-rw-r--r--articles/2021-09-26_noise_and_ray_marching.org254
-rw-r--r--default.nix32
-rw-r--r--flake.lock97
-rw-r--r--flake.nix28
l---------tags/development/2021-09-26_noise_and_ray_marching.org1
l---------tags/english/2021-09-26_noise_and_ray_marching.org1
l---------tags/math/2021-09-26_noise_and_ray_marching.org1
9 files changed, 390 insertions, 33 deletions
diff --git a/README.md b/README.md
index 9205238..206c262 100644
--- a/README.md
+++ b/README.md
@@ -3,3 +3,10 @@
This repository contains all textual contents of "blog.kummerlaender.eu" and is included into the [appropriate repository](https://github.com/KnairdA/blog.kummerlaender.eu/) as a subtree.
All articles and pages are freely available under the terms of the Creative Commons [CC-BY-SA](http://creativecommons.org/licenses/by-sa/3.0/) license. If you should require further permissions please feel free to contact me.
+
+## Build
+
+```sh
+nix flake clone git+https://code.kummerlaender.eu/blog_content --dest source
+nix build source/
+``` \ No newline at end of file
diff --git a/articles/2020-05-26_lambda_tuple_swallowing.md b/articles/2020-05-26_lambda_tuple_swallowing.md
index f64ea52..b33c046 100644
--- a/articles/2020-05-26_lambda_tuple_swallowing.md
+++ b/articles/2020-05-26_lambda_tuple_swallowing.md
@@ -1,6 +1,6 @@
# Working with tuples using swallowing and generic lambdas
-Suppose you have some kind of list of types. Such a list by can by itself be [used](/article/using_scheme_as_a_metaphor_for_template_metaprogramming/) to perform any compile time computation one might come up with. So let us suppose that you additionally want to construct a tuple from something that is based on this list. i.e. you want to connect the compile time only type list to a run time object. In such a case you might run into new question such as: How do I call constructors for each of my tuple values? How do I offer access to the tuple values using only the type as a reference? How do I call a function for each value in the tuple while preserving the connection to the compile time list? If such questions are of interest to you, this article might possibly also be.
+Suppose you have some kind of list of types. Such a list can by itself be [used](/article/using_scheme_as_a_metaphor_for_template_metaprogramming/) to perform any compile time computation one might come up with. So let us suppose that you additionally want to construct a tuple from something that is based on this list. i.e. you want to connect the compile time only type list to a run time object. In such a case you might run into new question such as: How do I call constructors for each of my tuple values? How do I offer access to the tuple values using only the type as a reference? How do I call a function for each value in the tuple while preserving the connection to the compile time list? If such questions are of interest to you, this article might possibly also be.
While the standard's tuple template is part of the C++ subset I use in basically all of my developments[^0] I recently had to revisit some of these questions while reworking OpenLB's core data structure using its [_meta descriptor_](/article/meta_descriptor/) concept. The starting point for this was a class template called `FieldArrayD` to store an array of instances of a single field in a SIMD vectorization friendly _structure of arrays_ layout. As a LBM lattice in practice stores not just one such field type but multiple of them (all declared in the central _descriptor_ structure) I then wanted a `MultiFieldArrayD` class template that does just that. i.e. a simple wrapper that accepts a list of fields as a variadic template parameter pack and instantiates a `FieldArrayD` for each of them. A sensible place for storing these instances is of course our trusty `std::tuple`:
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
diff --git a/default.nix b/default.nix
deleted file mode 100644
index a54c326..0000000
--- a/default.nix
+++ /dev/null
@@ -1,32 +0,0 @@
-{ pkgs ? import <nixpkgs> { }, mypkgs ? import <mypkgs> { }, ... }:
-
-pkgs.stdenv.mkDerivation {
- name = "blog.kummerlaender.eu";
-
- src = pkgs.fetchFromGitHub {
- owner = "KnairdA";
- repo = "blog.kummerlaender.eu";
- rev = "a043d5dd1933e4fa9cfa2b10a7fdfa05c6c4d0eb";
- sha256 = "0ykprjw97125miw8pqih3pd8hk2sdc1cginakg0p944svs0p6811";
- };
-
- LANG = "en_US.UTF-8";
-
- buildInputs = [
- pkgs.pandoc
- pkgs.highlight
- mypkgs.katex-wrapper
- mypkgs.make-xslt
- ];
-
- installPhase = ''
- mkdir source/00_content
- cp -r ${./articles} source/00_content/articles
- cp -r ${./tags} source/00_content/tags
- cp ${./meta.xml} source/00_content/meta.xml
-
- make-xslt
- mkdir $out
- cp -Lr target/99_result/* $out
- '';
-}
diff --git a/flake.lock b/flake.lock
new file mode 100644
index 0000000..4015285
--- /dev/null
+++ b/flake.lock
@@ -0,0 +1,97 @@
+{
+ "nodes": {
+ "cms": {
+ "inputs": {
+ "nixpkgs": "nixpkgs",
+ "pkgs-personal": "pkgs-personal"
+ },
+ "locked": {
+ "lastModified": 1632431502,
+ "narHash": "sha256-ZTauOeZwYIzKuQJoM4xmofwdNg9derVFz9vc2Lifh1E=",
+ "ref": "master",
+ "rev": "3fe48fff2bb624400a015e734bcecee245df5949",
+ "revCount": 250,
+ "type": "git",
+ "url": "https://code.kummerlaender.eu/blog.kummerlaender.eu"
+ },
+ "original": {
+ "type": "git",
+ "url": "https://code.kummerlaender.eu/blog.kummerlaender.eu"
+ }
+ },
+ "nixpkgs": {
+ "locked": {
+ "lastModified": 1629379628,
+ "narHash": "sha256-dI8wpEo7wIVWoTUk2oyWFUnlVHNKLs+ren1TqITN1mI=",
+ "owner": "NixOS",
+ "repo": "nixpkgs",
+ "rev": "a1007637cea374bd1bafd754cfd5388894c49129",
+ "type": "github"
+ },
+ "original": {
+ "owner": "NixOS",
+ "ref": "nixos-21.05",
+ "repo": "nixpkgs",
+ "type": "github"
+ }
+ },
+ "nixpkgs_2": {
+ "locked": {
+ "lastModified": 1629271619,
+ "narHash": "sha256-by9D3OkEKk4rOzJIMbC0uP2wP3Bt81auP5xmbmPg2a8=",
+ "owner": "NixOS",
+ "repo": "nixpkgs",
+ "rev": "7bbca9877caed472c6b5866ea09302cfcdce3dbf",
+ "type": "github"
+ },
+ "original": {
+ "owner": "NixOS",
+ "ref": "nixos-21.05",
+ "repo": "nixpkgs",
+ "type": "github"
+ }
+ },
+ "pkgs-personal": {
+ "inputs": {
+ "nixpkgs": "nixpkgs_2"
+ },
+ "locked": {
+ "lastModified": 1629652608,
+ "narHash": "sha256-eNcsdqMyK/Q3P0Tj16uDNcQzKIFf4CJkM7qTq3BdtF0=",
+ "ref": "master",
+ "rev": "fb63603b5eec859c84464e1a7f6f14931303f679",
+ "revCount": 52,
+ "type": "git",
+ "url": "https://code.kummerlaender.eu/pkgs"
+ },
+ "original": {
+ "type": "git",
+ "url": "https://code.kummerlaender.eu/pkgs"
+ }
+ },
+ "root": {
+ "inputs": {
+ "cms": "cms",
+ "stable": "stable"
+ }
+ },
+ "stable": {
+ "locked": {
+ "lastModified": 1632291606,
+ "narHash": "sha256-oEN24XJYAFK9tsD13TzLEizpgQigEfgC6i9x1b/1pVU=",
+ "owner": "NixOS",
+ "repo": "nixpkgs",
+ "rev": "83413f47809790e4ca012e314e7782adeae36cf2",
+ "type": "github"
+ },
+ "original": {
+ "owner": "NixOS",
+ "ref": "nixos-21.05",
+ "repo": "nixpkgs",
+ "type": "github"
+ }
+ }
+ },
+ "root": "root",
+ "version": 7
+}
diff --git a/flake.nix b/flake.nix
new file mode 100644
index 0000000..5241ae6
--- /dev/null
+++ b/flake.nix
@@ -0,0 +1,28 @@
+{
+ description = "blog.kummerlaender.eu";
+
+ inputs = {
+ stable.url = github:NixOS/nixpkgs/nixos-21.05;
+ cms.url = git+https://code.kummerlaender.eu/blog.kummerlaender.eu;
+ };
+
+ outputs = { self, stable, cms, ... }: {
+ defaultPackage.x86_64-linux = cms.generate ./.;
+
+ defaultApp.x86_64-linux = let
+ system = "x86_64-linux";
+
+ pkgs = import stable { inherit system; };
+
+ serve = pkgs.writeScriptBin "serve" ''
+ #!/bin/sh
+ pushd ${cms.generate ./.}
+ ${pkgs.gatling}/bin/gatling
+ popd
+ '';
+ in {
+ type = "app";
+ program = "${serve}/bin/serve";
+ };
+ };
+}
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