Triangulation de delaunay

Version imprimable

Bonjour messieurs,

Dans un elan de bonte, et ayant besoin du bouzin, je me suis mis a traduire cet algorithme de triangulation de delaunay du lua vers le C++.

Voici le bout de code en question sur lequel je bloque.

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
 
--- Triangulates a set of given vertices
-- @param ... a `vargarg` list of objects of type `Point`
-- @return a set of objects of type `Triangle`
-- @usage
-- local Delaunay = require 'Delaunay'
-- local Point    = Delaunay.Point 
-- local p1, p2, p3, p4 = Point(), Point(2,0), Point(1,1), Point(1,-1)
-- local triangles = Delaunay.triangulate(p1, p2, p3, p4)
-- for i = 1, #triangles do
--   print(triangles[i])
-- end
--
function Delaunay.triangulate(...)
  local vertices = {...}
  local nvertices = #vertices
  assert(nvertices > 2, "Cannot triangulate, needs more than 3 vertices")
  if nvertices == 3 then
    return {Triangle(unpack(vertices))}
  end
 
  local trmax = nvertices * 4
 
  local minX, minY = vertices[1].x, vertices[1].y
  local maxX, maxY = minX, minY
 
  for i = 1, #vertices do
    local vertex = vertices[i]
    vertex.id = i
    if vertex.x < minX then minX = vertex.x end
    if vertex.y < minY then minY = vertex.y end
    if vertex.x > maxX then maxX = vertex.x end
    if vertex.y > maxY then maxY = vertex.y end
  end
 
  local dx, dy = maxX - minX, maxY - minY
  local deltaMax = max(dx, dy)
  local midx, midy = (minX + maxX) * 0.5, (minY + maxY) * 0.5
 
  local p1 = Point(midx - 2 * deltaMax, midy - deltaMax)
  local p2 = Point(midx, midy + 2 * deltaMax)
  local p3 = Point(midx + 2 * deltaMax, midy - deltaMax)
  p1.id, p2.id, p3.id = nvertices + 1, nvertices + 2, nvertices + 3
  vertices[p1.id] = p1
  vertices[p2.id] = p2
  vertices[p3.id] = p3
 
  local triangles = {}
  triangles[#triangles + 1] = Triangle(vertices[nvertices + 1],
                                       vertices[nvertices + 2],
                                       vertices[nvertices + 3]
                              )
 
  for i = 1, nvertices do
 
    local edges = {}
    local ntriangles = #triangles
 
    for j = #triangles, 1, -1 do
      local curTriangle = triangles[j]
      if curTriangle:inCircumCircle(vertices[i]) then
        edges[#edges + 1] = curTriangle.e1
        edges[#edges + 1] = curTriangle.e2
        edges[#edges + 1] = curTriangle.e3
        remove(triangles, j)
      end
    end
 
    for j = #edges - 1, 1, -1 do
      for k = #edges, j + 1, -1 do
        if edges[j] and edges[k] and edges[j]:same(edges[k]) then
          remove(edges, j)
          remove(edges, k-1)
        end
      end
    end
 
    for j = 1, #edges do
      local n = #triangles
      assert(n <= trmax, "Generated more than needed triangles")
      triangles[n + 1] = Triangle(edges[j].p1, edges[j].p2, vertices[i])
    end
 
  end
 
  for i = #triangles, 1, -1 do
    local triangle = triangles[i]
    if (triangle.p1.id > nvertices or 
        triangle.p2.id > nvertices or 
        triangle.p3.id > nvertices) then
      remove(triangles, i)
    end
  end
 
  for _ = 1,3 do remove(vertices) end
 
  return triangles
 
end
 
return Delaunay

Et voici mon code en C++

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
 
std::vector<Triangle> Delaunay::triangulate(std::vector<Vec2f> points)
{
	std::size_t npoints = points.size();
	int trmax = npoints * 4;
 
	float minX = points[0].getX();
	float minY = points[0].getY();
	float maxX = minX;
	float maxY = minY;
 
	for(std::size_t i = 0; i < points.size(); ++i) {
		points[i].id = i;
		if (points[i].getX() < minX) 
			minX = points[i].getX();
    	if (points[i].getY() < minY)
			minY = points[i].getY();
    	if (points[i].getX() > maxX)
			maxX = points[i].getX();
    	if (points[i].getY() > maxY)
			maxY = points[i].getY();
	}
 
	float dx = maxX - minX;
	float dy = maxY - minY;
  	float deltaMax = std::max(dx, dy);
	float midx = (minX + maxX) * 0.5;
	float midy = (minY + maxY) * 0.5;
 
 
	Vec2f p1(midx - 2 * deltaMax, midy - deltaMax);
  	Vec2f p2(midx, midy + 2 * deltaMax);
  	Vec2f p3(midx + 2 * deltaMax, midy - deltaMax);
 
	p1.id = npoints + 1;
	p2.id = npoints + 2;
	p3.id = npoints + 3;
 
	points.push_back(p1);
	points.push_back(p2);
	points.push_back(p3);
 
	std::vector<Triangle> triangles;
	triangles.push_back(Triangle(points[points.size() - 1], points[points.size() - 2], points[points.size() - 3]));
 
	for(int i = 0; i < static_cast<int>(npoints); ++i) {
		std::vector<Edge> edges;
 
		for(int j = static_cast<int>(triangles.size()); j-- > 0; ) {
			Triangle curTriangle = triangles[j];
 
			if(curTriangle.inCircumCircle(points[i])) {
				edges.push_back(curTriangle.getE1());
				edges.push_back(curTriangle.getE2());
				edges.push_back(curTriangle.getE3());
				triangles.erase(triangles.begin() + j);
			}
		}
 
		for(int j = static_cast<int>(edges.size()) - 1; --j > 0; ) {
			for(int k = edges.size(); --k > j + 1; ) {
				if(edges[j].equals(edges[k])) {
					edges.erase(edges.begin() + j);	
					edges.erase(edges.begin() + k - 1);
					--k;
				}
			}	
		}
 
		for(int j = 0; j < static_cast<int>(edges.size()); ++j) {
			int n = triangles.size();
			std::cout << n << " " << trmax << std::endl;
			assert(n <= trmax && "Generated more than needed triangles");	
			triangles.push_back(Triangle(edges[j].getP1(), edges[j].getP2(), points[i]));
		}
	}
 
	for(int i = static_cast<int>(triangles.size()); --i > 0; ) {
		Triangle triangle = triangles[i];
		if(triangle.getP1().id > static_cast<int>(npoints) ||
			triangle.getP2().id > static_cast<int>(npoints) ||
			triangle.getP3().id > static_cast<int>(npoints)) {
			triangles.erase(triangles.begin() + i);
		}
	}
 
	return triangles;
}

Et voici ce qui cloche

Code:

1
2
3
 
pudding: delaunay.cpp:80: static std::vector<Triangle> Delaunay::triangulate(std::vector<Vec2f>): Assertion `n <= trmax && "Generated more than needed triangles"' failed.
[1]    11925 abort (core dumped)  ./pudding

D'apres le code lua, je suis sense generer 12 triangles pour un set de 10 points. Sauf que mon code C++ (en virant le assert), monte jusqu'a 60 n, et genere 1 miserable triangle.
Et j'ai aucune idee de pourquoi ni comment ca cloche.

Tous les fichier sont dans ce beau zip. Je me replonge dans mon code pour truver ce qui cloche.

Merci d'avance!

08/10/2015, 07h05
darkman19320

Salut,

Je pense qu'il y a déjà un problème quand tu parcours ta liste de triangles (c'est peut-être un mauvais copier-coller...):
Tu as écrit :

Code:

for(int j = static_cast<int>(triangles.size()); j-- > 0; )

Et il faudrait plutot:

Code:

for(int j = static_cast<int>(triangles.size() - 1); j-- > 0; )

Voir:

Code:

for(size_t j = triangle.size() - 1; j >= 0; --j)

Pour le reste, je n'ai pas encore regardé
08/10/2015, 07h34
Ehonn

Citation:

Envoyé par darkman19320

Code:

for(size_t j = triangle.size() - 1; j >= 0; --j)

j est toujours >= à 0 car il est non signé (size_t). Il y a une boucle infinie. Dans le même sens, il faudrait s'assurer que triangle.size() - 1 soit valide (triangle.size() doit être >= à 1).

Une solution pour régler ce souci : std::for_each et un reverse iterator.
08/10/2015, 08h17
darkman19320

Bien vu, j'ai le cerveau encore un peu endormi ce matin!

Tu penses qu'avec le std::for_each on peut directement supprimer un élément d'une liste?
08/10/2015, 08h24
Iradrille

Citation:

Envoyé par Bl4ckb0ne

Dans un elan de bonte, et ayant besoin du bouzin

Çà ne répond pas directement à ta question, mais regarde du coté de boost.
Passer d'un diagramme de Voronoy à une triangulation de Delaunay est simple et se fait en O(n) (il me semble).

@darkman19320, non avec un std::for_each / range-based loop tu ne peux pas supprimer des éléments (ni faire d'opérations qui invalide les itérateurs).
08/10/2015, 08h30
Ehonn

Oui on pourrait, mais il faudrait faire très attention à la possible invalidation des itérateurs lors de la suppression. (Edit: @Iradrille : on peut déplacer les éléments à supprimer à la fin et faire un .resize() après le std::for_each).
Le plus simple reste d'utiliser l'idiome .erase(...) avec std::remove(...)) ou std::remove_if(...))

L'autre solution (que je n'aime pas trop) pour continuer à utiliser les indices non signés avec des boucles inverses est de la décaler de 1 :

Code:

for (size_t i = v.size(); i >= 1; --i) { v[i-1]; }

Tout marche bien avec le fix de darkman. Merci pour le std::for_each et le reverse iterator, je connaissais pas du tout mais ca a l'air epic.

Maintenant resoudre ma couille algorithmique, voyes plutot :

- Execution du code c++

Code:

1
2
3
4
5
 
Triangl:
 Point x: 22.5 y: -4.5
 Point x: 4.5 y: 22.5
 Point x: -13.5 y: -4.5

- Execution du code lua

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
 
Triangle: 
  Point (2) x: 1.00 y: 3.00
  Point (3) x: 2.00 y: 5.00
  Point (4) x: 3.00 y: 4.00
Triangle: 
  Point (4) x: 3.00 y: 4.00
  Point (3) x: 2.00 y: 5.00
  Point (5) x: 9.00 y: 5.00
Triangle: 
  Point (1) x: 0.00 y: 0.00
  Point (6) x: 0.00 y: 2.00
  Point (7) x: 1.00 y: 0.00
Triangle: 
  Point (6) x: 0.00 y: 2.00
  Point (2) x: 1.00 y: 3.00
  Point (7) x: 1.00 y: 0.00
Triangle: 
  Point (2) x: 1.00 y: 3.00
  Point (4) x: 3.00 y: 4.00
  Point (8) x: 5.00 y: 3.00
Triangle: 
  Point (7) x: 1.00 y: 0.00
  Point (2) x: 1.00 y: 3.00
  Point (8) x: 5.00 y: 3.00
Triangle: 
  Point (4) x: 3.00 y: 4.00
  Point (5) x: 9.00 y: 5.00
  Point (8) x: 5.00 y: 3.00
Triangle: 
  Point (3) x: 2.00 y: 5.00
  Point (2) x: 1.00 y: 3.00
  Point (9) x: 1.00 y: 9.00
Triangle: 
  Point (2) x: 1.00 y: 3.00
  Point (6) x: 0.00 y: 2.00
  Point (9) x: 1.00 y: 9.00
Triangle: 
  Point (5) x: 9.00 y: 5.00
  Point (3) x: 2.00 y: 5.00
  Point (9) x: 1.00 y: 9.00
Triangle: 
  Point (7) x: 1.00 y: 0.00
  Point (8) x: 5.00 y: 3.00
  Point (10) x: 7.00 y: 3.00
Triangle: 
  Point (8) x: 5.00 y: 3.00
  Point (5) x: 9.00 y: 5.00
  Point (10) x: 7.00 y: 3.00