An automated workflow for the anatomo-functional mapping of the barrel cortex An automated workﬂow for the anatomo-functional mapping of the barrel cortex

Highlights: Background: The rodent barrel cortex is a widely used model to study the cortical processing of tactile sensory information. It is notable by the cytoarchitecture of its layer IV, which contains distinguishable structural units called barrels that can be considered as anatomical landmarks of the functional columnar organization of the cerebral cortex. To study sensory integration in the barrel cortex it is therefore essential to map recorded functional data onto the underlying barrel topography, which can be reconstructed from the post hoc alignment of tangential brain slices stained for cytochrome oxidase. New Method: This article presents an automated workﬂow to perform the registration of histological slices of the barrel cortex followed by the 2-D reconstruction of the barrel map from the registered slices. The registration of two successive slices is obtained by computing a rigid transformation to align sets of detected blood vessel cross-sections. This is achieved by using a robust variant of the classical iterative closest point method. A single fused image of the barrel ﬁeld is then generated by computing a nonlinear merging of the gradients from the registered images. Comparison with Existing Methods: This novel anatomo-functional mapping tool leads to a substantial gain in time and precision compared to conventional manual methods. It provides a ﬂexible interface for the user with only a few parameters to tune. Conclu-sions: We demonstrate here the usefulness of the method for voltage sensitive dye imaging of the mouse barrel cortex. The method could also beneﬁt other experimental approaches and model species.


Introduction
The rodent primary somatosensory cortex is a very convenient model for studying the cortical processing of sensory information because of its well defined structural and functional layout that is invariant from animal to animal [1,2,3]. In its layer IV, neurons are gathered into clusters called barrels that respect the same topology as the whiskers on the snout of the animal [4]. Each barrel is dedicated primarily to the processing of the input coming from its corresponding whisker (Figure 1A,B). When studying sensory processing in the barrel cortex either with electrophysiological or imaging methods, it is therefore of great interest to superimpose the recorded activ-ity onto the underlying barrel topography, which can be reconstructed from the post hoc alignment of tangential brain slices stained for cytochrome oxidase. In order to optimize this anatomo-functional mapping, which is usually accomplished manually, we developed an automated workflow for the registration of the histological slices of the barrel cortex and the 2-D barrel map reconstruction.
Here we focus our attention on voltage sensitive dye imaging (VSDI) of the mouse barrel cortex to illustrate the usefulness of the approach. However, the method can be extended to the study of the rat barrel cortex and applied to other techniques such as 2-photon calcium imaging.
The traditional first step to recover the map of the barrel cortex after imaging experiments is: brain fixation by perfusing the animal with a solution of paraformaldehyde, followed by the cutting of tangential slices (∼100 µm thick, with or without previous flattening of the cortex), which are subsequently stained for cytochrome oxidase using classical histological procedures that reveal the barrel arrangement in layer IV ( Figure 1C,D [5]).  Figure 1: Cytochrome oxidase staining of tangential sections from the mouse primary somatosensory cortex reveals the structural organization of layer 4 barrels that mirrors the arrangement of the vibrissae on the snout. A. Following the registration of tangential histological slices and the reconstruction of the barrel map, one can see that the spatial organization of the layer 4 barrels matches the layout of the vibrissae on the snout of the animal (B). C. Drawing of a coronal section of the left hemisphere of the mouse brain illustrating the position of layer 4 barrels within the primary somatosensory area of the cortex. After in vivo imaging of barrel cortex activity, sections are cut tangentially to reconstruct the layer 4 barrel map (cutting plane indicated by the red line). D. A series of tangential histological slices stained for cytochrome oxidase. On the first slice one can see superficial blood vessels. On the other slices, one can see white circular to elliptic spots that correspond to sections of plunging blood vessels. Depending on the exact axis of the cut, barrels can be spread over several slices.
Next, using digital microphotographs of the slices, it is necessary to: 1. register the slices ; 2. fuse the registered slices to define a reconstructed barrel image.
In this article we provide an automated solution for these two steps which are the most time-consuming tasks of the workflow when using conventional manual methods. After completing these steps, it is then relatively simple to define the barrel map by segmenting the reconstructed barrel cortex image. The superimposition of the map with the functional data can be finally achieved by using the superficial blood vessels as anatomical landmarks ( Figure 2). The proposed anatomo-functional mapping tool significantly speeds up the overall process and provides more accurate anatomo-functional mapping.

