Pietro Gori

Jean-Marc Thiery

Damien Rohmer Je

Je remercie les membres du jury pour avoir pris le temps de lire ce manuscrit malgré des délais restreints, ainsi que pour les nombreux retours positifs durant la soutenance.

We propose to progressively simplify tractograms by grouping similar bers into generalized cylinders that we call meta-bers. This generates a ne-grained multiresolution model that provides progressive and real-time navigation through different levels of detail. This model preserves the overall structure of the tractogram and can be built upon di erent similarity measures. A Delaunay tetrahedralization is used to infer neighboring connectivity, thus reducing the complexity of our algorithm. For high numbers of bers, our approach is faster and geometrically more accurate than usual simpli cation algorithms.

There is a need for a speci c geometric representation for bundles and metabers. In this work, we use generalized cylinders, Moving Least Square (MLS) projections, and implicit surfaces. The global MLS projection being identi ed as a very promising approach, we propose a method to reduce its complexity, while allowing for large kernel usage, making it possible to use on highly-challenging data e ciently. This opens a way toward global kernels analysis or modeling, which was previously unthinkable because of the complexity. We then investigate its use for the speci c case of supporting the various geometric operations required for building the previously presented multi-resolution representation.

In order to allow for individual analysis of each bundle, we introduce a segmentation technique based on the multi-resolution representation, combined with anatomical knowledge modeled in the framework of fuzzy sets, to de ne an anatomical coherence score. This score can be thresholded interactively, for ne-tuning of the nal bundle of interest. This technique manages better reproducibility than previous approaches.

Tractogram les can take a huge amount of space, making them di cult to transfer or even store. Using state of the art unit vector compression, combined with bers properties, we present a compression algorithm both fast and e cient, with a 10 times le size reduction, while preserving satisfying accuracy criteria. This algorithm can even be adapted for on-the-y GPU decompression of the data. Using this breakthrough compression, it is now possible to visualize datasets that would normally not t in memory.

We obtain some interesting results for geometric representations for tractograms that would bene t from some more investigations, especially for a better adaptation of MLS to bers, as well as for visualization techniques.

Résumé

L'Imagerie de di usion par Résonance Magnétique (IRM), associée à des algorithmes de tractographie, est la seule technique d'imagerie non-invasive capable de reconstruire les connexions du cerveau humain à l'intérieur de la matière blanche. Il en résulte de nombreuses polylignes en 3D, appelées bres, qui peuvent être utilisées dans des environnements cliniques ou en recherche. Ces bres sont des estimations des milliards de neurones présents dans le cerveau. Les méthodes actuelles sont incapables de reconstruire les neurones individuels. En e et, ils mesurent quelques micromètres de large, là où les images IRM ont des résolutions typiquement millimétriques. Cependant, les tractogrammes obtenus peuvent être composés de millions de bres, rendant la visualisation et l'interprétation extrêmement di ciles.

Dans cette thèse, nous nous e orçons de rendre la tractographie plus facile à utiliser et analyser, en termes de simpli cation, de représentation géométrique, d'espace de stockage et de visualisation. Les di érents travaux présentés ici se positionnent à l'intersection entre l'informatique graphique (en particuler, la modélisation géométrique et le rendu) et l'analyse de données médicales, et ont été développés pour répondre aux dé s mentionnés précédemment.

Nous proposons de simpli er progressivement les tractogrammes en groupant les bres similaires en des cylindres généralisés que nous appelons méta-bres. Cette approche construit un modèle multi-résolution à grain n qui permet une navigation en temps réel à di érents niveaux de détails. Ce modèle préserve la structure globale du tractogramme et peut être généré à partir de di érentes mesures de similarité. Une tétraédrisation de Delaunay est utilisée pour déterminer des connectivités voisines, réduisant la complexité de notre algorithme. Pour de grandes quantité de bres, notre approche est plus rapide et géométriquement plus précise que les algorithmes de simpli cation existants.

Une géométrie de représentation spéci que pour les faisceaux et les méta-bres est nécessaire. Dans ce travail, nous utilisons les cylindres généralisés, les déplacements aux moindres carrés (MLS) ainsi que les surfaces implicites. La projection MLS globale étant identi ée rapidement comme étant une approche très prometteuse, nous proposons une méthode pour réduire sa complexité, permettant l'utilisation de larges noyaux, rendant possible son utilisation e cace sur des données très complexes. Cela ouvre la voie à l'analyse et la modélisation de noyaux globaux, ce qui était impensable précédemment en raison de la complexité. Nous explorons ensuite son utilisation pour le cas spéci que du support de diverses opérations géométriques nécessaires à la construction de la représentation multi-résolution.

A n de pouvoir analyser individuellement chaque faisceau de bres, nous introduisons une technique de segmentation à partir de la représentation multi-résolution, combinée à des connaissances anatomiques modélisées dans la théorie des ensembles ous, a n de dé nir un score de cohérence anatomique. Ce score peut être seuillé de manière interactive, a n d'a ner le faisceau d'intérêt. Cette technique permet une meilleure reproducibilité comparée aux approches précédentes.

Les chiers de tractogrammes peuvent prendre une quantité importante d'espace, les rendant di ciles à transférer, voire à stocker. En utilisant l'état de l'art en compression de vecteurs unitaires, combiné aux propriétés des bres, nous présentons un algorithme de compression à la fois rapide et e cace, proposant une réduction par 10 de la taille des chiers, tout en préservant une précision satisfaisante. Cet algorithme peut même être adapté à une décompression des données à la volée sur le GPU. Cela permet de visualiser des jeux de données qui ne pourraient normalement pas rentrer en mémoire.

Nous obtenons des résultats intéressants pour la représentation géométrique de tractogrammes qui béné cieraient de recherches supplémentaires, particulièrement pour une meilleure adaptation des projections MLS aux bres, ainsi que pour les techniques de visualisation.

I Introduction

I.1 Medical context

Di usion Magnetic Resonance Imaging (dMRI) is an MRI technique allowing for the measurement of water di usion inside a living human brain. This allows for the reconstruction through tractography algorithms of some sets of connections inside the white matter [START_REF] Tournier | Di usion tensor imaging and beyond[END_REF], Soares et al., 2013]. This is explained with more details in Section III.1. These sets of connections are called bers or streamlines and represent collections of axons of the neurons. These axons connect neurons together through the white matter. In this thesis, we use the word ber to refer to these obtained sets of axons. This kind of data is really useful to help understand the brain. Moreover, it is used during pre-operative surgery for tumor removals for instance [START_REF] Soares | A hitchhiker's guide to di usion tensor imaging[END_REF]. Indeed, knowing where these connections are positioned in the patient's brain helps to prevent some complications, as the path used to remove the tumor can be optimized based on the bers localization.

These sets of bers are called tractograms. They are challenging to use because of their complexity, as they are commonly reconstructed with millions of bers, each one being a 3D polyline containing hundreds of points. Consequently, their usage without any prior adaptation is di cult, or even impossible sometimes. Indeed, the information they convey is numerous, making it di cult to apprehend with simple approaches. Their size of several gigabytes is even preventing their display in some cases, requiring some prior modi cations or simpli cations.

Fortunately, some redundancy is present in tractograms, as lots of similar bers can be obtained. It opens a way towards simpli cation approaches, on which we base the work of this thesis, which focuses on creating some useful algorithms to make tractography easier to use. These algorithms are aiming at simplifying tractograms, nding useful tracts of interest, or compressing the bers for lower memory requirements.

I.2 Computer graphics context

In order to simplify, compress, or visualize the obtained bers, some tools and techniques are needed. We will take inspiration from existing methods from the computer graphics domain for mesh simpli cations, geometric representations of data, or visualization techniques. These tools will be adapted and improved for our speci c data: brain tractogram.

The main elements that will be necessary in the context of this thesis are simplication approaches for sets of curves, which will need e cient and correct ways of comparing di erent curves with some constraints dictated by the nature of the data used. Simplifying sets of similar bers requires speci c geometric representations, some of which will be explored in this thesis. We call meta-bers these sets of similar bers.

Algebraic points set surfaces (APSS) methods are an interesting approach to building these meta-bers, thanks to the projection operator they provide. Consequently, we will improve these APSS methods in a generic way before applying it to tractograms to obtain a new geometric representation for meta-bers.

I.3 Objectives and framework

This research aims at providing some useful methods and insights for research as well as clinical usages. The objective is to improve the usability of tractograms, in terms of manipulation and understanding. Indeed, simplifying bers by grouping them into meta-bers with encompassing and adapted geometry can help both for manipulating the data and making out the important tracts. This is further improved with the segmentation approach we introduce as it allows for the selection of a meaningful subset of the tractograms. Compressing the data also goes in the same direction, as more tracts could be retrieved from dMRI for the same disk usage. Last but not least, visualization improvements present valuable bene ts for an easier understanding of the data.

The scope of these methods will be their use on a single computer, let it be a powerful laptop or desktop, but not on a remote server. This is needed as we would like to have interactive usage with visualization of the geometry created as well as the originally computed bers. It means that the codes and the methods proposed will rely on this kind of hardware and not on servers or supercomputers. As a result, a clinician might be able to try and use the results exposed here without the need for more than a personal computer.

We aim at obtaining interactive approaches for simpli cation, segmentation, or visualization. This requires fast algorithms with the use of e cient languages like C++, CUDA or OpenGL. For the user interface, the choice has been made to use Qt for the clean appearance it provides as well as its powerful functionalities, allowing for OpenGL integration, dialogs for various input-outputs, and all its numerous widgets. Some pre-computation will be necessary, but should ideally be the shortest possible, as it will be discussed in the following sections.

I.4 Organization and contributions

We organize this thesis in four parts, starting with a state of the art of previous approaches, both in the medical eld and the computer graphics eld (Chapter II). In Chapter III, we discuss the data we are using -tractograms -, how they are obtained, their usage, and their limitations. Then, we introduce a new multiresolution representation for tractograms, whose objective is to ease visualization. This algorithm requires some speci c geometric representation for tting sets of bers. Chapter IV presents di erent possible geometries that can be used for this purpose, with their advantages and drawbacks depending on the targeted use. This chapter also details our new and fast algorithm for Moving Least Squares (MLS) projections using Algebraic Point Set Surfaces (APSS). Chapter V exposes various contributions based on the previously presented algorithms. First, a segmentation method is presented, followed by a compression algorithm for tractograms, that is then used for visualization of massive sets of data. Finally, we conclude and present some perspectives for other research opened by this work.

Overall, we propose algorithms and methods for a multi-resolution representation, as well as a segmentation and a compression all suited to tractograms. We also introduce an e cient approach to reduce APSS complexity in its global form. The work presented in this thesis has led to the following publications: The code will be publicly available.

All algorithms are illustrated using the Human Connectome Project (HCP) [Van Essen et al., 2012] dataset for the parts focused on tractography. For the geometry part, we use some 3D point sets from [START_REF] Ohtake | Multi-level partition of unity implicits[END_REF] as well as from the EPFL Computer Graphics and Geometry Laboratory 1 , the Stanford Computer Graphics Laboratory 2 , and http://threedscans.com.

II

State of the art

This chapter aims at presenting the state of the art of the di erent methods upon which this thesis is built. These approaches are also our starting and comparison points.

II.1 Tractography from Di usion Magnetic Resonance Imaging

Di usion Magnetic Resonance Imaging (dMRI) is the only medical imaging technique able to non-invasively explore the white matter architecture of the brain. This method is able to measure the di usion of water molecules inside the brain. The brain itself is composed of gray matter and white matter, with the interface between the two being called the cortical surface. Functional regions of the brain are localized in the gray matter, which includes the actual neurons. White matter is mainly composed of nerve connections (the axons of the neurons), resulting in a highly anisotropic medium for water molecules. Indeed, the water is blocked by the axons, and cannot di use past them1 . These axons are the physical wires of the brain and transfer electrical signals between di erent brain regions. As such, zones in which water di uses anisotropically characterize places where a set of axons all have the same direction. As in the gray matter there are no such structural elements as the axons, the di usion measured is usually more isotropic and does not give a lot of information.

This di usion measuring technique is used on several consecutive planes, in order to ll information in a grid of voxels. Di usion is measured in each voxel, with the rate as well as the main directions of propagation, and is then recorded using one of several representations [START_REF] Basser | Mr di usion tensor spectroscopy and imaging[END_REF], Tuch et al., 2002, Tuch, 2004]. Tensors are a common representation, even though "crossing bers" are ill-represented with this approach [START_REF] Westin | Image processing for di usion tensor magnetic resonance imaging[END_REF], Tournier et al., 2011, Soares et al., 2013].

The images obtained with the use of tensors are called Di usion Tensor Imaging (DTI). A single value is often used to characterize the isotropy per voxel. It is called the Fractional Anisotropy (FA) [START_REF] Westin | Image processing for di usion tensor magnetic resonance imaging[END_REF].

Using these di usion images, retrieving the brain connections inside the white matter can be performed using tractography algorithms such as FACT ( ber assignment by continuous tracking)-algorithms [START_REF] Mori | Three-dimensional tracking of axonal projections in the brain by magnetic resonance imaging[END_REF], Mori and Van Zijl, 2002, Tournier et al., 2011]. More details on this process are provided in Section III.1.1. Brain connections that are obtained are usually called bers or streamlines, and constitute a tractogram. As a result, tractograms obtained from dMRI are bundles of 3D polylines, which are estimates of the trajectories of large groups of neural tracts. Tractography has been proven to be an invaluable tool for clinicians and researchers. It is nowadays used on a daily basis by neurosurgeons for preoperating planning and during surgical operations [START_REF] Jeurissen | Di usion MRI ber tractography of the brain[END_REF]. It also o ers important information for studying pathological processes in neurological diseases [START_REF] Ciccarelli | Di usion-based tractography in neurological disorders: concepts, applications, and future developments[END_REF].

The quality and accuracy of the bers depend on the precision of the dMRI used. Typically, 3 Tesla (3T) dMRI scanners can reach voxel sizes of around 1

× 1 × 1
, and this can improve with higher magnetic elds as 7 Tesla scanners or even more are emerging, allowing for smaller voxels and better accuracy. Figure II.1 -Visualization of images from MRI, Fractional Anisotropy and a tractogram using 3DSlicer [START_REF] Pieper | The na-mic kit: Itk, vtk, pipelines, grids and 3d slicer as an open platform for the medical image computing community[END_REF] and TrackVis [START_REF] Wang | Di usion toolkit: a software package for di usion imaging data processing and tractography[END_REF].

II.2 Simpli cation methods

Recent tractography methods produce up to several millions of bers [START_REF] Tournier | Di usion tensor imaging and beyond[END_REF]. This can complicate the rendering, visualization, and interpretation of tractograms, thus limiting the aforementioned clinical applications. Furthermore, the considerable number of bers can make computationally intractable processes such as non-linear registration or atlas construction [START_REF] Gori | Parsimonious approximation of streamline trajectories in white matter ber bundles[END_REF], which are important for research purposes. Many bers might have a similar trajectory and connectivity, making the tractogram redundant. For this reason, several authors have proposed new geometric representations and visualization techniques to simplify tractograms. One of the most popular approaches consists in grouping similar bers into clusters [START_REF] Zhang | Hierarchical Clustering of Streamtubes[END_REF], Maddah et al., 2007, Wassermann et al., 2010, Guevara et al., 2011, Garyfallidis et al., 2012, Demir and Cetingül, 2015, Gori et al., 2016, Alexandroni et al., 2017, Zhang et al., 2018, Siless et al., 2018] in order to remove statistical redundancy. These clusters are then approximated by one representative ber, usually called a prototype [START_REF] Guevara | Robust clustering of massive tractography datasets[END_REF], or with a new ber created for this purpose [START_REF] Garyfallidis | QuickBundles, a Method for Tractography Simplication[END_REF]. They are built using metrics speci cally designed for brain streamlines, thus preserving important properties such as extremities and path [START_REF] Olivetti | Comparison of distances for supervised segmentation of white matter tractography[END_REF], Siless et al., 2018]. Other authors have also proposed to represent the spatial extent of the clusters using an encompassing geometry [START_REF] Maddah | Probabilistic clustering and quantitative analysis of white matter ber tracts[END_REF]. These methods are usually controlled by one parameter, e.g. a threshold [START_REF] Zhang | An anatomically curated ber clustering white matter atlas for consistent white matter tract parcellation across the lifespan[END_REF], thus presenting only one level of resolution at a time. Furthermore, information such as the number of bers or the spatial extent (i.e. the volume) of the cluster might be lost in the process.

Progressive simpli cation methods, which focus on decimating complex 3D objects while preserving both important geometric features and topological relationships, have so far not been explored to approximate large scale tractograms. Zhang and Laidlaw [START_REF] Zhang | Hierarchical Clustering of Streamtubes[END_REF] propose to use a hierarchical clustering algorithm to progressively group together similar bers. They only applied it to bundles composed of a small number of bers.

II.3 Geometric methods

In this thesis, we study di erent geometries to represent bers and sets of bers.

Clusters methods can use prototypes or single bers to represent sets of simpli ed bers [START_REF] Guevara | Robust clustering of massive tractography datasets[END_REF], Garyfallidis et al., 2012]. A geometric representation of grouped bers can also be used [START_REF] Maddah | Probabilistic clustering and quantitative analysis of white matter ber tracts[END_REF], such as generalized cylinders [START_REF] Petrovic | Visualizing whole-brain dti tractography with gpu-based tuboids and lod management[END_REF] so that the spatial extent of the merged bers are kept at each resolution. To improve the visualization quality and e ciency, the geometric models used to represent the bers are often computed directly on the GPU [START_REF] Enders | Visualization of white matter tracts with wrapped streamlines[END_REF], Reina et al., 2006, Petrovic et al., 2007]. Using these approaches, the data will be easier to understand and visualize, or even to analyze according to some given criteria (general shapes, connections, etc.).

However, the spatial extent of groups of bers is usually lost in the process as the geometry used is not able to accurately represent all the information. For this reason, we focus in this thesis on point set surfaces and more speci cally on its algebraic version.

Point Set Surfaces Approaches constructing surfaces from point sets with normals can also be an interesting representation for bers, as it will be further explained in Section IV.3.

There has been an extremely large body of works for de ning and extracting smooth surfaces from point samples, as well as for processing and ltering point sets. We will only cover the most relevant works to our technique presented in Section IV.2, and we refer the interested reader to the recent survey of Berger et al. [START_REF] Berger | A survey of surface reconstruction from point clouds[END_REF] as well as to a previously published benchmark [START_REF] Berger | A benchmark for surface reconstruction[END_REF].

Moving Least Square surface approximation is a tool often used for its interesting properties [Levin, 1998]. It allows for local projections on point set surfaces [START_REF] Alexa | Point set surfaces[END_REF], Levin, 2004, Amenta and Kil, 2004] and presents some good denoising capabilities [Levin, 2004, Reuter et al., 2005]. On top of that, it is a fast method, which is useful for fast reconstructions or rendering [START_REF] Alexa | Point set surfaces[END_REF], Guennebaud et al., 2008]. It has been used for projection on di erent geometric primitives, such as hyperplanes [Levin, 2004], spheres [Guennebaud andGross, 2007, Guennebaud et al., 2008], Hermite data [Alexa and Adamson, 2009], or parabolic-cylinders [START_REF] Ridel | Parabolic-cylindrical moving least squares surfaces[END_REF]. Each of these projection versions has its own advantages and drawbacks depending on the point set to reconstruct. One common issue that MLS projections face is the reproduction of sharp features, because of their tendency to smooth the surface locally (which is why it is great for denoising). Some articles focus on these sharp features in order to preserve them while still removing the noise on other parts of the models [START_REF] Reuter | Point set surfaces with sharp features[END_REF], Öztireli et al., 2009, Guillemot et al., 2012]. Improving the quality of the resulting MLS surface was also done by using local features [START_REF] Dey | An adaptive mls surface for reconstruction with guarantees[END_REF] or gradient elds [START_REF] Chen | Non-oriented MLS Gradient Fields[END_REF]. The MLS framework is quite versatile and has also been used for other purposes, such as the introduction of opening and closing operators on 3D point sets [START_REF] Calderon | Point morphology[END_REF]. In this thesis, we will focus on algebraic point set surfaces (APSS) [Guennebaud andGross, 2007, Guennebaud et al., 2008], which consists of iteratively nding a local sphere tting the model before projecting on it. Obtaining this sphere is quite slow when considering all points from the input point set. Consequently, only the k-nearest neighbors are often used to reduce complexity.

Multi-level Partition of Unity Given an input point set with oriented normals, the Multi-level Partition of Unity (MPU) [START_REF] Ohtake | Multi-level partition of unity implicits[END_REF] aims at reconstructing an implicit surface from locally estimated quadratic functions, using an octreebased partitioning. During the top-down octree traversal, the points inside a cell are approximated by a quadric. If the residual error is above a xed tolerance, the cell is subdivided and the same process is done on the children cells until the tolerance is satis ed. This results in a set of locally-de ned primitives, and the global implicit function is obtained through a smooth blending using partitions of unity.

In MPU, as in point set surfaces methods, primitives are locally estimated from the input point set. However, opposite to point set surfaces, the primitive functions are attached to the cells and not to individual query points. As a consequence, there is no simple projection operator onto the reconstructed surface and further polygonalization has to be done for a visualization purpose. Moreover, the intrinsic local aspect of MPU makes it very sensitive to the quality of the input set, especially on the input normal vectors.

While it is possible to in uence the output by tuning the tolerance parameter, this technique is not appropriate for ltering: a higher tolerance value will result in a less detailed output but not in a smoothly-varying surface (i.e., which would be the result of a low-pass lter).

Xiao [Xiao, 2011] proposed to combine MPU with APSS by computing an implicit approximation of the input sample set before de ning the projective and implicit evaluations based on this approximation.

Global Surface Reconstruction Several methods exist that aim at reconstructing a surface globally from an input point set, with or without normals. The state of the art in terms of quality of reconstruction is the Poisson reconstruction [START_REF] Kazhdan | Poisson surface reconstruction[END_REF], which interprets the input normals as the gradient of the solid's indicator function. It then creates a vector eld in 3D and solves a Poisson equation to obtain the scalar function, allowing for the extraction of the surface. This method was further re ned with the addition of a screening term, reducing the smoothing of the nal surface that could occur previously [START_REF] Kazhdan | Screened poisson surface reconstruction[END_REF]. Another competing method, with results of equivalent quality, is [START_REF] Calakli | Ssd: Smooth signed distance surface reconstruction[END_REF]. Their approach is similar to the Poisson surface reconstruction, but instead of forcing the implicit function to approximate the indicator function of the volume, they force the implicit function to be a smooth approximation of the signed distance function to the surface. This removes some limitations that the Poisson reconstruction has. More recently, another approach using Deep Learning was introduced by [START_REF] Hanocka | Point2mesh: A self-prior for deformable meshes[END_REF]. They use self-prior computed using a neural network to iteratively shrink-wrap an initial mesh which is the convex hull. Other approaches are interesting, such as [START_REF] Mullen | Signing the unsigned: Robust surface reconstruction from raw pointsets[END_REF] that does not require normals to perform its reconstruction, or more recently [START_REF] Barill | Fast winding numbers for soups and clouds[END_REF].

These methods perform well on point sets with evenly distributed points. However, they do not always ll-in holes correctly, and they have no parameters for controlling the output surface, let it be for ltering, or adapting the tting surface locally.

Variational Implicit Point Set Surfaces Huang et al. [START_REF] Huang | Variational implicit point set surfaces[END_REF] de ne a variational problem to reconstruct a global implicit surface from an input point set. By minimizing a compromise between the Duchon's energy [Duchon, 1977] and the data tting energy, the normal vectors and the implicit function are jointly estimated. While this global approach allows for a robust reconstruction and does not require input normals, it requires inverting a fully-dense matrix as well as estimating its eigenvalues, which makes the complexity cubic with respect to the number of input points. This work is highly related to ours, as the regularization energy enforces closeto-triharmonic surfaces, which results in low variations of the curvature in the output, similar to what is achieved by the APSS approach [Guennebaud andGross, 2007, Guennebaud et al., 2008] when tting smoothly-varying spheres to the input point set. This follows most recent hole-lling (and alike) techniques [START_REF] Botsch | An intuitive framework for real-time freeform modeling[END_REF], Jacobson et al., 2010, Lescoat et al., 2019], which advocate the use of triharmonic functions for smooth surface modeling.

Similarly to other global approaches that we have previously mentioned, while it is possible to in uence the output by tuning the constant balancing the regularization and tting energies, this technique is not adequate for ltering input point sets.

Implicit Surfaces In this thesis, implicit surfaces are used as a visualization tool. They correspond to the isovalue of a scalar function F(x, y, z). Modeling approaches were initially using points, called the skeleton, to parameterize a scalar eld such as the blobs, metaball or soft objects [Blinn, 1982, Nishimura, 1985, Wyvill et al., 1986]. The nal function F being computed as a summation of all scalar elds related to all points. Extending this idea to curved skeletons can be achieved through the use of convolution surfaces [START_REF] Bloomenthal | Convolution surfaces[END_REF], where the eld is obtained as the continuous integration of a smooth kernel function along the skeleton. The rendering of the iso-surface can be obtained from the generation of an approximating triangular mesh [Lorensen andCline, 1987, Ju et al., 2002], or the use of e cient ray-casting able to evaluate the closest distance to the surface [START_REF] Mccormack | Creating and rendering convolution surfaces[END_REF]] that ts GPU capabilities.

II.4 Segmentation methods

It is of interest for clinicians and researchers to isolate precise white matter tracts of the brain for surgical planning or to test hypotheses related to a certain condition. Segmenting white matter tractograms into anatomically relevant and reproducible tracts is a di cult task due to the huge number of bers -up to millions -the presence of spurious bers and the high inter-subject variability. Furthermore, the anatomical de nition of a tract is usually qualitative and vague since it is based on relative spatial relationships (e.g. "anterior to", "close to") with respect to other anatomical structures.

