3.1. Projection pattern reconstruction

First, the camera image is converted from RGB to grayscale. Then, a Gaussian and a median filter are applied to reduce the noise. The result is converted into a binary image using a threshold obtained with the Otsu binarization method [Reference Otsu39]. Afterward, the ellipse fitted is used to reconstruct the shape of the spotlight projection. An example for each step is depicted in Fig. 3.

Figure 3.

The original image (a) is converted to a grayscale image (b). A Gaussian and median filter (c) are used reduce noise. The result is converted into a binary image (d) that is used for the contour detection. A closeup of the detected contour (green) on top of the original image can be seen in (e).

The camera projection of the intraocular surface onto the image plane can be described using the pinhole camera model. Based on the camera model and the surface shape (simplified to be perfectly spherical), we can reconstruct the three-dimensional projection directly from the microscope image. The setting for the reconstruction is depicted in Fig. 4.

Figure 4.

Cross-section along the optical axis. Used for the reconstruction of the projection.

Using a point $p_c$ on the camera sensor and the focal point, we can define a line $l$ that intersects the surface of the sphere at the point $p_s$ . By using a cross-section containing $p_s$ , the focal point (F), and the center of the sphere, the problem can be simplified to an intersection between $l$ (yellow in Fig. 4) and a circle. The line $l$ is given by Eq. (1) and the circle by Eq. (2), where $f$ is the focal length, $r$ is the radius of the sphere, and $d_0$ is the distance between the focal point and the bottom of the sphere. $y_1$ is the Euclidean distance between the center of the camera sensor and $p_c$ . Here, the coordinate system is defined with the center of the sphere as shown in Fig. 4.

(1) \begin{equation} y=\frac{f}{y_1}x+(d_0-r) \end{equation}

(2) \begin{equation} x^2+y^2=r^2 \end{equation}

The said intersection allows to calculate the distance ( $d$ ) between the point on the sphere surface $p_s$ represented by $p_1$ and the optical axis.

The resulting Eq. (3) is based on the quadratic formula used to derive the intersection between the line and the circle.

(3) \begin{equation} d=\frac{\frac{2 d_0 f}{y_1}-\frac{2 r f}{y_1}+\sqrt{\left(\frac{2 r f}{y_1}+\frac{2 d_0 f}{y_1}\right)^2-4 \left(1+\frac{f^2}{y_1^2}\right) (d_0^2-2 d_0 r)}}{2 \left(1+\frac{f^2}{y_1^2}\right)} \end{equation}

Knowing the distance ( $d$ ) between the optical axis and the point $p_s$ , we can obtain the corresponding height ( $h$ ) using Eq. (4). The height is defined as the distance between $p_s$ and the bottom of the sphere along axis $Y$ , as depicted in Fig. 4.

(4) \begin{equation} h=r-\sqrt{r^2-d^2} \end{equation}

Given the position $p_c=(x_c,y_c)$ and the distance $d$ , we can calculate the estimated position $p_s=(x_s,y_s,z_s)$ , given by Eqs. (5), (6), and (7). Here, $s$ is the physical size of the camera sensor in mm and $p$ is the resolution of the image sensor.

