Algorithms and tools for the automatic scanning of mid-complexity 3D objects

Algorithms and tools for the automatic scanning of mid-complexity 3D objects Master Thesis Author: Alex Carruesco Llorens Director: Carlos And´ ujar Gran Department of Computer Science (CS) Advisor: Pere Brunet Crosa Department of Computer Science (CS) Master in Innovation and Research in Informatics Computer Graphics and Virtual Reality specialization July 4th, 2016 Abstract (English) In recent years, interest in 3D applications has substantially grown and this tendency seems to increase in the future. There is also an increasing interest in 3D printing, thus it would be interesting to have a way to easily create 3D models from real objects. In this project we designed and implemented a complete scanning system which acquires 3D models of mid-complexity objects in a few minutes. The scanning system consists of a turntable and a 3D scanner, and our automatized system requires minimum user interaction. Moreover, the user is free to configure several parameters of the process. We also implemented a simulator for our 3D scanning system in order to discover which kind of situations may be problematic for the scanning process. This simulator allows us to detect some problems such as the under-scanned areas of an object, or configurations of the scanner location which hamper the scanning of the object. Finally, we tested both the 3D scanning system and the 3D scanner simulator with different configurations in order to evaluate their behaviour under different circumstances and we discussed the results. Abstract (Catalan) Als ´ultims anys, l’inter`es en aplicacions 3D ha crescut substancialment i aquesta tend`encia sembla que continuar`a augmentant en el futur. Tamb´e ha hagut un inter`es important en la impressi´o 3D, aix´ı que seria interessant disposar d’una forma senzilla de crear models 3D a partir d’objectes reals. En aquest projecte hem dissenyat i implementat un sistema d’escaneig complet que genera models 3D d’objectes de complexitat mitja en pocs minuts. El sistema d’escaneig est`a format per una taula girat`oria i un esc`aner 3D, i requereix molt poca interacci´o per part de l’usuari. A m´es, l’usuari ´es lliure de configurar diversos par`ametres del proc´es. Tamb´e hem implementat un simulador per al nostre sistema d’escaneig 3D per a descobrir quin tipus de situacions poden causar problemes al proc´es d’escaneig. El simulador ens permet detectar problemes com zones de l’objecte sense mostres o amb poca densitat de punts, o configuracions de l’esc`aner que dificulten l’escaneig. Finalment, hem provat el sistema d’escaneig 3D i el simulador d’esc`aner 3D amb diferents configuracions amb l’objectiu d’avaluar el seu comportament sota diferents circumst`ancies i hem comentat els resultats obtinguts. Abstract (Spanish) En los ´ultimos a˜nos, el inter´es en aplicaciones 3D ha crecido substancialmente y esta tendencia parece que seguir´a incrementando en el futuro. Tambi´en ha habido un inter´es importante en la impresi´on 3D, as´ı que ser´ıa interesante disponer de una forma sencilla de crear modelos 3D a partir de objetos reales. En este proyecto hemos dise˜nado e implementado un sistema de escaneado completo que genera modelos 3D de objetos de complejidad media en pocos minutos. El sistema est´a formado por una mesa giratoria y un esc´aner 3D, y requiere muy poca interacci´on por parte del usuario. Adem´as, el usuario es libre de configurar diversos par´ametros del proceso. Tambi´en hemos implementado un simulador para nuestro sistema de escaneo 3D con el objetivo de descubrir qu´e tipo de situaciones pueden ser problem´aticas para el proceso de escaneado. El simulador nos permite detectar problemas como zonas del objeto sin muestras o con poca densidad de puntos, o configuraciones del esc´aner que dificultan el escaneo. Finalmente, hemos probado el sistema de escaneo 3D y el simulador de esc´aner 3D con diferentes configuraciones con el objetivo de valorar su comportamiento bajo diferentes circunstancias y hemos comentado los resultados obtenidos. Contents 1 Introduction 9 1.1 Motivation ................................. 9 1.2 Goals .................................... 10 2 Previous work 11 2.1 Point cloud registration .......................... 11 2.1.1 The ICP algorithm ........................ 12 2.1.2 ICP based techniques ....................... 14 2.1.3 Solving the initial alignment ................... 16 2.1.4 Other registration techniques .................. 20 2.2 Surface reconstruction from 3D point clouds .............. 22 3 Scanning system 25 3.1 Scanning system setup .......................... 25 3.2 System pipeline .............................. 26 3.3 Data acquisition .............................. 27 3.3.1 Calibration ............................ 27 3.3.2 Real-time filtering ......................... 30 3.3.3 Kinect data issues ......................... 32 3.4 Point cloud processing .......................... 35 3.4.1 Merge of individual scanner frames ............... 35 3.4.2 Outliers removal ......................... 36 3.4.3 Selection of candidate points ................... 39 3.4.4 Mesh reconstruction ....................... 41 3.5 System pipeline revisited ......................... 42 1 3.6 Discarded methods ............................ 43 3.6.1 ICP refinement for the scan merge step ............. 43 3.6.2 ICP refinement for the global point cloud ............ 44 4 Scanner simulator 45 4.1 The simulator setup ........................... 45 4.1.1 Line-based scanner ........................ 46 4.1.2 Matrix-based scanner ....................... 47 4.2 Working pipeline ............................. 48 5 Software implementation 49 5.1 Dependencies ............................... 49 5.2 Implementation details .......................... 50 5.3 GUI .................................... 51 5.3.1 3D Scanner ............................ 51 5.3.2 Model Explorer .......................... 53 5.3.3 3D Scanner Simulator ...................... 55 6 Results 56 6.1 3D scanning system ............................ 56 6.1.1 Voxelization ............................ 57 6.1.2 Mesh reconstruction ....................... 59 6.1.3 More than a single object at the same time ........... 64 6.1.4 Our method compared to Kinect Fusion ............ 65 6.1.5 Final results ............................ 66 6.2 3D scanner simulator ........................... 70 6.2.1 Line-based ............................. 70 6.2.2 Matrix-based ........................... 72 7 Conclusions & Future work 74 7.1 Conclusions ................................ 74 7.2 Future work ................................ 75 A Kinect 2 specifications 77 2 List of Figures 2.1 a) The initial point clouds. b) Registration using ICPBS. c) Registration using sPICP. d) and e) are local amplifications of ICPBS and sPICP respectively. ............................ 16 2.2 Top: 10 input scans. Bottom left: registration after applying Robust Global Registration [GMGP05]. Bottom right: after ICP refinement. 17 2.3 A: Input point clouds. B: The erroneous alignment by direct correlation of histograms. C: Correct alignment using constellation images. D & E: The corresponding histograms for A. F: The histogram of the correct alignment. ................................ 18 2.4 Registration by 4-PCS dealing with outliers. The percentages are the proportion of outliers respect the original amount of data points. . . . 20 2.5 Intuitive illustration of Poisson reconstruction in 2D. ......... 23 2.6 Work of [GCSA13]. Left: raw point cloud with additional noise on the top half part. Middle: point cloud and reconstruction. Right: reconstruction only. ............................ 24 3.1 Scanning system running. The upper-right red ellipsoid shows the Kinect 2 sensor placed in a tripod and the bottom-left ellipsoid shows the turntable. ............................... 26 3.2 The scanning system pipeline. Top: The data acquisition process. Bottom: The point cloud processing. .................. 26 3.3 A screenshot of the manual point selection step. ............ 28 3.4 The circumcircle of three points A,Band C, and the corresponding circumcentre O............................... 29 3 3.5 Left: A complete raw data frame from Kinect 2. Right: Raw data with the real-time filtering. .......................... 31 3.6 Left: The real wooden A object. Right: A raw point cloud from the wooden A, captured with Kinect 2. ................... 32 3.7 Screenshot of wooden A raw data, showing the invalid points captured by Kinect 2 next to the object’s contour. ................ 34 3.8 Lateral view of the wooden A raw data. ................. 34 3.9 Left: Raw merged point cloud of the wooden letter A. Right: Merged point cloud of the wooden A after applying the statistical filter at each scan. .................................... 36 3.10 A voxelization example of size 1503and threshold 100 applied to the wooden letter A. ............................. 38 3.11 Left: a voxelization without filling the base of the object. Right: the same voxelization, but with a filled base. ................ 38 3.12 Example of a 2D case of the Sub-voxel Counting technique. The colours are used to mark the contribution of each pair of point-scanner to each traversed voxel. Thus the selected candidate is the brighter one, which has 3 contributions. ............................ 40 3.13 The scanning system pipeline revisited, with the user interaction coloured in green. .................................. 42 3.14 Merged point cloud of the wooden A, applying ICP between scans. . 43 4.1 A line-based scanner, from the LMI3D Gocator Series. ........ 46 4.2 The Armadillo model in the scanner simulator, using a line-based scanner. 47 4.3 The Armadillo model in the scanner simulator, using a matrix-based scanner. .................................. 47 5.1 Screenshot of the 3D Scanner mode GUI. ................ 52 5.2 Screenshot of the Model Explorer mode GUI. ............. 54 5.3 Screenshot of the 3D Scanner Simulator mode GUI. .......... 55 6.1 Test objects: Wooden A, Doll, Cylinder, Action Figure, Virgin. . . . 56 4 6.2 The merged scans for the non-centred Doll. ............... 57 6.3 Doll with different voxelizations. From left to right (size, threshold): (50, 600), (100, 300), (150, 120), (200, 30). .............. 58 6.4 Doll with different voxelizations. From left to right and up to down (dimensions, threshold): (128, 20), (128, 100), (128, 200), (128, 300). 59 6.5 Top: the centroids method. Bottom: the sub-voxel counting method. Left: the point cloud of selected candidates. Right: the reconstructed mesh. .................................... 60 6.6 Poisson reconstruction with different values of tree depth. From left to right: 4, 6, 8 and 10. ........................... 61 6.7 Poisson reconstruction with different values of samples per node. From left to right: 1, 2, 4, 8 and 16. ...................... 61 6.8 Poisson reconstruction with different values of scale. From left to right: 1.1, 1.5, 2.0, 4.0 and 8.0. ......................... 62 6.9 Zoom to the area where the Poisson reconstruction failed. ....... 63 6.10 Results varying the radius neighbourhood of the normals’ estimation method. From left to right: 0.1, 0.05, 0.02, 0.01, 0.005. ........ 63 6.11 Mesh reconstruction with the concave hull method. .......... 64 6.12 Scanning system results with 2 objects at the same time. ....... 64 6.13 Left: Results of our scanning system. Right: Results of Kinect Fusion. 65 6.14 Results obtained with the Poisson reconstruction (part 1). From left to right: voxelization, centroids candidates, centroids mesh, sub-voxel counting candidates, sub-voxel counting mesh. ............. 67 6.15 Results obtained with the Poisson reconstruction (part 2). From left to right: voxelization, centroids candidates, centroids mesh, sub-voxel counting candidates, sub-voxel counting mesh. ............. 68 6.16 Two different views of the first experiment for the line-based scanner. 71 6.17 The samples obtained by the scanner with the proposed extra degree of freedom for line-based scanners. .................... 71 6.18 Results of the simulation rotating the scanner about its Z axis. . . . . 72 6.19 Two different views of the first experiment for the matrix-based scanner. 73 5 where N=|P|,pi∈P,qi∈Q,c(i)denotes the pair correspondence of pi, and T(·)is an affine 3D transformation. We can distinguish different tasks in this problem: the selection of relevant points (ideally, we want to consider only the points which overlap and not the complete clouds); the correspondence between points of both point clouds (i.e. calculate the above mentioned indices c(i)); and the computation of the affine transform which minimizes the error metric. The state of the art shows that the majority of the existing methods are based on the ICP algorithm, either being a variant of this technique or using it to refine the obtained solution. Therefore, the following section will explain the basis of the ICP and its limitations; then we will visit some ICP interesting variants; afterwards we will see some solutions to the ICP initialization issue; and finally we will explore some other methods for point cloud registration. 2.1.1 The ICP algorithm The Iterative Closest Point (ICP) algorithm was formally introduced by Besl and McKay [BM92], but there was a previous work by Chen and Medioni which was closely related [CM91]. Since ICP showed quite good results, it has been the groundwork for a number of publications. The ICP algorithm for finding the transformation which best aligns Qto Pis detailed in Algorithm 1. Chen and Medioni [CM91] proposed a very similar algorithm, but instead of minimizing the distance between points, their algorithm minimized the distance of the points belonging to Qtowards the tangent plane at the points of the reference point cloud P(this method is also known as point-to-plane ICP). 12 Algorithm 1 ICP base algorithm for the 3D case 1: Data: Reference point cloud Pand point cloud Q. 2: Result: A 3D rigid transformation T. 3: procedure ICP(P,Q) 4: v= [1,0,0,0,0,0,0]T.The rotation quaternion and the 3D translation 5: τ= variation threshold 6: e0= inf 7: k=0 8: repeat 9: k=k+1 10: Qk= apply transform(Q,v) 11: Pairs = compute nearest point in Pfor each point from Qk 12: (v,ek) = minimize euclidean distance by least squares(Pairs) 13: until ek−1−ek>τand k≤maxIterations 14: return v The convergence if the ICP algorithm has been proved, and it is based on two key ideas: the least squares registration generically reduces the average distance between corresponding points during each iteration and the closest point selection reduces the distance for each point individually. The complete demonstration can be found in [BM92]. The main problem of this methodology is that it may easily fall in local minima, thus it will align correctly two point clouds only when they are relatively closer. This issue can be solved by estimating an initial rough transformation, and then applying ICP using the computed transformation as initial state. This technique achieves good results in a moderate amount of time (it depends on the convergence threshold and the amount of points, but it generally works fast if a good nearest neighbours algorithm is applied, such a kd-tree). But the algorithm does not take into account the noise and the outliers, and is also very sensitive to the real overlap between point clouds. All these issues opened the door for further improvements. 13 2.1.2 ICP based techniques In the literature, we can find a number of research publications based on the ICP algorithm. We present here some interesting ICP based techniques which improved the performance of the ICP algorithm or enhanced certain properties of the algorithm: 1) The PLICP (point-to-line ICP) [Cen08]: this ICP variant was intended for 2D problems and it aligns two point clouds by computing a polyline in the reference cloud. Then, it tries to align the second cloud to that polyline, thus it substitutes the point-to-point correspondence of the original ICP by a point-to-line correspondence and the linear minimization problem by a non-linear one. The method exhibits a quadratic convergence property (which implies a faster convergence) and the authors introduced a closed-form solution for solving the non-linear minimization problem in the 2D case, which allows a faster performance. However, its main drawback is the need of a good first guess for the ICP initialization to achieve accurate results. 2) Generalized ICP [SHT09]: this technique was intended to unify the ICP and the point-to-plane ICP into a single probabilistic framework. It replaces the common least squares minimization step by a new optimization function which takes into account the covariance matrices of the point clouds, which are assumed to be generated by a Normal distribution. Since this method introduces dispersion measures (i.e. covariance matrices), it allows improving the results of the basic ICP in presence of noise. 3) Colour point cloud registration with 4D ICP algorithm [MGP11]: Men et al. presented a new technique which uses the colour information of a point cloud to enhance the registration of the ICP. To integrate this information, the colour is codified using the hue calculated from the RGB data, and it is set as the 4th dimension of the point cloud. Then, the hue is weighted and added in the closest point search when computing the point pairs, achieving more precise correspondences. To allow an easier weighting of the 3D position and the hue, all the values are normalized: the 3D position 14 is normalized using the maximum range of the scanner and the hue is normalized between 0 and 1. The hue is only used in the point correspondence step, thus the minimization is not modified with respect to the original ICP. The results exhibit a faster convergence with a similar accuracy, thus this method improves a bit the ICP performance. 4) sPICP (Probability ICP with bounded scale) [DLB+16]: the goal of this method is to improve the performance of the ICP in the presence of noise, also considering an isotropic scale factor. This technique is based on a previous work called ICPBS (ICP with bounded scale, [DZYW07]), which only takes into account the scale factor in the optimization step. This new method introduces a Gaussian distribution in order to model the noise of the point clouds. The variance of the Gaussian is initially big, resembling a uniform distribution, but it is updated at each iteration using the estimated variance (which is computed using the data and the calculated 3D transformation at that iteration). Thus, the algorithm reduces the chance of falling into local minima by adjusting the variance step by step. The main difference with respect to the original ICP is the conversion of the optimization problem into a likelihood maximization using the Gaussian probability. The proposed algorithm is similar to the EM-ICP [GP02]: •The E-step establishes the one-to-one correspondence between Xand Y, using the previous (Rk−1,tk−1,sk−1), as ICPBS does (i.e. finding the nearest point). This step is accelerated by nearest neighbours search methods, such as Delaunay triangulation or k-d trees. •The M-step computes the rotation matrix, the translation vector, the scale factor and updates the variance. This technique showed a good accuracy (see Figure 2.1) and a fast performance, but it would be interesting to test it with other datasets with lower overlap between point clouds. Furthermore, there are some other interesting ICP based approaches: the LM-ICP [Fit03], which proposed a non-linear method for solving the optimization problem; the 15 Figure 2.1: a) The initial point clouds. b) Registration using ICPBS. c) Registration using sPICP. d) and e) are local amplifications of ICPBS and sPICP respectively. work of Jost and Hugli [JH03], which developed an heuristic approach for enhancing the ICP performance; Minguez et al. [MLM05] introduced a new distance metric for the minimization problem; or the Sparse ICP [BTP13], proposed by Bouaziz et al. in order to cope with the outliers and missing data of the point clouds. 2.1.3 Solving the initial alignment As mentioned before, one of the main problems of the ICP algorithm is the importance of the initialization. Gelfand et al. [GMGP05] proposed a method for estimating a coarse alignment, which is refined applying ICP-like methods later on. The main strengths of this technique are that it does not assume anything about the initial position of the point clouds and it is quite resistant to noise. The algorithm of Gelfand et al. [GMGP05] is based on finding correspondences for a small amount of feature points from the reference point cloud P with points of another point cloud Q, in order to acquire a coarse alignment between both point clouds, which is refined with some ICP iterations later on. The main property of the feature selection step is that feature points should come from regions with rare descriptor values, thus the potential correspondences are minimized. The algorithm works as follows: 16 1. Computation of histograms: compute a histogram of descriptor values for all points of P; and for all points of Q. 2. Selection of feature points: Select feature points in both point clouds, identifying the least populated bins of the histogram. To avoid selecting points belonging to the same feature, when a point is selected, all the neighbouring points that fall into a ball of a certain radius are marked as unavailable for selection. 3. Matching of points: match feature points from Q with their best correspondence in P according to a certain metric by applying a branch and bound method. 4. Computation of the transformation: calculate the 3D affine transformation which minimizes the root-mean-square error (RMSE) of the feature point pairs. 5. Refinement: refine the solution with some ICP iterations, starting from the previously obtained transformation. Figure 2.2: Top: 10 input scans. Bottom left: registration after applying Robust Global Registration [GMGP05]. Bottom right: after ICP refinement. 17 As feature descriptor, the authors recommended to use the local volume because it is an integral descriptor from the point neighbourhood that is easy to compute and is quite insensitive to noise. This algorithm showed promising results when there are strong features in both point clouds. But it becomes very computationally expensive when there are none of those remarkable features. Figure 2.2 demonstrates an example of 3D point cloud registration obtained with this algorithm. Makadia et al. [MPD06] introduced an algorithm which aligns two point clouds, independently from the initial poses. The authors claimed that the algorithm also works even with low percentages of overlap between the point clouds (as Figure 2.3 illustrates). Figure 2.3: A: Input point clouds. B: The erroneous alignment by direct correlation of histograms. C: Correct alignment using constellation images. D & E: The corresponding histograms for A. F: The histogram of the correct alignment. Given two input point clouds, the steps of the algorithm are: 1. Compute surface normals fields for both point clouds. 2. Create an Extended Gaussian Image (EGI, [Hor84], which is a mapping of surface normals of an object onto the unit sphere) approximated by an orientation histogram for each point cloud. 18 3. Estimate the rotation by correlating orientation histograms in Fourier space, since it simplifies the problem. 4. Estimate the translation by correlating the already rotation-aligned point clouds. 5. Accept the transformation only if it passes the verification step, which consists in checking the consistency of surface normals in the overlapping regions. 6. If the transformation does not pass the verification, convert histograms to constellation images (i.e. retain the local maxima of the histogram and suppress the remaining bins) and correlate them, selecting the best transformation. 7. Refine the solution with ICP. The main properties of the EGI are that it is invariant from translation and the rotation of an object causes an equal rotation of the EGI, thus it captures the orientation of the underlying surface. On the other hand, EGI is only unique for convex objects, which is not the general situation. However it was shown in [MSD04] that signal correlation provides a reliable measure for the rotational alignment, thus the authors use this fact in the Fourier space to simplify the problem. Similarly to the two previous papers, Aiger et al. designed the 4-PCS (4 points congruent set) technique [AMCO08], aiming to establish a coarse alignment between point clouds which is refined later on with ICP-based methods. Given two arbitrary point sets Pand Q, the 4-PCS algorithm consists in a RANSAC loop (introduced by Fischler et al. in [FB81]), where each iteration works as follows: 1. Select a coplanar base Bfrom P, formed by 4 coplanar points. 2. Find congruent bases Uto Bin Q, with a certain tolerance. 3. For each base Ui∈U, compute the best rigid transformation Ti(by solving the least squares problem) and assign it a score depending on the amount of aligned points below a certain threshold. 19 4. Update the optimal transformation matrix Topt with the best obtained Tiin this iteration only if it improves the solution (i.e. it aligns a larger amount of points). After this process, the algorithm achieves a rough alignment between both point clouds, which can be refined by applying an ICP method. The authors showed that 4-PCS is independent from the initial poses of the point clouds and it works correctly even when there is a small overlap between the point clouds. It also works well in the presence of outliers, as Figure 2.4 shows. Figure 2.4: Registration by 4-PCS dealing with outliers. The percentages are the proportion of outliers respect the original amount of data points. As 4-PCS had some performance issues (such as quadratic complexity), Super 4-PCS [MAM14] was developed to enhance its performance. Mellado et al. proposed an algorithm which reduced the O(n2)complexity of the 4-PCS to linear complexity with respect to the amount of points, achieving a more applicable approach. 2.1.4 Other registration techniques Aside from the ICP algorithm and its variants, there are other techniques which solve the problem of point cloud registration, achieving good results. Myoronenko et al. designed the Coherent Point Drift algorithm [MS10], which considers the alignment of two point clouds as a probability density estimation problem, 20 where one point cloud represents the GMM (Gaussian Mixture Model) centroids and the other one represents the data points. Then, the likelihood is maximized in order to achieve the best alignment. As an advantage, this technique works for rigid and non-rigid registration. The rigid case is solved by a closed-form solution for the EM algorithm and imposing spatial coherence constraints in the GMM centroids. In the non-rigid case, the coherence constraint is imposed by regularizing the distance field and using variational calculus to derive the optimal transformation. This technique requires quite complex mathematical algorithms, thus the authors proposed fast approximations to achieve a better time performance. Despite these efforts, the algorithm remains very slow compared to other state of the art solutions, but the obtained results exhibit a good accuracy and it works with some noise and outliers. Another interesting technique is the one proposed by Zeming et al. [ZB11], where the curvature feature is used for the task of registration. The algorithm follows these steps: 1. Compute the Gaussian and mean curvatures Kand Hof each point using a quadratic surface fitting method, and use them to compute the principal curvatures k1and k2. 2. Extract feature points according to the maximum local changes of curvature. 3. Obtain the initial match points by computing the Hausdorff distance of curvature. 4. Acquire the accurate match points by the circumference shape feature of the local surface. 5. Estimate the rotation and translation matrices from the quaternion, and the registration accuracy is improved by an iterative method. The proposed iterative method consists in finding a new match point set using the error ϵobtained by the quaternion method and computing the transformation again with that method until the error is small enough. 21 Hence, the calibration process consists in the following steps: 1. Place the sensor in a stable position, looking towards the turntable. 2. Acquire a single point cloud from the sensor, which should be obtained by averaging the information of some consecutive frames in order to reduce the impact of the noise. 3. Manually select Npoints of the turntable’s perimeter with the provided GUI, where Nshould be at least 3 (illustrated in Figure 3.3). Figure 3.3: A screenshot of the manual point selection step. 4. Compute the circumcircle and the normal for all the possible triangles from the selected points (henceforth the set of triangles T). Given 3 points, the circumcentre is the intersection of the perpendicular bisectors of the triangle’s edges, and the radius of the circumcircle is the distance from that centre to each of the points (Figure 3.4). As usual, the normal of the triangle is computed by normalizing the cross-product of 2 edges of the triangle’s edges. 5. Compute the turntable’s centre as the average of the circumcentres: O=P tri∈T circumcentre(tri) |T|(3.1) 28 6. Compute the turntable’s radius as the average of the circumcircles’ radius: r=P tri∈T circumcircle radius(tri) P tri∈T area(tri)(3.2) 7. Compute the turntable’s rotation axis as the weighted average of the triangles’ normals: −→ RA =P tri∈T area(tri)·normal(tri) P tri∈T area(tri)(3.3) Since normals of thin triangles can be very sensitive to noise, each normal is weighted by the area of its corresponding triangle in order to obtain a more robust normal estimation (instead of a simple average). Figure 3.4: The circumcircle of three points A,Band C, and the corresponding circumcentre O. Thus, in the calibration process the user only needs to select some points of the turntable’s perimeter; all the rest is automatically computed. Moreover, an error measure for the turntable centre is provided in order to have an idea of how precise the calibration is. The error measure is the mean of the difference 29 between the distance of each point towards the circumcentre and the radius, computed as follows: ϵO= N P i=0 k(pi−O)k − r N(3.4) where piis the i−th selected point, Ois the turntable centre, and ris the turntable radius. This error explains the quality of the turntable centre estimation using the selected points, thus it is necessary to select proper points at the turntable boundary. For instance, points forming an almost perfect circumference will report a near 0 error, but if the circumference is not centred at the turntable centre, the centres of the circumcircles will not coincide with the turntable centre. Therefore, we need the user to select points of circumferences centred at the turntable centre (the most easy-to- select circumference is the turntable’s boundary). 3.3.2 Real-time filtering From the previously obtained calibration, we know the position Owhere the turntable is placed, its rotation axis −→ RA and its radius r, thus we can use this information to filter points in real-time. More precisely, we know that we are only interested in those points inside the cylinder centred at Owith axis −→ RA and radius rand above the plane defined by Oand the normal −→ RA. Thus we can store only those points and discard the rest of the acquired point cloud. This filtering step could also be done as a post-process, but we decided to make it at acquisition time in order to avoid storing data which we certainly know that is completely outside of what we want to scan. This reduces the required physical storage memory and bandwidth, allowing the system to work at interactive frame rates. Figure 3.5 shows the big benefits of this procedure. 30 Figure 3.5: Left: A complete raw data frame from Kinect 2. Right: Raw data with the real-time filtering. The filtering consists in selecting only those points of the point cloud which fulfil two requirements: the point should be over the turntable plane, and its distance to the rotation axis should be less than the turntable’s radius. The first condition is checked using the implicit equation of a plane: −→ n·(X−X0)=0 (3.5) where nis the plane normal and X0is a supporting point of the plane. We evaluate this equation with −→ n=−→ RA,X0=Oand Xbeing a point of the point cloud. If the left side of the equation is greater than 0, the point is above the plane. The second condition is checked by computing the distance from a point Xto the line formed by Oand RA, which should be less than the turntable’s radius. Using trigonometric identities, we can easily deduce the following formula for computing the distance: dist =qk−−→ OX k2−(k−−→ OX k · cos(φ))2(3.6) where φis the angle between −−→ OX and −→ RA. Thus, cos(φ)can be computed using a dot product and, assuming that −→ RA is a normalized vector, the previous equation becomes: dist = v u tk−−→ OX k2−* ,k−−→ OX k · −−→ OX ·−→ RA k−−→ OX kk−→ RAk + - 2 =qk−−→ OX k2−(−−→ OX ·−→ RA)2(3.7) The pseudo-code of the real-time filter is given in Algorithm 2. 31 Algorithm 2 Real-time filter 1: Input: A raw point cloud P, the centre of the turntable O, the turntable’s radius r and the turntable’s rotation axis RA. 2: Output: The filtered point cloud. 3: procedure RealtimeFilter(P, O, r, RA) 4: filteredPoints ←{} 5: for p∈Pdo 6: pointDir ←p - O 7: dot ←dotprod(RA, pointDir) 8: if dot >0then .If the point is above the cylinder base 9: distToCylinderAxisSq = norm(pointDir)2- dot2 10: if distToCylinderAxisSq <r2then 11: filteredPoints.add(p) 12: return filteredPoints 3.3.3 Kinect data issues The acquired data by the Kinect sensor (see full specs in Appendix A) is not as good as we would desire, thus some extra problems arise when acquiring 3D points. In this section we discuss the problems of the raw data using a simple object, shown in Figure 3.6. The object is a wooden letter A, of dimensions 17.4 x 19.4 x 1.6 cm. Figure 3.6: Left: The real wooden A object. Right: A raw point cloud from the wooden A, captured with Kinect 2. 32 The Kinect 2 works at 30 FPS and it retrieves a depth image of 512x424 pixels, thus we can capture at most 217,088 points per frame. Its vertical field of view is 60◦. Since we scan medium-sized objects at approximately 0.9m of distance, we can only use around 30% of the field of view, and therefore around 9% of the image area. In practice, we obtained between 3,000 and 13,000 points per scan depending on the object dimensions. The first problem we faced was that the Kinect 2 scanning technology does not work properly with reflective surfaces, and our turntable’s surface was made of metallic materials. We solved this problem by covering the turntable with paper, but this fact added, as well, a constraint on the scanning system: it does not work properly when scanning objects with reflective materials. We should not think about this constraint as a limitation of the scanning system, but as a limitation of the 3D scanner. Replacing the scanner by another device with different technology may solve this issue. We found another problem when looking at the first real scans we captured: the sensor has some difficulties to capture object contours. As Figure 3.7 illustrates, the sensor captured a considerable amount of invalid points next to the contour of the object. We can see in the left part of the image a sequence of points going from the wooden A contour to the turntable platform, but the thickness of A is merely 16mm, and therefore only a small part of those points are correct. We treated this issue via two techniques: we applied an statistical filter for discarding outliers (explained in Section 3.4.2) and we later on used a voxelization of the point cloud with a certain threshold, which works as a second filter. We encountered as well another common problem: the variability of the depth information in the point cloud. The Kinect sensor captures a matrix of depth values per frame, which allows us to recover a 3D point cloud. But the same position of the matrix will contain quite different depth values in consecutive frames, producing a noisy point cloud. In Figure 3.8 we can see a lateral view of the wooden A raw data that shows a thick point cloud, while the wooden A is a complete flat surface (thus 33 Figure 3.7: Screenshot of wooden A raw data, showing the invalid points captured by Kinect 2 next to the object’s contour. ideally the point cloud should be a plane of points). Hence this thickness illustrates the variability across time of the point clouds acquired by the Kinect 2 sensor. Figure 3.8: Lateral view of the wooden A raw data. It is important to point out that Kinect 2 was not designed for scanning small objects; it was focused on human recognition and skeletal tracking. Therefore, the aforementioned problems are a trade-off for having a cheap 3D scanning device. 34 3.4 Point cloud processing 3.4.1 Merge of individual scanner frames At this point, we have a set of independent point clouds which were taken while the turntable was spinning, but all of them are relative to the same origin and reference system, located at the scanning sensor. In order to obtain a coherent global point cloud, we need to recover the rotation experienced by the object when the data was captured by the sensor. The turntable spins at a reasonably constant speed, thus we can compute the rotation angle using the time difference between the current scan and the first one: φi=ω·(ti−t0)(3.8) where ωis the turntable’s angular speed, t0is the time when the first scan was captured, and tiis the time of the i-th scan. We assume a very small latency with respect to the turntable’s angular speed, which was found to be ω=10.4478◦/s= 0.1823rad/s. Since the Kinect can work at 30 FPS, we can take around 1,000 scans during a spin. Therefore, the error in the angle calculation should be small. Since Kinect frames are separated by about 0.36◦, missing one Kinect frame inadvertently could result in an offset of up to 0.3mm (for an object fitting on a cylinder with 10cm of diameter). Using the rotation angle of the turntable and the calibration information, we can apply the necessary transformations to recover the rotation of the object at each instant. More formally, we apply the following transformation to all the points of the i-th scanned point set Si: p0 i=R(φi,−→ RA)·T(−O)·pi∀pi∈Si(3.9) where T(·)is a translation and R(·)is a rotation about an axis. As the equation shows, we first translate the point by −Oto convert the origin of the turntable into the origin of the reference system. Then, we apply the rotation associated to the i-th scan. 35 At this point all the scans will conform a coherent point cloud. But depending on the sensor position, this point cloud can have a vertical tilt, and for posterior processing steps we require the base of the object to be aligned with the XZ-plane (being the Y axis the vertical one). Hence, we apply a change of coordinates to align the rotation axis of the turntable with the vertical axis: p00 i= * . . . . . . . . . , −→ Xx −→ RAx −→ Zx0 −→ Xy −→ RAy −→ Zy0 −→ Xz −→ RAz −→ Zz0 0 0 0 1 + / / / / / / / / / - ·p0 i(3.10) where −→ RA is the turntable’s rotation axis, −→ Xis the cross product of −→ RA and the turntable’s origin O, and −→ Zis the cross product of −→ Xand −→ RA. As stated in Section 3.3.3, there are many problems with the acquired data. If we do not perform any kind of filtering, we will obtain a very noisy point cloud with a large amount of outliers, as seen in the left image of Figure 3.9. Thus we proposed to apply a statistical filter to each scan in order to remove some undesired points of the point cloud (as seen in the right image of Figure 3.9). Figure 3.9: Left: Raw merged point cloud of the wooden letter A. Right: Merged point cloud of the wooden A after applying the statistical filter at each scan. 3.4.2 Outliers removal There are two filtering steps in our scanning system: the first one is a statistical filter to remove outliers (used in the same way as in Rusu et al. [RMB+08]), and the second 36 one is a voxelization of the complete point cloud to reduce noise and outliers. Statistical filter The algorithm for the statistical filter iterates through the input point cloud twice. During the first iteration, it computes the average distance of each point towards its nearest kneighbours. From those values, both mean and standard deviation are computed, and they are used to determine a distance threshold: distth =mean +λ·sd (3.11) where distth is the distance threshold, and λis a multiplier used to adjust the threshold. All the points with an average neighbour distance above the threshold distth are discarded in the second iteration. As mentioned before, this statistical filter step is applied once to each of the captured scans in the merging phase (Section 3.4.1). Depending on the parameters of the filter (i.e. kand λ), we can remove the majority of the outliers. To avoid discarding feature points, which can be removed if we select strict parameters, we experimentally determined that a neighbourhood of k=10 points and a multiplier λ=1.0 are good enough to discard a number of invalid points, while maintaining as much relevant information as possible. Nevertheless, our framework allows the user to manually adjust these parameters (with the mentioned default values). Voxelization filter After the scans merging phase, we obtain a single huge point cloud, but it still contains a considerable amount of noise and outliers. Thus we need to further process this point cloud. The proposed voxelization filter consists in creating a 3D regular grid inside the bounding box of the point cloud, and count how many points each of the voxels has. After that, the points inside voxels with less than Mpoints are discarded. Voxelization therefore has two configurable parameters: the dimensions of the 3D grid and the threshold M. Figure 3.10 shows an example of such voxelization. 37 3.6.2 ICP refinement for the global point cloud Since the ICP between consecutive scans failed, we thought that it would be interesting to apply it in some way to the global point cloud. We expected that using ICP with all the scans would provide more reliable local information when registering each scan, thus the optimization could work in the proper direction. Hence we iterated through the scans using each scan as the point cloud which must be fitted and the rest of the scans as the target point cloud. Unfortunately, we also discarded this step because it was very time consuming, and there were almost no differences before and after applying this technique. 44 Chapter 4 Scanner simulator One of the goals of this project was to design and implement a 3D scanner simulator to reproduce the data acquisition phase of our proposed scanning system (or other envisioned systems). It would be a very useful tool for analysing the expected behaviour of the system before the real implementation. It also would allow to discover under-scanned regions of the target object and characterize which kind of shapes can be scanned with no significant under-sampled parts. The simulator fundamentals are explained in this chapter. 4.1 The simulator setup The scanner simulator consists in a virtual scanner which rotates about the Y axis of the target mesh. This simulates the behaviour of our proposed scanning system, but instead of making the turntable spin, we rotate the scanner (generating exactly the same result) in order to achieve a faster execution. In order to obtain a realistic and configurable simulator, we take into account some hardware specifications of real scanners, which are: •Near: the minimum distance at which the sensor is able to capture points. •Far: the maximum distance at which the sensor is able to capture points. •Fov-Y : the vertical field of view, which can be used together with Near and Far 45 to determine the scanning volume of the scanner. •Distance between emitter and receiver: the Euclidean distance between the scanner component which emits (emitter) the light signal and the scanner component which receives (receiver) the emitted signal. In addition to these parameters, we can classify the 3D scanners by its acquisition pattern in two types: line-based and matrix-based scanners. Both types of scanner are detailed in the sections below. 4.1.1 Line-based scanner This type of scanners is formed by those scanners which capture a line of points with depth information (see a real example in Figure 4.1). Thus, in addition to the previous parameters, we also need the resolution of the scanner in the line direction, which is used to determine how many point samples per scan will acquire the scanner at its nearest distance. Figure 4.1: A line-based scanner, from the LMI3D Gocator Series. Figure 4.2 shows the configuration of the simulator using a line-based scanner. The green lines represent the view frustum of the scanner. 46 Figure 4.2: The Armadillo model in the scanner simulator, using a line-based scanner. 4.1.2 Matrix-based scanner This type of scanners comprises the scanners which capture a 2D matrix of points with depth information. Hence, we also need the aspect ratio (i.e. width/heiдht) to determine the scanning volume, and also the amount of samples in each dimension of the 2D matrix. The amount of samples of the matrix is considered at the nearest distance. Figure 4.3 illustrates the simulator using a matrix-based scanner, where the green lines delimit the volume of the view frustum of the 3D scanner. Figure 4.3: The Armadillo model in the scanner simulator, using a matrix-based scanner. 47 4.2 Working pipeline The usage of the scanner simulator is quite simple: the user only needs to define the scanner parameters and its location and orientation. Then, the virtual scanning process starts and the scanner rotates about the centre of the turntable, which is the centre of the white circle (by default, the turntable is located at the base of the model’s bounding box, but this can be changed by the user). Therefore, the simulator produces a point cloud around the desired mesh, simulating the data acquisition of the real scanning system. That point cloud is visualized over the mesh in order to ease the inspection of under-scanned areas. In Chapter 6we evaluate the simulator behaviour with different configurations. 48 Chapter 5 Software implementation The software was developed in C++, using Microsoft Visual Studio 2013, and it has some dependencies, explained in Section 5.1. We explain some implementation details in Section 5.2, and the usage of the software in Section 5.3. 5.1 Dependencies These are the external library dependencies of the software: •GLM v0.9.7.2: it is a Mathematics library for graphics software based on the GLSL specifications. GLM is licensed under The Happy Bunny License and MIT License. •Embree v2.9.0: it is a collection of high-performance ray tracing kernels, developed at Intel. Embree is released as Open Source under the Apache 2.0 license. •Kinect SDK v2.0: it is the official SDK for controlling Kinect devices, developed by Microsoft. The software is licensed, not sold, which means it only gives you some rights to use the software. Microsoft reserves all other rights. •Qt v5.5: it is a very extended graphical user interface library. Qt is available under the GNU Lesser General Public License version 3. 49 •Point Cloud Library (PCL) v1.8.0: it is a standalone, large scale, open project for 2D/3D image and point cloud processing. PCL is released under the terms of the BSD license, and thus free for commercial and research use. •Eigen v3.2.8: it is a library for linear algebra. Eigen is licensed under the MPL2 license. •FLANN v1.8.4: it is a library for performing fast approximate nearest neighbour searches in high dimensional spaces. FLANN is distributed under the terms of the BSD License. •Boost v1.59.0: it is a very extended C++ library with a big amount of different functionalities. Boost is released under its own license, which adds no restrictions. •Visualization Toolkit (VTK) v7.0.0: it is an open-source, freely available software system for 3D computer graphics, image processing, and visualization. VTK is licensed under the BSD license. •Qhull v2015.2 : it is a collection of geometric algorithms, such as convex hull or Delaunay triangulation. Qhull has its own license, which states that the software can be freely copied, modified and distributed under certain conditions. 5.2 Implementation details We used some methods that were already implemented to ease the development of the scanning system: •The statistical filter used is the pcl::StatisticalOutlierRemoval, from the PCL library. •For the ICP we used the class pcl::IterativeClosestPoint from PCL, but we also implemented an original version of ICP. •The Poisson reconstruction was the pcl::Poisson implementation, and the normals’ estimation was the pcl::NormalEstimationOMP. 50 •The concave hull was implemented by the class pcl::ConcaveHull from PCL. •The ray tracing was executed using the Embree library. The rest of the scanning system and the scanner simulator was implemented from scratch. 5.3 GUI The graphical user interface of the software provides 3 independent modes of execution, according to the different purposes that the user may have. Those modes are: •3D Scanner: this mode is for calibrating the scanner with respect to the turntable and capturing real data. •Model Explorer: this mode allows exploring and processing the captured data. •3D Scanner Simulator: this mode is to be used for simulating the 3D scanner data acquisition of our scanning system. All three modes have some common UI components at the top of the right panel (see Figures 5.1,5.2 and 5.3): the first option restarts the original location of the camera; the second one defines the camera speed; the third one determines the size of the points in the different visualization modes; and the fourth one shows the bounding box of the rendered models. 5.3.1 3D Scanner The 3D scanner mode provides the calibration of the system and the data acquisition mode. Figure 5.1 shows a screenshot of the whole window of the software in the scanner mode, with the corresponding widgets for this mode on the right side. The proper usage of this mode should follow these steps: 51 1. Place the sensor looking towards the turntable, which should be empty. 2. Get an ”Averaged scan” to obtain a more reliable screenshot of the environment. 3. Press on ”Calibrate”, select the captured averaged scan, and follow the calibration process (i.e. manually select some points on the turntable’s boundary). 4. Press on ”Play” and select the ”Filter points” checkbox to see the captured data from the scanner in real-time already filtered. 5. Select the output folder where we want to store the scans. 6. Place the target object on the turntable, make it spin and press ”Rec” to start recording. Figure 5.1: Screenshot of the 3D Scanner mode GUI. Ideally, the rotation of the turntable should be also controlled from the software. However, as explained in Section 3.1, our turntable was outdated and had no usable SDK, and therefore we needed to control it manually. If we leave both turntable and sensor in the same poses, we can repeat steps 5 and 6 as many times as objects we want to scan, reusing the same calibration. The software 52 saves all user defined parameters, thus the calibration will not change until the user explicitly makes a new one. The calibration of the system can also be saved and loaded from a text file. It is interesting to state that we can also calibrate the system with an object on the turntable since we only need to select some points on the boundary. In any case, it is recommended to do so without any obstacle in order to capture the whole turntable, having more points to select. 5.3.2 Model Explorer This mode allows the user to load different models and apply them the point cloud processing steps in order to obtain the final mesh of the real object. Figure 5.2 shows the GUI for this mode. The user can load a single model or create a merge of scans by selecting the corresponding option in the ”File” menu. Then, the loaded model appears in the list view and the user can freely navigate around it. Below the list of loaded models, there is a section with the configurable parameters for each of the processing steps: •Statistical filter: the K-neighbours defines the amount of neighbours to visit per point when applying the statistical filter, and the SD-multiplier is the λ in the distance threshold (see Section 3.4.2). The user can decide to apply an extra statistical filter to the complete point cloud, but there is no need since it is already applied to each individual scan in the merge phase. •Voxelization: Grid dimensions define the dimensions of the 3D regular grid and Voxel threhsold determines the threshold for discarding voxels during the process. The user can also select the method for the candidates’ selection step using a combo box and, if desired, can save the obtained candidate points after the voxelization. •Mesh reconstruction: there is a combo box to select which mesh reconstruction method should be applied. If Poisson is selected, then the user can configure its 53 sponding mesh reconstruction with the same parameters. We can see that the point cloud formed by the centroids looks very regular, while the point cloud obtained with the sub-voxel counting algorithm seems more natural. Regarding to the obtained meshes, the one produced with the clustering is slightly smoother than the other, while the one generated from the sub-voxel counting seems to preserve slightly better the features of the object. Thus, we let the user choose which method fits better his goals. Figure 6.5: Top: the centroids method. Bottom: the sub-voxel counting method. Left: the point cloud of selected candidates. Right: the reconstructed mesh. Since centroids provide a more regular grid, we used them to evaluate the Poisson parameters in order to minimize the impact of the method in the evaluation. We believed that the sub-voxel counting technique might affect the assessment more than the centroids due to its irregularities. We discuss the following parameters: depth of the tree, samples per node, and scale. The first experiment with the Poisson parameters was varying the depth of the tree (i.e. octree) used in the isosurface extraction, while leaving the samples per node and the scale with values 5 and 1.1 respectively. Figure 6.6 illustrates the impact of the 60 parameter in the mesh reconstruction. For low values, we obtain a very coarse mesh which approximates the target object shape, but loses all the details. Increasing the tree depth, we achieve more detailed meshes, but going beyond certain values (such as tree depth 10) does not provide any substantial gain, while it exponentially increases the execution time. Figure 6.6: Poisson reconstruction with different values of tree depth. From left to right: 4, 6, 8 and 10. Figure 6.7: Poisson reconstruction with different values of samples per node. From left to right: 1, 2, 4, 8 and 16. The second experiment for the Poisson parameters consisted in changing the samples per node used in the isosurface extraction, setting the tree depth to 8 and the scale to 1.1. Figure 6.7 shows the importance of this parameter over the resulting mesh. As we can see, a single sample per node leads to a very irregular mesh. With 2 samples it improves, but it is still quite noisy, and with 4 and 8 samples it looks considerably better. If we set too many samples, the mesh becomes over-smoothed and loses the 61 surface details. Thus values between 4 and 8 seem more suitable for the kind of point clouds we deal with. The third experiment on the Poisson parameters consisted in varying the scale, fixing the tree depth to 8 and the samples per node to 5. Figure 6.8 illustrates the influence of the parameter on the reconstructed mesh. The default value of the PCL implementation is 1.1, thus we tested values from 1 to 8, and we observed that lower values are the ones that work better. With values between 1 and 2 we obtain good results, but a scale factor of 4 starts adding noise and with value 8 the mesh becomes too smooth. Figure 6.8: Poisson reconstruction with different values of scale. From left to right: 1.1, 1.5, 2.0, 4.0 and 8.0. Poisson surface reconstruction works generally well, but in all the previous examples we can see a blob in the feet of the doll (more detailed in Figure 6.9). This case also happens to other point clouds, and it seems to be caused by two facts: it occurs in under-scanned areas (e.g. the bottom part of the doll’s dress) which leads the point cloud to have holes; and it happens in regions with poor normals. Both facts seem to be closely related because the normals are poorly estimated in regions with holes. Therefore, we also studied the influence of the radius neighbourhood in the normals’ estimation step, which is executed just before the Poisson reconstruction. Figure 6.10 shows the results when varying the radius of the neighbourhood. It seems that high values make those blobs disappear, but the resulting mesh is too smooth. For lower values, the results seem better, thus this parameter needs to be carefully adjusted. If 62 Figure 6.9: Zoom to the area where the Poisson reconstruction failed. the value is too low, there are no neighbours, thus normals cannot be estimated. Figure 6.10: Results varying the radius neighbourhood of the normals’ estimation method. From left to right: 0.1, 0.05, 0.02, 0.01, 0.005. The other mesh reconstruction method that we considered was the concave hull, implemented in the PCL library. We can see the results of the method for α=0.005 in Figure 6.11. At a first glance, it may seem that the triangles are correctly computed, but some may have flipped normals. Looking carefully at the geometry we discover some edges with more than 2 incident faces and other geometric irregularities, which made very difficult the correction of the normals, as the mesh does not separate the interior from the exterior. 63 Figure 6.11: Mesh reconstruction with the concave hull method. 6.1.3 More than a single object at the same time Our scanning system is also able to process more than a single object per scan. After the voxelization, we search for different connected components by placing seeds in solid voxels and propagating them (considering the 26-neighbourhood of a voxel in 3D), and we treat each component separately. Figure 6.12: Scanning system results with 2 objects at the same time. Figure 6.12 shows the results of the Wooden A and the Action Figure reconstructed at the same time. As we can see, both meshes were independently reconstructed, obtaining proper results. 64 6.1.4 Our method compared to Kinect Fusion Microsoft provides a software for mesh reconstruction as part of the Kinect SDK, called Kinect Fusion. This programs applies an iterative alignment between frames to reconstruct the mesh. Since it uses all of the captured points, it requires to rotate the sensor instead of rotating the object; otherwise the alignment between frames would fail. We tested the software in order to compare it with our scanning system. Since both techniques apply a voxelization over the 3D domain, we tried to adjust parameters to provide similar voxel resolutions in both methods. As Figure 6.13 shows, the quality of our results and the ones obtained using Kinect Fusion is similar. Depending on the region of the object, one of the mesh reconstructions is better than the other. Considering the region between the hair and the hat of the Doll, our software obtained a better result, since the contour of the hair is more detailed, being more similar to the real object. On the other hand, Kinect Fusion reconstructed better the bottom part of the Doll’s dress, probably due to the fact that we have no points below the dress in our captured data, and thus the normals’ estimation fails in that region. Figure 6.13: Left: Results of our scanning system. Right: Results of Kinect Fusion. 65 One advantage of our system is the fact that it works automatically in an enclosed small space: the sensor is fixed and the turntable does not move, it only spins. On the other hand, Kinect Fusion needs to move the sensor around the object, which requires a larger space in order to automatize the process. Since we have no system for moving the Kinect around the object, we needed to move it manually. Moreover, Kinect Fusion does not separate the object from the environment; it scans everything at once. If we want to acquire a 3D model of a certain object, we need to separate it from its environment. In addition, our system ensures a watertight mesh as output, while Kinect Fusion provides a partially reconstructed mesh without any guaranteed properties. 6.1.5 Final results Aside from all of the aforementioned experimentation, we also captured data from the different objects presented in Figure 6.1 and we reconstructed their meshes. Figures 6.14 and 6.15 show the best obtained results for each one of the real objects. Table 6.1 shows the mesh characteristics, the corresponding parameters, and the time spent in the mesh processing phase (i.e. the scan merge, the voxelization, the candidates’ selection, and the mesh reconstruction steps). It is important to state that the source code has not been carefully optimized, thus the reported timings could be improved. 66 Figure 6.14: Results obtained with the Poisson reconstruction (part 1). From left to right: voxelization, centroids candidates, centroids mesh, sub-voxel counting candidates, sub-voxel counting mesh. 67 Figure 6.15: Results obtained with the Poisson reconstruction (part 2). From left to right: voxelization, centroids candidates, centroids mesh, sub-voxel counting candidates, sub-voxel counting mesh. As seen in Figures 6.14 and 6.15, the obtained meshes approximate quite well the real objects, preserving a considerable amount of features, but losing those which were small, such as the details of the doll’s head (they could probably be preserved with more precise scanners). Furthermore, this results confirm that the sub-voxel counting technique obtains more detailed meshes than the centroids method. 68 Wooden A Doll Cylinder Action figure Virgin Amount of scans 920 935 930 933 910 Number of points 2.28M 5.73M 2.84M 2.96M 10.53M Number of candidates 11,217 / 14,157 24,171 / 30,342 11,590 / 14,056 8,245 / 9,881 29,785 / 35,428 Output mesh vertices 3,130 / 6,904 20,498 / 21,300 8,366 / 8,580 5,032 / 7,330 25,472 / 28,019 Output mesh triangles 6,260 / 13,799 40,992 / 42,596 16,728 / 17,156 10,059 / 14,649 50,931 / 56,024 Voxelization size 150 150 150 150 150 Voxelization threshold 90 90 100 120 120 Normals’ estimation neighbourhood radius 0.05 0.06 0.05 0.05 0.18 Poisson tree depth 8 8 8 8 8 Poisson samples per node 5 4 3 4 4 Poisson scale 1.6 2.5 1.8 1.9 1.5 Scans merge time 18.2898 46.5359 23.1975 24.9143 84.2766 Voxelization and candidates’ selection time 0.5131 / 1.1006 1.2487 / 2.8199 0.6421 / 1.4292 0.6471 / 1.2978 1.6182 / 4.6286 Mesh reconstruction time 1.4737 / 1.9642 4.0760 / 4.9556 1.9722 / 2.1963 1.5747 / 1.9492 7.7220 / 10.4728 Total time 20.2766 / 21.3546 50.8606 / 53.3114 25.8118 / 26.8230 27.1361 / 28.1613 93.6168 / 99.3780 Table 6.1: The results’ information. For X / Y cells: X means centroids method, and Y means sub-voxel counting. Time in seconds. The tests where done using a computer with an Intel i7-4770k and 8GB of RAM. As Table 6.1 shows, our scanning pipeline substantially reduced the size of the point cloud: the mesh reconstruction step receives some thousands of points instead of the millions of points of the input. We can also observe that the parameters vary from 69 •Repeat the scanning process for the non-scanned parts of the object and apply registration algorithms to merge all the point clouds, such as 4-PCS. •Add more scanners to the scanning system. With scanners in different positions we could reduce the under-scanned regions (together, they will see a larger portion of the object surface). The overall process would be the same, but instead of calibrating a single scanner, the user should calibrate all of them. •Add more scanners to the 3D scanner simulator (in order to evaluate previous suggestion in our virtual system). •Add different types of noise to the simulator. This would allow us to experiment with more realistic data. 76 Appendix A Kinect 2 specifications The Kinect 2 or Kinect for Xbox One is a device which provides colour and depth information in real-time. The main improvement with respect to the previous version of Kinect is the change in the scanning technology, which has changed from structuredlight to time-of-flight. The TOF technology measures the distance based on the speed of light, while structured light infers the distance by the way a structured pattern of light is deformed over a surface. Figure A.1: The Kinect for Xbox One device. 77 This device is also capable of tracking the body skeleton of up to 6 people, where each skeleton is formed by 25 joints. Table A.1 shows the technical hardware specifications of the device. Attribute Value Color image resolution 1920x1080 Depth image resolution 512x424 IR image resolution 512x424 Minimum distance 0.5m Maximum distance 4.5m Vertical field of view 60◦ Horizontal field of view 70◦ Frame rate 30 FPS Audio format 16kHz, 24-bit PCM Audio input 4 microphones, with 24-bit ADC converter, echo cancellation and noise supression. Connectivity USB 3.0 Table A.1: Kinect 2 technical specifications. 78 Bibliography [AMCO08] Dror Aiger, Niloy J. Mitra, and Daniel Cohen-Or. 4-Points Congruent Sets for Robust Pairwise Surface Registration. In ACM SIGGRAPH 2008 Papers, pages 85:1–85:10. ACM, 2008. [BM92] Paul J. Besl and Neil D. McKay. A Method for Registration of 3-D Shapes.IEEE Transactions on Pattern Analysis and Machine Intelligence, 14(2):239–256, February 1992. [BTP13] Sofien Bouaziz, Andrea Tagliasacchi, and Mark Pauly. Sparse Iterative Closest Point. In Proceedings of the Eleventh Eurographics/ACMSIGGRAPH Symposium on Geometry Processing, SGP ’13, pages 113–123, Aire-la-Ville, Switzerland, Switzerland, 2013. Eurographics Association. [BTS+16] Matthew Berger, Andrea Tagliasacchi, Lee M. Seversky, Pierre Alliez, Ga¨el Guennebaud, Joshua A. Levine, Andrei Sharf, and Claudio T. Silva. A Survey of Surface Reconstruction from Point Clouds.Computer Graphics Forum, 2016. [Cen08] A. Censi. An ICP variant using a point-to-line metric. In Robotics and Automation, 2008. ICRA 2008. IEEE International Conference on, pages 19–25, May 2008. [CM91] Y. Chen and G. Medioni. Object modeling by registration of multiple range images. In Robotics and Automation, 1991. Proceedings., 1991 79 IEEE International Conference on, volume 3, pages 2724–2729, April 1991. [DLB+16] Shaoyi Du, Juan Liu, Bo Bi Bi, Jihua Zhu, and Jianru Xue. New iterative closest point algorithmfor isotropic scaling registration ofpoint sets with noise.Journal on Visualization and Image Representation, 38:207–216, July 2016. [DZYW07] S. Du, N. Zheng, S. Ying, and J. Wei. ICP with Bounded Scale for Registration of M-D Point Sets. In 2007 IEEE International Conference on Multimedia and Expo, pages 1291–1294, July 2007. [EKS83] H. Edelsbrunner, D. Kirkpatrick, and R. Seidel. On the shape of a set of points in the plane.IEEE Transactions on Information Theory, 29(4):551–559, July 1983. [FB81] Martin A. Fischler and Robert C. Bolles. Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography.Commun. ACM, 24(6):381–395, June 1981. [Fit03] A. W. Fitzgibbon. Robust Registration of 2D and 3D Point Sets. In Proceedings of the British Machine Vision Conference, volume 21, pages 1145–1153, January 2003. [GCSA13] Simon Giraudot, David Cohen-Steiner, and Pierre Alliez. Noise-adaptive Shape Reconstruction from Raw Point Sets. In Proceedings of the Eleventh Eurographics/ACMSIGGRAPH Symposium on Geometry Processing, SGP ’13, pages 229–238, Aire-la-Ville, Switzerland, Switzerland, 2013. Eurographics Association. [GMGP05] Natasha Gelfand, Niloy J. Mitra, Leonidas J. Guibas, and Helmut Pottmann. Robust Global Registration. In Proceedings of the Third Eurographics Symposium on Geometry Processing, SGP ’05. Eurographics Association, 2005. 80 [GP02] S´ebastien Granger and Xavier Pennec. Multi-scale EM-ICP: A Fast and Robust Approach for Surface Registration. In Proceedings of the 7th European Conference on Computer Vision-Part IV, ECCV ’02, pages 418–432, London, UK, UK, 2002. Springer-Verlag. [Hor84] B. K. P. Horn. Extended Gaussian images.Proceedings of the IEEE, 72(12):1671–1686, December 1984. [JH03] T. Jost and H. Hugli. A multi-resolution ICP with heuristic closest point search for fast and robust 3D registration of range images. In 3-D Digital Imaging and Modeling, 2003. 3DIM 2003. Proceedings. Fourth International Conference on, pages 427–433, October 2003. [KBH06] Michael Kazhdan, Matthew Bolitho, and Hugues Hoppe. Poisson Surface Reconstruction. In Proceedings of the Fourth Eurographics Symposium on Geometry Processing, SGP ’06, pages 61–70, Aire-la-Ville, Switzerland, Switzerland, 2006. Eurographics Association. [LC87] William E. Lorensen and Harvey E. Cline. Marching Cubes: A High Resolution 3D Surface Construction Algorithm. In Proceedings of the 14th Annual Conference on Computer Graphics and Interactive Techniques, SIGGRAPH ’87, pages 163–169, New York, NY, USA, 1987. ACM. [MAM14] Nicolas Mellado, Dror Aiger, and Niloy J. Mitra. Super 4PCS Fast Global Pointcloud Registration via Smart Indexing.Computer Graphics Forum, 33(5):205–215, 2014. [MGP11] H. Men, B. Gebre, and K. Pochiraju. Color point cloud registration with 4D ICP algorithm. In Robotics and Automation (ICRA), 2011 IEEE International Conference on, pages 1511–1516, May 2011. [MLM05] J. Minguez, F. Lamiraux, and L. Montesano. Metric-Based Scan Matching Algorithms for Mobile Robot Displacement Estimation. In Proceedings of the 2005 IEEE International Conference on Robotics and Automation, pages 3557–3563, April 2005. 81 [MPD06] Ameesh Makadia, Alexander IV Patterson, and Kostas Daniilidis. Fully Automatic Registration of 3D Point Clouds. In Proceedings of the 2006 IEEE Computer Society Conference on Computer Vision and Pattern Recognition - Volume 1, CVPR ’06, pages 1297–1304, Washington, DC, USA, 2006. IEEE Computer Society. [MS10] A. Myronenko and X. Song. Point Set Registration: Coherent Point Drift.IEEE Transactions on Pattern Analysis and Machine Intelligence, 32(12):2262–2275, December 2010. [MSD04] A. Makadia, L. Sorgi, and K. Daniilidis. Rotation estimation from spherical images. In Pattern Recognition, 2004. ICPR 2004. Proceedings of the 17th International Conference on, volume 3, pages 590–593, August 2004. [RMB+08] R. B. Rusu, Z. C. Marton, N. Blodow, M. Dolha, and M. Towards 3D Point cloud based object maps for household environments.Robotics and Autonomous Systems, 56(11):927–941, 2008. [RY07] Laurent Rineau and Mariette Yvinec. A Generic Software Design for Delaunay Refinement Meshing.Comput. Geom. Theory Appl., 38(1-2):100– 110, September 2007. [SHT09] A. Segal, D. Haehnel, and S. Thrun. Generalized-ICP. In Proceedings of Robotics: Science and Systems, Seattle, USA, June 2009. [ZB11] L. Zeming and H. Bingwei. A Curvature-based Automatic Registration Algorithm for the Scattered Points. In 2011 Third International Conference on Measuring Technology and Mechatronics Automation, volume 1, pages 28–31, Jan 2011. 82