Triangulación de una nube de puntos (una forma de calcularla) Parte 2

Aplicando Autolisp...

La siguiente funcion calcula el área de un triángulo proporcionándole sus tres vértices, (en estricto el área es la mitad del valor que devuelve esta función):

(defun signarea (pa pb pc)

   (+ (* (car pa) (- (cadr pb) (cadr pc)))

      (* (car pb) (- (cadr pc) (cadr pa)))

      (* (car pc) (- (cadr pa) (cadr pb))))
)

Una Implementación inicial del algoritmo:

Un algoritmo de triangulación consistente en encontrar el triángulo que contiene cada nuevo punto y, una vez hallado, formar tres nuevos triángulos, los que deben pasarse por la verificación de la condición de Delaunay, se presenta en la siguiente forma:

la variable: listadepuntos es la que contiene a la nube de puntos.

(defun tin ()
  ;;(obtener_xmax_xmin_ymax_ymin) ß para no distraernos de lo principal…
  ;; lo siguiente determina el triángulo mayor
  ;; o sea el que contiene a todos los puntos
  (setq pa (list xmin (+ ymax (* 2.0 (- ymax ymin))) 0)
        pb (list (- xmin (* 2.5 (- xmax xmin)))
                 (- ymin (* 2.5 (- ymax ymin))) 0)
        pc (list (+ xmax (* 2.5 (- xmax xmin))) ymin 0) )
  (setq listadepuntos (append (list pa pb pc) listadepuntos))

;; primeros 3 triángulos
  (setq tr1 '(3 0 1) tr2 '(3 1 2) tr3 '(3 2 0) )
  (setq grupo_triang (list tr1 tr2 tr3))
  (setq ii 4

        uu (length listadepuntos))
  (while (< ii uu)
    (setq eli (halla_tri))
    (setq grupo_triang (eliminar eli grupo_triang))
    (setq tr1 (list ii (car eli) (cadr eli))
          tr2 (list ii (cadr eli) (caddr eli))
          tr3 (list ii (caddr eli) (car eli)) )
    (setq gruponuevo (list tr1 tr2 tr3))
    (legalizar gruponuevo)

;; la siguiente línea añade los nuevos triángulos legalizados a
;; la triangulación. En el proceso de legalización se forman

;; nuevos triángulos en lugar de otros por la acción del flip.
    (setq grupo_triang (append grupo_triang gruponuevo))
    (eliminar el triangulo que contiene el punto nuevo)
    (setq ii (+ ii 1))
    )
  (dibuja_red_tin)
)

;; lo siguiente devuelve el triángulo que contiene al nuevo punto
(defun halla_tri ( / tr00 jj encontrado)
  (setq jj 0) 
  (while (not encontrado)
    (setq tr00 (nth jj grupo_triang)
          encontrado
           (and (> (signarea (list ii (car tr00) (cadr tr00)))  0.0)
                (> (signarea (list ii (cadr tr00) (caddr tr00))) 0.0)
                (> (signarea (list ii (caddr tr00) (car tr00)))) 0.0) )
          jj (1+ jj)
       )
  )
  tr00
) 

(defun signarea (a b c)
  (setq pa (nth a listadepuntos)
        pb (nth b listadepuntos)
        pc (nth c listadepuntos)
   )
  (+ (* (car pa) (- (cadr pb) (cadr pc)))
     (* (car pb) (- (cadr pc) (cadr pa)))
     (* (car pc) (- (cadr pa) (cadr pb))) ) 
 )

El código mostrado  lineas arriba complementado con las operaciones de legalización, eliminación de cada triangulo que contiene un punto nuevo y dibujado, nos produce la triangulación buscada.

El procedimiento que se aplica corresponde al algoritmo incremental aleatorio, en primer lugar se define un triangulo grande que contiene a toda la nube y luego se van añadiendo los puntos uno por uno, en cada incorporación de puntos se verifica que los triangulos formados cumplan con la condición Delaunay. Al terminar se hace el dibujado de los lados de la triangulación (líneas TIN) cuidando que no sean las que conectan los vértices del supertriángulo y que no se repitan.

De todas las operaciones a efectuar con este código; la de encontrar el triángulo que contiene a un nuevo punto se convierte en crucial si se quieren procesar miles de puntos. La función de búsqueda presentada líneas arriba es demasiado ordinaria. Es necesario implementar funciones más eficientes, una mejor opción es aprovechar una de las funciones de Visual Lisp como lo siguiente:

(defun halla_tri ( / x ) 
  (car (vl-member-if '(lambda (x)
              (and (>= (signarea (list ii (car x) (cadr x)))  0.0)
                   (>= (signarea (list ii (cadr x) (caddr x))) 0.0)
                   (>= (signarea (list ii (caddr x) (car x)))) 0.0) ) )
                      grupo_triang) )
 )

A esta solución podemos añadir un preordenamiento de la lista de puntos según las coordenadas X, tambien ayuda que los triangulos nuevos se añadan al comienzo de la lista: 

(setq grupo_triang (append gruponuevo grupo_triang) )

Sobre el "legalizado"

Una triangulación será correcta si se cumple con la condición de Delaunay, entonces cómo lo verificamos?.
Si los vértices de cada triángulo conforman una circunferencia dentro de la cual no se encuentra ningún otro punto de la nube, se dirá que se trata de una triamgulación Delaunay.
Conozco 2 maneras de efectuar la verificación:

1. Hallar el centro de la circunferencia circunscrita, por consiguiente el radio de dicha circunferencia y verificar la distancia del punto a evaluar hasta el circuncentro, si la distancia es menor que el radio, el punto está dentro de la circunferencia y debe efectuarse la modificación de la triangulación.

2. Resolver la determinante haciendo la verificación:



donde A, B, C son los vértices del triángulo y D es el punto que se desea verificar. Si resulta mayor que cero, el punto D está dentro de la circunferencia

También aquí el procedimiento debe ser eficiente. Si consideramos que no conocemos de antemano los puntos candidatos a verificar por cada triangulo, tendremos que verificar todos los lados de la triangulación lograda hasta encontrar la del triangulo que tiene el lado comun. Conforme crece la triangulación el proceso se hace más lento y se agrava si la cantidad de puntos de la nube es grande (> 2000).

Una versión del algoritmo en Autolisp se puede descargar de aqui:

triangulacion en lisp

Esta versión contiene un mejoramiento de los procedimientos y permite obtener triangulaciones más rápidas que con los procedimientos no refinados mostrados líneas arriba.

Es posible plantearse procedimientos más elaborados que mejoren la eficiencia del algoritmo. Se debe considerar la necesidad de que por cada triángulo se almacenen más registros que los 3 vértices que lo conforman. Mas adelante explicaré un planteamiento alternativo:

Cada registro de triángulo tiene los vértices que lo conforman y los "vecinos" por cada lado.

(continuará...)