E Cohen
email: cohen.emm@gmail.com

T Deffieux
email: thomas.deffieux@espci.fr

C Demené

L D Cohen
email: cohen@ceremade.dauphine.fr

M Tanter

4D point cloud registration for tumor vascular networks monitoring from ultrasensitive Doppler images

. In this study, we propose new computational tools to monitor the growth of a tumor vascular network by registering in time and space this new highly sensitive temporal data. We first quantify the acquired data using the minimal-path based framework we introduced in [4]; the vascular network paths around the tumor are segmented from images obtained for four days of observation; local geometrical parameters such as diameters are also estimated. Then, using a point cloud representation of the segmented vascular networks, we develop point cloud registration algorithms that automatically align similar vascular structures, thus allowing a better visualization of the growth and the evolution of the tumor vascular network. A rigid registration model is first considered by manually selecting similar features from two temporal different observations of the tumor. More accurate results are then obtained by automatically extracting invariant vascular patterns. Finally, combining rigid transformations to non-linear deformation models produce a very accurate time matching between invariant vascular structures.

1 Introduction

Context

Medical ultrasound imaging counts among the most used clinical imaging techniques, very appreciated for its portability, real time working and low-cost. It provides an anatomical imaging of high quality, usually associated to a Doppler examination for blood flows observation and quantification. In the last few years, ultrafast ultrasound imaging [START_REF] Tanter | Ultrafast imaging in biomedical ultrasound[END_REF] revolutionized the capacity to observe very fast physiological variations, such as shear waves propagation in the body. In fact, such waves propagate at typical speeds of a few tens of m.s -1 , involving very high frame rates ∼1 kHz (1000 images/sec) to be detected. Instead of using focused ultrasounds and forming the final image column by column, ultrafast imaging sends an ultrasonic plane wave that propagates and reflects into the whole image plan. Thus, echoes of all the structures back-propagate in the same time, and the final image can be reconstructed in only one ultrasonic emission. Ultrasensitive Doppler is one of these new ultrasensitive techniques which allows high sensitive acquisition of small vascular features without contrast agent.

Beyond anatomical vascular imaging, ultrasensitive Doppler allows for the first time to perform functional imaging with ultrasound (functional ultrasound or fUS). Mace et al. [START_REF] Macé | Functional ultrasound imaging of the brain[END_REF] demonstrated a first proof of concept by observing significant variations of blood flux in the rat brain, following a stimulation of the rat whiskers or an epileptiform seizures.

Being tens of times more sensitive than conventional Doppler [START_REF] Mace | Functional ultrasound imaging of the brain: theory and basic principles[END_REF], ultrasensitive Doppler allows very fine imaging of blood flows. Using a 15 MHz linear ultrasound probe, typical 2D spatial resolution obtained in the plane of imaging is 100 µm 2 . Yet, resolution in the out-of-plane direction remains higher, typically equal to 500 µm for rodent brain imaging. Thus, Demene et al. [START_REF] Demené | 4d microvascular imaging based on ultrafast doppler tomography[END_REF] proposed an ultrafast Doppler tomography based system (UFD-T) able to produce 3D microvascular maps of rodent brain, with high spatio-temporal resolution of 100 µm 3 × 10 ms. The acquisition is achieved in vivo via a simple tomographic mechanical set-up. UFD-T Fig. 1 Ultrasensitive Doppler imaging of the rat brain. Left: 2D ultrasensitive Doppler coronal plane. Right: 3D tomographic reconstruction UFD-T [START_REF] Demené | 4d microvascular imaging based on ultrafast doppler tomography[END_REF].

scanning times, ∼20 minutes for the whole rodent brain surface, are slightly below MRI and CT performances. For real time applications, other approaches using probes with multidimensional arrays of transducers could be more advantageous (e.g. [START_REF] Provost | 3d ultrafast ultrasound imaging in vivo[END_REF]). However, due to high temporal resolution of ultrafast imaging, 3D reconstructed volumes by UFD-T contain a large quantity of temporal information, revealing how promising is the technique in the context of 4D imaging. Nowadays, with the new design of matrix-array ultrasonic probes, real 3D ultrafast imaging can be achieved in vivo within a single ultrafast acquisition. It was demonstrated for 3D shear-wave imaging, 3D ultrafast Doppler imaging [START_REF] Provost | 3d ultrafast ultrasound imaging in vivo[END_REF] or functionnal imaging in rodents [START_REF] Rabut | Full 4d functional ultrasound imaging in rodents using a matrix array[END_REF].

