# Extending Text-to-Speech Synthesis with Articulatory Movement Prediction
using Ultrasound Tongue Imaging

###### Abstract

In this paper, we present our first experiments in text-to-articulation
prediction, using ultrasound tongue image targets. We extend a traditional
(vocoder-based) DNN-TTS framework with predicting PCA-compressed ultrasound
images, of which the continuous tongue motion can be reconstructed in
synchrony with synthesized speech. We use the data of eight speakers, train
fully connected and recurrent neural networks, and show that FC-DNNs are more
suitable for the prediction of sequential data than LSTMs, in case of limited
training data. Objective experiments and visualized predictions show that the
proposed solution is feasible and the generated ultrasound videos are close to
natural tongue movement. Articulatory movement prediction from text input can
be useful for audiovisual speech synthesis or computer-assisted pronunciation
training.

Index Terms: DNN-TTS, audiovisual synthesis, ultrasound

##  1 Introduction

Statistical parametric methods are frequently used in text-to-speech (TTS)
synthesis, with two main machine learning techniques: hidden Markov-models
(HMM, [1]) and deep neural networks (DNN, [2]). Recently, the focus of TTS
research has moved to end-to-end solutions, applying neural vocoders (e.g.
WaveNet [3] and WaveGlow [4]) and sequence-to-sequence models using attention
(e.g. Tacotron2 [5]). Still, traditional (non-end-to-end, vocoder-based) DNN-
TTS systems are useful in limited data scenarios, when there is few training
data available, for example with speech and biosignal recordings, or in case
of audiovisual speech synthesis.

###  1.1 Audiovisual speech synthesis

Audiovisual speech synthesis is a subfield of the more general areas of speech
synthesis and computer facial animation [6]. The goal of the visible speech
synthesis is typically to obtain a mask with realistic motions, not to
duplicate the musculature of the face to control this mask.

The field of visual speech synthesis is fairly well established and a variety
of approaches have been developed (including rule-based [7] and data-driven
methods [8]). Rule-based systems include models for speech sequence planning,
for muscle mechanisms and for the physical speech production apparatus. Within
the biomechanical model of the vocal tract, the tongue can be represented as a
finite element mesh [7] and complex biomechanical simulations are necessary to
estimate the internal muscle stresses during the movement of human
articulators [9]. In the context of data-driven approaches and HMM-based
synthesis, there are two main categories: image-based systems are supposed to
look like a video of a real person, while motion capture based systems derive
features from facial points tracked over time [8]. For HMM-based audiovisual
synthesis, a synchronous corpus of parametrized facial motion data and
acoustic speech data is necessary. Schabus et al. showed that in combined HMM-
based text-to-speech synthesis and facial animation, joint audiovisual models
perform better than training separate acoustic and visual models [8].

###  1.2 Predicting articulatory movement from text

Another type of TTS extension is when the target is to predict articulatory
motion (e.g. lip or tongue movement) and not just the face of the speaker,
besides the speech output. This requires special biosignals to be recorded,
which can track the movement of the articulatory organs (e.g. EMA, X-ray,
vocal tract MRI, and ultrasound tongue imaging). With such a system, by giving
an arbitrary input text, one is able to hear the speech and, in synchrony with
it, see how to move the tongue in 3D to produce target speech sounds. This
visual feedback can make a big difference for pronunciation training in L2
learning, especially when the target language contains speech sounds which are
difficult to articulate.

Most earlier studies in this context were using point-tracking devices, like
electromagnetic articulography (EMA) [10, 11, 12, 13, 14, 15]. Ling and his
colleagues proposed a HMM-based text-to-articulatory movement prediction
system, i.e. which can synthesize the speaker’s mouth from text [10]. Here,
the durations were not modeled, but in a subsequent study, they also
investigated the timing aspects and analyzed the critical articulators [11].
Wei et al. used DNNs for the text-to-EMA prediction and confirmed that stacked
bottleneck features are effective for this purpose [12]. Steiner and his
colleagues experimented similarly with text-to-EMA prediction using HMMs (with
synchronous text-to-speech), and they also included a geometric 3D tongue
model as the target [13]. Next, they compared HMMs and DNNs for the text-to-
tongue model prediction [14]. It was found that with less than 2 hours of
data, DNNs outperformed HMMs. Yu and her colleagues predicted 3D articulatory
movement, from text and audio inputs, therefore combining the text-to-speech
and acoustic-to-articulatory inversion fields [15]. For the machine learning
approach, they used a bottleneck long-term recurrent convolutional neural
network. They showed that the text information complements well the acoustic
features during the prediction of EMA-based articulation. The final output of
the system is speech synchronized with 3D articulatory animation, using a
facial mesh model [15].

As shown above, there have been several studies investigating text-to-
articulatory motion with HMMs or DNNs, but all of them are using point-
tracking equipment (electromagnetic articulography). Medical imaging target,
like ultrasound or MRI, have not been used before in this context.

Figure 1: Block diagram of the proposed approach.

###  1.3 Ultrasound tongue imaging

Ultrasound tongue imaging (UTI) is a technique suitable for the acquisition of
articulatory data. Phonetic research has employed 2D ultrasound for a number
of years for investigating tongue movements during speech [16]. Stone
summarized the typical methodology of investigating speech production using
ultrasound [17]. Usually, when the subject is speaking, the ultrasound
transducer is placed below the chin, resulting in mid-sagittal images of the
tongue movement. The typical result of 2D ultrasound recordings is a series of
gray-scale images in which the tongue surface contour has a greater brightness
than the surrounding tissue and air. Compared to other articulatory
acquisition methods (e.g. EMA, X-ray, XRMB, and vocal tract MRI), UTI has the
advantage that the tongue surface is fully visible, and ultrasound can be
recorded in a non-invasive way [17, 18, 19]. An ultrasound device is easy to
handle and move, since it is small and light, and thus it is suitable for
fieldworks, as well. Besides, it is a significantly less expensive piece of
equipment than the above mentioned devices.

In our earlier studies, we have shown that ultrasound is a feasible solution
for articulatory-to-acoustic mapping [18, 20] and acoustic-to-articulatory
inversion [21]. However, text-to-ultrasound synthesis has not been
investigated before.

###  1.4 Contributions of this paper

The goal of this paper is to extend DNN-TTS with articulatory movement
prediction, using ultrasound images of the tongue. We show on the data of
several speakers that the combined TTS and synthesized articulatory motion is
feasible and can result in acceptable articulatory movement video. Text-to-
articulatory movement prediction might be useful for computer-assisted
pronunciation training (CAPT) applications and articulatory visualization.

