Search tips
Search criteria 


Logo of digimagwww.springer.comThis JournalToc AlertsSubmit OnlineOpen Choice
J Digit Imaging. 2011 October; 24(5): 926–942.
Published online 2010 September 30. doi:  10.1007/s10278-010-9336-z
PMCID: PMC3180553

Development of Subject-Specific Geometric Spine Model through Use of Automated Active Contour Segmentation and Kinematic Constraint-Limited Registration


This paper describes the development of a patient-specific spine model through use of active contour segmentation and registration of intraoperative imaging of porcine vertebra augmented with kinematic constraints. The geometric active contours are fully automated and lead to a discrete representation of the image segmentation results. After determining errors within the segmentations, application of reliability theory allows the selection of active contour parameters to obtain best-fit segmentations from a stack of 2D images. The segmented images are then used in conjunction with C-arm fluoroscope images to simulate the result of intraoperative patient-specific model registration including patient and/or structure motion between preoperative and intraoperative scans. The results are validated through comparison of the error within the patient-specific model generated through use of the C-arm images with a model acquired directly from MRI images of the spine after motion. The results are applicable to the development of a wide variety of patient-specific geometric and biomechanical models.

Key words: Segmentation, template geometry, patient-specific model, spine, active contour, fluoroscopy, kinematic constraints, imaging, image error, registration


Segmentation and registration of magnetic resonance-derived images is of vital importance to the development of patient-specific computational models for visualization techniques. Generating an automatic method for segmentation and registration is expected to reduce surgical planning time and decrease the number of misdiagnoses. It also serves as an important first step toward providing surgical resident training simulators and has a potential application that includes integrating haptic or simulated force, feedback into surgical robots. It has the potential to vastly improve the current clinical visualization and feedback systems, which are presently provided by only a single C-arm fluoroscope image and cameras mounted on endoscopic tools1. Automated segmentation and registration of images from multiple imaging modalities can also lead to a direct measure of structure morphology and enable one to distinguish automatically between anatomical structures.

Differences from patient-to-patient exist and can affect surgical procedures. It is currently tedious and time-consuming for a manual observer to obtain patient-specific models from medical images. In addition, the resultant geometries are highly dependent on the expertise of the person selecting the edges and can vary significantly between two different observers26. Thus, a calibrated, predictable, and automated method of segmenting images is important because it can save time and provide a consistent segmentation result.

A standard process must be applied to generate patient-specific models. Active contours, or snakes, have been used extensively to develop curve evolutions via energy minimization techniques for segmentation and classification of images711. The active contour methods are processor intensive. Registration of an initial model with a partial model, individual segmented image, or key anatomical points provides for an updated model with much less computer processing time. In previous methods of automatic segmentation and registration, the mathematical derivations and corresponding segmented images were said to be correct if they resulted in contours that appeared to a human observer to be close to the image object but did not provide objective measures of the error (see7,1115). While some previous methods for registration provided error metrics, they were lacking either a consistent method to select the parameters or a numerical result for error calculation between the manually defined ‘truth’ contour and the automatically segmented contour1618. Though it was unclear how exactly the segmentation and registration parameters were selected, a former method measured the distance between the expected and registered locations of anatomical points to provide a very thorough understanding of the error19.

This paper seeks to provide an objective measure and document the performance of a novel registration algorithm that includes kinematic constraints. The registration process has been designed to expand on previous works2024 by combining the processes and applying them to develop a patient-specific model. The method handles anatomical damage and data abnormalities. This method will provide a reference to which future algorithms developed may be compared. Active contour methods as they relate to the fully automated segmentation and classification of the spinal porcine vertebrae imaged through standard magnetic resonance imaging (MRI) are compared and contrasted.

A segmentation and registration process was defined in an effort to provide clinicians with a three-dimensional model of the real-time surgical environment. Those model parameters to which the segmentation was most sensitive were identified and adjusted through iterative coordinate descent. Finally, kinematic constraint limitations were used to restrict the spinal motion and provide local reconstruction intraoperatively.


The parametric active contour segmentation model10,25,26 provides a basis to develop and discuss why the level-set active contour segmentation model7,2527 is appropriate for the segmentation applications applied herein. Validation is performed by applying an error metric, Hamming distance28, to identify the most sensitive parameters through Cotter statistical screening design29, and Latin hypercube sampling30. Next, iterative coordinate descent optimization31 selects the specific parameter values that are used to build the initial model and validation model. Registration is performed between the initial model and intraoperative C-arm fluoroscope images through use of anatomical landmarks and knowledge of kinematic constraints.

The images for the work described herein were all obtained through 3D Fast Spoiled Gradient-Recalled-Echo (3D FSPGR) MRI. The resultant images were 512 pixels with a scanning resolution of 1.0 mm in both the x- and y-directions. After data acquisition, the data was converted to .tif images and were treated as signal intensities throughout the parameterization, segmentation, and registration processes.

The complete process applied in the development of the patient-specific spine model is depicted by the flowchart shown in Figure 1. After the images were obtained and converted to .tif, a region of interest was determined and basic contrast stretching, or two-level threshholding, was performed. This was followed by the parameterization process, which consisted of four primary parts. A human observer first selected the spine from the images. The object selected by the human observer was assigned to be the correct shape for the following three parameterization steps. The next three parameterization steps, described in “Model Parameterization” section, consisted of Cotter statistical screening, Latin hypercube sampling, and iterative coordinate descent, respectively. The parameters obtained through the parameterization process were used as inputs to the segmentation process, developed in “Parametric Active Contour Model” and “Level-Set Active Contours” sections and accomplished through Eq. 27. At the termination of the segmentation process, the template geometry was constructed. The porcine vertebra samples were manipulated and secured into place so that they would not move between the intraoperative C-arm fluoroscope images and full MRI scans taken for validation purposes. A nurse selected anatomical points from the fluoroscope images that were observable on both the initial model and the fluoroscope images. These points were registered and used to independently reposition each vertebra with kinematic constraint limitations. The process, described in “Incorporation of Kinematic Constraints in Registration” section, led to the final registered model. Finally, the physical distance between each landmark location on the registered model was compared to the corresponding landmark location on the validation model.

