A extra dimension in triangulation Deaunay / Una dimensión más en la triangulación de Deaunay

Durante la carrera, debo reconocer que no presté demasiada atención a los temas relativos a la topografía, ni al trazado de carreteras. Pero el destino es cruel y mi mayor experiencia laboral es finalmente en temas de carreteras, ferrocarril, y sus habituales subconjuntos, trazados, calidad, estructuras... No hay nada como no desear algo para que se cumpla. He aprendido a odiar en bajito para que el destino no me oiga y me la vuelva a jugar. Aún así, al empezar a trabajar tuve que toparme con el concepto de triangulación. Este concepto no ha vuelto a abandonar. Como si me persiguiera. No se tarda mucho en aprender en que consite uno de los métodos más conocidos y veteranos de triangulación (http://es.wikipedia.org/wiki/Triangulación de Delaunay, http://en.wikipedia.org/wiki/Delaunay_triangulation). La tentación era demasiado fuerte como para no intentar realizar un programa.

During studies of the faculty, I must admit I did not pay much attention to issues relating to the topography, or the layout of roads. But the destiny is cruel and my older work experience is finally in the areas of road, rail, and their usual subsets, layouts, quality structures ... Nothing like not wanting something to be fulfilled. I have learned to hate in short, so the destiny can not hear me. Still, when I started working I had to run into the concept of triangulation. This concept has not been abandoned. As if that persecute me. It does not take long to learn that consists in one of the best known and veterans methods triangulation (Delaunay http://es.wikipedia.org/wiki/Triangulación, http://en.wikipedia.org/wiki/Delaunay_triangulation) . The temptation was too strong to not attempt to make a program.

Texto del programa denaunay.c:
Text inside delaunay.c program:

#include <stdio.h>

#include <string.h>

#include <math.h>

#include <stdlib.h>

int main(int argc,char *argv[]){

FILE *fi=fopen("/.../entrada.txt", "r");

FILE *fs=fopen("/.../salida.php", "w");

char cadena[100];

char xs[10],ys[10],zs[10];

char pana[10];//sobra pero si no da error raro

double x[12000],y[12000],z[12000],v[12000];

double martir[1000];//sobra pero si no da error raro

int i=0;

int j=0;

int k=0;

int m=0;

int u=0;//sobra pero si no da error raro

char wc[1];

char *t;

char je[2]="\n";

void concate(char s1[],char s2[]);

strcpy(xs,"");strcpy(ys,"");strcpy(zs,"");

  while (fgets(wc,2,fi)!=NULL){

if (strcmp(wc,je)==0) {m=0; t=fgets(wc,2,fi);x[i]=atof(xs);y[i]=atof(ys);z[i]=atof(zs);i++;

// fprintf(fs,"x%f\n",x[i]); fprintf(fs,"y%f\n",y[i]); fprintf(fs,"z%f\n",z[i]);

strcpy(xs,"");strcpy(ys,"");strcpy(zs,"");}

if (strcmp(wc,",")==0) {m=m+1;t=fgets(wc,2,fi);}

if (m==0){concate(xs,wc);}

if (m==1){concate(ys,wc);}

if (m==2){concate(zs,wc);}

  }

// ya esta leido y llega hasta i

// for (j=0;j<i;j++) {printf("x%f\n",x[j]);printf("y%f\n",y[j]);printf("z%f\n",z[j]);}

//-----------si se le mete el recorte de la distancia se convierte en una bala

m=i;

int pluf;

int n,orden=0,desorden;

double l,d,u1,u2,v1,v2,t1,c1,c2,cp,r2,r3;

fprintf(fs,"<?php\n");

if (argv[1]==NULL) {l=100;} else {l=pow(atof(argv[1]),2)/2;}

for(i=0;i<m;i++){v[i]=(x[i]*x[i]+y[i]*y[i])/2;//fprintf(fs,"%f %f %f\n",x[i],y[i],z[i]);

}

for(i=0;i<m-2;i++){//printf("%i %i %f\n",i,orden,l);

  for(j=i+1;j<m-1;j++){d=pow(x[i]-x[j],2)+pow(y[j]-y[i],2);if (d<l) {

    for(k=j+1;k<m;k++){d=pow(x[i]-x[k],2)+pow(y[k]-y[i],2);if (d<l) {

d=0;

 if (((j!=i) && (k!=j)) && (i!=k)) {pluf=1;

  u2=x[i]-x[j];u1=y[j]-y[i];   v2=x[i]-x[k];v1=y[k]-y[i];

  t1=((x[j]-x[k])*v2-(y[j]-y[k])*v1)/(u2*v1-u1*v2)/2;

   c1=(x[i]+x[j])/2+u1*t1;    c2=(y[i]+y[j])/2+u2*t1;

  r2=(pow(c1-x[i],2)+pow(c2-y[i],2))/2;

if (r2<=l) {

  cp=(pow(c1,2)+pow(c2,2))/2;

  

     for(n=0; n<m; n++){

       if (((n!=i) && (n!=j)) && (n!=k)){

             d=(v[n]+cp-(x[n]*c1+y[n]*c2)-r2);pluf=0;

  //     fprintf(fs,": %i,%i,%i,%f.\n",i,j,k,v[n]);

  if ( (v[n]+cp-(x[n]*c1+y[n]*c2)-r2)<0 ){pluf=1;break;}

  }}

}//if r2<d

  

   if (pluf==0){fprintf(fs,"$triang[]='%i %i %i';\n",i,j,k);orden++;

   

//   fprintf(fs,"%f\n",d);

   }

}

}}}}}

fprintf(fs,"\n?>");

printf("Numero de triangulos: %i.\n",orden);

//------------------------

    return 0;

}

void concate(char s1[],char s2[]){

register int i,j;

for (i=0;s1[i];i++);

for (j=0;s1[i]=s2[j];i++,j++);

}

En este ejemplo el archivo de entrada era: "/.../entrada.txt" de la forma:
In this example the file was: "/.../data-in.txt" of this way:

42.877,8.873,477

42.826,8.952,525

42.728,8.925,687

43.211,7.347,615

43.162,7.310,690

43.147,7.258,829

43.209,7.175,898

...

y a cambio soltaba un archivo de salida: "/.../salida.php" de la forma:
And the out-file: "/.../salida.php" was of the way:

$triang[]='0 1 92';

$triang[]='0 1 668';

$triang[]='0 38 56';

$triang[]='0 38 92';

$triang[]='0 48 93';

...

El código original es del año 1993. En aquella época no había los fabulosos programas de trazado que hay ahora y, a veces, había que leer los datos en un formato determinado y exportarl los triángulos en otro formato para ser leídos por otro programa de trazado. El programa original está perdido o en algún disquette de 31/2 ó 51/4. De todas maneras, ese programa lo fui manipulando y arriba escrito es una de las versiones que tenía preparadas para ser incrustado dentro de un código de php. Mi intención es realizar un programa de trazado web. Cuando tenga tiempo.

The original code is from the year 1993. At that time there was the fabulous layout programsthat there are now and sometimes we had to read the data in a certain format and export triangles in another format to be read by another draw-program . The original program is lost or some floppy 31/2 or 51/4. Anyway, that show what was written above manipulating and is one of the versions that I had prepared to be embedded inside a php code. My intention is to make a web layout program... When I have time.

Lógicamente la sintaxis variará según el lenguaje de programación y la finalidad del programa. Para el resto de la explicación que se realiza a continuación este código ha sido variado para ser coherente con la sintaxis de VB6

Hasta aquí es importante entender el concepto de cómo se hace la triangulación de Delaunay. Hay muy buenos artículos en la wikipedia y en el resto de Internet. Y lo segundo es visualizar los bucles  for que tanto he enfatizado. Esa redundancia del cálculo nos perseguirá por todo el artículo.


Logically the syntax will change depending on programming language. and the purpose of the program. For the rest of the explanation this code has been changed to be consistent with the syntax of VB6

So here it is important to understand the concept of how the Delaunay triangulation is done. There are very good articles on Wikipedia and the rest of the Internet. And the second is to display the "for" loops that I have emphasized. The redundancy of the calculation will pursue us throughout the article.

Ampliando dimensiones / Expanding dimensions

Siempre me ha gustado pensar que si hay un proceso que vale para n=1, valdrá para n=2 y asi sucesivamente hasta n=m. Si hay un proceso para un valor entero de una variable debe haberlo para un valor real y si hay algo que valga para dos dimensiones, éstas puedan aumentarse.

I always liked to think that if a process is true for n = 1, it will be valid for n = 2 and so on up to n = m. If there is a certain process to an integer value, should also having value for real and if there is anything that is valid for two dimensions, they can be increased.

En el caso de la triangulación el mejor método en dos dimensiones es el obtenido por el método de Delaunay: En nuestro artículo "http://carreteras-laser-escaner.blogspot.com.es/2014/08/geotecnia-fractales-y-laser-escaner.html" nos volcábamos con este procedimiento para tridimensionalizar áreas de una nube de puntos.

In the case of triangulation, the best method in two dimensions is obtained by the method of Delaunay: In our article "http://carreteras-laser-escaner.blogspot.com.es/2014/08/geotecnia-fractales- and-laser-escaner.html "we chose this procedure to areas of a 3D point cloud.

A partir de una realidad:
From a reality:

Tras obtener su nube de puntos:
It obtained its point cloud:

Se triangulaba una zona:
It was triangulated a zone:

Está claro que no se puede hacer sin más el método de Delaunay pensado para dos dimensiones. Tiene que haber un truco. En este caso dos. El primer truco consistía en escoger un número limitado de puntos. Este, que es un ejemplo sencillo consta sólo de 4000 puntos. Y el área afectada contenía unos 100.

Deunay method, designed for 2D, can not be done without further. There must be a trick. In this case, two. The first trick was to choose a limited number of points. This, which is a simple example only 4000 points. And the affected area contained about 100.

La realidad es que tenemos dos problemas:

• Un elevado número de puntos
• La realidad es 3D

The reality is that we have two problems: 

• A high number of points 
• The reality is 3D

¿Cómo lo hacíamos? / How we did it?

Al seleccionar un zona eliminábamos de nuestro cómputo los puntos fuera de la zona y al girar el modelo atendíamos a su proyección frontal. De esta manera no sólo reducíamos los puntos sino que además, partiendo de la proyección, ya teníamos un modelo 2D y aplicábamos Delaunay. La proyección hay que dirigirla con cabeza ya que una vista errónea podría mandarnos al traste un solución lógica.

When we selected an area, we eliminated our external points from the counting and rotate the model were attending only to its front projection. In this way we reduced not only points but also, based on the screening, we had a model and we was able to applied 2D Delaunay. The projection must direct with the mind, because a wrong view might send ruin one logical solution.

¿Como podremos hacerlo sin trucos? / How can we do it without tricks?

De todas maneras nuestro problema es grave si nos damos cuenta como se realiza internamente el cálculo. La triangulación de Delaunay consiste en obtener las ternas de puntos de una colección de puntos en la que ponemos la condición de que no tenga ningún otro punto de la colección dentro del círculo definido.

Our problem is serious if we realize it performs the calculation internally. Delaunay triangulation is to obtain the lists of points of a collection of points in which we put the condition that they have no other collection point within the defined circle.

Cualquier programa que realice este cálculo tendrá los siguientes bucles:
Any program that will make this calculation, it'll get following loops:

Dada una colección de puntos [x_i, y_i], i Î {1..max}

for i=1 .. max -2

   for j=1 .. max -1

      for k=1 .. max

         Se halla el círculo → R,  [Cx_i, Cy_i]

            for n=1 .. max (sin i, j, ni k)

           Se compueba que  [x_n, y_n] Ï {R,  [Cx_i, Cy_i]}  

           →  se guarda la terna [i, j, k] 

              como terna de índices del triángulo obtenido.

El proceso hace, aproximadamente, max⁴ operaciones. Con unos cientos de puntos, estamos hablando de unas cien mil operaciones. Para un ordenador de sobremesa normal, serían unos pocos segundos. Pero, si pudiéramos meter todo la colección de puntos (4.000) el número de operaciones dentro del bucle sería de 256.000.000.000.000 operaciones. La cosa empieza a ser un poco complicada.

Al margen de meternos en mejoras al algoritmo anterior (que se puede hacer y se hace), nos damos cuenta de que una nube de puntos de varios millones de éstos harían el tema cuanto menos, frustrante. Y, sin meternos en consideraciones de 3 dimensiones.

The process makes about max⁴ operations. With a few hundred points, we're talking about a hundred thousand operations. For a standard desktop computer, would be a few seconds. But, if we could put all the collection of points (4000) the number of operations within the loop would be 256,000,000,000,000 operations. Things start to get a little complicated. 
Aside from getting into improvements to the previous algorithm (which can be done and done), we realize that a scatter plot of several million of these would make the issue less is frustrating. And, without getting into 3-dimensional considerations.

Paso a paso / Step by step

Una primera aproximación sería no triangular. Buscaríamos símplemente los puntos más cercanos y los uniríamos entre sí. Haríamos una colección de segmentos. A cada punto [x_i, y_i] le corresponde un punto [x_j, y_j] que resulta que es el punto más cercano.

A first approach would not triangular. Simply would search the closest points and would unite together. We would do a collection of segments. At each point [x_i, y_i] corresponds to a point [x_j, y_j] will be the closest point.

SEGMENTOS / SEGMENTS

Esta colección de segmentos si la exportamos en formato DXF tendría esta apariencia:

This collection of segments exported in DXF format would look as this:

Obviamente no es una triangulación, pero es muy rápida con max² operaciones. no hay que despreciarla sino que hay que buscarle una utilidad.

It is obviously not a triangulation, but is very quick with max²  operations. This should not despise but must seek a profit.

UNOS CUANTOS TRIÁNGULOS (FRACTAL) / FEW TRIANGLES (FRACTAL)

La idea es sólo tomar los triángulos más fáciles. No es una triangulación completa, pero es muy rápida (max² operaciones).

The idea is to just take the easiest triangles. It is not a complete triangulation, but is very fast (max² operations).

Y los triángulos que se obtienen son los siguientes:

And the triangles obtained are follows:

Es una triangulación incompleta, pero es muy rápida ((2 · max)² operaciones). Para el caso que estudiábamos de la geotécnia, si bien no se tienen todos los triángulos con tantos puntos al menos tenemos una buena cantidad de ellos. es lo que denominábamos triangulación fractal. No nos vale como triangulación pero sí nos da información de cuales son las orientaciones de los planos y nos permiten hacer estadística.

It is an incomplete triangulation, but is very fast  ((2 · max)² operaciones). In the case we studied the geotechnical, although not all triangles are many points at least we have a good amount of them. it is what we called fractal triangulation. It is not a triangulation but it gives us information which the orientations of the planes and allow us to do statistics.

Continua en: http://carreteras-laser-escaner.blogspot.com/2014/11/a-extra-dimension-in-triangulation_13.html
Continued on: http://carreteras-laser-escaner.blogspot.com/2014/11/a-extra-dimension-in-triangulation_13.html

Geometric road runoff estimation from laser mobile mapping data

Artículo patrocinado por Extraco, Misturas, Lógica, Enmacosa e Ingeniería InSitu, dentro del proyecto SITEGI, cofinanciado por el CDTI. (2012). 

Article sponsored by Extraco, Misturas, Lógica, Enmacosa and Ingeniería Insitu inside the SITEGI project, cofinanced by the CDTI. (2012)

Para más información o para solicitar el documento original, contacte con nosotros:
If you would like more information or request the original document, contact with us:

http://carreteras-laser-escaner.blogspot.com.es/p/contacte-con-nosotros.html

Autores / Authors : Jinhu Wanga,(c), ,Higinio González-Jorge(b), Roderik Lindenbergh(a), Pedro Arias-Sánchez (b) and Massimo Menenti(a)

a) Dept. of Geoscience and Remote Sensing, Delft University of Technology Building 23, Stevinweg 1, Post Box 5048, 2628CN Delft, The Netherlands (jinhu.wang, r.c.lindenbergh, m.menenti)@tudelft.nl

b) Dept. of Natural Resources and Environmental Engineering, School of Mining Engineering, University of Vigo, E-36310 Vigo, Spain (higiniog, parias)@uvigo.es