##  2 Methods

###  2.1 Data

For the data, we used the UltraSuite-TaL80 database [22]
(https://ultrasuite.github.io/data/tal_corpus/). We chose four English male
(03mn, 04me, 05ms, 07me) and four female speakers (01fi, 02fe, 06fe, and
09fe). In parallel with speech, the tongue movement was recorded in
midsagittal orientation using the “Micro” ultrasound system of Articulate
Instruments Ltd. at 81.5 fps. Lip video was also recorded in UltraSuite-TaL80,
but we did not use that information in the current study. The ultrasound data
and the audio signals were synchronized using the tools provided by Articulate
Instruments Ltd. Each speaker read roughly 200 sentences – the duration of the
recordings was about 15 minutes, which was partitioned into training,
validation and test sets in a 85-10-5 ratio.

###  2.2 Processing the ultrasound data

In our experiments, articulatory features estimated from the raw scanline data
of the ultrasound were used as the additional target of the networks (see Fig.
1). The 64$\times$842 pixel images were resized to 64$\times$128 pixels using
bicubic interpolation, and we calculated PCA coefficients, similarly to
EigenTongues [23]. While calculating the PCA, we aimed at keeping the 70% of
the variance of the original images, thus having 128 coefficients. An example
for the PCA eigenvectors can be seen in Fig. 2, and the result of PCA is
presented in Fig. 4. To be in synchrony with the acoustic features (frame
shift of 5 ms), the ultrasound data was resampled to 200 Hz.

###  2.3 DNN-TTS framework

Fig. 1 illustrates the proposed approach, i.e. the combined acoustic and
articulatory feature prediction using a DNN from text input. The experiments
were conducted in the Merlin DNN-TTS framework [24] (https://github.com/CSTR-
Edinburgh/merlin). Textual / phonetic parameters are first converted to a
sequence of linguistic features as input (based on a decision tree). Next,
neural networks are employed to predict acoustic (60-dimensional MGC,
5-dimensional BAP, and 1-dimensional LF0, with delta and delta-delta features)
and articulatory features (128-dimensional ULT-PCA, with delta and delta-
delta) as output for synthesizing speech, at a 5 ms frame step with the WORLD
vocoder. From the predicted 128-dimensional articulatory features, the
64$\times$128 image is reconstructed using the PCA matrix, and bicubic
interpolation is applied to resize the image to 64$\times$842 pixels, to be
comparable with the original data. For visualization purposes, we transformed
this raw scanline data to wedge format, which shows how the real aspect ratio
of the tongue surface (for an example, see Fig. 4). The transformation was
done with ’ultrasuite-tools’ (https://github.com/UltraSuite/ultrasuite-tools)

####  2.3.1 FC-DNN

The DNN used here is a fully-connected feed-forward multilayer perceptron
architecture (FC-DNN, six hidden layers, 1024 neurons in each). We applied
tangent hyperbolic activation function, SGD optimizer, and a batch size of
256. The input features had min-max normalization, while output acoustic
features had mean-variance normalization. We trained the networks for 25
epochs with a warm-up of 10 epochs, applying early stopping, and a learning
rate of 0.002 after that with exponential decay. We only trained both a
duration model and an acoustic model, the latter also containing the
articulatory features.

####  2.3.2 LSTM

Recurrent networks are typically more suitable to process sequential data.
Therefore, we also trained an LSTM network following the Merlin recipe (four
FF layers with 1024 neurons each, and one LSTM layer with 512 neurons). To
ensure a longer training with the recurrent network, we used ADAM optimizer,
and a warm-up of 30 epochs with early stopping. The other parameters were the
same as for the FF-DDN. We trained both duration and acoustic models.

All neural network trainings were done individually with each speakers’ data,
without average voice training or adaptation.

##  3 Experimental results

After training the above models, we synthesized sentences from the test parts
of the ultrasound datasets. To measure the validation and test error, we
calculated both spectral prediction error (Mel-Cepstral Distortion, MCD), and
an articulatory feature related error (ULT-PCA/RMSE, calculated on the
normalized PCA values). We trained both duration and acoustic models, but for
the error calculations, we synthesized the test sentences with their original
timing. This way, warping the features in time was not necessary for
calculating the error measures.

###  3.1 Demonstration samples

An example for the PCA eigenvectors are in Fig. 2, and the predicted
articulatory feature sequence can be seen in Fig. 3 (speaker ’01fi’, sentence
’201_aud’). As lower dimensional PCA vectors contain more information, the
visualization was done in an exponential way and only 8 dimensions are plotted
out of 128. Clearly, both the FC-DNN and the LSTM are following the tendencies
found in the original data (e.g. in case of PCA-1, PCA-2, PCA-4), but the fine
details are not modeled well. This type of oversmoothing is a frequent
phenomena in statistical parametric speech synthesis. The higher dimensions
(e.g. PCA-64 and PCA-128) are almost constant; i.e. they could not be modeled
well with neural networks.

To visualize the individual ultrasound images, we plotted several ultrasound
frames from the original videos and from the predicted ones, as a function of
time, in Fig. 4. The reason why we are plotting every third frame is that for
the 5 ms frame step of the Merlin toolkit, the 81.5 fps ultrasound video was
interpolated to 200 Hz, and therefore, in the predicted data, roughly every
3rd frame contains visible tongue motion. In case of speaker ’01fi’, we can
see in the top row (original ultrasound images after PCA) that there is a
significant tongue movement, i.e. the tongue tip (on the right) goes higher,
as the time passes. Both the predictions with the DNN and LSTM network follow
the articulatory movement, but the images are smoothed – again, resulting from
the statistical training. For speaker ’03mn’, similar tendencies can be
observed: the movement of the tongue is changing its curvature as a function
of time, but in the DNN-predicted and LSTM-predicted images, the tongue
surface is not as clear as in the original data.

As the synthesized motion of the tongue is more visible in real-time videos,
we made available several samples at
http://smartlab.tmit.bme.hu/ssw11_txt2ult.

Figure 2: PCA eigenvectors, from speaker ’01fi’. Top: raw, scanline data
(resized to 64$\times$128 pixels). Bottom: wedge orientation. Figure 3:
Original and predicted articulatory features, from speaker ’01fi’. Sentence:
’”I leave it to nobody,” said Shakespeare, sulkily.’

Figure 4: Original and predicted articulatory feature sequences. Top: speaker
’01fi’, bottom: speaker ’03mn’. The numbers at the bottom-left correspond to
the frame number of the video.

###  3.2 Objective measures

Table 1 summarizes the MCD results. The MCD values of the test sentences with
the FC-DNN are between 5.8–7.0 dB (average: 6.228 dB), whereas with LSTM they
are between 6.0–7.5 dB (average: 6.593 dB), indicating that the recurrent
neural network was not helpful in estimating the acoustic features. The reason
for this might be that we have limited articulatory-acoustic databases
(roughly 200 sentences for each speaker), which is too small for training an
LSTM model.

The results of the RMSE calculated on the articulatory feature are summarized
in Table 2. The lowest error was achieved with the data of speaker 09fe: with
FC-DNN, the test error is 2.9, while with LSTM, the test error is 3.1. The
tendency is similar to the case of MCD: the LSTM network was not helpful in
predicting the articulatory features, probably due to the small size of the
data.

Table 1: MCD errors on the dev/test set. | MCD
---|---
Speaker | FC-DNN | LSTM
01fi | 6.995 / 6.971 | 6.647 / 6.588
02fe | 6.095 / 5.803 | 6.486 / 6.259
03mn | 5.781 / 5.785 | 5.977 / 5.948
04me | 5.896 / 6.024 | 6.318 / 6.312
05ms | 6.244 / 6.256 | 7.235 / 7.083
06fe | 5.758 / 5.582 | 6.444 / 6.330
07me | 6.589 / 6.562 | 6.831 / 6.749
09fe | 6.516 / 6.844 | 7.197 / 7.472
average | 6.234 / 6.228 | 6.642 / 6.593
Table 2: ULTPCA/RMSE errors on the dev/test set. | ULTPCA128/RMSE
---|---
Speaker | FC-DNN | LSTM
01fi | 3.292 / 3.223 | 3.319 / 3.208
02fe | 3.533 / 3.732 | 3.753 / 3.904
03mn | 3.147 / 3.660 | 3.289 / 3.680
04me | 3.849 / 3.985 | 4.031 / 4.033
05ms | 3.133 / 3.233 | 3.249 / 3.405
06fe | 3.439 / 3.250 | 3.743 / 3.451
07me | 3.544 / 3.595 | 3.498 / 3.461
09fe | 3.022 / 2.864 | 3.234 / 3.133
average | 3.370 / 3.443 | 3.515 / 3.534

##  4 Discussion and conclusions

We have shown above that text-to-ultrasound video prediction is feasible as an
extension to traditional DNN-based text-to-speech synthesis, despite the
relatively small amount of training data. Although the synchrony between
visual and speech output is not enforced by the model, the tied acoustic and
articulatory features during the DNN training ensure that the audio and visual
features are in synchrony, i.e. that in the generated ultrasound videos, the
tongue is moving according to the synthesized speech. To objectively check
this, SyncNet, part of Wav2Lip could be applied to assess synchrony [25]. Our
paper found that the joint learning of both acoustic and articulatory features
has advantages, but this is not substantiated – a comparison of the joint
model against two separate models remains future work.

Although there have been several earlier attempts for extending text-to-speech
synthesis with articulatory data, all of these studies were using EMA, being a
point tracking equipment [10, 11, 12, 13, 14, 15], and containing less spatial
information about the tongue than ultrasound. The advantage of ultrasound in
this context is that the resulting video shows a larger portion of the tongue,
compared to EMA.

Articulatory movement prediction from text input can be useful for audiovisual
speech synthesis. A specific application is computer-assisted pronunciation
training / computer-aided language learning [26, 27, 28], which can be
beneficial for learners of second languages. With such a combined TTS and
text-to-articulatory prediction system, by giving an arbitrary input text, one
is able to hear the speech and, in synchrony with it, see how to move the
tongue in 2D or 3D to produce target speech sounds. This visual feedback can
be helpful for pronunciation training in L2 learning, especially when the
target language contains speech sounds which are difficult to articulate.

In the future, we plan to investigate speaker adaptation and speaker-
independent training. For this, a common articulatory space has to be defined,
as the currently used PCA representation is specific for each individual
speaker. Also, multi-task learning might be useful in this context: a system
could potentially be pre-trained on speech-only material which is easier to
acquire, and the UTI be trained only in addition. Besides, we plan to
investigate the effect of the misalignments in the ultrasound transducer
position [29, 30] on the text-to-ultrasound prediction results.

The code is accessible at https://github.com/BME-SmartLab/txt2ult.

##  5 Acknowledgements

The author was partly funded by the National Research, Development and
Innovation Office of Hungary (FK 124584 and PD 127915 grants). We would like
to thank CSTR for providing the Merlin toolkit and the UltraSuite-TaL
articulatory database.

## References

  * [1] H. Zen, T. Nose, J. Yamagishi, S. Sako, T. Masuko, and A. Black, “The HMM-based speech synthesis system version 2.0,” in _Proc. ISCA SSW6_ , Bonn, Germany, 2007, pp. 294–299.
  * [2] H. Zen, A. Senior, and M. Schuster, “Statistical parametric speech synthesis using deep neural networks,” in _Proc. ICASSP_ , Vancouver, Canada, 2013, pp. 7962–7966.
  * [3] A. van den Oord, S. Dieleman, H. Zen, K. Simonyan, O. Vinyals, A. Graves, N. Kalchbrenner, A. W. Senior, and K. Kavukcuoglu, “WaveNet: A Generative Model for Raw Audio,” _CoRR_ , vol. abs/1609.0, 2016.
  * [4] R. Prenger, R. Valle, and B. Catanzaro, “Waveglow: A Flow-based Generative Network for Speech Synthesis,” in _Proc. ICASSP_ , Brighton, UK, 2019, pp. 3617–3621.
  * [5] J. Shen, R. Pang, R. J. Weiss, M. Schuster, N. Jaitly, Z. Yang, Z. Chen, Y. Zhang, Y. Wang, R. Skerrv-Ryan, R. A. Saurous, Y. Agiomvrgiannakis, and Y. Wu, “Natural TTS Synthesis by Conditioning Wavenet on MEL Spectrogram Predictions,” in _Proc. ICASSP_ , Calgary, Canada, 2018, pp. 4779–4783.
  * [6] D. W. Massaro, M. M. Cohen, M. Tabain, J. Beskow, and R. Clark, “Animated speech: research progress and applications,” in _Audiovisual Speech Processing_ , G. Bailly, P. Perrier, and E. Vatikiotis, Eds. Camridge, UK: Cambridge University Press, 2012, pp. 309–345.
  * [7] P. Perrier, “”GEPPETO”: A target-based model of speech production including optimal planning and physical modeling,” in _Adventures in Speech Science_ , Tokyo, Japan, 2014.
  * [8] D. Schabus, M. Pucher, and G. Hofer, “Joint audiovisual Hidden Semi-Markov Model-based speech synthesis,” _IEEE Journal on Selected Topics in Signal Processing_ , vol. 8, no. 2, pp. 336–347, 2014.
  * [9] I. Stavness, M. A. Nazari, F. Cormac, P. Perrier, Y. Payan, J. Lloyd, and S. Fels, “Coupled Biomechanical Modeling of the Face, Jaw, Skull, Tongue, and Hyoid Bone,” in _3D Multiscale Physiological Human_ , R. O. C. H. F. E. Magnenat-Thalmann Nadia, Ed. Springer London, 2014, pp. 253–274.
  * [10] Z.-H. Ling, K. Richmond, and J. Yamagishi, “An Analysis of HMM-based prediction of articulatory movements,” _Speech Communication_ , vol. 52, no. 10, pp. 834–846, oct 2010.
  * [11] ——, “HMM-Based Text-to-Articulatory-Movement Prediction and Analysis of Critical Articulators,” in _Proc. Interspeech_ , Makuhari, Japan, 2010, pp. 2194–2197.
  * [12] Z. Wei, Z. Wu, and L. Xie, “Predicting articulatory movement from text using deep architecture with stacked bottleneck features,” in _Proc. APSIPA_ , Jeju, South Korea, 2016, pp. 1–6.
  * [13] I. Steiner, S. Le Maguer, and A. Hewer, “Synthesis of Tongue Motion and Acoustics from Text Using a Multimodal Articulatory Database,” _IEEE/ACM Transactions on Audio Speech and Language Processing_ , vol. 25, no. 12, pp. 2351–2361, 2017.
  * [14] S. Le Maguer, I. Steiner, and A. Hewer, “An HMM/DNN comparison for synchronized text-to-speech and tongue motion synthesis,” in _Proc. Interspeech_. Stockholm, Sweden: International Speech Communication Association, 2017, pp. 239–243.
  * [15] L. Yu, J. Yu, and Q. Ling, “BLTRCNN Based 3D Articulatory Movement Prediction: Learning Articulatory Synchronicity From Both Text and Audio Inputs,” _IEEE Transactions on Multimedia_ , vol. PP, no. c, p. 1, 2018\.
  * [16] M. Stone, B. Sonies, T. Shawker, G. Weiss, and L. Nadel, “Analysis of real-time ultrasound images of tongue configuration using a grid-digitizing system,” _Journal of Phonetics_ , vol. 11, pp. 207–218, 1983.
  * [17] M. Stone, “A guide to analysing tongue motion from ultrasound images,” _Clinical Linguistics and Phonetics_ , vol. 19, no. 6-7, pp. 455–501, jan 2005.
  * [18] T. G. Csapó, T. Grósz, G. Gosztolya, L. Tóth, and A. Markó, “DNN-Based Ultrasound-to-Speech Conversion for a Silent Speech Interface,” in _Proc. Interspeech_ , Stockholm, Sweden, 2017, pp. 3672–3676.
  * [19] V. Ramanarayanan, S. Tilsen, M. Proctor, J. Töger, L. Goldstein, K. S. Nayak, and S. Narayanan, “Analysis of speech production real-time MRI,” _Computer Speech and Language_ , vol. 52, pp. 1–22, 2018.
  * [20] T. G. Csapó, C. Zainkó, L. Tóth, G. Gosztolya, and A. Markó, “Ultrasound-based Articulatory-to-Acoustic Mapping with WaveGlow Speech Synthesis,” in _Proc. Interspeech_ , 2020, pp. 2727–2731.
  * [21] D. Porras, A. Sepúlveda-Sepúlveda, and T. G. Csapó, “DNN-based Acoustic-to-Articulatory Inversion using Ultrasound Tongue Imaging,” in _International Joint Conference on Neural Networks_ , Budapest, Hungary, 2019, pp. N–19 221.
  * [22] M. S. Ribeiro, J. Sanger, J.-X. X. Zhang, A. Eshky, A. Wrench, K. Richmond, and S. Renals, “TaL: a synchronised multi-speaker corpus of ultrasound tongue imaging, audio, and lip videos,” in _2021 IEEE Spoken Language Technology Workshop (SLT)_ , Shenzhen, China, 2021, pp. 1109–1116.
  * [23] T. Hueber, G. Aversano, G. Chollet, B. Denby, G. Dreyfus, Y. Oussar, P. Roussel, and M. Stone, “Eigentongue feature extraction for an ultrasound-based silent speech interface,” in _Proc. ICASSP_ , Honolulu, HI, USA, 2007, pp. 1245–1248.
  * [24] Z. Wu, O. Watts, and S. King, “Merlin: An Open Source Neural Network Speech Synthesis System,” in _9th ISCA Speech Synthesis Workshop_. Sunnyvale, CA, USA: ISCA, sep 2016, pp. 202–207.
  * [25] K. R. Prajwal, R. Mukhopadhyay, V. P. Namboodiri, and C. V. Jawahar, “A Lip Sync Expert Is All You Need for Speech to Lip Generation In the Wild,” in _Proceedings of the 28th ACM International Conference on Multimedia_ , New York, NY, USA, 2020, pp. 484–492.
  * [26] W. Katz, T. Campbell, J. Wang, E. Farrar, J. Eubanks, A. Balasubramanian, B. Prabhakaran, and R. Rennaker, “Opti-speech: A real-time, 3D visual feedback system for speech training,” in _Proc. Interspeech_ , Singapore, Singapore, 2014, pp. 1174–1178.
  * [27] D. Jones, “Development of Kinematic Templates for Automatic Pronunciation Assessment Using Acoustic-to-Articulatory Inversion,” _Master’s Thesis_ , jul 2017.
  * [28] C. Agarwal and P. Chakraborty, “A review of tools and techniques for computer aided pronunciation training (CAPT) in English,” _Education and Information Technologies_ , vol. 24, no. 6, pp. 3731–3743, nov 2019.
  * [29] T. G. Csapó and K. Xu, “Quantification of Transducer Misalignment in Ultrasound Tongue Imaging,” in _Proc. Interspeech_ , online, 2020, pp. 3735–3739.
  * [30] T. G. Csapó, K. Xu, A. Deme, T. E. Gráczi, and A. Markó, “Transducer Misalignment in Ultrasound Tongue Imaging,” in _12th International Seminar on Speech Production_ , 2020.


capbtabboxtable[][] 11institutetext: 1 LIVIA, École de technologie supérieure
(ÉTS), Montreal, Canada
2 Aramis project-team, Inria Paris, Sorbonne Universités, UPMC Univ Paris 06,
Inserm, CNRS, Institut du cerveau et la moelle (ICM) - Hôpital Pitié-
Salpêtrière, Boulevard de l′hôpital, F-75013, Paris, France
3 AP-HP, Departments of Neurology and Neuroradiology, Hôpital Pitié-
Salpêtrière, 75013, Paris, France
4 LTCI lab - IMAGES group, Télécom ParisTech, Paris, France
5 Université de Montpellier, France

# White matter fiber segmentation using functional varifolds

Kuldeep Kumar1,2 Pietro Gori2,4 Benjamin Charlier2,5 Stanley Durrleman2
Olivier Colliot2,3 Christian Desrosiers1

###### Abstract

The extraction of fibers from dMRI data typically produces a large number of
fibers, it is common to group fibers into bundles. To this end, many
specialized distance measures, such as MCP, have been used for fiber
similarity. However, these distance based approaches require point-wise
correspondence and focus only on the geometry of the fibers. Recent
publications have highlighted that using microstructure measures along fibers
improves tractography analysis. Also, many neurodegenerative diseases
impacting white matter require the study of microstructure measures as well as
the white matter geometry. Motivated by these, we propose to use a novel
computational model for fibers, called functional varifolds, characterized by
a metric that considers both the geometry and microstructure measure (e.g.
GFA) along the fiber pathway. We use it to cluster fibers with a dictionary
learning and sparse coding-based framework, and present a preliminary analysis
using HCP data.

##  1 Introduction

Recent advances in diffusion magnetic resonance imaging (dMRI) analysis have
led to the development of powerful techniques for the non-invasive
investigation of white matter connectivity in the human brain. By measuring
the diffusion of water molecules along white matter fibers, dMRI can help
identify connection pathways in the brain and better understand neurological
diseases related to white matter [6]. Since the extraction of fibers from dMRI
data, known as tractography, typically produces a large number of fibers, it
is common to group these fibers into larger clusters called _bundles_.
Clustering fibers is also essential for the creation of white matter atlases,
visualization, and statistical analysis of microstructure measures along
tracts [12].

Most fiber clustering methods use specialized distance measures, such as Mean
Closest Points (MCP) distance [4, 11]. However, these distance-based
approaches require point-wise correspondence between fibers and only consider
fiber geometry. Another important aspect for white matter characterization is
the statistical analysis of microstructure measures. As highlighted in recent
publications, using microstructure measures along fibers improves
tractographic analysis [3, 10, 12, 15, 16, 17]. Motivated by these, we propose
to use a novel computational model for fibers, called functional varifolds,
characterized by a metric that considers both the geometry and microstructure
measure (e.g. generalized fractional anisotropy) along fiber pathways.

Motivation for this work comes from the fact that the integrity of white
matter is an important factor underlying many cognitive and neurological
disorders. In vivo, tissue properties may vary along each tract for several
reasons: different populations of axons enter and exit the tract, and disease
can strike at local positions within the tract. Hence, understanding diffusion
measures along each fiber tract (i.e., tract profile) may reveal new insights
into white matter organization, function, and disease that are not obvious
from mean measures of that tract or from the tract geometry alone [3, 17].
Recently, many approaches have been proposed for tract based morphometry [12],
which perform statistical analysis of microstructure measures along major
tracts after establishing fiber correspondences. While studies highlight the
importance of microstructure measures, most approaches either consider the
geometry or signal along tracts, but not both. The intuitive approach would be
to consider microstructure signal during clustering also. However, this has
been elusive due to lack of appropriate framework.

As a potential solution, we explore a novel computational model for fibers,
called functional varifolds [1], which is a generalization of the varifolds
framework [2]. The advantages of using functional varifolds are as follows.
First, functional varifolds can model the fiber geometry as well as signal
along the fibers. Also, it does not require pointwise correspondences between
fibers. Lastly, fibers do not need to have the same orientation as in the
framework of currents [5]. We test the impact of this new computational model
on a fiber clustering task, and compare its performance against existing
approaches for this task.

As clustering method, we reformulate the dictionary learning and sparse coding
based framework proposed in [8, 7, 9]. This choice of framework is driven by
its ability to describe the entire data-set of fibers in a compact dictionary
of prototypes. Bundles are encoded as sparse non-negative combinations of
multiple dictionary prototypes. This alleviates the need for explicit
representation of a bundle centroid, which may not be defined or may not
represent an actual object. Also, sparse coding allows assigning single fibers
to multiple bundles, thus providing a soft clustering.

The contributions of this paper are threefold: 1) a novel computational model
for modeling both fiber geometry and signal along fibers, 2) a generalized
clustering framework, based on dictionary learning and sparse coding, adapted
to the computational models, and 3) a comprehensive comparison of fully-
unsupervised models for clustering fibers.