Fig 1.
The steps required to obtain a patient-specific geometric spine model include parameterization, segmentation, template geometry construction, spine manipulation, fluoroscope imaging, and registration with kinematic constraints. For validation purposes, ...

The complete process was performed on six independent three-vertebra porcine samples to provide a basis for the reproducibility of the methods described herein. For simplicity, the data used to describe the development of the methods are from a single sample, but the average error across the entire six sample (18 vertebra) data set is also presented in “Performance Assessment” as part of the overall performance assessment.

Parametric Active Contour Model

A snake model is based on a parameterized contour v(s) defined by (x, y) coordinates in R2 such that

equation M1

The contour defined in Eq. 1 is said to possess an energy, Esnake, defined as the sum of the internal energy, Eint, and the external energy, Eext, according to

equation M2

where the internal energy, Eint, is comprised of an elastic and a bending term8. The first internal energy term, Eel, allows the parametric curve to be treated as a string with elastic potential energy. To control bending and to decrease the likelihood of stretching, it introduces a tension such that

equation M3

where the weight α(s) provides control of the elastic energy along different parts of the contour and the curve v′(s) is simply the first order derivative of the initial parametric curve.

The second internal energy term, Ebend, is dependent on the second derivative of the parametric curve v(s) and a constant β(s) such that

equation M4

where v″(s) is the second order derivative of the initial parametric curve.

The external energy, Eext, is the other term critical to defining the overall snake progression. Eext serves to pull the contour toward the image edges and is defined to be

equation M5

where Eimg is a scalar potential function defined on the image plane32.

In the work described herein, Eimg is assumed to be dependent on the gradient of the intensity, equation M6, of the initial gray scale image. Many different combinations for Eimg have been investigated in previous literature (see10,25,26), but the work described herein primarily applies an Eimg defined over the whole image that is dependent on the (x, y) location and has been adapted from7 such that

equation M7

where I(x, y) is the initial gray-scale image, Gσ is a Gaussian smoothing filter characterized with size hsize and standard deviation σ. The [nabla] symbol represents the gradient operator.

The quantity Eimg is approximately one when there is no edge present and approximately zero when there is an edge present. The local minima of Eimg thus attracts the snake. This attraction means that the selection of Eimg is crucial to defining the final curve that results from the active contour evolution. Note that an increase in the value of σ will lead to a more blurry resultant boundary. However, the selection of large σ values is often necessary in order to increase the capture range of the active contour7.

The final energy functional curve will result when Esnake is (locally) minimized. The global energy minimization is based on how the initial contour is specified, how edges are defined, how Eint and Eext are defined, and the process by which the minimum of Esnake is found. The important considerations in the minimization of Esnake are that the result is accurate (resulting in edges close to the true boundary), precise (repeatable), and efficient so that it requires low amounts of operator interaction33.

In order to minimize Esnake, the Euler–Lagrange equation must be satisfied34. This requires

equation M8

Note that the function v(s) in Eq. 7 is treated as v(s, τ) where τ represents the ‘level’ of the snake and will be described in “Level-Set Active Contours” section. This allows the equation to be treated dynamically. The dependence on τ will vanish once the solution to v(s) stabilizes. Subsequently, the partial derivative of v(s, τ) with respect to τ is set equal to the left-hand side of Eq. 7 such that

equation M9

where the solution of Eq. 8 evolves to the desired segmented image. Recall that v′(s) and v′′(s) were given earlier in this section to be the first and second derivatives, respectively, of v(s). Likewise, v′″(s)and v″″(s) are taken to be the third and fourth derivatives, respectively, of v(s).

The solution to Eq. 8 is nontrivial because the sampling rules and initial conditions affect the final segmentation result. Error can be introduced as a result of the discretization of the higher-order derivatives. In addition, concave regions can lead to poor convergence because the external energy described by Eq. 5 leads to opposing forces. The result is that the contour is not pulled down into the concave region7. Also, as the value of σ in Eq. 6 is increased, the effects of sharp boundaries will be maintained for points further away, but the boundary becomes blurry and less accurate. The level-set method described in “Level-Set Active Contours” section adds an additional dimension to the problem, which both aids in the numerical implementation and overcomes the aforementioned difficulties32.

Level-Set Active Contours

In the level-set model, the contour v(s) from Eq. 1 is defined as the level set of an implicit functionequation M10. The third dimension, τ, represents the ‘level’ of the function v(s).

This notation has been adopted because it is parameter free and intrinsic. The notation also covers topological discontinuities because different topologies of the zero-level set do not imply different topologies of ϕ since topology changes are connected via the fourth dimension. This means that the model covers the special case of joining and separating contour lines in R2 as the model iterates through different contours on its way to the final shape.

The question of minimizing Eq. 2 now becomes an issue of minimizing ϕ(x, y, τ) because v(s) is incorporated within the ϕ(x, y, τ) contour. The values of v(s) can be obtained at any time from the ϕ function by simply taking the intersection of ϕ and a given ‘level’ of τ.

First consider the general case of curve propagation such that

equation M11

where κ is the curvature of the zero set at a particular (x, y) location, λ is a scalar value that limits the curve length, F(κ) represents a function of the curvature of the propagating set, N is the normal vector to the propagating curve, and the level set v(s) represents a 2D surface mapped from a 3D implicit contour ϕ(x, y, τ) as

equation M12

Here, v(s) is derived from the initial contour ϕ(x, y, τ) according to the condition

equation M13

Note that ϕ in Eq. 11 is taken to be positive if the point (x, y) is outside the zero-level contour given by ϕ(τ = 0). Conversely, the negative sign is chosen if the point (x, y) is inside the initial zero-level contour given by ϕ(τ = 0)27.

