A Three-dimensional Deformable Brain Atlas for DBS Targeting. I. Methodology for Atlas Creation and Artifact Reduction

Atlas Digitization and Interpolation

The first step in the methodology was to scan photographs of each atlas slice for computer processing. Atlas slices in the sagittal plane were used. A high quality scanner (Epson^TM Expression 1600) was employed to produce TIFF image files. The scanner resolution was set to 600 dots per inch resulting in a resolution of 23.6 pixels/mm, scans of the photographs were later converted to high quality JPEG images. While the resolution of the scans may limit the accuracy of the depiction of the atlas structures, there are a number of other considerations that will practically limit its accuracy such as: overall print atlas 3D consistency (slice to slice consistency) [7], voxel resolution of the final reconstructed atlas (discussed in later sections), and size of the smallest feature to be modeled in the print atlas. While subcortical structures tend to be small, all were larger than around 2 mm along its most narrow axis. The resolution chosen for modeling of the structures within the digital atlas were therefore based off the size of the smallest relevant feature within these structures. The smallest features that we chose to model were about 1mm in size (about the width of the lamina within the globus pallidus interna) which would require a voxelized atlas with voxels of 0.5 mm or less (in each dimension).

Each section had a predefined coordinate system based on the anterior commissure and posterior commissure. For sagittal sections, each section had a set of predefined axes. The horizontal axis defined the commissural line, while the vertical axis defined the mid-commissural plane. When the sections with corresponding contours were physically scanned, the resulting images were not consistently straight or centered. To correlate these contours with each slice (through the coordinate system), it was necessary to reorient each scanned image. This reorientation required that the origin be defined as the intersection of the axes, and that the horizontal and vertical axes were rotated for proper alignment. A program was created that would read the image files and allow for selection of the origin, and selection of the furthest left and furthest right points on the horizontal axis. The horizontal axis was translated and rotated such that the horizontal axis was straight. The center of the image’s matrix was then defined as the origin, and the horizontal line was aligned along a single row on the image matrix. The required translation and rotation was calculated and a new centered and reoriented image was then saved to an updated image file.

Once each image was scanned and properly reoriented, the structures of interest were extracted. For example, the subthalamic nucleus (STN) would be extracted from each slice in the sagittal plane. Since each structure in each plane was defined by a contour, and there were many structures in each plane. The structures of interest had to be manually extracted from structures found in its surroundings. Each file was read into Adobe Photoshop^TM (San Jose, CA), and the structure of interest was identified and its contour selected and all other undesirable objects removed from the binary image. This method created an image that ultimately included only the contour of the structure of interest. The resulting image was saved as an updated file. Matlab™ (Nantick, MA) was used to contour each structure for each plane, and the resulting set of points was stored in a matrix that was indexed by the corresponding plane.

The set of points along one contour were connected to the set of points along the neighboring contour(s). In 2-D, this process would involve connecting points with a series of triangles in a method referred to as triangulation. A generalization of this method to 3-D is tetrahedization. Delaunay tetrahedization [8] was used to determine the set of lines that connect a point to its natural neighbors. Applying this algorithm across the set of contours of a structure results in a linearly interpolated surface. The surfaces of these structures typically are composed of a complex geometry and thus require a robust method for surface interpolation. Other groups have solved this problem with differing techniques and variable results [9-11], involving either linear or non-linear interpolation. For the interpolation of the surface and volume of the anatomic structure of interest we employed a linear interpolation method. Even though our atlas surfaces were based off of a linear interpolation method since we employ smoothing in a later step of our atlas creation process, the resulting atlas structures had a surface that was not as linear as the initially created surface. In general, higher order (non-linear) models are preferred to linear methods for data modeling. However in the special case of modeling print atlas data, there are known inconsistencies in the 3D spatial alignment of the atlas slices [7, 12, 13]. As such, the accurate creation of a surface to fit the print atlas slices may not be as important of an issue as the ability to reduce the effect of inconsistencies in the digital atlas which was accomplished in this work through smoothing.

