Proceedings of the 2005 IEEE
International Conference on Robotics and Automation
Barcelona, Spain, April 2005

Fast 2D-3D Registration for Navigation System ofSurgical Robot
Yumi Iwashita, Ryo Kurazume
and Tsutomu Hasegawa

Kozo Konishi, Masahiko Nakamoto
and Makoto Hashizume

Graduate School of Information Science and Electrical Engineering
Kyushu University
6-10-1, Hakozaki, Higashi-ku, Fukuoka, Japanyumi@irvs.is.kyushu-u.ac.jp

Graduate School of Medical Sciences
Kyushu University
3-1-1, Maidashi, Higashi-ku, Fukuoka, Japan

Abstract— This paper presents a new registration algorithm of
2D color images and 3D geometric models for navigation systemof surgical robot. A 2D-3D registration procedure is used tosuperimpose a tumor model on an endoscopic image precisely,
and is therefore indispensable for the surgical navigation system.
And thus, the performance of the 2D-3D registration procedure
influences directly the usability of the surgical robot operating
system. One of the typical techniques developed so far is the use of
external markers, but the accuracy of this method is getting worse
due to the breath, the heart beat, or other unknown factors. For
precise registration of 3D models and 2D images without externalmarkers or special measurement devices, the new registration
method is proposed which utilizes the 2D images and theirdistance maps created by the Fast Marching Method. We show
some results of fundamental experiments using simulated models
and actual images of the endoscopic operation.

I. I NTRODUCTION
Development work on a practical surgical robot system isproceeding aiming at a minimally invasive surgical operation
and a favourable postoperative process. Fig.1 shows the surgical robot system named Intuitive Surgical’s “da Vinci”.

(a) surgeon's consoleFig. 1.

(b) surgical cart

The surgical robot system “da Vinci”

In general, the surgical robot system is composed of two
main components, a manipulation system and a navigation
system.
The manipulation system of the “da Vinci” is so-called
the master-slave manipulator system. The slave manipulator
consists of 3 arms; two 7 degrees-of-freedom manipulators andone 5 degrees-of-freedom camera holder. The manipulators are

0-7803-8914-X/05/$20.00 ©2005 IEEE.

909