Consider the derivative of the curve ϕ(x, y, τ) with respect to τ. The flow function that leads to the largest gradient change is then given to be35

equation M14

where, Ψ is a function that will be small near the edge and will force the evolution to stop.

Now the goal is to determine how ϕ(x, y, τ) should progress so that the evolving surface v(s) will be maintained as the zero-level set within ϕ(x, y, τ).

Differentiating ϕ(x, y, τ) with respect to τ yields

equation M15

Likewise, taking the directional derivative of ϕ in Eq. 11 yields36

equation M16

Note that taking the directional derivative of ϕ as was done in Eq. 14 gives the arc length of the final segmented curve and will result in no change because the gradient will be zero on the boundary.

The relationships given in Eqs. 13 and 14 can now be used in conjunction with the definition of an orthonormal curve to show that [nabla]ϕ is orthonormal to v(s) when z is constant and the normal vector is given by36

equation M17

In Eq. 15, the left-hand side is related to the surface v(s), while the right-hand term is based off the terms of the surface ϕ(x, y, z, τ).

Now the relationship given by Eq. 15 can be rearranged to yield

equation M18

and then substituted into Eq. 13 to yield

equation M19

Combination of the information provided by Eqs. 9 and 17 gives

equation M20

The surface properties that allow the transition from Eqs. 17 to 18 need to be clarified via differential geometry. Specifically, at a particular point on the surface, the properties such as orthogonality that remain invariant are given by the curvature37, where curvature is defined as the divergence of the normal from Eq. 159. This means that the function of the curvature of the propagating set, F(κ) from Eq. 18 can be rewritten as

equation M21

Finally, substitution of Eq. 19 into Eq. 18 yields

equation M22

where Eq. 20 indicates how the implicit embedding function ϕ should deform for a given embedded level set evolution v(s).

Solving the geodesic curve evolution problem is now equivalent to searching for the steady state solution equation M23 that corresponds to the greatest gradient change26. The flow function that leads to the largest gradient change is9,35

equation M24

In the mapping equation M25, there exists near T = 0 an inverse mapping function ϕ defined by equation M26. The curvature function F(κ) that holds for Eq. 19 is36

equation M27

In Eq. 22, F(κ) represents a function of the curvature of the propagating set. The function ϕ refers to the implicit function that contains the level set function v(s) as discussed throughout “Level-Set Active Contours” section. The terms ϕx and ϕy refer to the partial derivatives of ϕ with respect to the x or y dimension. Likewise, a term ϕij refers to a second derivative of ϕ with respect to the i and j dimensions.

Now using the product rule for divergence, Eq. 21 can be rewritten as

equation M28

where, Ψ is a function that will be small near the edge and will force the evolution to stop. Several possibilities exist for the selection of Ψ (see7,25,26), but the function is applied as in10 such that

equation M29

where Eimg was given in Eq. 6.

One would expect that the evolution given by Eq. 12 would lead to the desired edge contour. However, to actually achieve the goal of attracting the contour to the bottom of the potential well, a constant inflation term, v, needs to be added26. Then Eq. 12 becomes

equation M30

Equation 25 acts on level sets of the image, treating each level as an individual surface under evolution of its own constraints. The term equation M31 is the sum of the two principle curvatures of the level sets of ϕ. The variable Ψ represents a scalar field function defined over the range of the model that represents the edges within the images. The final term in Eq. 25 corresponds to a stopping force for a particular image feature. In the event of an edge, this term will be close to zero.

In Eq. 25, the level set function ϕ can develop sharp and/or flat shapes during the active contour evolution. This makes further computation extremely inaccurate. Thus, it is important to maintain the evolving level set function as an approximate signed distance function, especially in the neighborhood of the zero-level set. Many previous schemes have approached this issue by performing a relatively ad hoc reinitialization process to periodically redefine the level set function as a signed distance function during the evolution79.

A recent study has integrated the reinitialization into the level set formulation so that the signed distance function is maintained38. The right-hand side of Eq. 25 is maintained in the neighborhood of the zero level set by multiplying it by a univariate Dirac function δ(ϕ) such that

equation M32

An additional term was also added to penalize the deviation of ϕ from a signed distance function such that, after applying the product rule for divergence, Eq. 25 can be rewritten as

equation M33

where Eq. 27 is the final level set formulation without reinitialization that is implemented for the work described herein.

Note that there are six parameters included in Eq. 14 that will need to be adjusted for the active contour. The first two parameters, σ and hsize, are related to the Gaussian filter used to find the edge-indicator function (Ψ). The value of σ is the standard deviation of the filter. The value of hsize represents the number of rows and columns in the filter. The next parameter, τ, is the level-set step size. The fourth and fifth parameters, λ and v, provide a limit to the curve length and elasticity, respectively. The final parameter, ε, gives the width of the Dirac function as was indicated in Eq. 26.

Incorporation of Kinematic Constraints in Registration

The origin of the coordinate system used for kinematic constraint application is set at the vertebral body centroid. As shown in Figure 2, the x-axis is defined to originate at the origin and point in the anterior direction. The y-axis is directed superiorly from the origin at the centroid along the line connecting the superior and inferior endplates of the vertebral body. The z-axis extends laterally from the centroid.

Fig 2.
The vertebral coordinate system has been defined to ensure an orthogonal basis for landmarks47. Points 1–4 on this image identify the anatomical points used to validate the transformed registration model, as described in “Anatomical Landmark ...

To find the z-axis shown in Figure 2, a preliminary z-axis, equation M34, is defined to be parallel to the line connecting the anatomical landmarks of the left and right pedicles. Then the x-axis, which points anteriorly, is calculated according to

equation M35

The z-axis can now be found such that

equation M36

where use of Eqs. 28 and 29 to find the z-axis ensures that the three coordinate axes are perpendicular.

