\begin{equation}
\TVA(u) = \int_\Omega \sqrt{\nabla u(x)^T A(x) \nabla u(x)} \, dx
\end{equation}
-for all $u \in C_c^1(\Omega)$. We see that if $A(x)$ is the identity
+for all $u \in C^1(\Omega)$. We see that if $A(x)$ is the identity
matrix we get the usual total variation from \fixme{ref}. Our goal
is to control $A(x)$ such that $\nabla u(x)$ is weighted down across
edges, and thus the smoothing across the edges will be reduced.
-\fixme{Why compact support? They assume this in the anisotropic coarea
-formula paper, but I think it is enough to assume it for the test
-functions?}
As we will not always be working with differentiable images, we extend
the definition of the total variation functional to the space
&= \sup_{\eta^T A^{-1} \eta \leq 1} \int_\Omega u \diver \eta \,
dx,
\end{align}
-where $\xi$ and $\eta$ are in $C_c^1(\Omega, \mathbb{R^2})$ the $C^1$
-vector fields with compact support. \fixme{Normally it's $C_c^\infty$ for
-the test functions, but maybe that can be fixed with some density
-argument?}
+where $\xi$ and $\eta = \Ahalf \xi$ are in $C_c^\infty(\Omega,
+\mathbb{R}^2)$, the space of $C^1$ vector fields with compact support.
If we let $\norm{\xi}_A = \sqrt{\xi^T A \xi}$ and $\norm{\eta}_A^* =
\sqrt{\eta^T A^{-1} \eta}$, then our new extended definition reads
\begin{equation}
\label{eq:extended_tv}
\end{equation}
-\fixme{Energy function here?}
+\fixme{Include the whole functional somewhere.}
\fixme{Describe how this is similar to a Euler time step in a related
equation \cite{grasmair2010anisotropic}. This is only without the square
root.}
-\fixme{Extend the functional}
\subsection{Anisotropy tensor}
shown. Some problems will be worked around, while other will be skipped
with a reference to further theory.
-The basic things we ask of our functional are lower semicontinuity
+The basic things we ask of our functional
+\begin{equation}
+ F(u) = \int_\Omega (u - v)^2 + \beta \, \TVA(u)
+\end{equation}
+are lower semicontinuity
and coercivity for existence and convexity for uniqueness. We restrict
ourself to $L^2(\Omega)$, and leave the extension to $\BV(\Omega)$ to
someone else \fixme{ref}.
\begin{equation}
\int_\Omega (u - v)^2 \, dx
\end{equation}
-is obviously strictly convex, as it is $0$ for $u = v$ and positive
-elsewhere. The anisotropic total variation
+is obviously strictly convex. It can be shown by expanding and
+rearranging the strict convexity condition
+\begin{equation}
+ \int_\Omega (\lambda x + (1-\lambda)y - v)^2 \, dx < \lambda
+ \int_\Omega (x - v)^2 \, dx + (1 - \lambda) \int_\Omega (y - v)^2 \,
+ dx
+\end{equation}
+to obtain that it is equivalent to
+\begin{equation}
+ - \lambda(1 - \lambda) \int_\Omega (x - y)^2 \, dx < 0
+\end{equation}
+which is true for $0 < \lambda < 1$ and $x \neq y$.
+
+The anisotropic total variation
\begin{equation}
\TVA(u) = \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega u \diver \xi \,
dx
\end{equation}
-can be thought of as a norm, and is therefore convex. The sum of the two
-is thus strictly convex, which, given the existence of a minimizer,
-implies uniqueness.
+can be thought of as, and has the properties of a norm, and is therefore
+convex. The sum of the two is thus strictly convex, which, given the
+existence of a minimizer, implies uniqueness.
\subsection{Coercivity}
We need coercivity to show that we cannot go further and further away to
-obtain a better and better solution. This means that as $\norm{u}_{L^2}
-\to \infty$ should imply that $F(u) \to \infty$, which is obvious from
-the fidelity term for some fixed $v \in L^2(\Omega)$.
+obtain a better and better solution. This means that $\norm{u}_{L^2} \to
+\infty$ should imply that $F(u) \to \infty$, which is obvious from the
+fidelity term for some fixed $v \in L^2(\Omega)$.
From the coercivity we can conclude that we should be able to find some
-near-minimal solutions somewhere in $L^2(\Omega)$ without going too far
-away.
+near-minimal solutions somewhere in $L^2(\Omega)$ without ``going too far
+away.''
\begin{figure}
\input{fig/lower_semicont}
\end{figure}
-\subsection{Lower semicontinuity}
+\subsection{Lower semi-continuity}
The lower semicontinuity is the most tricky part, and this is where we
will take some shortcuts. Lower semicontinuity for a functional $F$ at a
For a function $f : \mathbb{R} \to \mathbb{R}$ this can be visualized as
in Figure \ref{fig:lower_semicont}.
+Since our space $L^2(\Omega)$ is of infinite dimensions things become a
+little bit problematic here. The problem lies in the fact that a
+functional which is continuous with respect to sequences is not
+neccessarily continuous with respect to the underlying topology. In
+other words, in these spaces, there is a difference between sequencial
+continuity and topological continuity. Topological continuity implies
+sequential continuity, but not the other way. One way to get around this
+would be to consider topological \emph{nets}, an extension of sequences,
+but we will stick to proving sequential lower semi-continuity and
+refering to further theory. For further reading on the theory of
+sequential versus topological continuity see for example Megginson's
+book on Banach space theory \cite{megginson}.
+
\fixme{We consider the weak topology. Because it is convenient?}
We say that a sequence $f_n$ in $L^2(\Omega)$ converges weakly to $f$ if
\end{equation}
for all $\xi \in L^2(\Omega)$ and we write $f_n \rightharpoonup f$. The
weak topology is characterized by the fact that all weakly convergent
-sequences also converge in the weak topology. Thus the mapping $u
-\mapsto \int_\Omega \xi u \, dx$ is weakly continuous for all $\xi \in
-L^2(\Omega)$. \fixme{some nets and stuff.}
-
-Before arguing that our own functional is weakly lower semi-continous, we
-present a much needed result.
+sequences also converge in the weak topology. Also, the mapping $u
+\mapsto \int_\Omega u \, \xi \, dx$ is weakly continuous for all $\xi
+\in L^2(\Omega)$. Note that when we write weakly continuous it is not a
+weaker version of continuity, but rather continuity in the weak
+topology, and the same goes for lower semi-continuity.
+
+Before arguing that our own functional is sequentially weakly lower
+semi-continous, we present a much needed result.
\begin{lemma}
Assume that the functional $F : L^2(\Omega) \to \mathbb{R}$ is
defined by
\begin{equation}
F = \sup_i F_i
\end{equation}
- where all the $F_i$ are weakly lower semi-coninuous, then $F$ is
- weakly lower semi-continuous, meaning that for any sequence $u_k
- \rightharpoonup u$ we have $F(u) \leq \liminf_k F(u_k)$.
+ where all the $F_i$ are sequentially weakly lower semi-coninuous, then $F$
+ is sequentially weakly lower semi-continuous, meaning that for any sequence
+ $u_k \rightharpoonup u$ we have $F(u) \leq \liminf_k F(u_k)$.
\end{lemma}
\begin{proof}
For any sequence $u_k \rightharpoonup u$ in $L^2(\Omega)$ we have
\begin{equation}
F(u) = \sup_i F_i(u) \leq \sup_i \liminf_{k \to \infty} F_i(u_k)
\end{equation}
- from the weak lower semi-continuity of $F_i$. Using that $\liminf_{k
- \to \infty} = \sup_k \inf_{k \leq l}$, we obtain
+ from the sequential weak lower semi-continuity of $F_i$. Using that
+ $\liminf_{k \to \infty} = \sup_k \inf_{k \leq l}$, we obtain
\begin{equation}
\begin{aligned}
F(u) &\leq \sup_i \sup_k \inf_{k \leq l} F_i(u_k) \\
&= \liminf_{k \to \infty} F(u_k)
\end{aligned}
\end{equation}
- which proves that $F$ is weakly lower semi-continuous.
+ which proves that $F$ is sequentially weakly lower semi-continuous.
\end{proof}
-\fixme{This proof is not technically correct, as we have ignored the fact that
-in the weak topology, considering sequences is not really enough, and
-one should consider \emph{nets}. But we accept it for now and refer to
-(?) for further theory.}
-For our functional in \fixme{ref}, we first consider the fidelity term
+From our functional in \fixme{ref}, we first consider the fidelity term
\begin{equation}
- \int_\Omega (u - v)^2 \, dx = \sup_{\substack{\xi \in L^2(\Omega) \\
- \norm{\xi}_{L^2} \leq \norm{u-v}_{L^2}}} \int_\Omega (u - v) \, \xi
- \, dx.
+ \int_\Omega (u - v)^2 \, dx
+ %= \sup_{\substack{\xi \in L^2(\Omega) \\ \norm{\xi}_{L^2} \leq
+ %\norm{u-v}_{L^2}}} \int_\Omega (u - v) \, \xi \, dx.
+ = \sup \left\{\int_\Omega (u - v) \, \xi \, dx : \xi \in
+ L^2(\Omega), \norm{\xi}_{L^2} \leq \norm{u-v}_{L^2} \right\}
\end{equation}
As the map $u \mapsto \int_\Omega (u - v) \xi\, dx$ is continous in the
weak topology, the fidelity term is then a supremum of weakly continous
-functionals, and is thus weakly lower semi-continuous.
+functionals, and is thus sequentially lower semi-continuous.
For the regularization term the approach is similar. With our extended
definition from \eqref{eq:extended_tv}, we have
\begin{equation}
- \TVA(u) = \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega u \, \diver \xi \,
- dx
+ %\TVA(u) = \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega u \, \diver \xi \,
+ %dx
+ \TVA(u) = \sup \left\{\int_\Omega u \, \diver \xi \, dx : \xi \in
+ C_c^\infty(\Omega, \mathbb{R}^2), \norm{\xi}_A^* \leq 1 \right\}
\end{equation}
-where $\xi \in C_c^\infty(\Omega, \mathbb{R}^2)$. This is again is a
+This is again is a
supremum of weakly continuous functionals. Thus the regularization term
-is also weakly lower semi-continuous.
+is also sequentially weakly lower semi-continuous.
+\fixme{Could consider writing the supremums as a supremum over a set
+with $\{ \}$ instead of all the subscripts!}
+\fixme{It would probably be most correct to say: weakly sequentially
+lower semi-continuous?}
-The usual ways of going from coercivity and lower semicontinuity do not
-work in infinite dimensions. See for example \fixme{Theorem 5.1 in
-\cite{scherzer2008variational}}. \fixme{And probably some other
-references to fix all the small inaccuracies}
+The usual ways of going from coercivity and lower semicontinuity to
+existence do not work in infinite dimensions. But with our coercivity
+and sequential lower semi-continuity we can apply \fixme{Theorem 5.1 in
+\cite{scherzer2008variational}}.
\section{Anisotropic coarea formula}
u = \int_0^\infty u^s \, ds.
\label{eq:positive_int}
\end{equation}
-Note that it only holds for positive image, something which has to be
-worked around in the proof of the coarea formula. For the proof we will
-avoid measure theory and follow a proof given in
-\cite{olsson2009extending}.
+Note that \eqref{eq:positive_int} only holds for a non-negative image,
+something which has to be worked around in the proof of the coarea
+formula. For the proof we will avoid measure theory and follow a proof
+given in \cite{olsson2009extending}.
\begin{figure}
\input{fig/eta_r}
Assume that $u \geq 0$ such that the integral in
\eqref{eq:positive_int} holds, then inserting into the extended
total variation definition in \eqref{eq:extended_tv} gives
- \begin{align}
- \TVA(u) &= \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega \left(
- \int_0^\infty u^s ds \right) \diver \xi \, dx
- = \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega \int_0^\infty u^s
- \diver \xi \, ds \, dx \\
- &\leq \int_0^\infty \left( \sup_{\norm{\xi}_A^* \leq 1}
- \int_\Omega u^s \diver \xi \, dx \right) \, ds
- = \int_0^\infty \TVA(u^s) \, ds.
- \end{align}
+ \begin{equation}
+ \begin{aligned}
+ \TVA(u) &= \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega
+ \left( \int_0^\infty u^s ds \right) \diver \xi \, dx
+ = \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega \int_0^\infty u^s
+ \diver \xi \, ds \, dx \\
+ &\leq \int_0^\infty
+ \left( \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega u^s \diver
+ \xi \, dx \right) \, ds
+ = \int_0^\infty \TVA(u^s) \, ds.
+ \end{aligned}
+ \label{eq:tv_u_pos}
+ \end{equation}
For $u \leq 0$ we use that $\TVA(-v) = \TVA(v)$ and that $\TVA(c + v) =
\TVA(v)$ for any constant $c$. Note that $-u \geq 0$ and that its
- thresholded image $(-u)^s$ will be exactly the opposite of $u^{-s}$.
- This allows us to show that
- \begin{align}
- \TVA(u) &= \TVA(-u) \leq \int_0^\infty \TVA \big( (-u)^r \big) \,
- dr \\
- &= \int_0^\infty \TVA(1 - u^{-r}) \, dr = \int_0^\infty
- \TVA(u^{-r}) \, dr = \int_{-\infty}^0 \TVA(u^s) \, ds.
- \end{align}
+ thresholded image $(-u)^s$ will be exactly the opposite of $u^{-s}$,
+ that is, $(-u)^s = 1 - u^{-s}$. This allows us to show that
+ \begin{equation}
+ \begin{aligned}
+ \TVA(u) &= \TVA(-u) \leq \int_0^\infty \TVA \big( (-u)^r
+ \big) \, dr \\
+ &= \int_0^\infty \TVA(1 - u^{-r}) \, dr = \int_0^\infty
+ \TVA(u^{-r}) \, dr = \int_{-\infty}^0 \TVA(u^s) \, ds.
+ \end{aligned}
+ \label{eq:tv_u_neg}
+ \end{equation}
+ Following from the supremum definition of the anisotropic total
+ variation in \eqref{eq:extended_tv}, we obtain the inequality
+ \begin{equation}
+ \begin{aligned}
+ \TVA(u_1 + u_2) &= \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega (u_1 +
+ u_2) \diver \xi \, dx \\
+ &\leq \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega u_1
+ \diver \xi \, dx + \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega u_2
+ \diver \xi \, dx\\
+ &= \TVA(u_1) + \TVA(u_2).
+ \end{aligned}
+ \end{equation}
Next, we write a general $u$ as a difference between two positive
- functions $u = u_+ - u_-$ and conclude that
- \begin{align}
- \TVA(u) &\leq \TVA(u_-) + \TVA(u_+) = \TVA(-u_-) + \TVA(u_+) \\
- &\leq \int_{-\infty}^0 \TVA\big((-u_-)^s\big) \, ds +
- \int_0^\infty \TVA(u_+^s) \, ds \\
- &= \int_{-\infty}^0 \TVA(u^s) \, ds + \int_0^\infty \TVA(u^s) \,
- ds = \int_{-\infty}^\infty \TVA(u^s) \, ds.
- \end{align}
+ functions $u = u_+ - u_-$ where $u = u_+$ when $u \geq 0$ and $u =
+ -u_-$ when $u \leq 0$. Inserting \eqref{eq:tv_u_pos} and
+ \eqref{eq:tv_u_neg} we obtain
+ \begin{equation}
+ \begin{aligned}
+ \TVA(u) &\leq \TVA(u_-) + \TVA(u_+) = \TVA(-u_-) + \TVA(u_+) \\
+ &\leq \int_{-\infty}^0 \TVA\big((-u_-)^s\big) \, ds +
+ \int_0^\infty \TVA(u_+^s) \, ds \\
+ &= \int_{-\infty}^0 \TVA(u^s) \, ds + \int_0^\infty \TVA(u^s) \,
+ ds = \int_{-\infty}^\infty \TVA(u^s) \, ds.
+ \end{aligned}
+ \end{equation}
Note that $u_+$ and $u_-$ will not be differentiable everywhere, but
we did not use the differentiability of $u$ in this part of the
proof.
\, dx,
\end{equation}
and note that $m(\infty) = \TVA(u)$. Since $m(t)$ is non-decreasing
- with $t$, we can apply the existence theorem of Lebesgue \fixme{ref}
- to conclude that $m'(t)$ exists almost everywhere and that
- Lebesgue's inequality holds:
+ with $t$, we can apply the existence theorems of Lebesgue \cite[Thm.\
+ 17.12, 18.14]{hewstrom} to conclude that $m'(t)$ exists almost
+ everywhere and that the following inequality holds:
\begin{equation}
\int_{-\infty}^\infty m'(t)\, dt \leq m(\infty) - m(-\infty) =
\TVA(u).
+ \label{eq:tva_geq_mder}
\end{equation}
Next, fix an $s \in \mathbb{R}$ and define the function
\begin{equation}
\end{equation}
for all vector fields $\xi \in C_c^\infty(\Omega, \mathbb{R}^2)$.
Assuming that $\norm{\xi}_A^* \leq 1$ we have
- \begin{align}
- \frac{m(s+r) - m(s)}{r}
- &= \frac{1}{r} \int_{\{ s \leq u \leq s+r \}} \norm{\nabla u}_A
- \, dx \\
- &\geq \frac{1}{r} \int_{\{ s \leq u \leq s + r\}} \nabla u \cdot
- \xi \, dx
- = \int_\Omega -\eta_r(u) \diver \xi \, dx.
- \end{align}
+ \begin{equation}
+ \begin{aligned}
+ \frac{m(s+r) - m(s)}{r}
+ &= \frac{1}{r} \int_{\{ s \leq u \leq s+r \}} \norm{\nabla u}_A
+ \, dx \\
+ &\geq \frac{1}{r} \int_{\{ s \leq u \leq s + r\}} \nabla u \cdot
+ \xi \, dx
+ = \int_\Omega -\eta_r(u) \diver \xi \, dx.
+ \end{aligned}
+ \end{equation}
As the limit of the left-hand side when $r \to 0$ exists almost
everywhere, suppose it exists at $s \in \mathbb{R}$, then
\begin{equation}
since $\eta_r(u) \to u^s$ when $r \to 0$. As this holds for any
$\norm{\xi}_A^* \leq 1$, we get from the extended total variation
definition in \eqref{eq:extended_tv} that $m'(s) \geq \TVA(u_s)$
- almost everywhere and conclude that
+ almost everywhere and conclude using \eqref{eq:tva_geq_mder} that
\begin{equation}
\TVA(u) \geq \int_{-\infty}^\infty m'(t) \, dt \geq
\int_{-\infty}^\infty \TVA(u^s) \, ds.
\cdot \Ahalf \eta \, dt \\
&= \int_{\partial \{ u > s\} } \sqrt{\nu_s A \nu_s} \, dt.
\end{aligned}
+ \label{eq:perimeter_calc}
\end{equation}
Here, $\nu_s$ is the unit normal of the level set $\{ u > s \}$ and by
applying the divergence theorem we have assumed that the boundary is
\input{fig/line_param}
\end{figure}
-If we parametrize straight lines as shown in Figure \ref{fig:line_param}, we can
-define the set of all staight lines as $\mathcal{L} = \{ (\phi, \rho) :
-\phi \in [0, 2\pi), \rho \geq 0 \}$. We will write a line $l_{\phi,
-\rho} = l_{\nu, \rho}$ where $\nu$ is a unit vector along the line,
-i.e.\ $\nu = (-\sin \phi, \cos \phi)^T$. By defining the
+\fixme{Something something curious geometric measure theory interesting
+integral formulas etc, introduction one or two sentences}
+
+We write $\mathcal{L}$ for the set of all straight lines in the plane,
+and parametrize them as shown in Figure \ref{fig:line_param}. Thus a
+line is parametrized by the angle $\phi$ of the normal going to the
+origin, and the distance $\rho$ from origin to the line. Sometimes it is
+more convenient to consider a unit vector $\nu$ giving the direction of
+the line instead of the angle parameter $\phi$. We will write a line
+$\ell_{\phi, \rho} = \ell_{\nu, \rho}$ where $\nu$ is a unit vector
+along the line, i.e.\ $\nu = (-\sin \phi, \cos \phi)^T$. By defining the
measure on this set $d\mathcal{L} = \dpdr$ we are ready to introduce the
Cauchy--Crofton formula, which gives us a way to calculate the length of
a curve by looking at the measure of the set of lines that intersect the
this curve $\abs{C}$ is related to the set of lines $\mathcal{L}$ as
follows
\begin{equation}
- \int_\mathcal{L} \#(l_{\phi, \rho} \cap C) \,
- d\mathcal{L}(l_{\phi, \rho}) = 2 \abs{C},
+ \int_\mathcal{L} \#(\ell_{\phi, \rho} \cap C) \,
+ d\mathcal{L}(\ell_{\phi, \rho}) = 2 \abs{C},
\end{equation}
- where $\#(l_{\phi, \rho} \cap C)$ is the number of times the line
- $l_{\phi, \rho}$ intersects the curve $C$.
+ where $\#(\ell_{\phi, \rho} \cap C)$ is the number of times the line
+ $\ell_{\phi, \rho}$ intersects the curve $C$.
\label{thm:euclidean_cauchy_crofton}
\end{theorem}
This elegant formula is very useful when we later will discretize our
The Cauchy--Crofton formula for a differentiable curve $C$ of finite
length then becomes
\begin{equation}
- \abs{C}_M = \int_\mathcal{L} \sum_{x \in l_{\nu, \rho} \cap C}
+ \abs{C}_M = \int_\mathcal{L} \sum_{x \in \ell_{\nu, \rho} \cap C}
\, \frac{\det M(x)}{2 \left( \nu^T \cdot M(x) \cdot \nu
\right)^{\sfrac{3}{2}}}
- \, d\mathcal{L}(l_{\nu, \rho}).
+ \, d\mathcal{L}(\ell_{\nu, \rho}).
\label{eq:riemannian_cauchy_crofton}
\end{equation}
\label{thm:riemannian_cauchy_crofton}
\abs{C}_M &= \int_C \sqrt{\langle \dot{C}, M \, \dot{C} \rangle} \, dt
= \int_C \sqrt{\langle \Mhalf \dot{C}, \Mhalf
\dot{C} \rangle} = \abs{M^{\sfrac{1}{2}}C} \\
- &= \int_\mathcal{L} \#(l_{\phi, \rho} \cap M^{\sfrac{1}{2}}C) \,
- d\mathcal{L}(l_{\phi, \rho}) \\
- &= \int_\mathcal{L} \#(M^{-\sfrac{1}{2}}l_{\phi, \rho} \cap C) \,
- d\mathcal{L}(l_{\phi, \rho}) \\
- &= \int_\mathcal{L} \#(m_{\phi, \rho} \cap C) \abs{J_M(l_{\phi, \rho})}\,
+ &= \int_\mathcal{L} \#(\ell_{\phi, \rho} \cap M^{\sfrac{1}{2}}C) \,
+ d\mathcal{L}(\ell_{\phi, \rho}) \\
+ &= \int_\mathcal{L} \#(M^{-\sfrac{1}{2}}\ell_{\phi, \rho} \cap C) \,
+ d\mathcal{L}(\ell_{\phi, \rho}) \\
+ &= \int_\mathcal{L} \#(m_{\phi, \rho} \cap C) \abs{J_M(\ell_{\phi, \rho})}\,
d\mathcal{L}(m_{\phi, \rho})
\end{align}
- where $J_M(l_{\phi, \rho})$ is the Jacobian of our coordinate
+ where $J_M(\ell_{\phi, \rho})$ is the Jacobian of our coordinate
transformation $F : \mathcal{L} \to \mathcal{L}$, which maps
- $l_{\phi, \rho} \mapsto \Mhalf l_{\phi, \rho}$.
+ $\ell_{\phi, \rho} \mapsto \Mhalf \ell_{\phi, \rho}$.
- As $M\in \mathbb{R}^{2\times2}$ is symmetric, so is
- $M^{\sfrac{1}{2}}$, and it admits an eigendecomposition $\Mhalf
- = U\Sigma U^T$ where the components correspond to the following
- coordinate transformations
+ We will now compute the Jacobian $J_M(\ell_{\phi, \rho})$. As $M\in
+ \mathbb{R}^{2\times2}$ is symmetric, so is $M^{\sfrac{1}{2}}$, and
+ it admits an eigendecomposition $\Mhalf = U\Sigma U^T$ where the
+ components correspond to the following coordinate transformations
\begin{align}
%U(\phi,\rho) &= (\phi + \xi, \rho), \\
%U^T(\phi,\rho) &= (\phi - \xi, \rho), \\
- U(l_{\nu,\rho}) &= l_{\phi + \xi, \rho} = l_{U\nu, \rho} \\
- U^T(l_{\nu,\rho}) &= l_{\phi - \xi, \rho} = l_{U^T\nu, \rho} \\
+ U(\ell_{\nu,\rho}) &= \ell_{\phi + \xi, \rho} = \ell_{U\nu, \rho} \\
+ U^T(\ell_{\nu,\rho}) &= \ell_{\phi - \xi, \rho} = \ell_{U^T\nu, \rho} \\
\Sigma &= \begin{pmatrix}
\sigma_1 & 0 \\
0 & \sigma_2
\end{align}
As $U$ and $U^T$ correspond to rotations they do not contribute to
the Jacobian, except for changing the input angle of the operator
- $\Lambda$ such that $J_M(l_{\phi, \rho}) = J_{\Sigma^2}(U^T l_{\phi,
+ $\Lambda$ such that $J_M(\ell_{\phi, \rho}) = J_{\Sigma^2}(U^T \ell_{\phi,
\rho})$. Given a line
\begin{equation}
- l_{\phi, \rho} = \begin{pmatrix}
+ \ell_{\phi, \rho} = \begin{pmatrix}
\rho\cdot \cos \phi \\
\rho\cdot \sin \phi
\end{pmatrix}
\end{equation}
the operator $\Sigma$ transforms it into
\begin{equation}
- \Sigma l_{\phi, \rho} = \begin{pmatrix}
+ \Sigma \ell_{\phi, \rho} = \begin{pmatrix}
\sigma_1 \rho\cdot \cos \phi \\
\sigma_2 \rho\cdot \sin \phi
\end{pmatrix}
\sigma_2 \cos \phi
\end{pmatrix},
\end{equation}
- which is the line $l_{\theta, \eta}$ with
+ which is the line $\ell_{\theta, \eta}$ with
\begin{align}
\theta &= \arctan \left( \frac{\sigma_1}{\sigma_2} \tan \phi
\right) \\
+ \sigma_2 \rho \cdot \sin \phi \cdot \sin \theta.
\end{align}
As $\partial_\rho \theta = 0$, the Jacobian becomes
- $\abs{J_{\Sigma^2}(l_{\phi, \rho})} = \partial_\phi\theta \cdot
+ $\abs{J_{\Sigma^2}(\ell_{\phi, \rho})} = \partial_\phi\theta \cdot
\partial_\rho \eta$. Differentiation yields
\begin{align}
\partial_\phi \theta &= \frac{ \frac{\sigma_1}{\sigma_2}
\phi \, \frac{\sigma_1}{\sigma_2} \tan \phi}
{\sqrt{1 + \frac{\sigma_1^2}{\sigma_2^2}
\tan^2 \phi}}
- = \frac{\sigma_1 \sigma_2 \cos^2 \phi + \sigma_1 \sigma_2
- \sin^2 \phi}{\sqrt{\sigma_1^2 \sin^2 \phi + \sigma_2^2
- \cos^2 \phi}}.
+ = \frac{\sigma_1 \sigma_2}{\sqrt{\sigma_1^2 \sin^2 \phi +
+ \sigma_2^2 \cos^2 \phi}}.
\end{equation}
If $\nu = (\nu_x, \nu_y)^T$ is a unit vector along the line
- $l_{\phi, \rho} = l_{\nu, \rho}$ then
+ $\ell_{\phi, \rho} = \ell_{\nu, \rho}$ then
\begin{equation}
- \abs{J_{\Sigma^2}(l_{\nu, \rho})}
+ \abs{J_{\Sigma^2}(\ell_{\nu, \rho})}
= \frac{\sigma_1^2 \sigma_2^2}
{\left(\sigma_1^2 \sin^2 \phi + \sigma_2^2 \cos^2 \phi
\right)^{\sfrac{3}{2}}}
\nu)^{\sfrac{3}{2}}}.
\end{equation}
We are interested in the Jacobian of the whole transformation, so
- all that is left to do is insert $U^T l_{\nu, \rho}$ to obtain
+ all that is left to do is insert $U^T \ell_{\nu, \rho}$ to obtain
\begin{equation}
- \abs{J_M(l_{\nu, \rho})} = \abs{J_{\Sigma^2}(U^T l_{\nu, \rho})}
+ \abs{J_M(\ell_{\nu, \rho})} = \abs{J_{\Sigma^2}(U^T \ell_{\nu, \rho})}
= \frac{\det M}{\left(\nu^T U \cdot \Sigma^2 \cdot U^T
\nu\right)^{\sfrac{3}{2}}}
= \frac{\det M}{\left(\nu^T \cdot M \cdot \nu\right)^{\sfrac{3}{2}}}
\begin{equation}
\abs{C}_M = \int_C \sqrt{\langle \dot{C}, M \dot{C} \rangle} \,
dt
- = \int_\mathcal{L} \# (l_{\nu, \rho} \cap C)
+ = \int_\mathcal{L} \# (\ell_{\nu, \rho} \cap C)
\frac{\det M}{\left(\nu^T \cdot M \cdot \nu
- \right)^{\sfrac{3}{2}}} \, d\mathcal{L} (l_{\nu, \rho})
+ \right)^{\sfrac{3}{2}}} \, d\mathcal{L} (\ell_{\nu, \rho})
\label{eq:riemannian_const_m}
\end{equation}
We now argue that the similar formula in
then $M_\pi(x) = M(x_i)$ for some fixed $x_i \in U_i$. We then
approximate \eqref{eq:riemannian_const_m} by
\begin{equation}
- \abs{C}_{M_\pi} = \sum_i \int_\mathcal{L} \#(l_{\nu, \rho} \cap
- C \cap U_i) \, w_i(\nu) \, d\mathcal{L}(l_{\nu, \rho})
+ \abs{C}_{M_\pi} = \sum_i \int_\mathcal{L} \#(\ell_{\nu, \rho} \cap
+ C \cap U_i) \, w_i(\nu) \, d\mathcal{L}(\ell_{\nu, \rho})
\label{eq:mpi_approx}
\end{equation}
where $w_i$ is the weight-function used in the set $U_i$
\end{equation}
Using this weight in \eqref{eq:mpi_approx} we can get rid of the sum
over the partition $i$ and form a sum of all intersection point of
- $C$ and the line $l_{\nu, \rho}$ currently being integrated over.
+ $C$ and the line $\ell_{\nu, \rho}$ currently being integrated over.
The approximation becomes
\begin{equation}
\begin{aligned}
\abs{C}_{M_\pi} &=
- \sum_i \int_\mathcal{L} \sum_{x \in l_{\nu, \rho} \cap C \cap
- U_i} w_\pi(\nu, x) \, d\mathcal{L}(l_{\nu, \rho}) \\
- &= \int_\mathcal{L} \sum_{x \in l_{\nu, \rho} \cap C}
- w_\pi(\nu, x) \, d\mathcal{L}(l_{\nu, \rho}).
+ \sum_i \int_\mathcal{L} \sum_{x \in \ell_{\nu, \rho} \cap C \cap
+ U_i} w_\pi(\nu, x) \, d\mathcal{L}(\ell_{\nu, \rho}) \\
+ &= \int_\mathcal{L} \sum_{x \in \ell_{\nu, \rho} \cap C}
+ w_\pi(\nu, x) \, d\mathcal{L}(\ell_{\nu, \rho}).
\end{aligned}
\label{eq:riemannian_segments}
\end{equation}
= \int_C \abs{\dot{C}(t)}_{M_\pi} \, dt
= \int_C \sqrt{ \dot{C}(t)^T M_\pi\big( C(t) \big) \dot{C}(t) }
\, dt.
+ \label{eq:piecewise_curve_length}
\end{equation}
We know that $M_\pi(x)$ converges pointwise to $M(x)$, and thus
$\abs{\dot{C}(t)}_{M_\pi}$ converges pointwise to
$\abs{\dot{C}(t)}_M$. We have assumed bounds on the eigenvalues
- of $M(x)$, and we can parametrize $C$ by its arc length such that
- $\abs{\dot{C}(t)} = 1$, and thus the integrand is bounded and we can
- apply Lebesgue's dominated convergence theorem to see that
+ of $M(x)$ such that according to the Rayleigh principle
+ \begin{equation}
+ K \geq \lambda_1 = \max_\xi \frac{\xi^T M_\pi(x) \xi}{ \xi^T \xi
+ }
+ \implies \xi^T M_\pi(x) \xi \leq K \, \xi^T \xi \quad \forall \xi
+ \end{equation}
+ Thus the integrand of \eqref{eq:piecewise_curve_length} is bounded
+ by $g(t) = \sqrt{K \cdot \dot{C}(t)^T \dot{C}(t)}$. We know that
+ $g(t)$ is integrable as its integral is exactly $\sqrt{K} \abs{C}$
+ and we have assumed that the curve is of finite length. This means
+ we can apply Lebesgue's dominated convergence theorem to see that
$\abs{C}_{M_\pi} \to \abs{C}_M$.
+ \fixme{ugly line spacing}
+
We apply the same theorem to show that the right hand side of
\eqref{eq:riemannian_segments} converges. Recall the definition of
$w_\pi$ in \eqref{eq:wpi_def}. The numerator is equal to $\sigma_1^2
Next we need to bound $\nu^T M_\pi(x) \nu$ away from zero. According
to the Rayleigh principle
\begin{equation}
- \lambda_2 = \min_{\norm{\xi} = 1} \sqrt{\xi^T M_\pi(x) \xi}
+ \lambda_2 = \min_{\norm{\xi} = 1} \xi^T M_\pi(x) \xi
+ \end{equation}
+ and thus $\nu^T M_\pi(x) \nu \geq \lambda_2 \geq k$. The weight
+ function $w_\pi$ is then bounded such that
+ \begin{equation}
+ \sum_{x \in \ell_{\nu, \rho} \cap C} w_\pi(\nu, x)
+ \leq
+ \sum_{x \in \ell_{\nu, \rho} \cap C}
+ \frac{K^2}{k^{\sfrac{3}{2}}}
+ =
+ \frac{K^2}{k^{\sfrac{3}{2}}} \cdot \# (\ell_{\nu, \rho} \cap C)
+ = g(\ell_{\nu, \rho}).
\end{equation}
- and thus $\nu^T M_\pi(x) \nu \geq \lambda_2^2 \geq k^2$. The weight
- function $w_\pi$ is then bounded, but not neccesarily the sum
+ This is integrable following from the Euclidean Cauchy--Crofton
+ formula in Theorem \ref{thm:euclidean_cauchy_crofton} and the fact
+ that we assumed $C$ to be of finite length
\begin{equation}
- \sum_{\mathclap{x \in l_{\nu, \rho} \cap C}} w_\pi(\nu, x).
+ \int_\mathcal{L} g(\ell_{\nu, \rho}) \, d\mathcal{L}(\ell_{\nu,
+ \rho}) = \frac{K^2}{k^{\sfrac{3}{2}}} \abs{C} < \infty.
\end{equation}
- However, if the $l_{\nu, \rho} \cap C$ is infinite for a set of
- lines with measure greater than zero, the Euclidean Cauchy--Crofton
- formula in Theorem \ref{thm:euclidean_cauchy_crofton} implies that
- the curve is infinitely long which we have assumed it is not. So as
- the sum is bounded almost everywhere we can apply Lebesgue's
- dominated convergence theorem again and conclude that
+ Thus we can apply the dominated convergence theorem again and
+ conclude that
\begin{equation}
- \int_\mathcal{L} \sum_{x \in l_{\nu, \rho} \cap C}
- w_\pi(\nu, x) \, d\mathcal{L}(l_{\nu, \rho})
+ \int_\mathcal{L} \sum_{x \in \ell_{\nu, \rho} \cap C}
+ w_\pi(\nu, x) \, d\mathcal{L}(\ell_{\nu, \rho})
\to
- \int_\mathcal{L} \sum_{x \in l_{\nu, \rho} \cap C}
- w(\nu, x) \, d\mathcal{L}(l_{\nu, \rho})
+ \int_\mathcal{L} \sum_{x \in \ell_{\nu, \rho} \cap C}
+ w(\nu, x) \, d\mathcal{L}(\ell_{\nu, \rho})
\end{equation}
which--as both sides of the equality has been shown to
converge--leaves us with what we wanted to prove
\begin{equation}
\abs{C}_M =
- \int_\mathcal{L} \sum_{x \in l_{\nu, \rho} \cap C}
+ \int_\mathcal{L} \sum_{x \in \ell_{\nu, \rho} \cap C}
\frac{\det M(x)}{2 \left( \nu^T \cdot M(x) \cdot \nu
\right)^{\sfrac{3}{2}}} \,
- d\mathcal{L}(l_{\nu, \rho}).
+ d\mathcal{L}(\ell_{\nu, \rho}).
\end{equation}
\end{proof}
setting the anisotropy tensor $A$ equal to the identity matrix $I$, the
perimeter coincides nicely with the length of the boundary curve,
assuming some regularity for the boundary. In the general case we need
-to be more careful. As can be seen in \fixme{ref}, the anisotropic
-perimeter is calculated by integrating the norm of the \emph{normal
-vector} around the boundary, while the anisotropic curve length in
-\eqref{eq:riemannian_length} is the integral of the norm of the
-\emph{tangent vector} of the curve. If $P$ is a 90\textdegree{} rotation
-matrix and we assume an arclength parametrisation of our curve we have
+to be more careful. As can be seen in \eqref{eq:perimeter_calc}, the
+anisotropic perimeter is calculated by integrating the norm of the
+\emph{normal vector} around the boundary, while the anisotropic curve
+length in \eqref{eq:riemannian_length} is the integral of the norm of
+the \emph{tangent vector} of the curve. If $P$ is a 90\textdegree{}
+rotation matrix and we assume an arclength parametrisation of our curve
+we have
\begin{equation}
\begin{aligned}
\PerA(U; \Omega) &= \int_{\partial U} \sqrt{ \langle
expression for the continuous anisotropic total variation
\begin{equation}
\TVA(u) = \int_{-\infty}^{\infty} \int_\mathcal{L} \sum_{x \in
- l_{\nu, \rho} \cap \gamma_s } \frac{\det M(x)}
+ \ell_{\nu, \rho} \cap \gamma_s } \frac{\det M(x)}
{2 \left( \nu^T \cdot M(x) \cdot \nu \right)^{\sfrac{3}{2}}}
- \, d\mathcal{L}(l_{\nu, \rho}) \, ds,
+ \, d\mathcal{L}(\ell_{\nu, \rho}) \, ds,
\end{equation}
where $\gamma_s = \partial \{ u > s \}$ is the boundary of the level set
for level $s$. Note that $P$ does not affect the determinant, i.e.\
the discrete levels to get
\begin{equation}
\int_{-\infty}^\infty \PerA( \{ u > \lambda \}; \Omega) \, d\lambda
- \approx \sum_{\lambda = 0}^{L-2} \PerA( \{ u > \lambda \}; \Omega).
+ \approx \sum_{\lambda = 0}^{L-2} \PerA( \{ u > \lambda \}; \Omega)
+ \, \Delta \lambda.
\end{equation}
Note that we do not include any $\Delta \lambda$ difference, as we can
just absorb this into the $\beta$ parameter. The perimeter is then
\subsubsection{Discrete Riemannian Cauchy--Crofton formula}
\begin{figure}
- \centering
- \includegraphics[width=0.7\textwidth]{fig/neigh_lines.png}
- \caption{
- The set of lines $\mathcal{L}$ is discretized to $\mathcal{L}_D$
- where each line belongs to a family given by the angle $\phi$
- parameter. In the figure to the left we see a visualization of
- \fixme{almost} all the families as a neighborhood. \fixme{Note
- that the angle $\phi$ is normally taken as a parameter
- giving the normal to the actual line, but since we only
- consider the difference $\Delta\phi$, this does not matter?}
- To the right we see all the lines of one family.
- }
- \label{fig:neigh_lines}
+ \input{fig/line_disc}
\end{figure}
By approximating the integral in \fixme{REF} by a discrete sum we obtain
the approximation
\begin{align}
- \abs{C}_R &= \int_\mathcal{L} n_C \, \frac{\det M(p)}{2\left(u_L^T
- \cdot M(p) \cdot u_L\right)^{\sfrac{3}{2}}} \, d\mathcal{L} \\
- &\approx \sum_{\mathcal{L}_D} n_C \, \frac{\det M(p)}{2\left(u_L^T
- \cdot M(p) \cdot u_L\right)^{\sfrac{3}{2}}} \, \Delta\phi \,
- \Delta\rho \\
+ \abs{C}_M &= \int_\mathcal{L} \sum_{x \in \ell_{\nu, \rho} \cap C}
+ \, \frac{\det M(x)}{2\left(\nu^T
+ \cdot M(x) \cdot \nu\right)^{\sfrac{3}{2}}} \,
+ d\mathcal{L}(\ell_{\nu, \rho}) \\
+ &\approx \sum_{\ell_{\nu, \rho} \in \mathcal{L}_D}
+ \sum_{x \in \ell_{\nu, \rho} \cap C} \,
+ \frac{\det M(x)}{2\left(\nu^T
+ \cdot M(x) \cdot \nu\right)^{\sfrac{3}{2}}} \, \Delta\rho \,
+ \Delta\phi \\
+ &\approx \sum_\nu \sum_\rho \sum_{x \in \ell_{\nu, \rho} \cap C} \,
+ \frac{\det M(x)}{2\left(\nu^T
+ \cdot M(x) \cdot \nu\right)^{\sfrac{3}{2}}} \, \Delta\rho \,
+ \Delta\phi.
+ \label{eq:cauchy_crofton_approx1}
\end{align}
-where $\mathcal{L}_D$ is a discretization of the set of lines in the
-plane $\mathcal{L}$ \fixme{figure}. Since the final goal is to work with
-digital images, it makes sense to discretize our domain $\Omega$ as a
-regular (?) lattice? $\mathcal{G}$. The set of lines $\mathcal{L}_D$ can
-then be made up of lines through points on this lattice as shown in
-Figure \fixme{some figure with lines}. Further working towards our graph
-cut representation of the problem later, each line is made up of
-\emph{edges} going between the points in the lattice. We will denote an
-edge by $e$. The curve length can then be approximated by
+The set of lines $\mathcal{L}$ has been discretized to the lines
+$\mathcal{L}_D$. Note that we are approximating the length of the
+\emph{differentiable} curve $C$ in $\Omega$. Further we need to
+discretize the domain $\Omega$, which will lead to discrete level sets
+$\{ u > \lambda\}$ and ``discrete'' curves.
+
+Since the final goal is to work with digital images, it makes sense to
+discretize our domain $\Omega$ as a regular grid $\mathcal{G}$. Our
+image is then reduced to a function $u : \mathcal{G} \to \mathcal{L}$.
+\fixme{mathcal L is now two things.} Moreover, the level sets $\{ u >
+\lambda\}$ will be functions taking the value of 0 or 1 on this grid, as
+shown in Figure \fixme{ref}. Thus the perimeter of these sets will be
+piecewise linear lines following the boundaries of the pixels in
+$\mathcal{G}$.
+
+The choice of our discrete set of lines $\mathcal{L}_D$ is important, as
+it will decide the accuracy of our approximation in
+\eqref{eq:cauchy_crofton_approx1}. We will only consider lines going
+through points in in our grid $\mathcal{G}$, and for now we will
+consider a discretization which is uniform throughout the domain. The
+set of lines can then be represented by the neighborhood of a pixel as
+shown in Figure \ref{fig:line_neigh}. Extending the edges shown in the
+Figure gives all lines going through the point considered. Figure
+\ref{fig:line_family} shows all lines of a given family, i.e.\ lines
+having the same angle parameter $\phi$.
+
+Thus not only have we discretized the set of lines, but each line is
+made up of \emph{edges} going from one grid point to the next. We will
+denote such an edge by $e$ or $e_{ab}$ when its endpoints are $a, b \in
+\mathcal{G}$. Thus we rewrite the discretization of
+\eqref{eq:cauchy_crofton_approx1}, and sum over all the edges in the
+discretization $\mathcal{L}_D$ to obtain
\begin{equation}
- \abs{C}_R \approx \sum_{e} n_C(e) \, \frac{\det M(e) \norm{e}^3}{2
- \left(e^T \cdot M(e) \cdot e\right)^{\sfrac{3}{2}}} \, \Delta\phi \,
- \Delta\rho.
+ \abs{C}_M \approx \sum_{e} \sum_{x \in e \cap C} \, \frac{\det M(x)
+ \norm{e}^3}{2 \left(e^T \cdot M(x) \cdot e\right)^{\sfrac{3}{2}}} \,
+ \Delta\phi \, \Delta\rho.
\end{equation}
-As our metric tensor $M$ is only known in the grid points, we
-approximate the tensor in the intersection point $p$ (between the curve
-$C$ and the edge $e$) by
+This is beginning to look like something we can calculate. One
+difficulty is finding the intersections $e \cap C$. The exact
+calculations of these points will not fit into our graph cut framework
+later, and thus for an edge $e$ we will consider only the question of
+``did $e$ cross $C$ or not?'' This amounts to checking whether the
+terminals of $e$ lie on each side of the perimeter $C$, and the
+approximation is exact for zero or one intersection points, but will, as
+we see in Figure \fixme{ref}, not be entirely correct when we have more.
+
+\fixme{curve or perimeter here, maybe perimeter because then we know it
+follows the boundaries of the pixels.}
+
+The second difficulty is that in the discrete setting, we will only have
+an approximation of the metric tensor $M(x)$ for each point $x \in
+\mathcal{G}$, and it is thus not available for arbitrary intersection
+points in $\Omega$. For an edge $e$ we will utilize the average of the
+tensor in the two endpoints of the edge. Thus for an intersection point
+$x$ somewhere on the edge $e_{ab}$, we approximate the metric tensor by
\begin{equation}
- M(p) \approx M(e) = \frac{M\big(p_1(e)\big) + M\big(p_2(e)\big)}{2},
+ M(x) \approx M(e_{ab}) = \frac{M(a) + M(b)}{2},
\end{equation}
the component-wise average of the tensors in the two end points of the
-edge.
+edge. \fixme{really? componentwise? will that not mess up the
+eigenvalues?}
+
\begin{figure}
- \centering
- \includegraphics[width=0.7\textwidth]{fig/area_argument.png}
- \caption{
- We see that for any reasonable family of lines (i.e.\ we do
- not leave any unneccessary gaps), the area covered by a blue
- square, $\delta^2$ is equal to the area covered by a red
- rectangle $\Delta \rho \norm{e_k}$, as there will be the same
- number of blue squares and red rectangles if we cover the whole
- plane.
- }
- \label{fig:area_argument}
+ \input{fig/area_proof}
\end{figure}
-Assuming some regularity on our curve, and that our discretization is
-somehow fine, we approximate $n_C(e)$ by a function which is $1$ if $C$
-crosses $e$ an odd number of times, and $0$ if not, meaning that we
-ignore multiple passes across $e$ and only care about if $C$ goes from
-one side of $e$ to the other. Further, as can be seen in Figure
-\fixme{some figure with lines etc.}, the distance between lines in line
-family $k$ is $\Delta \rho_k = \delta^2 / \norm{e_k}$. Thus the curve
-length is approximated by
+Further, as can be seen in Figure \ref{fig:area_proof}
+the distance between lines in a given line family is $\Delta \rho =
+\delta^2 / \norm{e}$. Thus the curve length is approximated by
\begin{equation}
- \abs{C}_R \approx \sum_{e\cap C} \frac{\det M(e) \norm{e}^2
+ \abs{C}_M \approx \sum_{e \cap C} \frac{\det M(e) \norm{e}^2
\, \delta^2 \, \Delta\phi}{2 \left(e^T \cdot M(e) \cdot
e\right)^{\sfrac{3}{2}}}.
\end{equation}
+Going back to the perimeter of the level set $\{ u > \lambda \}$ we get
+\begin{equation}
+ \PerA(\{u > \lambda\}; \mathcal{G}) = \sum_{e_{ab}} \abs{
+ u^\lambda(a) - u^\lambda(b)} \frac{\det M(e_{ab}) \norm{e_{ab}}^2
+ \, \delta^2 \, \Delta\phi}{2 \left(e_{ab}^T \cdot M(e_{ab}) \cdot
+ e_{ab}\right)^{\sfrac{3}{2}}}.
+\end{equation}
+The absolute value $\abs{u^\lambda(a) - u^\lambda(b)}$ is one if one of
+$a$ and $b$ lie inside the level set and the other lies outside, and
+zero otherwise. In other words if the absolute value is one if $e_{ab}$
+crosses the perimeter of $\{ u > \lambda\}$ an odd number of times, and
+zero otherwise.
\fixme{We now have the problem that $C$ is a curve, but later we want it
to be a cut...}
-\subsubsection{Consistency}
+\subsubsection{Consistency, stability and convergence}
+
+As for any discretization these are important properties to consider.
+Consistency relates to how well a solution to the continuous problem
+fits in the discretized equation. This would optimally be controlled by
+the resolution of the discretization in such a way that if we refine the
+discretization, the
Consider a square centered around grid point $p$ with side lengths
$\sqrt{\delta}$. The size of this square will go to zero, while the