Motivations

Ultrafast ultrasound has introduced a new way of high resolution vascular imaging that benefits of all the convenient advantages of ultrasonic devices (real time, portability, low cost). Numerous of applications are currently in development with the production of new available data. This imaging data could be then analyzed by the implementation of new algorithms in the field of image analysis and computer vision.

In the past, Folkman [START_REF] Folkman | Tumor angiogenesis: therapeutic implications[END_REF] has already explained that tumor growth is deeply connected to angiogenesis, the creation of new blood vessels responsible for the provision in nutriments and oxygen to the tumor. Many cancers are treated with antiangiogenetic strategies acting on tumor vasculature [START_REF] Jain | Normalization of tumor vasculature: an emerging concept in antiangiogenic therapy[END_REF]. Monitoring of tumors and their surrounding vascular networks is thus essential to understand tumor growth and architecture.

Vascular networks can be very complex structures with many different type of vessels. In particular, tumors are characterized by highly tortuous vessels that can appear spontaneously during the tumor growth. The variability in time of those networks makes also their analysis very challenging.

UFD-T enables 3D imaging of tumors microvasculature. Demene et al. [START_REF] Demene | Comparison of tumor microvasculature assessment via ultrafast doppler tomography and dynamic contrast enhanced ultrasound[END_REF] have imaged the temporal evolution of a tumor implanted in the back of mice. Such impressive data need rigorous numerical analysis to be quantified. Thus, we used our minimal path based framework [START_REF] Cohen | An isotropic minimal path based framework for segmentation and quantification of vascular networks[END_REF] to segment vascular networks surrounding the tumor and estimate local geometrical parameters (diameter, curvature, etc.). Such extraction of the vasculature allowed us to perform 3D registration of the temporal data, leading to more visual understanding of the tumor evolution and a better description of local variations. We implemented a 3D reconstruction and rendering algorithm showing the registered tumor data and its evolution (see figure 2). We also investigated non-linear models to match similar tumor vascular features and try to detect temporal repeatability despite the disordered evolution of a tumor vascular network.

Fig. 2 Tumor growth reconstruction and visualization over four days of observations after implantation in a mouse back using our segmentation and registration pipeline.

2 Material and mathematical background

Material

The data used for this study was acquired by UFD-T on four mice, in the flank of which were implanted subcutaneously 10 mm 3 tumor fragments, coming from mouse Lewis lung carcinoma cells. More details can be found in [START_REF] Demene | Comparison of tumor microvasculature assessment via ultrafast doppler tomography and dynamic contrast enhanced ultrasound[END_REF]. The data has been jointly acquired on the four mice on days 8, 12, 16 and 20 after the implantation of the tumor.

We would like to understand how the tumor evolves in time from days 8 to 20. The growth of the tumor causes a global proliferation of the vascular network that is well observable, but a local comparison of the networks is hard without the use of other methods.

Rigid point cloud registration

Point cloud registration is a well-known issue in the field of robotics, computer vision and graphics, and the medical field e.g. for image-guided interventions [START_REF] Peters | Image-guided interventions: technology and applications[END_REF]. The underlying mathematical problem is to find the best geometric transformation describing the relationship between two sets of points or point clouds. The first question to answer before solving the problem is how far the two point clouds are from each other; in other words, how to model the geometric transformation between them? A first approach is to use a rigid, or linear, displacement model involving in general three kinds of operations: translation, rotation and scale. Let us note that there may be no exact solution to this problem, in particular when data are complicated involving non-linear deformation. Yet, it is still interesting to obtain a preliminary approximation that globally match the points, and further envisage a more refined model.

The Procrustes problem