It is important to note how anatomical movements are shown in this coordinate system. The x-axis represents lateral bending and anterior–posterior translation. Lateral bending to the right corresponds to positive rotation about the x-axis, while lateral bending to the left corresponds to negative rotation about the x-axis. The y-axis represents distraction–compression translation and axial rotation. Clockwise axial twisting corresponds to a negative rotation about the y-axis, while counterclockwise axial twisting corresponds to positive rotation about the y-axis. Medial–lateral translation and flexion/extension rotation occur along and around the z-axis. Flexion, or bending the body forward, is negative rotation about the z-axis. Similarly, extension, or bending the body backward, is positive rotation about the z-axis.

The transformation of the position and orientation of a collection of points in the local xyz frame can be expressed in terms of the fixed (global) XYZ frame through a rotation about the origin, O followed by a translation. The transformation matrix, [T] gives this conversion through a modified affine transformation where39

equation M37

In Eq. 30, the variables SxSy and Sz represent the position of the point in the xyz frame that is to be transformed to the origin, O. The cosIj indicates the cosine of an angle between the I-th axis in the global XYZ frame and the j-th axis in the local xyz frame. The transformation matrix [T] is used to convert the local vector equation M38 to the global frame vector equation M39 according to

equation M40

For the specific application described herein, the transformation given in Eq. 31 must be defined on external anatomical landmarks that are identifiable on C-arm fluoroscope images. In order to register each of the vertebra within a three vertebra column, the relationship provided by Eq. 31 had to be applied three times. The requirement that the relationship was applied once for each vertebra was necessary because the individual segments were registered as successive independent units with a defined relationship to adjacent vertebral segments. A human observer identified three landmarks per vertebra common to both the initial model and the C-arm fluoroscope images based on landmarks identified previously39. The spinal motion must also not exceed kinematics constraints. For this work, the transformation between initial spinal position (preoperative) and final spinal (intraoperative) position was first determined. Next, the system confirmed that kinematics constraints were satisfied. In the event that the kinematics constraints were exceeded, the motion was redefined so that it followed the previously defined maximum motion between lumbar vertebra39,40. The limits that were applied for the work described herein were obtained from a previous study. For completeness, the vertebra-to-vertebra limits of the human lumbar region obtained from the previous study are provided in Table 140. The cosines in Eq. 30 refer to the cosines of the values provided in Table 1.

Table 1.
Kinematics limits to lumbar spine range of motion40


Hamming Distance

The Hamming distance28 was originally defined as the number of bits that were different between two bit vectors. Extending this concept to the two-dimensional case indicates

equation M41

where H is the Hamming distance28. The symbol [plus sign in circle] represents a logical exclusive or function indicating that a given point within the union of the contours is deemed inside either the manually defined contour or the automatically segmented contour but not both. For the application of Hamming distance to the parameterization process, the image A(i, j) is a logical mask taking on values of 1 inside the automatically defined contour and 0 outside the automatically segmented contour. Likewise, the image M(i, j) is a logical mask taking on values of 1 inside the manually defined contour and 0 outside the manually defined ‘truth’ contour. In the work presented herein, the truth contour was determined by having a human observer select all the points along the bone edges. However, for the comparison between the transformed registration model after registration and the validation model obtained directly, A(i, j, k) was the transformed registration model after registration while the ’truth’ image M(i, j, k) was taken to be the validation model, where k was the image layer and the summation was performed over i, j, and k.

Anatomical Landmark Location Measurement

In order to validate that the registration worked properly, a validation model was obtained by directly segmenting the MRI data of the simulated intraoperative spinal positioning. The comparison of key anatomical landmarks on the final transformed registration model (intraoperative spinal position) and the validation model provided knowledge of the transformed registration model accuracy. The registration model was generated by first obtaining an MRI of the natural spinal positioning and then transforming it using the points that were selected on the C-arm fluoroscope images. The four anatomical landmarks identified in Figure 2 were found on both the transformed registration model and validation model. This was done while keeping in mind that the MRI scanning resolution was 1.0 mm. For each point, the distance between the respective landmark on the transformed registration model and the validated model were calculated.

Model Parameterization

Cotter Statistical Screening Design

The Cotter statistical screening design identifies the parameters with the greatest impact on error by estimating the combined effects of a parameter p for equation M42 as29

equation M43


equation M44


equation M45

In Eqs. 33, 34, and 35, H is the Hamming distance discussed in “Hamming Distance” section and H0 is the baseline case when all parameters are set at their minimum allowable values. In total, 2k + 2 trials had to be run in order to obtain all the Hamming distance values required to determine the Cotter statistical sensitivity levels for the k total parameters under investigation.

An observation of the correlation, rp1p2, between any two pairs of parameters for a design k  4 indicates that41

equation M46

Notice that the correlation approaches 1 as k  ∞. This means that the Cotter statistical design can identify the significant parameters, but it cannot specify the extent to which each parameter is independently significant because interactions among various parameters can generate an arbitrarily large significance value for a given parameter. Thus, Cotter statistical design was primarily used to identify the most significant parameters. These parameters were then evaluated using the Latin hypercube sampling method and measured via the sensitivity indicator discussed in “Latin Hypercube Sampling and Sensitivity Indicator” section.

Latin Hypercube Sampling and Sensitivity Indicator

The statistical method of Latin hypercube sampling was developed as a variance reduction technique designed to generate a distribution of the possible collections of parameter values from a multidimensional distribution30. It was viewed as a screening technique in which sample values were highly controlled but still allowed to vary. The Latin hypercube sampling technique has been widely used in engineering for deterministic simulation and risk analysis because it leads to an unbiased sample set with a variance less than that obtained through simple sampling30,42,43.