##  2 White matter fiber segmentation using functional varifolds

###  2.1 Modeling fibers using functional varifolds

In the framework of functional varifolds [1, 2], a fiber $X$ is assumed to be
a polygonal line of $P$ segments described by their center point
$x_{p}\in\mathbb{R}^{3}$ and tangent vector
$\boldsymbol{\beta_{p}}\in\mathbb{R}^{3}$ centered at $x_{p}$ and of length
$c_{p}$ (respectively, $y_{q}\in\mathbb{R}^{3}$,
$\boldsymbol{\gamma_{q}}\in\mathbb{R}^{3}$ and $d_{q}$ for a fiber $Y$ with
$Q$ segments). Let $f_{p}$ and $g_{p}$ be the signal values at center points
$x_{p}$ and $y_{q}$ respectively, and $\omega$ the vector field belonging to a
reproducing kernel Hilbert space (RKHS) $W^{*}$. Then the fibers $X$ and $Y$
can be modeled based on functional varifolds as:
$V_{(X,f)}(\omega)\approx\sum_{p=1}^{P}\omega(x_{p},\boldsymbol{\beta_{p}},f_{p})c_{p}$
and
$V_{(Y,g)}(\omega)\approx\sum_{q=1}^{Q}\omega(y_{q},\boldsymbol{\gamma_{q}},g_{p})d_{q}$.
More details can be found in [1].

