Multichannel segmentation of planar point clouds using evolving curves

The paper deals with the multichannel curve segmentation of planar point clouds, scanned by a terrestrial laser scanner (TLS), which usually represent walls of buildings. Information derived from the point cloud can be used for quality check of a building (e.g. wall flatness or mutual perpendicularity of walls) or creation of virtual models of buildings. Therefore, segmentation of the planar point cloud into subsets representing the surface of the wall and other structures (e.g. door, electrical outlets) is desired. We describe a mathematical model of evolving planar curves, which is used for the segmentation of the point cloud. The evolution of curves is controlled by the properties of the scanned walls, such as colour and intensity. We discretize the model using the finite-volume method and the semi-implicit Inflow-Implicit/Outflow-Explicit (IIOE) scheme. We demonstrate the functionality of segmentation on real data.


Introduction
One of the most effective methods for collecting spatial data is terrestrial laser scanning (TLS).The resulting data are provided as a 3D point cloud, possibly including radiometric information, namely colour and intensity.In general, point cloud data can serve as a geometric basis for modelling objects and phenomena.Point clouds are used for visualisation, measurement, and many other surveying tasks such as 3D model creation or surface reconstruction (Previtali et al. 2014), digitalisation of cultural heritage (Adamopoulos and Rinaudo 2021;Haličková and Mikula 2016;Remondino and Rizzi 2010), deformation analysis of structures (Holst et al. 2015), flood modelling or landslide surveying (Giussani and Scaioni 2012).
In the present paper, we work with point clouds resulting from TLS scans of buildings.Nowadays, virtual models are created for managing the whole lifecycle of a building, from the design, construction, operation to facility management and demolition using Building Information Modelling (BIM).A BIM model is an object-oriented virtual model defined as a combination of graphic and non-graphic data.The graphic data in most cases have the form of a 3D model.The use of information from the 3D point cloud and the information derived from the BIM model can be used for the monitoring and quality check of a given structure during the construction process or after its completion (Bariczová et al. 2021).
Our point cloud data include 3D coordinates of the points, colour channels captured by a digital camera, and the intensity channel containing information related to the strength of the reflected laser signal.Our goal is to find planar subsets of the original point cloud containing points that represent the surface of walls.These subsets could be used, e.g. for automation of the creation of as-built models of buildings.
In our typical workflow, we perform a two-phase segmentation.Primary segmentation is performed by plane fitting.Planar subsets of the entire point cloud are selected.Each subset contains points lying within a specified distance from the fitting plane and will be referred to as planar point cloud.The planar point cloud corresponding to a wall often contains objects that do not represent the surface of the wall (e.g.door and electrical outlets) and may affect the result of the quality check (Bariczová et al. 2021).The role of the secondary segmentation based on evolving curves is to segment each planar point cloud into subsets representing the wall surface and other objects, e.g.door, whiteboard, and the remaining unsegmented points (waste).
For the primary segmentation, many methods have been developed, see the review by Grilli et al. (2017).We usually use a modification of the Random Sample Consensus (RANSAC) algorithm (Fischler and Bolles 1981) called M-estimated Sample Consensus (MSAC) (Torr and Zisserman 2000), possibly followed by some filters, e.g. a normal filter which discards points whose local normals deviate from the normal of the fitting plane by more than a given threshold angle.In addition to the MSAC method, we also use the approach developed in Honti et al. (2022).
For the purpose of this paper, we consider a planar point cloud as the input data for the secondary segmentation, which will be described in detail.As we will see, the basic idea of our approach is to create a set of images representing the point cloud data1 and subsequently segment these images.In the field of image processing, many segmentation methods based on curve and surface evolution models have been developed.We distinguish two main approaches to handle evolution problems, the so-called direct (or Lagrangian) approach (Dziuk 1999;Mikula and Ševčovič 2004;Balažovjech and Mikula 2011) and level set (or Eulerian) approach (Sethian 1999;Osher et al. 2004;Caselles et al. 1997;Zhao et al. 2000;Corsaro et al. 2006).Our secondary segmentation is based on a Lagrangian curve evolution model, which is computationally fast because it solves only 1D curve evolution problem.The Eulerian approach treats the segmentation curve as the zero-level set of the level set function and the whole function is evolving.Therefore, one has to solve 2D evolution problem which is computationally more expensive.
The evolution of segmentation curves is driven by information from multiple channels of the point cloud data, e.g.colour, intensity, or distance from the fitting plane.The quality of the discretisation mesh is crucial for stability and accuracy in the Lagrangian approach.Various techniques for tangential redistribution of points improving the mesh quality have been developed (Hou et al. 1994;Kimura 1997; Ševčovič and Mikula 2001;Barrett et al. 2011;Ševčovič and Yazaki 2011;Benninghoff and Garcke 2014), and we adopt asymptotically uniform redistribution of points from (Mikula and Ševčovič 2004).Moreover, an algorithm for the detection and treatment of the topological changes is needed (Benninghoff and Garcke 2014;Frei and Garcke 2016).We use an efficient O(n) algorithm developed in (Mikula and Urbán 2012;Balažovjech et al. 2012;Ambroz et al. 2019).
The mathematical model for curve evolution is an extension of the model from the paper (Mikula et al. 2021a) which deals with one channel segmentation of forests from satellite images.Contributions of the present paper are mainly the fully automatic methodology for segmentation of the point clouds (meaning two segmentation phases, creation of the representative images, curve segmentation and creation of the final segments of the point cloud), the multichannel approach using a novel construction of the aggregated edge detector, and simultaneous segmentation of multiple regions (requiring a different detection of merging of curves).Minor contributions include construction of the initialisation mask, a simpler approximation of the signed curvature (without using arccos as in Mikula et al. (2021a)) and a local weight for the curvature regularisation.To make the text self-content, we included some ideas from Mikula et al. (2021a) in this paper.A preliminary version of the curve segmentation of planar point clouds is outlined in the paper (Bariczová et al. 2021) which deals with quality check of walls in buildings.In the paper (Bariczová et al. 2021), only a brief intuitive description of the model is presented, and the functionality of the segmentation is shown only on one example.In this paper, we present the final version of the methodology with the full mathematical description of the model and all computational details of our approach.The functionality of the model is demonstrated on multiple representative examples.
The text of this paper is organised into three main sections.In Sect.2, we introduce the main ideas of our segmentation methodology, mainly the creation of the representative planar point cloud images (Sect.2.1), mathematical model (Sect.2.2), construction of the edge detectors (Sect.2.3), and our approach to the multichannel segmentation (Sect.2.7).In Sect.3, we present the numerical discretisation of the mathematical model.The Sect. 4 shows several numerical experiments with real data that illustrate the functionality of our approach.

