IPSJ Transactions on Computer Vision and Applications

Vol. 1

209–219 (Sep. 2009)

Research Paper

A Fast Registration Method Using IP and Its
Application to Ultrasound Image Registration
Bo Zheng,†1 Ryo Ishikawa,†2 Takeshi Oishi,†1
Jun Takamatsu†3 and Katsushi Ikeuchi†1This paper presents a fast registration method based on solving an energy
minimization problem derived by implicit polynomials (IPs). Once a target
object is encoded by an IP, it will be driven fast towards a corresponding source
object along the IP’s gradient ﬂow without using point-wise correspondences.
This registration process is accelerated by a new IP transformation method.
Instead of the time-consuming transformation to a large discrete data set, the
new method can transform the polynomial coeﬃcients to maintain the same
Euclidean transformation. Its computational eﬃciency enables us to improve a
new application for real-time Ultrasound (US) pose estimation. The reported
experimental results demonstrate the capabilities of our method in overcoming
the limitations of a noisy, unconstrained, and freehand US image, resulting infast and robust registration.

1. Introduction
The task of registration aims to build a transformation relationship between agiven (source) object and a target object. It is often an intermediate but crucial
step inside the whole computer vision complexity.
Considerable eﬀorts have been made in this area. One of the most popular approaches is the Iterative Closest Point (ICP) method 2) , whose original form minimized the distances between the closest point pairs. Zhang 24) proposed using k-D
trees to speed up the closest point search. Several groups 4)–7),12)–14),17) improved
the original version by replacing the Euclidean distance or ﬁnding more robust
correspondences. However, the ICP-based methods inevitably require extra computation to ﬁnd the correspondences that are sensitive to local correspondence.
†1 Institute of Industrial Science, The University of Tokyo†2 Medical Imaging Project, Canon Inc.
†3 Nara Institute of Science and Technology

209

Avoiding the computation for searching for point-wise correspondences, thecoarse (global) registration methods 18),22) solve the problem by single (noniterative) computation, using the central and oriented information of two objects
being extracted from their polynomial coeﬃcients. These methods are very computationally eﬃcient but diﬃcult to deal with in cases where the two objects are
partially overlapped, so they are limited to some speciﬁc applications 16) .
Recently, consideration has been given to constructing a distance ﬁeld forachieving the registration, such as the work in Refs. 9), 11). These methods
spend much extra memory to preserve the distance ﬁeld for the models and thenregister them by the evolution generated in the distance ﬁeld. The methods also
avoid correspondence searching and thus speed up the convergence, but they are
liable to take too much memory space, especially for a dense model or a 3D
object, and registration is limited in the regions where distance ﬁeld has been
generated.
In this paper, we ﬁrst use an implicit polynomial (IP) to model the target object, and then we register this IP model to a discrete source object for quickly
ﬁnding the transformation relationship. The advantages of our method over theprior methods are that: i) unlike the ICP-based methods, it avoids the extra
computation for point-wise correspondences; ii) unlike the coarse registration
methods, it totally supports partial-overlapping registration; iii) unlike the registration methods of preserving a discrete distance ﬁeld, it needs very little memory
space for preserving a few IP coeﬃcients, and the algebraic model can generate
an inﬁnite distance ﬁeld to support registration in a wider space.
The remainder of this paper is organized as follows: After introducing IP modeling technique in Section 2, we present our registration method using the IPmodel in Section 3. Then in Section 4, in order to further accelerate our algorithm to be suitable for large data set cases we propose an IP transformation
method. Sections 5 and 6 present reports on experimental results and extend the
proposed method to Ultrasound image application. Conclusion and future work
follow in Section 7.
2. IP ModelingFirst we concentrate on how to model the target object with an implicit poly-c 2009 Information Processing Society of Japan210

A Fast Registration Method Using IP and Its Application to Ultrasound Image

