# Sparse Representation of Deformable 3D Organs with Spherical Harmonics and Structured Dictionary

Authors: Dan Wang; Ahmed H. Tewfik; Yingchun Zhang; Yunhe Shen
Journal: International Journal of Biomedical Imaging
 (2011)

---

## Abstract

This paper proposed a novel algorithm to sparsely represent a deformable surface (SRDS) with low dimensionality based on spherical harmonic decomposition (SHD) and orthogonal subspace pursuit (OSP). The key idea in SRDS method is to identify the subspaces from a training data set in the transformed spherical harmonic domain and then cluster each deformation into the best-fit subspace for fast and accurate representation. This algorithm is also generalized into applications of organs with both interior and exterior surfaces. To test the feasibility, we first use the computer models to demonstrate that the proposed approach matches the accuracy of complex mathematical modeling techniques and then both ex vivo and in vivo experiments are conducted using 3D magnetic resonance imaging (MRI) scans for verification in practical settings. All results demonstrated that the proposed algorithm features sparse representation of deformable surfaces with low dimensionality and high accuracy. Specifically, the precision evaluated as maximum error distance between the reconstructed surface and the MRI ground truth is better than 3 mm in real MRI experiments.


---

## Body

## 1. Introduction

Organ deformation during operations has imposed substantial challenges for performing precise diagnosis and surgery in minimum invasive surgery (MIS), such as natural orifice transluminal endoscopic (NOTES) [1]. The deformation markedly decreases the precision of the prior surgical plan that is based on the preoperative images (e.g. computed tomography (CT) or magnetic resonance imaging (MRI)), so it must be effectively compensated to lower surgical risks. However, this is not a trivial task due to the high degree of freedom of the 3D nonrigid deformation and limited field of view [1] for observation during MIS. To recovery, the 3D deformation with high resolution in real time, a critical issue is to seek an efficient representation of deformable surface, according to which the sampling and surface recovery strategy can be designed for updating the 3D visualization. This paper focuses on the topic of block sparse representation of deformable surfaces. The later topic of real time tracking an deformable organ with limited access to the organ is explored further in [2–4].Various techniques have been proposed for surface description, and each has its own advantage and disadvantage according to the application requirements. Broadly speaking, there are two major categories of surface representation methods: local feature-based models and global or parametric models. The work in [5] based on geometric partial differential equations (PDE) belongs to the former category which derives Euler-Lagrange equation and then a geometric evolution equation (or geometric flow) to describe the surfaces. Similarly, the method in [6] treats the whole surface as a union of localized patches. Global surface representation methods [7–10], particularly, decomposing surfaces into other primitive shapes, are more appropriate for shape analysis and classification due to the lower dimensionality of the parameter space. This paper falls into the category of parametric global surface description.Parametric surface representation describes a surface in a single functional form, such that the surface is fully characterized by a set of parameters in a particular domain. Surface harmonics such as spheroidal harmonics, cylindrical harmonics, and spherical harmonics (SH) [7] are widely used as building blocks for global surface description. Each harmonic does not bear localized features but contributes to the entire shape description. Among those different types of harmonics, a well-known approach is the spherical harmonic decomposition (SHD), which has advantages of smoothness and high accuracy [7, 11]. With proper parameterization [11, 12], any genus-0 surface (The genus of a connected surface is an integer representing the maximum number of cuttings along nonintersecting closed simple curves without rendering the resultant manifold disconnected. It is equal to the number of handles on it, so a sphere is genus 0 and a torus is genus 1.) can be analyzed in the harmonic domain with reduced data dimensions. SHD has been widely used in applications related to surface description, including static modeling of kidney [8] and brain [9, 10], as well as spatial-temporal modeling of left-ventricular with known motion period. In [13], hemisphere is also applied to open surfaces.Sparse signal representation has steadily gained attention over the years in the signal processing community. The aim is to find a representation which is sparse, or compact, such that most of the energy of a signal can be captured with only a few nonzero coefficients in a given dictionary. The first widely applied methods to seek sparse representation are greedy approaches, including matching pursuit (MP) [14], orthogonal matching pursuit (OMP) [15], and orthogonal least square (OLS) [16]. Those methods iteratively first select the most correlated element from dictionary and then remove the contribution of that element with decorrelation, before finding the next atom. Iteration terminates when any stopping criteria is met. The second type is global optimization algorithms, such as Basis Pursuit (BP) [17], FOCUSS [18] and Iterative Thresholding [19]. Global Optimization, in the approximate sense, relaxes the sparseness constrain, and its sparsity is a side-effect of the optimization. For example, the basis pursuit (BP) method approximates the l0-norm sparsity constraint with an l1-norm criteria, which effectively converts the problem into a convex optimization one, solved globally with linear programming. The orthogonal subspace pursuit (OSP) method [20] used in our paper belongs to the greedy category, which does not require prior knowledge of the dimension of the subspaces and combines the learned subspaces to produce a data-driven dictionary with good sparseness and generalizability.Besides sparse decomposition algorithms as mentioned above, an equally important issue for sparse representation is how to select a dictionary for an application. The two main groups for dictionary design methods are structured dictionaries built out of common bases, and trained dictionaries that are inferred from the training data. For the common bases, it is well known that the wavelet transform can be used to generate sparse multiscale representations of images, the short-time Fourier transform (STFT) generates sparse time-frequency representation of speech signals, and the DCT is another transform that has been used for compression in audio coding algorithms due to its good compaction property. For dictionary learning, the applicable approaches include maximum likelihood estimation (MLE) [21], method of optimal directions (MOD) [22], maximum a-posteriori (MAP) [23], and so forth. Those methods attempt to generalize the type of considered signal with the basis identified from a representative training data set. The proposed approach is based on trained dictionary, since, to the best of our knowledge, there is no common basis in which random surfaces can be sparsely represented.Although sparse representation has been widely applied in the fields of signal compression, image denoising, blind source separation, and compressed sensing, there is still very limited application in 3D surfaces [24]. In fact, the main statistically motivated surface modeling methods are based on principle component analysis which is not sparse [25–29]. Those methods first compute the mean shape and then build the model by establishing legal variations learned from a set of training data for a given type of images, such as bone [25]. With PCA, the major variations of the shape populations are described by the first few basis vectors, such that any surface of that shape population can be projected into an orthogonal subspace spanned by the retained vectors. More advanced techniques, such as multiresolution deformable model [27], are provided to improve the accuracy considering limited sample size.Most of the previous surface modeling/representation works are designed for either static models [8–10], or particular deformable organs with known physical properties (such as motion cycle) [30]. Further, computation bottleneck caused by large spherical harmonics basis hampered the applicability for real-time application. For PCA-based modeling methods, the resulting space that captures the variation in the population is either a super subspace including all training data or a truncated subspace with sacrificed generalization. In addition, PCA tends to be computational expensive when performing eigenvector decomposition as training data dimension increases, and it does not lead to any structure in the representation.To bring the demonstrated merits of sparse coding to 3D surface representation, we propose a generally applicable algorithm of parametric sparse representation of deformable surfaces based on SHD and orthogonal subspace pursuit. The main contributions of this paper include the following.(i)
Propose an algorithm of sparse representation of deformable surfaces.(ii)
Generalize the representation approach for organs involving both interior and exterior surfaces.(iii)
Present evaluation results conducted using computer models, ex vivo experiments based on 3D MRI scans of freshly excised porcine kidneys, and in vivo cardiac MRI scans of real patients.This paper is organized as follows: in Section2, we describe the proposed algorithm of sparse representation of deformable surface, denoted as SRDS thenceforth. Section 3 presents some experimental results using finite element model (FEM) data, ex vivo and in vivo experimental data. Finally, in Section 4, we finish with a few conclusions.


## 2. The SRDS Algorithm

The SRDS algorithm consists of three main steps to achieve sparse representations of deformable surfaces, as outlined in Figure1. Initially, SHD is performed to depict the deformable surfaces in the training set in the harmonic domain. Then OSP is applied in the transformed domain to identify the subspaces in which the SH coefficient vectors of the deformations can be linearly represented. Finally, each deformation is clustered to the proper subspace and represented with the corresponding coefficient vector with block sparsity. This representation method is also extended to organ deformations occurred on both interior and exterior layers as described in Section 2.4. Furthermore, as a practical issue, pixel-wise surface alignment among all the 3D surfaces is also addressed in Section 2.5.Figure 1
Flowchart of SRDS algorithm consists of three steps specified in the solid rectangulars; ellipsoids denote data input and output.


### 2.1. Step 1: SHD

Spherical harmonics are solutions to Laplace's equation expressed in the spherical coordinate system, defined as(1)Ylm(θ,φ)=(-1)m2l+14π(l-m)!(l+m)!Plm(cos⁡θ)eimφ,
where θ is the polar angle within [0, π], φ is the azimuthal angel within [0, 2π), l is the harmonic degree within [0, +∞], and m is the harmonic order varying in [-l, l]. Plm is the associated Legendre function. After proper parameterization [11, 12], a 3D surface x with finite energy can be expanded with SH series as(2)x(θ,φ)=∑l=0∞∑m=-l+lflmYlm(θ,φ).
Each harmonic coefficient flm is calculated using the inner product of the function x(θ,φ) and basis Ylm(θ,φ)(3)flm=∫φ=02π∫θ=0πx(θ,φ)Ylm(θ,φ)sinθdφdθ.Assume that harmonics up to levelL (l≤L) are involved in the transformation. Let Y denote the matrix composed of all (L+1)2 discretized harmonics, so Y has the following formation:(4)Y=|Y⃗0,0Y⃗1,-1Y⃗1,0Y⃗1,1⋯Y⃗L,L|N×(L+1)2.
Then a surface can be represented in the matrix format as(5)x=Yf⃗,
where x stands for a surface with N samples and f⃗=[f0,0f1,-1f1,0f1,1⋯fL,L]T is the harmonic coefficient vector. Notice that this equation is not exactly equal but approximately. For simplicity, we still use equal sign with least square estimation in this paper. The linear problem in (5) can be solved with the least square (LS) constraints outputting f⃗(6)f⃗=(YTY)-1YTx.Perform SHD for each of theK training deformed surfaces in X={xk}k=1K, so the group of deformations can be described by matrix F(7)F=|f⃗1f⃗2⋯f⃗K|(L+1)2×K.
Consequently, the training set of deformations can be totally characterized by columns in F as(8)X=YF.


### 2.2. Step 2: Subspace Identification with OSP

The aim of the second step is to explore the structures in those training deformations in the transformed harmonic domain and recognize the inherent subspaces in which the SHD coefficient vectors of the training deformations can be projected with high accuracy. The newly developed OSP algorithm [20] is adopted since it features better generalization and less computational cost compared to the gold standard K-SVD algorithm [31]. OSP is an iterative process that terminates when one of the predefined criteria is met. In this paper, we specify the following two stopping criteria: (1) an error threshold for ε subspaces selection and (2) a maximum number of iterations Emax for both controlling the subspace dimensions and avoiding deadlock searching. Further, the threshold for vector clustering is denoted as η, that is, we declare that a vector lives in a subspace if it can be projected to that subspace with error (l2 distance) no larger than η.A. Subspace Pursuit
Initially, each vectorf⃗k of length (L+1)2 in F is normalized by l2 norm. For convenience, we still use F to denote SH coefficient matrix even after normalization. The algorithm first identifies a subspace from F based on the stop criteria, then finds all the vectors in F that can be represented in that subspace with error level below η and remove those vectors from F to prepare for the next subspace pursuit. The process can be generalized as follows, in which A⊗B means that elements from B are excluded from A, and A⊕B stands for inclusion.
(1)
Initialization:i=0, D=∅, F0=F.(2)
Subspace searching and clustering
(i)
i=i+1; choose a vector f⃗i from Fi (e.g. first column of Fi) and let Fi=Fi-1⊗f⃗i to remove f⃗i.(ii)
Findni vectors from F′=F⊗f⃗i for representing f⃗i with error no larger than ɛ within Emax iterations and the ni vectors form Si.(iii)
Perform SVD decomposition onSi: UΣVT=Si; let Ai contain the first ni vectors of U; update D=D⊕Ai.(iv)
Select vectors fromFi that can be represented by Ai with error no larger than η, and then remove them from Fi.(v)
Repeat above steps until all the vectors are clustered.B. Subspace Pruning
One disadvantage of the traditional OSP algorithm is the presence of “spurious” or redundant subspaces especially as the dimension of the training data set increases. Those subspaces identified in the earlier iterations actually can be better represented by the later identified subspaces. Therefore, a postprocessing step is used to identify and then discard the redundancy among the subspaces without decreasing the performance. This is implemented by repartitioning the training data among the initial subspaces and then eliminating subspaces in which very few or no vector is clustered. In some cases, where the subspace size is limited to some constraint, an optimization step can be applied in conjunction with the pruning step. The details of the subspace optimization design is described in [20].C. MatrixF Factorization
After identifying the inherent subspaces, the coefficient matrixF of training set can be partitioned into two-part union of subspaces and corresponding coefficients via the following procedures.
Since each vectorf⃗k has been clustered into the belonging subspace during the subspace identification process, the corresponding coefficients for each f⃗k can be obtained accordingly. Suppose that f⃗k lives in subspace Ai which is spanned by ni orthogonal basis, so its corresponding coefficients c⃗k can be calculated using LS estimator(9)c⃗k=(AiTAi)-1AiTf⃗k.
Then f⃗k can be characterized by c⃗k in its subspace
(10)f⃗k=Aic⃗k.
If there are totallyJ subspaces identified from F, a structured dictionary constructed by concatenating all deformation subspaces is established as D=⋃i=1J{Ai}, with dimension I=∑i=1Jni. Since each vector f⃗k lies in one of the subspaces, f⃗k can also be represented in the structured dictionary with a block sparse vector {c⃗̃k}, which is obtained via extending the coefficients {c⃗k}k=1K by zero padding in positions corresponding to other subspaces in D. Figure 2 provides an example of 3 subspaces to illustrate the sparsity of coefficient vector c⃗̃1. If f⃗1 lies in subspace A2 which are spanned by the 5th, 6th, and 7th columns in D, then c⃗̃1 has nonzero values only at index of 5, 6, and 7. Consequently, F can be factorized as(11)F=DC,
where C=|c⃗̃1c⃗̃2⋯c⃗̃K|I×K is the corresponding coefficient matrix with block sparsity.Figure 2
Sparsity of coefficient vectorc⃗̃1.


### 2.3. Step 3: Structured Sparse Surface Representation

A. Sparse Representation of Training Surfaces
Integrating the subspace pursuing results in the harmonic domain in (11) with the initial SHD process in (8), the training deformations X can be sparsely represented in the original spatial domain as(12)X=YDC=GC,
where G=YD with size of N×I is the desired structured dictionary in the spatial domain. Since D=⋃i=1J{Ai}, G is inherently structured by subspaces of G=⋃i=1J{Gi} with Gi=YAi of size N×ni.
Up to this point, with matrixG that captures the deformation features in the considered population, each training deformation xk in X can be fully characterized by a sparse coefficient c⃗̃k as(13)xk=Gc⃗̃k.
The sparsity of c⃗̃k has already been illustrated in Figure 2.B. Sparse Representation of Testing Surfaces
For the testing deformations beyond the training set, we utilized the fact that the dictionary identified from an extensive training data features good generalization such that any deformation of that particular population can be represented in the subspaces with high accuracy. This is justified because organs only deform in limited ways due to their mechanical properties, so the deformation variations can be fully learned from a training data set. This applied structure allows fast deformation representation in subspaces of low dimensionality.
The testing set is denoted asH={hm}m=1M, where M is the number of deformations to be represented. The straightforward strategy is to find a best-fit subspace for hm by projecting it to every subspace {Gi}i=1J and choose the subspace with minimal projection error. Since the number of subspaces J and the dimension of each subspace {ni}i=1J are both small with the postprocessing of subspace pruning, this best-fit strategy still results in low computational cost. However, when the number of subspaces is too large, an alternative threshold approach can be applied by finding a subspace Gi in which hm can be represented with an error level smaller than η. The former best-fit method is implemented in Section 3 for performance validation.
Suppose thatGi* is the chosen subspace, so coefficient vector c⃗m can be estimated with LS as(14)c⃗m=(Gi*TGi*)-1Gi*Thm.
Then block sparse vector
(15)c⃗̃m=[00⋯c⃗m⋯000]
is obtained according to the rules described in Section 2.2. Further, the sparsity of c⃗̃k or c⃗̃m can be increased by trimming off nonzero elements with absolute value lower than a given threshold δ.
It is worth noting that, different from the traditional learning approaches relying on orthogonal least square (OLS) [32] or matching pursuit (MP) [33] algorithms which select atoms from the training set and recombine them for representing each surface in the testing set sparsely, our SRDS algorithm avoids this heavy overload caused by reshuffling all the atoms. Instead, we apply the block structure of the dictionary learned from a representative training data set. This essentially enables the representation of each deformable surface compactly and sparsely with high accuracy and low computational cost.


### 2.4. Extended Sparse Surface Representation

For an organ with both interior and exterior surfaces, such as bladder, deformations can take place on both layers. The above theory can be extended to achieve sparse surface representation for deformations occurred on both interior and exterior wall of the object.Initially, spherical parameterization is conducted on interior and exterior parts separately. We denotex⃗kin and x⃗kex(1≤k≤K) as the corresponding interior and exterior of each training surface x⃗k. Then each pair of x⃗kin (with N1 vertices) and x⃗kex (with N2 vertices) can be approximated by spherical harmonic basis as(16)x⃗k=[x⃗kin̂x⃗kex̂]=[YinOOYex]⋅[f⃗kinf⃗kex]
where Yin of size N1×(L+1)2 and Yex of size N2×(L+1)2 denote the spherical harmonic basis for inner and outer surfaces, respectively. L is the highest degree of harmonics included. f⃗kin and f⃗kex are the corresponding harmonic coefficient vectors. Therefore, each deformation is represented by vector f⃗k=[f⃗kinf⃗kex], and all K training frames can be characterized by F=⋃{f⃗k}k=1K as(17)X=YF=[YinOOYex][f⃗1in⋯f⃗Kinf⃗1ex⋯f⃗Kex].The following procedures of subspace identification and sparse surface representation as described in Sections2.2 and 2.3 can be applied straightforwardly. After identifying J subspaces D=⋃i=1J{Ai} from SH coefficient matrix F, each training deformation can be sparsely represented with block sparse coefficient vector c⃗̃k as:(18)x⃗k=[x⃗kin̂x⃗kex̂]=[YinOOYex]Dc⃗̃k=Gc⃗̃k=Gic⃗k,
where Gi=[YinOOYex]Ai is the subspace with size of (N1+N2)×ni, and G=⋃i=1J{Gi} is the desired structured dictionary. Accordingly, c⃗̃k is the block sparse vector, and c⃗k is the nonzero coefficient values in the selected subspace.


### 2.5. Surface Correspondence

Similar to other surface modeling methods [25, 27, 28], the proposed approach requires point-wise correspondence across difference surfaces besides rigid registration [34]. Specifically, this point-to-point alignment is established, such that an identical spherical parameterization can be applied in the SHD procedure. Figure 3 illustrates the goal of surface correspondence. Same colored vertices on deformation 1 and 2 indicate a matched pair. After established correspondence of the point pairs over the two deformations, vertices on deformation 2 can be numbered in the same order as deformation 1.Surface correspondence: (a) vertex index on deformation 1, (b) vertex index on deformation 2 before correspondance, (c) corresponding vertex index on deformation 2.
(a)(b)(c)Different correspondence methods have been proposed, such as minimum description length [35], SH-coefficient alignment [36], and so forth. We applied the SH based method [36] in this paper as well as ray-casting method for simple surfaces. The former SH-based method is based on the underlying fact that two points with the same parameter pair when mapped to a sphere are considered to be a corresponding pair. Therefore, it fixes parameterization of the template and rotate the other to optimize the surface correspondence by minimizing the root mean squared distance of the two SH coefficient vectors. The latter ray-casting method is introduced in the following section.


#### 2.5.1. Point Correspondence with Ray Casting

For surfaces, if unique intersection exists between a ray starting from its object center and the surface, a ray-casting method can be applied to obtain sample pairs across all the deformations. For illustration, Figures4(a) and 4(b) depict two different cases of ray-surface intersection in the simplified 2D space. In Figure 4(a), there is only one intersected point (p1 for ray r1⃗, and p2 for ray r2⃗) between each ray and the deformed surface S1. By contrast, Figure 4(b) gives an example when multiple intersections (p1, p1′, p1′′ for ray r1⃗ and p2 for ray r2⃗) are involved between rays and the surface S2.Ray casting sampling: (a) unique ray-surface intersection, (b) multiple ray-surface intersections, and (c) ray casting on two deformed surfaces.
(a)(b)(c)If the condition of single ray-surface intersection applies, deformations can be resampled through the following steps to achieve point correspondence.(i)
Construct an icosahedron ofW vertices with radius large enough to embrace the largest deformation volume among those under consideration; larger W results in denser samples to maintain the local details but incurs more computational cost.(ii)
Align the center of the 3D surfaces to the origin of icosahedron such that rays casting from the origin can intersect with the surface.(iii)
For each ray segment originated from the center to a vertex on the icosahedron, find the triangle on the surface mesh that intersects with the segment and use that intersected point as new surface sample.Figure4(c) illustrates the desired sample pairs as (p1,1,p2,1) and (p1,2,p2,2) over two deformations S1 and S2. This resampling process also establishes a one-to-one map between a point on the object and a point on the sphere (icosahedron), which naturally meets the purpose of spherical parameterization. As a result, point-wise correspondence can be achieved across all the resampled surfaces, and a uniform spherical harmonic matrix Y can be applied. As an example, Figures 5(a) and 5(b) compare an original kidney surface with the corresponding resampled surface. We can see that ray-casting procedure well maintains the shape.Resample with ray casting (a) original surface of kidney, (b) resampled surface of kidney.
(a)(b)


## 2.1. Step 1: SHD

Spherical harmonics are solutions to Laplace's equation expressed in the spherical coordinate system, defined as(1)Ylm(θ,φ)=(-1)m2l+14π(l-m)!(l+m)!Plm(cos⁡θ)eimφ,
where θ is the polar angle within [0, π], φ is the azimuthal angel within [0, 2π), l is the harmonic degree within [0, +∞], and m is the harmonic order varying in [-l, l]. Plm is the associated Legendre function. After proper parameterization [11, 12], a 3D surface x with finite energy can be expanded with SH series as(2)x(θ,φ)=∑l=0∞∑m=-l+lflmYlm(θ,φ).
Each harmonic coefficient flm is calculated using the inner product of the function x(θ,φ) and basis Ylm(θ,φ)(3)flm=∫φ=02π∫θ=0πx(θ,φ)Ylm(θ,φ)sinθdφdθ.Assume that harmonics up to levelL (l≤L) are involved in the transformation. Let Y denote the matrix composed of all (L+1)2 discretized harmonics, so Y has the following formation:(4)Y=|Y⃗0,0Y⃗1,-1Y⃗1,0Y⃗1,1⋯Y⃗L,L|N×(L+1)2.
Then a surface can be represented in the matrix format as(5)x=Yf⃗,
where x stands for a surface with N samples and f⃗=[f0,0f1,-1f1,0f1,1⋯fL,L]T is the harmonic coefficient vector. Notice that this equation is not exactly equal but approximately. For simplicity, we still use equal sign with least square estimation in this paper. The linear problem in (5) can be solved with the least square (LS) constraints outputting f⃗(6)f⃗=(YTY)-1YTx.Perform SHD for each of theK training deformed surfaces in X={xk}k=1K, so the group of deformations can be described by matrix F(7)F=|f⃗1f⃗2⋯f⃗K|(L+1)2×K.
Consequently, the training set of deformations can be totally characterized by columns in F as(8)X=YF.


## 2.2. Step 2: Subspace Identification with OSP