The inner product metric between $X$ and $Y$ is defined as:

$\langle V_{(X,f)},V_{(Y,g)}\rangle_{W^{*}}\ =\
\sum_{p=1}^{P}\sum_{q=1}^{Q}\kappa_{f}(f_{p},g_{q})\kappa_{x}(x_{p},y_{q})\kappa_{\beta}(\boldsymbol{\beta_{p}},\boldsymbol{\gamma_{q}})c_{p}d_{q}$
(1)

where $\kappa_{f}$ and $\kappa_{x}$ are Gaussian kernels and $\kappa_{\beta}$
is a Cauchy-Binet kernel. This can be re-written as:

$\langle
V_{(X,f)},V_{(Y,g)}\rangle_{W^{*}}\,=\,\sum_{p=1}^{P}\sum_{q=1}^{Q}\\!\exp\Big{(}\frac{-\|f_{p}-g_{q}\|^{2}}{\lambda_{M}^{2}}\Big{)}\exp\Big{(}\frac{-\|x_{p}-y_{q}\|^{2}}{\lambda_{W}^{2}}\Big{)}{\Big{(}\frac{\boldsymbol{\beta_{p}}^{T}\boldsymbol{\gamma_{q}}}{c_{p}\,d_{q}}\Big{)}}^{2}c_{p}\,d_{q}$
(2)