The most common technique for identifying a tract is to manually delineate Regions of Interest (ROIs) and select (or exclude) the bers that pass through them [START_REF] Wakana | Reproducibility of quantitative tractography methods applied to cerebral white matter[END_REF]. This method is tedious, time-consuming, it requires extensive training and it is not very reproducible for tracts with convoluted trajectories [START_REF] Zhang | Atlas-guided tract reconstruction for automated and comprehensive examination of the white matter anatomy[END_REF]. Another technique is based on the transfer of manually segmented ROIs from one (or multiple) training images to test subjects via non-linear deformations [START_REF] Zhang | Atlas-guided tract reconstruction for automated and comprehensive examination of the white matter anatomy[END_REF]. The resulting segmentation highly depends on the quality of the registration which might not be accurate when training and test images do not share the same anatomical topology (e.g. due to a tumor or illness). A third class of methods is based on clustering algorithms where clusters are usually classi ed into anatomical tracts using a labeled atlas or training subjects [O' Donnell and Westin, 2007, Guevara et al., 2012, Prasad et al., 2014, Garyfallidis et al., 2018]. These methods are usually computationally demanding, depend on user-tuned parameters (e.g. the size or the number of clusters) and the atlases are usually estimated from a small population of healthy subjects. With a di erent perspective, in [START_REF] Wassermann | The white matter query language: a novel approach for describing human white matter anatomy[END_REF] the authors proposed a query language (WMQL) to interactively de ne and segment anatomical tracts employing simple spatial relationships and logical operations. This method is fast and easy to use but the shape of the resulting tracts can vary among subjects since WMQL is based on binary relations and parallelepiped bounding boxes.

II.5 Compression methods

The several millions of bers that can be built using tractography algorithms can result in tens of gigabytes (GB) of data [START_REF] Rheault | Visualization, Interaction and Tractometry: Dealing with Millions of Streamlines from Di usion MRI Tractography[END_REF]. For instance, a le containing 1 million bers and obtained with a constant stepsize of 0.1 can weigh up to 8.7 . This massive amount of data complicates visualization, processing, sharing, or storage.

Existing methods, that propose a solution to this data size problem, can be divided into two di erent categories. They either compress the whole tractogram or the representation of each ber.

Tractogram-level compression

This kind of approach consists in reducing the number of bers of a tractogram. A simple way of doing so is to randomly select a subset of bers from the original dataset [START_REF] Gori | Parsimonious approximation of streamline trajectories in white matter ber bundles[END_REF]. This is a pragmatic approach, however, there is no guarantee that meaningful bers are not removed. Other approaches correspond to tractogram simpli cation (Section II.2), including cluster methods and progressive groupings. Nevertheless, depending on the application, keeping all the original bers might be required or preferred.

Fiber-level compression When trying to keep as much information as possible from the original tractogram, per-ber compression algorithms appear to be a good alternative. These methods can, in most cases, be combined with a tractogram-level compression.

It is possible to use general compression algorithms directly on the data. Whereas lossless compression algorithms compress data by identifying and removing statistical redundancy -and therefore do not lose any information from the original data -lossy algorithms are going further by removing unnecessary data or introducing approximations. In the case of brain tractograms, the data are composed of a set of 3D points expressed as oating-point values. A pragmatic approach is to compress them using a general oating-point data compression algorithm such as the one of [Lindstrom, 2014]. However, this approach cannot compete with more data-aware methods, which are more e cient as they take advantage of the speci cities of the data. For instance, a classical approach in neuroscience is to reduce the number of points, by using algorithms such as a linearization that removes collinear (or almost collinear) points [START_REF] Presseau | A new compression format for ber tracking datasets[END_REF]. This can be combined with a general compression algorithm as data are still oating-point values, improving the compression ratio.

E cient compression algorithms are data-aware, which means that they take advantage of the speci cities of the data, such as their representation, meaning, or usage.

One possible approach for a data-aware compression is to build a dictionary in a dedicated space of representation [START_REF] Presseau | A new compression format for ber tracking datasets[END_REF], Kumar and Desrosiers, 2016, Moreno et al., 2017]. Presseau's algorithm [START_REF] Presseau | A new compression format for ber tracking datasets[END_REF] uses a three steps process, starting with (i) linearization, then (ii) quantization, and eventually (iii) encoding the points using a dictionary. The method is similar to the zip compression, and the algorithm is named zfib. This method achieves good compression ratios (around 90%) but at the cost of high compression times (more than ten minutes for tractograms with 1M bers). Some points from the original bers are removed due to the linearization step, therefore reducing the performance of point-based algorithms [START_REF] Soares | A hitchhiker's guide to di usion tensor imaging[END_REF]. The authors also tried to apply some transformations -including wavelet-based ones -without success as they reduced the compression ratio.

Another compression algorithm method dedicated to tractography uses a sparse representation of bers [START_REF] Chung | E cient parametric encoding scheme for white matter ber bundles[END_REF], Kumar and Desrosiers, 2016, Moreno et al., 2017]. In [START_REF] Chung | E cient parametric encoding scheme for white matter ber bundles[END_REF], the bers are represented using cosine functions, parametrized with 60 parameters, reducing the memory cost. The average error introduced in this case is about 0.26

. [START_REF] Kumar | A sparse coding approach for the e cient representation and segmentation of white matter bers[END_REF]'s approach segments and clusters bers using a dictionary built with a sparse coding of the bers. In that case, the goal is not to reduce the memory usage but to ease some heavy computations on bers. In [START_REF] Alexandroni | The ber-density-coreset for redundancy reduction in huge ber-sets[END_REF] the idea is to use a dictionary combined with the ber-density-coreset method. As a result, it removes high-frequency data (sudden changes in the bers path are smoothed) by only retaining a few non-zero coe cients for each ber. The typical average error is around 2 on data with a resolution of 1.25 3 . In [START_REF] Caiafa | Multidimensional encoding of brain connectomes[END_REF], a sparse encoding is performed on the tensor representation. This results in a size in memory which is about the same as the one with the usual 3D points representation. Moreover, points are moved at the center of the voxels of the original grid, resulting in possible important errors, and point order is lost in the process.

Considering that current tractography methods are usually applied on di usion images with a resolution of 1 3 , the maximum error when compressing bers should not exceed the voxel size, and should ideally be an order of magnitude smaller. This should be the case not only with clinical DWI, whose voxels size is around 2

3 , but also for research protocols (e.g. 1.25 3 in [Van Essen et al., 2012]). Moreover, these dictionary-based methods do not allow for the compression and decompression of speci c fascicles or tracts, which could drastically reduce the memory load for visualization or processing.

II.6 Visualization techniques

The 3D visualization of tractograms is necessary for their study, no matter the context. Consequently, some tools exist and are widely used by surgeons and searchers [START_REF] Soares | A hitchhiker's guide to di usion tensor imaging[END_REF]. This tools include VTK [START_REF] Schroeder | The visualization toolkit: an object-oriented approach to 3D graphics[END_REF], paraview [START_REF] Ahrens | Paraview: An end-user tool for large data visualization[END_REF], 3DSlicer [START_REF] Pieper | The na-mic kit: Itk, vtk, pipelines, grids and 3d slicer as an open platform for the medical image computing community[END_REF], trackvis [START_REF] Wang | Di usion toolkit: a software package for di usion imaging data processing and tractography[END_REF], Fiberweb [START_REF] Ledoux | Fiberweb: di usion visualization and processing in the browser[END_REF] , FiberVis [START_REF] Osorio | Fibervis: a tool for a fast ber tractography visualization and segmentation[END_REF] and many others.

Visualizing the tractograms, with also the possibility to visualize other elements including the original MRI images is really important, however, it does not improve the visualization itself. Plenty of possibilities have been tried to either improve visualization quality or even add some more information onto the tractograms [Isenberg, 2015]. For instance, the authors in [START_REF] Everts | Exploration of the Brain's White Matter Structure through Visual Abstraction and Multi-Scale Local Fiber Tract Contraction[END_REF] try to ease visualization by bringing similar bers closer to each other, e ectively moving the bers. This approach is controversial as it does not improve visualization so much, at the expense of destroying some pieces of information. To improve readability, some labels with functional information can be added on the bers geometric representation [START_REF] Petrovic | Visualizing whole-brain dti tractography with gpu-based tuboids and lod management[END_REF]. This is done using some speci c shading and colors, allowing for easy visual separation of bers. Some geometric representations are also used to improve visualization, and mainly reduce the number of individual bers to display, as every single one of them might not be relevant [START_REF] Enders | Visualization of white matter tracts with wrapped streamlines[END_REF], Reina et al., 2006, Petrovic et al., 2007]. To display additional information, some colors, shapes or even distortions [START_REF] Reina | GPU-based hyperstreamlines for di usion tensor imaging[END_REF]] can be used. Rendering 3D lines while retaining a depth perception is achieved with geometric representations, however, these techniques can be heavy for large datasets and do not always help visualize tight line bundles. To this end, some methods were proposed to render dense line data by adding some lighting shading on lines [START_REF] Mallo | Illuminated lines revisited[END_REF], adding a halo e ect around bers depending on their respective depth [START_REF] Everts | Depth-dependent halos: Illustrative rendering of dense line data[END_REF] or add some ambient occlusion speci cally to lines [START_REF] Eichelbaum | Lineao-improved three-dimensional line rendering[END_REF].

All these techniques improve visualization over a basic display of the lines. However, they are not working in all conditions, especially with lots of lines, and might su er from performance issues.

On top of that, some other techniques from computer graphics could be extended to visualizing tractograms as sets of curves or meshes, with for instance hair rendering techniques [START_REF] Ward | A survey on hair modeling: Styling, simulation, and rendering[END_REF] or vector visualization like is encountered on road maps [START_REF] Amiraghdam | Localis: Locally-adaptive line simpli cation for gpu-based geographic vector data visualization[END_REF].

II.7 Motivations

There is a need for tractogram representations that would be able to simplify them while keeping their important properties: spatial extent, connectivity, and functions. To achieve this, it is necessary to investigate ways of grouping lines, as well as geometric representations that are adapted to the data speci cities. One possible approach is to use point sets de ning sets of bers. However, this requires being able to ll in potential gaps while controlling the nal surface for good spatial extent conservation. Then, with the help of this representation, we also want to explore ways of segmenting known bundles. This needs to be reproducible and should also take the natural relation de nitions as input to improve its usability. Tractogram data being redundant, it should also be feasible to reduce the necessary storage by considering the bers construction process. The visualization of the tracts is an important usage of the data. As a result, some visualization techniques are worth inquiring about.

III

Tractography and simplification

In this chapter, we present tractography, its limitations and usage, as well as their important properties. Then, we introduce a new multiresolution approach for simplifying bers datasets.

III.1 Tractography

III.1.1 Creation of the tractogram

Starting from the di usion image of the dMRI, several algorithms exist to reconstruct the sets of neural connections. These algorithms are either probabilistic or deterministic. They both build connections as 3D lines composed of points joint together. These lines are called bers, and each of these represents a large number of actual connections. Indeed, if the voxel size of the image from the scanner is 1 × 1 × 1 , the diameter of an axon is around a micrometer ( ) with length starting at a few centimeters and reaching about twenty centimeters. In this thesis, Mrtrix3 was used to obtain these bers, using di erent algorithms among those proposed.

Most of the time, the di usion sequence is acquired on top of an anatomical MRI. The white matter is obtained directly from the anatomical MRI, as well as the cortical surface. To create a possible ber, FACT ( ber assignment by continuous tracking)-algorithms [START_REF] Mori | Three-dimensional tracking of axonal projections in the brain by magnetic resonance imaging[END_REF], Mori and Van Zijl, 2002, Tournier et al., 2011] start by randomly placing a seed inside the white matter (Figure III.1a). Then, following the main direction dictated by the current voxel and its neighbors (each voxel contains a main direction of di usion, which corresponds to the greatest eigenvector for DTI), a new point is created at either a xed distance or a varying distance depending on various parameters such as the local anisotropy for instance (Figure III.1b). This process is iterated until the cortical surface is reached (Figure III.1c). Starting from the seed, the same process is repeated in the opposite direction (Figure III.1d). During this step-by-step process, the angle between consecutive segments has to be less than a predetermined threshold, usually chosen by the user (Figure III.1c). If this threshold is exceeded, the ber is considered invalid and discarded. The algorithms continue creating bers until as much valid bers as asked are obtained. The two crossed out segments in Figure III.1c are presented for illustration purposes and cannot occur at the same time, as only one possible segment is evaluated at each step. Some information represented by scalar values can be kept at each point of the bers, such as the functional anisotropy (FA) for instance [START_REF] Westin | Image processing for di usion tensor magnetic resonance imaging[END_REF], Brun et al., 2003, Yeatman et al., 2012]. Tractography algorithms can be either deterministic (as the one described above) or probabilistic. Probabilistic algorithms are based on the same principle. However, the next direction during the ber construction is chosen probabilistically with respect to the local tensors. This means that starting from the exact same seed, di erent bers can be obtained. Although it tends to obtain too much false positive and appears to be less precise than deterministic methods [START_REF] Sarwar | Mapping connectomes with di usion mri: deterministic or probabilistic tractography?[END_REF], it is better for some speci c reconstructions [START_REF] Schlaier | Probabilistic vs. deterministic ber tracking and the in uence of di erent seed regions to delineate cerebellar-thalamic bers in deep brain stimulation[END_REF]. The FACT-algorithms [START_REF] Mori | Three-dimensional tracking of axonal projections in the brain by magnetic resonance imaging[END_REF], Mori and Van Zijl, 2002, Tournier et al., 2011] are used for both deterministic and probabilistic approaches and are still the basic principle for various recent tracking algorithms such as iFOD1, SD_STREAM [START_REF] Tournier | MRtrix: Di usion tractography in crossing ber regions[END_REF] or iFOD2 [START_REF] Tournier | Improved probabilistic streamlines tractography by 2nd order integration over bre orientation distributions[END_REF]] that we are using in this thesis.

Usually, with the classical di usion images, a million bers is used, however, nothing prevents from obtaining more of them. The set of bers obtained is called a tractogram, and the process to obtain it is the tractography.

III.1.2 Tractogram usage

Tractography is guided by di erent objectives. The rst is a clinical one, as knowing the position of the important bundles (set of bers) is helpful for tumor ablation, for pre-operation planning, or even during the operation [START_REF] Smits | Incorporating functional mr imaging into di usion tensor tractography in the preoperative assessment of the corticospinal tract in patients with brain tumors[END_REF], Yamada et al., 2009, Soares et al., 2013]. The surgeon can use the tractogram to determine the best path to remove the tumor, avoiding as much as possible the connections considered the most important. The second objective is research, as the knowledge on these connections helps understand the way the brain operates and the anatomical changes induced by some pathologies.

III.1.3 Tractogram validity

Current tractography algorithms are able to retrieve almost every existing connection in the white matter at the available resolution. This was veri ed in two di erent manners, rst by realizing the whole tractography process on the brain of a dead person before dissecting it for comparison [START_REF] Bürgel | White matter ber tracts of the human brain: threedimensional mapping at microscopic resolution, topography and intersubject variability[END_REF], and secondly using phantoms (a set of synthetic bers) and by realizing some tractography simulation on these data [START_REF] Poupon | New di usion phantoms dedicated to the study and vali-dation of high-angular-resolution di usion imaging (hardi) models[END_REF], Fillard et al., 2011]. This is extremely positive and promising, however, a downside is that some constructed bers do not exist in the original brain. About two-thirds of the bers are false positive [START_REF] Maier-Hein | The challenge of mapping the human connectome based on di usion tractography[END_REF]. This is intrinsically due to the relatively high size of the voxels compared to the size of the actual axons, resulting in an approximation of the directions that can result in false bers being obtained. Indeed, when bers are crossing, there are a few possible reconstructions, some of them being incorrect (see Figure III.2). Some techniques exist to try to remove these non-existent bers, however, they cannot remove every false positive.

III.1.4 Important properties of brain bers

To be able to simplify tractograms, it is important to rst apprehend what are the important properties of brain bers that should be preserved. Fibers are a representation of a possible path for electrical neural signals. This means that the function of each ber is important. This function depends on the regions the ber connect and can be retrieved using its extremities. As such, the endpoints of each ber should be considered preferentially when comparing or simplifying bers. The path of each ber is also important as it allows knowing its position inside the brain, which is helpful for surgery for instance. At the level of the tractogram, this means that connections from one functional region to another need to be kept as well as the spatial extent of these sets of connections. However, the individual bers obtained through tractography algorithms may not all be important. They can be redundant by construction as some close-by seeds might often lead to similar bers. It makes simpli cation relevant as it is in this case possible to summarize the important information providing we have a way to extract it.

III.2 Multi-resolution III.2.1 Introduction

Current tractography methods generate tractograms composed of millions of 3D polylines, making visualization and interpretation extremely challenging, thus complexifying the use of this technique in a clinical environment. Consequently, simplifying tractograms is often necessary. Usual methods consist in choosing some representative bers [START_REF] Garyfallidis | QuickBundles, a Method for Tractography Simplication[END_REF]. Some parameters will determine the nal number of bers obtained, meaning that a recomputation is needed if the result is too scattered or too dense, which may happen quite often as they are not always intuitive. We do not focus here on single resolution approximations but propose a multiresolution representation that allows the users to navigate among the resolutions depending on their objective.

Taking inspiration from error-driven surface mesh simpli cation [START_REF] Garland | Surface Simpli cation Using Quadric Error Metrics[END_REF], we introduce a progressive merging strategy for grouping bers into meta-bers. The proposed method reduces the redundancy of the tractogram, producing a multi-resolution structure, which is organized into a nested hierarchy of levels of detail. Every fusion of bers represents a new level of resolution. Once the entire multi-resolution representation is computed, the user can navigate through di erent levels of detail in a continuous fashion and in real-time, while maintaining the overall structure of the original tractogram (see Figure III.3). Furthermore, we also propose an e cient implementation based on a Delaunay tetrahedralization which makes our method usable even on large tractograms containing millions of bers. We illustrate our method using the HCP dataset [Van Essen et al., 2012].

Figure III.3 -Thalamocortical bundle presented with three di erent resolutions. 50% and 25% resolutions are obtained using our multi-resolution. The color code depends on the local orientation of the ber: red for left-right, blue for inferior-superior, and green for anteroposterior.

This work was published as a conference paper in VCBM, a Eurographics workshop [START_REF] Mercier | Progressive and E cient Multi-Resolution Representations for Brain Tractograms[END_REF].

III.2.2 Method

III.2.2.1 Fiber grouping

We propose a tractogram simpli cation method based on the progressive mesh methodology [Garland andHeckbert, 1997, Hoppe et al., 1993]. Given a tractogram with bers, we rst look for the two most similar bers based on a similarity measure. Once detected, the couple is collapsed into a new meta-ber. A metaber is the name we will give to a geometric structure representing a set of bers. It is to be noted that all input bers can be seen as a meta-ber, which means that the geometry needs to be de ned in this speci c case. The similarity of this new meta-ber is compared against its neighbors, which are the union of the neighbors of the underlying bers. The process is then iterated until obtaining a single large meta-ber (i.e. -1 iterations) or using a stopping criterion to prevent over-simpli cation (see Section III.2.3.3).

The proposed "progressive brain tractograms" algorithm is general and can be used with any similarity measure. It is inherently multi-resolution, where every level of detail corresponds to the fusion of two meta-bers. Once the sequence of fusions is computed, we can visualize the tractogram at any resolution and switch among levels in real-time.

Fibers are considered similar when their pathway as well as their ending points are similar. In this work, we do not consider the possible functional signals describing the brain activity that could be attached to them. Fibers are inherently oriented as their points are stored in a given order, with a rst and a last point. This has no anatomical meaning as they already represent a set of axons and not individual nerve connections. Therefore, the comparison tools need to be invariant with respect to the orientation. We want to be able to compare bers together and classify them from the most to the least similar. This can be achieved using dissimilarities derived from distances. In the following, we will use three standard similarity measures for bers comparisons. Let = {x i , = 1... } and = {y j , = 1... } be two bers composed of and points respectively. Their dissimilarity can be expressed as:

• Mean of Closest distances (MC) [START_REF] Guevara | Robust clustering of massive tractography datasets[END_REF].

( , ) = ( ( , ) + ( , ))/2 (III.1)
where ( , ) = 1 =1 min y j ∈ ||x iy j || 2 .

• Minimum average direct-ip (MDF) [START_REF] Garyfallidis | QuickBundles, a Method for Tractography Simplication[END_REF]. It is de ned in a similar way than MC, assuming bers contain the same number of points: = = (constant).

( ,

) = min( direct ( , ), ipped ( , )) (III.2) with direct ( , ) = 1 =1 ||x i -y i || 2 and ipped ( , ) = 1 =1 ||x i -y L-i+1 || 2 .
• Weighted Currents (WC) [START_REF] Gori | Parsimonious approximation of streamline trajectories in white matter ber bundles[END_REF]. The dissimilarity measure of WC is de ned as:

( , ) = ( f a -t a 2 ) ( f b -t b 2 ) -1 =1 -1 =1 i ( c i -d j 2 ) j (III.3)
where c i and i (respectively d j and j ) are the centers and tangent vectors of (respectively ), f a , f b and t a , t b are the corresponding endpoints of and respectively, and , and are three Gaussian kernels parametrized by , and respectively.

The MDF constraint of having the same number of points is satis ed by resampling the bers. We use 100 points per ber in our implementation. It is also done for the other dissimilarities as it increases the computation speed by reducing the number of computations required for each comparison. The WC is one of the few distance measures that takes the extremities into account, thus, it should yields better comparisons, especially with our neighborhood simpli cation that is presented next.

III.2.2.2 Delaunay tetrahedralization

Comparing all bers to each other leads to quadratic complexity and intractable computations. Furthermore, most of the computations would be useless since similar close-by bers should be merged in priority, thus indicating that only neighboring bers need to be compared and considered as candidates for merging. As introduced previously, one major element of the bers is their connectivity, i. e., where they begin and end, making bers endpoints really meaningful. We exploit this idea and this property by computing once and for all an adjacency relationship among bers, relying on the Delaunay tetrahedralization of their extremities. This tetrahedralization uses the Euclidean distance between all the endpoints of the bers and is computed using the fast implementation of TetGen [Si, 2015]. The adjacency links are given by the edges of the tetrahedral mesh. By exploiting this tetrahedralization, we drastically reduce the number of computations since we now compare every ber with only the adjacent (linked) ones. When merging two meta-bers, the adjacency links are updated accordingly: a neighborhood relationship is built between the new meta-ber and all the neighbors of its two original meta-bers, while previous neighborhood relationships involving these two original bers are suppressed. Relying on the Delaunay tetrahedralization provides a geometrically well-distributed set of adjacency links, which is bene cial to our simpli cation. Furthermore, results obtained with this technique are not biased, namely equivalent to the ones obtained without using it.

The entire pipeline is summarized in Figure III.4.

III.2.2.3 Geometric representation

We only consider so far a description of every ber based on its geometry and connectivity (we do not take into account any function or scalar along the bers). The chosen geometric representation should preserve bers properties. To this end, we propose to use generalized cylinders with an elliptical basis as a rst geometric representation for meta-bers. Generalized cylinders are objects obtained by sweeping a 2D contour along a 3D curve. In our case, the 3D curve is our polyline, and we de ne the 2D contour as ellipses. We use one ellipse per vertex of the centercurve such that the resulting cylinder incorporates the trajectory and endpoints of the meta-bers of the previous level of resolution.

Given two cylinders 1 and 2 , Figure III.5 schematically presents how they are merged into cylinder 3 . We rst nd for each extremity of 1 and 2 the closest point on the other cylinder. If the closest point P c is not an extremity itself, we consider all the points past P c until the actual endpoint to be part of the extremity region, whereas all other points are in the body region (Figure III.5a). For the body, we use as reference the center-curve (or ber) with less points ( 1 in Figure III.5). For each point x 1 of 1 , we look for the closest point x 2 in 2 . The corresponding point x 3 in 3 is computed as:

x 3 = x 1 1 +x 2 2
1 + 2 , where the weights 1 and 2 correspond to the number of original bers that 1 and 2 represent. We set will have circular pro les with a minimum radius to start with.

Body

It is important to note that we use the Euclidean distance between the endpoints of 1 and 2 to compute the extremities of 3 instead of a more accurate (but computationally expensive) geodesic distance. On top of that, to make sure that the extremities of 3 would actually lie at the border with the gray matter (the cortical surface), they should be projected onto it.

III.2.2.4 Implementation

We implemented our method in C++ using Qt for the graphical user interface and OpenGL for the rendering. CUDA is also used to improve the rst similarity computation between neighbors. Algorithm 1 presents the steps used for obtaining our multi-resolution representation. The code is publicly available on github1 .

Algorithm We do not use CUDA in the second step of this algorithm, as the parallelism is not important enough to bene t from a GPU version, that will eventually be slower than the multi-core CPU version.

To improve e ciency and memory usage while navigating between resolutions, the generalized cylinder geometry is computed on the GPU, using the hardware tessellation unit to synthesize our visual approximations on-the-y. This is performed by recording, for every center curve, the parameters required for an ellipse while computing the multi-resolution: the two axes as well as the orientation vector.

III.2.3 Results and comparison

III.2.3.1 Dataset

We conducted our experiments on the HCP dataset [Van Essen et al., 2012]. Tractograms are obtained using the SD_STREAM deterministic tracking algorithm of MRTrix3 [START_REF] Tournier | MRtrix: Di usion tractography in crossing ber regions[END_REF]. All tracts employed in the experiments (Uncinate Fasciculus (UF), Inferior Fronto-Occupital Fasciculus (IFOF) and Thalamocortical bundle) are extracted using either the white matter query language (WMQL) [START_REF] Wassermann | The white matter query language: a novel approach for describing human white matter anatomy[END_REF] or manual segmentation. Please note that our method is general and could be used with any streamline tractography algorithm. In Figure III.6, we compare our method with QuickBundles (QB) [START_REF] Garyfallidis | QuickBundles, a Method for Tractography Simplication[END_REF], a well-known approximation algorithm for white matter brain tractograms based on prototypes. We rst executed QB on a bundle composed of 19, 782 bers. We used di erent thresholds in the range of 5-10 as suggested in [START_REF] Garyfallidis | QuickBundles, a Method for Tractography Simplication[END_REF]. We chose the middle threshold of 7 for our comparison. It resulted in 275 prototypes. We then executed our algorithm (once) with the MDF metric, used by QB, and then interactively changed the resolution to obtain III.1), creates an encompassing representation that well approximates the original bundle at each level of resolution. By contrast, QB produces prototypes that, depending on a user-de ned threshold, might not preserve the overall structure and volume of the original bundle. Moreover, one might need to execute QB several times before nding an optimal threshold (for a given clinical/research application). The proposed method results in a multi-resolution representation which approximates the original tractogram with a decreasing precision and is not intended to produce anatomically reproducible clusters among subjects. These results suggest that the overall structure of the bundles is preserved across the multi-resolution representation, even at low resolutions, e.g. 15%. The graphs on the right of Figure III.7 illustrate the average and maximum valences across the resolutions. The valence of a meta-ber is de ned as its number of neighbors. During initialization, the number of neighbors is obtained from the tetrahedralization, each connection creating a neighborhood. Then, it is sequentially updated when meta-bers are merged: a neighborhood relationship is built between the new meta-ber and all the neighbors of its two original meta-bers, while the previous neighborhood relationships involving these two original meta-bers are removed. The evolution of the valence depends on the merging order of the bers, and therefore on the employed distance/similarity measure. A measure that favors geometrically welldistributed cylinders -composed of a similar number of original bers-preserves a bounded valence. In this case, the size of the ordered list (see Algorithm 1) remains linear in the number of original bers , and our strategy exhibits a total time complexity of ( log( )). These complexities were observed in our experiments (see Table III.1). Computation times are given in Table III.1 and were obtained on an Intel Xeon E5-1650V4 (6 cores, 12 threads @3.60GHz), combined with an RTX 2080Ti. In this table and in Figure III.7, we notice that the average valence remains relatively stable using these three similarity measures, even with one million of bers (Figure III.7d). This suggests that our method, with the proposed dissimilarities, should preserve the overall structure of the brain since bers are aggregated all around the tractogram and not only in some speci c areas. For a relatively small number of bers, our approach is slower than QuickBundles as can be seen in Figures III.8 to III.10. This is the case even though QB approach is single-core while our algorithm can be multi-core, and take advantage of CUDA for the rst similarity computations between neighboring bers. We also see here the lack of parallelism mentioned previously, as our multi-core version is at best 6 times faster than the single-core version, even though the CPU used has 12 threads. Most of the processing, including ordering the list, and creating the geometry of the new meta-ber, is inherently single-core.