Registration of histological slices
A typical example of a series of images obtained after the histological process is shown in Figure 1D. Depending on their depth, the histological sections present different properties: in the first section, which corresponds to the surface of the cortex, some large superficial blood vessels are visible together with plunging blood vessels. This section is crucial for the whole process since it contains most of the superficial blood vessels that will be used for the final alignment of the histological data with the VSDI data ( Figure 2). The intermediate sections usually contain orthogonal blood vessels (white dots on the slices shown in Figure 1D) that can be used to align subsequent sections. Barrels start appearing on section 3 or 4 and can be visualized on up to 5 sections.
We do not use the deeper sections from subgranular layers as they do not contain any useful information for the reconstruction.
Image registration is a classical problem and its solutions find many applications in medical image analysis. Depending on the imaging modality and the specific prior knowledge of the object to register, a wide variety of methods have been considered in the literature (for an overview see [6,7]). Only a few previous studies explicitly deal with the problem of registering rodent brain histological sections, usually in order to reconstruct a whole brain in 3D [8,9].
In order to exploit the microscopic-scale information of the histological data, these applications require the precise registration of a large number of histological slices from the whole brain with a method that accounts for global but also local nonlinear deformations due to tissue shrinkage and tearing after histological preparation. Ourselin et al. presented in [8] a block matching strategy to compute local similarities and then estimate the rigid transformation that matches the maximum of similar regions in a robust way. Alternatively Ju et al. [9] used a method based on pairwise elastic imaging warps, with the specificity to compute the deformation of each section by considering not only two neighbors for each section, but an extended neighborhood including a group of images.
The specific problem of registering histological sections of Figure 2: Alignment of the barrel map with VSDI data using the superficial blood vessels as blueprints.
A. Barrel map reconstructed from the histological slices shown in figure 1D, using our anatomo-functional mapping method. B. Superimposition of the registered histological slices, allowing the delineation of the superficial blood vessels (in blue). C. The superficial blood vessels also appear in vivo on the fluorescent images taken during the VSDI session. The VSDI data can therefore be aligned with the underlying structural map of the barrel cortex using the blood vessels as anatomical landmarks. On the right, the cortical activity is shown imaged at 10ms following a single C2 whisker deflection with the voltage sensitive dye RH1691 under urethane anesthesia. The barrels outlined from the reconstruction in A are superimposed on the image as white lines.
the rodent barrel cortex has been rarely addressed in the literature. Egger et al. [3] proposed a tool for the 3-D reconstruction and standardization of the rat barrel cortex for the precise registration of single neuron morphology. Seventy micrometer thick tangential sections of rat barrel cortex are aligned pairwise by finding the rotation and translation that best superimpose blood vessels of two adjacent sections either manually or automatically, using a tool originally developed for the reconstruction of neuronal processes [10]. Finally the barrels are segmented by using a semi-automated method.
Here we developed a tool to compute the rigid transformations to align sets of detected blood vessels ( Figure 3A), and we decided to focus on a 2-D image fusion using an automated method. The usual approach to perform point cloud registration is the iterative closest point (ICP) algorithm introduced by Besl [11]. While the initial formulation is not robust to outliers, several approaches explicitly deal with this issue [12,13,14,15,16]. These approaches are related to reweighting least squares methods and we propose in Appendix B a unifying presentation and convergence analysis of a robust ICP method.

Fusion of histological slices
In our framework, the 2-D reconstruction of the barrel maps amounts to perform a fusion of the registered sections (Figure 3B). The goal is to reconstruct the edges of the barrels that are spread over the slices. Image fusion is classically addressed in the field of image processing and computer vision, to perform for instance image editing and stitching. To reconstruct sharp edges, it is necessary to use non-linear methods, the most popular one being based on gradient-domain blending [17,18]. The goal is to generate a novel image by locally keeping the content of the image having the highest local frequency. In this article, we focus on the use of gradient-domain methods, which have the advantage of being simpler and can be easily tuned for the purpose of barrel map fusion.

