Flying Laser Range Sensor for Large-Scale
Site-Modeling and Its Applications in Bayon
Digital Archival ProjectAtsuhiko Banno, Tomohito Masuda, Takeshi Oishi and Katsushi Ikeuchi
Abstract
We have been conducting a project to digitize the Bayon temple, located at the center of Angkor-Thom in the kingdom of Cambodia. Thisis a huge structure, more than 150 meters long on all sides and up to 45
meters high. Digitizing such a large-scale object in ﬁne detail requires developing new types of sensors for obtaining data of various kinds related
to irregular positions such as the very high parts of the structure occluded
from the ground. In this article, we present a sensing system with a moving platform, referred to as the Flying Laser Range Sensor (FLRS), forobtaining data related to these high structures from above them. The
FLRS, suspended beneath a balloon, can be maneuvered freely in the sky
and can measure structures invisible from the ground. The obtained data,
however, have some distortion due to the movement of the sensor during
the scanning process. In order to remedy this issue, we have developed
several new rectiﬁcation algorithms for the FLRS. One method is an extension of the 3D alignment algorithm to estimate not only rotation andtranslation but also motion parameters. This algorithm compares range
data of overlapping regions from ground-based sensors and our FLRS.
Another method accurately estimates the FLRS’s position by combining
range data and image sequences from a video camera mounted on the
FLRS. We evaluate these algorithms using a hardware-based method andverify that both methods achieve much higher accuracy than previous
methods.

1

Introduction

Preserving cultural heritage objects is a very important mission for us because
they are deteriorating or being destroyed. Digital preservation, measuring suchobjects in three dimensions and representing them in digital forms, is one of thebest ways to accomplish this preservation.
For the past seven years, we have been conducting a project entitled ”theGreat Buddha Project[25][14]” to develop sensors and software for digital preservation of large-scale outdoor structures of cultural interest. Along with thisproject, we have also archived several cultural heritage objects such as the
Kamakura Great Buddha in Kamakura, Japan, and the Atchana Buddha in
Skhotai, Thailand. We are currently extending these projects by attempting
the challenging project of digitizing the Bayon temple in the Angkor ruin in
Cambodia.
1

The Bayon Temple(Fig.1) is located at the center of Angkor-Thom in theKingdom of Cambodia and unites the outlook and traditions of ancient India
and the Khmer. The temple, constructed around the end of the 12th century, iswell known for its size (more than 150 meters long and up to 45 meters high),
its 51 towers, its 173 calm, smiling faces carved on the towers, and its double
corridors carved in beautiful and vivid bas-relief.

Figure 1: The Bayon Temple

We have been digitizing the Bayon Temple with the cooperation of theJapanese government team for Safeguarding Angkor (JSA) since 2002. The motivation for this digitization project is twofold: the historical value of digitallypreserving the Bayon Temple in the current form and the technical challenge ofapplying current technologies to large-scale cultural heritage objects. Regarding
the technical side, we have encountered various issues and derived solutions tomeet three goals:
1. To extend the geometric pipeline to determine the shape of extremely
large-scale objects
2. To extend the photometric pipeline to be able to texture the site in thereal color against strong and variable environmental illuminations
3. To develop new types of sensors that can measure regions invisible fromconventional ground-based laser range sensors
For the ﬁrst goal, alighnment, which determines relative relations between
range data sets, is one of the main components of our geometric pipeline. Traditional alignment algorithms[2][39][4] cannot be applied to our numerous range
data sets. We developed two kinds of alignment algorithms for large-scale sitemodeling. One is a rapid alignment algorithm using graphics hardware[29]. The
other is an accurate simultaneous registration algorithm running on a PC cluster
reducing the inﬂuence of error accumulation [27][28].
With respect to the second goal, two issues, illumination and size, need
to be considered in handling a large outdoor object such as the Bayon Temple. Compensation of illumination eﬀects are necessary[17] because illumination2

conditions vary during the picture-taking session. The second issue is the size:aligning more than 1000 color images over one single geometric model needs a
global bundle adjustment[18].
With respect to the third goal, digitizing large-scale structures requires developing new types of sensors for obtaining range data of various kinds situatedin irregular positions. One class of such irregular positions is the very high portions of the structures occluded from the ground. Figure 2, for instance, showsa part of the 3D model taken by conventional laser range sensors, such as commercial sensors, which require a stable base. There are several regions without
range data at the roofs, which are invisible from the ground. For scanning highpositions, temporary scaﬀolds are often built around structures. However, thisscaﬀold method is impractical for large cultural heritage objects because they
detract from the appearance of these priceless objects and also shaking occurs
at the top of the high scaﬀolds.

Figure 2: Lack of data for regions invisible from the groundTo overcome this diﬃculty, we proposed a novel 3D measurement system,a Flying Laser Range Sensor (FLRS)[13], suspended from the underside of aballoon platform (Fig.3). Several helicopter-based sensors[36][24] have been
proposed for large-scale site modeling. We avoided using a helicopter as theplatform for the sensor for three reasons. First, the high frequency vibrations
from the helicopter engine make it diﬃcult to design rectiﬁcation algorithms.
Second, even though helicopters have high maneuverability, it is still not entirely
safe to use them in proximity to very important cultural heritage structures.
Third, helicopters are considered as strategic, and we cannot import and export
them freely. Thus, our sensor was designed with a balloon platform.
Measurement by the FLRS has a unique problem: the FLRS moves during
its scanning and thus provides distorted range data, as shown in Fig.4. There
are several works dealing with shape rectiﬁcation and alignment of non-rigid
objects [8][35][10][6][16]. However, because we deal with deformed range data
due to movements of the range sensor during the scanning process, conventional
rectiﬁcation and alignment methods for non-rigid models are not appropriate
for our objective. We therefore developed new rectiﬁcation algorithms suitable
for the distorted range data taken by the moving range sensor. In this article,
we describe two methods to rectify 3D range data.
• 3D alignment (3DA) based method• Structure from motion (SFM) based method3

Figure 3: The Flying Laser Range Sensor

