Fast Simultaneous Alignment of Multiple Range Images
Using Index Images
Takeshi Oishi
Institute of Industrial Science,
The University of Tokyo

Atsushi Nakazawa
Osaka University
nakazawa@ime.cmc.osaka-u.ac.jp

oishi@cvl.iis.u-tokyo.ac.jp

Ryo Kurazume
Kyushu University
kurazume@is.kyushu-u.ac.jp

Katsushi Ikeuchi
Institute of Industrial Science,
The University of Tokyo
ki@cvl.iis.u-tokyo.ac.jp

Abstract
This paper describes a fast, simultaneous alignment
method for a large number of range images. Generally themost time-consuming task in aligning range images issearching corresponding points. The fastest searching
method is the “Inverse Calibration” method. However,this method requires pre-computed look-up tables andprecise sensor parameters. We propose a fast searching
method using “index images,” which work as look-up
tables and are rapidly created without any sensor
parameters by using graphics hardware. To accelerate thecomputation to estimate rigid transformations, weemployed a linear error evaluation method. When thenumber of range images increases, the computation time
for solving the linear equations becomes too long becauseof the large size of the coefficient matrix. On the other
hand, the coefficient matrix has the characteristic tobecome sparser as the number of range images increases.
Thus, we applied the Incomplete Cholesky Conjugate
Gradient (ICCG) method to solve the equations and foundthat the ICCG greatly accelerates the matrix operation by
pre-conditioningthecoefficientmatrix.
Some
experimental results in which a large number of range
images are aligned demonstrate the effectiveness of our
method.

1. Introduction
The modeling surface geometry of a real object is an
important issue in the fields of computer vision andcomputer graphics. The theory of modeling from reality
was proposed few years ago, and related technologies
have been highly developed in recent years [1]. Preserving
cultural heritage objects and art works in the world is oneof the important applications of these technologies, as

discussed in [2, 3]. In these projects, laser range sensors
are utilized for digitizing the surface geometry of theobjects because of the high accuracy of the sensors.
While the modeling technologies have been greatly
improved, there are still some problems in creating a 3Dmodel from a large number of range images. Because alaser range sensor can measure only the visible surface, itis necessary to take range images from many different
directions. Once the scanning has been completed, allrange images have to be aligned into a common
coordinate system. If an object is small enough to be put
on a turntable, it is easy to obtain the relative positions of
the range images. But in cases where an object is a large
statue, for example, it is difficult to record the accurate
position and direction of the laser range sensor. Therefore,
computation to obtain the relative positions of the range
images is required.
Many methods of aligning range images have been
proposed. These algorithms are based on the iterative
closest point (ICP) proposed by Besl [4] and are adaptedfrom the method proposed by Chen [5]. With ICP,
corresponding points are searched for as the closest points
between two range images, and a transformation matrix is
computed so that the mean square error of thecorresponding points is minimized. The computation isiterated until the mean square error falls below thethreshold value. In Chen’s method, the relative positionsof range images are calculated so that the distance
between vertices and the corresponding patches isminimized. In addition, there is a method to search forcorrespondences by projecting the points along with the
ray direction [6, 7]. Since the ICP algorithm tends to beaffected by false matching and noise, Masuda et al.
proposed a robust method that uses random sampling andthe Least Median Squares Estimation method (LMedS)[8].

