The Beltrami Flow

Let us denote by $(\Sigma,g)$ the image manifold and its metric and by $(M,h)$ the space-feature manifold and its metric, then the map ${X}:\Sigma\to M$ has the following measure:

$$S[X^i,g_{\mu\nu},h_{ij}]=\int d^m\sigma \sqrt{g}g^{\mu\nu}\partial_{\mu}X^i \partial_{\nu}X^j h_{ij}({X})$$ (1)

called the Polyakov action [16], where $m$ is the dimension of $\Sigma$, $g$ is the determinant of the image metric, $g^{\mu\nu}$ is the inverse of the image metric, $\mu,\nu \in [1,\dots ,\dim\Sigma]$, $i,j \in [1,\dots ,\dim M]$, and $h_{ij}$ is the metric of the embedding space.

This is a natural generalization of the $L2$ norm to manifolds.

As an example, a grey level volume of images can be treated as a $3D$ manifold embedded in $R^4$, i.e. a mapping

$$X : (x,y,z) \rightarrow \left( X^1 = x, X^2 = y, X^3 = z, X^4 = u(x,y,z) \right)$$ (2)

Many scale-space methods, linear and non-linear can be shown to be a gradient descent flows of this functional with appropriately chosen metric of the image manifold.

The gradient descent equation is

$$X^i_t = - \frac{1}{\sqrt{g}} \frac{\delta S}{\delta X^i}$$ (3)

As shown in [3], minimizing the area action in equation (1), with respect to the feature coordinate $u$, is equivalent to solving the following equation,

$$u_t = \displaystyle{\frac{u_{xx}(1+{u_y}^2+{u_z}^2) + u_{yy}(1+{u_x}^2+{u_z}^2) + u_{zz}(1+{u_x}^2+{u_y}^2)}{( 1 + {u_x}^2 + {u_y}^2 + {u_z}^2)^2}} +$$

$$\displaystyle{\frac{- 2 u_x u_y u_{xy} - 2 u_x u_z u_{xz} - 2 u_y u_z u_{yz}} {( 1 + {u_x}^2 + {u_y}^2 + {u_z}^2)^2}}$$ (4)

which is called the Beltrami flow

Implicit Scheme for Beltrami Flow

As noted by Weickert in [1], the classical nonlinear diffusion equation can be re-written as the following reaction-diffusion partial differential equation

$$u_t = \nabla \cdot \left (\frac{\nabla u}{g}\right ) = \fra... ...ght) + \frac{\partial}{\partial z} \left(\frac{u_z}{g} \right)$$ (5)

with $g(x,y) = (1+{u_x}^2 + {u_y}^2 + {u_z}^2)$. The Beltrami equation (4) may be reduced to a similar reaction-diffusion form, namely

$$u_t = \nabla \cdot \left(\frac{\nabla u}{2g} \right) + \frac... ...bla \frac{1}{g}\right) \cdot \nabla u + \frac{1}{g} \nabla^2 u$$ (6)

where $g$ is the edge indicator function. In this form, the Beltrami flow equation is not a ``pure'' diffusion equation. It has both a parabolic edge-preserving and a hyperbolic edge-sharpening terms.

Another modification is to replace $\frac{\partial}{\partial x}\frac{1}{g}\frac{\partial u}{\partial x} = \frac{\pa... ...partial u}{\partial x}\right) - \frac{1}{g} \frac{\partial^2 u}{{\partial x}^2}$ in equation (6). This enable to replace the discretization of $\left(\frac{\partial}{\partial x}\frac{1}{g}\right)\frac{\partial u}{\partial x}$ by a second order derivative term. Therefore the resulting linear system will be a tri-diagonal matrix, thus more stable than a system with no value on the diagonal.

In addition, the reaction-diffusion form of equation (6) hides the mixed derivatives, thereby making it conducive to the AOS approach. In other words, the equation can be rearranged into the form $u_t =(\mathbf{A}_x + \mathbf{A}_y + \mathbf{A}_z) u$, where $\mathbf{A}_x, \mathbf{A}_y$ and $\mathbf{A}_z$ are the following differential operators:

$$\left\{ \begin{array}{c} \displaystyle{\mathbf{A}_x} = \frac{1}{2... ...ght) + \frac{1}{g} \frac{\partial^2}{{\partial z}^2} \right) \end{array}\right.$$ (7)

Applying the backward difference formula to the equation (7) we get

$$\frac{u^{n+1} - u^n}{\Delta t } = (\mathbf{A}_x+\mathbf{A}_y+\mathbf{A}_z) u^{n+1}$$ (8)

