1

Heat Kernel Smoothing in Irregular Image Domains
Moo K. Chung1 , Yanli Wang2 , Gurong Wu3
1

2

Institute of Applied Physics and Computational Mathematics, Beijing, China
Email: wangyanliwyl@gmail.com
3

arXiv:1710.07849v1 [stat.ME] 21 Oct 2017

University of Wisconsin, Madison, USA
Email: mkchung@wisc.edu

University of North Carolina, Chapel Hill, USA
Email: guorong_wu@med.unc.edu
Abstract

We present the discrete version of heat kernel smoothing on graph data structure. The method is used to smooth
data in an irregularly shaped domains in 3D images. New statistical properties are derived. As an application, we
show how to filter out data in the lung blood vessel trees obtained from computed tomography. The method can
be further used in representing the complex vessel trees parametrically and extracting the skeleton representation
of the trees.

I. I NTRODUCTION
Heat kernel smoothing was originally introduced in the context of filtering out surface data defined on
mesh vertices obtained from 3D medical images [1], [2]. The formulation uses the tangent space projection
in approximating the heat kernel by iteratively applying Gaussian kernel with smaller bandwidth. Recently
proposed spectral formulation to heat kernel smoothing [3] constructs the heat kernel analytically using
the eigenfunctions of the Laplace-Beltrami (LB) operator, avoiding the need for the linear approximation
used in [1], [4].
In this paper, we present the discrete version of heat kernel smoothing on graphs. Instead of LaplaceBeltrami operator, graph Laplacian is used to construct the discrete version of heat kernel smoothing.
New statistical properties are derived for kernel smoothing that utilizes the fact heat kernel is a probability
distribution. Heat kernel smoothing is used to smooth out data defined on irregularly shaped domains in
3D images.
The connection between the eigenfunctions of continuous and discrete Laplacians has been well established by several studies [5], [6]. Although many have introduced the discrete version of heat kernel in
computer vision and machine learning, they mainly used the heat kernels to compute shape descriptors
or to define a multi-scale metric [7], [8], [9], [10]. These studies did not use the heat kernels in filtering
out data on graphs. There have been significant developments in kernel methods in the machine learning
community [11], [12], [13], [14], [15]. However, to the best of our knowledge, the heat kernel has never
been used in such frameworks. Most kernel methods in machine learning deal with the linear combination
of kernels as a solution to penalized regressions. On the other hand, our kernel method does not have a
penalized cost function.
As a demonstration, the method is applied in irregularly shaped lung blood vessel trees obtained from
computed tomography (CT). 3D image volumes are represented as a large 3D graph by connecting
neighboring voxels in the vessel trees. Since heat kernel smoothing is analytically represented, it is possible
to use the technique resample lung vessel trees in a higher resolution and obtain the skeleton representation
of the trees for further shape analysis [16], [17].

2

II. P RELIMINARY
Let G = {V, E} be a graph with vertex set V and edge set E. We will simply index the node set as
V = {1, 2, · · · , p}. If two nodes i and j form an edge, we denote it as i ∼ j. Let W = (wij ) be the edge
wight. The adjacency matrix of G is often used as the edge weight. Various forms of graph Laplacian
have been proposed [18] but the most often used standard form L = (lij ) is given by

i∼j
P−wij ,

i=j
lij =
i6=j wij ,
0,
otherwise
The graph Laplacian L can then be written as
L = D − W,
P
where D = (dij ) is the diagonal matrix with dii = nj=1 wij . For this study, we will simply use the
adjacency matrix so that the edge weights wij are either 0 or 1.
Unlike the continuous Laplace-Beltrami operators that may have possibly infinite number of eigenfunctions, we have up to p number of eigenvectors ψ1 , ψ2 , · · · , ψp satisfying
Lψj = λj ψj

(1)

with
0 = λ1 < λ2 ≤ · · · ≤ λp .
The eigenvectors are orthonormal, i.e., ψiT ψj = δij , the Kroneker’s delta. The first eigenvector is trivially
√
given as ψ1 = 1/ p with 1 = (1, 1, · · · , 1)T .
All other higher order eigenvalues and eigenvectors are unknown analytically and have to be computed
numerically. Using the eigenvalues and eigenvectors, the graph Laplacian can be decomposed spectrally.
From (1),
LΨ = ΨΛ,

(2)

where Ψ = [ψ1 , · · · , ψp ] and Λ is the diagonal matrix with entries λ1 , · · · , λp . Since Ψ is an orthogonal
matrix,
p
X
T
T
ΨΨ = Ψ Ψ =
ψj ψjT = Ip ,
j=1

the identify matrix of size p. Then (2) is written as
T

L = ΨΛΨ =

p
X

λj ψj ψjT .

j=1

This is the restatement of the singular value decomposition (SVD) for Laplacian.
For measurement vector f = (f1 , · · · , fp )0 observed at the p nodes, the discrete Fourier series expansion
is given by
n
X
f=
f˜j ψj ,
j=1

where f˜j = f T ψj = ψjT f are Fourier coefficients.

3

Fig. 1.

First few eigenvector of the Laplacian in a L-shaped domain. The first eigenvector ψ0 = 0

III. H EAT KERNEL ON GRAPHS
The discrete heat kernel Kσ is a positive definite symmetric matrix of size p × p given by
Kσ =

p
X

e−λj σ ψj ψjT ,

(3)

j=1

where σ is called the bandwidth of the kernel. Alternately, we can write (3) as
Kσ = Ψe−σΛ ΨT ,
where e−σΛ is the matrix logarithm of Λ. To see positive definiteness of the kernel, for any nonzero
x ∈ Rp ,
T

x Kσ x =

p
X

e−λj σ xT ψj ψjT x

j=1
p

=

X

e−λj σ (ψjT x)2 > 0.

j=1

When σ = 0, K0 = Ip , identity matrix. When σ = ∞, by interchanging the sum and the limit, we obtain
K∞ = ψ1 ψ1T = 11T /p.
K∞ is a degenerate case and the kernel is no longer positive definite. Other than these specific cases, the
heat kernel is not analytically known in arbitrary graphs.
Heat kernel is doubly-stochastic [18] so that
Kσ 1 = 1, 1T Kσ = 1T .
Thus, Kσ is a probability distribution along columns or rows.

4

Fig. 2. Heat kernels with bandwidth σ = 0.01, 0.1. We have used degree 70 expansions but the shape is almost identical if we use higher
degree expansions. The heat kernel is a probability distribution that follows the shape of domain.

Just like the continuous counterpart, the discrete heat kernel is also multiscale and has the scale-space
property. Note
Kσ2

=

p
X

e−(λi +λj )σ ψi ψiT ψj ψjT

i,j=1
p

=

X

e−2λj σ ψj ψjT = K2σ .

j=1

We used the orthonormality of eigenvectors. Subsequently, we have
Kσn = Knσ .
IV. H EAT KERNEL SMOOTHING ON GRAPHS
Discrete heat kernel smoothing of measurement vector f is then defined as
Kσ ∗ f = Kσ f =

p
X

e−λj σ f˜j ψj ,

(4)

j=0

This is the discrete analogue of heat kernel smoothing first defined in [1]. In discrete setting, the
convolution ∗ is simply a matrix multiplication. Then
K0 ∗ f = f
and

K∞ ∗ f = f¯1,

5

Fig. 3. From top left to right: 3D lung vessel tree. Gaussian noise is added to one of the coordinates. 3D graph constructed using 6-connected
neighbors.

P
where f¯ = pj=1 fj /p is the mean of signal f over every nodes. When the bandwidth is zero, we are not
smoothing data. As the bandwidth increases, the smoothed signal converges to the sample mean of all
values.
Define the l-norm of a vector f = (f1 , · · · , fp )T as
k f kl =

p
X

fj


l 1/l

.

j=1

The matrix ∞-norm is defined as
k f k∞ = max fj .
1≤j≤p

Theorem 1. Heat kernel smoothing is a contraction mapping with respect to the l-th norm, i.e.,
kKσ ∗ f kll ≤ kf kll .
Proof. Let kernel matrix Kσ = (kij ). Then we have inequality
kKσ ∗

f kll

=

p
p
X
X
i=1 j=1

l

|kij fj | ≤

p
X

|fj |l .

j=1

We used Jensen’s inequality and doubly-stochastic property of heat kernel.

6

V. S TATISTICAL PROPERTIES
Often observed noisy data f on graphs is smoothed with heat kernel Kσ to increase the signal-to-noise
ratio (SNR) and increases the statistical sensitivity [3]. We are interested in knowing how heat kernel
smoothing will have effects on the statistical properties of smoothed data.
Consider the following addictive noise model:
f = µ + e,

(5)

where µ is unknown signal and  is zero mean noise. Let e = (e1 , · · · , ep )T . Denote E as expectation and
V as covariance. It is natural to assume that the variability of noises at different nodes j is identical, i.e.,
Ee21 = Ee22 = · · · = Ee2p .

(6)

Further, we assume that data at two nodes i and j to have less correlation when the distance between the
nodes is large. So covariance matrix Re = Ve = E(eeT ) = (rij ) can be given by
rij = ρ(dij )

(7)

for some decreasing function ρ and geodesic distance dij between nodes i and j. Note rjj = ρ(0) with
the understanding that djj = 0 for all j. The off diagonal entries of Re are smaller than the diagonal.
Noise e can be further modeled as Gaussian white noise, i.e., Brownian motion or the generalized derivatives of Wiener process, whose covariance matrix elements are Dirac-delta. For the discrete counterpart,
rij = δij , where δij is Kroneker-delta with δij = 1 if i = j and 0 otherwise. Thus,
Re = E(eeT ) = Ip ,
the identity matrix of size p × p. Since δjj ≥ δij , Gaussian white noise is a special case of (7).
Once heat kernel smoothing is applied to (5), we have
Kσ ∗ f = Kσ ∗ µ + Kσ ∗ e.

(8)

We are interested in knowing how the statistical properties of data change from (5) to (8). For Re = Ip ,
the covariance matrix of smoothed noise is simply given as
RKσ ∗e = Kσ E(eeT )Kσ = Kσ2 = K2σ .
We used the scale-space property of heat kernel. In general, the covariance matrix of smoothed data Kσ ∗e
is given by
RKσ ∗e = Kσ E(eeT )Kσ = Kσ Re Kσ .
The variance of data will be often reduced after heat kernel smoothing in the following sense [1], [2].
This is formulated rigorously as follows.
Theorem 2. Heat kernel smoothing reduces variability, i.e.,
V(Kσ ∗ f )j ≤ Vfj
for all j. The subscript
Proof. Note

j

indicates the j-th element of the vector.
p
X
2
V(Kσ ∗ f )j = V(Kσ ∗ e)j = E
kij ei .
i=1

Since (kij ) is doubly-stochastic, after applying Jensen’s inequality, we obtain
p
p
X
2
X

E
kij ei ≤ E
kij e2i = Ee2i .
i=1

i=1

For the last equality, we used the equality of noise variability (6). Since Efj = Ee2i , we proved the
statement.

7