used to hold and operate many types of forceps, which canbe replaced according to an operation purpose. The camera
holder is used for controlling the position of an endoscope,
which is equipped with 2 CCD cameras. Captured stereo
images of an internal organ are displayed on two LCDs in asurgeon’s control console, and give stereognostic perception tothe operator. The operator is able to control the forceps and theendoscope quite intuitively using the stereognostic perception.
Our navigation system for the “da Vinci” is designed toassist an operator by presenting a proper operation procedure
to be taken and displaying patient’s information related to theoperation [1]. Another main function of the navigation system
is to display an invisible internal structure of an organ, locationof a tumor, or distribution of critical blood vessels on theendoscope image with high reality. The operator can obtain
useful information from these realistic images intuitively sothat the operation is performed smoothly.
The virtual images of the organ or the structure of blood vessels are created from 3D geometric models, which are obtainedusing medical imaging devices such as CT (Computerizing
Tomography) scanner or MRI (Magnetic Resonance Imaging)
beforehand. The obtained 3D geometric models are projectedon the image plane of the endoscope and superimposed on 2D
images of the target organs. The performance of this procedure
influences directly the usability of the surgical robot operating
system. Thus, the registration of 2D images and 3D geometricmodels has to be performed as precisely as possible.
Several registration techniques of 2D images and 3D geometric models have been proposed so far. One of the typicaltechniques is the use of external markers to calibrate the coordinate systems fixed on the patient’s body and the endoscope.
In this method, the position of several markers, pasted on theskin of the trunk and the endoscope respectively, are measured
by an external optical device. Then, the calibration calculationis performed to determine the relative pose between these twocoordinate systems. After the calibration, the 3D geometricmodels of organs are projected on the 2D images taken fromthe endoscope based on the standard position information of
the organ in the patient’s body. Therefore, the accuracy of thismethod is getting worse if the organ is away from the estimatedposition in the body due to the breath, the heart beat, or otherunknown factors.
This paper presents a new registration algorithm of 2D color
images and 3D geometric models without external markers
or special measurement devices. This method utilizes the 2D
images and their distance maps created by the Fast Marching
Method, and is able to determine a precise relative pose
between 2D images and 3D models. Since the registration
process can be executed at a quite high speed using actual
images, the registration accuracy is expected to be higher thanthe one of the conventional method using the external markers,
even if the position of the organ is disturbed by the breath,the heart beat, or other factors.
After an overview of previous works in Section II, we
will give a short introduction of the Fast Marching Method
in Section III, which is utilized for the fast constructionof a distance map. A detail description of the new 2D-3D
registration algorithm is presented in Section IV. In Section
V, we show some results of fundamental experiments using
actual images of the endoscopic operation.
II. R ELATED WORKS
For the alignment of a 2D color image and a 3D geometricmodel, Viola proposed a technique based on a statisticalmethod [2]. This method evaluated the mutual information
between the 2D image and the image of the 3D model basedon the distribution of intensity. Allen et al. [3] also proposed amethod using the intersection lines of several planes extracted
from the range data and the color edges. These methods work
well on the surface with little albedo variance. But theseare easily trapped in a local minimum on surfaces with richtextures.
Range sensors often provide reflectance images as sideproducts of range images. This 2D reflectance image is alignedwith the 3D range image because both images are obtained
through the same receiving optical device. Therefore, we
can utilize the 2D reflectance image for solving the 2D3D registration problem instead of the original 3D model.
By using the 2D reflectance image, the 2D-3D problem ismodified to the simple 2D-2D problem. Kurazume et al.
[4] proposed the registration algorithm using the reflectanceimage and the 2D texture image. In his method, a number ofphotometric edges extracted from both images were registeredand the relative pose is determined using a robust M-estimator.
Elstrom et al. [5] and Umeda et al. [6] also proposed themethods using the reflectance image and the texture image.
However, an organ doesn’t have a rich texture on its surface
and there are also differences of textures between individuals.
Therefore, the method using texture is quite difficult for the
2D-3D registration problem of the surgical navigation system.
On the other hand, some registration techniques using asilhouette image or a contour line have been proposed. Lensch[7],[8] proposed a silhouette-based approach. He used XOR
value of silhouette images of a 2D image and a 3D modelfor evaluating a registration error. The Downhill Simplex
Method was used for the convergent calculation of the 2D-3D
registration. In contour-based approach, the error is computed910

as the sum of distances between points on a contour line in a
2D image and on a projected contour line of a 3D model [9],[10]. Since a number of point correspondences between bothcontours has to be determined for calculating the registration
error, these algorithms are computationally expensive.
Our 2D-3D registration algorithm belongs to the contourbased approach. However, we adopt a distance map for evaluating a registration error instead of the point correspondence.
Therefore, once the distance map is created, our algorithm runs
faster than the conventional point-based approach[9], [10]. In
addition, since the distance map can be created quite rapidly
using the Level Set Method named the Fast Marching Method,
our method is able to track an object continuously even if theobject moves.
III. FAST M ARCHING M ETHOD FOR RAPIDCONSTRUCTION OF DISTANCE MAP