When the number of range images is very large, amethod that simultaneously aligns range images isrequired. The algorithms described above align two range
images; when using these algorithms, error accumulation
increases as the number of range images increases. In suchcases, a method that simultaneously aligns range images isuseful. Neugebauer et al. proposed a simultaneous
registration method that adopted projection search ofcorrespondences and point-plane error metric [9].
Benjamaa et al. extended the method proposed byBergevin et al. [10] and implemented a simultaneousalignment method while they accelerated the pair-wisealignment algorithm by using multi z-buffers [11].
Although various methods have been proposed, theproblem for every method is the computation cost ofcorrespondence search. If the number of vertices of two
range images is equally assumed to be N by the original
ICP, their complexity is O(N2) since correspondences aresearched for in all vertices. In order to accelerate ICP,
there are techniques [12, 13] that use Kd-trees and thatnarrow the search range by using data cache [14, 15, 16].
However, the complexity of Kd-tree search is O(NlogN).
That is, sufficient acceleration cannot be achieved bythese algorithms. The computational complexity of theinverse calibration method proposed by Blais is O(N) [6].
However, this method requires precise sensor parameters(intrinsic parameters of CCD camera, parameters ofscanning mechanism) and pre-computed look-up tables. In
addition, the creation of the look-up tables is very time
consuming because Euclidian distances between eachelement of a table and every ray of sampled points have to
be calculated.
Another problem in aligning a large number of range
images is the computation cost of matrix operations inwhich rigid transformations of range images are computed.
To directly solve a non-linear least squares problem isvery time consuming [13]. In this case, the linearizedalgorithm is effective in dealing with a large data set [7].
However, the computation time to solve the linear
equations with conventional solvers (SVD, Cholesky
decomposition, etc.) rapidly increases as the number ofrange images increases because the coefficient matrix
becomes very large.
We propose a fast method to align a large number ofrange images simultaneously. Our method has threecharacteristics. 1) The process of searching corresponding
points is accelerated by using index images, which are
rapidly created without sensor parameters. 2) The method
employs the point-plane error metric and linearized error
evaluation. 3) An iterative solver (incomplete Cholesky
conjugate gradient method) is applied in order toaccelerate the computation of the rigid transformations. In
Section 2, the details of our algorithm are described. Someexperimental results that demonstrate the effectiveness of
our method are shown in Section 3. Our conclusions aredescribed in Section 4.

2. Alignment algorithm
In this section, the details of our alignment algorithm
are explained. We assume that all range images have been
converted to mesh models. The algorithm is applied in thefollowing steps:
1. To compute, for all pairs of partial meshes,(a) to search all correspondence of vertices(b) to evaluate error terms of all correspondence pairs
2. To compute transformation matrices of all pairs for
immunizing all errors
3. To iterate steps 1 and 2 until the termination
condition is satisfied
First, we explain the fast method to search corresponding
points. Then, the details of error evaluation and thecomputation of rigid transformations are described.

2.1. Correspondence search
Our algorithm employs points and planes to evaluate
relative distance as the Chen and Medioni method [5]. The
corresponding pairs are searched along the line of sight(Fig. 1). Here, the line of sight is defined as the opticalaxis of a range sensor. Let us denote one mesh as the basemesh and its corresponding mesh as the target mesh. An
extension of the line of sight, from a vertex of the basemesh, crosses a triangle patch of the target mesh andcreates the intersecting point. In order to eliminate false
correspondences, if the distance between the vertex andthe corresponding point is larger than a certain thresholdvalue, the correspondence is removed. This
correspondence search is computed for all pairs of meshmodels.
Though the threshold distance is given empirically as
lgiven, the smaller value is selected as lth compared with theaverage distance of all corresponding points r̂ .
(if : l given < rˆ)⎧ll th = ⎨ given(2)ˆr(otherwize)
⎩
rˆ =

1
N

N

∑ y −x
i

(3)

i

i

N is the number of vertices included in the base mesh.
Target image
Correspondenceyx
Ray direction ofscene imagenynx

Base image
Range finder

Figure 1. Searching corresponding points

To search correspondences quickly, our method uses
index images. Though the complexity of this process is
O(N,), the same as that of the inverse calibration method,
sensor parameters are not required. Furthermore, thesearching process can be accelerated by graphics
hardware. The details of correspondence search using
index images are described below.

the range images that are measured by sensors that adopt a
method like light sectioning. On the other hand, in thecase of range images taken by sensors with scanningmechanisms using mirrors, the spherical projection isbetter because the angles between sampled points areequal (or near) to each other.