III.2.3.2 Comparison with QuickBundles

III.2.3.3 Multi-resolution

For larger datasets, as we have in Figure III.11, our approach is faster than QB even when we use our single-core version, meaning that we scale better than QB does.

Another important point is that we are computing not only one resolution, but all possible discrete resolutions from the original bers to one single meta-ber, while QB only computes one xed resolution, depending on the user's parameters.

Our algorithm also includes the creation of the geometry of the new meta-ber, while QB just uses a prototype.

Obtaining a single meta-ber at the end is not as interesting as it might sound, as the geometry used here is unable to keep all details of quite di erent bers, resulting in a bad looking appearance as well as imprecise representation of the whole dataset (see Figure III.12). This is the reason why, for similarity measures de ned on inner product spaces, we also propose an automatic stopping criterion to prevent over-simpli cation (e.g. a single cylinder). Two cylinders are not merged if they are almost orthogonal to each other. This threshold needs to be set up manually as a maximum angle value between two bers. This will not allow grouping bers with an angle greater than the threshold. Furthermore, all bers not merged at the end of the algorithm (and cylinders representing very few bers) can be considered as outliers and thus discarded. Indeed, they were really di erent from all the other bers of the tractogram, so they are likely to be invalid. Figure III.13 presents the Thalamocortical with and without the outliers obtained when using a maximum angle of 70 • .

III.2.3.4 Dissimilarities comparisons

Figures III.14 and III.15 show the IFOF and the Uncinate Fasciculus, respectively at 5% and 3% resolution, with the three metrics we used in our experiments. These resolutions are too low to be correct, however, to compare the merging order of the bers, it is interesting to see what remains or not depending on the metric at these resolutions. earlier than MDF, and WC is also grouping some bers on the upper part, resulting in a huge blue meta-ber that is not what we would like there.

On the other hand, in Figure III.15, both MC and MDF merge some bers on top before WC. We can also see that MDF has grouped more bers on the bottom part.

III.3 Conclusion III.3.1 Summary

We introduced a new multi-resolution representation for brain tractograms that reduces the redundancy, easing the visualization and interpretation. It supports any similarity measure between bers. We also proposed a method to determine adequate candidates for e cient mergings of groups of bers. It is initialized with a Delaunay tetrahedralization of the extremities, which means that we inherently assume that two bers are similar if their extremities are close to each other. The employed dissimilarities, like the ones proposed here, should, therefore, be consistent with this assumption. Our progressive method makes it possible to display the input tractogram at any resolution in a continuous and real-time way. From a technical point of view, our two main contributions are a multi-resolution representation for tractograms based on a progressive decimation algorithm, and a combinatorial strategy based on a Delaunay tetrahedralization to make it computationally tractable. Visualizing groups of similar bers as single generalized cylinders and being able to easily change the level of resolution may be very useful for clinicians. For instance, it can help neurosurgeons identify relevant anatomical tracts that should not be severed during surgery, thus reducing postoperative complications and improving the clinical outcome. It is to note that we received very positive feedback from our neurosurgeons' colleagues of the Ste Anne hospital in Paris.

III.3.2 Limitations

The entire multi-resolution algorithm relies on similarity measures between bers or meta-bers. All the metrics used here are line-based, which is great when comparing bers obtained with tractography algorithms. However, in our case, we are creating some volumic data with our meta-bers, and using only their skeleton to compare them means that we are not taking advantage of our representation. The use of a volumic similarity between those meta-bers would help to increase the accuracy of the merges. Such a similarity would also require a more complex geometry than the generalized cylinders presented here, as it would bene t from more precision. Indeed, in Figure III.5, we see that parts of the meta-ber 2 are not included into the new meta-ber 3 . This geometry is a rst approach that needs inquiring for better representations. Based on our dissimilarities comparison (Section III.2.3.4), it appears that using di erent metrics depending on the bundle of interest could be an interesting experiment. Other elements also have to be considered, like the mathematical validity of the metric, its anatomical meaning, or the parametrization of some metrics such as WC.

Another limitation in this simpli cation is the use of the Euclidean distance between bers extremities for the Delaunay tetrahedralization. The correct distance to use would be the geodesic distance on the cortical surface, as it is the anatomically meaningful distance between those extremities. This computation is quite expensive with usual methods, but comparing the results is an interesting point that should be explored. Another possible improvement could be to use the closest bers to build the neighborhood relations (with > 1) instead of relying only on direct connections between extremities, however, performances would su er from such an approach.

IV

Geometric representations

As presented in Section III.2, a geometric representation for bundles of bers and meta-bers is required to improve the visualization quality as well as the understanding of those complex tractograms. In this chapter, we investigate a few possibilities that could be used for this purpose, analyzing their advantages and drawbacks.

IV.1 Lines and generalized cylinders

IV.1.1 Lines

Fibers represent the brain connectivity, in the form of synapses of neurons, allowing for electrical signals to navigate through the brain. Using lines is the natural -and easiest -representation for these connections. It is natural as synapses connect only two neurons together, meaning that there is no intersection between connections at the level of the axons. However, even with the highest sensibility scanner, we are unable to retrieve such thin connections. The reconstruction performed using dMRI is only capable of obtaining possible paths for sets of connections, and not the individual connections directly. Thus, what we call bers are in fact sets of bers, provided they exist (see Section III.1.3).

Lines are the most used representation for brain tractograms in most visualization software (3DSlicer [START_REF] Pieper | The na-mic kit: Itk, vtk, pipelines, grids and 3d slicer as an open platform for the medical image computing community[END_REF], trackvis [START_REF] Wang | Di usion toolkit: a software package for di usion imaging data processing and tractography[END_REF]...) and in a lot of scienti c papers. This is for two major reasons:

• the format in which bers are obtained is line-based, with a list of points and their connections;

• lines are easy to use, visualize, compare...

There are also some drawbacks of using lines, the rst one being visualization. Brain connections are inherently three dimensional, while lines are one dimensional, thus not providing visual depth information when rendered. They also tend to look like noise when visualizing a lot of them. As they do not have a volume, they are also not suited to represent sets of connections. This is usually not an issue when using directly all the created bers, but when simpli cation or clustering is performed, this tends not to encompass the real placement of the connections.

IV.1.2 Generalized cylinders

A possibility to circumvent the limitations of simple lines is to create a volume around them, making them the skeleton of a geometric representation. This can be achieved for instance using cylinders with varying sections, called generalized cylinders. For the original bers from the tractograms, the section can simply be a circle with a small radius. For bundles -i. e., sets of several bers-the section can then be adapted to best t the spatial extent of the underlying individual bers. This relatively simple geometry improves a lot on several points compared to using lines, mainly on perceived depth. In Figure IV.1, both images are taken from the same point of view of the same tractogram. The top image is rendered using lines, whereas the bottom one uses generalized cylinders. The depth is clearly perceivable with cylinders, while the lines appear to be all at the same depth, making the interpretation of the 3D position inside the tractogram very complicated. On top of that, cylinders can also encompass sets of bers. This representation presents advantages in terms of GPU visualization as it is possible to only create the sections pro le on the CPU, letting the tessellation unit creating the nal mesh directly on the GPU, thus reducing GPU memory (VRAM) usage and allowing re ning or coarsening the entirety or parts of the mesh on-the-y when needed.

This geometry also presents several drawbacks. Guaranteeing the continuity of the pro le along the curve can be di cult. Indeed, it means that the pro le at each consecutive section needs to be connected correctly to avoid torsions of the nal geometry. More importantly, all possible con gurations cannot be represented with only this geometry. For instance, when a bundle is composed of two possible termination positions on one side, the cylinder representation is unable to create such a separation, and it will ll the gap between the two di erent parts, which can be an issue.

In our use of such geometry in Section III.2.2.3, we are using elliptical pro les as they are easy to manipulate. They do not encompass all previous cylinders or bers really well, sometimes removing parts of the bundle from the geometry as in

IV.2 FastAPSS: Fast approximate point set surfaces

Algebraic Point Set Surfaces (APSS) [START_REF] Guennebaud | Algebraic Point Set Surfaces[END_REF] aim at creating a surface, with oriented point clouds as input (points with normals). This could be used on our tractograms as a way of merging meta-bers. Indeed, we could use at rst a representation of the initial bers as cylinders with a circular pro le, create points on this cylindrical pro le, and use them for creating a surface. However, APSS is most commonly used in a simpli ed version, which has some limitations, mainly the fact that it struggles to ll-in huge missing parts. In order to investigate possible improvements to this approach, we focus rst on the general case, before investigating its application to tractography in the next Section (IV.3). This work is currently under submission in an international journal.

IV.2.1 Introduction

Moving Least Squares (MLS) projection operators [Levin, 1998, Levin, 2004] o er an elegant uni ed framework to model Point Set Surfaces (PSS) from unorganized point clouds [START_REF] Alexa | Point set surfaces[END_REF], Amenta and Kil, 2004, Dey and Sun, 2005, Guennebaud et al., 2008, Alexa and Adamson, 2009, Öztireli et al., 2009], with immediate application to noise ltering and mesh reconstruction. At their core, MLS operators take the form of iterative projection procedures from R 3 to R 3 , whose stationary set is the PSS. At each iteration, a geometric primitive (for instance a plane) is optimized to t the point set under a weighting kernel centered at the evaluation point, before moving this point onto the primitive. This simple scheme is very exible: one can use di erent weighting kernels to control the spatial in uence of the input samples at the projection location, and leverage di erent kinds of primitives as long as they are e cient to t to a set of point. Moreover, its meshless nature quali es it for immediate processing of acquired 3D data such as laser scans or photogrammetry captures. In terms of quality, the resulting PSS is guaranteed to be smooth as long as smooth global kernels are used.

Unfortunately, such kernels induce untractable computations, even on inputs of moderate size, as the projection of a single point has a ( ) complexity, with the number of input samples. Consequently, it is common to localize the support of the MLS projection. This is typically achieved with two alternative strategies. The rst one consists in using only samples located in a ball centered at the evaluation location and adjusted to the size of the local weighting kernel -this induces that the MLS projection is not de ned everywhere in space, resulting in a PSS which cannot ll holes larger than the ball [START_REF] Alexa | Point set surfaces[END_REF]. The second one consists in gathering the nearest neighboring (kNN) samples to the evaluation location, which brings down the complexity to ( ( )) but breaks the smoothness guarantee [START_REF] Pauly | E cient simpli cation of point-sampled surfaces[END_REF]. The latter remains the standard choice in practice, with kept small to preserve good performances, assuming that a suitable data structure equips the input point cloud, such as a kD-tree in 3D.

Localized MLS projections act as a frequency lter, with for instance large kernels and a proportionally large value of resulting in a low pass geometric ltering. By contrast, if only a small number of nearest neighbors contributes to the tting of the geometry, high frequencies always dominate as only local sample points are taken into account (see Figure IV.3). Moreover, the underlying assumption that the kNN is a relatively smoothly-varying spatial function does not hold anymore when input point sets miss large parts, resulting in poor performance far from the input samples (see Figure IV.3). These shortcomings have led users to abandon MLS methods in favor of global reconstruction methods [START_REF] Kazhdan | Poisson surface reconstruction[END_REF][START_REF] Kazhdan | [END_REF], which often require a volumetric solver, and do not provide any simple point-wise projective form for all-purpose ltering, like denoising, analysis, or local/progressive mesh extraction. Our objective is to design a point-wise operator that retains the simplicity of MLS but scales e ciently to non-local kernels and reaches a robustness which is competitive with global optimizations [Kazhdan andHoppe, 2013, Huang et al., 2019] when it comes to large missing parts in the input point set. To do so, we propose a new MLS approximation scheme designed for non-compact kernels and use it to model Algebraic Point Set Surfaces (APSS) [START_REF] Guennebaud | Dynamic Sampling and Rendering of Algebraic Point Set Surfaces[END_REF]. In particular, we make the following contributions:

• we present a simple hierarchical scheme for the approximate convolution of any d-dimensional point data set with a smoothly-varying non-compact kernel, running in (log( ));

• we provide a detailed CUDA parallel implementation allowing for interactive projection of thousands of points in parallel;

• we introduce an e cient world-space sampling scheme for precise, fast, and memory-e cient surface mesh extraction.

We demonstrate its use for the APSS ltering and reconstruction of point sets presenting large-magnitude geometric noise as well as large missing parts.

IV.2.2 Smooth approximation of non-compact APSSurfaces

We consider as input an oriented point set P = {p i , n i , } (point, normal, estimated area), with ∈ [1, |P |], as well as a smooth averaging kernel H (•, •) : R 3 ×R 3 → R + .

For any point q in space, APSS aims at nding the best algebraic sphere on which to project it, so that it locally ts P. Weights describing the in uence of points p i over q are rst computed using the smooth kernel H as = H (q, p i ). Then, we construct an algebraic sphere

(X) = ( 0 , u 123 , 4 ) • (1, X, X 2 ), with 0 , 4 ∈ R and u 123 = { 1 , 2 , 3 } ∈ R 3 .
) n i , with - → (X) = (u 123 + 2 4 X)).
The reason for tting an algebraic sphere (with parameters ( 0 , u 123 , 4 )) instead of a geometric sphere (with parameters a center c and a radius ) is that the former degenerates gracefully towards a plane (with

c = - u 123 2 4 = (0, ||c|| 2 - 0 4
)

Following [START_REF] Guennebaud | Dynamic Sampling and Rendering of Algebraic Point Set Surfaces[END_REF], we rst align the gradients onto the normals, before adjusting the level-set to t the positions:

( ) = ( 0 , u 123 , 4 ) • (1, X, X 2 ) (u 123 , 4 ) = argmin (u 123 , 4 ) - → (p i ) -n i 2 0 = argmin 0 (p i ) 2
where are the estimated areas covered by each point in the point set (see Section IV.2.2.1 for details).

This results in the following closed form expression for :

         4 = 1 2 p i •n i - p i • n i / p i 2 - p i 2 / u 123 = ( n i -2 4 p i ) / 0 = - p i • u 123 + 4 p i 2 / . (IV.1)
With in hand, one can obtain the projection q by projecting q onto , and de ne the normal n q by taking the gradient of at q:

q = Project(q, ) n q = - → ( q)/ - → ( q) . (IV.2)
Finally, one can extract a surface by contouring the implicit function (q) = (qq) • n q, by using the Marchings Cubes technique [START_REF] Lorensen | Marching cubes: A high resolution 3d surface construction algorithm[END_REF] on a uniform grid or the Dual Contouring technique [START_REF] Ju | Dual contouring of hermite data[END_REF] on an adaptive octree, among other methods.

IV.2.2.1 Estimated areas

To cope with variable density in the input point set and avoid over-tting regions that are highly-sampled, Guennebaud et al. [START_REF] Guennebaud | Algebraic Point Set Surfaces[END_REF] use an estimated local density to parameterize the kernel locally, meaning that each input point in uences the space around it di erently. Our formulation is slightly di erent, as we use an estimated local area instead for our summation to best approximate a spatial integral over the point set, without changing the kernel.

To estimate , we descend in an octree until no more than N points remain per node (or maximum depth D has been reached). We then estimate the plane best tting the points in the cell and intersect it with the node bounding box to compute the node total area, which we nally divide equally between all points in the node. Other reasonable strategies already introduced in the literature could be used instead. For example, Guennebaud et al. [START_REF] Guennebaud | Algebraic Point Set Surfaces[END_REF] simply query the nearest input point sample to estimate the radius of a local disk patch, and [START_REF] Barill | Fast winding numbers for soups and clouds[END_REF] t a local 2D Delaunay triangulation to the neighborhood of the point once projected on its tangent space to estimate the local surface area around it.

IV.2.2.2 Smooth approximation: overview

Looking at Equation IV.1, we can factorize the weights for points that in uence similarly q (for example ∈ p i

∈ p i if ∀ ∈ ).
This trivial observation points to a simple algorithm for an approximate computation of the algebraic sphere tting the weighted point set (w, P).

In a preprocessing, we organize P into an octree, where we store in each node the following 9D statistics aggregating the geometry of all contained input points { ∈ }:

node ( ) =                    = ∈ ∈ R p = ∈ p i ∈ R 3 n = ∈ n i ∈ R 3 = ∈ p i 2 ∈ R = ∈ p i • n i ∈ R (IV.3)
When computing the algebraic sphere at query point q, we traverse the octree in a top-down fashion (starting from the root of the octree, descending into it iteratively), and aggregate the statistics of the nodes that are deemed "uniform enough" from the point of view of q. Following common strategies [Greengard andRokhlin, 1987, Barill et al., 2018], we estimate whether this is the case by looking if q is su ciently far from the bounding sphere of the nodes. Finally, given the approximate statistics, we compute an estimation of .

Weight estimation Considering the set of points inside a given sub-tree corresponding to a node , various techniques can be designed to estimate the weight associated with all these points. In our case, we simply consider the average position of the points inside the sub-tree (p = p ( )/ ( )) and estimate the weight as = H (q, p ) (IV.4)

More complex strategies involving more detailed statistics about the point distribution inside the sub-tree (for instance approximating it with a normal distribution or a Gaussian Mixture Model) could be used, but it would require integrating the kernel H against arbitrarily complex distributions, which is not always feasible in practice. We found in our experiments that our simple strategy o ered good results, even for challenging input point sets.

IV.2.2.3 Algorithm and CPU implementation

A naive traversal strategy could consist in simply stopping at a node if the distance to its bounding sphere is larger than a given predetermined threshold, and rely in this case solely on the statistics stored at the node to represent its entire sub-tree. Similar binary decisions are for example used in [START_REF] Barill | Fast winding numbers for soups and clouds[END_REF] to estimate the winding number of a surface around query points, which is harmless in this case because the output of the algorithm is, in the end, an integer value, and the resulting approximation is below the nal quantization error.

In our case however, this simple strategy cannot be used, as it results in non-smooth statistics, and therefore in a non-smooth surface (for example the top of the face in inset). We describe next a simple hierarchical scheme that is adapted to the blending of any n-dimensional data, and which we use to blend the statistics in a smooth manner.

Smooth hierarchical partition of unity Looking at Figure IV.4 (left), we make the trivial observation that the bounding sphere of a node contains the bounding spheres of its child nodes entirely (these are tangent at some corner of the child's bounding box). By scaling these bounding spheres by the same factor > 1, we construct protection spheres {S( )} for all nodes { } that have the property that the sphere of the father strictly contains the spheres of its children. We rely on these protection spheres to design our hierarchical partition of unity.

Protection spheres (solid)

Node weights (black) and sub-trees weights (colored):

Figure IV.4 -Node/subtrees local partition of unity illustrated on a 2d example. { } (colored) denote the set of subtrees of the node (black). First: Bounding spheres (dotted) and protection spheres (solid: scaled bounding spheres) of the father node (black) and one of its child nodes (red). Others: The respective distances to the scaled bounding spheres induce a hierarchical partition of unity. Note that these weights do not sum up to 1, because the quantities that are averaged are integrated (cumulated) functions over di erent integration domains, and those domains depend on the level at which they are de ned.

When traversing the octree, we stop in a sub-tree (corresponding to a node and its child subtrees) if the query point q is outside the protection sphere S( ) of the node. Otherwise, if q is inside S( ), we blend the statistics of the node with the ones of its child sub-trees ℎ( ), based on the respective distance to the protection spheres. In this more complex case, we compute the following statistics for the entire sub-tree: tree ( , q) = ∈ ℎ ( )

(1 -(q)) × tree ( , q) (IV.5)

+ ∈ ℎ ( ) / × × (q) × node ( ),
(q) describing the two-sphere local blending (illustrated in Figure IV.4). The signed distance from the point q to a sphere of center c and radius being (q, ) = qc -:

(q) = Ω( (q, S( )), (q, S( ))), with (IV.6)

Ω( , ) =          0 if ≤ 0 (i. e., q ∈ S ) 1 if ≥ 0 (i. e., q ∉ S ) Ω ( /( -)) otherwise, with Ω ( ) = (-(1/( -1) ) / 2 )
Note that this scheme induces an in nitely di erentiable local partition of unity between the interior of the child's sphere S and the exterior of the parent's sphere S , because Ω : [0, 1] → [0, 1] is a strictly increasing function with null derivatives (of all orders) at 0 and 1 (see inset). The proof of this claim is given in Appendix A. Note that we also used Ω ( ) = (1cos( ))/2, which results in a 1 surface only (as the derivatives of higher orders are not null on the boundary), and we could not observe visual artifacts in practice.

We summarize in Algorithm 2 a pseudo-code of our projection operator as described in this section.

Algorithm 2: Computation of ( q, n q) = Project(Point q) = tree ( root , q) // see Algorithm 3 compute 0 , u 123 , 4 from // see Equation IV.1 q = Project(q, algebraic sphere ( 0 , u 123 , 4 )) : if ( 4 == 0) then // q: project q on plane 0 + u 123 • X = 0 else c = -u 123 /(2 4 ) = max(0, c 2 -0 / 4 ) // q: project q on sphere of center c and radius end n = 2 4 q + u 123 n q = n/ n Algorithm 3: Computation of tree (Node , Point q) // returns the cumulative stats for tting an algebraic sphere to the tree descending from node if is leaf then // accumulate statistics over points in the leaf return ∈ H (q, p i ) × node ( ) else if q ∉ S( ) then // q is su ciently far from the node return H (q, p ) × node ( ) else // blend statistics between and its children // (q) : blending function between and (EquationIV.6)

return ∈ ℎ ( ) (1 -(q)) × tree ( , q) + ∈ ℎ ( ) / × H (q, p ) × (q) × node ( ) end end
Analysis We stress that in Equation IV.5, one has to distinguish the node statistics node (computed according to Equation IV.3 and stored in the node) from the tree statistics tree that require a recursive evaluation in this formula. If q ∈ ( ) \ ( ) with ∈ ℎ( ), then tree ( , q) = H (q, p ) × node ( ), and in this case (1 -(q)) tree ( , q) + (q) node ( ) represents the blending between the statistics of and at the depth level of . However, if q ∈ ( ), then tree (q, ) blends further the statistics of 's children, with the recursive formula of Equation IV.5.

Note that the factor / appears in Equation IV.5 to account for the fact that while the node has to be blended with its more detailed representation (i. e., the entire set of its children), the local partition of unity if performed for each child independently, for which we can expect an amount of information given by / . For example, in the extreme case where all points are located inside the same sub-tree (i. e., 7 of the 8 children are empty), the summation is performed over this single sub-tree only, and the partition of unity is truly performed between the node and its only child.

Another way of understanding why this scheme is necessary is that the blended quantities are integrated (cumulated) functions over di erent integration domains (the domain being restricted to the cell where it is stored). Doing so, we equivalently compute in ne the integral over the entire domain of the blending between the approximations of each function over the various levels of detail of the hierarchy.

Note also that Equation IV.5 describes the case where q is outside the protection sphere of adequately as well, since in this case (q) = 1, ∀ ∈ ℎ( ), and / = 1, which leads to tree ( , q) = H (q, p ) × node ( ) in the end.

IV.2.2.4 GPU implementation

To further improve performance, the octree traversal (Algorithm 3), as well as the projection operation are also implemented in GPU, using CUDA. The octree is traversed on the GPU in an unrolled recursion using a while loop. Traversal is done both ways as values computed from the bottom of the tree are needed for the evaluation of the statistics. Each thread computes the statistics corresponding to one point to project, and then compute the projection of the point. This is presented in Algorithm 4. An example of a traversal is shown on a binary tree (instead of an octree) in Figure IV.5. Each node and arrow color corresponds to a line in the Algorithm 4. When reaching a stopping condition (a leaf, a point outside an external bounding box, or in-between two bounding boxes), necessary elements of the tree are evaluated, and then results are added at every level until coming back to the main node, which will contain the correct evaluation for the point being projected. This means that each thread needs to store the statistics at each possible depth.

The cumulated statistics and the last child traversed have to be known at each depth (except the child for the leaves) for the traversal. These values are illustrated in is taken on local memory. The maximum depth is xed at 12, which seems to be a good compromise in terms of precision and speed. This means that a single thread needs enough registers for 12 times the statistics and 11 integers for the children. These statistics are necessary at each level. Indeed, the example described in Figure IV.5 needs to store the statistics at depth 1 and 2 when going through nodes 3, 5, and 6 as the statistics of node 2 were previously computed and added into Statistics[1] (statistics at depth 1).

The octree is stored in global memory in breadth-rst order, even though traversal occurs in depth-rst order. This is done so that as many points are evaluated at the same time, the rst threads at a given depth will be loading from global memory values from neighboring nodes that will be likely to be used by other threads. To improve memory latency -which tends to be the most limiting factor in our implementation -the kernel parameters and the rst depth levels are stored in shared memory. It does not result in a major improvement but we still project points a little faster by about 0.7%.

IV.2.2.5 Performances

The CUDA implementation reduces by one order of magnitude the computation time to project one query point, compared with the CPU implementation (see Table IV.1). We used a GTX 1080Ti and a Xeon E5-1650 v4 (3.6 , 12 threads) for all our tests.

When projecting a signi cant amount of points in parallel, we rst sort them using the Morton code [Morton, 1966], to favor traversal coherency among concurrent threads. It consists in interleaving the bits of the 3D points rounded in a grid of a xed size. Points being close together in the 3D space obtain a close Morton code value, making it feasible to order them by localization in space. As a result, the nodes encountered in the octree will most likely be similar for a large set of collocated points, resulting in improved memory access as well as improved synchronization. We observe in practice a consistent gain in performance (see

IV.2.2.6 Practical projection strategies

Given an input query point q, we have seen so far how to compute an algebraic sphere and how to project q on it. The point q can further be projected onto the APSSurface by advecting it along the following vector eld:

f (t) = Project(f (t)) -f (t) (IV.7) f (0) = q
One simple strategy commonly used in the literature to approximate the ow presented in Equation IV.7 is to simply replace q by (q) iteratively (i.e. discretizing Equation IV.7 with a xed time step Δ = 1) and stop after a xed number of iterations or when the displacement is smaller than a xed threshold.

This strategy works well in practice for points that are located around the APSSurface (for example, the input points themselves in a denoising scenario). However, when projecting points from afar (which is made possible by the use of a non-compact kernel), we observe in practice that it is important to bound the displacement at each iteration by a reasonable threshold. In our implementation, we use min( , /10), with being the diagonal of the point set's bounding box. The main reason for it is that, when using a kernel with low-decreasing tails (which can be necessary in practice because of numerical instabilities), the point set is seen as a single sphere for all points su ciently far away (i.e., the weights (q) are almost constant for q ).

IV.2.2.7 Practical non-compact kernels

While most related works advocate the use of Gaussian kernels, we remark that this type of kernels is not appropriate for surface de nition in the presence of missing data (see Figure IV.7). The Gaussian has a tail that vanishes too quickly, and it is in practice often impossible to nd adequate parameters using this kernel that provide the desired intent. Figure IV.7 illustrates this point: it is for example impossible, on this input point set, to set a support allowing for smooth surface completion far away from the input samples while retaining high-frequency geometric details near them.

Approximating Interpolating

Input Gaussian Gaussian mixture Rational Rational singular

Figure IV.7 -As they vanish too quickly, Gaussian kernels are not suitable when approximating a surface with missing parts (left). Instead, Gaussian mixtures (center) or rational kernels (right) are more robust to holes. Whereas all these kernels are approximating, rational kernels can also be interpolating.

We advocate the use of either Gaussian Mixtures or Rational kernels:

H GM { } (q, p) = ∈ -3 -q-p 2 2 2
H Rat , (q, p) = ( qp 2 + ) -/2 (IV.8)

Note that setting to 0 when using a rational kernel results in an interpolatory surface, as the kernel becomes singular on the input samples. Interestingly, it allows the user to navigate continuously between interpolating ( = 0) and approximating ( > 0) surfaces. In practice, however, we cannot have < 10 -3 , as it might result in some weights being evaluated to 0 or almost 0, yielding to division by 0 (Equation IV.1) and impossible projection.

A simple and e ective strategy to set appropriate parameters for a Gaussian Mixture kernel is to take = 0 . We use this setup in all examples featuring APSSurfaces parameterized by this type of kernels in this thesis.

IV.2.3 Surface extraction

Having a dense sampling that closely follows the algebraic surface is key for multiple applications. It is for example necessary for the rendering of the surface using splats.

Perhaps more importantly, it allows for a simple, fast, memory-bounded mesh extraction of the surface. Contouring the implicit surface is commonly done using the Dual Contouring approach on an adaptive grid (an octree in 3D) that has been prepared in a preprocessing in order to o er good precision near the surface to be contoured, and less detail further away from where the surface is expected to lie. It allows contouring the surface using a volumetric structure with time and memory complexities that are linear in the size of the output extracted surface ( (1/ℎ 2 ) instead of (1/ℎ 3 ), ℎ denoting the size of the leaves of the structure [START_REF] Ju | Dual contouring of hermite data[END_REF]). Since we target point-based modeling of surfaces arbitrarily far away from input samples, we cannot simply prepare the octree using those input samples, as is often done in the literature. Figure IV.8 illustrates the type of artifacts that can be observed when contouring the implicit function on an octree fed with the input samples only. Large parts of the surface can be entirely missed (top) and the precision/connectivity of the resulting surface can be largely a ected away from the input data (bottom).

We propose a simple, yet e ective strategy to sample densely the APSSurface inside an adaptive octree that we construct iteratively, following these two simple observations:

• an octree cell in which a 3D point is projected contains most probably part of the algebraic surface: we can re ne the sampling inside this cell in a hierarchical manner by simply subdividing it and further investigating its children.

• an octree cell containing part of the surface is likely to have neighboring cells containing parts of the surface, as the implicit functions we construct are smooth.

As the surface might expand way outside the bounding box of P, we need to be exible in the generation of the octree; in particular, we allow nodes to exist anywhere in space, and index nodes via ( , , , ), with the depth of a node and ( , , ) its index in the grid at depth . The initial root is (0, 0, 0, 0), and spatially represents the bounding box of P.

Input Classical Ours

Figure IV.8 -The classical way to derive an octree from a point cloud focuses only on the existing points and completely misses the holes. When meshing the APSSurface on such an octree (via Dual-Contouring [START_REF] Ju | Dual contouring of hermite data[END_REF]), this yields a coarse or invalid surface. Instead, we propose a new way to generate an octree that will encompass the full APSSurface.

With this de nition, our goal is to generate nodes at depth following the surface, corresponding to potentially multiple nodes at depth 0. The key idea is to project nodes, by projecting their center and retrieving the node corresponding to the resulting position.

We propose the following re nement process: starting from nodes N 0 at depth 0 (initially, only (0, 0, 0, 0)), we obtain the nodes N by alternating node projection and node subdivision (Figure IV.9).

projection subdivision projection subdivision

Figure IV.9 -We re ne an octree by alternating 2 steps: projection, where we project the center of nodes and keep only nodes containing these projected centers, and subdivision, where we subdivide all retained nodes. This re nement process is repeated until a target depth is reached.

Some nodes in N might not have an ancestor in N 0 (they are outside the starting nodes), notably in the case of large holes in P. Hence, we add these unprocessed ancestors at depth 0 to N 0 and re ne them; this process is repeated until no new nodes need to be added to N 0 .

Since we want the nal octree to have depth , we need to adjust the re nement depth . Initially, = , but we adjust it in function of N 0 , as

= -log 2 max(|N 0 | , |N 0 | , |N 0 | ) .
This also means that we traverse less and less nodes when extending N 0 , limiting the memory and performance costs. Finally, we re-index the resulting nodes to form an octree (only 1 root) of depth that closely follows the APSSurface. To upsample P, we project points from the leaves. To create a mesh, we use Dual-Contouring [START_REF] Ju | Dual contouring of hermite data[END_REF] on either this octree (our experimental setup) or an octree derived from the upsampled point cloud.

Algorithm 5 is a pseudo-code summarizing the approach.

Algorithm 5: Surface-following octree.

Function P N (n): c ← center of p ← project c ← node at p at depth(n) + 1 return Function N O S (N , D): for ∈ [1, ] do M ← { P N (n) | ∈ N } N ← { ℎ ( ) | ∈ M } N ← { (N ) | ℎ ( + 1))} end return M ; Function C O (D): N ← ∅ ; N 0 ← ∅ ; M 0 ← { root } ; ← while |M 0 | > 0 do N ← N ∪ N O S (M 0 , ) ← (1, -log 2 max(|N 0 | , |N 0 | , |N 0 | ) ) A 0 ← { 0 ( ) | ∈ N } A 0 ← (A 0 ) N 0 ← N 0 ∪ M 0 ; M 0 ← A 0 \ N 0 end
Re-index N , and generate ancestors return N

IV.2.4 Evaluation

Our method was designed with the intent of i) preserving the main advantages of MLS methods (suitable for point set ltering using simple averaging kernels, useful to project any 3D point on the MLS surface independently, the MLS surface can be de ned either as a collection of splats using the projection operator or as a mesh by contouring the associated implicit function, depending on the target application requirements), while ii) reaching the robustness and stability of global methods such as Poisson reconstruction in presence of large missing data.

