Parallel Alignment of a Large Number of Range Images

Takeshi Oishi
Institute of Industrial Science,
The University of Tokyo
Atsushi Nakazawa
Osaka University

Ryusuke Sagawa
Osaka University

Ryo Kurazume
Kyushu University

Katsushi Ikeuchi
Institute of Industrial Science, The University of Tokyo

Abstract

This paper describes a method for parallel alignment ofmultiple range images. There are problems ofcomputational time and memory space in aligning a large
number of range images simultaneously. We developed aparallel method to address the problems. Searching forcorresponding points between two range images istime-consuming and requires considerable memory spacewhen performed independently. However, this process canbe preformed in parallel, with each corresponding pair ofrange images assigned to a node. Because thecomputation time is approximately proportional to thenumber of vertices, by assigning the pairs so that the
number of vertices computed is equal on each node, theload on each node is effectively distributed. In order toreduce the amount of memory required on each node, a
hypergraph that represents the correspondences of range
images is created, and heuristic graph partitioning
algorithms are applied to determine the optimal
assignment of the pairs. Moreover, by rejecting redundantdependencies, it becomes possible to acceleratecomputation time and reduce the amount of memory
required on each node. The method was tested on a
16-processor PC cluster, where it demonstrated highextendibility and improved performance.

1. Introduction
In recent years, the technology of modeling from reality
using laser range sensors has been highly developed [1].
Using this technology has enabled the digital preservationof precious cultural heritage objects around the world [2,3]. Laser range sensors can measure 3D surface geometry
with accuracy within a millimeter for a short distance. Fora long distance (50m to 100m), a laser range sensor thatapplies a time of flight can obtain surface geometry with a
1cm degree of accuracy.
While the performance of laser range sensors has been
greatly improved, there are still some problems in creatinga 3D model from a large number of range images.
Because a laser range sensor can measure only the visible
surface, it is necessary to take range images from manydifferent directions. Once the scanning has been
completed, all range images have to be aligned into acommon coordinate system. If an object is small enoughto 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 [8] and are adaptedfrom the method proposed by Chen [9]. 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 [10].
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 [12]. (Thisis the fundamental algorithm used in this study.)
Benjamaa et al. extended the method proposed byBergevin et al. [13] and implemented a simultaneousalignment method while they accelerated the pair-wisealignment algorithm by using multi z-buffers [14].
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 [15] that use kd-trees and that narrow
the search range by using data cache [16], as well as theparallel ICP algorithm proposed by Langis et al., which is
implemented on a PC cluster [17].
Despite the many alignment algorithms, it is difficult toalign the large number of range images that our activities
involve. The computation time for the pairs of range
images increases, and it is necessary to read all range
images into memory when such algorithms are used
(Parallel ICP [17] does not consider the amount ofmemory used). It is thought, moreover, that the amount ofdata will increase along with the development of
measurement technology. Therefore, we need a method inwhich the calculation time is short, the amount of memory
used is small, and the extendibility is high.
Thus, we propose a parallel simultaneous alignment
method that is implemented on a PC cluster because thismethod is cheap and highly extendible. In Section 2, thefundamental alignment algorithm is described. In Section
3, we present the algorithm of parallel computation.
Sections 4 and 5 contain the evaluations of this algorithm
and the alignment results of a large number of range
images, respectively. Our conclusions are presented in
Section 6.