nomial. In this section, we ﬁrst introduce the mathematical formulation of IP,
and then present a modeling technique.
2.1 Formulation of IP
IP is an implicit function deﬁned in a multivariable polynomial form. For
example, a 3D IP of degree n can be denoted by:
fn (x) =aijk xi y j z k = (1. . . z n)(a000 a100 . . . a00n )T , (1) x 
0≤i,j,k;i+j+k≤nm(x)Tawhere variable x = (x y z) denotes the Cartesian coordinate of a data point andan IP can be represented as an inner product between the monomial vector andthe coeﬃcient vector as m(x)T a. We employ the inverse lexicographical order 22)
for monomial indices {i, j, k}. The homogeneous binary polynomial of degree rin x, y, and z, i+j+k=r aijk xi y j z k , is called the r-th degree form of the IP.
2.2 Distance Field Constraints
Although in the past two decades, IP was attractive to shape representation for
object recognition purpose 3),10),19)–21),26) , in this paper we try to take advantage
of IP to model the distance ﬁeld in order to solve the registration problem. To dothis, we deﬁne optimization constraints that let a polynomial be a given signed
distance D at the data point x as:
fn (x) = D,(2)
where D is deﬁned as:⎧⎪0, x ∈ M⎨
D(x) =(3)+d(x, M), x ∈ Minner ,⎪⎩ −d(x, M), x ∈ M
outerwhere we suppose M represents a closed surface of target object; Minner andMouter represent its ﬁnite inner and outer regions respectively; and d(x, M) isa the minimum Euclidean distance: d(x, M) = minxi ∈M ||x − xi ||.
For constructing the signed distance ﬁeld (SDF), we employ the approximation
method proposed in Ref. 15) that generates two level sets of SDF along the surface
normals. The diﬀerence to their method is just that we generate SDF aroundthe surface more densely and in a wider region.

IPSJ Transactions on Computer Vision and Applications

Vol. 1

209–219 (Sep. 2009)

Fig. 1 (a) noisy data set; (b) ﬁtting result of 4-degree IP; (c) data set with missing data; (d) ﬁtting result of 4-degree IP; (e) resulting gradient ﬁeld of IP model of the same object.

2.3 Fitting
Using the optimization constraints (2), the IP coeﬃcients can be solved by thelinear least-squares method formulated as:
a = (M T M )−1 M T d,(4)
where matrix M consists of the monomial vectors, i.e., the i-th row of M is m(xi )(see Eq. (1)), and d is the right-hand vector whose ith element is D(xi ). For solving Eq. (4), we adopt our previous modeling method proposed in Ref. 26), whichallows us to adaptively determine the moderate degree for IP diﬀerent objects
according to the shape complexity. Figure 1 shows the modeling examples: forthe data set shown in Figs. 1 (a) and (c), the resulting ﬁts are shown in Figs. 1 (b)fand (d); Figure 1 (e) shows the resulting gradient ﬁeld generated by f · f . Wecan see that two important characteristics are: i) IP is robust against noises andocclusions, and ii) IP can provide a stable gradient ﬁeld around its zero-level set.
3. IP Registration
Having an IP model representing the target object, the next problem is how
to register this IP to the source object. To achieve that, we ﬁrst deﬁne theregistration problem as an energy minimization problem, and then, by minimizing
this energy function, the motion of the source object can be driven to the target
model (IP model).
3.1 Energy Functional
To achieve the registration, an energy functional E which will be minimized to

c 2009 Information Processing Society of Japan211

A Fast Registration Method Using IP and Its Application to Ultrasound Imageﬁnd the proper transformation parameters deﬁned as:
p = arg min E(p).

(5)

pwhere p is the transformation parameters (p ∈ R6 for rigid transformation). In
general, the energy functional E evaluates the registration by minimizing the
distance between the data set and IP, deﬁned as:dist(T (p, xi ), fn ), ∀xi ∈ Ω,
(6)
E=
i

where T (p, xi ) is the function R3 → R3 , which returns the transformed point ofxi by the rigid-transform operation respected to parameter p; dist(x, fn ) meansa certain distance from the data point x to the zero set of fn ; and Ω represents
the 3D region of the source model.
There are multiple choices for the distance dist(x, fn ) in Eq. (6), and the common one is to use L2 norm. Thus it becomes possible to form least-squares
regression with dist(x, fn ) = fn (x)2 . However, in this work, we choose moremeaningful approximation for the distance representation of IP (see Ref. 21)):fn (x)2
dist(x, fn ) =, ∀x ∈ Ω.
(7)
 fn (x) 2
3.2 Minimizing Energy Functional
To minimize the energy functional (6), the method employs the following two
steps. First, it minimizes the function without any constraint in the transformation. This means all points can move freely towards IP along their gradientsduring the ﬁrst minimization. Next, it determines the transformation parameters to maintain the Euclidean transformation. These two steps are repeatedalternately until the convergence. The eﬃciency of this minimization beneﬁts
from its success in avoiding time-consuming computation for ﬁnding point-wisecorrespondence.
3.2.1 First Step: Free Deformation
At the ﬁrst step, by calculus of variations (see Ref. 8)), the Gateaux derivative(ﬁrst variation) of the functional E in Eq. (6) to point x can be approximately
formulated as:∂E∂dist(x, f )fn (x)=≈ 2f (x)= 2g(x),(8)∂x∂x
 fn (x) 2

IPSJ Transactions on Computer Vision and Applications

Vol. 1

209–219 (Sep. 2009)

