From: Bjørn Rustad Date: Wed, 10 Dec 2014 10:30:22 +0000 (+0100) Subject: Writing and fixin X-Git-Url: http://git.rustad.me/?a=commitdiff_plain;h=7cc70d3f295ed56fddf51473378b8380cd6b2937;p=master Writing and fixin --- diff --git a/introduction.tex b/introduction.tex index 4642b62..aa0d2b1 100644 --- a/introduction.tex +++ b/introduction.tex @@ -115,19 +115,20 @@ The diffusivity is made directionally dependent by introducing a diffusion \emph{tensor} $A(u)$ such that the initial boundary value problem becomes \begin{equation} - \begin{cases} - \partial_t u &= \diver \big(A(u) \nabla u\big) - \text{ on } \Omega \times (0, \infty),\\ - u|_{t=0} &= f - \text{ on } \Omega,\\ - A(u) \nabla u \cdot \nu &= 0 - \text{ on } \partial \Omega \times (0, \infty). - \end{cases} + \left\{ + \begin{alignedat}{10} + &\partial_t u &{}={}& \diver \big(A(u) \nabla u\big) + &\quad \text{ on }& \Omega \times (0, \infty),\\ + &u|_{t=0} &{}={}& f + &\text{ on }& \Omega,\\ + &A(u) \nabla u \cdot \nu &{}={}& 0 + &\text{ on }& \partial \Omega \times (0, \infty), + \end{alignedat} + \right. \label{eq:aniso_diff} \end{equation} \nomenclature{$A(u)$}{Thermal diffusivity tensor, or anisotropy tensor.}% -\fixme{triple column alignment?} where $\nu$ is the outer normal of $\Omega$. The tensor $A(u)$ would normally be constructed such as to diminish the effect of $\nabla u$ across what we believe to be edges in the image. This way, there will diff --git a/results.tex b/results.tex index 3c51a53..5d32c5e 100644 --- a/results.tex +++ b/results.tex @@ -214,7 +214,7 @@ are much higher for the anisotropic total variation. \begin{figure} \centering - \input{fig/contrast_plot} + \input{fig/factory/finger/contrast_plot} \caption{Showing some contrast loss.} \label{fig:contrast_plot} \end{figure} @@ -317,19 +317,27 @@ overestimated as the sum of their three components. \begin{figure} \centering - \begin{subfigure}[b]{0.40\textwidth} + \begin{subfigure}[t]{0.30\textwidth} \centering - \includegraphics[width=\textwidth]{fig/circle100.png} + \includegraphics[width=\textwidth]{fig/factory/circle_deform/n_q0.png} \caption{ Circle. } \end{subfigure} ~ - \begin{subfigure}[b]{0.40\textwidth} + \begin{subfigure}[t]{0.30\textwidth} + \centering + \includegraphics[width=\textwidth]{fig/factory/circle_deform/r_p2_n8_b100000_g10000000000_r1_s1.png} + \caption{ + Size of neighborhood: 8. + } + \end{subfigure} + ~ + \begin{subfigure}[t]{0.30\textwidth} \centering - \includegraphics[width=\textwidth]{fig/neigh_artifacts.png} + \includegraphics[width=\textwidth]{fig/factory/circle_deform/r_p2_n32_b100000_g10000000000_r1_s1.png} \caption{ - Artifacts, and contrast loss. + Size of neighborhood: 32. } \end{subfigure} \caption{ diff --git a/theory.tex b/theory.tex index eb5d6b3..cc3cf49 100644 --- a/theory.tex +++ b/theory.tex @@ -13,8 +13,8 @@ important theorems, transform it into a form that is easy to discretize. \section{Anisotropic total variation} The method considered will build on the total variation regularization -method of Section \ref{sec:total_variation}. From anisotropic diffusion -in Section \ref{sec:anisotropic_diffusion} we borrow the idea of making +method of Section~\ref{sec:total_variation}. From anisotropic diffusion +in Section~\ref{sec:anisotropic_diffusion} we borrow the idea of making the regularization in each point directionally dependent. We introduce the anisotropic total variation \begin{equation} @@ -22,12 +22,12 @@ the anisotropic total variation \label{eq:aniso_tv_sqrt} \end{equation} for all $u \in C^1(\Omega)$. If $A(x)$ is the identity -matrix we get the regular total variation found in Definition -\ref{def:tv}. One problem with the regular total -variation method is that when reducing the total variation, it will also -try to reduce the variation over known edges in the image, which can -lead to contrast loss, especially in thin details. By controlling $A(x)$ -such that $\nabla u(x)$ is weighted down across known edges, we hope to +matrix we get the regular total variation found in +Definition~\ref{def:tv}. One problem with the regular total variation +method is that when reducing the total variation, it will also try to +reduce the variation over known edges in the image, which can lead to +contrast loss, especially in thin details. By controlling $A(x)$ such +that $\nabla u(x)$ is weighted down across known edges, we hope to retain the regularization properties of the method while reducing some of the negative effects. If the variation across an edge is ``ignored'' by the energy functional, there is no gain in reducing the height of the @@ -48,6 +48,7 @@ as $A(x) = \Ahalf(x) \Ahalf(x)$. We can then write &= \sup_{\eta^T A^{-1} \eta \leq 1} \int_\Omega u \diver \eta \, dx, \end{align} + 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^* = @@ -84,7 +85,7 @@ root.} \subsection{Anisotropy tensor} \label{sec:anisotropy_tensor} -\fixme{ +\fixme{% This tensor can be visualized, and can be used for edge detection. There is probably a bit to be said about the smoothing we do. We also discuss whether to use the noisy image (yes, probably), @@ -109,7 +110,7 @@ none, so we somehow want to be sure about the edges we find. Edges can be found in many different ways, but as suggested by Weickert in his book on Anisotropic Diffusion \cite{weickert1998anisotropic}, and -briefly mentioned in Section \ref{sec:anisotropic_diffusion}, a good +briefly mentioned in Section~\ref{sec:anisotropic_diffusion}, a good starting point is the \emph{edge detector} $\nabla \tilde{f}_\sigma$ and the structure tensor \begin{equation} @@ -125,9 +126,9 @@ control how sensitive the edge detector $\nabla \tilde{f}_\sigma$ is to noise and small details. \begin{figure} - \centering + \centering{} \includegraphics[width=0.4\textwidth]{fig/edges.png} - \caption{ + \caption{% The largest eigenvalue of the structure tensor is $\abs{\nabla \tilde{f}_\sigma(x)}^2$ and is as we see a good edge detector. } @@ -139,7 +140,7 @@ First consider the tensor $S_0(x) = \nabla \tilde{f}_\sigma(x) \otimes as the edge detector itself. Its eigenvalues will be 0 and $\labs{\nabla \tilde{f}_\sigma(x)}^2$ with corresponding eigenvectors $v_1$ and $v_2$ perpendicular and parallel to $\nabla \tilde{f}_\sigma(x)$ respectively. -Figure \ref{fig:edges} shows that the largest eigenvalue of the +Figure~\ref{fig:edges} shows that the largest eigenvalue of the structure tensor is a good edge detector. We are also interested in identifying and being sensitive to features in @@ -204,7 +205,7 @@ the norm across edges, we need $\sigma_1 \leq \sigma_2$, so we construct \begin{equation} \begin{aligned} \sigma_1 &= \left(1 + - \frac{(\lambda_1 - \lambda_2)^2}{\gamma^2}\right)^{-1}, \\ + \frac{(\lambda_1 - \lambda_2)^2}{\omega^2}\right)^{-1}, \\ \sigma_2 &= 1. \end{aligned} \label{eq:sigma_construction} @@ -218,7 +219,7 @@ Note that in smooth parts of the image we have $\sigma_1 \approx \sigma_2 = 1$ and the anisotropic total variation is close to the regular total variation. -\fixme{ +\fixme{% The numerical problems should be discussed somewhere but maybe not here. The eigenvalues are extracted such that } @@ -294,7 +295,7 @@ point $u$ means that at points $u_\epsilon$ close to $u$, the functional takes values either close to or above $F(u)$. More specifically, for a sequence $u_k$ converging to $u$, we have $F(u) \leq \liminf_k F(u_k)$. For a function $f : \mathbb{R} \to \mathbb{R}$ this can be visualized as -in Figure \ref{fig:lower_semicont}. +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 @@ -373,7 +374,7 @@ definition from \eqref{eq:extended_tv}, we have %\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\} + C_c^\infty(\Omega, \mathbb{R}^2), \norm{\xi}_A^* \leq 1 \right\} \end{equation} This is again is a supremum of weakly continuous functionals. Thus the regularization term @@ -513,8 +514,8 @@ given in \cite{olsson2009extending}. 1 & \text{if } t > s + r \end{cases} \end{equation} - visualized in Figure \ref{fig:eta_r} such that its derivative takes the - form shown in Figure \ref{fig:eta_r_diff}. By composing the function + visualized in Figure~\ref{fig:eta_r} such that its derivative takes the + form shown in Figure~\ref{fig:eta_r_diff}. By composing the function $\eta_r$ with our image $u$ and Green's identity we obtain \begin{equation} \int_\Omega - \eta_r(u) \diver \xi \, dx @@ -610,7 +611,7 @@ integration formula, and discretization, where an approximation of the perimeter will be computed using a graph cut machinery. \fixme{well, if we assumed differentiability, we wouldn't need the sup -definition of the TV...} +definition of the TV.} \section{Cauchy--Crofton formulas} @@ -627,7 +628,7 @@ counting the times it intersects line in the set of all lines. Intuitively, a long curve will intersect more lines. 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 +and parametrize them as shown in Figure~\ref{fig:line_param}. Thus a line is parametrized by the angle $\phi \in [0, 2\pi)$ of the normal going to the origin, and the distance $\rho \in [0, \infty)$ from origin to the line. Sometimes it is more convenient to consider a unit vector $\nu$ giving the direction of @@ -742,10 +743,11 @@ integration theory see for example \cite{bartle1995elements}. 0 & \sigma_2 \end{pmatrix}. \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(\ell_{\phi, \rho}) = J_{\Sigma^2}(U^T \ell_{\phi, - \rho})$. Given a line + As $U$ and $U^T$ correspond to rotations and our measure + $\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 \begin{equation} \ell_{\phi, \rho} = \begin{pmatrix} \rho\cdot \cos \phi \\ @@ -785,7 +787,7 @@ integration theory see for example \cite{bartle1995elements}. $\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} + \partial_\phi \theta &= \frac{\frac{\sigma_1}{\sigma_2} \sec^2 \phi}{1 + \frac{\sigma_1^2}{\sigma_2^2} \tan^2 \phi} = \frac{\sigma_1 \sigma_2}{\sigma_1^2 \sin^2 \phi + \sigma_2^2 @@ -811,10 +813,10 @@ integration theory see for example \cite{bartle1995elements}. = \frac{\sigma_1^2 \sigma_2^2} {\left(\sigma_1^2 \sin^2 \phi + \sigma_2^2 \cos^2 \phi \right)^{\sfrac{3}{2}}} - = \frac{\sigma_1^2 \sigma_2^2 + = \frac{\sigma_1^2 \sigma_2^2 }{\left(\sigma_1^2 \nu_x^2 + \sigma_2^2 \nu_y^2\right)^{\sfrac{3}{2}}} - = \frac{\det \Sigma^2}{(\nu^T \cdot \Sigma^2 \cdot + = \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 @@ -891,7 +893,7 @@ integration theory see for example \cite{bartle1995elements}. \begin{equation} \abs{C}_{M_\pi} = \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) } + = \int_C \sqrt{\dot{C}(t)^T M_\pi\big( C(t) \big) \dot{C}(t)} \, dt. \label{eq:piecewise_curve_length} \end{equation} @@ -900,7 +902,7 @@ integration theory see for example \cite{bartle1995elements}. $\labs{\dot{C}(t)}_M$. We have assumed bounds on the eigenvalues 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 + 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} @@ -928,7 +930,7 @@ integration theory see for example \cite{bartle1995elements}. function $w_\pi$ is then bounded such that \begin{equation} \sum_{x \in \ell_{\nu, \rho} \cap C} w_\pi(\nu, x) - \leq + \leq \sum_{x \in \ell_{\nu, \rho} \cap C} \frac{K^2}{k^{\sfrac{3}{2}}} = @@ -936,7 +938,7 @@ integration theory see for example \cite{bartle1995elements}. = 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 + formula in Theorem~\ref{thm:euclidean_cauchy_crofton} and the fact that we assumed $C$ to be of finite length \begin{equation} \int_\mathcal{L} g(\ell_{\nu, \rho}) \, d\mathcal{L}(\ell_{\nu, @@ -962,14 +964,14 @@ integration theory see for example \cite{bartle1995elements}. \end{equation} \end{proof} -With the anisotropic coarea formula in Theorem -\ref{thm:anisotropic_coarea} we have a way to calculate the anisotropic -total variation by integrating the anisotropic perimeter of each level -set of the image as introduced in \eqref{eq:functional_per}. In the -Euclidean case, which here would amount to 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. +With the anisotropic coarea formula in +Theorem~\ref{thm:anisotropic_coarea} we have a way to calculate the +anisotropic total variation by integrating the anisotropic perimeter of +each level set of the image as introduced in \eqref{eq:functional_per}. +In the Euclidean case, which here would amount to 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 \eqref{eq:perimeter_calc}, the @@ -981,9 +983,9 @@ are always perpendicular, so all we need is a 90\textdegree{} rotation. If $P$ is a 90\textdegree{} rotation matrix we have \begin{equation} \begin{aligned} - \PerA(U; \Omega) &= \int_{\partial U} \sqrt{ \langle + \PerA(U; \Omega) &= \int_{\partial U} \sqrt{\langle \nu_{\partial U}, A(x) \nu_{\partial U} \rangle} \, dt \\ - &= \int_{\partial U} \sqrt{ \langle P \nu_{\partial U}, P A(x) + &= \int_{\partial U} \sqrt{\langle P \nu_{\partial U}, P A(x) P^T \, P \nu_{\partial U} \rangle } \, dt. %&\stackrel{\mathclap{\gamma \hskip 0.07em = \hskip 0.07em %\partial U}}{=} @@ -998,7 +1000,7 @@ the boundary of $U$. Observe that a 90\textdegree{} rotation of the normal $\nu_{\partial U}$ gives us the tangent $\pm \dot{\gamma}$. Inserting this into \eqref{eq:per_to_length1} we get \begin{equation} - \PerA(U; \Omega) = \int_\gamma \sqrt{ \langle \dot{\gamma}, M(x) \, + \PerA(U; \Omega) = \int_\gamma \sqrt{\langle \dot{\gamma}, M(x) \, \dot{\gamma} \rangle }. \label{eq:per_to_length2} \end{equation} @@ -1007,19 +1009,19 @@ Note that since $\Omega$ is open, we assume that $\gamma = \partial U of $\Omega$. Now we make sure that all the assumptions of the Riemannian -Cauchy--Crofton formula in Theorem \ref{thm:riemannian_cauchy_crofton} +Cauchy--Crofton formula in Theorem~\ref{thm:riemannian_cauchy_crofton} are fulfilled so that it can be applied to the curve length integral we have constructed in \eqref{eq:per_to_length2}. -The structure tensor is constructed as described in Section -\ref{sec:anisotropy_tensor} +The structure tensor is constructed as described in +Section~\ref{sec:anisotropy_tensor} \begin{equation} S_\rho(x) = \left(K_{\rho} * \left( \nabla \tilde{f}_{\sigma} \otimes \nabla \tilde{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 -can see in \eqref{eq:s_eigenvalues}, the eigenvalues +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} @@ -1030,26 +1032,25 @@ are bounded from above and by the construction in our anisotropy tensor $A(x)$ is bounded away from zero as \begin{equation} \lambda_2 = \left(1 - + \frac{(s_1 - s_2)^2}{\gamma^2}\right)^{-1} \geq \left(1 + - \frac{s_1^2}{\gamma^2}\right)^{-1} \geq k > 0. -\end{equation} -\fixme{bah, now gamma is a curve...} + + \frac{(s_1 - s_2)^2}{\omega^2}\right)^{-1} \geq \left(1 + + \frac{s_1^2}{\omega^2}\right)^{-1} \geq k > 0. +\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 +\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 +Theorem~\ref{thm:anisotropic_coarea} obtain our final expression for the continuous anisotropic total variation \begin{equation} \TVA(u) = \int_{-\infty}^{\infty} \int_\mathcal{L} \sum_{x \in - \ell_{\nu, \rho} \cap \gamma_s } \frac{\det M(x)} + \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} -where $\gamma_s = \partial \{ u > s \}$ is the boundary of the level set +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 @@ -1138,7 +1139,7 @@ Note that we will later ignore the $\Delta \lambda$ difference, as we can just absorb it into the $\beta$ parameter of \eqref{eq:first_anisotropic_functional}. The perimeter is then calculated using a discretized version of the Cauchy--Crofton formula -introduced in Theorem \ref{thm:riemannian_cauchy_crofton}. Again, we +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 has zero perimeter. @@ -1147,8 +1148,8 @@ has zero perimeter. \begin{figure} \input{fig/line_disc} \end{figure} -By approximating the integral Theorem -\ref{thm:riemannian_cauchy_crofton} by a discrete sum we obtain +By approximating the integral +Theorem~\ref{thm:riemannian_cauchy_crofton} by a discrete sum we obtain the approximation \begin{equation} \begin{aligned} @@ -1171,12 +1172,12 @@ the approximation 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$. 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 be the -average of the distance to the two neighboring line families as shown in -Figure \ref{fig:line_neigh} +$\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 +be the average of the distance to the two neighboring line families as +shown in Figure~\ref{fig:line_neigh} In the discrete setting our domain $\Omega$ is discretized as a regular grid $\mathcal{G}$. Our image is then reduced to a function $u : @@ -1195,10 +1196,10 @@ possible lines of that family, meaning there are no grid points without a line of that family intersecting it. 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$. +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 @@ -1218,7 +1219,7 @@ 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 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 +we see in Figure~\ref{fig:curve_edge}, not be entirely correct when we have more. \begin{figure} @@ -1308,10 +1309,10 @@ 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 \begin{equation} - \PerA(\{u > \lambda\}; \Omega) \approx \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}}}. + \PerA(\{u > \lambda\}; \Omega) \approx \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}}}. \label{eq:per_approx2} \end{equation} The absolute value $\abs{u^\lambda_a - u^\lambda_b}$ is one if one of @@ -1371,7 +1372,7 @@ consistent with the continuous one. We will show that for an increasingly refined discrete domain $\mathcal{G}$, there exists a choice for $\mathcal{L}_D$ that leads to a consistent Cauchy--Crofton formula. For convenience we will use a neighborhood representation of $\mathcal{L}_D$ -similar to the one in Figure \ref{fig:line_neigh}. +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 @@ -1405,9 +1406,9 @@ in \eqref{eq:tensor_approx} converges to the tensor in the intersection point. For each $\phi$ parameter, our discretization in the $\rho$ dimension -can be regarded as a midpoint rule as shown in Figure -\ref{fig:line_midpoint}. Thus if $\sup \Delta \rho \to 0$, this part of -the discretization is fine. +can be regarded as a midpoint rule as shown in +Figure~\ref{fig:line_midpoint}. Thus if $\sup \Delta \rho \to 0$, this +part of the discretization is fine. \begin{figure} \input{fig/line_midpoint} @@ -1415,7 +1416,7 @@ 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 +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. @@ -1427,7 +1428,7 @@ $\sup \Delta \phi \to 0$, this discretization is also consistent. To show that all these properties can be fulfilled, we look at a particular neighborhood stencil construction. Consider a square centered around a grid point with side lengths -$\sqrt{\delta}$ as shown in Figure \ref{fig:square_cons}. As $\delta$ +$\sqrt{\delta}$ as shown in Figure~\ref{fig:square_cons}. As $\delta$ goes to zero, the size of this square will go to zero. Inside this square we can fit a square of $n^2 = \lfloor 1 / \sqrt{\delta} \rfloor^2$ grid points. This means that the number of grid points along the outer edge @@ -1441,9 +1442,9 @@ 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. -This construction can be seen in Figure \ref{fig:square_cons} for $n = +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} +vertical edge and its neighbors, shown in Figure~\ref{fig:square_cons} as angle $a$. These angles can be calculated to be \begin{equation} \sup \Delta \phi = \arctan \frac{1/n}{n/2} = \arctan @@ -1455,7 +1456,7 @@ diagonal of the square such that \begin{equation} \norm{e} \leq \sqrt{\delta / 2} \to 0. \end{equation} -And finally we know from Lemma \ref{lem:delta_rho} that for each line family +And finally we know from Lemma~\ref{lem:delta_rho} that for each line family $\delta^2 = \Delta \rho \norm{e}$ and the fact that $\norm{e} \geq \delta$. Thus for the inter-line distance $\Delta \rho$ we have \begin{equation} @@ -1468,7 +1469,7 @@ well-known, and consistent integral approximations, where the summand converges to the integrand, and the differences $\Delta \phi$ and $\Delta \rho$ go to zero. Thus the perimeter approximation in \eqref{eq:per_approx1} is consistent with the continuous formulation in -Theorem \ref{thm:riemannian_cauchy_crofton}. +Theorem~\ref{thm:riemannian_cauchy_crofton}. Note that as we will work with digital images with fixed resolutions, we do not really have the chance to refine our discretization. We do @@ -1614,7 +1615,7 @@ depend on. We will construct a graph in such a way that if a variable $u^\lambda_x$ ends up on the source side of the cut we set $u^\lambda_x = 0$, and if it ends up on the sink side we set $u^\lambda_x = 1$, as in -Definition \ref{def:graph_representable}. This is an arbitrary choice, +Definition~\ref{def:graph_representable}. This is an arbitrary choice, but still something we have to keep in mind through the rest of the section. @@ -1629,7 +1630,7 @@ term coming from our aim to minimize the total variation. $u^\lambda_x \in \{0,1\}$. The energy $E_\lambda^x(u^\lambda_x)$ and minimum cut obtaining this configuration is shown. The last two columns show the capacities of the cut for each of the two - graph constructions in Figure \ref{fig:norm_subgraph}. We + graph constructions in Figure~\ref{fig:norm_subgraph}. We verify that for each of the two graph constructions, the cut capacities are equal to the energies, plus a constant. } @@ -1655,7 +1656,7 @@ term coming from our aim to minimize the total variation. The fidelity term of our energy function in \eqref{eq:total_energy} simplifies to \begin{align} - E_\lambda^x(0) &= + E_\lambda^x(0) &= N_x(\lambda + 1) - N_x(\lambda) \label{eq:fidelity_energy0} \\ @@ -1665,18 +1666,18 @@ simplifies to where $E_\lambda^x(0)$ might be positive or negative depending on $\lambda$ and the pixel value $v_x$. -Figure \ref{fig:norm_subgraph} shows how graphs can be constructed to +Figure~\ref{fig:norm_subgraph} shows how graphs can be constructed to represent this part of the total energy. The construction differs -depending on whether $E_\lambda^x(0)$ is positive or negative. Table -\ref{tab:fid_energy} shows how the cuts correspond to the values of -$u^\lambda_x$ and we can easily verify that the constructed graph +depending on whether $E_\lambda^x(0)$ is positive or negative. +Table~\ref{tab:fid_energy} shows how the cuts correspond to the values +of $u^\lambda_x$ and we can easily verify that the constructed graph actually represents the fidelity term in the energy function. \subsubsection{Regularization term} For our neighboring relation in \eqref{eq:total_energy} of the form \begin{equation} - E^{x,y}(u^\lambda_x, u^\lambda_y) = + E^{x,y}(u^\lambda_x, u^\lambda_y) = w_{xy} \abs{u^\lambda_x - u^\lambda_y} \end{equation} we have @@ -1689,18 +1690,18 @@ we have \end{aligned} \label{eq:neigh_energies} \end{equation} -and by Theorem \ref{thm:graph_rep_id} our energy function is graph +and by Theorem~\ref{thm:graph_rep_id} our energy function is graph representable. In \cite{kolmogorov2002energy}, Kolmogorov and Zabih presents a way to construct a graph for any graph representable function -on the form shown in Theorem \ref{thm:graph_rep_id}. Since the +on the form shown in Theorem~\ref{thm:graph_rep_id}. Since the energies in \eqref{eq:neigh_energies} are especially simple, the construction and presentation is simplified. -Figure \ref{fig:neigh_subgraph} shows two different ways of how a graph -can be constructed to represent the regularization term. See Table -\ref{tab:neigh_energy} for an overview of how the two values of -$u^\lambda_x$ corresponds to cuts in the graph. +Figure~\ref{fig:neigh_subgraph} shows two different ways of how a graph +can be constructed to represent the regularization term. +See Table~\ref{tab:neigh_energy} for an overview of how the two values +of $u^\lambda_x$ corresponds to cuts in the graph. \begin{table}[t] \centering @@ -1709,7 +1710,7 @@ $u^\lambda_x$ corresponds to cuts in the graph. configuration the corresponding energy and the cut yielding this configuration is shown. The last two columns show the capacities of the cut in the two alternative graph constructions shown in - Figure \ref{fig:neigh_subgraph}. We verify that for each of the + Figure~\ref{fig:neigh_subgraph}. We verify that for each of the two graph constructions, the cut capacities are equal to the energies, plus a constant. } @@ -1729,7 +1730,7 @@ $u^\lambda_x$ corresponds to cuts in the graph. \label{tab:neigh_energy} \end{table} -Figure \ref{fig:big_graph} shows a visualization of how the final +Figure~\ref{fig:big_graph} shows a visualization of how the final graph might look with all its edges. The source will have a lot of outgoing edges, one for each pixel, while the sink has one incoming edge for each pixel. The vertices corresponding to the pixels are only @@ -1792,7 +1793,7 @@ T)$ as \label{eq:net_flow_cut} \end{equation} Note how this definition differs from the capacity of a cut $c(S,T)$ in -Definition \ref{def:s_t_cut}. While the capacity of a cut represents how +Definition~\ref{def:s_t_cut}. While the capacity of a cut represents how much flow it is maximally possible to send from $S$ to $T$, the net flow across a cut represents the net amount of flow going across the cut, counting negatively the flow that goes back from $T$ to $S$. @@ -1887,13 +1888,13 @@ max-flow min-cut theorem. \end{enumerate} \end{theorem} See \cite{cormen2009introduction} for a proof, and remark that because -of the inequality in \eqref{eq:flow_leq_cut}, the cut in Statement -\ref{it:max_flow_min_cut} is a minimum cut. The theorem does not tell us -how to find such a cut, and there are multiple ways. One possibility is -to take $S$ to be all vertices reachable from the source in the residual -graph and $T = V - S$. +of the inequality in \eqref{eq:flow_leq_cut}, the cut in +Statement~\ref{it:max_flow_min_cut} is a minimum cut. The theorem does +not tell us how to find such a cut, and there are multiple ways. One +possibility is to take $S$ to be all vertices reachable from the source +in the residual graph and $T = V - S$. -Figure \ref{fig:aug_flow} shows a simple graph which already has five +Figure~\ref{fig:aug_flow} shows a simple graph which already has five units flowing from $s$ to $t$. The marked path is a possible augmenting path, and note that it follows an edge in $E$ in the reverse direction, made possible by the construction of the residual graph just @@ -1906,7 +1907,7 @@ described. \subsection{Ford-Fulkerson} The Ford-Fulkerson algorithm is the most basic augmenting path algorithm, which can be extended to more advanced algorithms. It is -stated in pseudocode in Algorithm \ref{alg:ford_fulkerson}, and the idea +stated in pseudocode in Algorithm~\ref{alg:ford_fulkerson}, and the idea is to augment the flow along paths from $s$ to $t$ until it is no longer possible. \begin{algorithm} @@ -2028,7 +2029,7 @@ and a valid distance labeling $d$. The push procedure moves excess flow from an active vertex along an edge $(u, v) \in E_f$ for which $d(u) = d(v) + 1$, i.e.\ to a vertex with a smaller distance label. We call such edges \emph{admissible}. See -Algorithm \ref{alg:push} for a pseudocode implementation of the push +Algorithm~\ref{alg:push} for a pseudocode implementation of the push operation. Assuming that $f$ is a valid preflow, it is easy to verify that the @@ -2065,7 +2066,7 @@ procedure on some admissible edge $(u, v)$. %\State $e(v) \mathrel{+}= f_\text{aug}$ \EndFunction \end{algorithmic} -\caption{ +\caption{% \sf The push procedure of the Push-Relabel algorithm. } \label{alg:push} @@ -2076,7 +2077,7 @@ procedure on some admissible edge $(u, v)$. The relabel procedure is our tool for changing the distance labels of the vertices. It changes the label of a vertex to the greatest possible value, which is one more than the lowest label among its neighbors in -the residual graph. See Algorithm \ref{alg:relabel} for a pseudocode +the residual graph. See Algorithm~\ref{alg:relabel} for a pseudocode implementation. \begin{algorithm} \begin{algorithmic} @@ -2118,10 +2119,10 @@ vertices and a maximum preflow is obtained. A vertex $u$ can only be successfully relabeled to obtain a new label if the outgoing edges of $u$ in the residual graph have changed since the previous relabeling. This is why the push and relabel procedures often -are combined into a \emph{discharge} procedure as shown in Algorithm -\ref{alg:discharge}. When it is run on an active vertex $u$, we push as -much as possible of the excess flow to other vertices before the vertex -is relabeled. +are combined into a \emph{discharge} procedure as shown in +Algorithm~\ref{alg:discharge}. When it is run on an active vertex $u$, +we push as much as possible of the excess flow to other vertices before +the vertex is relabeled. \begin{algorithm} \begin{algorithmic} \Function{Discharge}{$u$} @@ -2258,13 +2259,13 @@ might perform differently on different kinds of graphs. The most used heuristics are the gap and global relabeling heuristics, both aiming to reduce the total number of relabeling steps. -The gap relabeling heuristic aims to find a label $k$ as in Theorem -\ref{thm:cut_identification} such that no vertex has that label. +The gap relabeling heuristic aims to find a label $k$ as in +Theorem~\ref{thm:cut_identification} such that no vertex has that label. From vertices $u$ with $d(u) > k$, there are no unsaturated edges going to vertices with smaller distance labels, so no more flow can ever find its way from these vertices to the sink. These vertices are therefore given the label $N$ and never considered again as they will never become -active. Algorithm \ref{alg:gap} shows a pseudocode representation of +active. Algorithm~\ref{alg:gap} shows a pseudocode representation of what is done once a gap $k$ is found. \begin{algorithm} @@ -2333,22 +2334,22 @@ valid initial labeling. When running the push-relabel algorithm with the gap heuristic, we can be sure that there will never be a vertex $u$ with label $d(u) = N-1$ at the end of the algorithm, i.e.\ we know that there will always be a gap -at label $N-1$. This can be seen using the same reasoning as in Theorem -\ref{thm:cut_identification}, because if there was a vertex with label +at label $N-1$. This can be seen using the same reasoning as in +Theorem~\ref{thm:cut_identification}, because if there was a vertex with label $N-1$, there would only be $N-3$ vertices possibly having labels in $\{1, \ldots, N-2\}$, so a gap must exist somewhere in that interval. When using the gap relabeling heuristic, such a gap can not exist, so we can conclude that there is no vertex with label $N-1$. -Using Theorem \ref{thm:cut_identification} we can then conclude that the +Using Theorem~\ref{thm:cut_identification} we can then conclude that the sets $S = \{ u \in V : d(u) \geq N\}$ and $T = V - S$ form a minimum cut of the graph. \subsection{Parametric push-relabel algorithm} Now we have an algorithm for finding a minimum $s$-$t$-cut in a graph, -so let's return to the graph constructed in Section -\ref{sec:graph_construction}. For every level $\lambda \in \{0, +so let's return to the graph constructed in +Section~\ref{sec:graph_construction}. For every level $\lambda \in \{0, \ldots, L\}$ we want to find a minimum $s$-$t$-cut which gives us the thresholded image $u^\lambda$. These can then hopefully be stacked together to form the final image $u$. @@ -2364,26 +2365,27 @@ the next. \input{fig/norm_evolution} \end{figure} -Going back to the graph representations in Figure -\ref{fig:norm_subgraph} and Figure \ref{fig:neigh_subgraph} we know that -only the capacity of edges from sub-graphs representing the fidelity -term depend on our level parameter $\lambda$. From the expression in -\eqref{eq:fidelity_energy0}, visualized in Figure -\ref{fig:norm_evolution}, we see that the energy term $E_\lambda^x(0)$ -increases monotonically with increasing $\lambda$ parameter. Let $u, v -\in V - \{s, t\}$. Since the edges in Figure \ref{fig:norm_subgraph} now -are the only ones depending on $\lambda$, the following is true for -\emph{decreasing} values of $\lambda$ +Going back to the graph representations in +Figure~\ref{fig:norm_subgraph} and Figure~\ref{fig:neigh_subgraph} we +know that only the capacity of edges from sub-graphs representing the +fidelity term depend on our level parameter $\lambda$. From the +expression in \eqref{eq:fidelity_energy0}, visualized in +Figure~\ref{fig:norm_evolution}, we see that the energy term +$E_\lambda^x(0)$ increases monotonically with increasing $\lambda$ +parameter. Let $u, v \in V - \{s, t\}$. Since the edges in +Figure~\ref{fig:norm_subgraph} now are the only ones depending on +$\lambda$, the following is true for \emph{decreasing} values of +$\lambda$ \begin{description} \item[Edges from $s$ to $u$\textmd{:}] - As seen in Figure \ref{fig:norm_subgraph_neg} the capacity of + As seen in Figure~\ref{fig:norm_subgraph_neg} the capacity of these edges will increase monotonically with decreasing $\lambda$ parameter. \item[Edges from $u$ to $v$\textmd{:}] These edges have no $\lambda$-dependence and will remain unchanged. \item[Edges from $v$ to $t$\textmd{:}] - As seen in Figure \ref{fig:norm_subgraph_pos} the capacity of + As seen in Figure~\ref{fig:norm_subgraph_pos} the capacity of these edges will decrease monotonically with decreasing $\lambda$ parameter. \end{description} @@ -2417,10 +2419,10 @@ we reset the flow and labels. \subsubsection{Output image construction} -We mentioned already in Section \ref{sec:total_energy} that in order to +We mentioned already in Section~\ref{sec:total_energy} that in order to be able to construct our output image $u$, the thresholded images $u^\lambda$ would have to stack one on top of the other as shown in -Figure \ref{fig:img_decomp}. Because of the reuse of the distance labels +Figure~\ref{fig:img_decomp}. Because of the reuse of the distance labels between the iterations of the push-relabel algorithm, we can guarantee that this is possible. @@ -2492,7 +2494,7 @@ Goldfarb and Yin \cite{goldfarb2009parametric} have found that the divide and conquer approach only yields improved performance when using the $L^2$ norm in the fidelity term. This has to do with the fact that the fidelity term for the $L^1$ norm only changes once for each pixel, -as we can see in Figure \ref{fig:norm_evolution}, reducing the amount of +as we can see in Figure~\ref{fig:norm_evolution}, reducing the amount of work that has to be done in each iteration of the regular parametric push-relabel algorithm. @@ -2501,9 +2503,9 @@ See \cite{gallo1989fast}, \cite{hochbaum2001efficient} and \subsection{Implementation} -A \cpp\ implementation can be found in Appendix -\ref{app:c++implementation}. It uses the open computer vision library -OpenCV \cite{opencv_library} to load and save image files. +A \cpp\ implementation can be found in +Appendix~\ref{app:c++implementation}. It uses the open computer vision +library OpenCV \cite{opencv_library} to load and save image files. Note that when implementing maximum flow algorithms it is not a good idea, memory- and performance-wise, to actually construct the residual