Calibration method for a breast intervention robot based on four-dimensional ultrasound image guidance

In breast interventional ultrasound therapy, it is difficult to directly diagnose the location of a tumor in 2-D ultrasound images. To assist surgeons in treatment more intuitively, a four-dimensional ultrasound image-guided breast intervention robot is proposed. The calibration approach of the ultrasonic image for the robot is one of the main contents of the research. This method is based on the establishment of a complete coordinate system conversion model, and it uses the ORB (oriented FAST and rotated BRIEF) feature extraction method to obtain and record the real-time image marker pixel positions, calculate the unknown parameters of the coordinate system conversion matrix, and establish a complete calibration system. This article demonstrates the feasibility of the calibration approach through experiments in our developed US-guided robotic system. Additional experimental and parametrical comparisons of the proposed method with state-of-the-art methods were conducted to thoroughly evaluate the outperformance of the proposed method.


Introduction
In recent years, with the continuous advancement of medical ultrasound (US) imaging technology and its outstanding advantages, US images have become the main diagnostic and treatment tool for breast lesions. US image-guided surgery has been used successfully in the auxiliary diagnosis and treatment of breast cancer, which is usually called interventional ultrasound (Solbiati & Tondolo, 2013). It is difficult for the surgeons to directly diagnose the location of a tumor in 2-D US images. Therefore it is necessary to reconstruct breast tumors to overcome the limitations of 2-D US imaging (Papenberg et al., 2008). Due to the shape of breast is prone to change in real-time during the surgery, the mapping of the US image frame and the patient frame changes in real-time, a reconstruction of the 4-D breast model, which shows the 3-D breast model in a moving state, is more suitable for clinical application of breast US intervention therapy.
To implement of breast US intervention treatment guided by 4-D US image, it is first necessary to convert the imaging frame to the surgical frame (original coordinate system), which is called the calibration of the US image. In the past ten years, scholars around the world have successfully researched calibration approaches and the corresponding optimization algorithms. Mostly, the US image calibration requires a phantom to provide a set of fiducials, e.g. the Nwire phantom method (Carbajal et al., 2013) and multi-wedge phantom method (Najafi et al., 2014), they are traditional US calibration methods. With the development of surgical robots, the traditional calibration methods no longer meet the clinical application, therefore some US image calibra-tion methods for robotic systems are developed. Huang et al. presented a robotic calibration method with a ball-shaped phantom (Huang et al., 2018). The calibration parameters were estimated with robotic translational motion, as while as rotational parameters were not taken into account in calibration. So this approach could not be applied in general US scanning from all angles of view. Aalamifar et al. designed an active echo phantom to transmit and receive US signals to a US probe for image calibration (Aalamifar et al., 2016). This approach designed a complex scanning trajectory to align the probe with the active echo phantom, without using an external tracking system. Wang et al. proposed a calibration method for the multiple cooperative robots (Wang et al., 2014). The calibration problem for multiple robot cooperation needs to solve a matrix equation with a priori knowledge of the temporal correspondence of the data, that are not always easy or possible to be obtained. All of the above are mechanism-based calibration methods. After integrating with the optical tracking system, the US system could interact with the robot arm and realize a fully automatic calibration procedure. Li et al. (2012) presented an automatic probe calibration approach based on rigid ultrasound image registration of 2D rotational slices for 3D ultrasound image reconstruction. The parameters was determined with an automatic image registration based approach. Xiong et al. (2021) presented a mechanism-image fusion-based calibration approach to a US-guided dual-arm robot. The pixel position in the US image of the needle-tip is calculated by manual annotation, which reduces the efficiency of registration. Due to the ultrasonic noises and artifacts, calibration accuracy and efficiency still have an apparent influence on navigation. Aiming the above influence from ultrasonic noise and artifacts, Wang et al. (2018) proposed an automatic probe calibration method based on the coherent point drift (CPD) point cloud registration. During calibration procedure, a large number of ultrasonic images are collected, each image frame has its own tracking information. It generates a lot of data which is time-consuming and results in a low computational efficiency.
In image-based calibrations mentioned above, different image processing pipelines were introduced to segment the target point in US images, the image processing is often complex and time-consuming. Moreover, these US calibration approaches ignored the deformation of soft tissue, they were all based on rigid US image registration. In our research, A US calibration approach of 4-D US-guided breast intervention robot with easy operability, low redundancy, and the ability to predict soft tissue deformation is proposed. A calibration phantom with special geometries and functions is designed. The markers can not only easily be fixed on the skin surface of the breast by the phantom; the spatial positions of the markers on the skin surface are but also calibrated by the phantom. A mathematical model of the coordinate system transformation between the various subsystems of the breast intervention robot system was established, so the calibration problem for 4-D US-guided breast intervention robot was formulated as a fundamental problem of solving a matrix equation. The calibration uses the ORB (oriented FAST and rotated BRIEF) (Rublee et al., 2011) feature extraction method to obtain and record the real-time image marker pixel positions, and calculates the parameters of conversion matrix by PSO (Particle Swarm Optimization) algorithm (Okwuand & Tartibu, 2021). Finally, the effectiveness of the proposed calibration approach was validated through experiments in our developed US-guided robotic system. Additional experimental and parametrical comparisons of the proposed method with state-of-the-art methods were conducted to thoroughly evaluate the outperformance of the proposed method.

