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 ) . 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   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 . the later topic of real time tracking an deformable organ with limited access to the organ is explored further in [ 24 ] . 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   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   treats the whole surface as a union of localized patches . global surface representation methods [ 710 ] , particularly , decomposing surfaces into other primitive shapes , are more appropriate for shape analysis and classification due to the lower dimensionality of the parameter space . 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 )   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 . ) shd has been widely used in applications related to surface description , including static modeling of kidney   and brain [ 9 , 10 ] , as well as spatial - temporal modeling of left - ventricular with known motion period . in , hemisphere 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 ) , orthogonal matching pursuit ( omp ) , and orthogonal least square ( ols ) . 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 . the second type is global optimization algorithms , such as basis pursuit ( bp ) , focuss   and iterative thresholding . 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   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 ) , method of optimal directions ( mod ) , maximum a - posteriori ( map ) , 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 . in fact , the main statistically motivated surface modeling methods are based on principle component analysis which is not sparse [ 2529 ] . 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 . 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 , are provided to improve the accuracy considering limited sample size . most of the previous surface modeling / representation works are designed for either static models [ 810 ] , or particular deformable organs with known physical properties ( such as motion cycle ) . 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 . 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 section 2 , 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 . the srds algorithm consists of three main steps to achieve sparse representations of deformable surfaces , as outlined in figure 1 . 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 . spherical harmonics are solutions to laplace 's equation expressed in the spherical coordinate system , defined as 
 ( 1)ylm(,)=(1)m2l+14(lm)!(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 ] . 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=0x(,)ylm(,)sin  d  d. 
   assume that harmonics up to level l ( l  l ) are involved in the transformation . let y denote the matrix composed of all ( l+1 ) discretized harmonics , so y has the following formation : 
 ( 4)y=|y0,0y1,1y1,0y1,1yl , 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,0  f1,-1  f1,0  f1,1    fl , l]t is the harmonic coefficient vector . , 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   perform shd for each of the k training deformed surfaces in x = { xk}k=1 , so the group of deformations can be described by matrix f 
 ( 7)f=|f1f2fk|(l+1)2k . consequently , the training set of deformations can be totally characterized by columns in f as 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   is adopted since it features better generalization and less computational cost compared to the gold standard k - svd algorithm . 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 pursuitinitially , each vector fk of length ( l+1 ) 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 =  , f = f.(2)subspace searching and clustering  

 i = i + 1 ; choose a vector fi from f ( e.g. first column of f ) and let fi = fi-1fi to remove fi . find ni vectors from f=ffi for representing fi with error no larger than  within emax iterations and the ni vectors form si . perform svd decomposition on si : uv = si ; let ai contain the first ni vectors of u ; update d = d  ai . select vectors from f that can be represented by ai with error no larger than  , and then remove them from f.repeat above steps until all the vectors are clustered . 


  initially , each vector fk of length ( l+1 ) 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 =  , f = f.(2)subspace searching and clustering  

 i = i + 1 ; choose a vector fi from f ( e.g. first column of f ) and let fi = fi-1fi to remove fi . find ni vectors from f=ffi for representing fi with error no larger than  within emax iterations and the ni vectors form si . perform svd decomposition on si : uv = si ; let ai contain the first ni vectors of u ; update d = d  ai . select vectors from f that can be represented by ai with error no larger than  , and then remove them from f.repeat above steps until all the vectors are clustered . initialization : i = 0 , d =  , f = f. subspace searching and clustering   

 i = i + 1 ; choose a vector fi from f ( e.g. first column of f ) and let fi = fi-1fi to remove fi . find ni vectors from f=ffi for representing fi with error no larger than  within emax iterations and the ni vectors form si . perform svd decomposition on si : uv = si ; let ai contain the first ni vectors of u ; update d = d  ai . select vectors from f that can be represented by ai with error no larger than  , and then remove them from f.repeat above steps until all the vectors are clustered . i = i + 1 ; choose a vector fi from f ( e.g. first column of f ) and let fi = fi-1fi to remove fi . find ni vectors from f=ffi for representing fi with error no larger than  within emax iterations and the ni vectors form si . perform svd decomposition on si : uv = si ; let ai contain the first ni vectors of u ; update d = d  ai . select vectors from f that can be represented by ai with error no larger than  , and then remove them from f. repeat above steps until all the vectors are clustered . b. subspace pruningone 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 . 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 . 
 c. matrix f factorizationafter identifying the inherent subspaces , the coefficient matrix f of training set can be partitioned into two - part union of subspaces and corresponding coefficients via the following procedures.since each vector fk has been clustered into the belonging subspace during the subspace identification process , the corresponding coefficients for each fk can be obtained accordingly . suppose that fk lives in subspace ai which is spanned by ni orthogonal basis , so its corresponding coefficients ck can be calculated using ls estimator 
 ( 9)ck=(aitai)1aitfk . 
								 then fk can be characterized by ck in its subspace 
									 ( 10)fk = aick.if there are totally j subspaces identified from f , a structured dictionary constructed by concatenating all deformation subspaces is established as d = i=1{ai } , with dimension i = i=1ni . since each vector fk lies in one of the subspaces , fk can also be represented in the structured dictionary with a block sparse vector { ck } , which is obtained via extending the coefficients { ck}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 c1 . if f1 lies in subspace a2 which are spanned by the 5th , 6 , and 7th columns in d , then c1 has nonzero values only at index of 5 , 6 , and 7 . consequently , f can be factorized as 
 ( 11)f = dc , 
								 where c=|c1  c2    ck|ik is the corresponding coefficient matrix with block sparsity . 
  after identifying the inherent subspaces , the coefficient matrix f of training set can be partitioned into two - part union of subspaces and corresponding coefficients via the following procedures . since each vector fk has been clustered into the belonging subspace during the subspace identification process , the corresponding coefficients for each fk can be obtained accordingly . suppose that fk lives in subspace ai which is spanned by ni orthogonal basis , so its corresponding coefficients ck can be calculated using ls estimator 
 ( 9)ck=(aitai)1aitfk . 
								 then fk can be characterized by ck in its subspace 
									 ( 10)fk = aick . if there are totally j subspaces identified from f , a structured dictionary constructed by concatenating all deformation subspaces is established as d = i=1{ai } , with dimension i = i=1ni . since each vector fk lies in one of the subspaces , fk can also be represented in the structured dictionary with a block sparse vector { ck } , which is obtained via extending the coefficients { ck}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 c1 . if f1 lies in subspace a2 which are spanned by the 5th , 6 , and 7th columns in d , then c1 has nonzero values only at index of 5 , 6 , and 7 . consequently , f can be factorized as 
 ( 11)f = dc , 
								 where c=|c1  c2    ck|ik is the corresponding coefficient matrix with block sparsity . a. sparse representation of training surfacesintegrating 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=1{ai } , g is inherently structured by subspaces of g = i=1{gi } with gi = yai of size n  ni.up to this point , with matrix g that captures the deformation features in the considered population , each training deformation xk in x can be fully characterized by a sparse coefficient ck as 
 ( 13)xk = gck . 
								 the sparsity of ck has already been illustrated in figure 2 . 
  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=1{ai } , g is inherently structured by subspaces of g = i=1{gi } with gi = yai of size n  ni . up to this point , with matrix g that captures the deformation features in the considered population , each training deformation xk in x can be fully characterized by a sparse coefficient ck as 
 ( 13)xk = gck . b. sparse representation of testing surfacesfor 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 as h = { hm}m=1 , 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=1 and choose the subspace with minimal projection error . since the number of subspaces j and the dimension of each subspace { ni}i=1 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 that gi * is the chosen subspace , so coefficient vector cm can be estimated with ls as 
 ( 14)cm=(gitgi)1githm . then block sparse vector 
									 ( 15)cm=[0  0    cm    0  0  0 ] 
								 is obtained according to the rules described in section 2.2 . further , the sparsity of ck or cm 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 )   or matching pursuit ( mp )   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 . 
  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 . the testing set is denoted as h = { hm}m=1 , 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=1 and choose the subspace with minimal projection error . since the number of subspaces j and the dimension of each subspace { ni}i=1 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 that gi * is the chosen subspace , so coefficient vector cm can be estimated with ls as 
 ( 14)cm=(gitgi)1githm . then block sparse vector 
									 ( 15)cm=[0  0    cm    0  0  0 ] 
								 is obtained according to the rules described in section 2.2 . further , the sparsity of ck or cm 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 )   or matching pursuit ( mp )   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 . for an organ with both interior and exterior surfaces , such as bladder the above theory can be extended to achieve sparse surface representation for deformations occurred on both interior and exterior wall of the object . we denote xkin and xkex(1  k  k ) as the corresponding interior and exterior of each training surface xk . then each pair of xkin ( with n1 vertices ) and xkex ( with n2 vertices ) can be approximated by spherical harmonic basis as 
 ( 16)xk=[xkin^xkex^]=[yinooyex][fkinfkex ] 
							 where y of size n1  ( l+1 ) and y of size n2  ( l+1 ) denote the spherical harmonic basis for inner and outer surfaces , respectively . therefore , each deformation is represented by vector fk=[fkinfkex ] , and all k training frames can be characterized by f={fk}k=1k as 
 ( 17)x = yf=[yinooyex][f1infkinf1exfkex ] . the following procedures of subspace identification and sparse surface representation as described in sections 2.2 and 2.3 can be applied straightforwardly . after identifying j subspaces d = i=1{ai } from sh coefficient matrix f , each training deformation can be sparsely represented with block sparse coefficient vector ck as : 
 ( 18)xk=[xkin^xkex^]=[yinooyex]dck = gck = gick , 
							 where gi=[yinooyex]ai is the subspace with size of ( n1 + n2 )  ni , and g = i=1{gi } is the desired structured dictionary . accordingly , ck is the block sparse vector , and ck is the nonzero coefficient values in the selected subspace . similar to other surface modeling methods [ 25 , 27 , 28 ] , the proposed approach requires point - wise correspondence across difference surfaces besides rigid registration . specifically , this point - to - point alignment is established , such that an identical spherical parameterization can be applied in the shd procedure . 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 . different correspondence methods have been proposed , such as minimum description length , sh - coefficient alignment , and so forth . we applied the sh based method   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 . 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 , figures 4(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 . if the condition of single ray - surface intersection applies , deformations can be resampled through the following steps to achieve point correspondence . construct an icosahedron of w 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 . align the center of the 3d surfaces to the origin of icosahedron such that rays casting from the origin can intersect with the surface .   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 . figure 4(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 . 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 . 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 . the initial 3d models of different organs are fed into a fem - based surgical simulation tool to generate deformation data for testing . for instance , figure 6 demonstrates two examples of shape distortions due to the endoscope poking and grasping one side of the gallbladder . table 1 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 , 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^kxk||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 . 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 cm . 
 a. training resultsduring 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 . table 2 shows that the subspace number j and dimensions of resulting dictionary ( dim  ( g ) = i ) are markedly small relative to n or l 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 . 
  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 . table 2 shows that the subspace number j and dimensions of resulting dictionary ( dim  ( g ) = i ) are markedly small relative to n or l 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 . b. sparsity and accuracy evaluationsparsity is examined in terms of ( / ) , where  is the average l0 norm of the coefficient vector ck ( training ) or cm ( 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 ) ) , 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 figure 7 . 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 . sparsity is examined in terms of ( / ) , where  is the average l0 norm of the coefficient vector ck ( training ) or cm ( 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 ) ) , 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 figure 7 . 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 . 
 c. effect of the performance of srds algorithm 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 . 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 . 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 . 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 . d. effect of the influence of coefficient truncation threshold  on the performance of srds algorithm is also tested while  is varied among [ 00.0001  0.0005  0.001  0.005  0.01  0.05  0.08  0.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 . the influence of coefficient truncation threshold  on the performance of srds algorithm is also tested while  is varied among [ 00.0001  0.0005  0.001  0.005  0.01  0.05  0.08  0.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 . 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 figure 15 . 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 . 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 surface x and x is defined as  
 ( 20)d(x , x)=max  px  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  px  ||p , p||2 
							 with ||||2 denoting the euclidean norm . 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 . table 4 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 . 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 . 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 . figure 17 visualizes the color - coded error distribution at all vertices on the represented surface with srds relative to the actual mri scans . 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 . 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 . .   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 . we can see that the sparsity and representation error resulting from srds method is very close to that using osp . figure 18 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 . 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 . figures 19 , 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 1921 indicates that errors in testing set is larger than that in training set and  ex1 the proposed approach is also tested over the in vivo cardiac mr images , 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 . 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 . table 8 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 . figure 22 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 , 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 . 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 particularly , in the 90th percentile evaluation in figure 23 , only 3% of all the point - wise errors are above 0.9  mm . to examine the efficiency of the proposed srds method quantitatively , the computational time to represent each surface in the testing set using srds method the results are summarized in table 10 , 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 . to 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 . 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 , srds algorithm has already been used in tracking organ deformations in minimum invasive surgery [ 24 ] . 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 .
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 .