Point cloud segmentation methodology
In this section, we describe our method for segmentation of the point clouds.We obtain our input point cloud data using the Trimble TX5 scanner.The point clouds usually contain millions of points with the average resolution about 2 − 6 mm.A typical example is a flat 1.Primary segmentation by plane fitting.Selection of planar subsets of the point cloud using the MSAC method, Fig. 1.Each subset consists of points located within a specified threshold distance from the fitting plane.The fitting can be possibly followed by a normal filter removing points whose local normals deviate from the normal of the fitting plane by more than a specified threshold angle.The resulting planar point clouds 3 serve as input data for secondary segmentation.2. Secondary segmentation by curve evolution.Segmentation of each planar subset of the point cloud using evolving curves (Fig. 2).
2 It might seem like a continuous surface, but it is an unorganised cloud of points without any structure.
Figure 1 shows a very dense point cloud with points visualised as tiny square planar patches.
3 For ease of expression, we call these subsets planar point clouds, although in reality they are only approximately planar.
(a) Preprocessing: Creation of representative images corresponding to the channels (such as colour or intensity) of the planar point cloud.(b) Evolution: Segmentation by evolving curves using representative images.(c) Postprocessing: Extraction of points in the point cloud corresponding to the segmented regions.
In the following subsections, we describe the steps of the secondary segmentation in detail.As we mentioned, the same procedure is performed for each planar point cloud.Therefore, in the following, we will work with a single planar point cloud.In Sects.2.2-2.6,we present the method for one-channel segmentation and in Sect.2.7, we propose a technique to use information from multiple channels.

Preprocessing: creation of the representative images from planar point cloud
To reduce the dimension of the problem, we transform (rotate and move) our planar point cloud to the new coordinate system in which the first two dimensions span the fitting plane.If we have a transformed point with coordinates x = (x 1 , x 2 , x 3 ), the coordinates x 1 and x 2 describe the position of the point in the plane and the x 3 axis is orthogonal to the fitting plane, i.e. the absolute value of x 3 is the distance from the fitting plane, denoted as D = |x 3 |.
Subsequently, we create a regular square mesh in the x 1 , x 2 plane.We denote the pixel size (edge length) as h and usually set it as a multiple of the resolution of the point cloud (e.g.3×resolution).Then, we represent the point cloud data by a set of bitmap images, one image for each available channel.Usually, we use the colour channels R, G, B, the intensity channel I, and the distance channel D. For example, we take the R channel of colour and compute the value in each pixel (square of the mesh) as the mean value of the R channel of the points which lie in the pixel, see Fig. 3. Finally, we rescale the image values to the range [0, 1].The domain of the images will be denoted by ⊂ R 2 .
The resulting images are pixel representations of the input planar point cloud and will be segmented using evolving curves.

Mathematical model of the curve evolution
In this section, we present a Lagrangian curve evolution model for segmentation of the planar point clouds obtained by TLS.The input data is an image I, which is one of the channels.The use of multiple channels will be described in the Sect.2.7.The mathematical model for curve evolution is an extension of the model from the paper (Mikula et al. 2021a) and we employ a similar notation.The segmentation curve is an evolving closed planar curve with a position vector of its points denoted by x(u, t), where u ∈ [0, 1] is a parameter that goes along the curve and t ∈ [0, t f ] is the time (with t f being a final time).The fact that the curve is closed means x(0, t) = x(1, t).
The evolution is driven by a suitably designed velocity field v(x, t), therefore, the basic evolution model is where ∂ t x := ∂x ∂t denotes the time derivative of the position vector, that is, the velocity of point x.Equation ( 1) is coupled with an initial condition x(u, 0) = x 0 (u), u ∈ [0, 1].The initial curve x 0 is a small circle (Fig. 4) automatically placed inside the segmented region (e.g. the small white circle in Fig. 2b).Details on automatic insertion of the initial curve are provided in Sect.2.8.We consider the velocity v(u, t) in the following form (see also Mikula et al. (2021a)): where N denotes the positively oriented unit normal vector defined by rotation of the unit tangential vector T = ∂x ∂s (also denoted as ) by π 2 in the positive direction, with s = u 0 ∂ u x du denoting the arc length parameter of the curve.Therefore, the normal can be expressed as 2) is the signed curvature and ∇ is the gradient operator.Functions E and B will be properly defined later.Now, we give a brief description of the terms in (2).The role of the first term B(u, t) N(u, t) is to expand the segmentation curve in the normal direction from its initial shape through the segmented region towards its border, and the 'blowing' function B(u, t) controls the speed of expansion.
The second term attracts the points of the curve towards the borders of the segmented region.Information about the borders is contained in the edge detector function E(x) ∈ [0, 1], x ∈ which has values close to 0 near the edges and close to 1 in homogeneous regions of the planar point cloud (see Fig. 5d).The negative gradient of E points towards the lower values and, therefore, is suitable to attract the segmentation curve towards the edges.
The time-dependent parameter λ(t) ∈ [0, 1] serves as a weight between the expansion and the edge attraction term.The standard approach is to set λ(t) at the beginning (t = 0) to a number 0 < λ 0 1 (i.e. the edge attraction does not dominate), keep it unchanged until the curve is close to the border (moving very slowly) and then switch λ(t) to 1, which turns Fig. 4 An evolving closed planar curve plotted at times t = 0, t 1 , t 2 , where 0 < t 1 < t 2 off the expansion term.The setting λ(t) = λ 0 will be called the expansion phase and the setting λ(t) = 1 we call the attraction phase.
The last term in (2) is called curvature regularisation and has a smoothing effect.We use it to deal with the noise and to smooth the sharp edges of the segmentation curve during the expansion phase.The importance of the curvature regularisation was studied, e.g. in Mikula et al. (2021b).The parameter δ(u, t) weights the influence of the term and will be properly defined in Sect.2.7.
The first two terms in (2) could be also called together as 'external velocity', because they substantially depend on external quantities which are properties of the segmented image, unlike the curvature regularisation term, which depends on the signed curvature k and normal N, which are properties of the curve.However, we may set the weight δ to be dependent also on the speed of the curve evolution, as we shall see later.