Before describing our registration algorithm, the Fast
Marching Method proposed by Sethian [11],[12] is briefly
introduced, which is utilized for the fast construction of a
distance map in our method.
The Fast Marching Method [11],[12] was initially proposed as a fast numerical solution of the Eikonal equations(|∇T (p)| F = 1), where T (p) is arrival time of a front
at a point p and F is a speed function. In general, since
this equation is solved using a convergent calculation, it
takes a long time to obtain a proper solution. However, byadding restriction such that a sign of the speed function F
is invariant, the Fast Marching Method solves the Eikonal
equation straightforwardly, and thus, quite rapidly. In thismethod, the arrival time of a front at each point is determinedin the order from the old to the new.
At first, the Fast Marching Method modifies the Eikonal
equation into the next difference equation.
−x+x(max(Dij
T, −DijT, 0)2 +1−y+yT, −DijT, 0)2 ) 2 = 1/Fij
max(Dij(1)

Next, since the arrival time is propagated in one direction
from the old to the new, a point which holds the oldest arrival
time in a whole region is chosen, and the arrival time of theboundary at this point is determined from Eq.(1).
The concrete procedure is as follows:
Step.1 (Initialization) The whole region is divided into a
number of grid points with a proper grid width.
Before starting the calculation, all the grid points are
categorized into three categories (known, trial, far)according to the following procedure.
1) The grid points belonging to the initial front(denote the boundary hereafter) are added tothe category of known and the arrival time ofthese grid points is set to 0 (T = 0).
2) Among the 4 neighboring grid points of a grid
point belonging to the category of known, theone which doesn’t belong to the category ofknown is categorized into to the category oftrial, and the arrival time of this grid pointis calculated temporarily from the equation,
Tij = 1/Fij . In addition, the arrival time ofthese grid points is stored in heap data structureusing the heap sort algorithm in ascending orderof T .
3) The grid points except for the above grid points
are all categorized into the category of far, andthe arrival time is set to infinity (T = ∞).
Step.2 Choose the grid point (imin , jmin ) which is placed at
the top of the heap structure. This grid point has thesmallest arrival time in the category of trial. Remove
this grid point from the category of trial and the heap
structure, and categorize it as the category of known.
Run “downheap” algorithm to reconstruct the heap
structure.
Step.3 Among the 4 neighbouring grid points((imin − 1, jmin ),(imin + 1, jmin ), (imin , jmin −1),(imin , jmin + 1)) of the selected grid point(imin , jmin ), the one which belongs to the categoryof far is changed to the category of trial.
Step.4 Among the 4 neighbouring grid points of the selected
grid point (imin , jmin ), the arrival time of the one
which belongs to the category of trial is calculatedusing Eq.(1). Then, run “upheap” algorithm to reconstruct the heap area.
Step.5 If there is a grid point which belongs to the categoryof trial, then go to Step 2. Otherwise, the process is
terminated.
This Fast Marching Method constructs a distance map of
a whole region quite rapidly, which indicates the distancefrom a boundary (an initial position of the front) to a certain
point. To construct the distance map, the speed function Fij
in Eq.(1) is set to 1 at first. Next, the arrival time of the front
T is determined using the above procedure. Since the speedfunction is 1, T suggests the distance from the boundary to
the point. Fig.2 shows an example of the calculated distancemap.

IV. A N EW 2D-3D REGISTRATION ALGORITHM BASED ON
THE DISTANCE MAP

This section describes our fast 2D-3D registration algorithm
using the distance map in detail. We assume that a 3D
geometric model of the target, such as a structure of an organ
or a tumor, has been constructed by CT or MRI images
beforehand, and is represented by a number of triangular
patches. A 2D color image is also taken from a camera
attached to an endoscope.
The proposed method is summarized as follows:
1) At first, a boundary of a target on the 2D color image is
detected using an active contour model (ex. Snakes or
the Level Set Method[13]).
2) Next, the distance map from the detected boundaryon the 2D image plane is constructed using the Fast
Marching Method.
3) The 3D geometric model of the target is placed at anarbitrary position and projected on the 2D image plane.
4) Contour points of the projected image and their corresponding patches on the 3D model are identified.
5) Force is applied to the selected patch of the 3D model in
3D space according to the distance value obtained fromthe distance map at each contour point.
6) Total force and moment around the COG (center ofgravity) of the 3D model is determined using the robust
M-estimator.
7) The pose of the 3D model is changed according to thetotal force and moment.
8) Repeat from step 1 to 7 until the projected image of the
3D model and the 2D image coincide each other.
The above procedure is explained in more details with someexamples in the following sections.
A. Construction of distance map
Fig.3 shows the calculated distance map. At first, a boundary is extracted from the color image using the Level Set
Method[13]. Next, the distance map from this boundary is
calculated using the Fast Marching Method explained in
Section III.
boudaryboundary

contour lines of
distance map(a) Object(b) Obtained distance map+20 pixels

Fig. 3.

+100 pixels

Detected boundary and distance map

B. Fast detection of triangular patches of contour points in
3D model+20 pixels

Fig. 2. An example of the distance map calculation using the Fast Marching
Method911

