Finding the rigid transformation that aligns two 3D models is a fundamental problem in computer vision with applications in multiple fields, ranging from robotics
Finding the rigid transformation that aligns two 3D models is a fundamental problem in computer vision with applications in multiple fields, ranging from robotics[Pomerleau] to medicine [almhdie], and passing by augmented reality [Wu2013]. This article is motivated by medical applications in general and surgical navigation in orthopaedics in particular [Mezger]. The workflow of surgical navigation is usually a two step process. First, the surgeon uses a pre-operative 3D image of the targeted anatomy, e.g. a CT or MRI, to plan the procedure. Second, an intra-operative system performs optical tracking of fiducial markers attached to instruments and bones for determining their 3D pose in real-time, that are used to guide the surgical execution according to what was established in advance. For this purpose the system must overlay the plan with the actual patient’s anatomy in the OR, which passes by using a registration algorithm to align intra-operative 3D data with the pre-operative CT or MRI (Fig. 1).
Thus, 3D registration is a crucial step in surgical navigation and a common approach consists in asking the surgeon to pin-point a number of recognisable anatomical landmarks with a tracked touch-probe. These landmarks are reconstructed in 3D providing explicit point correspondences with the pre-operative image that enable alignment with classical methods [Horn88]. Unfortunately the approach is not effective in practice because the surgeon often struggles to recognise and access the landmarks, and invariably fails to touch their exact location, which strongly affects the accuracy and robustness of registration. A better alternative is to use as intra-operative 3D data the curves reconstructed by randomly grasping the bone surface with the tracked probe (Fig. 1). However, and to the best of our knowledge, there are no methods in literature to perform the global alignment of a curve and a surface. Stryker [Stryker] and Exactech [Exactech] have recently introduced navigation systems that employ randomly reconstructed curves, but they are exclusively used to refine registration using a local ICP variant [Gruen05, Besl92], and touching anatomic landmarks is still mandatory to perform initial alignment.
This article presents for the first time a global method for registering 3D curves with 3D surfaces without requiring a coarse initial alignment. The algorithm works with pairs of points augmented with local differential information that define the so-called 2-tuples point+vector with the vector being the tangent at the point in case of curves, or the normal at the point in case of surfaces. It is shown that the rigid transformation that aligns curve and surface can be determined in closed-form from a single pair of matching 2-tuples for which the two points correspond. In addition, a 2-tuple point+vector can be described in an invariant manner by a 4-parameter descriptor from which it is possible to derive a set of necessary conditions for a pair of 2-tuples to match. These findings enabled to devise a fast search scheme to establish putative 2-tuple correspondences between curve and surface that are used in an hypothesise-and-test framework to accomplish global registration.
Comparative tests against plausible alternatives in the literature demonstrate that the proposed algorithm is the first effective solution for global alignment of curves and surfaces. The experiments also show that the approach is well tailored for solving the 3D registration problem in surgical navigation, with the method proving to be very fast and robust, being able to accomplish accurate alignment in situations of small overlap and large percentage of outliers.
As a final contribution, the framework is extended to the case of registration of curve-vs-curve, which is also a largely unsolved problem, and surface-vs-surface [Raposo17].
1.1 Related work
Despite the vast literature in 3D registration, there is a limited number of authors explicitly addressing curve-vs-surface (Fig. 2(c)) and curve-vs-curve (Fig. 2(d)) alignments. Most works concern the alignment of sparse fiducial points (Fig. 2(a)) or dense point clouds that are referred to as surfaces (Fig. 2(b)).
There are several approaches for the local alignment of two point clouds of which Iterative Closest Point (ICP) is probably the most prominent and broadly disseminated one
There are several approaches for the local alignment of two point clouds of which Iterative Closest Point (ICP) is probably the most prominent and broadly disseminated one[Besl92]. There are numerous variants and modifications of ICP [almhdie, Pomerleau2013], many of which work well in situations of refinement of curves-vs-surface or curve-vs-curve registration [Gruen05]. However, all these methods are local and require proper initialisation to converge. In this article, we aim at global, fast alignment of curves and surfaces, with no prior assumptions about the initial displacement or amount of overlap.
If point correspondences are explicitly known, then the alignment can be accomplished by classical methods without the need of initialization [almhdie]. The recent Go-ICP [Yang16] and GOGMA [GOGMA] algorithms use Branch-and-Bound (BB) over the 6-dimensional space of euclidean motions to achieve global registration of Point Clouds without point correspondences. Since the complexity of BB is exponential in the dimension, these methods are slow, computationally expensive and, from our experiments in curve registration, they often diverge because of small overlap. Several authors propose to handle complexity by searching for rotation and translation separately [Makadia, Straub_2017_CVPR]. However, the search for the rotation is invariably performed in the space of the surface normals, which precludes the application to curve-vs-surface registration because of the lack of normals in the curve side.
The 3D surface feature-based algorithms [Rusu, FGR] involve extracting local features, obtaining matches between features in the two point clouds, and finally estimating the relative pose using RANSAC or other robust estimators. Since curves and surfaces have very different topologies, it is difficult in practice to detect common, coincident saliencies. Moreover many of these methods use feature description for matching which is typically designed for dense point clouds. We run comparative tests with the method of
involve extracting local features, obtaining matches between features in the two point clouds, and finally estimating the relative pose using RANSAC or other robust estimators. Since curves and surfaces have very different topologies, it is difficult in practice to detect common, coincident saliencies. Moreover many of these methods use feature description for matching which is typically designed for dense point clouds. We run comparative tests with the method of[FGR] and show that the approach is not amenable for curve-vs-surface registration.
The new family of algorithms 4PCS [amo_fpcs_sig_08, Mohamad15] replaces features by sets of 4 coplanar points whose relations define affine invariants that are preserved under rigid displacements. They work in hypothesize-and-test schemes by selecting a random base of 4 points in the source 3D model and finding all the 4-point sets in the target model that are approximately congruent with the base, i.e. related by a rigid transformation. Despite the search being in linear time, the approach is not suitable for performing curve vs surface alignment because of the very small overlap that dramatically increases search time, as shown by our experiments.
More closely related with our work is the article of [Gourdon94] that also uses local differential information for 3D registration of curves. Two methods are proposed: the first only requires a point correspondence between curves, which considerably decreases the complexity of search, but involves the computation of third order derivatives which is impractical in real, noisy data; the second uses two correspondences, leading to a dramatic increase in the complexity. We present a more clear mathematical formulation of the problem, and provide new insights that lead to an effective search scheme.
Notation: Matrices are represented by symbols in sans serif font, e.g. , vectors are represented by bold symbols, e.g. , and scalars are indicated by plain letters, e.g. . Normals and tangents are represented by lower case bold symbols and 3D points are written in upper case bold letters.
2 Curve vs surface registration using 2-tuples point+vector
This section presents a method for estimating the rigid transformation , with rotation and translation components and , respectively, that aligns a curve C with a surface S, as depicted in Fig. 3. For this purpose, we start by showing that it is possible to compute from a pair of corresponding points and , together with the information of their tangents on the curve side and their normals on the surface side. In the remainder of this paper, the pair of points with the corresponding tangents/normals will be referred to as a 2-tuple point+vector and all the tangents and normals in the mathematical derivations are assumed to be unitary.
This section also shows how a 2-tuple point+vector can be described in a compact, translation- and rotation-invariant manner by a 4-parameter descriptor , and provides the derivation of the necessary conditions for a 2-tuple point+tangent to be a match of a 2-tuple point+normal. These conditions are used in Section 3 to effectively establish putative matches that allow a fast 3D registration.
2.1 Closed-form solution for curve vs surface registration
Let and be two corresponding 2-tuples point+vector in curve C and surface S, respectively, and the rigid displacement that aligns C with S. Rotation can be determined independently of translation as the succession of two rotations: that aligns vectors and , and that places tangents in the planes defined by normals , respectively. This can be written as
where rotation is represented in angle-axis format by
with being the normal to the plane defined by vectors and , as illustrated in Fig. 4(a), and being given by , with .
Having vectors and aligned using rotation , a second rotation around by an angle (Fig. 4(b)) must be performed in order to make vectors and be orthogonal to and , i.e., must satisfy the following conditions
Using Rodrigues’ formula, can be written as
where is given by
Please note that matrix is not an arbitrary matrix. Its structure must be such that the first two values of its right-side null space are consistent sine and cosine values. This idea will be further explored in Section 2.3.
Given rotation , the translation can be determined in a straightforward manner using one of the point correspondences: .
2.2 Translation- and rotation-invariant descriptor of 2-tuples point+vector
At this point, it is possible to compute given matching 2-tuples between a curve and a surface. However, there is still the challenge of, given a 2-tuple in one side, finding potential correspondences on the other side. This section describes a compact description of a generic 2-tuple that will prove to be useful for carrying this search.
Let be two points and be the corresponding vectors that can either be tangents, in case belong to a curve, or normals, in case lie on a surface.
Consider a local reference frame with origin in , with the axis aligned with , and with the axis oriented such that it is coplanar with vector and points in the positive direction. This arrangement is depicted in Fig. 5(a), where , and .
The local cartesian coordinates can now be replaced by spherical coordinates which are particularly convenient to represent vectors. Choosing these coordinates such that the azimuth of vector is zero, it comes that the mapping from cartesian to spherical () coordinates is
where and .
The cartesian coordinates of vectors and in the local reference frame, expressed in terms of azimuth and elevation are
Equation 7 emphasizes an important fact that is that an appropriate choice of local frame allows to uniquely describe a 2-tuple point+vector up to translation and rotation using only 4 parameters, which are used to construct vector (Fig. 5(b)):
Further mathematical manipulation enables to directly move from a 2-tuple to its descriptor by applying the following vector formulas
where represents the signal function.
2.3 Necessary conditions for a 2-tuple in a curve to match a 2-tuple in a surface
Let and be 2-tuples in curve C and surface S with descriptors and , as defined in Equation 8. If the 2-tuples are not a match, the matrix equation 5 does not provide a solution with the desired format and rotation cannot be estimated. This section explores this fact to derive the necessary conditions for the pair of 2-tuples and to be a match by enforcing that Equation 5 has a consistent solution.
Let be defined by . The first condition for and to be a match is that . Another necessary condition is that there exists a rotation that simultaneously makes be orthogonal to . Since we are considering local reference frames for description such that and are coincident and aligned with a common axis, the system of equations 3 becomes
Since the cosine varies between and , the following must hold to enable the existence of an angle :
Manipulating the previous equations on the elevation angles of descriptors and , we obtain a set of inequalities that, together with the distance condition, are necessary conditions for the pair of 2-tuples to be a match:
A careful analysis of the inequalities shows that they are the conditions on the elevation angles for making the cone defined by rotating vector (or ) to intersect the plane defined by (or ). This is illustrated in Fig. 6, where the two possible situations of non-existence or existence of a valid rotation are represented. This figure also clarifies the fact that the orientation of tangents/normals is irrelevant since all the derivations are independent of such orientation.
The previous inequalities must be satisfied in order to exist a rotation such that becomes orthogonal to and to in separate. A condition on the azimuthal and elevation angles that makes the two pairs of vectors orthogonal in simultaneous can be obtained by manipulating Equation 11:
with and .
3 Method for fast curve vs surface registration
At this point, and given 2 corresponding 2-tuples, we are able to determine the rigid transformation . In addition, we proposed a way to describe each 2-tuple by a compact 4-parameter vector, with such description being invariant to translations and rotations, and derived the necessary conditions on these parameters for a 2-tuple in curve C to be a potential match of a 2-tuple in surface S. The current challenge is in quickly establishing the correspondences such that a fast alignment of the curve and the surface is obtained. This section proposes a solution to this problem.
A typical CAOS procedure has an offline stage for obtaining a 3D model of the targeted bone that occurs before the actual surgery is performed. Knowing this, we propose an offline stage for processing the bone model (surface) whose output is used in the online correspondence search scheme, allowing a very fast operation. The sequence of steps of this stage is shown in Fig. 7(a). The advantage of performing an offline processing of the data is that most of the computational effort of the algorithm is transferred from the online stage to the pre-processing, where computational time is irrelevant.
We propose to build a data tree structure that contains the relevant information for all pairs of points in order to facilitate and accelerate the online search stage. Firstly, all 2-combinations of points are extracted and their 4-parameter vectors are computed. Then, a 3-dimensional R-tree is created using all points and , to account for the switched point-wise correspondences. Each object of the tree also includes the value for and two indices that identify the pair of points in the point cloud.
Our proposed online search scheme (Fig. 7(b)) starts by extracting a random pair of points from the curve, and its tangents, and computing its descriptor . This pair is then used for querying the R-tree for selecting all pairs in the surface that simultaneously have a distance , where is a parameter to account for noise in the data, and satisfy the conditions in Equation 13. The obtained set of pairs is afterwards pruned by choosing only the ones that satisfy Equation 14. The obtained correspondences of pairs of points are then processed in a RANSAC scheme in order to find the rigid transformation that yields the highest number of inliers. If all the correspondences have been processed and the stopping criteria was not met, the algorithm repeats this process for a new random pair of points extracted from the curve.
4 Extensions to curve vs curve and surface vs surface
This section shows how to solve the global 3D registration problem for two 3D models of the same type: two curves or two surfaces. We will provide an explanation for the curve vs curve alignment, with the derivations being identical for the case of surface vs surface registration.
Consider two curves C and Ĉ and two corresponding 2-tuples point+tangent with descriptors and . The methods and derivations of Section 2 hold with two differences. First, the constraints for determining the angle of rotation (Equation 3) become
meaning that can be computed in closed form using two points and just one tangent. This is valid because in this case is the rotation that aligns the corresponding tangent vectors.
where the sign accounts for the fact that the tangents are in general non-oriented. Instead of being inequalities, as in the curve vs surface alignment, in this case the conditions for matching are equalities, enabling search mechanisms other than R-trees. This is validated in the experimental section, where the search is carried by extracting the pairs of points that satisfy the conditions in Equation 16 using the pair extraction scheme proposed in [Mellado14] that runs in time, with being the number of points in the target curve. Note that this scheme does not contain an offline processing stage as the search algorithm proposed in Section 4.
As stated before, the surface vs surface registration is similar to curve vs curve, having the difference that tangents are replaced by normals. We will not discuss this problem further because the resulting method is equivalent to the one that has been recently presented and tested in [Raposo17].
This section reports tests performed on synthetic and real data in order to assess the accuracy and speed of the proposed registration methods. The first experiments use synthetic data for which the ground truth rigid transformations are known, and compare our curve vs surface and curve vs curve registration methods with two state-of-the-art approaches for which there is public implementation available: Super4PCS [Mellado14] and Fast Global Registration (FGR) [FGR]. The last experiment attempts to mimic a common CAOS procedure, where 3D data on the surface of a bone is reconstructed and registered with a pre-operative virtual dense model of that bone.
The normals were computed using the PlanePCA algorithm [planePCA] with a neighbourhood of 30 points and the tangents were estimated using a standard algorithm for computing Frenet frames. Both registration algorithms were implemented in C++ and all tests were performed on a Intel Core i5-6200U CPU @ 2.30GHz with 8GB of RAM.
5.1 Curve vs surface registration using synthetic data
In this experiment we used the 4 synthetic models shown in Fig. 8 for evaluating the performance of the proposed curve vs surface algorithm. The sets of 6 segments shown in red in the figure represent the curves that were manually extracted from each model in order to create the curves for performing the registration.
The extracted sets of curves were used for generating smaller sets of curves with the intent of assessing the performance of the algorithms for different amounts of data.
This was done by randomly choosing 2 and 4 out of the 6 segments and then selecting a random set of contiguous points in each segment such that the total number of points is about and of the total number of points of the original sets, respectively.
This scheme was used for creating 25 different sets of each of the two sizes, to which random rigid transformations are applied. We also consider the full curve (containing of the points) in 25 different initial poses generated randomly. This procedure yields 75 different input curves for each model, to which random noise drawn from the standard normal distribution with mean 0 and standard deviation
of the points) in 25 different initial poses generated randomly. This procedure yields 75 different input curves for each model, to which random noise drawn from the standard normal distribution with mean 0 and standard deviationis added. We test both with noise-free data () and by adding noise with which represents of the diameter of the models, whose dimensions were previously adjusted by setting their diameters to . The noise was added to each point independently. This level of noise () causes the direction of tangents to vary with respect to by an average of . Also, we defined as stopping criteria for the search algorithms a maximum execution time of 5 sec or a percentage of inliers of 95.
Fig. 9 shows the results obtained with our proposed curve vs surface method and Super4PCS. We also tested with the FGR method but it failed in all cases because it is a feature-based approach and we are working with two different types of 3D data. Results are given as rotation and translation errors, computed as in [Raposo17], and computational times. The results are merged for all models such that each boxplot corresponds to 100 registrations. Regarding our method, the reported times only include the online search stage and information on the offline processing of the models can be found in Table 1.
The superiority of our approach w.r.t. Super4PCS both in terms of accuracy and speed is evident as it was able to provide proper alignments for all the different conditions of size of input data and noise. On the other hand, Super4PCS performed poorly for the sets of curves with and of the points of the original data, being only able to provide acceptable solutions for the complete sets. As expected, the accuracy of our method also increases with the amount of input data as a larger coverage of the surface is given. However, even for the smallest curve sizes, it provides good alignments, indicating that the method is able to work with very local information. Concerning computational times, our method is very fast, being able to perform the search in less than 1 sec in all cases. Even if we consider the time corresponding to the offline processing of the surface, which in the worst case is 1.9 sec, the total execution time of our method is still far below Super4PCS’s.
5.2 Curve vs curve registration using synthetic data
This experiment is performed similarly to the previous one, having the difference that we perform curve vs curve alignment by replacing the surfaces with the curves represented in Fig. 8. We tested with FGR but do not show the obtained results since it was only able to provide acceptable solutions in the noise-free case. When noise was added, the number of corresponding features became very low and the method performed poorly for all curve sizes.
Results in Fig. 10 show the superiority of our approach w.r.t. Super4PCS in the noise-free case. When noise is added, our method still outperforms Super4PCS when the amount of input data is very small, performing similarly when of the data is provided. This can be explained by the fact that tangents are more affected by noise than points, leading to a degradation in the performance of our approach. However, it still manages to accomplish a proper alignment of the curves in all situations in under 1 sec.
5.3 Experiments in CAOS using a dry knee model
The last experiment mimics a common procedure in CAOS where 3D data is reconstructed by touching the surface of the bone with an instrumented touch probe and subsequently registered with a pre-operative 3D virtual model of that bone. In order to simulate this, we used a dry knee model to work as the bone, as shown in Fig. 11(a), and acquired 3D curves on the surface of the bone using a state-of-the-art optical tracking system (the Optotrak Certus). We acquired 30 curves in the condyliar region, highlighted in yellow in Fig. 11(a), where 15 of them contained of outliers. An example of a curve with outliers is given in Fig. 11(b). Each curve was used for performing curve vs surface registration using our proposed approach and the obtained rigid transformations are used to represent the 23 control points illustrated in Fig. 11(a) in the virtual model reference frame. These control points had been previously reconstructed by carefully placing the touch probe in the small holes of the model (Fig. 11(c)) and their ground truth 3D coordinates in the virtual model reference frame are known. Fig. 11(d) shows the distributions of distances between the transformed and the ground truth points for three different regions of the knee, both for the curves without forced outliers (solid color) and with outliers (transparent).
As expected, the obtained distances are smaller in the region where the data was acquired (condyles). However, in the other regions the errors do not increase substantially, not even in the femoral shaft that is more than 10 cm away from the area of acquisition. This is an important result since it confirms the previous observation that our method is able to properly register large surfaces with very local information, being advantageous in CAOS procedures where the area of the bone that is exposed is often restricted. Another relevant observation is that our method is able to deal with large amounts of outliers, shown by the fact that there was not a significant degradation of the registration accuracy when using data with outliers. This demonstrates that besides being accurate and fast, the proposed method is robust and resilient to outliers, which highly improves its usability.
We present the first method for fast global registration of curves and surfaces that does not require an initial coarse alignment. The method makes use of pairs of points, augmented with their local differential information, not only to solve the rigid transformation estimation problem but also to establish correspondences of pairs of points in a very fast manner. Experiments demonstrate that the proposed method significantly advances the state-of-the-art by providing a fast and robust 3D registration algorithm that dramatically outperforms two plausible alternatives for global registration.
As future work, we intend to extend the rigid transformation estimation algorithm to a more general method for determining not only the rotation and translation but also the scale. This has applications in CAOS has it would, for instance, allow problems of difference in size between the virtual models and the respective bones to be overcome.
The authors thank the Portuguese Science Foundation and COMPETE2020 program for generous funding through project VisArthro (ref.: PTDC/EEIAUT/3024/2014). This paper was also funded by the European Union’s Horizon 2020 research and innovation programme under grant agreement No 766850.