Monte Carlo Simulation of Errors for N-localizer Systems in Stereotactic Neurosurgery: Novel Proposals for Improvements

Introduction: Frame-based stereotaxis has been widely utilized for precise neurosurgical procedures throughout the world for nearly 40 years. The N-localizer is an integral component of most of the extant systems. Analysis of targeting errors related to the N-localizer has not been carried out in sufficient detail. We highlight these potential errors and develop methods to reduce them. Methods: N-localizer systems comprising three and four N-localizers of various geometries were analyzed using Monte Carlo (MC) simulations. The simulations included native and altered geometric dimensions (Width [W] x Height [H]). Errors were computed using the MC simulations that included the x- and y-axes of vertically oriented rods, that altered the W/H ratio, and that added a fourth N-localizer to a three N-localizer system. Results: The inclusion of an overdetermined system of equations and the geometries of the N-localizer systems had significant effects on target errors. Root Mean Square Errors (RMS-e) computed via millions of MC iterations for each study demonstrated that errors were reduced by (1) inclusion of the x- and y-coordinates of the vertically oriented rods, (2) a greater triangular area enclosed by the diagonal fiducials of the N-localizer system (stereotactic triangle), (3) a larger W/H ratio, and (4) an N-localizer system that comprised four N-localizers. Conclusion: Monte Carlo simulations of Root Mean Square error (RMS-e) is a useful technique to understand targeting while using N-localizer systems in stereotactic neurosurgery. The application of vertical rod positions enhances computational accuracy and can be performed on any N-localizer system. Keeping the target point within the stereotactic triangle enclosed by the diagonal rods can also reduce errors. Additional optimizations of N-localizer geometry may also reduce potential targeting errors. Further analysis is needed to confirm these findings which may have clinical importance.


Introduction
Frame-based stereotaxis has revolutionized precise neurosurgical procedures throughout the world. Stereotactic space localization is a key step in the use of frame-based systems that provides coordinates for direct targeting in the brain. For this localization step, the N-localizer has been the most commonly used system for the past 40 years [1][2][3]. Depending on the manufacturer, the N-localizer has several dimensions (i.e., width and height) and configurations. Therefore, a careful and robust analysis of these systems may suggest improvements to the accuracy of the localization step, which may directly impact the targeting during the neurosurgical procedure.
In this study, we utilized Monte Carlo (MC) simulations to highlight potential targeting errors during the localization step by computing targets within the stereotactic surgical volume. We incorporated a mathematical optimization of overdetermined equations that exploits the existing geometry of N-localizer systems. Moreover, we analyzed various existing and theoretical N-localizer geometries.

Mathematics of the stereotactic localization
The mathematics of three N-Localizers used together for stereotactic neurosurgery have been published and extensively used [1][2][3][4][5]. As a standard approach, the method computes three three-dimensional (3D) points along the diagonal rods, given three corresponding two-dimensional (2D) points of intersection of the diagonal rods with an axial computed tomography (CT) image. Then a transformation matrix can be constructed to transform from 2D to 3D and from 3D to 2D. This classic method operates under the predicate that any stereotactic image should include all necessary and sufficient information for registration to its corresponding stereotactic system in order to accurately compute the target. Typically, this method applies a determined system of equations that includes three equations per axis, which is the minimum number of required equations to compute three unknowns.
Rather than utilizing a determined set of equations to compute the transformation matrix, the utilization of an overdetermined system of equations may minimize errors. The vertical fiducials represent a source of data that positions that may improve the computation and are not utilized directly in the determined system. In order to develop an overdetermined system of equations, we first need to set up a general transformation equation. The first step is to construct column vectors that convert 2D coordinates (u, v, t) to 3D coordinates (x, y, z) (Equations 1-3). In this transformation for imaging, the coordinate is generally set to any non-zero value such as 1 [1]. Then these equations can be assimilated into a general transformation matrix that easily accommodates an overdetermined system of equations (Equation 4). This overdetermined system of equations takes the optimized form that is solved by using a pseudoinverse technique that minimizes errors, such as with Singular Value Decomposition (SVD) or QR Decomposition [6].
The transformation matrix can then be rearranged into a 3x3 stereotactic transformation matrix ( ) (Equation 5) that can then be used to transform 2D to 3D coordinates (Equation 6) or to transform 3D to 2D coordinates, where denotes the inverse of (Equation 7). The determinant of should be computed prior to computing the inverse and rarely can be zero at a perfectly orthogonal zero plane. It should also be noted that the inverse of can also be computed by interchanging the (u, v, t) coordinates with the (x, y, z) in Equation 4 to avoid performing a matrix inverse operation.