Let X = (x i ) and Y = (y i ) be two 3D point clouds of respectively N x and N y points. We will consider that X is the moving data set and Y the static model, and will search for the best rigid transformation T such as Y ∼ T (X). A simple case, called Procrustes problem, is when the correspondences between points is known, i.e. when N x = N y = n and each point x i corresponds to the same indexed point y i . Schonemann et al. [START_REF] Schönemann | Fitting one matrix to another under choice of a central dilation and a rigid motion[END_REF] and Horn [START_REF] Horn | Closed-form solution of absolute orientation using unit quaternions[END_REF] have both proved least-square solutions of this problem. In what follows, we sum up the main theoretical steps. The idea is to find the transformation minimizing the sum of square errors between points. The least-square functional to be minimized is

F(t, R, s) = n ∑ i=1 y i -sRx i -t 2 , (1)
where the rigid model parameters that must be found are t, R and s, respectively denoting a translation vector, a rotation matrix and a scale scalar parameter. Let us define the centroids of each point cloud by

x = 1 n n ∑ i=1 x i , ȳ = 1 n N x ∑ i=1 y i , (2)
and denote new demeaned coordinates

x i = x i -x, y i = y i -ȳ. (3)
The minimization function can be modified as follows

F(t, R, s) = n ∑ i=1 y i -sRx i -t 2 , (4)
with t = tȳ -sR x. Expanding the last equation gives

n ∑ i=1 y i -sRx i 2 -2 t , n ∑ i=1 y i -sRx i + n t 2 (5)
By noting that ∑ n i=1 x i = ∑ n i=1 y i = 0, it comes that the middle term of this expression is zero. This simplification of F allows to easily conclude that t should be equal to zero to minimize F. Indeed, the first term does not depend on t and the third one is a square so is positive. Thus, we have already found the translation

t = ȳ + sR x (6)
To determine the scale, we keep on expanding

F n ∑ i=1 y i -sRx i 2 = n ∑ i=1 y i 2 -2s n ∑ i=1 y i , Rx i + s 2 n ∑ i=1 Rx i 2 (7)
Since the rotation R is an orthogonal operator, it preserves the euclidean norm, allowing to write Rx i = x i . We then rewrite the last expression as follows

S y -2sD + s 2 S x , (8)
where

S y = ∑ n i=1 y i 2 , S x = ∑ n i=1 x i 2
, and D = ∑ n i=1 y i , Rx i . Finally, by factorizing this quadratic equation in s as

s S x - D √ S x 2 + S y S x -D 2 S x , (9)
and noting that S y S x -D 2 ≥ 0 (CauchySchwarz inequality), the first term should be equal to zero to minimize F and we get the optimal value that s should take

s = D S x = ∑ n i=1 y i , Rx i ∑ n i=1 x i 2 . ( 10
)
The determination of the rotation is more complicated and a prove can be found in [START_REF] Schönemann | A generalized solution of the orthogonal procrustes problem[END_REF] using the singular value decomposition (SVD) method, or in [START_REF] Horn | Closed-form solution of absolute orientation using unit quaternions[END_REF] using the quaternion representation of rotations. We give here the SVD-based result. Defining the cross-covariance 3 matrix

B = n ∑ i=1 y i x i T , (11)
and realizing its SVD decomposition B = UΣV T , where U, V are two unitary matrices and Σ a diagonal matrix, the rotation is given by

R = Udiag([1, 1, det(UV )])V T . (12)
The operator diag transforms an input vector into a diagonal matrix whose diagonal is composed of this vector components. An other expression of s using Σ can also be found in [START_REF] Schönemann | Fitting one matrix to another under choice of a central dilation and a rigid motion[END_REF] s = trace(Σ )

∑ n i=1 x i 2 . ( 13
)
Considering that s is a scalar, the above problem is called the isotropic Procrustes problem. Several authors have proposed to solve the anisotropic Procrustes problem which can be express as the minimization of the following least-square functional

F(t, R, s) = n ∑ i=1 y i -RSx i -t 2 , ( 14
)
where S is a diagonal matrix modeling the scale anisotropy. The last expression refers to the pre-scaling problem where scaling S precedes the rotation R, while the post-scaling problem minimizes