Contributions
Our main contribution is a comprehensive pipeline for the 2-D reconstruction of the barrel cortex from tangential histological sections. This pipeline is composed of 2 successive modules that respectively perform: histological section registration and barrel image reconstruction using image fusion. As a side contribution, we propose to recast the registration problem using a robust ICP optimization method, and we show that a Majorize-Minimize framework can be applied to provably ensure the convergence of the method. These contributions originate from specific needs raised by studies of sensory integration in the barrel cortex, however the histological section registration tool proposed here might be helpful to reconstruct any anatomical tissue in which blood vessels penetrate predominantly orthogonally to the cutting plane of the histological slices.

Overview of the proposed framework
This Section details each step of the proposed method. The corresponding Matlab code to reproduce the figures of this article is available online 1 . This code is packaged as a graphical user interface (GUI) that is helpful to guide the user through the various processing steps, from image loading to the final reconstructed barrel map. Table 1 reviews the notation introduced in the paper. Width of debiasing kernel Eq (6) 10 Table 2: Parameters used in the paper.

Step 1: Pairwise Registration of Histological Sections
The input of the algorithm is Q raw images which are digital images with range normalized in [0, 1] (0 being black and 1 white). Figure 1D shows examples of such images.

Pre-processing
Depending on the configuration of the microscope used to acquire the images, the input images contain 2 or 3 regions: In the center of the image, a gray region corresponding to the histological section of the cortex. This is the foreground region.
A circular white region corresponding to the lens/plate of the microscope and that surrounds the cortex region. It is considered as background.
In some cases, a black region can also surround the two other regions. It is also considered as background. We extract the background using a K-means algorithm with K = 3 regions. The background is then replaced by the value 0 and the resulting images (the so-called "slices" in the following) are denoted {S 1 , . . . , S Q }.

Vessel Detection
Blood vessels which are approximately orthogonal to the cutting plane look like slightly elliptic spots that can be approxi-mated by small Gaussians. To be invariant to local contrast fluctuations, we detect these orthogonal vessels using normalized cross correlations with Gaussian templates of varying standard deviations, assumed to be smaller than σ max .
For each slice S m , for m ∈ {1, . . . , Q}, we compute its associated normalized cross correlation NCC(S m ) against the set of Gaussian templates, as detailed in Appendix A. Given a threshold τ ∈ [0, 1] and a maximum number N of detected vessels, we define the detected vessel centers X m to be the set of pixels x satisfying both NCC(S m )(x) > τ and NCC(S m )(x) is among the N largest values of NCC(S m ).

Slice Registration by Robust ICP
For each m, we now register the slice S m with the slice S m+1 (see Figure 3A). Registration is obtained by computing an optimal transformation T m which maps pixels in slice S m+1 to pixels in image S m . We restrict here the computation to rigid transformations, i.e. of the form T (x) = R(x) + t where R is a planar rotation and t ∈ R 2 is a translation vector.
Variational registration. This optimal deformation is obtained by exploiting the fact that detected orthogonal vessels should have approximately the same position in two consecutive frames. We denote X m+1 = {x i } i∈I and X m = {y j } j∈J the two sets of vessel positions. We cast this problem as the optimization of a non-convex functional measure of the goodness of fit between the transformation T (x i ) of each detected vessel x i in S m and its closest neighbor y j in S m+1 . T m is obtained by computing a local minimizer of Here ρ : R + → R + is a penalty function. The most common choice is a quadratic loss ρ(r) = r 2 , which assumes some sort of Gaussian distribution of the fitting errors. This choice poorly handles outliers in the detected vessels, which are likely to be present in our datasets. We choose here to use the following robust loss function which gives less weight to outliers (large values of r) than a quadratic loss. Small values of ε are used to cope with many outliers. Note that setting ε → +∞ recovers the quadratic loss r 2 which assumes no outliers. Note also that other loss functions could be used as well, as long as they satisfy the hypotheses exposed in Appendix B.
ICP iterations. A classical algorithm to minimize (1) is the Iterative Closest Point (ICP), introduced by [11] for the quadratic loss ρ(r) = r 2 . This algorithm has been extended by several authors to cope with robust loss (see Section 1.1 for more details). We use a similar approach here, and provide more details in Appendix B. The ICP algorithm iterates between two steps. In the first step, T is known and assumed to be fixed, and one computes a A. Snapshots extracted from the GUI illustrating the results of the automated detection of blood vessel cross sections on consecutive histological slices (yellow filled circles and red open circles, respectively -S 3 to S 6 previously shown in Figure 1D). The result of the registration obtained by robust ICP applied on the detected blood vessels is illustrated on the right for each pair of slices. Note that the user can delimit manually a region of interest (black dashed lines) in order to restrict blood vessel detection to the barrel field region. B. Images from the same experiment obtained by using our automated barrel map reconstruction tool before (left) or after (right) post-processing steps.
nearest neighbor z i = y j for each vessel x i , where the index j minimizes min In the second step, the optimal T is updated by solving min T i∈I For the quadratic loss ρ(r) = r 2 , this second step is solved in closed form as detailed in Appendix B.2. For a generic loss ρ, there is no such closed form. We detail in Appendix B.1 a Majorize-Minimize (MM) method to compute a local minimizer. To the best of our knowledge, this presentation, and the corresponding convergence analysis, is new.