c) Key Laboratory of Quantitative Remote Sensing Information Technology Academy of Opto-Electronics, Chinese Academy of Sciences No. 9 Deng Zhuang South Road, HaiDian District, 100094 Beijing, China

ABSTRACT:

Mountain roads are the lifelines of remote areas but are often situated in complicated settings and prone to landslides, rock fall, avalanches and damages due to surface water runoff. The impact and likelihood of these types of hazards can be partly assessed by a detailed geometric analysis of the road environment. Field measurements in remote areas are expensive however. A possible solution is the use of a Laser Mobile Mapping System (LMMS) which, at high measuring rate, captures dense and accurate point clouds.

This paper presents an automatic approach for the delineation of both the direct environment of a road and the road itself into local catchments starting from a LMMS point cloud. The results enable a user to assess where on the road most water from the surroundings will assemble, and how water will flow over the road after e.g. heavy snow melt or rainfall. To arrive at these results the following steps are performed. First outliers are removed and point cloud data is gridded at a uniform width. Local surface normal and gradient of each grid point are determined. The relative smoothness of the road is used as a criterion to identify the road’s outlines. The local gradients are input for running the so-called D8 method, which simply exploits that surface water follows the direction of steepest descent. This method first enables the identification of sinks on the roadside, i.e. the locations where water flow accumulates and potentially enters the road. Moreover, the method divides the road’s direct neighbourhood into catchments, each consisting of all grid cells having runoff to the same sink. In addition the method is used to analyse the surface flow over the road’s surface. The new method is demonstrated on a piece of 153 meters long Galician mountain road as sampled by MMS data.

