aboutsummaryrefslogtreecommitdiff log msg author committer range
blob: c5ae92118f39acd6e9c840693964c7fb0c12c405 (plain)
 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254  pre { line-height: 125%; } td.linenos .normal { color: inherit; background-color: transparent; padding-left: 5px; padding-right: 5px; } span.linenos { color: inherit; background-color: transparent; padding-left: 5px; padding-right: 5px; } td.linenos .special { color: #000000; background-color: #ffffc0; padding-left: 5px; padding-right: 5px; } span.linenos.special { color: #000000; background-color: #ffffc0; padding-left: 5px; padding-right: 5px; } .highlight .hll { background-color: #ffffcc } .highlight .c { color: #888888 } /* Comment */ .highlight .err { color: #a61717; background-color: #e3d2d2 } /* Error */ .highlight .k { color: #008800; font-weight: bold } /* Keyword */ .highlight .ch { color: #888888 } /* Comment.Hashbang */ .highlight .cm { color: #888888 } /* Comment.Multiline */ .highlight .cp { color: #cc0000; font-weight: bold } /* Comment.Preproc */ .highlight .cpf { color: #888888 } /* Comment.PreprocFile */ .highlight .c1 { color: #888888 } /* Comment.Single */ .highlight .cs { color: #cc0000; font-weight: bold; background-color: #fff0f0 } /* Comment.Special */ .highlight .gd { color: #000000; background-color: #ffdddd } /* Generic.Deleted */ .highlight .ge { font-style: italic } /* Generic.Emph */ .highlight .gr { color: #aa0000 } /* Generic.Error */ .highlight .gh { color: #333333 } /* Generic.Heading */ .highlight .gi { color: #000000; background-color: #ddffdd } /* Generic.Inserted */ .highlight .go { color: #888888 } /* Generic.Output */ .highlight .gp { color: #555555 } /* Generic.Prompt */ .highlight .gs { font-weight: bold } /* Generic.Strong */ .highlight .gu { color: #666666 } /* Generic.Subheading */ .highlight .gt { color: #aa0000 } /* Generic.Traceback */ .highlight .kc { color: #008800; font-weight: bold } /* Keyword.Constant */ .highlight .kd { color: #008800; font-weight: bold } /* Keyword.Declaration */ .highlight .kn { color: #008800; font-weight: bold } /* Keyword.Namespace */ .highlight .kp { color: #008800 } /* Keyword.Pseudo */ .highlight .kr { color: #008800; font-weight: bold } /* Keyword.Reserved */ .highlight .kt { color: #888888; font-weight: bold } /* Keyword.Type */ .highlight .m { color: #0000DD; font-weight: bold } /* Literal.Number */ .highlight .s { color: #dd2200; background-color: #fff0f0 } /* Literal.String */ .highlight .na { color: #336699 } /* Name.Attribute */ .highlight .nb { color: #003388 } /* Name.Builtin */ .highlight .nc { color: #bb0066; font-weight: bold } /* Name.Class */ .highlight .no { color: #003366; font-weight: bold } /* Name.Constant */ .highlight .nd { color: #555555 } /* Name.Decorator */ .highlight .ne { color: #bb0066; font-weight: bold } /* Name.Exception */ .highlight .nf { color: #0066bb; font-weight: bold } /* Name.Function */ .highlight .nl { color: #336699; font-style: italic } /* Name.Label */ .highlight .nn { color: #bb0066; font-weight: bold } /* Name.Namespace */ .highlight .py { color: #336699; font-weight: bold } /* Name.Property */ .highlight .nt { color: #bb0066; font-weight: bold } /* Name.Tag */ .highlight .nv { color: #336699 } /* Name.Variable */ .highlight .ow { color: #008800 } /* Operator.Word */ .highlight .w { color: #bbbbbb } /* Text.Whitespace */ .highlight .mb { color: #0000DD; font-weight: bold } /* Literal.Number.Bin */ .highlight .mf { color: #0000DD; font-weight: bold } /* Literal.Number.Float */ .highlight .mh { color: #0000DD; font-weight: bold } /* Literal.Number.Hex */ .highlight .mi { color: #0000DD; font-weight: bold } /* Literal.Number.Integer */ .highlight .mo { color: #0000DD; font-weight: bold } /* Literal.Number.Oct */ .highlight .sa { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Affix */ .highlight .sb { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Backtick */ .highlight .sc { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Char */ .highlight .dl { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Delimiter */ .highlight .sd { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Doc */ .highlight .s2 { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Double */ .highlight .se { color: #0044dd; background-color: #fff0f0 } /* Literal.String.Escape */ .highlight .sh { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Heredoc */ .highlight .si { color: #3333bb; background-color: #fff0f0 } /* Literal.String.Interpol */ .highlight .sx { color: #22bb22; background-color: #f0fff0 } /* Literal.String.Other */ .highlight .sr { color: #008800; background-color: #fff0ff } /* Literal.String.Regex */ .highlight .s1 { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Single */ .highlight .ss { color: #aa6600; background-color: #fff0f0 } /* Literal.String.Symbol */ .highlight .bp { color: #003388 } /* Name.Builtin.Pseudo */ .highlight .fm { color: #0066bb; font-weight: bold } /* Name.Function.Magic */ .highlight .vc { color: #336699 } /* Name.Variable.Class */ .highlight .vg { color: #dd7700 } /* Name.Variable.Global */ .highlight .vi { color: #3333bb } /* Name.Variable.Instance */ .highlight .vm { color: #336699 } /* Name.Variable.Magic */ .highlight .il { color: #0000DD; font-weight: bold } /* Literal.Number.Integer.Long */* 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 #+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
Q Criterion
Curl Norm
#+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
Blue noise texture
Fourier transformation
#+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
White noise texture
Fourier transformation
#+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
Simple ray marching
Ray marching with blue noise jittering
#+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
Blue noise
White noise
#+END_EXPORT #+BEGIN_EXPORT html #+END_EXPORT