where $\lambda_{M}$ and $\lambda_{W}$ are kernel bandwidth parameters. For
varifolds [2], a computational model using only fiber geometry and used for
comparison in the experiments, we drop the signal values at center points.
Thus, the varifolds-based representation of fibers will be:
$V_{X}(\omega)\approx\sum_{p=1}^{P}\omega(x_{p},\boldsymbol{\beta_{p}})c_{p}$
and
$V_{Y}(\omega)\approx\sum_{q=1}^{Q}\omega(y_{q},\boldsymbol{\gamma_{q}})d_{q}$.
Hence, the inner product is defined as:

$\langle
V_{X},V_{Y}\rangle_{W^{*}}\,=\,\sum_{p=1}^{P}\sum_{q=1}^{Q}\exp\Big{(}\frac{-\|x_{p}-y_{q}\|^{2}}{\lambda_{W}^{2}}\Big{)}{\Big{(}\frac{\boldsymbol{\beta_{p}}^{T}\boldsymbol{\gamma_{q}}}{c_{p}\,d_{q}}\Big{)}}^{2}c_{p}\,d_{q}.$
(3)

###  2.2 Fiber Clustering using Dictionary learning and sparse coding

For fiber clustering, we extend the dictionary learning and sparse coding
based framework presented in [8, 7, 9]. Let $V_{T}$ be the set of $n$ fibers
modeled using functional varifolds, $A\in\mathbb{R}_{+}^{n\times m}$ be the
atom matrix representing the dictionary coefficients for each fiber belonging
to one of the $m$ bundles, and $W\in\mathbb{R}_{+}^{m\times n}$ be the cluster
membership matrix containing the sparse codes for each fiber. Instead of
explicitly representing bundle prototypes, each bundle is expressed as a
linear combination of all fibers. The dictionary is then defined as $D\ =\
V_{T}A$. Since this operation is linear, it is defined for functional
varifolds.