The Latin hypercube sampling defines a Latin square grid in an arbitrary number of dimensions where each sample is the lone sample in each axis-aligned hyperplane containing the sample30. A square grid containing sample positions is a Latin square if and only if a single isolated sample exists in each row and each column44. The Latin hypercube sampling generalizes this rule to an arbitrary number of dimensions where each sample is the lone sample in each axis-aligned hyperplane containing the sample30. The vector of input variables can be represented as

equation M47

and a given output variable can be modeled as

equation M48

where H(V) is a deterministic but unknown function (i.e., Hamming distance) that is a function of the input parameters vp for equation M49 are the parameters that define the active contour (i.e., bending, elasticity, and damping). To sample N vectors of input variables, the range of each input variable can be divided into N intervals with one observation on each input variable selected at random in each interval. This means that there are N observations on each of the k input variables. One of the observations on v1 is randomly selected and paired with a randomly selected observation on v2. This process is repeated for all possible vi through vk and results in the first vector of input variables. Next, one of the remaining observations of v1 is matched randomly with one of the remaining observations on v2 and so on until the second vector of input variables is obtained. This procedure is repeated until all of the input variables and observations are used and N Latin hypercube samples are defined.

Once the test bins are defined based on Latin hypercube sampling, a standard metric must be calculated to express the sensitivity of a given parameter. In this case, a standard economic principle applied to each bin measures the normalized activity level with respect to changes in the parameter45,46. The sensitivity indicator S for the Pth parameter from bin b is then given as

equation M50

where Hc is the Hamming distance for the test case, H0 is the Hamming distance for the baseline case for bin b, vc is the parameter value for the test case, v0 is the parameter value for the baseline case of bin b, and R is the range of possible values for the parameter. In the work described herein, the value for vc in a given bin was v0 plus 5% of the range for the parameter being tested. Once the sensitivity indicator was determined for each variable in each of the N bins, the average sensitivity across all bins was calculated for each of the variables according to

equation M51

Note that this method requires N(k + 1) total trials because the baseline error needs to be established for a particular bin when no parameters have been adjusted.

Iterative Coordinate Descent Optimization

Iterative coordinate descent, or cyclic coordinate descent, is a parameterization process used in practice to minimize a cost function31. In this work, Hamming distance was the cost function minimized through the approach. Iterative coordinate descent cycles through the k parameters using each in turn as a search direction. For example, at the first iteration, fixing all variables except for the variable with the largest sensitivity indicator as given in Eq. 40 allows the selection of that variable that minimizes the Hamming distance given by Eq. 32. On the next iteration, the process was repeated by fixing all variables except for the second most sensitive parameter and once again minimized the Hamming distance. This process was continued until all parameters were adjusted.


Parameter Sensitivity Assessment

The Cotter statistical design that was described in “Cotter Statistical Screening Design” section requires all variable parameters to be defined in order to determine the significance of modifying any given parameter29. The method requires each parameter to be set to its respective maximum and minimum values for the sensitivity analysis. The parameters, their meaning, initial settings, and minimum and maximum values are shown in Table 2. The results of running the Cotter statistical design are presented in Table 3.

Table 2.
Level set active contour parameters to select
Table 3.
Cotter statistical design results

The normalized sensitivity, S(p), was greatest for the v, τ, and ε parameters, as indicated by Table 3. Thus, v, τ, and ε were the most sensitive parameters and modifying one of these three parameters would lead to the largest effect on the overall segmentation error. However, a high correlation between any pair of variables means that the results from the Cotter design could not indicate which of the three variables would generate the most impact independently.

The Latin hypercube method previously described in “Latin Hypercube Sampling and Sensitivity Indicator” section resolved the ambiguity between any two pairs of parameters. As seen in Table 4, the Latin hypercube method demonstrated that the Dirac function width, ε, had an independent significance of more than 1.5 times that of the curve elasticity coefficient, v, which was the next most important parameter. The curve elasticity parameter was just slightly more important than the Gaussian filter standard deviation, σ. The independent significance of τ was the next most independently significant parameter, though it was only four-fifths the significance of the leading three parameters. The independent significance of the ε, v, σ, and τ parameters were more than twice that of the remaining two parameters.

Table 4.
Latin hypercube results

From the data presented in Tables 3 and and4,4, the iterative coordinate descent parameterization order previously described in “Model Parameterization” section was identified. In all cases, the parameter selections were made by first adjusting ε followed by v, σ, and τ, respectively. In all cases, this adjustment order was followed to perform parameter selections. The initial settings for v, σ, and ε were always 1.5. The initial settings for τ and λ were always 5. The initial setting for hsize was always 15 for the MRI.

Five images from the set of spinal MRIs were parameterized. A sample image indicating the manually defined contour and the automatically segmented contour is shown in Figure 3. The minimum and maximum values as well as the mean and standard deviation of each best-fit variable are presented in Table 5. Note that all parameterizations were only taken to two decimal places. Note also that all cases resulted in a final value for τ of 5.00. This outcome, which provided for a standard deviation of zero, indicated that the τ parameter is not as significant as initially indicated and could potentially be removed from future parameterizations for complex morphological structures by simply using τ set at 5. The remaining images from the complete data set were segmented using the mean value for each of the six parameters. Note that the mean value of hsize was rounded to the nearest integer value of 12.

Fig 3.
The method described in “Level-Set Active Contours” section was used to automatically segment the porcine cervical spine from the initial baseline MRI. The image at left shows the manual contour selected by the human observer in green ...
Table 5.
Level set parameter value ranges from training data

Model Generation

The mean parameter values given in Table 5 were used to generate the initial three-dimensional rendering of the spinal model shown in Figure 4 and the validation model shown in Figure 6. Recall that the validation model was generated directly from MRI images acquired after the spinal motion. The C-arm fluoroscope points of interest and the final transformed intraoperative patient-specific model after simulated patient motion generated through use of the initial model and C-arm fluoroscope image registration is shown in Figure 5.

