Department of Science and Technology
Institutionen för teknik och naturvetenskap
Linköpings Universitet
Linköpings Universitet
Examensarbete
LITH-ITN-MT-EX--04/033--SE
Digitizing the Parthenon
using 3D Scanning –
Managing Huge Datasets
Therese Lundgren
LITH-ITN-MT-EX--04/033--SE
Digitizing the Parthenon
using 3D Scanning
-Managing Huge Datasets
Examensarbete utfört i Medieteknik
vid Linköpings Tekniska Högskola, Campus
Norrköping
Therese Lundgren
Handledare: Paul Debevec
Examinator: Mark Ollila
Rapporttyp Report category Examensarbete B-uppsats C-uppsats D-uppsats _ ________________ Språk Language Svenska/Swedish Engelska/English _ ________________ Titel Title
Digitizing the Parthenon using 3D Scanning - Managing Huge Datasets
Författare
Author
Therese Lundgren
Sammanfattning
Abstract
Digitizing objects and environments from real world has become an important part of creating realistic computer graphics. Through the use of structured lighting and laser time-of-flight measurements the capturing of geometric models is now a common process. The result are visualizations where viewers gain new possibilities for both visual and intellectual experiences.
This thesis presents the reconstruction of the Parthenon temple and its environment in Athens, Greece by using a 3D laser-scanning technique.
In order to reconstruct a realistic model using 3D scanning techniques there are various phases in which the acquired datasets have to be processed. The data has to be organized, registered and integrated in addition to pre and post processing. This thesis describes the development of a suitable and efficient data processing pipeline for the given data.
The approach differs from previous scanning projects considering digitizing this large scale object at very high resolution. In particular the issue managing and processing huge datasets is described.
Finally, the processing of the datasets in the different phases and the resulting 3D model of the Parthenon is presented and evaluated.
ISBN
_____________________________________________________ ISRN
LITH-ITN-MT-EX--04/033--SE
_________________________________________________________________
Serietitel och serienummer ISSN
Title of series, numbering ___________________________________
Nyckelord
Keyword
Datum
Date
2004-04-07
URL för elektronisk version
http://www.ep.liu.se/exjobb/itn/2004/mt/033/
Avdelning, Institution
Division, Department
Institutionen för teknik och naturvetenskap Department of Science and Technology
Digitizing the Parthenon using 3D Scanning
-Managing Huge Datasets
Therese Lundgren
USC Institute for Creative Technologies
and Department of Science and Technology, Linköping University
Abstract
Digitizing objects and environments from real world has become an important part of creating realistic computer graphics. Through the use of structured lighting and laser time-of-flight measurements the capturing of geometric models is now a common process. The result are visualizations where viewers gain new possibilities for both visual and intellectual experiences.
This thesis presents the reconstruction of the Parthenon temple and its environment in Athens, Greece by using a 3D laser-scanning technique.
In order to reconstruct a realistic model using 3D scanning techniques there are various phases in which the acquired datasets have to be processed. The data has to be organized, registered and integrated in addition to pre and post processing. This thesis describes the development of a suitable and efficient data processing pipeline for the given data.
The approach differs from previous scanning projects considering digitizing this large scale object at very high resolution. In particular the issue managing and process-ing huge datasets is described.
Finally, the processing of the datasets in the different phases and the resulting 3D model of the Parthenon is presented and evaluated.
Keywords
3D scanning, range image, huge datasets, mesh generation, data processing pipeline, filtering, cultural heritage.
Acknowledgments
This thesis is part of the Master of Science degree in Media Technology and Engineer-ing at the Department of Science and Technology at LinköpEngineer-ing University, Sweden. It presents research and implementation done at the Institute for Creative Technologies (ICT), University of Southern California, Marina del Rey, CA, USA from September 2003 until February 2004.
I would like to express my special thanks to my examiner Dr. Mark Ollila at the Institute of Science and Technology, Linköping University, for making this thesis work possible. Also, many thanks to my supervisor Dr. Paul Debevec at the Institute for Creative Technologies for giving me the opportunity the be a part of a very interesting research project in addition to support and help.
Moreover, thanks to Andrew Gardner and Nathan Yun and all the other staff at the ICT for ideas, help and support, with both research issues and time consuming manual work. Finally, thanks to Paolo Cignoni at CNR in Pisa, Italy, for software support.
Contents
1 Introduction 7 1.1 Motivation . . . 7 1.2 Aim . . . 7 1.3 Problem Description . . . 8 1.3.1 Central Issues . . . 9 1.3.2 Restrictions . . . 9 1.3.3 Requirements . . . 10 1.4 Outline of Thesis . . . 102 Background and Related Work 11 2.1 History of the Parthenon . . . 11
2.2 Data Acquisition . . . 12 2.2.1 Range Images . . . 12 2.3 Registration Phase . . . 13 2.4 Integration Phase . . . 13 2.5 Post-Processing Phase . . . 14 2.5.1 Hole filling . . . 14 2.5.2 Mesh Simplification . . . 14
3 Data Processing Pipeline 15 3.1 Data acquisition and 3D Scanning . . . 15
3.2 Pre-processing . . . 16
3.3 Registration and the Iterative Closest Point Algorithm . . . 17
3.3.1 MeshAlign . . . 18
3.3.2 Selection of Regions . . . 18
3.3.3 Global Registration . . . 19
3.3.4 Transformation . . . 20
3.4 Handling large Datasets . . . 20
3.4.1 Partitioning into Voxels . . . 20
3.4.2 The Triangulation Issue . . . 21
3.4.3 The PLY Polygon File Format . . . 23
3.5 Integration by Volumetric Merging . . . 23
3.5.1 Bounding Box Parameters . . . 24
3.6 Post-Processing Procedures . . . 24
3.6.1 Redundant Triangle Representation . . . 25
3.6.2 Mesh Editing . . . 26
4 Implementation and Software 28
4.1 Filtering . . . 28
4.2 Registration . . . 29
4.3 Voxelization . . . 31
4.3.1 Adjustment of the Coordinate System . . . 31
4.3.2 Partitioning . . . 32
4.3.3 Triangulation . . . 33
4.4 Integration . . . 35
4.5 Post-Processing . . . 35
4.5.1 Zipper used on the Voxels . . . 35
4.5.2 Mesh Editing in GSI Studio . . . 36
4.6 Automatization of the Pipeline Process . . . 37
4.7 The Surrounding Environment . . . 37
5 Results 39 5.1 Processing Pipeline . . . 39 5.2 Reconstructed Model . . . 40 5.2.1 Partial Results . . . 41 5.2.2 Final Results . . . 42 5.3 Project Integration . . . 43 6 Discussion 46 6.1 Process discussion . . . 46 6.2 Conclusion . . . 47 6.3 Future Work . . . 47 A List of Abbreviations 52
List of Figures
1.1 Registered point data (left), pieces of data in the registration phase (top
right) and integrated geometry (bottom right). . . . 9
2.1 The west facade of the Parthenon in Athens, Greece. . . 11
3.1 The Quantapoint laser scanner. . . 15
3.2 The acquired scan data, bottom figure visualizes the depth differences. 16 3.3 Gaussian blur applied n times. . . . 17
3.4 Categories of the range scans, 1 (top), 2 (middle) and 3 (bottom). . . . 19
3.5 Pair alignment; less accumulated error in the right than in the left graph. 20 3.6 Voxel overlap without seams, the right figure with one voxel shaded to emphasize the overlap. . . 21
3.7 A range image and its (i,j) coordinate system. . . . 22
3.8 Thresholds control the triangulation. . . 22
3.9 The result when triangulating one single range image. . . 23
3.10 Grid after merging; original (left) and using new parameters (right). . 25
3.11 Boundary edge case. . . 25
3.12 Black areas will create holes in the geometry . . . 26
3.13 Decimation of the surrounding environment; 5, 10 and 20 percent re-spectively of the original amount of polygons. . . 27
4.1 Pre-processing of the range images. . . 28
4.2 Selection of regions (left) and corresponding new PFM:s (right). . . . 29
4.3 The manual alignment method in MeshAlign . . . . 30
4.4 Registration phase. . . 30
4.5 Voxelization of the facades; east (top left), north (top right), west (bot-tom left) and south (bot(bot-tom right). . . . 31
4.6 Voxelization phase. . . 33
4.7 Triangulation scheme. . . 33
4.8 The voxel array. . . 34
4.9 Integration phase. . . 35
4.10 Removal of redundant triangles. . . 35
4.11 Hole filling; before (left) and after (right). . . . 36
4.12 Mesh editing and decimation. . . 36
4.13 Graphical user interface; voxelization (left) and merging (right). . . . 37
4.14 The reconstructed surrounding environment of the Parthenon. . . 38
5.1 Flowchart of the 3D data processing pipeline for producing the model. 39 5.2 Flowchart of the post-processing phase. . . 40
5.3 The east facade assembled from four range images (left) compared to
the same part from the completely processed final model (right). . . . 41
5.4 A rendering of the whole model before completion of the post-processing phase. . . 42
5.5 Scanner location map (seen from above). . . 43
5.6 The final model and environment. The white area indicates a voxel and the blue regions show its neighboring voxels. . . 43
5.7 A real photograph (left) compared to a rendering under the same illu-mination (right). . . . 44
5.8 Rendering from the larger research project; the east facade. . . 44
5.9 Rendering from the larger research project; the west facade. . . 45
List of Tables
4.1 Number of points to be voxelized. . . 32
4.2 Connectivity thresholds. . . 34
5.1 Description of the phases in the processing pipeline. . . 40
A.1 List of Abbreviations . . . 52
Chapter 1
Introduction
This chapter gives an introduction to this thesis work by presenting the motivation to the work, a problem description and finally an overview of the structure of this report.
1.1
Motivation
Digitizing objects and environments from real world has become an important part of creating realistic computer graphics. Capturing geometric models through the use of structured lighting, laser triangulation, and laser time-of-flight measurements is now a common process. The scene digitized in this thesis is the Parthenon and its environ-ment in Athens, Greece, done in collaboration with the ongoing Acropolis Restoration project.
What makes this thesis and the scanning of the Parthenon different from previous larger scanning projects [18] [20] is mainly the data acquisition of a structure having the size of the Parthenon (75m × 35m) at this high resolution resulting in huge datasets. The reproduction of the model can not simply follow the traditional data processing pipeline. Instead the central issue is to develop a data processing pipeline suitable for the acquired data consisting of both traditional and new phases. By using the developed pipeline an accurate 3D model of the Parthenon is to be constructed.
Moreover, this thesis probably differ from traditional ones by being a part of a larger research project1with several people involved, and where the digitizing of the Parthenon is a very significant component. The final model have to fulfill the require-ments for the larger research project in order to contribute to good results. The continu-ous integration into the larger project is important as well as meeting internal deadlines during the time for this thesis.
1.2
Aim
The aim of this thesis is firstly to find an efficient, accurate and not too time consuming data processing pipeline for the acquired 3D scan data. This by both implementing own algorithms and by using existing software.
Secondly, the aim is to assemble all the acquired data into an accurate geometry model of the Parthenon by using the developed data processing pipeline. The recon-struction of the final geometry model is further divided into three partial goals:
• A model of the east facade. • A model of the west facade.
• A complete 3D representation of the whole structure.
Finally, the purpose of the complete model is to be integrated into the larger re-search project at the ICT, with the aim to publish a scientific rere-search paper as well as film for the electronic theater at the SIGGRAPH 20042conference.
1.3
Problem Description
In order to produce the final model of the Parthenon a suitable processing pipeline has to be developed. The problem of finding an efficient 3D data processing pipeline involves following phases (see also figure 1.1):
Data acquisition First the data has to be automatically collected by the use of a
scan-ning sensor. Depending on how the data is acquired the need for certain phases varies.
Pre-processing Pre-processing of the acquired data involves converting the data to an
appropriate format, size and quality. There is a major need for pre-processing of the acquired data in terms of filtering since the laser scan measurements include noise.
Registration In the registration phase pieces of data is registered to each other in order
to form a complete and continuous data set. This means transformations between the coordinate systems of the different pieces of data.
Integration The integration phase addresses the problem of combining and merging
the different pieces of data from the registration phase into one consistent piece.
Post-processing Data acquired from 3D scanning is not always complete, e.g.
oc-cluded geometry create holes, and the resulting mesh might include too many polygons. Therefore a need for various post-processing methods are necessary. The approaches in the phases have to be able to solve the processing of the acquired large datasets. That implies adjustments of both datasets and implementations.
Figure 1.1: Registered point data (left), pieces of data in the registration phase (top
right) and integrated geometry (bottom right).
1.3.1
Central Issues
The central issues in the various phases presented in this thesis work involves follow-ing:
• The large amount of noise in the acquired data, i.e. the range scans have to be
decreased. Depending on the strength of the return, the depth accuracy varies between 0.5cm and 3cm.
• To register the range scans without a pre-alignment done on site.
• To merge a large scale model with the size of the Parthenon (75m×35m×20m)
at high resolution as well as its surrounding terrain (200m×200m) at low reso-lution.
• How to efficiently handle the large number of scans and their huge amount of
data. One single full resolution range scan contains 58 million points and there are in total 120 scans.
• The processing of the data pipeline has to integrate into other processes
simulta-neously and the file formats have to be compatible since this thesis is a part of a larger project.
• The processing pipeline have to be partially automatized to be used efficiently
without supervision, e.g. processes being run over night.
1.3.2
Restrictions
Restrictions in the data acquisition have impact on this thesis work and its possible results. When performing the 3D scanning of the Parthenon temple the focus was on capturing the two facades of the east and west side of the temple. Scaffolding and equipment around the structure and in the middle section of the temple made the scanning more difficult, especially on the north and south sides.
Additionally, most photos were taken of the east facade. Since photos are an es-sential requirement for further work in the larger research project, this affected the ordering of the partial goals for the reconstruction of the final model (see section 1.2).
1.3.3
Requirements
In this thesis work there are requirements in the data processing which have to be fulfilled. These involves quality, resolution and efficiency.
Firstly, due to the varying strength of the return of the scanner’s laser, noise was created in the range measurements. The existing noise have to be filtered intelligently in order to obtain smooth geometry without loosing details.
Since the purpose of the reconstructed model of the Parthenon is to be integrated into the larger research project the requirements of that project apply. Texturing the model is not part of this thesis but one of the central issues in the larger research project. To perform a photo realistic texturing, i.e. to estimate the reflectance properties for the surface of the reconstructed model, the geometry has to be accurate enough. Thus, the resolution of the model has to be high, with the geometry in a rendering of the recon-structed model being comparable to a real photograph taken from the same location as the rendering. Furthermore, the final model has to be available in various resolu-tions. This since a complete full resolution model will not be possible to load into the memory of an ordinary workstation at once.
For some of the software used in the process there are restrictions regarding a max-imum amount of polygons allowed in the files to be processed. This requires that the original data cannot simply be processed in the existing form.
There are also requirements regarding the time aspect. The processing times in the various phases of the pipeline have to be optimized since there are many time consuming processes and much manual work involved with data acquired from 3D scanning.
1.4
Outline of Thesis
This report proceeds as follows. Chapter two presents the background and previous work in the 3D scanning area. Chapter three describes the methods and algorithms in the developed data processing pipeline. In chapter four the implementation and the software are discussed. Chapter five presents the results and evaluates them. Finally, chapter six summarizes and discusses the report and examples of future work are pre-sented.
Target group
The reader of this thesis preferably has a good knowledge in computer graphics and a general knowledge about math.
Chapter 2
Background and Related Work
This chapter will introduce some important theories and previous work in the area of research within digitizing real world objects as well as the history of the Parthenon. It also explains the necessary phases when processing data acquired from 3D scanning.
2.1
History of the Parthenon
The Parthenon was erected during the 5th Century B.C., replacing an earlier project that was ransacked by invaders before it was finished. The temple was dedicated to the goddess Athena and came to be her most magnificent monumental offering [17]. Sculptures of historical events and Greek myths decorates the Parthenon, not all of them there today. Strategies for digitally capturing a large number of sculptures have recently been described [25].
Over the centuries, the Parthenon has changed greatly [26]. During the 3rd century B.C., the temple was destroyed by fire and then minimally restored. In the 6th century A.D., the temple was turned into a cathedral and high windows were carved through the walls.
Later on when Athens was a small city in the Turkish Empire the Parthenon was turned into a mosque. In 1687 a Turkish garrison used it for storing gunpowder when it was bombed by the Venetian fleet. A cannonball ignited the gunpowder and significant portions of the building were destroyed in the catastrophic explosion [17].
2.2
Data Acquisition
Various people, such as Levoy et al. [18], Ikeuchi [16] and Rushmeier [23] have been involved in 3D scanning and have shown that accurate and detailed geometric models can be acquired of real-world objects using these techniques. Various non contact rangefinding techniques can be used and different scanning devices exist, including both academic prototypes and commercial systems. A common classification divides them into passive and active systems [4].
The passive systems includes for example photogrammetry1which is based on the reconstruction of the geometry from a silhouette approach and returns a nearly com-plete model of the object but suffers from scarce accuracy. The active systems project a laser or a structured light pattern on the object and measure the geometry of the points they hit by for example a time-of-flight technique. In this technique the time it takes for the laser or light is measured and converted into a distance. A problem in the active systems is to capture a complete model without any missing data.
A novel technique combining passive and active systems has recently been pre-sented by Dias and Vaz [11]. The technique demonstrates the possibility to register and combine the advantages from both range and intensity based data.
One of the most recent scanning project is the Digital Buddha Project at the Uni-versity of Tokyo [20]. As with all scanning projects the purpose is to accurate measure 3D coordinates of a huge amount of points over a short period. The points will then be used to reconstruct and model the scene.
2.2.1
Range Images
The data obtained in the data acquisition is a set of range images or point clouds con-sisting of a collection of single (x,y,z) points, which represent the distance between the scanning sensor and the geometry. A range image is a two-dimensional i×j grid of these range points sampled from the part of the surface of the target object that the scanning sensor was able to reach during each session.
Range images are usually very large, storing millions of points, and the file size of the datasets tend to be huge. One significant challenge that the 3D scanning project
the Digital Michelangelo Project [18] encountered was the handling of large datasets.
Their model of Michelangelo’s David was two orders of magnitude larger than any ex-isting scanned model. For that purpose the PLY polygon file format [30] was invented in 1998 in order to try to create a single simple file format that is flexible enough for
1Photogrammetry; a science or art of obtaining reliable measurements or information from photographs
future needs and usable for a wide range of models. The format stores graphical ob-jects described as a collection of polygons, and it has both an ascii representation and a binary version for compact storage and rapid saving and loading (see appendix B).
2.3
Registration Phase
The first step in creating a single model from multiple range images is the registration or aligning phase. Normally the range image geometry is relative to the current sensor location, a local coordinate system, and therefore has to be transformed into a common coordinate space where all the range images are well aligned in 3D space, a global
coordinate system. After the alignment, the 2D surfaces in the range images that are
partially overlapping also will be geometrically overlapping. The aligning procedure is usually a manual and time consuming phase in the 3D scanning pipeline.
The research area of aligning range images involves the Iterative Closest Point (ICP) algorithm proposed by Besl and McKay [2] which have showed to be the stan-dard solution to the alignment problem. Given two 3D point sets, P and Q, the task is to find the rigid body transformation, which means the rotations and translations, that brings P into the best possible alignment with Q. The ICP algorithm requires an initial estimate for which points that should have the same 3D position. Different approaches for efficient [24], stable [13] and trimmed [5] versions of the ICP algorithm have been proposed.
In a recent publication, Allen et al. [1] focus on an automatically algorithm to assemble data acquired from a time-of-flight laser scanner. One of their test beds for their model-building methods is the Cathedral Saint-Pierre in Beauvais, France. The range images obtained there provides up to one million points, the range images of the Parthenon contains originally up to 58 million points each.
2.4
Integration Phase
Integration or merging follows when multiple range images are aligned in order to build a consistent non-redundant mesh. Within this research area many reconstruction methods have been proposed during the last few years. They are being divided into
volumetric and non-volumetric methods. Volumetric approaches immerse the range
map in a discrete voxel space and reduce the integration problem to that for each voxel cell reconstruct that part of the surface included in that particular cell. Most of these approaches requires discrete distance functions. Curless and Levoy [8] presented a volumetric method consisting of a cumulative weighted signed distance function where they added one range image at a time.
In a paper of Turk and Levoy [27] a non-volumetric method called zippering is described. The method is incremental and works with pairs of meshes. In the first step all redundant triangles are discarded by simultaneously eating back the boundaries of each mesh that lie directly on top of the other mesh. This is followed by the zippering where the triangles of one mesh are clipped to the other mesh and the vertices on the boundary are shared.
A volumetric merging method which uses a Marching Intersection (MI) algorithm derived from the Marching Cubes2 [19] method is described, implemented and im-proved in publications from the Visual Computing Group in Pisa [21] [22]. The
gorithm is based on dividing the volume to be merged into a discrete voxel grid and individually extract a surface from each voxel .
2.5
Post-Processing Phase
Depending on the methods used for aligning and merging in the 3D scanning pipeline there are different needs for post-processing. The post-processing phase involves sev-eral research areas and the main areas in this thesis are hole filling and mesh simplifi-cation.
2.5.1
Hole filling
Models created from 3D scanning tend to include holes due to either occluded geom-etry or dark areas of the surface which the scanner cannot register. Hole filling can be performed either as post-processing operation or integrated in the integration phase.
Cureless and Levoy [8] integrates hole filling into their volumetric merging algo-rithm which is based on meshes. Another mesh based reconstruction that includes hole filling have been proposed [27] [29] treating each range image as an ordered 2D array, which is triangulated to form a polygonal mesh. Hoppe et al. [15] describes a regis-tration method based on all range images treated as a large point cloud. With no initial connectivity between the points holes are filled during the surface reconstruction.
One widely used approach to fill simple holes is to triangulate the surrounding components of the boundary by creating a surface with the topology of a disk. For holes with more complicated topology this may result in self-intersecting geometry and the problem is harder to solve. Davis et al. [9] address the hole filling problem by using volumetric diffusion. The volume is assigned a distance function with the surface having zero value and with positive values on the inside and negative values on the outside. The goal is to extend the distance function over the holes so that finally the function is zero over the entire surface implying the desired hole-free geometry.
2.5.2
Mesh Simplification
Mesh simplification or decimation is mandatory when managing meshes produced with 3D scanning methods. The high sampling rate, up to 10 points per square millimeter of the scanner devices produces meshes consisting of a huge amount of polygons, about 20M-100M faces per mesh are common [4]. To be able to import the meshes into other applications the number of faces has to be significantly reduced.
A comparative publication is presented by Cignoni et al. [6] introducing various approaches of mesh simplification algorithms. Different implementations are tested and the computational complexity and approximation accuracy of the resulting meshes are compared.
Previous work on simplifying meshes while preserving surface properties such as texture coordinates and color is described in a paper by Garland and Heckbert [12]. Their mesh simplification algorithm uses iterative contractions of vertex pairs to sim-plify the geometry. The error between old and new vertices is measured using quadric matrices, which is a symmetric 4×4 matrix associated with each vertex. For a given vertex the sum of squared distances to an associated set of planes is measured. This results in a fast algorithm with high quality approximations which both preserves ge-ometry and topology.
Chapter 3
Data Processing Pipeline
This chapter describes the different phases in the chosen 3D data processing pipeline in order to produce the final model of the Parthenon.
3.1
Data acquisition and 3D Scanning
During the summer 2003 the Graphics Lab at the ICT acquired in total 120 range images in and around the Parthenon site over five days. The chosen time-of-flight panoramic range scanner is manufactured by Quantapoint, Inc1and uses a 950nm in-frared laser measurement component [14]. It can be classified as an active system, as described in section 2.2. In high resolution mode, the scanner acquires scans of 18,000×3200 (x,y,z) points in eight minutes, with a maximum scanning range of 40 meters and a field of view of 360 degrees horizontal and 74.5 degrees vertical. The main reason for choosing this scanner is the very high speed of measuring a large amount of points at extremely high resolution. Some scans from within the structure or at a very short distance were scanned in low-resolution, acquiring one-quarter the number of points.
Figure 3.1: The Quantapoint laser scanner.
The laser system scans a volume of space where millions of digital measurements define every shape and dimension. The massive dataset produced by the scanner makes
it possible to obtain models with photo quality, which, as mentioned in section 1.3.3 is necessary for the purpose of the model. The output data from the scanner is in the Quantapoint format (QPG). Besides the range image a 16-bit monochrome image of the infrared intensity returned to the sensor for each measurement is produced by the scanner.
During this scanning procedure, depending on the strength of the return, the depth accuracy varied between 0.5cm and 3cm. This varying depth accuracy resulted in very noticeably noise in the geometry which had to be reduced, see section 3.2.
Figure 3.2: The acquired scan data, bottom figure visualizes the depth differences. Scans from outside the structure were initially aligned during the scanning process through the use of checkerboard fiducial markers placed within the scene. The pre-alignment procedure was a relatively time consuming process and was therefore not possible to perform for all the scanning sequences. Additionally, the scanner was con-strained to scan horizontally, no tilting possible, and four markers were necessary for each range image.
The pre-alignment proceeds as follows. The first scanning session in the sequence is set as the reference scan. For each new scan in the sequence four markers have to be the same as for the previous scanning session. The heads of the markers are allowed to be rotated between the sessions to face the scanning sensor.
The registration on site consisted of marking corresponding fiducial markers in the available interface from Quantapoint directly after each session. Thereafter the scan is transformed into the local coordinate space for the previous scan. In this way the first step in the registration phase is achieved.
3.2
Pre-processing
Initially, in the pre-processing phase, the obtained point clouds in the Quantapoint format are converted into range images in the portable float map2format (PFM). Each pixel in a PFM range image has a r, g and b value which is equal to that corresponding point’s (x,y,z) location in 3D space. Thus, hereafter when referring to the range images’ (r,g,b) values the notation (x,y,z) will be used.
Due to the varying strength of the return of the laser, noise was created in the range measurements. Therefore, in the pre-processing phase the noise has to be filtered in order to produce smooth geometry. For determining the amount of noise and the
accuracy of the range images they are triangulated (see section 3.4.2 for a description how it is implemented).
A few attempts using various smoothing algorithms were tested and the one lever-aging the best result for the test geometry, considering smoothing originally flat sur-faces is the Gaussian blur filter.
A filter kernel consisting of a 3×3 matrix with weights as in equation 3.1 is applied on the range images. To sequentially applying small kernels instead of larger ones reduces the execution time [28]. In order to preserve edges and small details a threshold is used and added to the Gaussian blur filter. If the difference between a pixel and its neighbor (either the x, y, or z value) is greater than this threshold the weight is set to zero implying that the neighbor will not affect the current value of the pixel. The threshold helps preserving edges and details in the geometry. Following, for each value the weight is multiplied and then the average value over the 3x3 kernel replaces the current value.
Filter kernel and its weights =
12 24 12
1 2 1
(3.1)
Figure 3.3: Gaussian blur applied n times.
The filtering can be applied multiple times over the same image resulting in differ-ent levels of smoothness. The best result for the range images in this thesis showed to be applying the filter five times. A comparison is shown in figure 3.3. As visible in the figure the filtering also results in small holes being filled. This advantage is obvious already after applying the filter once.
Since the filtering also implies that original information is lost a subsampling algo-rithm could be implemented after the filtering without any noticeable visual difference. Subsampling is performed by, in the image plane ignoring every other pixel in both the
i and j dimension. Running the algorithm once decreases the file size by four without
jeopardizing the accuracy and makes the pipeline faster and the data easier to handle.
3.3
Registration and the Iterative Closest Point
Algo-rithm
The registration of the range images is solved in a partially manual and a partially automatic manner. The aligning process is one of the most time consuming phases in the whole 3D scanning pipeline since it demands user assistance.
3.3.1
MeshAlign
To register the range images the Italian software MeshAlign v.2 developed by the Visual Computing Group at CNR-Pisa [4] is partially used. MeshAlign provides an interface to estimate an initial alignment and has an aligning procedure based on the ICP algo-rithm [2] introduced in section 2.3.
The ICP algorithm includes three basic steps: (P and Q are two 3D point sets.)
1. Pair each point of P to the closest point in Q.
2. Compute the rigid motion that minimizes the mean square error3(MSE) between the paired points.
3. Apply the rigid motion to P and update the MSE.
The three steps are iterated and the iterations have been proven to converge in terms of the MSE. Since two range images never include the same data but are overlapping the ICP algorithm cannot be run using all their points. Instead a group of paired points is used, i.e. the sample points.
The registration algorithm in MeshAlign is implemented by using the standard semi-automatic range map alignment approach [4] which starts with following two procedures.
Initial Estimate: The first step is to pairwise locate all range images in a single
com-mon global coordinate system. Each pair consists of two adjacent and over-lapping range images. This will provide the algorithm with a first rough initial estimate for which points should have the same position after the alignment.
Refinement: After the first step, the initial alignment is refined using the iterative ICP
algorithm which minimizes the alignment error between each pair of meshes. The initial estimate procedure is user dependent and is performed by either marking point pairs or by interactive manual manipulation. For the pre-aligned range images the initial estimate is already solved for, resulting in only the refinement step being computed. The second step, the refinement, is the automatic part of the semi-automatic alignment approach.
3.3.2
Selection of Regions
MeshAlign cannot handle meshes having the size of the meshes obtained from the
Parthenon range images, it has a limit of a file size about 30Mb, corresponding to about 1.5M polygons. The meshes representing the Parthenon are, after the filtering and included size reduction of a size up to 140Mb each.
In order to be able to take advantage of the already implemented ICP algorithm in
MeshAlign, a modified aligning method is developed. The method is based on selecting
and cropping out small regions with obvious features in the image plane. The features selected from one range image shall easily be found in the corresponding range image it will be aligned to. This is a manual and user assisted procedure where four regions are selected (see figure 4.2).
Thereafter, the new range image is triangulated into a mesh which can be imported into MeshAlign. The selection of regions is performed carefully since the aim is to choose easily recognisably features and to select regions as far away from each other in terms of position in 3D space. Thus, obtaining a rigid motion as robust and correct as possible from the ICP algorithm.
3.3.3
Global Registration
Organizing the global registration, which is what range images should be aligned to each other and treated as pairs, is a complex task. How complex this task actually is showed to be very dependent on how the scanning on the site is performed. In the Parthenon case there are basically three categories of range images (see also figure 3.4). 1. The category with the most data consists of range images acquired from the scanner when located at about 40 meters from the structure. The scanner was placed at this distance all around the structure, digitizing all its facades. Pre-aligning of the range images was performed.
2. Next group consists of data produced from the scanner when being tilted and facing upwards. This in order to scan the ceiling and the backsides of columns. Pre-aligning was not possible to perform.
3. The last category are range images acquired when placing the scanner very close to certain parts of interest of the Parthenon. Due to time restrictions pre-aligning was only performed for some of the scanning sequences.
Figure 3.4: Categories of the range scans, 1 (top), 2 (middle) and 3 (bottom). The global registration is already solved for considering the pre-aligned range im-ages. Pre-aligned datasets have one range image which is the identity range image, and its coordinate system will be considered the global coordinate system. The identity range images will have no translation in the (x,y,z) directions. The other range images in the sequence have their translations relatively to the identity range image, stored as
vectors in a separate file. Rotation of the scanning sensor was not allowed within the se-quence, consequently only the translation is being stored. Having pre-aligned datasets makes the global registration much more efficient since this otherwise is a very time consuming phase.
When choosing pair of meshes for the not pre-aligned range images there are some things to consider. Since the pairwise registration sequentially performs the error min-imization, the error tends to accumulate [4]. This may result in artifacts after a number of pairwise steps. Thus, if possible, a range image in the middle of a sequence is cho-sen as the identity range image in order to decrease the number of pairwise steps where the error might accumulate. Figure 3.5 shows how the sequence of pairwise alignment can be decreased by choosing another identity range image.
ID
ID
Figure 3.5: Pair alignment; less accumulated error in the right than in the left graph. The result from the global registration is a 4×4 matrix describing the translations and the rotations, i.e. the rigid motion of the aligned mesh, stored in a alignment file (ALN file).
3.3.4
Transformation
Next step is to actually apply the matrices in the corresponding ALN files on the range images. This proceeds in the image plane simply by multiplying every non-zero (x,y,z) value with the matrix. Noteworthy is that the pre-alignment showed to leverage very accurate results, the additional transformations for these range images were very small in most cases.
3.4
Handling large Datasets
Most previous 3D scanning projects [18] [2] [25] have involved statues or other smaller things in the range of a couple of meters. The Parthenon temple measures a surface of 75 meters × 35 meters and the surrounding terrain about 40000 square meters. Having 120 range images, each of a size up to 18000 × 3200 pixels, and each storing three float point numbers the initial data set is extremely large. Following section describe how the problem of handling large datasets is solved.
3.4.1
Partitioning into Voxels
Partitioning into voxels means to divide the data into smaller three-dimensional vol-umes, e.g. voxels, which hereafter will be referred to as voxelization. The main reason for this is that the merging phase cannot process too large meshes. Another reason is that the final model has to be manageable in other implementations outside this thesis.
The aim is to sequentially merge the voxels separately, resulting in a consistent model when viewing them all together afterwards. Different approaches for voxeliza-tion were investigated, e.g. octree representavoxeliza-tion4used in [10] or a regular grid. In the chosen approach a regular grid is most suitable to use. The space containing the Parthenon is divided into equally sized voxels by splitting the x, y and z axis an appro-priate number of times. The determination of a required voxel size in order to obtain a sufficient resolution when merging was carried out by testing various cases of different size.
To avoid holes being created during the merging process, the voxels are created partially overlapping but perfectly connected (see figure 3.6). Thus, neighboring vox-els share several points in the boundary regions. These points are triangulated in all neighboring voxels preventing holes to occur. A resulting fallout is redundant poly-gons being created, section 3.6.1 will describe this in detail.
Figure 3.6: Voxel overlap without seams, the right figure with one voxel shaded to emphasize the overlap.
The result of the voxelization for each range images is, a (large) number of equally sized, partially overlapping voxels all containing a different amount of polygons. When a range image has been voxelized the triangulation of its points follows.
3.4.2
The Triangulation Issue
The aligned i×j range images are in the portable float format where each pixel in the image has a (x,y,z) value corresponding to the pixels’ location in 3D space. The trian-gulation issue is solved by triangulating a range image based on the (i,j) coordinates in the image plane (see figure 3.7). Similar approaches are presented in [27] [29].
By connecting neighbors fulfilling requirements based on different thresholds im-plemented in a connectivity function triangles are created (see section 4.3.3 for further details). Connecting neighbors means that four points in adjacent rows and columns in the range image possibly creates zero, one or two triangles. Figure 3.8 shows how the thresholds control and avoid connecting portions of the surface that are separated by depth discontinuities.
4Octree; subdivision of space recursively into eight cubes, each of which may be subdivided into eight
i
j
Figure 3.7: A range image and its (i,j) coordinate system.
distance > threshold, no connection made scanned geometry
scanner location
Figure 3.8: Thresholds control the triangulation.
Triangulating a 2D range image based on its (i,j) coordinates is elementary since the width and height of the image are known and can be used for finding the right neighbors to compute connectivity on. The implementation of the voxelization implied that only the non-zero (i,j) coordinates of the points are saved for each voxel. The (i,j) coordinates are stored together with the voxel index in a two-dimensional voxel array containing all voxels.
In order to use the same triangulation method as for a complete range image the voxel array first is being searched for its minimum and maximum i-values to obtain the corresponding width of the voxel in the image plane. Knowing the width makes it possible for the triangulation algorithm to still be based on the ordering of the points in the original range image. The result of triangulating one single range image is showed in figure 3.9 and implementation details are described in section 4.3.3. When all voxels have been triangulated they together create one large mesh representing the input range image.
Figure 3.9: The result when triangulating one single range image.
3.4.3
The PLY Polygon File Format
One of the main issues with the data pipeline when handling huge amount of large files is memory and speed. Moreover, the files have to be in a format that other software further on in the pipeline are able to read. Previous research [18] has shown that the PLY polygon file format introduced in section 2.2.1 is very suitable for this purpose. Therefore, the range images are triangulated into binary meshes in this format. A comparison between the ascii file format OBJ and the binary PLY format showed that the last one requires almost half the space, e.g. a piece of geometry represented as an OBJ file demands 14390kB compared to 7600kB as a PLY file.
3.5
Integration by Volumetric Merging
For integration this thesis will describe a volumetric merging method which uses a
Marching Intersection (MI) algorithm [22] derived from the Marching Cubes5method
[19] introduced in section 2.4.
A standard volumetric method assigns a discrete voxel grid consisting of tiny voxels to the input meshes. These voxels are not to confuse with the larger voxels being described in the voxelization algorithm in section 3.4.1. Each tiny voxel is assigned a distance function, D(x) and a weight W(x). The weight for each voxel is normally dependent on the direction of vertex normals and viewing direction. W(x) is increasing the more similar the different directions are. The integration of the different meshes into one consistent mesh is done by extracting an isosurface at D(x)=0.
The MI algorithm [21] [22] was developed to have good performance but with less memory occupancy. This is obtained by only storing the intersections of the meshes with the grid instead of a 3D distance field over the whole surface as in the standard volumetric method described above.
The algorithm is implemented in the tool PlyMC [4] which is used for the merging of the Parthenon geometry. There are different parameters available, the most
tant are voxel size which determines the resolution of the finally merged mesh. The smaller voxel size the more time and memory consuming merging process why it is important to determine the largest size without loosing resolution. To avoid loosing details in the current meshes the size of the tiny voxels has to be 1.2cm × 1.2cm. The most important steps in PlyMC are the following:
For each mesh to be produced:
1. Compute the bounding box for the different input meshes.
2. In the bounding box, make a voxel grid of tiny voxels with the chosen size. 3. Perform the volumetric merging by using the MI algorithm.
4. Create and save an output mesh.
3.5.1
Bounding Box Parameters
The larger voxels are obtained by regularly partition the Parthenon space (described in section 3.4.1) and thus have the same size but contain different geometry, which not always extends all the way to the boundaries of the voxels. If using the ordinary PlyMC the computed bounding boxes for the large voxels will not have the same size. When adjacent voxels are processed in PlyMC and divided into the tinier voxels, the resulting grid intersections in the overlap regions do not correspond to each other. Consequently, the same vertex but in overlapping adjacent voxels are not likely to have the same 3D position in both voxels after merging since the voxel grids might be slightly shifted.
The MI algorithm, as mentioned earlier, is based on the intersections in the voxel grid. It stores the intersections of the geometry and the grid and creates new vertices at every intersection. In order to have the final voxels connecting without seams the voxel grid used for merging also have to line up perfectly, i.e. the bounding boxes for the voxels have to be of the same size and divisible with the tiny voxel size. This is obtained by changing PlyMC not to compute the bounding box itself but to use given input parameters for it.
Providing PlyMC with the correct bounding box parameters based on the position of the large voxels in the Parthenon grid, the voxel grid used in the merging phase also connects seamlessly. This results in vertices in the overlapping regions having the same positions before as well as after merging (figure 3.10). An exception, though, is if the vertices are too close to a boundary edge of a voxel. Close to an edge the MI algorithm gives different results for the same vertex since it does not have the same neighboring data as input (figure 3.11).
The existing redundant triangles in the overlapping regions are considered and pro-cessed in the next step in the pipeline (see section 3.6.1). Moreover, the changes in PlyMC and the input bounding box parameters are explained further in section 4.4.
3.6
Post-Processing Procedures
The post-processing procedures of the final voxels in the pipeline involve removing redundant polygons, mesh editing and mesh simplification. First redundant triangles are removed followed by clean up and hole filling. Finally, mesh simplification is performed.
Figure 3.10: Grid after merging; original (left) and using new parameters (right).
Figure 3.11: Boundary edge case.
3.6.1
Redundant Triangle Representation
Turk and Levoy [27] have developed a program called zipper which central algorithm removes overlapping portions of meshes and then clips one mesh against another. To solve the problem with the redundant triangle representation in the overlap regions in the voxels zipper is modified.
First, between all adjacent voxels a boundary halfway between the two correspond-ing meshes, A and B, is determined. Next, redundant triangles are discarded by tessel-lating all triangles on mesh A. Triangles before and on the boundary are preserved, the ones positioned past the boundary are discarded. Continuing with mesh B, all tri-angles in its overlapping region that are on the boundary or redundant with the ones preserved in mesh A are removed. These operations result in all redundant triangles being discarded.
In some cases this procedure for removing redundant triangles will leave tiny, hardly visible, holes where the two meshes meet. Therefore, in some places addi-tional vertices are added and triangulated. Triangulation is easily solved since the two vertices next to the new vertex are known as well as the one these two vertices had in common in the already existing triangle. These four vertices construct two new triangles that make the two meshes join smoothly.
3.6.2
Mesh Editing
A clean up procedure is necessary since there during the scanning were scaffolding inside and at some locations around the Parthenon. Digitizing items like scaffolding is difficult since the geometry is very fine. Additionally, the scaffolding inside the Parthenon were not completely captured when most scans are from the outside. There-fore and also because of esthetic reasons all scaffolding geometry are removed.
The following necessary mesh editing step is hole filling. Models acquired from 3D scanning usually include holes since the scanning sensor usually cannot observe the entire surface of the object due to occlusion. In addition, there are black parts on some locations of the Parthenon surface (see figure 3.12) that could not be registered, as well as plants growing on the structure preventing the sensor from registering these regions.
Figure 3.12: Black areas will create holes in the geometry
The mesh editing is done in a semi-automatic process by using the commercial GSI
Studio software from Geometry System Inc.6. It is claimed to be a complete software
for creating polygon meshes from range data but that was not the case. Firstly, the datasets of the Parthenon are too large and could not be loaded into the memory of the software. Secondly, their registration implementation was not fully developed. It was not possible to manually change the view, i.e. the 3D position of the camera of a mesh in order to align it to another mesh.
Fortunately the data in the developed processing pipeline has been voxelized into smaller pieces resulting in GSI Studio being convenient for mesh editing purposes. Important is also that loading and saving of PLY files are possible. The approach in GSI Studio uses the method mentioned in section 2.5.1 where a hole is being filled by creating a surface with the topology of a simple disk.
Due to noise or non-accurate alignment there are tiny parts of geometry floating around in space in the merged voxels. In future processing (outside this thesis) of the model it is essential that these polygons are removed. In GSI Studio it is possible to manually select the geometry to be deleted.
A hole filling approach is integrated in the volumetric surface reconstruction algo-rithm in PlyMC, i.e. the software used for merging, but the result is not sufficient. In the software from GSI, holes are filled using a hole filling tool where the user selects
the boundary of a hole followed by an automatic triangulation of that hole (further described in section 4.5.2).
3.6.3
Decimation
The choice of decimation procedure has to suit the processes and requirements for the larger project the model is used in. In that project an important feature is the ability to preserve edges and texture coordinates of the meshes when decimated. The require-ments include low-resolution and high-resolution voxels connecting without seams and the texture map for a voxel only computed once. Preserving the texture mapping space allows the same texture map to be used for either high or low-resolution geometry. By using the QSlim7software based on the quadric error technique [12] introduced in sec-tion 2.5.2 these requirements are achieved. The model is decimated into three different resolutions as shown in figure 3.13.
Figure 3.13: Decimation of the surrounding environment; 5, 10 and 20 percent respec-tively of the original amount of polygons.
Chapter 4
Implementation and Software
This chapter describes the various implementations further as well as it goes into detail about how the software in the 3D data processing pipeline is used.
4.1
Filtering
The implementation of the Gaussian blur filtering program have to consider memory issues. When performing the filtering algorithm a temporary portable float image with the same size as the input portable float image is stored. Two full size range images in full resolution are not possible to store into the working memory on an ordinary work station (1 GB RAM). Thus, an out-of-core technique had to be implemented. By dividing the filtering process into two temporary parts, filtering first the left half of the image and saving it on the hard drive and then repeating the same with the right half the problem is solved. The edge in between them is treated with one pixels overlap resulting in an continuous filtering. For the resulting image the left and right halves are imported one by one and the values are copied. This implies in worst case 1.5 full size range image being stored into memory.
To filter one single range image is rather time consuming. For efficiency reasons the filtering algorithm is implemented taking multiple (no upper limit) range images as input. The size reduction, i.e. keeping every other point in both the i and j direction in the range images is performed in the same implementation. The result of this step in the pipeline is that all range images are saved as new filtered range images having a fourth of the original size (see figure 4.1).
4.2
Registration
The implementation of making the selection of regions possible provides an user inter-face. It is very important to carefully perform the selection of regions since this is all data the actual aligning algorithm will be based on. In the optimal case there are four pairs of regions where each pair includes the same part of geometry.
The complex task to perform the global registration, i.e. to determine which range images to be pairwise aligned to each other is solved for regarding the pre-aligned scans. Range images without pre-alignment have no initial estimate and it is therefore even more important to select appropriate regions. Thus, manual work demanding time and effort was spend finding the best solutions.
From the pre-processed range images the user selects the regions by dragging the mouse over the image, the first region sets the size limit since a new PFM with a size four times the first selection’s size is created. When four regions (see figure 4.2) are selected the new range image is triangulated and ready to be imported into MeshAlign.
Figure 4.2: Selection of regions (left) and corresponding new PFM:s (right). The manual part in the semi-automatic aligning process is the initial estimate which has two different approaches in MeshAlign. These are either a point based method or a manual method.
The point based method allows the user to point out four pairs of points between the mesh to be registered and its reference mesh. This method is not very efficient for the geometry of the Parthenon since it showed to be very difficult to point out specific feature points. Sculptures with more obvious details are more suitable for this approach [25].
The better but more manual method gives the user the ability to move and rotate the mesh freely until it is placed in the correct position. In both methods, the process implementing the ICP algorithm is performed. By comparison, the last method gives a more accurate alignment but is more time consuming than the point based method.
For some cases a solution with a two-step process had to be invented since
Me-shAlign has some restrictions. If the mesh which is to be registered is too far off from
the reference mesh it is very hard for the user to perform the initial estimate in an accurate way due to the way the 3D graphical user interface in MeshAlign is imple-mented. The total bounding box containing the geometry of the mesh pair is too large and prevents the user to zoom in on details.
Figure 4.3: The manual alignment method in MeshAlign
possible to import the whole corresponding mesh into MeshAlign (see figure 4.3). A first initial estimate is performed and the original range image is transformed by using the obtained ALN file. In the second step, regions are selected from this transformed range image. Thereafter, when the triangulated range image is imported into MeshAlign again an accurate alignment is much easier to perform since zooming is possible, and ICP can be run on high resolution geometry resulting in a second ALN file. The various steps in the registration phase are illustrated in figure 4.4.
4.3
Voxelization
The voxelization constitutes one of the main problems in this thesis. It makes the processing and management of large datasets possible and efficient, why this section receives much focus.
4.3.1
Adjustment of the Coordinate System
The coordinate system in which the scanner stores the data is not of standard type having its z-axis pointing up. Additionally, both the xy-plane and the xz-plane are not perpendicular to the Parthenon structure. This results in the voxels not being aligned to the actual structure. Some voxels can only include a small part of a column when other hold a lot of geometry. As a consequence a voxelization is not as efficient as it theoretically can be. Thus, the voxels needs to be adjusted to the Parthenon in terms of location in 3D space.
In order to change coordinate system the actual orientation of the Parthenon have to be found. First a plan of the desired voxelization is made, and for convenience the grid is chosen to align with the principal architectural features of the structure (see figure 4.5).
Figure 4.5: Voxelization of the facades; east (top left), north (top right), west (bottom
left) and south (bottom right).
Thereafter, a voxel size of 4.2 meters is determined by averaging virtually per-formed measurements. The axis direction vectors are determined, and the cross prod-ucts calculated to prove that the basis is perfectly orthonormal. The rotation matrix is build by inserting the x, y and z vectors. To be able to transform the original coordi-nate system to the new one an anchor point has to be measured. Finally, a function is implemented that rotates and translates a point in the original coordinate system into the new. The result is a coordinate system having its origo at an edge of a voxel, and which is adjusted to have all points on the Parthenon having positive (x,y,z) values.
The adjustment of the coordinate system changes the partitioning of space so that it becomes static, i.e. the grid containing the Parthenon structure is now always divided into 8×5×17 equally sized voxels. Previously, the grid was dynamic and the number of voxels could be changed. Even though, that in order to obtain a reasonable size of the voxels, the total volume (with the environment included) was partitioned into
40×50×20 voxels. Additionally, besides the adjustment of the coordinate system all points belonging to the surrounding terrain are neglected and processed separately.
This change in the voxelization phase made the following integration phase sig-nificantly more efficient. Since the surrounding environment is not voxelized and the number of voxels including geometry is less, the time needed for the merging decreased from 46 hours to 5.5 hours.
4.3.2
Partitioning
The idea behind the voxelization is to traverse through each range image and organize all non-zero (x,y,z) points into corresponding voxels (see figure 4.6). An overlapping region is used in order to get a continuous surface with connecting points over the boundaries. Consequently, a point possibly belongs to one up to eight voxels. The size of the overlapping region is a dynamic parameter, thus the overlapping region can be adjusted to suit different resolutions of a range image. The index for a voxel consists of a value representing how far away in the (x,y,z) directions from origo the voxel is located. In the Parthenon case the indices start at voxel_0_0_0 and continues up to
voxel_7_4_16.
The most time consuming phase is the traversal of the huge range images why each range image only is traversed once. The algorithm follows this procedure:
For each point in the range image:
1. Get the (x,y,z) value if non-zero.
2. Transform the point into the new coordinate system.
3. Calculate the 3D position in the current voxel grid to obtain a voxel index for the point.
4. Check if the point lies within a boundary region and should be treated as an over-lap point. If so, determine additional voxel index (up to seven possible additional cases).
5. Store the original (i,j) position from the 2D range image for the current point and the determined voxel index/indices in the voxel array.
The voxel array storing the data consisting of voxel indices and the (i,j) positions of the points in the voxels tend to be very large (see table 4.1). For efficiency and memory issues a resize function of the array is implemented.
Phase Number of points (maximum)
Original range image 18000×3200 = 57,600,000 Pre-processed range image 57,600,000÷4 = 14,400,000
Table 4.1: Number of points to be voxelized.
The way the voxel grid is determined was implemented in two different ways. The first version have a dynamic voxel grid depending on the input range images. The
maximum and minimum (x,y,z) extents considering all range images is determined and set as the final voxel grid. The user then adjusts the parameters deciding how many times each direction is divided. During the voxelization all range images are traversed twice. This version is preferred when the extents of the geometry is unknown.
The second implementation has a voxel grid adjusted to known measurements of the Parthenon coordinates (described in section 4.3.1). Thus, the time consuming grid calculation is avoided and the range images are only traversed once resulting in a more efficient voxelization. This is the implementation used for the final model.
Figure 4.6: Voxelization phase.
4.3.3
Triangulation
The kernel in the implementation of the triangulation algorithm is the connectivity function. Neighboring points from the two-dimensional range image are sent to the function and tested if they are to connect or not (figure 4.7). The connectivity is based on two thresholds based on the euclidian distance between points in 3D space.
p1
p4 p3
p2
Figure 4.7: Triangulation scheme.
Firstly, the distance threshold is simply the euclidian distance between the two points. Secondly, the scaled distance threshold is the quota between the distance
threshold and the average distance from the two points to the location of the
scan-ner. The location of the scanner for each scan is obtained from a separate file which includes all scanner locations for each session. Consequently, points that are further away from the scanner are connected even though the distance between the two points is greater than the distance threshold. This is necessary due to a sparser distribution of points the further away from the scanner location. By using suitable thresholds, accurate geometry without unnecessary holes and incorrect connections is obtained. Equations 4.1 and 4.2 show the thresholds deciding whether the points p1 and p2 will connect or not.
distance_threshold =
q
(p1x− p2x)2+ (p1
scaled_distance_threshold = distance_threshold
scanner_distance (4.2)
The thresholds (see table 4.2) in the connectivity function are easily adjusted to suit different resolutions of the range images. This is convenient since subsampling showed to be necessary for certain cases in the alignment phase (mentioned in section 4.2).
Range image Distance threshold Scaled distance threshold
Original (100%) 3.0 cm 0.0002 Pre-processed (25%) 6.0 cm 0.0004
Table 4.2: Connectivity thresholds.
The triangulation problem, i.e. to triangulate based on the image plane but from points stored in a voxel array (introduced in section 3.4.2 and illustrated in figure 4.8) is solved and implemented in following algorithm. The algorithm is possible since the points in the array appear in scanline order.
(difference and j-difference are the differences between adjacent points’ different i-and j-values)
For all voxels in each range image:
1. Get the extent in the i-direction. This indicates maximum number of searches allowed in order to find p3 and p4.
2. Loop through the voxel array, when between two adjacent points the i-difference=1 and j-difference=0 the points p1 and p2 are found. Then perform a search for p3 (i-difference=0 and j-difference=1 between p1 and p3) and finally a search for
p4 (i-difference=1 and j-difference=0 between p3 and p4). The connectivity for
the found points are tested by using the connectivity function. 3. Let p2 be p1 and go to step 2.
4. Repeat until the end of the voxel array. 5. Goto step 1 and compute the next voxel.
p1 p2 p3 p4
extent in i-direction
Figure 4.8: The voxel array.
After all voxels have been processed by the algorithm the vertices and faces (i.e. the indices to the vertices) are saved in the PLY file format and each voxel mesh gets a specific filename revealing which original range image it cames from and which voxel it is. The writing of the PLY file is not described in the algorithm above, only the solution to the triangulation problem.