The smoothing algorithm was implemented by convolving the data with a kernel. Two different kernels were utilized; the kernel used for all small to medium sized structures was a Gaussian kernel with a standard deviation of 0.5 mm. The kernel used for large sized structures (only striatum and thalamus) was a Gaussian kernel with a standard deviation of 0.25 mm. Since these structures are binary prior to smoothing, smoothing essentially operates upon the surface. As the larger structures have greater surface areas, less smoothing was necessary to achieve a similar effect. Upon smoothing, the mesh was no longer purely composed of binary intensities and a new surface was thus determined from an isosurface of an appropriate value. The smoothing created a surface that was no longer a simple linear interpolation of the initial data and was also able to reject some of artifacts in the data.

In Fig. (1A), we present an example of a triangulated surface which lacks smoothing and possesses edges that are not differentially continuous. Therefore, we subsequently smoothed the surfaces by first voxelizing the data and then convolving with a smoothing kernel. The smoothed surface of this figure lacks rough edges and may produce a more reasonable reconstruction (Fig. 1B). This process first requires a mesh of the structure of interest, a mesh with voxels of 0.25 mm³ were used. This voxel size roughly corresponds to 1/3 to 1/4 of the size of the smallest feature of any anatomical structure of interest. The individual voxels were created to possess binary intensity (prior to smoothing) and thus the intensity was set as 1 for a voxel that corresponded to the inside of a structure, and 0 for a voxel that corresponded to the outside region of a structure (Fig. 2A).

One of the issues with Delaunay tetrahedrization is the inherent inability to correctly interpolate concave objects. The determination of lines that connect the points of two convex contours to create a surface is a relatively simple process using tetrahedrization. However, if one of these surfaces is concave, there is ambiguity in how the points of the objects connect. The solution we employed was to consider each concavity within a convex structure as a series of convex objects. We then simply subtracted out the portion of the object necessary to recreate the original concave object. This process of creating structures from the union or elimination of convex sub-pieces was repeated across all subcortical structures that were modeled. The following atlas structures were created (from anterior most to posterior most): striatum, globus pallidus externa (GPe), anterior commissure (AC), globus pallidus interna (GPi), optic tract (OT), thalamus, zona incerta (ZI), subthalamic nucleus (STN), substantia nigra pars reticulata (SNr), and red nucleus. Five thalamic nuclei were also created using this approach: ventralis oralis anterior (Voa), ventralis oralis posterior (Vop), ventralis intermedius (Vim), ventralis caudalis (Vc), and centromedian (CM).

Intensity Cut-Off Determination

Due to the necessary voxelization of the data, interpolation, and smoothing used to create the digitized three-dimensional atlas, it was not ultimately a perfect representation of the print atlas. The voxelization of the structure can only approximate the data because the smallest feature that could be represented was limited by the size of our voxels (0.25 mm in each direction). In addition, the smoothing process made the binary digital voxelized atlas possess smoother boundaries than the original dataset. Two steps were undertaken to evaluate the digital atlas accuracy versus the print atlas. First, an intensity cut-off value was determined for each structure when the surface was smoothed, since each voxel could range in values between 1 and 0. The atlas was evaluated at various cut off values to determine how accurate the smoothed digitized atlas was compared to the respective planes in the print atlas.

The smoothing used in our methodology convolves the binary matrix of the structure with a Gaussian kernel. The cut-off value used to designate the boundary of a structure was determined using three methods: 1) an initial value was determined from analysis of a test cube, 2) a qualitative evaluation of the effect of the change of the cut-off value using custom in-house print atlas and digital atlas viewing software, and 3) a quantitative analysis of the effect the intensity cut-off value had on the resultant area of the structure. For the initial value determination, a 49x49x49 voxel cube was created and smoothed. Figure 2 shows an edge of a square taken from a middle cross-section of a cube and the effect of smoothing on the structure. For the qualitative evaluation, each atlas slice was displayed over the digitally reconstructed atlas and compared to the corresponding slices from the print atlas version. Finally, the GPi was used as an example for a medium sized structure to quantitatively determine the effect of the intensity cut-off value. Areas of intersection between the print atlas structures and digital atlas structures were calculated to quantitatively determine the difference between the reconstructed atlas and the print atlas. Three areas were calculated, the area of intersection, A ∩ B, the area in A that was not in the intersection of A and B, A ∩ B^C, and the area in B that was not in the intersection of A and B, B ∩ A^C, where A was a print atlas structure and B was a digital atlas structure in the same plane.