1 INTRODUCTION

Light Detection And Ranging (LiDAR) surveying techniques enable to quickly obtain 3D geometry. Notably a Laser Mobile Mapping System (LMMS) which integrates a Global Navigation Satellite System (GNSS), an Inertial Measuring Unit (IMU) and LiDAR profilers on a moving platform, enables efficient and complete 3D data collection (Vosselman and Maas, 2010). Most applications of LMMS data focus on cities, but there are also applications considering highway surveying. For example, in (Kukko et al., 2009), point cloud and image data acquired by a LMMS is used to classify and model the road environment in a fully automatic way. In (Bitenc et al., 2011) a LMMS is used to generate a Digital Terrain Model (DTM) of a sandy beach of 6 km long in a study evaluating the possibility of LMMS for beach erosion assessment.

LMMS point cloud data has also been combined with Airborne Laser Scanning (ALS) data to map curbstones (Zhou and Vosselman, 2012, Tao, 2000). And based on LMMS point cloud data, an automatic feature extraction approaches were developed to extract basic road structures, e.g. lamp poles, road signs, lanes and crosswalks (Jae-Seung et al., 2007, Foy et al., 2007, Mancini et al., 2012, Pu et al., 2011).

For mountainous rural areas, roads are lifelines to the civilians and the safety of the road and its environment is an important concern. The safety and condition of the roads need regular inspection and monitoring for security reasons. LMMS ranging provides the possibility to sample the geometry of a road and its surroundings in an almost continuous way. The resulting point cloud data contains information that can serve as input for flood hazard and landslide prediction. In (Poppenga et al., 2010), a Digital Elevation Model (DEM) is constructed from point cloud data to model surface flow, and was applied to flood inundation and erosion estimation. Also in (Kazuhiro et al., 2005, White et al., 2010, Ziegler and Giambelluca, 1997) high-resolution DEM data is generated to predict surface erosion and to estimate the amount of sediment drained by streams. Especially for mountainous roads, rocks on roadside hills could fall down and cause risks. Also, water flow may cause erosion at the side of the road, eventually resulting in road damage. Moreover, steep and unstable roadsides may cause landslides resulting in further road damage.

