Image Anal Stereol 2002;21:31-36 Original Research Paper
VIDEOGRAMMETRIC RECONSTRUCTION APPLIED TO VOLCANOLOGY: PERSPECTIVES FOR A NEW MEASUREMENT TECHNIQUE IN VOLCANO MONITORING
Emmanuelle Cecchi1, Jean-Marc Lavest2 and Benjamin Van wyk de Vries1
1Laboratoire Magmas et Volcans, LMV, UMR 6524 du CNRS, Blaise-Pascal University, Clermont-Ferrand, 2Laboratoire des Sciences et Mate´riaux pour l’Electronique, et d’Automatique, LASMEA, UMR 6602 du CNRS, Blaise-Pascal University, Clermont-Ferrand e-mail: cecchi@opgc.univ-bpclermont.fr (Accepted February 21, 2002)
ABSTRACT
This article deals with videogrammetric reconstruction of volcanic structures. As a first step, the method is tested in laboratory. The objective is to reconstruct small sand and plaster cones, analogous to volcanoes, that deform with time. The initial stage consists in modelling the sensor (internal parameters) and calculating its orientation and position in space, using a multi-view calibration method. In practice two sets of views are taken: a first one around a calibration target and a second one around the studied object. Both sets are combined in the calibration software to simultaneously compute the internal parameters modelling the sensor, and the external parameters giving the spatial location of each view around the cone. Following this first stage, a N-view reconstruction process is carried out. The principle is as follows: an initial 3D model of the cone is created and then iteratively deformed to fit the real object. The deformation of the meshed model is based on a texture coherence criterion. At present, this reconstruction method and its precision are being validated at laboratory scale. The objective will be then to follow analogue model deformation with time using successive reconstructions. In the future, the method will be applied to real volcanic structures. Modifications of the initial code will certainly be required, however excellent reconstruction accuracy, valuable simplicity and flexibility of the technique are expected, compared to classic stereophotogrammetric techniques used in volcanology.
Keywords: calibration, digital elevation model (DEM), videogrammetry, volcanology.
INTRODUCTION
Quantification of morphological changes on volcanoes constitutes an important aspect of monitoring and hazard research. At present, vertical stereophotogrammetric techniques, sometimes using highly advanced equi pment Gwinner et al. (2000), are used for 3D reconstruction. With such techniques it is difficult to measure a surface of complex geometry, especially steep features on cones and craters, or failure scars. Moreover, these methods are cumbersome owing to the intensive stereopreparation and picture acquisition constraints Villeneuve (2000); Zlotnicki et al. (1990). Recent studies in the field of vision have showed the possibility of using new reconstruction tools in volcanology. The reconstruction approach presented here is based on a multi-view sensor calibration stage followed by a 3D reconstruction process using n independent views. This reconstruction process consists in deforming a generic meshed model of the studied object. At the moment, the method is being tested in laboratory on analogue experiments dealing with volcano deformation.
SENSOR MODELLING AND ESTIMATION OF ITS SPATIAL ATTITUDE
Multi-view Calibration
The initial stage of the method consists of modelling the sensor and positioning external geometry. A mathematical relation is established between 3D point coordinates of a scene and the 2D coordinates of the same points projected in the image.
A multi-view calibration method using a 2D target and based on a photogrammetric approach has been chosen Lavest et al. (1998); Beyer (1992). This method allow to estimate accurately and rapidely the internal and external parameters of the sensor, which will be used in the reconstruction process. In laboratory, the high reconstruction accuracy expected for the studied objects (sand cones) implies to have an accurate estimation of these parameters.
This multi-view calibration approach offers the following main advantages:
31
Cecchi E et al: Multi-view Reconstruction applied to Volcano Monitoring
-     An accurate knowledge of the calibration target is not necessary: its geometry is reestimated during the calibration process.
-     The use of several views allows reliable estimation of distortion parameters.
-     Thanks to its simplicity of practical implementation, the method allows sensor modelling (internal parameters) in the field.
The projection model used for image formation process is referred to Pin-hole optical model (e.g Fig. 1).
Fig. 1. Pin-hole model, image geometry and coordinate system. P is a 3D point in the world referential RW — XYZ. pu v is the projection of P in the image. RC camera referential, with zi = f camera focal length.   (u0,v0): principal point
i.e.
coordinates in the pixel referential. Following notation will be used:
-     {Tx, Ty, Tz)T: translation vector,
-     R:   rotation  matrix  (composed  of 9   elements
(r11r
12'
''33
), defined by the three Euler angles: a rotation around the X axis, b around the Y axis, and g around the Z axis),
-     dx,dy: scale factors of the elementary pixel size,
-     (doxr,doyr) and (doxt,doyt): elements of the radial and tangential optical distortion according to x and
y,
-     (a 1ja 2a3) and (p1,p2): polynomial coefficients modelling the radial and tangential distortion respectively.
During calibration, the Euclidean geometry of the observed target, and some extrinsic and intrinsic parameters are estimated. This approach takes into account the image deformation induced by optical distortion phenomena on the lens surface. Let us define the colinearity equations:
u + ex
v + ey
, + {do xr + doxt) jdx +
^ dx' r-
X+r   Y+r„Z+Tx
,31 X+r 32 Y+r33Z+T z
v0 + {doyr + doyt 33 dy +
...    ...    „.„
\dyr
P(F)
i
(1)
where ex and ey can be seen as measurement errors related to x and y coordinates respectively. Calibrating the sensor is done by estimating the following parameter vector:
F 9+6m+3*n ~    u 0i v 0'a 1 'a 2' a 3 i1 ipi f xif y
X1Y1Z1---,XnYn,Zn,
Tx1Ty1Tz1a1b 1 g 1---,Tym,Tzm,am,bm,g m
where:
(2)
-     uo,vo,a1,a2,a3,p1,p2,fx,fy are the intrinsic parameters (altogether 9).
-     Tx 1,Ty 1Tz1a1b 1 g 1-,Tym,Tzm,am,bm,g m are the 6m extrinsic parameters, m being the picture number,
-    X1Y1Z1---,XnYn,Zn are the 3n target point coordinates, n being the target point number.
The problem, in a least square sense, is to find F that minimizes the error:
S=ååe 2 xij+e 2 yij-
i1j=1
ij            j ij
(3)
As -P(F) and Q(F) are non-linear functions in relation to the parameters of the vector F, the Levenberg-Marquardt routine is used for the criterion S optimization.
This approach is known under the name ”Bundle adjustments”. In order to have an over-determined system, the number of measurements (2.m.n) have to be greater or equal to the number of unknown parameters {9.6m.3n) (internal, external parameters and target point coordinates).
In practice, pictures are taken around a plate composed of retro-reflective circles. The sensor turns around the target so as to create a convergent view bundle to constrain the 3D position of the target points. Circles should cover the entire sensor image to have the best estimation of the distortion coefficients. Fig. 2 shows four pictures among the m views used for the sensor calibration. On each view, 25 circles are visible. A subpixellar detection is done Lavest et al. (1999) and the parameter vector F is calculated for each view of the calibration set.
32
Image Anal Stereol 2002;21:31-36
Fig. 2. Some views of the calibration set.
Sensor Position and Orientation
Fig. 3 shows a set of views around the cone to model, lying on the calibration target. Only a small number of white circles are visible on the pictures. These circles allow positioning the sensor in the model referential.
Fig. 3. Two views among height (2000 x 1030pixels) of the cone picture set and the surrounding target points.
As the internal parameters of the sensor and the (Xi,Yi,Zi) coordinates of the target points are accurately known, the parameter vector F to determine becomes:
F
Tx1Ty 1Tz1a1b 1 g 1
TxqTyqTzqaqb q g q
¦    (4)
However, it is preferable to combine information from the two view sets (calibration + reconstruction view sets) and optimize the whole system. m images containing 25 circles and q images with a few circles ensure the calculation of the sensor position and orientation. So the last 6q extrinsic values of the vector F at the solution give the different sensor attitudes around the cone.
RECONSTRUCTION
The reconstruction stage is tested on scaled models, analogous to volcanic structures. Such analogue models make it possible to study the gravitational deformation of a cone induced by a ductile core. Cones, made up of a sand and plaster mixture and containing a silicone inclusion represent an altered core volcanic edifice. As the silicone has a viscous behaviour, cones deform with time.
The position and orientation of the sensor have been given by the calibration stage ((RiTi) between the camera and the model for each view). So q views around the object are available, all being referenced in relation to the origin of the model referential. Several ways to reconstruct can be envisaged:
-    To erode an initial volume, and only keep the voxels for which the projected luminance is coherent in all views Kutulakos (2000).
-    To make a dense correlation of the points in the different images and triangulate to obtain a dense three-dimensional point cloud Koch etal. (1998).
-    To construct a initial generic model and make an iterative deformation so as the projected texture of each facet of the meshed model be coherent in all views.
The approach using a initial surface model seemed the best adapted for our application. It appeared to be the simplest computational and most direct method, as it deals with a surface.
Initialization of the model
As the studied objects in laboratory are cones, the initial model is a conic surface of variable resolution, and is initialized from the information contained in pictures. The diameter, the height and the position of the cone are estimated. As it is shown in Fig. 4, three points selected by the user at the base of the cone in the image give, by ray tracing, their intersection with the support plane of the model (i.e. P1, P2, P3 points). The circle fitting the best these three points is then estimated. So we can define the vector normal to the circle passing by its centre PS. A fourth point selected at the cone summit in the image gives the cone height by calculating the point on the cone normal that is the nearest to the vector ci*s (e.g. Fig. 4 and Fig. 5).
33
Cecchi E et al: Multi-view Reconstruction applied to Volcano Monitoring
Fig. 4. Initialization of the model with a conic surface.
Fig. 5. One of the real cone view and its projected generic model.
Reconstruction
The meshed model is composed of n facets and p 3-D points. The number of facets can be defined, so a multi-scale approach is possible. We postulates that the model surface to reconstruct can be represented by a function of z = f(x,y) type, which is generally the case for the envisaged applications. This approach allows to limit the number of unknown parameters to estimate. After its initialization, the 3D generic model is deformed to fit the real object. This deformation consists, for a given model resolution, in calculating the z coordinate of each p points of the model in order to have textured areas, resulting from the projection of the facets in the pictures, coherent in all the views (minimization of a correlation criterion).
First, for each facet of the meshed model, a selection of pictures in which the correlation will be done is carried out. This pictures selection is based on (1) the visibility of the facet in the pictures, (2) the angle between view points. At this stage, each facet can be correlated in a minimum of two views
As the objective is to reconstruct as accurate as possible the studied object, a criterion measuring the similarity between this object (the cone) and the generic model is defined. This criterion, based on texture coherence, has to be minimized in order to have the greater adequacy between the object
and the generic model. As we have opted for a surface parametrization of z = f(x,y) type, the (xi,yi) coordinates of each p point will not be modified during the optimization process. Unknown parameters of the vector F are the zi coordinates of each p point of the meshed model (summit facet).
f=   z1?j
,2, ...,^p
(5)
Let us define the optimization criterion:
n   rNi-1    Ni
Cf = å åå Vij -Vik2
i1 j=1k=j+1
(6)
Where:
Ni is the number of views among the picture set where the facet i is visible (Ni > 2).
Vij and Vik are the centred and normalized luminance vectors of a facet i in the two images where the measurement is done. When a luminance vector is centred, its mean equal 0 and normalized means that the standart deviation to 1. Centering and normalizing vectors allows to take into account affine variations of light illumination during picture acquisition, that is variations of the mean luminance value and the grey-level range between pictures. The Vij and Vik vectors are composed of sampling point luminances of a facet. Their size depends on the model resolution and automatically adapts so that the 3D sampling of a triangular facet corresponds to a step of one pixel in the image. Besides, to reduce noise effect, a filtering is done for each measurement of the vectors, the mask filtering size being linked to the model resolution.
Fig. 6. Surface Adaptation. l0l 1 • • • ls are the projected sampling point luminances of a facet i.
To summarize, for a given model resolution, the zi coordinates of the p points are calculated in order that C(F) be minimum (e.g. Fig 6). As the criterion is non linear in relation to its parameters, its minimization is iteratively done, using the Levenberg-Marquardt method.
34
Image Anal Stereol 2002;21:31-36
RESULTS AND PERSPECTIVES             CONCLUSION
In Fig. 7, results from the optimization for two different resolutions of the meshed model (respectively 2048 and 32768 facets) are shown. From an experimental point of view, the reconstructed surface is accurate enough to highlight the main morphological characteristics of the cone, that is to say two nested collapses on one of the flank. This feature can be seen even at a quite low resolution (res.4 = 2048 facets), while at a higher resolution (res.6 > 30000 facets) fine-scale textural irregularities due to grain flow are visible. Increasing the resolution model leads to get progressively closer to the real cone surface. However, the convergence time becomes very important due to the high number of measurements. As we use a non linear optimisation method (Levenberg-Marquardt), we need to give initial conditions close to the solution. This implies that the 3D model initialization is an important stage for the system convergence. However, the multi-scale approach (control of the model resolution) allow to keep a relative flexibility for the initialization. Finally, at the end of the optimization process, calculated camera positions could also be changed to refine the correlation.
Fig. 7. Reconstruction results for different meshed model resolutions. Resolution (res.)4 = 2048 facets, res.6 = 32768 facets.
The videogrammetric reconstruction technique presented here is new to volcanology. A first step will be to validate the method and test its precision at a laboratory scale. The objective will be then to follow the analogue model deformation with time for accurate quantitative analysis. The technique is also going to be applied to real volcanic structures. The transition between laboratory to natural textures will certainly require important modifications of the initial code, especially concerning the orientation and position sensor estimation. Several situations can be envisaged according to the avalaible input data in the field. This is now in development. Howeverthe method promises to provide excellent reconstruction accuracy, a valuable simplicity and flexibility compared to classical stereophotogrammetry techniques currently used. This should make the technique useful and easily applicable for those who need quick and accurate terrain models on rough ground, such as volcanologists engaged in monitoring.
REFERENCES
Beyer H (1992). Geometric and Radiometric Analysis of a CCD-Camera Based Photogrammetric Close-Range System . PhD thesis, Institut fu¨r Geoda¨sie und Photogrammetrie, Nr 51, ETH, Zu¨rich.
Gwinner K, Hauber E, Jaumann R, Neukum G (2000). Highresolution, digital photogrammetric mapping: A tool for earth science. EOS, 81(44).
Koch R, Pollefeys M, Gool LV (1998). Multi viewpoint stereo from uncalibrated video sequences. In Proc. of Fifth European Conference on Computer Vision ECCV98, Freiburg, Germany. 55-71.
Kutulakos K (2000). Approximate n-view stereo. In Proc. of Sixth European Converence on Computer Vision ECCV2000 , Dublin, Ireland. 67-83.
Lavest J, Viala M, Dhome M (1998). Do we really need an accurate calibration pattern to achieve a reliable camera calibration. In Proc. of Fifth European Conference on Computer Vision ECCV98, Freiburg, Germany. 158-74.
Lavest J, Viala M, Dhome M (1999). Quelle precision pour une mire d’e´talonnage? Trait Signal, 16(3):241-54.
Villeneuve N (2000). Apports multi-sources a` une meilleure comprehension de la mise en place des coulees de lave et des risques associes au Piton de la Fournaise. PhD thesis, Institut de Physique du Globe de Paris. Departement des Observatoires. Observatoire volcanologique du Piton de la Fournaise, La Reunion.
35
Cecchi E et al: Multi-view Reconstruction applied to Volcano Monitoring
Zlotnicki J, Ruegg J, Bachelery P, Blum P (1990). Eruptive mechanism on Piton de la Fournaise volcano associated with the December 4, 1983, and January 18,
1984 eruptions from ground deformation monitoring and photogrammetric surveys. J Volcanol Geoth Res, 40:197-217.
36