From aebac67a8627fdffb9cb7b398cf99cf8aa2d7398 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Bj=C3=B8rn=20Rustad?= Date: Fri, 12 Dec 2014 17:48:26 +0100 Subject: [PATCH] Fix and fix --- bib.bib | 17 ++- introduction.tex | 303 ++++++++++++++++++++++++------------- results.tex | 377 +++++++++++++++++++++++++++++------------------ theory.tex | 79 +++++----- 4 files changed, 489 insertions(+), 287 deletions(-) diff --git a/bib.bib b/bib.bib index 879eba1..037da76 100644 --- a/bib.bib +++ b/bib.bib @@ -1,5 +1,19 @@ +@article{esedogluaniso, + author = {Esedoḡlu, Selim and Osher, Stanley J.}, + title = {Decomposition of images by the anisotropic {R}udin-{O}sher-{F}atemi model}, + journal = {Communications on Pure and Applied Mathematics}, + volume = {57}, + number = {12}, + publisher = {Wiley Subscription Services, Inc., A Wiley Company}, + issn = {1097-0312}, + url = {http://dx.doi.org/10.1002/cpa.20045}, + doi = {10.1002/cpa.20045}, + pages = {1609--1626}, + year = {2004}, +} + @book{bartle1995elements, title={The {E}lements of {I}ntegration and {L}ebesgue {M}easure}, author={Bartle, Robert G}, @@ -324,7 +338,7 @@ URL = {http://dx.doi.org/10.1007/s00245-010-9105-x}, } -@INPROCEEDINGS{olsson2009extending, +@inproceedings{olsson2009extending, author={Olsson, C. and Byrod, M. and Overgaard, N.C. and Kahl, F.}, booktitle={Computer Vision, 2009 IEEE 12th International Conference on}, title={Extending continuous cuts: Anisotropic metrics and expansion moves}, @@ -335,6 +349,7 @@ ISSN={1550-5499}, } + @book{megginson, AUTHOR = {Megginson, Robert E.}, TITLE = {An introduction to {B}anach space theory}, diff --git a/introduction.tex b/introduction.tex index aa0d2b1..28324a4 100644 --- a/introduction.tex +++ b/introduction.tex @@ -1,12 +1,88 @@ \chapter{Introduction} -So here maybe we introduce the field of image restoration. Ideas can be -found in the introduction of the project. Maybe we also say something -about the different noise models, and common assumptions concerning the -noise. - -\fixme{introduce u and f as images} -\fixme{introduce Omega} +In everyday life cameras are used to capture a moment and save it for +eternity, but imaging technology technology can be used in many other +contexts, including medical and astronomical applications. +With the advent of computers, tasks previously reserved for the human +brain, like recognizing textures, detecting edges and inferring shapes +and motion, can now be performed algorithmically. The background of +these methods span several fields, including psychology and biology for +the study of human vision, statistics and analysis for the mathematical +foundations, and computer science for the implementation and performance +analysis. + +It is possible to roughly spread the tasks of image processing out along +a line spanning from the raw capturing of light coming into the camera, +a purely physical problem, to the computer vision +methods in the other end, semantically interpreting the scene, +recognizing objects, their position and their movement. In between we +have algorithms working on the captured image, without implying too much +about what the image contains. These methods are often categorized as +\emph{early vision} and includes but is not limited to image +restoration, segmentation, filtering and edge detection. + +In the process of capturing the image with our physical apparatus, +there is always some noise included. Some might come from the physical +nature of how light travels from the scene to the objective, while some +might result from inaccuracies in the construction of the capturing +apparatus. These noise-inducing processes can be studied, and modeled +mathematically, and the process of image restoration looks at how one +can remove some -- or optimally all -- of the noise in the captured +image. These methods often take into account how we would expect a +``normal'' image to look in the capturing conditions, and also what +kinds of noise we expect to be a part of the image. + +In the next section we will see a small overview of some of the most +popular methods in image restoration, but in the rest of the report we +will focus on a total variation based method. One of the strong points +of these methods are their ability to preserve edges, instead of +smoothing over them. Different approaches to total variation image +restoration exist, but we will formulate it as an energy minimization +problem. + +As our input images are digital images, the first step will be to +discretize the energy function. Next comes an important part of this +project which is to reformulate the minimization problem as a series of +minimum cut problems from graph theory. This means we have to carefully +construct graphs and verify that finding the minimum cuts actually +yield a global minimizer of the original energy function. + +A few different algorithms for finding these minumum cuts are +considered, but we will mainly focus on the push-relabel algorithm and +exploit its ability to re-use some of the results across the separate +subproblems. + +Towards the end we will look at some image restoration results showing +how the method performs for different kinds of noise, and for different +input parameters. + +\chapter{Methods in image restoration} + +There are numerous methods in image restoration, and since we do not +have time or space to discuss them all, we will focus on the ones +related to the anisotropic total variation method considered in this +thesis. + +In this Chapter, and also in the rest of the thesis we will assume that +we are given an image $f$ on a rectangular, open domain $\Omega$. The +assumptions on which space this $f$ resides in will vary, but since we +are talking about image restoration, we assume that it includes some +kind of noise. + +There are different kinds of noise models for different situations and +applications, but we will assume that the given image $f$ is a +combination of an actual image $u^*$, and some noise $\delta$. There is +additive noise where the assumption is that $f = u^* + \delta$, and +there is also multiplicative noise where $f = u^* \cdot \delta$. + +In the real world the noise might be more complex of course, and even +come from a combination of sources. Depending on the application, the +goal might even not be to find $u^*$, but rather an output which +fulfills certain smoothness or regularity properties. In any case, we +will denote the output of the methods $u$, and loosely discuss their +properties. + +\fixme{superbad} \section{Diffusion filtering} @@ -28,8 +104,8 @@ 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 $v$ is therefore equivalent to solving the heat equation -with $v$ as initial value, until some time $T > 0$ depending on +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. @@ -44,7 +120,7 @@ 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} +\subsection{Non-linear diffusion filtering} In the theory of the heat equation one can introduce a \emph{thermal diffusivity} $\alpha$ such that @@ -73,20 +149,12 @@ equation \end{equation} The thermal diffusivity $\alpha(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. Although noise will also -contribute to an increase in the size of $\abs{\nabla u}$, the general -assumption is that this is a good measure of where the details are in -the image. 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}. +norm of the gradient $\abs{\nabla u}$ grows. -Note that we follow Weickert's terminology when it comes to the -distinction between non-linear and anisotropic methods. The -Perona--Malik method, and other methods with varying scalar -diffusivities, are often called anisotropic, but we will call them -non-linear, and spare the anisotropic term for methods where the -diffusivity is a tensor, and thus both location and direction dependent. +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 @@ -96,17 +164,18 @@ flow which can be formulated as 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 variational total variation formulation presented later. One -discrete time-step in the solution of this partial differential equation -corresponds to the Euler--Lagrange equation of the variational +forward Euler time-step in the solution of this partial differential +equation corresponds to the Euler--Lagrange equation of the variational formulation. -Common for all of these are that we control the amount of diffusion done -in each point of the image, based somehow on the image itself. But we -can go further by not only making the diffusivity location dependent, -but also direction dependent. This allows us to have have more -diffusion in directions with little variation. As mentioned, we will -use the term anisotropic diffusion, meaning ``not uniform in all -directions.'' +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, +since the diffusion will be different in different directions. +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} @@ -129,8 +198,8 @@ becomes \end{equation} \nomenclature{$A(u)$}{Thermal diffusivity tensor, or anisotropy tensor.}% -where $\nu$ is the outer normal of $\Omega$. The tensor $A(u)$ would -normally be constructed such as to diminish the effect of $\nabla +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 suggests in \cite{weickert1998anisotropic} constructing $A(u)$ based on the @@ -142,27 +211,28 @@ edge estimator $\nabla u_\sigma$ where \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 for example symmetrically extending $u$ across the edges. +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, we introduce the \emph{structure tensor} +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.} +\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 +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 a extremum principle as done in -\cite{weickert1998anisotropic}. +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 make our method more conserving of edges, the solution of \eqref{eq:aniso_diff} will still be @@ -170,10 +240,9 @@ 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, just like its isotropic version, 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. +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 @@ -238,11 +307,11 @@ variation using the distributional derivative. \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$. This means we assume -no-flux conditions on 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. +$\varphi$ vanish on the boundary of $\Omega$. Thus no variation is +measured across 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. @@ -266,17 +335,50 @@ Our optimization problem has thus become and there are many different numerical methods to find or approximate a solution. -As with anisotropic diffusion, the total variation method also has its +As with any restoration method, the total variation method also 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 -over sharp discontinuities like the diffusion methods. This comes at a -cost though, and in some cases the method may introduce new edges that -were not in the original image, since flat sections of zero variation -are encouraged by the functional, an effect called the -\emph{stair-casing effect}. Thin objects and corners may also suffer -from contrast loss since bringing them closer to their surroundings in -level value reduces the total variation. +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{% + 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 in level value 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 @@ -291,7 +393,7 @@ level value reduces the total variation. \includegraphics[width=\textwidth]{fig/ftv.png} \caption{TV} \end{subfigure} - \caption{ + \caption{% Example of how the contrast of thin stuff might be reduced when regularization is too high. \fixme{FIX} } @@ -308,54 +410,43 @@ be found in \cite{caselles2011total}. \subsection{Graph cut approach} +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 set of edges such that the sum of their weight is minimal. It +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. -Using graph cuts is also the approach we will be taking later when -considering the anisotropic total variation minimization, and it is -therefore valuable to briefly look into how graph cuts are used in the -case of regular total variation. +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 \emph{level} of the -image, before solving them separately and combining the results. +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 +$\{ 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} -\nomenclature{$u^\lambda$}{Thresholded image at level $\lambda$.}% 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 can then be combined to +give us the thresholded image $u^\lambda$ which then can be combined to form the complete image $u$. -Using the \emph{coarea formula} which we will spend more space on later, -the total variation can then be rewritten into an integral over the -perimeter of the level sets -\begin{equation} - \TV(u) - = \int_\Omega \abs{D u} \, dx - = \int_\mathbb{R} \int_\Omega \abs{D u^\lambda} \, dx \, d\lambda - = \int_\mathbb{R} \Per(\{u > \lambda\},\Omega) \, d\lambda. - \label{eq:coarea_formula} -\end{equation} - -Our input images will be digital images given as pixels on a grid, and -we will represent them as functions $u : \mathcal{G} \to \mathcal{P}$, -where $\mathcal{G}$ is our regular discrete pixel grid, and -$\mathcal{P} = \{0, \hdots, L-1\}$ is the set of discrete level values. -Through careful manipulation of the functional in -\eqref{eq:second_min_presentation} we obtain a discrete energy -functional decomposed as a sum over all the level values +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) @@ -364,23 +455,27 @@ functional decomposed as a sum over all the level values =: \sum_{\lambda=0}^{L-2} F_\lambda(u^\lambda) \label{eq:total_energy} \end{equation} -where the sum over $(x, y)$ is over all pixel pairs $x$ and $y$ such -that are in a neighbor relation, i.e.\ are ``close'' to each other. +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 it can be useful to consider the intuition behind -the construction, and how it relates to a graph cut. - -The energy $F^{x,y}$ represents how strongly the pixels at $x$ and $y$ -are connected, and will depend on their spatial relation. Pixels close -to each other are more strongly related, and thus more likely to have -similar values. The energy $F^x_\lambda$ represents how strongly we -believe that the pixel at $x$ is connected to the level $\lambda$, and -it is related to the level value of the input image $f$ in that pixel. - -By finding the minimum cut, we could say that we find the minimum set of -relations that has to be broken in order to create a thresholded image -$u^\lambda$. - -\fixme{wow. such handwavy.} +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}. + +\fixme{Wow. Such handwavy.} diff --git a/results.tex b/results.tex index 5d32c5e..fde9ef4 100644 --- a/results.tex +++ b/results.tex @@ -6,14 +6,12 @@ the image restoration results. The anisotropic total variation is the main advance from my project work \cite{project}, and we would like to see what its effects are. -\section{Anisotropy} - \begin{figure} \centering \begin{subfigure}[b]{0.30\textwidth} \centering \includegraphics[width=\textwidth]{fig/circle100n.png} - \caption{ + \caption{% Noisy circle. } \end{subfigure} @@ -21,7 +19,7 @@ see what its effects are. \begin{subfigure}[b]{0.30\textwidth} \centering \includegraphics[width=\textwidth]{fig/anisotropy_circle.png} - \caption{ + \caption{% Uniform tensor. } \label{fig:tensor_experiment_unif} @@ -30,11 +28,11 @@ see what its effects are. \begin{subfigure}[b]{0.30\textwidth} \centering \includegraphics[width=\textwidth]{fig/circle100a.png} - \caption{ + \caption{% Normal tensor. } \end{subfigure} - \caption{ + \caption{% Funny caption. Parameters. } \label{fig:tensor_experiment} @@ -50,70 +48,181 @@ $\epsilon \ll 1$ is small. This would result in us down-weighting the size of $\nabla \tilde{f}_\sigma$ in the $y$-direction, and thus regularizing mostly in the $x$-direction. -The results of this experiment can be seen in Figure -\ref{fig:tensor_experiment}, where a noisy picture of a circle has been -restored in two different ways, first with a uniform anisotropy tensor -$A = \operatorname{diag}(1, \epsilon)$, and then with the tensor -described in Section \ref{sec:anisotropy_tensor}. We see in Figure -\ref{fig:tensor_experiment_unif} that the uniform tensor gives a strong -smoothing in the $x$-direction, while no apparent smoothing in the -$y$-direction. +The results of this experiment can be seen in +Figure~\ref{fig:tensor_experiment}, where a noisy picture of a circle +has been restored in two different ways, first with a uniform anisotropy +tensor $A = \operatorname{diag}(1, \epsilon)$, and then with the tensor +described in Section~\ref{sec:anisotropy_tensor}. We see in +Figure~\ref{fig:tensor_experiment_unif} that the uniform tensor gives a +strong smoothing in the $x$-direction, while no apparent smoothing in +the $y$-direction. + +\begin{figure} + \centering + \begin{subfigure}[t]{0.30\textwidth} + \centering + \includegraphics[width=\textwidth]{fig/factory/lena_eye/lena_eye.png} + \caption{% + Lena's eye. + } + \label{fig:lena_process_orig} + \end{subfigure} + ~ + \begin{subfigure}[t]{0.30\textwidth} + \centering + \includegraphics[width=\textwidth]{fig/factory/lena_eye/r_p2_n16_b8000_g60_r5_s3_blur.png} + \caption{% + Smoothed. + } + \label{fig:lena_process_blur} + \end{subfigure} + ~ + \begin{subfigure}[t]{0.30\textwidth} + \centering + \includegraphics[width=\textwidth]{fig/factory/lena_eye/r_p2_n16_b8000_g60_r5_s3_edge.png} + \caption{% + Edge detector. + } + \label{fig:lena_process_edge} + \end{subfigure} + + \begin{subfigure}[t]{0.30\textwidth} + \centering + \includegraphics[width=\textwidth]{fig/factory/lena_eye/r_p2_n16_b8000_g60_r5_s3_structure.png} + \caption{% + Anisotropy tensor. + } + \label{fig:lena_process_tensor} + \end{subfigure} + ~ + \begin{subfigure}[t]{0.30\textwidth} + \centering + \includegraphics[width=\textwidth]{fig/factory/lena_eye/r_p2_n16_b8000_g60_r5_s3_color.png} + \caption{% + Anisotropy tensor. + } + \label{fig:lena_process_color} + \end{subfigure} + ~ + \begin{subfigure}[t]{0.30\textwidth} + \centering + \includegraphics[width=\textwidth]{fig/factory/lena_eye/r_p2_n16_b8000_g60_r5_s3.png} + \caption{% + Restored Lena. + } + \label{fig:lena_process_restored} + \end{subfigure} + \caption{% + Funny caption. Parameters. + } + \label{fig:lena_process} +\end{figure} + +It is useful to look at intermediate results to better understand +what happens and how the parameters affect our end result. In +Figure~\ref{fig:lena_process} we show the 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 +$\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$, which +is shown squared in Figure~\ref{fig:lena_process_edge}, is not too +sensitive to noise in the image. The structure tensor is then +constructed using $\nabla f_\sigma$, and further transformed into the +anisotropy tensor, which is visualized in +Figure~\ref{fig:lena_process_tensor}. In a selection of points, the +tensor has been drawn as its two eigenvectors, with the corresponding +eigenvalues as their length, and we clearly see that the eigenvalues are +reduced across edges. Another way of visualizing the tensor, which makes +it possible to see it in every point, is using the color wheel in +Figure~\ref{fig:color_wheel}. The color is decided by the tensor angle, +while the brightness, or the radius in the wheel, is set to $1 / +\lambda_2$, the inverse of the smallest eigenvalue. Thus the stronger +the anisotropy, the brighter the color, while for smooth areas we expect +black. + +\begin{figure} + \centering + \includegraphics[width=0.4\textwidth]{fig/wheel.png} + \caption{% + Color wheel. + } + \label{fig:color_wheel} +\end{figure} + +Finally, Figure~\ref{fig:lena_process_restored} shows the restored +image, and we see that it has been heavily regularized. -In the construction of the tensor in Section -\ref{sec:anisotropy_tensor}, we have several parameters. The \emph{noise -scale} $\sigma$ controls the smoothing done before $\nabla -\tilde{f}_\sigma$ is calculated, and thus controls how sensitive we are -to noise in the input image $f$. All details of scale less than $\sigma$ -will be considered to be noise. +In the construction of the tensor in +Section~\ref{sec:anisotropy_tensor}, we have several parameters. The +\emph{noise scale} $\sigma$ controls the smoothing done before $\nabla +f_\sigma$ is calculated, and thus controls how sensitive the edge +detector is to noise in the input image $f$. All details of scale less +than $\sigma$ will be considered to be noise. \begin{figure} \centering \begin{subfigure}[t]{0.30\textwidth} \centering - \includegraphics[width=\textwidth]{fig/gradlow.png} - \caption{ - Low smoothing. + \includegraphics[width=\textwidth]{fig/factory/lines/lines.png} + \end{subfigure} + ~ + \begin{subfigure}[t]{0.30\textwidth} + \centering + \includegraphics[width=\textwidth]{fig/factory/lines/r_p2_n8_b5000_g100_r4_s2.png} + \end{subfigure} + ~ + \begin{subfigure}[t]{0.30\textwidth} + \centering + \includegraphics[width=\textwidth]{fig/factory/lines/r_p2_n8_b5000_g100_r15_s0.2.png} + \end{subfigure} + + \begin{subfigure}[t]{0.30\textwidth} + \centering + \includegraphics[width=\textwidth]{fig/factory/lines/r_p2_n8_b5000_g100_r15_s2_blur.png} + \caption{% + Original and blurred with $\sigma = 2$. } + \label{fig:scale_comp_orig} \end{subfigure} ~ \begin{subfigure}[t]{0.30\textwidth} \centering - \includegraphics[width=\textwidth]{fig/gradnscale.png} - \caption{ - Higher noise scale. + \includegraphics[width=\textwidth]{fig/factory/lines/r_p2_n8_b5000_g100_r4_s2_color.png} + \caption{% + Restored with relatively high noise scale of $\sigma = 2$. } \label{fig:scale_comp_high} \end{subfigure} ~ \begin{subfigure}[t]{0.30\textwidth} \centering - \includegraphics[width=\textwidth]{fig/gradiscale.png} - \caption{ - Low noise scale, but high integration scale. + \includegraphics[width=\textwidth]{fig/factory/lines/r_p2_n8_b5000_g100_r15_s0.2_color.png} + \caption{% + Restored with low noise scale (0.2), but large integration + scale $\rho = 15$. } \label{fig:scale_comp_int} \end{subfigure} - \caption{ + \caption{% Funny caption. Parameters. } \label{fig:scale_comp} \end{figure} -Figure \ref{fig:scale_comp} shows a constructed zebra pattern of -increasing width. The anisotropy should in theory help reduce contrast +Figure~\ref{fig:scale_comp} shows a constructed zebra pattern of +increasing width. The anisotropy introduced should in theory help reduce contrast loss in this situation, by reducing the regularization done in the -$x$-direction, but there is the question of at which scale it considers -the lines to be detail instead of noise. +$x$-direction, but there is the question of how the different scales +affect the regularization. -In Figure \ref{fig:scale_comp_high} the noise scale is increased and the -tensor no longer consider the left-most part of the image to be details, -but rather noise. +In Figure~\ref{fig:scale_comp_high} the noise scale is increased such +that the edge detector, and thus the anisotropy tensor, considers these +fine lines to be noise rather than details. Thus this part is +regularized to a smooth gray area. The \emph{integration scale} $\rho$ of \eqref{eq:s_def} represents the -size of the structures we want to affect our anisotropy tensor. Although -all edges and details larger than the noise scale will affect the edge -detector $\nabla \tilde{f}_\sigma$, but the smoothing after of scale -$\rho$ may ignore them depending on their size. In Figure -\ref{fig:scale_comp_int} a high value for $\rho$ means that the +size of the structures we want to detect with our anisotropy tensor. In +Figure~\ref{fig:scale_comp_int} a high value for $\rho$ means that the structure found in the inner square is almost completely ignored in favour of the larger, more coherent structure around. @@ -121,57 +230,95 @@ favour of the larger, more coherent structure around. \centering \begin{subfigure}[t]{0.30\textwidth} \centering - \includegraphics[width=\textwidth]{fig/lenaedge.png} - \caption{ - Largest structure tensor eigenvalue. + \includegraphics[width=\textwidth]{fig/factory/finger/n_q300.png} + \caption{% + Noisy fingerprint. } - \label{fig:lena_process_edge} \end{subfigure} ~ \begin{subfigure}[t]{0.30\textwidth} \centering - \includegraphics[width=\textwidth]{fig/lenatensor.png} - \caption{ - The anisotropy tensor eigenvector and eigenvalues. + \includegraphics[width=\textwidth]{fig/factory/finger/r_p2_n32_b9000_g80_r10_s3_color.png} + \caption{% + $\rho = 10$ } - \label{fig:lena_process_tensor} \end{subfigure} ~ \begin{subfigure}[t]{0.30\textwidth} \centering - \includegraphics[width=\textwidth]{fig/lenarest.png} - \caption{ - Restored Lena. + \includegraphics[width=\textwidth]{fig/factory/finger/r_p2_n32_b9000_g80_r20_s3_color.png} + \caption{% + $\rho = 20$ } - \label{fig:lena_process_restored} \end{subfigure} - \caption{ - Funny caption. Parameters. + \caption{% + Not so funny caption. Include some parameters. } - \label{fig:lena_process} + \label{fig:finger_tensor} \end{figure} -It is often useful to look at intermediate results to better understand -what happens and how the parameters affect our end result. In Figure -\ref{fig:lena_process} we show three stages of our image restoration -algorithm. Figure \ref{fig:lena_process_edge} shows the larges -eigenvalue of our structure tensor, and is thus an indicator of where -the edge structures are. +Another visualization of the anisotropy tensor can be seen for a +fingerprint image in Figure~\ref{fig:finger_tensor}. We see that the +gradual changes in the direction of the edges are captured by the +tensor. Even the singularities in the fingerprint can be identified by +finding the dark spots, as the changes in these areas do not tend to one +single direction. -Figure \ref{fig:lena_process_tensor} is a depiction of the anisotropy -tensor in some of the points in the image, shown as eigenvectors with -their corresponding eigenvalues as their length. We see that across -edges, \fixme{what!?} +\fixme{We show how the different neighborhoods work. Here we might make some +nice, very contrived examples, maybe. To show that both too small and +too big is bad.} -Finally, Figure \ref{fig:lena_process_restored} shows the final restored -image, and we see that heavy modifications have been made. +\begin{figure} + \centering + \begin{subfigure}[t]{0.30\textwidth} + \centering + \includegraphics[width=\textwidth]{fig/factory/circle_deform/n_q0.png} + \caption{% + Circle. + } + \end{subfigure} + ~ + \begin{subfigure}[t]{0.30\textwidth} + \centering + \includegraphics[width=\textwidth]{fig/factory/circle_deform/r_p2_n8_b100000_g10000000000_r1_s1.png} + \caption{% + Size of neighborhood: 8. + } + \end{subfigure} + ~ + \begin{subfigure}[t]{0.30\textwidth} + \centering + \includegraphics[width=\textwidth]{fig/factory/circle_deform/r_p2_n32_b100000_g10000000000_r1_s1.png} + \caption{% + Size of neighborhood: 32. + } + \end{subfigure} + \caption{% + Not so funny caption. Include some parameters. + } + \label{fig:circle_arti} +\end{figure} + +The neighborhoods were introduced as a way to describe the discrete set +of lines $\mathcal{L}_D$ in the discretization of the regularization +term in Section~\ref{sec:disc_regularization}. And as discussed, we want +the stencil to have many short edges, such that the angular +differences $\Delta \phi$, the inter-line distances $\Delta \rho$ and the +edge lengths $e$ are ``small.'' + +Figure~\ref{fig:circle_arti} shows a circle that has been strongly +regularized using a 8-point stencil. The circle has clearly been +deformed to look more like an octagon. This is due to the fact that the +8-point stencil measures horizontal, vertical and 45\textdegree{} lines +exactly in the Cauhcy--Crofton formula, while the length of other lines +overestimated as the sum of their three components. \begin{figure} \centering \begin{subfigure}[t]{0.30\textwidth} \centering \includegraphics[width=\textwidth]{fig/fingern.png} - \caption{ + \caption{% Noisy fingerprint. } \end{subfigure} @@ -179,7 +326,7 @@ image, and we see that heavy modifications have been made. \begin{subfigure}[t]{0.30\textwidth} \centering \includegraphics[width=\textwidth]{fig/fingertv.png} - \caption{ + \caption{% Restored without much anisotropy. } \label{fig:finger_contrast_tv} @@ -188,12 +335,12 @@ image, and we see that heavy modifications have been made. \begin{subfigure}[t]{0.30\textwidth} \centering \includegraphics[width=\textwidth]{fig/fingera.png} - \caption{ + \caption{% Restored with much anisotropy. } \label{fig:finger_contrast_atv} \end{subfigure} - \caption{ + \caption{% Funny caption. Parameters. } \label{fig:finger_contrast} @@ -201,20 +348,20 @@ image, and we see that heavy modifications have been made. We have claimed that the anisotropic total variation will reduce some of the contrast loss one can encounter with regular total variation -regularization. An example of this can be seen in Figure -\ref{fig:finger_contrast} where a portion of a fingerprint with +regularization. An example of this can be seen in +Figure~\ref{fig:finger_contrast} where a portion of a fingerprint with added Gaussian noise has been restored in two different ways. We see that when the anisotropy is increased, more of the contrast between dark and light portions of the fingerprint is retained. This can also be seen -in Figure \ref{fig:contrast_plot} where a one-dimensional slice has been -taken through Figure \ref{fig:finger_contrast_tv} and -\ref{fig:finger_contrast_atv}. Although much of the same structure is -found for this particular regularization parameter $\beta$, the peaks +in Figure~\ref{fig:contrast_plot} where a one-dimensional slice has been +taken through Figure~\ref{fig:finger_contrast_tv} +and~\ref{fig:finger_contrast_atv}. Although much of the same structure +is found for this particular regularization parameter $\beta$, the peaks are much higher for the anisotropic total variation. \begin{figure} \centering - \input{fig/factory/finger/contrast_plot} + \input{fig/factory/contrast} \caption{Showing some contrast loss.} \label{fig:contrast_plot} \end{figure} @@ -232,7 +379,7 @@ 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. -Figure \ref{fig:lena_method_noise} shows the results after restoring the +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 regular total variation method shows and outline of Lena, which @@ -260,7 +407,7 @@ anisotropically restored image, has less details. \begin{subfigure}[t]{0.30\textwidth} \centering \includegraphics[width=\textwidth]{fig/lena_method_tv.png} - \caption{ + \caption{% Regular total variation. } \end{subfigure} @@ -268,7 +415,7 @@ anisotropically restored image, has less details. \begin{subfigure}[t]{0.30\textwidth} \centering \includegraphics[width=\textwidth]{fig/lena_method_atv.png} - \caption{ + \caption{% Anisotropic total variation. } \end{subfigure} @@ -276,73 +423,17 @@ anisotropically restored image, has less details. \begin{subfigure}[t]{0.30\textwidth} \centering \includegraphics[width=\textwidth]{fig/lena_method_atv72.png} - \caption{ + \caption{% Size 72 neighborhood. } \end{subfigure} - \caption{ + \caption{% Funny caption. Parameters. } \label{fig:lena_method_noise} \end{figure} -\begin{figure} - \centering - \input{fig/factory/lena_gamma_seq/corr.tex} - \caption{ - The correlation between the method noise and the noisy image for - different anisotropy parameters $\gamma$. - } -\end{figure} +\section{Conclusion} -\section{Neighborhood} - -We show how the different neighborhoods work. Here we might make some -nice, very contrived examples, maybe. To show that both too small and -too big is bad. - -The neighborhoods were introduced as a way to describe the discrete set -of lines $\mathcal{L}_D$ in the discretization of the regularization -term in Section \ref{sec:disc_regularization}. And as discussed, we want -the stencil to have many short edges, such that the angular -differences $\Delta \phi$, the inter-line distances $\Delta \rho$ and the -edge lengths $e$ are ``small.'' - -Figure \ref{fig:circle_arti} shows a circle that has been strongly -regularized using a 8-point stencil. The circle has clearly been -deformed to look more like an octagon. This is due to the fact that the -8-point stencil measures horizontal, vertical and 45\textdegree{} lines -exactly in the Cauhcy--Crofton formula, while the length of other lines -overestimated as the sum of their three components. - -\begin{figure} - \centering - \begin{subfigure}[t]{0.30\textwidth} - \centering - \includegraphics[width=\textwidth]{fig/factory/circle_deform/n_q0.png} - \caption{ - Circle. - } - \end{subfigure} - ~ - \begin{subfigure}[t]{0.30\textwidth} - \centering - \includegraphics[width=\textwidth]{fig/factory/circle_deform/r_p2_n8_b100000_g10000000000_r1_s1.png} - \caption{ - Size of neighborhood: 8. - } - \end{subfigure} - ~ - \begin{subfigure}[t]{0.30\textwidth} - \centering - \includegraphics[width=\textwidth]{fig/factory/circle_deform/r_p2_n32_b100000_g10000000000_r1_s1.png} - \caption{ - Size of neighborhood: 32. - } - \end{subfigure} - \caption{ - Not so funny caption. Include some parameters. - } - \label{fig:circle_arti} -\end{figure} +Bleep bloop. Bad sides, good sides? diff --git a/theory.tex b/theory.tex index cc3cf49..cb02987 100644 --- a/theory.tex +++ b/theory.tex @@ -5,10 +5,17 @@ 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 anisotropic diffusion and total variation methods, while eliminating -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, briefly discuss its well-posedness, and through some -important theorems, transform it into a form that is easy to discretize. +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, its +different forms and briefly discuss its well-posedness. Through the +coarea formula, the anisotropic total variation is rewritten as an +integral of the perimeter of all the level sets. + +Later, the 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 following chapter. \section{Anisotropic total variation} @@ -48,11 +55,11 @@ as $A(x) = \Ahalf(x) \Ahalf(x)$. We can then write &= \sup_{\eta^T A^{-1} \eta \leq 1} \int_\Omega u \diver \eta \, dx, \end{align} - where $\xi$ and $\eta = \Ahalf \xi$ are in $C_c^\infty(\Omega, -\mathbb{R}^2)$, the space of $C^1$ vector fields with compact support. -If we let $\norm{\xi}_A = \sqrt{\xi^T A \xi}$ and $\norm{\eta}_A^* = -\sqrt{\eta^T A^{-1} \eta}$. +\mathbb{R}^2)$, the space of $C^\infty$ vector fields with compact support. +In the following we let $\norm{\xi}_A = \sqrt{\xi^T A \xi}$ and +$\norm{\eta}_A^* = \sqrt{\eta^T A^{-1} \eta}$, 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 symmetric positive definite tensor $A : \Omega \to \mathbb{R}^{2 \times 2}$ we define the @@ -71,30 +78,22 @@ anisotropy tensor $A(u)$.}% With this extended definition, we have arrived at a minimization problem where we seek to find a minimizer of the functional \begin{equation} - F(u) = \int_\Omega (u - f)^2 \, dx + \beta \, \TVA(u) + F(u) = \int_\Omega (u - f)^2 \, dx + \beta \, \TVA(u). \label{eq:first_anisotropic_functional} \end{equation} -The question is then how to construct this anisotropy tensor $A(x)$ to -get the improvements we hope for, and how this affects our numerical +Similar functionals have been considered in +\cite{grasmair2010anisotropic} and \cite{olsson2009extending}. The +question is now how to construct this anisotropy tensor $A(x)$ to get +the improvements we hope for, and how this affects our numerical solution methods. -\fixme{Describe how this is similar to a Euler time step in a related -equation \cite{grasmair2010anisotropic}. This is only without the square -root.} - \subsection{Anisotropy tensor} \label{sec:anisotropy_tensor} \fixme{% - This tensor can be visualized, and can be used for edge - detection. There is probably a bit to be said about the smoothing we - do. We also discuss whether to use the noisy image (yes, probably), + We also discuss whether to use the noisy image (yes, probably), or the smoothed image (implicitly, complicated, possible - iteratively) in the structure tensor. We then describe how the - structure tensor is decomposed and transformed into our anisotropy - tensor. Then there are some problems related to stability which we - can discuss here, or maybe we should leave it for a later, more - implementation-focused chapter. + iteratively) in the structure tensor. } There are many possible choices for the anisotropy tensor $A(x)$. Our @@ -111,35 +110,37 @@ none, so we somehow want to be sure about the edges we find. Edges can be found in many different ways, but as suggested by Weickert in his book on Anisotropic Diffusion \cite{weickert1998anisotropic}, and briefly mentioned in Section~\ref{sec:anisotropic_diffusion}, a good -starting point is the \emph{edge detector} $\nabla \tilde{f}_\sigma$ and +starting point is the \emph{edge detector} $\nabla f_\sigma$, +where $f_\sigma = K_\sigma * \tilde{f}$ and $\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$. + +As seen in Figure~\ref{fig:edges}, the edge detector is fine for +detecting edges, but it can not give us information about larger +structures, like corners and textures, which is why we introduce the structure tensor \begin{equation} - S_\rho(x) := K_\rho * \left( \nabla \tilde{f}_\sigma(x) \otimes - \nabla \tilde{f}_\sigma(x) \right), + S_\rho(x) := K_\rho * \left( \nabla f_\sigma(x) \otimes + \nabla f_\sigma(x) \right). \label{eq:s_def} \end{equation} -where $\tilde{f}$ is the symmetric extension of the initial image $f$ in -$\mathbb{R}^2$. The $\sigma$ signifies a Gaussian smoothing -$\tilde{f}_\sigma = K_\sigma * \tilde{f}$. This initial smoothing -parameter $\sigma$, called the \emph{noise scale} by Weickert, will -control how sensitive the edge detector $\nabla \tilde{f}_\sigma$ is to -noise and small details. \begin{figure} \centering{} \includegraphics[width=0.4\textwidth]{fig/edges.png} \caption{% The largest eigenvalue of the structure tensor is $\abs{\nabla - \tilde{f}_\sigma(x)}^2$ and is as we see a good edge detector. + f_\sigma(x)}^2$ and can as we see function as an edge detector. } \label{fig:edges} \end{figure} -First consider the tensor $S_0(x) = \nabla \tilde{f}_\sigma(x) \otimes -\nabla \tilde{f}_\sigma(x)$. It obviously contains the same information +First consider the tensor $S_0(x) = \nabla f_\sigma(x) \otimes +\nabla f_\sigma(x)$. It obviously contains the same information as the edge detector itself. Its eigenvalues will be 0 and $\labs{\nabla -\tilde{f}_\sigma(x)}^2$ with corresponding eigenvectors $v_1$ and $v_2$ -perpendicular and parallel to $\nabla \tilde{f}_\sigma(x)$ respectively. +f_\sigma(x)}^2$ with corresponding eigenvectors $v_1$ and $v_2$ +perpendicular and parallel to $\nabla f_\sigma(x)$ respectively. Figure~\ref{fig:edges} shows that the largest eigenvalue of the structure tensor is a good edge detector. @@ -1016,8 +1017,8 @@ have constructed in \eqref{eq:per_to_length2}. The structure tensor is constructed as described in Section~\ref{sec:anisotropy_tensor} \begin{equation} - S_\rho(x) = \left(K_{\rho} * \left( \nabla \tilde{f}_{\sigma} - \otimes \nabla \tilde{f}_{\sigma} \right)\right)(x), + S_\rho(x) = \left(K_{\rho} * \left( \nabla f_{\sigma} + \otimes \nabla f_{\sigma} \right)\right)(x), \end{equation} Because of the convolutions with the Gaussian function, this is a smooth continuous map from $\bar{\Omega}$ to $\mathbb{R}^{2 \times 2}$. As we -- 2.47.3