Edge detector
The (smoothed) edge detector function E(x) is constructed as follows (Mikula et al. 2021a).First, we prefilter the original image I, for example, using the Gaussian filter to get I σ 0 = G σ 0 * I (the convolution with the Gaussian kernel G σ 0 with variance σ 0 ) and then compute the norm of the gradient of the filtered image ∇I σ 0 (x) .If the original image is very noisy and strong filtration is needed, we can use the anisotropic diffusion (Perona-Malik) filter (Perona and Malik 1990) instead.However, in some situations, the prefiltration is not necessary because some filtration is actually done in the preprocessing step (Sect.2.1) when computing pixel representations of the point cloud data, which is akin to a mean filter.In such cases, one may directly compute the norm of the gradient of the original image ∇I(x) .
Next, we find edges in the image I by computing an edge detector function (3) The gradients are large near the edges, and therefore, g(x) is close to zero in such regions.
In homogeneous areas, the gradients are small, and thus g(x) is close to 1.The parameter μ > 0 controls the sensitivity of edge detection.If it is small, the value of g(x) is close to zero (detects an edge) only for very sharp edges (i.e.very large gradients).If the parameter μ is large, the function g(x) has values close to zero even for smaller gradients, which means that g(x) is more sensitive to edges.The reasonable values of μ are from the interval (0, 10].Finally, we apply the Gaussian filter (convolve the function g(x) with the Gaussian kernel with variance σ ) to obtain the smoothed edge detector function (4) Smoothing makes the edges 'wider' and causes the edge attraction term −∇ E(x) in (2) to point in the right direction (towards the edge) in a larger neighbourhood of the actual edge, which helps the edge attraction.The edges become too wide and blurry for large σ , so small values of σ are reasonable, we use σ = 0.5.In Fig. 5a, we can see a point cloud of a wall with two pictures.Image I for the red channel is in Fig. 5c and the corresponding smoothed edge detector function E(x) is plotted in Fig. 5d.We note that if image I contains holes (as in Fig. 7b), i.e. not a number (NaN) values corresponding to empty pixels (containing no points of the point cloud), the gradients and filtration need to be computed carefully.For example, if we carelessly calculate the derivative ∂ x I in the (i, j)-th pixel by the finite difference (∂ x I) i j ≈ (I i+1, j −I i−1, j )/2 h and the value  I i−1, j or I i+1, j would be NaN, the result would be NaN.Therefore, we first interpolate (or extrapolate) empty pixels, then compute the finite difference, and eventually reset the NaN values to preserve the information about the holes (if I i j was NaN, then reset (∂ x I) i j to NaN).If a point x(u, t) of the curve lies in an empty (NaN) pixel, we set both the expansion and edge attraction term in (2) to zero.

Blowing function
The blowing function B used in the expansion term in ( 2) is defined as (see also Mikula et al. (2021a)) where H is called the homogeneity function and is defined as where d 0 is a threshold value and d(u, t) = I(x(u, t)) − I t is the difference function which computes the (absolute value of the) difference between the value I(x(u, t)) and the average value I dx in the region V (t), which denotes the set of points x ∈ visited by the curve since the beginning of the evolution, together with points inside the initial curve (i.e.V (0) is the interior of the initial curve).If the value of the image I at a point x(u, t) of the curve is similar to the points x ∈ already segmented by the curve, then the difference d(u, t) is small and the point x(u, t) should move, which is provided by homogeneity H (u, t) = 1 for small d(u, t).On the contrary, if the difference d(u, t) greater than threshold d 0 , the definitions ( 5), ( 6) give zero expansion of the curve at the point x(u, t).
A natural definition of the blowing function could simply be B = H .However, then the expanding term could sometimes drive the points of the curve across the edges.This is the reason why we multiply the homogeneity function by the edge detector function E (which is close to zero in the neighbourhood of edges).