Monte Carlo simulations for error analysis
Monte Carlo (MC) simulations were performed to introduce random errors into the computation of the stereotactic transformation matrix ( ). This process utilized 35x1024x1024 repeated random errors up to 1 millimeter (mm) per study set at 0.5 mm/pixel, which were isotropic in horizontal and vertical pixel widths.
The utilization of MC simulations has been previously published for single localizer comparisons [4], but here we utilize MC simulations for the computation of the stereotactic transformation matrix ( ), which is defined on the left-hand side of Equation 4, using three or more localizers. MC simulation computes the maximum deviation of all cases. The ideal transformation matrix ( ) and the perturbed MC computed transformation matrix ( ) were used with Equation 6 to compute a 3D target position, ( ) and ( ), respectively. Finally, the Root Mean Square Error (RMS-e) for target positions was computed as Equation 8 .

Overdetermined system of equations to compute the stereotactic transformation matrix
In the analysis of Equation 4, an overdetermined system of equations can be used to compute the stereotactic transformation matrix ( ). Then, analyzing the N-localizer system, one can identify several points via the vertical rods whose x-(lateral or LAT axis) and y-(antero-posterior or AP axis) coordinates are established by the manufactured geometry. Therefore, these vertical rod positions can be utilized in the same matrix computation for to minimize errors related to the x-and y-axes. Some localizer systems also contain a fixed x-or y-axis along the diagonal rods, but this geometry is not universal and was not utilized in the analyses. We analyzed the localizer frames (LF) of two stereotactic systems:   The Root Mean Square error (RMS-e) for the average target positions within 50 mm from the center of the volume in each direction (center, lateral, anterior, and posterior) is shown in Table 1 for four different Nlocalizer configurations having various geometries, where each configuration comprises three Nlocalizers. RMS-e is clearly reduced in every configuration by the inclusion of more data to specify an overdetermined system of equations. For three rods, the equations are determined as the standard approach, whereas the use of six (three diagonals plus three verticals) and nine rods (three diagonals plus six vertical) creates an overdetermined system of equations for computing the lateral (x-axis) and anteroposterior (yaxis) components of the transformation matrix ( ), as in Equations 1 and 2. The overdetermined system does not affect the vertical (z-axis), which is determined by three points.

Localizer Frame Width (mm) Height (mm) Width/Height 3 rods (RMS-e) 6 rods (RMS-e) 9 rods (RMS-e)
Leksell-G 120 120  Table 1 demonstrates that the BRW LF is the most accurate and stable three N-localizer system among the four systems analyzed. A prior study on the accuracy of the BRW LF was also reported by Grunert in 1999 [7]. One explanation for this finding is that the overall volume that the N-localizer system encloses is larger, not only in the vertical dimension but also within each axial CT image. Enclosing more pixels within the image leads to more pixels per millimeter, minimizing errors via greater precision. Another important finding is the standard deviations of the targeting errors in each system.

Localization areas
To better understand where the RMS errors are occurring, we fix points 50 mm in each direction along the m M M M center, lateral, anterior, and posterior aspects of the entire volume. The results of MC simulations at these points are plotted for both the BRW LF and LG LF ( Figure 2) and demonstrate that the BRW LF is relatively stable throughout the 50 mm space, but the LG LF has variable accuracy depending on the target location. The highest errors for LG LF are located at the inferior and posterior regions. The apparent differences as depicted in Figure 2 between the BRW LF and LG LF may be related to an area of a triangle formed by the three-diagonal fiducials for each image plane. The idea of keeping the target point inside of this triangle was first proposed by Perry, et al. [8] and subsequently tested by Grunert, et al. [7] and Brown [3]. For the BRW LF, the area of such a triangle is fairly constant along the height (z) but for the LG LF, the area of the triangle varies significantly depending on the height ( Figure 3). The largest predicted RMS-e for the LG LF occurs when a target lies posteriorly and inferiorly because for an inferior CT image, the triangle lies anteriorly and has a smaller area. Under these conditions, the target point has a higher probability of lying outside of the triangle. Heat maps can be generated using barycentric interpolation computed via Equation 36 in the Appendix, which are observed in Figure 3 for the BRW LF and LG LF [9][10]. For a more detailed analysis of these presumed errors related to targeting and the area of the triangle, please see the Appendix where two proofs demonstrate similar outcomes and introduce the term "stereotactic triangle". Given the above analysis of the areas of a triangle formed by the three diagonal rods, a simple modification to the two lateral N-localizer plates of the LG LF appears to be desirable. Specifically, the diagonal rods of these two N-localizer plates descend from posterior to anterior parallel to each other, as depicted in Figure  1. If instead, the right lateral N-localizer were inverted such that its diagonal rod descended from anterior to posterior, then all three diagonal rods would descend in a clockwise manner, similar to how the three diagonal rods of the BRW LF ascend in a clockwise manner. This modification of the LG LF would render the triangle area less sensitive to the height of the CT image, as is the case for the BRW LF. To explore the advantage of this inverted or antiparallel configuration, we evaluated a rectilinear 140x140 mm N-localizer system that comprises three N-localizers arranged in either standard (parallel) or inverted (antiparallel) configuration. The RMS-e for an overdetermined system of equations that exploits nine rods is 0.584 (+/-0.123) mm for the standard configuration and 0.488 (+/-0.093) mm for the inverted configuration.
The effect of the antiparallel design of the three N-Localizer system on the area of the stereotactic triangle was next analyzed. The parallel N-localizer system, such as with the Leksell G (LG) Localizer Frame (LF) (Leksell G, Leksell Stereotactic Frame, Elekta, Stockholm, Sweden) exhibits a steady decline in the areas of triangles normal to the vertical axis (z) as presented in Figure 3. Indeed, the antiparallel arrangement creates more uniformity of the areas of triangles normal to the vertical axis ( Figure 4).