The aim of the second step is to explore the structures in those training deformations in the transformed harmonic domain and recognize the inherent subspaces in which the SHD coefficient vectors of the training deformations can be projected with high accuracy. The newly developed OSP algorithm [20] is adopted since it features better generalization and less computational cost compared to the gold standard K-SVD algorithm [31]. OSP is an iterative process that terminates when one of the predefined criteria is met. In this paper, we specify the following two stopping criteria: (1) an error threshold for ε subspaces selection and (2) a maximum number of iterations Emax for both controlling the subspace dimensions and avoiding deadlock searching. Further, the threshold for vector clustering is denoted as η, that is, we declare that a vector lives in a subspace if it can be projected to that subspace with error (l2 distance) no larger than η.A. Subspace Pursuit
Initially, each vectorf⃗k of length (L+1)2 in F is normalized by l2 norm. For convenience, we still use F to denote SH coefficient matrix even after normalization. The algorithm first identifies a subspace from F based on the stop criteria, then finds all the vectors in F that can be represented in that subspace with error level below η and remove those vectors from F to prepare for the next subspace pursuit. The process can be generalized as follows, in which A⊗B means that elements from B are excluded from A, and A⊕B stands for inclusion.
(1)
Initialization:i=0, D=∅, F0=F.(2)
Subspace searching and clustering
(i)
i=i+1; choose a vector f⃗i from Fi (e.g. first column of Fi) and let Fi=Fi-1⊗f⃗i to remove f⃗i.(ii)
Findni vectors from F′=F⊗f⃗i for representing f⃗i with error no larger than ɛ within Emax iterations and the ni vectors form Si.(iii)
Perform SVD decomposition onSi: UΣVT=Si; let Ai contain the first ni vectors of U; update D=D⊕Ai.(iv)
Select vectors fromFi that can be represented by Ai with error no larger than η, and then remove them from Fi.(v)
Repeat above steps until all the vectors are clustered.B. Subspace Pruning
One disadvantage of the traditional OSP algorithm is the presence of “spurious” or redundant subspaces especially as the dimension of the training data set increases. Those subspaces identified in the earlier iterations actually can be better represented by the later identified subspaces. Therefore, a postprocessing step is used to identify and then discard the redundancy among the subspaces without decreasing the performance. This is implemented by repartitioning the training data among the initial subspaces and then eliminating subspaces in which very few or no vector is clustered. In some cases, where the subspace size is limited to some constraint, an optimization step can be applied in conjunction with the pruning step. The details of the subspace optimization design is described in [20].C. MatrixF Factorization
After identifying the inherent subspaces, the coefficient matrixF of training set can be partitioned into two-part union of subspaces and corresponding coefficients via the following procedures.
Since each vectorf⃗k has been clustered into the belonging subspace during the subspace identification process, the corresponding coefficients for each f⃗k can be obtained accordingly. Suppose that f⃗k lives in subspace Ai which is spanned by ni orthogonal basis, so its corresponding coefficients c⃗k can be calculated using LS estimator(9)c⃗k=(AiTAi)-1AiTf⃗k.
Then f⃗k can be characterized by c⃗k in its subspace
(10)f⃗k=Aic⃗k.
If there are totallyJ subspaces identified from F, a structured dictionary constructed by concatenating all deformation subspaces is established as D=⋃i=1J{Ai}, with dimension I=∑i=1Jni. Since each vector f⃗k lies in one of the subspaces, f⃗k can also be represented in the structured dictionary with a block sparse vector {c⃗̃k}, which is obtained via extending the coefficients {c⃗k}k=1K by zero padding in positions corresponding to other subspaces in D. Figure 2 provides an example of 3 subspaces to illustrate the sparsity of coefficient vector c⃗̃1. If f⃗1 lies in subspace A2 which are spanned by the 5th, 6th, and 7th columns in D, then c⃗̃1 has nonzero values only at index of 5, 6, and 7. Consequently, F can be factorized as(11)F=DC,
where C=|c⃗̃1c⃗̃2⋯c⃗̃K|I×K is the corresponding coefficient matrix with block sparsity.Figure 2
Sparsity of coefficient vectorc⃗̃1.


## 2.3. Step 3: Structured Sparse Surface Representation

A. Sparse Representation of Training Surfaces
Integrating the subspace pursuing results in the harmonic domain in (11) with the initial SHD process in (8), the training deformations X can be sparsely represented in the original spatial domain as(12)X=YDC=GC,
where G=YD with size of N×I is the desired structured dictionary in the spatial domain. Since D=⋃i=1J{Ai}, G is inherently structured by subspaces of G=⋃i=1J{Gi} with Gi=YAi of size N×ni.
Up to this point, with matrixG that captures the deformation features in the considered population, each training deformation xk in X can be fully characterized by a sparse coefficient c⃗̃k as(13)xk=Gc⃗̃k.
The sparsity of c⃗̃k has already been illustrated in Figure 2.B. Sparse Representation of Testing Surfaces
For the testing deformations beyond the training set, we utilized the fact that the dictionary identified from an extensive training data features good generalization such that any deformation of that particular population can be represented in the subspaces with high accuracy. This is justified because organs only deform in limited ways due to their mechanical properties, so the deformation variations can be fully learned from a training data set. This applied structure allows fast deformation representation in subspaces of low dimensionality.
The testing set is denoted asH={hm}m=1M, where M is the number of deformations to be represented. The straightforward strategy is to find a best-fit subspace for hm by projecting it to every subspace {Gi}i=1J and choose the subspace with minimal projection error. Since the number of subspaces J and the dimension of each subspace {ni}i=1J are both small with the postprocessing of subspace pruning, this best-fit strategy still results in low computational cost. However, when the number of subspaces is too large, an alternative threshold approach can be applied by finding a subspace Gi in which hm can be represented with an error level smaller than η. The former best-fit method is implemented in Section 3 for performance validation.
Suppose thatGi* is the chosen subspace, so coefficient vector c⃗m can be estimated with LS as(14)c⃗m=(Gi*TGi*)-1Gi*Thm.
Then block sparse vector
(15)c⃗̃m=[00⋯c⃗m⋯000]
is obtained according to the rules described in Section 2.2. Further, the sparsity of c⃗̃k or c⃗̃m can be increased by trimming off nonzero elements with absolute value lower than a given threshold δ.
It is worth noting that, different from the traditional learning approaches relying on orthogonal least square (OLS) [32] or matching pursuit (MP) [33] algorithms which select atoms from the training set and recombine them for representing each surface in the testing set sparsely, our SRDS algorithm avoids this heavy overload caused by reshuffling all the atoms. Instead, we apply the block structure of the dictionary learned from a representative training data set. This essentially enables the representation of each deformable surface compactly and sparsely with high accuracy and low computational cost.


## 2.4. Extended Sparse Surface Representation

For an organ with both interior and exterior surfaces, such as bladder, deformations can take place on both layers. The above theory can be extended to achieve sparse surface representation for deformations occurred on both interior and exterior wall of the object.Initially, spherical parameterization is conducted on interior and exterior parts separately. We denotex⃗kin and x⃗kex(1≤k≤K) as the corresponding interior and exterior of each training surface x⃗k. Then each pair of x⃗kin (with N1 vertices) and x⃗kex (with N2 vertices) can be approximated by spherical harmonic basis as(16)x⃗k=[x⃗kin̂x⃗kex̂]=[YinOOYex]⋅[f⃗kinf⃗kex]
where Yin of size N1×(L+1)2 and Yex of size N2×(L+1)2 denote the spherical harmonic basis for inner and outer surfaces, respectively. L is the highest degree of harmonics included. f⃗kin and f⃗kex are the corresponding harmonic coefficient vectors. Therefore, each deformation is represented by vector f⃗k=[f⃗kinf⃗kex], and all K training frames can be characterized by F=⋃{f⃗k}k=1K as(17)X=YF=[YinOOYex][f⃗1in⋯f⃗Kinf⃗1ex⋯f⃗Kex].The following procedures of subspace identification and sparse surface representation as described in Sections2.2 and 2.3 can be applied straightforwardly. After identifying J subspaces D=⋃i=1J{Ai} from SH coefficient matrix F, each training deformation can be sparsely represented with block sparse coefficient vector c⃗̃k as:(18)x⃗k=[x⃗kin̂x⃗kex̂]=[YinOOYex]Dc⃗̃k=Gc⃗̃k=Gic⃗k,
where Gi=[YinOOYex]Ai is the subspace with size of (N1+N2)×ni, and G=⋃i=1J{Gi} is the desired structured dictionary. Accordingly, c⃗̃k is the block sparse vector, and c⃗k is the nonzero coefficient values in the selected subspace.


## 2.5. Surface Correspondence

Similar to other surface modeling methods [25, 27, 28], the proposed approach requires point-wise correspondence across difference surfaces besides rigid registration [34]. Specifically, this point-to-point alignment is established, such that an identical spherical parameterization can be applied in the SHD procedure. Figure 3 illustrates the goal of surface correspondence. Same colored vertices on deformation 1 and 2 indicate a matched pair. After established correspondence of the point pairs over the two deformations, vertices on deformation 2 can be numbered in the same order as deformation 1.Surface correspondence: (a) vertex index on deformation 1, (b) vertex index on deformation 2 before correspondance, (c) corresponding vertex index on deformation 2.
(a)(b)(c)Different correspondence methods have been proposed, such as minimum description length [35], SH-coefficient alignment [36], and so forth. We applied the SH based method [36] in this paper as well as ray-casting method for simple surfaces. The former SH-based method is based on the underlying fact that two points with the same parameter pair when mapped to a sphere are considered to be a corresponding pair. Therefore, it fixes parameterization of the template and rotate the other to optimize the surface correspondence by minimizing the root mean squared distance of the two SH coefficient vectors. The latter ray-casting method is introduced in the following section.


### 2.5.1. Point Correspondence with Ray Casting

For surfaces, if unique intersection exists between a ray starting from its object center and the surface, a ray-casting method can be applied to obtain sample pairs across all the deformations. For illustration, Figures4(a) and 4(b) depict two different cases of ray-surface intersection in the simplified 2D space. In Figure 4(a), there is only one intersected point (p1 for ray r1⃗, and p2 for ray r2⃗) between each ray and the deformed surface S1. By contrast, Figure 4(b) gives an example when multiple intersections (p1, p1′, p1′′ for ray r1⃗ and p2 for ray r2⃗) are involved between rays and the surface S2.Ray casting sampling: (a) unique ray-surface intersection, (b) multiple ray-surface intersections, and (c) ray casting on two deformed surfaces.
(a)(b)(c)If the condition of single ray-surface intersection applies, deformations can be resampled through the following steps to achieve point correspondence.(i)
Construct an icosahedron ofW vertices with radius large enough to embrace the largest deformation volume among those under consideration; larger W results in denser samples to maintain the local details but incurs more computational cost.(ii)
Align the center of the 3D surfaces to the origin of icosahedron such that rays casting from the origin can intersect with the surface.(iii)
For each ray segment originated from the center to a vertex on the icosahedron, find the triangle on the surface mesh that intersects with the segment and use that intersected point as new surface sample.Figure4(c) illustrates the desired sample pairs as (p1,1,p2,1) and (p1,2,p2,2) over two deformations S1 and S2. This resampling process also establishes a one-to-one map between a point on the object and a point on the sphere (icosahedron), which naturally meets the purpose of spherical parameterization. As a result, point-wise correspondence can be achieved across all the resampled surfaces, and a uniform spherical harmonic matrix Y can be applied. As an example, Figures 5(a) and 5(b) compare an original kidney surface with the corresponding resampled surface. We can see that ray-casting procedure well maintains the shape.Resample with ray casting (a) original surface of kidney, (b) resampled surface of kidney.
(a)(b)


## 2.5.1. Point Correspondence with Ray Casting

For surfaces, if unique intersection exists between a ray starting from its object center and the surface, a ray-casting method can be applied to obtain sample pairs across all the deformations. For illustration, Figures4(a) and 4(b) depict two different cases of ray-surface intersection in the simplified 2D space. In Figure 4(a), there is only one intersected point (p1 for ray r1⃗, and p2 for ray r2⃗) between each ray and the deformed surface S1. By contrast, Figure 4(b) gives an example when multiple intersections (p1, p1′, p1′′ for ray r1⃗ and p2 for ray r2⃗) are involved between rays and the surface S2.Ray casting sampling: (a) unique ray-surface intersection, (b) multiple ray-surface intersections, and (c) ray casting on two deformed surfaces.
(a)(b)(c)If the condition of single ray-surface intersection applies, deformations can be resampled through the following steps to achieve point correspondence.(i)
Construct an icosahedron ofW vertices with radius large enough to embrace the largest deformation volume among those under consideration; larger W results in denser samples to maintain the local details but incurs more computational cost.(ii)
Align the center of the 3D surfaces to the origin of icosahedron such that rays casting from the origin can intersect with the surface.(iii)
For each ray segment originated from the center to a vertex on the icosahedron, find the triangle on the surface mesh that intersects with the segment and use that intersected point as new surface sample.Figure4(c) illustrates the desired sample pairs as (p1,1,p2,1) and (p1,2,p2,2) over two deformations S1 and S2. This resampling process also establishes a one-to-one map between a point on the object and a point on the sphere (icosahedron), which naturally meets the purpose of spherical parameterization. As a result, point-wise correspondence can be achieved across all the resampled surfaces, and a uniform spherical harmonic matrix Y can be applied. As an example, Figures 5(a) and 5(b) compare an original kidney surface with the corresponding resampled surface. We can see that ray-casting procedure well maintains the shape.Resample with ray casting (a) original surface of kidney, (b) resampled surface of kidney.
(a)(b)


## 3. Experiments

Three types of experiments are conducted to demonstrate the feasibility of the proposed SRDS algorithm. The computer-generated FEM data is first used to demonstrate that the SRDS approach matches the accuracy of complex mathematical modeling techniques, then an ex vivo experiment is conducted using 3D MRI scans of porcine kidneys for evaluation in practical settings, and finally in vivo experiment is carried over dynamic cardiac MRI scans for evaluation in real patients.


### 3.1. Experiment with FEM Data

Three representative organs are employed in this FEM experiment: 3D cortical mesh as an example of complicated shapes, gallbladder as an instance with geometrically simple shape, and bladder consisting of both interior and exterior walls.


#### 3.1.1. Computer Model Setup

The initial 3D models of different organs are fed into a FEM-based surgical simulation tool to generate deformation data for testing. For instance, Figure6 demonstrates two examples of shape distortions due to the endoscope poking and grasping one side of the gallbladder.Gallbladder deformation generated by surgical simulator: (a) and (b) are distortions caused when the endoscope pokes and grasps the gallbladder.
(a)(b)Table1 lists the FEM experimental setup of the three organs including number of vertices N, SH level L, number of deformations for training K, and testing M. “GBL” stands for gallbladder in all the tables. The maximum SH level used for brain model is chosen according to [9], and the levels for gallbladder and bladder are determined when the SHD representation error is below 0.1% (EOF). The complex brain structure requires more vertices and higher SH level for surface representation to achieve sufficient accuracy. To evaluate the representation precision qualitatively, an evaluation parameter EOF is defined as the normalized Euclidean distance between the original surface and the reconstructed surface(19)EOF=‖x̂k-xk‖2‖xk‖2.
All surfaces are centered to the origin of the coordinate system so that EOF will not be heavily affected by the denominator.Table 1
FEM Model setup.
VerticesNSH levelLTrainingKTestingMBrain40962803535GBL303825250114BladderN1 = 4434,307446N2 = 4274


#### 3.1.2. Results

With the FEM data, the proposed SRDS algorithm is evaluated from three perspectives: (1) subspace dimensionality, (2) sparsity and accuracy of representations, and (3) the effect of subspace pursuit thresholdɛ and coefficient truncation threshold δ on the performance. The sparsity is defined as the l0 norm of the coefficient vector c⃗̃m.A. Training Results
During training stage, we setɛ=0.005 for subspace detection, η=0.01 for clustering, Emax=50 as the maximum iteration times, and δ=0.005 for coefficient truncation. Subspaces on X, Y, and Z axis are identified separately. Table 2 shows that the subspace number J and dimensions of resulting dictionary (dim(G)=I) are markedly small relative to N or L2 in all three tests. We notice that the subspace dimensions of brain are relatively smaller than the other two. This is because of smaller training data size and minor extent of deformation considered in the brain experiment, which results in smaller dictionary size to capture the deformation features.Table 2
Dimension of dictionary (J/I).
BrainGBLBladderSubspace onX (J/I)2/132/391/37Subspace onY (J/I)2/62/481/30Subspace onZ (J/I)1/31/421/34B. Sparsity and Accuracy Evaluation
Sparsity is examined in terms of (μ/σ), where μ is the average l0 norm of the coefficient vector c⃗̃k (training) or c⃗̃m (testing) and σ is the corresponding standard deviation. To verify that whether our method achieves equivalent sparsity and precision when applying the structure of the dictionary, we also test the case without relying on any structure learned from training set, during which sparse representation of each deformation in the testing set is repursued from the training set using OSP approach. In the following tables, we use “OSP” to refer to the results obtained using such a repursuing process. Table 3 summarizes the sparsity of the SRDS representation of three organs for both training and testing set. It illustrates that the number of atoms needed for representing the complex deformations is much smaller than the dimension of spherical harmonic vectors ((L+1)2), and particularly the sparsity and accuracy via SRDS is very close to that from OSP repursuit for the testing deformations, which indicates the good generalization of the structured dictionary.
The reconstruction error in terms of EOF is further compared with that from standard SHD method, as shown in Figure7. In general, the accuracy of SRDS is equivalent to that of SHD method. Specifically, it shows that the SRDS method achieves average EOF of 1.32% (brain) and 0.14% (gallbladder) versus 1.29% (brain) and 0.13% (gallbladder) with SHD method. For bladder model with deformations on multiple layers, the overall representation error with SRDS is 0.07%, very close to 0.06% with SHD. Figures 11 and 12 show typical reconstructed deformations of the testing data for the three organs with SHD and SRDS methods. Figures 12(e) and 12(f) demonstrate the interior and exterior representation of the bladder at a same time instance. From those results, we can see that the SRDS algorithm achieves the accuracy equivalent to complex mathematical modeling techniques while significantly lowers the representation dimensionality.Table 3
Sparsity (μ/σ) and accuracy evaluation.
X (μ/σ)Y (μ/σ)Z (μ/σ)EOF (%)Brain train (SRDS)3.1/0.23.4/0.93.0/0.21.31Brain test (SRDS)3.5/0.53.9/0.43.0/0.01.32Brain test (OSP)4.1/0.64.5/0.53.9/0.81.30GBL train (SRDS)25.8/12.24.8/1.832.3/7.30.13GBL test (SRDS)33.2/1.044.1/1.140.8/1.30.15GBL test (OSP)44.4/4.048.3/2.343.9/4.50.13Bladder train (SRDS)24.1/3.720.8/2.924.8/4.00.076Bladder test (SRDS)24.6/2.721.0/2.525.6/3.20.073Bladder test (OSP)22.5/4.316.3/6.523.0/5.20.070EOF of brain, gallbladder, and bladder reconstruction with SHD and SRDS methods. The left pair is for training set; the right pair is for testing set.
(a)(b)(c)C. Effect ofɛ
The performance of SRDS algorithm is examined as the subspace pursuit thresholdɛ varies. Specifically, we study the effect of ɛ on the dimensionality (I) of the structured dictionary G, sparsity and accuracy of the surface representation. Figure 8 shows how the subspace dimensions on three axis change during the training stage as ɛ increases from 0.001 to 0.01. Figure 9 displays the influence of ɛ on the average sparsity μ of the surface representation results for both training and testing data sets. In general, smaller ɛ leads to larger subspace size and less description sparsity, since lower ɛ usually leads to more recruited atoms to meet the desired representation accuracy. Therefore, there is a tradeoff between representation accuracy and desired sparsity. Figure 10 reveals the representation EOF as a function of ɛ. Not surprisingly, the reconstruction error is increased as ɛ becomes larger. An empirical point can be chosen according to the training curve when space dimension I expands significantly but only trivial EOF improvement is gained, that is, ɛ=0.005  is a preferred value in this test according to Figure 10.Figure 8
DimensionI of trained dictionary G decreases on all three axis as ɛ increases.Sparsity (μ) of surface representation changes with different ɛ on three axis: (a) is for training set; (b) is for testing set.
(a)(b)Figure 10
EOF of gallbladder reconstruction increases asɛ gets larger for both training and testing sets.Reconstructed brain deformations using SHD approach: (a) is the initial brain shape and the circle marks one typical area under deformation, (b)–(e) are the reconstructed brain deformations.
(a)(b)(c)(d)(e)(f)Reconstructed brain deformations using SRDS approach: (a) is the initial brain shape, (b)–(e) are the reconstructed brain deformations.
(a)(b)(c)(d)(e)(f)D. Effect ofδ
The influence of coefficient truncation thresholdδ on the performance of SRDS algorithm is also tested while δ is varied among [00.00010.00050.0010.0050.010.050.080.1]. Figure 13 shows the effect of δ on the average sparsity μ of the surface representation results. We can see that, as the truncation threshold δ enlarges, the sparsity of the representation is increased for both training and testing data sets at the price of decreased representation error as shown in Figure 14, so there is tradeoff between sparsity and accuracy. Empirically, one can choose the δ value when the representation precision remarkably deteriorates while the sparsity is still increasing. Therefore, according to Figures 13 and 14, an appropriate value for δ is between 0.005 and 0.01.Sparsity (μ) of surface representation changes with different δ on three axis: (a) is for training set; (b) is for testing set.
(a)(b)Figure 14
EOF of gallbladder reconstruction increases asδ gets larger for both training and testing sets.


### 3.2. Ex vivo Experiment Using MRI

To evaluate the proposed algorithm in real applications, an ex vivo experiment using three porcine kidneys were conducted at the Center for Interdisciplinary Applications in Magnetic Resonance (CIA-MR) of University of Minnesota. Deformations imposed to each kidney were controlled and maintained still during imaging by a customized nonmagnetic mechanical device as shown in Figure15. Each deformed kidney shape was scanned in 3D MRI mode with spatial resolution of 1.2 mm to generate both training set and testing set. The SH degree L of the organ representation is set to be 20, and each 3D kidney mesh after surface correspondence has N=4002 vertices. Different from computer-generated deformations where surface correspondence is intrinsically established, the shapes from MRI scans are rendered independently, so the method described in Section 2.5 is applied to achieve point-wise correspondence.Images of three porcine kidneys for experiment.
(a)(b)(c)(d)(e)Both intramodel and intermodel experiments are conducted. The former uses training and testing deformations from the same kidney; the later utilizes two out of the three kidneys for training and the third one for testing in a cross-evaluation fashion. Besides sparsity and EOF, Hausdorff distance between the represented shape and corresponding MRI surface is also examined as a physical measurement of error. The Hausdorff distance between surfacex and x′ is defined as (20)d(x,x′)=max⁡p∈x⁡d(p,x′),
where d(p,x′) is defined as the distance between a point p on surface x and the closest point on surface x′, that is,(21)d(p,x′)=min⁡p′∈x′⁡‖p,p′‖2
with ∥·∥2 denoting the Euclidean norm.


#### 3.2.1. Intramodel Test

In the intramodel experiment, 31 deformations of the same kidney were generated and scanned by the MRI machine, among which 20 frames were randomly selected as training set, and the other 11 were applied for testing the generalization of the learned subspaces.Table4 shows the trained subspace dimensions (J as number of subspace, I as dictionary size of G), the sparsity of the descriptors in each axis for both training set and testing set, and the corresponding errors in terms of EOF and Hausdorff distance. Similar to the FEM experiment, the sparsity is also evaluated with OSP repursuit process in the testing set for comparison. The table shows that the sparsity and the accuracy achieved with SRDS is very close to that from OSP repursuing process. However, the SRDS method features delay-free surface representation by applying the structure in the identified dictionary. Further results about computational efficiency are shown in Section 3.4. One may notice that the size of training data in the MRI experiment is smaller than that in FEM test due to the less availability of 3D MRI images. As a rule of thumb, larger training set carries richer deformation information and thus leads to better generalization of the dictionary. However, given the size of training data and extent of deformation involved in the ex vivo experiment, high representation precision is still achieved.Table 4
Subspace dimension and sparsity (μ/σ) for intramodel experiment.
X (μ/σ)Y (μ/σ)Z (μ/σ)EOF (%)Haus (mm)J/I1/172/312/31Train (SRDS)16.9/0.415.5/1.515.8/1.50.240.55Test (SRDS)17.0/0.016.9/0.316.7/0.90.640.87Test (OSP)16.6/1.516.3/0.917.1/0.80.600.85Figure16 illustrates the accuracy of the surface representation in the intramodel test. The average EOF in Figure 16(a) for training set is 0.24% and 0.64% for testing set, with maximum rate less than 1%. Further, error as Hausdorff distance (shown in Figure 16(b)) is 0.55±0.23 mm with 95th percentile error of 0.86 mm for the training set, and 0.87±0.10 mm with 95th percentile error of 0.96 mm for the testing set. This intramodel experiment demonstrated that the SRDS algorithm identifies subspaces generalizable enough to accurately represent deformations beyond the training set for the same object.Boxplots of representation error in intramodel experiment: (a) EOF of training and testing set, (b) Hausdorff distance of training and testing sets.
(a)(b)Figure17 visualizes the color-coded error distribution at all vertices on the represented surface with SRDS relative to the actual MRI scans. Figure 17(a) illustrates the error range for different colors. Figures 17(b) and 17(d) show the error distribution for a typical reconstruction in the training and testing set, respectively. Figures 17(c) and 17(e) show maximum 90% level reconstruction errors, that is, 90% of all deformations in the training or testing set have representation point errors less than the values shown in the figures. Consistent with the EOF and Hausdorff distance results, the color diffusion in Figure 17 indicates that the precision in the testing group is relatively lower than that in the training group. However, among all the pixel-wise errors shown in Figure 17(e), less than 3% of all the surface points have error distance larger than 0.5 mm.Representation error in intramodel experiment: (a) color-coded scales, (b) example representation error in training set, (c) 90% representation error in training set, (d) example of representation error in testing set, and (e) 90% representation error of testing set.
(a)(b)(c)(d)(e)


#### 3.2.2. Intermodel Test

Three intermodel experiments are performed to further validate the proposed SRDS method applied to organs from different subjects. In the following context, “Ex1” stands for the experiment training on Kidney 2 and 3 plus one initial shape of Kidney 1 while testing on deformations of Kidney 1, and the like for “Ex2” and “Ex3”. In each experiment, both sparsity and accuracy are examined.The number of subspaces (J) and dimensions (I) of the identified dictionary are listed in Table 5. The training results vary among the three experiments but all features low subspace dimensions. Table 6 shows the sparsity of the intermodel experiments using the SRDS algorithm, and the error level is evaluated in terms of EOF and Hausdorff distance. Each testing deformation is also sparsely retrained using OSP for comparison. We can see that the sparsity and representation error resulting from SRDS method is very close to that using OSP.Table 5
Subspace dimension (J/I) in intermodel experiments.
X axisY axisZ axisEx1 (J/I)1/361/371/37Ex2 (J/I)2/872/862/87Ex3 (J/I)3/792/883/83Table 6
Sparsity (μ/σ) and accuracy evaluation for intermodel experiment with kidneys.
X (μ/σ)Y (μ/σ)Z (μ/σ)EOF (%)Haus (mm)Ex1 train (SRDS)34.7/5.535.6/5.536.4/1.60.320.69Ex1 test (SRDS)35.7/1.336.7/0.936.7/1.82.012.08Ex1 test (OSP)36.7/2.436.0/2.435.9/2.61.941.94Ex2 train (SRDS)42.8/4.542.8/2.243.5/3.70.260.59Ex2 test (SRDS)47.0/3.144.5/1.846.3/2.41.021.27Ex2 test (OSP)44.0/2.241.9/12.145.9/2.00.951.22Ex3 train (SRDS)37.1/9.143.3/2.240.8/1.30.250.59Ex3 test (SRDS)41.4/6.545.1/1.740.4/1.30.901.18Ex3 test (OSP)41.2/8.442.7/9.739.3/10.60.831.16Figure18 shows the representation accuracy using SRDS algorithm in training and testing sets for the three tests. In general, the error in testing set is larger than that in the training set. Particularly, as for EOF evaluation, “Ex1” leads to the largest EOF error relative to “Ex2” and “Ex3”, but the average error rate is still as low as 0.3% for training set and 2.0% for testing set. Table 7 lists the specific Hausdorff measurements corresponding to boxplots in Figures 18(c) and 18(d), including minimum, 95th percentile and mean. We can see that the 95th percentile Hausdorff distance across all experiments is below 3 mm, and the mean is belong 2.1 mm. Comparing those error levels with the intramodel test, one can see that the homology existing among the training and testing deformations contributes to better dictionary generalization and, thus, leads to higher representation accuracy.Table 7
Hausdorff distance for intermodel experiment with kidneys.
Min (mm)95th (mm)Mean (mm)Ex1 train0.231.140.69±0.42Ex1 test0.332.812.08±0.61Ex2 train0.301.030.59±0.23Ex2 test0.491.991.27±0.49Ex3 train0.231.030.59±0.39Ex3 test0.541.741.18±0.44Boxplots of representation error in intermodel experiments: (a) EOF of training set, (b) EOF of testing set, (c) Hausdorff distance of training set, (d) Hausdorff distance of testing set.
(a)(b)(c)(d)Figures19, 20 and 21 show the color-coded error fields of a typical representation and at the maximum 90% level for the three intermodel experiments. In either training set or testing set, it is observed that large errors are mostly distributed around the edge area where local details are rich. Consistent with the previous EOF and Hausdorff distance measurements, the color diffusion in Figures 19–21 indicates that errors in testing set is larger than that in training set and “Ex1” generates relatively larger error comparing to “Ex2” or “Ex3.”Representation error in intermodel Ex1: (a) color-coded scales, (b) example of color-coded point error in training set, (c) 90% color-coded point error in training set, (d) example of color-coded point error in testing set, and (e) 90% color-coded point error of testing set.
(a)(b)(c)(d)(e)Representation error in intermodel Ex2: (a) color-coded scales, (b) example of color-coded point error in training set, (c) 90% color-coded point error in training set, (d) example of color-coded point error in testing set, and (e) 90% color-coded point error of testing set.
(a)(b)(c)(d)(e)Representation error in intermodel Ex3: (a) color-coded scales, (b) example of color-coded point error in training set, (c) 90% color-coded point error in training set, (d) example of color-coded point error in testing set, and (e) 90% color-coded point error of testing set.
(a)(b)(c)(d)(e)


### 3.3. In Vivo Experiment Using MRI

The proposed approach is also tested over the in vivo cardiac MR images [37], consisting of automatically segmented images from volumetric MRI scans of a diastole-systole-diastole cycle. For each patient, there are around 22 phases in a cardiac cycle. Surface correspondence of LV shapes within and across patients are accomplished using the approach described in Section 2.5. Since the generated surfaces from automatic segmentation software are quite rough, we use the spherical harmonic representation as a filter to smooth out those surface noises and then apply the smoothed surfaces as training and testing data. Therefore, the demonstrated error in this section is relative to the SHD surfaces, not to the original raw surfaces. The iter-patient results are reported as follows.Similar to the ex vivo test, we use the segmented left ventricles (LV) of 2 different patients plus an initial LV surface for the third patient as training data, and the remaining LV shapes in a beating cycle of the third patient are used to test the generalization of the identified subspaces. The formulated tests are noted as “Ex1,” “Ex2,” and “Ex3.” Table8 lists the sparsity test results of the three cross validations for both training and testing sets. We can see that the sparsity in the training set is close to that in the testing set, but the former achieves much higher accuracy. This is because that the identified subspaces generalize perfectly for those elected atoms among the training set after spherical harmonic smoothing. Consistent with the previous experiments, the representation of testing surfaces using SRDS is also compared with that using repursuing OSP. According to the results, SRDS achieves performance slightly worse than but close to that of OSP. However, as demonstrated in Section 3.4, without relying on the structured dictionary learned from the training population, OSP is a computational expensive task, since for each testing surface, it requires to research for atoms from the training set to achieve sparse representation.Table 8
Sparsity (μ/σ) along X, Y, Z axis and accuracy of in vivo LV tests.
X (μ/σ)Y (μ/σ)Z (μ/σ)EOF (%)Haus (mm)Ex1 train (SRDS)37.0/0.035.4/0.531.7/9.00.130.08Ex1 test (SRDS)37.0/0.036.0/0.034.6/1.53.211.07Ex1 test (OSP)33.3/10.536.8/2.633.8/9.22.920.98Ex2 train (SRDS)35.4/7.432.8/14.132.4/3.30.170.10Ex2 test (SRDS)38.0/0.037.0/0.034.0/0.02.311.15Ex2 test (OSP)32.1/6.324.0/10.532.3/6.51.931.07Ex3 train (SRDS)36.7/3.238.8/1.538.4/0.50.120.07Ex3 test (SRDS)39.6/1.240.0/0.038.8/0.42.331.09Ex3 test (OSP)32.7/11.533.4/11.330.1/12.21.980.99Figure22 provides boxplots for the representation accuracy of the testing set in terms of EOF and Hausdorff distance. Table 9 provides the minimum, 95th percentile, and mean Hausdorff measurements corresponding to Figure 22(b). In coincidence, “Ex1” leads to slightly larger errors than the other two tests, with average EOF of 3.2% (“Ex1”), and mean Hausdorff distance of 1.67±0.39 mm. The 95th percentile Hausdorff distance across all experiments is below 2.2 mm. Figures 23(b) and 23(c) show the color-coded error field of a typical representation and at the maximum 90% level for the testing set in one interpatient experiment (c). As indicated by the color distribution, majority of the point errors are below 0.9 mm. Particularly, in the 90th percentile evaluation in Figure 23, only 3% of all the point-wise errors are above 0.9 mm.Table 9
Hausdorff distance for in vivo experiment with LV.
Min (mm)95th (mm)Mean (mm)Ex1 train00.500.08±0.18Ex1 test0.992.161.67±0.39Ex2 train00.630.10±0.24Ex2 test0.751.371.14±0.18Ex3 train00.600.07±0.21Ex3 test0.721.461.09±0.24Boxplots of representation error in LV interpatient experiments: (a) EOF, (b) Hausdorff distance.
(a)(b)Representation error in in vivo experiment: (a) color-coded scales, (b) example representation error in testing set, (c) 90% representation error of testing set.
(a)(b)(c)


### 3.4. Efficiency

To examine the efficiency of the proposed SRDS method quantitatively, the computational time to represent each surface in the testing set using SRDS method is compared with that resulting from OSP repursuing approach for the above five organs. The results are summarized in Table10, including training set size K, maximum SH level L, average time (in seconds) required with SRDS (t1) and OSP (t2), respectively, and the ratio between the two. As shown in Table 10, the time consumption for seeking sparse representation of the testing surfaces using the SRDS is at least 10 times lower than that using the original OSP method which does not rely on the dictionary structure learned from the training data set. The advantage is more pronounced when the training data size K or the SH level L is large. For example, in the brain model, the high SH level L leads to substantial computational delay during the search for proper atoms for representation, such that the SRDS achieves a speed orders of magnitude faster than the OSP method without training. On the other hand, for the case of gallbladder, the large training size also increases the time used by repursuing OSP, so it runs 65 times slower than SRDS.Table 10
Computational time of SRDS and OSP.
TrainingKSH levelLSRDS (sec)t1OSP (sec)t2t2/t1Brain35800.351.2170.7GBL350250.639.165.2Bladder74302.557.222.9Kidney52200.89.211.5LV51250.56.913.8To summarize, considering the test results for sparsity, accuracy, and efficiency given in this experiment section, we can see that the proposed SRDS method achieves sparse surface representation with high computational efficiency and accuracy.


## 3.1. Experiment with FEM Data

Three representative organs are employed in this FEM experiment: 3D cortical mesh as an example of complicated shapes, gallbladder as an instance with geometrically simple shape, and bladder consisting of both interior and exterior walls.


### 3.1.1. Computer Model Setup

The initial 3D models of different organs are fed into a FEM-based surgical simulation tool to generate deformation data for testing. For instance, Figure6 demonstrates two examples of shape distortions due to the endoscope poking and grasping one side of the gallbladder.Gallbladder deformation generated by surgical simulator: (a) and (b) are distortions caused when the endoscope pokes and grasps the gallbladder.
(a)(b)Table1 lists the FEM experimental setup of the three organs including number of vertices N, SH level L, number of deformations for training K, and testing M. “GBL” stands for gallbladder in all the tables. The maximum SH level used for brain model is chosen according to [9], and the levels for gallbladder and bladder are determined when the SHD representation error is below 0.1% (EOF). The complex brain structure requires more vertices and higher SH level for surface representation to achieve sufficient accuracy. To evaluate the representation precision qualitatively, an evaluation parameter EOF is defined as the normalized Euclidean distance between the original surface and the reconstructed surface(19)EOF=‖x̂k-xk‖2‖xk‖2.
All surfaces are centered to the origin of the coordinate system so that EOF will not be heavily affected by the denominator.Table 1
FEM Model setup.
VerticesNSH levelLTrainingKTestingMBrain40962803535GBL303825250114BladderN1 = 4434,307446N2 = 4274


### 3.1.2. Results

With the FEM data, the proposed SRDS algorithm is evaluated from three perspectives: (1) subspace dimensionality, (2) sparsity and accuracy of representations, and (3) the effect of subspace pursuit thresholdɛ and coefficient truncation threshold δ on the performance. The sparsity is defined as the l0 norm of the coefficient vector c⃗̃m.A. Training Results
During training stage, we setɛ=0.005 for subspace detection, η=0.01 for clustering, Emax=50 as the maximum iteration times, and δ=0.005 for coefficient truncation. Subspaces on X, Y, and Z axis are identified separately. Table 2 shows that the subspace number J and dimensions of resulting dictionary (dim(G)=I) are markedly small relative to N or L2 in all three tests. We notice that the subspace dimensions of brain are relatively smaller than the other two. This is because of smaller training data size and minor extent of deformation considered in the brain experiment, which results in smaller dictionary size to capture the deformation features.Table 2
Dimension of dictionary (J/I).
BrainGBLBladderSubspace onX (J/I)2/132/391/37Subspace onY (J/I)2/62/481/30Subspace onZ (J/I)1/31/421/34B. Sparsity and Accuracy Evaluation
Sparsity is examined in terms of (μ/σ), where μ is the average l0 norm of the coefficient vector c⃗̃k (training) or c⃗̃m (testing) and σ is the corresponding standard deviation. To verify that whether our method achieves equivalent sparsity and precision when applying the structure of the dictionary, we also test the case without relying on any structure learned from training set, during which sparse representation of each deformation in the testing set is repursued from the training set using OSP approach. In the following tables, we use “OSP” to refer to the results obtained using such a repursuing process. Table 3 summarizes the sparsity of the SRDS representation of three organs for both training and testing set. It illustrates that the number of atoms needed for representing the complex deformations is much smaller than the dimension of spherical harmonic vectors ((L+1)2), and particularly the sparsity and accuracy via SRDS is very close to that from OSP repursuit for the testing deformations, which indicates the good generalization of the structured dictionary.
The reconstruction error in terms of EOF is further compared with that from standard SHD method, as shown in Figure7. In general, the accuracy of SRDS is equivalent to that of SHD method. Specifically, it shows that the SRDS method achieves average EOF of 1.32% (brain) and 0.14% (gallbladder) versus 1.29% (brain) and 0.13% (gallbladder) with SHD method. For bladder model with deformations on multiple layers, the overall representation error with SRDS is 0.07%, very close to 0.06% with SHD. Figures 11 and 12 show typical reconstructed deformations of the testing data for the three organs with SHD and SRDS methods. Figures 12(e) and 12(f) demonstrate the interior and exterior representation of the bladder at a same time instance. From those results, we can see that the SRDS algorithm achieves the accuracy equivalent to complex mathematical modeling techniques while significantly lowers the representation dimensionality.Table 3
Sparsity (μ/σ) and accuracy evaluation.
X (μ/σ)Y (μ/σ)Z (μ/σ)EOF (%)Brain train (SRDS)3.1/0.23.4/0.93.0/0.21.31Brain test (SRDS)3.5/0.53.9/0.43.0/0.01.32Brain test (OSP)4.1/0.64.5/0.53.9/0.81.30GBL train (SRDS)25.8/12.24.8/1.832.3/7.30.13GBL test (SRDS)33.2/1.044.1/1.140.8/1.30.15GBL test (OSP)44.4/4.048.3/2.343.9/4.50.13Bladder train (SRDS)24.1/3.720.8/2.924.8/4.00.076Bladder test (SRDS)24.6/2.721.0/2.525.6/3.20.073Bladder test (OSP)22.5/4.316.3/6.523.0/5.20.070EOF of brain, gallbladder, and bladder reconstruction with SHD and SRDS methods. The left pair is for training set; the right pair is for testing set.
(a)(b)(c)C. Effect ofɛ
The performance of SRDS algorithm is examined as the subspace pursuit thresholdɛ varies. Specifically, we study the effect of ɛ on the dimensionality (I) of the structured dictionary G, sparsity and accuracy of the surface representation. Figure 8 shows how the subspace dimensions on three axis change during the training stage as ɛ increases from 0.001 to 0.01. Figure 9 displays the influence of ɛ on the average sparsity μ of the surface representation results for both training and testing data sets. In general, smaller ɛ leads to larger subspace size and less description sparsity, since lower ɛ usually leads to more recruited atoms to meet the desired representation accuracy. Therefore, there is a tradeoff between representation accuracy and desired sparsity. Figure 10 reveals the representation EOF as a function of ɛ. Not surprisingly, the reconstruction error is increased as ɛ becomes larger. An empirical point can be chosen according to the training curve when space dimension I expands significantly but only trivial EOF improvement is gained, that is, ɛ=0.005  is a preferred value in this test according to Figure 10.Figure 8
DimensionI of trained dictionary G decreases on all three axis as ɛ increases.Sparsity (μ) of surface representation changes with different ɛ on three axis: (a) is for training set; (b) is for testing set.
(a)(b)Figure 10
EOF of gallbladder reconstruction increases asɛ gets larger for both training and testing sets.Reconstructed brain deformations using SHD approach: (a) is the initial brain shape and the circle marks one typical area under deformation, (b)–(e) are the reconstructed brain deformations.
(a)(b)(c)(d)(e)(f)Reconstructed brain deformations using SRDS approach: (a) is the initial brain shape, (b)–(e) are the reconstructed brain deformations.
(a)(b)(c)(d)(e)(f)D. Effect ofδ
The influence of coefficient truncation thresholdδ on the performance of SRDS algorithm is also tested while δ is varied among [00.00010.00050.0010.0050.010.050.080.1]. Figure 13 shows the effect of δ on the average sparsity μ of the surface representation results. We can see that, as the truncation threshold δ enlarges, the sparsity of the representation is increased for both training and testing data sets at the price of decreased representation error as shown in Figure 14, so there is tradeoff between sparsity and accuracy. Empirically, one can choose the δ value when the representation precision remarkably deteriorates while the sparsity is still increasing. Therefore, according to Figures 13 and 14, an appropriate value for δ is between 0.005 and 0.01.Sparsity (μ) of surface representation changes with different δ on three axis: (a) is for training set; (b) is for testing set.
(a)(b)Figure 14
EOF of gallbladder reconstruction increases asδ gets larger for both training and testing sets.


## 3.1.1. Computer Model Setup

The initial 3D models of different organs are fed into a FEM-based surgical simulation tool to generate deformation data for testing. For instance, Figure6 demonstrates two examples of shape distortions due to the endoscope poking and grasping one side of the gallbladder.Gallbladder deformation generated by surgical simulator: (a) and (b) are distortions caused when the endoscope pokes and grasps the gallbladder.
(a)(b)Table1 lists the FEM experimental setup of the three organs including number of vertices N, SH level L, number of deformations for training K, and testing M. “GBL” stands for gallbladder in all the tables. The maximum SH level used for brain model is chosen according to [9], and the levels for gallbladder and bladder are determined when the SHD representation error is below 0.1% (EOF). The complex brain structure requires more vertices and higher SH level for surface representation to achieve sufficient accuracy. To evaluate the representation precision qualitatively, an evaluation parameter EOF is defined as the normalized Euclidean distance between the original surface and the reconstructed surface(19)EOF=‖x̂k-xk‖2‖xk‖2.
All surfaces are centered to the origin of the coordinate system so that EOF will not be heavily affected by the denominator.Table 1
FEM Model setup.
VerticesNSH levelLTrainingKTestingMBrain40962803535GBL303825250114BladderN1 = 4434,307446N2 = 4274


## 3.1.2. Results

With the FEM data, the proposed SRDS algorithm is evaluated from three perspectives: (1) subspace dimensionality, (2) sparsity and accuracy of representations, and (3) the effect of subspace pursuit thresholdɛ and coefficient truncation threshold δ on the performance. The sparsity is defined as the l0 norm of the coefficient vector c⃗̃m.A. Training Results
During training stage, we setɛ=0.005 for subspace detection, η=0.01 for clustering, Emax=50 as the maximum iteration times, and δ=0.005 for coefficient truncation. Subspaces on X, Y, and Z axis are identified separately. Table 2 shows that the subspace number J and dimensions of resulting dictionary (dim(G)=I) are markedly small relative to N or L2 in all three tests. We notice that the subspace dimensions of brain are relatively smaller than the other two. This is because of smaller training data size and minor extent of deformation considered in the brain experiment, which results in smaller dictionary size to capture the deformation features.Table 2
Dimension of dictionary (J/I).
BrainGBLBladderSubspace onX (J/I)2/132/391/37Subspace onY (J/I)2/62/481/30Subspace onZ (J/I)1/31/421/34B. Sparsity and Accuracy Evaluation
Sparsity is examined in terms of (μ/σ), where μ is the average l0 norm of the coefficient vector c⃗̃k (training) or c⃗̃m (testing) and σ is the corresponding standard deviation. To verify that whether our method achieves equivalent sparsity and precision when applying the structure of the dictionary, we also test the case without relying on any structure learned from training set, during which sparse representation of each deformation in the testing set is repursued from the training set using OSP approach. In the following tables, we use “OSP” to refer to the results obtained using such a repursuing process. Table 3 summarizes the sparsity of the SRDS representation of three organs for both training and testing set. It illustrates that the number of atoms needed for representing the complex deformations is much smaller than the dimension of spherical harmonic vectors ((L+1)2), and particularly the sparsity and accuracy via SRDS is very close to that from OSP repursuit for the testing deformations, which indicates the good generalization of the structured dictionary.
The reconstruction error in terms of EOF is further compared with that from standard SHD method, as shown in Figure7. In general, the accuracy of SRDS is equivalent to that of SHD method. Specifically, it shows that the SRDS method achieves average EOF of 1.32% (brain) and 0.14% (gallbladder) versus 1.29% (brain) and 0.13% (gallbladder) with SHD method. For bladder model with deformations on multiple layers, the overall representation error with SRDS is 0.07%, very close to 0.06% with SHD. Figures 11 and 12 show typical reconstructed deformations of the testing data for the three organs with SHD and SRDS methods. Figures 12(e) and 12(f) demonstrate the interior and exterior representation of the bladder at a same time instance. From those results, we can see that the SRDS algorithm achieves the accuracy equivalent to complex mathematical modeling techniques while significantly lowers the representation dimensionality.Table 3
Sparsity (μ/σ) and accuracy evaluation.
X (μ/σ)Y (μ/σ)Z (μ/σ)EOF (%)Brain train (SRDS)3.1/0.23.4/0.93.0/0.21.31Brain test (SRDS)3.5/0.53.9/0.43.0/0.01.32Brain test (OSP)4.1/0.64.5/0.53.9/0.81.30GBL train (SRDS)25.8/12.24.8/1.832.3/7.30.13GBL test (SRDS)33.2/1.044.1/1.140.8/1.30.15GBL test (OSP)44.4/4.048.3/2.343.9/4.50.13Bladder train (SRDS)24.1/3.720.8/2.924.8/4.00.076Bladder test (SRDS)24.6/2.721.0/2.525.6/3.20.073Bladder test (OSP)22.5/4.316.3/6.523.0/5.20.070EOF of brain, gallbladder, and bladder reconstruction with SHD and SRDS methods. The left pair is for training set; the right pair is for testing set.
(a)(b)(c)C. Effect ofɛ
The performance of SRDS algorithm is examined as the subspace pursuit thresholdɛ varies. Specifically, we study the effect of ɛ on the dimensionality (I) of the structured dictionary G, sparsity and accuracy of the surface representation. Figure 8 shows how the subspace dimensions on three axis change during the training stage as ɛ increases from 0.001 to 0.01. Figure 9 displays the influence of ɛ on the average sparsity μ of the surface representation results for both training and testing data sets. In general, smaller ɛ leads to larger subspace size and less description sparsity, since lower ɛ usually leads to more recruited atoms to meet the desired representation accuracy. Therefore, there is a tradeoff between representation accuracy and desired sparsity. Figure 10 reveals the representation EOF as a function of ɛ. Not surprisingly, the reconstruction error is increased as ɛ becomes larger. An empirical point can be chosen according to the training curve when space dimension I expands significantly but only trivial EOF improvement is gained, that is, ɛ=0.005  is a preferred value in this test according to Figure 10.Figure 8
DimensionI of trained dictionary G decreases on all three axis as ɛ increases.Sparsity (μ) of surface representation changes with different ɛ on three axis: (a) is for training set; (b) is for testing set.
(a)(b)Figure 10
EOF of gallbladder reconstruction increases asɛ gets larger for both training and testing sets.Reconstructed brain deformations using SHD approach: (a) is the initial brain shape and the circle marks one typical area under deformation, (b)–(e) are the reconstructed brain deformations.
(a)(b)(c)(d)(e)(f)Reconstructed brain deformations using SRDS approach: (a) is the initial brain shape, (b)–(e) are the reconstructed brain deformations.
(a)(b)(c)(d)(e)(f)D. Effect ofδ
The influence of coefficient truncation thresholdδ on the performance of SRDS algorithm is also tested while δ is varied among [00.00010.00050.0010.0050.010.050.080.1]. Figure 13 shows the effect of δ on the average sparsity μ of the surface representation results. We can see that, as the truncation threshold δ enlarges, the sparsity of the representation is increased for both training and testing data sets at the price of decreased representation error as shown in Figure 14, so there is tradeoff between sparsity and accuracy. Empirically, one can choose the δ value when the representation precision remarkably deteriorates while the sparsity is still increasing. Therefore, according to Figures 13 and 14, an appropriate value for δ is between 0.005 and 0.01.Sparsity (μ) of surface representation changes with different δ on three axis: (a) is for training set; (b) is for testing set.
(a)(b)Figure 14
EOF of gallbladder reconstruction increases asδ gets larger for both training and testing sets.


## 3.2. Ex vivo Experiment Using MRI

To evaluate the proposed algorithm in real applications, an ex vivo experiment using three porcine kidneys were conducted at the Center for Interdisciplinary Applications in Magnetic Resonance (CIA-MR) of University of Minnesota. Deformations imposed to each kidney were controlled and maintained still during imaging by a customized nonmagnetic mechanical device as shown in Figure15. Each deformed kidney shape was scanned in 3D MRI mode with spatial resolution of 1.2 mm to generate both training set and testing set. The SH degree L of the organ representation is set to be 20, and each 3D kidney mesh after surface correspondence has N=4002 vertices. Different from computer-generated deformations where surface correspondence is intrinsically established, the shapes from MRI scans are rendered independently, so the method described in Section 2.5 is applied to achieve point-wise correspondence.Images of three porcine kidneys for experiment.
(a)(b)(c)(d)(e)Both intramodel and intermodel experiments are conducted. The former uses training and testing deformations from the same kidney; the later utilizes two out of the three kidneys for training and the third one for testing in a cross-evaluation fashion. Besides sparsity and EOF, Hausdorff distance between the represented shape and corresponding MRI surface is also examined as a physical measurement of error. The Hausdorff distance between surfacex and x′ is defined as (20)d(x,x′)=max⁡p∈x⁡d(p,x′),
where d(p,x′) is defined as the distance between a point p on surface x and the closest point on surface x′, that is,(21)d(p,x′)=min⁡p′∈x′⁡‖p,p′‖2
with ∥·∥2 denoting the Euclidean norm.


### 3.2.1. Intramodel Test

In the intramodel experiment, 31 deformations of the same kidney were generated and scanned by the MRI machine, among which 20 frames were randomly selected as training set, and the other 11 were applied for testing the generalization of the learned subspaces.Table4 shows the trained subspace dimensions (J as number of subspace, I as dictionary size of G), the sparsity of the descriptors in each axis for both training set and testing set, and the corresponding errors in terms of EOF and Hausdorff distance. Similar to the FEM experiment, the sparsity is also evaluated with OSP repursuit process in the testing set for comparison. The table shows that the sparsity and the accuracy achieved with SRDS is very close to that from OSP repursuing process. However, the SRDS method features delay-free surface representation by applying the structure in the identified dictionary. Further results about computational efficiency are shown in Section 3.4. One may notice that the size of training data in the MRI experiment is smaller than that in FEM test due to the less availability of 3D MRI images. As a rule of thumb, larger training set carries richer deformation information and thus leads to better generalization of the dictionary. However, given the size of training data and extent of deformation involved in the ex vivo experiment, high representation precision is still achieved.Table 4
Subspace dimension and sparsity (μ/σ) for intramodel experiment.
X (μ/σ)Y (μ/σ)Z (μ/σ)EOF (%)Haus (mm)J/I1/172/312/31Train (SRDS)16.9/0.415.5/1.515.8/1.50.240.55Test (SRDS)17.0/0.016.9/0.316.7/0.90.640.87Test (OSP)16.6/1.516.3/0.917.1/0.80.600.85Figure16 illustrates the accuracy of the surface representation in the intramodel test. The average EOF in Figure 16(a) for training set is 0.24% and 0.64% for testing set, with maximum rate less than 1%. Further, error as Hausdorff distance (shown in Figure 16(b)) is 0.55±0.23 mm with 95th percentile error of 0.86 mm for the training set, and 0.87±0.10 mm with 95th percentile error of 0.96 mm for the testing set. This intramodel experiment demonstrated that the SRDS algorithm identifies subspaces generalizable enough to accurately represent deformations beyond the training set for the same object.Boxplots of representation error in intramodel experiment: (a) EOF of training and testing set, (b) Hausdorff distance of training and testing sets.
(a)(b)Figure17 visualizes the color-coded error distribution at all vertices on the represented surface with SRDS relative to the actual MRI scans. Figure 17(a) illustrates the error range for different colors. Figures 17(b) and 17(d) show the error distribution for a typical reconstruction in the training and testing set, respectively. Figures 17(c) and 17(e) show maximum 90% level reconstruction errors, that is, 90% of all deformations in the training or testing set have representation point errors less than the values shown in the figures. Consistent with the EOF and Hausdorff distance results, the color diffusion in Figure 17 indicates that the precision in the testing group is relatively lower than that in the training group. However, among all the pixel-wise errors shown in Figure 17(e), less than 3% of all the surface points have error distance larger than 0.5 mm.Representation error in intramodel experiment: (a) color-coded scales, (b) example representation error in training set, (c) 90% representation error in training set, (d) example of representation error in testing set, and (e) 90% representation error of testing set.
(a)(b)(c)(d)(e)


### 3.2.2. Intermodel Test

Three intermodel experiments are performed to further validate the proposed SRDS method applied to organs from different subjects. In the following context, “Ex1” stands for the experiment training on Kidney 2 and 3 plus one initial shape of Kidney 1 while testing on deformations of Kidney 1, and the like for “Ex2” and “Ex3”. In each experiment, both sparsity and accuracy are examined.The number of subspaces (J) and dimensions (I) of the identified dictionary are listed in Table 5. The training results vary among the three experiments but all features low subspace dimensions. Table 6 shows the sparsity of the intermodel experiments using the SRDS algorithm, and the error level is evaluated in terms of EOF and Hausdorff distance. Each testing deformation is also sparsely retrained using OSP for comparison. We can see that the sparsity and representation error resulting from SRDS method is very close to that using OSP.Table 5
Subspace dimension (J/I) in intermodel experiments.
X axisY axisZ axisEx1 (J/I)1/361/371/37Ex2 (J/I)2/872/862/87Ex3 (J/I)3/792/883/83Table 6
Sparsity (μ/σ) and accuracy evaluation for intermodel experiment with kidneys.
X (μ/σ)Y (μ/σ)Z (μ/σ)EOF (%)Haus (mm)Ex1 train (SRDS)34.7/5.535.6/5.536.4/1.60.320.69Ex1 test (SRDS)35.7/1.336.7/0.936.7/1.82.012.08Ex1 test (OSP)36.7/2.436.0/2.435.9/2.61.941.94Ex2 train (SRDS)42.8/4.542.8/2.243.5/3.70.260.59Ex2 test (SRDS)47.0/3.144.5/1.846.3/2.41.021.27Ex2 test (OSP)44.0/2.241.9/12.145.9/2.00.951.22Ex3 train (SRDS)37.1/9.143.3/2.240.8/1.30.250.59Ex3 test (SRDS)41.4/6.545.1/1.740.4/1.30.901.18Ex3 test (OSP)41.2/8.442.7/9.739.3/10.60.831.16Figure18 shows the representation accuracy using SRDS algorithm in training and testing sets for the three tests. In general, the error in testing set is larger than that in the training set. Particularly, as for EOF evaluation, “Ex1” leads to the largest EOF error relative to “Ex2” and “Ex3”, but the average error rate is still as low as 0.3% for training set and 2.0% for testing set. Table 7 lists the specific Hausdorff measurements corresponding to boxplots in Figures 18(c) and 18(d), including minimum, 95th percentile and mean. We can see that the 95th percentile Hausdorff distance across all experiments is below 3 mm, and the mean is belong 2.1 mm. Comparing those error levels with the intramodel test, one can see that the homology existing among the training and testing deformations contributes to better dictionary generalization and, thus, leads to higher representation accuracy.Table 7
Hausdorff distance for intermodel experiment with kidneys.
Min (mm)95th (mm)Mean (mm)Ex1 train0.231.140.69±0.42Ex1 test0.332.812.08±0.61Ex2 train0.301.030.59±0.23Ex2 test0.491.991.27±0.49Ex3 train0.231.030.59±0.39Ex3 test0.541.741.18±0.44Boxplots of representation error in intermodel experiments: (a) EOF of training set, (b) EOF of testing set, (c) Hausdorff distance of training set, (d) Hausdorff distance of testing set.
(a)(b)(c)(d)Figures19, 20 and 21 show the color-coded error fields of a typical representation and at the maximum 90% level for the three intermodel experiments. In either training set or testing set, it is observed that large errors are mostly distributed around the edge area where local details are rich. Consistent with the previous EOF and Hausdorff distance measurements, the color diffusion in Figures 19–21 indicates that errors in testing set is larger than that in training set and “Ex1” generates relatively larger error comparing to “Ex2” or “Ex3.”Representation error in intermodel Ex1: (a) color-coded scales, (b) example of color-coded point error in training set, (c) 90% color-coded point error in training set, (d) example of color-coded point error in testing set, and (e) 90% color-coded point error of testing set.
(a)(b)(c)(d)(e)Representation error in intermodel Ex2: (a) color-coded scales, (b) example of color-coded point error in training set, (c) 90% color-coded point error in training set, (d) example of color-coded point error in testing set, and (e) 90% color-coded point error of testing set.
(a)(b)(c)(d)(e)Representation error in intermodel Ex3: (a) color-coded scales, (b) example of color-coded point error in training set, (c) 90% color-coded point error in training set, (d) example of color-coded point error in testing set, and (e) 90% color-coded point error of testing set.
(a)(b)(c)(d)(e)


## 3.2.1. Intramodel Test

In the intramodel experiment, 31 deformations of the same kidney were generated and scanned by the MRI machine, among which 20 frames were randomly selected as training set, and the other 11 were applied for testing the generalization of the learned subspaces.Table4 shows the trained subspace dimensions (J as number of subspace, I as dictionary size of G), the sparsity of the descriptors in each axis for both training set and testing set, and the corresponding errors in terms of EOF and Hausdorff distance. Similar to the FEM experiment, the sparsity is also evaluated with OSP repursuit process in the testing set for comparison. The table shows that the sparsity and the accuracy achieved with SRDS is very close to that from OSP repursuing process. However, the SRDS method features delay-free surface representation by applying the structure in the identified dictionary. Further results about computational efficiency are shown in Section 3.4. One may notice that the size of training data in the MRI experiment is smaller than that in FEM test due to the less availability of 3D MRI images. As a rule of thumb, larger training set carries richer deformation information and thus leads to better generalization of the dictionary. However, given the size of training data and extent of deformation involved in the ex vivo experiment, high representation precision is still achieved.Table 4
Subspace dimension and sparsity (μ/σ) for intramodel experiment.
X (μ/σ)Y (μ/σ)Z (μ/σ)EOF (%)Haus (mm)J/I1/172/312/31Train (SRDS)16.9/0.415.5/1.515.8/1.50.240.55Test (SRDS)17.0/0.016.9/0.316.7/0.90.640.87Test (OSP)16.6/1.516.3/0.917.1/0.80.600.85Figure16 illustrates the accuracy of the surface representation in the intramodel test. The average EOF in Figure 16(a) for training set is 0.24% and 0.64% for testing set, with maximum rate less than 1%. Further, error as Hausdorff distance (shown in Figure 16(b)) is 0.55±0.23 mm with 95th percentile error of 0.86 mm for the training set, and 0.87±0.10 mm with 95th percentile error of 0.96 mm for the testing set. This intramodel experiment demonstrated that the SRDS algorithm identifies subspaces generalizable enough to accurately represent deformations beyond the training set for the same object.Boxplots of representation error in intramodel experiment: (a) EOF of training and testing set, (b) Hausdorff distance of training and testing sets.
(a)(b)Figure17 visualizes the color-coded error distribution at all vertices on the represented surface with SRDS relative to the actual MRI scans. Figure 17(a) illustrates the error range for different colors. Figures 17(b) and 17(d) show the error distribution for a typical reconstruction in the training and testing set, respectively. Figures 17(c) and 17(e) show maximum 90% level reconstruction errors, that is, 90% of all deformations in the training or testing set have representation point errors less than the values shown in the figures. Consistent with the EOF and Hausdorff distance results, the color diffusion in Figure 17 indicates that the precision in the testing group is relatively lower than that in the training group. However, among all the pixel-wise errors shown in Figure 17(e), less than 3% of all the surface points have error distance larger than 0.5 mm.Representation error in intramodel experiment: (a) color-coded scales, (b) example representation error in training set, (c) 90% representation error in training set, (d) example of representation error in testing set, and (e) 90% representation error of testing set.
(a)(b)(c)(d)(e)