Fig 4.
Initial spine model. The method described in “Level-Set Active Contours” section was used to automatically segment the porcine cervical spine from the initial baseline MRI. This model takes approximately 8 h to generate. This image ...
Fig 5.
Transformed Intraoperative Spine Model with Kinematics Constraints. The method described in “Incorporation of Kinematic Constraints in Registration” section was used to register the initial model shown in Figure 3 to the C-arm ...
Fig 6.
Validation spine model. The method described in “Level-Set Active Contours” section was used to automatically segment the MRI images that were obtained directly from the spine after motion. This model was generated for purposes of validating ...

The model shown in Figure 5 was generated through application of the registration methods described in “Incorporation of Kinematic Constraints in Registration” section, which required the initial model shown in Figure 4. It takes less than 5 min to register the initial patient model with the C-arm images obtained after motion (intraoperative scans). The five-minute estimate includes both the time necessary to have a nurse select points of interest and the time required to generate the registered model.

The validation model shown in Figure 6 was generated through the same manner as the initial spine model shown in Figure 4. This model was generated directly from MRI data of the spine after motion. This model was generated only for purposes of validating the transformed registration model through error measurements.

Performance Assessment

The Hamming distance associated with each of the three vertebra shown were calculated to be 7.68%, 8.44%, and 9.66%, for the structures of the top, middle, and bottom vertebra, respectively. As stated in the discussion on Hamming distance from “Hamming Distance” section, these numbers were calculated by assigning A(i, j) as the transformed registration model after registration while the ’truth’ image M(i, j) was taken to be the validation model with the percent error based on the union of the two surfaces.

Of particular note, the top vertebra in the spinal structure shown in Figure 4 had a fractured inferior articular process. When the registration was performed without application of the kinematic constraints, the observed error within the registration for that vertebra rose threefold to 24.2%. The initial porcine cervical spine model shown in Figure 4 was generated by segmenting the baseline MRI through the methods described in “Level-Set Active Contours” section. This model takes approximately 8 h to generate, which was a primary motivator for performing the much faster registration method described herein.

Final locations of the anatomical landmarks identified in Figure 2 were measured. The physical distance (or Hausdorff distance19) between the locations of the anatomical landmarks on the transformed registration model and the validation model are given in Table 6. The average distance between the anatomical landmarks was 2.42 mm for the model shown in Figure 2. Four landmarks were identified and measured on each of three vertebra from six total spine samples processed as part of the study described herein. For the 72 total landmarks measured, the average Hausdorff distance was 2.45 mm.

Table 6.
Distance between the anatomical landmarks identified in Figure 1 on the transformed model and the validation model


Active contours have been used extensively to develop curve evolutions because MRI-derived geometries are of vital importance to the development of patient-specific computational models for visualization techniques. Generation of such models leads to a reduction in surgical planning time, direct measurement of joint morphology, better mechanical modeling, and estimation of blast injuries. It could also provide an important component in development of haptic-feedback for surgical robot-assisted interventions. A standard process was applied to select the parameters required in the generation of a consistent model for determining differences in morphology from patient to patient. Hamming distance-based error calculations, statistical analysis of parameter importance, iterative coordinate descent, and affine transformations were combined to develop a novel registration procedure augmented with kinematic constraints against which future registration algorithms may be compared.

The parameterization process selected appropriate parameters that were subsequently used to build a three-vertebra porcine geometric template model from MRI images. Afterward, registration of the model with automatically segmented C-arm fluoroscope images simulated fully automated segmentation and registration of preoperative MR scans with intraoperative images obtained from the fluoroscope. This step was accomplished through use of kinematic constraints that indicated restrictions to spinal motion. The kinematic constraints accounted for anatomical damage and data abnormalities by using intraoperative images to update the initial patient-specific models. The results were validated through comparison of the patient-specific model generated through the initial MRI scan in conjunction with the C-arm images and a model developed from directly acquired MRI images of the spine after motion.

The initial goal was to generate the same resultant model through the C-arm fluoroscope method (Fig. 5) as would be obtained directly through MRI segmentations (Fig. 6) while saving overall processing time during an operation. It takes approximately 8 h to produce the models directly from the MRI images. This time may be decreased significantly by processing the collection of 2D images via parallel processing techniques. In the final model described herein, it takes less than 5 min to register an initial patient scan with the C-arm images obtained after motion (intraoperative scans). The five-minute estimate includes both the time necessary to have a nurse select points of interest and the time required to generate the registered model. It should be noted that all points were selected by a single nurse for the results described herein. Variances between observing nurses should be investigated in a future study in order to identify the learning curve that will be present in training observers while maintaining consistent and reproducible results.

The particular sensitivity levels of each parameter may be different for different imaging modalities or different morphological structures. For instance, if one were to obtain the images with computed tomography (CT) or micro-CT rather than through an MRI, or with an MRI data acquisition sequence other than the 3D FSPGR utilized herein, differences in the noise levels within the images could affect the importance of the Dirac function width. Second, the anatomical shape of the object being segmented may also influence the importance of the stiffness parameters and level step size. Before using the software in clinical applications, it will be important to repeat this work with multiple spines and the spinal column regions of interest (i.e., cervical, thoracic, lumbar), including mechanical measurement of surgical fiducial placement, to validate the accuracy and reproducibility of the model generation.

The method defined and applied in the work described herein is significant for four reasons. First, the method of kinematics constraint-limited motion registration handles differences in patient position between preoperative and intraoperative scans. Second, intraoperative scans account for changes in patient condition between preoperative and intraoperative scans (i.e., repositioning of individual vertebra as a result of worsening of bulging disks). Third, this method requires only an initial MRI scan and a C-arm fluoroscope. The MRI scans are common for disease diagnosis and surgical planning. The C-arm fluoroscopes are typically readily available in an operating room and often used in spinal surgeries to ensure proper tool placement. The models could be used to enhance the speed with which surgeons can ultimately perform interventions through haptic interfaces without adding any additional cost to the surgery. Finally, the additional preparation time to generate a better visualization method for the surgeon is minimal. This is important to minimize the amount of time the patient must spend under anesthesia.


