Pierre-Yves Gousenbourger, Applied Mathematics

Pierre-Yves Gousenbourger

Teaching assistant at Université catholique de Louvain

T.A. at Université catholique de Louvain

T.A. at UCLouvain

Academic curriculum in short

I started my PhD at Université catholique de Louvain in September 2014 under the supervision of Pierre-Antoine Absil and Laurent Jacques. I conduct my thesis at the ICTEAM research institute, in Louvain-la-Neuve.

I graduated from the Université catholique de Louvain as an engineer in Applied Mathematics and conducted a master thesis under the supervision of Pierre-Antoine Absil.

My research interests concern interpolation, optimization and imaging on manifolds, mainly conducted by Bézier functions.

Position

I am a teaching assistant at École polytechnique de Louvain (EPL) and very interested in all pedagogical aspects of my job. More precisely, I supervise students for classes in the Electrical Engineering (ELEN) department of the ICTEAM institute.

Miscellaneous

I recently shared my office with Vincent Schellekens (phD), proper substitute to Valerio.
Two former office mates are Dr. Valerio Cambareri, a gentleman now gone in the industry to share his great mind, and Dr. François Rottenberg, another gentleman, creator of WaveComBox and also known as The Brain Of The Lab.

Unidimensional Bézier fitting on manifolds

The natural next step after developping interpolation techniques on manifolds is to relax the interpolation constraint. I'm now interested in fast Bézier techniques to fit a curve $\gamma \in \mathcal{M}$ to manifold-valued data points $d_0,\dots,d_n \in \mathcal{M}$ associated to time-stamps $t_0,\dots,t_n$ where $t_i = i \in \mathbb{Z}$.

An application of this project concerns wind field estimation for unmanned aerial vehicule (UAV) control (drones, if you prefer). In that application the wind field can be characterized by a covariance matrix of rank r, represented in the manifold of SDP matrices of size p and rank r ($\mathcal{S}_+(p,r)$). However, computing such a matrix for a given prevalent wind orientation $\theta_i$ requires time-consumable numerical simulations. Therefore, give a bench of $n$ covariance matrices $C(\theta(i))$, therefore, there is an interrest in finding a fitting a curve $\gamma(\theta(t)): [0,n-1] \to \mathcal{S}_+(p,r)$ with a fast a numerically efficient method.

A first approach was presented at the ESANN2017 conference, as a joint work with Estelle Massart (UCL, Belgium) and Antoni Musolas (MIT, USA). That approach modifies the previous work on interpolation and thus still provides elegant and fast solutions in closed form.

Learn more/less.

In the solution we propose, the curve $\gamma$ takes the form of a composite cubic Bézier curve $$B: [0,n] \to \mathcal{M}: B(t) = \beta_i(t-i;p_i,b_i^{+}, b_{i+1}^{-}, p_{i+1}), \quad \text{with } i = \lfloor t \rfloor,$$ where $\beta_{i}$ is a cubic Bézier function $p_k, b_k^{\pm} \in \mathcal{M}$ are respectively the extreme control points of $\beta_i$ and the inner control points.

In the following paper, we relax the interpolation constraint from Arnould et al (2015) and we search the curve $B(t)$ such that $$ \int_0^n \left\| \frac{D^2 B(t)}{dt^2} \right\|^2_\mathcal{M} \mathrm{d} t + \lambda \sum_{i=0}^n \mathrm{d}_\mathcal{M}(d_i,B(t_i))$$ is minimized on the Euclidean space.

The most recent method to achieve this goal is called the Blended cubic spline method. The method consists in building polynomial pieces by solving the above-mentioned optimization problem in various (Euclidean) tangent spaces, and then blending together corresponding pieces by means of carefully chosen weights in order to satisfy all the following properties:

The curve $B(t)$ is differentiable on its domain;
The curve $B(t)$ interpolates the data points at their associated parameter when $\lambda$ goes to infinity;
The curve $B(t)$ is the natural smoothing spline on the Euclidean space (i.e., the optimization problem is solved);
The algorithm is only based on two operators: the Riemannian exponential and logarithm;
The algorithm builds the curve $B(t)$ at a given $t$ in $\mathcal O(1)$ operations;
The curve $B(t)$ is stored with only $\mathcal O(n)$ vectors, where $n$ is the number of data points.

This method is summarized in the following paper:

Pierre-Yves Gousenbourger, Estelle Massart, P.-A. Absil
Data fitting on manifolds with composite Bézier-like curves and blended cubic splines
Journal on Mathematical Imaging and Vision, Springer.
DOI: 10.1007/s10851-018-0865-2
2018.

Another approach consists in generalizing the optimality conditions to manifolds. The curve $B(t)$ is computed (in closed form) with very few elements of Riemannian geometry (the exp and the log, only !) and returns the natural Bézier spline on the Euclidean space and is differentiable everywhere.

Pierre-Yves Gousenbourger, Estelle Massart, Antoni Musolas, P.-A. Absil, Laurent Jacques, Julien M. Hendrickx, Yousef Marzouk
Piecewise-Bézier $C^1$ smoothing on manifolds with application to wind field estimation.
ESANN2017, Springer.
2017.

The efficiency of our fast solution is analysed in the following papers. We compare the quality of our reconstruction with the solution obtained by directly optimizing the objective function.

In the following paper, we derive a variational model in which we approximate the acceleration by discretizing the squared second order derivative along the curve. We derive a closed-form, numerically stable and efficient algorithm to compute the gradient of a Bézier curve on manifolds with respect to its control points, expressed as a concatenation of so-called adjoint Jacobi fields. This is joint work with Dr. Ronny Bergmann, from TU Chemnitz.

Ronny Bergmann, Pierre-Yves Gousenbourger
A Variational Model for Data Fitting on Manifolds by Minimizing the Acceleration of a Bézier Curve
Frontiers in Applied Mathematics and Statistics, 4(59), 1-16.
DOI: 10.3389/fams.2018.00059
2018.

The paper hereunder is the preliminary study that lead to the document above.

Pierre-Yves Gousenbourger, Laurent Jacques, P.-A. Absil
Fast Method to Fit a C1 Piecewise-Bézier Function to Manifold-Valued Data Points: How Suboptimal is the Curve Obtained on the Sphere S2?
Geometric Science of Information: 3rd International Conference 2017, France.
Proceedings (GSI2017), Springer.
2017.