Width versus height geometry of N-localizer systems
Monte Carlo analysis for RMS-e of the target can be performed using various sizes of N-localizer systems. For the BRW LF, we varied the vertical height subtended by the three N-localizers from 189 mm to 120 mm and used a nine-rod overdetermined system of equations to understand how height may affect the RMS-e. As expected, the average RMS-e improved from 0.51 (+/-0.013) mm to 0.37 (+/-0.008) mm as the height decreased ( Table 2).   Conceptually, this effect can be explained by considering the interpolation along the N-localizer's diagonal rod [4]. The steeper the diagonal rod, the greater the extent to which random noise will affect the interpolation in an interpolant. Therefore, the greater the width/height ratio of the N-localizer, the lower the RMS-e within that space. However, it should be noted that a considerable degree of height is necessary for localization throughout the entire stereotactic surgical space. Therefore, a balance between the practical limits to optimizing the width/height ratio and the encompassing space need to be considered.

N-localizer systems with four N-localizers
An overdetermined system of equations minimizes errors in computing the transformation matrix ( ). Previously, this overdetermined system had an influence on the x-and y-axes. However, the addition of another N-localizer adds more data for the x-axis, the y-axis, and, more importantly, the z-axis. Hence, incorporation of a fourth N-localizer to an N-localizer system that comprises three N-localizers may afford further improvement in RMS-e [3]. To explore this possibility, we used MC to analyze a 140x140 mm (WxH) rectilinear N-localizer system arranged in a standard configuration without inverting a diagonal rod. An incremental reduction in error is observed which is promoted by the incorporation of an overdetermined system of equations, here starting with four rods but reaching a total of twelve rod positions ( Table 3). These errors may be compared to the error shown in Table 1