The work described herein enhances previous medical image segmentation and registration techniques by combining reliability theory, segmentation error calculations, application of kinematic constraints, and validation with objective error measurements through a model obtained directly from MRI segmentation. It provides intraoperative local registration through application of kinematic constraints and handles anatomical damage and data abnormalities. This work applies all of these techniques in a combined package with a graphical user interface that provides an environment to easily identify the set of anatomical landmarks that will be used for registration purposes.


The authors would like to acknowledge the National Science Foundation (NSF) for providing funding support to the team members via the Graduate Research Fellowship (GRF) Program.

Contributor Information

Catherine G. Strickland, ude.eudrup@gcretals.

Daniel E. Aguiar, ude.eudrup@raiugad.

Eric A. Nauman, Phone: +1-765-4948602, Fax: +1-765-4967537, ude.eudrup@namuane.

Thomas M. Talavage, ude.eudrup@tmt.


1. Kwartowitz DM, Herrell SD, Galloway RL. Toward image-guided robotic surgery: determining intrinsic accuracy of the da Vinci robot. Int J Comput Assist Radiol Surg. 2006;1(3):157–165. doi: 10.1007/s11548-006-0047-3. [Cross Ref]
2. Elmore JG, Wells CK, Lee CH, Howard DH, Feinstein AR. Variability in radiologists’ interpretations of mammograms. N Engl J Med. 1994;331:1493–1499. doi: 10.1056/NEJM199412013312206. [PubMed] [Cross Ref]
3. Kundel H. Reader error, object recognition, and visual search. In Proc 2004 SPIE Conference. Bellingham, WA: SPIE; 2004. pp. 1–9.
4. Krupinski EA. The future of image perception in radiology: synergy between human and computers. Acad Radiol. 2003;10:1–3. doi: 10.1016/S1076-6332(03)80781-X. [PubMed] [Cross Ref]
5. Manning DJ, Gale A, Krupinski EA. Perception research in medical imaging. Br J Radiol. 2005;78:683–685. doi: 10.1259/bjr/72087985. [PubMed] [Cross Ref]
6. Abbey CK, Barrett HH. Linear iterative reconstruction algorithms: study of observer performance. In Inf Process Med Imaging. SPIE Bellingham. WA, Dordrecht: Klower Academic; 1995. pp. 65–76.
7. Xu C, Prince JL. Snakes, shapes, and gradient vector flow. IEEE Trans Image Process. 1998;7(3):359–369. doi: 10.1109/83.661186. [PubMed] [Cross Ref]
8. Kass M, Witkin A, Terzopoulos D. Snakes: active contour models. Int J Comput Vis. 1987;1(4):321–331. doi: 10.1007/BF00133570. [Cross Ref]
9. Osher S, Paragios N. Geometric level set methods in imaging vision and graphics. 1. New York, NY, USA: Springer; 2003.
10. Malladi R, Sethian JA, Vemuri BC. Shape modeling with front propagation: a level set approach. IEEE Trans Pattern Anal Mach Intell. 1995;17(2):158–175. doi: 10.1109/34.368173. [Cross Ref]
11. Tsai A, Yezzi A, Wells W, Tempany C, Tucker D, Fan A, Grimson WE, Willsky A. Shape-based approach to the segmentation of medical imagery using level sets. IEEE Trans Med Imaging. 2003;22(2):137–154. doi: 10.1109/TMI.2002.808355. [PubMed] [Cross Ref]
12. Kirbas C, Quek F. A review of vessel extraction techniques and algorithms. ACM Comput Surv. 2004;36(2):81–121. doi: 10.1145/1031120.1031121. [Cross Ref]
13. Suri JS, Setarehdan SK, Singh S, editors. Advanced algorithmic approaches to medical image segmentation: state-of-the-art applications in cardiology, neurology, mammography, and pathology, ch 3. London, England, UK: Springer-Verlag London; 2002. pp. 148–161.
14. Suri JS. Advanced algorithmic approaches to medical image segmentation: state-of-the-art applications in cardiology, neurology, mammography, and pathology, ch 8, pp 416–421. London, UK: Springer; 2002.
15. Pham DL, Xu C, Prince JL. Current methods in medical image segmentation. Annu Rev Biomed Eng. 2000;2(1):315–337. doi: 10.1146/annurev.bioeng.2.1.315. [PubMed] [Cross Ref]
16. Hill DLG, Batchelor PG, Holden M, Hawkes DJ. Medical image registration. Phys Med Biol. 2001;46:R1–R45. doi: 10.1088/0031-9155/46/3/201. [PubMed] [Cross Ref]
17. Hemler PF, Napel S, Sumanaweera TS, Pichumani R, VanDenElsen PA, Martin D, Drace J, Adler JR. Registration error quantification of a surface-based multimodality image fusion system. Med Phys. 1995;22(7):1049–1056. doi: 10.1118/1.597591. [PubMed] [Cross Ref]
18. Maintz JBA, Viergever MA. A survey of medical image registration. Med Image Anal. 1998;2(1):1–36. doi: 10.1016/S1361-8415(01)80026-8. [PubMed] [Cross Ref]
19. Banik S, Rangayyan RM, Boag GS. Automatic segmentation of the ribs, the vertebral column, and the spinal canal in pediatric computed tomographic images. J Dig Img. 2010;23(3):301–322. doi: 10.1007/s10278-009-9176-x. [PMC free article] [PubMed] [Cross Ref]
20. Siewerdsen JH, Moseley DJ, Burch S, Bisland SK, Bogaards A, Wilson BC, Jaffray DA. Volume ct with a flat-panel detector on a mobile, isocentric C-arm: pre-clinical investigation in guidance of minimally invasive surgery. Med Phys. 2005;32(1):241–254. doi: 10.1118/1.1836331. [PubMed] [Cross Ref]
21. Fahrig R, Butts K, Wen Z, Saunders R, Kee ST, Sze DY, Daniel BL, Laerum F, Pelc NJ. Truly hybrid interventional MR/X-Ray system: investigation of in vivo applications. Acad Radiol. 2001;8:1200–1207. doi: 10.1016/S1076-6332(03)80702-X. [PubMed] [Cross Ref]
22. Froelich JJ, El-Sheik M, Wagner HJ, Schenbach S, Scherf C, Klose KJ. Feasibility of C-arm-supported ct fluoroscopy in percutaneous abscess drainage procedures. Cardiovasc Intervent Radiol. 2000;23:423–430. doi: 10.1007/s002700010099. [PubMed] [Cross Ref]
23. Mitton D, Zhao K, Bertrand S, Zhao C, Laporte S, Yang C, An KN, Skalli W. 3D reconstruction of the ribs from lateral and frontal X-ray s in comparison to 3D CT-scan reconstruction. J Biomech. 2008;41:706–710. doi: 10.1016/j.jbiomech.2007.09.034. [PubMed] [Cross Ref]
24. Templeton A, Cody D, Liebschner M. Updating a 3-D vertebral body finite element model using 2-D images. Med Eng Phys. 2004;26(4):329–333. doi: 10.1016/j.medengphy.2004.01.004. [PubMed] [Cross Ref]
25. Geman S, McClure DE: Bayesian image analysis: an application to single photon emission tomography. In Proc Amer Stat Assoc Statistical Computing Section. 1985, pp 12–18
26. Caselles V, Kimmel R, Sapiro G. Geodesic active contours. Int J Comput Vis. 1997;22(1):61–79. doi: 10.1023/A:1007979827043. [Cross Ref]
27. Kichenassamy S, Kumar A, Olver P, Tannenbaum A, Yezzi A. Gradient flows and geometric active contour models. In Proceedings of the Fifth International Conference on Computer Vision. Washington, DC, USA: IEEE Computer Society; 1995. pp. 810–815.
28. Hamming RW. Error detecting and error correcting codes. Bell Syst Tech J. 1950;26(2):147–160.
29. Cotter SA. A screening design for factorial experiments with interactions. Biometrika. 1979;66(2):317–320. doi: 10.1093/biomet/66.2.317. [Cross Ref]
30. McKay MD, Beckman RJ, Conover WJ. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics. 1979;21(2):239–245. doi: 10.2307/1268522. [Cross Ref]
31. Nocedal J, Wright SJ. Numerical optimization. Springer series in operations research. New York, NY, USA: Springer; 1999. pp. 53–55.
32. Yezzi A, Kichenassamy S, Kumar A, Olver P, Tannenbaum A. A geometric snake model for segmentation of medical imagery. IEEE Trans Med Imaging. 1997;16(2):199–209. doi: 10.1109/42.563665. [PubMed] [Cross Ref]
33. Udupa JK, Herman GT, editors. 3D imaging in medicine, ch 1. 2. Boca Raton, FL, USA: CRC Press LLC; 2002. pp. 1–67.
34. Sapiro G. Geometric partial differential equations and image analysis, ch 1. New York, NY, USA: Cambridge University Press; 2001. pp. 44–63.
35. Epstein CL, Gage M. Wave motion: theory, modeling, and computation, v 7 of Math Sci Res Inst Publ, ch 2. New York, NY, USA: Springer; 1987. pp. 15–59.
36. Osher SJ, Sethian JA. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations. J Comput Phys. 1988;79(1):12–49. doi: 10.1016/0021-9991(88)90002-2. [Cross Ref]
37. Boyer CB: A history of mathematics, 2nd edition. New York, NY, USA: Wiley, 1991, Revised by Uta C. Merzbach
38. Li C, Xu C, Gui C, Fox MD. Level set evolution without re-initialization: a new variational approach. In Proceedings of the 2005 IEEE Computer Society International Conference on Computer Vision and Pattern Recognition. Washington, DC, USA: IEEE Computer Society; 2005. pp. 430–436.
39. Wade JA. An investigation of ovine lumbar kinematics using the Purdue spine simulator. Master’s thesis. West Lafayette: School of Mechanical Engineering, Purdue University; 2005.
40. Wilke HJ, Kettler A, Claes LE. Are sheep spines a valid biomechanical model for human spines? Spine. 1997;22(20):2365–2374. doi: 10.1097/00007632-199710150-00009. [PubMed] [Cross Ref]
41. Ryan TP, editor. Modern experimental design. Wiley series in probability and statistics. Hoboken, NJ, USA: Wiley-Interscience; 2007.
42. Owen AB. A central limit theorem for latin hypercube sampling. J R Stat Soc. 1992;54:541–551.
43. Loh WL. On latin hypercube sampling. The Annals of Statistics. 1996;24:2058–2080. doi: 10.1214/aos/1069362310. [Cross Ref]
44. Iman RL, Helton JC. The repeatability of uncertainty and sensitivity analyses for complex probabilistic risk assessments. Risk Anal. 1991;11:591–606. doi: 10.1111/j.1539-6924.1991.tb00649.x. [Cross Ref]
45. Zurada JM, Malinowski A, Usui S: Perturbation method for deleting redundant inputs of perceptron networks, 1997
46. Pannell DJ. Sensitivity analysis of normative economic models: theoretical framework and practical strategies. Agric Econ. 1997;16(2):139–152. doi: 10.1016/S0169-5150(96)01217-0. [Cross Ref]
47. Galle B. Development of a mechanical spine simulator and determination of lumbar kinematics. Master’s thesis. West Lafayette: School of Mechanical Engineering, Purdue University; 2005.

Articles from Journal of Digital Imaging are provided here courtesy of Springer