These methods do not rely on motion sensors: our FLRS combined with software
rectiﬁcation can achieve higher accuracy than hardware solutions. The two
methods have slightly diﬀerent purposes: one emphasizes accuracy and the otheremphasizes ease of use. The ﬁrst method, 3DA-based, can be applied easily, but
the results are less accurate, especially under strong wind. The second method,
SFM-based, can achieve high accuracy, but it requires preparation for scanningsuch as the calibration and synchronization of sensors.

Figure 4: A sample snapshot compared with the distorted range data obtainedby the FLRS

The ﬁrst method, the 3DA-based method, is an extension of our 3D alignment algorithm, 3DTM[37]. This algorithm assumes that the FLRS is usuallyused in conjunction with a ground-based sensor. It applies to data related tostructures that are both visible and invisible from the ground-based sensor, andto data from overlapping regions. The 3DTM algorithm determines rotationand translation parameters and compares data in the overlapping regions of4

two range images, while the 3DA algorithm determines motion parameters in
addition to rotation and translation parameters[22].
The second method, the SFM-based method, uses distorted range data gathered by a moving range sensor and image sequences gathered by a video camera
mounted on the FLRS. Combining image sequences and range data, we estimate the sensor motion parameters. This method can obtain camera motion
parameters of much greater accuracy than those obtained by a full-perspective
factorization[11], with beginning factorization results as initial estimations. Finally, using reﬁned camera motion parameters, the distorted range data arerectiﬁed[1].
We evaluate these algorithms using a hardware-based method and verify thatthe proposed methods achieve high accuracy. Then we compare two proposed
methods.
This article is organized as follows: Section 2 explains the design issues of
our FLRS system. The ﬁrst method, the 3DA-based method, is presented in
Section 3. In Section 4, we describe the SFM-based rectiﬁcation method, which
utilizes image sequence and range data. Experimental results by both methods
are described in Section 5, in which we compare the two methods as well as the
hardware-based method. Finally, in Section 6 we present our conclusions as well
as some digital photographs of the Bayon temple created from our 3D data.

2

Flying Laser Range Sensor Hardware

The FLRS (Flying Laser Range Sensor) was developed to measure large-scale
objects from the air by using a balloon as a base rather than constructing scaffolds. With respect to the measurement principle, the passive stereopsis method
could capture images without the inﬂuence of balloon motion. However, its results would lack the accuracy needed for the preservation and repair of cultural
heritage objects. The laser radar method is suitable for outdoor measurementof large objects, and we therefore adopted a laser range sensor that uses the”time-of-ﬂight” principle. Figure 5 shows a close-up of the FLRS and a scene
showing the measurement process being carried out in the Bayon temple. A
video camera is mounted on the platform near the range sensor (Fig.5(a)) tocheck the scanning area and to stock the image sequences.
We have two types of FLRSs, which have diﬀerent range capacities. EachFLRS is composed of a scanner unit, a controller and a personal computer (PC).
These three units are suspended beneath a balloon.
The scanner unit includes a laser range ﬁnder specially designed to be suspended from a balloon. Figure 6 shows the interior of the scanner unit. It
consists of a spot laser radar unit and two mirrors. We chose the LARA25200
and LARA53500 supplied by Zoller+Fröhlich GmbH[40] as laser radar units because of their high sampling rate. Each laser radar unit is mounted on an FLRS
scanner unit. The two systems equipped with LARA25200 and LARA53500 are
respectively referred to as the ”25m sensor” and the ”50m sensor”.
There are two mirrors inside each unit to give direction to the laser beam.
One is a polygonal mirror with four reﬂection surfaces, which determines theazimuth of the beam. In normal use, this mirror, which rotates rapidly, controls
the horizontal direction of the laser beam. The other is a plane mirror (swing
mirror) that determines the elevation of the laser beam. The plane mirror swings

5

(a)

(b)

Figure 5: FLRS. (a) Close-up. (b) Measurement scene.
Table 1: The speciﬁcations of the 25m and 50m Sensors
25m Sensor
50m Sensor
Ambiguity interval
25.2 m53.5 mMinimum range
1.0 m
1.0 m
Resolution
1.0 mm
1.0 mmSampling rate ≤ 625,000 pixels/sec. ≤ 500,000 pixels/sec.
Linearity error≤ 3 mm≤ 5mm
Range noise at 10m≥ 1.0 mm
≥ 1.5mm
Range noise at 25m≥ 1.8 mm
≥ 2.7mm
Laser output power
23 mW
32mW
Laser wavelength
780nm
780nm
Angle Resolution Horizontal
0.05◦
0.05◦
◦
Vertical
0.02
0.02◦
◦
Horizontal ﬁeld
≤ 90
≤ 90◦
◦
Vertical ﬁeld
≤ 30
≤ 30◦
Scanning period/range image≤ 15 sec.
≤ 1 sec.

slowly to control the vertical direction of the laser beam. Combining two mirrors,
both sensors take a range image in a one-second measurement with more pixels
along the horizontal direction (900 pixels) than vertical direction (160 pixels).
The view angle of this image is 45◦ horizontal and 30◦ vertical. In addition, the
25m sensor takes a range image of 900 × 800 pixels with the same view angle ina ﬁve-second measurement. The speciﬁcations of two units are shown in Table
1.
Figure 6 shows the interior view of the 25m sensor. The laser beam emitted
from the LARA is ﬁrst directed to a surface of the polygonal mirror. Then thepolygonal mirror reﬂects the laser beam into the plane mirror. The plane mirror
also reﬂects the beam into the outside of the unit (lower part of Fig.6).

6

Figure 6: Interior of a scanner unit (25m sensor)

2.1

Controller Unit and PC

The controller is composed of a signal processing unit, an interface unit, amirror controller, and a power supply unit. The signal processing unit receives
the signals from the PC and performs actual control of rotation angles of themirrors and the laser radar unit. The range data obtained by the laser radar andthe angle data obtained by the mirror encoders are synchronized and combinedin the interface board.
The PC includes a CPU board, an image capture board, and several PCI
boards. Main commands for the mirror operations and the laser on/oﬀ aresent through the LAN cable between the PC and the controller unit. The
synchronized range data and encoder data set are transmitted to the PC via
one of the PCI boards. The data set is stored in the PC and converted into
3-dimentional shape data. The PC on the balloon platform is operated remotely
via another laptop PC on the ground through a LAN cable.