F(t, R, s) = n ∑ i=1 y i -SRx i -t 2 . ( 15
)
To solve the pre-scaling problem, Gower et al. [START_REF] Gower | Procrustes problems[END_REF] used an iterative procedure called block relaxation (BR). Each iteration has two steps. First, for fixed R the computation of the three diagonal coefficients of S

λ i = (B T R) ii ∑ n i=1 x i x i T ii , (16)
then for fixed S, the computation of R

R = UV T , ( 17
)
where UΣV is now the singular value decomposition of BS. The algorithm iterates until the change in least-square error falls below a certain threshold. The authors used a similar BR algorithm for the post-scaling problem. Dosse et al. [START_REF] Dosse | Anisotropic orthogonal procrustes analysis[END_REF] pointed out that the BR procedure can lead to local minima corresponding to incorrect solutions. They proposed an other solution less sensitive to this kind of problem. In this study, we will always refer to the pre-scaling anisotropic problem.

Iterative Closest Point algorithm

Let us now assume that the correspondence between points is unknown. The problem becomes much more complicated. Indeed, obtaining an exact solution by solving a least-square problem is not possible anymore because we do not know the pairs of points whose square errors should be minimized. Thus, Besl et al. [START_REF] Besl | Method for registration of 3-d shapes[END_REF] introduced a new iterative method, called iterative closest point (ICP), to register 3D point clouds which attempts to match the points in X and Y by using a closest point approach.

Let us describe an iteration k of the algorithm. We denote by X k the transformed point cloud issued from X at the k-th iteration. Firstly, the algorithm determines the closest points ỹi in Y to each x i ∈ X k , by minimizing the euclidean distance ỹi = arg min y∈Y yx i [START_REF] Provost | 3d ultrafast ultrasound imaging in vivo[END_REF] Note that, as soon as a point in Y is designated as a closest point, it cannot be selected again. As a consequence, it should be imposed that N y ≥ N x , for instance by always choosing the moving point cloud X to be the one with the fewest points. Then, the algorithm solves a Procrustes problem to register X k on the set of closest points Ỹk . The Procrustes problem can be isotropic or anisotropic as proposed by [START_REF] Chen | Registration of 3d shapes under anisotropic scaling[END_REF] for the registration of medical ultrasound images. The rigid transformation found is finally applied to X k to obtain the next iteration X k+1 . The algorithm stops the iterations when the change in least-square error falls below a threshold parameter.

Deformable models

When dealing with real data, especially with tumor data where vascular networks may greatly vary in time, geometrical transformations that map corresponding points can be very complicated. Therefore, deformable models are usually preferred to get a more precise estimation of the transformation. A deformable model aims to find a local transformation. Minimizing an energy functional, a set of parameters is estimated locally for each pair of corresponding points, allowing to map the points with many degrees of freedom and making the transformation more realistic.

A very famous deformable model is the Thin-Plate-Spline (TPS) transformation introduced by [START_REF] Bookstein | Principal warps: Thin-plate splines and the decomposition of deformations[END_REF]. Given two sets of n corresponding points X = (X i ) and X = (X i ), the method tries to find the best mapping function f : X → X minimizing the following energy

R 2 ∂ 2 f i ∂ x 2 2 + 2 ∂ 2 f i ∂ x∂ y 2 + ∂ 2 f i ∂ y 2 2 dxdy, (19)
in the case of 2D points, where i ∈ {1, ..., 2}, or

R 2 ∂ 2 f i ∂ x 2 2 + ∂ 2 f i ∂ y 2 2 + ∂ 2 f i ∂ z 2 2 + • • • 2 ∂ 2 f i ∂ x∂ y 2 + 2 ∂ 2 f i ∂ y∂ z 2 + 2 ∂ 2 f i ∂ x∂ z 2 dxdydz, (20)
in 3D, where i ∈ {1, ..., 3} and f i denotes the Cartesian coordinates of f . It can be shown that minimizers exist and take the following respective 2D and 3D forms