Fig.4 shows an example of the contour detection of the
projected 3D geometric model. The contour detection andidentifying triangular patches on the 3D model correspondingto points on the contour line are computationally expensive andtime consuming. In our implementation, we utilize the highspeed rendering function of the OpenGL hardware accelerator
and execute these procedures quite rapidly.
The detailed algorithm is as follows: Initially, we assign
different colors for all the triangular patches in the 3D modeland draw the projected image of the 3D model on the image
buffer using the OpenGL hardware accelerator. The contour
points of the 3D model are detected by raster scanning of
the image buffer. By reading colors of the detected contour
points, we can identify the corresponding triangular patcheson the 3D geometric model. Fig.5 shows the color image ofthe projected 3D model with the patches of different colors.

where ri is a vector from the COG to the triangular patch iand ρ(z) is a particular estimate function.
In practice, a part of the target organ is occasionally hiddenby the other organs in the endoscopic image. In this case,the obtained boundary does not coincide with the projected
contour of the 3D model and the correct distance map cannot
be obtained. Therefore, the robust M-estimator is used in oursystem to ignore contour points with a large amount of errors.
Let’s consider the force fi and the moment ri × fi as an error
zi . Then, we modify Eqs.(2) and (3) as Fρ(zi )(4)
E(P ) ==
Miwhere P is the pose of the 3D geometric model.
The pose P which minimizes the E(P ) is obtained as thefollowing equation.

detected contour ∂ρ(zi ) ∂zi
∂E
=
=0
∂P
∂zi ∂Pi

Fig. 4.

(5)Here, we define the weight function w(z) as the following
equation in order to evaluate the error term.

Contour detection of 3D geometric modelw(z) =

1 ∂ρ
z ∂z

(6)

From the above equation, we obtain the following weighted
least squares method.

∂E
∂zi
=
=0
w(zi )zi∂P∂Pi(a) 3D modelFig. 5.

In our implementation, we adopt the Lorentzian function
for the estimation function ρ(z), and gradually minimize theerror E in Eq.(4) by the steepest descent method.

(b) Object image with different
color patches

Detection of triangular patches of the contour points

ρ(z) =
C. Force and moment calculation using the M-estimator anddetermination of relative pose
After obtaining the distance map on the 2D image and thelist of the triangular patches of the 3D model corresponding to
the contour points, the force fi is applied to all the triangular
patches of the contour points (Fig. 6). The force fi is the vector
perpendicular to the line of sight, and the projection of the fi
onto the 2D image plane coincides with the fDMi , which isthe vector toward the direction of the steepest descent of the
distance map(Fig. 7(a)). We assume that the magnitude of the
fDMi is proportional to the value on the distance map at the
projected contour point of the 3D model. Then, total force andmoment around the COG are calculated using the following
equations as shown in Fig.7.
F =ρ(fi )
(2)
i

M

=



ρ(ri × fi )

(7)

(3)

i

912

w(z)

=

σ2
log(1 + (z/σ)2 )2
1
1 + (z/σ)2

(8)
(9)

The P which minimizes the error in Eq.(4) is the estimatedrelative pose between the 2D image and the 3D geometricmodel.

fi

Fig. 6.
points.

Applying the force f to all the triangular patches of the contourfi
f DM ifv

fi

Objectr
M
F

Projected
contourImage planef DM

Steepest descentdirection

Projected
contour pointv
Focul point

Boundary(a)Fig. 7.

(b)Force and moment around COG

D. Coarse-to-fine strategy and the “distance band”
Though the 2D distance map can be constructed efficiently
by the Fast Marching Method, computing the distance value
in whole region of the color image is wasteful. Therefore we
adopt coarse-to-fine strategy and a band-shaped distance map
named a distance band. The distance band is a narrow band
along the image boundary and the precise distance map isconstructed in this region using the Fast Marching Method.
In other region except the distance band, the distance value
is roughly estimated as the distance from the center of theimage boundary. In the beginning of the registration process,we use a coarse image and a wide distance band in order toalign the 2D image and the 3D model roughly. After several
iterations, we change the image to the fine one and computethe distance band with narrower width. In experiments, weused the distance band with 10 pixels and three images withdifference resolution, that is, 160x120, 320x240, 640x480.
Computation time of the distance band in each image is shownin Table I.
TABLE IC OMPUTATION TIME OF THE DISTANCE BAND BY THE FAST M ARCHING
M ETHOD .
Image size
160x120320x240
640x480