2.2

Operation

The balloon is ﬁlled with helium gas; the diameter of the balloon is 5.0 meter andits maximum buoyancy is about 60kgf . The balloon is made of a particularly
ﬂexible chloroethene, which avoids rapid expansion of a hole in an emergency.
Floating in the air, the balloon is controlled by several hands on the groundwith four peaces of rope. To scan an objective region, the FLRS is broadly
moved to a nearby site by humans. The facing direction of the FLRS is controlled by the pan and tilt mechanism, which can point the scanner unit fromthe horizontal direction to the area directly below, with a scope of 180 degrees

7from side to side. Monitoring the video images, we can adjust the direction.
The combination of the polygonal mirror’s rotation and the plane mirror’sswing makes up a range image that alternates between a raster scan order andits reverse order. This mechanism can record the depth value and the timingof each pixel in the range image because the sampling rate of the laser radar
is known. In addition, the video camera on the FLRS is synchronized withthe laser radar; the image capture board on the PC begins to stock the image
sequence in 30fps simultaneously with the start of the scanning process.

3

3D Alignment with Deformation Parameters

In the case of measuring large-scale objects, we can utilize some range data setstaken by another range sensor ﬁxed on the ground. The ground-base sensors canmeasure many of the regions; our FLRS was originally devised to complement
the measurement for regions that are invisible from the ground. Some parts ofa range image taken by the FLRS are also taken by another range sensor ﬁxed
on the ground. Based on these overlapping regions, we derive an algorithm that
rectiﬁes the distorted range data obtained by the FLRS, and aligns them to therange data taken by the ground-based sensors.
In the ﬁrst part of this section, we describe a robust estimate method of thetranslation and rotation parameters for a conventional non-rigid body alignment. The second part deals with a preliminary alignment between a groundbased sensor’s data and deformed data. In the last part, we extend the method
to the case of the FLRS. The extended method rectiﬁes the distorted FLRS
data and aligns them to the other range data of the ground-based sensor.

3.1

Robust Determination of Translation and Rotation
Parameters

Several alignment algorithms that estimate translation and rotation parameters
for the registration of rigid bodies have been proposed. Besl and McKey designed an alignment algorithm, referred to as the Iterative Closed Point (ICP)
algorithm, of a pair of range data by ﬁnding corresponding point pairs and determining rotation and translation parameters to reduce the sum of squared distances (SSD) among these pairs [2]. However, this algorithm only handles cases
when one data set is a subset of another data set. Zhang extended the alignment
algorithm to handle general cases by introducing statistical measures[39]. Chen
and Medioni proposed to use the distance between a point and plane instead ofa point to point pair[4]. Wheeler and Ikeuchi employed the Lorentzian distanceinstead of SSD to avoid the eﬀect of outliers[37]. A recent survey of alignment
algorithms can be found in [32].
The original algorithm aligns the closest points together with two data sets
so as to minimize the distance summation of point-to-point distance[2]. The
minimization of the error function is represented as follows[37]:1 ∑ρ(zi (p)),
N i
N

E(p) =

8(1)

wherep = (T, q),zi (p) = ||R(q)xi + T − yi ||2 ,1
ρ(zi (p)) = log(1 + zi (p)),
2
T
R(q)
xi
yi
N

:
:
:
:
:

(2)
(3)
(4)

translation vector (sensor position),rotation matrix corresponding to quaternion q,
ith point in the data set of interest,
the corresponding point of xi in the measured data,
the number of data point,

The range images are aligned iteratively by moving (translating/rotating) the
measured data according to the estimated parameters. The movement is determined such that the total distance sum between the corresponding points is
minimized. As for the rotation matrix, we use the quaternion representation of
3 degrees of freedom. In the direct square summation of error function Eq.(3),considerable noise leads to imprecise alignment. In this algorithm, M-estimation(Lorentz function) is used for noise elimination by considering the probability
distribution of the error[37]. Using this error metric E(p), we compute theparameters p that fulﬁll the following equation:
popt = arg min E(p).
p(5)

For the gradient-based solution of non-linear optimization[31], the descentgradient is computed as follows:
∂E
∂p

=

1 ∑ ∂ρ(zi ) ∂zi
·
N i
∂zi
∂p

=

1 ∑
∂zi
w(zi )zi
,
N i
∂p

N

N

where w(zi ) =(6)

1 ∂ρ(zi )
·
.
zi
∂zi

If we evaluate ∂zi /∂p by identifying quaternion qI , we can represent ∂zi /∂p as¯∂(R(q)xi + T − yi ) ¯¯∂zi (p)= 2(R(q)xi + T − yi )
¯
∂p
∂p
qI
[
]
2(xi + T − yi )=.
(7)−4xi × (T − yi )The detailed explanation of the derivation with respect to the quaternion isshown in [38].
Our goal is to simultaneously determine deformation, translation, and rotation parameters by comparing the target data to transform with its corresponding data. The translation and rotation parameters are determined in aminimization paradigm. If we ﬁx these parameters, the determination of thedeformation parameter becomes an iterative shape-matching problem.
93.2

Preliminary Alignment for Deformed Data

We extend the parameter estimation of the alignment formulation for the distorted range data by the moving range sensor. Here, zi (p) in Eq.(3) is transformed into:∑zi (p) =||R(q)g(xi , k) + T − yi ||2 ,
(8)
i

where

p=
g(xi , k) :(T, q, k),deformation function of point xi
with respect to parameter k.

And the gradient described in Eq.(7) is extended as:¯∂zi (p)∂(R(q)g(xi , k) + T − yi ) ¯¯
= 2(R(q)g(xi , k) + T − yi )
¯
∂p
∂p
qI



2(g(xi , k) + t − yi )
−4g(xi , k) × (T − yi )=∂(g(xi ,k))2(g(xi , k) + T − yi ) ∂k

