Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2863144

Formats

Article sections

- Abstract
- 1 Introduction
- 2 Fourier Surface and Volume in 3D
- 3 Volumetric Parameterization through the Sweeping Plane
- 4 Experimental Results
- 5 Conclusion
- References

Authors

Related links

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

NIHMSID: NIHMS178047

Hon Pong Ho,^{1} Xenophon Papademetris,^{1,}^{2} Fei Wang,^{3} Hilary P. Blumberg,^{3} and Lawrence H. Staib^{1,}^{2}

See other articles in PMC that cite the published article.

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.

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.

(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.

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.

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.

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.

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 *x*(β* _{s}*),

$$\left[\begin{array}{c}\hfill x({\alpha}_{s},{\beta}_{s})\hfill \\ \hfill y({\alpha}_{s},{\beta}_{s})\hfill \\ \hfill z({\alpha}_{s},{\beta}_{s})\hfill \end{array}\right]={\displaystyle \sum _{{k}_{1}=0}^{{K}_{1}}{\displaystyle \sum _{{k}_{2}=0}^{{K}_{2}}[\begin{array}{c}\hfill {a}_{11}({k}_{1},{k}_{2}){a}_{14}({k}_{1},{k}_{2})\hfill & \hfill {a}_{21}({k}_{1},{k}_{2}){a}_{24}({k}_{1},{k}_{2})\hfill & \hfill {a}_{31}({k}_{1},{k}_{2}){a}_{34}({k}_{1},{k}_{2})\hfill \\ ]\end{array}\text{[cos(k1,\alpha s)cos(k2,\beta s)cos(k1,\alpha s)sin(k2,\beta s)sin(k1,\alpha s)cos(k2,\beta s)sin(k1,\alpha s)sin(k2,\beta s)]}}}$$

(1)

*K*_{1},*K*_{2} are the number of harmonics used to represent the long-axis and circumferential variations. To compute the coefficients *a _{ij}* for all frequencies, we perform Fourier Transforms:

$${a}_{11}({k}_{1},{k}_{2})={\displaystyle {\int}_{{\alpha}_{s}=0}^{2\pi}{\displaystyle {\int}_{{\beta}_{s}=0}^{2\pi}x({\alpha}_{s},{\beta}_{s})\phantom{\rule{thinmathspace}{0ex}}\text{cos}({k}_{1},{\alpha}_{s})\phantom{\rule{thinmathspace}{0ex}}\text{cos}({k}_{2},{\beta}_{s})d{\alpha}_{s}d{\beta}_{s}}}$$

(2)

Although we define the contours’ α* _{s}* parameter up to π, the

After computing the surface coefficients *a _{ij}* (

$$\left[\begin{array}{c}\hfill x(\alpha ,\beta ,\gamma )\hfill \\ \hfill y(\alpha ,\beta ,\gamma )\hfill \\ \hfill z(\alpha ,\beta ,\gamma )\hfill \end{array}\right]={\displaystyle \sum _{{k}_{1}=0}^{{K}_{1}}{\displaystyle \sum _{{k}_{2}=0}^{{K}_{2}}{\displaystyle \sum _{{k}_{3}=0}^{{K}_{3}}[\begin{array}{c}\hfill {a}_{11}({k}_{1},{k}_{2},{k}_{3}){a}_{18}\hfill & \hfill {a}_{21}({k}_{1},{k}_{2},{k}_{3}){a}_{28}\hfill & \hfill {a}_{31}({k}_{1},{k}_{2},{k}_{3}){a}_{38}\hfill \\ ]\end{array}}\text{[cos(k1\alpha )cos(k2\beta )cos(k3\gamma )cos(k1\alpha )sin(k2\beta )cos(k3\gamma )sin(k1\alpha )cos(k2\beta )cos(k3\gamma )sin(k1\alpha )sin(k2\beta )cos(k3\gamma )cos(k1\alpha )cos(k2\beta )sin(k3\gamma )cos(k1\alpha )sin(k2\beta )sin(k3\gamma )sin(k1\alpha )cos(k2\beta )sin(k3\gamma )sin(k1\alpha )sin(k2\beta )sin(k3\gamma )]}}}$$

(3)

$$\text{where}\phantom{\rule{thinmathspace}{0ex}}V(\alpha ,\beta ,0)=S({\alpha}_{s},{\beta}_{s})\phantom{\rule{thinmathspace}{0ex}}\text{are the boundary conditions}$$

(4)

which is an extension of the Fourier surface with an additional parameter γ representing the boundary-to-axis variation of the volume. *K*_{3} 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 *S*(α* _{s}*, β

$$\frac{\alpha V(\alpha ,\beta ,\gamma )=\mathbf{v}(D),\text{}\mathbf{v}:{6}^{\to}}{}$$

(5)

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:

(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

$$\frac{\gamma V({\alpha}_{0},\beta ,\gamma )=(1-\epsilon K(\beta ,\gamma ))\phantom{\rule{thinmathspace}{0ex}}N({\alpha}_{0},\beta ,\gamma ),\text{}\beta [0,2\pi ],\gamma [0,\pi ]}{}$$

(6)

Here, *K*(β, γ) is the curvature of the evolving front and ε is a small constant. Let *V* (α_{0}, β, γ) be the first sweeping plane, see Fig. 2 (left). Let *x*¯ (α_{0}), *y*¯ (α_{0}), *z*¯ (α_{0})) = *V* (α_{0}, β, π) 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\prime ({\alpha}_{0})=\phantom{\rule{thinmathspace}{0ex}}{\displaystyle \underset{}{(x,y,z)\Omega \phantom{\rule{thinmathspace}{0ex}}{W}_{\alpha}\phantom{\rule{thinmathspace}{0ex}}(}x-\overline{x}({\alpha}_{0}),y-\overline{y}({\alpha}_{0}),z-\overline{z}({\alpha}_{0}))\phantom{\rule{thinmathspace}{0ex}}\mathbf{v}(D(x,y,z))\mathit{\text{dxdydz}}}$$

(7)

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

$${W}_{\alpha}({d}_{1},{d}_{2},{d}_{3})={e}^{-}\sqrt{{({d}_{1})}^{2}+{({d}_{2})}^{2}+{({d}_{3})}^{2}}/d(\alpha )$$

(8)

$$d(\alpha )=\underset{\alpha}{\mathit{\text{max dist}}}(V(\alpha ,\beta ,0),(\overline{x}(\alpha ),\overline{y}(\alpha ),\overline{z}(\alpha ))$$

(9)

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:

$${\alpha}_{1}={\alpha}_{0}+\mathrm{\Delta}\alpha $$

(10)

$$V({\alpha}_{1},\beta ,\pi )=V({\alpha}_{0},\beta ,\pi )+\Delta sf\prime ({\alpha}_{0})/\left|f\prime ({\alpha}_{0})\right|$$

(11)

Before going back to (6), we need *V* (α_{1}, β, 0) which is not *S*(α_{1}, β* _{s}*) in general except when α = α

$$V({\alpha}_{1},\beta ,0)=\{S({\alpha}_{s},{\beta}_{s})|(S({\alpha}_{s},{\beta}_{s})-V({\alpha}_{1},\beta ,\pi ))\phantom{\rule{thinmathspace}{0ex}}\xb7\phantom{\rule{thinmathspace}{0ex}}f\prime ({\alpha}_{1})=0\}$$

(12)

After looping through Equations (6)–(12) for all α [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 *K*_{3} 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 2*K*_{2} 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* × 2*K*_{2} × *K*_{3} 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.

As a validation, we synthesize a tract with orientations by sampling the tangents of a portion of a parametric conical spiral [20]. *K*_{1} = 10 is the number of input contours. *K*_{2} = 41 and *K*_{3} = 6 are user-defined. *K*_{1} and *K*_{2} are the same in *S*(α* _{s}*, β

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)

Examples of reconstructed tract surfaces *V* (α, β, 0) (top row) and the corresponding orientations $\frac{\alpha V(\alpha ,\beta ,\gamma )}{}$ (bottom row)

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 *S*(α* _{s}*, β

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;

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |