Diffusive smoothing of 3D segmented medical data.

This paper proposes an accurate, computationally efficient, and spectrum-free formulation of the heat diffusion smoothing on 3D shapes, represented as triangle meshes. The idea behind our approach is to apply a [Formula: see text]-degree Padé-Chebyshev rational approximation to the solution of the heat diffusion equation. The proposed formulation is equivalent to solve r sparse, symmetric linear systems, is free of user-defined parameters, and is robust to surface discretization. We also discuss a simple criterion to select the time parameter that provides the best compromise between approximation accuracy and smoothness of the solution. Finally, our experiments on anatomical data show that the spectrum-free approach greatly reduces the computational cost and guarantees a higher approximation accuracy than previous work.

Introduction

In medical applications, the heat kernel is central in diffusion filtering and smoothing of images [1-6], 3D shapes [7,8], and anatomical surfaces [9,10]. However, the computational cost for the evaluation of the heat kernel is the main bottleneck for processing both surfaces and volumetric data; in fact, it takes from to time on a data set sampled with n points, according to the sparsity of the Laplacian matrix. This aspect becomes more evident for medical data, which are nowadays acquired by PET, MRI systems and whose resolution is constantly increasing with the improvement of the underlying imaging protocols and hardware.

To overcome the time-consuming computation of the Laplacian spectrum on large data sets (Section ‘Previous work’), the heat kernel has been approximated by prolongating its values evaluated on a sub-sampling of the input surface [11-13]; applying multi-resolution decompositions [14] or a rational approximation of the exponential representation of the heat kernel [15]; and considering the contribution of the eigenvectors related to smaller eigenvalues. The heat equation has been solved through explicit [16] or backward [17,18] Euler methods, whose solution no more satisfies the diffusion problem. Further approaches apply a Krylov subspace projection [19], which becomes computationally expensive when the dimension of the Krylov space increases, still remaining much lower than n.

This paper proposes an accurate, computationally efficient, and spectrum-free evaluation of the diffusive smoothing on 3D shapes, represented as polygonal meshes. The idea behind our approach (Section ‘Discrete heat diffusion smoothing’) is to apply the -degree Padé–Chebyshev rational polynomial approximation of the exponential map to the solution of the heat equation. This spectrum-free formulation converts the heat equation to a set of sparse, symmetric linear systems and the resulting computational scheme is independent of the evaluation of the Laplacian spectrum, the selection of a specific subset of eigenpairs, and multi-resolutive prolongation operators. Our approach has a linear computational cost, is free of user-defined parameters, and works with sparse, symmetric, well-conditioned matrices. Since the computation is mainly based on numerical linear algebra, our method can be applied to any class of Laplacian weights and any data representation (e.g., 3D shapes, multi-dimensional data), thus overcoming the ambiguous definition of multi-resolutive and prolongation operators on point-sampled or non-manifold surfaces. Bypassing the computation of the eigenvectors related to small eigenvalues, which are necessary to correctly recover local features of the input shape or signal, the spectrum-free computation is robust with respect to data discretization. As a result, it properly encodes local and global features of the input data in the heat diffusion kernel. For any data representation and Laplacian weights, the accuracy of the heat smoothing computed through the Padé–Chebyshev approximation is lower than , where is the degree of the rational polynomial, and can be further reduced by slightly increasing r. Finally (Section ‘Results and discussion’), our experiments on surfaces and volumes representing anatomical data show that the spectrum-free approach greatly reduces the computational cost (from 32 up to 164 times) and guarantees a higher approximation accuracy than previous work.

Previous work

Let us consider the heat equation , , on a closed, connected manifold of , where defines the initial condition on . The solution to the heat equation , , is computed as the convolution between the initial condition f and the heat kernel . Here, is the Laplacian eigensystem , .

The heat equation is solved through its FEM formulation [20] on a discrete surface (e.g., triangle mesh, point set) of . Indicating with the Laplacian matrix, which discretizes the Laplace–Beltrami operator on , the “power” method applies the identity , where m is chosen in such a way that is sufficiently small to guarantee that the approximation is accurate. Here, is the identity matrix. However, the selection of m and its effect on the approximation accuracy cannot be estimated a-priori. In [17,18], the solution to the heat equation is computed through the Euler backward method , . The resulting functions are over-smoothed and converge to a constant map, as . Krylov subspace projection [19], which replaces the Laplacian matrix with a full coefficient matrix of smaller size, has computational and memory bottlenecks when the dimension k of the Krylov space increases, still remaining much lower than n (e.g., ).

Once the Laplacian matrix has been computed, we evaluate its spectrum and approximate the heat kernel by considering the contribution of the Laplacian eigenvectors related to smaller eigenvalues, which are computed in superlinear time [21]. Such an approximation is accurate only if the exponential filter decays fast (e.g., large values of time). Otherwise, a larger number of eigenpairs is needed and the resulting computational cost varies from to time, according to the sparsity of the Laplacian matrix. Furthermore, the number of eigenpairs is heuristically selected and its effect on the resulting approximation accuracy cannot be estimated without computing the whole spectrum. Finally, we can apply multi-resolution prolongation operators [13] and numerical schemes based on the Padé–Chebyshev polynomial [22,15]. However, previous work has not addressed this extension, convergence results, and the selection of the optimal scale.

Discrete heat diffusion smoothing

Let us discretize the input shape as a triangle mesh , with vertices , which is the output of a 3D scanning device or a segmentation of a MRI acquisition of an anatomical structure. Let be the Laplacian matrix, where is a symmetric, positive semi-definite matrix and is a symmetric and positive definite matrix. On triangle meshes, is the Laplacian matrix with cotangent weights [23,24] or associated with the Gaussian kernel [25], and is the mass matrix of the Voronoi [18] or triangle [26] areas. For any class of weights, the Laplacian matrix is uniquely defined by the couple and is associated to the generalized eigensystem such thatwhere and Λ are the eigenvectors’ and eigenvalues’ matrices. From the relation (1), we get the identities andThen, the spectral representation of the heat kernel isFor a signal , , sampled at , the solution , , to the heat equation , , is achieved by multiplying the heat kernel matrix with the initial condition . Applying the Padé–Chebyshev approximation to the exponential of the Laplacian matrix in Eq. (3), we getand the vector is the sum of the solutions of r sparse linear systemsWe briefly recall that the weights and nodes of the Padé–Chebyshev approximation (4) are precomputed for any polynomial degree [27]. Each vector is calculated as a minimum norm residual solution [28], without pre-factorizing the matrices and . Algorithm 1 summarizes the main steps of the proposed computation.

Spectrum-free heat kernel smoothing.

According to Varga [29], the approximation error between the exponential map and its rational polynomial approximation is bounded by the uniform rational Chebyshev constant , which is independent of t and lower than . Assuming exact arithmetic, the approximation error is bounded asin particular, selecting the degree in Eq. (6) provides an error lower than , which is satisfactory for the approximation of on 3D shapes. Iterative solvers of sparse linear systems are generally efficient and accurate for the computation of the diffusion smoothing; for several values of t, a factorization (e.g., LU) of the coefficient matrix of the linear systems can be precomputed and used for their solution in linear time.

Optimal time parameter

Among the possible time parameters, we select a value that provides a small residual and a low value of the penalty term , which controls the smoothness of the solution. Rewriting these two functions in terms of the Laplacian spectrum asthe residual and penalty terms are increasing and decreasing maps with respect to t, respectively. If t tends to zero, then the residual becomes null and the smoothness term converges to the energy . If t becomes large, then the residual tends to and the solution norm converges to . Indeed, the plot of is L-shaped [30] and its corner provides the optimal regularization parameter, which is the best compromise between approximation accuracy and smoothness (Fig. 1a).

Fig. 1

L-curve and discrepancy. (a) Optimal parameter and corresponding diffusion smoothing (upper part, rightPadé–Chebyshev approximation of degree ) on the noisy 3D shapes of the teeth. (b) Error (y-axis) between the Padé–Chebyshev approximation () of and its truncated spectral approximation with k eigenpairs (, x-axis), and different values of t.

In previous work, the evaluation of the L-curve is computationally expensive, as it generally involves the evaluation of the Laplacian spectrum and/or the solution of a linear system with slowly converging iterative solvers. Through the Padé–Chebyshev approximation, we have an efficient way to evaluate the map for several values of t, thus precisely estimating the optimal time parameter. In fact, the terms in Eq. (7) are evaluated by applying the Padé–Chebyshev approximation of and computing and . In this way, we avoid the evaluation of the spectral representations (7) through the computation of the Laplacian spectrum.

Results and discussion

We consider the solution to the heat diffusion process, whose initial condition takes value 1 at the anchor point and 0 otherwise. For our tests on triangle meshes, we have selected the linear FEM weights [26,21]. In this case [15], the discretization of the inner product is induced by the matrix , which is intrinsic to the surface and is adapted to the local sampling through the variation of the triangles’ or Voronoi areas. In the paper examples, the level-sets are associated with iso-values uniformly sampled in the range of the solution to the heat equation, whose minimum and maximum are depicted in blue and red, respectively. Furthermore, the color coding represents the same scale of values for multiple shapes. Noisy examples have been achieved by adding a 20% Gaussian noise to the input shapes.

Truncated spectral and Padé–Chebyshev approximations

For the truncated spectral approximation of the solution to the heat equation, the number k of eigenpairs must be selected by the user and the approximation accuracy cannot be estimated without extracting the whole spectrum. The different accuracy (Fig. 1b) of the truncated spectral approximation and the Padé–Chebyshev method of the heat kernel is analyzed by measuring the approximation error (y-axis) between the spectral representation of the heat kernel , computed using a different number k (x-axis) of eigenfunctions, and the corresponding Padé–Chebyshev approximation. For small values of t, the partial spectral representation requires a large number k of Laplacian eigenvectors to recover local details. For instance, selecting eigenpairs the approximation error remains higher than ; in fact, local shape features encoded by are recovered for a small t using the eigenvectors associated with high frequencies, thus requiring the computation of a large part of the Laplacian spectrum. For large values of t, increasing k strongly reduces the approximation error until it becomes almost constant and close to zero. In this case, the behavior of the heat kernel is mainly influenced by the Laplacian eigenvectors related to the eigenvalues of smaller magnitude. Indeed, the spectral representation generally requires a high number of eigenpairs without achieving an accuracy of the same order of the spectrum-free approximation through the Padé–Chebyshev method.

Robustness to noise and sampling

Figs. 2–4 compare the diffusion smoothing of a noisy data set computed with the Padé–Chebyshev approximation of degree and the truncated approximation with k Laplacian eigenpairs. A low number of eigenpairs does not preserve shape details; increasing k reconstructs the surface noise. For both examples, the error between (a) and the smooth approximation of (b) is lower than 1% for (c) the Padé–Chebyshev method and (d) varies from 12% () up to 13% for the truncated spectral approximation.

On irregularly-sampled and noisy shapes (Figs. 5 and 6), the spectrum-free computation provides smooth level sets, which are well-distributed around the anchor point and remain almost unchanged and coherent with respect to the original shape. A higher resolution of improves the quality of the level-sets of the canonical basis function, which are always uniformly distributed around the anchor (black dot). Finally, an increase of the noise magnitude does not affect the shape and distribution of the level sets.

We also compare the accuracy of the heat kernel on the unitary sphere and computed with (i) the proposed approach; (ii) the spectral representation of the heat kernel , with k eigenpairs; (iii) the Euler backward method; and (iv) the power method (Section ‘Previous work’). For all the scales (Fig. 7), the approximation accuracy of the Padé–Chebyshev method is higher than the truncated Laplacian spectrum with k eigenpairs, , the Euler backward method, and the power method. Reducing the scale, the accuracy of the Padé–Chebyshev remains almost unchanged while the other methods are affected by a larger discrepancy and tend to have an analogous behavior . Finally, the Euler backward method generally over-smooths the solution, which converges to a constant map as , and the selection of m with respect to the shape details is guided by heuristic criteria.