Normal and tangential speed
The basic evolution model (1) moves the points on the curve both in the normal and tangential directions.From a continuous (analytical) point of view, only the normal component β = v•N of the velocity (2) has an effect on the shape of the curve.Therefore, we could use the evolution model ∂ t x = βN.However, in the discrete setting, it is convenient to enrich this model with a suitable tangential speed α which can be designed to control the distribution of the discretisation points along the curve.This becomes crucial in numerical computations because inappropriate positions of the points can lead to unacceptable errors.Therefore, instead of the basic model ( 1), we use the model where β = v • N is the normal speed calculated from (2) as where is the so-called external normal speed.
The tangential speed α is constructed to achieve a uniform distribution of the points in the discrete setting.We follow the approach from Mikula and Ševčovič (2004).The length of the curve segment between the points x( û, t) and x( û + u, t) is û+ u û g du, where g := ∂ u x .We want this length to be asymptotically (for t → ∞) same for every û (each pair of points), since this corresponds to a uniform distribution of grid points on the discretised curve.However, the total length of the curve L = 1 0 g du changes during evolution, so we have to study the relative length of the segment (with respect to the total length).For an asymptotically uniform redistribution of points, the relative length of the segment approaches a constant This is satisfied if g L → c, where c > 0 is a constant.Since 1 0 g L du = 1, we have to set c = 1 and the condition is lim t→∞ g L = 1.This is fulfilled if the ratio g L obeys a relaxation equation where the relaxation parameter ω > 0 controls the speed of relaxation.On the left-hand side, we need a formula for the derivative of g and L. The first one is where we used: 1. Derivative of the norm Integrating (10), we obtain the equation for the derivative of the curve length where kβ := 1 L L 0 kβ ds denotes the average value of kβ over the curve.In the above calculation we used ds = ∂ u s du = g du and The application of the quotient rule in ( 9) and the use of formulas ( 10), ( 11) give After multiplication by L g and rearrangement, we obtain the equation for the tangential speed α To ensure a unique solution, we can simply fix the value of α in one point (e.g.α(0, t) = 0).However, this can lead to a larger tangential motion than is needed for uniform redistribution.Therefore, to minimise the tangential motion, we ensure uniqueness by requiring zero average tangential motion α = 1 L L 0 α ds = 0 for each t ∈ [0, t f ], see also Ambroz et al. (2019).

Curvature regularisation
The weight δ(u, t) of the curvature regularisation term in ( 2) and ( 8) can be simply a constant δ(u, t) = δ k (usually a small number, e.g.0.01).However, it is convenient to use a more sophisticated setting that considers the phase of the evolution and the relative speed of the point x(u, t) where the local weight δ loc (u, t) is defined as follows (Fig. 6): (u, t) .Thus, in the expansion phase, the proposed weight δ(u, t) is a function of the relative external normal speed w(u, t) of the point x(u, t).For slow points, we apply a small regularisation4 weight, and for fast points, we set a high regularisation. 5In the attraction phase, we turn off curvature regularisation, which enables sharp corners of the final segmentation curve.

Multichannel segmentation
So far, we have only discussed single-channel segmentation.Now, we propose a technique to use information from multiple channels and demonstrate it in a simple but illustrative example in Fig. 7.There is a white wall with three plastic plates: blue, red, and white.There are holes (NaN values, black) in Fig. 7b in parts of the planar point cloud in which there are no point cloud points within the threshold distance from the fitting plane.Let us denote the images corresponding to the channels R, G, B, I and D as I l , l = 1, . . ., 5, respectively, see the left column in Fig. 8. First, we will discuss the construction of the smoothed edge detector function E and then the computation of the blowing function B. In the right column of Fig. 8, we can see the edge detector functions g l (x) corresponding to the images I l computed using (3).In practice, we can set different sensitivity parameter μ l for each channel.Notice that some edges are discovered in one edge detector g l (x) but not in another, e.g. the red plate (with (R,G,B) ≈ (1, 0, 0)) is 'invisible' for the R channel (because the wall is white, i.e. (R,G,B) ≈ (1, 1, 1), see first row of Fig. 8), but it is clearly distinguished by the B channel (and vice versa for the blue panel).The intensity edge detector (fourth row in Fig. 8) revealed multiple edges that are not seen in the colour (R, G, B) channels6 : the edge of the white plate, electrical outlets, light switch, and electrical box covers.These objects have different reflexivity than the wall, and this is captured by the intensity channel which is related to the strength of the reflected laser signal.The electrical outlets, light switch and box covers were also clearly detected by the D channel, see the last row of Fig. 8.
A natural idea is to use the data from all available channels R, G, B, I, D (sometimes the colour data might be missing) to construct one edge detector.Therefore, we aggregate the edge detector functions g l (x), l = 1, . . ., 5 into a single edge detector function which contains aggregated information about all edges in all channels, see Fig. 9. Indeed, if at a point x any (even one) of the edge detectors has a value g l (x) close to zero, then the value g(x) will be close to zero.Moreover, if multiple channels have values close to zero at a point x, then the resulting value g(x) will be much closer to zero.This approach is also convenient, in the situation if multiple edge detectors agree on some mild edge, the resulting aggregated edge will be more distinct.Finally, we smooth the edge detector function ( 14) using ( 4) to obtain the final smoothed aggregated edge detector function E(x) and use it in the curve segmentation.
Since the range of the edge detector function ( 3) is (0, 1], the aggregation ( 14) can be interpreted as fuzzy logic aggregation using the product aggregation operator A : van Krieken et al. (2022).Fuzzy logic in general is a many-valued logic in which the truth values of variables x i are real numbers from the interval [0, 1], i.e. there is a concept of partial truth while x i = 0 represents completely false and x i = 1 is completely true.The product is a generalisation of the AND operator from Boolean logic: the conjunction x 1 ∧ x 2 with x 1 , x 2 ∈ {0, 1} generalises to Now, we shall briefly discuss the aggregated blowing function B. We compute it using (5) as in the single-channel case.However, now E is the smoothed aggregated edge detector function described above and H is the aggregated homogeneity function computed as where H l (u, t) ∈ {0, 1} denotes the homogeneity function computed using (6) for the image I l .Equation ( 15) guarantees that the value B(u, t) of the blowing function is nonzero only if the difference d l (u, t) = I l (x(u, t)) − I l t is smaller than the threshold d l 0 for each channel.If d l (u, t) ≥ d l 0 for some channel,7 then the blowing function (expansion term) vanishes at the point x(u, t) of the curve.The threshold value d l 0 can be set differently for each channel.However, in this paper, we will use d l 0 = 0.5 for all channels.