Overall frame of the calibration approach
In a 4-D US image reconstruction system, the robot-holding US probe calibration system is as shown in Fig. 1. It consists of B-model US image equipment, an optical navigation system (an NDI optical navigator, passive rigid body, and passive 4-point probe), a 4-D dynamic calibration phantom (a calibration phantom holder and markers), a robot holds an US probe, a patient and an operating table. There is no direct connection among them in the breast intervention robot system. The NDI optical navigation is applied to the US probe calibration of breast intervention robot. The patient frame and the US probe frame in the calibration system are associated through NDI optical navigation, which lays the foundation for the real-time guidance of the US probe during robot operation and the dense reconstruction of the 3-D breast tissue.  For the proposed calibration approach, a 4-D dynamic calibration phantom was designed and applied to the US image calibration for a breast intervention robot. The functions and components of the 4-D dynamic calibration phantom bracket include a marker positioning bracket with marking holes, a coordinate system conversion bracket, an adjustable screw and nut, and a double-sided suction cup, as shown in Fig. 2. The coordinate system of the patient's breast is defined as the virtual coordinate system in the calibration approach, and the center of the nipple is set as the origin of the virtual coordinate system. The positions of the marking holes on the positioning bracket are evenly arranged on two flat arc curves, and the marking holes are processed by precision instruments. The parameters of the coordinate conversion between the conversion bracket and the virtual coordinate system are determined by the size and structure of the conversion bracket. The red double arrow indicates the direction in which the dynamic phantom can be moved and adjusted, and the two red curves indicate the locations of the markers.
Before calibration, the registration of US probe and operating table in the optical navigation system is completed, and the mappings between US probe, operating table and optical navigation system are established. The specific procedures for registration are as follows: (1) The passive rigid body is fixed at the center of US probe, and the optical positioner tracks its position and posture. The coordinate system of US probe is set as O S , the coordinates origin is the center of US probe, and the XOY plane is parallel to US imaging plane. The coordinate system of US image is set as O I , the coordinates origin is at the upper left corner of US image, and the UOV plane coincides with US image plane. Therefore, there is only translation without rotation between US image and US probe. (2) The NDI navigation system detects the coordinate origin of coordinate transformation bracket fixed on operating table through passive 4-point probe and records its position information. The operating table is placed horizontally, and the patient lying prone on it. Therefore, the XOY plane of operating table is parallel to horizontal plane, and there is only translation between operating table and optical navigator.
In the overall frame of US calibration, five coordinate systems are defined: the US image, the actual and virtual coordinates of the coordinate conversion bracket, the US probe, and the NDI optical navigator, which are presented as I , M , M, S, and R, respectively. The four corresponding coordinate transformation matrices are T S

Skin surface marker
The patient's breast surface is used as the carrier of the markers. Before the operation, the doctor adjusts the conversion bracket fit to the patient's position and shape so that the intersection of the curve on the marker positioning bracket coincides with the origin of the virtual coordinate system, then uses the marker positioning bracket to mark the skin surface of patient's breast, as shown in Fig. 3. The marker is a small steel ball with a diameter of 2 mm. A total of eight markers are divided into two groups, four in each group, and fixed on breast skin surface with glue. The sticking position of marker on breast surface is determined through graduation hole on the marker positioning bracket. Curves C 1 and C 2 intersect at right angles in space, and intersection point is O M . Intersection of two planes is set as the Z axis of virtual coordinate system. The virtual dynamic coordinate system changes in real time with the position and shape of breast during operation. In the US probe calibration, US probe scans these two cross-sections located by the two curves, as well as NDI optical navigator detects the position of US probe.
We construct a cylindrical coordinate system in the dynamic calibration phantom: O M is the coordinates origin, and Z axis is the center of cylindrical coordinate system. The coordinates of a point in the dynamic virtual rectangular Fig. 3 Markers on the skin surface of breast coordinate system of the breast can be expressed as P M = x m y m z m T . We convert P M from cylindrical coordinates to rectangular coordinates: P M = r cos ψ r sin ψ z T , and the initial angular coordinates ψ of markers on curves C 1 and C 2 are ψ 1 = 0 • and ψ 2 = 90 • , respectively. The coordinates of markers on curves C 1 and C 2 can be obtained as P M,1i = (r 1i , 0, z 1i ) and P M,2i = (0, r 2i , z 2i ), i = 0, 1, 2, 3.