(9)This straightforward extension causes unstable convergence of the deformation alignment because the obtained translation, rotation, and deformation parameters overreach their optimum if every parameter is applied simultaneously
to the deformation: every parameter interferes with the others.
In order to prevent this interference, we design our extended formulation
again to remove the translation and rotation eﬀect caused only by deformation.
The basic idea is to recover the position and posture that change due to deformation. This is implemented by the ”preliminary” rigid-body transformationthat determines only the deformation parameter. First, every parameter is acquired by Eqs.(6) and (9). Then the preliminary rigid-body transformation isdetermined only by the deformation parameter as follows:
g0 (xi , k) = Ro g(xi , k) + To ,

(10)

where
(Ro , To ) = (R(qo ), To ),such that(qo , To ) = arg min
q,T

N
∑

||R(q)g(xi , k) + T − xi ||2 .

(11)

i

Ro , To can be derived from the following equation:
∂

∑N 2
N
∑
∂²i
i ²i
=
2²i ·
= 0,
∂Po
∂Po
i

where
²i = g(xi , k) + To − xi ,Po = (qo , To ).
10

(12)

This is a conventional alignment problem, but it is unnecessary to strictly solve
the above equation. In fact, Ro does not aﬀect the stable convergence so muchas To . If Ro is ignored, To in Eq.(12) is concretely derived as follows:N∑(g(xi , k) + T − xi ) = 0i∑N∴ To = −i (g(xi , k) − xi )N.

Finally, Eq.(8) is replaced with:∑zi (p) =||R(q){Ro g(xi , k) + To } + T − yi ||2 .

(13)

(14)

i

The acquisition of a good initial parameter is signiﬁcant for the optimalalignment result. In our implementation, the initial transformation parameteris set manually by GUI with accuracy enough to reach the optimum.

3.3

Simultaneous Determination of Deformation Parameters for FLRSIn addition to the alignment of inter-scanning transformation parameters, wehave to determine parameters due to intra-scanning movement. This movement
causes obtained data to be distorted. In the case of a short-term measurement,
we can suppose that pose and motion change of the FLRS is generated by• Initial velocity• Initial angular velocity• Acceleration generated by external force• Angular acceleration generated by external moment
We can ignore the inﬂuence of translation and angular acceleration because
the FLRS needs only one second to scan, and insigniﬁcant rotation can beapproximated to translation movement. Hence we consider only constant velocity movement. Under this assumption, we set up the deformation equation inEq.(8).
In this case, the geometrical function g(xi , k) is represented only by theconstant velocity vector v of the FLRS movement, and Eq.(8) is replaced with:g(xi , v) = xi − τi v,(15)

where τi is the ith point’s captured time passed since the start of the scanning.
The descent gradient is represented in this case as follows:2(g(xi , v) + T − yi )∂zi (p) ,−4g(xi , v) × (T − yi )
=
(16)
∂p
∂g(xi ,v)2(g(xi , v) + T − yi ) ∂v

11

wherep=(T, q, v),

τi 0 0∂g(xi ,v)
=  0 τi 0  .
∂v
0 0 τi

While translation T and rotation q are due to the sensor movement among
multiple views, shape deformation alignment is due to the sensor movement
during one scan. Thus, we refer to this alignment as ”inter-and-intra scanningalignment.”

4

Rectiﬁcation by Combining Image Sequenceand Range Data

The second method, the SFM-based method, utilizes range data and an image sequence obtained by a video camera mounted on the FLRS as follows:First, the motion of the FLRS is roughly estimated only by the obtained images; we use a full-perspective factorization[5][11]. These approximate motion
parameters are utilized as an initial solution. Then the motion parameters forrectiﬁcation are estimated based on an optimization imposing some constraints,
which include information derived from the image sequence and the distorted
range data. For the sake of an accurate estimation of motion parameters, theoptimization consists of three constraints: range data, balloon motion smoothness, and tracking of interest points. The camera motion parameters can befound through the minimization of a global functional. Finally, by using thecamera motion parameters, the distorted range data are rectiﬁed.

4.1

Problem Setting

We assume that F images are given by the video camera while obtaining a single
range image; the period is deﬁned as Ts . Then P interest points are tracked
over all the frames; each 3D interest point p is expressed as Sp in the world
coordinate system.
The 3D point p is projected on the image plane at (up (t), vp (t))t at time t.
We can obtain these projected points at the discrete time tf , (1 ≤ f ≤ F ) byusing an interest point detector. The relation between the discrete time tf andthe frame number f is described as t = f /30 because we use a video cameraof 30 fps. The positions of these projected points as a continuous time aredetermined using an interpolation technique.
By using the intrinsic camera parameter matrix A, the image projection is
described as follows:up (t)
i(t)t
 vp (t)  ∝ A  j(t)t  (Sp − T(t))
(17)
1
k(t)t

12

t
f
p
T(t)
Sp
i(t), j(t), k(t)

:
:
:
:
:
:

scanning time 0 ≤ t ≤ Ts ,frame number 1 ≤ f ≤ F ,the number of the interest points 1 ≤ p ≤ P ,camera position at t,
the 3D position of p-th interest point,
the x, y and z axis of the camera coordinate system at t,

The range sensor outputs a cloud of 3D points; each point is converted intot
the camera coordinate description x(t) = (x(t), y(t), z(t)) because the range
sensor and the camera are calibrated.

4.2

Motion Estimate by Combining Range Data and Image Sequence

We adopt an optimization strategy to estimate motion parameters for the correct
rectiﬁcation. With respect to set the cost function to minimize, we impose threeconstraints that combine the image sequence and the range data. The ﬁrst is arange data constraint, which utilizes temporally geometric relations between 3D
points and the camera. The second is a smoothness constraint, which guarantees
a smooth motion of the balloon. The third is a tracking constraint, which is
known as the bundle adjustment.
4.2.1

Range Data Constraint

The range sensor outputs a 3D point x(t) at t. Since the range sensor and thevideo camera are calibrated, we can determine the position of the projected
point of x(t) on the image plane. Let us deﬁne the 2D projected point as u(t).
If the point u(t) matches any interest point (up (t), vp (t))t on the image at time
tp , we can judge that the range sensor scans the 3D interest point p at time tp .
There are few points that are scanned by the range sensor but invisible from thevideo camera because the video camera is located near the laser radar (about
10cm) and these points are more than 5 meters away from the sensor in actual
case.
Therefore, we can derive the ﬁrst constraint as follows:
FA =