The superscript $n$ is related to the present and $n+1$ to the next time step. The subscripts $i,j$ index the discrete pixel location; $u_{i,j}^n$ are known values, and $u_{i,j}^{n+1}$ are to be found. Using $u^{n+1}$ on the right side of equation (8) makes the integration scheme implicit and unconditionally stable, namely

$$[\mathbf{I}- \Delta t (\mathbf{A}_x+\mathbf{A}_y+\mathbf{A}_z)]u^{n+1} = u^n$$ (9)

where $\mathbf{I}$ is the identity matrix. Before proceeding in time, we calculate the values of the edge indicator function $g$, using the known values of $u^n$. Thus, the scheme is only semi-implicit. Although $g$ depends on the gradient of $u$, we treat it like a given function of $(x,y)$, making the governing PDE ``quasi-linear''.

Note that equation (9) includes a large bandwidth matrix, because all equations, related to new pixel values $u^{n+1}$ are coupled. Our aim is to decouple the set in equation (9) so that each row and each column of pixels can be handled separately. For this, we re-arrange the equations into the following form:

$$u^{n+1} = [\mathbf{I}-\Delta t (\mathbf{A}_x+\mathbf{A}_y+\mathbf{A}_z)]^{-1}u^n$$ (10)

Of course, we do not intend to invert the matrix to solve the linear set. This is only a symbolic form used for further derivation. For a small value of $\Delta t$, the matrix in the brackets on the right side of Equation (10) is close to the identity $\mathbf{I}$. Thus, its inverse can be expanded into the Taylor series in the proximity of $\mathbf{I}$: $[\mathbf{I}-\Delta t (\mathbf{A}_x+\mathbf{A}_y+\mathbf{A}_z)]^{-1} \approx \mathbf{I}+ \Delta t (\mathbf{A}_x+\mathbf{A}_y+\mathbf{A}_z)$, where the linear term is retained and the high order terms are neglected. Introducing this form into equation (10), we get

$$3 u^{n+1} = (\mathbf{I}+3\Delta t \mathbf{A}_x) u^n + (\math... ... t \mathbf{A}_y) u^n + (\mathbf{I}+3\Delta t \mathbf{A}_z) u^n$$ (11)

Introducing the notations $V=(\mathbf{I}+ 3 \Delta t \mathbf{A}_x ) u^n$, $W = (\mathbf{I}+ 3 \Delta t \mathbf{A}_y ) u^n$ and $Z = (\mathbf{I}+ 3 \Delta t \mathbf{A}_z ) u^n$ the solution is simply

$$u^{n+1} = \frac{V + W + Z}{3}$$ (12)

In order to get an implicit scheme, we apply the differential matrix operators $\mathbf{A}_x$ and $\mathbf{A}_y$ to $u^{n+1}$ (and not to $u^n$), namely

$$(\mathbf{I}+3\Delta t \mathbf{A}_x)^{-1} V = u^n \hbox{, }\h... ...}\hspace{4mm} (\mathbf{I}+3\Delta t \mathbf{A}_z)^{-1} Z = u^n$$ (13)

Following the procedure of expanding the matrix inverses into Taylor series and applying the linearization for small $\Delta t$, we finally obtain the equation sets for $V$ and $W$ as follows:

$$(\mathbf{I}-3\Delta t \mathbf{A}_x) V = u^n \hbox{, }\hspace... ...ox{, }\hspace{4mm} (\mathbf{I}-3\Delta t \mathbf{A}_z) Z = u^n$$ (14)

This leads for $V$ to the following numerical scheme

$$- \alpha_m \Delta t V_{i-1} + [1+(\alpha_m+\alpha_p)\Delta t]V_i - \alpha_p \Delta t V_{i+1} = u_i^n$$ (15)

where

$$\alpha_m = \frac{3}{4}\frac{h_{i-1} + 3 h_i}{{\Delta x}^2} \... ...mm} \alpha_p = \frac{3}{4}\frac{h_{i+1} + 3 h_i}{{\Delta x}^2}$$ (16)

and a similar equation in $y$, and $z$.

These equations can be solved with either the Dirichlet or Neumann boundary conditions; these and other details are described in [17].

Simulation Results for Beltrami Flow

1st Test: Varying the time-step:

We ran a series of numerical simulations to demonstrate the performance of the implicit scheme for the Beltrami flow on a 3D medical dataset. The original dataset, represented in figure 1 is a 3D X-Ray Angiography of the head of a patient.

Figure 1: Three orthogonal slices of a CTA scanner in the Circle-of-Willis

