Non-iterative rigid 2D/3D point-set registration using semidefinite programming

7.4 MB·English

by  Yuehaw Khoo

1
Non-iterative rigid 2D/3D point-set registration
using semidefinite programming
Yuehaw Khoo∗† and Ankur Kapoor∗
∗Imaging and Computer Vision, Siemens Research Coorporation, Princeton, NJ 08540, USA
†Department of Physics, Princeton University, Princeton, NJ 08540 USA
Abstract—We describe a convex programming framework for 2D/3D point-set registration problem described in Section III
6 pose estimation in 2D/3D point-set registration with unknown with guarantees using our semidefinite programs conreg1
1 point correspondences. We give two mixed-integer nonlinear and conreg2. We further demonstrate their usefulness in the
program (MINLP) formulations of the 2D/3D registration prob-
0 application of coronary vessel imaging.
lem when there are multiple 2D images, and propose convex
2
relaxations for both of the MINLPs to semidefinite programs From a broader perspective, 2D/3D point-set registration

r (SDP) that can be solved efficiently by interior point methods. problem arises in numerous medical imaging applications in
p
Ourapproachtothe2D/3Dregistrationproblemisnon-iterative fields such as neurology [1], orthopaedics [2], and cardiol-
A
in nature as we jointly solve for pose and correspondence.
ogy[3].Theassociatedbodyofliteratureon2D/3Dregistration
Furthermore, these convex programs can readily incorporate
6 is expanding rapidly, as is apparent in the thorough review of
feature descriptors of points to enhance registration results. We
prove that the convex programs exactly recover the solution to techniques recently published by Markelj et al. [4]. Amongst
]
V the MINLPs under certain noiseless condition. We apply these thevastliteratureonthe2D/3Dregistrationproblem,webriefly
formulationstotheregistrationof3Dmodelsofcoronaryvessels discussthetechniquesmostcloselyrelatedtotheonepresented
C
to their 2D projections obtained from multiple intra-operative
in this paper, namely, methods for registration of point-sets.
. fluoroscopic images. For this application, we experimentally
s As mentioned earlier, a typical point-set registration prob-
c corroborate the exact recovery property in the absence of noise
[ and further demonstrate robustness of the convex programs in lem consists of two mutually interlocked subproblems, pose
the presence of noise. estimation and point correspondence. The key idea in seminal

3 Iterative Closest Point (ICP) [5] algorithm is to alternatively
v Index Terms—Rigid registration, 2D/3D registration, iterative-
0 closest-point, convex relaxation, semidefinite programming. solve the two subproblems starting from an initial estimate of
3 pose. For the correspondence subproblem, it uses closeness in
6 terms of Euclidean distance between two points to determine
0 I. INTRODUCTION whether they correspond to each other. There are many
0 RIGID registration of two point sets is a classical problem variants proposed to enhance the original ICP, either to
.
1 in computer vision and medical imaging of finding a make fast correspondence [6] or to include more features to
0 transformation that aligns these two sets. Typically, a point-set obtain high quality correspondence [7]. Variants for 2D/3D
5 registration problem consists of two intertwined subproblems, registration include [8]–[11]. Though popular, these methods
1 pose estimation and point correspondence, where solving one suffer for common drawbacks, namely, they are all local
:
v is often the pre-condition for solving the other. A canonical methods and do not guarantee global convergence. Their
Xi formulation of the rigid point-set registration problem for two performances all rely on a good initialization and the spatial
point clouds is the following. Let X,Y ∈Rd×m be two point configuration (distribution) of 3D points. In 2D/3D point-set
r
a sets in dimension d, we want to solve: registration, for this type of local optimizers several strategies
min (cid:107)RX−YP(cid:107)2, (1) have been proposed to increase the capture range by the use
F
P∈Πm×m,R∈SO(d) of multi-resolution pyramids [12], use of stochastic global
d
optimization methods such as differential evolution [9] or
where Πm×m is the set of m×m permutation matrices, and
d direct search wherein the parameter space is systematically
R∈SO(d) is the rotation matrix from the special orthogonal
and deterministically searched by dividing it into smaller and
group in d dimension. Finding a solution to this problem is
smaller hyperrectangles [9]. In these cases, except in direct
difficult as it is a nonconvex, mixed integer nonlinear problem.
search, the guarantees on finding the correct global minima
Anothercloserelativeofthisproblem,namelythe2D/3Dpoint-
are very weak. In the case of direct search one requires the
set registration, assumes a 3D point-set and 2D projections of
parameter space to be finite.
the 3D point-set. The objective is to find out the pose of the
Another line of work [13]–[15] in point-set registration
3D model that gives rise to the 2D degenerate point-set upon
focus on using soft or inexact correspondences to enhance
projections. This adds an additional complication to that of
the search for global optimizers. In a recent work [16], the
the regular point-set registration problem, namely, the loss of
Oriented Gaussian Mixture Model (OGMM) method is proposed
depth information. In this paper, we focus on solving forms of
to extend the successful Gaussian mixture based nonrigid
registrationalgorithm(gmmreg)[15]tothe2D/3Dcase.These
Y. Khoo (email: [email protected]), A. Kapoor (email:
[email protected]). methods model the point configuration as Gaussian mixtures,
2
and they intend to find a transformation that maximizes the convex programs conreg2 gives solution to a variant of
overlap between these distributions. The structure of Gaussian 2D/3D registration problem where we simultaneously estimate
distribution encapsulates the idea of soft correspondence and the point correspondences between multiple 2D images while
enables fast implementation. Empirically, as more fuzziness respecting the epipolar constraint. This could be utilized in
is allowed for point correspondences (larger variance for the finding the corresponding image points for 3D reconstruction
Gaussians), the target function of optimization is smoother and purposes [18]. Another natural extension of these programs is
hence it is less likely to find a local optimum. Nevertheless, a the incorporation of local descriptors of points as additional
good initialization is still crucial for these types of algorithms terms in the objective of both the programs. For our clinical
due to the nonconvex nature of the cost. application, we use tangency of the vessel in the local neigh-
Other than local optimization approaches to the registration borhood of the point to illustrate the use of point descriptors.
problem, efforts have been put forth to exhaustively search Exact recovery analysis:Weproveexactnessof theconvex
for the global optimum using methods such as branch and programs,thatis,undercertainconditionstheproposedrelaxed
bound [17]. While such techniques can be fast in practice, in convex programs will give a solution to the original MINLPs.
principleregistrationproblemsthatrequireexponentialrunning We prove that under noiseless situation the relaxed convex
time can exist. In this paper, we introduce a different approach programs are in fact able to exactly recover the rotation
to jointly and globally solve the pose and correspondence and permutation matrices that match the projected points to
problem in 2D/3D point-set registration by relaxing the hard the 3D model. Our simulations show that algorithms’ global
registrationproblemtoconvexprograms.Unlikelocalmethods, convergence results also hold with noise when the error in
the solutions of convex programs do not depend on the recovery grows nearly proportionally to the added noise. Real-
initialization. We prove for certain noiseless situations the data examples corroborate the theoretical results, and suggest
global optimizer of the surrogate convex problem coincides potential applications in coronary tree matching.
with the global optimizer of the original problem. Therefore Here we outline the organization of this paper. In Section II,
through convex programming we are guaranteed to attain the we summarize the notation used in the paper. In Section III,
global optimum in polynomial time. In the noisy regime where we describe the type of 2D/3D registration problems we intend
the solution of the convex programs closely resembles the to solve and present two MINLPs formulations when there
optimizer of the original problem, the proposed method will are multiple images. In Section IV, we present the convexly
converge to the neighborhood of the original optimizer. In relaxed versions of the 2D/3D registration problem in terms of
these situations, it is practically advantageous to use convex tractable semidefinite programs. In Section V, we prove that
optimization instead of greedy methods or global optimization achieving global optimality is possible under certain situations.
that may run in exponential time. From a theoretical point In Section VI-A, we mention how to incorporate additional
of view, the convex relaxation framework can be used to features to the convex programs to enhance registration results,
characterize the polynomial time solvable cases in 2D/3D in particular in the application of coronary vessel imaging.
registration with unknown correspondences, through analyzing Lastly, in Section VII, we empirically verify the exact recovery
the exact recovery condition mathematically. property, and demonstrate the robustness of the algorithm on
Before proceeding to the rest of the sections, we note simulated and real medical datasets for the registration of 2D
that the correspondence estimation sub-problem, which is coronary angiogram with the 3D model of the coronaries.
typically treated as the quadratic assignment problem, is NP-
hard on its own. Therefore the solution to our proposed convex
programs may not always be close to the solution of the II. NOTATION
joint pose and correspondence estimation problem. There are
We use upper case letters such as A to denote matrices, and
certainly cases for which the use greedy approaches is the
lower case letters such as t for vectors. We use I to denote
only resolution, though not globally optimal. However, we d
the identity matrix of size d×d. We denote the diagonal of a
hopethepromisingresultsdemonstratedbytheproposednovel
matrix A by diag(A). We use An for integer n≥1 to denote
non-iterative approach towards 2D/3D registration will attract
the multiplication of A with itself n times. We will frequently
proposals with tighter convex relaxations. The shortcomings
use block matrices built from smaller matrices, in particular
of our approach are discussed in the concluding remarks.
when we deal with problem (REG2). For some block matrix
A. Our contributions A, we will use Aij to denote its (i,j)-th block, and A(p,q)
to denote its (p,q)-th entry. We also use A to denote the i-th
i
The main contributions of this paper can be summarized as
column of A. When we have an index set s = {s ,...,s }
1 n
the following:
where each s is an integer, we use A to denote the matrix
Algorithm: The original 2D/3D point-set registration of a i s
[A ,...,A ]. We use A (cid:23) 0 to mean that A is positive
3D model and multiple 2D images is formulated as nonconvex s1 sn
semidefinite. We use (cid:107)x(cid:107) to denote the Euclidean norm of
2
MINLPs that are difficult to solve. We propose two convex
x∈Rn. We denote the trace of a square matrix A by Tr(A).
relaxations for these MINLPs. The programs jointly optimize
We use the following matrix norms. The Frobenius, mixed
the correspondences and transformation, and the convex
(cid:96) /(cid:96) and entry-wise (cid:96) norms are defined as:
2 1 1
nature of these programs enable efficient search of global
optimum regardless of initialization using standard off-the- (cid:107)A(cid:107) =Tr(ATA)1/2, (cid:107)A(cid:107) =(cid:88)(cid:107)A (cid:107) , and
F 2,1 i 2
shelf conic programming software. Furthermore, one of the
i
3
(cid:107)A(cid:107) =(cid:88)(cid:88)|A(i,j)|. (2) of points between I and I be denoted by Pˆ , where
1 a b ab
i j Pˆab(i,j)=1 if point i in Ia corresponds to point j in Ib and
zero otherwise. We do not require one-to-one correspondence
The Kronecker product between matrices A and B is denoted
between points such as each row or column of Pˆ can
by A⊗B. The all-ones vector is denoted by 1. We use |s| ab
have more than one nonzero entries. In the presence of such
to denote the number of elements in a set s. For a set s we
correspondence information, R can be obtained as the solution
use conv(s) to denote the convex hull of s. We introduce the
to the following optimization problem:
following sets,
N N N
Πa×b ≡{A∈{0,1}a×b :(cid:88)a A(i,j)=1,(cid:88)b A(i,j)=1}, (REG2) R∈mSOin(3), (cid:88)(cid:107)ΨiRX−IiPi(cid:107)2,1+(cid:88)(cid:88)(cid:107)Pab−Pˆab(cid:107)1
d (Pi,Pab)∈Si=1 a=1b>a
i=1 j=1
a b where Pab is the permutation matrix that relates the points
(cid:88) (cid:88)
Πla×b ≡{A∈{0,1}a×b : A(i,j)=1, A(i,j)≤1}, of Ia and Ib for N(N −1)/2 image pairs, and (cid:107)·(cid:107)1 is the
i=1 j=1 entry-wise (cid:96)1 norm. The domain of optimization for the S will
a b be explained in the next subsection. We note that solving this
(cid:88) (cid:88)
Πa×b ≡{A∈{0,1}a×b : A(i,j)≤1, A(i,j)≤1}, problem could also be useful if we are interested in estimating
i=1 j=1 the correspondence between the N(N − 1)/2 image pairs
and we frequently call them the permutation, left permutation directly for reconstruction purpose.
and sub-permutation matrices. Fortheremainderofthispaperweonlyconsidertwoimages
to simplify notation. It should be obvious the solution we
propose in the subsequent sections can be easily extended
III. PROBLEMSTATEMENT
to the multiple images case. Problems (REG1) and (REG2)
A 3D centerline representation of a coronary artery tree, have nonconvex domains which consist of integer and rotation
segmentedfromapreoperativeCTAvolume,istoberegistered matrices, as thus are very difficult to solve. In subsequent
to N fluoroscopic images. Our observed 3D model with m sections, we formally define these domains and present convex
points is described by matrix X ∈ R3×m. The projection relaxation that can recover the exact solution under certain
operatorsbetweenthe3Dspaceandthei-thimage,represented conditions.
by Ψ , are known from the calibration of the X-ray apparatus.
i
The projection operator maps the 3D model to a degenerate A. Domains for permutation matrices
point cloud Ii ∈R3×ni that represents i-th projection image. A more formal way to understand the problems (REG1),
By degenerate we mean the affine rank of I is two. We will
i (REG2) is the following: Suppose there is a ground truth 3D
assume that ni ≥ m, 1 ≤ i ≤ N. The 2D/3D registration model with m(cid:48) points which is described by the coordinate
problem is to find an alignment matrix R, in the special matrix X ∈R3×m(cid:48). Our observed 3D model is described by
orthogonalgroupSO(3),thatmatchessomepermutationofthe matrix Xgt∈ R3×m. Assume m ≤ n ,n ≤ m(cid:48). In this case,
1 2
observed 3D model with each of the degenerate projections.
we have
We propose to solve the multiple images 2D/3D registration
problem as: I1 =Ψ1RXgtQ1, I2 =Ψ2RXgtQ2, (3)
X =X Q , (4)
N gt 3
(cid:88)
(REG1) min (cid:107)Ψ RX−I P (cid:107) ,
i i i 2,1 where R, Ψ ,Ψ ,I ,I are matrices as introduced as before,
R∈SO(3), 1 2 1 2
Pi∈Πnii×mi=1 Q1 ∈ Πml (cid:48)×n1,Q2 ∈ Πml (cid:48)×n2 and Q3 ∈ Πml (cid:48)×m. One can
intuitively regard Q ,Q ,Q as operators that generate the
where (cid:107)·(cid:107) is the mixed (cid:96) /(cid:96) norm. We do not know a 1 2 3
2,1 2 1
imagesandtheobservedmodelbysub-samplingandpermuting
priori the correspondences between points in the 3D model
the points in ground truth 3D model X . For example, if the
and projection image (which are encoded in P for each of N gt
i ith row of Q is zero, then it means the i-th point of X (or
given images), and the rotation R. These are found by solving 1 gt
more precisely, i-th column of X ) is not selected to be in
(REG1). Intuitively, the minimization of the cost in (REG1) gt
I .
simply ensures by subsampling and permuting I , the image 1
i Herewemakeanassumption,thatifapointiinX (thei-th
points should as close as possible to the projections of the 3D gt
columnofX )iscontainedinX,thentheprojectionsofpoint
model X posed by some rotation R. We employ (cid:107)·(cid:107) norm gt
2,1 imustcorrespondtosomecolumnsofI ,I .Looselyspeaking,
inordertoalleviatethecostsduetotheoutliers.Ifweconsider 1 2
it means the observed 3D model X is a subset of the point
a maximum likelihood estimation framework with Gaussian
clouds I ,I . In this case we know P = QTQ ∈ Πn1×m
type noise, a squared Frobenius norm (cid:107)·(cid:107)2 could be replaced 1 2 1 1 3 l
F and P = QTQ ∈ Πn2×m. Furthermore, if a point i of
instead. 2 2 3 l (cid:80)
X is not selected by Q (meaning Q (i,j) = 0), then
Next we consider a variant of problem (REG1) where the gt 1 j 1
Q QT ∈ Rn×n is almost an identity matrix I , except it is
correspondences of the points between the two images are 1 1 n
zero for the ith diagonal entry (similarly for Q QT). We use
available. Such correspondences could come from feature 2 2
these facts after multiplying the equations in (3) from the right
matching, or by exploiting epipolar constraints. This basically
by QTQ ,QTQ to get
means P need to be optimized dependently to preserve 1 3 2 3
i
correspondences between the images. Let the correspondence Ψ RX =I P , Ψ RX =I P . (5)
1 1 1 2 2 2
4
When there is no noise in the image, the equations in (5) have lay out the details of the fully relaxed problem in the next
tobesatisfied.Ifnotweturntosolvetheoptimizationproblem two subsections. We present the relaxations in the case of two
(REG1). images,andthegeneralizationtomorethantwoimagesshould
TosimultaneouslyestimatethecorrespondenceP ,weusea be obvious to the reader.
12
similar construction as [19] and [20]. Let P =QTQ . When
12 1 2
having exact correspondence matrix Pˆ without ambiguity,
12
we require
P12 =Pˆ12. (6) A. Relaxing (REG1)
We relate P ,P ,P through a new variable
1 2 12
Weremindreadersthatinproblem(REG1)wefindarotation
P ≡[Q1 Q2 Q3]T[Q1 Q2 Q3]∈R(n1+n2+m)×(n1+n2+m). to register the 3D model with the images. The points in the 3D
(7) model X and I ,I should be matched after a permutation
1 2
In this context, the variables P1,P2 used in (REG1), along operationthroughthevariablesP1,P2.P1,P2 havetheirentries
with P12 are simply the off diagonal blocks of the variable P. onthedomain{0,1}.Itisdifficulttooptimizeonsuchdiscrete
Then the domain of (P1,P2,P12) is simply domain. A typical way to bypass such difficulty is by the
following discrete-to-continuous relaxation,
S = {(P ,P ,P ):
1 2 12
P (cid:23)0,rank(P)≤n,P =I ,
ii ni
P1 ∈Πnl1×m,P2 ∈Πnl2×m,P12 ∈Πn1×n2}. (8) Pi ∈{0,1}ni×m → Pi ∈[0,1]ni×m (9)
P should be understood as the diagonal blocks of P.
ii
Again, we solve the optimization (REG2) when images or fori=1,2.TherelaxedmatricesP1 andP2 canbeinterpreted
correspondences are noisy. This formulation ensures that the as probabilistic map between the points in the image and
left permutation and sub-permutation matrices P ,P ,P are the points in the 3D model. For example, the size of an
1 2 12
cycle-consistent. This basically means the relative permutation entry P1(i,j) resembles the probability of point i in I1 being
matricesP ,P ,P allcanbeconstructedbysomeunderlying corresponded to point j in the 3D model. Such relaxation has
1 2 12
global “sampling” matrices Q ,Q ,Q for each image and the been used for problems such as graph isomorphism [25].
1 2 3
3D model (see [20] for a detail description). Wenowintroducetheconvexlyrelaxedproblemconreg1,
IV. CONVEXRELAXATION
(conreg1)
In this section, we propose two programs, which are the
min (cid:107)Ψ RX−I P (cid:107) +(cid:107)Ψ RX−I P (cid:107)
convexrelaxationsofproblem(REG1)and(REG2).Wereplace R,P1,P2 1 1 1 2,1 2 2 2 2,1
the non-convex constraints that prevents these problems from s.t R∈conv(SO(3)),
being convex to looser convex constraints. These relaxations P ∈[0,1]ni×m,
i
enables the efficient search of global optimum using standard
P 1≤1,1TP =1,for i∈{1,2}.
convex programming packages. To solve the convex problems i i
we use CVX, a package based on interior point methods [21]
for specifying and solving convex programs [22], [23]. Before
The last constraint simply follows from the definition of
looking at each of the two problems, we introduce a result
left permutation matrix in the notation section. In case when
from [24] to deal with the nonconvexity of SO(3) manifold. there are N >2 images, there should be N matrix variables
In [24] the authors give a characterization of the convex hull P ,...,P for i∈{1,··· ,N} with constraints as indicated
of SO(d) for any dimension d, denoted as conv(SO(d)), in 1 N
above.
terms of positive semidefinite matrix. The d=3 version will
be restated in the following theorem. Due to space limitation,
for this theorem we locally define the (i,j)-th entry of a matrix
X as X instead of X(i,j).
ij B. Relaxing (REG2)
Theorem 1 (Thm 1.3 of [24]). conv(SO(3))={X ∈R3×3 :
In problem (REG2) we include the estimation of P ,
12
(cid:34) (cid:35) the correspondence between image I ,I , on top of pose
1−XX111−3+XX223+1X331+XX111−3+XX223−1X33 XX1223−−XX2312 XX2132++XX3221 estimation. In this case, we relax the SO1(3)1and the integrality
X12−X21 X23−X32 1+X11+X22+X33 X31−X13
X23+X32 X12+X21 X31−X13 1−X11+X22−X33 constraints just like the case of problem (REG1). However,
is positive semidefinite }.
different from (REG1), the domain of optimization for (REG2)
includes the set S in (8), which has an extra rank requirement
Such linear matrix inequality is convex and can be included that prevents the convexity of (REG2). A popular way of
in a semidefinite program. For both problems (REG1) and dealing with the rank constraint is by simply dropping it. Such
(REG2), we relax the domain of R, which is SO(3) into relaxation is common in many works, for example in the
conv(SO(3)). However, this relaxation alone is not sufficient seminal paper of low rank matrix completion [26]. We now
to turn (REG1) and (REG2) into convex problems, and we state the convex problem conreg2
5
simply by finding the closest orthogonal matrix to R˜ in terms
of Frobenius norm, namely
(conreg2)
R(cid:63) =argmin(cid:107)R−R˜(cid:107)2 (12)
min (cid:107)Ψ RX−I P (cid:107) +(cid:107)Ψ RX−I P (cid:107) F
R,P1,P2,P12 1 1 1 2,1 2 2 2 2,1 R∈O(3)
+(cid:107)P12−Pˆ12(cid:107)1 where O(3) is the orthogonal group in 3D. For completeness
s.t R∈conv(SO(3)), we provide the analytical expression for R(cid:63). Let the singular
P ∈[0,1]ni×m, for i∈{1,2}, value decomposition (SVD) of R˜ be of the form
i
P12 ∈[0,1]n1×n2, R˜ =UΣVT. (13)
P 1≤1, 1TP =1, for i∈{1,2},
i i Then
P121≤1, 1TP12 ≤1, R(cid:63) =UVT. (14)
 
I P P
n1 12 1 The reader maybe curious that the special orthogonality
PT I P (cid:23)0.
 12 n2 2 constraint is not enforced in the rounding procedure. This
P1T P2T Im is indeed unnecessary since an element R˜ in conv(SO(3)) has
positivedeterminant[24],henceUVT obtainedfromtheabove
procedure necessarily is in SO(3).
Again, conreg2 inherits the equality and inequality con-
As we see in the case when the solutions of conreg1,
straints for P ’s and P from the the property of left permu-
i 12 conreg2arenotfeasiblefor(REG1)and(REG2),arounding
tation and sub-permutation matrices. In case when there are
procedure is required to project R(cid:63) to a nearest point in
N >2images,thereshouldbeN matrixvariablesP ,...,P
1 N SO(3). Although it is possible to prove (cid:107)R(cid:63) − R (cid:107) being
and N(N −1)/2 variables P . The relevant constraint that 0
ij small (as in other applications of convex relaxation, see [27]
encodes cycle consistency in such case should be
for example), R(cid:63) in general is still suboptimal in terms of the
 
I P ··· P P costoftheoriginalnonconvexproblem(REG1)and(REG2),as
n1 12 1N 1
P1T2 In2 ··· P2N P2 the rounding procedure does not consider optimizing any cost.
 ... ... ... ... ... (cid:23)0. (10) It is thus quite typical to use the approximate solution from
  convex relaxation as an initialization and search locally for a
P1TN P2TN ··· InN PN more optimal solution, for example in [28]. In this paper, we
PT PT ··· PT I
1 2 N m use the OGMM algorithm, which searches locally by optimizing
a nonconvex cost using gradient-descent based method.
C. Connection to classical registration problem in 3D
It is obvious that similar relaxation can be derived for the V. EXACTRECOVERY
classical rigid registration problem of (1). For the sake of Inthissection,wewillgiveproofthatthetworelaxedconvex
completeness, we state the convex formulation to the 3D programs exactly recover the solution to the original problems
registration problem in (1). (REG1) and (REG2) under certain situations. We will show
that the solutions returned by our algorithm conreg1 and
min (cid:107)RX−YP(cid:107)2
R,P F conreg2 lie in the domain of (REG1) and (REG2), hence
s.t R∈conv(SO(3)), showing the attainment of global optimizers of (REG1) and
P ∈[0,1]m×m,P1=1,1TP =1. (11) (REG2), as the domain of (REG1) and (REG2) are properly
contained in the domain of conreg1 and conreg2. By this
Again, in applications where robustness is a requirement, a we mean exact recovery. The proofs of exact recovery for
mixed (cid:96) /(cid:96) can always be replaced instead. The discussion conreg1 and conreg2 are almost the same. We will go
2 1
of this program is out of the scope of this paper and will be through proving exact recovery for conreg1 in detail, and
included in a future publication of the authors. just state the results for exact recovery of conreg2.
D. Projection and local optimization A. Our assumptions
Although in section V we prove that under certain situations We summarize the assumptions for our proof.
the two convex programs give solution to (REG1) and (REG2), (1) We consider the situation when n1 = n2 = m = m(cid:48),
this will not always be the case as the domain of search has which means that all the image points I1,I2 come from the
been made less restrictive. As we are concerned about posing projection of the 3D model.
the 3D model in this paper, we describe a strategy in this (2) We consider the noiseless situation when the equations
section to obtain an approximate solution from conreg1 and I P¯ =Ψ R¯X, I P¯ =Ψ R¯X, (15)
conreg2fortherotationingeneralcases.Denotetheoptimal 1 1 1 2 2 2
“rotation” in the domain conv(SO(3)) of both programs as can be satisfied, where P¯ ,P¯ ∈Πm×m and R¯ ∈SO(3). We
1 2 d
R˜. As R˜ is not necessarily in SO(3), we need a strategy to will call these variables the ground truth. Notice that under
project R˜ to a point R(cid:63) on SO(3). The strategy we use is assumption (1), P¯ ,P¯ are permutation matrices instead of
1 2
6
left-permutation matrices. In the case of conreg2, we further are in one-to-one correspondence with solutions to conreg1
assume with 3D model R¯X and image I P¯ ,I P¯ . Therefore without
1 1 2 2
Pˆ =P¯ . (16) loss of generality, using the change of variable introduced in
12 12
(17), we can let the 3D model being X(cid:48) = R¯X, the images
In order to show exact recovery, we want to show that the being I(cid:48) = I P¯ = Ψ R¯X = Ψ X(cid:48) and I(cid:48) = Ψ X(cid:48). We
solution R˜,P˜ ,P˜ to conreg1, is in fact R¯,P¯ ,P¯ . For 1 1 1 1 1 2 2
1 2 1 2 thus consider solving the problem conreg1 in the following
conreg2 we in addition show P˜ =P¯ .
12 12 form
(3) Without lost of generality we can consider R¯ =I and
3
P¯1,P¯2 =Im. Such consideration is backed by lemma 1. min (cid:107)Ψ1RX(cid:48)−Ψ1X(cid:48)P1(cid:107)2,1+(cid:107)Ψ2RX(cid:48)−Ψ2X(cid:48)P2(cid:107)2,1
(4) The m points are in generic positions. A set S = R,P1,P2
s.t R∈conv(SO(3)),
{s ,s ,...,s } of n real numbers are called generic if there
1 2 n
is no polynomial h(x1,...,xn) with integer coefficients such Pi ∈[0,1]ni×m,
that h(s1,...,sn) = 0 [29]. When we say the points are in Pi1≤1,1TPi =1,for i∈{1,2}.
genericpositions,wemeanthesetof3mrealnumbersfromX
(21)
is generic. This assumption excludes the possibilities that three
points can lie on a line or four points can lie on a plane. It In this case, instead of showing R˜ =R¯,P˜ =P¯ ,P˜ =P¯ ,
1 1 2 2
also implies the centroid of any subset of the points is nonzero. it suffices to show R˜ =I and P˜ ,P˜ =I being the unique
3 1 2 m
Furthermore, no two subsets can share the same centroid. solution to conreg1.
(5) There are at least four points, i.e. m≥4.
Lemma 2.
(6) Without lost of generality, we consider Ψ [0 0 1]T =
1
Ψ2[0 0 1]T = [0 0 1]T (the directions of projection are Ψ1R˜X(cid:48) = Ψ1X(cid:48)P˜1 (22)
perpendicular to z-axis). Ψ R˜X(cid:48) = Ψ X(cid:48)P˜ (23)
2 2 2
We note that although we only give proof of exact recovery
under these assumptions, in practice it is possible for our and
algorithm to recover the pose exactly when n ,n ≥m.
1 2 Ψ R˜nX(cid:48) = Ψ X(cid:48)P˜n (24)
1 1 1
Ψ R˜nX(cid:48) = Ψ X(cid:48)P˜n (25)
B. Exact recovery property of conreg1 2 2 2
Lettheoutcomeoftheoptimizationconreg1beR˜,P˜ ,P˜ . for any n≥1.
1 2
Theorem 2. Under our assumptions, the solution R˜,P˜ ,P˜ to Proof. See appendix A
1 2
conreg1 is unique, and R˜ =R¯,P˜1 =P¯1,P˜2 =P¯2. Lemma3. P˜ ,P˜ aredoublystochasticmatricesanditsuffices
1 2
When we have exact recovery, R(cid:63) obtained from the to consider (cid:107)P˜1−Im(cid:107)1 <1 and (cid:107)P˜2−Im(cid:107)1 <1. Then
projection of R˜ to SO(3) is simply R˜, as R˜ is already in lim P˜n = lim P˜n =A, (26)
SO(3). To prove the theorem, we first state a couple lemmas. n→∞ 1 n→∞ 2
where A is partitioned by index sets a ,...,a and
Lemma 1. If we make the following change of data 1 k
X →X(cid:48), I1 →I1(cid:48), I2 →I2(cid:48), (17) A(ai,ai)=(1/|ai|)11T, A(ai,aj)=0 for i(cid:54)=j. (27)
where Proof. See appendix B
X(cid:48) =R¯X, I(cid:48) =I P¯ , I(cid:48) =I P¯ , (18)
1 1 1 2 2 2 We are now ready to prove theorem 2.
then the solutions of conreg1 are changed according to the Proof. We first want to show R˜ = I . Taking limits of (24),
3
invertible maps
(25) and using lemma 3, we get
R˜ →R˜R¯T, P˜1 →P¯1TP˜1, P˜2 →P¯2TP˜2. (19) Ψ1BX(cid:48) =Ψ1X(cid:48)A, Ψ2BX(cid:48) =Ψ2X(cid:48)A, (28)
Proof. whereB =lim R˜n andA=lim P˜n =lim P˜n
n→∞ n→∞ 1 n→∞ 2
(cid:107)I P˜ −Ψ R˜X(cid:107) +(cid:107)I P˜ −Ψ R˜X(cid:107) . We saw in lemma 3 A has multiple irreducible components,
1 1 1 2,1 2 2 2 2,1
= (cid:107)I P¯ P¯TP˜ −Ψ R˜R¯TR¯X(cid:107) + each being an averaging operator for the points of relevant
1 1 1 1 1 2,1
(cid:107)I P¯ P¯TP˜ −Ψ R˜R¯TR¯X(cid:107) indices. The equations lead to
2 2 2 2 2 2,1
= (cid:107)I1(cid:48)P¯1TP˜1−Ψ1R˜R¯TX(cid:48)(cid:107)2,1+ BX(cid:48) =X(cid:48)A. (29)
(cid:107)I(cid:48)P¯TP˜ −Ψ R˜R¯TX(cid:48)(cid:107) . (20)
2 2 2 2 2,1
Since B =lim R˜n =R˜lim R˜n−1 =R˜B, multiply-
The last second equality is due to the fact that P¯1P¯1T = ing (29) from tnh→e∞left by R˜ we gne→t ∞
I ,P¯ P¯T = I ,R¯TR¯ = I , since permutation matrices are
m 2 2 m 3
necessarily orthogonal matrices. Therefore, a solution after the R˜X(cid:48)A=X(cid:48)A (30)
change of data in (17) is P¯1TP˜1, P¯2TP˜2 and R˜R¯T. Now we show R˜ =I for each of the following cases.
3
Since the map presented in (19) is invertible, we know the Case 1: Suppose A has three or more than three irreducible
solutions to conreg1 with 3D model X and image I ,I components. Let the centroid of any three components be
1 2
7
[c1,c2,c3]∈R3×3.Byequation(30)R˜[c1,c2,c3]=[c1,c2,c3]. VI. CLINICALAPPLICATIONANDADDITIONALFEATURES
By the assumption that columns of X(cid:48) are generic, [c ,c ,c ]
1 2 3
has full rank hence invertible. Thus R˜ =I .
3 Thoughtheproblemof2D/3Dregistrationappearsinseveral
Case 2: Suppose A has two irreducible components. This
medical imaging applications, this paper is motivated by the
means R˜[c ,c ]=[c ,c ] where c and c are independent by
1 2 1 2 1 2 clinicalsettinginwhichtwoormorefluoroscopicimagesofthe
the assumption that points are generic. Let cˆ and cˆ denote
1 2 coronaries are used for intervention guidance. In a procedure
the normalized version of c ,c . R˜ =Ψ +αcˆ ∧cˆ , where
1 2 c1c2 1 2 called percutaneous coronary intervention (PCI), commonly
Ψ denotes the projection operator to the plane containing
c1c2 referred to as angioplasty, a cardiac interventionalist introduces
c ,c , ∧ denotes the vector cross product and α some scalar
1 2 a thin, flexible tube with a device on the end via arterial
constant with |α| ≤ 1. We know |α| ≤ 1 since the largest
access (usually femoral) into the patient’s coronary arteries
singular value of R˜ is less than 1, as
which may be blocked by calcification. To guide and position
conv(SO(3))⊂conv(O(3))={O ∈R3×3|OTO (cid:22)I } the device accurately, the cardiac interventionalist uses X-ray
3
fluoroscopy with intermittent contrast injections (angiography).
(31)
whereconv(O(3))istheconvexhulloftheorthogonalgroupin The resultant fluoroscopic images contain projected outlines of
thecoronaryvesselsandblockagesifany.However,theseintra-
3D. Such characterization of the convex hull of the orthogonal
operative fluoroscopic images do not have the performance
group can be found in [30]. If |α|=1, then we indeed have
proven R˜ =I . If |α|<1, we will arrive at a contradiction. In characteristic of pre-operative images such as Computed
3
this situation B = lim R˜n = Ψ . By the assumption Tomography Angiography (CTA) with blood pool contrast
n→∞ c1c2 injection that enables the occlusions causing calcifications to
that m ≥ 4 and equation (29), there exists four points with
coordinates Y =[Y ,Y ,Y ,Y ]∈R3×4 such that either be clearly distinguished. The benefit of using pre-operative
1 2 3 4
Coronary CTA images are well documented [32]–[35] and
Ψ Y =[c ,c ,c ,c ], or Ψ Y =[c ,c ,c ,c ]. (32) are becoming standard of care. Coronary CTA may also
c1c2 1 1 1 2 c1c2 1 1 2 2
allowbetterdeterminationofcalcification,lesionlength,stump
ThismeanseitherY1,Y2,Y3 formaline,orY1,Y2,Y3,Y4 form morphology, definition of post-CABG anatomy, and presence
a plane and both cases violate the assumption that points are of side branches compared to fluoroscopy or angiography
in generic positions. alone [33].
Case 3: Suppose A has a single irreducible component. In Toaidthecardiacinterventionalistperformingtheprocedure
this case BX(cid:48) = [c1,...,c1], where c1 is the center of the on the coronaries, a strategy to provide additional information
point cloud. It is easy to show that B = Ψc1, where Ψc1 is byfusingpre-operativeimageswithintra-operativefluoroscopic
the projection onto the line spanned by c1. However, Ψc1X(cid:48) = has been developed. However, finding a transformation to
[c1,...,c1]impliescolumnsofX(cid:48) formsaplanewithc1 being align pre-operative and intra-operative images is challenging
thenormalvector.Againthisviolatestheassumptionofgeneric since finding inter modal correspondences is a nontrivial
positions. problem. Typically to achieve this multi-modality fusion,
Now since R˜ =I3, from equations (22) and (23) we have vessel centerlines segmented from both modalities are used
as landmarks for inter modal correspondences. A 3D model
[0 0 1]X(cid:48) =[0 0 1]X(cid:48)P˜ , [0 0 1]X(cid:48) =[0 0 1]X(cid:48)P˜ . (33)
1 2 of the vessels extracted from pre-procedure CTA is aligned
Again from the generic positions assumption, [0 0 1]X(cid:48) ∈Rm with the 2D fluoroscopic views, and overlaid on top of the
live fluoroscopy images, thereby augmenting the interventional
isavectorwithdistinctentries.Itisknownthatinthissituation,
the two doubly stochastic matrices P˜ ,P˜ =I (see corollary images.
1 2 m
8 of [31]). In this paper, we process the CTA images using the
segmentation algorithm presented in [36] and fluoroscopic
images using [37], [38] to get the respective point-sets. Thus,
C. Exact recovery property of conreg2 the clinical alignment transformation between pre-operative
Let the solution of conreg2 be R˜, P˜ , P˜ , P˜ . and intra-operative images is reduced to (REG1) or (REG2).
1 2 12
Though this paper was motivated by the clinical scenario
Theorem3. Underourassumptions,thesolutionR˜,P˜ ,P˜ ,P˜
1 2 2 of CTO, our algorithms are generalized to handle other
toconreg2isunique,andR˜ =R¯,P˜ =P¯ ,P˜ =P¯ ,P˜ =
1 1 2 2 12 clinical situations where natural structures can be obtained by
P¯ .
12 image segmentation methods. For instance, vessels from brain,
Proof. The recovery of R¯,P¯ ,P can be proven in the exact liver, kidney or peripheral images can be used by processing
1 2
mannerasinprevioussubsection.IfPˆ =P¯ ,(cid:107)P˜ −Pˆ (cid:107)> them with relevant segmentation algorithms. As illustrated in
12 12 12 12
0 unless P˜ =P¯ . We therefore conclude our proof. Section ??, the algorithm is not restricted to point-sets from
12 12
vessel like structures only. Other structures, notably surfaces,
Although it does not seem like the additional constraints in extracted by means of segmentation algorithm can be used
conreg2 helps our proof at this point, in practice in the case to obtain desired point sets required by methods (REG1)
when n ,n ≥ n , conreg2 demonstrates slight advantage or (REG2). Typically, the structure representations extracted
1 2 3
over conreg1 in simulation in terms of stability (Fig. 4b) from segmentation methods need to be re-sampled uniformly
and also in real data (Fig. 5). to result in respective point-sets.
8
A. Incorporating additional features
Inthecontextoffindinganalignmenttransformationbetween
two point clouds, point descriptors maybe valuable by favoring
transformation that matches the descriptors as well as the
coordinates. Our proposed formulations easily allow the incor-
poration of descriptors by adding additional terms in the cost
of (REG1), (REG2) that encourage matching of transformation
invariantfeatures,orfeaturestransformingaccordingtorotation.
To match transformation invariant features, we simply add the
Fig.1. Anillustrationofsyntheticdatausedforsimulations.Solidbluelines
following terms: representthecenterlineson2Danddashedredlinesrepresenttheprojection
of3Dmodelontheimageplane.Theleftfigureillustratesthecasewhenthe
(cid:107)XF −IFP (cid:107) +(cid:107)XF −IFP (cid:107) , (34) 3Dmodelisapropersubset.
1 1 2,1 2 2 2,1
toourcostin(REG1)and(REG2).XF ∈Rd×m,IF ∈Rd×n1,
and I2F ∈Rd×n2 denote some d dimensional feat1ure vectors m) 100 m) 100
m m
forthepointsin3Dmodelandtheimages,respectively.Todeal ( (
with features that transform with rotation, such as intensity MSD 10−3 MSD 10−3
R R
gradient, the following terms can be added to (REG1) and
(REG2): 10−6 10−6
[0,10][10,20][20,40][40,60][60,90] [0,10][10,20][20,40][40,60][60,90]
(cid:107)Ψ RXF −IFP (cid:107) +(cid:107)Ψ RXF −IFP (cid:107) . (35) InitialPose(deg) InitialPose(deg)
1 1 1 2,1 2 2 2 2,1
(a) (b)
??conreg1 ??conreg2 ??OGMM
B. Constructing features for coronary vessel Fig.2. SimulationresultsofthetwoconvexprogramsandOGMM(a)when
theimagesaretheprojectionofthecomplete3Dmodel.and(b)whenthe
imagescontaintheprojectionsofthe3Dmodelasapropersubset.Herewe
To incorporate features for the specific application of
onlyplotthecaseswhenOGMMfindsasolutionwithin25mmRMSD.Iffor
registering coronary vessel, we will leverage the idea of patch aparticularrotationrangethereisnoboxplotitmeansnosolutionhasless
[39]. We note that the typical features for vessels, such as the than25mmRMSDforthatrange.
tangent of the vessel, are in general extracted from the local
neighborhoodoftheparticularpoint[40].Therefore,wesimply
VII. EXPERIMENTS
combine the coordinates of neighborhood points around each
point, which we denote as patch, as the descriptor for each In this section we evaluate the performance of our algorithm
point. For the coordinate of point i, we concatenate it with the using synthetic data of the heart vessels, and real CT data of
coordinateofsubsequentpointsi+1,i+2,...,i+n −1toform six different patients. In particular, we demonstrate through
p
a patch of size n . We denote the result of such concatenation simulation that regardless of the initial pose of the 3D model,
p
as XP, IP and IP , and we store it as the i-th column of the our algorithm exactly recovers the ground truth pose when
i 1 i 2 i
matrices XP ∈ R3np×m and I1 ∈ R3np×n1, I2 ∈ R3np×n2. thereisnonoise.Inthepresenceofimagenoise,ouralgorithm
Now instead of registering point, our goal is then to register is able to return a solution close to the ground truth. These
these patches of the 3D model to the images. The cost in observations are also verified when working with real data. For
(REG1) when considering matching patches become these experiments we solve conreg1 and conreg2 using
the convex programming package CVX in MATLAB using
(cid:107)(I ⊗Ψ R)XP−IPP (cid:107) +(cid:107)(I ⊗Ψ R)XP−IPP (cid:107) , interiorpointmethod.Ingeneralthetimecomplexityofinterior-
np 1 1 1 2,1 np 2 2 2 2,1
(36) point-method for conic programming is O(N3.5) where N
and for (REG2) the cost is denotes the total number of variables [41]. The main source of
complexity in our program comes from the stochastic matrices
(cid:107)(Inp ⊗Ψ1R)XP −I1PP1(cid:107)2,1+ P1 and P2. Assuming there are m points, the time complexity
(cid:107)(I ⊗Ψ R)XP −IPP (cid:107) +(cid:107)P −Pˆ (cid:107) . (37) is then O(N3.5) ∼ O(m7). Since the interior-point-method
np 2 2 2 2,1 12 12 1
is a second order method requiring the storage of Hessian,
In our experiments, we found n = 3 suffices to improve the space complexity is O(N2) ∼ O(m4). While this looks
p
the solution without introducing too many variables into the intimidating, as reported in the upcoming sections it is quite
optimization. A practical issue is that one needs to take care fast in practice due to the fast implementation of the solvers
of the sampling density of the points on the centerline. If the in CVX. To measure the quality of registration, we use the
point spacing is very different between X and I , I , the root-mean-square-distance
1 2
patches cannot be well matched. To deal with this issue, we 1
RMSD= √ ((cid:107)Ψ R(cid:63)X−Ψ R¯X(cid:107) +(cid:107)Ψ R(cid:63)X−Ψ R¯X(cid:107) ).
first project the randomly posed 3D model using Ψ1 and Ψ2, 2 m 1 1 F 2 2 F
calculatetheaveragespacingbetweentwosubsequentpointsin (38)
the projection, and sample the image points I , I according As a reminder, R(cid:63) is the solution after rounding from
1 2
to this spacing. conreg1,conreg2, and R¯ is the ground truth rotation. R¯
9
2 0.025 10 10
conreg1 m) m)
0.02 conreg2 (m (m
1 MSD0.015 MSD 5 MSD 5
R 0.01 R R
0.005
0 0
−01 0 1 00 0.005 0.ε01 0.015 0.02 ∞ 21.6 15.6 12.1 9.6 ∞ 21.6 15.6 12.1 9.6
Noise(dB) Noise(dB)
(a)Bunnyε=0 (b)Stability (a) (b)
conreg1 conreg2
Fig. 3. (a) A simulation using the Stanford bunny dataset when there is
no noise. The red points are taken from I1, and the blue points are taken Fig.4. Simulationresultsofthetwoconvexprograms(a)whentheimages
fromthemodelX posedbythesolutionfromconreg1.Theposeisexactly
aretheprojectionofthecomplete3Dmodeland(b)whentheimagescontain
recoveredwhenthereisnonoise.(b)RMSDv.s.noiselevelforthebunny
theprojectionsofthe3Dmodelasapropersubset.Foreachnoiselevel,we
dataset.Eachdatapointisaveragedover20noiseandposerealizations.The
do50experimentswithdifferentposesofthe3Dmodel,withrotationranging
averagespacingbetweenthepointsisabout0.04. from0◦ to90◦.TheRMSDgraduallyincreasesasthenoiseincreasesfrom
nonoise(∞dB)to2mm(9.6dB).
is known in the case of synthetic data, and for the case of
real data it is manually given by medical experts. The unit of
RMSD will be in millimeter.
there are extra centerlines in the 2D images. In the latter case,
the 3D model is a proper subset. In our simulations, we first
A. Synthetic data
rotatethe3Dmodelofcoronaryvesselsintosomearbitrarypos-
We use synthetic data to demonstrate the ability of our tureandprojectitintwofixedorthogonaldirections.Wetested
algorithm to exactly recover the pose of the 3D model when our algorithm systematically when rotation R in (39) is within
there is no noise in the image. We also add bounded noise to each of the ranges [0◦,10◦],[10◦,20◦],[20◦,40◦],[40◦,60◦],
each point in the image in the following way and [60◦,90◦], using 10 different noise realizations in each
case.
I =Ψ (RX +ε ), I =Ψ (RX +ε ),
1i 1 i 1i 2j 2 j 2j
i∈1,...,n , j ∈1,...,n , (39) Weobservedexactrecoveryinbothprogramswhenthenoise
1 2
level is zero, which confirms our proof in Section V. Since our
where ε , ε ∼ U[−ε,ε]3 and U[−ε,ε]3 is the uniform
1i 2j proposed methods are convex programs, the solution, in this
distribution in the cube [−ε,ε]3.
casethegroundtruthrotation,doesnotdependoninitialization.
In order to run conreg2, we simply get Pˆ from the
12 Fig. 2 shows the results of both programs with no noise added
epipolarconstraints.Forapointinimage1,denotedasx∈R3
to the points with respect to the initial pose of the 3D model.
in homogeneous pixel coordinate, the epipolar line l in image
In our experiments, we set the tolerance level for the convex
2 can be computed as l=Fx, where F ∈R3×3 denotes the solverto10−3.Fig.4showstheresultsforbothprogramswith
fundamental matrix. If point x(cid:48) in image 2 satisfies |lTx(cid:48)|=
bounded noise added to each point in the image. In this case
|xTFx(cid:48)|≤δ where δ is some pre-defined threshold, we then
the RMSD increases proportionally as the noise ε increases
set the entry in Pˆ corresponding to the points x and x(cid:48) to be
12 from 0 to 2 mm for both input conditions. We observe that the
1. We pick δ to be 2 to 3 times the average spacing between
latter input condition of extra centerlines in the image may be
the neighboring points.
a better approximation of a typical clinical scenario as often
We first start with the Stanford bunny dataset to illustrate
non vascular objects such as catheters are mis-segmented as
the exact recovery of pose from conreg1 and conreg2.
vessels.
We first project the bunny in x and y-directions to obtain I
1
and I2, respectively. We then pose the model bunny X with Wecompareouralgorithmwiththerecent2D/3Dregistration
arbitrary rotation and register it to its projections. In Fig. 3a, algorithm OGMM [16], [42]. In our experiments, we use our
we show the registration results of conreg1 in one viewing MATLAB implementation of [16] for two images with each
direction with m = 473. When there is no noise, the exact individualgaussiandistributionbeingisotropic.Ifwearetouse
recovery of pose is verified in Fig. 3. While we do not give OGMM, using the identity as an initialization, the recovery of
guarantees on the approximate recovery of pose when there the pose is not guaranteed unless we are in the [0◦,10◦] range.
is noise, we also show in Fig. 3b that our proposed methods In fact for most experiment we cannot find solution within 25
are stable when perturbed by noise. We note that for Fig. mm RMSD. We show the results from OGMM in Fig. 2, where
3b, we sparsify the number of points to m = 60 and the we plot the RMSD against different 3D poses. We note that
average distance between neighboring points is 0.04. Although there are many ways one can improve on our implementation
conreg1 can solve registration problem with m = 473 in of OGMM in terms of radius of convergence, such as annealing
less than 10 seconds, the computational cost of conreg2 is to zero the variance of Gaussian distributions starting from a
high and we have to limit the number of points. large value as proposed in [15]. Nevertheless, it still remains
For the registration of coronary vessels, we perform simu- that registration methods based on Gaussian mixtures method
lations for two different cases (Fig. 1), namely, when the 2D [15] can get stuck in a local optimum due to the nature of the
images are exactly the projections of the 3D model and when cost.
10
60
MSD(mm) 1200 2400 MSD(mm) 1200 40
R R 20
0 1 2 3 4 5 6 0 1 2 3 4 5 6 0 [0,10][10,20][20,40][40,60][60,90] 0 [0,10][10,20][20,40][40,60][60,90]
Patient Patient
InitialPose(deg) InitialPose(deg)
(a) (b)
(a) (b)
Fig.5. (a)Resultsusingthetwoconvexprogramsconreg1(Red),conreg2
(Green),andOGMM(Blue)fordataof6patients.(b)Numberofinstancesfor Fig.6. (a)Resultsusingthetwoconvexprogramsconreg1(Red),conreg2
whichtheRMSDerrorwaslessthan10mm(solidcolor),15mm(northeast (Green),andOGMM(Blue)fordataof6patientswithrespecttoangleofpose.
lines)and25mm(northwestlines),respectively.Foreachofthe6patients, (b)NumberofinstancesforwhichtheRMSDerrorwaslessthan10mm(solid
foreachrotationrangewegenerate10differentinstancesresultinginatotal color),15mm(northeastlines)and25mm(northwestlines),respectively
of50instancesperpatient. with respect to angle of pose. For each of the 6 patients, for each rotation
rangewegenerate10differentinstancesresultinginatotalof50instances
perpatient.
B. Real data
To illustrate the application of our algorithm on clinical be resolved adequately using conreg1 without enforcing
data, we have used pre-processed images from Chronic Total explicit matching between images. In such cases, the point
Occlusion (CTO) cases treated using PCI. Both our methods consistency terms in conreg2 may aid in resolution of such
(REG1)and(REG2)operateonpoint-setsandhenceanyim- ambiguities. We present the images of the 3D model after the
ages must be pre-processed to segment relevant structures such 2D/3D registration in Fig. 7.
asvesselsororgansurfaces.IncaseofCTOcases,therelevant In Fig. 7 we not only perform registration using conreg1
structures to align are the centerlines of the coronaries as seen and conreg2 but also further refine their solutions using
in pre-operative CTA images and intra-operative fluoroscopic gmmreg. Typically for convex relaxation techniques such as
images.ToextractthecenterlinesinCTAimages,wehaveused the one adopted in this paper, it is not unusual to prove the
the segmentation algorithm presented in [36]. The resultant convex relaxed solutions to the noisy data lie in the vicinity of
centerlines also provide topological information allowing us to the ground truth [27], [43], [44]. We, in fact, demonstrate
computegeometricfeaturessuchastangents.Thesecenterlines via simulations with synthetic data that the solutions of
were re-sampled at uniform interval to result in the 3D point- algorithmsconreg1andcongreg2liewithincertainradius
set for registration methods. The fluoroscopic images were (proportional to the noise magnitude) of the ground truth (Fig.
processed using 2D vessel segmentation algorithms [37], [38] 4).Thus,localoptimizerssuchasICPorGMMbasedmethods
toproviderespectivecenterlinesintheprojectionimages.These after conreg1 and conreg2 could be used to further refine
were re-sampled to provide the 2D point-sets. The projection the results. As local search based registration methods often
matrix was estimated from the C-arm parameters stored in the fail to reach the global optimum due to non-convex nature of
DICOM header of the image. the cost, this combined approach ensures the ground truth is
We test our algorithm on six sets of clinical data in which within the basin of convergence of any local methods. It is
a medical expert has aligned the 3D model to 2D fluoroscopic clear that this strategy brings improvement in row 1,3 and 4
images. We consider this the ground truth to compare against in Fig. 7 especially when registering using conreg1. While
our results. To characterize the recovery of our algorithm, adding a local refinement step brings little changes visually
we perturb the aligned 3D model by an arbitrary rotation in row 2 and 5 and even seems counterproductive for row
within the ranges [0◦,10◦],[10◦,20◦],[20◦,40◦],[40◦,60◦], 6 in Fig. 7, the RMSD when comparing with the ground
and[60◦,90◦].Foreachrotationrange,wegenerate10different truth given by the clinician is in fact lowered. For example,
instances. The results are detailed in figure 5 and 6. Again it is shown in Fig. 5a that conreg1 and conreg2 fail to
for comparison we use results from OGMM. We see that we achieve 15 mm RMSD for images of patient 6. However after
are able to recover the pose to within 10 mm for most cases. a local refinement all test cases for 6 patients achieve 15 mm
Further, the recovery of pose in consistent irrespective of the RMSD (Table. I). The deterioration of registration quality for
initial pose. Whereas, in general we cannot obtain RMSD patient 6 in image plane A stems from the fact that the local
within 25 mm with OGMM. On average, there are 80 points refinement tries to fit the model to the vessels in both of the
in each model and image. The average running time for image planes. It may bring improvement to registration in one
conreg1 and conreg2 are 0.4s and 15s, respectively on an of the image plane while sacrificing the registration quality for
Intel(cid:13)R CORETM i7 2.2GHz running MATLAB(cid:13)R. conreg2 theotherimageplaneinordertolowertheoverallcost.Results
is significantly slower than conreg1 due to the additional combining algorithms conreg1 and congreg2 with OGMM
positive semidefinite matrix that encodes cycle consistency. and a variant of LM-ICP [45] as local search methods are
The effect of additional point consistency terms between the presented in Table I. We include two different local methods
images in conreg2 is more prominent in real patient data. to illustrate that the final results are not dependent on the
We note that clinical data typically may have more number of local refinement method utilized. Though the desired accuracy
ambiguous 3D model to projection point matches that cannot of registration is highly dependent on application with some