In this work both the road itself and the roadside environment are considered. This paper computes the roadside environment catchments and estimates where and how water would flow over the surface. To some extent, rock fall is expected to follow the water flow direction as well. Firstly LMMS point cloud data sampling a mountainous road was downsampled and the outliers and noisy points were removed. Based on the data, normal vectors, as well as the 2D slope were estimated at every point. Then, an automatic iterative floating window approach is introduced that takes advantage of point height, normal vector and slope to identify points sampling the road surface. After that, the D8 algorithm is used to estimate the water flow direction on the roadside. Based on these directions, the road environment is divided in runoff sections.

2 METHODOLOGY

The methodology described in this section aims at estimating water runoff of roadsides based on a LMMS point cloud dataset. The method consists of four steps. (1) Point cloud pre-processing. 

The original point cloud data have very high density and is downsampled before processing. Outliers are also removed as well. (2) Local surface normal and 2D slope estimation. (3) Road delineation, which allows to decompose the point cloud into road and roadside points. (4) Runoff estimation. The processing flowchart is illustrated in figure 1.

Figure 1: Method for estimation runoff from LMMS point cloud data

The original point cloud has very high point cloud density, thus for efficient processing purpose, a downsampling procedure is performed using a uniform voxel size. Also, the outliers were removed before processing by using a neighbourhood point cloud mean density criterion. Details on these procedures are given in (Wang et al., 2013). A grid point is estimated from the point cloud points within the voxel using inverse distance interpolation with power 2.

