International Journal of Computer Vision, 75(1), 189–208, 2007

The Great Buddha Project: Digitally Archiving, Restoring, and Analyzing
Cultural Heritage Objects
KATSUSHI IKEUCHI, TAKESHI OISHI, JUN TAKAMATSU, RYUSUKE SAGAWA, ATSUSHI NAKAZAWA,
RYO KURAZUME, KO NISHINO, MAWO KAMAKURA AND YASUHIDE OKAMOTO
University of Tokyo, Japan

Abstruct.
This paper presents an overview of our research project on digital preservation of cultural heritage objects anddigital restoration of the original appearance of these objects. As an example of these objects, this project focuses on thepreservation and restoration of the Great Buddhas. These are relatively large objects existing outdoors and providing various
technical challenges. Geometric models of the great Buddhas are digitally achieved through a pipeline, consisting of acquiring
data, aligning multiple range images, and merging these images. We have developed two alignment algorithms: a rapid
simultaneous algorithm, based on graphics hardware, for quick data checking on site, and a parallel alignment algorithm, basedon a PC cluster, for precise adjustment at the university. We have also designed a parallel voxel-based merging algorithm forconnecting all aligned range images. On the geometric models created, we aligned texture images acquired from color cameras.
We also developed two texture mapping methods. In an attempt to restore the original appearance of historical objects, we havesynthesized several buildings and statues using scanned data and a literature survey with advice from experts.
Keywords: cultural heritage, range data, alignment, merging, texturing

1. Introduction
Currently, a large number of cultural heritage objects
around the world are deteriorating or being destroyed
because of natural weathering, disasters, and civil wars.
Among them, Japanese cultural heritage objects, inparticular, are vulnerable to fires and other natural
disasters because most of them were constructed of woodand paper.
One of the best ways to prevent these objects from loss
and deterioration is to digitally preserve them. Digital dataof heritage objects can be obtained by using modern
computer vision techniques. Once these data have beenacquired, they can be preserved permanently, and thensafely passed down to future generations. In addition, such
digital data can be used for many applications that aim to
restore real objects through digital simulation and to plan
restoration projects through precise measures given bysuch digital models. And by creating multi-media content
from digital data, a user can view digital contents throughthe Internet from anywhere in the world, without moving
the objects or visiting the sites.
One of the origins of this line of research is Kanade’s
virtualized reality project [1]. The basic idea of thisproject was to create 3D virtual reality models throughobservation of the real objects. Kanade and his students
constructed an experimental room equipped with multiple
TV cameras, which were focused in the center of the room.
By observing some event at the center of the room, suchas a player playing basketball, they obtained sequences of
2D images captured through synchronized TV cameras.
By applying a multi-baseline stereo algorithm or its
variation to these image sequences, they created a seriesof 3D models of an object, a moving 3D model. By
combining this moving 3D model with a background 3Dmodel and rendering the resulting synthesized 3D model,
they generated a series of 2D images from many

International Journal of Computer Vision, 75(1), 189–208, 2007

Figure 1.

Three steps in geometric modeling.

directions. This line of research has been extensively
explored [2-5] and a couple of representative results areincluded in this special issue.
Another line of research has been directed towardmodeling of various cultural heritage objects in a very
precise manner, using laser range sensors [6]. Recently,
technologies of laser sensing have drastically advanced,
and laser range sensors provide very accurate 3D range
images of an object. However, those range images, given
by a range sensor, are partial mesh models of an entire
object, obtained from arbitrary directions. Thus, theresearch issues include how to align, or to determine
relative relations, among such partial meshes, and how tomerge these aligned partial mesh models into a unified
mesh model of the object. Representative examples
include Stanford University’s Michelangelo Project [7],IBM’s Pieta Project [8], and Columbia University’s
French cathedral project [9], to name a few.
We have been working to develop digital archival
methods by using laser range sensors [6]. Our project has
a number of unique features; among them is its focus on
digitizing large outdoor objects such as the Kamakura
great Buddha and Cambodia’s Bayon temple. Theselarge-scale objects present several challenges inprocessing range data into a unified mesh model. Our
project emphasizes not only geometric modeling but also
photometric and environmental modeling, particularly foroutdoor objects.
The remainder of this paper is organized as follows.
Section 2 describes the outline of the geometric pipeline
developed, a rapid alignment algorithm based on graphics
hardware, a parallel simultaneous alignment algorithm
based on a PC cluster, and a parallel voxel-based merging
algorithm, solving the issues given under digitizing large
objects. Section 3 describes methods to align observed
textures from a digital camera with range data for

texturing large outdoor objects. Section 4 reports ourefforts to restore the original appearance of these objects
using acquired digital data and a literature survey. Section
5 summarizes this paper.
2. Geometric Modeling
2.1.

Overview

Figure 1 shows an overview of three steps in geometricmodeling: data acquisition, alignment, and merging.
Several computer vision techniques, such as traditional
shape-from-X and binocular stereo, or modern range
sensors, provide 3D data points. In this paper, we mainly
utilize laser range sensors as the input device of the 3D
data points. These 3D data are represented as a set of 3D
data points with connecting triangular arcs. This data
representation is referred to as a (triangular) mesh modelof an object. By repeating the data acquisition process so
as to cover the entire object’s surface, we have a set ofpartial mesh models, overlapping each other and covering
the entire object surface as the combination of thosepartial mesh models.
The second step in geometric modeling is to align thesepartial mesh models. Since each sensor is located at anarbitrary position on data acquisition, we have todetermine relative relations of these partial mesh models,referred to as alignment, by considering resemblances inthe data set. When we handle a large object, even if we
can assume that the sensor itself maintains the same
degree of data accuracy as one for a smaller object, eachscan covers a smaller portion of the object; many scans
are necessary to cover the entire surface of the object. As
the result, a long alignment sequence is formed. It is
important to avoid accumulation of errors along thisInternational Journal of Computer Vision, 75(1), 189–208, 2007
sequence of alignment.
The third step is to integrate the aligned multiple meshmodels into a complete mesh model, representing anentire surface of an object. This step is referred to as'merging.' The procedure can be considered asdetermining one surface position from multiple
overlapping surface observations. In the merging
procedure, it is important to make the integration
framework robust against any noise that may occur when
scanning the range images or that may be inherited fromthe registration procedure.
2.2.

Alignment

2.2.1. A Rapid Alignment Based on Ggraphics Hardware.
Our rapid algorithm employs points and planes to evaluate
relative distance as the Chen and Medioni and Gagnon etal. methods [13,15]. Scanning a large architectural object
requires using different types of range sensors, due tocomplex scanning conditions at sites, whose resolutions
may be different from each other. An alignment algorithm,
based on point correspondence such as ICP [10], does notwork well on such range data due to inequality in point
resolutions. We can mitigate this situation with analgorithm based on point-face correspondence, because it
creates virtual points on face patches, and sets up
correspondence between the real points and virtual points,
even if there is unbalance in point distributions.
Our algorithm uses an M-estimator, in particular, theLorentzian function, as the error measure, to avoid theeffects from the outliers [18]. It is well known that the L-2
norm is susceptive to noise. Surfaces of cultural heritage
objects are often covered by foreign materials such asmolls or water. Sometimes returned range values are quite
noisy from such surfaces. The M-estimator effectively
removes outliers. Since we prefer a smooth differentiable
function for effective minimization process, we chose theindex image(scene mesh)

Lorentian function as the evaluation function.
In order to increase computational efficiency, ouralgorithm is designed to utilize graphics processing
hardware. The corresponding pairs are searched along theline of sight for the usage of the graphics hardware. Here,the line of sight is defined as the optical axis of a range
sensor. While this does not guarantee that all thecorresponding pairs are correct, using the M-estimator
removes the effect of such mismatching. It is also true thatafter several iterations of the minimization process, theprocess eventually converges into the correct
corresponding pairs. But for the sake of rapid computation,
we decided to use the line of sight searching method.
Let us denote one mesh as the model mesh and its
corresponding mesh as the scene mesh. One vertex in themodel mesh is depicted as the point in Fig. 2. Eachtriangular mesh in the scene model is assigned one
particular color, and then the color map is depicted on anindex plane using the painting capability of graphics
hardware. An extension of the line of sight, from a vertexof the model mesh, crosses a triangular mesh of the scenemesh and creates the intersecting point. This operation is
implemented using the Z-buffer capability of the graphics
hardware. In order to eliminate false correspondences, if
the distance between the vertex and the corresponding
point is larger than a certain threshold value, thecorrespondence is removed. This correspondence search iscomputed for every combination of mesh models.
The error measure between corresponding points is the
cosine distance between the point and the plane. Let thevertex of the model mesh and the corresponding crossingrr
point in the scene mesh be x and y , respectively. The
error measure between the pairs is written as(1)n ⋅ ( y − x)rwhere n is the normal of x defined around thevertex.
The transformation matrices of the model and scene

Vertex of model meshLine of sightcolor Ætriangular index

Color＝0x0000ff
Figure 2.

Search procedure.

Triangle Index＝256
Corresponding point yInternational Journal of Computer Vision, 75(1), 189–208, 2007
meshes are computed so that this error measure isminimized. The error evaluation function is rewritten as
r
r r
r r
(2)
RM n ⋅ {( RS y + tS ) − ( RM x + tM )}Here, the rotation matrix and the translation
vector ofrthe model and scene mesh are RM , Rs , t M , t s respectively.
The distance between the model and the scene mesh isexpressed as
r
r r
r r
2
(3)
ε 2 = min ∑ (RM n ⋅ {( RS y + t S ) − ( RM x + t M )})R ,ts ≠ m,kIf it is assumed that the angles of rotation are small, therotation matrix R can be approximated as⎛ 1⎜
R = ⎜ c3
⎜− c
⎝ 2

− c3
1
c1

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

The translation vector is expressed as
t = (t x t y t z )(5)After some algebraic manipulations [17], Equation (3)is rewritten as

ε 2 = min ∑ Amsk Δ − γ msk
δ

2

(6)

