1 -- Bowyer and Watson algorithm
3 function BowyerWatson (listPoints,bbox)
4 local triangulation = {}
5 local lgth = #listPoints
6 -- add four points to listPoints to have a bounding box
7 listPoints = buildBoundingBox(listPoints)
9 triangulation[1] = {lgth+1, lgth+2, lgth+3}
10 -- the second triangle
11 triangulation[2] = {lgth+1, lgth+3, lgth+4}
12 -- add points one by one
14 -- find the triangles which the circumcircle contained the point to add
15 badTriangles = buildBadTriangles(listPoints[i],triangulation)
16 -- build the polygon of the cavity containing the point to add
17 polygon = buildCavity(badTriangles, triangulation)
18 -- remove the bad triangles
19 for j=1,#badTriangles do
20 table.remove(triangulation,badTriangles[j]-(j-1))
22 -- build the new triangles and add them to triangulation
24 table.insert(triangulation,{polygon[j][1],polygon[j][2],i})
26 end -- end adding points of the listPoints
27 -- remove bounding box
28 if(bbox ~= "bbox") then
29 triangulation = removeBoundingBox(triangulation,lgth)
30 table.remove(listPoints,lgth+1)
31 table.remove(listPoints,lgth+1)
32 table.remove(listPoints,lgth+1)
33 table.remove(listPoints,lgth+1)
39 function buildBoundingBox(listPoints)
40 -- listPoints : list of points
41 -- epsV : parameter for the distance of the bounding box
42 local xmin, xmax, ymin, ymax, eps
47 for i=1,#listPoints do
48 if (listPoints[i].x < xmin) then
49 xmin = listPoints[i].x
51 if (listPoints[i].x > xmax) then
52 xmax = listPoints[i].x
54 if (listPoints[i].y < ymin) then
55 ymin = listPoints[i].y
57 if (listPoints[i].y > ymax) then
58 ymax = listPoints[i].y
61 eps = math.max(math.abs(xmax-xmin),math.abs(ymax-ymin))*0.15
66 -- add points of the bounding box in last positions
67 table.insert(listPoints,{x=xmin,y=ymin})
68 table.insert(listPoints,{x=xmin,y=ymax})
69 table.insert(listPoints,{x=xmax,y=ymax})
70 table.insert(listPoints,{x=xmax,y=ymin})
74 function removeBoundingBox(triangulation,lgth)
75 -- build the four bounding box edge
82 for i=1,#triangulation do
83 boolE1 = pointInTriangle(point1,triangulation[i])
84 boolE2 = pointInTriangle(point2,triangulation[i])
85 boolE3 = pointInTriangle(point3,triangulation[i])
86 boolE4 = pointInTriangle(point4,triangulation[i])
87 if((not boolE1) and (not boolE2) and (not boolE3) and (not boolE4)) then
88 table.insert(newTriangulation,triangulation[i])
91 return newTriangulation
95 function buildBadTriangles(point, triangulation)
97 for j=1,#triangulation do -- for all triangles
98 A = listPoints[triangulation[j][1]]
99 B = listPoints[triangulation[j][2]]
100 C = listPoints[triangulation[j][3]]
101 center, radius = circoncircle(A,B,C)
102 CP = Vector(center,point)
103 if(VectorNorm(CP)<radius) then -- the point belongs to the circoncirle
104 table.insert(badTriangles,j)
110 -- construction of the cavity composed by the bad triangles around the point to add
111 function buildCavity(badTriangles, triangulation)
113 for j=1,#badTriangles do -- for all bad triangles
114 ind = badTriangles[j]
115 for k=1,3 do -- for all edges
116 edge = {triangulation[ind][k],triangulation[ind][k%3+1]}
118 for l = 1,#badTriangles do -- for all badtriangles
119 badInd = badTriangles[l]
120 if(badInd ~= ind) then -- if not the current one
121 edgeBord = edgeBord or edgeInTriangle(edge,triangulation[badInd])
124 -- if the edge does not belong to another bad triangle
125 if(edgeBord == false) then
126 -- insert the edge to the cavity
127 table.insert(polygon,edge)
134 function edgeInTriangle(e,t)
149 function pointInTriangle(e,t)
161 local out = {x = B.x - A.x, y = B.y - A.y}
165 function VectorNorm(NP)
166 return math.sqrt(NP.x*NP.x +NP.y*NP.y)
170 function circoncircle(M, N, P)
171 -- Compute center and radius of the circoncircle of the triangle M N P
173 -- return : (center [Point],radius [float])
175 local MN = Vector(M,N)
176 local NP = Vector(N,P)
177 local PM = Vector(P,M)
178 m = VectorNorm(NP) -- |NP|
179 n = VectorNorm(PM) -- |PM|
180 p = VectorNorm(MN) -- |MN|
182 d = (m + n + p) * (-m + n + p) * (m - n + p) * (m + n - p)
184 rad = m * n * p / math.sqrt(d)
188 d = -2 * (M.x * NP.y + N.x * PM.y + P.x * MN.y)
193 om2 = math.pow(VectorNorm(OM),2) -- |OM|**2
194 on2 = math.pow(VectorNorm(ON),2) -- |ON|**2
195 op2 = math.pow(VectorNorm(OP),2) -- |OP|**2
196 x0 = -(om2 * NP.y + on2 * PM.y + op2 * MN.y) / d
197 y0 = (om2 * NP.x + on2 * PM.x + op2 * MN.x) / d
203 return Out, rad -- (center [Point], R [float])
207 -------------------------- TeX
208 -- build the list of points
209 function buildList(chaine, mode)
210 -- if mode = int : the list is given in the chaine string (x1,y1);(x2,y2);...;(xn,yn)
211 -- if mode = ext : the list is given in a file line by line with space separation
213 if mode == "int" then
214 local points = string.explode(chaine, ";")
217 Sx,Sy=string.match(points[i],"%((.+),(.+)%)")
218 listPoints[i]={x=tonumber(Sx),y=tonumber(Sy)}
220 elseif mode == "ext" then
221 io.input(chaine) -- open the file
223 lines=string.explode(text,"\n+") -- all the lines
226 xy=string.explode(lines[i]," +")
227 listPoints[i]={x=tonumber(xy[1]),y=tonumber(xy[2])}
230 print("Non existing mode")
237 function rectangleList(a,b,nbrA,nbrB)
244 listPoints[k] = {x = (i-1)*stepA, y=(j-1)*stepB}
251 -- trace a triangulation with TikZ
252 function traceMeshTikZ(listPoints, triangulation,points,color)
254 for i=1,#triangulation do
255 PointI = listPoints[triangulation[i][1]]
256 PointJ = listPoints[triangulation[i][2]]
257 PointK = listPoints[triangulation[i][3]]
258 output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
260 if(points=="points") then
261 for i=1,#listPoints do
262 output = output .. "\\draw[color=".."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
269 -- trace a triangulation with MP
270 function traceMeshMP(listPoints, triangulation,points,color)
272 output = output .. " pair MeshPoints[];"
273 for i=1,#listPoints do
274 output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
277 for i=1,#triangulation do
278 PointI = listPoints[triangulation[i][1]]
279 PointJ = listPoints[triangulation[i][2]]
280 PointK = listPoints[triangulation[i][3]]
281 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\mpcolor{"..color.. "};"
283 if(points=="points") then
284 for i=1,#listPoints do
285 output = output .. "dotlabel.llft( btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u) withcolor \\mpcolor{"..color.. "};"
293 -- buildMesh with TikZ
294 function buildMeshTikZ(chaine,mode,points,bbox,full,scale,color)
295 listPoints = buildList(chaine, mode)
296 triangulation = BowyerWatson(listPoints,bbox)
297 output = traceMeshTikZ(listPoints, triangulation,points,color)
298 if(full=="full") then
299 output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" .. output .."\\end{tikzpicture}"
306 function buildMeshMP(chaine,mode,points,bbox,full,scale,color)
307 listPoints = buildList(chaine, mode)
308 triangulation = BowyerWatson(listPoints,bbox)
309 output = traceMeshMP(listPoints, triangulation,points,color)
310 if(full=="full") then
311 output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale.. ";" .. output .."endfig;\\end{mplibcode}"
313 output = "u:="..scale.. ";".. output
319 function buildRect(largeur,a,b,nbrA, nbrB)
320 listPoints = rectangleList(a,b,nbrA,nbrB)
321 triangulation = BowyerWatson(listPoints,"none")
322 traceTikZ(listPoints, triangulation,largeur,"none")
326 -- function give a real polygon without repeting points
327 function cleanPoly(polygon)
331 table.insert(polyNew, e1)
332 table.insert(polyNew, e2)
335 bool1 = (polygon[i][1] == polyNew[j])
336 bool2 = (polygon[i][2] == polyNew[j])
337 if(bool1 or bool2) then -- the edge has a common point with polyNew[j]
339 table.insert(polyNew, polygon[i][1])
341 elseif(not bool2) then
342 table.insert(polyNew, polygon[i][2])
351 function TeXaddOnePointTikZ(chaine,point,step,color,colorBack, colorNew, colorCircle)
352 Sx,Sy=string.match(point,"%((.+),(.+)%)")
355 listPoints = buildList(chaine, "int")
356 -- build the triangulation
357 triangulation = BowyerWatson(listPoints,"none")
358 badTriangles = buildBadTriangles(P,triangulation)
359 if(step == "badT") then
361 for i=1,#triangulation do
362 PointI = listPoints[triangulation[i][1]]
363 PointJ = listPoints[triangulation[i][2]]
364 PointK = listPoints[triangulation[i][3]]
365 output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
367 -- draw and fill the bad triangle
368 for i=1,#badTriangles do
369 PointI = listPoints[triangulation[badTriangles[i]][1]]
370 PointJ = listPoints[triangulation[badTriangles[i]][2]]
371 PointK = listPoints[triangulation[badTriangles[i]][3]]
372 output = output .. "\\draw[fill="..colorBack.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
374 -- draw the circoncircle
375 for i=1,#badTriangles do
376 PointI = listPoints[triangulation[badTriangles[i]][1]]
377 PointJ = listPoints[triangulation[badTriangles[i]][2]]
378 PointK = listPoints[triangulation[badTriangles[i]][3]]
379 center, radius = circoncircle(PointI, PointJ, PointK)
380 output = output .. "\\draw[dashed, color="..colorCircle.."] ("..center.x .. "," .. center.y .. ") circle ("..radius ..");"
383 for i=1,#listPoints do
384 output = output .. "\\draw[color ="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
386 -- mark the point to add
387 output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
388 elseif(step == "cavity") then
389 polygon = buildCavity(badTriangles, triangulation)
390 polyNew = cleanPoly(polygon)
391 -- remove the bad triangles
392 for j=1,#badTriangles do
393 table.remove(triangulation,badTriangles[j]-(j-1))
395 -- draw the triangles
396 for i=1,#triangulation do
397 PointI = listPoints[triangulation[i][1]]
398 PointJ = listPoints[triangulation[i][2]]
399 PointK = listPoints[triangulation[i][3]]
400 output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
402 -- fill and draw the cavity
405 PointI = listPoints[polyNew[i]]
406 path = path .. "(".. PointI.x ..",".. PointI.y ..")--"
408 output = output .. "\\draw[color="..colorNew..",fill ="..colorBack..", thick] " .. path .. "cycle;"
409 -- mark the points of the mesh
410 for i=1,#listPoints do
411 output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
413 -- mark the adding point
414 output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
415 elseif(step == "newT") then
416 polygon = buildCavity(badTriangles, triangulation)
417 polyNew = cleanPoly(polygon)
418 -- remove the bad triangles
419 for j=1,#badTriangles do
420 table.remove(triangulation,badTriangles[j]-(j-1))
422 -- draw the triangle of the triangulation
423 for i=1,#triangulation do
424 PointI = listPoints[triangulation[i][1]]
425 PointJ = listPoints[triangulation[i][2]]
426 PointK = listPoints[triangulation[i][3]]
427 output = output .. "\\draw[color ="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
429 -- fill and draw the cavity
432 PointI = listPoints[polyNew[i]]
433 path = path .. "(".. PointI.x ..",".. PointI.y ..")--"
435 output = output .. "\\draw[color="..colorNew..",fill ="..colorBack..", thick] " .. path .. "cycle;"
436 -- draw the new triangles composed by the edges of the polygon and the added point
438 output = output .. "\\draw[color=TeXCluaMeshNewTikZ, thick]".."(".. listPoints[polygon[i][1]].x .. "," .. listPoints[polygon[i][1]].y .. ") -- (" .. listPoints[polygon[i][2]].x .. "," .. listPoints[polygon[i][2]].y ..");"
439 output = output .. "\\draw[color="..colorNew..", thick]".."(".. listPoints[polygon[i][1]].x .. "," .. listPoints[polygon[i][1]].y .. ") -- (" .. P.x .. "," .. P.y ..");"
440 output = output .. "\\draw[color="..colorNew..", thick]".."(".. listPoints[polygon[i][2]].x .. "," .. listPoints[polygon[i][2]].y .. ") -- (" .. P.x .. "," .. P.y ..");"
443 for i=1,#listPoints do
444 output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
446 -- mark the added point
447 output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
452 function TeXaddOnePointMP(listPoints,P,step,color,colorBack, colorNew, colorCircle)
454 output = output .. "pair MeshPoints[];"
455 -- build the triangulation
456 triangulation = BowyerWatson(listPoints,"none")
457 badTriangles = buildBadTriangles(P,triangulation)
458 for i=1,#listPoints do
459 output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
461 if(step == "badT") then
463 for i=1,#triangulation do
464 PointI = listPoints[triangulation[i][1]]
465 PointJ = listPoints[triangulation[i][2]]
466 PointK = listPoints[triangulation[i][3]]
467 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\mpcolor{" .. color .."};"
469 -- draw and fill the bad triangle
470 for i=1,#badTriangles do
471 PointI = listPoints[triangulation[badTriangles[i]][1]]
472 PointJ = listPoints[triangulation[badTriangles[i]][2]]
473 PointK = listPoints[triangulation[badTriangles[i]][3]]
474 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\mpcolor{" .. color .."};"
475 output = output .. "fill (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\mpcolor{" .. colorBack .."};"
477 -- draw the circoncircle
478 for i=1,#badTriangles do
479 PointI = listPoints[triangulation[badTriangles[i]][1]]
480 PointJ = listPoints[triangulation[badTriangles[i]][2]]
481 PointK = listPoints[triangulation[badTriangles[i]][3]]
482 center, radius = circoncircle(PointI, PointJ, PointK)
483 output = output .. "draw fullcircle scaled ("..radius .."*2u) shifted ("..center.x .. "*u," .. center.y .. "*u) dashed evenly withcolor \\mpcolor{" .. colorCircle .."};"
486 for i=1,#listPoints do
487 output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\mpcolor{" .. color .."};"
489 -- mark the point to add
490 output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\mpcolor{" .. colorNew .."};"
491 elseif(step == "cavity") then
492 polygon = buildCavity(badTriangles, triangulation)
493 polyNew = cleanPoly(polygon)
494 -- remove the bad triangles
495 for j=1,#badTriangles do
496 table.remove(triangulation,badTriangles[j]-(j-1))
498 -- draw the triangles
499 for i=1,#triangulation do
500 PointI = listPoints[triangulation[i][1]]
501 PointJ = listPoints[triangulation[i][2]]
502 PointK = listPoints[triangulation[i][3]]
503 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\mpcolor{" .. color .."};"
505 -- fill and draw the cavity
508 PointI = listPoints[polyNew[i]]
509 path = path .. "(".. PointI.x ..",".. PointI.y ..")*u--"
511 output = output .. "fill " .. path .. "cycle withcolor \\mpcolor{" .. colorBack .."};"
512 output = output .. "draw " .. path .. "cycle withcolor \\mpcolor{" .. colorNew .."} withpen pencircle scaled 1pt;"
513 -- mark the points of the mesh
514 for i=1,#listPoints do
515 output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\mpcolor{" .. color .."};"
517 -- mark the adding point
518 output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\mpcolor{" .. colorNew .."};"
519 elseif(step == "newT") then
520 polygon = buildCavity(badTriangles, triangulation)
521 polyNew = cleanPoly(polygon)
522 -- remove the bad triangles
523 for j=1,#badTriangles do
524 table.remove(triangulation,badTriangles[j]-(j-1))
526 -- draw the triangle of the triangulation
527 for i=1,#triangulation do
528 PointI = listPoints[triangulation[i][1]]
529 PointJ = listPoints[triangulation[i][2]]
530 PointK = listPoints[triangulation[i][3]]
531 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\mpcolor{" .. color .."};"
536 PointI = listPoints[polyNew[i]]
537 path = path .. "(".. PointI.x ..",".. PointI.y ..")*u--"
539 output = output .. "fill " .. path .. "cycle withcolor \\mpcolor{" .. colorBack .."};"
540 -- draw the new triangles composed by the edges of the polygon and the added point
542 output = output .. "draw".."(".. listPoints[polygon[i][1]].x .. "," .. listPoints[polygon[i][1]].y .. ")*u -- (" .. listPoints[polygon[i][2]].x .. "," .. listPoints[polygon[i][2]].y ..")*u withcolor \\mpcolor{" .. colorNew .."} withpen pencircle scaled 1pt;"
543 output = output .. "draw".."(".. listPoints[polygon[i][1]].x .. "," .. listPoints[polygon[i][1]].y .. ")*u -- (" .. P.x .. "," .. P.y ..")*u withcolor \\mpcolor{" .. colorNew .."} withpen pencircle scaled 1pt;"
544 output = output .. "draw".."(".. listPoints[polygon[i][2]].x .. "," .. listPoints[polygon[i][2]].y .. ")*u -- (" .. P.x .. "," .. P.y ..")*u withcolor \\mpcolor{" .. colorNew .."} withpen pencircle scaled 1pt;"
547 for i=1,#listPoints do
548 output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\mpcolor{" .. color .."};"
550 -- mark the added point
551 output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\mpcolor{ " .. colorNew .."};"
557 -- build the list of points extern and stop at nbr
558 function buildListExt(chaine, stop)
560 io.input(chaine) -- open the file
562 lines=string.explode(text,"\n+") -- all the lines
563 for i=1,tonumber(stop) do
564 xy=string.explode(lines[i]," +")
565 table.insert(listPoints,{x=tonumber(xy[1]),y=tonumber(xy[2])})
567 xy=string.explode(lines[stop+1]," +")
568 point={x=tonumber(xy[1]),y=tonumber(xy[2])}
569 return point, listPoints
573 function TeXOnePointTikZ(chaine,point,step,color,colorBack,colorNew,colorCircle,scale)
574 output = TeXaddOnePointTikZ(chaine,point,step,color,colorBack,colorNew,colorCircle)
575 output = "\\noindent\\begin{tikzpicture}[x="..scale..",y="..scale.."]".. output .. "\\end{tikzpicture}"
579 function TeXOnePointMP(chaine,point,step,color,colorBack,colorNew,colorCircle,scale,mode,picture)
581 Sx,Sy=string.match(point,"%((.+),(.+)%)")
583 listPoints = buildList(chaine, mode)
586 P, listPoints = buildListExt(chaine,tonumber(point))
588 output = TeXaddOnePointMP(listPoints,P,step,color,colorBack,colorNew,colorCircle)
589 if(picture=="full") then
590 output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale..";".. output .. "endfig;\\end{mplibcode}"
592 output = "u:="..scale..";".. output