Whole region
[ms]
6.9
32.8
189.1

is detected by applying the Level Set Method[13] to thebackground subtraction image.
The computation time of the force fi for one patch on the
3D model is 0.30 [µs] for Pentium IV processor, 3.06 GHz.
The total processing time is 9.6 [ms] for one update period
including projected contour detection, calculation of force andmoment, and the steepest descent method. The average andthe variance of the registration error between the boundary
and the projected contour are 1.19[pixels] and 0.90[pixels2 ]respectively after 60 iterative calculations. Here, the numberof points on the boundary is 652, image size is 640 × 480
pixels, and the 3D model is composed of 73,192 meshes with
219,576 vertices. On the other hand, the computation timeof the force fi for the conventional point-based method [4]is 1.15 [µs], which uses the k-D tree structure to search thenearest points between the contour points of the 2D image andthe projected 3D model.
Table II show the comparison of the processing time forseveral cases in which different numbers of points on theboundary were used. From these results, it is clear that theprocessing time of the proposed method is almost constanteven if the number of points on the boundary increases. On
the other hand, the point-based approach requires 3 to 7 times
longer processing time than the proposed method. Moreover,the larger the number of points becomes, the longer processing
time is required for the point-based approach. Therefore, it isverified that the proposed 2D-3D registration algorithm hasan advantage comparing with the conventional point based
approach in views of execution time, especially when there
are a large number of points on the boundary.
Next example is the tracking experiment of the moving
object. Fig.9 shows the tracking results of the doll whichmoves on the image plane. As seen in this example, the 2D3D registration can be performed even if the object moves up
to 20.5 pixels/second. This result suggests that it is possible
to overcome the registration error due to unknown disturbancesuch as the breath or the heart beat in endoscopic images.
TABLE IIC OMPARISON OF PROCESSING TIME FOR ONE PATCH

Distance band
(10 pixels) [ms]
0.48
1.1
2.7

Number of points
on boundary
628
1265
1868
2490

V. E XPERIMENTS
In this section, we show some fundamental results of theregistration and tracking experiments using simulated images
and actual images of the endoscopic operation.
A. 2D-3D registration of simulated images
Firstly, Fig.8 shows how the proposed 2D-3D registration
algorithm works with a model of a small doll (a height of5cm). The boundary of the object on the 2D color image

913

Proposedmethod [µs]
0.30
0.30
0.30
0.30

Point-based
method [4] [µs]
1.15
1.50
1.70
2.22

B. 2D-3D registration using actual images of the endoscopic
operation
We carried out fundamental experiments for a navigation
system of the endoscopic operation.
Firstly, we tried the 2D-3D registration using actual imagesof the endoscopic operation of the cirrhosis of the liver. Fig.10
shows the initial results of the experiments. Fig.10(a) shows
the actual image of the endoscopic operation of the cirrhosis of
the liver. In this image, the liver is partially hidden by the other2D image 3D model(a) Initial position

120 pixels(a) t = 0 s(b) 1 iteration (0.01s)

color image

3D model(b) t = 7.7 s(c) 2 iterations (0.02s)(d) 3 iterations (0.03s)(c) t = 12.8 s(e) 30 iterations (0.3s)
Fig. 8.

(f) 60 iterations (0.6s)

2D-3D registration of simulation images

organs. Fig.10(b) and Fig.10(c) show the 3D model of the liverobtained from the CT images and some experimental results
for the 2D-3D registration, respectively. The registration is
successfully executed if the initial position of the 3D modelis close enough to the image of the cirrhosis of the liver.
Fig. 11 shows the tracking results of the gallbladder in video
images taken by the endoscope. Boundaries of the gallbladderon the every 2D video images are detected using by Snakes. As
seen in this result, our algorithm is able to track the gallbladdereven if the gallbladder image moves due to the change of theposition of the endoscope.

(d) t = 18.1 sFig. 9.

An example of moving object tracking

is successfully performed, then the algorithm can track amoving object stably and continuously. Therefore the proposed
algorithm is quite useful for developing a practical navigation
system since it overcomes the registration error caused byorgan’s movement due to unknown disturbance such as thebreath or the heart beat.

