Multimodal image registration of the scoliotic torso for surgical planning
© Harmouche et al.; licensee BioMed Central Ltd. 2013
Received: 26 May 2012
Accepted: 26 November 2012
Published: 4 January 2013
Skip to main content
© Harmouche et al.; licensee BioMed Central Ltd. 2013
Received: 26 May 2012
Accepted: 26 November 2012
Published: 4 January 2013
This paper presents a method that registers MRIs acquired in prone position, with surface topography (TP) and X-ray reconstructions acquired in standing position, in order to obtain a 3D representation of a human torso incorporating the external surface, bone structures, and soft tissues.
TP and X-ray data are registered using landmarks. Bone structures are used to register each MRI slice using an articulated model, and the soft tissue is confined to the volume delimited by the trunk and bone surfaces using a constrained thin-plate spline.
The method is tested on 3 pre-surgical patients with scoliosis and shows a significant improvement, qualitatively and using the Dice similarity coefficient, in fitting the MRI into the standing patient model when compared to rigid and articulated model registration. The determinant of the Jacobian of the registration deformation shows higher variations in the deformation in areas closer to the surface of the torso.
The novel, resulting 3D full torso model can provide a more complete representation of patient geometry to be incorporated in surgical simulators under development that aim at predicting the effect of scoliosis surgery on the external appearance of the patient’s torso.
Idiopathic scoliosis is a disease characterized by a complex three-dimensional curvature of the spine and the rib cage; these internal curvatures are externally manifest as a lateral trunk asymmetry and/or a rib hump. Such external deformations are often aesthetically undesirable for patients and can cause psychological problems, and in more severe cases, chronic back problems or pulmonary problems . Treatments aim at slowing down the curvature progression or at correcting some of the undesired curvature. They include a brace in less severe cases and surgery in the form of vertebral fusion in more severe cases. Surgeons rely on their experience and intuition in order to establish the adequate instrumentation that would lead to the desirable post-operative external trunk appearance. However, the effects of the brace or surgical instrumentation on the external shape of the trunk cannot be reliably predicted prior to completion of treatment. Our group has recently developed a simulator with the aim of predicting the effect of scoliosis surgery on the spine and the torso . Interesting preliminary results were obtained. However, the simulation outcome lacked generalizability across patients. Current research aims to integrate both bone and soft tissue in the 3D model of the trunk in order to verify whether surgical simulators that model soft tissue information could be useful in improving the prediction of the effects of surgery on the external appearance of the patient’s trunk. Such a model would require fusion of soft tissue information, typically obtained from MRIs in prone position, spine and rib cage data, typically obtained from X-rays in standing position, and a representation of the external surface of the trunk, obtained from an active vision system. This registration task is not trivial because of a mixture of both rigid and non rigid deformations that occurs between the acquisition of the different image modalities, in particular between the MRI and the remaining modalities, since the patient is lying down during MRI acquisition. Adding to the difficulty is the lack of correspondences between the tissue visible on the MRI and the structures visible on the remaining modalities. The aim of the present work is to register MRI data acquired in prone position with X-ray and TP data acquired in standing position, all while compensating for the postural changes that occur between the acquisitions and while preserving the rigidity of bone structures, in order to obtain a 3D representation (incorporating bone structures, trunk surface, and soft tissue information) of the torso of a patient with scoliosis.
Multimodal medical image registration has been applied to several types of images such as MRI/X-ray, MRI/CT, CT/X-ray, MRI/ultrasound, etc. Registration techniques can either be rigid, affine, semi-rigid, or deformable. Rigid or affine registration techniques apply either a rigid or an affine transformation to the entire source image being registered with a target image. Early work used rigid registration techniques in order to register MRI/CT and MRI/PET data of the brain . In terms of non-rigid methods, several techiniques have been used in medical image registration. Those consist of thin-plate spline , free-form deformations using B-splines [5–8], elastic models , fluid models [10–13], and Markov random field approaches [14–16].
Most of the work on MRI/X-ray registration was not focused on the spine and consisted mainly of rigid registration methods [17–25] to the exception of  which use perspective transformations and , which use free-form deformations. A review of these techniques can be found in  and .
Rigid registration techniques cannot capture the complex deformation that occurs in the shape of the spine between the standing and prone positions in which the different image modalities are acquired. Furthermore, with the knowledge that the vertebrae themselves are rigid structures, traditional purely non-rigid registration algorithms are also not appropriate for the task at hand. For example, Skerl et al.  perform registration of vertebrae using MR and CT images, but do not model the rigidity of these structures. As a compromise, semi-rigid techniques have been used for the registration of spine data. For example, vertebral structures extracted from MRI data have previously been modeled as rigid bodies for registration purposes . Soft tissue was registered using modified thin-plate splines that allow segmented vertebral structures to be constrained to rigid deformations. Similar work was done by Rohr , requiring only a few correspondence points instead of a full segmentation of the rigid structures. Huesman et al.  register CT and MR images of the torso using thin-plate spline transformations. Vertebral rigidity constraints were incorporated into the thin-plate spline approximation parameter. However, acquiring CT data for the entire torso is not possible in clinical settings due to radiation issues. More importantly, CT and MR images are both acquired in prone position, which implies that non-rigid deformations due to posture changes are not taken into account. Loeckx et al.  use B-splines in order to register PET/PET and CT angiography/CT angiography while incorporating a rigidity constraint in the cost function. Wang et al.  use triangular B-splines to deform 3 sagittal 2D MR images of vertebrae by placing knots at the boundaries of segmented vertebrae in order to insure their rigidity. B-spline functions have the advantage of allowing multiresolution registration. However, they require a regular grid of image data, which is not available in the image modalities that we wish to register. None of the works above register MRI and X-ray data, and, since all the images are acquired in prone position, nor do they model the spinal deformation between postures. Some works that have registered MRI to X-ray spine data have either required fiducials on cadaveric data  or required 9 X-ray images per patient for adequate X-ray/ MRI registration ; both strategies which are not possible in clinical settings. In later work, Van de Kraats et al.  generate CT-like data from MRI in order to perform rigid registration of vertebral bodies from X-ray images of cadavers. Markelj et al.  use 2 2D X-rays in order to rigidly register MRI and X-ray images by matching 3D gradients of 3D images to 3D gradients reconstructed from the 2D X-ray images. Tang et al.  apply rigid registration to fit a model of the vertebrae only onto Sagittal MRI slices in order to measure spinal curvature. Harmouche et al.  used an articulated model representation of the spine in order to register the spine extracted from MRI to that obtained from X-ray data of the same patient. The remaining tissues on the MRI were not registered. Finally, to our knowledge, no previous works address the simultaneous registration of MRI, X-ray and TP data.
This paper proposes a method that registers MRI to both TP and X-ray data in order to construct a 3D model covering the thoracic and lumbar vertebral levels of the torso of a scoliotic patient. MRI slices are acquired in prone position and X-ray and TP data are acquired in standing position. External markers are placed on the patient’s skin prior to X-ray and TP acquisition in order to register the two modalities using thin-plate splines. The MRI is then registered to the TP and X-ray data using a two step process: first, bone structures guide an articulated model representing the 3D spinal deformations that occur between the two postures and serves as an initial transformation of the MRI data. This transformation is then refined by confining the soft tissue to the volume delineated by the trunk surface and the surface of the vertebrae using weighted thin-plate splines with constraints for bone rigidity and soft tissue elasticity. The availability of the surface contours allows us to model the rigidity constraint in a simpler fashion when compared to previous works. To our knowledge, this is the first work that combines bone, soft tissue, and surface topography information for a model of the human torso.
This article is separated as follows: The Methods section describes the data acquisition process and the articulated model based registration. Qualitative and quantitative results of the proposed method are compared with both rigid and simple articulated model registration in the Results and discussion section. The key conclusions and proposed future works are presented in the Conclusions section.
In order to obtain the 3D representation of a patient incorporating bone, soft tissue and surface information, MRI, X-ray images and TP data are obtained for scoliosis patients. The 3D spine models are then extracted from the MRI and X-ray images and the 3D positions of external markers are obtained from the TP and X-ray data. The goal is to then align MRIs of the torso with reconstructed thoracic and lumbar vertebrae from biplanar X-ray images and surface topography data (TP). Thus, we are searching for the independent transformations T TP−Xray and T MRI−Xray , which would transform MRI and TP data into X-ray space, respectively.
MRI, X-ray images and TP data are obtained for three adolescent patients with scoliosis at the Sainte-Justine Hospital in Montreal. Patient 1 is 153.4 cm tall, weighs 48.8 kg, and has a right thoracic curvature with a 50° Cobb angle. Patient 2 is 146.1 cm tall, weighs 41.7 kg, and has a left thoraco-lumbar curvature with a 57° Cobb angle. Patient 3 is 160.5 cm tall, weighs 46.2 kg, and has a right thoracic curvature with a 49° Cobb angle. The present research involving human subjects has been approved by the ethics committee of the CHU Sainte-Justine Research Center (Comité d’Éthique de la Recherche du CHU Sainte-Justine). Written consent is obtained from the participants and a parent or guardian. These 3 image modalities are obtained during the same day but at different times for each patient.
Postero-anterior and lateral radiographs are then obtained for each of the patients while they are in standing position, legs slightly appart, and arms to the front and bent upwards. An explicit calibration method  is then used in order to calculate the 3D position of vertebral landmarks manually identified by experts on both X-rays. The obtained landmarks allow the mapping of generic vertebral surface meshes onto the patient space; resulting in a patient-specific 3D geometric model of the spine from the X-ray data (Figure 2(a)). In this case, six of those landmarks are used per vertebra in order to generate the X-ray articulated model. These are placed on the centers of the superior and inferior plates of each of the 17 thoracic and lumbar the vertebral bodies, and below and above the left and right pedicles. The chosen landmarks are regularly used for vertebral reconstruction in our clinical setting since they are clearly visible on the X-rays. In addition, prior to image acquisition, radio-opaque markers which are clearly identifiable on the obtained X-rays are placed on the surface of the patient’s torso and affixed on top of the adhesive markers used for TP data acquisition. The 3D position of the markers is then calculated using the same methodology as for the vertebral landmarks. These external markers are used for TP/X-ray registration.
Sagittal MR images of the spine are routinely acquired in our clinical setting for scoliosis patients awaiting surgery (1.5 Tesla, TR/TE = 771/15, 704x704, 350 FOV, with a 0.5 mm by 0.5 mm in-plane resolution and 3 mm thickness, with a 3.6 mm separation between slices). In order to visualize the entire torso, T1-weighted axial MRI slices of 1 mm by 1 mm in-plane resolution, 2 mm thickness and 12 mm spacing between slices and covering the entire torso are acquired (Siemens Symphony system 1.5 Tesla, TR/TE = 650/12, 704x704) (Figure 2(c)). Since acquisition of the sagittal slices already requires a total of 30 minutes inside the MRI, the larger spacing between the axial slices allows for the entire trunk to be imaged in a reasonable amount of time for patient comfort. The time inside the MRI thus totals around 50 minutes including rest time between acquisitions. Since the sagittal and axial slices are acquired during the same scan, their location with respect to each other is known and they are thus pre-registered. The sagittal slices are thus used in order to situate the axial MRI acquisitions to the appropriate vertebral level and consequently to calculate the transformations required. This is required because the distance between the axial slices is large and the contours of the vertebrae are not all visible on these acquisitions. It must be noted that in cases where sagittal slices are not required for clinical purposes, we would be able to use the 50 minute total scan time in order to acquire the axial slices only, thus allowing for a significantly better resolution and better visibility of the vertebrae. In order to generate a 3D geometric model of the spine from MRI data, the 3D shape of the seventeen thoracic and lumbar vertebrae is manually segmented from the sagittal slices using TomoVision’s SliceOmatic software. Landmarks are then manually placed on the left and right edges of the posterior, anterior, inferior and superior ends of the vertebral body for all thoracic and lumbar vertebrae also using TomoVision’s SliceOmatic software. These landmarks are used in order to calculate the articulated model for MRI/X-ray registration.
where P is the vector of source points obtained from the TP data, V is the vector of target points obtained from the X-ray data, and K is the matrix containing U(r).
Since the spacing between consecutive axial MRI slices is large, they are treated as independent 2D slices and the continuity between successive slices is not adressed. Thus, each slice is registered separately to the remaining model using T MRI−Xray , which a composition of an articulated-model transformation (T MRI−bone ) and a thin-plate spline transformation (T MRI−soft ). The goal is to transform each MRI voxel into the space of the 3D X-ray model, taking into account the non-rigid deformation following the posture change subject to the following constraints: First, the spine extracted from the MR images has to be aligned with the X-ray spine model. Second, MRI data of the torso has to be contained within the TP volume, such that the contour of the torso on the MRI corresponds to the surface topography.
where D vertebra is the distance between point p and the border of its closest vertebra obtained from the MRI data and D surface is the distance between point p and the border of the segmented MRI torso. Since the meshes representing the X-ray vertebrae were originally obtained from cadaveric data and registered to patient space using a few interest points, the use of the MRI vertebral meshes were favoured in order to calculate this distance. The MRI vertebrae were segmented from real patient data. The weighting applied to the thin-plate spline at each voxel is such that the further we are from the vertebra, the more weight is given to the non-rigid deformation. The proposed transformation was chosen as a simple and sufficient way to model the deformation between the vertebrae and the torso surface, as there is a lack of correspondences between the different image modalities in that area. In addition, the present goal is to simply fill the volume contained within the surface topography with approximate soft tissue information that would allow for the creation of a 3D patient model for surgical simulations.
Following the desciption of the validation methods, qualitative and quantitative results for our registration method will be presented.
Higher Dice values signify better overlap between the 2 surfaces, with perfect overlap corresponding to a Dice value of 1. The assessment of Dice values is somewhat dependent on the size of the areas being compared. That is, larger areas have higher area to perimeter ratio resulting in inherently higher Dice values. Since in the present case the different methods are comparing the same area being deformed, the Dice value is adequate for the comparison. In addition, the Dice values calculated assume that the MRI slices are in the proper orientation following registration and correspond anatomically to the surface topography slice at the same location and orientation. We make this assumption as the MRI slices have been aligned using the articulated model. We tolerate any bias obtained following this assumption given the fact that our main aim is to verify whether the soft tissues fit within the space encompassed within the surface topography. The Dice measure is able to verify this aim. The mean and variance of the Dice values is calculated for all patients, in addition to testing for statistically significant differences between the Dice values of the different methods.
In order to provide a more precise Dice measure that does not depend on the orientation of the slices, a 3D Dice analysis is preformed. This analysis reduces the bias present in the 2D Dice calculations. In order to obtain a volume from the registered MRI slices, we first build a surface mesh over the entire torso by triangulating points sampled from the MRI contours using a Delaunay triangulation. Since the MRI slices have a large distance between them, a more detailed mesh is then obtained using linear subdivision. We then create a volume with a resolution of 1 mm x 1 mm x 3 mm and having slices parallel to the x-y plane, and set all voxels inside the surface to be equal to 1, otherwise 0. This resolution is chosen since it is only slightly larger than the MRI slice thickness, and such a voxel size has been successfully used for MRI segmentation applications in the past. A similar volume with the same resolution is obtained from the surface topography mesh. We then compare every voxel of the obtained MRI and TP volumes using the Dice measure.
In order to assess the quality of the deformation field resulting from our proposed registration, and to examine the areas that display the highest amount of deformation, the determinant of the Jacobian of the deformation field is calculated both with and without the weighting applied on the thin-plate spline transformation. The Jacobian of the deformation field J is defined as the matrix of partial derivatives of the deformation field. It has been used in previous works on image registration in order to ensure a smooth and topology preserving transformation . For each voxel, the transformation vectors are calculated using the methods described previously. The discrete partial derivatives making up the Jacobian are then approximated using the first order derivative of a Gaussian (sigma = 1). The determinant of the Jacobian is obtained from the partial derivatives and indicates the voxel-wise relative volume change. For example, a determinant value of 3 indicates that the original volume has decreased by a factor of 3. In order to insure volume preservation of the tissues and organs during registration, a determinant value of 1 is required. The determinant of the Jacobian for the deformation of all voxels of MRI slices are shown and areas of highest deformation are addressed. The mean and variance values for all slices are then calculated.
Figures 5(b) and 5(c) show the Dice values for patients 2 and 3, respectively. A significant improvement in Dice values is also observed in these cases as well when our method is compared with rigid and articulated registrations. For the case of patient 2, the average Dice values are 0.844±0.028, 0.830±0.031 and 0.974±0.008 for rigid, articulated, and proposed registrations, respectively. A statistically significant improvement is also observed in the case of patient 2 when the proposed method is used compared to both the rigid and the articulated model registration (p<0.01 in both cases). However, an improvement was not seen when rigid and articulated model Dice values were compared.
In the case of patient 3, the average Dice values obtained are 0.878 ± 0.052, 0.871 ± 0.010, and 0.976 ± 0.010 for rigid, articulated, and proposed registrations, respectively. As in the case of patient 2, a statistically significant improvement is observed when the proposed method is used compared to both the rigid and the articulated model registrations (p < 0.01 in both cases), and no improvement is observed when rigid and articulated model registrations are compared.
The 3D Dice measure is obtained for the 3 patients using the results from the proposed method. In this case, values of 0.939, 0.938, and 0.940 were obtained for patients 1, 2, and 3, respectively. These values are only slightly lower (by an average of 3.66%) than those obtained using the 2D Dice measure. This is to be expected since, unlike the 2D Dice measure, the 3D measure compares voxels within the patient’s surface that were not used in the registration process. More importantly, these errors are inherent to the fact that an interpolation, more so a linear interpolation, was used between consecutive MRI slices due to the low MRI resolution and not to the method itself. Nonetheless, the values obtained using the 3D measure are still considered excellent in the literature, and are still considerably higher (by an average of 11.52%) than the 2D values obtained for rigid and articulated registration. We can thus conclude that, although the spacing between consecutive MRI slices is high, our method still provides a volume of soft tissue information contained within the surface topography.
One obvious caveat with the Dice similarity coefficient to validate our work is that it does not measure the anatomical correctness with which our soft tissues were deformed. However, for the desired precision, given that the MRI resolution is low in the z-direction thus leading to a less precise overall model, that a measurable difference in trunk shape changes due to treatment is in the order of a tenth of a cm, and that no anatomical correspondences are present in the space contained between the vertebrae and the patient’s trunk surface, the current results are deemed satisfactory. Our main concern is to obtain a model with which the vertebrae maintain their rigid characteristics and the space contained within the surface of the patient is filled with soft tissue information. However, we must note that the presented registration framework does not take the physical characteristics of the tissues into account when modeling the deformations.
The results above show that the proposed method is able to register the MRI, X-ray and TP data of a human torso with satisfactory precision, and doing so while still compensating for the deformations that occur between images due to differences in posture in which these images are acquired. The residual registration errors still present in the results may be due to several factors. For example, the precision of the manual intervention required for the localization landmarks on all images has inherent limits. X-ray and MRI landmark localisation errors have been studied in previous works to be 2.1±1.5 mm and 3.17±3.3 mm. The MRI landmark localisation errors have been shown to significantly decrease to 1.57±1.13 mm when the centroids are compared, as is the case for our proposed registration method. In the case of the surface topography, where the resolution of the equipment is 1.1 mm, the landmark localisation error is assumed to be equivalent to the radius of the adhesive markers, thus 2.5 mm. The registration error between the surface topography and the X-ray data has been previously established at 2.7 mm. Ongoing work within our group is aimed at automating the landmark extraction process. Furthermore, unlike existing methods, our proposed articulated model method does not require actual correspondence points for registration. This is due to the fact that the center and the orientation of the vertebral bodies are used in order to calculate the transformation between vertebrae. Thus, the need for landmark extraction can be eliminated and replaced with the use of higher order primitives. This has the potential to reduce registration variability and to improve precision. In terms of the correspondence points used to drive the thin-plate spline registration between the MRI and surface topography, the accuracy is not studied. These correspondences might suffer from a lack of accuracy due to the fact that deformations in the z axis between the 2 modalities that are due to gravity have not been taken into account by the articulated model. However, a lack of anatomically significant correspondences makes it difficult to measure the accuracy. Since the landmark selection is automatic, we assume that the precision is mostly dependent on the resolution of the images. In this case, the image with the lower resolution is the MRI. The MRI resolution is 1 mm by 1 mm in-plane and has 2 mm thickness leading to a maximum distance of 2.45 mm between voxels.
Further refinement of the method can be made by incorporating tissue-specific elasticity constraints, which would require tissue-specific biomechanical analysis resulting in a more complex model. This would provide for more realistic anatomical deformations. It should also be noted that tissues contained within the boundaries of the ribcage are believed to deform differently from tissues outside of this boundary. Modeling the transformation of these tissues differently has the potential of improving registration precision. In addition, a thorough biomechanical analysis would allow us to model the effect of gravity on the non-rigid deformations resulting from a posture change. The implication is that gravity would have a different effect on the various anatomical structures being registered. However, the 3D voxel deformations resulting from gravitational forces are very likely to be outside of the plane of the acquired MRI slices. Thus, incorporating these deformations would require a higher resolution in the z plane, which is currently infeasible in a clinical setting due to prohibitive acquisition times.
A method to register MRI, X-ray and surface topography data was proposed. This method first registered surface topography and X-ray data using thin-plate splines, and then fit the MRI data onto the model by taking into account the non-rigid deformations that are due to the postural difference between acquisitions. An articulated model was used in order to approximate the vertebral deformations in the MRI, and the remainder of the soft tissues was deformed using thin-plate splines with vertebral rigidity and the surface topography as constraints. Visual results as well as 2D and 3D Dice values measuring the fit between the surface topography and MRI data were obtained for real data of 3 pre-operative patients with scoliosis in order to validate the proposed method. Both qualitative and quantitative results showed a significant improvement in fitting the MRI data with the X-ray and surface topography data when compared to both rigid and simple articulated model registration. The determinant of the Jacobian of the deformation field was obtained with and without the rigidity constraint, showing higher variations in the deformation in the anterior part of the slices, and showing lower variations closer to the vertebrae in the case of the proposed method. Future work aims at incorporating tissue-specific elasticity constraints to the registration process, and at using automatically extracted higher order primitives for the articulated model registration. The precision of the obtained registration results allows us to build a complete 3D model of a patient’s trunk including soft tissue, vertebral, and trunk surface information which can be incorporated in a surgical simulator under development in order to potentially better predict the outcome of scoliosis treatments.
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.