Mathematical model for dynamic US calibration approach
The space mapping of a point from the US image to the patient's breast is expressed as follows: Matrix T M M is the transformation matrix from operating table frame to breast frame. The translation parameters are determined by the structure and size of coordinate conversion bracket. Suppose the translation parameters are (x 1 , y 1 , z 1 ) = (0, 100, 30) and the units are mm. Assuming that quaternion parameter of the rotation transformation is (w 1 , a 1 , b 1 , c 1 ) in T M M , the specific expression is as follows (Yi et al., 2009): The rotation parameters [w 1 a 1 b 1 c 1 ] T change in real time in the 4-D dynamic calibration phantom. According to conventional and intuitive method of expressing rotation, we introduce a rotation vector − → k = [x y z] T (t) and a rotation angle θ (t) , which represents the rotation angle of coordinate system M around axis vector − → k at time t . Finally, the rotation of breast coordinate system is represented as, In the process of calibration, the rotation angle θ changes in real time with breast movement and deformation, and the real-time dynamic rotation parameters θ , x, y, and z will be solved by PSO algorithm. Transformation matrix T M R represents translation transformation between NDI optical navigator frame and template support frame. We use passive 4-marker probe to detect the position of coordinate origin of phantom holder to obtain the translation x m y m z m 1 T . Then, the translation transformation matrix is, Conversion matrix T R S represents the coordinate system conversion between NDI optical navigator frame and US probe frame. The posture of passive rigid body recorded by the NDI is in form of quaternion; therefore, T R S has a total of 7 parameters, four of which are rotation parameters and three of which are displacement parameters, which are w 2 a 2 b 2 c 2 x 2 y 2 z 2 .
In the process of calibration, US probe scans the crosssections where the two curves are located. Therefore, NDI optical navigator records two sets of pose information of US probe, corresponding conversion matrices in these two states are represented as T R(1) S and T R(2) S . The transformation equation from US image frame to US probe frame is, Substituting P I and P M into (1), we obtain, Expanding the zero elements on the left side of the two linear matrix equations, we obtain two equations below, which represent the calibration models of two poses of US probe; f 1 and f 2 represent the calibration equations when US probe scans curves C 1 and C 2 , respectively.
The parameters of transformation matrix in the US calibration approach are θ , x, y, and z, which will be solved in our research.
If the cycle of compound movement of patient's breast, which moves with breathing, heart beating and artery motion, is T, then 0 ≤ t ≤ T . US image sequences I 1 and I 2 are acquired in this compound motion cycle to obtain original image data set of markers. ORB method is used to extract the feature points of markers from the two sets of image sequence and obtain the pixel positions of markers. At the same time, NDI optical navigator records the pose information T R(1) S and T R(2) S of two states of US probe and calculates the conversion matrix between US image and patient's breast frame through above mathematical model.

Solution of initial rotation parameters
In the proposed US image calibration, PSO algorithm is used to solve the dynamic rotation parameters. The main contents of this algorithm is summarized as follows: 1. The parameters θ , x, y, and z are defined as the dimensional variables x ji of PSO, where j = 1, 2, · · · , N ; i = 0, 1, 2, 3. The ranges of particle positions are 0 ≤ Initialize the single extremum p ji = x ji (0) of the particle, and global extremum p gi is determined by one particle which has the smallest fitness value among the N particles. 3. Evolution equations for particle velocity and position are, v ji (s + 1) = ωv ji (s) + c 1 r 1 (s) p ji (s) − x ji (s) x ji (s + 1) = x ji (s) + v ji (s + 1) 4. The fitness function of PSO algorithm for the US calibration is defined as, 5. Calculate the fitness F(x) of each particle, update the particle position with PSO iterative algorithms, select and update the global optimal position. 6. When does the iteration terminate is determined by ξ = F(x) −F(x) = |F(x)| ≤ 0.001.