We compare our method mainly with local APSS with 20 nearest-neighbors and global APSS (see Figure IV.11), but also with robust methods such as MPU [START_REF] Ohtake | Multi-level partition of unity implicits[END_REF], Screened Poisson reconstruction [START_REF] Kazhdan | Screened poisson surface reconstruction[END_REF] and Point2Mesh [START_REF] Hanocka | Point2mesh: A self-prior for deformable meshes[END_REF] (see Figure IV.12). VIPSS [START_REF] Huang | Variational implicit point set surfaces[END_REF] is also considered when needed, with its limitation being the number of input points. In Figure IV.12, the input models for this method were all reduced to 7000 points, which is the maximum our machine equipped with 64 of RAM can handle.

IV.2.4.1 Accuracy

As shown in Figure IV.10, our method o ers a smooth approximation of global APSS. While the exponential blending function Ω is slightly more accurate than with the cosinus Ω , the e ect is limited, whereas the sphere scaling factor has a major impact on the quality, as it determines the blending between levels.

Increasing provides results closer to the ground truth, but with an increased computational cost. Moreover, we observe that the number of traversed nodes depends mostly on the sphere scaling factor , as well as the location of the points to project. Additional experiments (for up to 4) suggest that the number of traversed nodes follows a complexity of ( 2 ). We nd that = 2 gives good results in all our experiments, with no perceptible improvements in the output quality when increasing , although this is dependent on the type of kernel that is used.

As

IV.2.4.3 Resilience to missing data

We analyze our approach on synthetic point sets and real-life scans featuring uneven sampling and large missing parts, notably by arti cially removing parts of the samples on some examples. On most extreme tests, we only keep input points distributed along 3D planes, sketch-like, similar to what Huang et al. [START_REF] Huang | Variational implicit point set surfaces[END_REF] show in their VIPSS article.

We compare our results to two robust methods in Figure IV.15: Screened Poisson reconstruction [START_REF] Kazhdan | Screened poisson surface reconstruction[END_REF] and MPU [START_REF] Ohtake | Multi-level partition of unity implicits[END_REF]. We also compare to local APSS for completeness, and to VIPSS [START_REF] Huang | Variational implicit point set surfaces[END_REF] on subsampled inputs as their complexity -(|P | 3 ) -made execution impossible on large models.

We observe that only the Poisson reconstruction and our method cope nicely with missing data, but we o er a more natural completion on holes such as missing back of the head (second line in Figure IV.15). Note that while P contains both positions and normals, VIPSS only uses positions as it optimizes for the normals, and this can skew the results. A normal estimation step has been presented in the original article of Guennebaud & Gross [START_REF] Guennebaud | Algebraic Point Set Surfaces[END_REF], but it remains to be demonstrated that such an approach could work on input point sets featuring this type of degeneracies.

IV.2.4.4 Robustness to noise

To study how resilient to noise our method is, we arti cially degraded artifactfree inputs by introducing noise both on the positions and on the normals of the samples (Figure IV.16), and by removing parts of the shapes (Figure IV.17). Our experiments suggest that our method is more sensitive to position noise than to normal noise, especially in the presence of large holes in the data. This sensitivity remains rather modest however, even far away from the input samples, in contrast to what has been previously reported in the literature for APSS [START_REF] Berger | A benchmark for surface reconstruction[END_REF]. This di erence can be explained simply: while it is common to rely on a small number of nearest neighbors, in which scenario the APSSurface is indeed extremely subject to noisy samples, we rely on all samples to t the algebraic spheres.

Input larger kernel no noise position noise normal noise

Figure IV.16 -Our approach allows for noise removal by increasing the kernel extent. Here, we show the result with a Gaussian noise applied on the positions as well as on the normals orientations.

IV.2.4.5 Reliability of upsampling and application to surface mesh reconstruction

For practical applications, we need the surface contouring and the point cloud upsampling to be robust to missing data. We can see in Figure IV.18 that even with very large missing parts, our sampling scheme allows us to reliably reconstruct the APSSurface, the latter also being robust to missing parts. In this scenario, the initial estimation of the extent of the model (initialized with the bounding box of P) can be very di erent from the extent of the original input. Apart from point Figure IV.17 -Our approach allows for noise removal by increasing the kernel extent. Here, we show the result on data with missing parts with a Gaussian noise applied on the positions as well as on the normals orientations.

cloud upsampling, the construction of our octree also allows contoured meshes to keep an acceptable vertex density even inside missing parts (Figure IV.8). We observe in practice that the number of octree nodes we process is somehow similar for all experiments (see Table IV.2), even on inputs whose input geometry varies very strongly (Part Sphere 3's bounding box volume is roughly 1000000 times smaller than the one of Sphere's bounding box in Figure IV.18), validating that our approach allows investigating the 3D space e ciently. This is critical, as we can in practice make no assumption on the size of the output surface. The target level of detail has to be therefore re-evaluated on-the-y during octree construction, as the estimation of the extent of the output surface is re ned, in order to avoid over-consumption of computer memory.

Sphere

IV.2.4.6 Application to point set ltering

Finally, while we have focused so far on extracting surfaces with the highest amount of detail, we stress that our approach is suitable for ltering the input point set, an application for which compared methods are not suitable (Figure IV.19).

For the sake of completeness, we also compare to the standard -nearest neighbors based APSS approach and demonstrate that ultimately, to allow for low-pass ltering, it is not enough to use a large kernel (while keeping the number of neighbors low) and it is essential to consider a large region of the input point set around the evaluation point. Filtering can happen with kNN, however, it comes at the expense of looking for more nearest-neighbors, leading to slower projections overall when our approach becomes even faster when using larger kernels as we do not need to descend the octree as much. Such e cient low-pass lterings can allow for applications that leverage APSS in their machinery, such as multi-resolution point set deformations [START_REF] Nader | Adaptive Multi-scale Analysis for Point-based Surface Editing[END_REF] or point multi-resolution descriptors in the context of data analysis [START_REF] Mellado | Growing least squares for the analysis of manifolds in scale-space[END_REF].

We expect that improving upon the APSS method by allowing for e cient ltering of input point sets featuring artifacts typical of real-life scans will have an impact on both applications.

IV.2.5 Discussion, limitations, and future work

In order to o er instant feedback to users editing the input weighting kernel, our approach was designed with the constraint of remaining entirely agnostic to it, both for the construction of the structure and for its adaptive traversal. We believe for this reason that our approach can be used for a variety of geometry processing tasks requiring the e cient approximation of the convolution of large amounts of 3D data against a smooth kernel. Note that our approach was introduced for the 3D case, but that it makes no real assumption on the dimensionality of the ambient space and that it can trivially be extended to any n-dimensional Euclidean space. However, our approach could be improved in a number of ways, that merit future investigation.

First, our traversal criterion relies on the geometry of the protection spheres and the location of the evaluation point only. This criterion could be re ned by taking the distribution of the points inside the nodes into account. Indeed, even for a coarse node with a large bounding box, if all input points are located at the same position, there is no need to further traverse the tree, as the weight given to the mean point is the same for all points, and no approximation is made at all. Taking another con guration, if all points are distributed mainly along a single direction, assuming we use a radial kernel, the approximation that is made is largest if the gradient of the kernel is aligned with the main direction of the distribution. The squared ( rst order) error that is made on the weights { } = {H (p i , q)} when approximating all weights by the one H ( p, q) de ned at the centroid p is:

( ) = (H (p i , q) -H ( p, q)) 2 1 order 1 (q) = - → p H ( p, q) • (p i -p) 2 = - → p H ( p, q) • ((p i -p) • (p i -p) ) • - → p H ( p, q) .
By storing in each node the covariance matrix = p i • p i , the centered covariance matrix ¯ = ((p ip) • (p ip) ) can be recovered as ¯ = -(

) p•p . The squared ( rst order) error can then be evaluated in (1)

as 1 ( ) = - → p H ( p, q) ¯ - → p H ( p, q)
. Since the kernels we use are smooth, this function is smooth with respect to q. This opens an interesting avenue for anisotropic approximation: by making the size of the protection spheres smoothly dependent (between > 1 and ) on this tting error, it could be possible to exploit the large anisotropy of the input data to improve the adaptive traversal e ciency. This extension of our work is however non-trivial: one would have to design an adaptive scheme ensuring by construction that the child's protection sphere remains strictly inside the father's protection sphere, to ensure the validity of our partition of unity. This imposes bounds on the spatially-varying, node-dependent (q, ) (2 (q, ) > 1 + (q, child ) ∀q, ), which could be veri ed by construction in a preprocessing analysis of the octree. While we believe that little bene t would be observed for our application scenario in practice, it may reveal important for cases where the intrinsic dimension of the data (here, 2D) is much smaller than the dimension of the ambient space (here, 3D).

Secondly, while we o er instant feedback to users editing the kernel, it can be tedious to search appropriate parameters to obtain the desired output. We believe that the design of complex, spatially-varying kernels (see inset for a simple example) with the sole application of real-time point-based modeling would bene t the Shape Modeling community and that our framework is appropriate for its illustration.

Finally, while we have studied how large support kernels can be used for the ltering of point sets and their robust surface extraction (in particular, far away from the input samples), we believe that our methodology can be used for the parametric design of smooth skeletons. A number of existing techniques extract the medial axis of point sets using growing spheres tting the data (e.g., [START_REF] Rebain | Lsmat least squares medial axis transform[END_REF]) or clustering input samples advected inside the surface (e.g., [START_REF] Cao | Point cloud skeletons via laplacian-based contraction[END_REF], Huang et al., 2013]), and we believe that our technique can be used for the smooth, parametric tting of maximally-inscribed spheres as well, thus allowing for interactive shape modeling applications such as skeleton design and interactive quad mesh modeling [START_REF] Ji | B-mesh: A modeling system for base meshes of 3d articulated shapes[END_REF].

IV.2.6 Conclusion

We have presented a smooth scheme for approximating e ciently the convolution of multi-dimensional data against a smooth spatial kernel, and have used our construction to compute e cient approximations of Algebraic Point Set Surfaces. We have demonstrated that the use of ad-hoc kernels that are simple to manipulate, such as Gaussian Mixtures and Rational Singular kernels, allows for intuitive continuous point set ltering as well as high-quality surface extraction of input point sets featuring large missing data, which are typical of real-life scanned data. Our approach competes with all recent global surface extraction techniques while retaining all advantages of MLS local techniques such as the ability to project any 3D point on the surface in a smooth manner independently. Being entirely agnostic to the input weighting kernel, our approach opens the route to real-time point-based modeling by casting the problem to the design of parametric kernels guiding the APSSurface.

IV.3 FastAPSS, application to bers

In this section, we apply the FastAPSS projection operator to create a geometrical representation of bers and meta-bers. We will then use it directly inside the multi-resolution representation introduced in Section III.2, thanks to a new metric we propose which uses the volume of meta-bers.

The bers obtained from tractography algorithms are simple 3D polylines consisting of linked points. In order to create a point set with normals corresponding to these bers, points are sampled on a cylinder with a circular pro le following the 3D polyline. From our experiments, we noticed that better reconstructions happen when points are almost equidistant to each other. To reach this, we suppose that we either have bers with already equidistant points (which is the case for some tractography algorithms) or that we are able to resample the bers. Then, the radius of the circular pro le is chosen equal to the distance between the points ( ), and six points are created on each pro le, resulting in almost equidistant points on the pro le (they are not exactly equidistant due to the angles between consecutive segments). The pro le at a given point p is created along a plane which is the average of the two perpendicular planes corresponding to both segments originating from p. This allows us to easily de ne the normal at each point as the normalized vector from the original point p of the ber to the newly created pro le point. On top of that, the extremities are kept as such, without creating a pro le, and using as normal the direction of the segment starting at the extremity. We will call these points pro le points.

IV.3.1 Volumetric considerations for the comparison of bers

Using the points and normals created above, it is possible to apply the FastAPSS algorithm to build a surface corresponding to the bers. For this purpose, a kernel needs to be chosen. Experimentally, we choose a multiple Gaussians kernel (see Equation IV.8), containing four standard deviations, the rst one being the same as the distance between points . The factor used between the four standard deviations is 1.2. This kernel is de ned as:

H ber { 0 = } (q, p) = 3 =0 -3 -q-p 2 2 2 | ∀ ∈ [1, 3], = 1.2 × -1 (IV.9)
Figure IV.20 shows an example of a ber, its sampling, and the mesh reconstructed using FastAPSS.

IV.3.1.1 Building meta-bers using FastAPSS

To de ne a meta-ber using FastAPSS, the only requirement is to obtain pro le points, so that a corresponding surface can be built, and a volume computed. One really interesting property of the Moving Least Square projections is the ability to quickly project points onto the nal surface. Let us consider the merge of two meta-bers ( 1 and 2 ) into another meta-ber ( 3 ). We note the surface corresponding to meta-ber . The pro le points P of 3 can be de ned as:

P ∈

3 ⇔ (P ∈ 1 and P ⊄ 2 ) or (P ∈ 2 and P ⊄ 1 ) (IV.10) This means that we can project points of 1 onto 2 and vice-versa, and retain only the points that were not inside the other meta-ber (see Figure IV.21). These points are by de nition (Equation IV.10) the points belonging to 3 . Extremities are considered separately, as they remain a core part of brain connections:

EP ∈ 3 ⇔ EP ∈ 1 or EP ∈ 2 (IV.11)
This means that a meta-ber keeps all previous extremities of the original bers it represents. These extremities are divided into two di erent groups, that we arbitrarily call beginning and ending . Each beginning and each ending point with respect to their original ber are always forced to not be in the same group. We compute

( 1 , 2 ) = (p 1 ∈ 1 ,p 2 ∈ 2 ) (p 1 , p 2 ). If: ( 1 , 2 ) + ( 1 , 2 ) ≤ ( 1 , 2 ) + ( 1 , 2 ) (IV.12) then 3 = 1 ∪ 2 and 3 = 1 ∪ 2 , otherwise 3 = 1 ∪ 2 and 3 = 1 ∪ 2 .
Figure IV.21 -Pro le scheme of the merging of meta-bers 1 (blue) and 2 (green) into 3 (red). Points belonging to one meta-ber surface and inside the other meta-ber volume are in gray and are discarded (Equation IV.10).

IV.3.1.2 Volume aware similarity

To improve meta-bers comparisons, we said in Section III.2 that it would be bene cial to consider meta-bers volume and not only their skeletons, as is done with usual line-based similarities. The obtained triangular mesh allows for an easy volume computation. Indeed, for a mesh containing triangles, with the coordinates of the vertices (a, b, c) of triangle being p ka = ( , , ) (and similarly for b and c), the volume can be computed as [START_REF] Zhang | E cient feature extraction for 2d/3d objects in mesh representation[END_REF]]:

= -1 =0 + + 6 - - - - (IV.13)
This formula remains valid as long as the mesh is closed and contains only valid triangles (no intersections between triangles for instance).

When merging meta-bers, the objective is to combine the most similar ones. We already explained that it means they should connect the same functional regions and pass through the same local regions inside the white matter. Ensuring this proximity when having an encompassing and accurate representative geometry is possible by comparing the volume variation between the two merged meta-bers and the new representative meta-ber. Indeed, the volume of the new meta-ber will not be a lot bigger than the volume of each of the merged meta-ber if they are close to each other. This needs to be combined with a measure of the distances at the extremities. Taking inspiration from the Weighted Current structure [START_REF] Gori | Parsimonious approximation of streamline trajectories in white matter ber bundles[END_REF], we propose here to use the variation of volume to create a new volume-based similarity:

( 1 , 2 ) = =1 =1 ( ( M 1 b a -M 2 b c ) ( M 1 e a -M 2 e c )) ( ( 1 ∪ 2 ) ( 1 ) + ( 2 ) ) (IV.14)
where 1 (respectively 2 ) is a meta-ber representing (respectively ) original bers, M x b y and M x e y are respectively the beginning and ending extremity of meta-ber , and , , and are three Gaussian kernels parametrized by , , and , respectively. 1 ∪ 2 corresponds to the newly created meta-ber 3 .

Using this similarity instead of line-based one should allow the removal of some unwanted early merges, as were exposed in Section IV.1. (c) Mesh corresponding to the other original ber.

(d) Meta-ber corresponding to their merge.

Figure IV.22 -Example of a meta-ber obtained using FastAPSS.

However, the merge fails in some cases (Figure IV.23). This bad case happens mostly when we are merging two meta-bers that are really close to each other. As shown in Figure IV.24, this kind of result occurs when on the pro le, most points are considered inside the other meta-ber, resulting in them being discarded. This can be mitigated using a threshold on the distance inside of a meta-ber, however, it did not solve all the issues, as the kernel also plays an important role in the creation of a correct mesh. With the kernels we tested in FastAPSS (Gaussian, rational, and multiple Gaussian), we have sometimes strange reconstructions happening for the bers and the meta-bers. The resulting geometry (Figure IV.23) has many degenerated triangles as well as missing faces. This makes the volume computation meaningless, leading sometimes to negative values. These values being used in the ordering of the merges, it breaks the whole chain and cannot be used in its current state. Manually adjusting the kernel -even using only multiple Gaussian kernels -could result in a better result. However, this would not be compatible with an automated multi-resolution construction.

(a) Mesh corresponding to one original ber.

(b) Mesh corresponding to the other original ber.

(c) Meta-ber corresponding to their merge.

IV.3.2 Discussion

IV.3.2.1 Summary

FastAPSS projections are an interesting approach to obtain a geometric representation for meta-bers. It presents some advantages, mainly in the improvement it brings in terms of variety of possible geometries, allowing for better progressive grouping. Indeed, it is capable of grouping geometries seamlessly as can be seen in Figure IV.22. The mesh created also allows for an easy volume computation, providing another element to compare meta-bers. The proposed metric should improve the merging order with respect to the representative geometry, even though we were unable to test it.

IV.3.2.2 Limitations

Some work and investigations remain to better t fastAPSS approaches to brain bers. However, some limitations already appear with these rst attempts. As far as visualization is concerned, a few drawbacks emerge. The rst one is memory consumption. Creating a mesh for every single meta-ber is costly, as lots of triangles are needed to obtain an accurate surface. Typical computed tractograms easily weigh a few gigabytes for recording only the points linked together as lines ( bers). Here, we would have to store many more points as well as their connectivity. What is more, such high triangle count data is always a challenge to display with good performances, especially when some additional rendering techniques are bene cial (see Section V.3). Another big limitation is obviously the kernels, which prevented us from even being able to obtain some complete rst results. This limitation is the same as for FastAPSS, and there is a need for other kernels that would adapt to the subjacent geometry. Concerning the similarity from Equation IV.14, it would also be important to show if it is actually a distance measure and if the volume is to be considered inside a Gaussian kernel or not. Providing that we can manage all these limitations, there is another limit that appears here: intersections between di erent meta-bers. Indeed, meta-bers can become quite large, while being close to other larger meta-bers, which can lead to intersections between them, especially at the extremities. To prevent this issue, one could think of processing the whole set of pro le points from all bers and meta-bers as just one single geometrical element. Finding a kernel able to accurately t all this geometry would be challenging, providing that we are even able to build an octree on so many points that would t in GPU memory. This sums up in a simple observation: FastAPSS projection should be great to improve the merging order, but another geometry is required for better visualization.

IV.4 Implicit surfaces

Implicit surfaces are mathematically de ned by an equation ( , , ) = 0. This function can be a distance function to a skeleton or a mathematical de nition of a geometry, like a sphere for instance. It gives a value to any point in space, relatively to the surface, meaning that evaluating ( , , ) -= 0 represents a way to shrink or expand the surface.

One interesting property of implicit surfaces is the operators that can be used between implicit functions, creating di erent kinds of interactions between implicit geometries. Indeed, depending on the operator used, the di erent elds can remain independent or be mixed with various degrees. This o ers an easy way to avoid and control the intersections of the geometric representations for meta-bers.

IV.4.1 Applying implicit surfaces to tractograms

Following the previous approaches we used in Section IV.3, we want to de ne the bers as an implicit distance function based on a skeleton: the ber polyline. This can be done using cylinders or capsules on each segment and combining them. The segments will not connect very well with such an approach. For a better representation, we use some Cauchy kernels, allowing for a continuous distance function on the whole polyline. Let [p a , p b ] be a non degenerated segment, and p a point in space where we want to know the eld value. Let:

               = ||p b -p a || = ||p -p a || ℎ = 1 < p -p a , p b -p a > 2 = 1 + 2 ( 2 -ℎ 2 )
with a parameter controlling the decrement radius. The analytical expression of the Cauchy kernel integrated on the segment is de ned as [START_REF] Mccormack | Creating and rendering convolution surfaces[END_REF]]:

[p a ,p b ] (p) = 1 2 2 ( ℎ 2 + 2 ℎ 2 + -ℎ 2 ( -ℎ) 2 + 2 ) + 1 2 3 ( ( ℎ ) + ( ( -ℎ) )) (IV.15)
This expression corresponds to the convolution of a kernel function (H ( ) = 1 (1+ 2 2 ) 2 with the distance to the arbitrary point of evaluation) and a distance function to the segment.