2.2 Surface normal estimation

Surface normals are estimated during the iterative filtering of the road points. The normal at a certain discrete point is defined as a vector perpendicular to the tangential plane of the local surface at that point (Th¨urmer and W¨uthrich, 1997, Dey et al., 2005). 

In this work, surface normals are estimated from neighbouring points. For each point in the point cloud, radius neighbourhood searching was performed to acquire the points within a preset radius of the query point. From the points found, the local normals are computed as described in (Wang et al., 2013).

2.3 2D slope computation

The 2D slope, also known as 2D gradient, is a vector field of a surface. The vector direction points to the greatest change in height, and the vector magnitude is the rate of change. The vector direction is also referred to as aspect or local surface orientation. 

A first approximation of the 2D slope in a regular grid is obtained by selecting the largest 1D slope in one of the eight neighbouring directions. Suppose the grid size is w, then the 1D slope Si in the direction of each of the eight neighbouring grid cells is given by:

where Hi is the elevation of the i − th neighbour of the query point, and hq is the elevation of the query point itself, while di is the horizontal distance from the query point to the i − th neighbour. Note that di is √(2ω) in diagonal direction.

2.4 Roadside points segmentation

In this work, roadsides points are segmented from the original point cloud by an automatic iterative point cloud segmentation approach, which takes the normals and 2D slopes as obtained in the previous steps as input. A moving window was used to obtain for each point its local normal and relative height. Next, a procedure is followed to segment points into road and off-road points.

