From: Bjørn Rustad Date: Tue, 3 Feb 2015 18:59:48 +0000 (+0100) Subject: FIX EVERYTHING X-Git-Url: http://git.rustad.me/?a=commitdiff_plain;h=e89f8f503005ed63bcbbb3f77791ee6b6fbdd344;p=master FIX EVERYTHING --- diff --git a/commands.tex b/commands.tex index 2a7a689..cb1018b 100644 --- a/commands.tex +++ b/commands.tex @@ -4,7 +4,7 @@ \newcommand{\lnorm}[1]{\lVert #1 \rVert} \newcommand{\diver}{\mathop{\mathrm{div}}\nolimits} \newcommand{\TV}{\mathrm{TV}} -\newcommand{\TVA}{\mathrm{TV}_A} +\newcommand{\TVA}{\mathrm{TV}_{\!A}} \newcommand{\BV}{\mathrm{BV}} \newcommand{\Per}{\mathrm{Per}} \newcommand{\PerA}{\mathrm{Per}_A} diff --git a/continuous.tex b/continuous.tex index ff44c6e..8c3e80e 100644 --- a/continuous.tex +++ b/continuous.tex @@ -1,11 +1,11 @@ \chapter{Continuous formulation} -From the introduction we see that there are many different approaches to +In the previous chapter we saw that there are many different approaches to the image restoration problem, all with their own strengths and weaknesses. The method considered in this thesis is an anisotropic total -variation formulation, and the hope is that we keep the strengths of the +variation formulation, and the aim is to keep the strengths of the anisotropic diffusion and total variation methods, while eliminating -some of their respective weaknesses. +some of their weaknesses. This chapter will be devoted to the continuous formulation of the method. We will look at the functional we want to minimize and its @@ -13,7 +13,7 @@ different forms, and briefly discuss its well-posedness. Through the anisotropic coarea formula, the anisotropic total variation is rewritten as an integral of the perimeter of all the level sets of the image. -Later, the anisotropic Cauchy--Crofton formula is introduced to make it +Following that, the anisotropic Cauchy--Crofton formula is introduced to make it feasible to calculate the perimeter of these level sets. All of this leads up to the discretization of our functional in the next chapter. @@ -46,8 +46,8 @@ functional, there is no gain in reducing the height of the edge as before. Note that $u(x)$ and $A(x)$ are always dependent on the position in the -image $x$, but we will sometimes drop the $x$, when no misunderstandings -are possible. +image $x$, but we will sometimes drop the $x$, when the intended meaning +is clear. As we will not always be working with differentiable images, we extend the definition of the total variation functional. Being symmetric @@ -67,9 +67,10 @@ matrices as $A = \Ahalf \Ahalf$. We can then write \end{aligned} \end{equation} where $\xi$ and $\eta = \Ahalf \xi$ are in $C_c^\infty(\Omega, -\mathbb{R}^2)$, the space of $C^\infty$ vector fields with compact support. -In the following we let $\norm{\xi}_A = \sup_x \sqrt{\xi^T A \xi}$ and -$\norm{\eta}_A^* = \sup_x \sqrt{\eta^T A^{-1} \eta}$, and with that we +\mathbb{R}^2)$, the space of smooth vector fields with compact support. +In the following we define the norms $\lnorm{\xi}_A = \sup_x (\xi^T A +\xi)^{\sfrac{1}{2}}$ and +$\lnorm{\eta}_A^* = \sup_x (\eta^T A^{-1} \eta)^{\sfrac{1}{2}}$, and with that we present the formal definition of the anisotropic total variation. \begin{definition}[Anisotropic total variation] For a function $u \in BV(\Omega)$ and a continuous symmetric @@ -122,8 +123,8 @@ The image is smoothed by a Gaussian filter as described in Section~\ref{sec:diffusion_filtering}: $f_\sigma = K_\sigma * \tilde{f}$, where $\tilde{f}$ is the symmetric extension of the initial image $f$ in $\mathbb{R}^2$. The smoothing parameter $\sigma$ is called the -\emph{noise scale}, and with this we aim to avoid being too sensitive to -noise in $f$. +\emph{noise scale}, and it controls the scale at which details are +considered to be noise. As seen in Figure~\ref{fig:edges}, the edge detector is fine for detecting edges, but it can not give us information about larger @@ -136,10 +137,9 @@ Its eigenvalues are $\lambda_1 = \labs{\nabla f_\sigma(x)}^2$ and $\lambda_2 = 0$ with corresponding eigenvectors $v_1$ and $v_2$ parallel and perpendicular to $\nabla f_\sigma(x)$ respectively. -We are also interested in identifying and being sensitive to features in -a neighborhood around the point $x$, such as corners or curved edges and -\fixme{coherent structures}. This is why we introduce the component-wise -convolution with $K_\rho$ such that +To detect features in a neighborhood around the point $x$, such as +corners, curved edges and coherent structures we introduce the +component-wise convolution with $K_\rho$ such that \begin{equation} S_\rho(x) := K_\rho * \big( \nabla f_\sigma(x) \otimes \nabla f_\sigma(x) \big)(x). @@ -147,21 +147,33 @@ convolution with $K_\rho$ such that \end{equation} The parameter $\rho$, called the \emph{integration scale}, controls the size of the neighborhood -which affects the structure tensor. It controls the size of the +which affects the structure tensor. Thus it defines the size of the structures we want our anisotropy tensor to be sensitive to. \begin{figure} \centering{} - \includegraphics[width=0.4\textwidth]{fig/edges.png} + \begin{subfigure}[b]{0.45\textwidth} + \includegraphics[width=\textwidth]{fig/factory/finger/n_q200.png} + \end{subfigure} + ~ + \begin{subfigure}[b]{0.45\textwidth} + \includegraphics[width=\textwidth]{fig/factory/finger/r_p2_n32_b15000_g150_r20_s3_edge.png} + \end{subfigure} \caption[The edge detector $\abs{\nabla f_\sigma(x)}$]{% - The largest eigenvalue of the structure tensor is $\abs{\nabla - f_\sigma(x)}^2$ and can as we see function as an edge detector. + A noisy fingerprint on the left, and the largest eigenvalue of + the structure tensor is $\abs{\nabla f_\sigma(x)}^2$ on the + left, which---as we can see---functions as an edge + detector. } \label{fig:edges} \end{figure} The smoothed tensor $S_\rho(x)$ can easily be verified to be symmetric -positive semi-definite, just like $S_0(x)$. We order the two real +positive semi-definite, just like $S_0(x)$. In addition, when $\rho > +0$, the elements of $S_\rho$ are smooth maps from $\Omega$ to +$\mathbb{R}$. + +We order the two real eigenvalues $\lambda_1 \geq \lambda_2$ with corresponding eigenvectors $v_1$ and $v_2$. From the characteristic polynomial of $S_\rho(x) = (\begin{smallmatrix} s_{11} & s_{12} \\ s_{12} & s_{22} @@ -175,8 +187,8 @@ expression for the eigenvalues The vector $v_1$ will then indicate the direction of most variation in the -neighborhood. From this we see that an edge would give $\lambda_1 \gg -\lambda_2 \approx 0$, while smooth areas would give $\lambda_1 \approx +neighborhood. An edge will give $\lambda_1 \gg +\lambda_2 \approx 0$, while smooth areas will give $\lambda_1 \approx \lambda_2 \approx 0$. In corners we have variation in the direction of $v_1$ but also perpendicular to $v_1$, so we will have $\lambda_1 \approx \lambda_2 \gg 0$. Thus the quantity $(\lambda_1 - \lambda_2)^2$ @@ -219,9 +231,14 @@ eigenvalues are different. A visualization of the two tensors can be seen in Figure~\ref{fig:tensor_viz} where the two tensors are shown at an edge in the image. +\begin{figure} + \input{fig/tensor_viz} +\end{figure} + In smooth areas, $\sigma_1 \approx 1$ and -$A(x)$ will be close to the identity matrix, while $A(x)$ will reduce -the effect of $\nabla u$ across edges. +$A(x)$ will be close to the identity matrix. At or around edges, +$\sigma_1$, which corresponds to the eigenvector perpendicular to the +edge, will be small. Around corners $A(x)$ will be close to the identity matrix, which gives regularization similar to smooth areas. This is one possible down-side @@ -229,8 +246,8 @@ of this tensor choice. The parameter $\omega$ controls the amount of anisotropy in the method, such that if it is very large we are left with the identity matrix. Note -that changing the parameter $\omega$ implicitly affects how much -restoration is done. For an image $u$, decreasing $\omega$ will, all +that changing the parameter $\omega$ implicitly affects the amount of +regularization applied. For an image $u$, decreasing $\omega$ will, all else being equal, decrease the lowest eigenvalue of $A(x)$ and in turn decrease the anisotropic total variation $\TVA(u)$. @@ -239,16 +256,18 @@ decomposition is not well-defined. This is not a problem though, since $\Sigma(x)$ will be the identity matrix, so any orthogonal matrix will suffice for $U(x)$. +Note that the eigenvalues of $S_\rho$ are continuous, and so are the +eigenvectors (ignoring their sign) except possibly when $\lambda_1 = +\lambda_2$. Thus $A$ is also continuous except possibly in these points. +When $\lambda_1 = \lambda_2$ however, the eigenvalues of $A$ will both +be 1, and $A$ is the identity matrix. Thus we can argue that if +$S_\rho(x) \to \lambda I$ then $A(x) \to I$ and $A$ is continuous in all +of $\Omega$. + See \cite{weickert1999coherence} for a different tensor construction, made to enhance flow structures in the image, relevant in for example fingerprint analysis. -\fixme{argue that $A$ is continuous?} - -\begin{figure} - \input{fig/tensor_viz} -\end{figure} - \section{Well-posedness} \fixme{Rating: 7/10} @@ -306,25 +325,25 @@ The anisotropic total variation dx \end{equation} can be thought of as---and has the properties of---a semi-norm, and is -therefore convex. The sum of the two is thus strictly convex, which, -given the existence of a minimizer, implies uniqueness. +therefore convex. The sum of the fidelity and regularization terms is +thus strictly convex, which, given the existence of a minimizer, implies +uniqueness. \subsection{Coercivity} Coercivity relates to how the functional behaves when the norm of the image $u$ tends to infinity. What we need in order to conclude with existence -is sequential coercivity. Thus we need all level sets $F^\alpha = \{ u \in -L^2(\Omega) : F(u) \leq \alpha \}$ to be \emph{sequentially +is weak sequential coercivity. Thus we need all level sets $F^\alpha = \{ u \in +L^2(\Omega) : F(u) \leq \alpha \}$ to be \emph{weakly sequentially pre-compact}, meaning that all sequences in the set contain a -subsequence converging to an element of the closure of the set. +subsequence weakly converging to an element of the closure of the set. It is obvious from the fidelity term that for some fixed $f \in L^2(\Omega)$, if $\lnorm{u}_{L^2} \to \infty$ then $F(u) \to \infty$. This implies that all the level sets $F^\alpha$ are bounded. -Since $L^2(\Omega)$ is a Hilbert space, all bounded sequences contain a +Since $L^2(\Omega)$ is a Hilbert space all bounded sequences contain a weakly convergent subsequence. Thus all the level sets $F^\alpha$ are -weakly sequentially pre-compact, since all sequences in these sets have -a subsequence weakly converging to a point in the closure of the set. +weakly sequentially pre-compact. \begin{figure} \input{fig/lower_semicont} @@ -550,7 +569,7 @@ image $u \geq 0$ as an integral over all the layers \label{eq:positive_int} \end{equation} Note that \eqref{eq:positive_int} only holds for non-negative images, -which complicates the following proof a little. +which complicates the proof of the anisotropic coarea formula a little. \begin{figure} \input{fig/eta_r} @@ -685,7 +704,7 @@ integration theory see for example \cite{bartle1995elements}. \end{cases} \end{aligned} \end{equation} - visualized in Figure~\ref{fig:eta_r} and \ref{fig:eta_r_diff}. By + visualized in Figure~\ref{fig:eta_r_both}. By composing the function $\eta_r$ with our image $u$ and using Green's formula, for example from \cite[Corollary 9.32]{grasmair2010anisotropic} we obtain @@ -699,8 +718,8 @@ integration theory see for example \cite{bartle1995elements}. for all vector fields $\xi \in C_c^\infty(\Omega, \mathbb{R}^2)$. The measure of $\{ x : u(x) = \lambda \text{ and } \nabla u(x) \neq 0\}$ is zero for all $\lambda$ following from \cite[Corollary \rom{1}, - Section 3.1.2]{evans1991measure}. Thus we can ignore the sets $\{ u - = s\}$ and $\{ u = s + r \}$. + Section 3.1.2]{evans1991measure}, thus we can ignore the sets $\{ u + = s\}$ and $\{ u = s + r \}$ in the integrals. Assuming that $\norm{\xi}_A^* \leq 1$ we obtain from \eqref{eq:mdef} and \eqref{eq:eta_greens} that \begin{equation} @@ -717,11 +736,14 @@ integration theory see for example \cite{bartle1995elements}. As the limit when $r \to 0$ of the left-hand side exists almost everywhere, suppose it exists at $s \in \mathbb{R}$. The integrand on the right-hand side $-\eta_r(u) \diver \xi$ approaches $-u^s - \diver \xi$ pointwise almost everywhere. We apply - Lebesgue's dominated convergence theorem - using that $\abs{\eta_r(u) \diver \xi} \leq - \abs{u^s \diver \xi}$ and that $\xi \in C^\infty_c(\Omega, - \mathbb{R}^2)$ is bounded by the extreme value theorem. Thus + \diver \xi$ pointwise almost everywhere. + We also have the bound + $\abs{\eta_r(u) \diver \xi} \leq + \abs{u^s \diver \xi}$ and know that + $\xi \in + C^\infty_c(\Omega, \mathbb{R}^2)$ is bounded following from the extreme value + theorem. Thus we can apply + Lebesgue's dominated convergence theorem, giving that $\abs{u^s \diver \xi}$ is integrable and \begin{equation} \int_\Omega -\eta_r(u) @@ -733,7 +755,7 @@ integration theory see for example \cite{bartle1995elements}. \end{equation} As this holds for any $\norm{\xi}_A^* \leq 1$, we get from the extended total variation - definition in \eqref{eq:extended_tv} that $m'(s) \geq \TVA(u_s)$ + definition in \eqref{eq:extended_tv} that $m'(s) \geq \TVA(u^s)$ almost everywhere and conclude using \eqref{eq:tva_geq_mder} that \begin{equation} \TVA(u) \geq \int_{-\infty}^\infty m'(t) \, dt \geq @@ -742,7 +764,9 @@ integration theory see for example \cite{bartle1995elements}. Combining the upper and lower bounds just proved, we have equality. \end{proof} This coarea formula is our first step in transforming the anisotropic -total variation into an easily discretizable expression. +total variation into an easily discretizable expression. It allows us to +consider each level $\lambda$ separately when calculating the total +variation. The anisotropic total variation of the thresholded images occurring in the anisotropic coarea formula is very much related to the size of the @@ -794,8 +818,8 @@ functional \PerA(\{u > \lambda\}; \Omega) \, d\lambda. \label{eq:functional_per} \end{equation} -The transformation is motivated by our upcoming Cauchy--Crofton -integration formula, and discretization, where an approximation of the +The transformation is motivated by our upcoming anisotropic Cauchy--Crofton +integration formula, and the discretization, where an approximation of the perimeter will be computed using a graph cut machinery. \section{Cauchy--Crofton formulas} @@ -812,9 +836,9 @@ a number of interesting integral formulas. Several of them fall in a category often referred to as \emph{Cauchy--Crofton style formulas}, and give ways to measure geometric objects using the set of all lines in the plane. The formulas presented here will give a way to measure the length -of a curve by counting the times it intersects line in the set of all +of a curve by counting the times it intersects lines in the set of all lines. The first formula will be for the isotropic case, and we will -use this to prove the anisotropic formula following it. +use it to prove the anisotropic formula following it. We write $\mathcal{L}$ for the set of all lines in the plane, and parametrize them as shown in Figure~\ref{fig:line_param}. A @@ -846,7 +870,7 @@ and the distance to origin $\rho$.}% \label{thm:euclidean_cauchy_crofton} \end{theorem} \begin{proof} - See \cite[Theorem 3, Section 1-7.]{do1976differential}. + See \cite[Theorem 3, Section 1-7]{do1976differential}. \end{proof} 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 @@ -867,7 +891,7 @@ 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] +\begin{theorem}[The anisotropic Cauchy--Crofton formula] Assume that our space $\Omega$ is equipped with a continuous positive definite metric tensor $M(x)$, whose eigenvalues are bounded by $0 < k \leq @@ -883,7 +907,7 @@ by a sum over all these lines. \end{equation} \label{thm:riemannian_cauchy_crofton} \end{theorem} -\begin{proof}[Proof of the Riemannian Cauchy--Crofton formula] +\begin{proof}[Proof of the anisotropic Cauchy--Crofton formula] 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 @@ -1002,7 +1026,7 @@ by a sum over all these lines. \nu)^{\sfrac{3}{2}}}. \end{equation} We are interested in the Jacobian of the whole transformation - $J_{\Sigma^2}(U^T \ell_{\theta, \rho})$, so + $J_{\Sigma^2}(U^T \ell_{\nu, \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})} @@ -1021,7 +1045,7 @@ by a sum over all these lines. \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 + Further we argue that the similar formula in \eqref{eq:riemannian_cauchy_crofton} holds for a non-constant but continuous metric tensor $M(x)$. By partitioning the domain into disjoint sets $U_i$ such that $\Omega = \cup_i U_i$, we make a @@ -1095,7 +1119,7 @@ by a sum over all these lines. \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 + by $g(t) = (K \cdot \dot{C}(t)^T \dot{C}(t))^{\sfrac{1}{2}}$. We know that $g(t)$ is integrable as its integral is exactly $\sqrt{K} \abs{C}$ and we have assumed that the curve is of finite length. This means we can apply Lebesgue's dominated convergence theorem to see that @@ -1154,7 +1178,9 @@ by a sum over all these lines. 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}. +each level set of the image. We will now see how the anisotropic +Cauchy--Crofton formula can help us calculate the perimeters of the +level sets. 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 @@ -1166,7 +1192,7 @@ 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. Thus a 90\textdegree{} rotation -is needed. +separates the two. If $P$ is a 90\textdegree{} rotation matrix we have \begin{equation} @@ -1191,7 +1217,7 @@ the boundary of $U$ that does not overlap with the boundary of $\Omega$ \label{eq:per_to_length2} \end{equation} -Now we make sure that all the assumptions of the Riemannian +Now we make sure that all the assumptions of the anisotropic 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}. @@ -1200,7 +1226,7 @@ The structure tensor is constructed as described in Section~\ref{sec:anisotropy_tensor} \begin{equation} S_\rho(x) = \left(K_{\rho} * \left( \nabla f_{\sigma} - \otimes \nabla f_{\sigma} \right)\right)(x), + \otimes \nabla f_{\sigma} \right)\right)(x). \end{equation} Because of the convolutions with the Gaussian function, this is a smooth map from $\bar{\Omega}$ to $\mathbb{R}^{2 \times 2}$. As we @@ -1210,17 +1236,17 @@ elements in 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 +are bounded above. Moreover, by the construction in \eqref{eq:sigma_construction}, there exists uniform bound $k$ such that -the smallest eigenvalue $\sigma_2$ of the anisotropy tensor $A(x)$ is bounded away +the smallest eigenvalue $\sigma_1$ of the anisotropy tensor $A(x)$ is bounded away from zero, as \begin{equation} - \sigma_2 = \left(1 + \sigma_1 = \left(1 + \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$ is continuous and positive -definite and thus the curve length calculation in +definite with bounded eigenvalues and thus the curve length calculation in \eqref{eq:per_to_length2} fulfills all the assumptions of the anisotropic Cauchy--Crofton formula in Theorem~\ref{thm:riemannian_cauchy_crofton}. @@ -1261,7 +1287,7 @@ and $\gamma_s = \partial \{u > s\} \cap \Omega$. 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} +corresponding eigenvector perpendicular to the edge. As mentioned previously, these transformations are motivated by the discretization that follows in the next chapter. diff --git a/discrete.tex b/discrete.tex index fd98fbf..257fb6d 100644 --- a/discrete.tex +++ b/discrete.tex @@ -15,9 +15,8 @@ a global minimizer in polynomial time. 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. - +L-1\}$, such that $u : \mathcal{G} \to \mathcal{P}$. 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. @@ -123,8 +122,9 @@ the approximation \end{aligned} \end{equation} 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 +$\mathcal{L}_D$. Note that $C$ is still a differentiable curve, not yet +discretized. +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 @@ -133,25 +133,21 @@ 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}, and thus also depends on $\phi$. -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 -value of 0 or 1 on this grid, as shown in Figure \fixme{ref}. - The choice of our discrete set of lines $\mathcal{L}_D$ is important, as it will decide the accuracy of our approximation. We need some sensible restrictions on the set $\mathcal{L}_D$ to simplify the further discussion. All lines intersect at least two grid points, and from the periodicity of our grid they thus intersect an infinite number of grid points. This puts some -restrictions on the angles we can choose. For each angle, we include all +restrictions on the angles we can choose. For each angle included, we include all 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. +shown in Figure~\ref{fig:line_neigh}. We write $\mathcal{N}(x)$ for the +neighborhood of grid point $x$. Thus if there is an edge $e_{ab}$ in our +discretization $\mathcal{L}_D$, then $a \in \mathcal{N}(b)$. Extending the edges shown in +Figure~\ref{fig:line_neigh} 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$. @@ -172,7 +168,7 @@ 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 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 +terminals of $e$ lie on different sides of the curve $C$. This approximation is exact for zero or one intersection points, but will, as we see in Figure~\ref{fig:curve_edge}, be wrong when we have more. @@ -217,8 +213,9 @@ will be provided by the following lemma. \label{lem:delta_rho} \end{lemma} \begin{proof} - Consider a line $\ell$ intersecting the point $(p, q) \in - \mathbb{Z}^2$ in the grid. + Consider a line $\ell$ intersecting the point in the grid given by + the indices $(p, q) \in + \mathbb{Z}^2$. 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 grid points. @@ -258,9 +255,9 @@ approximation of \eqref{eq:cauchy_crofton_approx2} we obtain \end{equation} where the sum is over all edges crossing the curve. -The curve length we initially wanted to calculate was the perimeter +The curve length we initially wanted to calculate was the level set perimeter $\PerA(\{u > \lambda\}; \Omega)$ in \eqref{eq:per_approx1}. To find -edges that crosses the curve, we identify the edges that have one +edges that crosses this boundary 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} @@ -273,7 +270,6 @@ Thus we rewrite the sum over $e \cap C$ such that 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 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. @@ -300,7 +296,7 @@ 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 \fixme{(?)} in increasing level values, i.e.\ +be monotonically decreasing 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}. @@ -310,8 +306,6 @@ Later we will present a graph cut algorithm that find the thresholded images minimizing each level, \emph{while guaranteeing that they meet this requirement.} -\fixme{maybe with a $w_{xy}$ definition} - \subsubsection{Consistency} \fixme{Rating: 6/10} @@ -365,10 +359,11 @@ Further, if $\sup \norm{e} \to 0$ it is obvious that the tensor average 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 +Consider now the discretization in $\rho$. 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. +part of the discretization is consistent. \begin{figure} \input{fig/line_midpoint} @@ -406,7 +401,7 @@ of this square $n$ goes to infinity. \end{figure} 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 +for multiple points that lie on the same line from the origin. 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 @@ -418,11 +413,12 @@ 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 +largest $\Delta \phi_k$ will be when $\phi_k = \sfrac{m \pi}{2}$ for $m +\in \mathbb{Z}$, so +around vertical and horizontal edges. +The supremum can then be calculated to be \begin{equation} - \sup \Delta \phi_k = \frac{2}{2} \arctan \frac{1/n}{n/2} = \arctan + \sup \Delta \phi_k = 2 \cdot \sup \frac{\phi_{k+1} - \phi_k}{2} = \arctan \frac{1/n}{n/2} = \arctan \frac{2}{n^2} \to 0. \end{equation} @@ -432,14 +428,14 @@ diagonal of the square such that \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 -$\delta^2 = \Delta \rho \norm{e}$ and the fact that $\norm{e} \geq -\delta$. Thus for the inter-line distance $\Delta \rho$ we have +$\delta^2 = \Delta \rho \norm{e}$ and $\norm{e} \geq +\delta$. Thus for the inter-line distance $\Delta \rho_k$ we have \begin{equation} - \sup \Delta \rho = \sup \frac{\delta^2}{\norm{e}} \leq + \sup \Delta \rho_k = \sup \frac{\delta^2}{\norm{e}} \leq \frac{\delta^2}{\delta} = \delta \to 0. \end{equation} -And thus the approximation has been shown to be equivalent to +Hence the approximation has been shown to be equivalent to 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 @@ -461,7 +457,7 @@ are then combined to form the final restored image $u$. In this section we will look at how these graphs are constructed such that their minimum cuts correspond to the minimizers of the functional -$F^\lambda$. The description is taken from my project work +$F^\lambda$. The description is taken with some small adjustments from my project work \cite{project}, and is included here for completeness. \subsection{Graphs} @@ -480,8 +476,8 @@ capacity $c(u, v) = q > 0$ and it means that it is possible to send a will let $c(u, v) = 0$ for any pair $(u, v) \not\in E$, and we do not allow self-loops in our graph. When a directed graph $G$ is equipped with capacity function $c$, one might call it a capacitated graph or a -graph, but as all our graphs will be capacitated from this point, we -will just call them graphs and we write $G = (V, E, c)$. +flow network, but as all our graphs will be capacitated from this point, +we will just call them graphs and we write $G = (V, E, c)$. There are two special vertices in the graph, the source $s$ and the sink $t$. Contrary to other vertices, which can neither produce nor @@ -594,7 +590,7 @@ depend on. We will construct a graph in such a way that if a vertex $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 += 1$, and if it ends up on the sink side we set $u^\lambda_x = 0$, as in 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. @@ -607,8 +603,8 @@ term coming from our aim to minimize the total variation. \begin{table}[b] \centering - \caption[Different combinations for the term - $F_\lambda^x(u^\lambda_x)$]{Each row represents one of the two possible values of + \caption[Graph construction for $F_\lambda^x(u^\lambda_x)$]{Each row + represents one of the two possible values of $u^\lambda_x \in \{0,1\}$. The functional $F_\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 @@ -622,8 +618,10 @@ term coming from our aim to minimize the total variation. & Graph \subref{fig:norm_subgraph_pos} cut cap. & Graph \subref{fig:norm_subgraph_neg} cut cap. \\ \hline - $0$ & $0$ & $(\{s, u^\lambda_x\}, \{t\})$ & $0$ & $-F_\lambda^x(1)$ \\ - $1$ & $F_\lambda^x(1)$ & $(\{s\}, \{u^\lambda_x, t\})$ & $F_\lambda^x(1)$ & $0$ + $0$ & $0$ & $(\{s\}, \{u^\lambda_x, t\})$ & $0$ & + $-F_\lambda^x(1)$ \\ + $1$ & $F_\lambda^x(1)$ & $(\{s, u^\lambda_x\}, \{t\})$ & + $F_\lambda^x(1)$ & $0$ \end{tabular} \label{tab:fid_energy} \end{table} @@ -646,7 +644,7 @@ simplifies to \label{eq:fidelity_energy1} \end{align} where $F_\lambda^x(1)$ might be positive or negative depending on -$\lambda$ and the pixel value $v_x$. +$\lambda$ and the pixel value $f_x$. Figure~\ref{fig:norm_subgraph} shows how graphs can be constructed to represent this part of functional. The construction differs @@ -681,12 +679,12 @@ 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. +See Table~\ref{tab:neigh_energy} for an overview of how the values +of $(u^\lambda_x, u^\lambda_y)$ corresponds to cuts in the graph. \begin{table}[t] \centering - \caption[Different combinations for the term $F^{x,y}(u^\lambda_x, + \caption[Graph construction for $F^{x,y}(u^\lambda_x, u^\lambda_y)$]{An overview of the possible configurations of the variables in the term $F^{x,y}(u^\lambda_x, u^\lambda_y)$. For each configuration the corresponding functional value and the cut yielding this @@ -704,10 +702,11 @@ of $u^\lambda_x$ corresponds to cuts in the graph. & Alt.\ \subref{fig:neigh_subgraph_alt1} cut cap. & Alt.\ \subref{fig:neigh_subgraph_alt2} cut cap. \\ \hline - $(0, 0)$ & $0$ & $(\{s, u^\lambda_x, u^\lambda_y\}, \{t\})$ & $w_{xy}$ & $0$ \\ - $(0, 1)$ & $w_{xy}$ & $(\{s, u^\lambda_x\}, \{u^\lambda_y, t\})$ & $2w_{xy}$ & $w_{xy}$ \\ - $(1, 0)$ & $w_{xy}$ & $(\{s, u^\lambda_y\}, \{u^\lambda_x, t\})$ & $2w_{xy}$ & $w_{xy}$ \\ - $(1, 1)$ & $0$ & $(\{s\}, \{u^\lambda_x, u^\lambda_y, t\})$ & $w_{xy}$ & $0$ + $(0, 0)$ & $0$ & $(\{s\}, \{u^\lambda_x, u^\lambda_y, t\})$ + & $w_{xy}$ & $0$ \\ + $(0, 1)$ & $w_{xy}$ & $(\{s, u^\lambda_y\}, \{u^\lambda_x, t\})$ & $2w_{xy}$ & $w_{xy}$ \\ + $(1, 0)$ & $w_{xy}$ & $(\{s, u^\lambda_x\}, \{u^\lambda_y, t\})$ & $2w_{xy}$ & $w_{xy}$ \\ + $(1, 1)$ & $0$ & $(\{s, u^\lambda_x, u^\lambda_y\}, \{t\})$ & $w_{xy}$ & $0$ \end{tabular} \label{tab:neigh_energy} \end{table} diff --git a/fig/area_proof.tex b/fig/area_proof.tex index d19c6e1..096c378 100644 --- a/fig/area_proof.tex +++ b/fig/area_proof.tex @@ -6,62 +6,63 @@ } } - \draw[blue] (0,0) rectangle (1,1); - \draw[blue] (1,0) rectangle (2,1); - \draw[blue] (2,0) rectangle (3,1); - \draw[blue] (3,0) rectangle (4,1); - \draw[blue] (4,0) rectangle (5,1); + \draw[blue, thick] (0,0) rectangle (1,1); + \draw[blue, thick] (1,0) rectangle (2,1); + \draw[blue, thick] (2,0) rectangle (3,1); + \draw[blue, thick] (3,0) rectangle (4,1); + \draw[blue, thick] (4,0) rectangle (5,1); - \draw[blue] (0,1) rectangle (1,2); - \draw[blue] (1,1) rectangle (2,2); - \draw[blue] (2,1) rectangle (3,2); - \draw[blue] (3,1) rectangle (4,2); - \draw[blue] (4,1) rectangle (5,2); + \draw[blue, thick] (0,1) rectangle (1,2); + \draw[blue, thick] (1,1) rectangle (2,2); + \draw[blue, thick] (2,1) rectangle (3,2); + \draw[blue, thick] (3,1) rectangle (4,2); + \draw[blue, thick] (4,1) rectangle (5,2); - \draw[blue] (0,2) rectangle (1,3); - \draw[blue] (1,2) rectangle (2,3); - \draw[blue] (2,2) rectangle (3,3); - \draw[blue] (3,2) rectangle (4,3); - \draw[blue] (4,2) rectangle (5,3); + \draw[blue, thick] (0,2) rectangle (1,3); + \draw[blue, thick] (1,2) rectangle (2,3); + \draw[blue, thick] (2,2) rectangle (3,3); + \draw[blue, thick] (3,2) rectangle (4,3); + \draw[blue, thick] (4,2) rectangle (5,3); \clip (0,0) rectangle (5,3); - \path[nedge=4cm,red] (0,3) -- (2,2); - \path[nedge=4cm,red] (0,2) -- (2,1); - \path[nedge=4cm,red] (0,1) -- (2,0); + \path[nedge=4cm,red, thick] (0,3) -- (2,2); + \path[nedge=4cm,red, thick] (0,2) -- (2,1); + \path[nedge=4cm,red, thick] (0,1) -- (2,0); - \path[nedge=4cm,red] (-1,1) -- (1,0); + \path[nedge=4cm,red, thick] (-1,1) -- (1,0); - \path[nedge=4cm,red] (1,3) -- (3,2); - \path[nedge=4cm,red] (1,2) -- (3,1); - \path[nedge=4cm,red] (1,1) -- (3,0); + \path[nedge=4cm,red, thick] (1,3) -- (3,2); + \path[nedge=4cm,red, thick] (1,2) -- (3,1); + \path[nedge=4cm,red, thick] (1,1) -- (3,0); - \path[nedge=4cm,red] (2,3) -- (4,2); - \path[nedge=4cm,red] (3,3) -- (5,2); - \path[nedge=4cm,red] (4,3) -- (6,2); + \path[nedge=4cm,red, thick] (2,3) -- (4,2); + \path[nedge=4cm,red, thick] (3,3) -- (5,2); + \path[nedge=4cm,red, thick] (4,3) -- (6,2); - \path (5,1) [edge,-,red] -- ($(6,0)!(5,1)!(4,1)$); - \path (4,1) [edge,-,red] -- ($(5,0)!(4,1)!(3,1)$); - \path (3,1) [edge,-,red] -- ($(4,0)!(3,1)!(2,1)$); - \path (2,1) [edge,-,red] -- ($(3,0)!(2,1)!(1,1)$); - \path (1,1) [edge,-,red] -- ($(2,0)!(1,1)!(0,1)$); + \path (5,1) [edge,-,red, thick] -- ($(6,0)!(5,1)!(4,1)$); + \path (4,1) [edge,-,red, thick] -- ($(5,0)!(4,1)!(3,1)$); + \path (3,1) [edge,-,red, thick] -- ($(4,0)!(3,1)!(2,1)$); + \path (2,1) [edge,-,red, thick] -- ($(3,0)!(2,1)!(1,1)$); + \path (1,1) [edge,-,red, thick] -- ($(2,0)!(1,1)!(0,1)$); - \path (5,2) [edge,-,red] -- ($(6,1)!(5,2)!(4,2)$); - \path (4,2) [edge,-,red] -- ($(5,1)!(4,2)!(3,2)$); - \path (3,2) [edge,-,red] -- ($(4,1)!(3,2)!(2,2)$); - \path (2,2) [edge,-,red] -- ($(3,1)!(2,2)!(1,2)$); - \path (1,2) [edge,-,red] -- ($(2,1)!(1,2)!(0,2)$); + \path (5,2) [edge,-,red, thick] -- ($(6,1)!(5,2)!(4,2)$); + \path (4,2) [edge,-,red, thick] -- ($(5,1)!(4,2)!(3,2)$); + \path (3,2) [edge,-,red, thick] -- ($(4,1)!(3,2)!(2,2)$); + \path (2,2) [edge,-,red, thick] -- ($(3,1)!(2,2)!(1,2)$); + \path (1,2) [edge,-,red, thick] -- ($(2,1)!(1,2)!(0,2)$); - \path (5,3) [edge,-,red] -- ($(6,2)!(5,3)!(4,3)$); - \path (4,3) [edge,-,red] -- ($(5,2)!(4,3)!(3,3)$); - \path (3,3) [edge,-,red] -- ($(4,2)!(3,3)!(2,3)$); - \path (2,3) [edge,-,red] -- ($(3,2)!(2,3)!(1,3)$); - \path (1,3) [edge,-,red] -- ($(2,2)!(1,3)!(0,3)$); + \path (5,3) [edge,-,red, thick] -- ($(6,2)!(5,3)!(4,3)$); + \path (4,3) [edge,-,red, thick] -- ($(5,2)!(4,3)!(3,3)$); + \path (3,3) [edge,-,red, thick] -- ($(4,2)!(3,3)!(2,3)$); + \path (2,3) [edge,-,red, thick] -- ($(3,2)!(2,3)!(1,3)$); + \path (1,3) [edge,-,red, thick] -- ($(2,2)!(1,3)!(0,3)$); \end{tikzpicture} \caption[A visual argument for $\delta^2 = \Delta \rho \norm{e}$]{ A visual argument showing that $\delta^2 = \Delta \rho \norm{e}$. If extended to the whole plane, there will be the same amount of blue - squares as red rectangles, so their areas must be equal. + squares as red rectangles as each grid point is the upper left + corner of both a blue and a red rectangle. Thus their areas must be equal. } \label{fig:area_proof} diff --git a/fig/bk_norm_subgraph.tex b/fig/bk_norm_subgraph.tex index 8237a3c..90fefb8 100644 --- a/fig/bk_norm_subgraph.tex +++ b/fig/bk_norm_subgraph.tex @@ -4,8 +4,8 @@ \node[vertex] (u) at (0, -1) {$u^\lambda_x$}; \node[vertex] (t) at (0, -2) {t}; - \path[edge] (s) -- node[weight noslope] {$\max\{F_{L-1}^x(1), 0\}$} (u); - \path[edge] (u) -- node[weight noslope] {$\max\{F_{L-1}^x(1), 0\} - F_\lambda^x(1)$} (t); + \path[edge] (s) -- node[weight noslope] {$\max\{F_{L-1}^x(1), 0\} - F_\lambda^x(1)$} (u); + \path[edge] (u) -- node[weight noslope] {$\max\{F_{L-1}^x(1), 0\}$} (t); \end{tikzpicture} \caption[Modified fidelity subgraph for the Boykov--Kolmogorov algorithm]{% diff --git a/fig/factory/contrast.tex b/fig/factory/contrast.tex index 7a6daae..915d697 100644 --- a/fig/factory/contrast.tex +++ b/fig/factory/contrast.tex @@ -1,7 +1,7 @@ % Title: glps_renderer figure % Creator: GL2PS 1.3.8, (C) 1999-2012 C. Geuzaine % For: Octave -% CreationDate: Fri Jan 30 17:04:57 2015 +% CreationDate: Tue Feb 3 19:56:32 2015 \begin{pgfpicture} \pgfsetlinewidth{0.01pt} \color[rgb]{1.000000,1.000000,1.000000} diff --git a/fig/line_disc.tex b/fig/line_disc.tex index 6358bbb..ebbd653 100644 --- a/fig/line_disc.tex +++ b/fig/line_disc.tex @@ -7,15 +7,24 @@ \node[tiny vertex] (\x\y) at (\x, \y) {}; } } - \node[tiny vertex,draw] (a) at (1, 2) {}; - \node[tiny vertex,draw] (b) at (1, 3) {}; - \node[tiny vertex,draw] (c) at (2, 2) {}; - \node[tiny vertex,draw] (d) at (2, 3) {}; - \node[tiny vertex,draw] (e) at (2, 4) {}; - \node[tiny vertex,draw] (f) at (3, 3) {}; - \node[tiny vertex,draw] (g) at (3, 1) {}; - \node[tiny vertex,draw] (h) at (2, 1) {}; - \node[tiny vertex,draw] (i) at (2, 0) {}; + \node[tiny vertex,draw] (a) at (2, 2) {}; + \node[tiny vertex,draw] (b) at (2, 3) {}; + \node[tiny vertex,draw] (c) at (3, 2) {}; + \node[tiny vertex,draw] (d) at (3, 3) {}; + \node[tiny vertex,draw] (e) at (3, 4) {}; + \node[tiny vertex,draw] (f) at (4, 3) {}; + \node[tiny vertex,draw] (g) at (4, 1) {}; + \node[tiny vertex,draw] (h) at (3, 1) {}; + \node[tiny vertex,draw] (i) at (3, 0) {}; + + \node[tiny vertex,draw] (j) at (2, 1) {}; + \node[tiny vertex,draw] (k) at (1, 2) {}; + \node[tiny vertex,draw] (l) at (1, 3) {}; + \node[tiny vertex,draw] (m) at (1, 4) {}; + \node[tiny vertex,draw] (n) at (1, 0) {}; + \node[tiny vertex,draw] (o) at (1, 1) {}; + \node[tiny vertex,draw] (p) at (0, 3) {}; + \node[tiny vertex,draw] (q) at (0, 1) {}; \path[edge,-] (a) -- (b); \path[edge,-] (a) -- (c); @@ -26,11 +35,20 @@ \path[edge,-] (a) -- (h); \path[edge,-] (a) -- (i); - \path[dash] (1, 2) -- +(35.785:4cm); - \path[dash] (1, 2) -- +(54.215:4cm); + \path[edge,-] (a) -- (j); + \path[edge,-] (a) -- (k); + \path[edge,-] (a) -- (l); + \path[edge,-] (a) -- (m); + \path[edge,-] (a) -- (n); + \path[edge,-] (a) -- (o); + \path[edge,-] (a) -- (p); + \path[edge,-] (a) -- (q); - \draw (1,2) +(35.785:3.1cm) arc (35.785:54.215:3.1cm); - \path (1,2) ++(45:3.4cm) node {$\Delta \phi$}; + \path[dash] (2, 2) -- +(35.785:4cm); + \path[dash] (2, 2) -- +(54.215:4cm); + + \draw (2, 2) +(35.785:3.1cm) arc (35.785:54.215:3.1cm); + \path (2, 2) ++(45:3.4cm) node {$\Delta \phi$}; %\path (00) [edge, <->, anchor=center] -- node[anchor=south] {h} (10) ; diff --git a/fig/line_midpoint.tex b/fig/line_midpoint.tex index 995a483..0cad458 100644 --- a/fig/line_midpoint.tex +++ b/fig/line_midpoint.tex @@ -9,29 +9,29 @@ \begin{scope} \clip (0, 0) rectangle (3, 2); - \path[nedge=4cm] (0,5) -- (2,4); - \path[nedge=4cm] (0,4) -- (2,3); - \path[nedge=4cm] (0,3) -- (2,2); - \path[nedge=4cm] (0,2) -- (2,1); - \path[nedge=4cm] (0,1) -- (2,0); + \path[nedge=4cm, thick] (0,5) -- (2,4); + \path[nedge=4cm, thick] (0,4) -- (2,3); + \path[nedge=4cm, thick] (0,3) -- (2,2); + \path[nedge=4cm, thick] (0,2) -- (2,1); + \path[nedge=4cm, thick] (0,1) -- (2,0); - \path[nedge=4cm] (-1,1) -- (1,0); + \path[nedge=4cm, thick] (-1,1) -- (1,0); - \path[nedge=4cm] (1,5) -- (3,4); - \path[nedge=4cm] (1,4) -- (3,3); - \path[nedge=4cm] (1,3) -- (3,2); - \path[nedge=4cm] (1,2) -- (3,1); - \path[nedge=4cm] (1,1) -- (3,0); + \path[nedge=4cm, thick] (1,5) -- (3,4); + \path[nedge=4cm, thick] (1,4) -- (3,3); + \path[nedge=4cm, thick] (1,3) -- (3,2); + \path[nedge=4cm, thick] (1,2) -- (3,1); + \path[nedge=4cm, thick] (1,1) -- (3,0); - \path[nedge=4cm] (2,5) -- (4,4); - \path[nedge=4cm] (3,5) -- (5,4); - \path[nedge=4cm] (4,5) -- (6,4); + \path[nedge=4cm, thick] (2,5) -- (4,4); + \path[nedge=4cm, thick] (3,5) -- (5,4); + \path[nedge=4cm, thick] (4,5) -- (6,4); \end{scope} - \path[dash] (2.5,1.5) -- (3.5,1); - \path[dash] ($(3.5,0.5)!(2.5,1.5)!(1.5,1.5)$) -- ($(4.5,0.0)!(3.5,1.0)!(2.5,1.0)$); + \path[dash, thick] (2.5,1.5) -- (3.5,1); + \path[dash, thick] ($(3.5,0.5)!(2.5,1.5)!(1.5,1.5)$) -- ($(4.5,0.0)!(3.5,1.0)!(2.5,1.0)$); - \path (3.5,1.0) [edge, <->, anchor=center] -- node[anchor=north + \path (3.5,1.0) [edge, thick, <->, anchor=center] -- node[anchor=north west,yshift=+3.0pt,xshift=+1pt, fill=white, font=\small, inner sep=1pt] {$\Delta \rho$} ($(4.5,0.0)!(3.5,1.0)!(2.5,1.0)$) ; diff --git a/fig/norm_evolution.tex b/fig/norm_evolution.tex index ed85765..704e44d 100644 --- a/fig/norm_evolution.tex +++ b/fig/norm_evolution.tex @@ -1,55 +1,24 @@ \centering -\begin{subfigure}[t]{0.4\textwidth} - \centering - \begin{tikzpicture} - \draw[->] (0,0) -- (4,0) node[right] {$\lambda$}; - \draw[->] (0,-2) -- (0,2) node[above] {$F_\lambda^x(1)$}; - \draw[line,domain=0:1] plot ({\x},{-1}); - \draw[line,domain=1:2,dashed] plot ({\x},{2*\x-3}); - \draw[line,domain=2:4] plot ({\x},{1}); +\begin{tikzpicture} + \draw[->] (0,0) -- (4,0) node[right] {$\lambda$}; + \draw[->] (0,-2) -- (0,2) node[above] {$F_\lambda^x(1)$}; + \draw[line,domain=0:4] plot ({\x},{.8*\x - 1.6}); - \draw (1,2pt) -- (1,-2pt) node[ - font=\small, - anchor=south - ] - {$f_x-1$}; + \draw (2.5,2pt) -- (2.5,-2pt) node[ + font=\small, + anchor=north + ] + {$f_x$}; - \draw (2,2pt) -- (2,-2pt) node[ - font=\small, - anchor=north - ] - {$f_x$}; + \draw (1.5,2pt) -- (1.5,-2pt) node[ + font=\small, + anchor=south + ] + {$f_x-1$}; - \end{tikzpicture} - \caption{$L^1$ fidelity term.} - \label{fig:l1_norm_evolution} -\end{subfigure} -~ -\begin{subfigure}[t]{0.4\textwidth} - \centering - \begin{tikzpicture} - \draw[->] (0,0) -- (4,0) node[right] {$\lambda$}; - \draw[->] (0,-2) -- (0,2) node[above] {$F_\lambda^x(1)$}; - \draw[line,domain=0:4] plot ({\x},{.8*\x - 1.6}); - - \draw (2.5,2pt) -- (2.5,-2pt) node[ - font=\small, - anchor=north - ] - {$f_x$}; - - \draw (1.5,2pt) -- (1.5,-2pt) node[ - font=\small, - anchor=south - ] - {$f_x-1$}; - - \end{tikzpicture} - \caption{$L^2$ fidelity term.} - \label{fig:l2_norm_evolution} -\end{subfigure} +\end{tikzpicture} \caption[Fidelity energy $F_\lambda^x(1)$ as a function of $\lambda$]{% - Two figures showing how the fidelity energy term $F_\lambda^x(1)$ + A sketch of how the fidelity energy term $F_\lambda^x(1)$ in \eqref{eq:final_discretization} increases monotonically with $\lambda$. } \label{fig:norm_evolution} diff --git a/fig/norm_subgraph.tex b/fig/norm_subgraph.tex index dbb9461..4ee8a91 100644 --- a/fig/norm_subgraph.tex +++ b/fig/norm_subgraph.tex @@ -6,7 +6,7 @@ \node[vertex] (u) at (0, -1) {$u^\lambda_x$}; \node[vertex] (t) at (0, -2) {t}; - \path[edge] (s) -- node[weight noslope] {$F_\lambda^x(1)$} (u); + \path[edge] (u) -- node[weight noslope] {$F_\lambda^x(1)$} (t); \end{tikzpicture} \caption{% The graph when $F_\lambda^x(1) > 0$, with constant equal to 0. @@ -21,7 +21,7 @@ \node[vertex] (u) at (0, -1) {$u^\lambda_x$}; \node[vertex] (t) at (0, -2) {t}; - \path[edge] (u) -- node[weight noslope] {$-F_\lambda^x(1)$} (t); + \path[edge] (s) -- node[weight noslope] {$-F_\lambda^x(1)$} (u); \end{tikzpicture} \caption{Graph when $F_\lambda^x(1) < 0$, with constant equal to $-F_\lambda^x(1)$. diff --git a/introduction.tex b/introduction.tex index e020ab8..4ac94fc 100644 --- a/introduction.tex +++ b/introduction.tex @@ -125,443 +125,3 @@ affect the restoration. %the result section will focus on looking at how these parameters affect %the restoration performance. -\chapter{Methods in image restoration} - -\fixme{Rating: 7/10} - -There are numerous methods in image restoration, and since we do not -have time nor space to discuss them all, this short overview will focus on the ones -related to the anisotropic total variation method considered later in -this thesis. See \cite{jahne} and \cite{aubert2006proc} for more -background on image processing in general. - -In this chapter, and also in the rest of the thesis we will assume that -we are given an image $f : \Omega \to \mathbb{R}$ where $\Omega$ is a -rectangular, open domain. Because of limitations in the numerical method -used, the codomain is $\mathbb{R}$ and we are thus restricted to -monochrome, or grayscale images. Such images are produced in large -numbers by for example ultrasound, X-ray and MRI machines. - -The space in which the image $f$ resides in will vary, but since we are -looking at image restoration methods, we assume that it includes some -kind of noise. Depending on the application and how the image is -obtained, one might construct different models describing different -types of noise. - -We will assume that the given image $f$ is a combination of an -underlying, actual image $u^*$, and some noise $\delta$. The simplest -model is additive noise where the assumption is that $f = u^* + \delta$. -There is also multiplicative noise where $f = u^* \cdot \delta$. An -other much seen noise type is salt and pepper noise, which is when black -and white pixels randomly appear in the image, independently of the -actual image $u^*$. - -These are only models, and in the real world the noise might be more -complex, and even come from a combination of sources. Depending on the -application, the goal might not even be to recover $u^*$, but rather to -obtain an output which fulfills certain smoothness or regularity -properties. In any case, we will continue using $f$ as the input and $u$ -as the output in the description of our methods. - -\section{Diffusion filtering} -\label{sec:diffusion_filtering} - -\fixme{Rating: 8/10} - -Diffusion filtering is a broad group of filtering and restoration methods -based on physical diffusion processes. The basic idea is to take the -noisy image as the initial value of some diffusion process, and then let it -evolve for some time. The best known method is probably the Gaussian -filter or Gaussian blur, in which one convolves the image with the -Gaussian function -\begin{equation} - K_\sigma(x,y) = \frac{1}{2\pi\sigma^2} \exp \left( - \frac{x^2 + - y^2}{2\sigma^2} \right). - \label{eq:gaussian_function} -\end{equation} -\nomenclature{$K_\sigma(x,y)$}{The Gaussian kernel.}% -In the discrete setting where the image consists of a grid of pixels, -the Gaussian blur amounts to calculating each pixel in the output image -as a weighted average of its neighboring pixels in the input image. - -The Gaussian function happens to be the fundamental solution of the heat -equation $\partial_t u = \Delta u$. Convolving $K_\sigma(x,y)$ with the -original image $f$ is therefore equivalent to solving the heat equation -with $f$ as initial value, until some time $T > 0$ depending on -$\sigma$. Boundary conditions have to be specified of course, and one -common choice is to symmetrically extend the image in $x$ and $y$ -directions, which corresponds to zero flux boundary conditions. - -By basic Fourier analysis it is possible to show that the Gaussian -filter is a low-pass filter which attenuates high frequencies. Further -theory can be found in Weickert's book on anisotropic diffusion -\cite{weickert1998anisotropic}. - -The main concern with the Gaussian filter is that it will, in addition -to smoothing out possible noise, remove details from the image, which -motivates the next set of methods, where the amount of diffusion can -vary for different parts of the image. - -\subsection{Non-linear diffusion filtering} - -\fixme{Rating: 8/10} - -In the theory of the heat equation one can introduce a \emph{thermal -diffusivity} $\alpha$ such that the equation becomes -\begin{equation} - \left\{ - \begin{aligned} - \partial_t u &= \diver \big( \alpha(\nabla u) \nabla u\big), \\ - u |_{t=0} &= f. - \end{aligned} - \right. -\end{equation} -\nomenclature{$\alpha(\nabla u)$}{Scalar thermal diffusivity.}% -The thermal diffusivity $\alpha(\nabla u) = \alpha(x, \nabla u)$ is material dependent, and can also -vary throughout the object. It specifies how well heat travels through -the specific point in the object. We can make use of this in the image -restoration context by controlling the diffusivity in different parts of -the image, where our goal is to reduce noise without loosing image -detail. Optimally, we would like there to be a lot of diffusion in -smooth parts of the image, and not so much in areas with a lot of -details. - -One much-studied non-linear diffusion equation is the Perona--Malik -equation -\begin{equation} - \partial_t u = \diver \left( \frac{\nabla u}{1 + - \frac{\abs{\nabla u}^2}{\lambda^2}} \right). - \label{eq:perona_malik} -\end{equation} -The thermal diffusivity $\alpha(\nabla u) = (1 + \sfrac{\abs{\nabla u}^2}{ -\lambda^2})^{-1}$ varies from $1$ in smooth areas to $0$ as the -norm of the gradient $\abs{\nabla u}$ grows. - -This particular form of the thermal diffusivity has been shown to be -related to how brightness is perceived by the human visual system. The -model has some theoretical problems related to well-posedness, for more -information see \cite{weickert1998anisotropic}. - -A different kind of non-linear diffusion model is the total variation -flow which can be formulated as -\begin{equation} - \partial_t u = \diver \frac{\nabla u}{\abs{\nabla u}}, -\end{equation} -where the diffusivity has a similar effect of reducing the diffusion in -areas of high variation. As the name suggests this model can be related -to the total variation formulation presented later. One forward Euler -time-step in the solution of this partial differential equation -corresponds to the Euler--Lagrange equation of the variational -formulation. - -Note that we follow Weickert's terminology when it comes to the -distinction between non-linear and anisotropic diffusion methods. The -Perona--Malik equation, and other diffusion equations with -non-homogenous diffusivities, are often by others called anisotropic, as -the diffusivity depends on the location. We will name these methods -non-linear and spare the anisotropy term for the ``real'' anisotropic -methods. These are methods where the diffusivity is a tensor, and thus -both location and direction dependent. - -\subsection{Anisotropic diffusion} -\label{sec:anisotropic_diffusion} - -The diffusivity is made directionally dependent by introducing a -diffusion \emph{tensor} $A(u)$ such that the initial boundary value -problem becomes -\begin{equation} - \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.}% -where $\nu$ is the outer normal of $\Omega$. The tensor $A(u)$ -is constructed such as to diminish the effect of $\nabla -u$ across what we believe to be edges in the image. This way, there will -also be less diffusion through these edges. Weickert -\cite{weickert1998anisotropic} suggests constructing $A(u)$ based on the -edge estimator $\nabla u_\sigma$ where -\nomenclature{$\nabla u_\sigma$}{Edge estimator.}% -\begin{equation} - u_\sigma := K_\sigma * \tilde{u} -\end{equation} -\nomenclature{$u_\sigma$}{Image $u$ smoothed with a Gaussian of -parameter $\sigma$.}% -and $\tilde{u}$ is an extension of $u$ from $\Omega$ to $\mathbb{R}^2$ -made by symmetrically extending $u$ across the edges. -Assuming we are at an edge in the image, the -direction of $\nabla u_\sigma$ should be perpendicular to the edge, -while its magnitude will provide information on the steepness of the -edge. - -To extract this information, and also to identify features on a larger -scale, the \emph{structure tensor} is introduced -\begin{equation} - S_\rho(x) := K_\rho * (\nabla u_\sigma \otimes \nabla u_\sigma), -\end{equation} -\nomenclature{$S_\rho(x)$}{Structure tensor.}% -where the convolution with the Gaussian function $K_\rho$ is done -component-wise. The anisotropy tensor $A(u)$ can then be constructed -based on the eigenvectors and eigenvalues of $S_\rho(x)$. -The structure tensor and its properties will be discussed further when -we introduce our anisotropic total variation functional. - -Assuming some -smoothness, symmetry and uniform positive definiteness on $A(u)$ one can -prove well-posedness, regularity and an extremum principle of the -problem \eqref{eq:aniso_diff} as done in \cite{weickert1998anisotropic}. - -Even if the diffusivity tensor was introduced to reduce the amount of -smoothing across edges, the solution of \eqref{eq:aniso_diff} will still be -infinitely differentiable \cite{weickert1998anisotropic}, i.e.\ $u(T) \in -C^\infty(\Omega)$ for $T > 0$. Thus we have no real discontinuities, and -no real edges in our solution. - -Further, the anisotropic diffusion may introduce structure based on -noise, when there really was no structure to begin with. This is a -problem we hope to avoid in our anisotropic total variation method. - -\begin{figure} - \centering - \begin{subfigure}[t]{0.4\textwidth} - \centering - \includegraphics[width=\textwidth]{fig/fnoise.png} - \caption{Noisy, $\sigma > 9000$} - \end{subfigure} - ~ - \begin{subfigure}[t]{0.4\textwidth} - \centering - \includegraphics[width=\textwidth]{fig/fanisodiff.png} - \caption{Regularized, parameters: bleep bloop} - \end{subfigure} - \caption[bleapofsijdosij]{\fixme{UPDATE THIS}} - \label{fig:ansi_diff} -\end{figure} - -\section{Total variation} -\label{sec:total_variation} - -\fixme{Rating: 8/10, numerical methods: 5/10} - -The total variation image restoration method is usually formulated as a -minimization problem -\begin{equation} - \begin{gathered} - \min_{u \in L^p(\Omega)} F(u), \\ - F(u) = \underbrace{\int_\Omega \abs{u - f}^p \, dx}_{\text{fidelity - term}} {}+{} \beta \underbrace{\int_\Omega - \abs{\nabla u} \, dx}_{\mathclap{\text{regularization - term}}}, - \label{eq:first_min_presentation} - \end{gathered} -\end{equation} -where $p$ is normally taken to be 1 or 2. -The fidelity term penalizes images $u$ that are far from the original -image $f$, and thus controls the \emph{fidelity} of our solution. The -regularization term is the total variation, and minimizing it will -reduce the variation and thus regularize the image. The $\beta$ parameter -controls the strength of the regularization. - -Since we do not only want to consider differentiable images $u \in -C^1(\Omega)$ for which the gradient exists, we introduce the total -variation using the distributional derivative. -\begin{definition}[Total variation] - Given a function $u \in L^1(\Omega)$, the total variation of $u$, - often written $\int_\Omega \abs{Du}\, dx$, where the $D$ is the - gradient taken in the distributional sense, is - \begin{equation} - \TV(u) - = \int_\Omega \abs{Du} \, dx - = \sup \left\{ \int_\Omega u \cdot \diver \varphi \, dx : - \varphi \in C^\infty_c\left(\Omega, \mathbb{R}^2\right), - \norm{\varphi}_{L^\infty(\Omega)} \leq 1 - \right\}. - \end{equation} - The test functions $\varphi$ are taken from - $C^\infty_c\left(\Omega, \mathbb{R}^2\right)$, the space of smooth - functions from $\Omega$ to $\mathbb{R}^2$ with compact support in - $\Omega$. - \label{def:tv} -\end{definition} -\nomenclature{$\TV(u)$}{Total variation of image $u$.}% - -Note that since $\Omega$ is open and bounded, the test functions -$\varphi$ vanish on the boundary of $\Omega$. Thus no variation is -measured at the boundary. As this restoration method is the one -which will be extended later in this thesis, we will look a little bit -more deeply into the background and the numerical methods relating to -it. - -As we are searching for an image with low total variation, it is useful -to introduce the space of functions of bounded variation. -\begin{definition}[Functions of bounded variation] - The space of functions of bounded variation $\BV(\Omega)$ is the - space of functions $u \in L^1(\Omega)$ for which the total variation - is finite, i.e., - \begin{equation} - \BV(\Omega) = \left\{ u \in L^1(\Omega) : \TV(u) < \infty - \right\}. - \end{equation} -\end{definition} -\nomenclature{$\BV(\Omega)$}{Functions of bounded variation in -$\Omega$.}% -Our optimization problem has thus become -\begin{equation} - \min_{u \in \BV(\Omega)} \int_\Omega \abs{u - v}^p \, dx + \beta \, - \TV(u). - \label{eq:second_min_presentation} -\end{equation} - -As with any restoration method, the total variation method has its -strengths and weaknesses. Its main strength is its ability to recover -edges in the input image. The total variation of a section only takes -the absolute change into account, and does not favor gradual changes -like the diffusion methods. - -There is also a theoretical result stating that the set of edges in the -solution $u$ is contained in the set of edges in the original image $f$, -thus no new edges are created \cite{caselles2011total}. However, in the -presence of noise, the method may introduce or rather ``find'' new edges -that were not in the original image, since flat sections of zero -variation are encouraged by the functional. This effect is called the -\emph{stair-casing effect}, and can be seen in Figure~\ref{fig:grad_tv} -where a noisy gradient has been restored using total variation -restoration. - -\begin{figure} - \centering - \begin{subfigure}[t]{0.4\textwidth} - \centering - \includegraphics[width=\textwidth]{fig/factory/grad/n_q40.png} - \caption{Noisy gradient} - \end{subfigure} - ~ - \begin{subfigure}[t]{0.4\textwidth} - \centering - \includegraphics[width=\textwidth]{fig/factory/grad/r_p2_n16_b5000_g1000000000_r5_s3.png} - \caption{Total variation restoration} - \end{subfigure} - \caption[The stair-casing effect of regular total variation - restoration]{% - Although the original gradient was smooth, the total variation - method manages to find structure in the noise, and create edges - in the restored image. - } - \label{fig:grad_tv} -\end{figure} - -Thin objects and corners may also suffer from contrast loss since -bringing them closer to their surroundings reduces the -total variation. An example of this is shown in -Figure~\ref{fig:tv_example}, where a not particularly noisy fingerprint -image has been strongly regularized. The original black and white levels -have been brought closer to yield a lower total variation in the -regularized image. - -\begin{figure} - \centering - \begin{subfigure}[t]{0.4\textwidth} - \centering - \includegraphics[width=\textwidth]{fig/finger.png} - \end{subfigure} - ~ - \begin{subfigure}[t]{0.4\textwidth} - \centering - \includegraphics[width=\textwidth]{fig/ftv.png} - \end{subfigure} - \caption[Heavily regularized fingerprint showing contrast loss]{% - A fingerprint heavily regularized using the total variation - method. The originally white and black ridges have been brought - closer in value, to reduce the total variation. - } - \label{fig:tv_example} -\end{figure} - -\subsection{Numerical methods} - -See \cite{chan2011numerical} for an overview of some of the numerical -methods relating to total variation image restoration. Amongst others it -describes some dual and primal-dual methods, as well as the graph cut -approach we take in this thesis. - -\subsubsection{Graph cut approach} - -\fixme{Could probably be shortened a lot, but I thought it could be -useful to at least give an idea how it is done.} - -Using graph cuts is the approach we will be taking later when -considering the anisotropic total variation regularization, and it is -therefore valuable to briefly look into how graph cuts are used in the -case of regular total variation. - -In the theory of graphs, a graph cut is a set of edges that when -removed will separate the graph into two disconnected parts. A minimum -cut is a cut such that the sum of the weights is minimal. It -has been shown that using special graph constructions, the graph cuts -can be used to minimize a special type of energy functions. - -In the discrete setting our image consists of pixels, and is represented -by a function $u : \mathcal{G} \to \mathcal{P}$ where $\mathcal{G}$ is a -regular grid over $\Omega$, and $\mathcal{P} = \{0, \hdots, L-1\}$ is -the discrete set of pixel values, or \emph{levels}. - -The idea of the graph cut approach is to decompose the minimization -problem into one minimization problem for each level of the image, and -then solve them separately before combining the results. - -For an image $u$ and a level $\lambda$ we denote the \emph{level set} by -$\{ u > \lambda\}$, defined as the set $\{ x \in \Omega : u_x > \lambda -\}$. The thresholded image $u^\lambda$, an indicator function, is then -defined as -\begin{equation} - u^\lambda = \idfun_{u > \lambda}. -\end{equation} -Here, $\idfun_E$ signifies the characteristic function of the set $E$, -\nomenclature{$\idfun_E$}{Characteristic function of the set $E$.}% -the function which is equal to one in every point in $E$, and zero -elsewhere. The graph cut we find later will for each level $\lambda$ -give us the thresholded image $u^\lambda$ which then can be combined to -form the complete image $u$. - -Through careful manipulation of the continuous functional in -\eqref{eq:second_min_presentation} it is possible to obtain a discrete -functional decomposed as a sum over all the level values on the form -\begin{equation} - F(u) = - \sum_{\lambda=0}^{L-2} \sum_x F^x_\lambda(u^\lambda_x) - + \beta \sum_{\lambda = 0}^{L-2} \sum_{(x, y)} F^{x,y}(u^\lambda_x, - u^\lambda_y) - =: \sum_{\lambda=0}^{L-2} F_\lambda(u^\lambda) - \label{eq:old_total_energy} -\end{equation} -where the sum over $(x, y)$ is over all pixel pairs $(x, y)$ in a neighbor -relation, i.e.\ where the pixels are ``close'' to each other. - -The actual form of the functional, and the steps to construct it will be -presented later, but hopefully, it does not seem too infeasible. And for -those with some knowledge of graph theory, it might be valuable to -consider the intuition behind the graph cut approach. - -When constructing the graph used to find the thresholded image -$u^\lambda$, we have two special vertices, one representing the set $\{ -u > \lambda\}$, while the other represents the set $\{ u \leq -\lambda\}$. The pixels are then connected to these vertices with a -weight representing how strongly they are related to the corresponding -set. This weight will be based on the value of $F^x_\lambda$. - -Then there are connections between pixels in a neighborhood relation, -representing the energy $F^{x,y}$. Thus when finding a cut, we partition -the pixels into the sets $\{ u > \lambda \}$ and $\{ u \leq \lambda \}$. -And if in addition the cut is minimal, we know that the edges cut have -minimal weight, and can prove that the $u^\lambda$ produced minimizes -the functional in \eqref{eq:total_energy}. - diff --git a/main.tex b/main.tex index e08df0b..80bc52a 100644 --- a/main.tex +++ b/main.tex @@ -75,10 +75,11 @@ \usepackage{algpseudocode} \newtheorem{theorem}{Theorem} -\newtheorem{proposition}{Proposition} -\newtheorem{lemma}{Lemma} -\newtheorem{definition}{Definition} - +\numberwithin{theorem}{chapter} +\newtheorem{proposition}[theorem]{Proposition} +\newtheorem{lemma}[theorem]{Lemma} +\newtheorem{definition}[theorem]{Definition} + % Indeksgenerering er kommentert ut her. Ta bort prosenttegnene % hvis du vil ha en indeks: %\usepackage{makeidx} @@ -200,6 +201,7 @@ Bjørn Rustad, \today. % Inkluder alle de andre kildefilene \input{introduction} +\input{methods} \input{continuous} \input{discrete} \input{maxflow} diff --git a/maxflow.tex b/maxflow.tex index 378aa8c..e54e92d 100644 --- a/maxflow.tex +++ b/maxflow.tex @@ -6,8 +6,8 @@ minimizing the functional for one level value $\lambda$. We will in this and the next section see how such a minimum cut can be found by sending flow through the graph and trying to identify the ``bottleneck''. This chapter, except for the description of the -Boykov--Kolmogorov algorithm is taken from my project work -\cite{project} and is included here for completeness. +Boykov--Kolmogorov algorithm is taken with some adjustments from my +project work \cite{project} and is included here for completeness. \section{Flow graphs} @@ -89,12 +89,13 @@ cuts are minimum cuts. %graph flow theory, stating that the capacity of a minimum $s$-$t$-cut %in a graph, is equal to the maximum flow from the source to the sink. -But how does this help us? We know that if we know the maximum flow +But how does this help us? If we know the maximum flow value, and we have an $s$-$t$-cut with capacity equal to the maximum flow, we actually have a minimum cut. The question is then, how do we find a maximum flow, and how do we find a minimum cut? \section{Augmenting path algorithms} + The family of augmenting flow algorithms represent a popular approach to the maximum flow problem. The idea is simply to look for paths from the source to the sink that can carry additional flow, so-called augmenting @@ -103,6 +104,7 @@ path. When no such path exists anymore, no more flow can be sent from the source to the sink, and a maximum flow has been reached. \subsection{Residual graph} + When further discussing approaches to solving the maximum flow problem we will need the notion of a residual graph $G_f = (V_f, E_f, c_f)$, which is derived from the original graph $G$ and contains the edges @@ -251,7 +253,7 @@ As in most of the cited push-relabel literature, we define $N = \begin{equation} e(u) = \sum_{v \in V} f(v, u) - \sum_{v \in V} f(u, v), \end{equation} -which represents the amount of flow which \emph{disappears} in vertex +which represents the amount of flow which disappears in vertex $u$. Equivalent to the preflow conservation constraint is stating that $e(u) \geq 0$ for all vertices $u \in V - \{s,t\}$. @@ -495,7 +497,7 @@ the end of the queue. The original article of Goldberg and Tarjan complexity of $O(\abs{V}^3)$, which can be very good if you have a dense graph. -Using highest level selection rule one always discharges the active +Using the highest level selection rule one always discharges the active vertex with the largest distance label. Goldberg and Tarjan state that this rule also gives a complexity of $O(\abs{V}^3)$ while this bound is improved to $O(\abs{V}^2 \sqrt{\abs{E}})$ in an article by Cheriyan and @@ -541,7 +543,7 @@ The preflow and capacity constraints are still valid after the gap relabeling procedure, as only the labels $d$ are changed. For the labeling $d$ one has to verify that $d(u) \leq d(v) + 1$ for every edge $(u,v) \in E_f$ in the residual network. If none, both, or only $v$ is -relabeled, this is trivial. It is not possible that only $v$ would be +relabeled, this is trivial. It is not possible that only $u$ would be relabeled, as this would imply that $d(u) \geq d(v) + 2$ which is not a valid initial labeling. @@ -670,7 +672,7 @@ we do not change the labeling $d$, so the labeling constraints are also satisfied in the new graph. Through this procedure we have easily created the graph for $\lambda = -k-1$, and the distance labels remain the same. As these labels always +k-1$, and the distance labels remain the same. As the distance labels always increase monotonically, we have a head start compared to the case where we reset the flow and labels. @@ -693,10 +695,10 @@ inclusion \{ u \in V : d^\lambda(u) \geq N \} \subseteq \{ u \in V : d^{\lambda-1}(u) \geq N \}. \end{equation} -For a pixel $x \in S$ we will set $u^\lambda_x = 0$, which together with +For a pixel $x \in S$ we will set $u^\lambda_x = 1$, which together with the previous inclusion property implies \begin{equation} - u^\lambda_x \geq u^{\lambda-1}_x + u^\lambda_x \leq u^{\lambda-1}_x \end{equation} for all $x \in \mathcal{G}$. Being equivalent with the inequality in \eqref{eq:decreasing_sets}, this means our algorithm produces stackable @@ -704,13 +706,13 @@ thresholded images $u^\lambda$. We then construct our output image $u$ by giving each pixel the value \begin{equation} - u_x = \min \{ \lambda \in \{0, \ldots, L-1\} : u^\lambda_x = 1 \}. + u_x = \min \{ \lambda \in \{0, \ldots, L-1\} : u^\lambda_x = 0 \}. \end{equation} This marks the end of the description of the implemented algorithm, but -we will further discuss some possible improvements, and also look at the -results when using the method on different kinds of noisy images. +we will further discuss some possible improvements. \subsection{Divide and conquer} + The possibility of re-using the graph between separate level is a very nice property of the push-relabel algorithm, but there are further room for improvements. Consider one pixel $x$ with value $u_x$, @@ -984,7 +986,7 @@ A proof of correctness can be found in \fixme{kolmogorov-phd}. \subsection{Graph reuse} -As in the push-relabel algorithm described earlier, it is not necessary +As with 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. @@ -1005,18 +1007,43 @@ graph in Figure~\ref{fig:norm_subgraph_neg}. As $\lambda$ goes from $L-1$ to $0$ stay non-negative. The edge from $u_\lambda^x$ to $t$ is non-decreasing with decreasing $\lambda$ parameter. +\begin{table}[b] + \centering + \caption[Graph construction for $F_\lambda^x(u^\lambda_x)$ in the + Boykov--Kolmogorov algorithm]{ + Each row represents one of the two possible values of + $u^\lambda_x \in \{0,1\}$. The functional $F_\lambda^x(u^\lambda_x)$ + and minimum cut obtaining this configuration is shown. The last + column shows the capacities of the cut for the graph + construction in Figure~\ref{fig:bk_norm_subgraph}. We verify + that the cut capacities are equal to the functional value, plus + a constant. + } + \begin{tabular}{cccl} + \hline + $u^\lambda_x$ & $F_\lambda^x(u^\lambda_x)$ & Min.\ cut $(S, T)$ + & Graph cut cap. \\ + \hline + $0$ & $0$ & $(\{s\}, \{u^\lambda_x, t\})$ & $\max \{ + F_{L-1}^x(1), 0\} - F_\lambda^x(1)$ \\ + $1$ & $F_\lambda^x(1)$ & $(\{s, u^\lambda_x\}, \{t\})$ & + $\max \{F_{L-1}^x(1), 0\}$ + \end{tabular} + \label{tab:fid_energy} +\end{table} + This construction allows us to update the capacities of the edges, while retaining the flow and the trees $S$ and $T$. A nice property that -follows is that the partition $(V - T, T)$ is after each run a minimum -cut since there is no more residual capacity between the trees $S$ -and $T$. Further, no vertex already in the tree $T$ will leave $T$ when +follows is that the partition $(S, V-S)$ is after each run a minimum +cut since the tree $S$ has been grown as much as possible. +Further, no vertex already in the tree $S$ will leave $S$ when the trees are updated and algorithm is run for lower $\lambda$ values. -This can be seen from the fact that if $u$ is in $T$ and the algorithm +This can be seen from the fact that if $u$ is in $S$ and the algorithm has terminated, all paths from $s$ to $t$ going -through $u$ have a saturated edge somewhere before $u$. For all these -paths, the capacity will only change for edges after $u$, as we only -change the capacity of edges $(v, t)$. Thus these saturated edges will -stay saturated, and $u$ will stay in $T$. \fixme{shaky} +through $u$ have a saturated edge somewhere after $u$. For all these +paths, the capacity will only change for edges before $u$, as we only +change the capacity of edges $(s, v)$. Thus these saturated edges will +stay saturated, and $u$ will stay in $S$. \subsection{Performance improvements} diff --git a/methods.tex b/methods.tex new file mode 100644 index 0000000..8524664 --- /dev/null +++ b/methods.tex @@ -0,0 +1,446 @@ +\chapter{Methods in image restoration} + +\fixme{Rating: 7/10} + +There are numerous methods in image restoration, and since we do not +have time nor space to discuss them all, this short overview will focus on the ones +related to the anisotropic total variation method considered later in +this thesis. See \cite{jahne} and \cite{aubert2006proc} for more +background on image processing in general. + +In this chapter, and also in the rest of the thesis we will assume that +we are given an image $f : \Omega \to \mathbb{R}$ where $\Omega$ is a +rectangular, open domain. Because of limitations in the numerical method +used, the codomain is $\mathbb{R}$ and we are thus restricted to +monochrome, or grayscale images. Such images are produced in large +numbers by for example ultrasound, X-ray and MRI machines. + +The space in which the image $f$ resides in will vary, but since we are +looking at image restoration methods, we assume that it includes some +kind of noise. Depending on the application and how the image is +obtained, one might construct different models describing different +types of noise. + +We will assume that the given image $f$ is a combination of an +underlying, actual image $u^*$, and some noise $\delta$. The simplest +model is additive noise where the assumption is that $f = u^* + \delta$. +There is also multiplicative noise where $f = u^* \cdot \delta$. An +other much seen noise type is salt and pepper noise, which is when black +and white pixels randomly appear in the image. + +These are only models, and in the real world the noise might be more +complex, and even come from a combination of sources. Depending on the +application, the goal might not even be to recover $u^*$, but rather to +obtain an output which fulfills certain smoothness or regularity +properties. In any case, we will continue using $f$ for the noisy input +image and use $u$ for the output image in the description of the +restoration methods. + +\section{Diffusion filtering} +\label{sec:diffusion_filtering} + +\fixme{Rating: 8/10} + +Diffusion filtering is a broad group of filtering and restoration methods +based on physical diffusion processes. The basic idea is to take the +noisy image as the initial value of some diffusion process, and then let it +evolve for some time. The best known method is probably the Gaussian +filter or Gaussian blur, in which one convolves the image with the +Gaussian function +\begin{equation} + K_\sigma(x,y) = \frac{1}{2\pi\sigma^2} \exp \left( - \frac{x^2 + + y^2}{2\sigma^2} \right). + \label{eq:gaussian_function} +\end{equation} +\nomenclature{$K_\sigma(x,y)$}{The Gaussian kernel.}% +In the discrete setting where the image consists of a grid of pixels, +the Gaussian blur amounts to calculating each pixel in the output image +as a weighted average of its neighboring pixels in the input image. + +The Gaussian function happens to be the fundamental solution of the heat +equation $\partial_t u = \Delta u$. Convolving $K_\sigma(x,y)$ with the +original image $f$ is therefore equivalent to solving the heat equation +with $f$ as initial value, until some time $T > 0$ depending on +$\sigma$. Boundary conditions have to be specified of course, and one +common choice is to symmetrically extend the image in $x$ and $y$ +directions, which corresponds to zero flux boundary conditions. + +By basic Fourier analysis it is possible to show that the Gaussian +filter is a low-pass filter which attenuates high frequencies. Further +theory can be found in Weickert's book on anisotropic diffusion +\cite{weickert1998anisotropic}. + +The main concern with the Gaussian filter is that it will, in addition +to smoothing out possible noise, remove details from the image, which +motivates the next set of methods, where the amount of diffusion can +vary between different parts of the image. + +\subsection{Non-linear diffusion filtering} + +\fixme{Rating: 8/10} + +In the theory of the heat equation one can introduce a \emph{thermal +diffusivity} $\alpha$ such that the equation becomes +\begin{equation} + \left\{ + \begin{aligned} + \partial_t u &= \diver \big( \alpha(\nabla u) \nabla u\big), \\ + u |_{t=0} &= f. + \end{aligned} + \right. +\end{equation} +\nomenclature{$\alpha(\nabla u)$}{Scalar thermal diffusivity.}% +The thermal diffusivity $\alpha(\nabla u) = \alpha(x, \nabla u)$ is material dependent, and can also +vary throughout the object. It specifies how well heat travels through +the specific point in the object. We can make use of this in the image +restoration context specifying different diffusivity in different parts of +the image, in an effort to reduce noise without loosing image +detail. Optimally, we would like there to be a lot of diffusion in +smooth parts of the image, and not so much in areas with a lot of +details. + +One much-studied non-linear diffusion equation is the Perona--Malik +equation +\begin{equation} + \partial_t u = \diver \left( \frac{\nabla u}{1 + + \frac{\abs{\nabla u}^2}{\lambda^2}} \right). + \label{eq:perona_malik} +\end{equation} +The thermal diffusivity $\alpha(\nabla u) = (1 + \sfrac{\abs{\nabla u}^2}{ +\lambda^2})^{-1}$ varies from $1$ in smooth areas to $0$ as the +norm of the gradient $\abs{\nabla u}$ grows. + +This particular form of the thermal diffusivity has been shown to be +related to how brightness is perceived by the human visual system. The +model has some theoretical problems related to well-posedness, for more +information see \cite{weickert1998anisotropic}. + +A different kind of non-linear diffusion model is the total variation +flow which can be formulated as +\begin{equation} + \partial_t u = \diver \frac{\nabla u}{\abs{\nabla u}}, +\end{equation} +where the diffusivity has a similar effect of reducing the diffusion in +areas of high variation. As the name suggests this model can be related +to the total variation formulation presented later. One forward Euler +time-step in the solution of this partial differential equation +corresponds to the Euler--Lagrange equation of the variational +formulation. + +Note that we follow Weickert's terminology when it comes to the +distinction between non-linear and anisotropic diffusion methods. The +Perona--Malik equation, and other diffusion equations with +non-homogenous diffusivities, are often by others called anisotropic, as +the diffusivity depends on the location. We will name these methods +non-linear and spare the anisotropy term for the ``real'' anisotropic +methods. These are methods where the diffusivity is a tensor, and thus +both location and direction dependent. + +\subsection{Anisotropic diffusion} +\label{sec:anisotropic_diffusion} + +The diffusivity is made directionally dependent by introducing a +diffusion \emph{tensor} $A(u)$ such that the initial boundary value +problem becomes +\begin{equation} + \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.}% +where $\nu$ is the outer normal of $\Omega$. The tensor $A(u)$ +is constructed such as to diminish the effect of $\nabla +u$ across what we believe to be edges in the image. This way, there will +also be less diffusion through these edges. Weickert +\cite{weickert1998anisotropic} suggests constructing $A(u)$ based on the +edge estimator $\nabla u_\sigma$ where +\nomenclature{$\nabla u_\sigma$}{Edge estimator.}% +\begin{equation} + u_\sigma := K_\sigma * \tilde{u} +\end{equation} +\nomenclature{$u_\sigma$}{Image $u$ smoothed with a Gaussian of +parameter $\sigma$.}% +and $\tilde{u}$ is an extension of $u$ from $\Omega$ to $\mathbb{R}^2$ +made by symmetrically extending $u$ across the edges. +Assuming we are at an edge in the image, the +direction of $\nabla u_\sigma$ should be perpendicular to the edge, +while its magnitude will provide information on the steepness of the +edge. + +To extract this information, and also to identify features on a larger +scale, the \emph{structure tensor} is introduced +\begin{equation} + S_\rho(x) := K_\rho * (\nabla u_\sigma \otimes \nabla u_\sigma), +\end{equation} +\nomenclature{$S_\rho(x)$}{Structure tensor.}% +where the convolution with the Gaussian function $K_\rho$ is done +component-wise. The anisotropy tensor $A(u)$ can then be constructed +based on the eigenvectors and eigenvalues of $S_\rho(x)$. +The structure tensor and its properties will be discussed further when +we introduce our anisotropic total variation functional. + +Assuming some +smoothness, symmetry and uniform positive definiteness on $A(u)$ one can +prove well-posedness, regularity and an extremum principle of the +problem \eqref{eq:aniso_diff} as done in \cite{weickert1998anisotropic}. + +Even if the diffusivity tensor was introduced to reduce the amount of +smoothing across edges, the solution of \eqref{eq:aniso_diff} will still be +infinitely differentiable \cite{weickert1998anisotropic}, i.e.\ $u(T) \in +C^\infty(\Omega)$ for $T > 0$. Thus there are no real discontinuities, and +no real edges in the solution. + +Further, the anisotropic diffusion may introduce structure based on +noise, when there really was no structure to begin with. This is a +problem we aim to avoid in our anisotropic total variation method. + +%\begin{figure} +% \centering +% \begin{subfigure}[t]{0.4\textwidth} +% \centering +% \includegraphics[width=\textwidth]{fig/fnoise.png} +% \caption{Noisy, $\sigma > 9000$} +% \end{subfigure} +% ~ +% \begin{subfigure}[t]{0.4\textwidth} +% \centering +% \includegraphics[width=\textwidth]{fig/fanisodiff.png} +% \caption{Regularized, parameters: bleep bloop} +% \end{subfigure} +% \caption[bleapofsijdosij]{\fixme{UPDATE THIS}} +% \label{fig:ansi_diff} +%\end{figure} + +\section{Total variation} +\label{sec:total_variation} + +\fixme{Rating: 8/10, numerical methods: 5/10} + +Total variation was initially introduced to the field of image +restoration by Rudin, Osher and Fatemi in \cite{rudin1992nonlinear} and +is usually formulated as a minimization problem +\begin{equation} + \begin{gathered} + \min_{u \in L^p(\Omega)} F(u), \\ + F(u) = \underbrace{\int_\Omega \abs{u - f}^p \, dx}_{\text{fidelity + term}} {}+{} \beta \underbrace{\int_\Omega + \abs{\nabla u} \, dx}_{\mathclap{\text{regularization + term}}}, + \label{eq:first_min_presentation} + \end{gathered} +\end{equation} +where $p$ is normally taken to be 1 or 2. +The fidelity term penalizes images $u$ that are far from the original +image $f$. The +regularization term is the total variation of the image, and minimizing it will +reduce the variation and thus regularize the image. The $\beta$ parameter +controls the strength of the regularization. Note that $u = f$ is a +minimizer of the fidelity term, while a constant image $u = c$ is a +minimizer of the regularization term. + +As this restoration method is the one +which will be extended later in this thesis, we will look a little bit +more deeply into the background and the numerical methods relating to +it. + +Since we do not only want to consider differentiable images $u \in +C^1(\Omega)$ for which the gradient exists, we introduce the total +variation using the distributional derivative. +\begin{definition}[Total variation] + Given a function $u \in L^1(\Omega)$, the total variation of $u$, + often written $\int_\Omega \abs{Du}\, dx$, where the $D$ is the + gradient taken in the distributional sense, is + \begin{equation} + \TV(u) + = \int_\Omega \abs{Du} \, dx + = \sup \left\{ \int_\Omega u \cdot \diver \varphi \, dx : + \varphi \in C^\infty_c\left(\Omega, \mathbb{R}^2\right), + \norm{\varphi}_{L^\infty(\Omega)} \leq 1 + \right\}. + \end{equation} + The test functions $\varphi$ are taken from + $C^\infty_c\left(\Omega, \mathbb{R}^2\right)$, the space of smooth + functions from $\Omega$ to $\mathbb{R}^2$ with compact support in + $\Omega$. + \label{def:tv} +\end{definition} +\nomenclature{$\TV(u)$}{Total variation of image $u$.}% + +Note that since $\Omega$ is open and bounded, the test functions +$\varphi$ vanish on the boundary of $\Omega$. Thus no variation is +measured at the boundary. + +As we are searching for an image with low total variation, it is useful +to introduce the space of functions of bounded variation. +\begin{definition}[Functions of bounded variation] + The space of functions of bounded variation $\BV(\Omega)$ is the + space of functions $u \in L^1(\Omega)$ for which the total variation + is finite, i.e., + \begin{equation} + \BV(\Omega) = \left\{ u \in L^1(\Omega) : \TV(u) < \infty + \right\}. + \end{equation} +\end{definition} +\nomenclature{$\BV(\Omega)$}{Functions of bounded variation in +$\Omega$.}% +Our optimization problem has thus become +\begin{equation} + \min_{u \in \BV(\Omega)} \int_\Omega \abs{u - v}^p \, dx + \beta \, + \TV(u). + \label{eq:second_min_presentation} +\end{equation} + +As with any restoration method, the total variation method has its +strengths and weaknesses. Its main strength is its ability to recover +edges in the input image. The total variation of a section only takes +the absolute change into account, and does not favor gradual changes +like the diffusion methods. + +There is also a theoretical result stating that the set of edges in the +solution $u$ is contained in the set of edges in the original image $f$, +thus no new edges are created \cite{caselles2011total}. However, in the +presence of noise, the method may introduce or rather ``find'' new edges +that were not in the original image, since flat sections of zero +variation are encouraged by the functional. This effect is called the +\emph{stair-casing effect}, and can be seen in Figure~\ref{fig:grad_tv} +where a noisy gradient has been restored using total variation +restoration. + +\begin{figure} + \centering + \begin{subfigure}[t]{0.4\textwidth} + \centering + \includegraphics[width=\textwidth]{fig/factory/grad/n_q40.png} + \caption{Noisy gradient} + \end{subfigure} + ~ + \begin{subfigure}[t]{0.4\textwidth} + \centering + \includegraphics[width=\textwidth]{fig/factory/grad/r_p2_n16_b5000_g1000000000_r5_s3.png} + \caption{Total variation restoration} + \end{subfigure} + \caption[The stair-casing effect of regular total variation + restoration]{% + Although the original gradient was smooth, the total variation + method manages to find structure in the noise, and create edges + in the restored image. + } + \label{fig:grad_tv} +\end{figure} + +Fine details, thin objects and corners may suffer from contrast loss since +bringing them closer to their surroundings reduces the +total variation. An example of this is shown in +Figure~\ref{fig:tv_example}, where a not particularly noisy fingerprint +image has been strongly regularized. The original black and white levels +have been brought closer to yield a lower total variation in the +regularized image. + +\begin{figure} + \centering + \begin{subfigure}[t]{0.4\textwidth} + \centering + \includegraphics[width=\textwidth]{fig/finger.png} + \end{subfigure} + ~ + \begin{subfigure}[t]{0.4\textwidth} + \centering + \includegraphics[width=\textwidth]{fig/ftv.png} + \end{subfigure} + \caption[Heavily regularized fingerprint showing contrast loss]{% + A fingerprint heavily regularized using the total variation + method. The originally white and black ridges have been brought + closer in value, to reduce the total variation. + } + \label{fig:tv_example} +\end{figure} + +\subsection{Numerical methods} + +See \cite{chan2011numerical} for an overview of some of the numerical +methods relating to total variation image restoration. Amongst others it +describes some dual and primal-dual methods, as well as the graph cut +approach we take in this thesis. + +\subsubsection{Graph cut approach} + +\fixme{Could probably be shortened a lot, but I thought it could be +useful to at least give an idea how it is done.} + +Using graph cuts is the approach we will be taking later when +considering the anisotropic total variation regularization, and it is +therefore valuable to briefly look into how graph cuts are used in the +case of regular total variation. + +In the theory of graphs, a graph cut is a set of edges that when +removed will separate the graph into two disconnected parts. A minimum +cut is a cut such that the sum of the weight of the edges cut is minimal. It +has been shown that for some energy functionals, it is possible to +construct graphs for which the minimum cuts correspond to minimizers of +the functional. + +In the discrete setting our image consists of pixels, and is represented +by a function $u : \mathcal{G} \to \mathcal{P}$ where $\mathcal{G}$ is a +regular grid over $\Omega$, and $\mathcal{P} = \{0, \hdots, L-1\}$ is +the discrete set of pixel values, or \emph{levels}. We denote the value +in pixel $x$ as $u(x) = u_x$. + +For an image $u$ and a level $\lambda$ we denote the \emph{level set} by +$\{ u > \lambda\}$, defined as the set $\{ x \in \Omega : u_x > \lambda +\}$. The thresholded image $u^\lambda$, an indicator function, is then +defined as +\begin{equation} + u^\lambda = \idfun_{u > \lambda}. +\end{equation} +Here, $\idfun_E$ signifies the characteristic function of the set $E$, +\nomenclature{$\idfun_E$}{Characteristic function of the set $E$.}% +the function which is equal to one in every point in $E$, and zero +elsewhere. + +The idea of the graph cut approach is to decompose the minimization +problem into one minimization problem for each level of the image, and +then solve them separately before combining the results. + +Through careful manipulation of the continuous functional in +\eqref{eq:second_min_presentation} it is possible to obtain a discrete +functional decomposed as a sum over all the level values on the form +\begin{equation} + F(u) = + \sum_{\lambda=0}^{L-2} \sum_x F^x_\lambda(u^\lambda_x) + + \beta \sum_{\lambda = 0}^{L-2} \sum_{(x, y)} F^{x,y}(u^\lambda_x, + u^\lambda_y) + =: \sum_{\lambda=0}^{L-2} F_\lambda(u^\lambda) + \label{eq:old_total_energy} +\end{equation} +where the sum over $(x, y)$ is over all pixel pairs $(x, y)$ in a neighbor +relation, i.e.\ where the pixels are ``close'' to each other. +The actual form of the functional, and the steps to construct it will be +presented later. + +The graph cut we find will for each level $\lambda$ give us the +thresholded image $u^\lambda$, and they can then be combined to form the +complete image $u$. + +When constructing the graph used to find the thresholded image +$u^\lambda$, we have two special vertices, one representing the set $\{ +u > \lambda\}$, while the other represents the set $\{ u \leq +\lambda\}$. The pixels are then connected to these vertices with a +weight representing how strongly they are related to the corresponding +set. This weight will be based on the value of $F^x_\lambda$. + +Additionally there are connections between pixels in a neighborhood relation, +representing the energy $F^{x,y}$. Thus when finding a cut, we partition +the pixels into the sets $\{ u > \lambda \}$ and $\{ u \leq \lambda \}$. +And if in addition the cut is minimal, we know that the edges cut have +minimal weight, and can prove that the $u^\lambda$ found minimizes +the functional in \eqref{eq:old_total_energy}. + diff --git a/results.tex b/results.tex index 0e2635e..14d45d9 100644 --- a/results.tex +++ b/results.tex @@ -134,7 +134,7 @@ what happens and how the parameters affect our end result. In Figure~\ref{fig:lena_process} we show the different stages of our image restoration algorithm. We start out with a section of the much-used Lena test image in Figure~\ref{fig:lena_process_orig}. The first step in the -construction of the structure tensor is to blur the image with parameter +construction of the structure tensor is to apply a Gaussian blur with parameter $\sigma$, and the result is shown in Figure~\ref{fig:lena_process_blur}. The blurring is done so that the edge detector $\nabla f_\sigma$, visualized as $\abs{\nabla f_\sigma}^2$ in Figure~\ref{fig:lena_process_edge}, is not too @@ -147,11 +147,11 @@ eigenvectors of the tensor have been drawn. The length of the vectors have been scaled by the corresponding eigenvalue. Another way of visualizing the tensor is shown in -Figure~\ref{fig:lena_process_color}. The color in each point is decided -by the size of the smallest eigenvalue, and the direction of its -corresponding eigenvector, using the color wheel in +Figure~\ref{fig:lena_process_color}. The brightness and color in each point is decided +by the size of the smallest eigenvalue and the direction of its +corresponding eigenvector using the color wheel in Figure~\ref{fig:color_wheel}. The direction decides the color, -while the brightness, here the radius, is set to $1 / \sigma_1 - 1$, +while the brightness, which is the radius in the color wheel, is set to $1 / \sigma_1 - 1$, where $\sigma_1$ is the smallest eigenvalue. Thus the stronger the anisotropy, the brighter the color, while we expect uniform areas in the original image to be black in the @@ -203,7 +203,8 @@ anisotropy affected the regularization is not obvious however. \end{overpic} %\includegraphics[width=\textwidth]{fig/factory/lines/r_p2_n8_b10000_g100_r4_s2_color.png} \caption{% - Restored with relatively high noise scale of $\sigma = 2$. + High noise scale of $\sigma = 2$ means thinnest lines are + regularized. } \label{fig:scale_comp_high} \end{subfigure} @@ -215,8 +216,8 @@ anisotropy affected the regularization is not obvious however. \end{overpic} %\includegraphics[width=\textwidth]{fig/factory/lines/r_p2_n8_b10000_g100_r15_s0.2_color.png} \caption{% - Restored with low noise scale $\sigma = 0.2$, but large integration - scale $\rho = 15$. + Low noise scale $\sigma = 0.2$, but large integration + scale $\rho = 15$ means inner structure is ignored. } \label{fig:scale_comp_int} \end{subfigure} @@ -225,12 +226,10 @@ anisotropy affected the regularization is not obvious however. An example constructed to show the effects of the parameters $\sigma$ and $\rho$ in the anisotropy tensor. The noise scale $\sigma$ controls the smoothing done before edge - detection, and $\sigma=2$ means the thinnest lines are - not considered to be edges by the tensor, and they are thus - regularized. + detection. The integration scale $\rho$ controls the size of the structures - considered by the tensor, and does in this case allow us to - ignore the smaller structures. Parameters: $\abs{\mathcal{N}} = + considered by the tensor. + Parameters: $\abs{\mathcal{N}} = 8$, $\beta = 10000$, $\omega = 100$. } \label{fig:scale_comp} @@ -570,7 +569,7 @@ approximately equal to the noise removed when using regular total variation in Figure~\ref{fig:finger_contrast_tv}. By visual inspection, the last image seems to have somewhat higher contrast. This can be confirmed by calculating the standard deviation of the images, one of -many possible contrast measures, which gives 22.9, 37.9 and 26.6 +many possible contrast measures, which gives 22.9, 37.9 and 27.2 respectively. To further compare the three restored images, a single row has been @@ -580,7 +579,8 @@ Figure~\ref{fig:finger_contrast_atv2} look promising compared to Figure~\ref{fig:finger_contrast_tv}, some details in the singularity in the upper left is actually lost. The metric tensor is approximately the identity matrix there, and thus an increased restoration parameter -$\beta$ leads to increased smoothing. +$\beta$ leads to increased smoothing compared to the regular total +variation of Figure~\ref{fig:finger_contrast_tv}. \begin{figure} \centering @@ -597,9 +597,9 @@ $\beta$ leads to increased smoothing. So far we have just inspected the restored image visually to judge the results of the method. Depending on the application, we might have different wishes for the results. In many cases we just want to remove the -noise $\delta$ in \fixme{ref} where the captured image $f$ is assumed to -consist of an actual image $u$ and an additional noise term $\delta$ -such that $f = u + \delta$. The noise removed by the restoration method, +noise $\delta$ in $f = u^* + \delta$ where the captured image $f$ is assumed to +consist of an actual image $u^*$ and an additional noise term $\delta$. +The noise removed by the restoration method, $\hat{\delta} = f - u^*$ is called the method noise. Any assumptions of properties we have on $\delta$, we optimally want $\hat{\delta}$ to fulfill as well. Thus if we assumed independently distributed Gaussian @@ -607,6 +607,8 @@ noise, we want to remove that kind of noise. An indication of problems can in that case be if the method noise contains a lot of the details of the image. +\fixme{fuck} + %Figure~\ref{fig:lena_method_noise} shows the results after restoring the %noisy Lena picture in Figure \fixme{noisylena}, with different %parameters, and their method noise. We note that the method noise of the