Initialization.
A major difficulty to solve (2) is that it is highly non-convex, and the ICP algorithm is likely to converge to a local minimizer T . To improve the quality of the result, it is important to test several initializations to obtain a good registration.
Registered images. Once the registration transforms {T 1 , . . . , T Q−1 } have been computed, they can be cascaded to warp the input slice images to obtain the sequence {S 1 , . . . ,S Q } of sections, all registered with respect to the initial one S 1 =S 1 as follow Step 2: reconstruction of the barrel image We now have a set {S 1 , . . . ,S Q } of registered slices. We fuse them in a single image S which gathers the edge information of the relevant images to reconstruct the barrel map.

Pre-processing
In order to avoid the amplification of artifacts during the gradient fusion process detailed next, we inpaint (i.e. remove) the orthogonal vessel traces and denoise the resulting image. The inpainting method is detailed in Appendix C. The output of the inpainting is then denoised using a median filter on 3 × 3 patches. We use this non-linear filter to reduce salt-and-pepper noise instead of a convolution, as it removes noise while preserving edges. We denote the output of this pre-processingS m .

Slice gradient fusion
We denote m ∈ M the set of relevant slices containing partial barrel information. To reconstruct a sharp image, we fuse together the gradient of the input slices {S m } m∈M by keeping at each pixel the gradient with the largest magnitude. This method is partly inspired by some recent works in computational photography, such as [17,18]. The details of this method are given in Appendix D. We denote S the output of the fusion process.

Drift removal
The histological sections often present variations in intensity across the barrel cortex that might be due to anatomical reasons. Usually, the anterior lateral barrel subfield (small barrels corresponding to small vibrissae) appears darker than the posterior medial barrel subfield (large barrels corresponding to large vibrissae). This drift is enhanced by the gradient fusion operation that is applied on each pair. As a consequence, the merged image S obtained by the procedure explained above exhibits a strong drift in intensity. We thus filter the merged image with a high-pass filter to remove this low frequency component where h σ is a low-frequency gaussian filter of standard deviation σ , and is the discrete 2-D convolution.

In vivo VSD imaging and DiI staining
Animal preparation and VSDI setup. Experiments were performed in conformity with the French (authorization number: 2012-0068) and European (2010/63/UE) legislations relative to the protection of animals used for experimental and other scientific purposes. VSDI of the cortical activity evoked by single whisker deflections was performed on 6-12 week-old C57Bl6 mice under urethane anesthesia (1.7 mg/g), essentially as previously described in [19]. Briefly, the left barrel cortex was exposed and stained for 1 hour with the VSD RH1691 (1mg/ml, in Ringer's solution containing [in mM]: 135 NaCl, 5 KCl, 5 HEPES, 1.8 CaCl2, 1 MgCl2). After removal of the unbound dye, the cortex was covered with agarose (0.5-1% in Ringer's) and a coverslip. Cortical imaging was performed through a tandem-lens fluorescence microscope (Sci-Media), equipped with one Leica PlanApo 5x (objective side) and one Leica PlanApo 1x (condensing side), a 630 nm excitation filter, a 650 nm dichroic mirror, and a long pass 665 nm emission filter. The field of view was 2.5 × 2.5 mm, resulting in a pixel resolution of 25 × 25 µm.
Whisker stimulation. Individual deflections of the right 24 posterior macrovibrissae of the mice were performed using a custom built multi-whisker stimulator based on a matrix of 24 multidirectional piezoelectric benders [20]. The whiskers were inserted in 27G stainless steel tubes attached to the benders, leaving 2 mm between the tip of the tube and the whisker base. The 24 whiskers were stimulated individually, in the 4 cardinal directions, at 0.1 Hz within pseudo randomized sequences containing extra blank trials (each stimulation being repeated 10 times). Each whisker deflection consisted of a 100 µm displacement (measured at the tip of the tube), with a 2 ms rising time, a 2 ms plateau and a 2 ms fall (specific filters were used to correct for the mechanical ringing of the stimulators).
Image analysis. Acquisition and data preprocessing were done using in-house software (Elphy, G. Sadoc, UNIC-CNRS), further analyses were made using custom written routines in Ig-orPro (Wavemetrics). Subtraction of the averaged unstimulated blank trials was used to correct for bleaching artifacts. For each whisker, data corresponding to the 4 directions of deflection were averaged.
DiI in vivo staining. DiI stain (Molecular Probes, LifeTechnologies) was deposited on the shanks of silicon electrodes [21] that were inserted in the barrel cortex perpendicularly with a microcontroller (Luigs & Neumann).

Histological Procedures
Following the experiments and the administration of an overdose of urethane, mice were perfused with saline followed by paraformaldehyde (4% in 0.1 M phosphate buffer). After an overnight post-fixation in paraformaldehyde, the brains were cut in 100 µm thick tangential sections and stained for cytochrome oxidase.

Results
Starting from tangential slices stained for cytochrome oxidase, the reconstructed 2-D barrel map can be obtained in a few clicks by using the provided GUI. The whole procedure takes 7-9 minutes including 3-5 minutes of computation that do not require the intervention of the user. The same process using traditional manual methods (with the help of a raster graphics editor such as Adobe Photoshop) takes, for a well-trained person, 16 to 28 minutes. Most importantly, our automated approach prevents user dependent variability in the obtained barrel map. In order to assess the accuracy of our registration method we first carried out histological control experiments, and then used VSD imaging to functionally evaluate the precision of the barrel map obtained following the full reconstruction procedure.

Histological Validation of the Registration Method
In order to assess the efficiency of our registration method based on automated blood vessel detection and robust ICP, we perpendicularly inserted DiI coated electrodes in the barrel cortex of urethane anesthetized mice (n = 5 experiments, 3 to 6 electrode penetrations per experiments), before processing the brain for histological staining of the cytochrome oxidase following our standard procedures. The electrodes being flat, they did not leave any round or elliptic white marks in the tissue and therefore did not interfere with our registration method (Figure 4). Because the penetration of DiI coated electrodes was perpendicular to the cortical surface, the DiI staining should appear aligned on consecutive cortical sections following proper registration of the slices. In order to control this alignment, the location of DiI spots was reported for each slice on the final fused barrel map image for each section ( Figure 4C). The calculated mean distance between DiI spots from the same electrode penetration (34.41 ± 18.93 µm, mean ± SD, n = 5) revealed the subcolumnar resolution of our registration method. These control experiments were further used to evaluate the eventual shrinkage of the cortical tissue due to brain fixation and histological procedures. The distances of 250 µm and 500 µm separating the electrode penetration sites in vivo, were compared to the distances measured between DiI spots on the slices following histological procedures. Over our 5 control experiments, the observed tissue shrinkage within the x-y plane was minimal (< 1.5%).