To obtain the eld function for a ber , we add the Cauchy kernels for all its -1 segments:

(p) = -1 =1 [p i ,p i+1 ] (p) (IV.16)
Then, the overall function for the whole bundle or tractogram can be obtained by combining the functions. Here, we choose to add those functions together, resulting in a merging of all bers that are in the same vicinity.

IV.4.2 Visualizing the implicit meta-bers

Di erent approaches can be used to visualize implicit surfaces. This visualization always needs the evaluation of the eld function in order to nd the isovalue at which we want to create the surface.

One possible approach is to create a mesh corresponding to the isovalue of the implicit surface. Two main methods exist to achieve this goal: Marching Cubes [START_REF] Lorensen | Marching cubes: A high resolution 3d surface construction algorithm[END_REF] and Dual Contouring [START_REF] Ju | Dual contouring of hermite data[END_REF]. We used the Dual Contouring with FastAPSS (Section IV.2) because it allows for a higher quality mesh construction compared to Marching Cubes, and we wanted to obtain a mesh as the objective of the technique is to create a surface corresponding to the input point set. When visualizing a tractogram using FastAPSS or Cauchy kernels, we could use this technique, however, it requires an octree (or a grid) thin enough to retrieve all details. This is an issue we have when using FastAPSS with bers (Section IV.3). The precision of the mesh depends heavily on the resolution of this octree/grid, resulting in a massive mesh for a whole tractogram. This is the reason why, for the use of Cauchy kernels, we experiment with another method called ray marching [Tuy andTuy, 1984, Perlin andHo ert, 1989]. Its principle is similar to ray tracing, in the sense that we start by launching rays from a camera and nd intersections with the geometry. At each step of a ray, the eld function (Equation IV.16) is evaluated, leading to either another step parametrized by this evaluation or evaluation of the color with lighting in case the ray has reached the surface. As long as the stepping process is ne enough, this leads to a more precise representation than the reconstruction of a mesh. This comes at the expense of the complexity of the algorithm, which can require a long time for a single image computation as many eld evaluations are required. We use this ray marching approach, based on https://www.shadertoy.com/view/ Xds3zN.

We visualize three di erent bundles using Cauchy kernels added for all bers and ray marching. Although it is not a high number of bers, this gure took a few hours to be computed because of the complexity of the Cauchy kernel, combined with the ine ciency of the basic ray marching approach we use. However, the result is very promising, with elegant mergings between bers occurring.

IV.4.3 Perspectives

Some issues that can happen when using implicit surfaces are linked to the di erence in sizes between segments or bers. Indeed, the implicit functions will behave di erently based on the size di erence of the elements, which does not correspond to our intent, as we would prefer segments of di erent sizes to behave in the same way. The Scalis approach [START_REF] Zanni | Scale-invariant integral surfaces[END_REF] prevents this from happening by working in a di erent space where the scale does not matter. Our bers would bene t from such an approach when using tractograms with segments of di erent lengths for instance.

As we previously mentioned, implicit functions can be merged with di erent operators, with various e ects. This could be used in conjunction with the multiresolution representation (Section III.2), using di erent operators for bers that are part of a same meta-ber or between di erent meta-bers at a given level of resolution. In this case, it means that the only requirement from the multiresolution representation is an ordered list of progressive mergings of the individual bers as the geometry will be implicitly built based on these merging operators.

One major limitation of our implementation is the complexity and the time required to obtain a single image. This could be improved in several ways, among which is the possibility to use more intelligent ray marching approaches, that would try to maximize the steps of the ray marching while guaranteeing that the surface is reached, like for instance this recent article [START_REF] Galin | Segment tracing using local lipschitz bounds[END_REF]. Other possibilities would be to try reducing the number of Cauchy kernels computations required, by for instance removing bers far behind the visualized model, that are not participating in the nal rendering.

IV.4.4 Summary

This quick exploration of implicit surfaces for visualizing tractograms is showing that it is indeed an interesting representation. It solves most issues exposed in Section IV.1.2 with generalized cylinders, at the expense of the complexity. This complexity could be lowered by rst using state-of-the-art techniques, and then by investigating some new simpli cations adapted to tractograms.

IV.5 Conclusion

IV.5.1 Summary

In this chapter, we investigated di erent methods for representing bers as well as sets of bers that we call meta-bers. First, we experimented using lines and generalized cylinders with elliptical bases. These representations have the advantage of being quite easy to implement, and o er good visualization performances, with lines being better in this regard. However, it lacks precision in the diversity of con gurations it can represent. For instance, a meta-bers with two regions on one extremity cannot be represented with such a geometry. This is why we also explored algebraic point set surfaces (APSS). The global approach was interesting for us but su ered from its complexity. Indeed, the evaluation of the sphere on which to project any point q in space requires to consider all input points with a weight depending on the distance to q, even though this in uence is limited for faraway points. We take advantage of this by grouping points in an octree and considering the average point per leaf when su ciently far away from the leaf.

The results obtained with this simpli cation, that we call fastAPSS, are extremely close to the standard version, with reduced complexity. We then experimented the use of fastAPSS on brain tractograms, showing that the geometric quality obtained is promising. One limitation remains, however: the intersections between the geometries of di erent meta-bers. It can be solved using implicit surfaces, at the expense of either creating a mesh for visualization or using some ray launching techniques, both being high complexity approaches.

IV.5.2 Limitations

We experimented di erent geometric representations for brain tractograms, with each its own advantages and drawbacks. Some improvements remain to be done for fastAPSS to better work with bers without the exposed issues. Other kernels could be explored, as well as some speci c handling of points, for determining with better accuracy the points that we need to keep for the new meta-bers. Solving these issues will then allow further experiments on the proposed volume aware similarity.

Our implicit surfaces visualization would also bene t from some optimizations and the reduction of the complexity with the use of more advanced techniques. On top of that, one interesting work would be to use di erent geometric representations depending on the objective. For instance, fastAPSS meta-bers could be used to rst compute a merging order for the multi-resolution representation. Then, real-time visualization of this order could be performed with generalized cylinders for their lower rendering cost. The users can then choose a level of resolution, on which a better visualization could be performed with the implicit surfaces, using di erent combination strategies depending on a ber belonging to a meta-ber or not. This approach would bring the best of each tested method, and try to minimize their drawbacks.

interpretation. We illustrate our method using the HCP dataset [Van Essen et al., 2012] and compare it to other well-known white matter segmentation techniques.

This segmentation was presented at OHBM by Alessandro Delmonte [START_REF] Delmonte | Segmentation of white matter tractograms using fuzzy spatial relations[END_REF]] and accepted at ISBI as a conference paper [START_REF] Delmonte | White matter multi-resolution segmentation using fuzzy set theory[END_REF] where it was also presented as a poster. This is a joint work with Alessandro Delmonte.

V.1.2 Method

V.1.2.1 Modeling spatial relations

Qualitative anatomical de nitions of white matter tracts are mostly based on spatial relations with respect to anatomical structures (e.g. "anterior to the hippocampus").

Translating these de nitions into operational algorithms calls for mathematical models representing the intrinsic vagueness of the de nitions. To this end, we rely on the fuzzy set theory and the spatial relation models proposed in [Bloch, 2005]. The general idea is to de ne for each point in space the degree to which it satis es a given relation with respect to a reference object (i.e. an anatomical structure). As an example, Figure V.1 presents a fuzzy cone, corresponding to the relation "anterior to the hippocampus". This is formalized next. V.1.2.1.1 Directionality First, we explain how we model directional relations, such as "in direction with respect to a reference structure ". Let ì be a unit vector pointing in direction . We de ne a fuzzy spatial region as a membership nal fuzzy region satisfying the combination of relations. Every point P in space (i.e. every voxel) is thus assigned a membership value * describing the degree of satisfaction of the combined relations. Every ber of a tractogram is then assigned a fuzzy score FS which is computed as the weighted average of the membership values * of the voxels through which the ber passes. The fuzzy score FS of a ber is computed as:

( ) = 1 ℎ( ) ∈ * ( ) × ℎ( ∩ ) (V.3) with = { | ∩ ≠ ∅}
being the set of voxels traversed by the ber . The score weight is computed as the proportion of the length of the ber within all voxels of ( ℎ( ∩ )).

V.1.2.1.2 Connectivity In addition to the relative directions, we also model another common anatomical de nition about the location of the tract terminations (e.g. " bers terminate in temporal lobe"). Let f a be one of the endpoints of a ber and the region of the ending region (Figure V.3), we de ne the degree of connectivity as:

(f a ) = max r∈ ( - ||f a -r|| 2 2 2 ). (V.4)
where is a parameter that can be tuned to adapt to the considered region. If ∃r ∈ such that f a = r (i. e., f a ∈ ), then (f a ) = 1 , and (f a ) decreases when the distance from f a to increases. When the de nition involves only one ending region , f a is the closest endpoint to . Otherwise, when using two ending regions, the ber orientation is the one minimizing the sum of the distances between the ending points and the regions (each extremity being linked with a di erent region). Let 1 and 2 be two regions.

If (f a , 1 ) + (f b , 2 ) ≤ (f a , 2 ) + (f b , 1 ), then = 1 and = 2 . Otherwise, = 2 and = 1 . Then, we compute an endpoint score for the extremities f a and f b of the ber:

( ) = (f a ) × (f b ) (V.5)
In case only one extremity f a is considered, Equation V.5 becomes ( ) = (f a ). Note that, in this work, bers have been cut at the boundary between white matter and gray matter.

V.1.2.1.3 Anatomical coherence score Eventually, all relations are combined in a conjunctive way, resulting in one score per ber, that we name the Anatomical Coherence Score (ACS):

( ) = ( ) × ( ) (V.6)
This score corresponds to the level at which a ber is conforming to the de nitions used. Its value is included in [0, 1] by construction, with 1 meaning that the ber completely satis es the description, and 0 meaning that it is not at all part of the tract of interest.

V.1.2.2 Multi-resolution representation

The segmentation process involves choosing the ACS threshold to select only bers with higher ACS values. It is time-consuming and can be quite hard, especially when considering a whole-brain tractogram composed of millions of bers. This is particularly the case for applications demanding high accuracy such as surgical planning. For instance, Figure V.4 presents a half brain whose bers have been colored according to the ACS values for the UF.

We take advantage of the previously introduced simpli cation method (Section III.2).

Based on this model, we create a GUI where the user can navigate in real-time through di erent levels of detail and display the bers/cylinders with an ACS higher than a selected threshold. As it was exposed in Section III.2.2.1, the multiresolution approach can use any dissimilarity metric. Consequently, we propose to use an extension of the computational model of Weighted Currents ( ) [START_REF] Gori | Parsimonious approximation of streamline trajectories in white matter ber bundles[END_REF], that will use this ACS to improve the merging order based on the tract of interest. Let = {c i , i } and = {d j , j } be two bers composed of and segments respectively, where c i and i (respectively d j and j ) are the centers and tangent vectors of (respectively ). The proposed dissimilarity measure is de ned as: The proposed technique simpli es the geometric representation, preserving at the same time the overall structure of the original tractogram (i.e. shape, connectivity, and ACS). The same Delaunay tetrahedralization of bers extremities is used to reduce the complexity. This provides a geometrically well-distributed set of adjacency links. The similarity measure in Equation V.7 is only computed between adjacent (linked) bers. Fuzzy scores for cylinders are computed as the weighted average of the membership values * of the voxels the cylinders cross. As for bers, weights are computed as the proportion of the length of the cylinder inside the voxels:

( , ) = (|ACS( ) -ACS( )|) ( f a -t a 2 ) ( f b -t b 2 ) -1 =1 -1 =1 i ( c i -d j 2 ) j (V.7)
( ) = 1 ( ) ∈ * ( ) × ( ∩ ) (V.8) with = { | ∩ ≠ ∅}
being the set of voxels traversed by the cylinder . In practice, this computation is done by considering a pattern of segments in each volume section of the cylinders, computing their average value, and using the length of the center-line of the cylinder as a proportion of the overall length of the cylinder. The distance between cylinder terminations and ending regions is computed as the average distance between the ending surface of the cylinder and the region . We also propose an automatic stopping criterion for the multiresolution to prevent oversimpli cation (e.g. a single cylinder). We use the inner product of Equation V.7 to compute angles between bers and cylinder center-lines. Two cylinders/ bers are not merged if they are almost orthogonal. Fibers that were never merged in the whole process are then considered as outliers and discarded.

V.1.3 Results

V.1.3.1 Dataset and Numerical aspects

Experiments are based on 10 Uncinate Fasciculus (UF) and 10 Inferior Fronto-Occipital Fasciculus (IFOF) segmented from 5 healthy subjects of the HCP dataset [Van Essen et al., 2012]. Gray matter structures are extracted using FreeSurfer.1 Whole-brain tractograms of one million bers are estimated with MRTrix3 [START_REF] Tournier | MRtrix: Di usion tractography in crossing ber regions[END_REF] using a probabilistic algorithm (iFOD2) and the Constrained Spherical Deconvolution model. We separately analyze the two brain hemispheres, removing bers measuring less than 4cm. In Equation V.7, we used the value of 6 for , and . The de nitions of fuzzy relations are implemented in Python and the computational time is about 10-15s per de nition. Multi-resolution is instead implemented in C++, Qt and OpenGL. The computational time depends on the number of bers with a complexity of ( ( )) (Section III.2.3.3).

V.1.3.2 Tract modeling

Based on both imaging and dissection studies [START_REF] Ebeling | Topography of the uncinate fascicle and adjacent temporal ber tracts[END_REF], Wakana et al., 2007, Catani and Thiebaut de Schotten, 2008, Sarubbo et al., 2013] and with the help of an experienced neurosurgeon, we de ne the Uncinate Fasciculus bundle (UF) and the Inferior Fronto Occupital Fasciculus bundle (IFOF) considering the labels available in FreeSurfer. The employed de nition for the UF is: inferior of putamen AND inferior of claustrum AND anterior of amygdala. Moreover, one endpoint of the bers should be in the anterior temporal lobe. For the IFOF, we combine together the frontal, occipital and stem de nitions as: internal AND lateral of external capsule AND anterior of hippocampus OR anterior AND superior AND lateral of putamen OR posterior of putamen AND posterior AND superior of hippocampus. Furthermore, bers should have one endpoint in the occipital lobe and the other one in the orbital and inferior frontal gyrus.

We also propose to model the occipital lobe using the following rule: the cuneus, lateral-occipital section and pericalcarine regions are all considered part of the occipital lobe. Then, we nd a plane perpendicular to the principal direction of the hippocampus and we move it along the posterior direction until it reaches one of the three previous regions. All voxels composing the fusiform or lingual label and lying posterior to the plane de ne the complementary part of the occipital lobe.

These anatomical regions are illustrated in Figure V.6.

V.1.3.3 Experiments

In Figure V.7, we show the segmentation of the IFOF and UF for three subjects, with three di erent ACS thresholds. It can be noticed that results are very similar among subjects even when varying the thresholds. Reproducibility has also been observed at other thresholds.

In Figure V.8, di erent levels of resolution are shown for both IFOF and UF, using two thresholds for the ACS. The proposed parsimonious model based on generalized cylinders eases the interpretation of a tract by giving a summarized representation of the original bers. The structure of the original bundle is more visible and recognizable at a lower resolution (2nd and 3rd columns). Furthermore, the proposed GUI allows the user to navigate in real-time through di erent levels of resolution and display bers/cylinders with ACS higher than a selected threshold. Using it, our clinical colleagues were able to analyze the structures of the tracts at di erent resolutions and thresholds and choose an optimal one (ACS2 in In Figure V.9, we compare the results of the proposed method with the ones obtained using a manual segmentation, delineated according to [Wakana et al., 2007, Catani andThiebaut de Schotten, 2008], and the default WMQL queries [START_REF] Wassermann | The white matter query language: a novel approach for describing human white matter anatomy[END_REF]. These methods show a greater dispersion in the pathways of the bers, and thus less reproducible results, compared to the proposed method.

The percentage of white-matter voxels containing the IFOF bers of at least one subject decreases from the 15% of the manual segmentation, to 9% in WMQL, to the achieved 6%. This is mainly due to the fact that our method produces more compact segmentations where there are fewer voxels with bers belonging to only one subject.

Note that, for both hemispheres, the computational time for the IFOF (resp. UF) was ∼100s (resp. ∼25s) for the proposed method, ∼300s (resp. ∼240s) for WMQL and ∼900s (resp. ∼900s) for the manual segmentation.

V.1.4 Discussion

V.1.4.1 Summary

We proposed a segmentation algorithm for white matter tractograms that produces anatomically-relevant and reproducible tracts. Di erently from most segmentation algorithms, we directly model the qualitative anatomical de nitions of the tracts using the fuzzy set theory. This makes it possible to cope with the inherent imprecision of the anatomical de nitions and to combine them into a single quantitative value (ACS) which summarizes the degree of satisfaction of all relations. To ease interpretation when working with millions of bers, we combined this method with the multi-resolution representation from Section III.2. It produces a multi-resolution representation and provides a GUI where the user can navigate in real-time through di erent levels of resolution. This helped our clinical colleagues to better understand the structure of the bundles and nd an optimal ACS threshold Figure V.9 -Reproducibility analysis carried on ve subjects from the HCP dataset.

A colormap encodes, for each voxel of the image, the number of subjects that have at least one of the segmented IFOF bers passing through this voxel. The areas showing a smaller ber dispersion with our method are highlighted using a white arrow.

for the segmentation of IFOF and UF. Using 5 subjects of the HCP dataset, we showed that our results were reproducible and less scattered than using a manual segmentation or WMQL [START_REF] Wassermann | The white matter query language: a novel approach for describing human white matter anatomy[END_REF].

V.1.4.2 Limitations

One major limitation of this approach is the need for a segmentation of brain regions to be able to create the de nitions of a ber bundle in anatomical terms. This is however also the case for WMQL, and manual segmentation requires even more knowledge. Another issue is the de nitions themselves, for which there is not always a consensus. Indeed, even neurosurgeons do not always agree on the correct relations to use. This is the reason why we had to propose some speci c relations for the two bundles we segmented. An extension of the proposed technique to other tracts would be a great contribution, as well as its use in statistical analyses, taking advantage of its reproducibility. Also, a study on the e ect of the di erent parts of the ACS (Equation V.6) would be interesting to understand the impact of every de nition on the nal segmented bundle.

V.2 Compression V.2.1 Introduction

Di usion MRI ber tracking datasets can contain millions of 3D bers, and their representation can weight tens of gigabytes of memory. Their size makes them di cult to store, visualize, process, or exchange over the network. There is an important need for a new e cient and fast compression algorithm for brain tractography data. We propose a new compression algorithm well-suited for tractograms, by taking advantage of the way bers are obtained with usual tracking algorithms.

Our approach is based on unit vector quantization methods combined with a spatial transformation which results in low compression and decompression times, as well as a high compression ratio. For instance, an 11.5GB tractogram can be compressed to a 1.02GB le and decompressed in 11.3 seconds on an Intel Xeon E5-1650v4 processor (6cores, 12 threads, 3.6GHz). Moreover, our method allows for the compression and decompression of individual bers, reducing the need for a costly out-of-core algorithm with heavy datasets. Last, we open a way towards on-the-y compression and decompression for handling larger datasets without needing a load of RAM (i.e. in-core handling), for faster network exchanges and faster loading times for visualization or processing. To combine speed, compression ratio, the ability to compress independently each ber of the tractogram, and the scalability, we focus our work on the representation of individual bers and introduce a new compression and decompression algorithm based on this representation. This method has been presented at OHBM [START_REF] Rousseau | Q b: Fast and accurate compression of white matter tractograms[END_REF] and published in NeuroInformatics [START_REF] Mercier | Q b: Fast and e cient brain tractogram compression[END_REF]. It is a joint work with Sylvain Rousseau.

V.2.1.1 Data constraints for our representation

To compress data, it is necessary to nd a space of representation smaller than the original one. In our case, we take advantage of the way bers are acquired. In particular, we use two constraints: a constant stepsize ( ) and a maximum angle of deviation ( ).

Commonly used tractography algorithms use a constant stepsize while tracking bers. This property is sometimes used in some post-processing [START_REF] Soares | A hitchhiker's guide to di usion tensor imaging[END_REF]. Furthermore, termination of bers happens when the angle between two consecutive segments of the ber is higher than a maximum angle (see Section III.1.1). Using these two constraints allows us to consider a ber as a succession of unit vectors whose angles between consecutive vectors are bounded by the maximum angle.

103

V.2.1.2 Unit vectors quantization

When considering lossy compression of independent unit vectors, the literature has shown that unit vectors quantizations are the most e cient [START_REF] Cigolle | A Survey of E cient Representations for Independent Unit Vectors[END_REF]. They consider unit vectors as points on the surface of the unit sphere. A unit vector is compressed (or quantized) by being approximated with the closest point of a point set distributed onto the surface of the unit sphere, as shown in Figure V.10. Each point is then encoded using an identi er that is smaller in memory than its original 3D coordinates. To minimize the maximal error (angle between the original unit vector and the quantized unit vector) the chosen point set needs to have a distribution as uniform as possible onto the surface of the unit sphere. Among possible distributions, octahedral quantization [START_REF] Meyer | On Floating-Point Normal Vectors[END_REF]] is a good candidate for a fast compression and decompression of the unit vectors [START_REF] Cigolle | A Survey of E cient Representations for Independent Unit Vectors[END_REF]. More recently, a new inverse mapping for the spherical Fibonacci point set was introduced by [START_REF] Keinert | Spherical Fibonacci Mapping[END_REF]. This point set is known for the uniformity of its distribution over the surface of the unit sphere. This quantization will, therefore, decrease the quantization error, but at the cost of computational complexity. Both distributions make sense with our approach as we can either want a fast decompression for streaming algorithms (octahedral) or sacri ce a little bit of speed (spherical Fibonacci) to decrease the compression error that is directly linked to the compression ratio. These two methods will be used in the remainder of this chapter. The common point of all unit vector quantization techniques is their speed. These algorithms can usually decompress millions of unit vectors per second on a single processor [START_REF] Cigolle | A Survey of E cient Representations for Independent Unit Vectors[END_REF] and billions of them on a graphics processor unit (GPU).

V.2.1.3 Contribution

We introduce a new representation model for the bers of brain tractograms. We then propose qfib, a new compression and decompression algorithm based on the proposed representation. It handles each ber individually using two speci c constraints of tractography methods (constant stepsize and maximum angle between consecutive points), a mapping on the unit sphere, and a unit vector quantization. This method has the following properties:

• high compression ratio (between 80 and 90%);

• fast compression and decompression (few seconds);

• low compression error;

• ability to compress/decompress individual bers;

• conservation of the number of points of each ber;

• ability to access any random single ber from the compressed dataset.

Moreover, we provide a publicly available open-source implementation of our compression and decompression algorithm along with a docker le2 .

V.2.2 Method

We de ne as the constant stepsize and as the maximum angle in the following parts of this section.

V.2.2.1 Brain bers representation

A ber containing points is described as a 3D polyline = {p 1 . . . p N } of ordered points. Each ber has its own number of points as we use a constant stepsize . The ℎ point of the ber is de ned as

p i = p i-1 + ----→ -1 for ∈ [2, ].
Since two successive points along the ber are supposed to be at a constant distance, p i can also be written as:

p i = p i-1 + - ---→ -1 (V.9)
where -→ denotes the normalized vector -→ /|| -→ ||. This means that each 3D polyline from the dataset can be represented by its rst point and a set of unit vectors. The stepsize needs to be constant on a per-ber basis for this representation to work.

V.2.2.2 Unit vectors quantization

Using unit vectors already reduces the size required to store bers, as only two dimensions (spherical coordinates) are needed to encode each segment instead of the three required by the Cartesian coordinates. When encoding the rst two points of each ber, their size, and unit vectors, the compression ratio will be slightly under 33%, depending on the average number of points in the ber of the dataset. In order to further increase the compression ratio, we take advantage of the tractography resolution -linked to the dMRI resolution -and use a lossy compression. As introduced in Section V.1.1, when choosing a quantization method, we can prioritize the speed with the octahedral quantization [START_REF] Meyer | On Floating-Point Normal Vectors[END_REF], or the precision with the spherical Fibonacci quantization [START_REF] Keinert | Spherical Fibonacci Mapping[END_REF].

V.2.2.2.1 Octahedral quantization [START_REF] Meyer | On Floating-Point Normal Vectors[END_REF] It is a unit vector representation method that projects the vector [ , , ] de ned on the surface of the unit sphere to an octahedron by normalizing it using an L1-norm. This octahedron is then unwrapped onto a 2D unit square as shown in Figure V.11. The resulting 2D coordinates [ , ] are discretized prior to encoding. When encoding using an bits quantization, each of these two coordinates is discretized on /2 bits. The unquantized vector can be retrieved using the following formula [START_REF] Meyer | On Floating-Point Normal Vectors[END_REF]: [START_REF] Meyer | On Floating-Point Normal Vectors[END_REF] projects a unit vector to an octahedron, then, to a unit square to encode the discretization of the resulting 2D coordinates.