if we assume  fn (x)  is a constant and let the gradient vector g(x) =fn (x)f (x) f2.
n (x)
Therefore, we need to minimize this functional to satisfy the Euler-Lagrange
equation ∂E∂ x = 0. Thus the steepest descent process is executed in the following
gradient ﬂow for each point x:∂x≈ −2g(x),(9)∂twhere t denotes the time step. From the view of implementation, each point canbe updated as:
xk+1 = xk − 2g(xk ),(10)

where xk denotes the point at k-th time step.
Similarly, if we let X ∈ RN ×3 represent the integral data set consisting of thepoint coordinates as rows, and G ∈ RN ×3 be the integral gradients of X, thenX k at the k-th time step can be updated as the following recurrence:
X k+1 = X k − 2G(X).

(11)

This shows the fact that, under this operation, all points in the discrete data
set of source object will be moved towards the zero set of polynomials along their
gradients.
3.2.2 Second Step: Rigid Deformation
In the ﬁrst step, although every point in the source model can move close to
the IP surface, the transformation dose not maintain the rigid transformationthat is necessary for rigid parameter estimation. Therefore, in the second step,
we alter the outcome of the ﬁrst step so that the result is rigid transformation.
For achieving this, since from the Eq. (11) we can obtain data sets represented
by two matrices X k and X k+1 in two continuous steps, the covariant matrix ofthese two matrices can be calculated as:
A = (X k − X̄ k )T (X k+1 − X̄ k+1 ),(12)

where X̄ is a matrix in which each row consists of the mean value (mass center
point) of X. If A (∈ R3×3 ) is SVD-decomposed as A = U SV T , where S is the
diagonal matrix and U and V are the unitary matrices, then rigid transformation

c 2009 Information Processing Society of Japan

212

A Fast Registration Method Using IP and Its Application to Ultrasound Imageis obtained by:
R = U V T , t = X̄ k − X̄ k+1 RT ,(13)

where R and t are rotation and translation parameters respectively (see Ref. 23)).
Thus, with R and t, we can update the transformation in the ﬁrst step shown in
Eq. (11) to be rigid transformation as in:
X k+1 = X k RT + t.

(14)

Therefore, the above two steps can be alternately repeated until the discrete
data set of the source model is moved near the IP model with desired accuracy.
4. Acceleration Using IP Transformation
In this section, we present an acceleration method that makes registration more
eﬃcient when the number of source data points is much larger than the numberof IP coeﬃcients of the target model. Since in this case moving the large-scaledata sets to the IP model is obviously time-consuming, instead we inversely move
the IP model towards the discrete model.
To drive this motion of IP requires a transformation method through which
the IP’s coeﬃcients can be transformed to maintain an inverse transformation towhat was done in Eq. (14). Therefore the acceleration method can be described
as a process where, instead of updating the discrete data set in Eq. (14), we
update the coeﬃcients of IP at each step:
ak+1 = V (R, t)ak ,5. Numerical Experiments
In this section, we report results of experiments dealing with some synthetic
data sets to evaluate the method on computational performance. All the experiments were implemented in Matlab 7 with a PC having an Intel core 2 CPU,
2.4 GHz, and 2 GB memory.
5.1 Registration for Synthesized Data Set
Figure 2 shows the registrations in a simple case that registers two synthetic
data sets to target object “bunny” modeled by an IP of degree 8. Also their(15)

where we suppose ak+1 and ak are IP coeﬃcients at the (k + 1) and k-th time
step, and V is a square transformation matrix corresponding to the Euclideantransformation of rotation R and translation t. Thus the remained problem is
how to calculate the transformation matrix V in Eq. (15).
4.1 IP Transformation Method
To our knowledge, although prior literatures such as Ref. 22) have introduced
the existence of IP transformation matrix V from the aspect for extracting algebraic invariants of IP, unfortunately computational implementation of IP transformation was not explicitly described. Tarel, et al. 18) showed a tensor-based

IPSJ Transactions on Computer Vision and Applicationstransformation that represents 3D IP in 4D homogeneous form, and then transformed the tensor elements. The method is good for IP invariants extraction but
still suﬀers from heavy computational cost for tensor transformation, since fortransforming an r-degree 3D IP, it need to transform 4d tensor elements.
We describe our implementation based on symbolic computation in Appendix A.1. This solution can be seen as a kind of incremental calculationthat increases, beginning from the transformation for ﬁrst-degree terms, until
it obtains a transformation for n-th degree terms.

Vol. 1

209–219 (Sep. 2009)

(a)

(b)

Fig. 2 IP registration with 2D and 3D feature points: (a) registration process between“bunny” IP model and points (blue points: initial position; green points: selected
mid-steppositions; red points: ﬁnal position). (b) the mean absolute distance1|dist(x, f )| versus iteration number.
N