## 3.2.2. Intermodel Test

Three intermodel experiments are performed to further validate the proposed SRDS method applied to organs from different subjects. In the following context, “Ex1” stands for the experiment training on Kidney 2 and 3 plus one initial shape of Kidney 1 while testing on deformations of Kidney 1, and the like for “Ex2” and “Ex3”. In each experiment, both sparsity and accuracy are examined.The number of subspaces (J) and dimensions (I) of the identified dictionary are listed in Table 5. The training results vary among the three experiments but all features low subspace dimensions. Table 6 shows the sparsity of the intermodel experiments using the SRDS algorithm, and the error level is evaluated in terms of EOF and Hausdorff distance. Each testing deformation is also sparsely retrained using OSP for comparison. We can see that the sparsity and representation error resulting from SRDS method is very close to that using OSP.Table 5
Subspace dimension (J/I) in intermodel experiments.
X axisY axisZ axisEx1 (J/I)1/361/371/37Ex2 (J/I)2/872/862/87Ex3 (J/I)3/792/883/83Table 6
Sparsity (μ/σ) and accuracy evaluation for intermodel experiment with kidneys.
X (μ/σ)Y (μ/σ)Z (μ/σ)EOF (%)Haus (mm)Ex1 train (SRDS)34.7/5.535.6/5.536.4/1.60.320.69Ex1 test (SRDS)35.7/1.336.7/0.936.7/1.82.012.08Ex1 test (OSP)36.7/2.436.0/2.435.9/2.61.941.94Ex2 train (SRDS)42.8/4.542.8/2.243.5/3.70.260.59Ex2 test (SRDS)47.0/3.144.5/1.846.3/2.41.021.27Ex2 test (OSP)44.0/2.241.9/12.145.9/2.00.951.22Ex3 train (SRDS)37.1/9.143.3/2.240.8/1.30.250.59Ex3 test (SRDS)41.4/6.545.1/1.740.4/1.30.901.18Ex3 test (OSP)41.2/8.442.7/9.739.3/10.60.831.16Figure18 shows the representation accuracy using SRDS algorithm in training and testing sets for the three tests. In general, the error in testing set is larger than that in the training set. Particularly, as for EOF evaluation, “Ex1” leads to the largest EOF error relative to “Ex2” and “Ex3”, but the average error rate is still as low as 0.3% for training set and 2.0% for testing set. Table 7 lists the specific Hausdorff measurements corresponding to boxplots in Figures 18(c) and 18(d), including minimum, 95th percentile and mean. We can see that the 95th percentile Hausdorff distance across all experiments is below 3 mm, and the mean is belong 2.1 mm. Comparing those error levels with the intramodel test, one can see that the homology existing among the training and testing deformations contributes to better dictionary generalization and, thus, leads to higher representation accuracy.Table 7
Hausdorff distance for intermodel experiment with kidneys.
Min (mm)95th (mm)Mean (mm)Ex1 train0.231.140.69±0.42Ex1 test0.332.812.08±0.61Ex2 train0.301.030.59±0.23Ex2 test0.491.991.27±0.49Ex3 train0.231.030.59±0.39Ex3 test0.541.741.18±0.44Boxplots of representation error in intermodel experiments: (a) EOF of training set, (b) EOF of testing set, (c) Hausdorff distance of training set, (d) Hausdorff distance of testing set.
(a)(b)(c)(d)Figures19, 20 and 21 show the color-coded error fields of a typical representation and at the maximum 90% level for the three intermodel experiments. In either training set or testing set, it is observed that large errors are mostly distributed around the edge area where local details are rich. Consistent with the previous EOF and Hausdorff distance measurements, the color diffusion in Figures 19–21 indicates that errors in testing set is larger than that in training set and “Ex1” generates relatively larger error comparing to “Ex2” or “Ex3.”Representation error in intermodel Ex1: (a) color-coded scales, (b) example of color-coded point error in training set, (c) 90% color-coded point error in training set, (d) example of color-coded point error in testing set, and (e) 90% color-coded point error of testing set.
(a)(b)(c)(d)(e)Representation error in intermodel Ex2: (a) color-coded scales, (b) example of color-coded point error in training set, (c) 90% color-coded point error in training set, (d) example of color-coded point error in testing set, and (e) 90% color-coded point error of testing set.
(a)(b)(c)(d)(e)Representation error in intermodel Ex3: (a) color-coded scales, (b) example of color-coded point error in training set, (c) 90% color-coded point error in training set, (d) example of color-coded point error in testing set, and (e) 90% color-coded point error of testing set.
(a)(b)(c)(d)(e)


## 3.3. In Vivo Experiment Using MRI

The proposed approach is also tested over the in vivo cardiac MR images [37], consisting of automatically segmented images from volumetric MRI scans of a diastole-systole-diastole cycle. For each patient, there are around 22 phases in a cardiac cycle. Surface correspondence of LV shapes within and across patients are accomplished using the approach described in Section 2.5. Since the generated surfaces from automatic segmentation software are quite rough, we use the spherical harmonic representation as a filter to smooth out those surface noises and then apply the smoothed surfaces as training and testing data. Therefore, the demonstrated error in this section is relative to the SHD surfaces, not to the original raw surfaces. The iter-patient results are reported as follows.Similar to the ex vivo test, we use the segmented left ventricles (LV) of 2 different patients plus an initial LV surface for the third patient as training data, and the remaining LV shapes in a beating cycle of the third patient are used to test the generalization of the identified subspaces. The formulated tests are noted as “Ex1,” “Ex2,” and “Ex3.” Table8 lists the sparsity test results of the three cross validations for both training and testing sets. We can see that the sparsity in the training set is close to that in the testing set, but the former achieves much higher accuracy. This is because that the identified subspaces generalize perfectly for those elected atoms among the training set after spherical harmonic smoothing. Consistent with the previous experiments, the representation of testing surfaces using SRDS is also compared with that using repursuing OSP. According to the results, SRDS achieves performance slightly worse than but close to that of OSP. However, as demonstrated in Section 3.4, without relying on the structured dictionary learned from the training population, OSP is a computational expensive task, since for each testing surface, it requires to research for atoms from the training set to achieve sparse representation.Table 8
Sparsity (μ/σ) along X, Y, Z axis and accuracy of in vivo LV tests.
X (μ/σ)Y (μ/σ)Z (μ/σ)EOF (%)Haus (mm)Ex1 train (SRDS)37.0/0.035.4/0.531.7/9.00.130.08Ex1 test (SRDS)37.0/0.036.0/0.034.6/1.53.211.07Ex1 test (OSP)33.3/10.536.8/2.633.8/9.22.920.98Ex2 train (SRDS)35.4/7.432.8/14.132.4/3.30.170.10Ex2 test (SRDS)38.0/0.037.0/0.034.0/0.02.311.15Ex2 test (OSP)32.1/6.324.0/10.532.3/6.51.931.07Ex3 train (SRDS)36.7/3.238.8/1.538.4/0.50.120.07Ex3 test (SRDS)39.6/1.240.0/0.038.8/0.42.331.09Ex3 test (OSP)32.7/11.533.4/11.330.1/12.21.980.99Figure22 provides boxplots for the representation accuracy of the testing set in terms of EOF and Hausdorff distance. Table 9 provides the minimum, 95th percentile, and mean Hausdorff measurements corresponding to Figure 22(b). In coincidence, “Ex1” leads to slightly larger errors than the other two tests, with average EOF of 3.2% (“Ex1”), and mean Hausdorff distance of 1.67±0.39 mm. The 95th percentile Hausdorff distance across all experiments is below 2.2 mm. Figures 23(b) and 23(c) show the color-coded error field of a typical representation and at the maximum 90% level for the testing set in one interpatient experiment (c). As indicated by the color distribution, majority of the point errors are below 0.9 mm. Particularly, in the 90th percentile evaluation in Figure 23, only 3% of all the point-wise errors are above 0.9 mm.Table 9
Hausdorff distance for in vivo experiment with LV.
Min (mm)95th (mm)Mean (mm)Ex1 train00.500.08±0.18Ex1 test0.992.161.67±0.39Ex2 train00.630.10±0.24Ex2 test0.751.371.14±0.18Ex3 train00.600.07±0.21Ex3 test0.721.461.09±0.24Boxplots of representation error in LV interpatient experiments: (a) EOF, (b) Hausdorff distance.
(a)(b)Representation error in in vivo experiment: (a) color-coded scales, (b) example representation error in testing set, (c) 90% representation error of testing set.
(a)(b)(c)


## 3.4. Efficiency

To examine the efficiency of the proposed SRDS method quantitatively, the computational time to represent each surface in the testing set using SRDS method is compared with that resulting from OSP repursuing approach for the above five organs. The results are summarized in Table10, including training set size K, maximum SH level L, average time (in seconds) required with SRDS (t1) and OSP (t2), respectively, and the ratio between the two. As shown in Table 10, the time consumption for seeking sparse representation of the testing surfaces using the SRDS is at least 10 times lower than that using the original OSP method which does not rely on the dictionary structure learned from the training data set. The advantage is more pronounced when the training data size K or the SH level L is large. For example, in the brain model, the high SH level L leads to substantial computational delay during the search for proper atoms for representation, such that the SRDS achieves a speed orders of magnitude faster than the OSP method without training. On the other hand, for the case of gallbladder, the large training size also increases the time used by repursuing OSP, so it runs 65 times slower than SRDS.Table 10
Computational time of SRDS and OSP.
TrainingKSH levelLSRDS (sec)t1OSP (sec)t2t2/t1Brain35800.351.2170.7GBL350250.639.165.2Bladder74302.557.222.9Kidney52200.89.211.5LV51250.56.913.8To summarize, considering the test results for sparsity, accuracy, and efficiency given in this experiment section, we can see that the proposed SRDS method achieves sparse surface representation with high computational efficiency and accuracy.


## 4. Conclusions and Discussion

This paper introduced a new algorithm for block sparse representation of deformable organ surfaces with high accuracy. The proposed SRDS design first identifies the deformation subspaces from the training data set in the transformed spherical harmonic domain and then represents each deformed surface with a block sparse vector in the structured dictionary. SRDS is generalized to applications involving organs with multiple surface layers, such as bladder. The algorithm has been validated with FEM data and real 3D MRI scans under both ex vivo and in vivo conditions. The FEM test results demonstrate that SRDS achieves accuracy matching that of complex mathematical modeling techniques. Further, the maximum representation error in ex vivo experiment is below 1 mm for intramodel test and below 3 mm for intermodel test. For the in vivo experiment, the SRDS achieves an accuracy of better than 2.5 mm.SRDS algorithm has already been used in tracking organ deformations in minimum invasive surgery [2–4]. The structure introduced in the dictionary enables efficient surface recovery from limited samples. In addition, the merits of block sparse surface representation presented here can be applied to various medical organ modeling, shape classification, and similarity retrieval where reduced parameter dimension can potentially speed up the implementations.



---

*Source: 658930-2011-09-19.xml*