Inserting initial curves
At the beginning of the segmentation, we insert a small circle (or multiple circles8 ) into the image.We use the radius r 0 = 3h, with h denoting the pixel size (edge length).In order to find suitable candidates for positions of the initial circles, we construct an initialisation mask M, Fig. 10.First, we filter the smoothed edge detector Gaussian filter with a higher standard deviation (e.g.σ M = 2) and a larger kernel size (e.g. 7 × 7).Then, we threshold the function E σ M as follows: The region where M(x) = 1 will be referred to as segmentable region S ⊂ .We note that here we perform a careless filtration, which intentionally spreads NaN values to their neighbourhood (differently from the Sect.2.3).The centres of the initial circles are randomly placed in the segmentable region S (more precisely to random pixel centres in S ).The construction of mask M prevents the insertion of circles on the edges, holes (NaN values), and their neighbourhood.If all curves stop evolving, we remove the segmented region V (t) of each curve from the segmentable region by setting M(x) = 0 for all x ∈ V (t).Subsequently, we randomly add a new circle (or multiple circles) into the segmentable region.If we insert multiple circles (either at the beginning or during the process of segmentation), it is convenient to remove the interior and possibly some neighbourhood of a circle from the segmentable region after each curve insertion.This prevents the generation of intersecting or too close curves.When a segmentation target is reached, e.g.99 % of the initial segmentable area is segmented, we end the segmentation process.

Postprocessing: creation of the point cloud segments
When the curve segmentation is finished, the planar point cloud is divided into segments.We take the segmented region V (t f ) of each segmentation curve and select all points in the point cloud that lie in V (t f ).The last segment is a 'waste after segmentation' and corresponds to regions of which were not visited by any curve.

Numerical scheme
In this section, we discuss the numerical discretisation of our method in great detail.For the purpose of the derivation of the numerical scheme, let us rewrite the model ( 7), (8) using T = ∂ s x, N = ∂ s x ⊥ and the Frenet-Serret formula kN = ∂ s T = ∂ 2 s x.We obtain the following equation: with is actually a system of two PDEs, since x = (x 1 , x 2 ) is the unknown position vector.Equation ( 17) has a form of so-called intrinsic PDE.The tangential term −α∂ s x represents an intrinsic advection along the curve, with −α being the speed of the advection, and the curvature term δ∂ 2 s x can be regarded as an intrinsic diffusion with diffusion coefficient δ.Now, we will discretize equation ( 17).First, we perform the space discretisation which is based on the finite-volume method (Mikula et al. 2021a).Then, we present the time discretisation both by the standard semi-implicit scheme (Ševčovič and Mikula 2001) and semi-implicit scheme with Inflow-Implicit/Outflow-Explicit (IIOE) technique to treat the advection term (Mikula and Ohlberger 2011;Balažovjech et al. 2012;Mikula et al. 2014).

Space discretisation
The mesh will consist of m grid points x i , i = 1, . . ., m, where x i (t) ≈ x(u i , t), u i = (i −1) u and u = 1 m .The construction of the finite volume V i corresponding to the point x i consists of two line segments connecting the points , where h i := x i+1 −x i is the length of the segment {x i , x i+1 }.We integrate equation ( 17) over the finite volume V i and obtain Fig. 11 The notation on the curve mesh where we approximate some quantities by constants on V i , namely ∂ t x ≈ ∂ t x i , α ≈ α i , δ ≈ δ i and β e ≈ β e,i .After using the Newton-Leibniz formula, we have the following: which evaluates to Finally, we approximate the derivative ∂ s x in the midpoints , respectively.The space discretisation of the PDE (17) reads for each i = 1, . . ., m, with x 0 := x m and x m+1 := x 1 .

