Триангуляция полигона

Автор Тема: Триангуляция полигона  (Прочитано 10712 раз)

0 Пользователей и 1 Гость просматривают эту тему.

Оффлайн BearDyuginАвтор темы

  • ADN Club
  • ****
  • Сообщений: 270
  • Карма: 24
  • Геодезист
Триангуляция полигона
« : 04-11-2016, 13:44:55 »
Здравствуйте!

В общем столкнулись с тем что одна функция не всегда адекватно обрабатывает полигоны, особенно не выпуклые.
Было придумано решения разбить полигон на треугольники.

Нашёл описание алгоритма http://grafika.me/node/12
Цитировать
Алгоритм триангуляции:
1. Берем три вершины A1, A2, A3
2. Проверяем образуют ли вектора A1A3, A1A2 и их векторное произведение левую тройку векторов.
3. Проверяем нет ли внутри треугольника A1A2A3 какой-либо из оставшихся вершин многоугольника.
4. Если оба условия выполняются, то строим треугольник A1A2A3, а вершину A2 исключаем из многоугольника, не трогая вершину A1, сдвигаем вершины A2 (A2 на A3), A3 (A3 на A4)
5. Если хоть одно условие не выполняется, переходим к следующим трем вершинам.
6. Повторяем с 1 шага, пока не останется три вершины.

Рис.1 Алгоритм триангуляции невыпуклого многоугольника

На рисунке 1
• треугольник A1A2A3 удовлетворяет обоим условиям (п.2, п.3);
• треугольник A2A3A4 не удовлетворяет условию (п.2);
• треугольник A3A4A5 не удовлетворяет условию (п.3).

и даже код к нему на C# http://grafika.me/node/457
Код - C# [Выбрать]
  1. using System;
  2. using System.Drawing;
  3.  
  4. namespace RobotInLabuteny
  5. {
  6.     class Polygon
  7.     {
  8.         private PointF[] points; //вершины нашего многоугольника
  9.         private Triangle[] triangles; //треугольники, на которые разбит наш многоугольник
  10.         private bool[] taken; //была ли рассмотрена i-ая вершина многоугольника
  11.  
  12.         public Polygon(float[] points) //points - х и y координаты
  13.         {
  14.             if (points.Length % 2 == 1 || points.Length < 6)
  15.                 throw new Exception(); //ошибка, если не многоугольник
  16.  
  17.             this.points = new PointF[points.Length / 2]; //преобразуем координаты в вершины
  18.             for (int i = 0; i < points.Length; i += 2)
  19.                 this.points[i / 2] = new PointF(points[i], points[i + 1]);
  20.  
  21.             triangles = new Triangle[this.points.Length - 2];
  22.  
  23.             taken = new bool[this.points.Length];
  24.  
  25.             triangulate(); //триангуляция
  26.         }
  27.  
  28.         private void triangulate() //триангуляция
  29.         {
  30.             int trainPos = 0; //
  31.             int leftPoints = points.Length; //сколько осталось рассмотреть вершин
  32.  
  33.             //текущие вершины рассматриваемого треугольника
  34.             int ai = findNextNotTaken(0);
  35.             int bi = findNextNotTaken(ai + 1);
  36.             int ci = findNextNotTaken(bi + 1);
  37.  
  38.             int count = 0; //количество шагов
  39.  
  40.             while (leftPoints > 3) //пока не остался один треугольник
  41.             {
  42.                 if (isLeft(points[ai], points[bi], points[ci]) && canBuildTriangle(ai, bi, ci)) //если можно построить треугольник
  43.                 {
  44.                     triangles[trainPos++] = new Triangle(points[ai], points[bi], points[ci]); //новый треугольник
  45.                     taken[bi] = true; //исключаем вершину b
  46.                     leftPoints--;
  47.                     bi = ci;
  48.                     ci = findNextNotTaken(ci + 1); //берем следующую вершину
  49.                 }
  50.                 else { //берем следующие три вершины
  51.                     ai = findNextNotTaken(ai + 1);
  52.                     bi = findNextNotTaken(ai + 1);
  53.                     ci = findNextNotTaken(bi + 1);
  54.                 }
  55.  
  56.                 if (count > points.Length * points.Length) { //если по какой-либо причине (например, многоугольник задан по часовой стрелке) триангуляцию провести невозможно, выходим
  57.                     triangles = null;
  58.                     break;
  59.                 }
  60.  
  61.                 count++;
  62.             }
  63.  
  64.             if (triangles != null) //если триангуляция была проведена успешно
  65.                 triangles[trainPos] = new Triangle(points[ai], points[bi], points[ci]);
  66.         }
  67.  
  68.         private int findNextNotTaken(int startPos) //найти следущую нерассмотренную вершину
  69.         {
  70.             startPos %= points.Length;
  71.             if (!taken[startPos])
  72.                 return startPos;
  73.  
  74.             int i = (startPos + 1) % points.Length;
  75.             while (i != startPos)
  76.             {
  77.                 if (!taken[i])
  78.                     return i;
  79.                 i = (i + 1) % points.Length;
  80.             }
  81.  
  82.             return -1;
  83.         }
  84.  
  85.         private bool isLeft(PointF a, PointF b, PointF c) //левая ли тройка векторов
  86.         {
  87.             float abX = b.X - a.X;
  88.             float abY = b.Y - a.Y;
  89.             float acX = c.X - a.X;
  90.             float acY = c.Y - a.Y;
  91.  
  92.             return abX * acY - acX * abY < 0;
  93.         }
  94.  
  95.         private bool isPointInside(PointF a, PointF b, PointF c, PointF p) //находится ли точка p внутри треугольника abc
  96.         {
  97.             float ab = (a.X - p.X) * (b.Y - a.Y) - (b.X - a.X) * (a.Y - p.Y);
  98.             float bc = (b.X - p.X) * (c.Y - b.Y) - (c.X - b.X) * (b.Y - p.Y);
  99.             float ca = (c.X - p.X) * (a.Y - c.Y) - (a.X - c.X) * (c.Y - p.Y);
  100.  
  101.             return (ab >= 0 && bc >= 0 && ca >= 0) || (ab <= 0 && bc <= 0 && ca <= 0);
  102.         }
  103.  
  104.         private bool canBuildTriangle(int ai, int bi, int ci) //false - если внутри есть вершина
  105.         {
  106.             for (int i = 0; i < points.Length; i++) //рассмотрим все вершины многоугольника
  107.                 if (i != ai && i != bi && i != ci) //кроме троих вершин текущего треугольника
  108.                     if (isPointInside(points[ai], points[bi], points[ci], points[i]))
  109.                         return false;
  110.             return true;
  111.         }
  112.  
  113.         public PointF[] getPoints() //возвращает вершины
  114.         {
  115.             return points;
  116.         }
  117.  
  118.         public Triangle[] getTriangles() //возвращает треугольники
  119.         {
  120.             return triangles;
  121.         }
  122.  
  123.     }
  124.  
  125.     public class Triangle //треугольник
  126.     {
  127.         private PointF a, b, c;
  128.  
  129.         public Triangle(PointF a, PointF b, PointF c)
  130.         {
  131.             this.a = a;
  132.             this.b = b;
  133.             this.c = c;
  134.         }
  135.  
  136.         public PointF getA()
  137.         {
  138.             return a;
  139.         }
  140.  
  141.         public PointF getB()
  142.         {
  143.             return b;
  144.         }
  145.  
  146.         public PointF getC()
  147.         {
  148.             return c;
  149.         }
  150.     }
  151.  
  152.  
  153. }

