the process of aligning images so that the corresponding features can be related is called image registration . image registration methods have been discussed and 
 classified in books [ 14 ] and surveys [ 510 ] . a way to reduce the ad hoc nature of registration methods is to 
 require them to satisfy certain properties . researchers have realized the 
 importance of symmetry and transitivity of registration methods [ 1120 ] . in , ashburner et al . proposed an approximately symmetric 
 image registration method that uses symmetric priors . in , christensen and johnson 
 proposed a registration algorithm that approximately satisfies the symmetry 
 property . ( christensen and johnson used the term  inverse consistency  for 
 what we refer to as  symmetry . ) their idea is to estimate the forward and 
 reverse transformation simultaneously by minimizing an objective function 
 composed of terms that measure the similarity between the two images and the 
 consistency of the forward and reverse transformations . this approach requires 
 maintaining two transformations , computing their inverses and it has a tradeoff 
 among the terms in the objective function . in , rogelj and kovai proposed a 
 registration method that uses symmetrically designed forces that deform the two 
 images . the method is approximately symmetric , it requires maintaining forward 
 and backward transformation , but it does not use their inverses . in , krinjar and tagare proposed 
 an exactly symmetric registration method that is based on a symmetrically 
 designed objective function , but it requires the computation of the inverse 
 transformation . in , lorenzen et al . proposed an exactly symmetric 
 registration method that is based on a symmetrically designed fluid model . the 
 method uses two transformations whose compositions define the forward and 
 backward transformations in such a way that they are inverses of each other . beg and khan in 
 and avants et al . in   used an exactly symmetric registration method that 
 maintains two functions , which when composed appropriately give   forward and backward transformations that are exact inverses of each other . in , cachier and rey analyzed 
 the reasons behind the asymmetry in registration , proposed symmetrized 
 similarity and smoothing energies , but their implementation of the method was 
 not exactly symmetric . in , tagare et al . proposed an exactly symmetric 
 registration method that does not require to maintain both forward and reverse 
 transformations and compute their inverses . instead , the objective function , 
 which can be based on intensity differences ( e.g. , mean squared difference , 
 normalized cross - correlation ) or distributions ( e.g. , mutual information , 
 normalized mutual information ) is modified such that the method is symmetric 
 and only the forward transformation is needed . consequently , the objective 
 function has only one term , which avoids the problem of tradeoff among multiple 
 terms . their implementation is symmetric up to the machine precision . in , christensen and johnson 
 realized the importance of transitivity of image registration but did not provide a way to satisfy it . the above methods are either approximately or exactly 
 symmetric but none of them is transitive . in this paper , we propose a method to 
 modify any image registration algorithm such that it is provably symmetric and 
 transitive on an image sequence . symmetry and transitivity are especially 
 important in motion tracking applications ; they insure that a physical point is 
 tracked in the same way regardless of the order in which the images are 
 registered . if there are topological changes present in the image sequence , the 
 two properties can hold only over corresponding regions . motion within the body is present at 
 the system level , organ level , tissue level , cellular level , subcellular level , 
 and molecular level . in addition to the normal motion , pathology - induced motion 
 or changes can occur ( e.g. , osteoporosis , multiple sclerosis , and tumor 
 growth ) . in both normal and pathology - induced motion or changes , it is often 
 useful to compute the motion , that is , to register image sequences . such 
 information can improve our understanding of the normal function and diseases 
 as well as help develop better treatments . the presented approach is 
 illustrated on sequences of cardiac mr images , which if accurately registered 
 can provide clinically useful myocardial displacement and strain information . however , the same or similar approach can be used for the registration of any 
 other image sequence . let  denote a set of 
 real numbers and  an n - dimensional 
 metric space . an n - dimensional 
 intensity image is a function i :   . intensity images will be referred to as just images . 
 without loss of generality the set of all images is denoted as . an n - dimensional 
 geometric transformation is a function t :   . geometric transformations will be referred to as 
 just transformations . the set of all transformations is denoted as . let tid denote the 
 identity transformation , that is , tid(r ) = r ,     r  . let  denote the 
 composition of transformations . an image   registration operator is a function  :   . ideally , any image registration operator should have 
 the following three properties  i , j , k  . 

 identity . 
 an image registration operator , the result of the registration should not depend on the order of images , that 
 is , when an image registration operator is applied to two images , the obtained 
 transformation should be the inverse of the transformation obtained when the 
 order of images is reversed . the generated transformation from the second to the third 
 image composed with the generated transformation from the first to the second 
 image should be equal to the generated transformation from the first to the 
 third image . an image registration operator , when applied to two identical images should 
 generate the identity transformation . the result of the registration should not depend on the order of images , that 
 is , when an image registration operator is applied to two images , the obtained 
 transformation should be the inverse of the transformation obtained when the 
 order of images is reversed . the generated transformation from the second to the third 
 image composed with the generated transformation from the first to the second 
 image should be equal to the generated transformation from the first to the 
 third image . the proposed 
 approach is simple ; select a reference image and then perform the registration 
 of any two images from the sequence of images through the reference . the 
 reference can be an image from the sequence of images or an image similar to 
 the images in the sequence . let the reference be denoted as r and let  represent an 
 image registration operator that does not necessarily have any of the 
 properties from section 2.2 . a new image registration operator  is defined 
 as(4)(i , j)=[(j , r)]1(i , r),where i and j are any two 
 images from the sequence . it is assumed that the transformations generated by  from images in 
 the sequence to the reference image are invertable , that is , that [ (j , r ) ] always exists . proofidentity holds since(5)(i , i)=[(i , r)]1(i , r )  by  ( 4)=tid . symmetry holds 
 since(6)(j , i)(i , j)=[(i , r)]1(j , r)[(j , r)]1  (i , r )  by  ( 4)=[(i , r)]1(i , r)=tid  by  ( 4 ) . transitivity holds 
 since(7)(j , k)(i , j)=[(k , r)]1(j , r)[(j , r)]1  (i , r )  by  ( 4)=[(k , r)]1(i , r)=(i , k )  by  ( 4 ) . identity holds since(5)(i , i)=[(i , r)]1(i , r )  by  ( 4)=tid . symmetry holds 
 since(6)(j , i)(i , j)=[(i , r)]1(j , r)[(j , r)]1  (i , r )  by  ( 4)=[(i , r)]1(i , r)=tid  by  ( 4 ) . transitivity holds 
 since(7)(j , k)(i , j)=[(k , r)]1(j , r)[(j , r)]1  (i , r )  by  ( 4)=[(k , r)]1(i , r)=(i , k )  by  ( 4 ) . it should be noted that the only requirement for the 
 above result to hold is that [ (j , r ) ] exists . this 
 means that the obtained transformations and their inverses do not need to be 
 differentiable . however , a number of registration methods involve 
 regularization terms that use transformation derivatives , in which case the 
 registration operator needs to generate diffeomorphic transformations . if the jacobian of the transformation (j , r ) is positive 
 then the inverse transformation [ (j , r ) ] exists . if the 
 jacobian of the transformation 
					 (j , r ) is zero or 
 negative , the inverse transformation does not exist and the 
 reference - based registration operator given by ( 4 ) can not be used or it can be 
 used only over the part of the domain where the jacobian is positive . many 
 registration methods control the jacobian either directly [ 2226 ] or indirectly [ 12 , 2730 ] by using a smoothness term that penalizes extreme 
 warps to prevent singularities ( zero jacobian ) or folding of the space 
 ( negative jacobian ) , in which case the inverse transformation exists and the 
 reference - based registration can be used . while the 
 result of the previous section holds for any registration operator that 
 generates invertable transformations , here we illustrate the approach on two 
 sequences of cardiac magnetic resonance images using two nonrigid image 
 registration algorithms . we acquired a 
 3d anatomical cine mri scan together with a 3d tagged cine mri scan of a 
 healthy volunteer and then repeated the acquisitions four months later . the 
 volunteer was a 27-year - old male subject with no history of heart disease . the 
 purpose of the tagged scan was to validate the myocardial deformation recovered 
 from the anatomical scan . the scans were acquired using steady - state 
 free - precession short axis cine imaging ( flip angle   =   65 , tr   =   3.4 ms , te   =   1.7 ms ) covering the entire heart 
 on a 1.5 t clinical mri scanner ( intera , philips medical systems , best , the 
 netherlands ) . all the scans had contiguous short - axis slices with similar field 
 of view and phases covering the entire cardiac cycle and their parameters are 
 given in table 1 . the tags were applied immediately after the detection of the 
 r - wave from the ekg signal and the first frame was acquired at a delay of 15 
 milliseconds after the r - wave . two orthogonal sets of parallel planar tags with 
 about 9 mm plane separation were oriented orthogonal to the image planes . for both scans , for each acquired slice the scanner 
 recorded the rigid body transformation from the scanner coordinate system to 
 the slice . this allowed us to map all the slices to a common coordinate system , 
 that is , to spatially align the anatomical and tagged scans . similarly , the 
 scanner recorded the start time for each phase ( frame ) relative to the peak of 
 the r wave , which allowed us to temporally align the anatomical and tagged 
 scans . since the heart rate , that is , the duration of the cardiac cycle , was 
 not the same for the anatomical and tagged scans , we used relative time ( as a 
 percentage of the cardiac cycle ) for the temporal alignment . to recover the 
 myocardial deformation , we use thin plate splines ( tps )   to represent the 
 transformation between any two frames and then maximize the normalized mutual 
 information   to determine 
 the transformation parameters ( tps node positions ) . we use normalized mutual 
 information since it was shown to outperform several other images similarity 
 measures . since 
 myocardium is nearly incompressible and its volume does not change by more than 
 4% over the cardiac cycle , we constrain the optimization of the tps node 
 positions such that the jacobian of the transformation never deviates from 1 
 ( which corresponds to exact incompressibility ) by more than 4% . ( the purpose of this section is to illustrate the 
 approach of section 2.3 ( construction of symmetric and 
 transitive registration algorithms from nonsymmetric and nontransitive 
 registration algorithms ) , and instead of the method given in   we could have used any 
 other registration method . for this reason we did not present here all the 
 details of the used registration method and the interested reader this registration 
 algorithm we denote as 1 , while 2 represents its 
 unconstrained version . given that near incompressibility is a physical property 
 of the myocardium , 1 is expected to 
 be more accurate than 2 . each of the two operators was used to recover 
 myocardial deformation from the two cardiac mr image sequences in two ways : 
 sequential and reference - based . in the sequential approach , the deformation was 
 recovered from the first to the second frame , then from the second to the third 
 frame , and so on . in the reference - based approach , figure 2 shows 
 a short - axis section through the 3d image and 3d lv wall model surface for the 
 sequential and referenced - based recoveries using the two registration 
 operators .   to quantitatively evaluate the deformation recovery 
 accuracy we compared the cardiac deformation recovered from the anatomical cine 
 mri against the corresponding information from the tagged cine mri . we 
 developed a tool for interactive positioning of virtual tag planes in tagged mri 
 scans . the tag planes are modeled as splines that the user can position and 
 deform by moving control points . this allows the user to fit the virtual tag 
 planes to the tagged image as well as to compute tag plane intersections . once 
 the cardiac deformation is recovered from the anatomical cine mri using the 
 proposed method , it is applied to the virtual tag planes at ed and then 
 compared to the interactively positioned tag planes in other frames . in each 
 image slice , we compute the distances between the corresponding intersections 
 of orthogonal virtual tag lines ( in - slice cross - sections of the virtual tag 
 planes ) generated interactively and by the model . this allows for in - plane 
 ( short - axis ) deformation recovery validation . the out - of - plane ( long - axis ) deformation 
 is not evaluated with this procedure since the tag planes , being perpendicular 
 to the short - axis image slices , do not contain information about the 
 out - of - plane motion . virtual tag lines for the sequential and referenced - based 
 recoveries for the two operators are shown in figure 3 . table 2 contains the 
 distances between corresponding intersections of virtual and real tag lines at 
 end systole , which is the most deformed state relative to end diastole . let 1 and 2 represent the 
 symmetric and transitive modifications of 1 and 2 , respectively . operators 1 and 2 are defined by 
 ( 4 ) ( the end diastole frame is used for r ) , which 
 involves transformation inversion . since we use tps for transformation representation 
 and the inverse of a tps transformation can not be obtained analytically , we 
 invert the transformation numerically using the newton - raphson method for 
 solving nonlinear systems of equations . the numerical error of the computation of the 
 inverse transform is denoted as . if the jacobian of the transformation is positive 
 then the inverse transformation exists and  can be 
 specified to be an arbitrary small positive number , that is , the inverse 
 transformation can be computed with an arbitrary small error . while  can be set to 
 an arbitrary small positive value , in practical applications little is gained 
 if  is set to a 
 value smaller than two orders of magnitude smaller than the pixel size . the 
 reason for that is that the registration error is on the order of the pixel 
 size , and by setting  to one 
 hundredth the pixel size the error of the computation of the inverse 
 transformation is already negligible compared to the registration error , and 
 further reducing  does not 
 improve the registration accuracy . for a given image registration operator  , we define the following three errors . 

 identity error of image i is(8)eiden = maxr(i , i)(r)r. 
 symmetry 
 error of images i and j is(9)esym = maxr(j , i)((i , j)(r))r. 
 transitivity 
 error of images i , j , and k is(10)etran = maxr(j , k)((i , j)(r))(i , k)(r). 
 
 identity error of image i is(8)eiden = maxr(i , i)(r)r. 
 symmetry 
 error of images i and j is(9)esym = maxr(j , i)((i , j)(r))r. 
 transitivity 
 error of images i , j , and k is(10)etran = maxr(j , k)((i , j)(r))(i , k)(r). identity , symmetry and transitivity errors for 1 ,   2 ,   1 , and 2 are given in 
 tables 3 , 4 , and 5 , respectively . the errors for 1 and 2 in the three 
 tables were computed using    =   0.001 mm . the 
 dependence of eiden ,   esym , and etran on  is depicted in 
 figure 4 for 1 . theorem 1 says 
 that reference - based modification of any registration operator satisfies 
 identity , symmetry , and transitivity on an image sequence . while the theorem 
 always holds and the modified registration operator satisfies the three 
 properties exactly , section 3 demonstrates a practical application of the 
 reference - based registration using a relatively accurate registration algorithm 
 ( 1 ) and its less accurate version ( 2 ) . the method was applied to two cardiac cine mri 
 scans , which are periodic image sequences , but the same conclusions hold in the 
 case of linear image sequences . as expected , the constrained registration ( 1 ) outperformed its unconstrained version ( 2 ) , which can be seen in figures 2 and 3 and in table 
 2 for both sequential and reference - based approaches . the two figures and the 
 table also show that in this case the reference - based registration was more 
 accurate than the sequential registration . the advantage of the sequential 
 registration is that the difference between any two consecutive frames is 
 relatively small , while the reference - based registration deals with larger 
 deformations ( e.g. , from end diastole to end systole ) . however , the 
 disadvantage of the sequential registration is that there is accumulation of 
 error from frame to frame , which seems to overweigh the advantage of small 
 frame - to - frame difference . the reference - based constrained registration 
 ( column ( c ) in figures 2 and 3 ) was more accurate than the other three algorithms 
 ( table 2 ) . the difference in accuracy of the four algorithms is most 
 pronounced at end systole , and it can be seen as the different level of 
 agreement between the model contours and the underlying image in figure 2 and 
 between the virtual tag lines and the underlying image in figure 3 . 
			  since 1 and 2 use tid as the initial 
 transformation in searching for the transformation that maximizes the 
 normalized mutual information , both operators satisfy ( 1 ) exactly and 
 consequently have eiden = 0 for any image , 
 as it can be seen in table 3 . operators 1 and 2 involve 
 transformation inversion as defined in ( 4 ) , which is done numerically with an 
 accuracy of . this is why the identity errors for 1 and 2 in table 3 
 are approximately equal to . the symmetry errors ( table 4 ) for 1 and 2 are 
 approximately 1 mm , while for 1 and 2 they are 
 approximately 2. the reason for this is that the evaluation of ( 9 ) 
 involves a composition of the operators , each contributing approximately  to esym due to the 
 numerical inversion of transformation in ( 4 ) . similarly , the transitivity errors ( table 5 ) for 1 and 2 are 
 approximately 2 mm , while for 1 and 2 they are 
 approximately 2. the reason for this is that the evaluation of ( 10 ) 
 involves a composition of the operators , each contributing approximately  to etran due to the 
 numerical inversion of transformation in ( 4 ) . it should be noted that 1 and 2 have similar 
 values for eiden ,   esym , and etran . the reason for this is that the three errors depend 
 only on the accuracy of the computation of the inverse transformation , and not 
 on the registration accuracy ( 1 is more accurate than 2 ) . in fact , if 
 the inverse transformation could be computed exactly , the three errors would be 
 zero regardless of the registration operator . the three errors scale with  , as shown in figure 4 , and they can be made as small 
 as the machine precision . figure 4 also shows that eiden   ,   esym  2 , and etran  2. thus , the reference - based registration slightly 
 worsens the identity error and it significantly improves the symmetry and 
 transitivity errors over the sequential registration . for the two tested image sequences , the symmetric and 
 transitive registration methods 1 and 2 had smaller 
 registration errors than their nonsymmetric and nontransitive counterparts 1 and 2 , respectively ( table 2 ) . while all four registration 
 methods had either zero or nearly zero identity errors ( table 3 ) , 1 and 2 had very small 
 ( nearly zero ) symmetry and transitivity errors and 1 and 2 had these 
 errors on the order of the pixel size or even larger ( tables 4 and 5 ) . thus , 
 in this limited study , very small ( nearly zero ) symmetry and transitivity 
 errors were accompanied by reduced registration errors . however , to determine 
 if this is always the case one would need to prove it mathematically or at 
 least repeat the experiment on a large number of image sequences . to simplify the notation , it was assumed that all the 
 images had the same domain , but the same conclusions hold when the domains are 
 different . furthermore , the method involves the transformation inversion , which 
 is done only once ( after both images are registered to the reference ) , as 
 opposed to the methods proposed in [ 12 , 1416 ] that require computing the 
 inverse transformation in each iteration of the optimization . we used the normalized mutual information as the image 
 similarity measure for the registration algorithms 1 and 2 . we repeated all the experiments by using the mutual 
 information , mean square difference , and normalized cross - correlation as 
 alternative image similarity measures and the obtained results were nearly 
 identical to those reported in section 3 and for this reason they are not 
 included in the paper . these repeated experiments confirm that the conclusions 
 of the paper are not specific to the normalized mutual information . the 
 reference - based registration of image sequences is provably symmetric and 
 transitive . furthermore , a limited study showed that the reference - based 
 registration was more accurate than the sequential registration , although this 
 can not be guaranteed to always hold .