2. Simultaneous Alignment Algorithm
In this section, the outline of the fundamental alignment
algorithm is explained. We assume that all range images
have been converted to mesh models. The algorithm isapplied in the following 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
Our algorithm employs points and planes to evaluate
relative distance as the Chen and Medioni method [9]. The
corresponding pairs are searched along the line of sight.
Here, the line of sight is defined as the optical axis of arange sensor. Let us denote one mesh as the base meshand 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 every pairs ofmesh models.
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 between the pairs is written as(1)n ⋅ ( y − x)rrwhere n is the normal of x defined around thevertex.
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
(2)
ε = RB n ⋅ {( RT y + t T ) − ( RB x + t B )}Here, the rotation matrix and the translation
vector of
r r
the base and target mesh are RM , RS , t M , t S respectively.
The distance between the base and the target mesh isexpressed as
r
r 2
r
r
r
(3)
ε = minr ∑ (Ri nik ⋅ {( R j y ijk + t j ) − ( Ri x ik + 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
⎝ 2

− c3
1
c1

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

The translation vector is expressed asrT(5)t = (t x t y t z )After some algebraic manipulations [12], (3) is rewritten
as
r
r
2
(6)
ε = minr∑ Aijk ⋅ δ − sijkδi ≠ j ,krrr
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 ⎠rrr
δ = (m 0 ...m n−1 )TrTm i = (c1i c 2 i c 3i t xi t yi t zi )

(7)
T

(8)
(9)
(10)

(11)
r
where the number of mesh models is n. By (6) δ iswritten as

r r ⎞rr⎛⎜ ∑ A T ijk Aijk ⎟δ = ∑ A T ijk sijk⎜⎟
i , j ,k
⎝ i , j ,k⎠meshes are accurately estimated.

(12)
2.

The angle θ between ray directions of two meshmodels is less than a threshold value.
Two observation directions of the meshes are
relatively near. This condition also reduces thepossibility of false correspondences between frontand backside meshes, by setting the threshold, as
θ=90º. We could use a more accurately estimatedvalue for this threshold, but since this value is used asa constraint to reduce the possibility described above,
we use this θ=90º for the sake of safety andsimplicity.

3.

The overlapping area of two meshes is larger thana threshold value.
Overlapping area is expressed as the ratio of thenumber of vertices included in one mesh model andthe number of corresponding points between two
meshes. Corresponding points are searched for a few
vertices selected randomly. We used 10% of thevertices for this search. A pair whose overlapping area
is less than threshold value will be removed as weakdata dependency. We set the threshold value as 0.03
to 0.05. Since the computation of overlapping areas
can be performed independently and sequentially foreach pair, the computations are performed easily inparallel without the problem of memory usage.

4.

Two range images are adjacent to each other.
This condition removes non-adjacent relations
sequentially. For example, as shown in Figure 2, if
the length from I0 to I3 is larger than the length from
I1 to I3 (l01 < l03), the arc between I0 and I3 is removed.
Here, the distance is evaluated from the center of amesh model.

3. Parallel alignment based on a PC cluster
Among the simultaneous alignment operations
described in Chapter 2, 1(a) correspondence search and1(b) error evaluation require a large amount ofcomputational time. They also require data space to readin data of all vertices. On the other hand, these two
operations can be conducted independently in each pair ofpartial mesh models. Computation of transformation instep 2 does not require much computational time or
memory space. Thus, we designed correspondence searchand error evaluation in step 1 to be conducted in slave PCs
in a PC cluster, and computation of transformation in step
2 to be conducted in a master PC.

3.1. Graph simplification
We remove redundant or weak data dependency
relations of partial mesh models for the sake of efficiency
in parallel computation. Figure 1 shows overlapping
data-dependency relations. Each node in the graphrepresents one mesh model, and each arc represents anoverlapping dependency relation among mesh models.
The left graph shows the original state in which all themesh models overlap each other. If we conduct alignmentof one mesh as is, we would have to read into a PC’s
memory all the remaining mesh models. By removingsome of redundant overlapping dependencies, we cantransform the original graph into a simpler one as shownin the right figure. By using this simpler relational graph,
we only need adjacent data with respect to a vertex foralignment of a vertex, and we can reduce the necessary
memory space.
Range image
Pair

l01 < l02 < l03
Slave2

I3

I1

The bounding-boxes of two mesh models overlap
each other.
A sufficient overlapped region exists between twomesh models, provided that initial positions of two

I0

l02
l13

I1

if (l13 < l03 )

l01

I0

pair ( I 0 , I 3 )Figure 2. Non-adjacency relation

Figure 1. Data dependency relations

1.

l01

I3

rejectSlave1

We will remove the dependency relation between thetwo mesh models if any of the mesh pairs does not satisfy
any one of the following four conditions:

I2

l02

l03
l13

Slave0I2

3.2.

Parallelization

by

graph

partitioning

algorithms
The problem of load balancing with a minimum amountof required memory is an NP-hard problem [18]. It isdifficult to obtain an optimal solution in a reasonable time.
Alternatively, we employ an approximation method to
solve
this
problembyapplying
heuristic
graph-partitioning algorithms.

3.2.1 Pair-node hyper-graphFirst, we define the pair-node hyper-graph. The left
image of Figure 3 shows a graph that expresses therelations of partial meshes In. The graph is converted tothe hyper-graph in which each node expresses pairs Pi,j of
two partial meshes i and j, and networks represent meshes,
as shown in the right figure of Figure 3. We refer to it as a“pair-node hyper-graph.”
The weight of the network Wneti is defined as thenumber of vertices vi in the partial mesh, i; the weight
Wnodei,j of the node is defined as the sum of the number ofvertices vi and vj.
(13)
W net i = vi
(14)

W node i , j = vi + v j

A pair-node hyper-graph is partitioned so that the sum of
the node weights in each subset is roughly equal forcomputational load balance, and summation of all the
net-weight in each subset is minimized for efficiency ofmemory usage.

P2,3

I2

I3

I0

P1,3

P0, 2

P0,3

I1

P0,1

Figure 3. Pair node hyper-graph
It is necessary to consider both node weights and net
weights in optimization, even though they are related toeach other, and using them seems to be redundant.
Reducing the computational load requires each sub-groupto have equal values in the node-weights. On the other
hand, even when a hyper-graph is portioned equally interms of node-weight, depending on the method, each
sub-group has different memory usage. Let us consider
the example, shown in Figure 3, to divide the hyper-graphinto two sub-graphs. For the sake of simplicity, we assume
that all node-weights and net-weights are the same in all
the nodes and all the networks. When the hyper-graph isdivided into two groups, {P0,2, P1,3, P2,3} and {P0,1, P0,3},
the node balance is achieved in two sub-graphs. The first
sub-graph needs to load in all the data {I0, I1, I2, I3}. The
maximum value in sums of net-weights is four units.
When the hyper-graph is divided into two groups, {P0,2,
P0,3, P2,3} and {P1,3, P0,1}, each sub-group needs only toload in three data sets. The maximum value in the sum of
net-weights is three units. In these two cases, both
portioning methods have roughly equal load balance interms of node-weights, but have different memory usage.
When we divide the graph by considering only memory
usage, it is not guaranteed that each sub-graph has equalload balance. Thus, we will consider both node-weights
and net-weights in the optimization procedure.

3.2.2 Initial partitioning
The pair-node hyper-graph is initially partitioned so thatthe sum of the node-weights in each subset is roughly
equal. Spectral bisection methods [19, 20] that minimize
the edge-cut by using second eigenvector are widely
available, but it is difficult to apply the method to our
problem. Intelligent graph growth algorithm [21] canobtain a fairly optimal solution in a small computation
time. However, this method tends to be trapped in a poor
partitioning [22]. We used the random seeded breadth firstsearch method for initial partitioning. Since the sum of
net-weight included in each subset is greatly influenced by
the selection of the seed, we created initial partitions for
multiple seeds and adopted the partition in which the sumof net-weight included is minimized. In order to obtaink-way partitions, the recursive bisection method is used.
After logk phases, the hyper-graph is partitioned into k
sub-graphs [23].
3.2.3 Refinement of the partition
The partitioned graphs are refined so that the sum of
net-weights included in each subset graph is minimized.
We improved the KLFM algorithm, which is an iterativerefinement algorithm [24, 25]. The algorithm moves anode from one partition to another so that the operation
causes the greatest improvement in the cut-size. While theoriginal KLFM algorithm moves a node at one iteration,
our method moves a net at one iteration. That is, all nodes
connected to the net are moved at the same time. For
k-way refinement, the subset graph of which the sum of
net-weight is maximum weight is computed with all othersubsets. The refinement process is reiterated until there is
no more improvement.
The net gain is computed for all nets along theboundary of two subset graphs. Now, we consider the kth
net at the boundary between the subset graphs, Gi and Gj.
In the case the net N(i,j),k is moved to Gi, the gain gi,j,k isexpressed using two values, Dinti,j,k, the variation of thesum of net weight of Gi and Dexti,j k, the variation of thesum of net-weight of Gj as(15)
g i , j ,k = D ext i , j ,k − D int i , j ,k .
On the other hand, in the case where N(i,j),k is moved toGj, the gain gj,i,k is expressed in the similar way as(16)
g j ,i ,k = D ext j ,i ,k − D int j ,i ,k .
The two lists, Li, Lj, consisting of all gains of the allnets at the boundary, are created. The list with the larger
sum of the total node-weight (computational time) is
selected for consideration of the movement, and thecomponents, candidate nets in the list, are processed one
by one in descending order of the gain. At each movementof one net, all nets and nodes concerned with the net are
updated, and the moved net is locked in order to avoid
thrashing. The sum of the net-weight (memory usage) andthe moved net’s ID are also recorded at each movement.
After all nets are moved, the minimum value of the sum of
the net-weight (memory usage) is compared with thevalue at the starting stage. If the minimum value is smaller
than that of the starting state, the corresponding
movement-sequence is performed, and the next iteration
begins. If not, the refinement process is terminated. See
Figure 4 for the flow chart of the refinement process.

Start

ATijksijk in (12) independently, and sends the matrices tothe master program. The master program computes thetransformation matrices of all range images from thematrices ATijkAijk and ATijksijk received from the slave
programs. The results are applied to all master/slave data.
Each iteration process is continued until the error falls
below a certain threshold value.

Create two lists (Li, Lj) for gains (gi,j,k, gj,i,k)Sort these lists in descending order by gains
Select one list by following condition
if(ΣWnodei > ΣWnodej ) : select Lielse : select Lj
Move an unfixed net which contains minimum gain
Update the neighbors’ gains
Fix the moved net
No
Yes

Terminated?
Yes
Updated?

Consider Graphs Gi and Gj

No

Finish

4. Performance evaluation
This method was implemented on a PC cluster thatconsisted of 8 PCs. Each PC had dual AthlonMP2400+
processors and 4Gbytes of memory, and was connected by
100Base-TX ethernet. The range images used forevaluation were 50 images created artificially from thecomplete 3D model of the Great Buddha of Kamakura.
Figure 5 shows the original 3D model of the statue and thepartial mesh models created artificially. These meshmodels contain an average of 83,288 vertices and 158,376
patches.
In this section, our method is evaluated from theviewpoints of convergence and accuracy, computation
time, and memory usage.

Figure 4. Flowchart of refinement process

3.3. Implementation
We implemented our method as a master/slave system.
The procedures of the computation is are as follows
Algorithm Procedure of Parallel Alignment/* Check correspondence of all pairs of the partial meshes*/
Create-Pair-Table();/* Create the lists of the files for each processor */Create-File-Lists():
while(error > threshold){/* Slave Process*/for(i = 0; i < nmeshes; ++i)
for(j = 0; j < nmeshes; ++j)
Whether-i-and-j-overlap-each-other?{
Correspondence-Search(i, j);Calculation-Each-Matrix(i, j);}/* Master Process */CalculationMatrix(all);/* Master & Slave process */UpdatePosition();}The master program holds bounding-boxes andtransformation matrices from initial position to current
position of all partial meshes, checks all pairs, and creates
the list of computations for each node. The pairs list foreach slave is computed at the beginning of the entire
iteration process based on the relational table using thealgorithm described above. The slave programs receive
the lists and read the required partial meshes into memory.
Then, each slave computes the matrices ATijkAijk and(a) Original model (b) Created mesh models
Figure 5. Partial mesh models for evaluations

4.1. Convergence and accuracy
Because our method rejects redundant dependencies,
the influence of the rejection on convergence andaccuracy has to be evaluated. In this case, the number ofall pairs is 2,450, but it is reduced to 160 by the rejection
process. We needed to verify whether accurate
convergence is performed even when the number of pairs
becomes very small. Virtually created mesh models haveaccurate positions of measured points, so convergence andaccuracy can be evaluated by the distance between an
accurately aligned mesh model and the target mesh model.
The distance between two meshes is defined as an averageof the Euclidean length of all vertices. Each mesh model
added Gaussian noise along the line of sight at maximum
length 10mm. All mesh models were moved at random inthe maximum length of 100mm in the directions of x, y,
and z, respectively, and rotated at random in the maximum
angle of 0.05 radians to the x-axis, y-axis, and z-axis,
respectively.

processors. Computation time Tn is expressed as the ratio
to the computation time with one processor T0. This figure
shows that the computation time is linearly improved as
the number of processors increases. Moreover, our method
improves computation time in a predictable way unlike
the sequential method in which the mesh models are
assigned in arbitrary selected order. The actual
computation time with one processor averages 20560ms,
and the computation time with 16 processors averages
1784ms. Thus, the computation time with 16 processors isapproximately 11.5 times faster than that with 1 processor.

Average of error [m]1
Threshold distance : 0.1m
Threshold distance : 0.5m
Threshold distance : 1.0m
Threshold distance : 2.0m
Threshold distance : 5.0m

0.8

0.6

0.4

0.2

12
0
5

10

15

20

Number of iterations [times]

Figure 6. Convergence with original method

Threshold distance : 0.1m
Threshold distance : 0.5m
Threshold distance : 1.0m
Threshold distance : 2.0m
Threshold distance : 5.0m

0.8

0.6

Sequential method

8
6
4
2
0
0

0.4

2

4

6

8

10

12

14

Number of processors n

16

Figure 8. Computational Efficiency

0.2

4.3 Amount of required memory

0
0

5

10

15

20

Number of iterations [times]

Figure 7. Convergence with our method
The results of the original method and our method areshown in Figures 6 and 7, respectively. The thresholddistance for rejecting outliers while searching forcorrespondences is changed gradually. Although both theoriginal method and our method do not converge at thecorrect positions when the threshold distance is 5m, ourmethod converges at a better position than the position of
the original method. Although the convergence speed ofour method is slower than that of the original method, ourmethod tends to converge at a better position than theposition of the original method. It is thought that this isbecause the rejection of redundant pairs reduces false
correspondence of mesh models. A feature of thealignment algorithm that we used is that it tends to be
influenced by false correspondence and noise. Therefore,
by rejecting redundant pairs, transformations areaccurately estimated. When the threshold values are 0.1m
and 0.5m, the error converges at approximately 0. So we
see that accurate estimation is acquired by our method.

4.2. Computational efficiency
Here, the computation time is evaluated. Computation
time is defined as the time taken for one iteration, and an
average of time of all iterations is used for the evaluation.
Figure 8 shows the time ratio with the number ofNext, the evaluation of memory performance is shown.
The amount of memory usage is shown in Figure 9 with
the number of processors. Each value shows the ratio withthe amount of memory used with a single processor. It
appears that the amount of required memory decreases as
the number of processors increases. Compared with thesequential method, the performance is highly improved byour method. An actual maximum size of required memory
with a single processor is 269Mbytes and that with 16
processors is 48Mbytes. Therefore, our method could
reduce the amount of memory used by approximately 17%
for these mesh models.
7
Graph partitioning method

6

Required memory M1 /Mn

Average of error [m]1

Graph partitioning method10

Calculation time T1 /Tn

0Sequential method

5
4
3
2
1
0
0

2

4

6

8

10

12

14

Number of processors n

Figure 9. Required memory16

Figure 10. Alignment result (Nara Buddha)

5. Experimental results
In this section, we will show the result of parallelalignment of a large number of range images that could
not be aligned by one PC because of limitation of memory
space. We used the following two sets of partial meshmodels.
Model-1. 114 mesh models that measured the Great
Buddha of Nara by Cyrax2400 [18]. These models
contain an average of 327,470 vertices and 606,072
meshes.
Model-2. 210 mesh models that measured the Bayon
Temple in Cambodia by Cyrax2500. These models
contain an average of 433,785 vertices and
798,890 meshes.
Due to the limitation of memory space, the minimum
numbers of processors required for aligning these data sets
were 2 and 4 for Model-1 and Model-2, respectively.Thealignment results computed by the minimum andmaximum numbers of processors are shown in Table 1
and Table 2. These tables show the average computation
time, the maximum amount of memory usage, and theminimum amount of memory usage.
Table 1. Total performance (Nara)
Processors
2
16

Ave. Time(s)76
13.2

Max. Mem(MB)
1287
292

Min. Mem(MB)
1275
254

Table 2. Total performance (Bayon)
Processors
4
16

Ave. Time(s)
103.9
40.2

Max. Mem(MB)
1608
559

Min. Mem(MB)
1456
472

Figure 11. Alignment result (Bayon)

In the case of Model-1, the computation time with 16
processors is 5.75 times faster than that with 2 processors,
and the amount of required memory is reduced 22.6%. For
Model-2, the computation time with 16 processors is 2.58
times faster than that with 4 processors, and required
memory is reduced 34.8%. As for the size of required
memory, these results show an improvement better thanthat described in the previous section (30% for 2-16 and
47 % for 4-16). On the other hand, in the case of Model-2,
although the number of processors is increased 4 times,
the reason the computation time is not greatly improved(2.58 times) is that the time taken for calculation of thetransformation matrix, which is not parallelized and is
performed on the server program, is lengthened. An actual
computation time taken by the server program is anaverage of 14 seconds, and is 35% of the total time taken
for one iteration.
Figures 10 and 11 show the alignment results of theGreat Buddha of Nara and the Bayon Temple in Cambodia,
respectively. Alignment takes approximately 5 minutes for
20 iterations for Model-1 and approximately 15 minutes
for Model-2.

6. Conclusion
In this paper, we have proposed the parallel method for
simultaneous alignment of multiple range images. In
considering time performance and memory performance,
we parallelized the alignment algorithm. Then, we
implemented this method on a PC cluster, and showed its
validity by aligning a large number of range images
simultaneously. Future work will deal with accelerating
the computation of transformation matrices.

Acknowledgments
This work is supported by Ikeuchi CREST (Core Research
for Evolutional Science and Technology) of JST (Japan

Science and Technology Corporation). The authors would
like to thank the staffs of the Todaiji Temple in Nara,
Japan, and the Koutokuin Temple in Kamakura, Japan.
The Bayon Temple in Cambodia was digitized with thecooperation of the Japanese Government Team for
Safeguarding Angkor (JSA).

References[1] K. Ikeuchi and Y. Sato, Modeling from Reality, Kluwer
Academic Press, 2001.
[2] K. Ikeuchi, “Modeling from Reality,” Proc. Third
International Conference on 3D Digital Imaging andModeling (3DIM), Quebec City, Canada. May 28- June 1,
2001[3] M. Levoy et al. “The Digital Michelangelo Project, “Proc.
SIGGRAPH 2000, pp.131-144.
[4] Greg Turk and Marc Levoy. Zippered polygon meshes fromrange images. Proc. ACM SIGGRAPH '94, pp. 311-318,
1994.
[5] Curless, B., Levoy,M., A Volumetric Method for Building
Complex Models from Range Images, Proc. ACM
SIGGRAPH '96, pp. 303-312, 1996.
[6] M. Wheeler, Y. Sato, and K. Ikeuchi, Consensus surfaces
for modeling 3D objects from multiple range images, Proc.
IEEE International Conference on Computer Vision,
pp.917-923, 1997.
[7] R. Sagawa, K. Nishino, K. Ikeuchi, ”Robust and Adaptive
Integration of Multiple Range Images with Photometric
Attributes”, Proc. IEEE Computer Society Conference on
Computer Vision and Pattern Recognition, vol.2,
pp.172-179, 2001.
[8] P. J. Besl and N. D. McKay, “A method for registration of
3-D shapes,” IEEE Transactions on Pattern Analysis andMachine Intelligence, 14(2) 1992, 239-256.
[9] Y. Chen and G. Medioni, “Object modeling by registrationof multiple range images,” Image and Vision Computing
10(3), 1992, 145-155.
[10] Blais, G. and Levine, M. “Registering Multiview Range
Data to Create 3D Computer Objects,” IEEE Transactions
on Pattern Analysis and Machine Intelligence, Vol. 17, No.
8, 1995.
[11] 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.
[12] P. J. Neugebauer. “Reconstruction of Real-World Objects
via Simultaneous Registration and Robust Combination ofMultiple Range Images.” International Journal of Shape
Modeling, 3(1&2):71-90, 1997.
[13] R. Bergevin, M. Soucy, H. Gagnon, and D. Laurendeau.
“Towards a general multi-view registration technique,”IEEE Transactions on Pattern Analysis and MachineIntelligence, 18(5):540–547, May 1996.
[14] 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.
[15] Z. Zhang. Iterative point matching for registration of
free-form curves and surfaces. International Journal ofComputer Vision, 13(2):119–152, 1994.
[16] David A. Simon, Martial Hebert and Takeo Kanade.
Realtime 3-D pose estimation using a high-speed range
sensor. In IEEE Intl. Conf. Robotics and Automation, pages
2235-2241,San Diego, California, May 8-13 1994.
[17] C. Langis , M. Greenspan and G. Godin “The parallel
iterative closest point algorithm,” In Proc. International
Conference on 3D Digital Imaging and Modeling (3DIM),
Quebec City, Quebec. May 28- June 1, 2001.
[18] M. Garey and D. S. Johnson, Computers and Intractability:
A Guide to the Theory of NP-Completeness, San Francisco,
CA: Freeman, 1979.
[19] Alex Pothen, Horst D. Simon, and Kang-Pu Liou.
“Partitioning sparse matrices with eigenvectors of graphs,”
SIAM Journal of Matrix Analysis and Applications,11(3):430–452, 1990.
[20] L. Hagen, A. B. Kahng, “New Spectral Methods for Ratio
Cut Partitioning and Clustering.” IEEE Transactions on
Computer-Aided Design, Vol. 11, No. 9, pp. 1074-1085,
September, 1992.
[21] G. Karypis and V. Kumar. “A fast and high quality
multilevel scheme for partitioning irregular graphs,” SIAM
Journal on Scientific Computing, 20(1):359-392, 1998.
[22] S. Hauck and G. Borriello, "An Evaluation of Bipartitioning
Techniques", IEEE Transactions on ComputerAided Designof Integrated Circuits and Systems, Vol. 16, No. 8, pp.
849-866, August 1997.
[23] A. George and J. W.-H. Liu. Computer Solution of Large
Sparse Positive Definite Systems. Prentice-Hall,Englewood Cliffs, NJ, 1981.
[24] B. W. Kernighan and S. Lin, "An Efficient Heuristic
Procedure for Partitioning of Electrical Circuits", Bell
Systems Technical Journal, Vol. 49, No. 2, pp. 291- 307,1970.
[25] C. M. Fiduccia and R. M. Mattheyses, "A Linear-Time
Heuristic for Improved Network Partitions", Proc. Design
Automation Conference, pp. 241-247, 1982.
[26] B. Krishnamurthy, "An Improved Min-Cut Algorithm for
Partitioning VLSI Networks", IEEE Transactions on
Computers, Vol. C-33, No. 5, pp. 438-446, 1984.
[27] http://www.cyra.com.