Bidimensional Bézier interpolation on manifolds

In this project, we receive a set of data points evaluated on a Riemannian manifold $\mathcal{M}$ (for instance the space of matrices of rotation like in the picture; or the sphere). These points $p_{ij}$ are time-labeled, so associated to time-stamps $(t_1^i,t_2^j)$ organised on a regular rectangular grid $D$.

Based on these points, the goal of the game is to propose a smooth (i.e. $C^1$) surface of Bézier $$\mathfrak{B}: D \to \mathcal{M}: (t_1,t_2) \mapsto \mathfrak{B}(t_1,t_2)$$ such that the data points are interpolated by the surface, meaning $\mathfrak{B}(t_1^i, t_2^j) = p_{ij}$.

This task is very well understood when the space is Euclidean, but quite less when points live on a nonlinear manifold. It is also a generalization of my master thesis.

Learn more/less.

This work was done in collaboration with the university of Münster (Applied mathematics department) joinly with Benedikt Wirth and Paul Striewski. It has application in shape interpolation, matrix reconstruction or cartography.

The amont of possible functions interpolating the data points is infinite. We reduce the search set by constraining the surface to be a piecewise cubic Bézier spline, i.e. driven by control points. This family of function is of the form $$\beta(t_1,t_2,\mathbf{b})$$ where $\mathbf{b} \in \mathcal{M}$ is a set of 16 control points on the manifold, and $\beta$ corresponds to a weighted (Karcher) mean of these control points. This choice implies a very simple (geodesic) constrain on the control points to achieve the smoothness of the spline. Finally, to propose the best function possible, a solution is to minimize its mean square (Riemannian) acceleration over the control points.

P.-A. Absil, Pierre-Yves Gousenbourger, Paul Striewski, Benedikt Wirth
Differentiable piecewise-Bézier surfaces on Riemannian manifolds
SIAM Journal on Imaging Sciences, 9(4), 1788-1828.
DOI: 10.1137/16M1057978
2017.

As the control points generation is quite long in the above article, we developped an other method to compute those based on a B-spline representation of the curve. We increase drastically the speed of the generation by paying the price of a smallest accuracy.

P.-A. Absil, Pierre-Yves Gousenbourger, Paul Striewski, Benedikt Wirth
Differentiable piecewise-Bézier interpolation on Riemannian manifolds.
24th European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning (ESANN 2016).
2016.

We developped different algorithms and methods to compute these surfaces on a general Riemannian manifold. Mainly, the computation of control points and the reconstruction of the curve may vary. Check our Numerical Tour.

Unidimensional Bézier interpolation on manifolds

This project was my master thesis lead in 2013 under Pierre-Antoine Absil's supervision. We receive a set of data points $p_0,\dots,p_n \in \mathcal{M}$ associated to time-stamps $t_0,\dots,t_n$ where $t_i = i \in \mathbb{Z}$. The goal of the project is to find a curve $\gamma$ such that the curve is smooth and looks "nice". Furthermore, we constraint the curve to interpolate the data points as $\gamma(t_i) = p_i$. This problem finds place in medical imaging and in volume reconstruction based on MRI.

With Pierre-Antoine, we investigated on manifold-valued Bézier functions and on the usability of the De Casteljau algorithm in such situations. At the end of the day, we obtained that curves can be computed from a simple linear (Euclidean !) system resolution.

Learn more/less.

Research groups

I'm part of the following research group:

PhD thesis

P.-Y. Gousenbourger. Interpolation and fitting on Riemannian manifolds. PhD thesis, UCLouvain., 2020.
The access to constantly increasing computational capacities has revolutionized the way engineering is seen. We are now able to produce a large quantity of data, thanks to cheap sensors. However, processing such data remains costly both in computational time and in energy. One of the reasons is that the structure of the data is often omitted or unknown. The search space becomes so large that finding the solution to simple problems often turns out to be finding a needle in a haystack.
A classical problem in data processing is called the "fitting problem". It consists in fitting a d-dimensional curve to a set of data points associated to d parameters. The curve must pass sufficiently close to the data points while being regular enough. When the underlying structure of the data points (i.e., the manifold) is known, one can impose to the curve to preserve this structure (i.e., to remain on the manifold), such that the search space is drastically reduced.
The goal of this thesis is to develop methods to (approximately) solve this fitting problem; the bet is to require "less" (less computational capabilities, power, storage, time) by leveraging “more” knowledge on the search space. The objective is the following: provide a toolbox that produces a differentiable fitting curve to data points on manifolds, based on very few and simple geometric tools, at low computational cost and storage capacity, all this while maintaining an acceptable quality of the solutions.
The algorithms are applied to different illustrative problems. In 3D shape reconstruction, the data points are organs contours acquired via MRI, and the parameter is the acquisition depth; in wind fields estimation, the data points belong to the manifold of positive semi-definite matrices of given rank, and the parameters are the prevalent wind amplitudes and angles. We also show the performances of our algorithms in applications for parametric model order reduction.
```
@phdthesis{Gousenbourger2020thesis,
  author  = {Gousenbourger, Pierre-Yves},
  title	  = {Interpolation and fitting on Riemannian manifolds},
  school  = {UCLouvain, ICTEAM institute},
  year	  = {2020}
}
```

Journal papers

