Nonparametric Estimation of Cortical Thickness via ... - Semantic Scholar
10-th Annual Meeting of the Organization for Human Brain Mapping 2004
Nonparametric Estimation of Cortical Thickness via Heat Kernel Smoothing Moo K. Chung1,2,3, Shijie Tang1 1 Department of Statistics, 2Department of Biostatistics and Medical Informatics 3 Keck Laboratory for Functional Brain Imaging and Behavior Waisman Center, University of Wisconsin, Madison, WI 53706 June 16, 2004
Introduction
Since kernel Kσ is not a probability distribution on manifolds, we normalize kernel in a small geodesic ball Bp = {q ∈ ∂Ω : d(p, q) ≤ r} ⊂ ∂Ω : d2(p,q) exp − 2σ2 1Bp (q) e σ (p, q) = R (3) K d2(p,q) dµ(q) Bp exp − 2σ 2
For triangular mesh that represents the cortical thickness, it is natural to take a discrete measure µ in defining convolution. Let q1, · · · , qm be neighboring vertices of p = q0 (Figure 3). Then the discrete version of heat kernel is given by
The cerebral cortex has the topology of a 2D highly convoluted sheet. The cortical thickness has been used to characterize the brain shape [2]. The thickness measurements are always contaminated with kp−qik2 exp − 2σ2 noise. The noise may come from a scanner or due fσ (p, qi) = W kp−qj k2 Pm where indicator function 1 is defined as 1 (q) = to the partial volume effect. In order to increase Bp Bp j=0 exp − 2σ 2 the signal-to-noise ratio and smoothness of mea- 1 if q ∈ Bp and 1Bp (q) = 0 otherwise. surements on data defined on the cortex, diffusion and discrete convolution smoothing has been usually used before [1, 2]. We m X present a completely new smoothing technique that fσ (p, qi)Y (qi). fσ ∗ Y (p) = W W is simpler than diffusion smoothing for the estimai=0 tion of the cortical thickness.
This is the generalization of Nadaraya-Watson estimator [4] defined in Euclidean space to manifolds.
Results Heat kernel smoothing with large bandwidth σ = 5 mm is performed iteratively with smaller bandwidth Figure 2: Cortical thickness computed at the Figure 1: Left: original cortical thickness of the posterior right hemisphere of the autistic brain and σ = 0.5 mm. Figure 1 and 2 shows before and after right hemisphere of an autistic subject. Right: heat its iterated heat kernel smoothing with σ = 0.5 and heat kernel smoothing. Heat kernel smoothing dekernel smoothing with σ = 0.5 and k = 100 crease the variability (Theorem 3) while increasing k = 100 iterations. iterations. We present a couple of selected properties of a heat the Gaussianess (Figure 4) which would be useful kernel smoother. Other properties can be found in for random field based multiple comparison. [3]. Methods Theorem 1 Kσ ∗ Y is the unique solution of the folThe cortical surfaces are segmented from T1 lowing initial value problem after time t = σ 2/2: weighted magnetic resonance images using a de∂f formable surface algorithm (FreeSurfer) [5]. The = ∆f, f (p, 0) = Y (p) (4) ∂t cortical thickness is computed as the minimum distance between the outer and inner cortical surfaces. where ∆ is the Laplace-Beltrami operator. The thickness measurement Y is assumed to follow Theorem 2 Heat kernel smoother minimizes the sum Figure 4: QQ-plot of the original cortical thickness measure. Right: QQ-plot after heat kernel the additive model of true signal θ plus noise ǫ on the of weighted squared residuals Z smoothing with σ = 5 mm. cortex ∂Ω: 2 dµ(q). K (p, q) Y (q) − θ] σ Y (p) = θ(p) + ǫ(p), p ∈ ∂Ω. (1) ∂Ω QQ Plot of Sample Data versus Standard Normal
7
8
6
6
5
Quantiles of Input Sample
Quantiles of Input Sample
QQ Plot of Sample Data versus Standard Normal
10
4
2
0
−2 −5
4
3
2
−4
−3
−2
−1 0 1 Standard Normal Quantiles
2
3
4
5
1 −5
−4
−3
−2
−1 0 1 Standard Normal Quantiles
2
3
4
5
Theorem 3 If the covariance function of Y in (1) is Acknowledgements decreasing isotropic function of the form RY (p, q) = The authors wish to thank Bruce Fischl and Yasunari ρ(d(p, q)), then Tosa of the Harvard University for valuable advice Var[Kσ ∗ Y (p)] ≤ VarY (p) for each p ∈ ∂Ω. on segmentation using FreeSurfer. Hence heat kernel smoothing will reduce the variability of cortical thickness measurements. References Theorem 4 Heat kernel smoothing with large bandwidth can be decomposed into multiple kernel [1] Andrade, A., Kherif, F., Mangin, J., Worsley, K.J., Par∂Ω where µ(q) is a surface measure and heat kernel Kσ smoothing with smaller bandwidth via adis, A., Simon, O., Dehaene, S., Le Bihan, D., and Poline, J-B. 2001. Human Brain Mapping. 12:79–93. is defined in terms of the eigenvalues of the Laplace√ K · · ∗ Kσ} ∗f = K kσ ∗ f. σ ∗ ·{z | Beltrami operator. Based on the parametrix expan[2] Chung, M.K., Worsley, K.J., Robbins, S., Paus, P., Taylor, k times J., Giedd, J.N., Rapoport, J.L., Evans, A.C. 2003. Neusion [6], roImage. 18:198-213. h i 2 d (p, q) 1 exp − u0(p, q) Kσ (p, q) ≈ [3] Chung, M.K., Robbins, S., Dalton, K.M., Davidson, R.J. 2 1/2 2σ (2πσ) Evans, A.C. 2004. Cortical thickness Analysis in Autism.
We assume ǫ to be a zero mean Gaussian random field. To overcome the complexity of solving diffusion equations to estimate θ on the cortex [2], we have developed a much simpler method based on the heat kernel smoothing which generalizes Gaussian kernel smoothing in Euclidean space. We define heat kernel smoothing estimator of θ to be the convolution Z b = Kσ ∗ Y (p) = θ(p) Kσ (p, q)Y (q) dµ(q) (2)
where d(p, q) is the geodesic distance between x and y. The first term u0(p, q) → 1 as p → q. When the manifolds is flat, u0(p, q) = 1 and d(p, q) = kp − qk, the Euclidean distance between p and q so the heat kernel Kσ becomes Gaussian kernel h kp − qk2 i 1 Gσ (p, q) = exp − . 2 1/2 2σ (2πσ)
Technical Report, Department of Statistics, University of Wisconsin-Madison. [4] Fan, J., and Gijbels, I. 1996. Local Polynomial Modelling and Its Applications. Chapman & Hall/CRC. [5] Fischl, B. and Dale, A.M. 2000. PNAS 97:11050-11055.
Figure 3: Typical triangular mesh with m = 6 neighboring vertices around p = q0.
[6] Rosenberg, S. 1997. The Laplacian on a Riemannian Manifold, Cambridge University Press.
Nonparametric Estimation of Cortical Thickness via Heat Kernel Smoothing Moo K. Chung1,2,3, Shijie Tang1 1 Department of Statistics, 2Department of Biostatistics and Medical Informatics 3 Keck Laboratory for Functional Brain Imaging and Behavior Waisman Center, University of Wisconsin, Madison, WI 53706 June 16, 2004
Introduction
Since kernel Kσ is not a probability distribution on manifolds, we normalize kernel in a small geodesic ball Bp = {q ∈ ∂Ω : d(p, q) ≤ r} ⊂ ∂Ω : d2(p,q) exp − 2σ2 1Bp (q) e σ (p, q) = R (3) K d2(p,q) dµ(q) Bp exp − 2σ 2
For triangular mesh that represents the cortical thickness, it is natural to take a discrete measure µ in defining convolution. Let q1, · · · , qm be neighboring vertices of p = q0 (Figure 3). Then the discrete version of heat kernel is given by
The cerebral cortex has the topology of a 2D highly convoluted sheet. The cortical thickness has been used to characterize the brain shape [2]. The thickness measurements are always contaminated with kp−qik2 exp − 2σ2 noise. The noise may come from a scanner or due fσ (p, qi) = W kp−qj k2 Pm where indicator function 1 is defined as 1 (q) = to the partial volume effect. In order to increase Bp Bp j=0 exp − 2σ 2 the signal-to-noise ratio and smoothness of mea- 1 if q ∈ Bp and 1Bp (q) = 0 otherwise. surements on data defined on the cortex, diffusion and discrete convolution smoothing has been usually used before [1, 2]. We m X present a completely new smoothing technique that fσ (p, qi)Y (qi). fσ ∗ Y (p) = W W is simpler than diffusion smoothing for the estimai=0 tion of the cortical thickness.
This is the generalization of Nadaraya-Watson estimator [4] defined in Euclidean space to manifolds.
Results Heat kernel smoothing with large bandwidth σ = 5 mm is performed iteratively with smaller bandwidth Figure 2: Cortical thickness computed at the Figure 1: Left: original cortical thickness of the posterior right hemisphere of the autistic brain and σ = 0.5 mm. Figure 1 and 2 shows before and after right hemisphere of an autistic subject. Right: heat its iterated heat kernel smoothing with σ = 0.5 and heat kernel smoothing. Heat kernel smoothing dekernel smoothing with σ = 0.5 and k = 100 crease the variability (Theorem 3) while increasing k = 100 iterations. iterations. We present a couple of selected properties of a heat the Gaussianess (Figure 4) which would be useful kernel smoother. Other properties can be found in for random field based multiple comparison. [3]. Methods Theorem 1 Kσ ∗ Y is the unique solution of the folThe cortical surfaces are segmented from T1 lowing initial value problem after time t = σ 2/2: weighted magnetic resonance images using a de∂f formable surface algorithm (FreeSurfer) [5]. The = ∆f, f (p, 0) = Y (p) (4) ∂t cortical thickness is computed as the minimum distance between the outer and inner cortical surfaces. where ∆ is the Laplace-Beltrami operator. The thickness measurement Y is assumed to follow Theorem 2 Heat kernel smoother minimizes the sum Figure 4: QQ-plot of the original cortical thickness measure. Right: QQ-plot after heat kernel the additive model of true signal θ plus noise ǫ on the of weighted squared residuals Z smoothing with σ = 5 mm. cortex ∂Ω: 2 dµ(q). K (p, q) Y (q) − θ] σ Y (p) = θ(p) + ǫ(p), p ∈ ∂Ω. (1) ∂Ω QQ Plot of Sample Data versus Standard Normal
7
8
6
6
5
Quantiles of Input Sample
Quantiles of Input Sample
QQ Plot of Sample Data versus Standard Normal
10
4
2
0
−2 −5
4
3
2
−4
−3
−2
−1 0 1 Standard Normal Quantiles
2
3
4
5
1 −5
−4
−3
−2
−1 0 1 Standard Normal Quantiles
2
3
4
5
Theorem 3 If the covariance function of Y in (1) is Acknowledgements decreasing isotropic function of the form RY (p, q) = The authors wish to thank Bruce Fischl and Yasunari ρ(d(p, q)), then Tosa of the Harvard University for valuable advice Var[Kσ ∗ Y (p)] ≤ VarY (p) for each p ∈ ∂Ω. on segmentation using FreeSurfer. Hence heat kernel smoothing will reduce the variability of cortical thickness measurements. References Theorem 4 Heat kernel smoothing with large bandwidth can be decomposed into multiple kernel [1] Andrade, A., Kherif, F., Mangin, J., Worsley, K.J., Par∂Ω where µ(q) is a surface measure and heat kernel Kσ smoothing with smaller bandwidth via adis, A., Simon, O., Dehaene, S., Le Bihan, D., and Poline, J-B. 2001. Human Brain Mapping. 12:79–93. is defined in terms of the eigenvalues of the Laplace√ K · · ∗ Kσ} ∗f = K kσ ∗ f. σ ∗ ·{z | Beltrami operator. Based on the parametrix expan[2] Chung, M.K., Worsley, K.J., Robbins, S., Paus, P., Taylor, k times J., Giedd, J.N., Rapoport, J.L., Evans, A.C. 2003. Neusion [6], roImage. 18:198-213. h i 2 d (p, q) 1 exp − u0(p, q) Kσ (p, q) ≈ [3] Chung, M.K., Robbins, S., Dalton, K.M., Davidson, R.J. 2 1/2 2σ (2πσ) Evans, A.C. 2004. Cortical thickness Analysis in Autism.
We assume ǫ to be a zero mean Gaussian random field. To overcome the complexity of solving diffusion equations to estimate θ on the cortex [2], we have developed a much simpler method based on the heat kernel smoothing which generalizes Gaussian kernel smoothing in Euclidean space. We define heat kernel smoothing estimator of θ to be the convolution Z b = Kσ ∗ Y (p) = θ(p) Kσ (p, q)Y (q) dµ(q) (2)
where d(p, q) is the geodesic distance between x and y. The first term u0(p, q) → 1 as p → q. When the manifolds is flat, u0(p, q) = 1 and d(p, q) = kp − qk, the Euclidean distance between p and q so the heat kernel Kσ becomes Gaussian kernel h kp − qk2 i 1 Gσ (p, q) = exp − . 2 1/2 2σ (2πσ)
Technical Report, Department of Statistics, University of Wisconsin-Madison. [4] Fan, J., and Gijbels, I. 1996. Local Polynomial Modelling and Its Applications. Chapman & Hall/CRC. [5] Fischl, B. and Dale, A.M. 2000. PNAS 97:11050-11055.
Figure 3: Typical triangular mesh with m = 6 neighboring vertices around p = q0.
[6] Rosenberg, S. 1997. The Laplacian on a Riemannian Manifold, Cambridge University Press.
