Information-Theoretic Registration with Explicit Reorientation of Diffusion-Weighted Images
IIEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. X, NO. X, MAY 2018 1
Information-Theoretic Registration with ExplicitReorientation of Diffusion-Weighted Images
Henrik G. Jensen ∗ , Franc¸ois Lauze ∗ , and Sune Darkner ∗ Abstract —We present an information-theoretic approach to registration of DWI with explicit optimization over the orientational scale,with an additional focus on normalized mutual information as a robust information-theoretic similarity measure for DWI. The frameworkis an extension of the LOR-DWI density-based hierarchical scale-space model, that varies and optimizes over the integration, spatial,directional, and intensity scales. We extend the model to non-rigid deformations and show that the formulation provides intrinsicregularization through the orientational information. Our experiments illustrate that the proposed model deforms ODFs correctly and iscapable of handling the classic complex challenges in DWI-registrations, such as the registration of fiber-crossings along with kissing,fanning and interleaving fibers. Our results clearly illustrate a novel promising regularizing effect, that comes from the nonlinearorientation-based cost function. We illustrate the properties of the different image scales, and show that including orientationalinformation in our model make the model better at retrieving deformations compared to standard scalar-based registration.
Index Terms —Registration, Diffusion-Weighted MRI, Normalized Mutual Information, Explicit Reorientation, Analytical Gradients,Density-Based Formulation, Watson Density Function, Free-Form Deformation, Cubic B-Spline. (cid:70)
NTRODUCTION D IFFUSION
Weighted Images (DWI) is a non-invasiveMRI protocol that can be used to infer cytoarchitectureby tracking the movement of water molecules, otherwiseinvisible in structural MRI. However, the geometry of DWImakes it a challenge in image registration, a key tool forcomparing and analyzing medical image data. In addi-tion, DWI acquired on different scanners or with differentprotocols have a non-linear relationship and defining thesimilarity between two DWI is inherently difficult. [1].Our contribution is a scale-space formulation of densityestimation for image similarity for DWI in a nonrigid reg-istration setting. The model encompasses explicit reorienta-tion providing the computational framework for estimatingnon-linear similarity measures well-suited for DWI. Themodel allows for image similarities such as normalized mu-tual information (NMI) by extending the Locally OrderlessRegistration for Diffusion-Weighted Images (LOR-DWI) [2]to nonrigid registration. LOR-DWI includes orientationalinformation in the image similarity, which allows explicitoptimization over the reorientation of diffusion gradients.This extends from DTI to the raw high-angular resolutionscans (HARDI) or the topographically inverted OrientationDensity Functions (ODF).We demonstrate the properties of the framework, suchas the effect of the reorientation on the ODFs, the effectof the scale-space on optimization and regularization onsimulated DWI data and on a synthetic deformation of datafrom a subject of the Human Connectome Project (HCP) [3].We show that LOR-DWI formulation preserves the ODFsand produce excellent mappings for crossing, kissing and
Corresponding author: [email protected] ∗ Department of Computer Science, University of Copenhagen, DenmarkManuscript, 2018. curving fibers, while providing an inherent regularizingeffect.The density formulation also allows us to optimize overthe isocurves of the DWI signal I : R × P → R . Inprinciple, this framework could be extended to other non-trivial geometry of other image modalities where scale-space structures can be defined. In this work, we focus onDWI data where we use a classical scale-space structure onthis non-linear geometry. ELATED W ORK
This work addresses two major challenges in voxel-basedregistration of DWI: The reorientation of DWI in imageregistration, and the non-linear similarity between DWI.Image registration refers to a process that transformsdata into a shared coordinate system. For DWI, the commonway to register two images is to use scalar-based methodson quantitative measures, such as the fractional anisotropy(FA) or the mean diffusivity (MD) [4]. However, as suchmethods disregard most of the directional information inDWI, methods have been developed to also account forthe reorientation of the diffusion profile. Most of theseare created on top of scalar-based methods and iterativelyreorients the gradients after changing the deformation field.Some of the most popular methods can be found in [5], [6],[7], [8]. However, registration frameworks have also beendesigned with an objective function that explicitly optimizesover the reorientation of the gradients, such as
DTI-TK [9],
DT-REFineD [10], and the more recent
DR-TAMAS [11].These frameworks have been shown to generally outper-form scalar-based frameworks for DWI [12], [13], [14], [15].For scalar based approaches to image similarity, registra-tion is straightforward as any popular scalar-based measurecan be used, e.g. the sum-of-squares difference (SSD) [10]or mutual information (MI) [16], [15]. Explicit reorientation a r X i v : . [ c s . C V ] J un EEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. X, NO. X, MAY 2018 2 strategies inherently define the similarity over both positionand orientation and, once the full diffusion profile is partof the similarity, such measures can be defined in a well-suited way for the non-linear relationship between DWI.Both [9] and [10] used variations on SSD in the objectivefunction, while MI was used in [2] through an extensionof the LOR framework in [17]. As argued in [2], the in-variant and statistical properties of MI makes it a logicalchoice for DWI, where multiple factors result in a morestatistical relationship, such as variations in b -values, non-monoexponential behavior in biological tissue, and inter-scanner variability [1]. MI is often used in the standardregistration of complex modalities and, as such, in scalar-based registration of DWI [18] and in the pre-processing ofDWI [19]. MI and normalized MI (NMI) [20] provides a non-linear statistical measure for DWI. However, it is also likelythat more functional measures could be well-suited, suchas cross-correlation (CC) and normalized CC (NCC), bothof which can be defined from the LOR density-formulation[21], [22]. The density-based DWI comes from the general-ized way of estimating image similarity measures based onLocally Orderless Images (LOI) [23]. The first mention ofLOI in the context of image registration was in [24] wherea variational approach to image registration was presented.The LOR framework [25], [22] generalized a range of simi-larity measures as linear and non-linear functions of densityestimates for scalar-valued images. OCALLY O RDERLESS
DWI
We start by introducing a few necessary notations. Ω ⊂ R is the spatial domain of the images under consideration. Ascalar image is a function I : Ω → R . We assume that wecan extend it on the whole R , for instance by extending itwith 0. The projective space of directions of R is denotedby P , and the unit sphere of R by S . We will encounter spatio-directional images I : Ω × P → R , which we similarlyassume too be extendable to R × P . This is necessary inboth cases in order to define their spatial smoothing viaconvolution. We use the following elementary property: As P can naturally be identified as the quotient S / {± } bythe antipodal symmetry, a function f : P → R can belifted to an antipodal symmetric function ˜ f : S → R .Conversely, any antipodal symmetric function g : S → R factors through P . A spatio-directional image can (and will)be lifted to an antipodal symmetric image ˜ I : Ω × S : ˜ I ( x , − v ) = ˜ I ( x , v ) . We will denote by I both the spatio-directional image and its antipodal symmetric lifting in thefollowing. The LOR framework defines the density estimates overthree scales: The image scale σ , the intensity scale β , andthe integration scale α . In the context of scalar registration, for a transformation φ : R → R , the estimated histogram h and the corresponding density p is computed as h βασ ( i, j | φ, x ) = (1) (cid:90) R P β ( I σ ( φ ( x )) − i ) P β ( J σ ( x ) − j ) W α ( τ − x ) d τ p βασ ( i, j | φ, x ) (cid:39) h βασ ( i, j | φ, x ) (cid:82) Λ h βασ ( k, l | φ, x ) dk dl (2)where i, j ∈ [ a , a ] =: Λ are values in the image intensityrange, I σ ( φ ( x )) = ( I ∗ K σ )( φ ( x )) and J σ ( x ) = ( J ∗ K σ )( x ) are images convolved with the kernel K σ with standarddeviation σ , P β is a Parzen-window of scale β , and W α isa Gaussian integration window of scale α . The marginalsare trivially obtained by integration over the appropriatevariable. The LOR-approach to similarity lets us use a set ofgeneralized similarity measures, the linear and non-linear F lin = (cid:90) Λ f ( i, j ) p ( i, j ) di dj (3) F non − lin = (cid:90) Λ f ( p ( i, j )) di dj (4)where the linear measure f ( i, j ) includes e.g. sum ofsquared differences and Huber, and the non-linear f ( p ( i, j )) includes e.g. MI, normalized MI (NMI), see [22] for details. This work addresses the estimation of the image similarity F of DWI in the context of nonrigid registration an exten-sion of our previous work [2]. In this context, I and J arespatio-directional signals. Specifically, DWI MR attenuationsignals at location x , for a gradient direction v , are modeledby S ( x , v ) = S ( x ) e − b I ( x , v ) [26] and apparent diffusioncoefficients volumes are given by I ( x , v ) = − b log S ( x , v ) S ( x ) .Gradient directions v belong to S but the diffusion isorientation-free such that I ( x , v ) ≈ I ( x , − v ) which nat-urally defines a spatio-directional image Ω × P → R . Inorder to apply LOR-DWI, the histogram and density esti-mates Equation (1) and eq. (2) must be extended to spatio-directional data, and the action of the spatial transformationon the directions must be defined.We introduce a kernel on the sphere as an extension tothe density estimates of LOR to include directional infor-mation. This kernel accounts for directional smoothing anddefines our LOR-DWI framework. Thus, we extend spatialsmoothing to be spatio-directional, such that the directionalsmoothing preserves this symmetry and the projective struc-ture, via a symmetric kernel Γ κ ( ν , v ) on S . We define thesmoothed signal I σ,κ at scales ( σ, κ ) by I σκ ( x , v ) = (cid:90) S (cid:18)(cid:90) R I ( τ , ν ) K σ ( τ − x ) d τ (cid:19) Γ κ ( ν , v ) d ν = ( I ∗ ( K σ ⊗ Γ κ ))( x , v ) (5)where K σ ( x ) is a Gaussian kernel with σ standard devi-ation. We employ a symmetric Watson distribution [27] as Γ κ ( ν , v ) for directional smoothing on S , given by Γ κ ( ν , v ) = Ce κ ( ν (cid:62) v ) (6) C = M ( 12 , d , κ ) (7) EEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. X, NO. X, MAY 2018 3 where M ( , d , κ ) the confluent hypergeometric functionalso called the Kummer function ( d = 3 in our case) [28], ± v the center of the distribution, and κ the concentrationparameter, which is roughly inverse proportional to thevariance on the sphere. Because of the symmetry propertyof the Watson distribution and the antipodal symmetry of I ( x , v ) , it is clear that I σκ ( x , v ) is antipodal symmetrictoo. As one alternative, a symmetrized von Mises-Fisher[27] distribution or a symmetrized heat kernel could beconsidered. A transformation φ : Ω → Ω acts on directions via itsJacobian: at x ∈ Ω , J x φ sends R v to R J x φ v which iswell-defined as soon as det( J x φ ) (cid:54) = 0 which is the case ifwe assume that φ is diffeomorphic. Via the representation P (cid:39) S / {± } , we can write P x φ : {± v } (cid:55)→ {± J x v | J x v | } .We denote by ψ the mapping ( x , v ) (cid:55)→ J x v | J x v | . Becauseit will only appear inside an antipodal symmetric kernel, ψ x = ψ ( x , · ) represents P x φ without ambiguity . Because ofspatial convolutions, we assume that φ can be extented to R , assuming that it is the identity out of a compact set D containing Ω . In the following, we will often omit the parameter x . We set Φ = ( φ, ψ ) , and we write the joint histogram and densityfor similarity in image registration as h βασκ ( i, j | Φ , x ) = (8) (cid:90) R × S P β ( I σκ ( φ ( x ) , ψ ( v )) − i ) P β ( J σκ ( x , v ) − j ) W α ( τ − x ) d τ × d v p βασκ ( i, j | Φ , x ) = h βασκ ( i, j | Φ , x ) (cid:82) Λ h βασκ ( i, j | Φ , x ) dl dk . (9)Assume that W is a Gaussian kernel with standard devi-ation α . If we left α → ∞ , we would obtain full spatialintegration and can ignore the integration scale α . In [2], we assumed that φ was a global affine transformation,with ψ x = ψ being the projectivization of its linear part. Inthe present work, we assume φ to be a nonrigid transfor-mation, and we use Rueckert et al. framework [29] wherethe transformation φ is given as a hierarchical spline repre-sentation linearly parameterized by a spatial grid of controlpoints c . In the next subsection, we use its vectorized form,so that φ = Bc , Φ c = ( φ c , ψ c ) . We seek the transformation Φ ∗ c which maximizes the regularized NMI Φ ∗ c = arg max c M (Φ c ) = arg max c F ( I ◦ Φ c , J ) + S (Φ c ) (10)The dependency of M in c is complex. The followingdependency diagram (Fig. 1) for the mutual informationterm illustrates it. N M I ( I σκ ◦ Φ c , J σκ ; β ) H I σκ ◦ Φ c ; β H J σκ ; β H I σκ ◦ Φ c , J σκ ; β p β, I σκ ◦ Φ c p β, J σκ p β, I σκ ◦ Φ c , J σκ h β, I σκ ◦ Φ c , J σκ I σκ ( φ c , ψ c ) ψ c J σκ c φ c Fig. 1. Dependency graph of the nonrigid DWI registration betweenthe moving image I and the target image J , with normalized mutualinformation (NMI) as the similarity measure. The deformation is pa-rameterized by c so that any change in c will eventually affect thetotal similarity between the two images. The red frame contains MIregistration elements which are classical, while the blue frame containsthe part which extends the LOR framework to the LOR-DWI. We use quasi-Newton methods to compute an optimum ofEquation (10), in particular L-BFGS[] from [30], and thusneed to compute the gradient of the similarity with respectto control points vector c . Most of the required calculationsare provided in [29], [22], [2], they correspond to the redframe in Figure 1, a full formula for the gradient is longand not informative for the present work. We thereforeonly describe how spatio-directional smoothing contributesto the image similarity, i.e., the part of Figure 1 enclosedin the blue frame. Thanks to the LOR-DWI representa-tion, the present extension only appears within the spatio-directional smoothing of I . One complication originatesfrom choices in our implementation, where the Kummerfunction M ( , , κ ) as the normalization constant in theWatson kernel is replaced by an estimate of this factor froma discrete set of N directions ν , . . . , ν N at each voxel x . Itcan therefore no longer be considered constant. We rewritethe discrete spatio-directional smoothing as I σκ ( φ c ( x ) , ψ c ( v )) = N (cid:88) n =1 (cid:90) R I ( τ, ν n ) K σ ( φ c ( x ) − τ )¯Γ κ ( ν n , ψ c ( v )) dτ (11)where we have set ¯Γ κ ( ν n , ψ c ( v ) = e κ ( ν (cid:62) n ψ c ( v )) (cid:80) Ni =1 e κ ( ν (cid:62) i ψ c ( v )) . (12) EEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. X, NO. X, MAY 2018 4
We compute the Jacobian of the spatio-directional smooth-ing with respect to the control point parameter c : J c I σκ ( φ c ( x ) , ψ c ( v )) = N (cid:88) n =1 (cid:18)(cid:90) R I ( τ, ν n ) J c K σ ( φ c ( x ) − τ ) dτ (cid:19) ¯Γ κ ( ν n , ψ c ( v )+ N (cid:88) n =1 (cid:18)(cid:90) R I ( τ, ν n ) K σ ( φ c ( x ) − τ ) dτ (cid:19) J c ¯Γ κ ( ν n , ψ c ( v ) (13)From eq. (13), we have J c K σ ( φ c ( x ) − τ ) = − K σ ( φ c ( x ) − τ ) σ ( φ c ( x ) − τ ) (cid:62) B (14)because K σ is a Gaussian kernel with variance σ and φ c = Bc , its Jacobian J c φ c is simply B . The part involvingdirectional kernel is somewhat more complex, as it involvesthe Jacobian J x φ c and the normalizing factor in Equa-tion (12). Set f ( x ) = e κx and define f i = f ( ν i (cid:62) ψ c ( v )) .Since f (cid:48) ( x ) = 2 κxf ( x ) , we have J c f ( ν (cid:62) ψ c v ) = 2 κf ( ν (cid:62) ψ c ( v )) ψ c ( v ) (cid:62) νν (cid:62) J c ψ c ( v ) . (15)By a straightforward and careful calculation, we obtain J c ¯Γ κ ( ν n , ψ c ( v )) = (16) − κ ¯Γ κ ( ν n , ψ c ( v )) (cid:80) Ni =1 f i ψ c ( v ) (cid:62) (cid:88) i (cid:54) = n f i ν i ν (cid:62) i J c ψ c ( v ) (17)To understand J c ψ c , we need to deal with tensors . Thematrix B is built of cubic B-splines. With k control points, c has dimension k and B ( x ) ∈ R × k . Its spatial Jacobianis a 3D- tensor made of quadratic B-splines, J x B ∈ R × k × . J x φ c = J x B · c where the ’ · ’ operator represents the contrac-tion R × k × × R k → R × , ( r uvw , s v ) (cid:55)→ ( (cid:80) v r uvw s v ) uw .One has J c J x φ c ∈ R × × k (cid:54) = J x J c ( Bc ) = J x B ∈ R × k × , though the difference is a matter of swappingindices. J c ψ c is a 3D tensor of the same order as J c J x φ c .Another contraction comes from J c J x φ c ( v ) . This is a matrixin R × k and one can write J c J x φ c ( v ) = J x B • v where • is the contraction R × k × × R → R × k , ( r uvw , t w ) (cid:55)→ ( (cid:80) w r uvw t w ) uv . The differentiation of the inner product (cid:104) J x φ c v , J x φ c v (cid:105) is given by J c (cid:104) J x φ c v , J x φ c v (cid:105) = ( J x φ c v ) (cid:62) ( J x B • v ) (18)Denoting by V the vector J x φ c v , so that ψ c ( v ) = V | V | , theJacobian J c ψ c v is J c ψ c v = (cid:32) I − V V (cid:62) | V | (cid:33) J x B • v | V | (19)where I is the identity of R . Note also that V = J x φ c v =( J x B · c ) v and that I − V V (cid:62) | V | is the orthogonal projection π V ⊥ onto V ⊥ , the plane orthogonal to V . Putting things together, the Jacobian with respect to the control pointparameter c of the spatio-directional smoothing is given by J c I σκ ( φ c ( x ) , ψ c ( v )) = − N (cid:88) n =1 (cid:90) R I ( τ, ν n ) K σ ( φ c ( x ) − τ )¯Γ κ ( ν n , ψ c ( v )) × ( φ c ( x ) − τ ) (cid:62) B σ +2 κ V (cid:62) (cid:16)(cid:80) i (cid:54) = n ν i ν (cid:62) i (cid:17) | V | (cid:80) Ni =1 f i π V ⊥ (cid:18) J x B • v | V | (cid:19) dτ. (20) So far, we have not specified the form of the regularizer S (Φ c ) in Equation (10). The regularization has receivedlittle attention due to the inherent regularization from thesmooth kernels and the additional directional structure.However, we found that the last steps in the hierarchicaltransformation model required some regularization to keepthe deformation stable, due to the high resolution of thedeformation field. We used a simple regularization thatpenalizes a non-uniform grid by the squared differencebetween a point c i and its direct neighbours. The controlpoints are organized as a family of R grids, from coarseto fine resolution, and the regularizer S (Φ c ) is the sum (cid:80) Rr =1 S r (Φ c r ) at each resolution, where c r = ( c r , . . . c rp r ) T is the grid of control points at resolution level r , and wedenote by N r ( i ) the set of indices j such that control point c rj is neighbor to control point c rj and by | N r ( i ) | its cardinal.We set S r (Φ c ) = − λ r p r (cid:88) i =1 (cid:107) c ri − | N r ( i ) | (cid:88) j ∈ N r ( i ) c rj (cid:107) . (21) λ r is a strictly positive parameter controlling the degree ofsmoothness. The negative sign in Equation (21) comes fromthe fact that we maximize the objective function. In order tocompute the gradient of S r (Φ c r ) , we define series of linearmappings T ri : c (cid:55)→ | N r ( i ) | c ri − (cid:80) j ∈ N r ( i ) c rj . The regularizerin Equation (21) can be rewritten as S r (Φ c ) = − λ r p r (cid:88) i =1 | N r ( i ) | (cid:107) T ri c r (cid:107) (22)and by classical manipulation, we obtain that ∇ c S r (Φ c ) = − λ r p r (cid:88) i =1 | N r ( i ) | T r ∗ i T ri c r . (23) T r ∗ i is the adjoint of T ri . Operator − (cid:80) p r i =1 1 | N r ( i ) | T r ∗ i T ri is adiscrete Laplacian. Our sought gradient is ∇ c S (Φ c ) = − R (cid:88) r =1 λ r p r (cid:88) i =1 | N r ( i ) | T r ∗ i T ri c r . (24)In this paper, we chose λ = · · · = λ R . XPERIMENTS
To illustrate the properties of the LOR-DWI method, weconduct a series of experiments on simulated data andartificially warped real data.
EEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. X, NO. X, MAY 2018 5
Artificially generated samples enables visual inspection andvalidating of the framework in a highly controlled en-vironment. While artificial experiments can be found inmost DWI registration papers, the artificial data is rarelyavailable. To our knowledge, there are no open sources ofsimulated DWI data for comparing registration frameworks,although the
DIPY project may be a good source for gener-ating simulations [31]. To this end, we created our own sim-ulated DWI data, which is freely available. The generationof the artificial data is described in the following. The simulated DWI (HARDI) samples were created by de-forming a unit sphere of equally distributed directions [32]to a certain HARDI or ODF shape. Figure 2 shows twosimulated HARDI samples and their corresponding ODFs,where the first is a single fiber ODF, and the second is acrossing fiber ODF. DWI samples are antipodal symmetricand every ODF from the simplest to the most complex canbe constructed through a combination of single fiber ODFs.The simulated data is visualized using the regularized QBIalgorithm, which uses a linear combination of real sphericalharmonics to represent either the direct QBI sample ortransformed ODF [33]. These models can be combined to (a) HARDI sample (b) ODF of (a)(c) HARDI sample (d) ODF of (c)Fig. 2. Simulated DWI samples. The left column shows the raw DWIsignal. The right column shows the reconstructed diffusion ODFs thatfollow anisotropic diffusion. The red lines indicate fiber orientations. form simulated DWI tracts in various DWI shapes, suchas crossing fiber tracts Figure 3. A × grid is usedthroughout this section to create blueprints of fiber tractconstellations with the purpose of performing LOR-DWIregistrations. The images are coloured according to the gen-eralized FA (GFA) value where dark blue regions representfree isotropic diffusion. To simulate a more realistic DWIscenario, random uniform noise has been added to thesamples in Figure 3. The simulated voxels with unit densityare rescaled for the free diffusion to have a low density ormean diffusivity to resemble real b normalized data.
1. Contact us on [email protected] for the code or examples. (a) HARDI fiber tract crossing(b) ODFs of (a) showing the fibersFig. 3. A simulated DWI fiber tract crossing with uniform random noise.(a) HARDI fiber tract crossing (b) ODFs of (a) showing the fibersFig. 4. Simulated DWI fiber tract crossing shown form above. Theisotropic diffusion have now been normalized to have a lower meandiffusivity.
For consistency the same setup is used for all experimentson the simulated data, unless specifically stated otherwise.
Hierarchical mesh resolution.
In the B-spline deformation,the spacing between the control points is decreased inorder to iteratively increase the degrees of freedom inthe registration. We use Φ local = Φ δ =4 +Φ δ =3 . +Φ δ =3 +Φ δ =2 . 10 iterations is used for all resolutions except thelast, which terminates based on the optimal toleranceof (cid:15) = 1e − , or 90 iterations. Watson concentration, directional resolution.
The concen-tration parameter is set to κ = 15 , which is sufficientlysmooth to represent the 100 uniform directions used. Spatial resolution.
A full spatial resolution is used with aB-spline smoothing at a near-Gaussian variance of σ = EEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. X, NO. X, MAY 2018 6 . voxels. Histogram size.
Due to the relatively small data sample,a low number of bins is used for the histogram (i.e.intensity resolution) × . This allows for some larger,stable but less refined deformations [22]. HARDI registration, ODF visualization.
We register theraw HARDI models but we visualize the ODFs of thewarped data, based on the Funk-Radon transformed(FRT). We do this to illustrate that the warped raw datais correctly reoriented and do not suffer visibly fromaffine shearing.
Regularization.
We use regularization by a factor of λ =1e − . The regularization could have be omitted inthe simulated experiments and substituted by multiplelevels of resolution of the deformations (large to small),as the simulated data is very structured. As an aside,every example below could be highly improved by aexperiment-specific choice of parameters, but for con-sistency the same settings were use for all experiments. The first set of experiments are based on artificially gen-erated distributions of HARDI shells and correspondingODFs, each representing different fiber constellations. Theexperiment maps a straight fiber tract to a curved tract ofthe same width. It offers an opportunity to discuss some ofthe differences between a good reconstruction and a correctmapping.
Three experiments illustrating the different scenarios wherestraight fibers are registered to wavy fibers. These exper-iments illustrate the regularizing effects of using the fulldiffusion profile for registration. In the first experiment, theposition of the end of the fibers is constrained by addingintersecting fibers at the end and beginning. Figure 5 showsthe simulated moving (straight) and target (wavy) images.The results of the first experiment is shown in Figure 6, (a) Moving Image (b) Target ImageFig. 5. Experiment 1: Simulated fiber tract images with intersecting fibersat the boundaries. where 6a is the reconstructed warp, and 6b shows thefinal spatial mapping from the moving image, overlaidon the original target image. As the figures illustrate thestraight fibers are stretched in correspondence with thecurvature, and the reconstructed ODFs from the HARDI-based registration are rotated correctly, although smoothedby the interpolation. The second experiment is performed (a) Reconstructed warped image (b) Spatial mappingFig. 6. Experiment 1: Registration of single tract images of varying shapeand (given the boundary) length. without features on the boundaries. The results are shownin Figure 7, where we observe that the length is preservedas the straight fiber tract is mapped to a sub-part of thecurved tract. Intuitively, a correct registration in such acase should preserve the length of the straight fiber afterdeformation. The fact that this happens with no strong (a) Reconstructed warped image (b) Spatial mappingFig. 7. Experiment 1: Registration of without boundary fibers. The spatialmapping indicates a nice preservation of length of the straight tractmapped to a subsection of the curved tract. outside regularization, illustrate an inherent regularizationeffect of the cost function.In the third experiment, the proposed similarity is com-pared with the equivalent scalar-based registration by per-forming a mean diffusivity registration ( κ = 0 ) of the tractswith no signal on the boundaries. The mean diffusivitycarries no directional information. The result can be seenin Figure 8. This pure scalar-based registration is driven bythe edges of the simulated tracts only. The reconstructionappears to be correct, but the final spatial mapping indicatesa lack of regularization as the fibers are stretched unevenly. The second set of experiments are designed to investigatethe registration of crossing fiber tracts.
The first experiment examines the framework’s ability toregister two crossing tracts with a horizontal and verticalshift Figures 9a and 9b. Circular fibers have been added asfixed points in the image to illustrate the local shift of thecrossing tracts. The result of the registration is shown inFigures 9c and 9d. The final spatial mapping, shown with
EEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. X, NO. X, MAY 2018 7 (a) Reconstructed warped image (b) Spatial mappingFig. 8. Experiment 1: Scalar-based registration without boundary fibersand no directional directional information in the cost function duringoptimization. arrows, is accurate including the reconstruction. Note thatthe reconstruction is subject to the effect of the smoothing inthe interpolation.
The second experiment involving crossing fibers show threefiber tracts crossing under a varying amount of shear with afixed horizontal tract ( Figures 10a and 10b). The purpose isto investigate a relatively large deformation combined witha change to the complex crossing at the center. The resultsare shown in Figures 10c and 10d. As the figure illustratesthe structure of the complex center-crossing closely matchesthe orientations of 45 degree crossing fibers. We remindthe reader that the registration was not performed on theODFs but directly on the simulated HARDI models andsubsequently reconstructed.
The third set of experiments examines the most complexconfigurations namely the registration of kissing and fan-ning fibers. Note that although the artificial cases con-structed for these experiments have a 1-1 correspondence.
The first experiment consist of two DWI images, that sim-ulates both fanning (dispersing) and kissing (interleaving)fiber tracts. The moving image in Figures 11a and 11bcontains a crossing with a few fibers fanning in and outalong the vertical center line. The target image containstwo curved tracts fanning in and out, and merging at thecentral crossing. The results of the registration experimentare shown in Figures 11c and 11d, where both the recon-structed warp and spatial mapping show a registration thatfollows the lines of the original target image. The resultingregistration is able to move, shrink and turn the center-crossing to fit to fit the target image.
The second experiment involves two straight fiber tractsthat are registered to two curving and interleaving tracts(Figures 12a and 12b). The result, shown in Figures 12cand 12d, displays a large and difficult deformation that byall accounts appears to be a successful registration.
In this set of experiments, the registration framework isevaluated on real data obtained from the HCP [3]. By in-troducing a random synthetic warp on a subset of the braindata, a ground truth is obtained where the objective is toregister the warped image back to the original image. Theseexperiments illustrate the effect of the parameters of themodel in a realistic scenario in a guaranteed diffeomorphicscenario.
Fig. 13. Selected region of interest for synthetic warp experiments (blue)overlaid on the b image of HCP subject . The DWI data used in this experiment is shown in Figure 13,where the region of interest (ROI) is the blue overlay on thesubjects b image. An ROI of the brain is used to improvethe visualization of the results. The ROI was chosen at theedge of the corpus callosum (CC) in the left hemisphere dueto the characteristic C-shape of the CC and the intersectionwith other well-known structures e.g. the cingulum. Fur-thermore, the ROI is in an area with crossing fiber tractsand is near the cortex. A b = 1000 DWI volume is usedwith the ROI being × × at 1.25mm isotropic voxelswith 90 directions. Only the central sagittal slice along withcorresponding ODFs is visualised (Figure 14a), while thedeformation is applied to the whole ROI. The deformationfield for the central slice is shown in Figure 14b. The experimental setup provides a ground truth, which issimply the identity transformation. As the similarity mea-sure, NMI, does not reflect a unique point-wise correspon-dence, we introduce three measures for evaluating the re-sulting registration: (1) mean squared error (MSE) betweencoordinates, (2) curl of the final deformation field, and (3)divergence of the final deformation field [34].
Curl quantifiesthe amount of orientational change in the deformation field,also referred to as rotation or vorticity. We use the L -normof the curl vector, i.e. the velocity of rotation. A highermagnitude of the curl is expected for the results of a scalar-based registration compared to a registration over the fulldiffusion profile. Divergence quantifies the density of theoutward flux of a vector field which can be either positive ornegative, indicating an expansion or a contraction at a givenpoint. In combination, these three measures provide infor-mation about the magnitude of the voxel-wise distance andthe state of the final deformation field. Curl and divergencehas previously been used as regularization in a nonrigidregistration framework [35].In the following, the experiments are performed whilevarying the intensity scale, the orientation scale, and thespatial scale, as it is the intensity and orientational scale
EEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. X, NO. X, MAY 2018 8 (a) Moving Image (b) Target Image (c) Reconstructed warped image (d) Spatial mappingFig. 9. Experiment 2: Simulated crossing tracts (a) with a vertical shift (b). The registered result of the crossing under a both vertical and horizontalshift is reconstructed in (c), and shown with the final spatial mapping from the original position of the moving image in (d).(a) Moving Image (b) Target Image (c) Reconstructed warped image (d) Spatial mappingFig. 10. Experiment 2: Simulated crossing tracts with ∼
30 degree (a) to 45 degree shearing (b). The reconstruction of the registered result is shownin (c), and the corresponding spatial mapping in (d).(a) Moving Image (b) Target Image (c) Reconstructed warped image (d) Spatial mappingFig. 11. Experiment 3: (a) Simulated fanned opening in the moving image, and (b) an interleaving curved target image. The reconstruction is shownin (c), and the spatial mapping in (d).(a) Moving Image (b) Target Image (c) Reconstructed warped image (d) Spatial mappingFig. 12. Experiment 3: Simulated straight (a) and kissing (b) fiber tracts. The reconstruction is shown in (c), and the spatial mapping in (d).
EEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. X, NO. X, MAY 2018 9 (a) Sagittal slice at ROI center. (b) Deformation field.Fig. 14. Original central ROI slice with ODFs (a), and the deformation tobe applied (b). that sets this registration framework apart form other frame-works. The following is investigated:
1. Isointensity curves.
The density-based formulation al-lows us to smooth the image inverse proportionally tothe image gradients, such as the borders near the CSF,which can be strongly affected by partial volume ef-fects. We iteratively decrease the control point-spacingin the FFD model to get an increasingly localized reg-istration, and increase the intensity range accordingly.This is realized with an initial histogram with relativelyfew bins and a fixed degree of smoothing and thensuccessively increase the number of bins as the defor-mation resolution is increased of iterations.
2. Explicit reorientation.
The directional information is, ac-cording to our experiments, expected to provide a morestable and regularized transformation. To investigatethis hypothesis, different levels of directional smooth-ing are examined from κ = 30 a highly peaked ODF toa scalar-based mean diffusivity registration i.e. settingthe concentration parameter κ to . The experimentson simulated data indicated that the directional infor-mation results in direct regularized solutions throughthe similarity, and these experiments seek to uncover ifthese observations hold for real data.Performing multi-resolution registration (a continuationmethod) is a key element in most high resolution reg-istration frameworks, e.g. in the FFD model [29], ANTs [5],
Elastix [36],
FSL [7], etc. It provides stability whileallowing for large and small deformations while reducingthe computational complexity. While experiments regardingthe spatial part will be performed, its impact on registrationquality is fairly well-described in contrast to iso-intensity curves and explicit reorientation.
Hierarchical mesh resolution.
Similar to the simulated ex-periments, the spacing between the control points issuccessively decreased. To this end we use the com-position Φ local = Φ δ =10 ◦ Φ δ =5 ◦ Φ δ =3 . ◦ Φ δ =3 , wherethe control point spacing δ is scaled down according tothe dimensions of the image. The low initial resolutionof the deformation permits the generation of a largerandom and valid synthetic warp at around half of thecontrol point spacing ( ∼ . δ ) [37]. For the optimizationwe use a quasi-Newton L-BFGS optimizer which runsfor 50 iterations for all resolutions, unless the optimalitycondition of (cid:15) = 1e − is fulfilled. Note that we refer tothe result after each optimizer termination as a step , thus4 steps in total. Accumulated curl and divergence.
While the MSE is easyto calculate at each step, the curl and divergencedepend on the first-order derivatives of the spatialdeformation, which are accumulated over each step.Thus the final curl and divergence is defined by theproduct of Jacobians for each step
Curl (Φ δ =3 ) = J (Φ δ =10 ) J (Φ δ =5 ) J (Φ δ =3 . ) J (Φ δ =3 ) for all voxels. Change in resolution.
Changing spatial resolution isachieved by convolution.
Other fixed parameters.
The regularization is fixed to λ =1e − , and we interpolate and optimize over 30 interpo-lated orientations of the 90 in the HCP data.All the following experiments will show results of the fouraccumulating steps. All error measures are reported for theentire ROI and not just the slice visualized. This experiment illustrates the effect of changing the sizeof the smoothed joint histogram. Three experiments areperformed with (i) a fixed histogram of × bins, (ii) ahistogram with × bins, and (iii) gradually increasingthe histogram size in 50, 100, 200 and 500 bins. Figure 15shows the results in terms of MSE, curl and divergenceas mean value over all points along. It is evident how asmall histogram with wide isocurves results in an initiallyfaster convergence Figure 15 (blue line). However, wide iso-curves result in flat regions with small gradients and littlestructure which may cause the result to become sub-optimalor degenerate with increasing degrees of freedom in thetransformation. In contrast, starting with a high-resolutionof the histogram with thin iso-curves has the opposite effectin the initial steps (red line). However, by iteratively refiningthe joint histogram, we allow for a wide to thin movementin the image, which provides superior result (green line). In this experiment, we investigate if directional informationincreases the stability and improves the registration. Weemploy the size progression of the histogram from theprevious experiment (i.e. [50,100,200,500]), and vary theconcentration parameter κ from mean diffusivity at κ = 0 tosharp angular features at κ = 30 . The results are shown inFigure 16. The figure shows how the directional information EEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. X, NO. X, MAY 2018 10
Step M SE Bins = 50Bins = 500Bins = [50,100,200,500] (a) MSE
Step C u r l Bins = 50Bins = 500Bins = [50,100,200,500] (b) Curl
Step D i v e r gen c e Bins = 50Bins = 500Bins = [50,100,200,500] (c) (Absolute) DivergenceFig. 15. Results from testing the effect of the number of bins through the four step registration. The lines are the mean over all points. The bestresults are (a): green at . − , (b): green at . − , and (c): green at . − . Step M SE κ = 0 κ = 5 κ = 10 κ = 15 κ = 20 κ = 30 (a) MSE Step C u r l κ = 0 κ = 5 κ = 10 κ = 15 κ = 20 κ = 30 (b) Curl Step D i v e r gen c e κ = 0 κ = 5 κ = 10 κ = 15 κ = 20 κ = 30 (c) (Absolute) DivergenceFig. 16. Results from testing the effect of smoothness of the directional interpolation through the four step registration. The lines are the mean overall points. The best results are the yellow line at κ = 30 with (a): . − , (b): . − , and (c): . − . However, the overall difference between κ = 10 and κ = 30 is not more than around − . results in better registration, with significantly less curl thanthe scalar registration at κ = 0 and we observe that the bestvalue κ = 30 is stable in high directional resolution datasuch as the HCP. We suggest setting κ = 15 as this is shouldsuffice for both high and low resolution data, such as DTI. In the last experiment, we use κ = 15 , set the bins to [50 , , , , and investigate the effects of changingthe spatial scale. The spatial resolution is set to s =[4 , , , , which equivalent of smoothly interpolating every4th point, followed by every 3rd, etc. As with the controlpoint spacing of the deformation field, we scale this to fitthe image if the image is not equilateral. For instance, if theimage has the spatial dimension × × then spacebetween spatial interpolations for s = 3 will be [2 , , witha bound on no less than 1 (i.e. full resolution). The results areshown in Figure 17. The results are very similar in the finalstep. However, hierarchical resolution approach comparedto the full resolution gave a speedup of a factor 2.4. Finally, we visualize the registration to perform a qualitativeassessment of the warp, shape and orientation of the indi-vidual ODFs. All registrations are based on the raw, noise-corrected HARDI data, while we show the tomographicinversion (FRT) indicating the direction of the diffusionand likely fiber tract orientations. Unlike the simulated experiments, we have fitted a B-spline to the image priorto deforming and visualizing the result, in contrast to thesmoothing spline used in the artificial examples. The resultswhere obtained using an increasing control point resolu-tion with κ = 15 , bins = [50 , , , , and spatialresolution = [4 , , , and Figure 18 shows the result forthe central ROI slice, along with a zoomed in version inFigure 19. ONCLUSION
We have presented a scale-space formulation of density es-timation that extends LOR to spatio-directional data, explic-itly for the registration of DWI data with explicit reorienta-tion of the full diffusion profile. We have provided empiricalevidence showing that the underlying structure of the datais preserved during registration, while providing excellentregistration results through a number of classical artificialexamples known to be difficult for registration. In addition,we have shown that the formulation of the similarity itselfprovides regularization through the additional informationgiven by the orientational dimension, illustrated clearly inthe artificial examples. We have investigated the differentscales provided by the framework and shown how thedifferent parameters influences the registration results. Inconclusion, the LOR-DWI provides a smooth and well-matched deformation of DWI data, and the registrationresults are improved by the induced regularization obtainedby integrating the orientational information in the objective
EEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. X, NO. X, MAY 2018 11
Step M SE × -3 s = 1s = [4,3,2,1] (a) MSE Step C u r l s = 1s = [4,3,2,1] (b) Curl Step D i v e r gen c e s = 1s = [4,3,2,1] (c) (Absolute) DivergenceFig. 17. Results from testing the effect of changing the spatial resolution from low to full. The lines are the mean over all points. For the full (blue) andlow-to-full (red) resolution, the difference in results are (a): blue . − vs red . − , (b): blue . − vs red . − , and (c): blue . − vs red . − . (a) Ground truth (b) Warped image (c) Registered image reconstructionFig. 18. 2D visualization of Figure 14a, and the reconstructed warped image after applying the deformation field to the original image. We will beregistering (b) back to (a). (a) Ground truth(b) Registered image reconstructionFig. 19. Same as Figure 18, zoomed in on the left side (anterior) ofthe figure (around the Genu). The images are not a 100 percent thesame, but the difference is hard to notice, and the registration more thanadequate. function. Finally, we showed that the deformation have verylittle impact on the shape of the deformed ODF. A CKNOWLEDGMENTS
This research was supported by Centre for Stochastic Geom-etry and Advanced Bioimaging, funded by grant 8721 fromthe Villum Foundation. R EFERENCES [1] H. Johansen-Berg and T. E. Behrens,
Diffusion MRI: from quantita-tive measurement to in vivo neuroanatomy . Academic Press, 2013.[2] H. G. Jensen, F. Lauze, M. Nielsen, and S. Darkner, “Locallyorderless registration for diffusion weighted images,” in
Interna-tional Conference on Medical Image Computing and Computer-AssistedIntervention . Springer, 2015, pp. 305–312.[3] D. C. Van Essen, S. M. Smith, D. M. Barch, T. E. Behrens, E. Yacoub,K. Ugurbil, W.-M. H. Consortium et al. , “The wu-minn humanconnectome project: an overview,”
Neuroimage , vol. 80, pp. 62–79,2013.[4] L. J. O’Donnell, A. Daducci, D. Wassermann, and C. Lenglet,“Advances in computational and statistical diffusion mri,”
NMRin Biomedicine , 2017.[5] B. B. Avants, N. Tustison, and G. Song, “Advanced normalizationtools (ants),”
Insight j , vol. 2, pp. 1–35, 2009.[6] J. Tournier, F. Calamante, A. Connelly et al. , “Mrtrix: diffusiontractography in crossing fiber regions,”
International Journal ofImaging Systems and Technology , vol. 22, no. 1, pp. 53–66, 2012.[7] M. Jenkinson, C. F. Beckmann, T. E. Behrens, M. W. Woolrich, andS. M. Smith, “Fsl,”
Neuroimage , vol. 62, no. 2, pp. 782–790, 2012.[8] B. Fischl, “Freesurfer,”
Neuroimage , vol. 62, no. 2, pp. 774–781, 2012.[9] H. Zhang, P. A. Yushkevich, D. C. Alexander, and J. C. Gee, “De-formable registration of diffusion tensor mr images with explicitorientation optimization,”
Medical image analysis , vol. 10, no. 5, pp.764–785, 2006.[10] B. T. Yeo, T. Vercauteren, P. Fillard, J.-M. Peyrat, X. Pennec,P. Golland, N. Ayache, and O. Clatz, “Dt-refind: Diffusion tensorregistration with exact finite-strain differential,”
IEEE transactionson medical imaging , vol. 28, no. 12, pp. 1914–1928, 2009.
EEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. X, NO. X, MAY 2018 12 [11] M. O. Irfanoglu, A. Nayak, J. Jenkins, E. B. Hutchinson,N. Sadeghi, C. P. Thomas, and C. Pierpaoli, “Dr-tamas: Diffeo-morphic registration for tensor accurate alignment of anatomicalstructures,”
NeuroImage , vol. 132, pp. 439–454, 2016.[12] Y. Wang, A. Gupta, Z. Liu, H. Zhang, M. L. Escolar, J. H. Gilmore,S. Gouttard, P. Fillard, E. Maltbie, G. Gerig et al. , “Dti registrationin atlas based fiber analysis of infantile krabbe disease,”
Neuroim-age , vol. 55, no. 4, pp. 1577–1586, 2011.[13] P. Zhang, M. Niethammer, D. Shen, and P.-T. Yap, “Large defor-mation diffeomorphic registration of diffusion-weighted imagingdata,”
Medical image analysis , vol. 18, no. 8, pp. 1290–1298, 2014.[14] Y. Wang, Q. Yu, Z. Liu, T. Lei, Z. Guo, M. Qi, and Y. Fan,“Evaluation on diffusion tensor image registration algorithms,”
Multimedia Tools and Applications , vol. 75, no. 13, pp. 8105–8122,2016.[15] Y. Wang, Y. Shen, D. Liu, G. Li, Z. Guo, Y. Fan, and Y. Niu,“Evaluations of diffusion tensor image registration based on fibertractography,”
Biomedical engineering online , vol. 16, no. 1, p. 9,2017.[16] W. Wells, P. Viola, H. Atsumi, S. Nakajima, and R. Kikinis, “Multi-modal volume registration by maximization of mutual informa-tion,”
Medical Image Analysis , vol. 1, no. 1, pp. 35–51, 1996.[17] S. Darkner and J. Sporring, “Locally orderless registration,”
IEEEtransactions on pattern analysis and machine intelligence , vol. 35, no. 6,pp. 1437–1450, 2013.[18] W. Van Hecke, A. Leemans, E. D’Agostino, S. De Backer, E. Van-dervliet, P. M. Parizel, and J. Sijbers, “Nonrigid coregistration ofdiffusion tensor images using a viscous fluid model and mutualinformation,”
IEEE transactions on medical imaging , vol. 26, no. 11,pp. 1598–1612, 2007.[19] J. M. Treiber, N. S. White, T. C. Steed, H. Bartsch, D. Holland,N. Farid, C. R. McDonald, B. S. Carter, A. M. Dale, and C. C. Chen,“Characterization and correction of geometric distortions in 814diffusion weighted images,”
PloS one , vol. 11, no. 3, p. e0152472,2016.[20] C. Studholme, D. L. Hill, and D. J. Hawkes, “An overlap invariantentropy measure of 3d medical image alignment,”
Pattern recogni-tion , vol. 32, no. 1, pp. 71–86, 1999.[21] J. Sporring and S. Darkner, “Jacobians for lebesgue registration fora range of similarity measures,”
Department of Computer Science,University of Copenhagen, Tech. Rep , vol. 4, 2011.[22] S. Darkner and J. Sporring, “Locally Orderless Registration,”
IEEETransactions on Pattern Analysis and Machine Intelligence , vol. 35,no. 6, pp. 1437–1450, 2013.[23] J. J. Koenderink and A. J. Van Doorn, “The structure of locallyorderless images,”
International Journal of Computer Vision , vol. 31,no. 2, pp. 159–168, 1999.[24] G. Hermosillo, C. Chefd’Hotel, and O. Faugeras, “Variationalmethods for multimodal image matching,”
International Journal ofComputer Vision , vol. 50, no. 3, pp. 329–343, 2002.[25] S. Darkner and J. Sporring, “Generalized partial volume: An infe-rior density estimator to parzen windows for normalized mutualinformation,” in
IPMI , ser. LNCS, vol. 6801. Springer, 2011, pp.436–447.[26] X. Tao and J. V. Miller, “A method for registering diffusionweighted magnetic resonance images,” in
Medical Image Computingand Computer-Assisted Intervention–MICCAI 2006 . Springer, 2006,pp. 594–602.[27] P. Jupp and K. Mardia, “A unified view of the theory of directionalstatistics, 1975-1988,”
International Statistical Review/Revue Interna-tionale de Statistique , pp. 261–294, 1989.[28] M. Abramowitz and I. Stegun,
Handbook of Mathematical Functions,With Formulas, Graphs, and Mathematical Tables, . Dover, 1974.[29] D. Rueckert, L. I. Sonoda, C. Hayes, D. L. Hill, M. O. Leach, andD. J. Hawkes, “Nonrigid registration using free-form deforma-tions: application to breast mr images,”
IEEE transactions on medicalimaging , vol. 18, no. 8, pp. 712–721, 1999.[30] M. Schmidt, “minfunc: unconstrained differentiablemultivariate optimization in matlab,” , 2005.[31] E. Garyfallidis, M. Brett, B. Amirbekian, A. Rokem, S. VanDer Walt, M. Descoteaux, and I. Nimmo-Smith, “Dipy, a libraryfor the analysis of diffusion mri data,”
Frontiers in neuroinformatics ,vol. 8, p. 8, 2014.[32] A. Semechko, “Suite of functions to perform uni-form sampling of a sphere,” ,2012.[33] M. Descoteaux, E. Angelino, S. Fitzgibbons, and R. Deriche, “Reg-ularized, fast, and robust analytical q-ball imaging,”
Magneticresonance in medicine , vol. 58, no. 3, pp. 497–510, 2007.[34] H. Helmholtz, “Uber integrale der hydrodynamischen gleichun-gen, welcher der wirbelbewegungen entsprechen” (on integralsof the hydrodynamic equations which correspond to vortex mo-tions),”
Journal fr die reine und angewandte Mathematik , vol. 55, no. 8,pp. 25–55, 1858.[35] S. Riyahi-Alam, M. Peroni, G. Baroni, and M. Riboldi, “Regular-ization in deformable registration of biomedical images based ondivergence and curl operators,”
Methods of information in medicine ,vol. 53, no. 01, pp. 21–28, 2014.[36] S. Klein, M. Staring, K. Murphy, M. A. Viergever, and J. P. Pluim,“Elastix: a toolbox for intensity-based medical image registration,”
IEEE transactions on medical imaging , vol. 29, no. 1, pp. 196–205,2010.[37] D. Rueckert, P. Aljabar, R. A. Heckemann, J. V. Hajnal, andA. Hammers, “Diffeomorphic registration using b-splines,” in
International Conference on Medical Image Computing and Computer-Assisted Intervention . Springer, 2006, pp. 702–709.
Henrik G. Jensen received his Masters De-gree in Computer Science from the Universityof Copenhagen in 2014. Shortly after, he wasawarded an open grant and started as a PhDfellow in the same department, studying medi-cal image analysis and registration of diffusion-weighted MRI. His research interests also in-clude machine learning and data analysis, andbefore finishing his Bachelor’s Degree, he hadalready defended his first international confer-ence paper on automatic diagnosis of rare dis-eases. He has also worked as a research assistant analyzing aerialdrone footage with the purpose of automatic weed detection to reducethe use of agricultural herbicides.
Francois Lauze studied Mathematics in France,University of Nice Sophia Antipolis where re-ceived a PhD in Algebraic Geometry in 1994.He spent some years in Burkina Faso, wherehe taught Mathematics at the University of Oua-gadougou. He moved to Denmark and engagedin PhD n o
2, at the IT University of Copenhagen,which he received in 2004. He worked on Varia-tional Methods for Motion Compensated Inpaint-ing and motion recovery among other.He cur-rently holds a position as Associate Professorin Mathematical Image Analysis at the department of Computer Sci-ence, University of Copenhagen. He has mainly worked with variationalmethods, and geometry for Image Analysis - mainly Riemannian butalso some metric geometry. More recently with inverse problems inphotometric stereo and tomographic imaging.