Assessment of the barrel map reconstruction tool using
VSDI of cortical responses to individual whisker deflections To finally validate our 2-D barrel map reconstruction method, we confronted its resulting map with the functional organization of the barrel cortex established in vivo by real time imaging of cortical responses to individual whisker deflections under urethane anesthesia (n = 4 experiments). Using a mechanical multi-whisker stimulator [20], we deflected independently the 24 principal whiskers in a pseudo random order, and imaged the evoked cortical responses in the contralateral barrel cortex using the VSD RH1691. Figure 5 illustrates the results obtained from one experiment. As previously reported in similar conditions [19], the earliest responses to whisker deflections were localized to the corresponding barrel-related columns ( Figure 5A). When reporting the 90% contours of the early cortical responses onto the aligned barrel map, we observe a good anatomo-functional match ( Figure 5B). To quantify this match, the distance between the centroid of the anatomically defined barrels and the centroid of the early VSD response (area above a 90% threshold) was measured ( Figure 5C). The mean centroid-centroid distance over the 4 control experiments for all the barrels is 60.5 ± 21.2 µm (34.8-116.5 µm range across the barrel field). We observed slightly higher values for the columns located at the border of the map which might result from the curvature of the cortex, the maxima of cortical responses were located within the corresponding barrel area in the majority of cases (86.55%), attesting to the accuracy of the method.

Discussion
We have designed a comprehensive pipeline for the reconstruction of the 2-D barrel map from histological sections. This tool enables a fast reconstruction of a precise barrel map, thus saving a significant amount of time for the experimentalist. Indeed, most of the studies based on optical imaging which required a post-hoc anatomo-functional mapping of the barrel cortex relied on manual reconstruction and alignment of the barrel map with the functional images [22,23,19,24,25,26,27]. Other studies [28,29,30] used a method based on a warping algorithm described by [31]. However, the manual detection of the fiducial markers required with this approach remains time-consuming in comparison to the automatic detection of blood vessel cross sections proposed here.
In a recent article, Guy and collaborators [32] used a warping algorithm in order to align the layer IV barrel map, reconstructed manually from histological slices of flattened cortex, with the in vivo functional images using sets of fiducial points. Although the algorithm is not described in detail and probably involves manual selection of the fiducial points, it might be an interesting complement to our approach when working with flattened barrel cortex slices, since it allows a compensation for the curvature of the brain and distortion of the tissue. Finally, instead of using the superficial blood vessels as anatomical landmarks to align the barrel map on the functional images, an alternative approach is to use the images of early cortical responses to single whisker deflections as landmarks [33,34,35]. Using such a method as a standard requires the acquisition of several additional single whisker responses, which might be difficult to implement for instance when working with awake head fixed animals. Note that the tool we propose here to reconstruct the barrel map is valuable whatever solution is chosen in the end to realign the barrel map with the functional images.
On the methodological and mathematical sides, we mainly re-use a set of already existing tools (cross-correlation, ICP and gradient fusion). Our main contribution is to put them together in a coherent processing pipeline. A result of independent interest, that seems to be new, is to show how a family of robust ICP algorithms can be recasted as majorization-minimization algorithms. This in turn allows us to analyze the convergence of these methods.
Obviously the efficacy of the proposed anatomo-functional mapping tool depends upon the quality of the histological slices. Although the preparation of these slices relies on standard protocols which often belong to the daily routines of neurophysiology laboratories, two aspects are essential for the accuracy of the outcome: the quality of the perfusion, and the right thickness of the first (most superficial) slice. Indeed, on the one hand blood vessels have to appear as white circular or elliptical spots on the images to allow the proper registration of consecutive slices and, on the other hand, the superficial blood vessels should be intact to allow the final overlay of the obtained barrel map with the in vivo recordings. When cutting the brain, it is therefore important to set the zero position of the blade with care to ensure a 100 µm thickness to the first slice and thus preserve the integrity of the superficial blood vessels. Finally, although one could deplore that this overlay is a crucial step of the analysis that remains to be achieved manually, we propose here a solution that automatises the most time-consuming phases of the overall process and thus represents a substantial gain in time and precision. Although its use has been demonstrated successfully for VSDI of the adult mouse barrel cortex, it could be expanded to other experimental model species or to the developing brain. Furthermore, the histological section registration method described here might be helpful to reconstruct any tissue in which a majority of blood vessels are orthogonal to the cutting plane of the slices. We give here the details of an iterative algorithm to compute a local minimizer (in fact a stationary point) of (3), which reads This method is similar to re-weighting 2 methods often used for robust ICP (see for instance [36]), but we integrate it into a Majorize-Minimize framework, which ensures its convergence. Starting from an initial transform R (0) , we compute the iterations as whereẼ is a so-called surrogate function, which should satisfy Under these conditions, it can be shown that the iterations enjoy some good convergence properties. The sequence E(T ( ) ) is decaying and converges to some value E . If E is smooth (which is the case here), ||∇E(T ( ) )|| → 0. Since in our case, the energy E is coercive, the sequence T ( ) is bounded, and all its cluster points T are stationary (i.e. ∇E(T ) = 0) with same energy The main difficulty in general is to devise a "good" surrogate functionẼ, i.e. such that one can compute the iteration (B.2) in closed form. The following proposition shows that one can actually design such a surrogate function using a quadratic loss. Proposition 1. If ρ is C 1 (R) and w(r) = ρ (r) 2r is decreasing, there exists a constant C(T ) independent of T so that the func-tionalẼ is a majorizing functional for (B.1) and thus satisfies properties Proof. We rewriteẼ as where we defined Thanks to the separabilityẼ (it is a summation over i of functions involving independent variables) and the change of variables r = ||T (x i ) − z i ||, it thus suffices to prove thatρ is a surrogate functional for ρ on R + . Hypothesis (H 1 ) holds because ρ is C 1 , and one verifies thatρ(r , r ) = ρ(r ) so that (H 3 ) holds. For any r 0, we consider It satisfies h(r ) = h (r ) = 0 and h (r) = 2r(w(r ) − w(r)). Since w is decaying, r is the only point where h is vanishing on R + . This implies that h 0, hence (H 2 ).
The hypothesis that w is decreasing should be interpreted as the condition that ρ should penalize less than a quadratic loss, which makes sense for a robust penalization. Note that the loss (2) that we use in our method satisfies this condition, and that the weighting function satisfies w(r) = 1 ε 2 + r 2 .

