From: Bjørn Rustad Date: Mon, 19 Jan 2015 14:48:55 +0000 (+0100) Subject: Fix figure, theory, table and BK X-Git-Url: http://git.rustad.me/?a=commitdiff_plain;h=c9baaf970b49ed88493e783e279932238de79673;p=master Fix figure, theory, table and BK --- diff --git a/commands.tex b/commands.tex index 3ef03af..2a7a689 100644 --- a/commands.tex +++ b/commands.tex @@ -13,6 +13,7 @@ \newcommand{\Ahalf}{A^{\sfrac{1}{2}}} \newcommand{\Mhalf}{M^{\sfrac{1}{2}}} \newcommand{\rom}[1]{\uppercase\expandafter{\romannumeral #1\relax}} +\AtBeginDocument{\renewcommand{\setminus}{\mathbin{\backslash}}} \definecolor{fixme-color}{rgb}{0.7,0.2,0} \newcommand{\fixme}[1]{\textcolor{fixme-color}{FIXME: #1}} \def\cpp{C{}\texttt{++}} diff --git a/fig/circ_rule.tex b/fig/circ_rule.tex index d1ce2bf..db1df7c 100644 --- a/fig/circ_rule.tex +++ b/fig/circ_rule.tex @@ -1,5 +1,5 @@ \centering -\begin{tikzpicture}[scale=1.0] +\begin{tikzpicture}[scale=1.3] \node[tiny vertex] (o) at (1, 1) {}; @@ -10,17 +10,19 @@ \draw (1,1) +(0:1.5cm) arc (0:360:1.5cm); %\path (2,2) ++(22.5:3cm) node {$\phi_{xy}$}; - \draw[dash] (1,1) +(45:0cm) -- +(45:2.6cm); - \draw[dash] (1,1) +(63.3:0cm) -- +(63.43:2.6cm); + \draw[dash] (1,1) +(35.785:0cm) -- +(35.785:2.6cm); + \draw[dash] (1,1) +(54.215:0cm) -- +(54.215:2.6cm); + \draw (1,1) +(26.57:1.4cm) -- +(26.57:1.6cm) +(23:2.1cm) + node[font=\small, fill=white, inner sep=0.5pt] {$\phi_{k-1}$}; \draw (1,1) +(45:1.4cm) -- +(45:1.6cm) +(45:1.9cm) node[font=\small, fill=white, inner sep=0.5pt] {$\phi_k$}; - \draw (1,1) +(63.43:1.4cm) -- +(63.43:1.6cm) +(67.5:1.9cm) + \draw (1,1) +(63.43:1.4cm) -- +(63.43:1.6cm) +(67.4:1.9cm) node[font=\small, fill=white, inner sep=0.5pt] {$\phi_{k+1}$}; %\path (2,2) ++(67.5:3.8cm) node{$\Delta\phi_{xy}$}; - \draw (1,1) +(45:2.5cm) arc (45:63.43:2.5cm); - \path (1,1) ++(54.215:3cm) node {$\Delta \phi_k$}; + \draw (1,1) +(35.785:2.5cm) arc (35.785:54.215:2.5cm); + \path (1,1) ++(45:2.9cm) node {$\Delta \phi_k$}; \foreach \p in {0,90,180,270} { \draw (1,1) +(\p:1.4cm) -- +(\p:1.6cm); @@ -30,7 +32,7 @@ } \end{tikzpicture} -\caption{ +\caption{% We showed that the angular difference between subsequent angle parameters goes to zero. The discretization in the $\phi$ dimension can be viewed as a rectangle approximation rule of the integral, as diff --git a/main.tex b/main.tex index 9645236..67be508 100644 --- a/main.tex +++ b/main.tex @@ -6,10 +6,10 @@ \usepackage{mathtools} \usepackage{unicode-math} \usepackage{fontspec} -%\setmainfont[Ligatures=TeX]{Asana Math} +%\setmainfont[Mapping=tex-text,Ligatures=TeX]{Asana} %\setmainfont[Ligatures=TeX]{TG Termes Math} -\setmainfont[Mapping=tex-text,Ligatures=TeX]{TeX Gyre Pagella} -\setmathfont[Mapping=tex-text,Ligatures=TeX]{TeX Gyre Pagella Math} +%\setmainfont[Mapping=tex-text,Ligatures=TeX]{TeX Gyre Pagella} +%\setmathfont[Mapping=tex-text,Ligatures=TeX]{TeX Gyre Pagella Math} \usepackage{amsthm} \usepackage{graphicx} \usepackage[percent]{overpic} diff --git a/results.tex b/results.tex index 00c6317..71c08d5 100644 --- a/results.tex +++ b/results.tex @@ -390,6 +390,29 @@ overestimated as the sum of their three components. \input{fig/pixel_perimeter.tex} \end{figure} +\begin{table} + \caption{% + The circumference of different circles measured by the + discretized Cauchy--Crofton formula in \eqref{ref}, using different + neighborhood stencils. + } + \centering + \begin{tabular}{cccc} + \hline + $\abs{\mathcal{N}}$ & $r = 0.5$ & $r = 5.5$ & $r = 50.5$ \\ + \hline + 4 & 3.14 & 34.56 & 317.3 \\ + 8 & 2.68 & 33.94 & 317.5 \\ + 16 & 2.08 & 34.59 & 316.8 \\ + 32 & 1.63 & 34.44 & 317.2 \\ + 48 & 1.40 & 34.29 & 317.3 \\ + 72 & 1.21 & 33.95 & 317.2 \\ + \hline + $2\pi r$ & 3.14 & 34.56 & 317.3 \\ + \hline + \end{tabular} +\end{table} + As mentioned before, there is another discretization error which relates to the length of the edges in the neighborhood. Thus a larger neighborhood is not always better, even if it will reduced the artifacts diff --git a/theory.tex b/theory.tex index 7820bdf..1e3f2d1 100644 --- a/theory.tex +++ b/theory.tex @@ -212,6 +212,7 @@ U(x)^T$ where \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} @@ -681,10 +682,10 @@ calculated in the following way 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} @@ -766,11 +767,6 @@ and the distance to origin $\rho$.}% \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 @@ -786,8 +782,13 @@ the metric tensor $M$.}% \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 @@ -802,8 +803,8 @@ where our domain is equipped with a metric tensor in each point. \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} @@ -815,9 +816,9 @@ where our domain is equipped with a metric tensor in each point. &= \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 @@ -825,7 +826,7 @@ where our domain is equipped with a metric tensor in each point. 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), \\ @@ -835,13 +836,18 @@ where our domain is equipped with a metric tensor in each point. \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 \\ @@ -863,10 +869,10 @@ where our domain is equipped with a metric tensor in each point. \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 @@ -889,7 +895,8 @@ where our domain is equipped with a metric tensor in each point. \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} @@ -913,7 +920,8 @@ where our domain is equipped with a metric tensor in each point. = \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})} @@ -929,7 +937,7 @@ where our domain is equipped with a metric tensor in each point. 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 @@ -944,14 +952,14 @@ where our domain is equipped with a metric tensor in each point. 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}}}. @@ -971,18 +979,20 @@ where our domain is equipped with a metric tensor in each point. \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} @@ -994,11 +1004,14 @@ where our domain is equipped with a metric tensor in each point. 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 @@ -1007,14 +1020,12 @@ where our domain is equipped with a metric tensor in each point. 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} @@ -1029,11 +1040,11 @@ where our domain is equipped with a metric tensor in each point. \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. @@ -1047,8 +1058,9 @@ where our domain is equipped with a metric tensor in each point. \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} @@ -1072,11 +1084,10 @@ 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. 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 @@ -1092,9 +1103,7 @@ from an arc length parametrization of the boundary. If $P$ is a \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 }. @@ -1116,16 +1125,17 @@ Section~\ref{sec:anisotropy_tensor} \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 + @@ -1133,55 +1143,80 @@ our anisotropy tensor $A(x)$ is bounded away from zero as \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} @@ -1193,13 +1228,13 @@ the fidelity term can be discretized as in my project work 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} @@ -1218,14 +1253,14 @@ where $I(x)$ is the indicator function that takes the value 1 if $x$ is 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. @@ -1244,7 +1279,7 @@ the discrete levels to get \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 @@ -1276,17 +1311,18 @@ the approximation \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 @@ -1309,7 +1345,7 @@ 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 +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 @@ -1318,16 +1354,16 @@ discretization $\mathcal{L}_D$ to obtain \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} @@ -1349,10 +1385,10 @@ structure tensor in \eqref{eq:s_def} corresponding to the \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} @@ -1363,7 +1399,7 @@ we need a way to calculate the inter-line distance $\Delta \rho$ which 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} @@ -1374,47 +1410,48 @@ will be provided by the following lemma. \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}) @@ -1422,16 +1459,18 @@ the sum over $e \cap C$ such that \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) - @@ -1442,13 +1481,14 @@ Thus we have arrived at our final discretization, which takes the form 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}. @@ -1463,7 +1503,7 @@ this requirement.} \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. @@ -1482,7 +1522,7 @@ convenience we will use a neighborhood representation of $\mathcal{L}_D$ 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 &= @@ -1523,10 +1563,17 @@ part of the discretization is fine. 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} @@ -1545,16 +1592,24 @@ of this square $n$ goes to infinity. \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} @@ -2064,6 +2119,8 @@ to $t$ following edges in $E_f'$. After the blocking flow $f'$ has been 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 @@ -2625,8 +2682,234 @@ For the gap relabeling heuristic, we need to have a easy way of finding 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.