f 1 (x i , y i ) = a 0 + a 1 x i + a 2 y i + n ∑ j=1 F j r 2 i j ln(r 2 i j ) f 2 (x i , y i ) = b 0 + b 1 x i + b 2 y i + n ∑ j=1 G j r 2 i j ln(r 2 i j ), (21)
and

f 1 (x i , y i , z i ) = a 0 + a 1 x i + a 2 y i + a 3 z i + n ∑ j=1 F j r i j f 2 (x i , y i , z i ) = b 0 + b 1 x i + b 2 y i + b 3 z i + n ∑ j=1 G j r i j f 3 (x i , y i , z i ) = c 0 + c 1 x i + c 2 y i + c 3 z i + n ∑ j=1 H j r i j , (22)
for all points i ∈ {1, ..., n}, and where r i j = X i -X j . Inspired by [START_REF] Fitzpatrick | Image registration[END_REF], let us give some technical details in the case of 2D TPS about how the unknown parameters

a 0 , a 1 , a 2 , b 0 , b 1 , b 2 , F 1 , • • • , F n and G 1 , • • • , G n can be computed.
If those parameters are determined, we found a solution to the corresponding minimization problem. If we introduce the following matrices and vectors

P =    1 x 1 y 1 . . . . . . . . . 1 x n y n    , (23)
K =      0 r 2 12 ln(r 2 12 ) • • • r 2 1n ln(r 2 1n ) r 2 21 ln(r 2 21 ) 0 • • • r 2 2n ln(r 2 2n ) . . . . . . . . . . . . r 2 n1 ln(r 2 n1 ) r 2 n2 ln(r 2 n2 ) • • • 0      , (24)
L = P K 0 P T (25)
w 1 =          a 0 a 1 a 2 F 1 . . . F n          , w 2 =          b 0 b 1 b 2 G 1 . . . G n          , ( 26
) x =          x 1 . . . x n 0 0 0          , y =          y 1 . . . y n 0 0 0          , ( 27
)
the system of equations 21 can be expressed in terms of matrices as follows

x = Lw 1 y = Lw 2 . ( 28
)
Notice that the last three lines of the system correspond to 6 additional equations given by n ∑

1

F i = 0, n ∑ 1 G i = 0 n ∑ 1 x i F i = 0, n ∑ 1 x i G i = 0 n ∑ 1 y i F i = 0, n ∑ 1 y i G i = 0. ( 29
)
They can be interpreted as some boundary conditions that force f to be an affine transformation in the case x and y tend to infinity. The affine part of the transformation is parametrized by the a i , b i coefficients. Thus, by inverting the matrix L, we can derive all the unknown parameters and finally solve the problem

w 1 = L -1 x w 2 = L -1 y . ( 30
)
3 Methods and results

Once vascular networks have been segmented, they can be used to accurately register 3D images. It is clear that a sparse representation of vascular networks offers many more accurate possibilities of quantification and registration. Numerous works can be found in the literature dealing with vascular structures registration, e.g. [START_REF] Groher | 2d-3d registration of vascular images[END_REF][START_REF] Dufour | Segmentation et modélisation des structures vasculaires cérébrales en imagerie médicale 3d[END_REF][START_REF] Reinertsen | Validation of vessel-based registration for correction of brain shift[END_REF].

The model used to represent vascular networks is essential for the registration to be succesful. We proposed a representation in graph structure as described in [START_REF] Cohen | An isotropic minimal path based framework for segmentation and quantification of vascular networks[END_REF]. This graph representation is very useful to understand how the network is organized, and to decompose in a logically and organized order the segmented structures. However, this model is somewhat limited by the errors introduced by the algorithm of segmentation. Due to noise, complexity of data, or lack of precision, robustness in the method, there should be aberrant bifurcations or edges in the resulting segmentation. Therefore, a more general representation in point clouds should be associated to the graph-based model. For registration of point clouds, our work is based on a modified version of ICP and inspired by [START_REF] Reinertsen | Validation of vessel-based registration for correction of brain shift[END_REF].

The methods are presented and applied on 3D data acquired on rodent brains and tumors. They have been developed to solve several issues: the detection of vascular invariant structures in the brain for the production of cerebral atlas; the detection of invariant structures in the vasculature of tumors for a better understanding of their architecture; the visualization of tumors growth.

