\sigma_1 & 0 \\
0 & \sigma_2
\end{pmatrix}.
+ \label{eq:sigma_construction}
\end{equation}
and for $\sigma_1$ and $\sigma_2$ we choose
\begin{equation}
dx \\
&= \sup_{\norm{\xi}_A^* \leq 1} \int_{\partial \{ u > s\} } \nu_s
\cdot \xi \, dt \\
- &= \sup_{\abs{\eta} \leq 1} \int_{\partial \{ u > s\} } \nu_s
+ &= \sup_{\norm{\eta} \leq 1} \int_{\partial \{ u > s\} } \nu_s
\cdot \Ahalf \eta \, dt \\
&= \int_{\partial \{ u > s\} } \Ahalf \nu_s
- \cdot \frac{\Ahalf \nu_s}{\abs{\Ahalf \nu_s}} \, dt \\
+ \cdot \frac{\Ahalf \nu_s}{\norm{\Ahalf \nu_s}} \, dt \\
&= \int_{\partial \{ u > s\} } \sqrt{\nu_s A \nu_s} \, dt.
\end{aligned}
\label{eq:perimeter_calc}
\begin{proof}
See \cite[Theorem 3, Section 1-7.]{do1976differential}.
\end{proof}
-This elegant formula is very useful when we later will discretize our
-perimeter calculation. The set of lines $\mathcal{L}$ is then discretized
-in a reasonable way, and the length of the curve $C$ can be approximated
-by a sum over all these lines.
-
If our space is equipped with a metric tensor $M(x)$ such that the inner
product of two vectors $a$ and $b$ in a point $x$ is calculated as
$\langle a, b\rangle_M = \langle a, M(x) b \rangle$, then the length of
\end{equation}
We will now present and prove a Cauchy--Crofton formula in this case
where our domain is equipped with a metric tensor in each point.
+This elegant formula is very useful when we later will discretize our
+perimeter calculation. The set of lines $\mathcal{L}$ is then discretized
+in a reasonable way, and the length of the curve $C$ can be approximated
+by a sum over all these lines.
\begin{theorem}[The Riemannian Cauchy--Crofton formula]
- Assume that our space $\Omega$ is equipped with a continuous metric
+ Assume that our space $\Omega$ is equipped with a continuous
+ positive definite metric
tensor $M(x)$, whose eigenvalues are bounded by $0 < k \leq
\lambda_2 \leq \lambda_1 \leq K < \infty$ for all $x \in \Omega$.
The Cauchy--Crofton formula for a differentiable curve $C$ of finite
\label{thm:riemannian_cauchy_crofton}
\end{theorem}
\begin{proof}[Proof of the Riemannian Cauchy--Crofton formula]
- Assume first that our space is equipped with at constant metric
- tensor $M$. The length of our curve using this tensor can be
+ Assume first that our space is equipped with a constant metric
+ tensor $M$. The length of a curve in this space can be
calculated by transforming the curve and applying the Euclidean
Cauchy--Crofton formula
\begin{align}
&= \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})
+ d\mathcal{L}(m_{\phi, \rho}).
\end{align}
- where $J_M(\ell_{\phi, \rho})$ is the Jacobian of our coordinate
+ Here $J_M(\ell_{\phi, \rho})$ is the Jacobian of the coordinate
transformation $F : \mathcal{L} \to \mathcal{L}$, which maps
$\ell_{\phi, \rho} \mapsto \Mhalf \ell_{\phi, \rho}$.
\nomenclature{$J_M(\ell_{\phi, \rho})$}{Jacobian of the coordinate
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
+ it admits a decomposition $\Mhalf = U\Sigma U^T$ where the
components correspond to the following coordinate transformations
\begin{align}
%U(\phi,\rho) &= (\phi + \xi, \rho), \\
\Sigma &= \begin{pmatrix}
\sigma_1 & 0 \\
0 & \sigma_2
- \end{pmatrix}.
+ \end{pmatrix}
+ = \begin{pmatrix}
+ \sqrt{\lambda_1} & 0 \\
+ 0 & \sqrt{\lambda_2}
+ \end{pmatrix}
\end{align}
As $U$ and $U^T$ correspond to rotations and our measure
- $\mathcal{L}$ is invariant under rigid motions, $U$ and $U^T$ do not
+ $\mathcal{L}$ is invariant under rigid motions, $U$ and $U^T$ do not
have direct contributions to the Jacobian. They do however affect the
input angle of the operator $\Lambda$ such that $J_M(\ell_{\phi,
- \rho}) = J_{\Sigma^2}(U^T \ell_{\phi, \rho})$. Given a line
+ \rho}) = J_{\Sigma^2}(U^T \ell_{\phi, \rho})$. Thus we will now
+ compute $J_{\Sigma^2}(\ell_{\phi, \rho})$. Given a line
\begin{equation}
\ell_{\phi, \rho} = \begin{pmatrix}
\rho\cdot \cos \phi \\
\sigma_2 \cos \phi
\end{pmatrix},
\end{equation}
- which is the line $\ell_{\theta, \eta}$ with
+ which equals the line $\ell_{\theta, \eta}$ with
\begin{align}
\theta &= \arctan \left( \frac{\sigma_1}{\sigma_2} \tan \phi
- \right) \\
+ \right) \label{eq:theta}\\
\eta &= \left\langle \begin{pmatrix}
\sigma_1 \rho\cdot \cos \phi \\
\sigma_2 \rho\cdot \sin \phi
\partial_\rho \eta &= \sigma_1 \cos \phi \cdot \cos \theta +
\sigma_2 \sin \phi \cdot \sin \theta.
\end{align}
- In the expression for $\partial_\rho \eta$ we use that $\sin
+ In the expression for $\partial_\rho \eta$ we insert $\theta$ from
+ \eqref{eq:theta} and use that $\sin
(\arctan( x )) = x / \sqrt{1 + x^2}$ and that $\cos ( \arctan ( x ))
= 1 / \sqrt{1 + x^2}$ to obtain
\begin{equation}
= \frac{\det \Sigma^2}{(\nu^T \cdot \Sigma^2 \cdot
\nu)^{\sfrac{3}{2}}}.
\end{equation}
- We are interested in the Jacobian of the whole transformation, so
+ We are interested in the Jacobian of the whole transformation
+ $J_{\Sigma^2}(U^T \ell_{\theta, \rho})$, so
all that is left to do is insert $U^T \ell_{\nu, \rho}$ to obtain
\begin{equation}
\abs{J_M(\ell_{\nu, \rho})} = \abs{J_{\Sigma^2}(U^T \ell_{\nu, \rho})}
dt
= \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} (\ell_{\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
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$
+ where $w_i$ is the weight-function used in the set $U_i$, that is,
\begin{equation}
w_i(\nu) = \frac{\det M(x_i)}{\left(\nu^T \cdot M(x_i) \cdot \nu
\right)^{\sfrac{3}{2}}}.
\end{equation}
We further simplify the approximation by introducing the global
weight-function $w_\pi(\nu, x)$ which is equal to
- $w_i(\nu)$ when $x \in U_i$. It can be written
+ $w_i(\nu)$ when $x \in U_i$. It can be written as
\begin{equation}
w_\pi(\nu, x) = \frac{\det M_\pi(x)}{\left(\nu^T \cdot M_\pi(x)
\cdot \nu \right)^{\sfrac{3}{2}}}.
\end{aligned}
\label{eq:riemannian_segments}
\end{equation}
- As our partition $\pi$ is refined the weight $w_\pi(x)$ converges
+ Now it only remains
+ to show that the left- and right-hand side of
+ \eqref{eq:riemannian_segments} converges to the left- and right-hand
+ side of \eqref{eq:riemannian_cauchy_crofton}.
+
+ As our partition $\pi$ is refined, the weight $w_\pi(x)$ converges
pointwise to the continuously varying weight
\begin{equation}
w(\nu, x) = \frac{\det M(x)}{\left( \nu^T \cdot M(x) \cdot \nu
\right)^{\sfrac{3}{2}}}
\end{equation}
- found in \eqref{eq:riemannian_cauchy_crofton}. Now it only remains
- to show that the left and right hand side of
- \eqref{eq:riemannian_segments} converges to the left and right hand
- side of \eqref{eq:riemannian_cauchy_crofton}.
+ found in \eqref{eq:riemannian_cauchy_crofton}.
- Recall from \eqref{eq:riemannian_length} that the left hand side is
+ Recall from \eqref{eq:riemannian_length} that the left-hand side is
calculated as
\begin{equation}
\abs{C}_{M_\pi}
We know that $M_\pi(x)$ converges pointwise to $M(x)$, and thus
$\labs{\dot{C}(t)}_{M_\pi}$ converges pointwise to
$\labs{\dot{C}(t)}_M$. We have assumed bounds on the eigenvalues
- of $M(x)$ such that according to the Rayleigh principle
+ 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}
+ and therefore we have the bound
+ \begin{equation}
+ \xi^T M_\pi(x) \xi \leq K \, \norm{\xi}^2, \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
we can apply Lebesgue's dominated convergence theorem to see that
$\labs{C}_{M_\pi} \to \labs{C}_M$.
- We apply the same theorem to show that the right hand side of
+ 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
\sigma_2^2 = \lambda_1 \lambda_2$ and is by assumption bounded from
above by $K^2$.
- \fixme{sigma squared meh?}
-
Next we need to bound $\nu^T M_\pi(x) \nu$ away from zero. According
to the Rayleigh principle
\begin{equation}
\frac{K^2}{k^{\sfrac{3}{2}}}
=
\frac{K^2}{k^{\sfrac{3}{2}}} \cdot \# (\ell_{\nu, \rho} \cap C)
- = g(\ell_{\nu, \rho}).
+ =: g(\ell_{\nu, \rho}).
\end{equation}
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
+ that we assumed $C$ to be of finite length:
\begin{equation}
\int_\mathcal{L} g(\ell_{\nu, \rho}) \, d\mathcal{L}(\ell_{\nu,
\rho}) = \frac{K^2}{k^{\sfrac{3}{2}}} \abs{C} < \infty.
\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
+ which---as both sides of the equality in
+ \eqref{eq:riemannian_segments} have been shown to converge---leaves
+ us with what we wanted to prove
\begin{equation}
\abs{C}_M =
\int_\mathcal{L} \sum_{x \in \ell_{\nu, \rho} \cap C}
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. The normal and tangent vector
-are always perpendicular, so all we need is a 90\textdegree{} rotation,
-assuming that the tangent vector is unit length. It will be if it arises
-from an arc length parametrization of the boundary. If $P$ is a
-90\textdegree{} rotation matrix we have
+the \emph{tangent vector} of the curve. Thus a 90\textdegree{} rotation
+is needed.
+
+If $P$ is a 90\textdegree{} rotation matrix we have
\begin{equation}
\begin{aligned}
\PerA(U; \Omega) &= \int_{\partial U} \sqrt{\langle
\end{equation}
We simplify the equation by defining the metric tensor $M(x) = P A(x)
P^T$ and letting $\gamma = \partial U$ be an arclength parametrization of
-the boundary of $U$. Observe that a 90\textdegree{} rotation of the
-normal $\nu_{\partial U}$ gives us the tangent $\dot{\gamma}$.
-Inserting this into \eqref{eq:per_to_length1} we get
+the boundary of $U$ to obtain
\begin{equation}
\PerA(U; \Omega) = \int_\gamma \sqrt{\langle \dot{\gamma}, M(x) \,
\dot{\gamma} \rangle }.
\otimes \nabla f_{\sigma} \right)\right)(x),
\end{equation}
Because of the convolutions with the Gaussian function, this is a smooth
-continuous map from $\bar{\Omega}$ to $\mathbb{R}^{2 \times 2}$. As we
+map from $\bar{\Omega}$ to $\mathbb{R}^{2 \times 2}$. As we
can see in \eqref{eq:s_eigenvalues}, the eigenvalues
depend continuously on the coefficients of the
elements in
-the structure tensor $S_\rho(x)$. The extreme value theorem \fixme{ref}
-states that a continuous real-valued function on a nonempty compact
-space is bounded above. Thus the eigenvalues $\lambda_1$ and $\lambda_2$ of $S_\rho(x)$
+the structure tensor $S_\rho(x)$. The extreme value theorem states that
+a continuous real-valued function on a nonempty compact
+space is bounded above. Thus the eigenvalues $s_1 \geq s_2$ of $S_\rho(x)$
are bounded above and by the construction in
-\eqref{eq:sigma_construction}, the smallest eigenvalue of
-our anisotropy tensor $A(x)$ is bounded away from zero as
+\eqref{eq:sigma_construction}, there exists uniform bound $k$ such that
+the smallest eigenvalue of the anisotropy tensor $A(x)$ is bounded away
+from zero, as
\begin{equation}
\lambda_2 = \left(1
+ \frac{(s_1 - s_2)^2}{\omega^2}\right)^{-1} \geq \left(1 +
\end{equation}
Hence, our metric tensor $M(x) = P A(x) P^T$ and curve length
calculation in
-\eqref{eq:per_to_length2} fulfill all the assumptions of
-Theorem~\ref{thm:riemannian_cauchy_crofton}. \fixme{did we argue for
-continuous?} Thus we can apply the Cauchy--Crofton formula of
-\eqref{eq:riemannian_cauchy_crofton} to calculate the curve length in
-\eqref{eq:per_to_length2} and through the anisotropic coarea formula in
-Theorem~\ref{thm:anisotropic_coarea} obtain our final expression for the
-continuous anisotropic total variation
+\eqref{eq:per_to_length2} fulfill all the assumptions of the anisotropic
+Cauchy--Crofton formula in Theorem~\ref{thm:riemannian_cauchy_crofton}.
+\fixme{did we argue for
+continuous? and pos def?} Thus we can apply the formula to calculate the
+perimeter in \eqref{eq:per_to_length2} as
\begin{equation}
- \TVA(u) = \int_{-\infty}^{\infty} \int_\mathcal{L} \sum_{x \in
- \ell_{\nu, \rho} \cap C_s } \frac{\det M(x)}
+ \PerA(U; \Omega) = \int_\mathcal{L} \sum_{x \in
+ \ell_{\nu, \rho} \cap \gamma } \frac{\det M(x)}
{2 \left( \nu^T \cdot M(x) \cdot \nu \right)^{\sfrac{3}{2}}}
\, d\mathcal{L}(\ell_{\nu, \rho}) \, ds,
\end{equation}
-where $C_s = \partial \{ u > s \}$ is the boundary of the level set
-for level $s$. Note that $P$ does not affect the determinant, i.e.\
-$\det A = \det PAP^T = \det M$, and from our eigendecomposition in
-\fixme{ref} we see that the transformation $PAP^T \to M$ actually
-amounts to switching the two eigenvalues $\sigma_1$ and $\sigma_2$ in
-$\Sigma$.
+where $\gamma = \partial U$.
+Note that $P$ does not affect the determinant, i.e.\
+$\det A = \det PAP^T = \det M$, and from our decomposition in
+\eqref{eq:sigma_construction} we see that the transformation $PAP^T \to
+M$ actually amounts to switching the two eigenvalues $\sigma_1$ and
+$\sigma_2$ in $\Sigma$.
+
+This concludes the treatment of the continuous problem. We have seen how
+the anisotropic coarea formula in Theorem~\ref{thm:anisotropic_coarea}
+allows us to calculate the anisotropic total variation as an integral of
+the perimeter of all the level sets. Through the anisotropic
+Cauchy--Crofton formula in Theorem~\ref{thm:riemannian_cauchy_crofton}
+these perimeters are calculated by an integral over the set of all
+lines. We are then left with the functional
+\begin{equation}
+ F(u) = \int_\Omega (u - f)^2 + \beta \, \TVA(u),
+ \label{eq:final_cont_rep}
+\end{equation}
+where
+\begin{equation}
+ \TVA(u) = \int_{-\infty}^\infty \int_{\mathcal{L}} \sum_{x \in
+ \ell_{\nu, \rho} \cap C_s} \frac{\det M(x)}
+ {2 \left( \nu^T \cdot M(x) \cdot \nu \right)^{\sfrac{3}{2}}}
+ \, d\mathcal{L}(\ell_{\nu, \rho}) \, ds.
+\end{equation}
+\fixme{define $C_s$. $C_s$ vs $\gamma$.}
+Within the restrictions that these theorems put on the tensor $M(x)$, we
+have chosen a construction where one eigenvalue is always 1, while the
+other varies from 1 in smooth areas towards 0 around edges, with the
+corresponding eigenvector perpendicular to the edge. \fixme{meh}
-\fixme{Small summary, complete functional, tensor construction.}
+As mentioned previously, these transformations are motivated by the
+discretization that follows in the next chapter.
\chapter{Discrete formulation}
The whole transformation from the initial functional in
\eqref{eq:first_anisotropic_functional} through the anisotropic coarea
formula and the Cauchy--Crofton formula was motivated by the
-discretization which will be described here. The anisotropic coarea
-formula allows us to minimize the functional for each level separately,
-while the Cauchy--Crofton formula gives a feasible way of calculating
-the perimeter of each level set.
-
-The integrals we had in \fixme{ref} will be approximated by discrete
-sums, and we will be careful to discretize in a way which is consistent
-with the initial continuous formulation.
+discrete formulation which will be described here. After discretizing
+the functional, we will se how a graph cut approach can be used to find
+a global minimizer in polynomial time.
\section{Discretization}
\fixme{Rating: 4/10}
-\fixme{Sometimes I write uniform grid, and sometimes regular. Uniform is
-probably better.}
+Assume that our discrete images are given on a uniform grid $\mathcal{G}$,
+where each grid point is called a pixel. The image is a function giving
+each pixel a value in the set of levels $\mathcal{P} = \{0, \hdots,
+L-1\}$. This is a reasonable assumption for digital grayscale
+images.
+
+Thus, when discretizing the functional in \eqref{eq:final_cont_rep}, we
+have to consider that our images now have both discrete domain and
+co-domain.
-We assume that our discrete images are given on a uniform grid $\mathcal{G}$,
-where each discrete point is called a pixel and we further assume that
-each pixel takes a value in the set of levels $\mathcal{P} = \{0,
-\hdots, L-1\}$. This is a reasonable assumption for digital grayscale
-images. We now want to discretize the energy function in
-\eqref{eq:continuous_energy}.
+The integrals in \eqref{eq:final_cont_rep} will be approximated by
+discrete sums. First the fidelity term is discretized without too much
+trouble, while with the regularization term, there is more choice as to
+how to discretize the set of lines $\mathcal{L}$. In the end we will
+verify that our discretization is consistent with the continuous
+functional.
\subsection{Fidelity term}
level value $k \in
\mathcal{P}$, we define the following function
\begin{equation}
- N_x(k) = \abs{k - v_x}^p
+ N_x(k) = \abs{k - f_x}^2
\end{equation}
which is the value of the fidelity term if we were to give $u_x$ a value
of $k$. This allows us write
\begin{equation}
- \int_\Omega \abs{u - v}^p\, dx \approx \sum_{x \in \mathcal{G}}
- \abs{u_x - v_x}^p \Delta x = \sum_{x \in \mathcal{G}} N_x(u_x)
+ \int_\Omega \abs{u - f}^2\, dx \approx \sum_{x \in \mathcal{G}}
+ \abs{u_x - f_x}^2 \Delta x = \sum_{x \in \mathcal{G}} N_x(u_x)
\Delta x.
\label{eq:fidelity_approx_1}
\end{equation}
true, and 0 if $x$ is false. Since $I(\lambda < u_x) = u_x^\lambda$ we
rewrite \eqref{eq:fidelity_approx_1} and obtain
\begin{equation}
- \sum_{x \in \mathcal{G}} \abs{u_x - v_x}^p =
+ \sum_{x \in \mathcal{G}} \abs{u_x - f_x}^2 =
\sum_{x \in \mathcal{G}} N_x(u_x) = \sum_{\lambda = 0}^{L-2} \sum_{x
\in \mathcal{G}}
\big( N_x(\lambda + 1) - N_x(\lambda) \big) \, u_x^\lambda + N_x(0).
\end{equation}
As our domain is discretized uniformly, we drop the constant
$\Delta x$, and absorb it into our parameter $\beta$ of
-\eqref{eq:first_anisotropic_functional}. Note that since
+\eqref{eq:final_cont_rep}. Note that since
our image takes values in $\mathcal{P} = \{0, \hdots, L-1\}$, the
thresholded image $u^{L-1}$ is equal to zero everywhere.
\end{equation}
As with the $\Delta x$ difference, we can absorb the $\Delta \lambda$
difference into the $\beta$ parameter of
-\eqref{eq:first_anisotropic_functional}. The perimeter is then
+\eqref{eq:final_cont_rep}. The perimeter is then
calculated using a discretized version of the Cauchy--Crofton formula
introduced in Theorem~\ref{thm:riemannian_cauchy_crofton}. Again, we
stop the sum at $L-2$ since the level set $\{ u > L - 1\}$ is empty and
\label{eq:cauchy_crofton_approx1}
\end{aligned}
\end{equation}
-The set of lines $\mathcal{L}$ has been discretized to the lines
+The set of lines $\mathcal{L}$ has been discretized into the set
$\mathcal{L}_D$. Note that we are approximating the length of the
\emph{differentiable} curve $C$ in $\Omega$. Being a difference in the
$\rho$ parameter of our line discretization in
Figure~\ref{fig:line_param}, the difference $\Delta \rho$ represents the
distance from one line to the next in a line family as shown in
-Figure~\ref{fig:line_family}. The difference $\Delta \phi$ is taken to
+Figure~\ref{fig:line_family}, and thus depends on the angle $\phi$
+considered. The difference $\Delta \phi$ is taken to
be the average of the distance to the two neighboring line families as
-shown in Figure~\ref{fig:line_neigh}
+shown in Figure~\ref{fig:line_neigh}, and thus also depends on $\phi$.
-In the discrete setting our domain $\Omega$ is discretized as a regular
+In the discrete setting our domain $\Omega$ is discretized as a uniform
grid $\mathcal{G}$. Our image is then reduced to a function $u :
\mathcal{G} \to \mathcal{P}$.
Moreover, the level sets $\{ u > \lambda\}$ will be functions taking the
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
+made up of 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
\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.
+ \label{eq:cauchy_crofton_approx2}
\end{equation}
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
+later, and thus for an edge $e$ we will only consider the question of
``did $e$ cross $C$ or not?'' This amounts to checking whether the
terminals of $e$ lie on different sides of the curve $C$, and the
approximation is exact for zero or one intersection points, but will, as
-we see in Figure~\ref{fig:curve_edge}, not be entirely correct when we
-have more.
+we see in Figure~\ref{fig:curve_edge}, be wrong when we have more.
\begin{figure}
\input{fig/curve_edge}
\emph{integration scale} $\rho$, and thus we expect the tensors $M(a)$
and $M(b)$ to be similar for edges $e$ of reasonably short length.
-We also remark that using the min-max theorem, it is easy to conclude
-that the eigenvalues of the tensor approximation $M(e_{ab})$ is bounded
+We also remark that using the Rayleigh principle, it is easy to conclude
+that the eigenvalues of the tensor approximation $M(e_{ab})$ are bounded
below and above by the smallest and largest eigenvalues of $M(a)$ and
-$M(b)$. \fixme{Reference to min-max/Rayleigh}
+$M(b)$.
\begin{figure}
\input{fig/area_proof}
will be provided by the following lemma.
\begin{lemma}
For each family of lines given by an angle parameter $\phi$ in the
- regular grid of size $\delta$ we have the relation
+ uniform grid of size $\delta$ we have the relation
\begin{equation}
\delta^2 = \norm{e} \Delta \rho.
\end{equation}
\mathbb{Z}^2$ in the grid.
The distance $\Delta \rho$ from this line $\ell$ to the neighboring
lines can then be calculated as a minimum over the distance to all
- other lines.
+ other grid points.
The lines are split into edges $e = (\delta s, \delta t)^T$ where
$s, t \in \mathbb{Z}$ are coprime such that $e$ does not intersect
- any other points than its two endpoints.
+ any other grid points than its two endpoints.
- Let $(p\prime, q\prime)$ be an arbitrary point not on the line
- $\ell$. We can then calculate
+ We then calculate the minimal distance to a grid point not on the
+ line $\ell$ as
\begin{equation}
\begin{aligned}
- \Delta \rho &= \min_{(p\prime, q\prime)} \left\{
+ \Delta \rho &= \min_{(p\prime, q\prime) \in \mathcal{G}
+ \setminus \ell} \left\{
\left\langle \delta [p - p\prime, q - q\prime],
\frac{e^\perp}{\norm{e^\perp}} \right\rangle \right\} \\
- &= \min_{(p\prime, q\prime)} \left\{\delta^2 \cdot
- \frac{t(p-p\prime) - s(q - q\prime)}{\norm{e}} \right\}.
+ &= \min_{(p\prime, q\prime) \in \mathcal{G} \setminus \ell}
+ \left\{\delta^2 \cdot \frac{t(p-p\prime) - s(q -
+ q\prime)}{\norm{e}} \right\}.
\end{aligned}
\end{equation}
Since $s$ and $t$ are coprime, there exists $a, b \in \mathbb{Z}$
- such that $at - bs = 1$, and since $(p\prime, q\prime)$ is any point
- not on the same line as $(p, q)$ we obtain
+ such that $at - bs = 1$, and since the $\Delta \rho$ cannot be zero,
+ we obtain
\begin{equation}
\Delta \rho = \frac{\delta^2}{\norm{e}}.
\end{equation}
\end{proof}
Inserting $\Delta \rho = \delta^2 / \norm{e}$ and the tensor
approximation of \eqref{eq:tensor_approx} into the curve length
-approximation of \eqref{eq:cauchy_crofton_approx1} we obtain
+approximation of \eqref{eq:cauchy_crofton_approx2} we obtain
\begin{equation}
\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}}},
- \label{eq:cauchy_crofton_approx2}
+ \label{eq:cauchy_crofton_approx3}
\end{equation}
-where the sum is over all edges crossing the curve an odd number of
-times.
+where the sum is over all edges crossing the curve.
The curve length we initially wanted to calculate was the perimeter
-$\PerA(\{u > \lambda\}; \Omega)$ in \eqref{eq:per_approx1}. As this
-curve is closed, we know that every edge intersecting it must have one
-terminal in $\{ u > \lambda \}$ and the other outside. Thus we rewrite
-the sum over $e \cap C$ such that
+$\PerA(\{u > \lambda\}; \Omega)$ in \eqref{eq:per_approx1}. To find
+edges that crosses the curve, we identify the edges that have one
+terminal inside the level set, and the other outside.
+Thus we rewrite the sum over $e \cap C$ such that
\begin{equation}
\PerA(\{u > \lambda\}; \Omega) \approx \sum_{e_{ab}}
\abs{u^\lambda_a - u^\lambda_b} \frac{\det M(e_{ab})
\cdot M(e_{ab}) \cdot e_{ab}\right)^{\sfrac{3}{2}}}.
\label{eq:per_approx2}
\end{equation}
-The absolute value $\abs{u^\lambda_a - u^\lambda_b}$ is one if one of
+The absolute value $\abs{u^\lambda_a - u^\lambda_b}$ is 1 if one of
$a$ and $b$ lie inside the level set and the other lies outside, and
-zero otherwise. In other words the absolute value is one if $e_{ab}$
-crosses the perimeter of $\{ u > \lambda\}$ an odd number of times, and
-zero otherwise.
+0 otherwise.
+
+In other words the absolute value is one if $e_{ab}$ crosses the
+perimeter of $\{ u > \lambda\}$ an odd number of times, and zero
+otherwise.
Thus we have arrived at our final discretization, which takes the form
\begin{gather}
F(u) = \sum_\lambda \sum_x F_x^\lambda(u_x^\lambda) + \beta \sum_\lambda
- \sum_{(x, y)} F_{x,y}^\lambda(u_x^\lambda, u_y^\lambda) :=
+ \sum_{(x, y)} F_{x,y}^\lambda(u_x^\lambda, u_y^\lambda) =:
F^\lambda(u^\lambda), \\
\begin{aligned}
F_x^\lambda(u_x^\lambda) &= \big(N_x(\lambda + 1) -
M(e_{xy}) \cdot e_{xy} \right)^{\sfrac{3}{2}}}.
\end{aligned}
\end{gather}
+Recall that $N_x(\lambda) = \abs{ \lambda - f_x }^2$.
If we minimize $F_\lambda$ to obtain $u^\lambda$ for each level
separately, it is obvious that we will also minimize the sum over all
$F_\lambda$. However, it is not guaranteed that the obtained thresholded
images $u^\lambda$ can be combined to make an output image $u$. They
were defined as $u^\lambda = \idfun_{u > \lambda}$, so we need them to
-be monotonically decreasing (?) in increasing level values, i.e.\
+be monotonically decreasing \fixme{(?)} in increasing level values, i.e.\
\begin{equation}
u_x^\lambda \geq u_x^\mu, \quad \forall \lambda \leq \mu, \quad
\forall x \in \mathcal{G}.
\fixme{Rating: 6/10}
-Consistency relates to how well a solution to the continuous problem
+Consistency relates to whether a solution to the continuous problem
fits in the discretized equation, in other words, whether the
discretized equation approximates the continuous one.
similar to the one in Figure~\ref{fig:line_neigh}.
If we consider the edges $e$ of each family separately, the curve length
-approximation in \eqref{eq:cauchy_crofton_approx2} can be written
+approximation in \eqref{eq:cauchy_crofton_approx3} can be written
\begin{equation}
\begin{aligned}
\abs{C}_M &=
The discretization in the $\phi$ dimension can also be regarded as a
version of the \emph{rectangle method}, although not the midpoint rule.
-As shown in Figure~\ref{fig:circ_rule}, the summand is evaluated on the
-endpoint of the partition intervals $[\phi_k, \phi_{k+1}]$ and the
-difference is taken to be $\Delta \phi = \phi_{k+1} - \phi_k$. Thus if
-$\sup \Delta \phi \to 0$, this discretization is also consistent.
+As shown in Figure~\ref{fig:circ_rule}, the circle is split into
+intervals
+\begin{equation}
+ \left[
+ \frac{\phi_{k-1} + \phi_{k}}{2}, \frac{\phi_{k} + \phi_{k+1}}{2}
+ \right]
+\end{equation}
+of length $\Delta \phi_k = (\phi_{k+1} + \phi_{k-1}) / 2$. The summand
+is evaluated at $\phi_k$, somewhere inside the interval.
+Thus if
+$\sup \Delta \phi_k \to 0$, this discretization is also consistent.
\begin{figure}
\input{fig/circ_rule}
\input{fig/square_cons}
\end{figure}
-For each grid point along the outer edge of this square, we include in
-our neighborhood a grid point having the same angle $\phi$ to the
-$x$-axis. This means either including the actual grid point at the
-outer edge, or one having the same angle, just closer to the center.
+We include all grid points inside the square in our neighborhood, except
+for multiple points that lie on the same line from the origin. Thus if
+two or more grid points lie on the same line, we include only the one
+closest to the origin. This implies that for each grid point along the
+outer edge of this square, we include in our neighborhood a grid point
+having the same angle $\phi$ to the $x$-axis.
+
This construction can be seen in Figure~\ref{fig:square_cons} for $n =
-5$. The maximal $\Delta \phi$ will then be between the horizontal or
-vertical edge and its neighbors, shown in Figure~\ref{fig:square_cons}
-as angle $a$. These angles can be calculated to be
+5$.
+The maximal angle between two lines $\phi_k - \phi_{k-1}$ will always be
+when one of $\phi_k$ or $\phi_{k-1}$ is horizontal or vertical, shown in
+Figure~\ref{fig:square_cons} as the angle $a$.
+Thus the
+largest $\Delta \phi_k$ will then be around the horizontal or
+vertical edge.
+These angles can be calculated to be
\begin{equation}
- \sup \Delta \phi = \arctan \frac{1/n}{n/2} = \arctan
+ \sup \Delta \phi_k = \frac{2}{2} \arctan \frac{1/n}{n/2} = \arctan
\frac{2}{n^2} \to 0.
\end{equation}
found and added to $f$, the distance labels are recalculated, and the
label of the sink will be increased by at least one.
+
+\fixme{it's own section now}
Boykov and Kolmogorov present a variant of the augmenting path algorithm
in \cite{boykov2004experimental}, specialized for the kinds of graphs
occuring in graphical applications. It keeps two non-overlapping trees
when a gap occurs. This is done by keeping track of how many vertices
exist with each label.
+\fixme{move this since we also implemented boykov--kolmogorov}
+
\section{Boykov--Kolmogorov algorithm}
+\begin{algorithm}
+ \begin{algorithmic}
+ \Function{BoykovKolmogorov}{$G$, $s$, $t$}
+ \State $A := {s, t}$, $O := \emptyset$, $S = {s}$, $T = {t}$
+ \State $s\text{.color} = S$, $t\text{.color} = T$
+ \While{\texttt{True}}
+ \State $e \gets$ \Call{Grow}{$G$, $A$}
+ \If{not $e$}
+ \State break
+ \EndIf
+ \State \Call{Augment}{$G$, $s$, $t$, $e$, $O$}
+ \State \Call{Adopt}{$G$, $s$, $t$, $O$}
+ \EndWhile
+ \EndFunction
+ \end{algorithmic}
+ \caption{\sf The Boykov--Kolmogorov maximum flow algorithm}
+ \label{alg:boykov_kolmogorov}
+\end{algorithm}
+
+A maximum flow algorithm specialized for the type of graphs found in
+imaging applications is described by Boykov and Kolmogorov in
+\cite{boykov2004experimental}. Their algorithm is an augmenting path
+algorithm where the paths are found using trees that keep track of
+possible paths.
+
+There is one tree denoted $S$ with the source $s$ as its root, and
+another tree denoted $T$ with the sink $t$ as its root. These trees are
+non-overlapping, and all edges in the tree $S$ can carry flow towards
+the leafs of the tree, while all edges in $T$ can carry flow towards the
+sink $t$.
+
+A vertex can either be part of these trees, or be a free vertex, and we
+write
+\begin{equation}
+ u\text{.color} = \begin{cases}
+ S & \text{ if } u \in S \\
+ T & \text{ if } u \in T \\
+ \texttt{Free} & \text{ otherwise.}
+ \end{cases}
+\end{equation}
+The idea of the algorithm is then to grow these trees towards eachother,
+and when they touch, we have found an augmenting path from $s$ to $t$.
+As much flow as possible is then pushed along this path, and some edges
+will become saturated. If an edge in the tree $S$ or $T$ becomes
+saturated, the tree is split up into multiple trees, also called a
+forest. The parts no longer connected to the root are called orphans,
+and we write $u \in O$ for an orphan vertex $u$. These vertices are---if
+possible---joined back into the tree, or become free nodes.
+
+The vertices in the trees $S$ and $T$ are either active or passive, and
+we write $u \in A$ if $u$ is active. The active vertices are those at
+the boundary of the tree, which can possibly connect to other vertices
+to grow the tree. The passive vertices are internal in the tree, and
+edges to their neighbors are either a part of the tree, or completely
+saturated.
+
+The main loop of the algorithm is then as shown in
+Algorithm~\ref{alg:boykov_kolmogorov}. In the grow procedure, the trees are grown
+from their active vertices until a new path is found. The flow
+is then augmented along this path in the augment procedure, and orphan
+vertices might be created. These orphan vertices are then either joined
+back into their respective trees, or become free, in the adopt
+procedure.
+
+\begin{algorithm}
+\begin{algorithmic}
+ \Function{Grow}{$G$, $A$}
+ \While{$\abs{A} \neq 0$}
+ \State $u \gets $ one node from $A$
+ \ForAll{$v$ such that \Call{treeCap}{$u$, $v$} > 0}
+ \If{$v\text{.color} = \texttt{Free}$}
+ \State $v\text{.color} \gets u\text{.color}$
+ \State $v\text{.parent} \gets u$
+ \State $A := A \cup \{v\}$
+ \ElsIf{$v\text{.color} \neq u\text{.color}$}
+ \State \Return $e := (u, v)$
+ \EndIf
+ \EndFor
+ \State Remove $u$ from $A$
+ \EndWhile
+ \State \Return \texttt{False}
+ \EndFunction
+\end{algorithmic}
+\caption{\sf The grow procedure of the Boykov--Kolmogorov maximum flow
+algorithm}
+\label{alg:grow}
+\end{algorithm}
+
+The grow procedure is shown in Algorithm~\ref{alg:grow}. It goes through
+the set of active vertices $A$ to try to expand the trees $S$ and $T$.
+When considering an active vertex $u \in S$, we want to grow by finding
+vertices $v$ such that $c_f(u, v) > 0$, while when considering an active
+vertex $u \in T$, we want to grow by finding vertices $v$ such that
+$c_f(v, u) > 0$. This is why the \textsc{treeCap} is introduced which is
+defined as
+\begin{equation}
+ \textsc{treeCap}(u, v) = \begin{cases}
+ c_f(u, v) & \text{ if } u \in S, \\
+ c_f(v, u) & \text{ if } u \in T.
+ \end{cases}
+\end{equation}
+Thus if a non-saturated edge is found from an active vertex $u$ to a free
+vertex $v$, then $v$ is added to the tree of $u$. The tree structure is
+stored by keeping a parent attribute in each non-free node, such that
+for example $v\text{.parent} = u$.
+
+If a non-saturated edge is found from an active vertex $u$ to a vertex
+$v$ in the other tree, the two trees connect, and we can return the edge
+$e$ since we have an augmenting path. Note that if this happens, the
+vertex $u$ is still active. It only becomes passive when all neighbors
+are considered without yielding an augmenting path.
+
+\begin{algorithm}
+ \begin{algorithmic}
+ \Function{Augment}{$G$, $s$, $t$, $e$, $O$}
+ \State $p \gets $ path from $s$ to $t$ through $e$
+ \Comment through the trees $S$ and $T$
+ \State $\Delta f \gets \min \{c_f(u, v) : (u,v) \in p\}$
+ \ForAll{$(u, v) \in p$}
+ \State $f(u, v) \mathrel{+}= \Delta f$
+ \State $f(v, u) \mathrel{-}= \Delta f$
+ \If{$c_f(u, v) = 0$}
+ \If{$u\text{.color} = v\text{.color} = S$}
+ \State $v\text{.parent} = \texttt{Null}$
+ \State $O := O \cup \{v\}$
+ \ElsIf{$u\text{.color} = v\text{.color} = T$}
+ \State $u\text{.parent} = \texttt{Null}$
+ \State $O := O \cup \{u\}$
+ \EndIf
+ \EndIf
+ \EndFor
+ \EndFunction
+ \end{algorithmic}
+ \caption{\sf The augment procedure of the Boykov--Kolmogorov maximum
+ flow algorithm}
+ \label{alg:augment}
+\end{algorithm}
+
+The augment procedure finds the augmenting path going through the tree
+$S$, the edge $e$ and the tree $T$. The maximal possible flow is then
+sent along this path. At least one edge will then become saturated. If a
+saturated edge occurs in the trees $S$ or $T$, the edge terminal
+farthest from the tree root is marked as an orphan. Note that even if a
+whole subtree is disconnected from the main tree, only the root of this
+subtree is marked as an orphan. The adopt procedure will take care of
+either reconnecting, or freeing all the vertices in the disconnected
+subtree.
+
+\begin{algorithm}
+\begin{algorithmic}
+ \Function{Adopt}{$G$, $s$, $t$, $O$}
+ \While{$\abs{O} \neq 0$}
+ \State $u \gets $ one node from $O$
+ \State Remove $u$ from $O$
+ \State found $\gets \texttt{False}$
+ \ForAll{$v$ such that \Call{treeCap}{$v$, $u$} > 0}
+ \If{$u\text{.color} \neq v\text{.color}$}
+ \State continue
+ \EndIf
+ \State orig $\gets$ \Call{treeOrigin}{$v$}
+ \If{orig $\neq s$ and orig $\neq t$}
+ \State continue
+ \EndIf
+ \State found $\gets \texttt{True}$
+ \State $u\text{.parent} \gets v$
+ \State break
+ \EndFor
+ \If{found $\neq \texttt{True}$}
+ \ForAll{$v$ such that $u\text{.color} = v\text{.color}$}
+ \If{\Call{treeCap}{$v$, $u$} > 0}
+ \State $A := A \cup \{ v \}$
+ \EndIf
+ \If{$v\text{.parent} = u$}
+ \State $O := O \cup \{ v \}$
+ \State $v\text{.parent} \gets \texttt{Null}$
+ \EndIf
+ \State $u\text{.color} \gets \texttt{Free}$
+ \State $A := A - \{ u \}$
+ \EndFor
+ \EndIf
+ \EndWhile
+ \EndFunction
+\end{algorithmic}
+\caption{\sf The adopt procedure of the Boykov--Kolmogorov maximum flow
+algorithm}
+\label{alg:adopt}
+\end{algorithm}
+
+The adopt procedure processes all the vertices in the set of orphans
+$O$. These vertices are either single vertices, or roots of disconnected
+subtrees. For an orphaned vertex $u$ we look through its neighbors to
+find a possible parent vertex $v$. It should belong to the same tree as
+$u$, and also satisfy $\textsc{treeCap}(v, u) > 0$. In addition, $v$
+should be connected to one of the tree roots $s$ or $t$. This keeps us
+from connecting to other orphaned vertices or subtrees, and is checked
+by the procedure $\textsc{treeOrigin}(v)$, which follows the parent
+information until reaching either $s$, $t$ or an orphaned vertex.
+
+If a potential parent $v$ of $u$ is \emph{not} found, $u$ becomes a free
+vertex. All vertices that had $u$ as its parent are orphaned, and are thus
+treeted by the adopt procedure later. Vertices that are in one of the
+trees $S$ or $T$, and have a non-saturated edge to this newly freed
+vertex $u$, become active.
+
+When no more active vertices remain, there are no more possible
+augmenting paths, and the algorithms terminates with a maximum flow.
+A proof of correctness can be found in \fixme{kolmogorov-phd}.
+
+\subsection{Performance improvements}
+
+Adopt to vertex closest to root. And FIFO. Do two-edge pushes initially.
+
+\subsection{Graph reuse}
+
+As in the push-relabel algorithm described earlier, it is not necessary
+to completely restart the Boykov--Kolmogorov algorithm for every level
+of the image, but we can reuse the graph and the trees $S$ and $T$ in
+successive runs.
+
+Since this is algorithm does not deal with the relaxed preflow concept,
+the updates from one level to the next has to be done in a different
+way. We can no longer decrease the capacity of edges, as this could
+break the flow conservation constraint.
+
\section{Performance comparison}
Maybe.
projects / master / commitdiff