The problem of dictionary learning using sparse coding [8, 7] can be expressed
as finding the matrix $A$ of $m$ bundle prototypes and the fiber-to-bundle
assignment matrix $W$ that minimize the following cost function:

$\operatorname*{arg\,min}_{A,W}\ \
\frac{1}{2}||V_{T}-V_{T}AW||_{W^{*}}^{2},\quad\textrm{subject to:}\ \
||\boldsymbol{w}_{i}||_{0}\leq S_{\mathrm{max}}.$ (4)

Parameter $S_{\mathrm{max}}$ defines the maximum number of non-zero elements
in $\boldsymbol{w}_{i}$ (i.e., the sparsity level), and is provided by the
user as input to the clustering method.

An important advantage of using the above formulation is that the
reconstruction error term only requires inner product between the varifolds.
Let $Q\in\mathbb{R}^{n\times n}$ be the Gram matrix denoting inner product
between all pairs of training fibers, i.e., $Q_{ij}=\langle
V_{X_{i},f_{i}},V_{X_{j},f_{j}}\rangle_{W^{*}}$. Matrix $Q$ can be calculated
once and stored for further computations. The problem then reduces to linear
algebra operations involving matrix multiplications. The solution of Eq. (4)
is obtained by alternating between sparse coding and dictionary update [8].
The sparse codes of each fiber can be updated independently by solving the
following sub-problem:

$\operatorname*{arg\,min}_{\boldsymbol{w}_{i}\in\mathbb{R}_{+}^{m}}\ \
\frac{1}{2}||V_{X_{i}}-V_{T}A\boldsymbol{w}_{i}||_{W^{*}}^{2},\quad\textrm{subject
to:}\ \ ||\boldsymbol{w}_{i}||_{0}\leq S_{\mathrm{max}}.$ (5)

which can be re-written as:

$\operatorname*{arg\,min}_{\boldsymbol{w}_{i}\in\mathbb{R}_{+}^{m}}\ \
\frac{1}{2}\big{(}Q(i,i)+\boldsymbol{w}_{i}^{\top}A^{\top}QA\boldsymbol{w}_{i}-2Q(i,:)A\boldsymbol{w}_{i}\big{)},\quad\textrm{s.t.:}\
\ ||\boldsymbol{w}_{i}||_{0}\leq S_{\mathrm{max}}.$ (6)

The non-negative weights $\boldsymbol{w}_{i}$ can be obtained using the
kernelized Orthogonal Matching Pursuit (kOMP) approach proposed in [8], where
the most _positively_ correlated atom is selected at each iteration, and the
sparse weights $\boldsymbol{w}_{s}$ are obtained by solving a non-negative
regression problem. Note that, since the size of $\boldsymbol{w}_{s}$ is
bounded by $S_{\mathrm{max}}$, it can be otained rapidly. Also, in case of a
large number of fibers, the Nystrom method can be used for approximating the
Gram matrix [7]. For dictionary update, $A$ is recomputed by applying the
following update scheme, until convergence:

$A_{ij}\ \leftarrow\
A_{ij}\dfrac{\left(QW^{\top}\right)_{ij}}{\left(QAWW^{\top}\right)_{ij}},\quad
i=1,\ldots,n,\quad j=1,\ldots,m.$ (7)

##  3 Experiments

#### Data:

We evaluate different computational models on the dMRI data of 10 unrelated
subjects (6 females and 4 males, age 22-35) from the Human Connectome Project
(HCP) [14]. DSI Studio [18] was used for the signal reconstruction (in MNI
space, $1$mm), and streamline tracking employed to generate $50,000$ fibers
per subject (minimum length $50$ mm, maximum length $300$ mm). Generalized
Fractional Anisotropy (GFA), which extends standard fractional anisotropy to
orientation distribution functions, was considered as along-tract measure of
microstructure. While we report results obtained with GFA, any other along-
tract measure may have been used.

#### Parameter impact:

We performed k-means clustering and manually selected pairs of fibers from
clusters most similar to major bundles. We then modeled these fibers using
different computational models, and analyzed the impact of varying the kernel
bandwidth parameters. The range of these parameters were estimated by
observing the values of distance between centers of fiber segments and
difference between along tract GFA values for selected multiple pairs of
fibers. Figure 1 (top left) shows GFA color-coded fibers for $3$ pairs
corresponding to a) right Corticospinal tract – CST (R), b) Corpus Callosum –
CC, and c) right Inferior Fronto-Occipital Fasciculus – IFOF (R). Cosine
similarity (in degrees) is reported for the fiber pairs modeled using
varifolds (Var) and functional varifolds (fVar), for $\lambda_{W}$ = 7 mm and
$\lambda_{M}$ = 0.01.

|
---|---
|
|
Figure 1: Along-fiber GFA visualization and cosine similarity between pairs of
fibers from three prominent bundles: a) CST (R), b) CC, c) IFOF (R), using
framework of varifolds (Var) and functional varifolds (fVar) (top left), and
Comparing variation of cosine similarity for the select fiber pairs over
kernel bandwidth parameters $\lambda_{W}$ and $\lambda_{M}$ for the framework
of functional varifolds (top right: CST (R), middle left: CC, middle right:
IFOF (R)); Impact of $\lambda_{M}$ on clustering consistency (measured using
Average Silhouette) for $m=100,125,150$ for functional Varifolds vs Varifolds
(bottom left), and functional Varifolds vs GFA only (bottom right)

Figure 1 (top left) shows GFA color-coded fiber pairs. The color-coded
visualization reflect the variation of fiber geometry, microstructure measure
(i.e. GFA) along fiber, and difference in GFA along fiber for the select fiber
pairs. This visualization of variation and difference in GFA values along
fibers support our hypothesis that modeling along tract signal along with
geometry provides additional information. The change in cosine similarity for
CC from $45.8$ degrees (using varifolds) to $66.3$ degrees (using functional
varifolds) while for CST (R) from $45.6$ degrees to $72.4$ degrees, reflect
more drop in cosine similarity if along tract signal profiles are not similar.
This shows that functional varifolds imposes penalty for different along fiber
signal profiles.

Figure 1 also compares the impact of varying the kernel bandwidth parameters
for functional varifolds using similarity angle between pairs of these
selected fibers (top right: CST (R), bottom left: CC, bottom right: IFOF (R)).
We show variation over $\lambda_{W}$ = 3, 5, 7, 9 and 11 (mm) and
$\lambda_{M}$ = 0.001, 0.005, 0.01, 0.05, and 0.1.

|  Model |  _m_ =100 |  _m_ =125 |  _m_ =150
---|---|---|---
fVar |  0.3624 |  0.3451 |  0.3314
Var |  0.3356 |  0.3089 |  0.2905
GFA |  -0.0579 |  -0.0584 |  -0.0610
MCP |  0.3240 |  0.2888 |  0.2619
Figure 2: Mean silhouette obtained with Varifolds, Varifolds, GFA, and MCP,
computed for varying a number of clusters, over 10 subjects and 3 seed values
(_left)_. Detailed results obtained for 10 subjects using $m$=100 (_right_).

Comparing the parameter variation images in Figure 1 we observe that the
cosine similarity values over the parameter space show similar trends for all
$3$ pairs of fibers. This observation allows us to select a single pair of
parameter values for our experiments. We have used $\lambda_{W}$ = $7$ mm and
$\lambda_{M}$ = $0.01$ for our experiments based on the cosine similarity
values in Figure 1. The smaller values for $\lambda_{W}$ ($<7$mm) and
$\lambda_{M}$ ($<0.01$mm) will make the current fiber pairs orthogonal while
for larger values we lose the discriminative power as all fiber pairs will
have very high similarity.

#### Quantitative analysis:

We report a quantitative evaluation of clusterings obtained using as
functional varifolds (fVar), varifolds (var), MCP and GFA computational model.
The same dictionary learning and sparse coding framework is applied for all
computational models. For each of the $10$ HCP subjects, we compute the
Gramian matrix using $5,000$ fibers randomly sampled over the full brain for
$3$ seed values. The MCP distance $d_{ij}$ is calculated between each fiber
pair $(i,j)$, as described in [4], and the Gramian matrix obtained using a
radial basis function (RBF) kernel:
$k_{ij}=\exp\big{(}\\!-\\!\gamma\\!\cdot\\!d_{ij}^{2}\big{)}$. Parameter
$\gamma$ was set empirically to $0.007$ in our experiments.

fVar Var MCP

Figure 3: Full clustering visualization ($m$ = 100, top row), single cluster
visualization (mid row), and GFA based color coded visualization of the
selected single cluster (bottom row). Using following computational models for
fibers: functional varifolds (left column), varifolds (middle column), and MCP
distance (right column). Superior axial views. Note: (top row) each figure has
a unique color code.