Time discretisation by semi-implicit scheme
In this paper, we show two different time discretisations.In this section, we perform the time discretisation using a semi-implicit scheme (Ševčovič and Mikula 2001) and motivate the need for a more sophisticated semi-implicit IIOE scheme (Mikula and Ohlberger 2011;Balažovjech et al. 2012;Mikula et al. 2014), which is presented in the Sect.3.3 and used in our implementation.The motivation is given from a different point of view from those of (Mikula and Ohlberger 2011;Mikula et al. 2014;Balažovjech et al. 2012).We discretize the time interval [0, t f ] using discrete times t n , n = 0, . . ., n max with t 0 = 0 and t n max = t f .Let us denote the length of the time step by τ n = t n+1 −t n .The discretisation is not necessarily uniform in time, and we usually choose a smaller τ n in the attraction phase.
We approximate the time derivative in (18) by the forward finite difference ∂ t x i ≈ Then, we take the coefficients h i , α i , δ i , β e,i from the n-th (old) time step, the advection and diffusion terms are approximated implicitly, that is, we take x i from (n + 1)-th (new) time step, and the expansion term explicitly, meaning that x i is taken from the old time step.We obtain the following: which can be rearranged as for n = 0, . . ., n max − 1.We denoted The formula ( 19) is a system of m vector equations (2m scalar equations) with unknowns It is a cyclic tridiagonal system. 9Such systems can be solved by a modified Thomas algorithm (also known as tridiagonal matrix algorithm) (Press et al. 2007). 10If the system matrix is diagonally dominant, the Thomas algorithm is stable.The diagonal dominance in our case reads where we cancelled the absolute value on the LHS because all quantities are nonnegative.However, we cannot cancel it in the terms on the RHS, because α n i ∈ R. Therefore, the condition ( 22) is not always satisfied.We can achieve diagonal dominance with the appropriate choice of the time step length τ n .From ( 22), we have the following: which must hold for each i = 1, . . ., m, see also Balažovjech and Mikula (2011).Therefore, we can evaluate the RHS of ( 23) for each i and find the minimum that will give us the suitable time step length τ n .Such computation of τ n slows down the evolution algorithm and, moreover, if we obtain a very small τ n , the evolution algorithm would be even slower.However, there is a trick (Mikula et al. 2014) which will guarantee stability without any restriction of the time step length.Now, we give a motivation for the trick.Let us consider the case with no tangential speed, α n i = 0.Then, we can cancel the absolute values also on the RHS of ( 22) and obtain the condition ≥ 0 which holds (even strictly) for any τ n , thus for zero tangential speed the system matrix is strictly diagonally dominant.In the case of non-vanishing tangential speed, the idea will be the same, i.e. to obtain the system matrix for which we can cancel all absolute values in diagonal dominance condition and then cancel out all terms except 9 Since the subsystems for x 1 and x 2 coordinates are independent of each other, we can solve the subsystems independently.If we would not treat the expansion term explicitly, we would end up with a more complicated system with 2m by 2m matrix with 5 nonzero elements in each row. 10The presence of the coefficients A n 1 and C n m does not allow to use the Thomas algorithm directly.However, one can use a trick from Press et al. (2007) to obtain the standard tridiagonal matrix and then use the Thomas algorithm.

Time discretisation by semi-implicit IIOE scheme
In this section, we present the semi-implicit Inflow-Implicit/Outflow-Explicit (IIOE) scheme (Mikula and Ohlberger 2011;Balažovjech et al. 2012;Mikula et al. 2014).The goal is to obtain a diagonally dominant matrix with no restriction on the time step length also in the case with nonzero tangential speed.We have seen that the tangential term causes problems.We will take a term containing α n i implicitly (i.e. it will occur in the system matrix) only if it has a suitable sign; otherwise, we will take it explicitly (move it to the RHS).By a 'suitable sign', we mean the sign that will allow us to cancel all absolute values in the diagonal dominance condition.
First, we rearrange the tangential term in (18) as follows: where we introduced the notations To clarify the notation: if the advection speed −α i is positive, we have inflow into the finite volume V i at the point x i− 1 2 and outflow at x i+ 1 2 , if the advection speed −α i is negative, there is outflow at x i− 1 2 and inflow at x i+ 1 2 .As we shall see, if the sign of α n i corresponds to the inflow (either at the point ), the corresponding terms have a suitable sign and allow cancellation of absolute values in the diagonal dominance criterion.Therefore, to perform the time discretisation of the tangential term (24), we take the Inflow terms Implicitly and Outflow terms Explicitly The final system has the form (19), i.e.
for n = 0, . . ., n max − 1.The resulting coefficients A n i , B n i , C n i , and D n i are The diagonal dominance criterion All terms are nonnegative, and therefore, we can cancel all absolute values and the condition reduces to ≥ 0. This holds for any τ n and the solution of the system (26) by the modified Thomas algorithm is guaranteed.
Let us briefly discuss the stopping criterion for evolution.If the segmentation is finished and the curve does not evolve, we should break the time loop even before the final n max -th time step is computed by ( 26).If a large percentage (e.g.99 %) of the points do not move, we switch the evolution phase from expansion to attraction (setting λ n = 1) and perform 20 evolution steps in this phase.By not moving points, we naturally mean points which move very slowly (e.g.x n+1 i − x n i < 0.05h) but also the points jumping back and forth between two neighbouring pixels, which happens very often at an edge.

Discretisation of the normal and tangential speed
The normal speed is given by ( 8) and we approximate it at the point x n i as follows: The evaluation of δ n i by ( 13) is straightforward.The discretisation of the signed curvature k n i will be discussed in Sect.3.5.The blowing function is calculated as , with the difference d n i in the calculation of H n i (using ( 6)) approximated as p∈V n I p where V n is the set of pixels p that approximates the region V (t n ), |V n | denotes the number of pixels in V n and I p denotes the value of image I in the pixel p.The edge detector function E is computed using (3), (4) in the bitmap image I and the gradients ∇I and ∇ E are approximated by central differences.The positively oriented normal in (28) is calculated as . The tangential speed is the solution of the equation ( 12) and we discretize it as follows.
The derivative on the LHS is approximated by the finite difference, ∂ s α ≈ . We regard it as a central difference at the point x n i+ 1 2 and, therefore, approximate the RHS at that point.The average (integral) term is approximated by the Riemannian sum, ).The approximation of the signed curvature k n i+ 1 2 will be explained in the next section.The norm g = ∂ u x is approximated by Putting all together, we have for i = 1, . . ., m, which is a two-diagonal system for unknown α n i , i = 1, . . ., m.The system can be easily solved by setting α n 1 = 0 and successively computing α n 2 , . . ., α n m by Finally, we redefine the tangential speeds by subtracting their average to ensure zero average tangential motion.