VI. C ONCLUSIONSACKNOWLEDGMENT

This paper described a new registration algorithm of 2D
color images and 3D geometric models for navigation systemof surgical robot. This method utilizes the 2D images andtheir distance map created by the Fast Marching Method, anddetermines a precise relative pose between 2D images and 3Dmodels. The efficiency of the proposed algorithm was verified
through the fundamental experiments using simulation images
and actual images of the endoscopic operation of liver.
We have conducted experiments with various images and
3D models and obtained almost good results. Occasionallythe performance of the proposed algorithm becomes worse
especially in case that only the partial edge of an organ isvisible from the endoscope. However, once the registration

This research was partly supported by the 21st Century
COE Program ”Reconstruction of Social Infrastructure Related
to Information Science and Electrical Engineering”, and theMinsty of Public Management, Home Affairs, Posts andTelecommunications in Japan under Strategic Information andCommunications R&D Promotion Programme (SCOPE).

914

Initial position(a) Endoscpoic image ofthe cirrhosis of the liver(b) 3D model of the liver(c-1) Initail position(c-2) after 20 iterations

Fig. 11.

(c-3) after 36 iterations(c-4) after 72 iterations

Fig. 10. 2D-3D registration using actual images of the cirrhosis of the liver

R EFERENCES[1] M. Hashizume, M. Shimada, M. Tomikawa, Y. Ikeda, I. Takahashi,
R. Abe, F. Koga, N. Gotoh, K. Konishi, S. Maehara, and K. Sugimachi,“Early experiences of endscopic procedures in general surgery assisted
by a computer-enhanced surgical system,” in Surg Endosc, 2002, pp.
1187–1191.
[2] P. Viola and W. W. III, “Alignment by maximization of mutual information,” Int. J. of Computer Vision, vol. 24, no. 2, pp. 137–154, 1997.
[3] I. Stamos and P. K. Allen, “Integration of range and image sensing forphotorealistic 3d modeling,” in Proc. of the 2000 IEEE International
Conference on Robotics and Automation, 2000, pp. 1435–1440.
[4] R. Kurazume, K. Noshino, Z. Zhang, and K. Ikeuchi, “Simultaneous 2d
images and 3d geometric model registration for texture mapping utilizing
reflectance attribute,” in Proc. of Fifth Asian Conference on Computer
Vision (ACCV), 2002, pp. 99—106.
[5] M. D. Elstrom and P. W. Smith, “Stereo-based registration of multisensor imagery for enhanced visualization of remote environments,”
in Proc. of the 1999 IEEE International Conference on Robotics andAutomation, 1999, pp. 1948–1953.
[6] K. Umeda, G. Godin, and M. Rioux, “Registration of range and color
images using gradient constraints and range intensity images,” in Proc.
of 17th International Conference on Pattern Recognition, 2004, pp. 12–
15.
[7] H. Lensch, W. Heidrich, and H.-P. Seidel, “Automated texture registration and stitching for real world models,” in In Pacific Graphics ’00,
2000, pp. 317–326.
[8] H. Lensch, W. Heidrich, and H. P. Seidel, “Hardware-accelerated silhouette matching,” in In SIGGRAPH Sketches, 2000.
[9] K. Matsushita and T. Kaneko, “Efficient and handy texture mapping on
3d surfaces,” in Comput. Graphics Forum 18, 1999, pp. 349–358.
[10] P. J. Neugebauer and K. Klein, “Texturing 3d models of real world
objects from multiple unregistered photographic views,” in Computer
Graphics Forum 18, 1999, pp. 245–256.
[11] J. Sethian, “A fast marching level set method for monotonically advancing fronts,” in Proceedings of the National Academy of Science, vol. 93,
1996, pp. 1591–1595.

915

2D-3D registration in endoscopic video images of the gallbladder.

[12] ——, Level Set Methods and Fast Marching Methods, second edition.
UK: Cambridge University Press, 1999.
[13] Y. Iwashita, R. Kurazume, T. Tsuji, K. Hara, and T. Hasegawa, “Fastimplementation of level set method and its realtime applications,” inIEEE International Conference on Systems, Man and Cybernetics 2004,2004.