Appendix B.2. ICP Step 2 with Weighted Quadratic Loss
We consider the problem of solving min T i∈I where R is a rotation and t ∈ R 2 . This minimization appears in the MM iteration (B.2) when using the majorizing function (B.3). This problem has a closed form solution, as detailed for instance in [37]. For the sake of completeness, we recall the steps of the method. One first centers the points, for i ∈ Ĩ x i = x i − k∈I w k x k k∈I w k andz i = z i − k∈I w k z k k∈I w k .
The optimal rotation is obtained as R = VU T where (U, V) are the eigenvectors of the correlation matrix i∈I w ixiỹi T = UΛV T (B.4) (here Λ is the diagonal matrix of eigenvalues). The optimal translation is then computed as t = i∈I w i (z i − Rx i ) i∈I w i .

Appendix C. Inpainting
We consider a registered sliceS m and its associated vessel locations X m = {x i } i∈I . We denotex i = T m (x i ) the registered vessel locations, where the cumulative transform T m is defined in (5). We recall that the cross-correlation minimization (as detailed in Appendix A) outputs at each pixel x the index k(x) of the optimal Gaussian template at this location, which has a radius σ k(x) .
We define a mask Φ, which is the set of pixels that are at distance smaller than σ k(x i ) from the pointx i . It is thus a union of disks. Pixels in Φ should be discarded and inpainted. This is achieved using a quadratic minimization that seeks a smooth interpolation of missing data where ∇S is a finite difference approximation of the gradient of the image S . The solution of (C.1) corresponds to solving a Poisson equation ∆S = 0 on Φ with Dirichlet boundary conditions given by the constraints. This can be solved using a conjugate gradient method.

Appendix D. Gradient Domain Image Fusion
We consider a set {S m } m∈M of input images to fuse. At each pixel x, we denote the index of largest gradient magnitude as where ∇ is a finite differences approximation of the gradient operator.
We design a fused vector field as Since the vector field u is obtained by gluing together gradients from several different images, it is in general not anymore the gradient of an image. We thus reconstruct a valid fused image S using the minimal norm pseudo-inverse, i.e. by computing an image S whose gradient is as close as possible to u with adequate boundary conditions, where ∆ = div •∇ is the Laplacian operator and div = −∇ * is the divergence. When using periodic boundary conditions (which can be used in our case), one solves (D.1) in O(N log(N)) operations (where N is the number of pixels) using an FFT Poisson solver.