c 2009 Information Processing Society of Japan213

A Fast Registration Method Using IP and Its Application to Ultrasound Image

energy convergence behaviors is shown in Fig. 2 (b) respectively. From this result,
we can see that our method can totally support the partially overlapped objects.
5.2 Comparison between ICP Method and Our Method
We solved a registration task shown in which we initially rotated and translated the target model to the position shown as blue points in Fig. 3. Then
we tested three methods, standard ICP method 2) , ICP method with KD tree 24)
and our method, to let the process iterate until it satisfed the same accuracy forregistration.
Since they are set with the same stopping criterion (same threshold for registration accuracy), we only compared the computational time and memory inTables 1, 2 and 3, which show the registration tasks using 1k points, 2.5k points
and 10k points respectively for both target and source data sets.
It is worth noting that, in the factor of accuracy, IP makes it diﬃcult to achieve−→
Fig. 3 Illustration for registration: blue points for target object and red points for source
object.
Table 1 Registration using 1k points.
CPU Time (sec)MemoryICP
11.95
2K points

ICP with KD tree
0.91
2K points + KD treeIP
0.43
1K points + 165 coef.

Table 2 Registration using 2.5k points.
CPU Time (sec)MemoryICP
30.54
5K points

CPU Time (sec)MemoryICP
123.0120K points

ICP with KD tree
6.56
5K points + KD treeIP
0.51
2.5K points + 165 coef.

Table 3 Registration using 10k points.
ICP with KD tree
35.12
20K points + KD tree

IPSJ Transactions on Computer Vision and Applications

IP
1.07
10K points + 165 coef.

Vol. 1

209–219 (Sep. 2009)better registration accuracy, because IP often models the data set “coarsely” by
the global algebraic description. But, on the contrary, the IP model has few
coeﬃcients and little robustness against noises, which make it faster and morestable when facing noise-condition applications, such as the ultrasound image
application described in the following section.
6. Application to US Image Registration
Ultrasound (US) imaging is widely used for assistance with surgical operations
and clinical diagnosis. The relative pose estimation for a US image to anothermodality, such as MRI or CT, is desired and helpful for diagnostic guidance. But
such estimation suﬀers from diﬃculties due to poor image quality with speckle
noises, low signal-to-noise ratio, uniform brightness, and the fact that only crosssection images can be obtained. Fortunately, considering the characteristics ofour method as described above, we can confront this registration problem bymaking use of IP.
6.1 Registration for US Images
Registration for US images can be simply viewed as 2D-3D registration, if thesource and target models, 2D feature points, and a 3D organ model represented
by an IP are given.
Three-dimensional organ models are often obtained from a medical modality,
such as CT, MRI or PET, by 3D image segmentation. For example, segmentation using the Graph Cut or Level set method can extract 3D surface models
from volumetric data. Then the 3D surface models can be modeled by IPs withthe modeling method described in Section 2, and the modeling process can beachieved oﬄine. On the other hand, 2D feature points extraction must be done
online. The 2D segmentation method for a single image, such as the Level set
method, is not suitable to a real-time requirement.
In this paper, we employed our previous method 25) to register an IP model toa single US image. The process of this method can be simply viewed as: i) using
an edge detector, it extracts two kinds of feature points around the IP model’ssurface from the US image plane: boundary-like feature points and region-likefeature points, and ii) by our registration method described in Sections 3 and 4,
it drives the motion of boundary-like feature points to the IP model’s surface

c 2009 Information Processing Society of Japan214

A Fast Registration Method Using IP and Its Application to Ultrasound Image

iteration 1#

iteration 3#

iteration 6#

iteration 8#

Fig. 4 Pose estimation for a single US image.

(zero level set) and also drives the motion of region-like feature points into the
IP model’s inner region (negative level set), but the registration process is thesame.
6.2 Experimental Results
6.2.1 Pose Estimation for Phantom ATS514
Figure 4 shows pose estimation process for a US image, where the US images
are obtained by measuring phantom ATS514 1) consisting of the cylinder-likeshapes inside in diﬀerent materials. Then we use a cylinder-like IP model (2degree IP shown in red in Fig. 4) to register the US images. The initial positionis set as: the cylinder-like IP model crosses the US image plane with an angle of85◦ . For this image, 8 iterations are required, and in the process the motion of
the IP is driven to ﬁt to the hole-like region in the image. We show this process
in Fig. 4 with cross-section contour and relative pose between the IP model andthe US plane at each iteration.
Figure 5 shows three US image frames with the same initialization and in
the second row of Fig. 5 the iteration number required and ﬁnal relative pose areshown. The third row shows the cross-section contour between the IP and the
US plane at the ﬁnal iteration. For each case, the consumed CPU time is within
30 ms.
6.2.2 Pose Estimation for Duck Toy
Figure 7 shows another result of the US image pose estimation where the images
are obtained by scanning a duck toy made of rubber in the cistern shown in

