Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Med Image Comput Comput Assist Interv. Author manuscript; available in PMC 2010 May 4.
Published in final edited form as:
Med Image Comput Comput Assist Interv. 2009; 12(Pt 2): 18–25.
PMCID: PMC2863144

Volumetric Shape Model for Oriented Tubular Structure from DTI Data


In this paper, we describe methods for constructing shape priors using orientation information to model white matter tracts from magnetic resonance diffusion tensor images (DTI). Shape Normalization is needed for the construction of a shape prior using statistical methods. Moving beyond shape normalization using boundary-only or orientation-only information, our method combines the idea of sweeping and inverse-skeletonization to parameterize 3D volumetric shape, which provides point correspondence and orientations over the whole volume in a continuous fashion. Tangents from this continuous model can be treated as a de-noised reconstruction of the original structural orientation inside a shape. We demonstrate the accuracy of this technique by reconstructing synthetic data and the 3D cingulum tract from brain DTI data and manually drawn 2D contours for each tract. Our output can also serve as the input for subsequent boundary finding or shape analysis.

1 Introduction

Shape priors are powerful tools for image analysis [1,2,3,4,5,6,7,8]. In cases where the signal-to-noise ratio is low, using a prior model, typically generated by manual delineation, can substantially improve the accuracy of the end result. With Magnetic Resonance Diffusion Tensor Imaging (MR-DTI) [9,10,11], see Fig. 1 (left), one can now evaluate directional structural information within a volume in addition to the boundary of a shape. There is a natural demand for incorporating volumetric structural information into prior models for DTI data [12,13,14]. Such models are useful for quantifying and comparing white matter structure in neurologic and psychiatric disorders.

Fig. 1
(Left) Diffusion Tensor Image showing fractional anisotropy (FA) and main diffusion orientations. (Middle) Position and orientation around the cingulum tract within the white matter. (Right) One of the manual segmentations of the cingulum tract.

1.1 Shape Representation

Common shape representations include the boundary mesh [3,8], the distance map [1,4,5,7], and coefficients of harmonic functions [6]. In the context of shape sampling, we can manually specify boundary points and use straight lines in 2D or triangles in 3D to represents the actual boundary of a shape. A distance map can be computed based on the shortest distance from each point to the boundary Such implicit representations have the advantage of representing complex shapes using a single map. On the other hand, a projection from boundary points to harmonic functions transforms the representation from the spatial to a frequency domain. A frequency representation is generally shift and rotation invariant as well as continuous in the spatial domain and therefore easy to resample.

1.2 Point Correspondence in a Volume

Shape normalization [15] establishes a common ground for comparing shapes. Boundary samples on one training shape need to be at the same relative locations on all the others in order to provide a meaningful comparison. If point correspondence cannot be established manually, this could be a point-based registration problem. When using a labeled map, for example the FA map in Fig. 1 (left), instead of boundary points, this could be an intensity-based registration [16]. However, accurate non-rigid intensity-based inter-subject registration can be very difficult for small structures. Assuming intensity maps can be registered reasonably well, if no interior landmark points or features [3] are given, internal point correspondence established over these homogeneous regions through registration is purely controlled by the interpolation between boundaries and smoothing constraints, which are unrelated to the internal structure of the training shape. In cm-reps [17], a shape is parameterized based on an object map and a manual skeleton, which provides interior landmarks, using the idea of inverse-skeletonization. Sun et al. [15] further analyze the cm-rep approach targeted at the corpus callosum from 2D brain DTI.

1.3 Training Shape with Orientation Data

Our goal is to create a set of normalized 3D training shapes with meaningful interior point correspondence. In most cases, specifying internal landmark points over all shapes would be impractical. So, we restrict our input to a limited number of 2D contours, see Fig. 1 (right), which are drawn on several slices of an image volume. Unlike previous surface-driven [15,17] or line-based orientation-driven [12,14] normalization techniques, we represent the training shapes using combinations of continuous basis functions and incorporate orientation data into our normalization process. The analysis of normalized shapes and their underlying orientation patterns can be performed on either the frequency domain or at freely selected spatial locations without the problem of under-sampling or re-meshing.

2 Fourier Surface and Volume in 3D