Then the same process is repeated with a smaller window size until the smallest window size threshold is reached. The details are given in (Wang et al., 2013).

2.5 D8 algorithm

The D8 algorithm introduced by (O’Callaghan and Mark, 1984), is a grid based algorithm and is widely used due to its simplicity. For a given query grid point, the D8 algorithm approximates the primary flow direction by choosing the direction to the neighbour with maximal 2D gradient, as illustrated in Figure 2.

For example, the flow direction from the central pixel, with value 16, is downward, because the gradient towards the pixel directly below, with value 11, is maximal among the eight neighbours of the central pixel. In the next step of the D8 algorithm, the flow is followed. In Figure 2, all flow eventually terminates at the pixels in the bottom row.

Applying this method on all the roadside pixels results in a decomposition of the sampled roadside into different catchments.

Large catchments correspond to a large local water inflow at the sink of the catchment, as shown in Figure 3, and the area is defined as upstream catchment area. The flow direction is determined for each pixel and pixels that are flowing towards the same bottom pixel, the sink, are assigned in the same randomly allocated colour.

In this work, the original downsampled point cloud data is organized in a uniform grid, and the height assigned to a grid cell is the mean height of all points belonging to the grid cell. Each grid cell is potentially surrounded by eight neighbouring grid cells. 

The gradient for each of these eight directions is obtained using Equation 1. Then the D8 algorithm is applied to the gridded point cloud to compute the local flow directions and, by accumulating flow to consecutively compute catchments and sinks.

2.4 Roadside points segmentation

In this work, roadsides points are segmented from the original point cloud by an automatic iterative point cloud segmentation approach, which takes the normals and 2D slopes as obtained in the previous steps as input. A moving window was used to obtain for each point its local normal and relative height. Next, a procedure is followed to segment points into road and off-road points. 

Then the same process is repeated with a smaller window size until the smallest window size threshold is reached. The details are given in (Wang et al., 2013).

2.5 D8 algorithm

The D8 algorithm introduced by (O’Callaghan and Mark, 1984), is a grid based algorithm and is widely used due to its simplicity. 

For a given query grid point, the D8 algorithm approximates the primary flow direction by choosing the direction to the neighbour with maximal 2D gradient, as illustrated in Figure 2.

For example, the flow direction from the central pixel, with value 16, is downward, because the gradient towards the pixel directly below, with value 11, is maximal among the eight neighbours of the central pixel. In the next step of the D8 algorithm, the flow is followed. In Figure 2, all flow eventually terminates at the pixels in the bottom row.

Applying this method on all the roadside pixels results in a decomposition of the sampled roadside into different catchments.

Large catchments correspond to a large local water inflow at the sink of the catchment, as shown in Figure 3, and the area is defined as upstream catchment area. The flow direction is determined for each pixel and pixels that are flowing towards the same bottom pixel, the sink, are assigned in the same randomly allocated colour.

In this work, the original downsampled point cloud data is organized in a uniform grid, and the height assigned to a grid cell is the mean height of all points belonging to the grid cell. Each grid cell is potentially surrounded by eight neighbouring grid cells. 

The gradient for each of these eight directions is obtained using Equation 1. Then the D8 algorithm is applied to the gridded point cloud to compute the local flow directions and, by accumulating flow to consecutively compute catchments and sinks.

Figure 2: D8 algorithm flow directions

Figure 3: Upstream catchment area of grid cells

3.1 Data description and Implementation

The point cloud data analysed in this paper is acquired with the Laser Mobile Mapping System of the University of Vigo, Spain. The study area is a piece of 153 meters long mountainous road, as shown in Figure 4.