Discretisation of the signed curvature
In the discretisation of the normal and tangential speeds, we need a formula for approximate curvature at the grid points x n i and at the midpoints x n i+ 1 2 of the edges.We will discretize the formula It can be derived directly from the definition.Using the notation g = ∂ u x , we have where we used: 1. Definition of the unit tangent vector T = ∂ u x g and the chain rule ∂ s = 1 g ∂ u (see 2. under (10)).

Quotient rule and ∂
To approximate the signed curvature at the grid point x n i , we can use the central differences (see Fig. 11) which leads to the discretisation At the midpoint x n i+ 1 2 of the edge, we use the central differences (see Fig. 11) and the normal , which gives

Topological changes
When dealing with curve evolution using the Lagrangian approach, one has to treat topological changes that occur during the evolution.The evolving curve can split into multiple curves, or several initial evolving curves can merge.The detection and treatment of the topological changes can be highly time-consuming because the natural approach has computational complexity O(m 2 ), where m denotes the number of grid points (it includes computation of pairwise distances between all grid points).However, we adopt the method from Mikula and Urbán (2012); Balažovjech et al. (2012); Ambroz et al. (2019), which has complexity only O(m).This makes the Lagrangian method for curve evolution efficient and applicable in complex situations.The basic idea of Mikula and Urbán (2012); Balažovjech et al. (2012); Ambroz et al. ( 2019) is to construct a background mesh (we use the mesh described in Sect.2.1).To detect splitting, we walk through all curve grid points and label the underlying cells by the grid point numbers.If the current i-th point belongs to a cell already labelled by the j-th point and the difference i − j is large enough (e.g.i − j > 3), the splitting is detected.Since we pass the curve only once, the computational complexity is O(m).The merging detection is similar, we walk through the points of all curves and label underlying cells by the curve number.If the current point belongs to the cell already labelled by a point of another curve, the merging is detected.The most complete explanation of the algorithms (including pseudocodes) is given in Ambroz et al. (2019).
For the reliable topological changes detection and treatment, it is required to maintain the mean curve segment length close to a specific constant h d called desired curve segment length (desired distance between neighbouring grid points).We use setting h d = h/2, where h is the pixel size.In our method, every curve is asymptotically uniformly discretised (see Sect. 2.5), i.e. all segment lengths are close to their mean value.However, the mean value can differ between different curves, and it increases for expanding curves and decreases for shrinking ones.Such difference is undesirable, especially in merging of curves, so we maintain the mean curve segment length close to h d for all curves by adding or removing points, after which the asymptotically uniform redistribution is quickly obtained due to the tangential velocity α.The detailed algorithm for adding and removing grid points is in Ambroz et al. (2019).We note that small curves with less than m = 15 grid points are deleted.
The merging detection in this paper is slightly different.If two evolving curves numbered c 1 and c 2 get close to each other, we have to decide if they are segmenting regions with similar properties (curves should merge) or different properties (curves should not merge).First, let us discuss the single-channel case.We compute the average values I t c 1 and I t c 2 of image I in the regions V c 1 (t), V c 2 (t) already segmented by the curves c 1 and c 2 , respectively, see definitions in the Sect.2.4.If the difference d c 1 ,c 2 (t) := | I t c 1 − I t c 2 | is smaller than the threshold value d 0 , merging is detected.
For the correct topological changes detection, we use the time step length τ n < h/β n e,max , where β n e,max := max i∈{1,...,m} β n e,i denotes the maximal external normal speed in the n-th time step.Small time step length is crucial because it prevents the points of the curve from skipping a pixel one time step (and missing possible topological change).
In the case of multichannel segmentation, we compute the difference d l c 1 ,c 2 (t) for each channel I l , l = 1, . . ., 5 and compare it with the corresponding threshold d l 0 .If d l c 1 ,c 2 (t) < d l 0 for all l, then the curves c 1 and c 2 segment the regions V c 1 (t), V c 2 (t) with similar properties (in all channels) and the merging is detected.If d l c 1 ,c 2 (t) ≥ d l 0 for some l, merging is forbidden.In the first experiment, we perform the curve segmentation of the planar point cloud corresponding to the white wall with three plastic plates depicted in Fig. 7.The point cloud consists of approximately 1.2 million points, its resolution is about 3 × 3 mm and the pixel size was set to h = 12 mm.The images and edge detectors for individual channels were discussed in the Sect.2.7 (see Fig. 8).The smoothed aggregated edge detector is visualised in Fig. 9. Several steps of segmentation are shown in Fig. 12.The curves coloured in magenta are stopped.We start the segmentation with a single curve, and when all curves stop, we  12.The curve meets the first edges around the time step n = 400 and the edge detector correctly prevents the curve from crossing the edge.In the next subfigures, we also colour already segmented region V n .In the time step n = 1700, we encounter a surprising observation: the curve stops on the white plastic plate, which is not visible in the RGB image, but thanks to the intensity channel (fourth row of Fig. 8), it is clearly captured in the aggregated edge detector, Fig. 9.The right part of the curve is enclosing the electrical box cover and the curve splitting is going to occur.In the time step n = 4000, we see that the first segmentation curve correctly segmented the wall and excluded unwanted objects: the plates, but also the electrical elements (e.g.electrical outlets in the lower part of the wall).Another curve (green) was inserted onto the blue plate.It correctly segmented the plate and did not merge with the previous curve (see n = 4400).After it stopped, a new (black) curve was inserted onto the red plate.The last subfigure of Fig. 12 shows the segmentation just before the termination.After the dark green curve on the left stops, the segmentation target is reached and no more curves are inserted.The final segments of the planar point cloud are visualised in Fig. 13.The magenta points represent the waste after segmentation, i.e. the regions that were not segmented by any of the curves.
In the second experiment, we segment a wall from the interior of P.O.Hviezdoslav Theatre in Bratislava, Fig. 14.The planar point cloud has roughly 0.4 million points with an average resolution of about 6 × 6 mm.Pixel size was set to h = 18 mm.We can see two pictures on the wall, a sink with facing in the left corner.The hole above the facing corresponds to a mirror (from which the laser beam did not reflect directly to the scanner).The representative images I l and the corresponding edge detectors g l are shown in Fig. 15.We see the edges of the pictures in all channels.Some edges of these pictures are not complete in each channel, however, in the smoothed aggregated edge detector, Fig. 16, we see distinct and complete edges (mainly due to the intensity channel).The edge around the facing is mild in all edge detectors, but since the edge is present in multiple edge detectors and the aggregation is done by multiplication ( 14), the aggregated edge is strong in Fig. 16, lower left part.An interesting situation is in the top left area of the wall strongly illuminated by the lamp.There is a mild edge in all colour channels, but not precisely in the same region.As a result, the aggregated edge in Fig. 16 is mild and diffuse.Several steps of the segmentation process are shown in Fig. 17.In this experiment, we start with four random initial curves.In the time steps n = 300 and n = 352, we see that three of them merged into one.However, as we see in the time step n = 600, the green and blue curves did not merge, which is correct.The evolution stopped in the time step n = 924.There are two main segments (the first corresponding to the wall and the second to the facing) and the waste after segmentation, see Fig. 18.The curve segmentation correctly excluded pictures, electrical outlets, and switch (above the facing).Interestingly, the blue curve did not reach the facing edge in the area to the left of the mirror.The reason is that the wall is crooked (not flat) in this region, which is slightly captured by the edge detector corresponding to the D channel (last row of Fig. 15).The resulting edge is mild, so this is not the reason for stopping the curve.The curve stopped because the homogeneity function ( 6) is zero in this region due to the large D = |x 3 | compared to the average in the already segmented region V n .The detail of this region is shown in Fig. 19.
The last experiment is the segmentation of computer classroom B306 at the Faculty of Civil Engineering of the Slovak University of Technology in Bratislava.The input point cloud with approximately 10 million points and an average resolution of about 2×2 mm is in Fig. 20.The ceiling and one wall (with windows) were removed from the original data in order to make the images more readable.The point cloud was scanned from four scanner positions.First, we performed primary segmentation to find planar point clouds.Secondary segmentation was performed only for vertical planes corresponding to the largest walls visualised in Fig. 21, left column: front wall with two whiteboards, the right one is partially covered by a projection screen; Side wall consisting of two planes, the first with a door and the second with a long wall hanger and alcove with a sink; and finally, the back wall.We also added a normal filter to remove points whose local normals deviate from the normal of the fitting plane by more than 30 degrees (coloured blue).The red points form the resulting planar point clouds.In the right column of Fig. 21, we can see corresponding smoothed edge detector functions.Note that some unwanted 'artificial' edges are also present, e.g.long inclined edge on the back wall caused by uneven illumination by natural light, edges on the whiteboard of the same origin, or edges due to merging of multiple scans (mainly on the front wall and bottom of the back wall).The final segments are shown in Fig. 22. Dark grey points represent waste from the primary segmentation (points that do not belong to any of the segmented planes).Light grey

