1 require "luamesh-polygon"
4 -- Bowyer and Watson algorithm
6 function BowyerWatson (listPoints,bbox)
7 local triangulation = {}
8 local lgth = #listPoints
9 -- add four points to listPoints to have a bounding box
10 listPoints = buildBoundingBox(listPoints)
12 triangulation[1] = {lgth+1, lgth+2, lgth+3,type="bbox"}
13 -- the second triangle
14 triangulation[2] = {lgth+1, lgth+3, lgth+4,type="bbox"}
15 -- add points one by one
17 -- find the triangles which the circumcircle contained the point to add
18 badTriangles = buildBadTriangles(listPoints[i],triangulation)
19 -- build the polygon of the cavity containing the point to add
20 polygon = buildCavity(badTriangles, triangulation)
21 -- remove the bad triangles
22 for j=1,#badTriangles do
23 table.remove(triangulation,badTriangles[j]-(j-1))
25 -- build the new triangles and add them to triangulation
27 if((polygon[j][1]>lgth) or (polygon[j][2]>lgth) or (i>lgth)) then
28 table.insert(triangulation,{polygon[j][1],polygon[j][2],i,type="bbox"})
30 table.insert(triangulation,{polygon[j][1],polygon[j][2],i,type="in"})
33 end -- end adding points of the listPoints
34 -- remove bounding box
35 if(bbox ~= "bbox") then
36 triangulation = removeBoundingBox(triangulation,lgth)
37 table.remove(listPoints,lgth+1)
38 table.remove(listPoints,lgth+1)
39 table.remove(listPoints,lgth+1)
40 table.remove(listPoints,lgth+1)
46 function buildBoundingBox(listPoints)
47 -- listPoints : list of points
48 -- epsV : parameter for the distance of the bounding box
49 local xmin, xmax, ymin, ymax, eps
54 for i=1,#listPoints do
55 if (listPoints[i].x < xmin) then
56 xmin = listPoints[i].x
58 if (listPoints[i].x > xmax) then
59 xmax = listPoints[i].x
61 if (listPoints[i].y < ymin) then
62 ymin = listPoints[i].y
64 if (listPoints[i].y > ymax) then
65 ymax = listPoints[i].y
68 eps = math.max(math.abs(xmax-xmin),math.abs(ymax-ymin))*0.15
73 -- add points of the bounding box in last positions
74 table.insert(listPoints,{x=xmin,y=ymin,type="bbox"})
75 table.insert(listPoints,{x=xmin,y=ymax,type="bbox"})
76 table.insert(listPoints,{x=xmax,y=ymax,type="bbox"})
77 table.insert(listPoints,{x=xmax,y=ymin,type="bbox"})
81 function BoundingBox(listPoints)
82 -- listPoints : list of points
83 -- epsV : parameter for the distance of the bounding box
84 local xmin, xmax, ymin, ymax, eps
89 for i=1,#listPoints do
90 if (listPoints[i].x < xmin) then
91 xmin = listPoints[i].x
93 if (listPoints[i].x > xmax) then
94 xmax = listPoints[i].x
96 if (listPoints[i].y < ymin) then
97 ymin = listPoints[i].y
99 if (listPoints[i].y > ymax) then
100 ymax = listPoints[i].y
103 eps = math.max(math.abs(xmax-xmin),math.abs(ymax-ymin))*0.15
108 return xmin, xmax, ymin, ymax
111 function removeBoundingBox(triangulation,lgth)
112 -- build the four bounding box edge
118 local newTriangulation = {}
119 for i=1,#triangulation do
120 boolE1 = pointInTriangle(point1,triangulation[i])
121 boolE2 = pointInTriangle(point2,triangulation[i])
122 boolE3 = pointInTriangle(point3,triangulation[i])
123 boolE4 = pointInTriangle(point4,triangulation[i])
124 if((not boolE1) and (not boolE2) and (not boolE3) and (not boolE4)) then
125 table.insert(newTriangulation,triangulation[i])
128 return newTriangulation
132 function buildBadTriangles(point, triangulation)
133 local badTriangles = {}
134 for j=1,#triangulation do -- for all triangles
135 A = listPoints[triangulation[j][1]]
136 B = listPoints[triangulation[j][2]]
137 C = listPoints[triangulation[j][3]]
138 center, radius = circoncircle(A,B,C)
139 CP = Vector(center,point)
140 if(VectorNorm(CP)<radius) then -- the point belongs to the circoncirle
141 table.insert(badTriangles,j)
147 -- construction of the cavity composed by the bad triangles around the point to add
148 function buildCavity(badTriangles, triangulation)
150 for j=1,#badTriangles do -- for all bad triangles
151 ind = badTriangles[j]
152 for k=1,3 do -- for all edges
153 edge = {triangulation[ind][k],triangulation[ind][k%3+1]}
155 for l = 1,#badTriangles do -- for all badtriangles
156 badInd = badTriangles[l]
157 if(badInd ~= ind) then -- if not the current one
158 edgeBord = edgeBord or edgeInTriangle(edge,triangulation[badInd])
161 -- if the edge does not belong to another bad triangle
162 if(edgeBord == false) then
163 -- insert the edge to the cavity
164 table.insert(polygon,edge)
171 function edgeInTriangle(e,t)
186 function pointInTriangle(e,t)
198 local out = {x = B.x - A.x, y = B.y - A.y}
202 function VectorNorm(NP)
203 return math.sqrt(NP.x*NP.x +NP.y*NP.y)
207 function circoncircle(M, N, P)
208 -- Compute center and radius of the circoncircle of the triangle M N P
210 -- return : (center [Point],radius [float])
212 local MN = Vector(M,N)
213 local NP = Vector(N,P)
214 local PM = Vector(P,M)
215 m = VectorNorm(NP) -- |NP|
216 n = VectorNorm(PM) -- |PM|
217 p = VectorNorm(MN) -- |MN|
219 d = (m + n + p) * (-m + n + p) * (m - n + p) * (m + n - p)
221 rad = m * n * p / math.sqrt(d)
225 d = -2 * (M.x * NP.y + N.x * PM.y + P.x * MN.y)
230 om2 = math.pow(VectorNorm(OM),2) -- |OM|**2
231 on2 = math.pow(VectorNorm(ON),2) -- |ON|**2
232 op2 = math.pow(VectorNorm(OP),2) -- |OP|**2
233 x0 = -(om2 * NP.y + on2 * PM.y + op2 * MN.y) / d
234 y0 = (om2 * NP.x + on2 * PM.x + op2 * MN.x) / d
240 return Out, rad -- (center [Point], R [float])
243 -- compute the list of the circumcircle of a triangulation
244 function listCircumCenter(listPoints,triangulation)
246 for j=1,#triangulation do
247 A = listPoints[triangulation[j][1]]
248 B = listPoints[triangulation[j][2]]
249 C = listPoints[triangulation[j][3]]
250 center, radius = circoncircle(A,B,C)
251 table.insert(list,{x=center.x,y=center.y,r=radius})
256 -- find the three neighbour triangles of T
257 function findNeighbour(T,i,triangulation)
259 -- i : index of T in triangualation
263 -- define the three edge
267 for j=1,#triangulation do
269 if(edgeInTriangle(e1,triangulation[j])) then
272 if(edgeInTriangle(e2,triangulation[j])) then
275 if(edgeInTriangle(e3,triangulation[j])) then
283 -- test if edge are the same (reverse)
284 function equalEdge(e1,e2)
285 if(((e1[1] == e2[1]) and (e1[2] == e2[2])) or ((e1[1] == e2[2]) and (e1[2] == e2[1]))) then
292 -- test if the edge belongs to the list
293 function edgeInList(e,listE)
296 if(equalEdge(e,listE[i])) then
303 -- build the edges of the Voronoi diagram with a given triangulation
304 function buildVoronoi(listPoints, triangulation)
305 local listCircumCircle = listCircumCenter(listPoints, triangulation)
306 local listVoronoi = {}
307 for i=1,#listCircumCircle do
308 listN = findNeighbour(triangulation[i],i,triangulation)
311 if( not edgeInList(edge, listVoronoi)) then
312 table.insert(listVoronoi, edge)
319 -------------------------- TeX
320 -- build the list of points
321 function buildList(chaine, mode)
322 -- if mode = int : the list is given in the chaine string (x1,y1);(x2,y2);...;(xn,yn)
323 -- if mode = ext : the list is given in a file line by line with space separation
324 local listPoints = {}
325 if mode == "int" then
326 local points = string.explode(chaine, ";")
329 Sx,Sy=string.match(points[i],"%((.+),(.+)%)")
330 listPoints[i]={x=tonumber(Sx),y=tonumber(Sy)}
332 elseif mode == "ext" then
333 io.input(chaine) -- open the file
335 lines=string.explode(text,"\n+") -- all the lines
338 xy=string.explode(lines[i]," +")
339 listPoints[i]={x=tonumber(xy[1]),y=tonumber(xy[2])}
342 print("Non existing mode")
348 -- function to build a gridpoints from the bounding box
350 function buildGrid(listPoints,h)
351 -- listPoints : list of the points of the polygon, ordered
352 -- h : parameter for the grid
353 xmin, xmax, ymin, ymax = BoundingBox(listPoints)
355 local grid = rectangleList(xmin,xmax,ymin,ymax,h)
359 -- function to build the list of points in the rectangle
360 function rectangleList(xmin,xmax,ymin,ymax,h)
361 nbrX = math.floor(math.abs(xmax-xmin)/h)
362 nbrY = math.floor(math.abs(ymax-ymin)/h)
363 local listPoints = {}
367 listPoints[k] = {x = xmin+(i-1)*h, y=ymin+(j-1)*h}
375 -- function to add points from a grid to the interior of a polygon
376 function addGridPoints(polygon, grid)
377 local listPoints = polygon
380 --print(grid[i].x,grid[i].y)
381 --print(isInside(polygon,grid[i]))
382 if(isInside(polygon,grid[i])) then
384 listPoints[k] = grid[i]
391 local function shallowCopy(original)
393 for key, value in pairs(original) do
399 -- function give a real polygon without repeting points
400 function cleanPoly(polygon)
402 polyCopy = shallowCopy(polygon)
405 table.insert(polyNew, e1)
406 table.insert(polyNew, e2)
407 table.remove(polyCopy,1)
412 while (i<=#polyCopy and find==false) do
413 bool1 = (polyCopy[i][1] == polyNew[j])
414 bool2 = (polyCopy[i][2] == polyNew[j])
415 if(bool1 or bool2) then -- the edge has a common point with polyNew[j]
417 table.insert(polyNew, polyCopy[i][1])
419 table.remove(polyCopy,i)
421 elseif(not bool2) then
422 table.insert(polyNew, polyCopy[i][2])
424 table.remove(polyCopy,i)
436 -- build the list of points extern and stop at nbr
437 function buildListExt(chaine, stop)
439 io.input(chaine) -- open the file
441 lines=string.explode(text,"\n+") -- all the lines
442 for i=1,tonumber(stop) do
443 xy=string.explode(lines[i]," +")
444 table.insert(listPoints,{x=tonumber(xy[1]),y=tonumber(xy[2])})
446 xy=string.explode(lines[stop+1]," +")
447 point={x=tonumber(xy[1]),y=tonumber(xy[2])}
448 return point, listPoints
451 function split(pString, pPattern)
452 local Table = {} -- NOTE: use {n = 0} in Lua-5.0
453 local fpat = "(.-)" .. pPattern
455 local s, e, cap = pString:find(fpat, 1)
457 if s ~= 1 or cap ~= "" then
458 table.insert(Table,cap)
461 s, e, cap = pString:find(fpat, last_end)
463 if last_end <= #pString then
464 cap = pString:sub(last_end)
465 table.insert(Table, cap)
470 function readGmsh(file)
471 io.input(file) -- open the file
473 local lines = split(text,"\n+") -- all the lines
482 if(lines[i+J] == "$EndNodes") then
484 -- go to the next line
486 if(boolNodes) then -- we are in the Nodes environment
487 xy=split(lines[i+J]," +")
488 table.insert(listPoints,{x=tonumber(xy[2]),y=tonumber(xy[3])})
490 if(lines[i+J] == "$Nodes") then
492 -- go to the next line
495 if(lines[i+J] == "$EndElements") then
497 -- go to the next line
499 if(boolElements) then -- we are in the Nodes environment
500 xy=split(lines[i+J]," +")
501 if(tonumber(xy[2]) == 2) then -- if the element is a triangle
503 table.insert(triangulation,{tonumber(xy[2+nbrTags+1]),tonumber(xy[2+nbrTags+2]),tonumber(xy[2+nbrTags+3])})
506 if(lines[i+J] == "$Elements") then
508 -- go to the next line
512 return listPoints, triangulation