We begin with a set of 2D manual contours drawn on z-planes, see Fig. 1 (right). We represent each individual contour as a closed curve, with xs), ys) expressed in terms of the arc-length parameter βs [set membership] [0, 2π] and decompose as Fourier series. The Fourier coefficients of the contour at each z level are them-selves decomposed as Fourier series over the range αs [set membership] [0, π]. We construct the Fourier surface [6]:Ssβs) =

[x(αs,βs)y(αs,βs)z(αs,βs)]=k1=0K1k2=0K2[a11(k1,k2)(...)a14(k1,k2)a21(k1,k2)(...)a24(k1,k2)a31(k1,k2)(...)a34(k1,k2)] [cos(k1,αs)cos(k2,βs)cos(k1,αs)sin(k2,βs)sin(k1,αs)cos(k2,βs)sin(k1,αs)sin(k2,βs)]

K1,K2 are the number of harmonics used to represent the long-axis and circumferential variations. To compute the coefficients aij for all frequencies, we perform Fourier Transforms:


Although we define the contours’ αs parameter up to π, the xs, βs) are repeated in reverse z-direction such that we have a full period of samples, see [6] for details. The other aij can be computed similarly. To perform integration in (2) by discrete summations, Ss, βs) is preprocessed to be equally spaced.

After computing the surface coefficients aij (k1, k2) from the input contours, our goal is to solve for the unknown coefficients aij(k1, k2, k3), defined in Eq. (3), which also carry orientation information from the image data D(x, y, z) [set membership] R6. D is the symmetric diffusion tensor computed at each voxel from the DTI data. This computation can be performed by inverse Fourier Transform of V (α, β, γ) :

[x(α,β,γ)y(α,β,γ)z(α,β,γ)]=k1=0K1k2=0K2k3=0K3[a11(k1,k2,k3)(...)a18a21(k1,k2,k3)(...)a28a31(k1,k2,k3)(...)a38] [cos(k1α)cos(k2β)cos(k3γ)cos(k1α)sin(k2β)cos(k3γ)sin(k1α)cos(k2β)cos(k3γ)sin(k1α)sin(k2β)cos(k3γ)cos(k1α)cos(k2β)sin(k3γ)cos(k1α)sin(k2β)sin(k3γ)sin(k1α)cos(k2β)sin(k3γ)sin(k1α)sin(k2β)sin(k3γ)]

whereV(α,β,0)=S(αs,βs)are the boundary conditions

which is an extension of the Fourier surface with an additional parameter γ representing the boundary-to-axis variation of the volume. K3 defines the number of harmonics used for this variation. Generally, parameter α ≠ αs and β ≠ βs. Also, we do not have interior samples of V. In the next section, we will discuss our method of generating ordered samples of V (α, β, γ) from Ss, βs). The input structural orientation data is given by the eigen-decomposition of D and taking the main eigenvector v denoting the orientations of the underlying tissue. We will represent the orientations by the partial derivative of V (α, β, γ) w.r.t. α:

[partial differential][partial differential]αV(α,β,γ)=v(D),   v:R6R3

3 Volumetric Parameterization through the Sweeping Plane

Similar to the challenge in white matter fiber tracking [11,12,13], v(D) contains noise and irrelevant orientations belonging to other structures and thus we cannot normalize the shape solely depending on the data orientations. Our method combines the idea of sweeping [18,13] and inverse-skeletonization [15] to generate ordered samples of V (α, β, γ) which satisfies (4) as well as (5). Without specifying a skeleton, we iteratively track the tract’s axis V (α, β, π) fromthe orientation data. We use a sweeping plane to assist the estimation to reduce local tracking error by taking boundary information into account. With the estimated axis position and a 2D boundary extracted from the input surface S(α, β), we find the inverse-skeleton on the sweeping plane by solving a minimum energy front propagation equation using a 2D level set method [19], see Fig. 2 (right), as described below:

Fig. 2
(Left) A cross section of the parameterized shape with the sweeping plane. (Right) An example of minimum energy 2D front propagation on the sweeping plane.

Assume the first manual contour is the same in our final volume, i.e. V (0, β, 0) = S(0, βs). Let N0, β, γ) be the inward normal of the contour. We solve for V0, β, γ), the area enclosed by V (0, β, 0), according to:

[partial differential][partial differential]γV(α0,β,γ)=(1εK(β,γ))N(α0,β,γ),   β[set membership][0,2π],γ[set membership][0,π]