A dye product has been injected in an artery of the patient before acquisition. This product is radio-active and enhance the contrast between the arteries and the neighboring tissues of the brain. With this image, a clinician can visually inspect the region of the Circle-of-Willis where cerebral aneurysms occur. An aneurysm ([18]) is an abnormal dilatation of the artery wall. Under high pressures, it can burst, leading to a stroke. To study this pathology, the clinician can use the three orthogonal slices of the data represented in figure 1. He can also use a Maximum-Intensity-Projection (MIP) image, and we are going to filter out the noise inside this scanned data in order to visualize this projection image.

Considering the Courant-Friedrich-Levy (CFL) condition [19] for the explicit scheme, int the case of an anisotropic cubic grid the maximum time step allowed should be $\frac{1}{2}\left(\frac{1}{{\Delta x}^2+{\Delta y}^2 +{\Delta z}^2}\right)$. The unconditionally stable scheme that we are using here enables us to multiply this time step by an acceleration factor $f$. In the following test, we are varying this acceleration factor, and displaying the MIP image in each case. The results are shown in figure 2.

Figure 2: Results of Beltrami flow solved until a fixed scale 250 and with different acceleration factors

Original image

$f = 1$, $t = 1457$

$f = 2$, $t = 728$

$f = 5$, $t = 291$

$f = 10$, $t = 145$

$f = 20$, $t = 73$

$f = 50$, $t = 29$

$f = 100$, $t = 14$

$f = 200$, $t = 7$

The MIP images represent a volume of interest around the Circle-of-Willis. As we see, the implicit scheme is always stable, but for $f \gg 50$, the resulting accuracy may be insufficient for certain applications. At the same time it is important to notice that small vessels cleared for a small acceleration factor like $f=1,\ldots,10$ are enhanced for bigger $f$. The reason is that due to accumulated errors, the problem with $f = 200$ seems like it is run to a slightly lower scale than the same problem with $f=10,20$, thus giving different results. Finally, observation of this dataset clearly depicts a particularity in the region of the Circle-of-Willis: the upper part of the cerebral arteries and the lower part do not seem connected. This observation should help the clinician in proper visualization and interpretation of this data.

2nd Test: Varying the ratio between reaction and diffusion:

The next series of numerical simulations are carried out to study the edge enhancement effect on gray level images. We fix the acceleration factor to $100$ and, inspired by the decomposition of equation (5) which leads to the reaction-diffusion form of equation (6), we solve the following normalized reaction-diffusion equation in this series:

$$\displaystyle{\frac{\partial u}{\partial t}} = \cos \beta \nabla h \cdot \nabla u + \sin \beta h \nabla^2 u$$ (17)

The first term on the right side of equation (17) is a reaction term, while the second is a diffusion term and $\beta$ is a parameter controlling the relative contribution of these opposing effects. The reaction term is responsible for edge enhancement, while the diffusion term smooths the noise away from the edges. Notice that for $\beta = \arctan{2}$, this reaction-diffusion equation (17) is equivalent to the Beltrami flow (4), multiplied by a factor of $2\cos{(\arctan{2})}$. Results of varying $\beta$ between $\frac{\pi}{9}$ and $\frac{\pi}{2}$ are shown in figure 3.

Figure 3: Results of solving equation (17) until scale $=250$ with different values of $\beta$

Original image

$\beta = 20$, reaction

$\beta = \frac{\pi}{6}$

$\beta = \frac{\pi}{4}$

$\beta = \arctan{2}$, Beltrami flow

$\beta = \frac{\pi}{2}$, Pure diffusion

The first row in figure 3, presents the initial image and the results for $\beta = \frac{\pi}{9}$ (high reaction) and $\beta = \pi/6$. The reason why the image for $\beta = \frac{\pi}{9}$ looks smoother than for $\beta = \frac{\pi}{6}$ is rather simple: due to the instability of the pure reaction process, the image needs to be smoothed for small values of $\beta$. The second row presents the results for $\beta = \pi/4$ (a nonlinear diffusion flow equation), $\beta = \arctan 2$ (the Beltrami flow), and $\beta = \pi/2$ (scaled ``linear'' diffusion). According to equation (17), the edge enhancement effect should decay with increasing $\beta$. Indeed, we see that the edge enhancement is stronger for the nonlinear diffusion flow ($\beta = \pi/4$) than for the Beltrami flow ( $\beta = \arctan{2}$).

Notice that due to the additional component in $z$ in the denominator of the scheme 4, the $3D$ Beltrami flow is slower than the corresponding $2D$ flow studied in [20]. This means that to achieve the same degree of noise reduction and edge enhancement, larger scale values should be employed in $3D$.