P.-Y. Gousenbourger, R. Bergmann. A variational model for data fitting on manifolds by minimizing the acceleration of a Bézier curve. Frontiers in Applied Mathematics and Statistics, 4(59), 2018.
We derive a variational model to fit a composite Bézier curve to a set of data points on a Riemannian manifold. The resulting curve is obtained in such a way that its mean squared acceleration is minimal in addition to remaining close the data points. We approximate the acceleration by discretizing the squared second order derivative along the curve. We derive a closed-form, numerically stable and efficient algorithm to compute the gradient of a Bézier curve on manifolds with respect to its control points, expressed as a concatenation of so-called adjoint Jacobi fields. Several examples illustrate the capabilites and validity of this approach both for interpolation and approximation. The examples also illustrate that the approach outperforms previous works tackling this problem.
```
@article{Bergmann2018a,
  title   = {{A variational model for data fitting on manifolds by minimizing 
             the acceleration of a B\'ezier curve}},
  author  = {Bergmann, Ronny and Gousenbourger, Pierre-Yves},
  journal = {Frontiers in Applied Mathematics and Statistics},
  doi     = {10.3389/fams.2018.00059},
  volume  = {4},
  number  = {59},
  pages   = {1--16},
  year    = {2018}
}
```
P.-Y. Gousenbourger, E. Massart, P.-A. Absil. Data Fitting on Manifolds with Composite Bézier-Like Curves and Blended Cubic Splines. Journal of Mathematical Imaging and Vision, 61(5), pp. 645-671, 2018.
We propose several methods that address the problem of fitting a $C^1$ curve $\gamma$ to time-labeled data points on a manifold. The methods have a parameter, $\lambda$, to adjust the relative importance of the two goals that the curve should meet: being "straight enough" while fitting the data "closely enough". The methods are designed for ease of use: they only require to compute Riemannian exponentials and logarithms, they represent the curve by means of a number of tangent vectors that grows linearly with the number of data points, and, once the representation is computed, evaluating $\gamma(t)$ at any $t$ requires a small number of exponentials and logarithms (independent of the number of data points). Among the proposed methods, the blended cubic spline technique combines the additional properties of interpolating the data when $\lambda \to \infty$ and reducing to the well-known cubic smoothing spline when the manifold is Euclidean. The methods are illustrated on synthetic and real data.
```
@article{Gousenbourger2018,
  author  = {Gousenbourger, Pierre-Yves and Massart, Estelle and Absil, P.-A.},
  title   = {Data fitting on manifolds with composite {B\'e}zier-like curves and blended cubic splines},
  journal = {Journal of Mathematical Imaging and Vision},
  doi     = {10.1007/s10851-018-0865-2},
  volume  = {61},
  number  = {5},
  pages   = {645-671},
  year    = {2018}
}
```
P.-A. Absil, P.-Y. Gousenbourger, P. Striewski, B. Wirth. Differentiable Piecewise-Bézier Surfaces on Riemannian Manifolds. SIAM Journal on Imaging Sciences, 9(4), pp. 1788-1828, 2016.
We propose several methods that address the problem of fitting a $C^1$ curve $\gamma$ to time-labeled data points on a manifold. The methods have a parameter, $\lambda$, to adjust the relative importance of the two goals that the curve should meet: being "straight enough" while fitting the data "closely enough". The methods are designed for ease of use: they only require to compute Riemannian exponentials and logarithms, they represent the curve by means of a number of tangent vectors that grows linearly with the number of data points, and, once the representation is computed, evaluating $\gamma(t)$ at any $t$ requires a small number of exponentials and logarithms (independent of the number of data points). Among the proposed methods, the blended cubic spline technique combines the additional properties of interpolating the data when $\lambda \to \infty$ and reducing to the well-known cubic smoothing spline when the manifold is Euclidean. The methods are illustrated on synthetic and real data.
```
@article{Absil2016,
  author  = {Absil, P.-A. and Gousenbourger, Pierre-Yves and Striewski, Paul and Wirth, Benedikt},
  title   = {{Differentiable Piecewise-B\'ezier Surfaces on Riemannian Manifolds}},
  doi     = {10.1137/16M1057978},
  journal = {SIAM Journal on Imaging Sciences},
  volume  = {9},
  number  = {4},
  pages   = {1788--1828},
  year    = {2016}
}
```

Conference papers

