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 . 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.

Parametric Active Contour Model

A snake model is based on a parameterized contour

*v*(

*s*) defined by (

*x*,

*y*) coordinates in

^{2} such that

The contour defined in Eq.

^{1} is said to possess an energy,

*E*_{snake}, defined as the sum of the internal energy,

*E*_{int}, and the external energy,

*E*_{ext}, according to

where the internal energy,

*E*_{int}, is comprised of an elastic and a bending term

^{8}. The first internal energy term,

*E*_{el}, 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

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,

*E*_{bend}, is dependent on the second derivative of the parametric curve

*v*(

*s*) and a constant

*β*(

*s*) such that

where

*v*″(

*s*) is the second order derivative of the initial parametric curve.

The external energy,

*E*_{ext}, is the other term critical to defining the overall snake progression.

*E*_{ext} serves to pull the contour toward the image edges and is defined to be

where

*E*_{img} is a scalar potential function defined on the image plane

^{32}.

In the work described herein,

*E*_{img} is assumed to be dependent on the gradient of the intensity,

, of the initial gray scale image. Many different combinations for

*E*_{img} have been investigated in previous literature (see

^{10}^{,}^{25}^{,}^{26}), but the work described herein primarily applies an

*E*_{img} defined over the whole image that is dependent on the (

*x*,

*y*) location and has been adapted from

^{7} such that

where

*I*(

*x*,

*y*) is the initial gray-scale image,

*G*_{σ} is a Gaussian smoothing filter characterized with size

*h*_{size} and standard deviation

*σ*. The

symbol represents the gradient operator.

The quantity *E*_{img} is approximately one when there is no edge present and approximately zero when there is an edge present. The local minima of *E*_{img} thus attracts the snake. This attraction means that the selection of *E*_{img} 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 contour^{7}.

The final energy functional curve will result when *E*_{snake} is (locally) minimized. The global energy minimization is based on how the initial contour is specified, how edges are defined, how *E*_{int} and *E*_{ext} are defined, and the process by which the minimum of *E*_{snake} is found. The important considerations in the minimization of *E*_{snake} 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 interaction^{33}.

In order to minimize

*E*_{snake}, the Euler–Lagrange equation must be satisfied

^{34}. This requires

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

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 region

^{7}. 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 difficulties

^{32}.

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 function

. 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

^{2} 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

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,

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

Here,

*v*(

*s*) is derived from the initial contour

*ϕ*(

*x*,

*y*,

*τ*) according to the condition

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 be

^{35}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

Likewise, taking the directional derivative of

*ϕ* in Eq.

^{11} yields

^{36}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

*ϕ* is orthonormal to

*v*(

*s*) when

*z* is constant and the normal vector is given by

^{36}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

and then substituted into Eq.

^{13} to yield

Combination of the information provided by Eqs.

^{9} and

^{17} gives

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 curvature

^{37}, where curvature is defined as the divergence of the normal from Eq.

^{15}^{9}. This means that the function of the curvature of the propagating set,

*F*(

*κ*) from Eq.

^{18} can be rewritten as

Finally, substitution of Eq.

^{19} into Eq.

^{18} yields

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

that corresponds to the greatest gradient change

^{26}. The flow function that leads to the largest gradient change is

^{9}^{,}^{35}In the mapping

, there exists near

*T*=

0 an inverse mapping function

*ϕ* defined by

. The curvature function

*F*(

*κ*) that holds for Eq.

^{19} is

^{36}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

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 Ψ (see

^{7}^{,}^{25}^{,}^{26}), but the function is applied as in

^{10} such that

where

*E*_{img} 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 added

^{26}. Then Eq.

^{12} becomes

Equation

^{25} acts on level sets of the image, treating each level as an individual surface under evolution of its own constraints. The term

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 evolution^{7}^{–}^{9}.

A recent study has integrated the reinitialization into the level set formulation so that the signed distance function is maintained

^{38}. 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

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

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 *h*_{size}, 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 *h*_{size} 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 , 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.

To find the

*z*-axis shown in Figure , a preliminary

*z*-axis,

, 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

The

*z*-axis can now be found such that

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 where

^{39}In Eq.

^{30}, the variables

*S*_{x}*S*_{y} and

*S*_{z} represent the position of the point in the

*xyz* frame that is to be transformed to the origin,

*O*. The cos

_{Ij} 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

to the global frame vector

according to

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 previously^{39}. 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 vertebra^{39}^{,}^{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 ^{40}. The cosines in Eq. ^{30} refer to the cosines of the values provided in Table .

| **Table 1.**Kinematics limits to lumbar spine range of motion^{40} |