(5) \begin{equation} x_s= \left \{ \begin{matrix} \dfrac{d}{\sqrt{1+\frac{y_c^2}{x_c^2}}},& \text{if } x_c\geq 0\\ \\[-5pt] \dfrac{-d}{\sqrt{1+\frac{y_c^2}{x_c^2}}}, & \text{otherwise} \end{matrix} \right. \end{equation}

(6) \begin{equation} y_s=\dfrac{y_p}{x_p}x_s \end{equation}

(7) \begin{equation} z_s=h \end{equation}

This allows to fully reconstruct the three-dimensional contour of the intersection based on the shape on the camera sensor.

3.2. Cone-sphere intersection

The intersection between a cone and a sphere is a rather complicated three-dimensional curve that does not lie on a two-dimensional plane. The only exception is the special case, where the center of the sphere lies on the axis of the cone, producing a circle-shaped intersection.

A parametric equation for this curve can be derived using equations defining a sphere and a cone. A right circular cone with the vertex in the origin can be defined using Eq. (8), where $\beta$ is the opening angle of the spotlight. $\beta$ is defined as the angle between the axis of the cone and every line from the vertex to a point on its surface. The axis is equal to the $Z$ axis.

(8) \begin{equation} z^2=\frac{x^2+y^2}{\tan\!(\beta )^2} \end{equation}

A sphere can be defined using Eq. (9), where $r$ is the radius of the sphere and $(x_0, y_0, z_0)$ is the position of the center.

(9) \begin{equation} (x-x_0)^2+(y-y_0)^2+(z-z_0)^2=r^2 \end{equation}

To set $y_0 =0$ and simplify the equation of the intersection, we can rotate the coordinate system around the axis of the cone, so that the center of the sphere is in the plane defined by the Z and X-axis. This does not lead to a loss of generality, as the cone is not affected by the rotation. The resulting equation used for the sphere is given in Eq. (10).

(10) \begin{equation} (x-x_0)^2+y^2+(z-z_0)^2=r^2 \end{equation}

Figure 5.

(a) Intersection between a sphere and a cone. (b) The blue plane is used to show the similarity between the cone-sphere intersection and an ellipse.

We can then obtain an equation for the intersection by combining Eqs. (8) and (10). The resulting equation is parametric with $x_i=x$ as a parameter. The points $p_i=(x_i,y_i,z_i)$ of the intersection can be calculated with Eqs. (11) and (12), where $c=\tan\!(\beta )$ . The range of values for $x$ is given in Eq. (13) and can be calculated using Eqs. (8) and (9). The definitions for $x_1$ and $x_2$ are given in Eqs. (14) and (15).

(11) \begin{equation} z_i=\frac{z_0+\sqrt{z_0^2-(1+c^2)(x_0^2-r^2+z_0^2-2x_0x)}}{(1+c^2)} \end{equation}

(12) \begin{equation} y_i=\pm \sqrt{z_i^2c^2-x^2} \end{equation}

(13) \begin{equation} x = [x_1,x_2] \end{equation}

(14) \begin{equation} x_1=\frac{\left(x_0+\frac{z_0}{c}\right)+\sqrt{\left(x_0+\frac{z_0}{c}\right)^2-1-\frac{1}{c^2} (x_0^2+z_0^2-r^2)}}{1+\frac{1}{c^2}} \end{equation}

(15) \begin{equation} x_2=\frac{\left(x_0-\frac{z_0}{c}\right)-\sqrt{\left(x_0-\frac{z_0}{c}\right)^2-1-\frac{1}{c^2} (x_0^2+z_0^2-r^2)}}{1+\frac{1}{c^2}} \end{equation}

3.3. Cone-plane intersection

When inspecting the shape of a cone and plane intersection in three dimensions, it is very similar to an ellipse. An example is depicted in Fig. 5(a). This similarity motivates us to simplify the real intersection to the shape of an ellipse, as this can significantly reduce the localization effort. Because instead of trying to reconstruct the location of the light source based on a projection of a three-dimensional curve, we can reconstruct based on the projection of an ellipse. It is known that the intersection between a cone and a plane has the shape of an ellipse, if the angle between the axis of the cone and the plane is higher than the opening angle of the cone. To show the similarity between the cone-sphere intersection to an ellipse, we construct a plane $P$ that intersects the cone and therefore produces an ellipse-shaped intersection.

The cone-plane intersection should be close to the cone-sphere intersection. First, we take the two points (A ( $x_1$ , 0, $z_A$ ) and B ( $x_2$ , 0, $z_B$ )) on the cone-sphere intersection with the biggest distance between each other and connect them with a line. The resulting plane $P$ is made perpendicular to the $XOY$ plane and contains this line. The x-coordinates of these two points are the ends of the range for the values of $x$ and can be calculated using Eqs. (14) and (15). The Z-coordinates of A and B are calculated using the cone equation as shown in Eq. (16). Figure 5(b) depicts an example for the constructed plane $P$ (blue).

(16) \begin{equation} z_A=\frac{x_1}{c}, z_B=\frac{-x_2}{c} \end{equation}

The resulting plane $P$ is defined by Eq. (17).

(17) \begin{equation} P\;:\; z=\frac{z_A-z_B}{c(z_A+z_B)}x+z_1-\frac{x_1(z_A-z_B)}{c(z_A+z_B)} \end{equation}

Table I.

The maximum difference between the intersections.

Figure 6.

Triangle (green) used to derive the vertex position based on a given ellipse. $S$ denotes the spotlight source. $C$ denotes center of ellipse. $s_1$ is the distance between $S'$ and $C$ . $s_2$ is the height of the spotlight source $S$ with the $XOY$ plane.

This equation can be used in combination with Eq. (8) to obtain the cone-plane intersection points $p'_{\!\!i}=(x'_{\!\!i},y'_{\!\!i},z'_{\!\!i})$ as a parametric ( $x'_{\!\!i}=x$ ) equation. $y'_{\!\!i}$ and $z'_{\!\!i}$ are given by Eqs. (18) and (19). Due to the definition of the plane, the range of values for $x$ is also given by Eq. (13).

(18) \begin{equation} y'_{\!\!i}=\pm \sqrt{\left(\frac{z_A-z_B}{c(z_A+z_B)}x-z_A-\frac{x_1(z_A-z_B)}{c(z_A+z_B)}\right)^2c^2-x^2} \end{equation}

(19) \begin{equation} z'_{\!\!i}=\sqrt{\frac{x^2+y_i^{\prime 2}}{c^2}} \end{equation}

From the definition of the plane $P$ , we know that the intersections have the two points A and B in common. To find the maximum difference between these two intersections, we can use the points $(x_h,0,z_h)$ where the surfaces intersecting the cone have their biggest difference. $z_h$ and $x_h$ are defined in Eqs. (20) and (21).

(20) \begin{equation} z_h=\cos\!\left(\sin^{-1}\left(\cos\!(\beta )\frac{x_0}{r}\right)\right)r+z_0 \end{equation}

(21) \begin{equation} x_h= \left \{ \begin{matrix} \sqrt{r^2-(z_h-z_0)^2}+x_0,& \text{if } x_0\geq 0\\\\ -\sqrt{r^2-(z_h-z_0)^2}+x_0, & \text{otherwise} \end{matrix} \right. \end{equation}

This allows us to directly calculate the maximum difference between the two intersections (the real intersection and simplified ellipse) by using the parametric equations. For our use case, we define an area of interest which is shown in Fig. 3. It is a 10 mm × 10 mm × 5 mm range, mainly defined by the microscope view and surgical region. The opening angle $\beta$ of the cone is independent of the instrument location and depends only on the designed spotlight’s lens. Therefore, we can calculate the maximum difference for different values for $\beta$ . This gives guidance on which angles could be suitable in regard to an error tolerance. The resulting maximum differences are given in Table I.

For our use case, these differences are negligible for the listed opening angles and we can see the projection as an ellipse without introducing a significant error ( $\lt\!\lt$ 10 μm).

3.4. Ellipse to cone reconstruction

For the ellipse fitting, the reconstructed shape of the contour is rotated onto the $XOY$ plane as shown in Fig. 6. After reconstructing the vertex of the cone, the inverse rotations are applied. The ellipse can be defined using the position of its center, the length of the major axis $a$ , and minor axis $b$ . The size of the minor axis is related to the distance between the vertex of the cone and the plane. The relationship between $a$ and $b$ depends on the angle between the cone and the $XOY$ plane.

Table II.

Properties of the camera in the Blender simulation.

Table III.

Rotations applied to each spotlight to achieve the correct angles.

To find the vertex position, a right triangle is used as depicted in Fig. 6. One corner of the triangle is the vertex position, and one corner is in the center position. The side $s_2$ is perpendicular to the XOY plane, and the side $s_1$ follows the major axis. The length of the side $s_1$ and $s_2$ can be calculated using Eqs. (22), (23), and (24),

(22) \begin{equation} \alpha =\sin^{-1}\left(\sqrt{1-\frac{b^2}{a^2}}\cos\!(\beta )\right) \end{equation}

(23) \begin{equation} s_1=a\frac{\sin\!(2\alpha )}{\sin\!(2\beta )} \end{equation}

(24) \begin{equation} s_2=a\left(\frac{\cos\!(2\alpha )}{\sin\!(2\beta )}+\frac{1}{\tan\!(2\beta )}\right) \end{equation}

where $\alpha$ is the angle between $s_2$ and the hypotenuse of the triangle $SC$ . As the ellipse lies in the $XOY$ plane. The spotlight position in $XYZ$ is defined as $p_l=(x_l,y_l,z_l)$ shown in Fig. 6. The $x_l$ and $y_l$ can be calculated using the rotation of the ellipse and $s_1$ . $z_l$ equals to the length of $s_2$ . Due to the symmetry of the ellipse, two possible positions for the vertex exist. Knowing the rough position of the insertion point allows us to narrow it down to one position. To derive the final result, the inverse rotations have to be applied to the vertex position.

Figure 7.

(a) Path of the spotlight used during the first test (Box). (b) Path of the spotlight used during the second test (Helix).

Figure 8.

The error performance with the box trajectory for a single spotlight. $e_x$ , $e_y$ , and $e_z$ denote the error in $X_S$ , $Y_S$ , and $Z_S$ axis shown in Fig. 4, respectively. $\epsilon _x$ and $\epsilon _y$ denote the rotation error in $Y_S$ and $Z_S$ axis. (a) Positioning error in each direction $X_S$ , $Y_S$ , and $Z_S$ . (b) Orientation error in $Y_S$ and $Z_S$ . (c) Overall error in position. (d) Overall error of orientation.

Figure 9.

The error performance with the box trajectory for multiple spotlights. (a) Positioning error in each direction $X_S$ , $Y_S$ , and $Z_S$ . (b) Orientation error in $Y_S$ and $Z_S$ . (c) Overall error in position. (d) Overall error of orientation.

3.5. Multiple spotlights

As the single spotlight may be prone to errors, we further analyze a setup with multiple spotlights. To evaluate the performance of such a setup, an instrument with three attached spotlights is evaluated. For each projection, the possible vertex positions are reconstructed independently following the algorithm for the single spotlight. To choose the resulting position, all possible combinations including three positions are evaluated based on their spatial difference. The set with the lowest difference is selected. From these three positions, the median position is selected as the final result. The workflow is depicted in Fig. 1(e–g).

Table IV.

Error for different tests.

Figure 10.

The error performance with different light intensities. The whiskers show the minimum and maximum recorded distance changes. The start and end of the boxes denote the first and third quartile. The band, red dot, and cross represent the median, mean, and outliers of the recorded changes, respectively. (a) Positioning error with the square trajectory for a single spotlight. (b) Rotation error with the square trajectory for a single spotlight. (a) Positioning error with the helix trajectory for a single spotlight. (d) Rotation error with the helix trajectory for a single spotlight.

Figure 11.

The error performance with different light intensities. (a) Positioning error with the square trajectory for multiple spotlights. (b) Rotation error with the square trajectory for multiple spotlights. (a) Positioning error with the helix trajectory for multiple spotlights. (d) Rotation error with helix trajectory for multiple spotlights.