N. T. Son, P.-Y. Gousenbourger, E. Massart, P.-A. Absil. Online balanced truncation for linear time-varying systems using continuously differentiable interpolation on Grassmann manifold Proceedings of the 6th International Conference on Control, Decision and Information Technologies (CoDIT 2019), 2019, to appear.
We consider model order reduction of linear time-varying systems on a finite time interval using balanced truncation. A standard way to perform MOR is to first numerically integrate the associated pair of differential Lyapunov equations for the two gramians, then compute projection matrices using the square root method, and finally formulate the reduced systems at each time instant of a chosen grid. This approach is well-known for delivering good approximation, but rather costly in computation and storage requirement. Furthermore, if one needs to compute the reduced system for any new time instant that is not included in the chosen grid, the mentioned procedure must be performed again without explicitly making use of the already computed data. For dealing with such a situation, we propose to store the projection matrices corresponding to a simplified sparse time grid and to use them to recover the projection subspaces at any other time instant via curve interpolation on the Grassmann manifold. By doing this, we can avoid the repetition of solving the differential Lyapunov equations which is the most expensive step in the procedure and therefore, as shown in a numerical example, accelerate the online reduction process.
```
@inproceedings{Son2019,
  author  = {Son, Nguyen Thanh and Gousenbourger, Pierre-Yves and Massart, Estelle and Absil, P.-A.},
  title   = {Online balanced truncation for linear time-varying systems using continuously differentiable interpolation on Grassmann manifold},
  booktitle = {2019 6th International Conference on Control, Decision and Information Technologies (CoDIT)},
  pages   = {165--170},
  year    = {2019},
  organization = {IEEE}
}
```
E. Massart, P.-Y. Gousenbourger, N. T. Son, T. Stykel, P.-A. Absil. Interpolation on the manifold of fixed-rank positive-semidefinite matrices for parametric model order reduction: preliminary results. ESANN2019, submitted.
We present several interpolation schemes on the manifold of fixed-rank positive-semidefinite (PSD) matrices. We explain how these techniques can be used for model order reduction of parameterized linear dynamical systems, and obtain preliminary results on an application.
```
@inproceedings{Massart2019a,
  author  = {Massart, Estelle and Gousenbourger, Pierre-Yves and Son, Nguyen Thanh and Stykel, Tatjana and Absil, P.-A.},
  title   = {Interpolation on the manifold of fixed-rank positive-semidefinite matrices for parametric model order reduction: preliminary results},
  booktitle = {ESANN2019},
  pages   = {281--286},
  isbn    = {978-287-587-065-0},
  publisher = {Springer},
  year    = {2019}
}
```
P.-Y. Gousenbourger, L. Jacques, P.-A. Absil. Fast method to fit a $C^1$ piecewise-Bézier function to manifold-valued data points: how suboptimal is the curve obtained on the sphere $\mathbb S^2$? Nielsen, F. and Barbaresco F., editors, GSI2017, pp. 595-603, 2017.
We propose an analysis of the quality of the fitting method proposed in [7]. This method fits smooth paths to manifold-valued data points using $C^1$ piecewise-B ́ezier functions. This method is based on the principle of minimizing an objective function composed of a data-attachment term and a regularization term chosen as the mean squared acceleration of the path. However, the method strikes a tradeoff between speed and accuracy by following a strategy that is guaranteed to yield the optimal curve only when the manifold is linear. In this paper, we focus on the sphere $S^2$. We compare the quality of the path returned by the algorithms from [7] with the path obtained by minimizing, over the same search space of $C^1$ piecewise-B́ezier curves, a finite-difference approximation of the objective function by means of a derivative-free manifold-based optimization method.
```
@inproceedings{Gousenbourger2017a,
  author  = {Gousenbourger, Pierre-Yves and Jacques, Laurent and Absil, P.-A.},
  editor  = {Nielsen, Frank and Barbaresco, Fr{\'e}d{\'e}ric},
  title   = {Fast method to fit a {$C^1$} piecewise-{B\'{e}}zier function to manifold-valued data points: how suboptimal is the curve obtained on the sphere {$\mathbb S^2$}?},
  booktitle = {Geometric Science of Information},
  doi     = {10.1007/978-3-319-68445-1_69},
  publisher = {Springer},
  address = {Berlin, Heidelberg},
  series  = {Lecture Notes in Computer Sciences},
  volume  = {10589},
  pages   = {595--603},
  year    = {2017}
}
```
P.-Y. Gousenbourger, E. Massart, A. Musolas, P.-A. Absil, L. Jacques, J.M. Hendrickx, Y. Marzouk. Piecewise-Bézier $C^1$ smoothing on manifolds with application to wind field estimation. ESANN2017, Springer eds,, pp. 305-301, 2017.
We propose an algorithm for fitting $C^1$ piecewise-Bézier curves to (possibly corrupted) data points on manifolds. The curve is chosen as a compromise between proximity to data points and regularity. We apply our algorithm as an example to fit a curve to a set of low-rank covariance matrices, a task arising in wind field modeling. We show that our algorithm has denoising abilities for this application.
```
@inproceedings{Gousenbourger2017,
  author  = {Gousenbourger, Pierre-Yves and Massart, Estelle and Musolas, Antoni and Absil, P.-A. and Jacques, Laurent and Hendrickx, Julien M and Marzouk, Youssef},
  booktitle = {ESANN2017},
  pages   = {305--310},
  publisher = {Springer},
  title   = {Piecewise-{B}{\'e}zier {$C^1$} smoothing on manifolds with application to wind field estimation},
  year    = {2017}
}
```
P.-A. Absil, P.-Y. Gousenbourger, P. Striewski, B. Wirth. Differentiable piecewise-Bézier interpolation on Riemannian manifolds. ESANN2016, Springer eds,, pp. 95-100, 2016.
We propose a generalization of classical Euclidean piecewise-Bézier surfaces to manifolds, and we use this generalization to compute a $C^1$-surface interpolating a given set of manifold-valued data points associated to a regular 2D grid. We then propose an efficient algorithm to compute the control points defining the surface based on the Euclidean concept of natural $C^2$-splines and show examples on different manifolds.
```
@inproceedings{Absil2016a,
  author    = {Absil, P.-A. and Gousenbourger, Pierre-Yves and Striewski, Paul and Wirth, Benedikt},
  title     = {Differentiable piecewise-{B\'{e}}zier interpolation on {R}iemannian manifolds},
  booktitle = {ESANN2016},
  isbn      = {978--287587027--8},
  publisher = {Springer},
  pages     = {95--100},
  year      = {2016}
}
```
C. Samir, P.-Y. Gousenbourger, S. H. Joshi. Cylindrical Surface Reconstruction by Fitting Paths on Shape Space. EProceedings of the 1st International Workshop on DIFFerential Geometry in Computer Vision for Analysis of Shapes, Images and Trajectories (DIFF-CV 2015), 2015.
We present a differential geometric approach for cylindrical anatomical surface reconstruction from 3D volumetric data that may have missing slices or discontinuities. We extract planar boundaries from the 2D image slices, and parameterize them by an indexed set of curves. Under the SRVF framework, the curves are represented as invariant elements of a nonlinear shape space. Differently from standard approaches, we use tools such as exponential maps and geodesics from Riemannian geometry and solve the problem of surface reconstruction by fitting paths through the given curves. Experimental results show the surface reconstruction of smooth endometrial tissue shapes generated from MRI slices.
```
@inproceedings{Samir2015,
  author = {Samir, Chafik and Gousenbourger, Pierre-Yves and Joshi, Shantanu H.},
  booktitle = {DIFF-CV 2015},
  pages  = {1--10},
  title  = {{Cylindrical Surface Reconstruction by Fitting Paths on Shape Space}},
  year   = {2015}
}
```
A. Arnould, P.-Y. Gousenbourger, C. Samir, P.-A. Absil, M. Canis. Fitting smooth paths on Riemannian manifolds: Endometrial surface reconstruction and preoperative MRI-based navigation. Nielsen, F. and Barbaresco F., editors, GSI2015, pp. 491-498, 2015.
We present a new method to fit smooth paths to a given set of points on Riemannian manifolds using $C^1$ piecewise-Bézier functions. A property of the method is that, when the manifold reduces to a Euclidean space, the control points minimize the mean square acceleration of the path. As an application, we focus on data observations that evolve on certain nonlinear manifolds of importance in medical imaging: the shape manifold for endometrial surface reconstruction; the special orthogonal group $SO(3)$ and the special Euclidean group $SE(3)$ for preoperative MRI-based navigation. Results on real data show that our method succeeds in meeting the clinical goal: combining different modalities to improve the localization of the endometrial lesions.
```
@inproceedings{Arnould2015,
  author  = {Arnould, Antoine and Gousenbourger, Pierre-Yves and Samir, Chafik and Absil, P.-A. and Canis, Michel},
  title   = {{Fitting Smooth Paths on Riemannian Manifolds : Endometrial Surface Reconstruction and Preoperative MRI-Based Navigation}},
  booktitle = {Geometric Science of Information},
  doi     = {10.1007/978-3-319-25040-3_53},
  editor  = {Nielsen, Frank and Barbaresco, Fr{\'e}d{\'e}ric},
  publisher = {Springer},
  address = {Berlin, Heidelberg},
  series  = {Lecture Notes in Computer Sciences},
  volume  = {9389},
  pages   = {491--498},
  year    = {2015}
}
```
P.-Y. Gousenbourger, C. Samir, P.-A. Absil. Piecewise-Bézier $C^1$ interpolation on Riemannian manifolds with application to 2D shape morphing. International Conference on Pattern Recognition (ICPR), pp. 4086-4091, 2014.
We propose a generalization of classical Euclidean piecewise-Bézier surfaces to manifolds, and we use this generalization to compute a $C^1$-surface interpolating a given set of manifold-valued data points associated to a regular 2D grid. We then propose an efficient algorithm to compute the control points defining the surface based on the Euclidean concept of natural $C^2$-splines and show examples on different manifolds.
```
@inproceedings{Gousenbourger2014,
  author = {Gousenbourger, Pierre-Yves and Samir, Chafik and Absil, P.-A.},
  title  = {Piecewise-{B\'e}zier {$C^1$} interpolation on {R}iemannian manifolds with application to {2D} shape morphing},
  booktitle = {International Conference on Pattern Recognition (ICPR)},
  pages  = {4086--4091},
  year   = {2014},
  doi    = {10.1109/ICPR.2014.700},
  isbn   = {9781479952083},
  issn   = {10514651}
}
```