Rigid initialization

Our goal is to realize the temporal registration of the UFD-T 3D tumor images. Here, we will start finding only rigid transformations between the data, being an initial guess into our global registration process. For each day of observation, we apply the segmentation framework presented in [START_REF] Cohen | An isotropic minimal path based framework for segmentation and quantification of vascular networks[END_REF], and obtain a complete description of the vascular network including centrelines, bifurcations, and local diameters (figure 3). We also get the corresponding local blood flow intensities for each point of the segmented network, by interpolating the values on the original image.

Fig. 3 Segmentation of the tumor data for each day of observation. The color map corresponds to diameters, with large and small vessels respectively in red and blue.

Manual selection of feature points

Now visualizing the segmented networks, it is possible to recognize some vascular features which seem to repeat themselves in the different days of observation of the same tumor. They generally correspond to vessels with high diameters and high blood flows. Our assumption is that it should always exist such time-invariant vessels that feed the tumor and from which the tumor builds a new multitude of small and tortuous vessels. A first attempt to register two vascular networks is to select manually on each network few feature points characterizing the shape of the recognizable vessels, and try to find the rigid transformation between each point set. In fact, we select pairs of feature points corresponding to the same approximate positions on each vascular network. This allows us to find the translation, rotation and scale parameters by resolving a Procrustes problem with 4 pairs of points. Figure 4 shows the segmented networks of two days of observation of the same tumor, respectively in green (day 8) and green (day 12). Four points are manually selected on each data set around the shape of interest (in bold on the figure). Once the rigid transformation is computed, the transformation is then applied to all the moving network. We observe that the shapes in bold do not exactly overlap but globally the networks have been brought closer to each other.

If we apply the same process to each day of observation of the same tumor, trying now to register day 16 on 12 and day 20 on 16, we obtain a global visualization of the tumor growth shown by figure 5. The dark points are the feature points manually selected on each vascular network.

Automatic detection of features

The previous registration based on manual selection of feature points is clearly not enough accurate, and one must take benefit of more feature points or shapes. Thus, we turn to a point cloud representation of the segmented vascular networks. The major question to answer is how to reject non-relevant points and keep only repeatable vascular patterns that will make the registration successful. Besides, the correspon- dences between points are now unknown, so we should use the ICP algorithm to iteratively match the selected points and converge towards an acceptable solution.

Using ICP directly on the whole point clouds, the presence of outliers and vascular structures of high temporal variability may cause similar or even worse results than those of the previous manual selection model. Even if such result has the advantage to be an automatic procedure without manually selecting points, we must improve it in accuracy. Therefore, we pre-select the points used for ICP by choosing only those with high diameters and blood flow intensities. This is based on the assumption that tumors grow from an invariant set of large vessels. We order points according to diameter and blood flow intensity and only keep the union of points with the highest values. It is remarkable that considering vessel segments instead of points leads to better ICP performance, probably due to the continuity property of a segment allowing the algorithm to better identify vascular structures.

Thus, we apply the ICP algorithm on the pre-selected vessel segments with a nonrigid transformation adding a scale parameter to the model. The result, presented in figure 6, gives the great impression that the data could be finely aligned and proves that only four observations allow to describe the tumor growth dynamic. Some vascular structures, whose the repeatability is easily recognizable on the different days of observation of the tumor, are manually highlighted in bold. After registration, we observe that those structures are mainly situated in the same area, with around them the proliferation of a big network of small and more unstable vessels.

Automatic detection of local invariant patterns

The previous method remains particularly unstable with the possible detection of false local minima in the registration point cloud algorithm. The registration should be much more improved. In order to accurately align the repeatable structures, we propose a more local approach. In fact, the pre-selected segment features previously used from each tumor data can be partitioned into sets of connected segments, enabling to locally register smaller vascular shapes one on the other instead of performing the global registration of several shapes. To detect if two segments are connected, we take advantage of the graph representation used in [START_REF] Cohen | An isotropic minimal path based framework for segmentation and quantification of vascular networks[END_REF] for the segmentation of vascular networks. When two segments E 1 and E 2 have a node in common, we connect and include them into a subgraph