Discussion
The N-localizer is a critical component of most frame-based stereotactic neurosurgical procedures, which are often used for precise targeting such as for deep brain stimulation (DBS), brain biopsies, invasive electroencephalography (iEEG), and laser interstitial thermal therapy (LITT). The N-localizer allows every axial slice in the localizer imaging space to be localized precisely, which can then be used during the neurosurgical procedure. Despite the important use of the N-localizer, there are no focused publications that analyze various optimizations of the N-localizer that can be implemented. Herein, we report the first such approach to improve the technology.
Since every image will have some associated noise that relatively uniform, using a system of overdetermined equations minimizes some of the errors during the stereotactic localization step. Introduction of an error during this step, therefore, may propagate into the actual surgical procedure, which may also have its own errors related to image fusion, imprecision from manually setting frame coordinates, and displacements related to imperfect instruments and deflections along a trajectory. Therefore, these targeting errors can be additive. Further, in comparing target accuracy for the CRW system versus the Leksell G (LG), average errors reported are 1.94 (+/-0.41) versus 2.35 (+/-1.03), respectively [11]. The larger average error for the LG relative to the CRW correlates well with the larger RMS-e predicted via MC for the LG Localizer Frame (LF) relative to the BRW LF (see Figure 1). While there are several factors and different methods that may affect these reported results, optimizing precision during the localization step may improve targeting for the neurosurgical procedure. In addition, techniques for improving Stereotactic Intraoperative Localization (StIL) and navigating coordinate systems have been previously reported [12][13][14].
There are several forms of optimization that play a role in error minimization. One such example is the novel utilization of all the vertical rod positions for their associated x-(LAT) and y-(AP) axes. As such, an overdetermined set of equations significantly reduces final computation errors for a target. This method is particularly attractive as it can be immediately performed on all N-localizer systems.
Other improvements in localization include optimization of the geometry of the N-localizer. The area of a triangle formed by three diagonal rods relative to a target can also play a role in errors. Targeting within this stereotactic triangle may be key to minimizing the errors. This observation leads to relatively simple improvements such as inverting opposite diagonal rods (antiparallel) and considering larger overall structures. There will be CT/MRI limitations related to size and perhaps inhomogeneity further away from the center of the bore, which would require further study. Further, the slope of the diagonal bar can play a role in errors as is exemplified by observing the width versus height of the N-localizer system.
While many systems primarily utilize three N-localizer registration, higher precision appears to be gained by utilizing a four N-localizer system [3]. For square-or rectangular-shaped three N-localizer systems, generally, the plates used are the lateral and anterior ones to avoid the posterior plate, which may be distorted during a CT or MRI study, or have materials that impair proper rod localization. Therefore, routine use of systems to support/suspend the stereotactic frames and suppress movement artifact optimize imaging conditions. Care during image acquisition cannot be overemphasized ( Figure 5).

Conclusions
Monte Carlo simulations of Root Mean Square error (RMS-e) is a useful technique in understanding targeting while using N-localizer systems in stereotactic neurosurgery. The authors have described many novel techniques that could lead to immediate improvements in the localization step. Some of these improvements include mathematical techniques and others are geometrical alterations of the systems. Utilizing all the fiducials of the N-localizer system and keeping the target within a stereotactic triangle may have clinical importance. Further investigation is needed to confirm these findings.

Appendices The Stereotactic triangle
As described in the technical report, the 3 diagonal N-localizer fiducials encode three 3D positions in addition to three 2D positions. Unlike the vertical bars, these diagonal bars are the only locations in 3D space that provide a value for the z-axis (vertical, z). Therefore, study of errors related to the z is the focus of the two proofs, which support the evidence computed via Monte Carlo simulations provided in the technical report.
Let us consider an axial computed tomography (CT) slice as a Euclidean plane containing 3 N-localizer fiducials. From those 3 N-localizer fiducials, we have 3 points ( ) at the diagonals. From these data, a target point (TAR, ) will need to be computed. Computing the x (lateral, LAT) and y (anteroposterior, AP) coordinates of the target will not be discussed here, because there are better methods available, as discussed in the technical report. For each diagonal point ( ), their associated values for z carry random errors from noise. The three points ( ) form a triangle and the TAR in this case is presumed to lie outside of the triangle ( Figure 6). has an unknown z value. Taking the diagonal rod farthest away from the target point ( ), is the intersection point between the line -and the line -. The stereotactic target may be inside or outside of the triangle area ( Figure 6). We introduce the term, the "stereotactic triangle", to refer to this triangle for 3 N-localizer systems. We stipulate, therefore, that to ensure accurate targeting, a target should be within the area of this triangle and that errors of larger and larger magnitudes may arise as the target moves farther and farther outside of the triangle. , , For the two subsequent proofs, we define a general triangle using 3 points to develop a function using partial derivatives which will describe the error [15,16]. We also define a target and a point that intersects two lines formed by and . Next, let us define areas and vectors related to the triangle (Figure 7).

Proof 1: General interpretation of the relation between z error and a triangle area using Barycentric Coordinates and partial differential equations
We next setup two vectors from to (equation 9) and to (equation 10). The variable now relates the fractional distance of divided by (equation 11). The symbol := in the subsequent equations implies , , It can then be determined that there are subsections of the triangle that can be defined in fractions as , , (equations [12][13][14]. We can then compute the target position (equation 15). It then follows that the areas and the fractional areas can be rearranged and combined to form equations 16-25.
We now choose to parameterize further to develop a function . Taking equation 15, we now introduce equations 22, 23, and 25 into a combined form (equation 26). We now take partial differentials (equations 27-29). Rearranging equations 27-29, taking their squares, we can then combine them to form equation 30. We make the assumption that the noise in the function is homogenous (equation 31). Finally, defining a parameter with boundaries (equation 33-34), then we arrive at equation 35. The function which may describe an error is then defined in equation 36.