IPSJ Transactions on Computer Vision and Applications

Vol. 1

209–219 (Sep. 2009)

Fig. 5 Pose estimation for phantom ATS514. First row: initial position; Second row: ﬁnal
position; Third row: cross-section contour of IP and US image (white points).

Fig. 6 Left: photograph of measuring a duck toy. Right: IP model of duck object.

Fig. 6 left. We modeled the duck object with an 8-degree IP shown in Fig. 6
right. Figure 7 shows the relative position and the cross-section contour at each
iteration resulting from the registration. In this case the consumed CPU time isabout 180 ms.
6.2.3 Pose Estimation and Tracking for CT Data
The third example is to test by solving a registration problem between real CT
data and a US image. To do this, we ﬁrst segmented the CT data to obtain thedesired organ object by modern segmentation methods such as Graph Cut (seeFig. 8). Then we model the organ object by a 4-degree IP as shown in Fig. 8right. Finding the relative pose to CT data for a US image now becomes ﬁnding
the relative pose to the IP model.
The ﬁrst example is pose estimation for a single US image shown in Fig. 9 where
we show the registration process with relative pose and cross-section contour at

c 2009 Information Processing Society of Japan215

A Fast Registration Method Using IP and Its Application to Ultrasound Image

Fig. 9 US pose estimation for CT data with 4-degree IP.

Fig. 7 Pose estimation for duck toy. First row: relative poses of IP and US images at each
iteration. Second row: cross-section contour at each iteration.

Fig. 8 Segmenting kidney data from the CT scans and modeling with a 4-degree IP.

each selected iteration.
We also tried a tracking problem for a US image sequence by tracking theposition of the kidney object with the above IP model. Results are shown in
Fig. 10 where we show the cross-section contours of IP for selected frames.
7. Conclusion and Future WorkA fast registration method using the IP gradient ﬁeld has been proposed. It is

IPSJ Transactions on Computer Vision and Applications

Vol. 1

209–219 (Sep. 2009)

Fig. 10 Tracking for US image sequence. Cross-section contours of frames are shown inblack points.

computationally eﬃcient due to the adoption of a two-step minimization process
without point matching, and acceleration supported by IP transformation. We
demonstrated its high performance by applying it into real-time US image registration, which shows robustness against the presence of weak boundaries andstrong noise, and is possibility suitable for real-time implementation. In future
work, we are considering modeling multiple organ objects using IPs to enhance
the robustness of registration.
Acknowledgments This work was supported by Canon Inc., under the
project: Physics-based vision theories for next-generation medical image processing.
References
1) http://www.atslabs.com
2) Besl, P. and McKay, N.: A Method for Registration of 3-D Shapes, IEEE Trans.
Patt. Anal. Mach. Intell. (TPAMI ), Vol.14, No.2, pp.239–256 (1992).
3) Blane, M., Lei, Z.B. and Cooper, D.B.: The 3L Algorithm for Fitting Implicit
Polynomial Curves and Surfaces to Data, IEEE Trans. Patt. Anal. Mach. Intell.

c 2009 Information Processing Society of Japan216