В целом всё понятно и не вижу проблем перенести его на Lisp

Но что-то мне подсказывает, что алгоритм этот будет верным лишь для полигонов с обходом по часовой стрелки.

Вопросы:
1) Вдруг где-то видели готовое решение на Lisp, чтоб велосипед не изобретать
2) Как определить по часовой или против у нас список вершин полигона?


Оффлайн Александр Ривилис

  • Administrator
  • *****
  • Сообщений: 13886
  • Карма: 1788
  • Рыцарь ObjectARX
  • Skype: rivilis
Re: Триангуляция полигона
« Ответ #1 : 04-11-2016, 14:52:08 »
Я так понимаю, что тебя интересует триангуляция Дэлоне. Посмотри алгоритм Евгения Елпанова: http://elpanov.com/index.php?id=6
Не забывайте про правильное Форматирование кода на форуме
Создание и добавление Autodesk Screencast видео в сообщение на форуме
Если Вы задали вопрос и на форуме появился правильный ответ, то не забудьте про кнопку Решение

Оффлайн BearDyuginАвтор темы

  • ADN Club
  • ****
  • Сообщений: 270
  • Карма: 24
  • Геодезист
Re: Триангуляция полигона
« Ответ #2 : 04-11-2016, 15:23:06 »
Я так понимаю, что тебя интересует триангуляция Дэлоне. Посмотри алгоритм Евгения Елпанова: http://elpanov.com/index.php?id=6
Нет, не просто триангуляция в чистом виде, а триангуляция конкретно полигона.
Даже сделать просто триангуляции, а потом выкинуть все треугольники вне полигона не вариант, т.к. триангуляция штука не однозначная, она может построить треугольники которые будут пересекать периметр полигона :-(
Ну и второй момент, я смотрел тот код Евгения, ещё лет 10 назад... Но к своему стыду так и не сумел понять, не то что код, но и саму теорию Делона.
Но в данной ситуации мне достаточно того кода что я привёл. Осталось только разобраться с тем как мы получили вершины полигона, по часовой или против, чтобы к ним либо сразу применить эту функцию или сначала сделать списку вершин реверс.
Хотя подозреваю, что если у нас полигон обходит против часовой, данная функция тупо не сможет его решить :-)