G = i E i (31)
We continue to fill each subgraph until no more pairs of common nodes can be detected. Figure 7 shows on the left the pre-selected vascular segments from a oneday-observation tumor data. We observe the presence of several groups of connected segments, especially a large one looking like half of an ellipse. This large vascular structure corresponds to a large vessel repeated on each day of observation of the tumor. The other structures seem to be more variable and may not be so much useful for the alignment of the different tumor data. On the right of the same figure, the connected components are separated by individual colors and will be now analyzed separately. We propose to identify among those subsets of vessels the most invariant structures when we observe the tumor along its days of growth. For this purpose, we must find some features characterizing the components that can be best matched as for instance the observable semi-elliptic pattern. As a first assumption, we claim that the longer and the largest a pattern G is, the more invariant it will be. Therefore, we use as a first feature the total curve length of the connected components and as a second one, a simple estimator of the space extension of the pattern Max Distance(G) = max

(E i ,E j )∈G m(E i ) -m(E j ) , (32)
where m(E i ) is the geodesic middle of the segment E i . By plotting on a graph, figure 8, those two features for each of the four days of observation of the same tumor (one color corresponds to one day of observation, data 1 to 4), we observe that one specific component designated by a gold star can be linearly separated from the rest of the patterns. This specific component is exactly the half elliptic form mentioned above. This method allows us to detect the most proper component in each day of observation and use it to register all the data. The process of pre-selection of high diameters and blood volume and the detection of the best characterizing component is shown on figure 9. At the end, we obtain four components looking very similar. We then apply the ICP algorithm with a rigid transformation to register those four vascular patterns. The rigid transformation found by the algorithm is then applied to all the data points, leading finally to a global rigid registration of the four days of observation of the tumor (figure 10). The result is much more accurate than the previous attempts, and the growth of the tumor can be properly visualized.

Non-linear registration

Once initial rigid transformations have been found, the different temporal data can be visualized in the same coordinates system and the growth of the tumor is more clearly depicted (figure 10). It is now reasonable to think about a non-rigid registration step, in order to match similar vascular structures with higher accuracy. This will help us to understand how the tumor is changing locally during its growth, improving the detection of invariant or unstable vascular areas.

For this purpose, we looked for a TPS deformable transformation that would best model tumor evolution. The TPS operator outputs a non-linear mapping between the matching points of two different point clouds. Since the correspondence between matching points is a priori unknown, we must incorporate the TPS model into an ICP procedure. Thus, we apply a modified version of the ICP algorithm in which the Procrustes problem solving step is replaced by the computation of the TPS transformation.

The result is very satisfying as many vessels are accurately matched. Figure 11 shows a tumor vascular network at two different days of observation. On the left, the networks have only been aligned according to the initial rigid transformation computed in the last section, but they are not completely registered yet. On the right, the non-linear component of the transformation, including scaling and other local deformation, is found using the TPS modified ICP version. In black are represented the TPS transformed points of the green vascular network (the green is the moving point cloud while the red remains static) whose distance from the red network is lower than one pixel. For those points, the matching seems to be very accurate.

To be sure that vascular patterns are consistently matched, we plotted figure 12 a color map showing the similarity between the two networks found by TPS. This map can be easily obtained by plotting on the left graph the inverse coordinates (i.e. their original coordinates before TPS transformation) of the best matching points colored in black on figure 11, and on the right the corresponding matching points of the red network. We can observe that many vascular patterns owning the same color in both maps own also very similar shapes. This confirms at least qualitatively that the TPS model performs well.

3D reconstruction of tumor growth