A Fast Registration Method Using IP and Its Application to Ultrasound Image(TPAMI ), Vol.22, No.3, pp.298–313 (2000).
4) Boughorbel, F., Koschan, A., Abidi, B. and Abidi, M.: Gaussian Energy Functions for Registration without Correspondences, 17th International Conference on
Pattern Recognition (ICPR) (2004).
5) Champleboux, G., Lavallee, S., Szeliski, R. and Brunie, L.: From Accurate Range
Imaging Sensor Calibration to Accurate Model-Based 3D Object Localization, Proc.
IEEE Conf. Computer Vision and Patt. Rec. (CVPR), pp.83–89 (1992).
6) Chen, Y. and Medioni, G.: Surface Description of Complex Objects from Multiple
Range Images, Proc. IEEE Conf. Computer Vision and Patt. Rec. (CVPR), pp.153–
158 (1994).
7) Chetverikov, D., Svirko, D., Stepanov, D. and Krsek, P.: The trimmed iterative
closest point algorithm, International Conference on Pattern Recognition (ICPR)(2002).
8) Evans, L.: Partial Diﬀerential Equations, Providence: American Mathematical Society (1998).
9) Iwashita, Y., Kurazume, R., Konishi, K., Nakamoto, M., Aburaya, N., Sato,Y., Hashizume, M. and Hasegawa, T.: Fast Model-Image Registration using 2D
Distance Map for Surgical Navigation System, Advanced Robotics, Vol.21, No.7,
pp.751–770 (2007).
10) Keren, D. and Cooper, D.: Describing Complicated Objects by Implicit Polynomials, IEEE Trans. Patt. Anal. Mach. Intell. (TPAMI ), Vol.16, No.1, pp.38–53(1994).
11) Munim, H.E. and Farag, A.: Shape Representation and Registration using Vector
Distance Functions, Proc. IEEE Conf. Computer Vision and Patt. Rec. (CVPR),
pp.Y1–Y8 (2007).
12) Nguyen, V.-D., Nzomigni, V. and Stewart, C.V.: Fast and Robust Registration of
3D Surfaces Using Low Curvature Patches, 2nd International Conference on 3-D
Imaging and Modeling (3DIM ’99 ), p.201 (1999).
13) Oishi, T., Nakazawa, A., Kurazume, R. and Ikeuchi, K.: Fast Simultaneous Alignment of Multiple Range Images using Index Images, Proc. 5th International Conference on 3-D Digital Imaging and Modeling (3DIM 2005 ), pp.476–483 (2005).
14) Rusinkiewicz, S. and Levoy, M.: Eﬃcient Variants of the ICP Algorithm, International Conference on 3D Digital Imaging and Modeling (3DIM 2001 ) (2001).
15) Sahin, T. and Unel, M.: Fitting Globally Stabilized Algebraic Surfaces to Range
Data, Proc. IEEE Conf. Int. Conf. on Computer Vision (ICCV ), Vol.2, pp.1083–
1088 (2005).
16) Salvi, J., Matabosch, C., Foﬁ, D. and Forest, J.: A review of recent range image registration methods with accuracy evaluation, Image and Vision Computing,
Vol.25, No.5, pp.578–596 (2007).
17) Stewart, C.: Covariance-based registration, Technical Report RPI-CS-TR 02-8,
Department of Computer Science Rensselaer Polytechnic Institute (2002).

IPSJ Transactions on Computer Vision and Applications

Vol. 1

209–219 (Sep. 2009)18) Tarel, J.-P., Civi, H. and Cooper, D.B.: Pose Estimation of Free-Form 3D Objects
without Point Matching using Algebraic Surface Models, Proc. IEEE Workshop
Model Based 3D Image Analysis, pp.13–21 (1998).
19) Tarel, J.-P. and Cooper, D.B.: The Complex Representation of Algebraic Curves
and Its Simple Exploitation for Pose Estimation and Invariant Recognition, IEEE
Trans. Patt. Anal. Mach. Intell. (TPAMI ), Vol.22, No.7, pp.663–674 (2000).
20) Tasdizen, T., Tarel, J.-P. and Cooper, D.B.: Improving the Stability of Algebraic
Curves for Applications, IEEE Trans. Imag. Proc. (TIP ), Vol.9, No.3, pp.405–416(2000).
21) Taubin, G.: Estimation of Planar Curves, Surfaces and Nonplanar Space Curves
Deﬁned by Implicit Equations with Applications to Edge and Range Image Segmentation, IEEE Trans. Patt. Anal. Mach. Intell. (TPAMI ), Vol.13, No.11, pp.1115–
1138 (1991).
22) Taubin, G. and Cooper, D.: Symbolic and Numerical Computation for Artiﬁcial
Intelligence, chapter 6, Computational Mathematics and Applications. Academic
Press (1992).
23) Umeyama, S.: Least-Squares Estimation of Transformation Parameters Between
Two Point Patterns, IEEE Trans. Patt. Anal. Mach. Intell. (TPAMI ), Vol.13, No.4,
pp.376–380 (1991).
24) Zhang, Z.: Iterative point matching for registration of free-form curves and surfaces, Int. J. Computer Vision (IJCV ), Vol.13, No.2, pp.119–152 (1992).
25) Zheng, B., Ishikawa, R., Oishi, T., Takamatsu, J. and Ikeuchi, K.: 6-DOF Pose
Estimation form Single Ultrasound Image using 3D IP Models, Proc. IEEE Comput.
Vision Patt. Recog. Workshop 2008, pp.1–8 (2008).
26) Zheng, B., Takamatsu, J. and Ikeuchi, K.: Adaptively Determining Degrees of
Implicit Polynomial Curves and Surfaces, Proc. Conf. Asian Conf. on Computer
Vision (ACCV ), Vol.LNCS 4844, pp.289–300 (2007).

Appendix
A.1 Implementation for IP Transformation
A.1.1 Pure Rotation Transformation
First, let us explain some notations for polynomial operations which wereﬁrst introduced by Taubin and Cooper (see Ref. 22)). Let a coeﬃcientΦijk
aijk of an IP in Eq. (1) be presented as i!j!k!
, and a vector Φ[l] =√Φl00
l!0!0!

