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