m≠ s,kHere, m,s,k denotes one model mesh, one scene mesh,and one corresponding pair, respectively.
(7)
γ msk = n mk ⋅ ( x mk − ymsk )⎞⎛⎞ ⎛Amsk = ⎜ 0{...0C msk 0{...0 ⎟ + ⎜ 0{...0− C msk 0{...0 ⎟⎜ 6 m×1 { 6 ( N − m −1)×1 ⎟ ⎜ 6 s×1 123 6 ( M − s −1)×1 ⎟
6×1
⎝
⎠ ⎝
6×1
⎠

(8)

⎛ n × ymsk ⎞⎟Cmsk = ⎜ mk
⎟
⎜ −n
mk
⎠
⎝

(9)

Δ = (δ 0 Lδ N −1 )

(10)

(11)
δ m = (c1m c2 m c3 m t xm t ym t zm )where the number of range images is N . By (6) Δ iswritten as
−1

⎛
⎞
(12)
Δ = ⎜⎜ ∑ AT msk Amsk ⎟⎟ ∑ AT msk γ msk⎝ m ≠ s,k
⎠ m≠ s,k
To avoid accumulation of errors, we have developed asimultaneous alignment method. Traditional sequential
methods [10,13,14] such as the Iterative Closest Point(ICP) algorithm align these meshes one by one, andprogressively align a new partial mesh with previouslyaligned meshes. If a few partial meshes can cover theentire surface of an object, the accumulation of alignment
errors is relatively small and can be ignored; sequentialalignment works well for a small object. However, forlarge objects, sometimes, we may need more than ahundred partial mesh models. In such cases, the error
accumulation would be very large, if we employ
sequential alignment. Thus, we align all partial meshes so
as to reduce the errors among all the pairs simultaneously

in equation (6).
We have implemented a simultaneous alignment
algorithm, and verified the effectiveness of the algorithm,
with respect to the pair-wise alignment algorithm. Since
we need a true value of the aligned result, we generate
synthesized 49 data sets by segmenting a whole 3D meshmodel of the Great Buddha of Kamakura with small
overlapping boundaries. We add a Gaussian noise, whose
standard deviation is 3 mm with 1cm cut-off. This level ofnoise is typical in the data obtained by our sensors. These
data were perturbed along a random direction to simulate
initial alignment errors in the manual alignment process.
Figure 3 shows synthesized range data. The simultaneousand the pair-wise alignment programs are applied to these
data. Figure 4 shows the resulting converging process.
While the simultaneous algorithm provides a constanterror along the number of data sets, the pair-wisealgorithm accumulates errors from the true value along thenumber of data sets processed in the horizontal axis.
2.2.2. Parallel alignment based on a PC cluster. The
rapid
simultaneousalignment
algorithm,
described in the previous section, is mainly utilized
for rough alignment of the data acquired daily at
the site. Although the algorithm employs
simultaneous alignment for accuracy, thealgorithm, designed for a single notebook type PC,
without considering the aspect of the memory
capability, may cause memory overflow forlarge-scale data of the entire target. We havedesigned a parallel alignment algorithm based ona PC cluster for large-scale alignment of the entire
data set.
2.2.2.1. Overview of the Parallel Algorithm. The
simultaneous alignment 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 conditionis satisfied
Among these operations, 1(a) correspondence searchand 1(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.

International Journal of Computer Vision, 75(1), 189–208, 2007

Figure 3. Synthesized range data. (a) Original whole mesh model of Kamakura Buddha, (b) Synthesized range data with initial perturbation andGaussian noise to simulate alignment error and simulate sensor errors, respectively.

0.012

Pair-wise alignment
Simultaneous alignment

Average of Errors [m]

0.01

0.008

0.006

0.004

0.002

0
0

5

10

15

20

25

30

35

40

45

Data number
Figure 4. Converged result.

2.2.2.1. Graph simplification. We remove redundant or
weak data dependency relations of partial mesh models
for the sake of efficiency in parallel computation. Figure 5
shows overlapping data-dependency relations. Each nodein the graph represents one mesh model, and each arc
represents an overlapping dependency relation among
mesh models. The left graph shows the original state inwhich all the mesh models overlap each other. If we
conduct alignment of one mesh as is, we would have toread into a PC’s memory all the remaining mesh models.
By removing some of redundant overlapping

dependencies, we can transform the original graph into a
simpler one as shown in the right figure. By using thissimpler relational graph, we only need adjacent data withrespect to a vertex for alignment of a vertex, and we canreduce the necessary memory space.
We will remove the dependency relation between thetwo mesh models if any of the mesh pairs does not satisfy
any one of the following three conditions:International Journal of Computer Vision, 75(1), 189–208, 2007

Slave3l01 < l02 < l03

I2

I3

l02

l03
l13
I1

Slave1

l01

I2
I3

l02

I0

l13
I1

if (l13 < l03 )Slave2

rejectI0l01

pair ( I 0 , I 3 )Figure 5. Data dependency relation.

Figure 6.

1.

The bounding-boxes of two range images overlap
each other. A sufficient overlapped region exists
between two mesh models, provided that initial
positions of two meshes are accurately estimated.

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 the possibility offalse correspondences between front- and backsidemeshes, by setting the threshold, as θ = 90° . We
could use a more accurately estimated value for thisthreshold, but since this value is used as a constraint
to reduce the possibility described above, we use thisθ = 90° for the sake of safety and simplicity.

3.

Two range images are adjacent to each other. Thiscondition removes non-adjacent relations sequentially.
For example, as shown in Fig. 6, if the length fromI 0 to I3 is larger than the length from I1 to I 3(l01 < l03 ), the arc between I 0 and I3 is removed.
Here, the distance is evaluated from the center of amesh model.

A pair-node hyper-graph is partitioned so that the sumof 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.
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 Fig. 7, 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.

2.2.2.3 Parallelization by graph partitioning algorithms.
The problem of load balancing with a minimum amountof required memory is an NP-hard problem. It is difficult
to obtain an optimal solution in a reasonable time.
Alternatively, we employ an approximation method to
solve
this
problembyapplying
heuristic
graph-partitioning algorithms.
Pair-node hyper-graph. First, we define the pair-node
hyper-graph. The left image of Fig. 7 shows a graph that
expresses the relations of partial meshes In. The graph isconverted to the hyper-graph in which each node
expresses pairs Pi,j of two partial meshes i and j, andnetworks represent meshes, as shown in the right figure of
Fig. 7. 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
W node i , j = vi + v j(14)

Non-adjacency relation.

I2

I3

I0
I1
Figure 7.

P2, 3

I3
P1,3
I1

Pair-node hyper-graph.

I2
P0, 2

P0,3
P0,1

I0

International Journal of Computer Vision, 75(1), 189–208, 2007
Initial partitioning: The pair-node hyper-graph is initiallypartitioned so that the sum of the node-weights in eachsubset is roughly equal. We used the random seededbreadth first search 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 inwhich the sum of net-weight included is minimized. Inorder to obtain k-way partitions, the recursive bisection
method is used. After logk phases, the hyper-graph ispartitioned into k sub-graphs.
Refinement of the partition: The partitioned graphs arerefined so that the sum of net-weights included in eachsubset graph is minimized. We improved the KLFM
algorithm, which is an iterative refinement algorithm. The
algorithm moves a node from one partition to another sothat the operation causes the greatest improvement in the
cut-size. While the original KLFM algorithm moves anode 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 graphof which the sum of net-weight is maximum weight iscomputed with all other subsets. The refinement process isreiterated 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 , L j , 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 thevalue 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
Fig. 8 for the flow chart of the refinement process.
Implementation: We implemented our method as amaster/slave system. The procedures of the computation isare as follows

Algorithm Procedure of Parallel Alignment/* Check correspondence of all pairs of the range images*/
Create-Pair-Table;/* Create the lists of the files for each processor */Create-File-Lists:
while(error > threshold){/* Slave Process*/for(i = 0; i < nImage; ++i)
for(j = 0; j < nImage; ++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 range images intomemory.
T
Then,each slave computes the matrices A msk Amsk andTA msk γ msk in (12) independently, and sends the matrices tothe master program. The master program computes thetransformation
matrices of Tall range images from theTmatrices A msk Amsk and A msk γ msk received from theslave programs. The results are applied to all master/slave
data. Each iteration process is continued until the error
falls below a certain threshold value. Equation (12) isevaluated in all the data. The matrix is very sparse; theICCG provides computation 10 times faster than usual
SVD.
2.3. Merging
Once we can determine relative configurations among all
the mesh models through the alignment operations, we
can connect these aligned mesh models into a unified
mesh model representing an entire object’s surface. This
operation is referred to as merging. One of the simplest
methods is to directly connect each mesh, and remove
overlapping meshes. This direct method discards multiple
observations of overlapping regions, and creates a sharp
jump around the connecting boundaries. In order to avoid
these shortcomings, we employ a voxel-based merging
algorithm.

International Journal of Computer Vision, 75(1), 189–208, 2007

Start

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

FinishFigure 8. Flow chart of refinement Process.

2.3.1. Voxel-based Merging Algorithm. After all partial
mesh models have been aligned, a volumetricview-merging algorithm generates a consensus mesh of
the object from them. Our method, first, calculates a
volumetric implicit-surface representation, where eachvoxel in the volumetric space has an estimated signed
distance from the position of the estimated consensus
mesh. Unlike previous techniques based on
implicit-surface representations, our method estimates thesigned distance to the object’s surface by determining a
consensus of locally coherent observations of the surface[20-23]. This consensus method is effective, in particular,
for noisy data, often provided at the occluding boundaries.
We utilize octrees to represent volumetric implicit
surfaces, thereby effectively reducing the computation andmemory requirements of the volumetric representation
without sacrificing accuracy of the resulting surface. Then,
this signed distance representation is converted to asurface mesh by using a variant of the marching-cubes
algorithm [19].
We originally designed software for a single PC.
However, because recent input data is unpredictably huge,
we decided to build a PC cluster to run this merging
software; the cluster parallel-processes the merging
algorithm to save computation time and utilize the large
memory space of many PCs [24]. We produced one
integrated digital of the Great Buddha of Kamakura andanother of the Bayon Temple in Cambodia with thissoftware.
2.3.2. Results of merging process. We have digitallyarchived Japanese Buddhas, including Asuka (7th century),
Kamakura (13th century), and Nara (17th century). Here,the great Buddha of Nara, although the name indicates the
Nara period of the 8th century, was burned, and the current
one was rebuilt in the 17th century. Thus, it is the newest
among the three. It is interesting to note that the faces of

the Buddhas have become wider as the time passed, as
shown in Fig. 9.
2.4. Analysis of geometric modeling
We extended this effort toward foreign cultural heritage
objects and obtained all 173 Buddha faces of the Bayontemple in the Angkor ruin in Cambodia. Eighteen of thesefaces are shown in Fig. 10. According to the JSA(Japanese Government Team for Safeguarding Angkor)
research, those faces can be classified into three
categories: Deva, Devata, and Ashura. See Fig. 11 forexamples. However, some of them are quite difficult toclassify due to parallel work by different craftsmen anddifferent techniques employed. Others have deteriorated
due to weathering. And still others are unfinished.
In order to classify these faces in an accurate manner,
we examined the 173 faces using cluster analysis. First,
we converted all the faces into normalized depth images,
defined a standard face among 173 faces, and adjusted all
the remaining faces onto this standard face using aclassification procedure. We extracted three key points: apair of inner corners of the eyes, and the middle point
between the mouth and the nose. We obtained translation,
rotation, and scaling parameters to minimize the differences
between two sets of key points of the standard and current
faces. We applied those parameters to the range data of thecurrent face, and obtained a 64 by 64 range image of thecurrent face. We conducted two types of analyses:
supervised and unsupervised. Details are given in [25].
2.4.1 Supervised analysis: Linear discrimination function.
The purpose of supervised analysis is to clarify thedifferences among the given classes. JSA has alreadyclassified all faces into three types based on subjective
evaluations by an artist. Through such a supervised
analysis, we can verify correctness of the process, andthen objectively evaluate the differences using statisticalanalysis methods.
In this paper, we use a linear function f (x ) = n ⋅ x + d
as the classification function. The dimension of thesample space, that is, image size ( = 64 × 64 , in this case),
is much greater than the number of samples; there are only
173 faces in the Bayon Temple. It is preferred that thedimension and parameters of the function are small inorder to prevent a so-called “over-fitting” problem.

Table 1. 3D data of great Buddha.
HeightNumberNumber(m)of verticesof meshes

2.7

1.5 Million

2.9 Million

Kamakura(13 Century)

13

5.0 Million

9.8 Million

Nara (17th Century)

15

36.3 Million

69.1 Million

Asuka (7th Century)thInternational Journal of Computer Vision, 75(1), 189–208, 2007

Figure 9.

Three great Buddha: Asuka (7th century), Kamakura (13th century), and Nara(16th century) Buddha.

Figure 10.

Bayon face library: 18 faces out of 153 faces in Bayon temple in Angkor ruin, Cambodia.

Roughly speaking, the parameters of the function n, dcan be determined by maximizing S B ST , where S B andST are intraclass and interclass variances, respectively.
Because the matrix S is not full-rank matrix, we solve
the equation using the singular value decomposition(SVD) method while minimizing n . Figure 10 illustrates
the result of the analysis. The graph is obtained byprojecting all vector points of faces on the 2D flat surface
that is determined by two linear discriminant functions inorder to clearly express the classification result.
The former discriminant function classifies them into
Devata (female god) and Deva (male god) with separate
planes corresponding to the y-axis in Fig. 12. The rightside area is a female area, and the left side area is a malearea. Although Asura data (male devil) are not used fortraining the function, almost all of Asura data are in themale area.
The second function classifies them into Deva (god)and Asura (devil) with separate planes corresponding to
the diagonal line. The upper side area of the diagonal line
is a god area and the lower side is a devil area. Almost allof Devata (female god) are also in the god area, eventhough we did not use Devata data for training thisfunction.
2.4.2. Cluster analysis. The purpose of unsupervised
analysis is to discover new knowledge throughclassification of the faces without any a priori standards.

International Journal of Computer Vision, 75(1), 189–208, 2007

Figure 11.

Ashura, Dava, and Devata. The black parts in the figures indicate areas where data is missing due to occlusion.

30
51S
25

Devata
Deva
Asura

20
35N

15

10

5

50E

0
-25

-20

-15

-10

-5

0

5

10

15

20

-5

-10

Figure 12.

Linear discrimination function

Cluster analysis provides us with some classification ofsamples according to distances among them. We calculate
the distance based on the Ward method, and define thedistance as the Euclidean distance in the distance image
space. We employ agglomerative hierarchical cluster
analysis. This analysis begins with each sample being
considered as one cluster and then proceeds to combine
the nearest two clusters until all samples belong to one
cluster. By conducting this cluster analysis, we observed
spatial resemblance groups among faces, as shown in Fig.
13. We determined that there were four or fiveindependent worker groups working on these faces in aparallel manner.

3. Photometric modeling: Texturing
Surface color distribution is important in representing
the appearance of cultural heritage objects. Some clay
statues still retain the surface colors that were painted at
the time the statue was originally made. For such types ofcultural properties, we have to archive color information
as well as geometric information. For this purpose, wehave developed two kinds of texturing methods:
calibration-based and reflectance edge-based methods.

International Journal of Computer Vision, 75(1), 189–208, 2007

Figure 13.

Similarity group.

3.1. Calibration-based Method
The main issue in texturing is how to determine therelationship between image sensors and geometricalsensors. When short-distance range sensors can be used,
as shown in Fig. 14, the most promising method is tocalibrate the geometrical relationship between the image
sensor and the range sensor before scanning.
Assume that the coordinate system of the image sensor
is (xc,yc) and the corresponding point in the range image is(X,Y,Z); the relationship between them can be describedas:⎡X ⎤⎡hx c ⎤⎢ ⎥⎢hy ⎥ = c ⎢ Y ⎥
34
⎢ c⎥
⎢Z ⎥
⎢⎣ h ⎥⎦
⎢ ⎥
⎣1⎦
⎡k u
= ⎢⎢ 0⎢⎣ 0− k u cot θ
k v / sin θ0(17)
u0 ⎤ ⎡ fv 0 ⎥⎥ ⎢⎢ 0
1 ⎥⎦ ⎢⎣ 0

0
f
0

⎡X ⎤
0 0⎤
t⎤⎢Y ⎥⎡R
0 0⎥⎥ ⎢ 33 ⎥ ⎢ ⎥0 1⎦ ⎢ Z ⎥1 0⎥⎦ ⎣
⎢ ⎥⎣1⎦

The matrix C34 represents the relationship between theimage and the world coordinate, and it can be calculatedby scanning the calibration box. Inversely, when we map
the texture images onto the geometrical triangular mesh

X n = {( xn , y n , z n 1 ≤ n ≤ 3 } , the corresponding points inimage coordinate xc = {xc , y c } can be easily calculated
as:
xc =

c11 X n + c12Yn + c13 Z n + c14 ,c31 X n + c32Yn + c33 Z n + c34

yc =c 21 X n + c22Yn + c 23 Z n + c 24
c31 X n + c32Yn + c33 Z n + c34(18)

For the modeling of the Koumoku-Ten clay figure, weused 60 range images and color images that were taken at
the same time. Figure 15 shows a picture of Komoku-Ten(a), the geometric model (b), and the synthesis result
under a new lighting condition generated using the
texturing result (c).
ColorsensorRange
sensorRelative
relation

Calibration box
Figure 14.

Calibration method.

International Journal of Computer Vision, 75(1), 189–208, 2007

(a)

(b)

(c)

Figure 15. Synthesized Komoku ten. (a) Picture of Komoku Ten, (b) 3D model, (c) Synthesized result under a new lighting condition generatedusing the texturing result.

3.2. Reflectance Edge-based Method
As shown in the previous section, one solution for
determining the relationship between range and color
image is through calibration using a calibration fixture.
However, this method requires that the range and color
sensors be fixed on the fixture once the relationship iscalibrated. Further, the calibration-based method is
accurate only around the position occupied by thecalibration fixture. When a target object is very large, thismethod becomes unreliable due to the lens distortion.
Thus, we need a method that does not rely on calibration
for handling a large object.
Generally speaking, range sensors often provide
reflectance images as side products of range images. The
returned timing provides a depth measurement, while thereturned strength provides a reflectance measurement. A
reflectance image is a collection of the strength ofreturned laser energy at each pixel. This reflectance imageis aligned with the range image because both images are
obtained through the same optical receiving device.
Commonly available range sensors, including ERIM,
Preceptron, and our main sensor, CYRAX, provide thisreflectance image.
We employ this reflectance image as a vehicle for thealignment of range images with color images [9,26].
Reflectance images share characteristics similar to color
images due to the fact that both images are somehow
related with surface roughness, as shown in Fig. 16. Since
our CYRAX range scanner uses a green laser diode,reflectance edges can be observed along materialboundaries between two different reflectance ratios forthis wavelength. Since different materials are of different
colors, a discontinuity also appears in the color images.
Jump edges along small ranges in a range image also
appear as jump edges in a reflectance image as well as in acolor image. Occluding boundaries are observed both inreflectance images and in color images.

International Journal of Computer Vision, 75(1), 189–208, 2007

ReflectanceColor
Material
difference
Jump edge

Occlusion
Figure 16.

Reflectance and color images.

Prior to the alignments, we paste the necessaryreflectance edges onto the 3D geometric model. As
mentioned above, since occluding boundaries varydepending on the viewing direction, edges along theoccluding boundaries are first removed from thereflectance images. Edges along the current occluding
boundaries will be estimated from the 3D geometricmodel and the current viewing direction. Our algorithm
extracts them automatically, and uses them for alignment.
We align edges extracted from reflectance images withthose in color images so that the 3D position error of thoseedges is minimized by iterative calculation as shown in
Fig. 17. Extracted edges are represented as a collection ofpoints along them. The alignment is done between 3D
reflectance points on a 3D geometric model projected on
the image plane and 2D color edge points in the 2D
image.
To establish correspondence, the system finds thecolor image points that are nearest to the projectedreflectance points. This operation is similar to the ICP
operation.
To determine the relative pose that coincides with theposition of 2D color edges and projected 3D reflectanceedges, we use the M-estimator.

3D error zi

First, the distance between corresponding 2D color
edge points and 3D reflectance edge points is evaluated asshown in Fig. 17, where z i is a 3D error vector that is ona perpendicular line from a 3D reflectance edge point tothe stretched line between the optical center and a 2D
color edge point on the image plane.
(19)
ε i = Z i sin θ
where Z i is the distance between the optical centerand a 3D reflectance edge point, and θ is the angle
between the color edge point and the reflectance edge
point.
The system finds the configuration, P, which
minimizes the total error, E, where ρ is an error
function. The minimum of E ( p ) can be obtained by∂ρ (ε i ) ∂ε i
∂E
=∑
=0
∂P
∂ε i ∂Pi

We can consider ω (ε ) as a weight function toevaluate error terms.

ω (ε ) =

1 ∂ρ
ε ∂ε

(21)

By substituting Eq(21) into (20) we obtain∂ε∂E= ∑ ω (ε i )ε i i = 0∂ρ∂P
i

3D edge point

2

Geometric model

Nearest pointθProjection to the image plane

Zi
Figure 17. Texturing algorithm.

Image plane(22)

We choose the Lorentzian function for this function.
⎛ 1⎛ ε ⎞ ⎞ω (ε ) = ⎜1 + ⎜ ⎟ ⎟⎜ 2⎝σ ⎠ ⎟⎝⎠2D image edge

(20)

−1

(23)

By solving this equation using the conjugate gradientmethod, we can obtain the configuration P that minimizes
the error term and gives the relative relationship between
the camera and the range sensor. Figure 18 shows the
texture-mapped Kamakura Buddha. Since this methodminimizes a non-linear equation, we need an initialalignment. The initial alignment is given manually using
our GUI. For the current implementation, relatively
accurate alignment is necessary for rotation, but it is not
the case for translation.

International Journal of Computer Vision, 75(1), 189–208, 2007

Figure 18. Texturing result.

4. Restoring the Original Nara Great Buddha
One of the advantages of obtaining digital data ofcultural heritage objects is to modify those data anddisplay the original appearance of the object. In order todemonstrate this ability, after we obtain the precise
geometric and photometric information about cultural
heritage objects in their current state, we can modify thecurrent data into a hypothesized original state. In thissection, we describe one of the examples: the restorationof the Nara Great Buddha and its main hall.
The Nara Great Buddha is one of the most important
heritage objects in Japan. The Buddha statue is sitting in
the Buddha palace at the Toudaiji temple in Nara, Japan.
Originally, this temple and statue were constructed byorder of the Shomu emperor, in the 8th century. Here, theoriginal one is referred to as “Tempyou Big Buddha.” Theoriginal Tempyou Buddha was made of bronze andcovered with gold plate. Unfortunately, however, thepalace was burned and the statue was melted twice due tocivil war in Japan. The current Buddha and palace wererebuilt in the 17th and 18th centuries. Accordingly, theshape of the current Great Buddha is different from that of
the original one in the 8th century.
4.1. Restoring the Nara Great Buddha
As the first step, we acquired the complete 3D meshmodel of the Nara Great Buddha in its current state byusing the geometrical modeling techniques described in
Section 2. We collected more than 100 partial meshmodels using CYRAX sensors. Those partial mesh models
were aligned using the parallel alignment algorithm on a
PC cluster and merged into a unified mesh model with
70M polygon. Figure 19 shows the picture of the current
Buddha, and its 3D geometric model.
We synthesized the original state by morphing the 3D
mesh of the model from this mesh model. From someliterature inherited at various temples, we knew the sizesof various face parts such as the nose and mouth.

“Enryaku-so-rokubun,”“Daibutsuden-hibun,”“Hichidaiji-nikki,” and “Gokokuji-honnsyoji-enngishu”
are representative documents that contain those sizes.
Unfortunately, however, those numbers often contradict
each other. Some researchers investigated which number
is the most reliable one. We followed their method tocompare them and determined a common figure for each
part.
Table 2 shows the obtained estimated and the current
dimensions of various face parts. Here, all the documents
employ the unit called “shaku.” We interpreted shaku as
the tempyo shaku, and one shaku is assumed to be 0.2964
meters among the various interpretations of shaku. Noticethat relatively large differences exist in heightmeasurements.
Using these data, we designed a two-step morphingalgorithm. First, we globally changed the scale of thewhole portions (for example, Height when sitting, Face
Length, Nose Length); these are gradually modified. In
the second stage, vertices were moved one by one
iteratively, similarly to the constraint propagation
algorithm, using smoothness and uniform constraints. The
two-stage morphing enabled us to obtain the complete
model of the original Great Buddha. Figure 20 shows the
3D models of the current (a) and the original Great
Buddha (Tempyou Buddha) (b). We can easily recognizethat the original Buddha is larger and rather thin.
4.2. Restoring the Toudaiji Main Hall
The main hall of the Toudaiji Temple was built during
the same decades as those of the Great Buddha (8thcentury). It was also rebuilt twice: in the 12th and 17th
centuries. In the 12th century, Tenjiku architecture wasimported from China, and the main hall was rebuilt in a
totally different architecture style. The rebuilding in the
18th century followed the same new style. As a result, thestyle of the current main hall is entirely different from thatof the original building.

Table 2. Current and estimated dimensions of various face parts.
Parts Name

Current (m)Original (m)Height when sitting

14.98

15.85

Eye length

1.02

1.16

Face length

3.20

2.82

Ear length

2.54

2.52

Palm length

1.48

1.66

Foot length

3.74

3.56

Nose height

0.50

0.47

Mouth length

1.33

1.10

International Journal of Computer Vision, 75(1), 189–208, 2007(a)
Figure 19.

(b)

Nara Buddha (a) Picture of current Buddha, (b) 3D geometric model of current Buddha.

(a)Figure 20.

(b)Comparison in 3D models. (a) current buddha, (b) original Buddha.

Fortunately, the Toudaiji temple has been displaying a
miniature model of the original hall, constructed for theParis Expo in 1900, as shown in Fig. 21. We digitized itusing the Pulsteck TDS-1500 and scaled it up to theoriginal size as shown in Fig. 22(a). The TDS-1500 can
scan a range from 3.5 meters through 10 meters with theaccuracy of 0.5mm to 5mm and the spatial resolution of
420 X 280. We obtained 12 range images from various
observation directions. As shown in Fig. 22(b), due to thelimits of the sensor’s accuracy and constraints ofobservation directions, though the model provides roughdimensions of locations of columns and walls, it does notprovide a precise and accurate picture of the detailed
parts.
According to Prof. Keisuke Fujii, who is an architecture
professor at the University of Tokyo and one of theexperts on building style in the era, the Toudaiji andToushoudaiji temples share a similar format. The main
hall of the Toushoudaiji Temple was also built during thesame period (8th century). We have decided to combine
the detailed part model of the Toushoudaiji Temple withthe rough whole model of the Todaiji temple.
We digitized various key parts of the main hall at
Toushoudaiji. Using the suggestions of Prof. Fujii, we
chose 20 important parts of the main hall. Figure 23
shows 6 parts among 20 important parts. We employed
Cyrax 2004 and Pulsteck TDS-1500, which have a range
from 0.5 meter through 1meter, with resolution of 0.23
mm through 0.83 mm, to obtain 780 range images. Figure
24 shows the obtained range images of the detailed parts.
We pasted these partial range data of Touhoudaiji parts(Fig. 24) to the scaled-up range data of the Toudaiji (Fig.
22(b)), using as a scale the average size difference
between those temples, roughly 1 to 2.3. But each part

International Journal of Computer Vision, 75(1), 189–208, 2007

Figure 21.

Miniature model of Buddha palace.

(a) a cloud of points representation(b) mesh model (Closed up)
Figure 22.

3D model acquired from the Miniature model.

needed more precise scaling parameters. The traditionalalignment algorithm determined translation and rotation
parameters as six unknown parameters. However, we
designed an extended alignment algorithm that determined
not only translation and rotation parameters but also scaleparameters.

( r

r

r

r

r

)

ε 2 = min ∑ Rm n ⋅ {( Rs K ( y, k s ) + ts ) − ( Rm x + tm )}R ,ts ≠ m, k2(24)

where K (, ) is a scaling function to expand each arc
length in a mesh, and k s is an unknown scaling parameter.
By using this extended alignment algorithm, we
completed the 3D model of the original Buddha’s palace.

International Journal of Computer Vision, 75(1), 189–208, 2007

Figure 23. Key parts of the main hall at Tousho-daiji digitized.

Figure 24.

3D models of those key parts.

Figure 25(a) shows the original Buddha’s palace
digitally restored by our method. By combining theoriginal Buddha’s palace and the original Buddha, wecreated the virtual appearance of the Nara Buddha in the
8th century, as shown in Fig. 25(b-c). The virtualappearance of this and other historic objects can be used
for education about and promotion of our cultural
heritage.

4.3. Analysis
As one of the demonstrations of utilizing digital
restoration, we conducted an experiment to determine theamount of gold used to plate the surface of the Buddha. It
is well known that the original Buddha was golden due togold plating of its surface. However, several contradictory
numbersexist
indocuments.
For
example,“Daibutu-den-hibun” and “Enryaku-sorokubun” say it
required 5412 ryou and 4187 ryou of gold, respectively, tocover the body of the great Buddha. Moreover, there wereInternational Journal of Computer Vision, 75(1), 189–208, 2007(a) Outside appearance(b) Insider appearance(c) Inside appearance
Figure 25.

Restored nara Buddha.

International Journal of Computer Vision, 75(1), 189–208, 2007

Table 3.

Four interpretations determining the amount of gold used.

Document nameInterpretation

Written amount in the
document

Large ryou (42g)

Small ryou (14g)

Daibutu-den-hibun

5412 ryou

227 kg

76 kg

Enryaku-sorokubun

4187 ryou

176 kg59 kg

two interpretations of “ryou”; A large ryou was 42 g,
while a small ryou was 14 g. Thus, there are fourinterpretations determining the amount of gold required.
In order to disambiguate this discussion, we used our
restored digital model of the Tempyou great Buddha. The
surface area, 597 m 2 , is obtained from the restored digitalmodel by taking a summation of all surface areas oftriangular meshes. For comparison, the surface area of thecurrent Buddha is 556 m 2 . From the documents, it is
known that the amalgam method was used to put gold
over the Buddha’s surface. Usually, this method puts
6 ~ 10 mg / cm . This number was also confirmed byexamining the thickness of gold plate on various treasures
stored in Sho-so-in. By multiplying the surface area of
Tempyo and the current Buddha with this number, we
obtained the gold amount as 36kg ~60kg and 33kg ~ 56kg.
Those numbers indicate that the interpretation of“enryaku-sorokubun” with a small ryou is most likely.

Acknowledgments
This research is sponsored, in part, by JST under Ikeuchi
Crest program, and, in part, by Ministry of Education
under Leading Project. The Bayon temple in Cambodia
was digitized with the cooperation of the Japanese
Government Team for Safeguarding Angkor (JSA).
References

2[1] T. Kanade, P. Rander, and P.J. Narayanann, “Virtualized Reality:
Constructing Virtual Worlds from Real Scenes,” IEEE Multimedia,
4(2):34-47, Januray 1997.
[2] K.N. Kutulakos and S.M. Seitz, “A Theory of Shape by Space
Carving,” Intern. Journal Computer Vision, 38(3): 199-218, July
2000.
[3] E. Borovik and L. Davis, “A Distributed System for Real-timeVolume Reconstruction,” Proc. Of CAMP2000, pp.183-189, 2000.

5. Conclusion
This paper introduced our project to digitally archiveand restore cultural heritage objects. Our project’s main
goal was to develop software to create VR models ofheritage objects through observation of real heritage
objects. As the input sensor of geometrical information,
we used laser range sensors because of their accuracy.
Since only partial mesh models of an object are obtained
from such sensors, we have developed post-process
algorithms. These included a rapid alignment algorithm
based on graphics hardware, a parallel alignment
algorithm based on a PC cluster, and a parallel merging
algorithm based on a PC cluster. For texturing color
information onto geometric modes, we developed a
non-calibrated texturing method based on laser reflectancefeatures.
Digital restoration of lost cultural heritage objects has agreat advantage compared with other restoration methods
such as physical construction of actual temples, because
we can examine various hypotheses without any physicalchanges or long building periods. We demonstrated theeffectiveness of this method through the restoration of the
Nara Great Buddha and its main hall.

[4] I. Kitahara and Y. Ohta, “Scalable 3D Representation for 3D Videoin a Large-Scale Space,” Presence, 13(2): 164-177, 2004
[5] T. Matsuyama, X. Wu, T. Takai, and S. Nobuhara, Real-Time 3D
Shape Reconstruction, Dynamic 3D Mesh Deformation, and High
Fidelity Texture Mapping for 3D Video,” IEEE Trans. Circuits andSystems for Video Technology, Vol. CSVT-14(3): 357-369, March
2004.
[6] K. Ikeuchi and Y. Sato, Modeling from Reality, Kluwer Academic
Press, 2001.
[7] M. Levoy et. al., “The digital Michelangelo project,” SIGGRAPH
2000, pp.131-144, New Orleans.
[8] J. Wasserman, Michelangelo’s Florence Pieta, Princeton University
Press 2003.
[9] I. Stamos and P. Allen, “Automatic registration of 2-D with 3-Dimagery in urban environments,” ICCV2001, pp.731-737,
Vancouver.
[10] P.J. Besl and N.D. McKay, "A method for registration of 3-d
shapes," IEEE Trans. Patt. Anal. Machine Intell., 14(2):239-256,
1992.
[11] R. Benjemma and F. Schmitt, “Fast global registration of 3D
shample surfaces using a multiple-z-buffer technique,” Int. Conf on
Recent Advances in 3-D Digital Imaging and Modeling, pp.
113-120, May 1997.
[12] P. Neugebauer, “Geometrical cloning of 3D objects via
simultaneous registration of multiple range images,” Int. Conf on
Shape Modeling and Application, pp.130-139, March 1997.

International Journal of Computer Vision, 75(1), 189–208, 2007[13] Y. Chen and G. Medioni, “Object modeling by registeration ofmultiple range images, Image and Vision Computing,10(3):145-155, April 1992.
[14] S. Rusinkiewicz and M. Levoy, “Efficient variants of the IPC
algorithm,” Int. Conf 3-D Digital Imaging and Modeling,pp.145-152, May 2001.
[15] H. Gagnon, M. Soucy, R. Bergevin, and D. Laurendeau,“Registeration of multiple range views for automatic 3-D modelbuilding,” CVPR94, pp.581-586.
[16] K. Nishino and K. Ikeuchi, "Robust Simultaneous Registration ofMultiple Range Images", Fifth Asian Conference on Computer
Vision ACCV '02, pp454-461, 2002.
[17] T. Oishi, R. Sagawa, A. Nakazawa, R. Kurazume, and K. Ikeuchi,“Parallel Alignment of a Large Number of Range Images on PC
Cluster,” Int. Conf 3-D Digital Imaging and Modeling, pp. 195-202,
Oct 2003.
[18] M. D. Wheeler and K. Ikeuchi, "Sensor Modeling, Probablistic
Hypothesis
Generation,andRobustLocalizationfor
ObjectRecognition", IEEE PAMI, 17(3): 252-265, 1995.
[19] B. Curless and M. Levoy, “A volumetric method for building
complex models from range images,” SIGGRAPH 96, pp.303-312,
New Orleans, LA.
[20] M. Wheeler, Y. Sato, and K. Ikeuchi, “Consensus surfaces formodeling 3D object from multiple range images,” ICCV98,
pp.917-924.
[21] R. Sagawa, K. Nishino, M.D. Wheeler and K. Ikeuchi, "Parallel
Processing of Range Data Merging", IEEE/RSJ International
Conference on Intelligent Robots and Systems, Vol. 1, pp577-583,
2001[22] R. Sagawa, T. Masuda, and K. Ikeuchi, “Effective Nearest
Neighbor Search for Aligning and Merging Range Images,” Int.
Conf 3-D Digital Imaging and Modeling, pp. 79-86, Oct 2003[23] R. Sagawa and K. Ikeuchi, “Taking Consensus of Signed DistanceField for Complementing Unobservable Surface,” Int. Conf 3-D
Digital Imaging and Modeling, pp. 410-417, Oct 2003[24] R. Sagawa and K. Ikeuchi, “Adaptively Merging Large-Scale
Range Data with Reflectance Properties,” IEEE Trans. Patt. Anal.
Machine Intell., 27(3): 392-405, March 2005.
[25] M. Kamakura, T. Oishi, J. Takamatsu, and K. Ikeuchi,“Classification of Bayon faces using 3D models,” Virtual Systems
and Multimedia, pp.751-760, October 2005.
[26] R. Kurazume, K. Nishino, Z. Zhang, and K. Ikeuchi, "Simultaneous
2D images and 3D geometric model registration for texture
mapping utilizing reflectance attribute," Fifth Asian Conference on
Computer Vision, Vol. 1, pp. 99-106, 2002.
[27] http://www.cvl.iis.u-tokyo.ac.jp