2.1.1 Creation of index images
An index image works as a look-up table to retrieve theindex of corresponding patches. Here, we describe theprocedures for creating an index image as follows:
1. A unique index number is assigned to each triangle
patch of a target mesh.
2. Index numbers are converted to unique colors.
3. Triangle patches of the target mesh are rendered on animage plane with the index colors.
First, a unique integer value is assigned to each triangle
patch. Since the assigned value can be any integer number,
0 to n-1 are assigned sequentially, where n is the numberof triangle patches.
Next, the assigned index values are converted to uniquecolors. If the precision of index values is the same as thatof rendering colors, the index values are converted
directly. Assume that the precision of each color channel
is expressed by q bits. [0 q-1] bits of an index value areassigned to Red, and [q 2q-1] bits to Green, and the next
to Blue, and highest q bits are assigned to Alpha. If the
precision is not the same, the indices have to be converted
carefully.
All triangle patches are rendered onto an image plane
using the index colors (Fig. 2). The pixels in which thetriangle patches are not rendered are filled with an
exceptional color like white. The target mesh is assumed
to be described in its measured coordinate system.
Target range image

w
y

z

h

Index image

x
Figure 2. Rendering of index image

Projection method
Generally, perspective projection is used for rendering
the index images. Perspective projection works well for

Zoom image

Figure 3. Example of index images
View frustum
To obtain a sufficient number of corresponding points,
rendering areas have to be determined properly. Allvertices of a target mesh are projected onto the index
image plane. Then, the rectangular area (umin, vmin, umax,vmax) that involves all projected vertices is obtained. Theview frustum is computed so that all vertices are renderedin this area. Minimum and maximum depths are alsoacquired at the projection process.
Image resolution
The resolution of an index image (Lu, Lv) has to be
determined as the following conditions are fulfilled.
(4)
Lu ≥ 2 × w / Δwmin
Lv ≥ 2 × h / Δhmin(5)

w = u max − u min(6)

h = v max − v min(7)

Variables w and h represent the height and width of therendering area respectively (Fig. 2). Δwmin and Δhmin
represent the minimum height and minimum width of all
triangles projected onto the index image plane.
The image resolution (Lu, Lv) can be roughly
determined so that the conditions described in inequality 4
and inequality 5 are satisfied. Although the parameters (w,h, Δwmin , Δhmin ) are different in each partial mesh, aunique resolution that satisfies the conditions of all meshmodels works well for all index images.
The rendering process is accelerated by using graphics
hardware. The rendering time becomes small enough toignore even if the images are rendered at each iterative
step. A large memory space for storing look-up tables is
not required. The memory space for only one index image
can be shared by all mesh models. Figure 3 shows anexample of index images.

2.1.2 Acquisition of corresponding points
By using the index image, corresponding points are
rapidly searched. The procedure for this process is thereverse of that used to make an index image. Here, we
assume that all vertices of the base mesh are previously
converted to the local coordinate system of the target
mesh. The following steps are applied to all vertices of thebase mesh:
1. A vertex is projected onto the index image plane by
the same projection method as is used for the index
image.
2. A color is obtained from the projected pixel.
3. The obtained color is converted to the index valueof a patch of the target mesh.
4. The vertex is projected onto the corresponding
patch; a corresponding point is acquired.
The procedures are depicted in Fig. 4. A vertex of the baseimage is projected onto the index image plane. Then, acolor is acquired from the projected pixel and is converted
to the index of a corresponding patch. Since a correctindex value may not be obtained because of round-off
errors, the vertex is re-projected onto a corresponding
patch, and the crossing point is checked to see whether it
is inside the patch or not. If the crossing point is inside thepatch, the accurate corresponding point is computed. Until
the correct corresponding point is obtained, steps 2-4 areapplied to 3×3 pixels around the projected pixel. The
computational complexity of this process is also O(N).

errors, we used the average normal vectors for error
evaluation.
The transformation matrices of the base and target meshmodels are computed so that this error measure isminimized. The error evaluation function is rewritten as
r
r r
r r
(3)
ε = RB n ⋅ {( RT y + t T ) − ( RB x + t B )}Here, the rotation matrix and the translation vector ofr r
the base and target mesh are RM , RS , t M , t S respectively. Tormake the function simple, the average normal n isassumed to be rotated by the matrix of base range image R
Ｂ. The distance between the base and the target mesh isexpressed as
r
r 2
r
r
r
(4)
ε = minr ∑ (Ri nik ⋅ {( R j y ijk + t j ) − ( Ri xik + t i )})R ,t