Atlas Evaluation

First, a test cube was smoothed and examined along edges and corners to determine reasonable values for the intensity threshold. In the case of an unsmoothed structure (as in Fig. 2A), the linearly interpolated value at the boundary of the cube is 0.5. The interpolated intensity value at the boundary of either of the smoothed cubes was 0.5. Upon smoothing with a Gaussian kernel using a standard deviation equal to one voxel (0.25 mm), the intensity at a boundary along the edge of the cube varied between 0.3 and 0.7 within 0.125 mm of the actual boundary. It varied between 0.18 and 0.82 within 0.25 mm of the actual boundary (see Fig. 2B). Similarly smoothing with a Gaussian kernel with a standard deviation of two voxels (0.5 mm), intensities of 0.4 and 0.6 were within 0.125 mm of the actual boundary, and intensities of 0.31 and 0.69 were within 0.25 mm of the actual boundary (see Fig. 2C). For these relatively small amounts of smoothing used, the intensities at the boundaries dropped off quickly to zero, so most intensity cut-off values between zero and one may produce reasonable descriptions of the boundaries of the respective structures.

Second, the sagittal reconstruction of the GPi was examined to quantitatively determine the effect of the intensity cut off value on an actual (medium sized) atlas structure. Since the GPi was of an ellipsoidal type shape, the results were not the same when compared to a cubed shape. In the case of the cube, the intensity cut-off value that exactly represented the boundary of the cube was found to be 0.5. As can be seen from Fig. (3), the optimal intensity cut-off value for any plane is not necessarily 0.5. Optimally, the dotted line and the solid line for all colors should be zero as that would show that the intersection of the print atlas slice and the digital atlas slice were exactly equal (A = B), and thus A ∩ B^C and B ∩ A^C would be the null set. In all cases, these intersection metrics did not reach zero area. This result, while seemingly unwanted, may actually be useful in reducing the effect of known artifacts in the print atlas such as spatial 3D inconsistency. Since the exact location and extent of these artifacts are not well characterized, smoothing may provide a reasonable methodology for reduction of these artifacts. Each plane exhibited differing effects on accuracy following smoothing. This effect may be useful in “smoothing out” inconsistencies from 3D spatial slice alignment.

The effect of the cut-off value was examined quantitatively, as seen in Fig. (3), as well as qualitatively (as in Fig. 4). Values that were less than 0.5 were chosen as the intensity cut-off value (these values were found to best model and approximate the structures). As noted previously, only particularly small or large values had a marked effect on the accuracy of the boundary. As can be seen in Fig. (4), the intensity cut-off values in the range of 0.3 and 0.5 did not drastically alter the shape and size of the structures, but did show a noticeable effect on accuracy of the boundary. The voxelization of the atlas, in addition, limited the size of the smallest distinguishable component of the structure. An example of this effect can be viewed on sagittal atlas planes 21.5 mm, 23.0 mm, and 24.5 mm, where it is evident that the protrusions of the striatum were not present due to their small size.

Atlas Application

Using in-house created software, we have displayed example images showing the digital atlas on sample MRI and microelectrode data. The panels of Fig. (5) are a set of sagittal images showing a T1-weighted MRI (parameters from [14]) with an intraoperatively acquired electrophysiological map (green dots) corresponding to the thalamus overlay. Fig. (5A) provides an example of target visualization without the aid of the atlas, panel B shows the same image of panel A but with an undeformed atlas overlaid on top. In Fig. (5C), a best fit linear affine deformed atlas is shown overlaid on the sagittal T1-weighted MRI. Atlas deformation (Fig. 5C) qualitatively appears to improve the fit of the atlas to the subcortical anatomy as well as the microelectrode map (specifically of the thalamus) of this patient.