Once tumor data has been acquired at several different times, one difficult task is to adopt the proper display to understand the temporal dynamic of the data. The acquired images are a priori not properly registered, therefore it is hardly possible to identify and follow by eye the expansion and development of the tumor. Our rigid registration algorithm allows us to align the data in the same coordinate system and then to create a film displaying clearly the evolution of the tumor. Note that given the local diameter and the blood flow intensity at each point of the segmented vascular network, we can reconstruct the data in 3D using a color map corresponding to the blood flow. Figure 2 shows the 3D reconstruction of a tumor observed at 4 different days. Those images are screen shots of a film showing the data aligned and growing from a fixed point, allowing to clearly observe the expansion of the tumor, the proliferation of small unstable vessels, and the invariance of several big vascular patterns. The part of the tumor represented in blue is actually the real tumor tissue that was segmented from the same ultrasensitive Doppler acquisition using B-mode as explained in [START_REF] Demene | Comparison of tumor microvasculature assessment via ultrafast doppler tomography and dynamic contrast enhanced ultrasound[END_REF].

Conclusion

In this study, we have developed advanced numerical methods for the analysis of new medical ultrasensitive Doppler data. We used our segmentation framework [START_REF] Cohen | An isotropic minimal path based framework for segmentation and quantification of vascular networks[END_REF] as a basis for the development of tumor data monitoring in 3D. Extracted tumor networks are represented as point clouds that can be automatically registered using rigid and non-rigid transformations.

For the detection of vascular brain invariant structures and the possible creation of a brain vascular atlas, our framework of segmentation represents an important tool to extract the geometry parameters and get a fine description of vascular networks. Using unsupervised machine learning techniques could be a very interesting way to cluster groups of segmented vessels and automatically detect invariant patterns. Further, by acquiring a significant amount of data, one could build a statistical atlas by discriminating the invariant from the variable parts of the networks. The nonlinear registration methods used in this study are a good start to extract the invariant part of a network by selecting the best matching structures.

Finally, supervised machine and deep learning methods would be of great interest if we could build large labeled datasets from microvascular ultrasensitive Doppler data. Techniques of data augmentation could improve performances if the number of data is not sufficient and one could train a model to automatically detect repeatable vascular shapes on new data.

Fig. 4

 4 Fig. 4 Registration of tumor temporal data by manual selection of features. Two days of observations of the same tumor, days 8 and 12, are respectively in red and green. Left: four pairs of matching points are manually selected and linked with dot lines on the figure. Right: rigid registration of the red (moving) network onto the green (static) network, computing the transformation using the matching pairs of points in (a).

Fig. 5

 5 Fig. 5 Tumor growth visualization using rigid registration with manual selection of features. The features points are in dark. Observation days 8, 12, 16 and 20 of the tumor are respectively in green, red, cyan and blue.

Fig. 6

 6 Fig. 6 Non-rigid ICP registration with scale parameter using automatic pre-selection of features. From top to bottom: three different viewing angles of the registered tumor data. Observation days 8, 12, 16 and 20 of the tumor are respectively in green, red, cyan and blue. Some invariant vascular patterns were manually annotated in bold to simply evaluate the quality of the registration.

Fig. 7

 7 Fig.7Decomposition of selected features into connected sub-graphs of paths that can be processed separately. Left: selected features. Right: each connected sub-graphs of paths is plotted with a particular color.

Fig. 8

 8 Fig. 8 Classification of sub-graphs according to their max-distance values. Each temporal data has its own color. Each point in the graph corresponds to a connected sub-graph. The points annotated with a star are the invariant detected structures, presenting highest max-distance.

Fig. 9

 9 Fig. 9 The proposed registration pipeline. Observation days 8, 12, 16 and 20 of the tumor are respectively in green, red, cyan and blue

Fig. 10

 10 Fig. 10 Best invariant structures rigid registration. Left: similar connected components are automatically extracted from each days of observation of the tumor, and finally registered. Right: the transformation found is applied to all the point clouds in order to nicely visualize the tumor growth.

Fig. 11

 11 Fig. 11 Non-linear ICP registration with TPS. Left: unregistered tumor data. Right: registered tumor data using TPS. The green point cloud is the moving one. Black points are the transformed points whose distance from the red network is lower than one pixel. Those points most probably belong to vessel invariant structures.

Fig. 12

 12 Fig. 12 Matching of vascular structures from temporal tumor data. Left: observation day 8. Right: observation day 12. Matching structures are drawn with the same color on both graphs.