Talks

Data fitting on manifolds: applications, challenges and solutions
ISP Group, Louvain-la-Neuve, Belgium. Dec. 11, 2019.
Storm trajectories prediction, birds migrations follow-up, rigid rotations of 3D objects, wind field estimation, model order reduction of superlarge parameter-dependent dynamical systems, MRI 3D body volumes reconstruction... All these applications have two things in common: first, they have a geometrical data-structure, i.e., the data lives on a (generally) Riemannian manifold; second, they can benefit of parameter(s)-dependent fitting methods somewhere in the process. If data fitting is a basic problem in the Euclidean space (where natural cubic splines and thin plates splines are the superstars in the domain), it become more intricate when data structure constrains the problem.

This talk is an opportunity to present you an efficient, ready-to-use algorithm for data fitting on manifolds based on B\'ezier curves, applied to some of the aforementioned applications.
Data fitting on manifolds by minimizing the mean square acceleration of a Bézier curve
Benelux Meeting 2019, Lommel, Belgium. Mar. 19, 2019.
Curve fitting on manifolds with Bézier and blended curves
Seminar at Research Group "Numerical Mathematics (Partial Differential Equations)", from TU Chemnitz, Chemnitz, Germany. Jul. 18, 2018.
Fast method to fit a $C^1$ piecewise-Bézier function to manifold-valued data points: how suboptimal is the curve obtained on the sphere $\mathbb S^2$?
GSI2015, Paris, France. Nov. 7, 2017.
We propose an analysis of the quality of the fitting method proposed in [7]. This method fits smooth paths to manifold-valued data points using $C^1$ piecewise-B ́ezier functions. This method is based on the principle of minimizing an objective function composed of a data-attachment term and a regularization term chosen as the mean squared acceleration of the path. However, the method strikes a tradeoff between speed and accuracy by following a strategy that is guaranteed to yield the optimal curve only when the manifold is linear. In this paper, we focus on the sphere $S^2$. We compare the quality of the path returned by the algorithms from [7] with the path obtained by minimizing, over the same search space of $C^1$ piecewise-B́ezier curves, a finite-difference approximation of the objective function by means of a derivative-free manifold-based optimization method.
Wind field estimation via $C^1$ Bézier smoothing on manifolds
Workshop on "Regularized Inverse Problem Solving and High-Dimensional Learning Methods" (WIPS), Louvain-la-Neuve, Belgium. Aug. 30, 2017.
Unmanned aerial vehicle control is a hot topic in research and at the crossroad of a lot of disciplines. For instance, safe and reliable navigation of UAVs requires consideration of the surrounding environment, in particular, the external wind conditions. This external wind configuration is usually evaluated by computationnaly expensive efforts and depends on external meteorological parameters. Hopefully, it can be modelled as a Gaussian process characterized by a covariance matrix belonging to the space of PSD matrices of rank r. In this work, we both expoit the manifold structure of this specific space and also propose a method to fit a small set of pre-computed solutions. That way, for a new value of the external meteorological parameters, we are able to recover a sufficiently accurate wind field configuration in a computationnaly tractable effort. Our method is based on manifold-valued Bezier curves. Joint work with E.M. Massart, A. Musolas, P.-A. Absil, J.M. Hendrickx, L. Jacques and Y. Marzouk.
Wind field estimation via $C^1$ Bézier smoothing on manifolds
ESANN2017, Bruges, Belgium. Apr. 27, 2017.
We propose an algorithm for fitting $C^1$ piecewise-Bézier curves to (possibly corrupted) data points on manifolds. The curve is chosen as a compromise between proximity to data points and regularity. We apply our algorithm as an example to fit a curve to a set of low-rank covariance matrices, a task arising in wind field modeling. We show that our algorithm has denoising abilities for this application.
Wind field estimation via $C^1$ Bézier smoothing on manifolds
Benelux Meeting 2017, Spa, Belgium. Mar. 28, 2017.
Interpolation and fitting on manifolds with differentiable piecewise-Bézier functions
ISP Group, Louvain-la-Neuve, Belgium. Mar. 16, 2017.
Fitting and interpolation are well known topics on the Euclidean space. However, when it turns out that the data points are manifold valued (understand: when the data point belongs to a certain manifold), most of the algorithms are no more applicable because even the notion of distance is not as simple as on the Euclidean space.