Here, K(β, γ) is the curvature of the evolving front and ε is a small constant. Let V0, β, γ) be the first sweeping plane, see Fig. 2 (left). Let x¯ (α0), y¯ (α0), z¯ (α0)) = V0, β, π) be the estimated tract’s axis location relative to the first contour. We then compute the average orientation f′ (α0) from D over a neighborhood Ω of the most recent axis point defined by the maximum inscribed disc from that point:

f(α0)=[triple integral](x,y,z)[set membership]ΩWα(xx¯(α0),yy¯(α0),zz¯(α0))v(D(x,y,z))dxdydz

with weighting function and maximum axis point-to-boundary distance:


d(α)=max distα(V(α,β,0),(x¯(α),y¯(α),z¯(α))

Now we move to the next α value by advancing Δα step in the parametric space and Δs step in the image space along the average orientation:



Before going back to (6), we need V1, β, 0) which is not S1, βs) in general except when α = α0. However, we know that V1, β, 0) must be on Ss, βs). So we reconstruct V1, β, 0), which will later become part of the reconstructed surface, from parameterized data points by the intersection between the sweeping plane passing through V1, β, π) with the normal f′(α1) and the input surface Ss, βs). That is, we solve by gradient descent on a set of surface parameters {(αs, βs)} for:


After looping through Equations (6)–(12) for all α [set membership] [0, π], we obtain a set of P sweeping planes with nested level contours on each plane. P will increase if Δs decreases. From the nested contours on each plane, we sample K3 number of contours with uniform separation. The first contour on each plane must be the boundary and the last one must be the tract axis. Each level contour is also equally sampled, having 2K2 number of points, with the first point defined as the closest point to the previous contour. The first boundary point is defined to be the closest point to the first boundary point on the previous plane. In order to measure the distance of points between two sweeping planes, we apply a rigid transformation to align their tract centers and plane orientations. At the end, we have P × 2K2 × K3 number of ordered points of V (α, β, γ). We then project our sampling points V (α, β, γ) into Fourier space defined by (3) and take derivatives to obtain the orientation vectors.

4 Experimental Results

As a validation, we synthesize a tract with orientations by sampling the tangents of a portion of a parametric conical spiral [20]. K1 = 10 is the number of input contours. K2 = 41 and K3 = 6 are user-defined. K1 and K2 are the same in Ss, βs) so that V (α, β, γ) reconstructs Ss, βs) with very small numerical errors. Table 1 compares results on a synthetic image. Input contours are obtained from numerical intersections of vertical planes and the spirals at different locations. The table shows the range of results from 5 sets of contours, and there are 10 contours in each set. Synthetic tensors are formed from the spiral tangents with additional gaussian noise variance from 0.01 to 0.2. With orientation noise variance below 0.1, the basis function approach with sweeping performs the best in terms of average and maximum error as well as error variance. Using basis functions has a smoothing effect which slightly increases the minimum error. However, finite differences are sensitive to individual inaccurate input contours. When orientation data becomes unreliable, the average error of the first method will be similar to the finite difference schemes; however, it is still better in terms of maximum error and error variance. We also applied this method to a set of 18 normal human MR-DTI images with manually traced cingula. Fig. 3 shows our reconstructed orientations of the cingulum tract from real DTI data and visually compares to the original main diffusion direction. Our process filtered out unaligned diffusion vectors, mostly near the tract boundary, which are likely due to partial volume effects or were corrupted by noise. Additional reconstructed surfaces and orientations are given in Fig. 4. Our method does not change the input surface Ss, βs), but modifies its parameterization to become V (α, β, 0).

Fig. 3
Unaligned orientations (lines in the left figure) in a cingulum tract from DTI data and the reconstructed shape (transparent gray volume) and orientations (lines in the right figure)
Fig. 4
Examples of reconstructed tract surfaces V (α, β, 0) (top row) and the corresponding orientations [partial differential][partial differential]αV(α,β,γ) (bottom row)
Table 1
Reconstruction errors for the synthetic tract orientations. Example artificial contours shown on the right. Errors are the angles in degrees between the reconstructed tangents and ground truth orientations. Four methods compared: (1) our method, (2) fit ...

5 Conclusion