√Φl−1,1,0(l−1)!1!0!

00l. . . √Φ0!0!l!

T

be according to the l-th form of IP. An op-

eration on this vector is deﬁned as: Φ[l,m] = Φ[l]  ΦT[m] , where  represents theclassic matrix multiplication, with the diﬀerence that the individual elements

c 2009 Information Processing Society of Japan


217

√Φijk

i!j!k!

A Fast Registration Method Using IP and Its Application to Ultrasound ImageΦabc√and √Φa!b!c!
in each vector are multiplied as √i+a,j+b,k+c. For example,

Φ[1] =

i!j!k! a!b!c!

ΦΦΦ√ 100 √ 010 √ 001
1!0!0! 0!1!0! 0!0!1!

and

⎛

Φ200
⎜
T
Φ[1,1] = Φ[1]  Φ[1] = ⎝ Φ110
Φ101

T

Φ110
Φ020
Φ011

,

(16)
⎞
Φ101
⎟
Φ011 ⎠ .
Φ002

(17)

In Ref. 22), it was pointed out that under a non-singular coordinate transformation A, e.g., Φ[l] = A[l] Φ[l] , the transformed coeﬃcient matrix is given by−1Φ[l,m] = A−T[l] Φ[l,m] A[m] ,(18)

where A[l] is a non-singular hl × hl transformation matrix (assuming hl be thenumber of monomial in the l-th form), and A[m] is the same. From this equation,
we know that if A[l] and A[m] are given, then a linear relationship between Φ[l,m]and Φ[l,m] can be found; that is, the element-wise correspondence can be linearly
expressed asΦijk =
aαβ Φijk = aT(19)
α Φ[l+m] .
β;i,j,k≤l+mThen since all of the elements in vector Φ[l+m] are contained in matrix Φ[l,m] ,and similarly all of the elements in Φ[l+m] are contained in Φ[l,m] , a new linear
correspondence can be constructed between Φ[l+m] and Φ[l+m] by arranging thenecessary elements in the right order into Φ[l+m] asΦ[l+m] = (Φl+m,0,0 , Φl+m−1,1,0 , Φl+m−2,2,0 , . . .)T
T
T
T
= (aT
1 Φ[l+m] , a2 Φ[l+m] , a3 Φ[l+m] , . . .)= A[l+m] Φ[l+m] ,(20)

if the (l + m)-th transformation matrix is constructed as: A[l+m] =
T
T
T
[aT
1 , a2 , a3 , . . .] .
 
 
Φ020
Φ002 TAs a simple example, the elements of Φ[2] = Φ200
Φ
Φ
Φ
110
101
011
2
2
2
are contained in Φ[1,1] shown in Eq. (17), and the elements of Φ[2] are contained

IPSJ Transactions on Computer Vision and Applications

Vol. 1

209–219 (Sep. 2009)T(supposing R[1] is a purein Φ[1,1] . From the equation Φ[1,1] = R[1] Φ[1,1] R[1]rotation matrix), it is easy to ﬁnd out the relationship of a linear expansion for
the elements, e.g., the ﬁrst element of Φ[1,1] is expanded as: Φ200 = r1 Φ[1,1] rT
1 =
2
2
2
r11
Φ200 + 2r11 r12 Φ110 + 2r11 r13 Φ101 + r12
Φ020 + 2r12 r13 Φ011 + r13
Φ002 , where r1
is the ﬁrst row of R[1] and rij is an element of R[1] . Then since all the elementsin Φ[1,1] are contained in Φ[2] , we can ﬁnd a linear relationship 12 Φ200 = aT
1 Φ[2] ,22
2 T
where a1 = (r11 r11 r12 r11 r13 r12 r12 r13 r13 ) is viewed as the ﬁrst row of A[2] .
Therefore, transformation matrices for diﬀerent forms can be calculated in an
incremental manner as

Algorithm
1) Initialization: given A[1] ;−12) Φ[1,1] = A−T[1] Φ[1,1] A[1] −→ A[2] ;
3)
n)

−1
Φ[1,2] = A−T[1] Φ[1,2] A[2] −→ A[3] ;..
.
−1Φ[1,n−1] = A−TΦ[1,n−1] A[n−1] −→ A[n] .
[1]Here, −→ represents constructing the transformation matrix for the next step,after ﬁnding out the linear mapping relationship between Φ[1,l] and Φ[1,l] . Note if
the initial argument A[1] represents the pure rotation, then A[l] is an orthogonal
matrix (see Ref. 22)), and thus Eq. (18) degenerates toΦ[l,m] = A[l] Φ[l,m] AT[m] .
A.1.2 Aﬃne Transformation
By using homogeneous coordinates, an IP of n variables in Euclidean spacecan be described in projective space by a corresponding homogeneous IP of n + 1
variables. To convert a ternary (i.e., 3D) IP of degree dd(x, y, z) =aijk xi y j z k(21)
f3D
0≤i,j,k;i+j+k≤dinto its homogeneous representation, a new component t = 1 is added to the 3D
IP as