Manifold-valued data appear in various domains as MRI acquisition (diffusion tensors, segmented shapes,...), fluid mechanics (SDP matrices or values on the Grassmann manifold), matrix completion and many more. A particular example concern the representation of local wind fields in a given domain. These can be expressed as covariance matrices which are SDP and of fixed (small) rank. The wind field changes when the prevalent wind has different orientation or magnitude. The computation of this wind field, however, resort in solving complicated and time-consuming PDEs in fluid mechanics. In order to obtain approached solutions as fast as possible, a solution is to compute several wind field covariance matrices representations for a set of given orientation/magnitude of the prevalent wind and then interpolate (or fit) a manifold-valued curve in the space of fixed rand SDP matrices.

In this talk, I will propose a method to fit manifold-valued data points based on a generalization of Bezier curves. This method is motivated by (but not specific to) the wind field orientation problem. I will try to introduce the concepts of manifolds and Bezier curves as a tutorial to finally show interesting results for that specific application.
Interpolation and fitting on manifolds with differentiable piecewise-Bézier functions
Anuj Srivastava's working group meeting, Skype presentation, Tallahassee, Florida, USA. Mar. 10, 2017.
Differentiable Bézier interpolation on manifolds with B-splines
GAMM 2017, Weimar, Germany. Mar. 6, 2017.
Interpolation and fitting on manifolds with differentiable piecewise-Bézier functions
Workshop on Manifold-Valued Image Processing (MVIP) 2016, Kaiserslautern, Germany. Dec. 2, 2016.
Interpolation on manifolds with differentiable surfaces of Bézier
Benelux Meeting 2016, Soesterberg, The Netherlands. Mar. 24, 2016.
Bézier interpolation on Riemannian manifolds
ASCII Tutorial Seminars, Louvain-la-Neuve, Belgium. Nov. 27, 2015.
Mani-what you say? Manifolds are mathematical sets with a smooth geometry (such as spheres) used in applications like model reduction, optimization and many more. What do manifold-valued objects look like? Can we interpolate them? Is it costly? With Bézier functions, interpolation on manifolds becomes very easy. Bézier? You should remember this from your bachelor, isn't it? From a very few concepts and drawings, a variety of possibilities is open.
Endometriosis: MRI navigation and surface reconstruction on manifolds
GSI2015, Paris, France. Oct. 30, 2015.
We present a new method to fit smooth paths to a given set of points on Riemannian manifolds using $C^1$ piecewise-Bézier functions. A property of the method is that, when the manifold reduces to a Euclidean space, the control points minimize the mean square acceleration of the path. As an application, we focus on data observations that evolve on certain nonlinear manifolds of importance in medical imaging: the shape manifold for endometrial surface reconstruction; the special orthogonal group $SO(3)$ and the special Euclidean group $SE(3)$ for preoperative MRI-based navigation. Results on real data show that our method succeeds in meeting the clinical goal: combining different modalities to improve the localization of the endometrial lesions.
How to design interpolating surfaces on manifolds?
One-Day Workshop on manifolds optimization, Louvain-la-Neuve, Belgium. Sep. 25, 2015.
Interpolation on Riemannian manifolds with a $C^1$ piecewize-Bézier path.
GDR MIA 2014, Paris, France. Nov. 21, 2014.
Interpolation on Riemannian manifolds with a $C^1$ piecewize-Bézier path.
ISP Group, Louvain-la-Neuve, Belgium. Oct. 8, 2014.
Nowadays, more and more problems are solved through specific manifold formulation. This often allows important reduction of computation time and/or memory management compared to classical formulations on the classical Euclidean space (because of non-linear constraints like restricting the solutions to a certain subdomain of a larger ambiant space). Interpolation and optimization tools can be useful for solving some of these problems (like defining the optimal trajectory of a humanitory plane dropping supplies, or fitting two objects orientations). However, current procedures are only defined on the Euclidean space. In this presentation, I focus on interpolation methods and, more precisely, I propose a new general framework to fit a path through a finite set of data points lying on a Riemannian manifold. The path takes the form of a continuously-differentiable concatenation of Riemannian Bézier segments. This framework will be illustrated by results on the Euclidean space, the sphere, the orthogonal group and the shape manifold. The content of this presentation meets also very recent research carried out in this institute for providing novel efficient manifold-based optimization methods.

Posters

Interpolation on the manifold of fixed-rank positive-semidefinite matrices for parametric model order reduction: preliminary results.
Conference ESANN 2019, Bruges, Belgium, 2019.
We present several interpolation schemes on the manifold of fixed-rank positive-semidefinite (PSD) matrices. We explain how these techniques can be used for model order reduction of parameterized linear dynamical systems, and obtain preliminary results on an application.
Blended smoothing splines on Riemannian manifolds.
Conference iTwist 2018, Marseille, France, 2018.
We present a method to compute a fitting curve $B$ to a set of data points $d_0,\dots, d_m$ lying on a manifold $\mathcal M$. That curve is obtained by blending together Euclidean Bézier curves obtained on different tangent spaces. The method guarantees several properties among which $B$ is $C^1$ and is the natural cubic smoothing spline when $\mathcal M$ is the Euclidean space. We show examples on the sphere $S^2$ as a proof of concept.
Data fitting on manifolds with blended cubic splines.
ICML2018, workshop GiMLi, Stockholm, Sweden, 2018.
Interpolation on manifolds using Bézier functions.
Conference iTwist 2018, Aalborg, Denmark, 2016.
Given a set of data points lying on a smooth manifold, we present methods to interpolate those with piecewise Bézier splines. The spline is composed of Bézier curves (resp. surfaces) patched together such that the spline is continuous and differentiable at any point of its domain. The spline is optimized such that its mean square acceleration is minimized when the manifold is the Euclidean space. We show examples on the sphere $S^2$ and on the special orthogonal group $SO(3)$.
Differentiable piecewise-Bézier interpolation on Riemannian manifolds.
Conference ESANN2016, Bruges, Belgium, 2016.
We propose a generalization of classical Euclidean piecewise-Bézier surfaces to manifolds, and we use this generalization to compute a $C^1$-surface interpolating a given set of manifold-valued data points associated to a regular 2D grid. We then propose an efficient algorithm to compute the control points defining the surface based on the Euclidean concept of natural $C^2$-splines and show examples on different manifolds.
Piecewise-Bézier $C^1$ interpolation on Riemannian manifolds with application to 2D shape morphing.
Dysco Day, Gent, Belgium, 2014.
Piecewise-Bézier $C^1$ interpolation on Riemannian manifolds with application to 2D shape morphing.
International Conference on Pattern Recognition (ICPR), pp. 4086-4091, 2014.
We present a new framework to fit a path to a given finite set of data points on a Riemannian manifold. The path takes the form of a continuously-differentiable concatenation of Riemannian Bézier segments. The selection of the control points that define the Bézier segments is partly guided by the differentiability requirement and by a minimal mean squared acceleration objective. We illustrate our approach on specific manifolds: the Euclidean plane (for sanity check), the sphere (as a first nonlinear illustration), the special orthogonal group (with rigid body motion applications), and the shape manifold (with 2D shape morphing applications).