Figure 4: Photo of the studied road

The LMMS used was the Lynx Mobile Mapper from OPTECH. The Lynx contains two LiDAR profilers collecting LiDAR point cloud data at 500,000 measurements per second with 360 degree field of view (Puente et al., 2013b, Puente et al., 2013a). All data is geo-referenced by differential GPS post processing. 

Figure 5: Original LMMS point cloud data in a 3D side view

The coordinate system used is UTM-WGS84. The original point cloud dataset contains 5,838,794 points and has an average point density of 2,084 points per square meter. This particular location has suffered from rock fall and landslides along the roadside slope. Figure 5 depicts the original point cloud.of the studied road, the points are colourized by elevation. Figure 6 is a Triangulated Irregular Network (TIN) generated from the point cloud data of the studied road, which has steep cliffs on both sides. All the procedures are implemented using C++.

Figure 6: TIN generated from the road point cloud

3.2 Roadside points segmentation

Following the methods described in Section 2, the point cloudwas filtered and voxelized using a uniform width of 0.1 meter. Then the point cloud was segmented and decomposed into three parts: Road points, Northern roadside and Southern roadside points. This is illustrated in Figure 7. The points in blue are road points, while the points in red and green are the northern and the southern roadside points respectively.

3.3 Catchments estimation results

Before the application of the D8 method to obtain the catchments from the roadside slopes, a uniformed size grid was generated from the point cloud data. In this work, the grid size was preset to 2.0 m. The on-road water flow directions are estimated, as shown in Figure 8. In the figure, the flow directions are denoted by arrows. Using the D8 method, the road is divided in catchments, which are indicated in Figure 8 by different colours.

Dark cell have no outflow. After the flow directions on the road were determined, the flow directions for off-road point cloud data are also estimated, as shown in Figure 9. In this figure, the sinks are colourized gray and all cells eventually flowing to the same sink are colourized by the same colour. Each sink is labelled by a digit. The number of grid cells having runoff to each labelled

Figure 7: Roadside points segmentation result

Figure 8: Road water flow directions on each grid cells

sink in Figure 9 is given in Figure 10. There are 25 sinks on the southern roadside and 29 on the northern roadside respectively. The results shows for example that the sink labelled as No. 15 on the north roadside, has 44 contributing cells, which indicates that this sink has a lot of potential water inflow. Comparison to the original terrain model in Figure 5 shows that this sink is actually located directly below the landslide area also shown in Figure 4. The location of this sink is indicated in Figure 5 by a green ellipse. The shape of the terrain at this location is indeed such that more water is expected to accumulate. On the other hand, the sink labelled as 26 has only 5 contributing grid cells. 

Figure 9: Labelled roadside catchment area

Figure 10: Count of catchment area cells for both roadsides

And indeed, at this location, the roadside is very steep and water flows directly on the road. In Figure 11, the amount of saturation of the grid cells conrresponds to the flow accumulatation. That is, a cell with a high colour saturation collects water from many cells. This also denotes water flow direction.

4 CONCLUSION AND RECOMMENDATION

Since mountainous roads have complicated morphological environments and face threat from landslides and rock fall, there is a need for road and road environment safety inspection and monitoring. 

To meet this obligation, detailed and continuous road environment surface flow modelling has to be acquired. LMMS can acquire point clouds in an efficient way, both from a time and costs perspective. For this reason we have presented a method to estimate roadside properties, which are the gradient and slope, and then the catchments on the roadside slope are computed with the D8 algorithm. The number of cells in each catchments is a measure for the amount of water flow into the corresponding road surface location.

Figure 11: Accumulated inflow of each grid cell

In this work, the cell size was set to 2 meters only for the feasibility demonstration of the D8 method in the catchments and runoff estimation. But for practical and high quality purpose, the resolution could be much higher, e.g. up to 25 cm, as long as the point cloud density in the original data is high enough. Also, to validate the catchment estimation results, other data sets could be introduced, like airborne laser scanning data, total station surveying or GNSS profiling of the terrain.

To evalute the results, other Geopraphy Information System (GIS) software could be used to evaluate the flow direction and compare the results. A future work would be the monitoring of the sink locations, and to inspect if local road erosion is correlated with the size of the inflowing roadside catchment. Note that the D8 method as presented here, requires a non-trivial slope. That is, if the surface off or on the road is locally flat, the method would be stuck. A possible solution is to take the expected speed and direction of water flow into account.

