Mars 2017
Reconstruction of trees from 3D point clouds
Martin Stålberg
Teknisk- naturvetenskaplig fakultet UTH-enheten
Besöksadress:
Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0
Postadress:
Box 536 751 21 Uppsala
Telefon:
018 – 471 30 03
Telefax:
018 – 471 30 00
Hemsida:
http://www.teknat.uu.se/student
Reconstruction of trees from 3D point clouds
Martin Stålberg
The geometrical structure of a tree can consist of thousands, even millions, of branches, twigs and leaves in complex arrangements. The structure contains a lot of useful information and can be used for example to assess a tree's health or calculate parameters such as total wood volume or branch size distribution. Because of the complexity, capturing the structure of an entire tree used to be nearly impossible, but the increased availability and quality of particularly digital cameras and Light Detection and Ranging (LIDAR) instruments is making it increasingly possible.
A set of digital images of a tree, or a point cloud of a tree from a LIDAR scan, contains a lot of data, but the information about the tree structure has to be extracted from this data through analysis.
This work presents a method of reconstructing 3D models of trees from point clouds. The model is constructed from cylindrical segments which are added one by one. Bayesian inference is used to determine how to optimize the parameters of model segment candidates and whether or not to accept them as part of the model.
A Hough transform for finding cylinders in point clouds is presented, and used as a heuristic to guide the proposals of model segment candidates.
Previous related works have mainly focused on high density point clouds of sparse trees, whereas the objective of this work was to analyze low resolution point clouds of dense almond trees. The method is evaluated on artificial and real datasets and works rather well on high quality data, but performs poorly on low resolution data with gaps and occlusions.
ISSN: 1401-5757, UPTEC F 17011 Examinator: Tomas Nyberg Ämnesgranskare: Thomas Schön Handledare: James Underwood
Popul¨ arvetenskaplig sammanfattning
Tr¨ ad spelar en viktig roll i ekosystem v¨ arlden ¨ over, och de har utgjort en viktig resurs f¨ or m¨ anniskan genom m¨ ansklighetens historia. Anv¨ andningsomr˚ adena f¨ or tr¨ ad ¨ ar m˚ anga och olika, men en aspekt som ¨ ar viktig i m˚ anga av fallen ¨ ar tr¨ adens fysiska form.
Inom skogsbruk spelar tr¨ adens form en avg¨ orande roll. Avverkningsbeslut baseras p˚ a det och s˚ agverk kapar och klyver stockar enligt deras form. Inom jordbruk anv¨ ands tr¨ ad som producerar frukt och n¨ otter, och dessas former modifieras genom besk¨ arning f¨ or att
¨
oka sk¨ ord och f¨ orb¨ attra v¨ axtf¨ orh˚ allanden.
Ett tr¨ ads fysiska form utg¨ ors av tusentals eller miljontals stammar, grenar, kvistar och l¨ ov arrangerade i komplexa strukturer. Denna komplexitet g¨ or det sv˚ art att m¨ ata och representera den exakta strukturen av ett tr¨ ad. N¨ ar siffror p˚ a kvantitativa v¨ arden beh¨ ovs anv¨ ands i dagsl¨ aget d¨ arf¨ or ofta uppskattningar baserade p˚ a ¨ ogonm˚ att eller en- klare m¨ atningar, s˚ asom stamdiameter vid br¨ osth¨ ojd.
Om en 3D-modell av ett tr¨ ad av intresse fanns att tillg˚ a kunde m˚ att f˚ as snabbt, nog- grannt och automatiskt fr˚ an ber¨ akningar p˚ a modellen ist¨ allet f¨ or att m¨ ata p˚ a det fysiska tr¨ adet. Hur s˚ adana 3D-modeller kan skapas med noggrannhet och effektivitet ¨ ar ett p˚ ag˚ aende forskningsomr˚ ade.
I detta arbete presenteras en metod f¨ or att skapa 3D-modeller av tr¨ ad utifr˚ an data i form av s˚ a kallade punktmoln, vilka kan skapas med exempelvis LIDAR-instrument. LIDAR st˚ ar f¨ or Light Detection and Ranging, och simplifierat kan man s¨ aga att en LIDAR sveper en laserstr˚ ale fram och tillbaka ¨ over ytan p˚ a ett objekt f¨ or att noggrannt m¨ ata positionen av ett stort antal punkter p˚ a objektets yta. En m¨ angd s˚ adana punktpositionsm˚ att kallas f¨ or ett punktmoln. F¨ or att med h¨ og detaljniv˚ a f˚ anga strukturen av ett helt tr¨ ad kr¨ avs ett punktmoln med flera hundra tusen eller miljontals individuella punktpositionsm˚ att.
Fr˚ an ett s˚ adant punktmoln av ett visst tr¨ ad ˚ aterskapar den presenterade metoden en 3D-modell av tr¨ adet, best˚ aende av cylindriska segment. Rekonstruktionen utg˚ ar ifr˚ an ett enskilt segment som representerar en del av stammen vid tr¨ adets bas, och l¨ agger sedan till ett segment i taget p˚ a s˚ a s¨ att att modellen passar b¨ attre och b¨ attre med punktmolnet f¨ or varje segment som l¨ aggs till. Nya segment kan antingen l¨ aggas till vid
¨
andarna av existerande segment, som forts¨ attningar av stammar eller grenar, eller p˚ a sidorna som nya grenar.
iii
Exakt vilka segment modellen ska ut¨ okas med, var och med vilken storlek och orientering
¨
ar k¨ arnan i rekonstruktionsproblemet. Problemet f¨ orenklas av att inte behandla hela punktmolnet p˚ a en g˚ ang utan att ist¨ allet titta p˚ a mindre bitar av det i taget. P˚ a grund av det stora antalet punkter ¨ ar det mycket tids¨ odande att utf¨ ora ber¨ akningar p˚ a hela punktmolnet. Mindre delar har f¨ arre antal punkter som ofta ocks˚ a ¨ ar arrangerade i mindre komplexa strukturer.
Varje nytt segment som l¨ aggs till modellen ¨ ar ett resultat av vad som kan liknas vid en serie t¨ arningskast. Varje t¨ arningskast best¨ ammer en viss storlek, position och orienter- ing av det nya segmentet, och efter ett visst antal t¨ arningskast v¨ aljs den b¨ asta av de konfigurationer som genererats av t¨ arningskasten. Om modellen ¨ ar mer lik punktmol- net med det nya segmentet ¨ an utan det nya segmentet s˚ a l¨ aggs det nya segmentet till modellen, annars inte.
Huruvida modellen ¨ ar b¨ attre med eller utan ett visst segment best¨ ams av tv˚ a faktorer.
Den f¨ orsta faktorn ¨ ar hur v¨ al modellen st¨ ammer ¨ overens med punktmolnet, hur sanno- likt det ¨ ar att punktmolnet skulle se ut som det g¨ or, om tr¨ adet s˚ ag ut exakt s˚ a som modellen g¨ or. Den andra faktorn ¨ ar hur v¨ al modellen uppfyller f¨ orubest¨ amda antagan- den om hur tr¨ ad ¨ ar uppbyggda. Vissa s˚ adana antaganden ¨ ar simpla, s˚ asom att grenar inte kan genomkorsa varandra och att de generellt sett blir kortare n¨ ar de str¨ acker sig upp˚ at och ut˚ at. Andra antaganden ¨ ar mer subtila och specifika f¨ or vissa arter eller typer av tr¨ ad, s˚ asom antaganden om vilka vinklar som kan f¨ orv¨ antas mellan grenar vid f¨ orgreningspunkter.
Den presenterade metoden utv¨ arderades p˚ a punktmoln fr˚ an b˚ ade m¨ atningar p˚ a verkliga
tr¨ ad och simulerade m¨ atningar p˚ a artificiella tr¨ ad. P˚ a punktmoln med h¨ og uppl¨ osning
och t¨ ackning fungerade metoden v¨ al. P˚ a punktmoln med l˚ ag uppl¨ osning och d˚ alig
t¨ ackning, d¨ ar delar av tr¨ adet d¨ oljs bakom grenar och stammar, fungerade metoden
d˚ aligt.
Thank you James and Thomas for your help, guidance, immense patience and the op- portunity to do this thesis at the Australian Centre for Field Robotics. It was a fantastic experience, even though the project turned out to be even more difficult than expected.
My warmest thanks to the people of the ACFR and the Treehouse for your hospitality and friendliness.
Uppsala, February 2017 Martin St˚ alberg
v
Popul¨ arvetenskaplig sammanfattning iii
Acknowledgements v
Contents vi
1 Introduction 1
1.1 Background . . . . 2
1.2 Related Work . . . . 3
1.3 Contributions . . . . 4
1.4 Outline . . . . 5
2 Reconstruction algorithm 6 2.1 Premises and assumptions . . . . 6
2.1.1 Unstructured point cloud . . . . 6
2.1.2 Color and laser return intensity . . . . 6
2.1.3 Initialization . . . . 7
2.1.4 Computational performance of implementation . . . . 7
2.2 Motivation . . . . 7
2.3 Adaptation to point clouds . . . . 8
2.4 Inference method . . . . 8
2.5 Model . . . . 8
2.6 The reconstruction algorithm . . . 10
2.6.1 Inference conducted on subsets . . . 10
2.6.2 Assigned and unassigned data points . . . 11
2.6.3 Initialization of new model segments . . . 13
2.6.4 Proposal selection . . . 14
2.6.5 Backwards extension of model segments . . . 15
3 Proposals for model segment candidates 17 3.1 Partial proposals . . . 17
3.1.1 Radius change . . . 18
3.1.2 Length change . . . 18
3.2 Full proposals . . . 18
3.2.1 RANSAC . . . 18
3.2.2 Alignment . . . 20
vi
3.3 Local Hough transform for cylinders . . . 21
3.3.1 Spherical grid through icosahedron subdivision . . . 22
3.3.2 Algorithm . . . 23
3.3.3 Test results . . . 24
4 Evaluation of proposals 29 4.1 Likelihood . . . 29
4.1.1 Data points and model points . . . 30
4.1.2 Removing occluded model points . . . 30
4.1.3 Likelihood function . . . 31
4.1.4 Reduction of computational cost of distance searches . . . 32
4.2 Priors . . . 32
4.2.1 Angle . . . 33
4.2.2 Radius . . . 34
4.2.3 Intersection . . . 34
4.2.4 Adjacency . . . 36
5 Experimental results 37 5.1 Generating artificial data . . . 37
5.1.1 Arbaro and POV-Ray . . . 37
5.1.2 Encoding and decoding coordinates . . . 38
5.2 Results . . . 40
5.2.1 Results on artificial data . . . 40
5.2.2 Results on real data . . . 42
6 Conclusion and future work 46 6.1 Conclusion . . . 46
6.2 Future work . . . 46
A Artificial data 49 A.1 Arbaro . . . 49
A.2 POV-Ray . . . 49
A.3 Matlab . . . 52
Bibliography 54
Trees play an important role in ecosystems all over the world and they have been an important resource for humans throughout human history. The various uses for trees are many and diverse but one aspect of trees that is important in many cases is their physical structure.
In forestry, tree geometry is crucial. Harvesting decisions and yield estimates are based on it and sawmills cut logs according to it. In agriculture, trees that produce fruit and nuts have their structures altered through the removal of branches, a process called pruning, to improve yield.
The complexity of trees make their structures difficult to measure and represent. There- fore, when quantitative qualities of the structure of a tree are of importance, they are often not measured directly. Instead, they are estimated, for example by human visual inspection or from statistical models with simple measurements as input, such as stem diameter at breast height.
If a 3D model of a tree of interest was available measurements could be taken accurately, quickly and automatically using the model instead of on the physical tree itself. How to accurately and efficiently construct such models is an ongoing research topic and what this work is concerned with.
This work was conducted at the Australian Centre for Field Robotics (ACFR) at the University of Sydney, where there are ongoing research efforts about the use of robotics in agriculture. One particular agricultural environment that is being investigated is orchards. If a method for reliably and accurately reconstrucing 3D models of trees is developed, one of the potential uses of such models is for calculating which branches to prune to optimise the structure, growth and yield of the trees. The pruning could then potentially be performed by robots, relieveing humans of a physically demanding task and reducing the need for human labour in orchards.
In Section 1.1 further background information about the topic is presented. Section 1.2 provides an overview of related work. The chapter ends with a summary of the contri- butions of this work and an outline of the rest of this document.
1
1.1 Background
The geometrical structure of a tree can consist of thousands, even millions, of branches, twigs and leaves in complex arrangements. The structure contains a lot of useful infor- mation and can be used for example to assess a tree’s health or calculate parameters such as total wood volume or branch size distribution. However, because of the immense complexity, capturing the structure of an entire tree used to be nearly impossible. Early attempts involved manual digitization, using local positioning devices placed at differ- ent positions on the surface of the tree in succession to map one 3D point at a time.
This was very time consuming and difficult to carry out in the field [1]. The increased availability of computational power and technologies including digital cameras and Light Detection and Ranging (LIDAR) instruments facilitate the task of capturing tree struc- ture. Such technologies have, in particular during the last decade, sparked an increase in the research efforts concerning digitization of tree geometry (see e.g. [2] and [3]).
Advantages of the photography approach are that high-resolution cameras provide dense, instant sampling including color information. However, no depth information is recorded.
Photos from different viewpoints can be used jointly to infer depth, possibly requiring painstaking parameter tuning, but the result may turn out to be very noisy because of the irregular, complex and self-similar tree structure. The main advantage of the laser scanning approach is the direct capture of 3D data. Many current laser scanning systems do not capture any color information, however, other than possibly the return intensity of the laser. Also, producing a dense scan can require significantly more time than taking a photo. The time requirement can cause noise in the data if wind moves the tree during the data acquisition. Because of self occlusion (parts of the tree hiding other parts of the tree from the sensor’s view) scans or photographs from multiple viewpoints are generally required to thoroughly capture the tree structure. The occlusion problem is amplified greatly by the presence of foliage, meaning that data should preferably be acquired during the winter season when foliage is at a minimum. Of course, this only applies to deciduous trees, those that shed their leaves during fall. Acquiring good data on evergreen trees is likely to be equally difficult during all seasons.
Both photography and laser scanning, especially the latter, produce large amounts of raw data which needs to be processed or analyzed to extract the structural information.
This involves identifying and characterizing branches and how they are connected. The result can be represented as a 3D model with small triangular faces closely matching the surface of the tree, or a collection of geometric primitives, such as cylinders and cones.
Yet another option is to use a skeleton, a small set of connected points describing the
topology of the tree.
1.2 Related Work
Most previous work focuses on the reconstruction of the trunk and major branches as data on twigs and leaves is generally too noisy to allow reconstruction with acceptable accuracy [2]. Some of the literature on the topic focuses on the production of 3D models for graphical purposes, e.g. for computer games or other virtual environments, and is more concerned with the visual appeal of the models than how accurately they correspond to the real world originals. Extrapolation or artificial synthesis is commonly used to replace parts of the tree that were missing or excessively noisy in the acquired data [4]. Several works use the fact that a pair of points on a tree that are close to each other are likely to be on the same branch.
In [4] a connectivity graph is created by letting each point in the LIDAR point cloud be represented by a vertex and connecting each vertex with the vertices whose points are among the k nearest neighbours of the first point. The edges are weighted with the distances between the points. A root point (the point chosen to represent the root of the tree) is then selected and a minimum weight spanning tree is constructed from the connectivity graph using Dijkstra’s algorithm [5]. The points are clustered based on the quantized path length through the spanning tree to the root node. The centroids of the clusters are connected to each other according to the neighbouring information of their constituent points to form a skeleton of the tree.
A similar approach to [4] is taken in [2], but the spanning tree is taken directly as the initial skeleton. The skeleton is refined by smoothing and contracting the point cloud, then removing redundant points deemed to contribute negligibly to the overall structure of the skeleton. The method described in [6] improves on the one in [4] by introducing a more sophisticated clustering method. The point cloud is repeatedly randomly clustered with a k-means algorithm and clusters that have a cylindrical enough shape are saved and their points removed from the point cloud before the next iteration of clustering.
The resulting cylindrical clusters are hopefully segments of branches and are connected to produce the tree skeleton.
This work is primarily based on the paper Modeling of Tree Branches by Bayesian Network Structure Inference by Ma et al. [7], where wide baseline stereo image pairs are used to infer tree structure, modeled as a hidden Bayesian network. The tree model is built in stages by generating hypotheses at each stage for the on-growing of current branches and generation of new branches, then accepting those hypotheses that best fit the data and prior assumptions about properties such as branch smoothness and radii.
The end result is a skeleton of the tree consisting of points called skeleton nodes, or
just nodes, and connections between them. The first stage is initialized manually by
specifying the skeleton root node with a position centered in the base of the tree, with
a growing direction and a radius. In each subsequent stage child nodes for each parent
node, one at a time, are hypothesized and evaluated according to the properties of the
parent node and how well they match the stereo images. The child nodes are inferred
using the Metropolis-Hastings algorithm [8], the driving proposals for which are based on the properties of the parent node.
In [9] the point cloud is rasterized and centroids of adjacent voxels are connected with an edge if the points contained in the respective voxels are distributed in a fashion that indicates a connection of the underlying structure. The resulting graph is then refined, e.g. by removing any loops, and finally taken as the skeleton of the tree.
The method in [10] produces a tree skeleton by contracting the point cloud using con- strained Laplacian smoothing. An initial skeleton is produced by connecting each point to its closest neighbours and the skeleton is then refined by removing loops and points that do not contribute much to the overall skeleton structure.
In [11] a method first developed to produce artifical trees, called a space colonization algorithm [12], is applied to the point cloud resulting from the scan of a real tree. The algorithm grows a skeleton of the tree in a bottom-up manner from a manually selected root node. Points function as attractors and the placement of new skeleton nodes is determined by the points within an area of influence from the current existing skeleton node. The method is similar to the one in [7] but less tunable.
The previous works often focus on sparse trees and rely on high density point clouds with little noise produced through long duration scans from static platforms and combining scans from several view points to mitigate occlusion (see e.g. [3]). In essence, they often rely on high quality data, whereas the data that is the target of this work is low density point clouds of dense trees, i.e. poor data.
1.3 Contributions
The ACFR has performed scans of an almond orchard with a LIDAR of the LMS5xx family from SICK Sensor Intelligence. The main objective of this work was to develop a method able to reconstruct accurate models of the almond trees in this data to use for further research towards autonomous pruning and other applications. Characteristic of this data is that the number and density of branches is high and the resolution of the data is low.
The developed method uses Bayesian inference to reconstruct models of trees from point clouds. The goal of producing accurate models from poor data was not achieved, but the method has some success on high and medium quality artificial datasets.
Additionally a Hough transform for finding cylinders in point clouds was developed, and
a method for producing simulated LIDAR scans of computer generated tree models. The
latter can be used during development to compare results on simulated point clouds with
the underlying ground truth data.
1.4 Outline
A method for reconstructing models of trees from point clouds is presented in Chapter 2.
The model consists of cylindrical segments which are added iteratively to the model.
Chapter 3 describes how the model segment candidates are initialized and modified prior
to being accepted or rejected as part of the model. Chapter 4 details how model segment
candidates are evaluated to determine which modifications to accept and whether to
accept them as part of the model. Results of the method on some datasets are presented
in Chapter 5. Finally conclusions are presented in Chapter 6.
This chapter gives a high level overview of the structure of the algorithm used for reconstruction. Many of the details are deferred to later chapters.
The chapter begins with Section 2.1, presenting some premises for this work and some important assumptions that were made. Section 2.2 presents the motivation for bas- ing this work on a particular related work, then Section 2.3 and Section 2.4 explain adaptations made according to the particulars of the problem investigated in this work.
Section 2.5 presents how the tree models are structured, and finally Section 2.6 gives an overview of the reconstruction algorithm.
2.1 Premises and assumptions
The premises for this work and some important assumptions are detailed in this section.
2.1.1 Unstructured point cloud
When dealing with structured three dimensional data, such as depth images or volume images, the structure of the data can be leveraged for various purposes, e.g. performance gains for certain algorithms. This work assumes no such particular structuring of the data to enable using input data that is a composite of several scans of the same sensor from different viewpoints or even from different sensors. An example of such data, relevant to making this assumption, is data from a side scanning LIDAR moving past the target tree.
2.1.2 Color and laser return intensity
This work makes no use of any color or intensity information associated with the input point cloud. If accurate color information was available it could be used to classify the points into classes such as leaves, flowers and bark, which could guide the reconstruction algorithm. There are instruments that combine LIDAR and photography to co-register range and color data. This can be achieved for example with a LIDAR and camera
6
mounted on the same rotating sensor head which enables capture of range scans and images from the same viewpoint with very small time discrepancy. However, the accu- racy of such color information when capturing tree crowns is low due to the high spatial frequency and the small size of twigs and leaves.
Many LIDARs record the return intensity of the laser pulses. Although it is not as straightforward as with color information, intensity can also be used to distinguish be- tween different parts of the tree. The decision was made not to make use of intensity information, to not widen the scope of this work too much.
However, if a classifier that could use color and intensity information to classify points was developed, it could easily be used to increase the performance of the algorithm described in this work. It could simply be used in a preprocessing step to remove all points not classified as part of trunk, branches or twigs, providing the reconstruction algorithm with clearer data.
2.1.3 Initialization
The reconstruction algorithm has to be initialized with a starting point, including a position, radius, growing direction etc. For trees with sufficiently simple geometry, and with good scanning data, these initialization parameters could conceivably be derived rather easily through an automatic analysis of the point cloud. Such an analysis was not developed in this work, instead initialization parameters were acquired manually.
2.1.4 Computational performance of implementation
In many numerical methods, quality of the result is improved with the allowed runtime or the number of iterations. This is also the case with the method presented in this work. Critical during the development of this method has been the tradeoff between allowing more runtime to get better results, and allowing less runtime to get quicker feedback on the effects of changes to the method. To get more of both, computational performance has been an important consideration in the implementation, which was done in Matlab
1. Some of the measures taken to improve computational performance was the vectorization of loops and the reuse of calculation results, where possible.
2.2 Motivation
The approach in [7] by Ma et. al. was chosen as a basis for this work. Although it uses stereo images rather than point clouds the framework seemed sound and the details could be adapted to also work for point clouds. A particularly important advantage of an approach based on Bayesian inference is the ability to incorporate existing knowledge
1More information about Matlab at https://mathworks.com/products/matlab.html
about tree morpohology and growth patterns into the prior probability functions (see Chapter 4). In cases where the type or species of the tree subject to reconstruction is known one could easily tune the parameters of the prior probability functions to the particular characteristics of that species or type.
2.3 Adaptation to point clouds
As stated previously, this work is primarily based on [7], where stereo images are used as the input for the structural reconstruction (for further details please refer to the paper itself or a short introduction in Section 1.2). This work is concerned with point clouds, and therefore certain adaptations and extensions to the approach of Ma et al. where necessary. These adaptations and extensions are detailed here.
The original likelihood function is based on pixel intensities in the stereo images. The point clouds considered here do not have any intensity information
The algorithm essentially tries to locally optimally fit cylinders to the data points.
Essential for this is good initialization. The initial proposals, and possibly also later proposals, are generated from the data. The details are described in Chapter 3.
2.4 Inference method
Whereas Ma et. al. [7] use the Metropolis-Hastings algorithm [8] in the inferene of the tree structure, this work choses a different approach. With 8 scalar parameters per model segment (two 3D positions and two radii), inferring multiple model segments simultaneously was deemed too computationally expensive. Instead model segments are inferred one at a time, sequentially. Essentially, it was deemed easier to find several local optima in separate lower dimensional spaces than a single global optimum in a higher dimensional space.
2.5 Model
The tree is modeled as a set of tapering or straight cylindrical segments. An example
of such a model segment is shown in Figure 2.1, where the most important geometric
attributes of a model segment are also described. The points labeled position 1 and
position 2 are called the first and second endpoints respectively. Many other derivative
properties that are repeatedly used in calculations, such as the distance between the
endpoints, are precalculated and stored with the model segment for increased compu-
tational efficiency. The model is constructed in such a way that in any pair of a parent
segment and child segment, the second endpoint of the parent is closer to the first end- point of the child than the first endpoint of the parent is to the second endpoint of the child.
Apart from the geometric attributes, each model segment also keeps a reference to its parent segment and a list of references to all its child segments. These references store the hierarchical information of the model. The model segments can be grouped into branches according to this hiearchy. If a model segment is the first child of its parent model segment, it belongs to the same branch as the parent segment. Otherwise it is the first model segment of a new branch. A depiction of the hierarchical aspect of the model is shown in Figure 2.2.
Figure 2.1: The tree is modeled as a set of straight or tapering cylindrical segments.
An example of such a segment is shown in the figure above. The legend shows the most important attributes of the model segment.
Figure 2.2: A graph representation of a tree model consisting of model segments A through H. Each edge goes from a parent to a child and is labeled with the one- based child index. Branches made up of multiple model segments are represented by rectangles. Remaining nodes constitute their own single segment branches. Each child
segment with index 1 is part of the same branch as its parent.
2.6 The reconstruction algorithm
This section gives a high level overview of the reconstruction algorithm, detailing the control flow and how different parts are used together.
The top level control flow is detailed in Algorithm 1. The input to the algorithm is a point cloud and an initial model segment, called the root segment, and the main output is a set of model segments. The algorithm is centered around two sets of model segments called current stage and next stage. Current stage holds references to model segments whose child segments are inferred during the current iteration. For any new segments that are added to the model, references are also added to next stage. At the end of each iteration, the contents of current stage are replaced with the contents of next stage, and next stage is overwritten with the empty set.
For each model segment in current stage, inference of new child segments are repeated until the inference was unable to find a child segment that improved the model. This may happen after 0 or more child segments were created for that model segment.
This is repeated until current stage is empty, a set maximum runtime has been reached or a maximum number of segments in the model has been reached.
A high level description of the inference of new model segments is given in Algorithm 2 and Algorithm 3, with varying level of abstraction. In Algorithm 3 a description is given of how a new model segment is initialized. After initialization Algorithm 2 describes how the new segment is subject to repeated modification proposals, which are accepted if they increase the posterior probability of the model. The evaluation of proposals is described in detail in Chapter 4. Eventually, the new model segment is returned if the posterior probability of the model is higher with the new model segment than without it.
To prevent the listings of Algorithm 2 and Algorithm 3 to become overly long the level of abstraction was set high and many details omitted. These details are described in the following subsections.
2.6.1 Inference conducted on subsets
For performance reasons the inference does not consider the whole point cloud at once,
instead it works on subsets. The inference of child segments for a model segment s
considers only points within a distance d from the second endpoint of s, called e
2. The
distance d is given by d = max(d
min, 4r
2) where r
2is the radius of s at e
2. The factor
4 was chosen empirically, and the value for d
minwas chosen empirically based on the
resolution of the point cloud to tune the behaviour of the algorithm with small model
segments. The sphere centered at e
2with radius d is called the region of inference. Like-
wise, only a subset of the model is considered at once. This subset of the model is called
the collision set, and consists of the model segment s itself and all other model segments
Algorithm 1 The reconstruction algorithm, high level
1:
procedure reconstruct(p, s
0)
2:
S
c← {s
0} . Add the root segment to the current stage
3:
M ← {s
0} . Add the root segment to the model
4:
while S
c6= ∅ do . While there are segments in the current stage
5:
S
n← ∅ . Overwrite c
nwith empty set
6:
for each segment s
iin S
cdo . For every segment in the current stage
7:
while s ← inferChildSegment(p, s
i) do
8:
S
n← S
n∪ {s} . Add the new segment to next stage
9:
M ← M ∪ {s} . Add the new segment to model
10:
end while
11:
end for
12:
S
c← S
n. Overwrite contents of S
cwith S
n13:
end while
14:
return M
15:
end procedure
whose bounding spheres are contained within, or intersect, the region of inference. This is further explained in Chapter 4.
2.6.2 Assigned and unassigned data points
To indicate which part of the data is already covered by the model, each data point is associated with a flag, i.e. a boolean variable. If the flag for a data point p is set, i.e. its value is true, the data point is said to be assigned. Assigned points are ignored when computing full proposals (see Section 3.2).
Initially the flag of each data point is set to false. The flag of a data point p is set to true when a new model segment is added to the model, and p is inside of or close to the new model segment.
Given a model segment s with endpoints ~ A and ~ B with radii r
Aand r
Brespectively, define a polar coordinate system with origin at ~ A and positive z-axis passing through B. Given a point p, its r and z coordinates in this polar coordinate system are given by ~
z = ~ v • ~ q
r = |~ q − z~ v| (2.1)
where ~ q = ~ p − ~ A and ~ v is a unit vector along the line from ~ A to ~ B given by
~ v =
B− ~~lAl = | ~ B − ~ A| (2.2)
~
p is defined to be inside of or close to s if it is true that
Algorithm 2 The reconstruction algorithm, inference of child segments
1:
procedure inferChildSegment(p, s)
2:
Calculate the current (initial) posterior probability
3:
s
0← initializeChildSegment(p, s)
4:
5:
for integers i in range [1, n] do
6:
if i > 1 then
7:
if the posterior increased in the previous iteration then
8:
Modify s
0with the alignment proposal
9:
else
10:
Modify s
0with a randomly chosen proposal
11:
end if
12:
end if
13:
Move the first endpoint of s
0backwards to touch existing segments
14:
Calculate the posterior probability with s
015:
if the posterior decreased then
16:
Revert s
0to its state at the end of the last iteration
17:
end if
18:
end for
19:
20:
if the last posterior is greater than the initial posterior then
21:
Mark all points inside or close to s
0as assigned
22:
return s
023:
else
24:
Delete s
025:
return nothing
26:
end if
27:
end procedure
z ≥ −c z ≤ l + c r ≤ r
z+ c
(2.3)
where c is the fuzziness factor, which defines exactly how close ~ p must be to s to be considered to be close to s. The value given to c was 0.2r
A. The value of r
s, given by (2.4), is the radius of s at z, or the radius of s at whichever of ~ A and ~ B is closest to ~ p.
r
s=
r
A, if z < 0 r
B, if z > l
z
l
(r
B− r
A) + r
B, otherwise
(2.4)
A visual representation of this definition is provided in Figure 2.3.
Figure 2.3: For a model segment s (blue) with endpoints ~A and ~B with radii rA and rB respectively, points on or within the thick solid black line are defined to be inside
of or close to s. The value of c is 0.2rA.
2.6.3 Initialization of new model segments
As described in Algorithm 3, the initial parameters for a new model segment can be set in a number of ways. The parameters are obtained either from RANSAC (see Section 3.2.1 for more about RANSAC) or from the best Hough transform result (see Section 3.3 for more about the Hough transform used). The parameters to be set are the first endpoint, A, its radius r ~
A, the second endpoint, ~ B and its radius r
B.
Central to the choice is a set of radii values, R. Let the parent segment be called s, and its radius at its second endpoint be called r
s2. If r
s2is less than or equal to r
min, R has only one element, namely r
min. Otherwise R has 15 elements, linearly spaced in the range [r
min, 1.5r
s2]. r
minis the minimum branch radius and its value is set based on the characteristics of the tree being modeled and on the resolution of the point cloud.
If the cylinder returned by RANSAC has a radius that is less than or equal to the maximum element of R, it is used to set the initial values for the new model segment.
The endpoints ~ A and ~ B are given by
A ~ = min(D)~ V B ~ = max(D)~ V
D = {(~ q − ~ C) • ~ V | ~ q ∈ Q}
(2.5)
where Q is the inlier set, ~ V is a unit vector along the axis of the fitted cylinder and ~ C is a point on the axis of the fitted cylinder. To get the correct orientation of the new model segment, if ~ A is closer to the second endpoint of s than what ~ B is, ~ A and ~ B are swapped.
If the radius of the cylinder fitted by RANSAC greater than max(R), the new model segment is instead initialized using the best result from the Hough transform. The Hough transform is computed with R at the second endpoint of s, and if s has at least on child segment already, the Hough transform is also computed at 7 points, p
H, linearly spaced between the end points of the child segment according to
p
H= {
6iA + (1 − ~
6i) ~ B | i ∈ Z, 0 ≤ i ≤ 6}. (2.6) The first child segment of a model segment s is considered as the continuation of the branch partially modeled by s. Computing the Hough transform along the continuation of s is done to detect branches that branch off of the continuation.
A Hough transform result can be seen as a quadruplet of a number of votes, a grid point, a center point and a radius. The hough transform with the highest number of votes is selected and its grid point is assigned to ~ A, its center point to ~ B and its radius to r
Aand r
B.
The results from the Hough transform are kept and used for full proposals during the in- ference, regardless of whether the new model segment was initialized with the RANSAC result or the best Hough transform result.
Regardless of how the child segment is initialized, the Hough transform results are stored and used for generating change proposals during the inference.
2.6.4 Proposal selection
Algorithm 2 describes how at one point a model segment candidate is modified with a randomly chosen change proposal. The proposal options and their probabilities are listed in Table 2.1. Most details of the proposals are described in Chapter 3. The Hough transform results are computed during the initialization of the new child segment, and stored for generating change proposals. The next Hough transform result is the one with the highest number of votes that has not already been used for a proposal.
Table 2.1: Proposal options and their probabilities
Proposal Probability
Radius change 0.4
Length change 0.4
Next Hough transform result 0.2
Algorithm 3 The reconstruction algorithm, initialization of child segments
1:
procedure initializeChildSegment(p, s)
2:
e
2← secondEndPoint(s)
3:
p
0← unassigned points in p within distance d from e
24:
if too few points in p
0then
5:
return nothing
6:
end if
7:
8:
R ← radii lineraly spaced around radius at second endpoint of s
9:
Calculate Hough transform: Hough(p
0, e
2, R)
10:
if segment s has children then
11:
for points q
ispaced linearly between the endpoints of the first child of s do
12:
Calculate Hough transform: Hough(p
0, q
i, R)
13:
end for
14:
end if
15:
16:
s
0← new segment with s as parent
17:
18:
C ← RANSAC(p
0)
19:
if Radius of C less than maximum of R then
20:
Assign properties to s
0from C
21:
else
22:
Assign properties to s
0according to the best Hough transform result
23:
end if
24:
25:
return s
026:
end procedure
2.6.5 Backwards extension of model segments
Algorithm 2 describes how at one point a model segment candidate’s first endpoint is moved backwards to touch existing segments. This procedure is called backwards extension and its purpose is to prevent gaps in the model, i.e. between adjacent model segments.
A set of points P are calculated as
P =
A + ~ −2i
n − 1 ( ~ B − A) | i ∈ Z, 0 ≤ i ≤ n
(2.7)
where ~ A and ~ B are the first and second endpoints of the model segment candidate,
respectively, and the value of n is 50. The points P are spaced linearly along the
extended axis of the model segment candidate, with the first point coinciding with ~ A
and the last point being situated at a distance of 2l from ~ A in the opposite direction of
B, where l is the length of the model segment candidate. ~
Let ~ p
c∈ P be the point closest to ~ A that is inside of a model segment in the collision
set. If ~ p
cexists, ~ A is moved to ~ p
c, otherwise ~ A is moved to the second endpoint of the
parent model segment. Note that the parent segment is part of the collision set.
candidates
As described in Chapter 2 a structural model of a tree is built by an iterative procedure.
The model is built of model segments shaped as tapering or straight cylinders. The segments each model a part of a stem or branch. As the procedure constructs the model it creates model segment candidates which might be added to the model if they fit the data well enough. The model segment candidates are modified by proposals which change one or several of their attributes. A proposal is accepted if it increases the goodness of fit between the model and the data and vice versa. Finally, a model segment candidate is accepted as part of the model if it after a set number of proposals has increased the goodness of fit compared to what it was in the absence of the model segment candidate.
This chapter details the proposals and Chapter 4 details how goodness of fit is evaluated.
Partial proposals change one or a few attributes of a model segment candidate whereas full proposals change most of, or all of, the attributes. Section 3.3 describes a local Hough transform for cylinders. This Hough transform and RANSAC (described in Section 3.2.1) are used as heuristics for getting good initial sets of attributes which can then be improved by the partial proposals and the alignment proposal (described in Section 3.2.2).
3.1 Partial proposals
The partial proposals change only one or a few of the attributes of a model segment by small random perturbations to their current values. Some of these proposals will nudge the state of the model segment, the set of all current values of its attributes, towards a local optimum in the posterior probability distribution.
Two types of partial proposals are used. One changes radii and the other changes length.
A third type of partial proposal that moved either of the endpoints a small distance in a random direction was also implemented. However, it was abandoned as it was assessed that its contribution was already covered by the full proposals in combination with the alignment proposal.
17
3.1.1 Radius change
This partial proposal changes one of the radii of the model segment candidate. One of the ends is chosen at random, with equal probability, and its radius is changed to a new value selected randomly from a uniform distribution on the interval [0.7r
0, 1.1r
0), where r
0is the current radius. Note that the interval is not centered on the current radius but rather the proposals are biased towards smaller radii. This is to model the tendency of branches to get progressively thinner towards their tips.
3.1.2 Length change
This partial proposal changes the length of the model segment candidate by moving the second endpoint along the axis of the model segment. The new length is selected randomly from a uniform distribution on the interval [0.7l
0, 1.3l
0), where l
0is the current length.
3.2 Full proposals
The full proposals change most of or all of the attributes of a model segment candidate.
Unlike the partial proposals they are not stochastic but rather data driven, i.e. they are based on an analysis of the local data.
Three types of full proposals are used. RANSAC (RANdom SAmple and Consensus) is used to produce the initial state of a model segment candidate by attempting to find the largest subset of the local data that is well modeled by a cylinder. A local Hough transform for cylinders assesses how cylindrical different subsets of the local data appear to be in different directions, information which is then used to form proposals. The alignment proposal utilizes distance and neighbour information from the evaluation of the previous proposal to move the model segment candidate closer to its closest local data points.
3.2.1 RANSAC
RANSAC is an iterative method for estimating parameters of a mathematical model from data that contains outliers. The data points that are not outliers are called inliers.
Each iteration consists of two steps. The first step is to randomly select n data points, where n is the smallest number of data points sufficient for calculating the parameters of the model in question. The second step is to evaluate the model with those parameters.
The underlying assumption is that if all of the points randomly selected in an iteration
are inliers, the resulting model is highly likely to fit the inlier set well. This assumes, of
course, that the selected mathematical model is suitable.
After a predetermined number of iterations the algorithm stops and the parameters that resulted in the best fit are returned. RANSAC is a non-deterministic method and is therefore not guaranteed to produce good results, but the probability to do so increases with the number of iterations. Several algorithm parameters can also be tuned to further increase the probability. More in-depth knowledge about RANSAC can be found in [13].
In this work a cylindrical model is used. In each iteration three data points are selected at random. In the plane which they define, the center and radius of their circumcircle (the circle intersecting all three points) are calculated. The radius of the cylinder is defined as the radius of the circumcircle and the axis of the cylinder is defined by the center of the circumcircle and a normal to the plane. The error function, which evaluates the model, computes the error, E, for each point in the point cloud according to
X ~
r= X − ~ ~ P ,
X ~
t= X ~
r− ( ~ X
r• ~ V ) · ~ V , E = (| ~ X
t| − R)
2,
(3.1)
where ~ X is the point itself, ~ P is the circumcenter, ~ V is a unit vector orthogonal to the circumcircle and R is the radius of the circumcircle (see Figure 3.1 for illustration). Those points whose errors are below a certain threshold are said to constitute the consensus set.
The variety of RANSAC used here is called MSAC, M-estimator SAmple and Consensus.
In MSAC the score of the model is the sum of the errors of the points in the consensus set, as opposed to RANSAC where the score is simply the cardinality of the consensus set. The RANSAC Toolbox
1for Matlab was used to run RANSAC.
Figure 3.1: The parameters of a cylinder (light blue) are determined by three ran- domly selected points (magenta). The errors for other points (black) are the squares of their shortest distances to the cylinder (red). The axis of the cylinder (purple) is per- pendicular to, and passes through the center of, the circumcirle of the selected points
(green).
Ideally the model should be constructed such that, if disregarding noise, whenever all n randomly selected points are part of the inlier set the optimal set of model parameters
1By M. Zuliani. Available at https://github.com/RANSAC/RANSAC-Toolbox (2014-11-12)
for modeling the inliers will result. To produce good parameters for the model used here, the three selected points must not only be part of the inlier set but also in, or close to, a plane perpendicular to the axis of the cylinder on whose surface they are located. There are more sophisticated methods for calculating cylinder parameters from a minimal set of points that would be better suited to the task (see e.g. [14]), but these were not looked further into due to time constraints.
3.2.2 Alignment
Whenever a proposal of any kind is accepted, the next proposal is always the alignment proposal. The alignment proposal uses distance and neighbour information from the evaluation of the previous proposal to move the model segment candidate closer to its closest local data points. This is done by calculating new end points for the model segment candidate, a calculation which makes use of the model segment’s model points.
The model points are points spaced evenly in a grid pattern on the curved part of the surface of a model segment. They are described in detail in Section 4.1.1. Each valid model point of the model segment candidate is assigned two weights, w
1and w
2, according to its position, as illustrated by Figure 3.2. The weights are designed so that model points close to one of the endpoints have a higher impact on its new position than model points further away.
Figure 3.2: The weights w1and w2 are assigned to each model point according to its longitudinal distance along the axis of the model segment candidate, measured from
the first endpoint. The length of the model segment is denoted by L.
New positions for both of the model segment candidate’s endpoints are proposed ac- cording to
P ~
i= ~ P
i0+ 2 P
Nj
w
ij· ~ V
jP
Nj
w
ij, i ∈ {1, 2}, (3.2)
where P
iand P
i0are the new and current position of endpoint i, j enumerates the model
points, N is the number of model points and ~ V
jis a vector from model point j to its
nearest neighbour data point. The factor 2 was found empirically to increase the rate of
convergence, i.e. decrease the number of alignment proposals required to move the model segment candidate a certain distance towards a cylindrical point cloud. The reason this works can be understood intuitively by considering two parallel but offset cylinders and realizing that the closest point on the second cylinder from a point on the first cylinder is (almost) always closer than the corresponding point on the second cylinder (see Figure 3.3). A multiplicative factor can provide some measure of compensation for this.
Figure 3.3: The distance from a point P on circle A to the corresponding point on circle B (red line) is almost always greater than the shortest distance from P to any point on B (blue line). The only exception is when P is located at Q in which case the
two distances are equal (green line).
3.3 Local Hough transform for cylinders
A local Hough transform for cylinders was developed to provide coarse but good initial guesses for configurations for model segment candidates. The result is used for full proposals that set both radii and both endpoints of a model segment candidate.
The original Hough transform was developed for detecting lines in images. Extensions and generalizations have since been employed for detecting various kinds of shapes in various kinds of data. The central idea of the Hough transform is that each data point casts votes for all possible configurations of a target shape that could cause that point to exist. The more points a certain configuration explains, the more votes it receives, and the more likely it is to be the cause of the data. To make computation feasible only a finite number of configurations can be considered, as computational and memory costs increase with more configurations. The configurations are described by points in a discretized parameter space and the data points cast votes for these parameter space points, the voting pattern depending on the shape. The main drawback of the Hough transform is its high memory and computation costs for high-dimensional parameter spaces. The costs increase exponentially with the dimensionality of the parameter space making naive implementations intractable for shapes with anything more than a few parameters, assuming a reasonably high resolution in the parameter space.
Motivated by the assumption that tree branches are locally approximately cylindrical
a Hough transform for cylinders was developed and used in this work for detecting
directions and radii of branch segments. It is designed to detect only cylinders with a specific radius whose axes pass through one specific point, called the transform center.
This limitation cuts computation and memory costs significantly and is not detrimental to the performance of the model construction. As the model construction is sequential it always aims to find connecting branch segments for one model segment at a time.
As those are necessarily connected physically to the branch represented by the current model segment their axes pass through it. The transform also detects cylinders whose axes pass close to the transform center, although the directional result for such cylinders is biased. Thus, by performing the transform at several locations along the axis of the current model segment, for many radii at each location, any connecting branches can potentially be detected.
3.3.1 Spherical grid through icosahedron subdivision
This Hough transform is designed for point clouds of the sort resulting from a LIDAR scan. This means each point in the point cloud is, disregarding noise, located on the surface of the scanned object. Each point votes for a set of points in a parameter space, which is discretized by a set of grid points. These grid points are approximately uniformly spaced on the surface of a sphere and describe cylinder axes together with the transform center. The grid is created by an iterative process starting with a regular unit icosahedron, a polyhedron with twenty equilaterally triangular faces and twelve vertices, each vertex at unit distance from the origin
2. The triangular faces of the icosahedron are each subdivided into four smaller triangles and the newly created vertices are projected outwards to unit distance from the origin. This procedure of subdivision and projection is repeated until a sufficiently high resolution, and number of grid points, has been reached. See Figure 3.4 for an illustration of the procedure. The vertices of the final shape are taken as the grid points.
Figure 3.4: An approximately isotropic spherical grid can be created by a procedure of subdivision and projection starting with a regular icosahedron.
The resulting grid is approximately uniform, i.e. the spacing between any point and its neighbours is almost the same for every point, regardless of its position on the sphere. A perfectly uniform grid is not necessary but desirable, as the computational cost increases with the number of grid points. A uniform grid therefore minimizes the computational
2Implementation by Kurt. Available at http://au.mathworks.com/matlabcentral/fileexchange/
28842-gridsphere (2014-11-13)
cost for a particular lowest local grid resolution. In this work a grid of 2562 points was used, which is achieved after four iterations of the subdivision procedure. With resolution defined as the angle between neighbouring grid points, as measured from the sphere’s center, this gives a resolution in the interval [4.0
◦, 4.7
◦] with an average of 4.3
◦. As the procedure for generating the grid is rather costly it is run once at the start of reconstruction to produce a grid of unit radius which is then stored and used in all Hough transform computations.
3.3.2 Algorithm
Considering the case of just one single point, offset from the transform center, the task of the algorithm is to find cylinders of a specific radius with axes passing through the transform center, that explain this point, i.e. with their surfaces constrained to the point. No point on the surface of any such cylinder is closer to the transform center than one cylinder radius. Therefore the point must be at least one cylinder radius away from the transform center. The axes of the cylinders fulfilling the criteria form a cone with apex at the transform center. The intersection of the cone and the unit sphere form a circle and all grid points close enough to that circle, called the voting circle, get one vote each. See Figure 3.5 for an illustration. If the point cloud the transform is applied to contains multiple points on the surface of a cylinder passing through the transform center with a radius close to the radius used in the transform, the voting circles of the points will overlap at the intersection of the unit sphere and the cylinder axes and grid points close to that location will receive a large number of votes.
The criterion stating that grid points need to be close enough to the voting circle to get votes is more specifically defined as within a certain threshold angle of the voting circle. The threshold angle used is one half of the maximum angle between any two neighbouring grid points, which with 2562 grid points is approximately 2.4
◦. This value ensures that the voted for grid points are as close to the voting circle as possible while still forming a gapless circle.
Strictly speaking the intersection of the axes of the described cylinders and the unit sphere form not one circle but two, on opposing sides of the sphere. It is advantageous, however, to only use one of these circles, the one closest to the point, as this gives a directionality to the result.
Pseudocode for the algorithm is shown in Algorithm 4 . The actual Matlab implementa- tion is entirely vectorized through the use of multidimensional arrays and executes much faster than what would the naive but more readable description in Algorithm 4 using for-loops. The Matlab implementation can be run for several radii simultaneously and in that case reuses some of the intermediate variables for increased efficiency.
The grid points are all unit vectors and the voting result gives likely directions of cylin-
ders through or close to the transform center. When applying this transform locally to
Figure 3.5: Four of the possible cylinders with axes constrained to the sphere’s center and surfaces constrained to some specific point (cyan). The intersection of the sphere and all such cylinders is a circle (red). Voted for grid points are shown in red, others
in grey.
point clouds of trees it is important that the size of the region of the point cloud supplied to the transform is big enough for relevant branches to be represented properly but not so big that it includes points from branches not immediately connected to the current region of interest.
3.3.3 Test results
In this section the properties of the transform are illustrated by results on artificial test data. Figure 3.6 shows the result on an ideal point cloud where all points are on the surface of a cylinder of the exact same radius as the transform and with the cylinder axes passing exactly through the transform center. As evident in the figure, the grid point closest to the cylinder axis has received the most votes. The same result can be seen from another angle in Figure 3.7 (F), together with results on the same point cloud with different transform radii and transform centers. This figure serves to illustrate that with an off-axis transform center and a mismatched transform radius the result is slightly biased and imprecise, but still useful for determining the general direction of the cylinder.
Figure 3.8 illustrates a situation in which the result is somewhat poor. The point cloud in this case is a section of a point cloud from a simulated scan of an artificial tree. The section is located where the vertical stem of the artificial tree splits into two sub-stems.
Gaussian noise with a standard deviation of one tenth of the stem radius has been added
to the data. Local maxima of grid point votes are located close to the axes of the stem
Algorithm 4 Local Hough transform for cylinders
1:
procedure cylinderHough(c, r, p)
2:
for i ← 1, ||p|| do . For every point in p
3:
p
i← p
i− c . Translate point by subtracting transform center
4:
d ← |p
i| . Calculate point’s distance from origin
5:
if (d < r) then . If point too close to transform center
6:
continue . Skip this point
7:
end if
8:
r
p← r/d . Projected cylinder radius
9:
θ ← arcsin r
p. Angle of vote circle
10:
α
l← max(0, θ − tol) . Lower angle bound with tolerance
11:
α
u← min(π/2, θ + tol) . Upper angle bound with tolerance
12:
for j ← 1, ||gridP oints|| do . For every grid point
13:
α ← ∠(p
i, gridP oints
j) . Angle between point i and grid point j
14:
if (α > α
l) ∧ (α < α
u) then . If angle between bounds
15:
votes
j← votes
j+ 1 . Increment votes of grid point
16:
end if
17:
end for
18:
end for
19:
end procedure
and the left sub-stem. As the transform center is rather far from the axis of the right
sub-stem the corresponding local maximum is similarly biased to the other side of the
axis. Several grid points around that local maximum also have similar amounts of votes,
making the local maximum less distinct.
Figure 3.6: Result on an ideal point cloud. The grid points are sized and colored according to their number of votes. See color bar for comparison. The grey sphere
shows the radius used in the transform.
(a) r = 0.7r0, x = r0 (b) r = 0.7r0, x = 0.5r0 (c) r = 0.7r0, x = 0
(d) r = r0, x = r0 (e) r = r0, x = 0.5r0 (f) r = r0, x = 0
(g) r = 1.3r0, x = r0 (h) r = 1.3r0, x = 0.5r0 (i) r = 1.3r0, x = 0 Figure 3.7: Results on the same point cloud with different transform radii and trans- form centers. Grid points are sized and colored according to their number of votes. The rows show different transform radii and the columns show different transform centers.
The cylinder radius is denoted by r0. The transform radius is denoted by r. The radial displacement of the transform center from the cylinder axis is denoted by x. The grey spheres and red crosses illustrate transform radii and transform centers respectively.
Figure 3.8: The point cloud in this case is a spherical selection from an artificial point cloud of a branching point scanned from only one side. The point cloud is of a vertical stem splitting in two slightly below the middle of the figure. Grid points are sized and colored according to their number of votes and grid points with few votes have been
removed for clarity.