Conclusion
In this paper, we presented a method for multichannel segmentation of planar point clouds using evolving curves.The planar point clouds were represented by a set of images, one for each available channel (colour and intensity channels and distance from the fitting plane).We discussed details of the mathematical model of the Lagrangian curve evolution, construction of the edge detectors using all channels, the design of the normal velocity driving the evolution.The model was equipped with a suitable tangential motion to improve the quality of the mesh.We derived the numerical discretisation of the mathematical model and used the fast O(n) approach to treat topological changes.The method was successfully tested on multiple examples of real data.The method could be used to automate the processing of point clouds, especially in the area of as-built documentation and geometry verification of building structures.Further development of the method can include automatic tuning of the model parameters; correction of the rough TLS intensity channel (Höfle and Pfeifer 2007;Xu et al. 2017

Fig. 1
Fig. 1 Primary segmentation of an apartment scanned by TLS.A planar subset corresponding to a wall was segmented (a) Transformed point cloud and mesh.(b) Corresponding bitmap image.

Fig. 3
Fig. 3 Creation of the bitmap images from the point cloud data visualised for R, G, B channels Smoothed edge detector function E(x).

Fig. 5
Fig. 5 Construction of the smoothed edge detector function for the red channel and the definition of the unit tangent vector T = ∂ s x. 4. Frenet-Serret formulas ∂ s T = kN, ∂ s N = −kT and orthonormality of the basis T, N (i.e.T • T = 1, N • N = 1 and T • N = 0).

Fig. 9
Fig. 9 The smoothed aggregated edge detector function E(x)

Fig. 12
Fig. 12 Curve segmentation of the wall with three plastic plates, visualisation of the selected time steps

Fig. 13
Fig. 13 The final segments of the planar point cloud

Fig. 14 A
Fig. 14 A wall in P.O.Hviezdoslav Theatre in Bratislava

Fig. 15 Fig. 16 Fig. 17 Fig. 18
Fig. 15 Left column: Bitmap images generated from the channels R, G, B, I, D of the point cloud.Right column: Corresponding edge detector functions

Fig. 20
Fig. 20 Point cloud of the classroom B306 at the Faculty of Civil Engineering of the Slovak University of Technology in Bratislava (ceiling and one wall removed)

Fig. 21
Fig. 21 Left column: primary segmentation of the point cloud.Red points belong to the resulting planar point clouds and the blue points were rejected by the normal filter.Right column: corresponding smoothed edge detector functions