VI. A PPLICATION
We applied heat kernel smoothing to the computed tomography (CT) of human lung obtained from
DIR-lab dataset (https://www.dir-lab.com) [19], [20]. Super-resolution procedure was used to
enhance the resolution from 1 × 1 × 2.5 mm3 to 1 × 1 × 1 mm3 resolution for CT images. The small part
(50 × 50 × 40 voxels) of the CT image was used to illustrate the method and display the results better.
The binary vessel segmentation was done using the multiscale Hessian filters at each voxel, as shown in
the top left of Figure 3 [21], [22], [23].
The binary segmentation was converted into a 3D graph by taking each voxel as a node and connecting
neighboring voxels. Using the 18-connected neighbor scheme, we connect two voxels only if they touch
each other on their faces or edges. If voxels are only touching at their corner vertices, they are not considered as connected. Although we used the 18-connected neighbor scheme in this study, for visualization
purpose only, Figure 3 displays the 6-connected neighbor scheme. This results in an adjacency matrix
and the 3D graph Laplacian. Figure 3 example generates a 3D graph with p = 6615 nodes. The center
of voxel is taken as the node coordinates. The large-scale eigenvector problem was subsequently solved
using an Implicitly Restarted Arnoldi Iteration method [24]. We used 6000 eigenvectors. Note we cannot
have more eigenvectors than the number of nodes. Numbers in Figure 3 are kernel bandwidths. At σ = 0,
heat kernel smoothing is equivalent to Fourier series expansion. Thus, we get almost identical results. As
the bandwidth increases, smoothing converges to the mean value. The regions that are different colors
even with σ = 10000 are regions that are disconnected. Thus, each disconnected regions are converging
to their own different mean values.
VII. D ISCUSSION
The proposed technique can be used to extract the skeleton representation of vessel trees. Here we
show the proof of concept. We perform heat kernel smoothing on node coordinates with σ = 1. Then
rounded off the smoothed coordinates to the nearest integers. The rounded off coordinates were used to
reconstruct the binary segmentation. This gives the thick trees in Figure 4 (top left). To obtain thinner
trees, the smoothed coordinates were scaled by the factor of 2,4 and 6 times before rounding off. This
had the effect of increasing the image size relative to the kernel bandwidth (Figure 4 clockwise from top
right). We hope that our method offers a new way of making the skeleton representation of the vessel
trees [16], [17].
ACKNOWLEDGEMENT
This work was supported by NIH research grants R01 EB022856. We would like to thank Jim Ramsay
of McGill University and Michelle Carey of University College Dublin for valuable discussion on solving
partial differential equations on irregular domains. We would like to thank Ruth Sullivan, Michael Johnson
and Michael Newton of University of Wisconsin-Madison for useful discussion on modeling vessel trees.
R EFERENCES
[1] M.K. Chung, S. Robbins, and A.C. Evans, “Unified statistical approach to cortical thickness analysis,” Information Processing in
Medical Imaging (IPMI), Lecture Notes in Computer Science, vol. 3565, pp. 627–638, 2005.
[2] M.K. Chung, S. Robbins, K.M. Dalton, R.J. Davidson, A.L. Alexander, and A.C. Evans, “Cortical thickness analysis in autism with
heat kernel smoothing,” NeuroImage, vol. 25, pp. 1256–1265, 2005.
[3] M.K. Chung, A. Qiu, S. Seo, and H.K. Vorperian, “Unified heat kernel regression for diffusion, kernel smoothing and wavelets on
manifolds and its application to mandible growth modeling in CT images,” Medical Image Analysis, vol. 22, pp. 63–76, 2015.
[4] X. Han, J. Jovicich, D. Salat, A. van der Kouwe, B. Quinn, S. Czanner, E. Busa, J. Pacheco, M. Albert, R. Killiany, et al., “Reliability
of MRI-derived measurements of human cerebral cortical thickness: The effects of field strength, scanner upgrade and manufacturer,”
NeuroImage, vol. 32, pp. 180–194, 2006.
[5] G.M.L. Gladwell and H. Zhu, “Courant’s nodal line theorem and its discrete counterparts,” The Quarterly Journal of Mechanics and
Applied Mathematics, vol. 55, pp. 1–15, 2002.
[6] T. Tlusty, “A relation between the multiplicity of the second eigenvalue of a graph laplacian, courants nodal line theorem and the
substantial dimension of tight polyhedral surfaces,” Electrnoic Journal of Linear Algebra, vol. 16, pp. 315–24, 2007.

8

Fig. 4. The skeleton representation of vessel trees. Using the heat kernel series expansion with bandwidth σ = 1 and 6000 basis, we
upsampled the binary segmentation at 2, 4, 6 times (clockwise from top right) larger than the original size (top left).

[7] Mikhail Belkin, Partha Niyogi, and Vikas Sindhwani, “Manifold regularization: A geometric framework for learning from labeled and
unlabeled examples,” The Journal of Machine Learning Research, vol. 7, pp. 2399–2434, 2006.
[8] J. Sun, M. Ovsjanikov, and L. J. Guibas, “A concise and provably informative multi-scale signature based on heat diffusion.,” Comput.
Graph. Forum, vol. 28, pp. 1383–1392, 2009.
[9] M. M. Bronstein and I. Kokkinos, “Scale-invariant heat kernel signatures for non-rigid shape recognition,” in IEEE Conference on
Computer Vision and Pattern Recognition (CVPR), 2010, pp. 1704–1711.
[10] F. de Goes, S. Goldenstein, and L. Velho, “A hierarchical segmentation of articulated bodies,” Computer Graphics Forum, vol. 27, pp.
1349–1356, 2008.
[11] Bernhard Schölkopf and Alex J Smola, Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond,
MIT Press, 2002.
[12] J. Nilsson, F. Sha, and M.I. Jordan, “Regression on manifolds using kernel dimension reduction,” in Proceedings of the 24th international
conference on Machine learning. ACM, 2007, pp. 697–704.
[13] J Shawe-Taylor and N Cristianini, Kernel methods for pattern analysis, Cambridge University Press, 2004.
[14] Florian Steinke and Matthias Hein, “Non-parametric regression between manifolds,” Advances in Neural Information Processing

9

Systems, 2008.
[15] F. Yger and A. Rakotomamonjy, “Wavelet kernel learning,” Pattern Recognition, vol. 44, no. 10-11, pp. 2614–2629, 2011.
[16] L. Lindvere, R. Janik, A. Dorr, D. Chartash, B. Sahota, J.G. Sled, and B. Stefanovic, “Cerebral microvascular network geometry
changes in response to functional stimulation,” Neuroimage, vol. 71, pp. 248–259, 2013.
[17] L. Cheng, J. De, X. Zhang, F. Lin, and H. Li, “Tracing retinal blood vessels by matrix-forest theorem of directed graphs,” in
International Conference on Medical Image Computing and Computer Assisted Intervention. Springer, 2014, pp. 626–633.
[18] F.R.K. Chung and S.T. Yau, “Eigenvalue inequalities for graphs and convex subgraphs,” Communications in Analysis and Geometry,
vol. 5, pp. 575–624, 1997.
[19] R. Castillo, E. Castillo, R. Guerra, V.E. Johnson, T. McPhail, A.K. Garg, and T. Guerrero, “A framework for evaluation of deformable
image registration spatial accuracy using large landmark point sets,” Physics in medicine and biology, vol. 54, pp. 1849–1870, 2009.
[20] G. Wu, Q. Wang, J. Lian, and D. Shen, “Estimating the 4d respiratory lung motion by spatiotemporal registration and super-resolution
image reconstruction,” Medical physics, vol. 40, no. 3, pp. 031710, 2013.
[21] A.F. Frangi, W.J. Niessen, K.L. Vincken, and M.A. Viergever, “Multiscale vessel enhancement filtering,” in International Conference
on Medical Image Computing and Computer-Assisted Intervention, 1998, vol. 1496, pp. 130–137.
[22] P.D. Korfiatis, C. Kalogeropoulou, A.N. Karahaliou, A.D. Kazantzi, and L.I. Costaridou, “Vessel tree segmentation in presence of
interstitial lung disease in MDCT,” IEEE Transactions on Information Technology in Biomedicine, vol. 15, pp. 214–220, 2011.
[23] Y. Shang, R. Deklerck, E. Nyssen, A. Markova, J. de Mey, X. Yang, and K. Sun, “Vascular active contour for vessel tree segmentation,”
IEEE Transactions on Biomedical Engineering, vol. 58, pp. 1023–1032, 2011.
[24] R.B. Lehoucq and D.C. Sorensen, “Deflation techniques for an implicitly restarted arnoldi iteration,” SIAM Journal on Matrix Analysis
and Applications, vol. 17, pp. 789–821, 1996.