Download the newest version of the code

Bidimensional interpolation with Bézier spline

We are given a set of data points in the Euclidean space associated with two timestamps (m,n) organised on a regular grid. These timestamps give the order in which points must be interpolated. In this problem, for visualization convenience and without loss of generality, the (x,y) value of the data points correspond to the regular grid (m,n) and the z-value is randomly chosen. The points are organized in a 2D-cell corresponding to the order of interpolation.

From these data points, the goal here is to draw a 2d-Bézier spline interpolating those in a $C^1$ way. For this, you can use our code here, general for manifolds. The possible manifolds are the Euclidean space, the sphere and the special orthogonal group SO(3).

Download the code.

Show/hide the code tutorial.

The code

Here is the code for interpolating data with a 2d-Bézier spline.

% Some random data points.
m = 6; n = 7;
data = cell(m,n);
for i=1:m
for j=1:n
  data{i,j} = [i-1,j-1,randn(1)];
end
end

% Some useful information...
manifold = 'euclidean';
set_path
global_variables
geo_functions

% Create the problem structure.
problem  = prepare_structure(data,manifold,10);

% Phase 1: Control points generation.
problem = control_points_simple_generation_2d(problem);

% Phase 2: Reconstruction of the curve.
problem = curve_reconstruction_double_bezier_c1(problem,1); % horizontal-vertical

% Plotting
figure;
draw_bezier_surface(problem);
title(['Bezier surface with method: ',problem.method,'.']);

Okay, let's dig a bit into that...

Data points

As mentionned, the data points must be stored in a 2d-cell. Let's construct a random set of data points where the (x,y)-values correspond to (m,n) and the z-value is random.

Hint: You may also call data_points_2d(manifold,type) to create some predefined test sets.

m = 6; n = 7;
data = cell(m,n);
for i=1:m
for j=1:n
  data{i,j} = [i-1,j-1,randn(1)];
end
end

Problem structure preparation

The system is based on a structure called successively by different subfunctions in matlab. Each subfunction will use elements of differential geometry (the exponential, the logarithm, the distance, etc.). This is why this preparation step is crucial in order to inform the computer on the paths needed (the script set_path depends on the manifold) and on the elements of differential geometry available for the code (global_variables and geo_functions depend on the manifold).

manifold = 'euclidean';

% define the paths where the manifold-tools are stored
set_path
% define the global variables, i.e. the elements of differential geometry.
global_variables
geo_functions

Create the problem structure.

% the structure stored the name of the manifold, the data points and the discretization of the surface when reconstructed.
problem  = prepare_structure(data,manifold,10);

The structure should look like

problem = 				
interp: {6x7 cell}
type: ''
manifold: 'euclidean'
   n: 6
   m: 7
dim1: 1
dim2: 3
   d: []
   t: 10
nint: 0

Bézier spline optimization and reconstruction

The interpolation problem works in two steps:

find the control points leading the Bézier curve
reconstruct the actual spline.

Phase 1: Control points generation

The subfunction takes the structure as entry and will update it by adding the computed control points.

problem = control_points_simple_generation_2d(problem);

There exist different methods to compute the control points:

control_points_generation(problem): By an exact technique minimizing the Euclidean acceleration of the path. Warning, this technique is costly in time when the log and exp are not trivial. The technique is documented in the technical report of 2015.
control_points_simple_generation_2d(problem) With an efficient but not perfect algorithm based on the B-spline representation of curves generalized to surfaces. This technique is documented in ESANN2016
control_points_double_tensorization(problem) With the same efficient algorithm but applied to 1D curves horizontally and then vertically. The technique is also documented in ESANN2016

In the Euclidean space, all these techniques are identical.

Phase 2: Reconstruction of the curve

Here also, the subfunction takes the structure as entry and will update it by adding the computed spline. If the structure does not contain the control points, the reconstruction is not possible.

problem = curve_reconstruction_double_bezier_c1(problem,1); % horizontal-vertical

There also exist different techniques to reconstruct the Bézier spline. The following ones return a C^1 surface. They are documented in technical report of 2015.

curve_reconstruction_mean_all_c1(problem) uses the altered defintion of Bezier surfaces based on averaging of all control points of the patch.
curve_reconstruction_double_bezier_c1(problem,k) uses the altered definition of Bezier curves horizontal-vertical (k=1) or vertical-horizontal (k=2) depending on the argument.

The two following ones return a C^0 surface and their limits are also discussed in the technical report of 2015.

curve_reconstruction_mean_all(problem) uses the basic definition of Bezier surfaces based on averaging of all the points of the patch.
curve_reconstruction_double_bezier(problem,k) uses the basic definition of Bezier curves horizontal-vertical (1) or vertical-horizontal (2) depending on the argument.

Plotting

We propose a method to plot the solution. The plotting method is chosen based on the manifold and on the paths set at the beginning.

figure;
draw_bezier_surface(problem);
title(['Bezier surface with method: ',problem.method,'.']);

Unidimensional interpolation with hybrid composite Bézier curves

We are given a set of data points in the Euclidean space associated with a timestamp $t_i \in \mathbb{Z}$. We would like to compute a composite Bézier curve that interpolates the data points $d_i$ at time $t_i$. In this problem, we try to interpolate four data points $d_i \in \mathbb{S}^2$, the unit sphere. These points are places such that they represent a triangle.

To fit the curve $\mathfrak{B}(t): [0,3] \to \mathbb{S}^2$, you can use our code here, general for manifolds. The possible manifolds are the Euclidean space, the sphere, the special orthogonal group SO(3),... actually all manifolds you can find in Manopt are operational. It is however mandatory to add manopt to your working path. The zip file hereunder contains an old version of Manopt for consistency. You can download the last version of Manopt here (recommended) : www.manopt.org.

Download the code.

Show/hide the code tutorial.

The code

Here is the code for interpolating data with a hybrid composite Bézier curve on the sphere.

% Add the paths for manopt and for the tools
cd manopt;
addpath(genpath(pwd));
cd ..;
addpath(genpath([pwd,'/tools']));
	
% Some data points on the sphere.
phi     = [pi/3 ; 4*pi/7 ; 4*pi/7 ; pi/3] - 2*pi/12; 
theta   = [pi/2 ; pi/4 ; 3*pi/4 ; pi/2] + 10*pi/12;
x       = cos(theta).*sin(phi); y = sin(theta).*sin(phi); z = cos(phi);
points  = [x y z];
data    = reshape(points',1,3,4);

% Problem parameters
n = size(data,3);
t = linspace(0,n-1,100);

% Structure creation
problem.M = spherefactory(3);
problem.data = data;

% Phase 1: Control points generation.
problem.control = cp_arnould(problem);

% Phase 2: Reconstruction of the curve.
problem = bezier_reconstruction(problem,t);

% Optional: velocity and acceleration computation.
problem.velocity = bezier_velocity(problem,t);
problem.acceleration = bezier_acceleration(problem,t);
	
% Plotting
y = squeeze(problem.curve)'; b = squeeze(problem.control)'; p = squeeze(problem.data)';	
a = problem.acceleration; v = problem.velocity;

[X,Y,Z] = sphere(); 
figure;
subplot(2,2,[1 3]); surf(0.95*X,0.95*Y,0.95*Z,...
    'FaceAlpha',0.7,'EdgeAlpha',1,...
    'FaceColor', [238 197 145]/255);
  hold on;
  plot3(y(:,1),y(:,2),y(:,3),'b','LineWidth',2);
  plot3(b(:,1),b(:,2),b(:,3),'.','Color',[0 0.7 0],'Markersize',30);
  plot3(p(:,1),p(:,2),p(:,3),'.r','Markersize',30);
  title('Bezier spline');
subplot(222); plot(t(2:end-1),v,'-b','LineWidth',2); title('velocity');
subplot(224); plot(t(2:end-1),a,'-b','LineWidth',2); title('acceleration');

Okay, let's dig a bit into that...

Data points

The data points must be stored in a (l,c,n) matrix, where (l,c) is the matrix size of the embedded space. Here, for the sphere, it would be (1,3), as $\mathbb{S}^2$ is embedded in $\mathbb{R}^3$. The points here are draw a triangle on the sphere, and the first and last points coincide.

Hint: You may also call data_points(manifold,type) to create some predefined test sets. Do not forget to add the path to it (addpath([pwd,'/data_points']);).

% Some data points on the sphere.
phi     = [pi/3 ; 4*pi/7 ; 4*pi/7 ; pi/3] - 2*pi/12; 
theta   = [pi/2 ; pi/4 ; 3*pi/4 ; pi/2] + 10*pi/12;
x       = cos(theta).*sin(phi); y = sin(theta).*sin(phi); z = cos(phi);
points  = [x y z];
data    = reshape(points',1,3,4);

Problem structure preparation

The system is based on a structure called successively by different subfunctions in matlab. Each subfunction will use the elements of differential geometry (the exponential, the logarithm, the distance, etc.) extracted from Manopt and stored in that structure. The names of the different elements of the structure will be very important as they are used by the rest of the code later on.

We prepare first the time parameter t to be able, later on, to compute the reconstruction spline. Typically, it is a discretization of the domain (here $t \in [0,3]$).

% Problem parameters
n = size(data,3);
t = linspace(0,n-1,100);

Now, we store the crucial information in a structure.

% Structure creation
problem.M = spherefactory(3);
problem.data = data;

The structure problem is supposed to be at the end composed of the following elements (in this specific problem):

problem = 				
            data: [1×3×4 double]
               M: [1×1 struct]
         control: [1×3×8 double]
           curve: [1×3×100 double]
        velocity: [98×1 double]
    acceleration: [98×1 double]

The data points are stored in the field data, the control points in control, the curve in curve and the manifold in M. It is important to respect these four names. The other names have less importance.

Bézier curve optimization and reconstruction

The interpolation problem works in two steps:

find the control points leading the Bézier curve
reconstruct the actual spline.

Phase 1: Control points generation

The subfunction takes the structure as entry and will return the control points. Store it in the structure for reconstruction.

problem.control = cp_arnould(problem);

There exist different methods to compute the control points:

cp_arnould(problem): technique proposed in the paper Arnould et al. 2015. This technique is recommended.
cp_gousenbourger(problem): suboptimal technique proposed in the paper Gousenbourger et al. 2014, where the velocity vectors must be specified in the beginning of the method. It works but... well, there is better, so why not using it ? ;-).

Phase 2: Reconstruction of the curve

Here also, the subfunction takes the structure as entry and will update it by adding the computed spline. If the structure does not contain the field control, filled with the control points returned by the control points generation method, the reconstruction is not possible.

problem = bezier_reconstruction(problem,t);

The method is a simple De Casteljau algorithm adapted to compute a hybrid composite Bézier curve.

Optional phase: Acceleration and velocity of the curve

Once the curve is reconstructed, we offer a generic code to evaluate the velocity and the acceleration of the curve. Just use the following piece of code:

a = bezier_velocity(problem,t);
v = bezier_acceleration(problem,t);

Be sure that the field curve is available in the structure problem.

The rest is just plotting the curve on a sphere. Try it :-) !

License

The code provided here is copyright and distributed under the terms of the GNU General Public License (GPL) version 3 (or later).

In short, this means that everyone is free to use the provided code, to modify it and to redistribute it on a free basis. However, this code is not in the public domain; it is copyrighted and there are restrictions on its distribution (see the license and the related frequently asked questions). For example, you cannot integrate this version of the code (in full or in parts) in any closed-source software you plan to distribute (commercially or not). Please contact me for more information.

This section was inspired by the very good Gabriel Peyré's Numerical Tours.