= 1 -| | -| | [ , ] = [ , ] if ≥ 0 [ ( ) -, ( ) -] if < 0 (V.
V.2.2.2.2 Spherical Fibonacci quantization [START_REF] Keinert | Spherical Fibonacci Mapping[END_REF] This method is based on the Spherical Fibonacci point set as it yields a nearly uniform point distribution on the surface of the unit sphere [START_REF] Keinert | Spherical Fibonacci Mapping[END_REF]. Using this point set, the spherical coordinates ( , ) of the ℎ point of the point set containing points are de ned as [González, 2010]:

= (1 -2 +1 ) = 2 ((3 - √ 5)/2) (V.11)
When compressing a vector using an bits quantization, the number of points of discretization used will be = 2 . The mapping proposed in [START_REF] Keinert | Spherical Fibonacci Mapping[END_REF] is used to nd the closest spherical Fibonacci point to the unit vector to encode (Figure V.10) and the identi er of this closest point is stored. The advantages of both quantization methods are discussed in Section V.2.3.

V.2.2.3 Mapping

We showed in Section V.2.2.2 that quantized unit vectors can be used to represent bers. This representation can be considered as a di erential description of the ber for which the decompression is performed by an integration using the rst point and the stepsize as constant terms. It provides good compression ratios, but it is possible to even push forward the compression by adding the constraint on the maximal angle between two consecutive segments. We take advantage of this property by expressing each unit vector with respect to the previous one. It can be seen as a second-order di erential representation (second derivative). The relative position between consecutive unit vectors is limited to a spherical cap, parametrized by the angle (in yellow in Figure V.12). As such, the quantization point set would only be partially used, which is not optimal. De ning a uniform point set on a spherical cap could solve this issue, however, none exists [START_REF] Rousseau | Fast Lossy Compression of 3D Unit Vector Sets[END_REF].

Another possibility is to de ne a mapping from the spherical cap to the whole unit sphere. In this way, vectors could be quantized on the whole unit sphere. To be retrieved, vectors will then need to be unquantized before being unmapped. To minimize the error, this mapping should preserve the uniformity of the point set distribution. Such a mapping was recently introduced by [START_REF] Rousseau | Fast Lossy Compression of 3D Unit Vector Sets[END_REF]. As illustrated in Figure V.12, the conservation of the point set distribution is ensured by keeping the equality between the ratios 1 / and 2 / ℎ where 1 is an in nitesimal ring on the surface of the spherical cap (orange), 2 is the in nitesimal ring on the whole sphere corresponding to 1 mapped to the whole sphere (blue), and and ℎ are respectively the surface of the spherical cap (yellow) and of the whole unit sphere (light blue). This method can be adapted to our problem by encoding each unit vector from its predecessor. A point x on a spherical cap is transformed into a point x on the whole unit sphere by:

x = P 0 + √ 1 -2 ℎ = 1 - 1 - ( 1 ) (V.12) where = 1-( ) 2
+ is the ratio between the length of the projection of the spherical cap on axis P 0 and the diameter of the unit sphere (Figure V.12). In practice, the unit vector ----→ -1 is mapped relative to - -----→ -2 -1 to nd the vector x , by using 0 = ------→ -2 -1 and x = ----→ -1 . As indicated in the original article [START_REF] Rousseau | Fast Lossy Compression of 3D Unit Vector Sets[END_REF], due to numerical imprecision, we add a small value to the ratio to ensure that all unit vectors lie on the spherical cap. Representing unit vectors relative to the previous ones can be seen as the second derivative, as the rst point is needed for the rst integration and the rst unit vector for the second one. and with a maximal angle , represented as the yellow area on the surface of the unit sphere, to another point x , de ned on the surface of the whole unit sphere, colored in light blue.

V.2.2.4 Propagation error reduction

Each time a unit vector is compressed, a small error is introduced due to the quantization process, meaning that starting from the point p i-1 , the new position will be p i instead of p i . Therefore, continuously compressing the unit vectors ----→ -1

can lead to an important error as shown with the orange curve in Figure V.13.

To avoid this error from being fully propagated, we substitute the points p i-1 for p i-1 . This point (p i-1 ) has already been mapped and quantized (compressed) at the previous step, and then unquantized and unmapped (decompressed). As a result, it is not at the exact same place as p i-1 , and instead of compressing the unit vectors ----→ -1 , we compress the unit vectors - ---→ -1 . This slight modi cation compensates for the accumulated error. Indeed, it allows the ber being compressed to always use a direction related to its current position, thus reducing the error propagation. This is presented on the green curve of Figure V.13. In Figure V.14, we show on a small ber example the di erence between a ber compressed without error propagation reduction (in orange) and with our error propagation reduction strategy (in green). The second one is closer to the original ber (in black).

V.2.2.5 Compression algorithm and encoding pattern

To compress the ber = {p 1 , ..., p N } into a ber for which the slightly changed decompressed coordinates will be {p 1 , ..., p N }, combining all the steps described in Section V.2.2, we obtain the following algorithm:

• we store the rst and second points as they are, using their Cartesian coordinates;

• for the other points, we use the mapping proposed in [START_REF] Rousseau | Fast Lossy Compression of 3D Unit Vector Sets[END_REF] on the unit vector - ---→ -1 where the axis of the spherical cap is - -----→ -2 -1 , and the ratio depends on the maximum angle ;

• we quantize the mapped vector using the octahedral [START_REF] Meyer | On Floating-Point Normal Vectors[END_REF] or the spherical Fibonacci quantization [START_REF] Keinert | Spherical Fibonacci Mapping[END_REF] depending on the We rst illustrate the results of our method on a few toy examples, then compare our algorithm with zfib [START_REF] Presseau | A new compression format for ber tracking datasets[END_REF] (see Section II.5 for details on their approach). We evaluate the compression ratios of the two approaches and their impact on mapped scalar values, using fractional anisotropy (FA) as an example.

We conclude by showing that our method can also run out-of-core.

V.2.3.1 Error

To better illustrate the error introduced using our compression method, we show in Figure V.16 four di erent bers (in black) and their compressed and decompressed version (in green). As we quantize the relative directions of each segment, a straight ber (Figure V.16(a)) remains totally straight, with no error introduced. However, when introducing some curvature in the bers, we notice that the error increases according to the angle as shown in Figures V.16(b), V.16(c), V.16(d), for which the maximum angles are respectively 40.3°, 45.7°and 73.7°. With the spiral curve (Figure V.16(d)), the angle is too important for the mapping to be e ective. In that case, the error is higher, and the di erence with the original curve is more visible than with the other curves. However, this ber is not anatomically meaningful as such a high curvature should not appear in a tractogram.

Our compression keeps the original number of points. As a result, we are able to measure pointwise errors. The errors are computed between the original dataset containing the bers = { ,1 . . . , , } and the compressed and decompressed bers = { ,1 , . . . , , } using the following formulas:

= =1 =1 || , - , || =1 (V.13) = {|| , - , || | ∈ [1, ], ∈ [1, ]}
where is the number of points in the ber , and the number of bers in the dataset. Table V.2 shows these errors for every le described in Table V.1. It emphasizes that increasing the precision of the quantization from 8 bits to 16 bits greatly reduces the error as there are more possible points on the unit sphere during quantization. However, this results in lower compression ratios (see Section V.2.3.2). Even though it is usually advised to use stepsizes equal to a tenth of the voxel size (so around 0.1 or 0.2 in our case), greater stepsizes can be used. In these cases, according to our errors, it is better to use a 16 bits quantization. The method used for quantization impacts the error. The spherical Fibonacci quantization introduces less error, but octahedral quantization remains interesting for its faster computation time with still relatively low error (see Section V.2.3.3). Figure V.17 evaluates the error depending on bers length on a histogram where each column represents a cluster of 10 length di erence. The tractograms were made using a probabilistic algorithm and a 0.1 stepsize. They were compressed and decompressed using an 8 bits octahedral quantization. The histogram shows that the maximum error seems to be mostly independent of the length, in contrast to the average error which increases with the length of the bers. Table V.3 shows the maximum and average errors obtained at the extremities of the bers. This point is important as it characterizes the connectivity of the ber. We notice that the error (including the maximum one) is really low for ≤ 0.2 with an 8 bits quantization. For = 0.5 , the error starts to be too important (same order of magnitude than the voxel size), especially with the octahedral quantization (with maximum errors of 0.5

), and a 16 bits quantization is preferable.

To visualize these errors, we color-mapped the compressed bers according to their distance from the original dataset (Figure V.18). The quantization used is the octahedral one. The color scale is going from 0 (dark blue) to 1.25 (dark red) of error, the voxel size of our data. We used the bers computed with a 1 stepsize using the deterministic method. This case is an illustration of the situations where it is important to use a 16 bits precision quantization. Indeed, the error obtained in this case (Figure V.18(b)) is almost not visible on the picture, compared to the 8 bits precision (Figure V.18(a)).

V.2.3.2 Compression Ratio

The compression ratio is computed using the following formula:

= 100 × (1 - ) (V.14)
Our compression ratio (Table V.4) is mostly stable around 90% using an 8 bits quantization and around 82% using a 16 bits quantization. We emphasize that our compression ratios are strictly better than zfib for the same error values with an 8 bits quantization and stepsizes under 0.5 , which correspond to the advised values for the stepsize. When the zfib error is set to 0.2 (their default value), we have either better compression ratios or smaller errors. For the bers generated using a probabilistic method, that are more tangled, we achieve better compression ratios with the same error value in every shown case.

V.2.3.3 Computation time

In Table V.5, we show the compression and decompression times of qfib and zfib. These gures were obtained using a computer with an Intel Xeon E5-1650v4 (6 same ★ : same error as qfib when using an 8 bit octahedral quantization (Table .V.2).

cores, 12 threads, 3.6GHz). In both methods, we only account for the compression and decompression time and not reading and writing from/to the hard drive. We set the error of zfib to 0.2 mm as it is the default value on the provided source code. We can see that, for both compression and decompression, our method is at least two orders of magnitude faster than zfib. Moreover, during our experiments using the source code provided by the authors, zfib failed to compress very large les whereas our method is more scalable.

V.2.3.4 Impact on Fractional Anisotropy

Scalar quantities such as Fractional Anisotropy (FA) are often mapped onto the bers. As such, it is necessary to verify that their values remain respected even when bers are compressed. We compute the FA map using a Bresenham-like integration similar to [START_REF] Presseau | A new compression format for ber tracking datasets[END_REF]. Results are presented in Table V.6 as an error percentage between the original bers and the compressed and decompressed ones. The error used for Presseau's algorithm was the default one of 0.2 . This table shows that we obtain a lower average FA error in all tested con gurations.

To go even further, we compute the FA error pointwise similarly to Equation V.13. We use a trilinear interpolation of the 8 voxels surrounding each point. This cannot be done for zfib as it does not keep the same amount of points as the original dataset. The result of this comparison is presented in Table V.7. We can see that our maximum errors are under 0.1 (so a 10% error) for 8 bits quantization when using stepsizes less than 0.5 . Moreover, the average error is really low, which was expected considering the low errors reported in Table V.2.

V.2.3.5 Out-of-Core version

To push forward the scalability of our method, we developed an out-of-core version of our algorithm. Compression times are presented in Table V.8. They are higher Table V.6 -Comparison of zfib and qfib for FA computation. We compute the average error in FA computation using a Bresenham-like integration [START_REF] Presseau | A new compression format for ber tracking datasets[END_REF]. than the in-core version. This is explained by the need to use the hard drive during compression and decompression in contrast to Table V.5. This algorithm requires a few MB of RAM to run, no matter the tractogram size, making it possible to be used on any computer. To demonstrate that our algorithm is not limited to powerful hardware or small les, we decided to compress a huge dataset using a single-board computer. We generated a dataset containing 10 millions of bers, with a stepsize of 0.1 using the probabilistic method, resulting in an 87.3GB tck le. For the hardware part, we used a Raspberry Pi 2 (≈ 35$) with 1GB of RAM, and a processor with 4 cores at 900 MHz. The data were stored in an external hard drive connected to the Raspberry Pi using the USB 2.0 interface. The compression resulted in a 7.49GB le with a compression ratio of 91.4% using an 8 bits octahedral quantization. It took around 4 hours, whereas the decompression took around 5.3 hours. These timings are consequent because of the USB 2.0 interface limitation. An important amount of time was indeed required to read and write all the data. This explains why the decompression took more time than the compression, as writing is slower than reading on a hard drive. The main limitation of our method is the requirement for a constant stepsize per ber. To overcome this issue, one could previously use a resampling method. Furthermore, we conceived our method especially for brain tractography, using speci c properties and constraints, thus other limitations could appear when trying to apply the same compression algorithm to other kinds of data. For instance, in the case of bers with a bigger maximum angle between consecutive segments, the mapping might not be as e cient, as illustrated in Figure V.16. In such a case, it could be necessary to use a 16 bit or higher precision quantization, thus reducing the resulting compression ratio.

V.2.3.7 Future work

Our method is pleasingly parallel, which means that it can be trivially parallelized at the ber level. Therefore, it would be straightforward to implement a massively parallel version on the GPU. This could improve a lot the compression and decompression times, however, the limiting factor will remain reading and writing the data to a disk. In the future, it could be interesting to integrate the proposed method in existing tractography algorithms to generate bers directly in this compressed format. To further improve the ber compression, a more parsimonious space of representation could be used, using, for instance, Bézier curves, B-splines, or NURBS. Another possibility would be to group bers depending on their maximum angle , so that we could use even fewer bits per point for almost straight bers while keeping a low error as highly curved bers could use 16 bits or higher precision.

V.2.4 Conclusion

We have presented a novel ber compression algorithm for brain tractograms, and its associated encoding format -qfib. We evaluated and validated it using a wide variety of brain tractography con gurations with di erent stepsizes and numbers of bers. This algorithm provides lower errors in general cases and better compression and decompression times than existing methods. The compression ratio is high, around 90%, with errors under the dMRI precision. In contrast to other methods, qfib does not remove any point from the original dataset, which is important when point-based algorithms need to be applied afterward. Moreover, the compression and decompression steps handle each ber individually. This implies that they can both be easily parallelized and that qfib is as e cient on small as on large databases in contrast to dictionary-based methods. This also means that one can decompress single bers from the compressed representation, allowing users to directly access speci c fascicles or tracts with negligible memory usage and loading time. Consequently, our approach opens a way towards on-the-y applications where an algorithm could work with the compressed representation of the bers in memory, and decompress only single bers or small bundles before applying some computation on them. This reduces the need for loads of RAM or VRAM. We also demonstrated that our algorithm works out-of-core. In this case, the size of the data to compress is only limited by the hard drive.

V.3 Visualization

We previously exposed the geometric aspect of visualizing tractograms, and which kind of geometry can be useful in di erent cases. However, we did not yet discuss visualization. Using a good geometry to represent bers is very important for the readability of the data, which can be enhanced by adding some visualization techniques. Another aspect is the possibility o ered by the hardware, to even be able to output an image of a whole huge tractogram. We explore some of these aspects in this section, starting with the visualization of huge datasets using a modi ed version of our compression algorithm qfib presented in Section V.2, before exploiting some computer graphics methods to help with their visualization, such as ambient occlusion.

V.3.1 Visualizing compressed bers

The use of a geometry to represent individual bers as well as meta-bers helps a lot for interpretability, at the cost of some more computation. As a result, it is extremely useful for bundle or simpli ed tractogram studies.

We presented qfib in Section V.2 with its capabilities in terms of compressing bers les to reduce disk space requirements or improve sharing. Verifying the validity of the retrieved bers when there are millions of them is challenging, as using the compressed bers requires to rst decompress the ones we want to use. For comparing bers together in order to obtain bundles, this is usually ne as metrics compare bers in pairs, requiring only these two bers to be decompressed. However, it becomes challenging for the visualization of the entire tractogram as all bers are required in this case, and some tractogram les can easily exceed RAM when decompressed, let alone VRAM! To circumvent this issue, we will take advantage of qfib decompression speed and realize it on-the-y on the GPU, requiring only enough RAM and VRAM to t the compressed bers. This is a joint work with Sylvain Rousseau (for the compression) and Jérémie Schertzer (for the optimizations).

V.3.1.1 Organizing data for GPU processing

Applying decompression directly is a possibility, however, some limitations may appear on the GPU depending on the graphics pipeline used. To overcome this, and also improve parallelism and optimization possibilities, we propose to slightly modify the structure of qfib compressed bers. Instead of compressing each ber individually, we will compress blocks containing at most 64 points. These blocks will weigh up to 80 bytes, decomposed as explained in Table V.9: This number of points will need to be chosen depending on the pipeline as we will explain later. The GPU aligns data at di erent levels for optimization purposes, meaning that the size of the smallest element it is able to access will determine the possible size of the block. The smallest element is an integer represented with 4 bytes, requiring each block to have a size multiple of 4 for the memory to work properly. This allows for any number of compressed points in each block, as long as a padding is used to make its byte size divisible by 4.

A new qfib version was made to obtain this compression format, which is using the 8-bit octahedral quantization for its higher compression capabilities and speed. The compressed blocks are then sent onto the GPU using Shader Storage Bu er Objects (SSBO). With this block format, two points are stored every 62 compressed points, reducing qfib error even more for long bers compared to Section V.2. We wanted to have 64 points in a block for some GPU limitations, so we optimized its compactness by saving the rst and second points as half-oat instead of standard oats, requiring only 2 bytes per coordinate. This reduces accuracy a little, but not signi cantly. Indeed, an average human brain is 150 in the largest dimension. With 16 bits, there are 65, 535 di erent possible values when data are normalized (which is the case here), meaning that the nal precision will be 150/65535 = 0.0023

. The maximum number of bers in the tractogram will be limited to 2 32 = 4, 294, 967, 296 because of the 4 bytes used for storing the ber number. Only one byte is used for the number of compressed points (which is at most 64 due to the block structure) and the number of the current block inside the ber. It means that a ber can contain at most 2 8 × 64 = 16, 384 points. This is not a real limitation as a ber could be divided into two di erent ones with one common point, in case more points are required. The block structure has a xed size of 80 bytes (or less) for the whole tractogram, no matter how many points it contains, which can be as low as one single point. This means that some memory is wasted, but this is in favor of a better parallelism as each GPU thread will be decompressing the exact same quantities of elements in parallel. We are storing 2 points as half-oats every block (so every 64 points), instead of every ber for qfib as it was presented in Section V.2, and we also store some IDs as well as beginning and ending information. Added with the wasted space that occurs for bers that are not multiples of 64 points, this means that this version of our compression will compress bers less than qfib does, with the added bene t of being able to display the data directly.

In OpenGL, there is a limitation for the size of an SSBO, which is limited to 2 of data. As a result, we rst evaluate the required size for the entirety of the tractogram we want to display and create the required number of SSBO to be able to t all data in memory.

V.3.1.2 Strategies for on-the-y decompression

Decompressing bers on the GPU means that we need a way to create the geometry directly on the GPU, instead of sending geometry from the GPU as it is done when rendering uncompressed bers. Figure V.19 presents the two graphics pipelines that are possible since the Turing architecture of Nvidia. The traditional pipeline o ers two possibilities for our on-the-y decompression: the geometry shader and the tessellation, which we will rst experiment with. We will then explore the task/mesh pipeline which was recently introduced, before presenting some possible optimizations and comparative results between pipelines. V.3.1.2.1 Decompression using the geometry shader Using the geometry shader is the most straightforward approach for this kind of application, as its primary goal is precisely to generate geometry on the GPU. One limitation of geometry shader is the number of points and connections that it can output. For a line-strip, which is what we are using here, the limit is 64 points and 63 lines between them. This means that the block format presented in Table V.9 is adapted, as long as we do not want to connect blocks between them. However, bers will be more than 64 points long, especially for datasets with a small stepsize , thus we will use here a maximum of 63 points per block instead of 64, as connecting the last point of a block to the rst one of the next block will require to emit one more vertex. The use of this pipeline -which is schematized in Figure V.20 -does not present any major di culty apart from data organization, as shown previously. We 123 anticipate the worst performances from the geometry shader as it is a well-known fact that it is not really well adapted to the parallel aspect of the GPU.

V.3.1.2.2 Tessellation approach The tessellation unit is capable of generating for each input patch a maximum of 64 output lines (so 65 linked points). This is the reason why we chose to divide each ber into blocks of 64 points, as it allows for inter-block connections for bers continuity. Tessellation works by rst generating a pattern in the tessellation control shader (TCS) in a given con guration (in our case a line composed of 0 to 65 points). These points are then placed at the correct 3D position in the tessellation evaluation shader (TES). This approach is perfect for subdividing geometries as, usually, knowing each point position inside the output points from the TCS, combined with the input patch information, is enough to be able to compute each position. Our compression algorithm being iterative, computing a point requires to know the two previous ones. This cannot be performed inside the TES, which is evaluated in parallel, without points knowing each other. As a result, we chose to decompress points inside the TCS, as we have here all the necessary information. Each patch sent to the TCS will contain a unique integer: the position of the block inside the SSBO.

Passing information from the TCS to the TES is possible but limited to 120 oats maximum. A set of 64 3D points requires 192 oats, which do not t, so we are using half-oats instead, thus requiring the space of only 96 oats. These values are grouped into unsigned integers (uint16) using the packUnorm2x16 function. Then, the TES unpacks the data into oats, for them to be drawn by the fragment shader. To link blocks between them when they are part of the same ber, we use the "Is ending" bit (Table V.9) to know whether we need to create a line to a 65 ℎ point (the rst point of the next block) or not. When using 39 or less points per block, half-oats are no longer required as 40 points need exactly 120 oats to be represented, however packing the data is still a requirement, as the GPU will add padding to the data so that each element is the size of a vec4.

Tessellation is typically faster than the geometry shader as the patterns are better adapted to GPUs.

V.3.1.2.3 Task/Mesh shader decompression The mesh shader pipeline is composed of two programmable shaders, the task (which is optional) and the mesh shader. The task shader function is to launch mesh shaders when required, while the mesh shaders will perform the actual geometry generation. Each of these shaders is using a structure close to compute shaders, allowing a user to control how many threads are running in parallel (with a maximum of 32 threads). The mesh shader being the one responsible for outputting geometry, it will do the decompression, meaning that its task will be iterative, and we cannot bene t from this thread control. At rst, we implemented a version containing only this mesh shader, which had about as good performances as the geometry shader. Without more optimizations, a task shader is not necessary, as it will just slow down the execution, resulting in slower performances as we noticed. However, for the strategy we propose just after, adding this shader will improve a lot the decompression speed. Indeed, using the task shader to launch only some mesh shaders is exactly its purpose, and it can be computed in parallel for 32 blocks, as the choice of rendering a block depends only on the block. Figure V.22 is demonstrating the pipeline using the task shader to early discard blocks that are not visible from the camera. A lot of them are hidden behind other blocks, especially with huge tractograms such as the ones we use. To prevent decompressing blocks that do not participate, we are taking advantage of the previous frame already computed. While rendering a frame, we are recording the depth of each fragment in a depth bu er texture that is then reprojected following the new camera position at the next frame. This depth bu er is moved a little bit away from the camera and is used to estimate the visibility of each block, using its rst and last points. This allows us to render only a portion of all blocks and improves rendering rates by a signi cant factor. When using tessellation, another possible improvement is to adapt the tessellation level (the number of displayed points) depending on the proximity of those points. Indeed, when a ber orientation results in some of its points ending in the same nal fragment, rendering all of them is wasting some resources. By computing this in the TCS, we are able to reduce the nal number of rendered points, however not by a signi cant margin.

V.3.1.3.2 Results

We ran our experiments on a computer equipped with a Core i7-6700k and an RTX 2080 super (8GB of VRAM) with 48GB of RAM. For our results, we are using two tractograms presented in Figure V.23. The rst one contains 500 thousand bers, 358 million points (358, 085, 380) and weigh 4.00 in the tck format, while the second one is composed of 3 million bers, for a total of 4.4 billion points (4,389,303,412) and 49.09 uncompressed (so more than even the RAM!). Figure V.24 presents the size of the two les with di erent numbers of points in the blocks. The highest number of points per block, the smallest the le size, which is expected as we store some more information than the compressed points in each block, requiring more memory. Whatever the number of points per block, qfib compresses more the tractogram as it requires a block per ber, meaning that long bers will require less memory overall. When decreasing the number of points per block, we notice that the le size increases, however, it is lower for 60 points than it is for 61, 62, and 63 points. This is due to the padding that occurs on the GPU, as it can only align memory on 4 bytes (32 bits), meaning that for 61 points, 3 unused bytes are needed for padding. Thus, each block has the same size as for 64 points, with more blocks necessary. In Figures V.23a and V.23b, there is not a real performance improvement when using these block sizes and as a result, we are only using aligned block sizes to look for the best performances we can obtain while maintaining small le sizes. Performances are presented in Figures V.25 and V.26, where the average framerate over 314 di erent points of view around the brain is plotted for di erent block sizes and graphic pipelines. Note that the timings presented in these gures also account for the Screen Space Ambient Occlusion (SSAO) presented in the next section. However, its in uence on timings is negligeable here, as most of the time is spent decompressing the bers. Setting 64 points per block rst appeared to be the ideal quantity in a rst implementation, as it is the highest number of points we can use in our tessellation pipeline. We can observe that the framerate is almost the same for 60 and 64 points per block, making 64 points per block the best value regarding tessellation. The geometry shader is unable to connect blocks together when 64 points per block are used, requiring to use 60 points per block maximum. It is to note that performances increase while the block size decreases, in contrast to tessellation, as the geometry bene ts from more and more parallelization when decreasing the number of emitted points. However, the geometry pipeline remains slower than tessellation. The task/mesh shader pipeline is the most interesting one. Without any optimization, it is the slowest approach, especially at 64 points per block. The di erence between 60 and 64 points can be explained by the limitation of points a single mesh shader can generate: 64. When using 64 points per block, 65 points are required to connect blocks together. This is possible with the mesh shader, in contrast to the geometry shader, but there is a cost as another mesh shader will be launched for this purpose, resulting in slower computing. Where mesh shaders shine however is when using our optimizations. The task shader will analyze 32 blocks in parallel and will create the tasks for the mesh shaders only if the block is visible. This approach being the main purpose of this pipeline, it is really e cient and allows for performances that are about 2.5 times better than tessellation using the same optimizations for the 3 million bers le. This means that the best con guration to optimize both performances and le size is 60 points per block, using the mesh shader pipeline with our optimizations. When this pipeline is not possible (only Nvidia Turing allows it at the moment of writing), the tessellation pipeline is the one o ering the best performances, with this time 64 points per block.

V.3.1.4 Conclusion and future work

To sum up, we introduced a slight modi cation of qfib (Section V.2) that allows for on-the-y decompression of huge tractograms in order to visualize them. This proves the speed of our compression approach and opens a way towards manipulating even bigger datasets as the hardware improves. We are able to obtain interactive visualization framerates with this approach, using only some simple optimizations, and we think that this could be improved even more by nding other ways of testing blocks visibility. Another possible improvement would be to take advantage of all the VRAM available, by decompressing a portion of the tractogram and rendering it with the rest being decompressed on-the-y. Using previous frames, the bers that are most likely to be seen could be selected to be the ones to store decompressed, improving performances even more when compressed tractograms do not take all the available memory.

The compression and visualization techniques presented here and in Section V.2 were also tested on pelvis tractography, with great results. Some surgeons from the Necker hospital tested the visualization tool and were impressed by what they were able to observe that was not possible for them otherwise. They usually extract some portions of the tractography (the nerves inside the pelvis) before visualizing, whereas with our prototype, they were able to see every single line built by tractography algorithms. They could observe the nerves as well as some muscles that are reconstructed in this area, allowing for a better overview of the whole area.

Limitations for this method are the same as the ones of qfib, mainly the fact that we do not compress or take into account other elements along bers, let it be Fractional Anisotropy, a geometry, or other elements. These elements could be easily added using the bers and blocks numbers to retrieve the necessary data from some other SSBOs. Our implementation is also imperfect because we have some ickering when rendering using the optimization. The depth plane is adapted every frame and is causing this issue. Bigger tractograms would bene t from an out-of-core algorithm, with some spatial ordering of blocks in order to minimize data transfers between RAM and VRAM.

V.3.2 Improving tractograms readability

One aspect that was not yet discussed is the readability of the di erent visualizations. Indeed, we presented some geometrical representations, and a visualization algorithm for compressed data, but without further details.

Obviously, the di erent geometric representations are an improvement over using only lines to represent bers, as was already presented in Section IV.1.2, with especially Figure IV.1 where we can compare the impact of using simple generalized cylinders. This volumetric representation is bene cial because it enables better depth perception. However, this is not the only way that depth can be enhanced.

A usual computer graphics method called Ambient Occlusion can be used to achieve good depth perception results [Mittring, 2007]. This consists in darkening parts of a 3D scene, in places where there are some close-by objects, like for instance the corners of a room, or recessed areas on objects. This e ect tries to mimic the self-shadowing e ect that will happen in such regions. It can be computed using several approaches, with a common one being Screen Space Ambient Occlusion (SSAO), which consists in getting the depth of every single pixel before darkening them depending on the depths of the surrounding pixels. Figure V.27 presents the e ect of adding ambient occlusion to the rendering of a 3D object (here a rhinoceros statue). On the left, a render is made in blender without ambient occlusion, compared to the middle where ambient occlusion is activated with an exaggerated factor to make it more visible. The right image represents the computed ambient occlusion at every pixel. The middle image is then a blend of the two outer images. Ambient occlusion helps in this case for sharper edges and contours. This ambient occlusion technique can also be used for tractogram rendering [START_REF] Eichelbaum | Lineao-improved three-dimensional line rendering[END_REF]. In this paper, entitled lineAO, the authors propose an adapted screen space ambient occlusion suited for dense line data, and they experiment it on brain tractograms. Their approach is multi-scale in the occlusion in uence around bers. This allows far away structures to occlude and be occluded while maintaining a good local occlusion as well. The bene ce of adding this e ect is visible in Figure V.28. It improves the volume and depth perception of the streamlines.

Here, we implemented a simpler approach with only local occlusion happening, yielding to satisfying results. However, it would be interesting to compare with lineAO to measure the improvements of their methods as well as the computation cost it might add.

Our simple SSAO requires (like classical ambient occlusion) a di ered rendering pipeline, in which the rst pass will compute a G bu er, consisting in our case of an image rendered without ambient occlusion and a depth image. Then, the second pass will use the depth image to infer an ambient occlusion value for every pixel, and it will combine the value of the rendered image with the ambient occlusion. This combination makes it possible to vary the quantity of ambient occlusion a user might want. This is illustrated in Figure V.29 where a whole-brain tractogram is displayed with three di erent color schemes and three di erent ambient occlusion levels: no ambient occlusion, 50% ambient occlusion, and 100% ambient occlusion. At 100%, it simply means that the image is multiplied by the ambient occlusion value.

For brain bers, the ambient occlusion is computed such that pixels representing bers that are further away from the camera than the surrounding bers will appear darker, allowing for better depth perception. The e ect is even more convincing when a single color is used to render bers as lines, like in images d, e, and f of Figure V.29. Indeed, when using the same colors for all bers without a light shading, it results in a solid color without ambient occlusion. However, as soon as it is added, individual bers become visible.

We also noticed that when rendering massive tractograms, or when zooming out, this darkening e ect reveals the cortical surface. This is shown in Figure V.30. This is an expected result, because when zooming out or displaying lots of bers, the depth di erence is actually between the extremities of the bers, and other parts of the bers, with extremities being on the cortical surface by de nition and construction. This simple approach yields interesting results and it would be interesting to add other elements such as light shading [START_REF] Mallo | Illuminated lines revisited[END_REF].

Transparency is another visualization feature that would help a lot in speci c applications, for instance visualizing a bundle or a fascicle inside another one, or even inside a whole-brain tractogram. In this case, being able to display the surrounding bers in transparency would help for the localization of the tract of interest. However, this presents some major challenges, as the rendering of bers as lines is usually fast enough because most of them are hidden behind the front bers, meaning that they are quickly discarded and do not require as much processing power. Transparency would require depth ordering of the bers as well as the rendering of all bers, which is costly. A possibility to explore would be to render separately both the whole-brain and the bundle with their depth, before combining the two images. The depth layers would then help to determine in which case to apply transparency or not.

For improving readability further, some inspiration could also be taken from vector rendering. For instance, in a recent article [START_REF] Amiraghdam | Localis: Locally-adaptive line simpli cation for gpu-based geographic vector data visualization[END_REF], lines representing roads on a map are simpli ed when rendered far away from the camera, with a guarantee for no intersections between roads. These kinds of algorithms could be quite interesting to adapt to bers while integrating their speci cities.

V.4 Conclusion V.4.1 Summary

In this chapter, we presented a segmentation method that mathematically de nes the inherent imprecision of relation de nitions of bundles using the theory of fuzzy sets. An anatomical coherence score is introduced, quantifying the probability of belonging to a bundle with respect to the de nition used. Consequently, it requires determining a threshold of this score to obtain the bundle. This is performed with the help of our multi-resolution representation, that even considers the score in the similarity between bers, leading to a better representation of the bundle. We also presented a compression algorithm for bers. It considers their segments as a succession of unit vectors, and only records them as quantized angles with respect to the previous unit vector. This results in high compression ratios, as well as a fast algorithm. Each ber is independent regarding the compression, meaning that it can be parallelized on the GPU. Some adjustments were required for a better adaptation to their speci cities, resulting in an on-the-y decompression for visualization of the compressed tractograms. Therefore, tractograms containing millions of bers can be visualized in interactive time, which was previously impossible. Visualizing so many lines requires some visualization techniques in order mainly to be able to perceive depth. We presented some images with a simple screen space ambient occlusion, which already drastically improves the 3D perception of the tractogram.

V.4.2 Discussion

Other bundles than the two presented here could be segmented using the same approach. For this purpose, the de nition of the bundle would need to be clari ed, so that specialists agree on their validity. Then, the reproducibility would need to be evaluated in the same way. The major limitation of the compression is its requirement of equally distant points. In order for this limitation to disappear and the compression to remain e cient, one should nd a way of compressing the varying distances. One possibility would be to quantize the variation of distance compared to a reference that could exist per-ber. The visualization optimization could also be further improved by discarding more occluded bers and solving the ickering that can happen in some cases with our implementation. We only scratched the surface in terms of tractogram visualization, and lots of other techniques would need to be explored, as well as interaction methods that play an important role in the user understanding of the data.

VI Conclusion

VI.1 Summary of the contributions

In this thesis, we studied brain tractography and provided new algorithms to help make them easier to use. We also improved the complexity of a well-known Moving Least Square (MLS) projection algorithm, opening the way for the study of the use of di erent kernels.

We rst presented a new multi-resolution representation suited to brain tractograms. This algorithm is inspired by the computer graphics simpli cation method QSlim [START_REF] Garland | Surface Simpli cation Using Quadric Error Metrics[END_REF], which progressively groups vertices minimizing a quadric error metric. We perform a similar grouping for bers, merging rst the pairs that are the closest ones according to a chosen metric. To make this technique usable, we also proposed to use a Delaunay tetrahedralization to create an adjacency relation between bers. It relies on the extremities of these bers, which are an important anatomical element. This allows us to create a whole multi-resolution representation for a tractogram in a reasonable time. This representation can be navigated in real-time, bringing an accessible resolution choice to the user depending on their needs and computer power.

The geometric representation of bers and meta-bers is done using generalized cylinders with elliptical pro les, which provide good performances, at the expense of an unideal geometric representation because of the crossing of elements and spatial extents that are not always accurate. Guided by this observation, we then studied possible geometric representations and focused on MLS projections using algebraic point set surfaces (APSS) [START_REF] Guennebaud | Algebraic Point Set Surfaces[END_REF]. It appeared that for our needs, the usual nearest neighbors projection was not suited for its inability to ll in holes. However, the global approach was really expensive to compute.

As a result, we proposed to reduce APSS complexity by considering groups of close points (nodes) instead of individual points when far away from those points. This is performed with a blending between those nodes and smaller nodes of individual points for continuity. We created a GPU version of this algorithm that is able to quickly project many points. Kernels are used to leverage the importance of points depending on their distance to the evaluation position, and their values can be changed for ltering or changing the nal appearance of the obtained surface. The hole-lling capabilities of the method required us to create a new sampling technique based on an adaptive octree that is created using the MLS operator. This results in a fast initialization close to the nal surface.

This operator was experimented on tractogram data and seems promising. However, some other problems appeared, mainly due to the shape of bers that might require some speci c kernels. Consequently, the proposed volumetric metric to compare bers was not tested with the multi-resolution representation. Provided it results in correct mergings, the introduced MLS operator solves the spatial extent issues we previously encountered, however, the crossing geometries still happen as long as we consider every ber and meta-ber independently. One solution to avoid those problems is to consider every ber and meta-ber as implicit surfaces and to combine them. This allows the construction of a single surface corresponding to the whole dataset, with implicit combinations of the geometries, parametrized by these operators. We showed some results obtained using a basic ray-marching approach.

Segmentation is an important aspect of tractograms studies as it is often required to extract some known fascicles. The challenge is to transform the vagueness of the de nitions of these fascicles into something usable by a computer. We did that using the fuzzy set theory, enabling the computation of an anatomical coherence score for every ber. This score corresponds to the con dence we can have that this ber belongs to the fascicle. We combined this score with the weighted current metric [START_REF] Gori | Parsimonious approximation of streamline trajectories in white matter ber bundles[END_REF] and used it to compute the multi-resolution representation. The result is an interactive multi-resolution visualization that can be used to nd the desired threshold to obtain the fascicle. We also showed that using the same threshold for a fascicle in di erent brains results in a better reproducibility than some state-of-the-art approaches.

Tractograms can contain lots of bers, leading to massive le sizes. Existing compression methods [START_REF] Presseau | A new compression format for ber tracking datasets[END_REF] do not scale well and are quite slow. Taking advantage of the ber construction, we introduced a new accurate and fast compression algorithm. It quantizes consecutive directions of individual bers using unit vector quantization combined with a recent mapping from a spherical cap to a unit sphere [START_REF] Rousseau | Fast Lossy Compression of 3D Unit Vector Sets[END_REF]. It reaches 80 to 90% compression with computation times in seconds on CPUs. Going even further, we implemented some GPU version of the decompression with some adaptations. This allows us to reach interactive visualizations of massive tractograms without the need for decompression. As a result, tractograms weighing more than the available VRAM of the GPU can be visualized in their entirety. This is performed using the graphics pipeline in some innovative ways. We nished by presenting a few computer graphics visualization techniques that bene t a lot to tractograms visualization.

To validate the work of this thesis, we asked surgeons for their opinion on the multi-resolution representation, the di erent de nitions of the Uncinate Fasciculus and the IFOF as well as the interactive visualization of compressed tractograms. Those feedbacks were really positive and helpful for guiding us into a path that would hopefully lead to useful improvements in their work.

VI.2 Discussion and future work

Our multi-resolution representation is, in its current state, limited to metrics using only the points of the bers. An extension to other information present along them, like fractional anisotropy or functional information, would be interesting. Some of the work presented in this thesis is still a work in progress. It is especially the case for the use of FastAPSS to improve the geometric representation of bers. Some issues appeared while testing this representation that are still present at the moment of writing this manuscript. A way of distinguishing cases when bers are intricated into one another, to make sure that we do not remove too many points for the new meta-ber, needs to be found. This could be reached by using a clever threshold or some neighboring check with points from both bers to merge. Then following this, testing the proposed comparison tool which uses the volume di erence induced by the merge into a meta-ber is something to study and to compare with other metrics. This geometric representation using fastAPSS could also be combined with the implicit surface representation. The volume aware similarity could rst compute a merging order of the individual meta-bers built with fastAPSS. For a given resolution, the rendering could then be performed using di erent operators for grouped meta-bers and separated ones, to obtain for instance natural merges for grouped meta-bers (like what was presented at the bundle level in Section IV.4) and no merges between di erent groups. Some approaches to improve ray marching performances could also be explored. This could be done using information from the bundles, and especially taking into account which parts are visible from the camera instead of computing everything. It could be something a little similar to the optimizations proposed for the interactive decompression (Section V.3.1). Another possibility would be using a dual-contouring to create a mesh, however, some optimizations would be needed to obtain a reasonably sized mesh.

Segmentation with our algorithm yields good results on the chosen fascicles. Other fascicles would be interesting to obtain using the same approach. However, relation de nitions for fascicles based on brain regions are not clearly existing, with variations between di erent surgeons or researchers, meaning that this procedure is time-consuming.

Concerning the new compression algorithm, an improvement would be to handle bers that have points at di erent intervals. This could be performed using a quantization of the variation of length between segments on top of the quantization on the directions. This would result in less compression and more importantly a reduced precision. However, considering only small chunks of bers-like we do for on-the-y decompression on the GPU -would help improve and control this precision. This is especially necessary as the available tools in the community for resampling bers do not equally spaced points as they are supposed to do. Indeed, they build a new ber by following the polyline, dividing its total length into the required number of segments. This is usually ne for increasing the number of points or for almost straight bers, however, when reducing the number of points, especially on a highly-curved ber, it tends to create some points that are closer than others. In the case of the use of our compression algorithm, it then results in a high error as the points will indeed be moved a lot compared to the original ones. A recent work in tractogram compression compares itself to us [START_REF] Haehn | Trako: E cient transmission of tractography data for visualization[END_REF]. They obtain high errors with our algorithm because of these resampling approaches. This is visible mainly in their use of our dataset (with equally distant points by construction), where we remain better than them.

Visualization of tractograms was only brie y addressed here and should be further investigated. As mentioned, transparency would improve a lot the quality of the perceived information, but other elements would also be really bene cial. Using some colors, textures, or even geometries to represent information along the bers have already been tried a few times, but are not yet adopted by surgeons. This is true for every work that remains to be done here, and even for other domains, but frequent exchanges, discussions, and demonstrations with surgeons are a necessity for reaching good tools that are bene cial to patients and doctors.

It is to be noted that tractography can also be obtained on other parts of the human body, such as the heart or the pelvis. It would be interesting to try the presented algorithms on these speci c tractographies. We already experimented the interactive visualization of compressed data on the pelvis tractography, and it seemed to function as intended. However, some further testings and investigations on the properties of these bers should be performed. Les données de tractographie (fibres) obtenues à partir d'IRM de diffusion sont difficiles d'utilisation. Dans cette thèse, nous proposons des méthodes et algorithmes pour la simplification, la visualisation et la manipulation de ces données. Nous introduisons une représentation multirésolution des tractogrammes, plus rapides et avec une meilleure précision géométrique que les approches de simplification existantes. Nous explorons aussi diverses représentations géométriques et nous nous concentrons sur les approches de projections aux moindres carrés (MLS) par l'intermédiaire des surfaces algébriques d'ensemble de points (APSS), pour lesquelles nous réduisons la complexité, permettant l'utilisation de noyaux globaux pour l'analyse et la modélisation. Une technique de segmentation utilisant la représentation multi-résolution et permettant une meilleure reproductibilité que d'autres approches est ensuite présentée. Les tractogrammes pouvant être volumineux, nous introduisons un algorithme de compression exploitant la manière d'obtenir les données à partir des IRM de diffusion. La vitesse de cet algorithme permet même son utilisation pour la visualisation de données compressées, la décompression se faisant à la volée sur le GPU. Ces travaux de recherche et les résultats obtenus se situent à l'intersection de l'informatique graphique et de l'analyse de données médicales, ouvrant de nombreuses perspectives.

Title: Geometrical modeling, simplification and visualization of brain white matter tractograms Keywords: shape modeling, medical imaging, geometry, simplification, visualization, segmentation Abstract: Tractography data (fibers) obtained from diffusion MRI present several challenges. In this thesis, we propose some useful methods and algorithms for simplification, visualization, and manipulation of these data. We introduce a new multi-resolution representation for tractograms, faster, and with higher geometric accuracy than existing simplification approaches. We also investigate various geometric representations and focus on moving least square (MLS) projection with algebraic point set surfaces (APSS), on which we reduce the complexity, allowing for the use of global kernels for analysis and modeling. A segmenta-tion technique using the multi-resolution representation is presented, achieving better reproducibility than other approaches. Tractograms being massive, we also introduce a compression algorithm taking advantage of data obtention from diffusion MRI. The algorithm speed even allows for the direct use of compressed data for visualization, as it can be decompressed on-the-fly on the GPU. This research and the obtained results lie at the intersection between Computer Graphics and Medical Data Analysis, paving the way for numerous perspectives.

Institut Polytechnique de Paris 91120 Palaiseau, France

Figure

  Figure II.1 presents an MRI image of a human brain, the fractional anisotropy measured in it, as well as a tractogram built from these data.

  (a) A random seed is placed in the white matter (white region).(b) Following the di usion direction, the trajectory is built.(c) Fibers stop at the cortical surface (black region), and the angle between consecutive segments is limited.(d) The process is iterated in the opposite direction.

Figure III. 1 -

 1 Figure III.1 -Construction of a ber.

( a )

 a Imaging data. (b) Reconstruction 1. (c) Reconstruction 2. (d) Reconstruction 3.

Figure III. 2 -

 2 Figure III.2 -Possible reconstructions of imaging data, some of them can be incorrect.

Figure

  Figure III.4 -Multi-resolution pipeline: a) original bundle tractogram, here the thalamocortical one, b) connections created by the Delaunay tetrahedralization based on the extremities of the bers, c) the similarity is computed for each ber, taking into account only the neighbors found in step b), d) couples of closest bers are progressively merged into generalized cylinders, e) the nal multi-resolution representation makes it possible to navigate through the di erent levels of detail in real-time. The percentage refers to the fraction of employed generalized cylinders compared to the original number of bers.

  First, the body and the extremity parts are found. Finally, 3 pro le is created.