REFERENCES

Bitenc, M., Lindenbergh, R., Khoshelham, K. and Van Waarden, A. P., 2011. Evaluation of a LIDAR Land-Based Mobile Mapping System for Monitoring Sandy Coasts. Remote Sensing 3(12), pp. 1472–1491.

Dey, T., Li, G. and Sun, J., 2005. Normal estimation for point clouds: a comparison study for a voronoi based method. In: Point-Based Graphics, 2005. Eurographics/IEEE VGTC Symposium Proceedings, pp. 39–46.

Foy, S., Deegan, C., Mulvihill, C., Fitzgerald, C., Markham, C. and McLoughlin, S., 2007. Road sign safety identification through the use of a mobile survey system. Vol. XXXVI-5Number C55, International Society of Photogrtammetry and Remote Sensing.

Jae-Seung, J., Jae-Min, P., Dong-Hun, J. and Byung-Guk, K., 2007. Automatic identification of road sign in mobile mapping system. Vol. XXXVI-5Number C55, International Society of Photogrtammetry and Remote Sensing.

Kazuhiro, A., John, S. and S., M. E., 2005. Forest road design with soil sediment evaluation using a high-resolution DEM. Journal of Forest Research 10(6), pp. 471–479.

Kukko, A., Jaakkola, A., Lehtomaki, M., Kaartinen, H. and Chen, Y., 2009. Mobile mapping system and computing methods for modelling of road environment. In: Urban Remote Sensing Event, 2009 Joint, pp. 1–6.

Mancini, A., Frontoni, E. and Zingaretti, P., 2012. Automatic road object extraction from mobile mapping systems. In: Mechatronics and Embedded Systems and Applications (MESA), 2012 IEEE/ASME International Conference on Advanced Intelligent Mechatronics (AIM), pp. 281–286.

O’Callaghan, J. F. and Mark, D. M., 1984. The extraction of drainage networks from digital elevation data. Computer Vision, Graphics, and Image Processing 28(3), pp. 323 – 344.

Poppenga, S. K., Worstell, B. B., Stoker, J. M. and Greenlee, S. K., 2010. Using selective drainage methods to extract continuous surface flow from 1-meter lidar-derived digital elevation data. Technical report, U.S. Geological Survey Scientific Investigations Report.

Pu, S., Rutzinger, M., Vosselman, G. and Elberink, S. O., 2011. Recognizing basic structures from mobile laser scanning data for road inventory studies. ISPRS Journal of Photogrammetry and Remote Sensing 66(6, Supplement), pp. S28 – S39.

Puente, I., González-Jorge, H., Riveiro, B. and Arias, P., 2013a. Accuracy verification of the Lynx Mobile Mapper system. Optics and Laser Technology 45, pp. 578–586.

Puente, I., González-Jorge, H., Martnez-Sánchez, J. and Arias, P., 2013b. Review of mobile mapping and surveying technologies. Measurement 46(7), pp. 2127 – 2145.

Tao, C. V., 2000. Mobile mapping technology for road network data acquisition. Journal of Geospatial Engineering 2, pp. 1–14.

Thürmer, G. and Wüthrich, C. A., 1997. Normal computation for discrete surfaces in 3D space. Computer Graphics Forum 16, pp. C14–C26.

Vosselman, G. and Maas, H.-G., 2010. Airborne and Terrestrial Laser Scanning. Vol. XXXVI, Whittles Publishing. 

Wang, J., González-Jorge, H., Lindenbergh, R., Arias-Sánchez, P. and Menenti, M., 2013. Automatic Estimation of Excavation Volume from Laser Mobile Mapping Data for Mountain Road Widening. Remote Sensing 5, pp. 4629–4651.

White, R. A., Dietterick, B. C., Mastin, T. and Strohman, R., 2010. Forest roads mapped using lidar in steep forested terrain. Remote Sensing 2(4), pp. 1120–1141.

YLLs Gip-Krop, A. Ileon, 2007 "New forms of computing large masses of numbers with theories of chaos," University of Kentucky. 
Navi Anait, Tat Siul 2011 "The chaos computation ."

Zhou, L. and Vosselman, G., 2012. Mapping curbstones in airborne and mobile laser scanning data. International Journal of Applied Earth Observation and Geoinformation 18, pp. 293–304. 

Ziegler, A. D. and Giambelluca, T.W., 1997. Importance of rural roads as source areas for runoff in mountainous areas of northern Thailand. Journal of Hydrology 196, pp. 204–229.