d
f4D
(x, y, z, t) =aijk xi y j z k tl .
(22)0≤i,j,k,l;i+j+k+l=dc 2009 Information Processing Society of Japan218

A Fast Registration Method Using IP and Its Application to Ultrasound Image

Therefore a homogeneous polynomial corresponding to a 3D IP of degree d is aform of degree d in a 4D IP. Thus the procedure mentioned in the last section
can be used to transform the homogeneous IP (4D IP). The diﬀerent is that theincremental algorithm starts from a 4 × 4 aﬃne transformation matrix A[1] , e.g.,in the Euclidean case
R3×3 t,A[1] =
01×3 1

Ryo Ishikawa received his M.E. and D.E. degrees in systemengineering from University of Electro-Communications, Japan,
in 2002 and 2005 respectively. He has been working in Canon Inc.
since 2005 and now is a researcher of Medical Imaging Project ofCanon. He is a member of IEEE, JSST and SOFT.

where R3×3 and t are a pure 3D rotation matrix and a translation vector. As
a result, a 3D IP of degree d can be aﬃne-transformed only by A[d] , once it isworked out. Furthermore, the incremental scheme can be modiﬁed for saving thecomputational cost, e.g., for calculating A[9] we can select the incremental orderas:A[1] → A[2] → A[4] → A[8] → A[9] .
(Received October 5, 2008)(Accepted February 24, 2009)(Released September 24, 2009)

Takeshi Oishi is a Project Lecturer at the Interfaculty Initiative in Information Studies, the University of Tokyo. He received
the B.Eng. degree in Electrical Engineering from Keio University
in 1999, and the Ph.D. degree in Interdisciplinary Information
Studies from the University of Tokyo in 2005. His research interests are in modeling from reality, digital archive and augmented
reality.

(Communicated by Ryo Kurazume)
Bo Zheng received his M.E. and Ph.D. in the Applied Physics
School of Engineering in 2005 and Information Science and Technology in 2008 respectively, from the University of Tokyo. He
worked for Fujitsu Co., Ltd. during the Summer of 2005. He iscurrently a project research associate in the Institute of Industrial
Science, the University of Tokyo. His research interests include algebraic representation, object recognition, and solving large-scaleand sparse linear systems.

IPSJ Transactions on Computer Vision and Applications

Vol. 1

209–219 (Sep. 2009)

Jun Takamatsu received the Ph.D. degree in Computer Science from the University of Tokyo, Japan, in 2003. After working
in the Institute of Industrial Science, the University of Tokyo, hejoined Nara Institute of Science and Technology (NAIST), Japan
in 2008. At present, he is an associate professor in Graduate
School of Information Science, NAIST. His research interests
are in 3D shape modeling and analysis, physics-based vision, androbotics including learning-from-observation and task planning using feasible motion analysis.

c 2009 Information Processing Society of Japan219

A Fast Registration Method Using IP and Its Application to Ultrasound Image

Katsushi Ikeuchi received the Ph.D. degree in Information
Engineering from the University of Tokyo in 1978. After working
at the Artiﬁcial Intelligence Laboratory, Massachusetts Instituteof Technology, Electrotechnical Laboratory of Japanese Ministryof International Trade and Industry, and the School of Computer
Science, Carnegie Mellon University, he joined the university in
1996. At present, he is a Professor at the Interfaculty Initiativein Information Studies, the University of Tokyo. His research interests include
modeling-from-reality, object recognition, and learning-from-observation. He hasreceived several awards, including the David Marr Prize in computational vision,and IEEE R&A K-S Fu memorial best transaction paper award. In addition, in
1992, his paper, “Numerical Shape from Shading and Occluding Boundaries” wasselected as one of the most inﬂuential papers to have appeared in Artiﬁcial Intelligence Journal within the past ten years. His IEEE activities include General
Chair, IROS95, ITSC00, IV01; Program Chair, CVPR96, ICCV03; Associate
Editor, IEEE TRA, IEEE TPAMI; distinguished lecture SPS (2000–2002), RAS(2004–2006). Dr. Ikeuchi was selected to IEEE Fellow in 1998. He is the EIC ofInternational Journal of Computer Vision.

IPSJ Transactions on Computer Vision and Applications

Vol. 1

209–219 (Sep. 2009)

c 2009 Information Processing Society of Japan