Figure III. 5 -

 5 Figure III.5 -Scheme of the merging of two cylinders 1 and 2 into 3 .

Figure

  Figure III.6a) Thalamocortical bundle with 19,782 bers, b) reduced to 275 prototypes with QuickBundles (Threshold=7mm), c), d) and e) reduced to 6925, 1422, and 275 cylinders with our method.

  Figure III.6 (c,d, and e). Our result, in Figure III.6e, has as many cylinders as prototypes in Figure III.6b. The proposed method, after a single and fast pre-computation (see Table

Figure

  Figure III.7 shows results on three di erent bundles: the Uncinate Fasciculus (a, b), the IFOF (c), and a whole-brain (d). The measures used for the experiments are MC (Equation III.1) for (a), WC (Equation III.3) for (b and d), with = 6 , = 6 and = 8 , and MDF (Equation III.2) for (c).

Figure III. 7 -

 7 Figure III.7 -Bundles at di erent resolutions: a) Uncinate Fasciculus, computed with the mean of closest distances; b) Uncinate Fasciculus, computed with the dissimilarity of the weighted currents; c) IFOF bundle, computed with the dissimilarity of the weighted currents; d) whole-brain tractogram, with 1 million bers, computed with the dissimilarity of the weighted currents. The graphs on the right represent the valence of the bers depending on the resolution. The maximum valence is in red, and the average valence is in blue.

Figure III. 8 -

 8 Figure III.8 -Computation time for the multi-resolution representation on the Uncinate Fasciculus (UF).

Figure III. 9 -

 9 Figure III.9 -Computation time for the multi-resolution representation on the Inferior Fronto-Occipital Fasciculus (IFOF) bundle.

Figure III. 10 -

 10 Figure III.10 -Computation time for the multi-resolution representation on the Thalamocortical bundle.

Figure III. 11 -

 11 Figure III.11 -Computation time for the multi-resolution representation on a 1 million bers tractogram of the whole brain.

  (a) 5% resolution. (b) 4% resolution.(c) 3% resolution.(d) 2% resolution.(e) 1% resolution.(f) 1 meta-ber.

Figure III. 12 -

 12 Figure III.12 -Example of oversimpli cation on the Uncinate Fasciculus.

( a )

 a With outliers. (b) Without ouliers.

Figure III. 13 -

 13 Figure III.13 -Thalamocortical bundle with and without outliers. The maximum angle was set to 70 • . The red arrows points to some places where outliers where removed.

Figure III. 14 -

 14 Figure III.14 -Comparison of the metrics used on the IFOF at 5% resolution.

Figure III. 15 -

 15 Figure III.15 -Comparison of the metrics used on the Uncinate Fasciculus at 3% resolution.

Figure IV. 1 -

 1 Figure IV.1 -Visualization from the same point of view of lines (top) and cylinders (bottom)

  Figure III.5. This geometry also creates quickly varying pro les, non-smooth geometries, and auto-intersections at several places, as can be observed in Figure IV.2.

  (a) Auto-intersections and non-smooth geometry can appear with the use of generalized cylinders for bundles.(b) The pro le can vary quickly in some cases, resulting in non-smooth geometries.(c) Auto-intersections are preventing a good visualization of this bundle.

Figure

  Figure IV.2 -Geometric limitations of the use of generalized cylinders.

  Figure IV.3 -Modeling Algebraic Point Set Surfaces using a xed number of input nearest points results in unacceptable approximations far from densely sampled regions. Our smooth approximation copes naturally with complex inputs featuring large missing parts and competes with global approaches while allowing for simple point set ltering and smooth point projection.

Figure

  Figure IV.5 -Example of a tree traversal. Nodes are traversed following the colored arrows. Colors on arrows and nodes correspond to the colors of the pseudo-code in Algorithm 4.

Figure

  Figure IV.6 -Projecting points, with and without ordering points following the Morton order. Even though ordering points using the Morton code is not free, doing so enables us to project points consistently faster.

FigureFigure

  Figure IV.10 -Our method is accurate with respect to global APSS, as shown here by the distance between them, relative to the bounding box diagonal.

  Figure IV.14 -While the time to project 100000 points is in O (|P |) for GlobalAPSS, it is in O (log(|P |)) for FastAPSS.

Figure

  Figure compares some results on the same kind of data, this time with 20 kNN and global APSS only. It shows the bene t of the global and fast APSS approaches, which leads to similar results, compared to the standard nearestneighbors approach. Our results have a massive speed bene t compared to global APSS, for equivalent results (see TableIV.1).

  Figure IV.18 -Even for very small portions of the unit sphere (see Part Sphere 3, 21 sample points), our octree generation algorithm still allows for a faithful point cloud upsampling leading to the expected dual-contoured surface.

  Figure IV.19 -Our approach allows for e cient large-scale continuous ltering of 3D point sets. The methods we compare to take as input precision-related parameters and are not suited for point set smoothing and ltering.

Figure IV. 20 -

 20 Figure IV.20 -Example of a ber, its sampling on a cylinder, and the mesh obtained using FastAPSS with the kernel from Equation IV.9.

Figure

  Figure IV.22 shows a merge resulting in a convincing meta-ber, in the sense that the obtained mesh corresponds to an expected geometry.

( a )

 a Two original bers. (b) Mesh corresponding to one original ber.

Figure IV. 23 -

 23 Figure IV.23 -Example of an incorrect reconstruction when using FastAPSS.

Figure IV. 24 -

 24 Figure IV.24 -Illustration of the reason of a merge failure, resulting in 3 = 1 ∪ 2 . The gray points are discarded because they are considered inside of the other metaber.

Figure IV. 25 -

 25 Figure IV.25 -Visualization of the Uncinate Fasciculus using implicit surfaces. The top gure is using a smaller parameter than the bottom one, merging less bers.

Figure

  Figure IV.26 -Visualization of the IFOF bundle using implicit surfaces. The top gure is using a smaller parameter than the bottom one, merging less bers.

Figure IV. 27 -

 27 Figure IV.27 -Visualization of the Thalamocortical bundle using implicit surfaces.

Figure V. 1 -

 1 Figure V.1 -Illustration of a fuzzy cone, corresponding to the relation "anterior to the hippocampus".

Figure V. 3 -

 3 Figure V.3 -Scheme with notations used for connectivity computation (EP).

Figure V. 4 -

 4 Figure V.4 -Half brain with bers colored according to the Anatomical Coherence Score values of the UF.

Figure V. 5 -

 5 Figure V.5 -Scheme with notations used for the extension of the weighted currents ( ) in Equation V.7.

Figure

Figure V. 6 -

 6 Figure V.6 -Anatomical regions used for segmentation. In beige: gyrus, in purple: hippocampus, in black: lingual, in yellow: fusiform, in red: cuneus, in blue: pericalcarine, and in green: lateraloccipital.

Figure V. 7 -

 7 Figure V.7 -Segmentation of IFOF and UF bundles of three subjects with three di erent ACS thresholds (0.5, 0.65, 0.7 for the IFOF and 0.7, 0.85, 0.9 for the UF respectively). Results are shown on the MNI152 T1w image.

Figure V. 8 -

 8 Figure V.8 -IFOF and UF bundles visualized at three di erent resolution levels and segmented with two di erent ACS thresholds (0.4, 0.57 for the IFOF and 0.6, 0.81 for the UF respectively).

Figure

  Figure V.10 -A unit vector (in green) is approximated by the closest quantization point (in orange).

  Figure V.11 -Octahedral quantization[START_REF] Meyer | On Floating-Point Normal Vectors[END_REF] projects a unit vector to an octahedron, then, to a unit square to encode the discretization of the resulting 2D coordinates.

  Figure V.12 -Notations used in Section V.2.2.3. The mapping introduced in[START_REF] Rousseau | Fast Lossy Compression of 3D Unit Vector Sets[END_REF] transforms a point , de ned on a spherical cap centered in P 0 and with a maximal angle , represented as the yellow area on the surface of the unit sphere, to another point x , de ned on the surface of the whole unit sphere, colored in light blue.

Figure V. 13 -

 13 Figure V.13 -Quantization with (in green) and without (in orange) propagation error reduction. The error computed per point is shown in red.

Figure V. 14 -

 14 Figure V.14 -Example on an S curve ber (black). The compressed and decompressed ber is shown without (in orange) and with error propagation reduction (in green).

Figure V. 16 -

 16 Figure V.16 -Original bers (black), compressed and decompressed bers (green) using an 8 bits octahedral quantization.

Figure V. 17 -

 17 Figure V.17 -Error in mm depending on bers length for an 8 bits octahedral quantization with = 0.1 .

Figure V. 18 -

 18 Figure V.18 -The gure on the left of each panel represents the worst bers extracted from the gure on the right hand side of the panel, showing the local compression error (see color code on the top of the gure) for the 8-bits (a) and 16-bits (b) octahedral quantization. The dataset was computed using the deterministic algorithm and = 1 .

Figure V. 19 -

 19 Figure V.19 -Graphics pipeline.

Figure V. 20 -

 20 Figure V.20 -Decompression scheme using the geometry shader.

  Figure V.21 illustrates this pipeline.

Figure V. 21 -

 21 Figure V.21 -Decompression scheme using the tessellation shader.

Figure V. 22 -

 22 Figure V.22 -Decompression scheme using the mesh shader.

Figure

  Figure V.23 -Visualization of compressed tractograms.

Figure V. 24 -

 24 Figure V.24 -Compressed les size depending on block size.

Figure V. 25 -

 25 Figure V.25 -Performances when displaying 500k bers.

Figure V. 26 -

 26 Figure V.26 -Performances when displaying 3M bers.

Figure V. 27 -

 27 Figure V.27 -From left to right, Rhinoceros rendered without Screen Space Ambient Occlusion (SSAO), rendered with SSAO, and SSAO layer.

( a )

 a Phong illuminated lines[START_REF] Mallo | Illuminated lines revisited[END_REF].(b) LineAO[START_REF] Eichelbaum | Lineao-improved three-dimensional line rendering[END_REF] with Phong illumination.

Figure

  Figure V.28 -Streamlines around the main vortices of a delta wing data set obtained from a uid dynamics simulation. Figure extracted from the article of lineAO.

  No ambient occlusion, grey colors.(e) 50% ambient occlusion, grey colors.(f) Full ambient occlusion, grey colors.(g) No ambient occlusion, colors by orientation.(h) 50% ambient occlusion, colors by orientation.(i) Full ambient occlusion, colors by orientation.

