In my master thesis I designed and implemented a multi-threaded C++ template library for total variation minimization of manifold-valued data. The library has the rather unwieldy name Manifold Total Variation Minimization Template Library (MTVMTL).
Hence, should anybody come up with something catchier which captures the main features I would love to know it. What I found particularly interesting about the topic is the combination of differential geometry, numerical analysis and software engineering. Furthermore, despite the rather abstract nature of a manifold, the results and effects of the algorithms can be illustrated quite visually, which I hope to achieve in the next couple of paragraphs.
A manifold is a mathematical object which locally looks like Euclidian, that means flat, space. A prime example for a manifold is the sphere, which globally has a non-negligible (constant) curvature, but looks flat locally. Since Euclidian space is in fact flat everywhere it is trivially also a manifold.
One does not need to look too far away to find data that can be modelled using a manifold: A color pixel in the RGB-color-model can be described as a tuple of three values describing the red, green and blue value. As such it can be interpreted as an element of Euclidian . It is, however, also possible to use a different non-linear color model like the Chromaticity-Brightness model. Here the brightness describes the intensity of the pixel and the chromaticity describes its normalized color such that the latter in fact has values on the sphere. More complex examples are , the group of rotations and , the cone of symmetric positive matrices or Grassmannian , which are usually matrix valued and can be used as model for a variety of data occuring in medical imaging, computer vision, machine learning, robotic and more.
Let's assume that a simple grayscale picture, like the one on the left which, is corrupted by noise. The defining characteristic of noise between the two pictures is the fact that formerly uniformly colored areas of the picture are now subject to small variations from one pixel to another. Hence, denoising a picture is equivalent to removing this small variations but at the same time keeping important features of image like edges, which are characterized by large variations when moving from one pixel to the next.
Since a variation in the case of a picture is given by the numerical difference of neighboring pixel values, an indicator of how noisy the picture is, is the sum of the absolute values of all those differences between neighboring pixels, which is called the total variation.
For the case of a of real-valued function , total variation is given by
where denotes the partition of the intervall into subintervals. For differentiable functions, this corresponds to the norm of its gradient. It is important to note that the absolute value makes the above expression for non-differentiable. To avoid that one could also use squared differences , which corresponds to the use of the norm. This, however, will favor smooth transitions in the minimization process, as one can see in the following plot.
We have three functions with steps which all have total variation equal to 1.0. The squared variations however are 1.0 for , 0.1 for and 0.01 for . If we performed a minimization of the sum of quadratic variations, clearly smoother transitions like the blue curve will lead to a smaller sum and will consequently be preferred. Hence, edges, which are characterized by a non-smooth transition, will become blurred.
The final minimization problem we have to deal with is given by the following functional
where is a function describing a 2D grayscale image, while denotes the original noisy image. The first sum is the so-called fidelity term which becomes large if the denoised image becomes very different from the original, noisy picture . Without it minimization just leads to a uniformly colored picture since then the total variation is zero.
Generalization to Manifolds
The above functional needs to be adjusted since the differences inside the absolute values require a vector space structure. Put differently, in non-flat space the shortest distance between two points is not given by a straight line any more. Also obtaining the vector connecting two points on the manifold cannot be obtained by subtracting those two points any more: The difference between north and south pole of the sphere is point which lies in the surrounding space outside the sphere. Without going to much into detail the vector space operations need to be reinterpreted in the following way
- Subtraction: becomes
- Addition: becomes
- Distance: becomes .
Here, and are special mappings between the manifold and its tangent space at a given point.
The (anisotropic) generalized functional then becomes
The problem that remains is that it is still non-differentiable such that standard methods from smooth optimization theory don't work. There are two algorithms implemented in MTVMTL to bypass that problem. The first is based on the proximal point algorithm, the second is an adaption of the Iteratively Reweighted Least Squares (IRLS) algorithm. The description of both is rather technical but to briefly sum it up, IRLS works by alternating between calculating weights of a regularized functional
and minimizing the regularized functional using smooth methods. For the minimization step the previously calculated steps are considered constant. In the case of MTVMTL the minimization step is realized using the Riemannian Newton Method.
Proximal point, on the other hand, is based on projection onto convex sets and simplifies in the differentiable case to a simple gradient descent iteration.
Architecture of the library
The above picture shows the overall structure of the library and its main classes. At the very bottom the Manifold class handles contains the the whole differential geometric structure of the given data along with methods to compute various relevant quantities such as distances or exponential mappings. The Data class deals with storage , I/O as well as different tool like edge detection while the Visualization class is used to render interactive 2D and 3D visualizations like
Finally, the Functional class contains all parameters and methods to define a functional and calculate its values as well as its derivatives and the TV_Minimizer class contains the actual implementation of the minimization algorithms.
Diffusion Tensor MRI