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,type="bbox"}
10 -- the second triangle
11 triangulation[2] = {lgth+1, lgth+3, lgth+4,type="bbox"}
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 if((polygon[j][1]>lgth) or (polygon[j][2]>lgth) or (i>lgth)) then
25 table.insert(triangulation,{polygon[j][1],polygon[j][2],i,type="bbox"})
27 table.insert(triangulation,{polygon[j][1],polygon[j][2],i,type="in"})
30 end -- end adding points of the listPoints
31 -- remove bounding box
32 if(bbox ~= "bbox") then
33 triangulation = removeBoundingBox(triangulation,lgth)
34 table.remove(listPoints,lgth+1)
35 table.remove(listPoints,lgth+1)
36 table.remove(listPoints,lgth+1)
37 table.remove(listPoints,lgth+1)
43 function buildBoundingBox(listPoints)
44 -- listPoints : list of points
45 -- epsV : parameter for the distance of the bounding box
46 local xmin, xmax, ymin, ymax, eps
51 for i=1,#listPoints do
52 if (listPoints[i].x < xmin) then
53 xmin = listPoints[i].x
55 if (listPoints[i].x > xmax) then
56 xmax = listPoints[i].x
58 if (listPoints[i].y < ymin) then
59 ymin = listPoints[i].y
61 if (listPoints[i].y > ymax) then
62 ymax = listPoints[i].y
65 eps = math.max(math.abs(xmax-xmin),math.abs(ymax-ymin))*0.15
70 -- add points of the bounding box in last positions
71 table.insert(listPoints,{x=xmin,y=ymin,type="bbox"})
72 table.insert(listPoints,{x=xmin,y=ymax,type="bbox"})
73 table.insert(listPoints,{x=xmax,y=ymax,type="bbox"})
74 table.insert(listPoints,{x=xmax,y=ymin,type="bbox"})
78 function removeBoundingBox(triangulation,lgth)
79 -- build the four bounding box edge
86 for i=1,#triangulation do
87 boolE1 = pointInTriangle(point1,triangulation[i])
88 boolE2 = pointInTriangle(point2,triangulation[i])
89 boolE3 = pointInTriangle(point3,triangulation[i])
90 boolE4 = pointInTriangle(point4,triangulation[i])
91 if((not boolE1) and (not boolE2) and (not boolE3) and (not boolE4)) then
92 table.insert(newTriangulation,triangulation[i])
95 return newTriangulation
99 function buildBadTriangles(point, triangulation)
101 for j=1,#triangulation do -- for all triangles
102 A = listPoints[triangulation[j][1]]
103 B = listPoints[triangulation[j][2]]
104 C = listPoints[triangulation[j][3]]
105 center, radius = circoncircle(A,B,C)
106 CP = Vector(center,point)
107 if(VectorNorm(CP)<radius) then -- the point belongs to the circoncirle
108 table.insert(badTriangles,j)
114 -- construction of the cavity composed by the bad triangles around the point to add
115 function buildCavity(badTriangles, triangulation)
117 for j=1,#badTriangles do -- for all bad triangles
118 ind = badTriangles[j]
119 for k=1,3 do -- for all edges
120 edge = {triangulation[ind][k],triangulation[ind][k%3+1]}
122 for l = 1,#badTriangles do -- for all badtriangles
123 badInd = badTriangles[l]
124 if(badInd ~= ind) then -- if not the current one
125 edgeBord = edgeBord or edgeInTriangle(edge,triangulation[badInd])
128 -- if the edge does not belong to another bad triangle
129 if(edgeBord == false) then
130 -- insert the edge to the cavity
131 table.insert(polygon,edge)
138 function edgeInTriangle(e,t)
153 function pointInTriangle(e,t)
165 local out = {x = B.x - A.x, y = B.y - A.y}
169 function VectorNorm(NP)
170 return math.sqrt(NP.x*NP.x +NP.y*NP.y)
174 function circoncircle(M, N, P)
175 -- Compute center and radius of the circoncircle of the triangle M N P
177 -- return : (center [Point],radius [float])
179 local MN = Vector(M,N)
180 local NP = Vector(N,P)
181 local PM = Vector(P,M)
182 m = VectorNorm(NP) -- |NP|
183 n = VectorNorm(PM) -- |PM|
184 p = VectorNorm(MN) -- |MN|
186 d = (m + n + p) * (-m + n + p) * (m - n + p) * (m + n - p)
188 rad = m * n * p / math.sqrt(d)
192 d = -2 * (M.x * NP.y + N.x * PM.y + P.x * MN.y)
197 om2 = math.pow(VectorNorm(OM),2) -- |OM|**2
198 on2 = math.pow(VectorNorm(ON),2) -- |ON|**2
199 op2 = math.pow(VectorNorm(OP),2) -- |OP|**2
200 x0 = -(om2 * NP.y + on2 * PM.y + op2 * MN.y) / d
201 y0 = (om2 * NP.x + on2 * PM.x + op2 * MN.x) / d
207 return Out, rad -- (center [Point], R [float])
210 -- compute the list of the circumcircle of a triangulation
211 function listCircumCenter(listPoints,triangulation)
213 for j=1,#triangulation do
214 A = listPoints[triangulation[j][1]]
215 B = listPoints[triangulation[j][2]]
216 C = listPoints[triangulation[j][3]]
217 center, radius = circoncircle(A,B,C)
218 table.insert(list,{x=center.x,y=center.y,r=radius})
223 -- find the three neighbour triangles of T
224 function findNeighbour(T,i,triangulation)
226 -- i : index of T in triangualation
230 -- define the three edge
234 for j=1,#triangulation do
236 if(edgeInTriangle(e1,triangulation[j])) then
239 if(edgeInTriangle(e2,triangulation[j])) then
242 if(edgeInTriangle(e3,triangulation[j])) then
250 -- test if edge are the same (reverse)
251 function equalEdge(e1,e2)
252 if(((e1[1] == e2[1]) and (e1[2] == e2[2])) or ((e1[1] == e2[2]) and (e1[2] == e2[1]))) then
259 -- test if the edge belongs to the list
260 function edgeInList(e,listE)
263 if(equalEdge(e,listE[i])) then
270 -- build the edges of the Voronoi diagram with a given triangulation
271 function buildVoronoi(listPoints, triangulation)
272 listCircumCircle = listCircumCenter(listPoints, triangulation)
274 for i=1,#listCircumCircle do
275 listN = findNeighbour(triangulation[i],i,triangulation)
278 if( not edgeInList(edge, listVoronoi)) then
279 table.insert(listVoronoi, edge)
286 -------------------------- TeX
287 -- build the list of points
288 function buildList(chaine, mode)
289 -- if mode = int : the list is given in the chaine string (x1,y1);(x2,y2);...;(xn,yn)
290 -- if mode = ext : the list is given in a file line by line with space separation
292 if mode == "int" then
293 local points = string.explode(chaine, ";")
296 Sx,Sy=string.match(points[i],"%((.+),(.+)%)")
297 listPoints[i]={x=tonumber(Sx),y=tonumber(Sy)}
299 elseif mode == "ext" then
300 io.input(chaine) -- open the file
302 lines=string.explode(text,"\n+") -- all the lines
305 xy=string.explode(lines[i]," +")
306 listPoints[i]={x=tonumber(xy[1]),y=tonumber(xy[2])}
309 print("Non existing mode")
316 function rectangleList(a,b,nbrA,nbrB)
323 listPoints[k] = {x = (i-1)*stepA, y=(j-1)*stepB}
331 -- trace Voronoi with MP
332 function traceVoronoiMP(listPoints, triangulation,listVoronoi, points, tri)
333 listCircumC = listCircumCenter(listPoints,triangulation)
335 output = output .. " pair MeshPoints[];"
336 for i=1,#listPoints do
337 output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ")*u;"
339 output = output .. " pair CircumCenters[];"
340 for i=1,#listCircumC do
341 output = output .. "CircumCenters[".. i .. "] = (" .. listCircumC[i].x .. "," .. listCircumC[i].y .. ")*u;"
344 for i=1,#triangulation do
345 PointI = listPoints[triangulation[i][1]]
346 PointJ = listPoints[triangulation[i][2]]
347 PointK = listPoints[triangulation[i][3]]
348 if(triangulation[i].type == "bbox") then
349 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBbox;"
351 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
355 for i=1,#listVoronoi do
356 PointI = listCircumC[listVoronoi[i][1]]
357 PointJ = listCircumC[listVoronoi[i][2]]
358 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u withcolor \\luameshmpcolorVoronoi;"
360 if(points=="points") then
362 for i=1,#listPoints do
363 if(listPoints[i].type == "bbox") then
364 output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{"..j.."}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
367 output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
370 for i=1,#listCircumC do
371 output = output .. "dotlabel.llft (btex $\\CircumPoint_{" .. i .. "}$ etex, (" .. listCircumC[i].x ..",".. listCircumC[i].y .. ")*u ) withcolor \\luameshmpcolorVoronoi ;"
375 for i=1,#listPoints do
376 if(listPoints[i].type == "bbox") then
377 output = output .. "drawdot (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u withcolor \\luameshmpcolorBbox withpen pencircle scaled 3;"
380 output = output .. "drawdot (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u withcolor \\luameshmpcolor withpen pencircle scaled 3;"
383 for i=1,#listCircumC do
384 output = output .. "drawdot (" .. listCircumC[i].x ..",".. listCircumC[i].y .. ")*u withcolor \\luameshmpcolorVoronoi withpen pencircle scaled 3;"
392 -- trace Voronoi with TikZ
393 function traceVoronoiTikZ(listPoints, triangulation,listVoronoi, points, tri,color,colorBbox,colorVoronoi)
394 listCircumC = listCircumCenter(listPoints,triangulation)
396 for i=1,#listPoints do
397 output = output .. "\\coordinate (MeshPoints".. i .. ") at (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
399 for i=1,#listCircumC do
400 output = output .. "\\coordinate (CircumPoints".. i .. ") at (" .. listCircumC[i].x .. "," .. listCircumC[i].y .. ");"
403 for i=1,#triangulation do
404 PointI = listPoints[triangulation[i][1]]
405 PointJ = listPoints[triangulation[i][2]]
406 PointK = listPoints[triangulation[i][3]]
407 if(triangulation[i].type == "bbox") then
408 output = output .. "\\draw[color="..colorBbox.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
410 output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
414 for i=1,#listVoronoi do
415 PointI = listCircumC[listVoronoi[i][1]]
416 PointJ = listCircumC[listVoronoi[i][2]]
417 output = output .. "\\draw[color="..colorVoronoi.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..");"
419 if(points=="points") then
421 for i=1,#listPoints do
422 if(listPoints[i].type == "bbox") then
423 if(listPoints[i].type == "bbox") then
424 output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint^*_{" .. j .. "}$};"
427 output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
431 for i=1,#listCircumC do
432 output = output .. "\\draw[color="..colorVoronoi.."] (" .. listCircumC[i].x ..",".. listCircumC[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\CircumPoint_{" .. i .. "}$};"
436 for i=1,#listPoints do
437 if(listPoints[i].type == "bbox") then
438 output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} ;"
441 output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} ;"
444 for i=1,#listCircumC do
445 output = output .. "\\draw[color="..colorVoronoi.."] (" .. listCircumC[i].x ..",".. listCircumC[i].y .. ") node {$\\bullet$};"
453 -- buildVoronoi with MP
454 function buildVoronoiMPBW(chaine,mode,points,bbox,scale,tri)
455 listPoints = buildList(chaine, mode)
456 triangulation = BowyerWatson(listPoints,bbox)
457 listVoronoi = buildVoronoi(listPoints, triangulation)
458 output = traceVoronoiMP(listPoints,triangulation,listVoronoi,points,tri)
459 output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale.. ";" .. output .."endfig;\\end{mplibcode}"
464 -- buildVoronoi with TikZ
465 function buildVoronoiTikZBW(chaine,mode,points,bbox,scale,tri,color,colorBbox,colorVoronoi)
466 listPoints = buildList(chaine, mode)
467 triangulation = BowyerWatson(listPoints,bbox)
468 listVoronoi = buildVoronoi(listPoints, triangulation)
469 output = traceVoronoiTikZ(listPoints,triangulation,listVoronoi,points,tri,color,colorBbox,colorVoronoi)
470 output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" .. output .."\\end{tikzpicture}" tex.sprint(output)
475 -- trace a triangulation with TikZ
476 function traceMeshTikZ(listPoints, triangulation,points,color,colorBbox)
478 for i=1,#listPoints do
479 output = output .. "\\coordinate (MeshPoints".. i .. ") at (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
481 for i=1,#triangulation do
482 PointI = listPoints[triangulation[i][1]]
483 PointJ = listPoints[triangulation[i][2]]
484 PointK = listPoints[triangulation[i][3]]
485 if(triangulation[i].type == "bbox") then
486 output = output .. "\\draw[color="..colorBbox.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
488 output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
491 if(points=="points") then
493 for i=1,#listPoints do
494 if(listPoints[i].type == "bbox") then
495 output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint^*_{" .. j .. "}$};"
498 output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
506 -- trace a triangulation with MP
507 function traceMeshMP(listPoints, triangulation,points)
509 output = output .. " pair MeshPoints[];"
510 for i=1,#listPoints do
511 output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ")*u;"
514 for i=1,#triangulation do
515 PointI = listPoints[triangulation[i][1]]
516 PointJ = listPoints[triangulation[i][2]]
517 PointK = listPoints[triangulation[i][3]]
518 if(triangulation[i].type == "bbox") then
519 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBbox;"
521 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
524 if(points=="points") then
526 for i=1,#listPoints do
527 if(listPoints[i].type == "bbox") then
528 output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{"..j.."}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
531 output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
540 function buildMeshMPBW(chaine,mode,points,bbox,scale)
541 listPoints = buildList(chaine, mode)
542 triangulation = BowyerWatson(listPoints,bbox)
543 output = traceMeshMP(listPoints, triangulation,points)
544 output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale.. ";" .. output .."endfig;\\end{mplibcode}"
548 -- buildMesh with MP include code
549 function buildMeshMPBWinc(chaine,beginning, ending,mode,points,bbox,scale)
550 listPoints = buildList(chaine, mode)
551 triangulation = BowyerWatson(listPoints,bbox)
552 output = traceMeshMP(listPoints, triangulation,points)
553 output = "\\leavevmode\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"
557 -- buildMesh with TikZ
558 function buildMeshTikZBW(chaine,mode,points,bbox,scale,color,colorBbox)
559 listPoints = buildList(chaine, mode)
560 triangulation = BowyerWatson(listPoints,bbox)
561 output = traceMeshTikZ(listPoints, triangulation,points,color,colorBbox)
562 output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" .. output .."\\end{tikzpicture}"
566 -- buildMesh with TikZ
567 function buildMeshTikZBWinc(chaine,beginning, ending,mode,points,bbox,scale,color,colorBbox)
568 listPoints = buildList(chaine, mode)
569 triangulation = BowyerWatson(listPoints,bbox)
570 output = traceMeshTikZ(listPoints, triangulation,points,color,colorBbox)
571 output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" ..beginning.. output..ending .."\\end{tikzpicture}"
576 -- print points of the mesh
577 function tracePointsMP(listPoints,points)
579 output = output .. " pair MeshPoints[];"
580 for i=1,#listPoints do
581 output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ")*u;"
583 if(points=="points") then
585 for i=1,#listPoints do
586 if(listPoints[i].type == "bbox") then
587 output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{" .. j .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
590 output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
594 for i=1,#listPoints do
595 if(listPoints[i].type == "bbox") then
596 output = output .. "drawdot (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u withcolor \\luameshmpcolorBbox withpen pencircle scaled 3;"
598 output = output .. "drawdot (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u withcolor \\luameshmpcolor withpen pencircle scaled 3;"
605 -- print points of the mesh
606 function tracePointsTikZ(listPoints,points,color,colorBbox)
608 for i=1,#listPoints do
609 output = output .. "\\coordinate (MeshPoints".. i .. ") at (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
611 if(points=="points") then
613 for i=1,#listPoints do
614 if(listPoints[i].type == "bbox") then
615 output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint^*_{" .. j .. "}$};"
618 output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
622 for i=1,#listPoints do
623 if(listPoints[i].type == "bbox") then
624 output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} ;"
626 output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} ;"
633 -- print points to mesh
634 function printPointsMP(chaine,mode,points,bbox,scale)
635 listPoints = buildList(chaine, mode)
636 if(bbox == "bbox" ) then
637 listPoints = buildBoundingBox(listPoints)
639 output = tracePointsMP(listPoints,points)
640 output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale.. ";" .. output .."endfig;\\end{mplibcode}"
645 -- print points to mesh
646 function printPointsMPinc(chaine,beginning, ending, mode,points,bbox,scale)
647 listPoints = buildList(chaine, mode)
648 if(bbox == "bbox" ) then
649 listPoints = buildBoundingBox(listPoints)
651 output = tracePointsMP(listPoints,points)
652 output = "\\leavevmode\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"
656 -- print points to mesh
657 function printPointsTikZ(chaine,mode,points,bbox,scale,color,colorBbox)
658 listPoints = buildList(chaine, mode)
659 if(bbox == "bbox" ) then
660 listPoints = buildBoundingBox(listPoints)
662 output = tracePointsTikZ(listPoints,points,color,colorBbox)
663 output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" .. output .."\\end{tikzpicture}"
668 -- print points to mesh
669 function printPointsTikZinc(chaine,beginning, ending, mode,points,bbox,scale,color,colorBbox)
670 listPoints = buildList(chaine, mode)
671 if(bbox == "bbox" ) then
672 listPoints = buildBoundingBox(listPoints)
674 output = tracePointsTikZ(listPoints,points,color,colorBbox)
675 output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" ..beginning.. output..ending .."\\end{tikzpicture}"
681 function buildRect(largeur,a,b,nbrA, nbrB)
682 listPoints = rectangleList(a,b,nbrA,nbrB)
683 triangulation = BowyerWatson(listPoints,"none")
684 traceTikZ(listPoints, triangulation,largeur,"none")
687 local function shallowCopy(original)
689 for key, value in pairs(original) do
695 -- function give a real polygon without repeting points
696 function cleanPoly(polygon)
698 polyCopy = shallowCopy(polygon)
701 table.insert(polyNew, e1)
702 table.insert(polyNew, e2)
703 table.remove(polyCopy,1)
708 while (i<=#polyCopy and find==false) do
709 bool1 = (polyCopy[i][1] == polyNew[j])
710 bool2 = (polyCopy[i][2] == polyNew[j])
711 if(bool1 or bool2) then -- the edge has a common point with polyNew[j]
713 table.insert(polyNew, polyCopy[i][1])
715 table.remove(polyCopy,i)
717 elseif(not bool2) then
718 table.insert(polyNew, polyCopy[i][2])
720 table.remove(polyCopy,i)
731 function TeXaddOnePointTikZ(listPoints,P,step,bbox,color,colorBack, colorNew, colorCircle,colorBbox)
733 -- build the triangulation
734 triangulation = BowyerWatson(listPoints,bbox)
735 badTriangles = buildBadTriangles(P,triangulation)
736 for i=1,#listPoints do
737 output = output .. "\\coordinate (MeshPoints".. i .. ") at (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
739 if(step == "badT") then
741 for i=1,#triangulation do
742 PointI = listPoints[triangulation[i][1]]
743 PointJ = listPoints[triangulation[i][2]]
744 PointK = listPoints[triangulation[i][3]]
745 if(triangulation[i].type == "bbox") then
746 output = output .. "\\draw[color="..colorBbox.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
748 output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
751 -- draw and fill the bad triangle
752 for i=1,#badTriangles do
753 PointI = listPoints[triangulation[badTriangles[i]][1]]
754 PointJ = listPoints[triangulation[badTriangles[i]][2]]
755 PointK = listPoints[triangulation[badTriangles[i]][3]]
756 output = output .. "\\draw[fill="..colorBack.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
758 -- draw the circoncircle
759 for i=1,#badTriangles do
760 PointI = listPoints[triangulation[badTriangles[i]][1]]
761 PointJ = listPoints[triangulation[badTriangles[i]][2]]
762 PointK = listPoints[triangulation[badTriangles[i]][3]]
763 center, radius = circoncircle(PointI, PointJ, PointK)
764 output = output .. "\\draw[dashed, color="..colorCircle.."] ("..center.x .. "," .. center.y .. ") circle ("..radius ..");"
768 for i=1,#listPoints do
769 if(listPoints[i].type == "bbox") then
770 output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint^*_{" .. j .. "}$};"
773 output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
776 -- mark the point to add
777 output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
778 elseif(step == "cavity") then
779 polygon = buildCavity(badTriangles, triangulation)
780 polyNew = cleanPoly(polygon)
781 -- remove the bad triangles
782 for j=1,#badTriangles do
783 table.remove(triangulation,badTriangles[j]-(j-1))
785 -- draw the triangles
786 for i=1,#triangulation do
787 PointI = listPoints[triangulation[i][1]]
788 PointJ = listPoints[triangulation[i][2]]
789 PointK = listPoints[triangulation[i][3]]
790 if(triangulation[i].type == "bbox") then
791 output = output .. "\\draw[color="..colorBbox.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
793 output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
796 -- fill and draw the cavity
799 PointI = listPoints[polyNew[i]]
800 path = path .. "(".. PointI.x ..",".. PointI.y ..")--"
802 output = output .. "\\draw[color="..colorNew..",fill ="..colorBack..", thick] " .. path .. "cycle;"
803 -- mark the points of the mesh
805 for i=1,#listPoints do
806 if(listPoints[i].type == "bbox") then
807 output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint^*_{" .. j .. "}$};"
810 output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
813 -- mark the adding point
814 output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
815 elseif(step == "newT") then
816 polygon = buildCavity(badTriangles, triangulation)
817 polyNew = cleanPoly(polygon)
818 -- remove the bad triangles
819 for j=1,#badTriangles do
820 table.remove(triangulation,badTriangles[j]-(j-1))
822 -- draw the triangle of the triangulation
823 for i=1,#triangulation do
824 PointI = listPoints[triangulation[i][1]]
825 PointJ = listPoints[triangulation[i][2]]
826 PointK = listPoints[triangulation[i][3]]
827 if(triangulation[i].type == "bbox") then
828 output = output .. "\\draw[color="..colorBbox.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
830 output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
833 -- fill and draw the cavity
836 PointI = listPoints[polyNew[i]]
837 path = path .. "(".. PointI.x ..",".. PointI.y ..")--"
839 output = output .. "\\draw[color="..colorNew..",fill ="..colorBack..", thick] " .. path .. "cycle;"
840 -- draw the new triangles composed by the edges of the polygon and the added point
842 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 ..");"
843 output = output .. "\\draw[color="..colorNew..", thick]".."(".. listPoints[polygon[i][1]].x .. "," .. listPoints[polygon[i][1]].y .. ") -- (" .. P.x .. "," .. P.y ..");"
844 output = output .. "\\draw[color="..colorNew..", thick]".."(".. listPoints[polygon[i][2]].x .. "," .. listPoints[polygon[i][2]].y .. ") -- (" .. P.x .. "," .. P.y ..");"
848 for i=1,#listPoints do
849 if(listPoints[i].type == "bbox") then
850 output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint^*_{" .. j .. "}$};"
853 output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
856 -- mark the added point
857 output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
862 function TeXaddOnePointMPBW(listPoints,P,step,bbox)
864 output = output .. "pair MeshPoints[];"
865 -- build the triangulation
866 triangulation = BowyerWatson(listPoints,bbox)
867 badTriangles = buildBadTriangles(P,triangulation)
868 for i=1,#listPoints do
869 output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ")*u;"
871 if(step == "badT") then
873 for i=1,#triangulation do
874 PointI = listPoints[triangulation[i][1]]
875 PointJ = listPoints[triangulation[i][2]]
876 PointK = listPoints[triangulation[i][3]]
877 if(triangulation[i].type == "bbox") then
878 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBbox;"
880 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
883 -- draw and fill the bad triangle
884 for i=1,#badTriangles do
885 PointI = listPoints[triangulation[badTriangles[i]][1]]
886 PointJ = listPoints[triangulation[badTriangles[i]][2]]
887 PointK = listPoints[triangulation[badTriangles[i]][3]]
888 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
889 output = output .. "fill (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBack;"
891 -- draw the circoncircle
892 for i=1,#badTriangles do
893 PointI = listPoints[triangulation[badTriangles[i]][1]]
894 PointJ = listPoints[triangulation[badTriangles[i]][2]]
895 PointK = listPoints[triangulation[badTriangles[i]][3]]
896 center, radius = circoncircle(PointI, PointJ, PointK)
897 output = output .. "draw fullcircle scaled ("..radius .."*2u) shifted ("..center.x .. "*u," .. center.y .. "*u) dashed evenly withcolor \\luameshmpcolorCircle;"
901 for i=1,#listPoints do
902 if(listPoints[i].type == "bbox") then
903 output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{" .. j .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
906 output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
909 -- mark the point to add
910 output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\luameshmpcolorNew;"
911 elseif(step == "cavity") then
912 polygon = buildCavity(badTriangles, triangulation)
913 polyNew = cleanPoly(polygon)
914 -- remove the bad triangles
915 for j=1,#badTriangles do
916 table.remove(triangulation,badTriangles[j]-(j-1))
918 -- draw the triangles
919 for i=1,#triangulation do
920 PointI = listPoints[triangulation[i][1]]
921 PointJ = listPoints[triangulation[i][2]]
922 PointK = listPoints[triangulation[i][3]]
923 if(triangulation[i].type == "bbox") then
924 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBbox;"
926 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
929 -- fill and draw the cavity
932 PointI = listPoints[polyNew[i]]
933 path = path .. "(".. PointI.x ..",".. PointI.y ..")*u--"
935 output = output .. "fill " .. path .. "cycle withcolor \\luameshmpcolorBack;"
936 output = output .. "draw " .. path .. "cycle withcolor \\luameshmpcolorNew withpen pencircle scaled 1pt;"
937 -- mark the points of the mesh
939 for i=1,#listPoints do
940 if(listPoints[i].type == "bbox") then
941 output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{" .. j .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
944 output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
947 -- mark the adding point
948 output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\luameshmpcolorNew ;"
949 elseif(step == "newT") then
950 polygon = buildCavity(badTriangles, triangulation)
951 polyNew = cleanPoly(polygon)
952 -- remove the bad triangles
953 for j=1,#badTriangles do
954 table.remove(triangulation,badTriangles[j]-(j-1))
956 -- draw the triangle of the triangulation
957 for i=1,#triangulation do
958 PointI = listPoints[triangulation[i][1]]
959 PointJ = listPoints[triangulation[i][2]]
960 PointK = listPoints[triangulation[i][3]]
961 if(triangulation[i].type == "bbox") then
962 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBbox;"
964 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
970 PointI = listPoints[polyNew[i]]
971 path = path .. "(".. PointI.x ..",".. PointI.y ..")*u--"
973 output = output .. "fill " .. path .. "cycle withcolor \\luameshmpcolorBack;"
974 -- draw the new triangles composed by the edges of the polygon and the added point
976 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 \\luameshmpcolorNew withpen pencircle scaled 1pt;"
977 output = output .. "draw".."(".. listPoints[polygon[i][1]].x .. "," .. listPoints[polygon[i][1]].y .. ")*u -- (" .. P.x .. "," .. P.y ..")*u withcolor \\luameshmpcolorNew withpen pencircle scaled 1pt;"
978 output = output .. "draw".."(".. listPoints[polygon[i][2]].x .. "," .. listPoints[polygon[i][2]].y .. ")*u -- (" .. P.x .. "," .. P.y ..")*u withcolor \\luameshmpcolorNew withpen pencircle scaled 1pt;"
982 for i=1,#listPoints do
983 if(listPoints[i].type == "bbox") then
984 output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{" .. j .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
987 output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
990 -- mark the added point
991 output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\luameshmpcolorNew ;"
997 -- build the list of points extern and stop at nbr
998 function buildListExt(chaine, stop)
1000 io.input(chaine) -- open the file
1001 text=io.read("*all")
1002 lines=string.explode(text,"\n+") -- all the lines
1003 for i=1,tonumber(stop) do
1004 xy=string.explode(lines[i]," +")
1005 table.insert(listPoints,{x=tonumber(xy[1]),y=tonumber(xy[2])})
1007 xy=string.explode(lines[stop+1]," +")
1008 point={x=tonumber(xy[1]),y=tonumber(xy[2])}
1009 return point, listPoints
1013 function TeXOnePointTikZBW(chaine,point,step,scale,mode,bbox,color,colorBack,colorNew,colorCircle,colorBbox)
1014 if(mode=="int") then
1015 Sx,Sy=string.match(point,"%((.+),(.+)%)")
1017 listPoints = buildList(chaine, mode)
1019 -- point is a number
1020 P, listPoints = buildListExt(chaine,tonumber(point))
1022 output = TeXaddOnePointTikZ(listPoints,P,step,bbox,color,colorBack,colorNew,colorCircle,colorBbox)
1023 output = "\\noindent\\begin{tikzpicture}[x="..scale..",y="..scale.."]".. output .. "\\end{tikzpicture}"
1027 function TeXOnePointTikZBWinc(chaine,point,beginning, ending,step,scale,mode,bbox,color,colorBack,colorNew,colorCircle,colorBbox)
1028 if(mode=="int") then
1029 Sx,Sy=string.match(point,"%((.+),(.+)%)")
1031 listPoints = buildList(chaine, mode)
1033 -- point is a number
1034 P, listPoints = buildListExt(chaine,tonumber(point))
1036 output = TeXaddOnePointTikZ(listPoints,P,step,bbox,color,colorBack,colorNew,colorCircle,colorBbox)
1037 output = "\\noindent\\begin{tikzpicture}[x="..scale..",y="..scale.."]".. beginning..output ..ending.. "\\end{tikzpicture}"
1041 function TeXOnePointMPBW(chaine,point,step,scale,mode,bbox)
1042 if(mode=="int") then
1043 Sx,Sy=string.match(point,"%((.+),(.+)%)")
1045 listPoints = buildList(chaine, mode)
1047 -- point is a number
1048 P, listPoints = buildListExt(chaine,tonumber(point))
1050 output = TeXaddOnePointMPBW(listPoints,P,step,bbox)
1051 output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale..";".. output .. "endfig;\\end{mplibcode}"
1055 function TeXOnePointMPBWinc(chaine,point,beginning,ending,step,scale,mode,bbox)
1056 if(mode=="int") then
1057 Sx,Sy=string.match(point,"%((.+),(.+)%)")
1059 listPoints = buildList(chaine, mode)
1061 -- point is a number
1062 P, listPoints = buildListExt(chaine,tonumber(point))
1064 output = TeXaddOnePointMPBW(listPoints,P,step,bbox)
1065 output = "\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"