P
∑

αp kx(tp ) − Rt (tp )(Sp − T(tp ))k2

(18)

p=1

{

where

1 if point p is scanned by the range sensor
0 otherwise()
R(t) = i(t) j(t) k(t) means the rotation matrix of the camera pose at t. As
described above, tp is the scan time when the 3D interest point p is scanned by
the range sensor. Each tp is estimated as follows:()
up (t)
tp = arg min ku(t) −k(19)
vp (t)tαp =It is not always possible to use all interest points in images for this constraint.
Several interest points in images are on 3D corners and 3D edges, which cause
13mismatches in the range image. In order to exclude these interest points on 3D
corners, we set αp = 0 in Eq.18 for the points whose neighboring points in therange image have a variance larger than a threshold. The interest points, which
are located out of the laser’s range, are also excluded.
4.2.2

Smoothness Constraint

One of the most signiﬁcant reasons for adopting a balloon platform is to be
free from the high frequency that occurs with a helicopter platform. A balloon
platform is only under the inﬂuence of low frequency: the balloon of our FLRS
is held with some wires swayed only by wind. This means that the movement of
the balloon is expected to be smooth. Certainly, the movement of the balloon
is free from rapid acceleration, rapid deceleration, or acute change of course.
Taking this fact into account, we consider the following function:)∫ Ts (° ∂ 2 q(t) °2° ∂ 2 T(t) °2
°
°
°
°
FB =
w1
+w2
dt
(20)
∂t2
∂t2
0
Here, T(t) denotes the position of the camera; w1 , w2 are weighted coeﬃcients; and q(t) is a unit quaternion that represents the camera pose. The ﬁrst
term of the above integrand represents smoothness with respect to the camera’stranslation while the second one represents smoothness with respect to the camera’s rotation. When the motion of the camera is smooth, the function FB has
a small value.
4.2.3

Tracking Constraint

Any interest point Sp must be projected near the observed interest points(up (tf ), vp (tf )) on each frame f . This constraint is well known as Bundle
Adjustment[3]. When the structure, motion, and shape have been roughly obtained, this technique is utilized to reﬁne them through an image sequence.
From Eq.17, the following relationship is conducted: 
i(tf )t · (Sp − T(tt ))up (tf )(21)
A−1  vp (tf )  ∝  j(tf )t · (Sp − T(tt )) 1
k(tf )t · (Sp − T(tt ))Supposing (uˆp (tf ), vˆp (tf ), 1) = A−1 (up (tf ), vp (tf ), 1) , we set the third constraint:
(
F ∑
P
t
(
∑
i(tf ) · (Sp − T(tf )) )2FC =uˆp (tf ) −tk(tf ) · (Sp − T(tf ))f =1 p=1)t(j(tf ) · (Sp − T(tf )) )2+ vˆp (tf ) −(22)tk(tf ) · (Sp − T(tf ))tt

The minimization of FC leads to the correct tracking of the interest points
by a moving camera.

144.2.4The Global Cost Function

Based on the above considerations, it will be found that the next cost function
should be minimized. Consequently, the weighted sum
F = wA FA + wB FB + wC FC(23)

leads to the global cost function. The coeﬃcients wA , wB and wC are determined
experimentally so that three terms take almost the same magnitude. Measuringx(t), Sp and T(t) in the metric system and (up (t), vp (t)) in pixel, we set threecoeﬃcients as (wA , wB , wC ) = (102 , 105 , 100 ) and two coeﬃcients in Eq.20 as(w1 , w2 ) = (100 , 101 ) in this study.
To minimize this function, we employ a conjugate gradient method[30][15][33][31]
to ﬁnd the search directions. Then, we use the golden section search to determine the minimum along the directions. For optimization, the LevenbergMarquardt method[21] is generally employed to minimize a functional value.
The Levenberg-Marquardt method is very eﬀective in estimating a function’s
parameters, especially to ﬁt a certain function. However, in our function, it
is not a parameter ﬁtting problem to minimize the value of FB , but a simple
decreasing of FB . Therefore we adopt the conjugate gradient method.

4.3

Initial Estimation

To minimize the value of Eq.23, we need a proper initial guess. With an inproper
initial value, the minimization will be trapped into local minimums. In thisstudy, we utilize the solution by full-perspective factorization[5][11] as the initialvalue to minimize the global cost function. Using the weak-perspective projection model, the full perspective factorization iteratively estimates the shape andthe camera motion under the perspective model.
For the factorization, we need P interest points tracked from F frames. There
are several methods to derive interest points of images[26][34][20]. Among them,
we adopt Harris operator[12] for derivation of interest points. Over 300 interest
points are derived from each frame; identiﬁed points are connected by using
window matching. In order to reject outliers, we impose the Epipolar constraint
for all adjacent frame pairs. Moreover, we adopt RANSAC[9] to estimate all
F-Matrices.

5

Experimental Results

We show the experimental results of two data sets; these data sets were obtained in the case of a balloon with a moderate motion and a wide motion. We
compared the rectiﬁed shapes resulting from our two proposed methods anda hardware-based method. Our FLRS system is equipped with a regular andan angular accelerometer. To demonstrate the advantage of our methods, we
utilized the rectiﬁed shape produced by these physical devices.
In all the experiments, it took one second to scan a range image of 900 × 160
pixels; the horizontal view angle is 45◦ and the vertical angle is 30◦ . The video
camera stocks 70 images with 30 fps including the one second of the scanning
period. For calibration of the range sensor and the video camera, we utilized
the range data obtained by the FLRS ﬁxed on the ground.

15

Case 1The data set of Case 1 was obtained when the balloon motion was moderate.
Figure 7 (a) shows a photograph of the scanned area. The dense model of Fig.7(b) is obtained by the Cyrax2500[19] ﬁxed on the ground. We treat this modelas the correct data hereafter.

(a)

(b)

Figure 7: Case1: (a) The object for the experiment. (b) The 3D model by theground-base sensor.
The original distorted range data and some rectiﬁed shapes are shown in
Fig. 8. Fig.8 (a) shows the original distorted range data obtained by the FLRS.
The shape of Fig.8 (b) is the rectiﬁed shape based on the measurements of
the acceleration meter and the angular accelerometer. Fig.8 (c) is the shaperectiﬁed by the 3DA-based method. To produce the rectiﬁed shape of Fig.8(c), we utilized the shape of Fig.7 (b). Fig.8 (d) is the shape rectiﬁed by the
SFM-based method of combining the range data and the image sequences of the
FLRS.
These shapes in Fig.8 do not clearly demonstrate the eﬀectiveness of therectiﬁcations due to a moderate balloon motion. To determine these, we evaluated the validity of the rectiﬁed shapes numerically, using geometric similarities
between the correct shape and each rectiﬁed shape.
The results are shown in Table 2. After aligning each rectiﬁed model to thecorrect one, the distances between the corresponding pairs are calculated. The
values in Table 2 show the percentages of the pairs with closer distances than
the thresholds among all pairs. We set up three thresholds as 1.0 cm, 5.0 cm
and 10.0 cm. When we scanned objects using the FLRS ﬁxed on the ground,the accuracy was estimated at about 1.0 cm. We think that a 3D model withaccuracy of less than 5.0 cm would be suﬃcient for modeling cultural assets.
More than half of the points in the original distorted range data have a 1.0
cm accuracy since this data set was obtained under moderate balloon motion.
Using three rectiﬁcation methods, the areas with 1.0 cm accuracy are increased.
The hardware-based method can rectify the distorted range data; the advantage

16

(a)

(b)

(c)

(d)

Figure 8: Case1: (a)The original distorted range data by FLRS. (b)The range
data rectiﬁed by the hardware-based method. (c)The range data rectiﬁed by
the 3DA-based method. (d)The range data rectiﬁed by the SFM-based method.
of this method is its rapidity because it does not need any iterative calculations.
In this data set, the hardware-based method shows high performance. The 3DAbased method does not increase the points with 1.0 cm accuracy; nevertheless,this method increases the points with 5.0 cm accuracy almost ninety percent of
the time. The SFM-based method demonstrates even greater eﬀectiveness thanthe hardware-based method: more than eighty percent of points have 1.0 cm
accuracy.
Table 2: The evaluation of the rectiﬁed models (Case1). The percentages of thepairs with closer than the thresholds.
Original Hardware-based 3DA-based SFM-based≤ 1.0cm (%)
58.08
67.23
60.20
80.42
≤ 5.0cm (%)
76.45
93.16
89.93
96.18
≤ 10.0cm (%)
79.28
97.74
97.13
98.49

The appearances of the regions with 1.0 cm accuracy are shown in Fig.9.
The white regions indicate the areas of high accuracy of 1.0 cm.
Figure 10 shows two cross-sections of four models: the correct model by
Cyrax2500, the one rectiﬁed by the hardware-based method, the one rectiﬁedby the 3DA-based method, and the one rectiﬁed by the SFM-method.

17

(a)

(b)

(c)

(d)

Figure 9: Case1: Comparisons with the range data obtained by the groundbased sensor. The white regions indicate where the distances of the corresponding pairs are less than 1.0 cm. (a)The original distorted range data by FLRS.
(b)The range data rectiﬁed by the hardware-based method. (c)The range rectiﬁed by the 3DA-based method. (d)The range data rectiﬁed by the SFM-based
method.

Case 2The data set of Case 2 was, on the other hand, obtained when the balloon moved
widely. Figure 11 (a) shows a photo picture of the scanned area. The dense
model of Fig.7 (b) is the correct data by the Cyrax2500.
Figure 12 (a) is the original distorted shape obtained by the FLRS. It isfound that the shape is widely deformed. In Fig.12 (b), the shape rectiﬁed by
the hardware-based method is shown. Unlike the previous case, this rectiﬁcation
method leads to a noticeable deformation. There is almost no variation inappearance of Fig.12 (c). On the other hand, a glance at Fig.12 (d) shows that
the SFM-based method seems to rectify the shape properly.
Table 3 shows the numerical evaluations in this data set. The wide balloon
motion caused a quarter of the area of 1.0 cm accuracy in the original range
data. This value clariﬁes numerically that the hardware-based method makes
rectiﬁcation less eﬃcient. In the case of a wide and rapid balloon motion, it
is found that the hardware-based method is improper to rectify the distorted
FLRS data. The 3DA-based method increases the area of 1.0 cm accuracy,
but it could not reach the level attained by the SFM-based method. Using theSFM-based method, more than half the area has a 1.0 cm accuracy, while morethan eighty percent of points have a 5.0 cm accuracy. We think the reason why
the 3DA-based method could not achieve the performance by the SFM-based
method is as follows: the assumption of the constant velocity movement of theballoon does not hold true in this case. On the other hand, the 3DA-based
method works out in the case of moderate balloon motion.
18

Cutting plane 1

Cross sections by the plane 1

Cutting plane 2

Cross sections by the plane 2

Figure 10: Case1: Cross sections of four models; ”CORRECT”:the correctmodel by Cyrax2500. ”HARD”: the range data rectiﬁed by the hardware-based
method. ”3DA”: the range data rectiﬁed by the 3DA-based method. ”SFM”:
the range data rectiﬁed by SFM-based method.

(a)

(b)

Figure 11: Case2: (a) The object for the experiment. (b) The 3D model by theground-base sensor.
The appearances of the regions of the 1.0cm accuracy in Case2 are shown in
Fig.13.
Figure 14 shows two cross sections of four models in Case2.

19

(a)

(b)

(c)

(d)

Figure 12: Case2: (a)The original distorted range data by FLRS. (b)The range
data rectiﬁed by the hardware-based method. (c)The range data rectiﬁed by
the 3DA-based method. (d)The range data rectiﬁed by the SFM-based method.
Table 3: The evaluation of the rectiﬁed models (Case2)
Original Hardware-based 3DA-based SFM-based≤ 1.0cm (%)
26.30
22.55
35.33
54.30
≤ 5.0cm (%)
55.97
45.47
65.18
86.54
≤ 10.0cm (%)
70.19
58.13
77.10
93.58

Additional Case (Case3)
We demonstrate an additional data set, in which it takes ﬁve seconds to scan arange image of 900 × 800 pixels; similarly the horizontal view angle is 45◦ andthe vertical angle is 30◦ . Since it was impossible to rectify the distorted data by
the hardware-based and 3DA-based methods, we show the range data rectiﬁed
only by the SFM-based method.
Figure 15 shows the scene of Case 3 and the correct range data by theground-base sensor. Although the original data are widely deformed due to theﬁve-seconds scanning, it is found that the SFM-based method rectiﬁes the data
properly.
The original distorted range data and the shape rectiﬁed by the SFM-based
method are shown in Fig.16.
The appearances of the regions of 1.0cm accuracy in Case 3 are shown in
Fig.17, which indicates the high performance of the SFM-based method.

20

(a)

(b)

(c)

(d)

Figure 13: Case2: Comparisons with the range data obtained by the groundbased sensor. The white regions indicate where the distances of the corresponding pairs are less than 1.0 cm. (a)The original distorted range data by FLRS.
(b)The range data rectiﬁed by the hardware-based method. (c)The range data
rectiﬁed by the 3DA-based method. (d)The range data rectiﬁed by the SFMbased method.

6

Conclusios

We have described our novel sensor, FLRS, for measuring large-scale objects.
The FLRS is a new-sprung system in our ongoing project to digitize the Bayontemple ar the Angkor ruin, in Cambodia. The FLRS can digitize objects fromthe air while suspended from the underside of the balloon platform. The advantage of the FLRS over conventional sensors is that the FLRS can easily measure
the regions invisible from the ground. We have also described two kinds ofrectiﬁcation methods to solve the unique problem for the FLRS: the distorted
range data due to the balloon s motion. One is the 3DA-based method, which
rectiﬁes the distorted range data by using another range data set obtained fromthe ground-based sensor; another is the SFM-based method, which estimates
the sensor motion combining the range data and the image sequences. We have
conﬁrmed the eﬀectiveness of these methods by comparing the hardware-based
rectiﬁcation method. Among them, the SFM-based method achieved higher
accuracy for the rectiﬁed range data.
Finally, ﬁgure 18 shows the measurement result reconstructed using all
range data obtained by all sensors, FLRSs, Cyrax, Z+F Imager and Climbing Sensor[23]. This 3D model consists of 20,000 range images obtained in all
missions and the total size is about 200 GByte.

21

Cutting plane 1

Cross sections by the plane 1

Cutting plane 2

Cross sections by the plane 2

Figure 14: Case2: Cross sections of four models; ”CORRECT”:the correctmodel by Cyrax2500. ”HARD”: the range data rectiﬁed by the hardware-based
method. ”3DA”: the range data rectiﬁed by the 3DA-based method. ”SFM”:
the range data rectiﬁed by the SFM-base method.

(a)

(b)

Figure 15: Case3: (a) The object for the experiment. (b) The 3D model by theground-base sensor.

Acknowledgement
The research described herein was supported, in part, by the Japanese Ministryof Education, Culture, Sports, Science and Technology, under the program, ”Development of High Fidelity Digitization Software for Large-Scale and Intangible
Cultural Assets.” Measurement missions in the Bayon temple were conductedin cooperation with the Japanese government team for Safeguarding Angkor

22

(a)

(b)

Figure 16: Case3: (a)The original distorted range data by FLRS. (b)The range
data rectiﬁed by the SFM-based method.

(a)

(b)

Figure 17: Case3: Comparisons with the range data obtained by the groundbased sensor. The white regions indicate where the distances of the corresponding pairs are less than 1.0 cm. (a)The original distorted range data by FLRS.
(b)The range data rectiﬁed by the SFM-based method.
(JSA).

References[1] A. Banno and K. Ikeuchi. Shape Recovery of 3D Data Obtained from a
Moving Range Sensor by using Image Sequences, Proc. of the International
Conference on Computer Vision (ICCV2005), Vol. 1, pp.792–799, 2005.
[2] P. J. Besl and N. D. McKay. A method for registration of 3-D shapes.
IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. 14,
pp.239–256, 1992.
[3] D. Brown. The bundle adjustment – progress and prospect. XIII Congressof the ISPRS, Helsinki, 1976.
[4] Y. Chen and G. Medioni. Object modeling by registration of multiple range
images. Image and Vision Computing, Vol. 10, No. 3, pp.145–155, 1992.
[5] S. Christy and R. Horaud. Euclidean shape and motion from multiple perspective views by aﬃne iterations. IEEE Transactions on Pattern Analysis
and Machine Intelligence, Vol. 18, No. 11, pp.1098–1104, 1996.

23

Figure 18: The Overview of the ”Digital Bayon”[6] H. Chui and A. Rangarajan. A new point matching algorithm for non-rigid
registration. Computer Vision and Image Understanding, Vol. 89, pp.114–
141, 2003.
[7] B. Curless and M. Levoy, A volumetric method for building complex models
from range images. Proc. of SIGGRAPH’96, pp.303-312, ACM, 1996.
[8] J. Feldmar and N. Ayache, Rigid, aﬃne and locally aﬃne registration24

of free-form surfaces. International Journal of Computer Vision, Vol. 18,
No. 2, pp.99–119, 1996.
[9] M. A. Fischler and R. C. Bolles. Random sample consensus: A paradigm
for model ﬁtting with applications to image analysis and automated cartography. Communications of the ACM, Vol.24, No.6, pp.381–395, 1981.
[10] D. Hähnel, S. Thrun and W. Burgard. An Extension of the ICP Algorithm
for Modeling Nonrigid Objects with Mobile Robots. Proc. of the International Joint Conference on Artiﬁcial Intelligence (IJCAI), 2003.
[11] M. Han and T. Kanade. Perspective factorization methods for euclidean
reconstruction. Technical Report:CMU–RI–TR–99–22, Robotics Institute,Carnegie Mellon University, 1999.
[12] C. Harris and M. Stephens. A combined corner and edge detector. Proc. of
Alvey Vision Conference, pp.147–152, 1988.
[13] Y. Hirota, T. Masuda, R. Kurazume, K. Ogawara, K. Hasegawa and K.
Ikeuchi. Designing a laser range ﬁnder which is suspended beneath a balloon. Proc. of the 6th Asian conference on Computer Vision (ACCV2004),
Vol. 2, pp.658–663, 2004.
[14] K. Ikeuchi, A. Nakazawa, K. Hasegawa and T. Ohishi. The Great Buddha
Project: Modeling Cultural Heritage for VR Systems through Observation,
Proc. the 2nd IEEE and ACM International Symposium on Mixed andAugmented Reality (ISMAR2003), 2003.
[15] D. A. Jacobs. The State of the Art in Numerical Analysis, Academic Press,
London, 1977.
[16] V. Jain, H. Zhang, O. van Kaick. Non-Rigid Spectral Correspondence of
Triangle Meshes. accepted to International Journal on Shape Modeling (viainvitation to Special Issue of SMI 2006).
[17] R. Kawakami, R.T.Tan and K. Ikeuchi, Consistent Surface Color for Texturing Large Objects in Outdoor Scenes, Proc. of the International Conference on Computer Vision(ICCV2005), Vol. 2, pp.1200–1207, 2005[18] R. Kurazume, K. Nishino, Z. Zhang and K. Ikeuchi, Simultaneous 2D images and 3D geometric model registration for texture mapping utilizing
reﬂectance attribute, Proc. of Fifth Asian Conference on Computer Vision(ACCV2002), Vol. 1, pp.99–106, 2002.
[19] Leica Geosystems. http://www.leica-geosystems.com/[20] D. G. Lowe. Distinctive image features from scale-invariant keypoints. International Journal of Computer Vision, Vol. 60, No. 2, pp.91–110, 2004.
[21] D. W. Marquardt. An algorithm for least-squares estimation of nonlinear
parameters. Journal of the Society for Industrial and Applied Mathematics,
Vol. 11, pp.431–441, 1963.

25

[22] T. Masuda, Y. Hirota, K. Nishino and K. Ikeuchi. Simultaneous determination of registration and deformation parameters among 3D range images.
Proc. of 5th International Conference on 3-D Digital Imaging and Modeling(3DIM2005), pp.369–376, 2005.
[23] K. Matsui, S. Ono and K. Ikeuchi. The Climbing Sensor: 3-D Modeling of a Narrow and Vertically Stalky Space by Using Spatio-Temporal
Range Image. International Conference on Intelligent Robots and Systems(IROS2005), 2005.
[24] R. Miller and O. Amidi. 3-D Site Mapping with the CMU Autonomous
Helicopter. The 5th International Conference on Intelligent Autonomous
Systems, 1998.
[25] D. Miyazaki, T. Oishi, T. Nishikawa, R. Sagawa, K. Nishino, T. Tomomatsu, Y. Yakase and K. Ikeuchi. The great buddha project: Modellingcultural heritage through observation. Proc. the 6th International Conference on Virtual Systems and Multimedia (VSMM2000), pp.138–145, 2000.
[26] H. P. Moravec. Towards automatic visual obstacle avoidance. Proc. of 5thInternational Joint Conference on Artiﬁcial Intelligence, p.584, 1977.
[27] K. Nishino and K. Ikeuchi, Robust Simultaneous Registration of Multiple Range Images, Proc. of 5th Asian Conference on Computer Vision(ACCV2002), pp454-461, 2002[28] T. Oishi, R. Sagawa, A. Nakazawa, R. Kurazume and K. Ikeuchi, Parallel
Alignment of a Large Number of Range Images, Proc. of 4th International
Conference on 3-D Digital Imaging and Modeling (3DIM2004), pp.195–202,
2003.
[29] T. Oishi, A. Nakazawa, R. Kurazume and K. Ikeuchi, Fast Simultaneous
Alignment of Multiple Range Images using Index Images, Proc. of 5th International Conference on 3-D Digital Imaging and Modeling (3DIM2005),
pp.476–483, 2005.
[30] E. Polak. Computational Methods in Optimization, Academic Press, New
York, 1971.
[31] W. H. Press, B. P. Flannery, S. A. Teukolsky and W. T. Vetterling. Numerical Recipes in C, Cambridge University Press, 1988.
[32] S. Rusinkiewicz and M. Levoy. Eﬃcient variant of the ICP algorithm. Proceedings of the 3rd International Conference on 3-D Digital Imaging andModeling (3DIM2001), pp.145–152, 2001.
[33] J. Stoer and R. Bulirsh. Introduction to Numerical Analysis, SpringerVerlag, New York, 1980.
[34] S. M. Smith and M. Brady. SUSAN - a new approach to low level image
processing, International Journal of Computer Vision, Vol.23, No.1, pp.45–
78, 1997.

26

[35] R. Szeliski and S. Lavallée. Matching 3-D Anatomical Surfaces with NonRigid Deformations using Octree-Splines. International Journal of Computer Vision, vol. 18, no. 2, pp.171–186, 1996.
[36] S. Thrun, M. Diel and D. Haehnel, Scan alignment and 3-D surface modeling with a helicopter platform, Proc, of the 4th International Conference
on Field and Service Robotics, 2003.
[37] M. D. Wheeler and K. Ikeuchi, Sensor Modeling, Probabilistic Hypothesis Generation, and Robust Localization for Object Recognition, IEEE
Transactions on Pattern Analysis and Machine Intelligence, Vol.17, No.3,
pp.252–265, 1995.
[38] M. D. Wheeler, Automatic modeling and localization for object recognition, PhD Thesis, School of Computer Science, Carnegie Mellon University, 1996.
[39] Z. Zhang, Iterative point matching for registration of free-form curves andsurfaces, International Journal of Computer Vision, Vol.13, pp.119–152,
1994.
[40] Zoller+Fröhlich GmbH, http://www.zf-lase.com/

27