Fig. 7

Comparison of the accuracy of different approximations of the heat kernel on the unitary sphere. error (y-axis) between the ground-truth diffusion smoothing on the cylinder, with a different sampling (x-axis) and scales. For different scales, the accuracy of the Padé–Chebyshev method (, orange line) remains almost unchanged and higher than the truncated approximation with 100 and 200 eigenpairs (red, blue line), the Euler backward (green line) and power (black line) methods.

Numerical stability

According to Section ‘Discrete heat diffusion smoothing’, the scale t influences the conditioning number of the coefficient matrices , , which are generally well-conditioned, as also confirmed by our experiments (Fig. 8). While previous work requires to heuristically tune the number of selected eigenpairs to the chosen scale, the Padé–Chebyshev approximation has a higher approximation accuracy, which is independent of the selected scale. Furthermore, those scales close to zero would require a larger number of eigenpairs, thus resulting in a larger computational cost for the truncated spectral approximation.

Fig. 8

Numerical stability of the Padé–Chebyshev approximation. With reference to Fig. 4, conditioning number (y-axis) of the matrices , for different time parameters t; the indices of the coefficients are reported on the x-axis.

Computational cost

Approximating the exponential map with a (rational) polynomial of degree r, the evaluation of the solution to the heat diffusion equation and the evaluation of the heat kernel at is reduced to solve r sparse, symmetric, linear systems (c.f., Eq. (5)), whose coefficient matrices have the same structure and sparsity of the adjacency matrix of the input triangle mesh. Applying an iterative and sparse linear solver (e.g., Gauss–Seidel method, conjugate gradient) [28] (Ch. 10), the computational cost for the evaluation of the heat kernel and the diffusion distance between two points is , where is the computational cost of the selected solver. Here, the function , which depends on the number n of shape samples and the sparsity of the coefficient matrix, typically varies from to . In fact, is the average computational cost of the aforementioned iterative solvers of sparse linear systems. Timings (Table 1) are also reduced from 32 up to 164 times with respect to the approximation based on a fixed number of Laplacian eigenpairs.

Table 1

Timings (in seconds) for the evaluation of the heat kernel on 3D shapes with n points, approximated with eigenpairs (Eigs) and the Padé–Chebyshev approximation (Cheb.). Column ‘’ indicates the number of times the computational cost is reduced. Tests have been performed on a 2.7 GHz Intel Core i7 Processor, with 8 GB memory.

Teeth surf. (Fig. 3)				Brain (Fig. 5)
n (K)	Eigs	Cheb.	×	n (K)	Eigs	Cheb.	×
10	39.01	0.32	122	20	99.77	0.61	164
50	154.13	2.50	62	50	189.02	2.08	91
80	188.21	4.12	46	100	299.20	4.98	60
100	307.03	6.21	49	200	658.11	11.20	59
200	450.21	10.03	45	400	850.11	18.21	47
500	670.31	21.11	32	500	1001.11	32.11	78

Conclusions and future work

We have presented an efficient computation of the diffusion soothing of medical data and the selection of the optimal scale, which provides the best compromise between approximation accuracy and smoothness of the solution. As future work, we foresee a specialization of the spectrum-free computation and the selection of the optimal time parameter for the analysis of brain structures and the smoothing of MRI images.

Conflict of interest

The author declares no conflict of interest.

Compliance with Ethics Requirements

This article does not contain any studies with human or animal subjects.

Require: A noisy map f:P→R, f≔(f(p_i))i=1n.

Ensure: A smooth approximation F(t)=K_tf of f.

1: Select the value of t (e.g., optimal value, Section ‘Discrete heat diffusion smoothing’).

2: fori=1,…,r-1do

3:  Compute g_i:(tL+θ_iB)g_i=-α_iBf.

4: end for

5: Approximate K_tf as α_0f+∑i=1rg_i.