Оффлайн Александр Ривилис

  • Administrator
  • *****
  • Сообщений: 13886
  • Карма: 1788
  • Рыцарь ObjectARX
  • Skype: rivilis
Re: Триангуляция полигона
« Ответ #3 : 04-11-2016, 15:35:58 »
У Евгения есть и функции, которые определяют направление обхода полилинии. Поищи.
Не забывайте про правильное Форматирование кода на форуме
Создание и добавление Autodesk Screencast видео в сообщение на форуме
Если Вы задали вопрос и на форуме появился правильный ответ, то не забудьте про кнопку Решение

Оффлайн BearDyuginАвтор темы

  • ADN Club
  • ****
  • Сообщений: 270
  • Карма: 24
  • Геодезист
Re: Триангуляция полигона
« Ответ #4 : 04-11-2016, 16:31:00 »
У Евгения есть и функции, которые определяют направление обхода полилинии. Поищи.
Нашёл на Caduser'e хороший бфл форум  :'(
от VovKa
Код - Auto/Visual Lisp [Выбрать]
  1. (defun vk_IsBoundaryClockwise (CoordsList / YList)
  2.   (setq YList (mapcar 'cadr CoordsList))
  3.   (minusp
  4.     (/ (apply '+
  5.          (mapcar (function (lambda (xi yi-1 yi+1) (* xi (- yi+1 yi-1))))
  6.             (mapcar 'car CoordsList)
  7.             (cons (last YList) YList)
  8.             (append (cdr YList) (list (car YList)))
  9.          )
  10.        )
  11.        2
  12.     )
  13.   )
  14. )
  15. (vk_IsBoundaryClockwise
  16.   (mapcar 'cdr
  17.      (vl-remove-if-not
  18.        '(lambda (x) (= (car x) 10))
  19.        (entget (car (entsel)))
  20.      )
  21.   )
  22. )
и от Евгения
Код - Auto/Visual Lisp [Выбрать]
  1. (defun lwcl( lw  / LST MAXP MINP)
  2. (if (= (type lw) 'ENAME)
  3.     (setq lw (vlax-ename->vla-object lw)))
  4. (vla-GetBoundingBox lw 'MinP 'MaxP)
  5. (setq
  6. minp(vlax-safearray->list minp)
  7. MaxP(vlax-safearray->list MaxP)
  8. lst(mapcar(function(lambda(x)
  9. (vlax-curve-getParamAtPoint lw
  10. (vlax-curve-getClosestPointTo lw x))))
  11. (list minp(list(car minp)(cadr MaxP))
  12. MaxP(list(car MaxP)(cadr minp)))))
  13. (if(or
  14. (<=(car lst)(cadr lst)(caddr lst)(cadddr lst))
  15. (<=(cadr lst)(caddr lst)(cadddr lst)(car lst))
  16. (<=(caddr lst)(cadddr lst)(car lst)(cadr lst))
  17. (<=(cadddr lst)(car lst)(cadr lst)(caddr lst))) t nil))
В качестве полигона у меня список вершин, думаю лучше подойдёт вариант от VovKa

Оффлайн BearDyuginАвтор темы

  • ADN Club
  • ****
  • Сообщений: 270
  • Карма: 24
  • Геодезист
Re: Триангуляция полигона
« Ответ #5 : 06-11-2016, 20:40:09 »
В общем получилась вот такая функция
Код - Auto/Visual Lisp [Выбрать]
  1. (defun polygon_triangul (point_list      /
  2. ;;;local functions
  3.                          vk_IsBoundaryClockwise
  4.                          is_left         is_point_Inside
  5. ;;;
  6.                          new_list        triangl_lict
  7.                          a               b
  8.                          c
  9.                         )
  10. ;;;
  11. ;;;vk_IsBoundaryClockwise
  12. ;;;
  13.   (defun vk_IsBoundaryClockwise (CoordsList / YList)
  14.     (setq YList (mapcar 'cadr CoordsList))
  15.     (minusp
  16.       (/
  17.         (apply
  18.           '+
  19.           (mapcar (function (lambda (xi yi-1 yi+1) (* xi (- yi+1 yi-1))))
  20.                   (mapcar 'car CoordsList)
  21.                   (cons (last YList) YList)
  22.                   (append (cdr YList) (list (car YList)))
  23.           )
  24.         )
  25.         2
  26.       )
  27.     )
  28.   )
  29. ;;;
  30. ;;;is_left
  31. ;;;
  32.   (defun is_left (a b c)
  33.     (minusp
  34.       (-
  35.         (* (- (car b) (car a)) (- (cadr c) (cadr a)))
  36.         (* (- (car c) (car a)) (- (cadr b) (cadr a)))
  37.       )
  38.     )
  39.   )
  40. ;;;
  41. ;;;is_point_Inside
  42. ;;;
  43.   (defun is_point_Inside (a b c p / ab bc ca)
  44.     (setq
  45.       ab (-
  46.            (*
  47.              (- (car a) (car p))
  48.              (- (cadr b) (cadr a))
  49.            )
  50.            (*
  51.              (- (car b) (car a))
  52.              (- (cadr a) (cadr p))
  53.            )
  54.          )
  55.       bc (-
  56.            (*
  57.              (- (car b) (car p))
  58.              (- (cadr c) (cadr b))
  59.            )
  60.            (*
  61.              (- (car c) (car b))
  62.              (- (cadr b) (cadr p))
  63.            )
  64.          )
  65.       ca (-
  66.            (*
  67.              (- (car c) (car p))
  68.              (- (cadr a) (cadr c))
  69.            )
  70.            (*
  71.              (- (car a) (car c))
  72.              (- (cadr c) (cadr p))
  73.            )
  74.          )
  75.     )
  76.     (or
  77.       (and
  78.         (<= 0 ab)
  79.         (<= 0 bc)
  80.         (<= 0 ca)
  81.       )
  82.       (and
  83.         (>= 0 ab)
  84.         (>= 0 bc)
  85.         (>= 0 ca)
  86.       )
  87.     )
  88.   )
  89.  
  90.  
  91.  
  92. ;;;;;;;;;;;;;;;;;;;;;;;;
  93.   (setq new_list
  94.          (if (not (vk_IsBoundaryClockwise point_list))
  95.            (reverse point_list)
  96.            point_list
  97.          )
  98.   )
  99.  
  100.   (while (> (length new_list) 3)
  101.     (setq a (car new_list)
  102.           b (cadr new_list)
  103.           c (caddr new_list)
  104.     )
  105.     (if (and
  106.           (is_left a b c)
  107.           (not (vl-remove-if
  108.                  'not
  109.                  (mapcar (function (lambda (x)
  110.                                      (is_point_Inside a b c x)
  111.                                    )
  112.                          )
  113.                          (vl-remove-if
  114.                            (function (lambda (y)
  115.                                        (or
  116.                                          (eq y a)
  117.                                          (eq y b)
  118.                                          (eq y c)
  119.                                        )
  120.                                      )
  121.                            )
  122.                            point_list
  123.                          )
  124.                  )
  125.                )
  126.           )
  127.         )
  128.       (setq
  129.         triangl_lict (cons (list a b c) triangl_lict)
  130.         new_list     (cons a (cddr new_list))
  131.       )
  132.       (setq
  133.         new_list
  134.          (append (cdr new_list) (list (car new_list)))
  135.       )
  136.     )
  137.   )
  138.   (cons new_list triangl_lict)
  139. )
может не очень изящная, но вроде работает  ;D
Команда для теста на замкнутых полилиниях
Код - Auto/Visual Lisp [Выбрать]
  1. (defun c:test_triangl (/ Make_LWPoly)
  2. ;;;
  3. ;;;Make_LWPoly
  4. ;;;
  5.   (defun Make_LWPoly (lst)
  6.     (entmakex
  7.       (append
  8.         (list (cons 0 "LWPOLYLINE")
  9.               (cons 100 "AcDbEntity")
  10.               (cons 100 "AcDbPolyline")
  11.               (cons 90 (length lst))
  12.               (cons 70 1)
  13.               (cons 210 '(0 0 1))
  14.         )
  15.         (mapcar (function (lambda (p) (cons 10 p)))
  16.                 lst
  17.         )
  18.       )
  19.     )
  20.   )
  21. ;;;;;;;;;;;
  22.   (mapcar 'Make_LWPoly
  23.           (polygon_triangul
  24.             (mapcar 'cdr
  25.                     (vl-remove-if-not
  26.                       '(lambda (x) (= (car x) 10))
  27.                       (entget (car (entsel)))
  28.                     )
  29.             )
  30.           )
  31.   )
  32. )