Materials and equipments
In order to verify the accuracy of dynamic US calibration, experiments are completed on breast intervention robotic system proposed by us, and a self-designed dynamic calibration phantom is used to calibrate US image. All structural parts of phantom are printed using a fused deposition modeling printer, polylactic acid (PLA) is used as the materials. Use the phantom to paste two groups of markers on patient's breast skin surface, each group contain 4 ones, a total of 8 markers, and the processing accuracy of positioning holes and markers reaches 0.1mm. Measure the angle values of two groups of markers in the phantom frame with a universal bevel protractor. The positions of each markers in the phantom frame is recorded in Table 1. The B-mode ultrasound, type LOGIQ P7, is produced by GE with a supporting image capture function. It can collect and transmit US images to the image workstation in real time through a video transmission port.
The NDI Polaris Spectra optical navigator, made in Canada, is used to track the pose of the US probe and transmit the signal of the pose data collected by the NDI optical navigator to the graphics workstation as the input information of the optimization algorithm. Obtain the position of coordinate conversion bracket by NDI optical navigator, it's 277.31 −235.02 −2217.08 1 T , then substituted it into equation (4) to get the transformation matrix T R M . NDI optical navigator record the postures of US probe in the two direction are . The Dell Graphics Workstation is used to extract the feature of markers in US images and to run the optimization program and calibration application program.
In the master-slave control system of UR5 robot, the effector of the robot holds US probe, the master site of it is the Omage.6.

Image acquisition of markers on breast surface
Human breathing cycle is approximately 3 ∼ 3.7 seconds, and the breathing cycle of patient in this experiment is T = 3 seconds. In two consecutive breathing cycles, Omage.6 as the master site is used to control UR5 robot holding US probe to scan the cross-section of patient's breast and obtain US image sequences I 1 and I. The equipments, materials and processes used in the experiment are shown in Fig. 4.
The US image acquisition frequency is set to 30, and then 90 frames of serial images are acquired in period T. The US image of markers is shown in Fig. 5. The metal markers cause artifacts in the image, they are highlighted in the US image.

Position extraction of markers
The ORB method is used to extract the features of markers in this experiment, as ORB feature is currently a very representative real-time image feature. ORB maintains the   invariance of the rotation scale of the feature and has a significant increase in speed (Gupta et al., 2019). To completely extract the 4 markers from each sequence image, the number of ORB features is set to 10 in ORB algorithm. The coordinates of the highlights in the sequence image were extracted, there may be multiple very close highlights at one marker, remove them appropriately and calculate their average value. The result are shown in red and green circle marks, and the fitted feature curve are shown in Fig. 6, and their pixel positions are recorded in Table 2. , and P I into (7) and (8), then according to (9), a equation set containing 4 target parameters is obtained, and a fitness function, which is introduced in PSO algorithm, is constructed with this equation set. The particle population size is set to N = 50, and the evolution parameters of particle velocity and position are set to ω = 0.8, c 1 = c 2 = 0.5. Solve the target parameters through PSO and obtain: (2) and (3)  , and T S I into (1) to get the mapping from US image frame to patient's breast frame on direction of curve 1(C1), in a similar way, to get the mapping on direction of curve 2(C2).

Quantitative analysis
The metric used to quantify the precision of a calibration method is the calibration reproducibility (CR) (Carbajal et al., 2013).
Target registration error (TRE)s used to evaluate the accuracy of calibration (Pagoulatos et al., 2001).
Where P B is the position of test markers in calibration phantom frame.
In the experiment, select a frame image from the image sequence as test image, calculate calibration matrix T using the remaining frame images. The markers in test image is used as the registration target for testing. Quantitative analysis of the dynamic US image calibration is as shown in Table 3.
Compared with other calibration approaches, the comparison results are shown in Table 4.

Conclusion
In the experiment, ORB feature extraction method can quickly obtain the accurate pixel of markers in US images and there was no need to denoise the original US image before extracting the marker position. The pixel positions of a total of eight markers could be obtained on the two curves while ensuring that US image was not distorted, which provided data for US probe calibration of the breast intervention robot. The initial parameters (θ , x, y, z) of each marker curve in the calibration can be calculated preliminarily, which provides a foundation for parameter solution of dynamic US probe calibration.
This article demonstrates the feasibility of US probe calibration for a breast intervention robot through experiments and quantity analysis. Specifically, the CR and TRE of the whole system in our proposed method are only 0.83 and 0.89mm, respectively. Those values are far less than the CR (1.97 mm) and TRE (2.06 mm) of the N-wire phantom calibration method and the CR (1.58 mm) and TRE (1.80 mm) of the multi-wedge phantom calibration method. Our proposed 4-D dynamic calibration method is significantly better than the other two image-based methods in terms of the calibration accuracy of the whole system.