Fluid Dynamics¶

Important: Please read the installation page for details about how to install the toolboxes. $\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\si}{\sigma}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$

This numerical tour explores fluid dynamics for image generation.

In [2]:

# options(warn=-1) # turns off warnings, to turn on: "options(warn=0)"

library(png)
library(pracma)
library(imager)

for (f in list.files(path="nt_toolbox/toolbox_general/", pattern="*.R")) {
    source(paste("nt_toolbox/toolbox_general/", f, sep=""))
}

for (f in list.files(path="nt_toolbox/toolbox_signal/", pattern="*.R")) {
    source(paste("nt_toolbox/toolbox_signal/", f, sep=""))
}

source("nt_toolbox/toolbox_wavelet_meshes/meshgrid.R")
options(repr.plot.width=5, repr.plot.height=5)

Velocity Flow Field¶

A velocity flow is simply a 2-D vector field $V = (V_i)_{i=1}^N \in \RR^{n \times n \times 2}$ where $V_i \in \RR^2$ is one of the $N=n \times n$ vectors at a position indexed by $i$.

It can be generated as a realization of Gaussian process. The blurring creates correlations in the flow.

In [2]:

n <- 128
V <- perform_blurring( array( rnorm(n*n*2), c(n,n,2) ), c(40), "per" )

[1] 127 127
[1] 127 127
[1] 127 127
[1] 127 127

Subsampling display operator.

In [3]:

myplot <- function(V){ plot_vf(V[seq(1,n,6),seq(1,n,6),]) }

We can display the vector field using arrow.

In [4]:

myplot(V)

We can renormalize the flow, which enhances the singularities. It defines $\tilde V$ as $\tilde V_i = V_i/\norm{V_i}$.

In [5]:

normalize <- function(V){ V/ array(rep( pmax( array(1e-9, dim(V)[1:2]), sqrt(apply(V**2, c(1,2), sum)) ), 2 ), dim(V)) }

Display.

In [6]:

myplot(normalize(V))

Incompressible Flows¶

An incompressible flow has a vanishing divergence. The set of vector incompressible flow defines a sub-space of $\RR^{n \times n \times 2}$ $$ \Ii = \enscond{V}{ \text{div}(V)=0 } \qwhereq \text{div}(V) = \pd{V}{x_1} + \pd{V}{x_2} \in \RR^{n \times n}. $$ Here $\pd{}{x_s}$ for $s=1,2$ are finite differences approximation of the horizontal and vertical derivative operators (we suppose here periodic boundary conditions).

The orthogonal projection $U = \text{Proj}_{\Ii}(V)$ on $\Ii$ is computed by solving a Poisson equation $$ U = V-\nabla A \qwhereq \Delta A = \text{div}(V). $$

This is especially simple for periodic boundary conditions since $A$ can be computed over the Fourier domain as $$ \forall \om \neq 0, \quad \hat A(\om) = \frac{\hat Y(\om)}{\mu(\om)} \qwhereq Y = \text{div}(V) \qandq \mu(\om_1,\om_2) = -4 \sin(\om_1 \pi / n)^2 -4 \sin(\om_2 \pi / n)^2 $$ and $\hat A(0)=0$.

Compute the kernel $\mu(\om)$.

In [7]:

grid <- meshgrid_2d(0:(n-1), 0:(n-1))
Y <- grid$X ; X <- grid$Y
mu <- sin(X*pi/n)**2
mu <- -4*(mu + t(mu))
mu[1,1] <- 1

Computation of $A$.

In [8]:

A <- function(V){ FFT <- fft( div(V[,,1], V[,,2], bound="per") )
                  return( Re( fft( FFT/mu, inverse=TRUE)/length(FFT) ) ) }

Projection on incompressible flows.

In [9]:

ProjI <- function(V){ V - grad(A(V), bound="per") }

Display $U=\text{Proj}_{\Ii}(V)$.

In [10]:

U <- ProjI(V)

myplot(U)

Display $W=U-V$ the irrotational component of $V$.

In [11]:

myplot(V-U)

Note that the decomposition $V=U+W$ is called the Hoge decomposition of the vector field.

Image Advection Along the Flow¶

A flow defines a warping operator that transport the content of an image along the streaming of the flow.

We load an image $f$.

In [12]:

f <- as.matrix(load_image("nt_toolbox/data/lena.png", 2*n))
f <- f[(n-round(n/2)+1):(n+round(n/2)),(n-round(n/2)+1):(n+round(n/2))]

Given some vector field $U$, the warping operator $f_1 = \Ww_U(f)$ along the flow is defined $$ f_1(x) = f(x+U(x)) $$ i.e. it advects the values of $f$ by the vector field $U$ to obtain the values of $f_1$.

We define $U$ as a scaled normalized incompressible flow.

In [13]:

U <- normalize(ProjI(V))

Helper function: enforce periodicity.

In [14]:

periodic <- function(P){ mod(P,n) }

Helper function: extend an image by 1 pixel to avoid boundary problems.

In [15]:

extend1 <- function(f){ ext <- array(0, c(dim(f)[1], dim(f)[2]+1))
                        ext[,1:dim(f)[2]] <- f
                        ext[,dim(ext)[2]] <- f[,1]
                        return(ext) }

extend <- function(f){ t(extend1(t(extend1(f)))) }

Helper function: bilinear interpolation on a grid.

In [16]:

myinterp <- function(f1, Pi){
    dim_f1 <- dim(f1)
    f1 <- as.cimg(t(f1))
    locations <- data.frame(x=as.vector(Pi[,,2]), y=as.vector(Pi[,,1]))
    return( as.matrix(as.cimg(imager::interp(f1, locations))) )
}

First we compute the initial and wraped grids.

In [17]:

grid <- meshgrid_2d(1:n, 1:n)
Y <- grid$X ; X <- grid$Y
P <- array(0, c(dim(X),2)) ; P[,,1] <- X ; P[,,2] <- Y

Defines the warping operator $\Ww_U$.

In [18]:

W <- function(f, U){ myinterp(extend(f), periodic(P - U)) }

Display a warped image $\Ww_{\rho U}(f)$ for some scaling $\rho$.

In [19]:

rho <- 2
options(repr.plot.width=5, repr.plot.height=5)
imageplot(W(f, rho*U))

Exercise 1

Display $\Ww_{\rho U}(f)$ for various values of $\rho$.

In [20]:

options(repr.plot.width=7, repr.plot.height=7)

source("nt_solutions/graphics_5_fluids/exo1.R")