i , j ,kIf it is assumed that the angles of rotation are minute,the rotation matrix R is written as⎛ 1⎜
R = ⎜ c3⎜−c⎝ 2Base image
1. Projection toimage plane

Ray direction(Target image)
2. Pick color
Color＝0x0000ff

3. Convert toindex value
Patch Index＝255

Figure 4. Searching a corresponding patch

2.2 Error metric
The error measure between corresponding points is the
cosine distance between the point and the plane. Let thevertex of the base mesh and the corresponding crossingrr
point in the target mesh be x and y , respectively. The
error measure of a pair k is written as
r r r
(1)
ek = n ⋅ ( y − x )r rr nx + n y(2)
n= r r ,nx + n yrr
where nr x and nr y are the normal vectors of x and ydefined around the vertices respectively. Since normal
vectors tend to be greatly influenced by measurement

1

c1

c2 ⎞
⎟
− c1 ⎟
1 ⎟⎠(5)

The translation vector is expressed asrTt = (t x t y t z )(6)After some algebraic manipulations [8], equation 4 isrewritten as
r
r
2
(7)
ε = minr ∑ Aijk ⋅ δ − s ijk
δPatch index image

− c3

i≠ j

r
r
r
s ijk = n ik ⋅ ( x ik − y ijk )

r
r
r
⎧⎪⎛
⎞⎫⎪
⎞ ⎛
Aijk = ⎨⎜⎜ 0{...0 C ijkT 0{...0 ⎟⎟ + ⎜⎜ 0{...0 − C ijkT 0{...0 ⎟⎟⎬⎪⎩⎝ 6i×1
6 ( l − i −1)×1 ⎠6 ( l − j −1)×1 ⎠ ⎪⎝ 6 j×1⎭rrr⎛n × y ⎞Cijk = ⎜⎜ ik r ijk ⎟⎟⎝ − nik ⎠rr r Tδ = (m0 ...mn−1 )rTmi = (c1i c 2 i c3i t xi t yi t zi ) ,

(8)
T

(9)

(10)
(11)

(12)
r
where the number of mesh models is n. By (7) δ iswritten asr r ⎞rr⎛(13)⎜ ∑ A T ijk Aijk ⎟δ = ∑ A T ijk sijk
⎜
⎟
≠
i
≠
j
i
j
⎝
⎠

2.3 Solving linear equationsr
From Eq. 13, δ is computed as the solution of linear
equations that include n×6 arguments. However,ambiguity remains in the equation. Then, the first meshmodel is assumed not to be moved. That is, as shown in
Fig. 5, the linear equations with ((n-1)×6) × ((n-1)×6)
coefficient matrix are solved. If all mesh models areconnected to the first mesh, the coefficient matrix issymmetric positive definite. Also, the matrix becomeslarger and sparser as the number of mesh models increases
and has 6×6 non-zero patterns as shown in Fig. 5.
n×6(n-1)×6∑Ai ≠ j ,k

T

ijk

Aijk

n×6
(n-1)×6

6×6 Non-Zero Area
Positive Definite
Symmetry Matrix

Figure 5 Characteristics of coefficient matrix
Since the computational complexity of direct solvers is
too high, we applied an iterative solver to this problem.
The
computational
complexityof
Choleskydecomposition that is the most popular direct solver for asymmetric positive definite matrix is O(n3). Then, we
employed the pre-conditioned conjugate gradient method(PCG). Though the complexity of PCG is O(n3), the sameas Cholesky decomposition, the number of iterations canbe drastically reduced by pre-conditioning. We employedincomplete
Choleskydecompositionasthepre-conditioner. Since it is known that the coefficientmatrix has 6×6 non-zero patterns, we implemented theincomplete Cholesky conjugate gradient method (ICCG)
specialized for the matrix pattern.
× ((n-1)×6)
Assume that the matrix M is the ((n-1)×6)
r
coefficient matrix shown in Fig. 5. β is the (n-1)×6rvector that is a part of the right side of Eq. 13. β does
not include the transformations of the first mesh model.
Equation 13 is re-written as follows.
r r(14)
Mα = βrr r
T(15)
α = (m1 ...mn −1 )A matrix C is assumed to be a regular ((n-1)×6)×((n-1)×6)
matrix. Equation 14 is written as follows:rr(16)
C −1 M (C T ) −1 C Tα = C −1 β~To simplify rthe equation, we define the matrix M andthe vector β ′ as follows.
~(17)
M = C −1 M (C T ) −1rr
(18)
β ′ = C −1 βEquation 16 is redefined by using the variables above.
r
~r
(19)
Mα ′ = β ′
r
T r
(20)
C α = α′~If the coefficient matrix M is near the identity matrix,
solving Eq. 19 is drastically accelerated by the conjugate
gradient method. That is, the CTC has to be nearly equal tothe original coefficient matrix M.
(21)
M ≅ C TC

But the computation cost to decompose the matrix is veryhigh. The matrix M is incompletely decomposed by the
Cholesky decomposition. In this process, only non-zeroareas are computed: Another element is filled with zero.
(22)
M = UD
Since the matrix M is sparse, the decomposition process isperformed very quickly. Then, matrix C is given as
follows:(23)
C = UD 1 / 2
Once matrix C has been computed, rigid transformations
are calculated from Eq. 19 and Eq. 20 by the conventional
conjugate gradient method.

3. Experimental results
In this section, the effectiveness of our method isdemonstrated by some experimental results. Two data
sets are used for the experiments. Target objects are theface of Deva in Cambodia (Fig. 6(a)) and the Nara Great
Buddha statue in Japan (Fig. 6(b)). The face of Deva wasmeasured by VIVID900 [19]. The resolution of
VIVID900 was fixed to 640×480, and the view angle
depended on mounted lenses. We used a wide lens forscanning the face of Deva. The Great Buddha statue wasmeasured by Cyrax2400 [20]. The resolution and view
angle of the sensor were flexible: users could change
them arbitrarily. In the scanning of Nara Great Buddha,
we adjusted the parameters according to measurement
environments. Generally, 800×800 was used as themeasurement resolution. The details of these data sets
are shown in Table 1 and Table 2 respectively.

Table 1. Data set 1: The face of Deva
Sensor
VIVID900
Images45
Vertices
Max: 76612, Min: 38190, Ave: 67674
Triangles Max: 150786, Min: 71526, Ave: 130437
Table 2. Data set 2: The Nara Great Buddha
SensorCyrax2500
Images114
Vertices
Max: 155434, Min: 11231, Ave: 81107
Triangles Max: 300920, Min: 18927, Ave: 148281
Vertices that measured outside the objects had been
removed previously. Obtained point crowds had been
converted into triangle mesh models. Since the original
data sets were too large to deal with using one PC, thesizes of the data were reduced to 1 / 4.
Our method is evaluated according to the following
three criteria:
1. Number of corresponding points with respect to theresolution of index images
2. Computation time of matrix operations withrespect to the number of mesh models
3. Computation time of alignment with respect to thewith the image resolution.
1

Number of correspondences(v'c/vc)

number of vertices
The PC used for the experiments had Athlon MP 2400+
processor, 2Gbyte memory, and a GeForce4Ti4600
graphics card.

0.8

0.6

0.4
Minimum(Model1)
Maximum(Model1)
Minimum(Model2)
Maximum(Model2)

0.2

0
0

(b)
Figure 6. Target objects(a : the face of Deva, b : the Nara Great Buddha)

400

600

800

1000

Figure 7. Number of corresponding points withresolution of index images

3.1 Number of corresponding points with theresolution of the index image

1200

1000

Face of Deva
Nara Great Buddha

800

Time(ms)

The number of corresponding points is evaluated withrespect to the resolution of the index image. As describedabove, if enough resolution cannot be assigned to theindex image, all triangle patches are not rendered: all thecorresponding points cannot be acquired. Here, therelation between the resolution of index images and thenumber of corresponding points is verified.
Corresponding points were searched for several pairs ofmesh models by gradually changing the resolution ofindex images. We selected a set of mesh models that haveminimum and maximum number of triangle patches fromeach data set as target meshes. Base meshes werearbitrarily selected. Experimental results are shown in Fig.
7. The vertical axis represents the ratio of the number ofcorresponding points acquired by our method v’c to theground truth vc. The resolutions of index images are
represented in the horizontal axis as the square root of
total pixels.
As shown in Fig. 7, when the resolution of the index
image becomes larger than a certain size (800×800),
almost all of the corresponding points are obtained.
Furthermore, the resolution required for obtaining enoughcorresponding points becomes larger as the number oftriangle patches of target mesh increases. In fact, insteadof the number of triangle patches, the measurement
resolution is concerned with the required index image
resolution. However, in this case, it can be said that the
number of triangle patches has the similar characteristics
with the measurement resolution because the sampled
points are distributed densely and uniformly in bothvertical and horizontal directions.
It is not required to estimate the image resolution foreach mesh model. Even if, due to any problems, forexample the limitation of graphics memory space, a large
enough resolution cannot be assigned to the index image,several corresponding points are acquired in compliance

200

Resolution of index images(pixel1/2)

(a)

600

400

200

0
0

25000

50000

75000 100000 125000 150000

Number of vertices

Figure 8. Time to solve linear equations

3.2 Time to solve the linear equations
The computation time to solve the linear equations isevaluated. As described above, the computation time tosolve the linear system is greatly influenced by the
number of mesh models. Thus, the relation between thecomputation time and the number of mesh models is
evaluated here. Data set 2 is used for this experiment. The
computation time of the matrix operations only issequentially measured by changing the number of meshmodels.
The experimental results are shown in Fig. 8. The
horizontal axis represents the number of mesh models,
and the vertical axis represents the computation time. The
results with usual Cholesky decomposition are also shownin this figure in comparison with our method. The
threshold value of the ICCG was set to 1.0 × 10-6.
In the case that the number of mesh models is lower
than 60, Cholesky decomposition is faster than our
method. On the other hand, the computation time of ourmethod increases at a slow rate and becomes smaller
than that of Cholesky decomposition when the number ofmesh models is higher than 70. Moreover, the differences
between Cholesky decomposition and our method becomelarger as the number of mesh models increases. That is, it
can be said that ICCG is effective for aligning a large
number of mesh models simultaneously.

alignment method to these models. Figure 10 shows thealignment results of data set 1. The total computation timewas 1738 seconds after 20 iterations. Figure 11 shows theresults of data set 2. The number of iterations was 20, thesame as for data set 1. Total computation time was 7832
seconds. The figures show that all mesh models werecorrectly aligned.

3.3. Computation time with number of vertices
The computation time is evaluated with respect to thenumber of vertices included in the mesh models. In thisexperiment, it is proven that the computational complexityof aligning a pair of mesh models is O(N), where N is thenumber of vertices. It can be also said that the complexityof searching corresponding points is O(N) because it canbe assumed that the computation time of another task issmall enough to ignore.
Each mesh model was aligned to itself so that the
number of vertices of base model and target model wereequal to each other. Since the number of corresponding
points also affects the computation time, all mesh models
were not moved. That is, the amount of movement isinfinitely zero; the number of corresponding points isnearly equal to the number of vertices. Index images wererendered at each iterative step. The image resolution was
fixed to 800×800. The computation time was evaluatedaccording to the average time taken for 20 iterations.
Experimental results are shown in Fig. 9. The
horizontal axis represents the number of vertices, and thevertical axis represents the computation time. It is clearthat the computation time is increasing linearly with thenumber of vertices. Moreover, there are no differences
between two data sets though these data were taken bydifferent sensors. That is, the efficiency of our method
does not depend on the sensors used for measurements.

Figure 10. Alignment results (the face of Deva)

1400
1200
Cholesky decomposition
6x6 ICCG

Time (ms)

1000
800
600
400
200
0
0

20

40

60

80

100

Number of range images

Figure 9. Computation time with number of vertices

3.4. Alignment results
The alignment results of these data sets are shown in
Fig. 10 and Fig. 11. Previously, we had aligned all meshmodels one by one. Then we applied the simultaneous

Figure 11. Alignment results (Nara Great Buddha)

4. Conclusion
In this paper, we proposed a fast, simultaneousalignment method for a large number of range images. In
order to accelerate the task of searching corresponding
points, we utilized index images that are rapidly renderedusing graphics hardware and are used as look-up tables.
Instead of sensor parameters, only an approximate
resolution of index images is required for this search. Inorder to accelerate the computation of rigid
transformations, we employed a linearized error function.
Since the computation time to solve the linear system
becomes large as the number of range images increases,
we applied the incomplete Cholesky conjugate gradient(ICCG) method. Experimental results showed theeffectiveness of our method. One of our future works will
be to improve the accuracy of alignment results.

Acknowledgment
This work is supported by Leading Project (Digitalarchive of cultural properties) of the Ministry of
Education, Culture, Sports, Science and Technology of theJapanese Government. The author would like to thank thestaffs of the Todaiji Temple in Nara, Japan. The Bayontemple in Cambodia was digitized with the cooperation of
JSA (Japanese Government Team for Safeguarding
Angkor).

via Simultaneous Registration and Robust Combination ofMultiple Range Images.” International Journal of Shape
Modeling, 3(1&2):71-90, 1997.
[10] R. Bergevin, M. Soucy, H. Gagnon, and D. Laurendeau.
To-wards a general multi-view registration technique. IEEE
Transactions on Pattern Analysis and Machine Intelligence,
18(5):540–547, May 1996.
[11] R. Benjemaa and F. Schmitt. “Fast global registration of 3d
sampled surfaces using a multi-z-buffer technique,” Proc.
Int. Conf. on Recent Advances in 3-D Digital Imaging andModeling, pages 113–120, May 1997.
[12] Z. Zhang, “Iterative point matching for registration of
free-form curves and surfaces,” International Journal ofComputer Vision, 13(2):119–152, 1994.
[13] K. Nishino and K. Ikeuchi, “Robust Simultaneous
Registration of Multiple Range Images,” Proc. Fifth Asian
Conference on Computer Vision, pp454-461, Jan., 2002.
[14] David A. Simon, Martial Hebert, and Takeo Kanade,“Realtime 3-D pose estimation using a high-speed range
sensor,” Proc. IEEE Intl. Conf. Robotics and Automation,
pages 2235-2241,San Diego, California, May 8-13 1994.
[15] M. Greenspan and G. Godin, “A Nearest Neighbor Method
for Efficient ICP,” Proc. International Conference on 3D
Digital Imaging and Modeling (3DIM), pp.161-168, 2001.
[16] R. Sagawa, T. Masuda and K. Ikeuchi, “Effective Nearest
Neighbor Search for Aligning and Merging Range Images,”
Proc. the Fourth International Conference on 3-D Digital
Imaging and Modeling (3DIM-03), pp.79-86, 2003.

References[17] http://www.minoltausa.com/vivid/[1] K. Ikeuchi and Y. Sato, Modeling from Reality, Kluwer
Academic Press, 2001.

[18] http://www.cyra.com.

[2] K. Ikeuchi, A. Nakazawa, K. Hasegawa and T. Oishi, “TheGreat Buddha Project: Modeling Cultural Heritage for VR
Systems through Observation,” IEEE ISMAR’03, Tokyo,
Japan. Nov., 2003.
[3] M. Levoy et al. “The Digital Michelangelo Project,” Proc.
SIGGRAPH 2000, pp. 131-144, 2000.
[4] P. J. Besl and N. D. McKay, “A method for registration of
3-D shapes,” IEEE Transactions on Pattern Analysis andMachine Intelligence, 14(2), 239-256, 1992.
[5] Y. Chen and G. Medioni, “Object modeling by registrationof multiple range images,” Image and Vision Computing,10(3), 145-155, 1992.
[6] G. Blais and M. Levine, “Registering Multiview Range
Data to Create 3D Computer Objects,” IEEE Transactions
on Pattern Analysis and Machine Intelligence, Vol. 17, No.
8, 1995.
[7] S. Rusinkiewicz, O. Hall-Holt and M. Levoy, “Real-Time
3D Model Acquisition,” ACM Transactions on Graphics.
21(3): 438-446, July 2002.
[8] T. Masuda, K. Sakaue and N. Yokoya, “Registration andIntegration of Multiple Range Images for 3-D Model
Construction,” Proc. Computer Society Conference on
Computer Vision and Pattern Recognition, 1996.
[9] P. J. Neugebauer. “Reconstruction of Real-World Objects