Figure

  Figure V.29 -E ect of ambient occlusion when visualizing 500k bers.

  Figure V.30 -The cortical surface becomes visible thanks to ambient occlusion when zooming out.

Titre:

  Modélisation géométrique, simplification et visualisation des fibres de la matière blanche du cerveau Mots clés: modélisation 3D, imagerie médicale, géométrie, simplification, visualisation, segmentation Résumé:

•

  Segmentation of White Matter Tractograms Using Fuzzy Spatial Relations, Alessandro Delmonte, IsabelleBloch, Dominique Hasboun, Corentin Mercier, Johan Pallud, Pietro Gori (OHBM, 2018)

	• Progressive and E cient Multi-Resolution Representations for Brain Tractograms, Corentin Mercier, Pietro Gori, Damien Rohmer, Marie-Paule Cani, Tamy Boubekeur, Jean-Marc Thiery, Isabelle Bloch (EG VCBM, Euro-graphics Workshop, 2018) https://github.com/CorentinMercier/neural-meta-tracts • White Matter Multi-resolution Segmentation using Fuzzy Theory, Corentin Mercier * , Alessandro Delmonte * , Johan Pallud, Isabelle Bloch, Pietro Gori (ISBI, 2019) https://github.com/CorentinMercier/FBTS • QFib: Fast and Accurate Compression of White Matter Tractograms, • QFib: Fast and E cient Brain Tractogram Compression, Corentin Mer-cier * , Sylvain Rousseau * , Pietro Gori, Isabelle Bloch, Tamy Boubekeur (Neu-roInformatics, 2020) https://github.com/syrousseau/qfib • FastAPSS: E cient Algebraic Point Set Surfaces using Non-compact Kernels, Corentin Mercier, Thibault Lescoat, Pierre Roussillon, Tamy Boube-Corentin Mercier Boubekeur (OHBM, 2019) keur, Jean-Marc Thiery (In submission)

* , Sylvain Rousseau * , Pietro Gori, Isabelle Bloch, Tamy

Table III

 III

			Number Maximum Average	Tet		Total time (s)
	Bundle	Metric
			of bers	valence	valence	time (s) single-core Multi-core	CUDA
		MC		85	25.4 ± 3.2	0.01	5.62	1.90	1.89
	Uncinate	WC		93	23.9 ± 3.2	0.01	9.37	2.79	2.17
	Fasciculus	MDF	940	85	25.1 ± 3.3	0.01	0.97	1.08	1.02
		QB (MDF)		-	-	-	0.02	-	-
		MC		494	29.8 ± 5.8	0.03	31.7	6.87	6.80
		WC		513	30.3 ± 6.9	0.03	64.6	10.7	9.71
	IFOF	MDF	1, 983	428	30.2 ± 3.8	0.03	5.55	2.63	2.60
		QB (MDF)		-	-	-	0.07	-	-
		MC		383	32.0 ± 5.0	0.28	163.6	43.0	38.3
	Thalamo-	WC		733	29.5 ± 2.6	0.28	333.2	70.4	61.5
	cortical	MDF	19, 782	226	30.2 ± 2.9	0.28	28.2	19.4	19.2
	Bundle
		QB (MDF)		-	-	-	4.03	-	-
		MC		973	31.0 ± 3.9	15.42	6, 426.3	2, 322.9	1, 944.9
	Whole-	WC		1, 018	28.1 ± 2.0	15.42	16, 827.4	3, 685.1	3, 222.8
	Brain	MDF	1, 000, 000	404	29.9 ± 2.6	15.42	1, 553.6	1, 075.7	1, 046.8
	Tractogram
		QB (MDF)		-	-	-	4, 719.3	-	-
	Figures III.8 to III.11 present the same computation times using some histograms
	for an easier comparison between QuickBundles and our multi-resolution repre-
	sentation.

.1 -Maximum and average valence with computation times for some bundles. The computation time for QuickBundles (QB) is presented for comparison. The threshold used for obtaining a single simpli cation was 7 .

  It should t the weighted point set (w, P) both in terms of positions (i. e., (p i )

	0) and normals (i. e.,	-→ (p i

  ± 4.35 17.20 ± 4.37 36.19 ± 3.92 146.31 ± 32.32 1850291.25 ± 1284.64 16.24 ± 0.85 45532.1 ± 4580.6 Table IV.1 -Average time to project one point (mean + standard deviation).

					CPU time ( )		GPU time ( )
	Model	Points
			20 kNN	40 kNN	100 kNN	FastAPSS	GlobalAPSS	FastAPSS	Global APSS
	Guitar		1.44 ± 0.07	2.76 ± 0.11	6.67 ± 0.12	32.24 ± 1.69	614.69 ± 24.45	2.91 ± 0.17	8.22 ± 0.19
	Face		2.78 ± 0.28	3.65 ± 0.46	5.79 ± 0.42	5.82 ± 0.31	851.87 ± 9.91	0.48 ± 0.02	25.25 ± 0.24
	Dinosaur		1.45 ± 0.10	2.77 ± 0.12	6.82 ± 0.16	28.62 ± 2.01	3049.67 ± 16.41	2.76 ± 0.16	50.04 ± 0.23
	Lord quasimodo		1.43 ± 1.05	2.47 ± 1.09	5.55 ± 1.21	18.13 ± 0.84	2259.44 ± 12.89	1.84 ± 0.07	41.23 ± 0.40
	Daratech		3.51 ± 0.08	6.60 ± 0.10 16.28 ± 0.12	71.19 ± 8.56	8719.35 ± 60.22	6.94 ± 0.69	113.17 ± 1.58
	Dancing children		1.12 ± 0.53	1.96 ± 0.58	4.48 ± 0.64	13.32 ± 0.61	2334.96 ± 28.60	1.40 ± 0.05	47.35 ± 0.41
	Anchor		0.66 ± 0.06	1.21 ± 0.09	2.77 ± 0.14	7.47 ± 0.50	1332.92 ± 52.77	0.85 ± 0.05	51.63 ± 0.42
	Gargoyle		1.11 ± 0.36	1.98 ± 0.42	4.56 ± 0.52	12.68 ± 0.48	3000.82 ± 1.72	1.40 ± 0.05	63.58 ± 0.48
	Suitcase		0.91 ± 0.32	1.59 ± 0.39	3.43 ± 0.50	9.93 ± 0.63	1958.61 ± 13.12	1.05 ± 0.05	64.22 ± 0.38
	Igea	134346	1.07 ± 0.28	1.93 ± 0.31	4.47 ± 0.40	12.98 ± 0.72	3978.08 ± 19.56	1.39 ± 0.06	93.70 ± 0.57
	Armadillo	172974	2.26 ± 0.12	5.30 ± 2.52 10.97 ± 0.16	39.43 ± 2.83	15971.80 ± 41.18	3.83 ± 0.18	220.64 ± 1.23
	African Statue	220317	2.05 ± 0.39	3.65 ± 0.39	8.28 ± 0.40	32.07 ± 2.34	11956.50 ± 21.60	3.03 ± 0.19	200.38 ± 1.03
	Napoleon	495000	1.29 ± 0.20	2.28 ± 0.22	5.00 ± 0.24	15.20 ± 1.34	12165.80 ± 23.29	1.64 ± 0.14	344.54 ± 2.02
	Bearded Man	499500	1.04 ± 0.26	1.77 ± 0.33	3.60 ± 0.36	10.07 ± 1.01	9223.72 ± 18.41	1.28 ± 0.14	320.86 ± 3.53
	Eisbar	781865	5.75 ± 1.09	8.53 ± 0.91 14.62 ± 0.46	24.56 ± 1.89	44060.40 ± 154.45	2.16 ± 0.16	696.79 ± 7.69
	Owl	1031960	5.40 ± 4.55	7.49 ± 4.56 13.15 ± 4.67	30.90 ± 2.42	67452.00 ± 168.53	2.46 ± 0.18	1076.94 ± 23.67
	Dragon	1180060	2.04 ± 1.90	3.00 ± 2.02	5.60 ± 2.40	37.32 ± 2.80	26116.10 ± 85.07	4.26 ± 0.35	776.88 ± 5.51
	Rhinoceros	1410356	6.28 ± 6.60	9.50 ± 6.45 18.05 ± 6.75	32.87 ± 2.79	86970.40 ± 137.95	3.39 ± 0.22	1418.64 ± 10.01
	Xyzrgb Dragon	3609455 5.06 ± 10.66 5.69 ± 10.96 7.21 ± 11.82	8.65 ± 1.90	17003.30 ± 642.83	1.77 ± 0.57	2078.57 ± 25.93
	Lucy	14027872 10.94
						,

Table IV

 IV

	Model	Points D T	(ms) U		(ms) # nodes
		6	980.3	189.1	18690
	Part Sphere3	21
		8	6295.8	1994.7
		6	321.1	200.1	23417
	Part Sphere2	2947
		8	4168.4	3391.7
		6	225.4	55.6	3245
	Part Sphere	30650
		8	2034.7	406.4	48967
		6	572.5	351.3	40387
	Sphere	100000
		8	8331.1	6710.8
		6	297.8	87.1	8180
	Face	40881
		8	3113.4	1325.2
		6	707.0	178.9	22281
	Bearded Man 499500
		8	7215.5	3700.8

.2 -Octree sampling computation time for the models in Figure IV.18, as well as two other models from Figure IV.12. D: maximum depth of the octree.

  The Uncinate Fasciculus (UF) and the Inferior Fronto-Occipital Fasciculus (IFOF) are shown in Figure IV.25 and IV.26 respectively, with two values of the parameter (see Equation IV.15). These fascicles contain 940 and 1983 bers, respectively. This parameter allows for a di erent radius around individual bers, changing the number of merges obtained through the addition of the implicit functions. The top image in both gures is obtained with a smaller parameter value than the bottom one, thus having a smaller encapsulating radius, making individual bers more visible.Figure IV.27 presents the Thalamocortical bundle, composed of 19782 individual bers.

  Table V.1 -Size of the les (tck format).

	Stepsize	0.1 mm	0.2 mm	0.5 mm
	Nb. of bers 500k 3M 500k 3M 500k 3M
		Size of input le (GB)
	Deterministic 3.80 22.8 1.92 11.5 0.84 5.01
	Probabilistic	4.36 26.2 2.35 14.1 1.26 7.53

Table V .

 V 2 -Maximum and average compression errors of our method depending on the dataset, precision in bits and quantization method.

	Stepsize	0.1 mm	0.2 mm	0.5 mm
	Nb. of bers	500k 3M 500k 3M 500k 3M
	Maximum error (×10 -2	)
	Quantization Precision	Deterministic
	8 bits	4.94 5.30 10.3 10.6 33.2 34.1
	Fibonacci
	16 bits	0.50 0.57 0.28 0.28 0.39 0.38
	8 bits	7.53 8.03 16.5 16.5 46.7 51.0
	Octahedral
	16 bits	0.50 0.56 0.27 0.29 0.50 0.52
	Quantization Precision	Probabilistic
	8 bits	3.04 2.95 5.86 5.91 19.7 20.6
	Fibonacci
	16 bits	0.13 0.14 0.14 0.16 0.69 0.72
	8 bits	4.79 4.90 8.55 9.38 29.8 31.7
	Octahedral
	16 bits	0.13 0.14 0.17 0.19 0.87 0.93
	Average error (×10 -2	)
	Quantization Precision	Deterministic
	8 bits	1.20 1.20 2.62 2.62 7.90 7.90
	Fibonacci
	16 bits	0.03 0.03 0.03 0.03 0.11 0.11
	8 bits	1.74 1.74 3.80 3.80 11.5 11.5
	Octahedral
	16 bits	0.03 0.03 0.03 0.03 0.12 0.12
	Quantization Precision	Probabilistic
	8 bits	1.33 1.34 2.14 2.15 5.77 5.78
	Fibonacci
	16 bits	0.02 0.02 0.04 0.04 0.21 0.21
	8 bits	2.10 2.11 3.37 3.39 8.47 8.49
	Octahedral
	16 bits	0.02 0.02 0.04 0.04 0.23 0.23

Table V .

 V 3 -Maximum and average compression errors on the endpoints. Compression ratios of qfib and zfib. The N/A values are the ones for which the algorithm was not able to perform the compression and decompression.

	Stepsize	0.1 mm	0.2 mm	0.5 mm
	Nb. of bers	500k 3M 500k 3M 500k 3M
		Maximum error (×10 -2	)
	Quantization Precision	Deterministic
		8 bits	4.91 5.28 10.3 10.6 32.9 33.5
	Fibonacci
		16 bits	0.50 0.57 0.28 0.28 0.36 0.38
		8 bits	7.53 8.03 15.9 16.5 46.1 50.7
	Octahedral
		16 bits	0.50 0.56 0.26 0.29 0.47 0.46
	Quantization Precision	Probabilistic
		8 bits	2.76 2.85 5.25 5.43 18.4 18.6
	Fibonacci
		16 bits	0.13 0.13 0.12 0.14 0.64 0.63
		8 bits	4.78 4.78 8.30 8.64 26.7 31.8
	Octahedral
		16 bits	0.13 0.13 0.14 0.15 0.79 0.80
		Average error (×10 -3	)
	Quantization Precision	Deterministic
		8 bits	19.9 19.9 42.2 42.2 109 109
	Fibonacci
		16 bits	0.44 0.44 0.32 0.32 1.17 1.18
		8 bits	29.0 29.0 61.4 61.4 161 161
	Octahedral
		16 bits	0.44 0.44 0.34 0.34 1.31 1.31
	Quantization Precision	Probabilistic
		8 bits	17.5 17.5 25.1 25.3 62.5 62.7
	Fibonacci
		16 bits	0.21 0.21 0.41 0.41 2.12 2.12
		8 bits	28.0 28.1 40.7 40.5 94.4 94.7
	Octahedral
		16 bits	0.22 0.22 0.46 0.46 2.36 2.39
	Table V.4 -Stepsize	0.1 mm	0.2 mm	0.5 mm
	Nb. of bers	500k 3M 500k 3M 500k 3M
		Compression ratios (in percentage)
	Method Parameter	Deterministic
		8 bits	91.4 91.4 91.1 91.1 90.4 90.4
	qfib	16 bits	83.1 83.1 82.8 82.8 82.3 82.2
		same ★	N/A N/A 78.4 N/A 96.6 96.8
	zfib	0.2 mm	98.1 98.1 95.9 96.0 87.5 87.5
	Method Parameter	Probabilistic
		8 bits	91.4 91.4 91.2 91.2 90.8 90.8
	qfib	16 bits	83.1 83.1 82.9 82.9 82.6 82.6
		same ★	N/A N/A 78.1 N/A 87.1 N/A
	zfib	0.2 mm	96.0 N/A 88.7 N/A 69.9 N/A

Table V .

 V 5 -Computation times of qfib and zfib. With zfib, we set the maximal error to 0.2 mm. N/A are the values for which the algorithm was not able to perform the full compression and decompression.

	Stepsize	0.1 mm	0.2 mm	0.5 mm
	Nb. of bers	500k 3M 500k 3M 500k 3M
		Compression time (s)
		qfib ( bo)	24.1	144	12.8 74.8 5.49 32.0
	Deterministic	qfib (octa) 7.83 46.5 3.81 22.6 1.67 9.76
		zfib	702 4243 387 2284 387 2373
		qfib ( bo)	27.9	167	15.3 90.6 8.27 49.7
	Probabilistic	qfib (octa) 8.61 54.7 4.86 29.0 2.53 15.5
		zfib	910 N/A 1052 N/A 1418 N/A
		Decompression time (s)
		qfib ( bo)	4.98 30.1 2.61 15.4 1.14 6.79
	Deterministic	qfib (octa) 3.56 20.7 1.90 11.3 0.88 5.30
		zfib	12.1 72.7 12.9 77.1 17.3	103
		qfib ( bo)	5.77 34.9 3.23 18.9 1.75 10.3
	Probabilistic	qfib (octa) 4.08 24.3 2.24 13.3 1.29 7.47
		zfib	28.5 N/A 43.0 N/A 60.8 N/A

  Table V.7 -Errors of FA computation between original bers and the compressed and decompressed ones using qfib.

	Stepsize		0.1 mm	0.2 mm	0.5 mm
	Nb. of bers	500k	3M	500k	3M	500k	3M
		Average error of FA (in percentage)
	Method	Bits		Deterministic
	qfib	8	0.004 0.004 0.008 0.009 0.023 0.023
	bonacci	16	0.000 0.000 0.000 0.000 0.000 0.000
	qfib	8	0.006 0.006 0.012 0.012 0.033 0.036
	octahedral 16	0.000 0.000 0.000 0.000 0.001 0.000
	zfib	-	2.506 2.527 2.572 2.553 2.200 2.219
	Method	Bits			Probabilistic
	qfib	8	0.002 0.002 0.004 0.002 0.003 0.002
	bonacci	16	0.000 0.000 0.000 0.000 0.000 0.000
	qfib	8	0.003 0.003 0.005 0.003 0.004 0.006
	octahedral 16	0.000 0.000 0.000 0.000 0.000 0.000
	zfib	-	0.575 N/A 0.347 N/A 0.077 N/A
	Stepsize		0.1 mm		0.2 mm	0.5 mm
	Nb. of bers		500k 3M 500k 3M 500k 3M
		Maximum error (absolute value ×10 -2 )
	Quantization Precision		Deterministic
		8 bits	3.91 3.93 7.92 8.32 22.6 22.9
	Fibonacci
		16 bits	0.29 0.29 0.16 0.15 0.19 0.23
		8 bits	5.89 6.23 12.0 12.4 32.7 32.5
	Octahedral
		16 bits	0.30 0.29 0.15 0.15 0.25 0.28
	Quantization Precision		Probabilistic
		8 bits	2.43 2.25 3.66 3.93 11.1 12.4
	Fibonacci
		16 bits	0.09 0.09 0.09 0.09 0.44 0.48
		8 bits	3.48 3.85 5.61 6.39 17.6 17.8
	Octahedral
		16 bits	0.09 0.09 0.10 0.10 0.48 0.57
		Average error (absolute value ×10 -3 )
	Quantization Precision		Deterministic
		8 bits	0.45 0.37 0.99 0.88 2.69 2.61
	Fibonacci
		16 bits	0.01 0.01 0.01 0.01 0.04 0.03
		8 bits	0.67 0.53 1.45 1.31 3.95 3.83
	Octahedral
		16 bits	0.01 0.01 0.01 0.01 0.04 0.04
	Quantization Precision		Probabilistic
		8 bits	0.44 0.35 0.70 0.63 1.87 1.78
	Fibonacci
		16 bits	0.01 0.00 0.01 0.01 0.07 0.06
		8 bits	0.69 0.54 1.12 0.99 2.77 2.63
	Octahedral
		16 bits	0.01 0.01 0.01 0.01 0.08 0.07

Table V .

 V 8 -Compression times of our out-of-core algorithm (qfib) with an 8 bits octahedral quantization.

	Stepsize	0.1 mm	0.2 mm	0.5 mm
	Nb. of bers 500k 3M 500k 3M 500k 3M
		Compression time (s)
	Deterministic 29.1 166 14.3 88.3 6.82 40.8
	Probabilistic	33.0 192 17.9 103 9.99 57.9
	V.2.3.6 Limitations

Table V .

 V 9qfib block format used on the GPU.

	Element	Size
	First and second points	12 bytes
	Fiber ID	4 bytes
	Is beginning (Boolean)	1 bit
	Is ending (Boolean)	1 bit
	Number of compressed points	6 bits
	Block ID	1 byte
	62 compressed points	62 bytes
	Total	80 bytes

For more details on the MRI process and the physics involved in measuring di usion, we advise the reader to consult the website mriquestions (https://mriquestions.com).

= 1 + 2 accordingly. Every point in the extremity region is connected to the extremity of the other cylinder with the same weighted process. This results in a skeleton for 3 as shown in Figure III.5b. Then, the pro le of 3 is built using ellipses that will encompass 1 and 2 locally(Figure III.5c). The previous pro le of 1 and 2 is considered (not shown in the gure). To avoid ickering, the thickness of rendered cylinders is clamped to a minimum, including for original bers that

https://github.com/CorentinMercier/neural-meta-tracts

→ 0) whereas the latter does not. When 4 ≠ 0, the center c as well as the radius of the sphere can be retrieved using:

Please note that any gray matter segmentation pipeline could be used.

https://github.com/syrousseau/qfib

Remerciements

V Applications

In this chapter, we present some applications based on the contributions from previous chapters: the multi-resolution representations, geometric representations, and the knowledge on the bers. First, we propose a segmentation approach for automatically extracting some bundles from the tractograms. Then, a compression algorithm suited to bers is presented, as well as its use for on-the-y decompression on a GPU, allowing for the visualization of massive tractograms.

V.1 Segmentation

V.1.1 Introduction

The neural architecture of the white matter of the brain, obtained using tractography algorithms, can be divided into di erent tracts. Their function might be a ected in some syndromes or conditions and is in many cases still an object of study. Obtaining a reproducible and correct segmentation is therefore crucial both in clinics and in research. However, it is di cult to obtain due to the huge number of bers and high inter-subject variability. We propose to segment and recognize tracts by using anatomical de nitions, which are usually based on relations between structures (e.g. "anterior to", "close to", with respect to other anatomical structures). Since these de nitions are mainly qualitative, we propose to model their intrinsic vagueness using fuzzy spatial relations [Bloch, 1999, Bloch, 2005], as illustrated in Figure V.1. All qualitative relations are eventually summarized into a single quantitative "anatomical coherence score" (ACS) which is assigned to every ber of the tractogram. To cope with the high redundancy of tractograms and ease interpretation, we also take advantage of the simpli cation scheme based on the multi-resolution representation presented in Section III.2. This allows for a real-time ne-tuning of the threshold for better tract segmentation and eases the function (P) which assigns to every point P in space the degree of being in direction with respect to . In 2D (see Figure V.2), for any point P, let ( ì , ì ) be the angle between ì and ì , where r is a point in , and is the minimum angle among all points r ( = r∈ ( ì , ì )). We de ne: This means that (P) continuously decreases from 1 to 0, being equal to 1 when P is exactly in direction with respect to and 0 when it is in the opposite direction. Note that other decreasing functions could be used. In 3D, the direction is represented by two angles, and computing from Equation V.1 might be expensive. A more e cient approach was proposed in [Bloch, 1999] by expressing as a fuzzy dilation:

where is a fuzzy dilation by the fuzzy structuring element (P) = max(0, ( ( ì , ì ))) where O is the origin of space. Using this approach, we model the directions anterior, posterior, superior, inferior, right, and left as well as the relations "lateral" and "medial", which are commonly used in the neuro-anatomical literature. For these two last relations, we use as reference the mid-sagittal plane which is automatically detected using the method described in [START_REF] Tuzikov | Evaluation of the Symmetry Plane in 3D MR Brain Images[END_REF]].

In the literature, a white matter tract is usually described as a logic combination of several relations, using operators such as AND and OR. The proposed fuzzy models of spatial relations are combined using fuzzy AND (using t-norms) and fuzzy OR (using t-conorms). Here, we use the minimum for AND and the maximum for OR, computed voxel-wise. The complete list of de nitions is represented as an abstract syntax tree, that is collapsed in a hierarchical order to eventually compute the application scenario.

This algorithm is applied independently on each ber of the dataset. As the rst unit vector is recomputed for each ber, it allows for a di erent for each ber, making qfib usable for resampled bers, as long as it guarantees constant stepsize along each ber. The C++ pseudo code for the compression and decompression is available in Algorithm 6. The pseudocode of the mapping and inverseMapping functions can be found in [START_REF] Rousseau | Fast Lossy Compression of 3D Unit Vector Sets[END_REF]. The quantize and unquantize functions are the octahedral quantization [START_REF] Meyer | On Floating-Point Normal Vectors[END_REF] or the spherical Fibonacci [START_REF] Keinert | Spherical Fibonacci Mapping[END_REF]. For more details on the implementation, the reader can refer to the source code available on GitHub (https://github. com/syrousseau/qfib). When used for le storage, we de ne the qfib format in Figure V.15. We evaluate it in Section V.2.3.

V.2.3 Results and discussion

We test our algorithm using bers traced from a randomly selected subject of the HCP dataset [Van Essen et al., 2012]. As our method requires a constant stepsize, we use the SD_STREAM (deterministic) and iFOD1 (probabilistic) algorithms from MRtrix3 [START_REF] Tournier | MRtrix: Di usion tractography in crossing ber regions[END_REF] to compute the tractograms. We test the e ectiveness of our algorithm on 12 bundles where we vary the stepsize between 0.1 and 1 , and the number of bers between 60 and 3 . The maximum angle is not speci ed when generating the bers, which means that the default formula of MRtrix3 [START_REF] Tournier | MRtrix: Di usion tractography in crossing ber regions[END_REF] is used: = 90 • × /

. The voxels size in the dMRI images is 1.25

3 . Fibers lengths are constrained between 40 and 256

and are forced to end on the brain cortical surface. All con gurations are presented in Table V.1 with le sizes (tck format).

Our technique is general and any di usion model or tractography algorithm could be used provided it is based on a per-ber constant stepsize ( ). We implemented a tck le reader but any le format for 3D points can be handled, assuming a suitable reader is provided to access a list of point coordinates. Although its application is demonstrated with the octahedral and spherical Fibonacci quantizations, our method is compatible with any actual or future unit vector quantization method.

A

Smoothness of the local partition of unity

This appendix is linked to Section IV.2.2.3 and is the proof of the smoothness of the proposed function for the local partition of unity. These proofs were made by Pierre Roussillon and Jean-Marc Thiery.

Let us show that Ω admits null derivatives of any order in 0 and 1 -.

It is clear, by recurrence, that derivatives of all orders will be composed of products of (-(1/( -1) ) / 2 ) / and

(1/( -1) )

( -1) .

For = 0, (-(1/( -1) ) / 2 ) dominates ( -∞ tends quickly to 0) over 1/ , and thus ( )

For = 1 -, (1/( -1)) dominates ( -∞ tends quickly to 0) over 1/( -1) , and thus

Let us now show that ( ) is smooth everywhere.

For a point outside and inside , we have that:

( ) = Ω( ( ), ( )), with Ω( , ) = Ω ( /( -)), ( ) (resp. ( )) denoting the signed distance from to the sphere (resp. ).

It follows that

By recurrence, it is straightforward to see that derivatives of any order of ( ) are products of ( ) Ω , which go to 0 as tends to ( ( ) Ω (0) = 0 ∀ > 0) or when tends to ( ( ) Ω (1 -) = 0 ∀ > 0). The derivatives vanish therefore quickly enough to match the null derivatives of ( ) at the boundary of the domain ( ( ) is constant for ∈ and for ∉ ).

Additionally, it is clear that the derivatives of the term on the right side are all well-posed: The only point where this scheme could result in a non-di erentiable partition of unity is at the center of S (where ( ) is continuous only), but by construction this point is strictly inside S and thus outside the domain of utility of Ω ( ( ) < 0 and ( ) = 0 here).

Finally, note that this derivation would remain valid for any other function Ω exhibiting the mentioned properties.

List of Figures

List of Algorithms