We demonstrated a technique to parameterize tracts with boundary and volumetric orientations. Using volumetric basis functions makes the computation of the derivative simple, flexible and precise, which is important for brain tractography [18] based on a prior model derived from these training shapes. Future work includes improving the manual tract surface Ss, βs) by simultaneously sweeping the FA data.


1. Chan T, Zhu W. Level set based shape prior segmentation. IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR); 2005. pp. 1164–1170.
2. Cootes TF, Edwards GJ, Taylor CJ. Active appearance models. IEEE Transactions on Pattern Analysis and Machine Intelligence. 2001;23(6):681–685.
3. Cootes TF, Taylor CJ, Cooper DH, Graham J. Active shape models-their training and application. Computer Vision and Image Understanding. 1995;61(1):38–59.
4. Huang X, Metaxas DN. Metamorphs: Deformable shape and appearance models. IEEE Transactions on Pattern Analysis and Machine Intelligence. 2008;30(8):1444–1459. [PubMed]
5. Leventon ME, Grimson WEL, Faugeras O. Statistical shape influence in geodesic active contours. IEEE Conference on Computer Vision and Pattern Recognition, Proceedings; 2000. pp. 316–323.
6. Staib LH, Duncan JS. Model-based deformable surface finding for medical images. IEEE Transactions on Medical Imaging. 1996;15(5):720–731. [PubMed]
7. Tsai A, Yezzi A, Wells W, Tempany C, Tucker D, Fan A, Grimson WE, Willsky A. A shape-based approach to the segmentation of medical imagery using level sets. IEEE Transactions on Medical Imaging. 2003;22(2):137–154. [PubMed]
8. Wang Y, Staib L. Boundary finding with prior shape and smoothness models. IEEE Transactions on Pattern Analysis and Machine Intelligence. 2000;22(7):738–743.
9. Basser PJ, Jones DK. Diffusion-tensor MRI: theory, experimental design and data analysis - a technical review. NMR Biomed. 2002;15(7–8):456–467. [PubMed]
10. Jellison BJ, Field AS, Medow J, Lazar M, Salamat SM, Alexander AL. Diffusion tensor imaging of cerebral white matter: A pictorial review of physics, fiber tract anatomy, and tumor imaging patterns. AJNR Am. J. Neuroradiol. 2004;25(3):356–369. [PubMed]
11. Mori S, van Zijl PC. Fiber tracking: principles and strategies - a technical review. NMR Biomed. 2002;15(7–8):468–480. [PubMed]
12. O’Donnell L, Westin CL. White matter tract clustering and correspondence in populations. Med. Image Comput. Comput. Assist. Interv. MICCAI. 2005 October [PubMed]
13. Parker GJM, Wheeler-Kingshott CAM, Barker GJ. Estimating distributed anatomical connectivity using fast marching methods and diffusion tensor imaging. IEEE Transactions on Medical Imaging. 2002;21(5):505–512. [PubMed]
14. Wakana S, Jiang H, Nagae-Poetscher LM, van Zijl PC, Mori S. Fiber tract-based atlas of human white matter anatomy. Radiology. 2004;230(1):77–87. [PubMed]
15. Sun H, Yushkevich P, Zhang H, Cook P, Duda J, Simon T, Gee J. Shape-based normalization of the corpus callosum for DTI connectivity analysis. IEEE Transactions on Medical Imaging. 2007;26(9):1166–1178. [PubMed]
16. Rueckert D, Sonoda LI, Hayes C, Hill DL, Leach MO, Hawkes DJ. Non-rigid registration using free-form deformations: application to breast mr images. IEEE Trans. Med. Imaging. 1999;18(8):712–721. [PubMed]
17. Yushkevich P, Zhang H, Gee J. Continuous medial representation for anatomical structures. IEEE Trans. on Medical Imaging. 2006;25(12):1547–1564. [PubMed]
18. Jackowski M, Kao CY, Qiu M, Constable TR, Staib LH. White matter tractography by anisotropic wavefront evolution and diffusion tensor imaging. Medical Image Analysis. 2005;9(5):427–440. [PMC free article] [PubMed]
19. Osher S, Sethian JA. Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations. J. of Comp. Physics. 1988;79:12–49.
20. von Seggern DH. CRC Standard Curves and Surfaces: A Mathematica Notebook. Boca Raton: CRC Press;