Since our evaluation is performed in an unsupervised setting, we use the
silhouette measure [11, 13] to assess and comparing clustering consistency.
Silhouette values, which range from $-1$ to $1$, measure how similar an object
is to its own cluster (cohesion) compared to other clusters (separation).
Figure 1 (bottom row) shows impact of $\lambda_{M}$ on clustering consistency
for functional Varifolds w.r.t Varifolds and GFA only. Figure 2 (right) gives
the average silhouette for $m$ = 100, 125, and 150 clusters, computed over
$10$ subjects and $3$ seed values. The impact of using both geometry and
microstructure measures along fibers is evaluated quantitatively by comparing
clusterings based on functional varifolds with those obtained using only
geometry (i.e., varifolds, MCP), and only along-fiber signal (i.e., GFA). As
can be seen, using GFA alone leads to poor clusterings, as reflected by the
negative silhouette values. Comparing functional varifolds with varifolds and
GFA, we observe a consistently improved performance for different numbers of
clusters. To further validate this hypothesis, we also report the average
silhouette (over $3$ seed values) obtained for $10$ subjects using $m$ = 100.
These results demonstrate that functional varifolds give consistently better
clustering, compared to other computational models using the same
framework111Silhouette analyzes only clustering consistency, not the along-
fiber signal profile..

#### Qualitative visualization:

Figure 3 (top row) shows the dictionary learned for a single subject ($m$ =
100) using functional varifolds (fVar), varifolds (Var), and MCP distance. For
visualization purposes, each fiber is assigned to a single cluster, which is
represented using a unique color. The second and third rows of the figure
depict a specific cluster and its corresponding GFA color-coded profiles. We
observe that all three computational models produce plausible clusterings.
From the GFA profiles of the selected cluster (with correspondence across
computational models), we observe that functional varifolds enforce both
geometric _as well as_ along-tract signal profile similarity. Moreover, the
clustering produced with varifolds or MCP (i.e., using only geometric
properties of fibers), are similar to one another and noticeably different
from that of functional varifolds.

##  4 Conclusion

A novel computational model, called functional varifolds, was proposed to
model both geometry and microstructure measure along fibers. We considered the
task of fiber clustering and integrated our functional varifolds model within
framework based on dictionary learning and sparse coding. The driving
hypothesis that combining along-fiber signal with fiber geometry helps
tractography analysis was validated quantitatively and qualitatively using
data from Human Connectome Project. Results show functional varifolds to yield
more consistent clusterings than GFA, varifolds and MCP. While this study
considered a fully unsupervised setting, further investigation would be
required to assess whether functional varifolds augment or aid the
reproducibility of results.

####

Acknowledgements: Data were provided by the Human Connectome Project.

## References

  * [1] Charlier, B., Charon, N., Trouvé, A.: The fshape framework for the variability analysis of functional shapes. Foundations of Computational Mathematics pp. 1–71 (2014)
  * [2] Charon, N., Trouvé, A.: The varifold representation of nonoriented shapes for diffeomorphic registration. SIAM Journal on Imaging Sciences 6(4), 2547–2580 (2013)
  * [3] Colby, J.B., Soderberg, L., Lebel, C., Dinov, I.D., Thompson, P.M., Sowell, E.R.: Along-tract statistics allow for enhanced tractography analysis. Neuroimage 59(4), 3227–3242 (2012)
  * [4] Corouge, I., Gouttard, S., Gerig, G.: Towards a shape model of white matter fiber bundles using diffusion tensor MRI. In: ISBI 2004. pp. 344–347. IEEE (2004)
  * [5] Gori, P., Colliot, O., Marrakchi-Kacem, L., Worbe, Y., Fallani, F.D.V., Chavez, M., Lecomte, S., Poupon, C., Hartmann, A., Ayache, N., et al.: A prototype representation to approximate white matter bundles with weighted currents. In: MICCAI. pp. 289–296. Springer (2014)
  * [6] Hagmann, P., Jonasson, L., Maeder, P., Thiran, J.P., Wedeen, V.J., Meuli, R.: Understanding diffusion mr imaging techniques: From scalar diffusion-weighted imaging to diffusion tensor imaging and beyond 1. Radiographics 26(suppl_1), S205–S223 (2006)
  * [7] Kumar, K., Desrosiers, C.: A sparse coding approach for the efficient representation and segmentation of white matter fibers. In: ISBI, 2016. pp. 915–919. IEEE (2016)
  * [8] Kumar, K., Desrosiers, C., Siddiqi, K.: Brain fiber clustering using non-negative kernelized matching pursuit. In: Machine Learning in Medical Imaging, LNCS, vol. 9352, pp. 144–152 (2015)
  * [9] Kumar, K., Desrosiers, C., Siddiqi, K., Colliot, O., Toews, M.: Fiberprint: A subject fingerprint based on sparse code pooling for white matter fiber analysis. NeuroImage 158, 242 – 259 (2017)
  * [10] Maddah, M., Grimson, W.E.L., Warfield, S.K., Wells, W.M.: A unified framework for clustering and quantitative analysis of white matter fiber tracts. Medical image analysis 12(2), 191–202 (2008)
  * [11] Moberts, B., Vilanova, A., van Wijk, J.J.: Evaluation of fiber clustering methods for diffusion tensor imaging. In: VIS 2005. pp. 65–72. IEEE (2005)
  * [12] O’Donnell, L.J., Westin, C.F., Golby, A.J.: Tract-based morphometry for white matter group analysis. Neuroimage 45(3), 832–844 (2009)
  * [13] Siless, V., Medina, S., Varoquaux, G., Thirion, B.: A comparison of metrics and algorithms for fiber clustering. In: PRNI, 2013. pp. 190–193. IEEE (2013)
  * [14] Van Essen, D.C., Smith, S.M., Barch, D.M., Behrens, T.E., Yacoub, E., Ugurbil, K., Consortium, W.M.H., et al.: The wu-minn human connectome project: an overview. Neuroimage 80, 62–79 (2013)
  * [15] Wang, Q., Yap, P.T., Wu, G., Shen, D.: Application of neuroanatomical features to tractography clustering. Human brain mapping 34(9), 2089–2102 (2013)
  * [16] Wassermann, D., Bloy, L., Kanterakis, E., Verma, R., Deriche, R.: Unsupervised white matter fiber clustering and tract probability map generation: Applications of a Gaussian process framework for white matter fibers. NeuroImage 51(1) (2010)
  * [17] Yeatman, J.D., Dougherty, R.F., Myall, N.J., Wandell, B.A., Feldman, H.M.: Tract profiles of white matter properties: automating fiber-tract quantification. PloS one 7(11), e49790 (2012)
  * [18] Yeh, F.C., Tseng, W.Y.I.: Ntu-90: a high angular resolution brain atlas constructed by q-space diffeomorphic reconstruction. Neuroimage 58(1), 91–99 (2011)