PDE-based registration
The basic premise of this work was to evaluate whether the equivalent potential energy distributions modeled by a partial differential equation (PDE) over an undeformed ("source") surface and a deformed ("target") surface could be used to determine correspondence between the two surfaces. For our specific method, finite element models of two classic PDEs (Laplace's equation and the diffusion equation) were solved independently over the source and target surfaces using the Galerkin method of weighted residuals with Lagrange polynomial interpolation [22]. Equivalent potential energy isocontours were calculated and matched on a closest-point basis, and from this correspondence the final displacements at all mesh nodes were interpolated.
Laplace's equation is most commonly used to describe potential flow problems in thermal, fluid, and electrostatic systems and is given by
where Φ represents the potential and σ describes the spatially varying conductivity. The diffusion equation, which utilizes a time-varying component, is given by
where Φ represents the potential and α is the diffusion coefficient.
To solve Laplace's equation (1), Dirichlet (Type I) boundary conditions were set to allow "flow" from a high- to low-potential area (Figure 1, Step 2). Specifically, nodes in the nipple and chest wall area were given boundary potential values of 1 and 0, respectively, and the conductivity σ was set to unity. The diffusion equation (2) was solved by temporally stepping the FEM solution using a fully explicit forward Euler scheme. For each data set utilized in this study, optimal ranges were empirically determined for time step ([1e-7, 8e-7]) and final time ([0.005, 0.01]). A pure Neumann (Type II, no flux) boundary condition was prescribed at the chest wall, and the potential field was allowed to propagate from a source located at the nipple. The diffusion coefficient α was set to unity. The diffusing front was stopped once the potential field reached the chest wall.
The solutions Φsource and Φtarget obtained from the PDEs were used to establish correspondence between the source and target nodes. This involved two distinct processes: finding point correspondence between isocontours of Φsource and Φtarget and interpolating the displacements at these isocontour points to all nodes in the mesh. In the first step, isocontours were extracted from Φsource and Φtarget for a set of selected isovalues (Figure 2, Step 4). The correspondence between the source and target isocontour points was determined by aligning the contours by their centroids and using the symmetric closest point (SCP) algorithm (Figure 2, Step 5. See Figure 3 for detailed description of SCP). In the second step, the displacement vectors at the source isocontours points were interpolated to all source nodes.
The method can be summarized in the following steps (Figures 1, 2):
-
1.
Obtain the undeformed source mesh and deformed target mesh that define a breast surface before and after deformation.
-
2.
Assign boundary conditions at nipple and/or chest wall nodes
-
3.
Solve the PDE (diffusion or Laplace) over the source and target meshes using FEM.
-
4.
Extract isocontours on the source and target surfaces.
-
5.
Determine point correspondence between source and target isocontours using SCP (Figure 3).
-
6.
Interpolate displacements at source isocontours to all source nodes.
TPS registration
As noted earlier, there are numerous methods of spline-based interpolation. TPS interpolation was chosen in part because it is a standard, well-characterized method in the literature [20] that has been successfully used in many non-rigid registration applications. Because it does not require a regular grid, the effects of changing a control point are relatively localized. The overall registration is achieved by the warping of a hypothetical thin sheet of metal using a series of radial basis functions based on a number of fixed control points. The global deformation field was then interpolated back to the surface node coordinates of the finite element mesh. The displacement vector at point (x, y, z) is therefore described by the following linear system:
where
with constraints
where X
i
, Y
i
, Z
i
are the coordinates of the control points, N is the number of control points, and a, b, c, and F are scalar weighting factors.
Simulation experiments
To assess the accuracy of the PDE-based and TPS methods described above, breast surfaces deformed using a biomechanical model were registered and the simulated displacements were used to calculate the registration error. A CT image volume of a pendant human breast (256 × 256 × 130, voxel size 0.6 mm3) and an MR image volume of a pendant breast from a different patient (256 × 256 × 98, voxel size 1.0 mm3, 3D T1-weighted fat-nulling inversion pulse sequence) were obtained. Both image volumes were segmented using ANALYZE 6.0 (Mayo Clinic, Rochester, MN), and triangular source meshes consisting of 6,313 and 3,942 nodes, respectively, were constructed.
Target surfaces were created by deforming the source surfaces using a finite element model of the breast as a linear elastic, Hookean solid under localized compression. The nodes along the chest wall were made to be fixed (Type I) and a Gaussian-shaped stress distribution (Type II) applied to the lateral surface of the breast. To provide a challenge to our registration method, the shape and magnitude of the applied stress was varied. First, the CT source surface was deformed using a circular contact area with a maximum displacement of approximately 33 mm. In the second simulation, a rectangular contact area was used with a maximum displacement of 13 mm. Finally, in the third simulation, the MR source surface was deformed assuming a rectangular contact area and a maximum displacement of 6 mm.
In addition to running the Laplace and diffusion registrations, TPS registration was performed for further comparison. Because the accuracy of the TPS method varies based on the number and distribution of control points, two different sets of registrations were performed. In the first analysis, a uniform distribution of control points over the breast surface was selected using a k-means algorithm. The number of control points was varied, and for a particular number of control points desired, 20 different configurations were selected to account for the initial random seeding of the k-means clustering algorithm. In the second analysis, a high number of control points in the deformed region (the part of the surface in contact with the simulated inflation bladder) and a lower number over the rest of the surface was used. Similarly, the error was calculated for a varying number of control points and configurations.
To assess the accuracy of the registration methods, the target registration error (TRE) was calculated as the Euclidean distance between the coordinates determined by non-rigid registration and the true target points. Because of the controlled model-based deformation, every node on the surface had a known correspondence from source to target surfaces, allowing for individual TRE determinations along with evaluation of mean and max error for the surface as a whole.
Phantom experiments
To test the registration methods on real-world data, a semi-anthropomorphic breast phantom was fabricated from an 8% w/v solution of polyvinyl alcohol (Flinn Scientific, Batavia, IL) that was frozen at -37°C in a plastic mold for 16 hours and thawed to ambient room temperature. Thirty-four 1-mm stainless steel ball bearings were implanted directly under the surface of the resulting tissue-mimicking cryogel to act as fiducials. The phantom was then placed inside a custom-built acrylic chamber designed to deliver compression by means of an air bladder positioned against the surface of the phantom (Figure 4).
CT images (512 × 512 × 174, 0.54 × 0.54 × 1 mm voxel size) were acquired with the phantom subjected to three different states of mechanical deformation (undeformed, 50% of maximum bladder pressure, and full inflation of approximately 200 mm Hg). Triangular surface meshes were obtained by segmentation of the image volumes using the surface extraction tools in ANALYZE, and the coordinates of the fiducial centroids localized. These meshes contained 8,127, 6,777, and 8,260 nodes, respectively. The Laplace, diffusion, and TPS methods were then used to register the phantom surface meshes as described above. Accuracy was assessed by calculating the TRE at the embedded fiducials. For the TPS method, 33 of the fiducials were used as control points in the interpolation and the remaining fiducial was reserved for calculating the TRE. To assess the error over the entire surface, the TPS registration was repeated in a "leave one out" scheme, each time using a different fiducial to calculate the error, with the final TRE representing the average over all trials.