X-Git-Url: https://melusine.eu.org/syracuse/G/git/?a=blobdiff_plain;f=luamesh.lua;h=179024daa2f26712cf60cdabb1a7d6548f13ee9c;hb=93a6d2fd77327d5abb1ee4a8324e5d0f79f29927;hp=83c2c257a34ea2ae9be4511885f314ddf3355ebc;hpb=8f69ee8fd8ba6b022d1aaae4c4f2b93e1694765f;p=delaunay.git

diff --git a/luamesh.lua b/luamesh.lua
index 83c2c25..179024d 100644
--- a/luamesh.lua
+++ b/luamesh.lua
@@ -207,6 +207,81 @@ function circoncircle(M, N, P)
    return Out, rad  -- (center [Point], R [float])
 end
 
+-- compute the list of the circumcircle of a triangulation
+function listCircumCenter(listPoints,triangulation)
+   list = {}
+   for j=1,#triangulation do
+      A = listPoints[triangulation[j][1]]
+      B = listPoints[triangulation[j][2]]
+      C = listPoints[triangulation[j][3]]
+      center, radius = circoncircle(A,B,C)
+      table.insert(list,{x=center.x,y=center.y,r=radius})
+   end
+   return list
+end
+
+-- find the three neighbour triangles of T
+function findNeighbour(T,i,triangulation)
+   -- T : triangle
+   -- i : index of T in triangualation
+   -- triangulation
+
+   list = {}
+   -- define the three edge
+   e1 = {T[1],T[2]}
+   e2 = {T[2],T[3]}
+   e3 = {T[3],T[1]}
+   for j=1,#triangulation do
+      if j~= i then
+         if(edgeInTriangle(e1,triangulation[j])) then
+            table.insert(list,j)
+         end
+         if(edgeInTriangle(e2,triangulation[j])) then
+            table.insert(list,j)
+         end
+         if(edgeInTriangle(e3,triangulation[j])) then
+            table.insert(list,j)
+         end
+      end
+   end
+   return list
+end
+
+-- test if edge are the same (reverse)
+function equalEdge(e1,e2)
+   if(((e1[1] == e2[1]) and (e1[2] == e2[2])) or ((e1[1] == e2[2]) and (e1[2] == e2[1]))) then
+      return true
+   else
+      return false
+   end
+end
+
+-- test if the edge belongs to the list
+function edgeInList(e,listE)
+   output = false
+   for i=1,#listE do
+      if(equalEdge(e,listE[i])) then
+         output = true
+      end
+   end
+   return output
+end
+
+-- build the edges of the Voronoi diagram with a given triangulation
+function buildVoronoi(listPoints, triangulation)
+   listCircumCircle = listCircumCenter(listPoints, triangulation)
+   listVoronoi = {}
+   for i=1,#listCircumCircle do
+      listN = findNeighbour(triangulation[i],i,triangulation)
+      for j=1,#listN do
+         edge = {i,listN[j]}
+         if( not edgeInList(edge, listVoronoi)) then
+            table.insert(listVoronoi, edge)
+         end
+      end
+   end
+   return listVoronoi
+end
 
 -------------------------- TeX
 -- build the list of points
@@ -252,9 +327,179 @@ function rectangleList(a,b,nbrA,nbrB)
    return listPoints
 end
 
+
+-- trace Voronoi with MP
+function traceVoronoiMP(listPoints, triangulation,listVoronoi, points, tri)
+   listCircumC = listCircumCenter(listPoints,triangulation)
+   output = "";
+   output = output .. " pair MeshPoints[];"
+   for i=1,#listPoints do
+      output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ")*u;"
+   end
+   output = output .. " pair CircumCenters[];"
+   for i=1,#listCircumC do
+      output = output .. "CircumCenters[".. i .. "] = (" .. listCircumC[i].x .. "," .. listCircumC[i].y .. ")*u;"
+   end
+   if(tri=="show") then
+      for i=1,#triangulation do
+         PointI = listPoints[triangulation[i][1]]
+         PointJ = listPoints[triangulation[i][2]]
+         PointK = listPoints[triangulation[i][3]]
+         if(triangulation[i].type == "bbox") then
+            output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBbox;"
+         else
+            output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
+         end
+      end
+   end
+   for i=1,#listVoronoi do
+      PointI = listCircumC[listVoronoi[i][1]]
+      PointJ = listCircumC[listVoronoi[i][2]]
+      output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u withcolor \\luameshmpcolorVoronoi;"
+   end
+   if(points=="points") then
+      j=1
+      for i=1,#listPoints do
+         if(listPoints[i].type == "bbox") then
+            output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{"..j.."}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
+            j=j+1
+         else
+            output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
+         end
+      end
+      for i=1,#listCircumC do
+         output = output .. "dotlabel.llft (btex $\\CircumPoint_{" .. i .. "}$ etex, (" .. listCircumC[i].x ..",".. listCircumC[i].y .. ")*u ) withcolor \\luameshmpcolorVoronoi ;"
+      end
+   else
+      j=1
+      for i=1,#listPoints do
+         if(listPoints[i].type == "bbox") then
+            output = output .. "drawdot  (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u  withcolor \\luameshmpcolorBbox withpen pencircle scaled 3;"
+            j=j+1
+         else
+            output = output .. "drawdot  (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u  withcolor \\luameshmpcolor withpen pencircle scaled 3;"
+         end
+      end
+      for i=1,#listCircumC do
+         output = output .. "drawdot  (" .. listCircumC[i].x ..",".. listCircumC[i].y .. ")*u  withcolor \\luameshmpcolorVoronoi withpen pencircle scaled 3;"
+      end
+   end
+
+   return output
+end
+
+
+-- trace Voronoi with TikZ
+function traceVoronoiTikZ(listPoints, triangulation,listVoronoi, points, tri,color,colorBbox,colorVoronoi)
+   listCircumC = listCircumCenter(listPoints,triangulation)
+    output = ""
+   for i=1,#listPoints do
+      output = output .. "\\coordinate (MeshPoints".. i .. ") at  (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
+   end
+   for i=1,#listCircumC do
+      output = output .. "\\coordinate (CircumPoints".. i .. ") at  (" .. listCircumC[i].x .. "," .. listCircumC[i].y .. ");"
+   end
+   if(tri=="show") then
+      for i=1,#triangulation do
+         PointI = listPoints[triangulation[i][1]]
+         PointJ = listPoints[triangulation[i][2]]
+         PointK = listPoints[triangulation[i][3]]
+         if(triangulation[i].type == "bbox") then
+            output = output .. "\\draw[color="..colorBbox.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
+         else
+            output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
+         end
+      end
+   end
+   for i=1,#listVoronoi do
+      PointI = listCircumC[listVoronoi[i][1]]
+      PointJ = listCircumC[listVoronoi[i][2]]
+      output = output .. "\\draw[color="..colorVoronoi.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..");"
+   end
+   if(points=="points") then
+      j=1
+      for i=1,#listPoints do
+         if(listPoints[i].type == "bbox") then
+            if(listPoints[i].type == "bbox") then
+               output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint^*_{" .. j .. "}$};"
+               j=j+1
+            else
+               output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
+            end
+         end
+      end
+      for i=1,#listCircumC do
+         output = output .. "\\draw[color="..colorVoronoi.."] (" .. listCircumC[i].x ..",".. listCircumC[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\CircumPoint_{" .. i .. "}$};"
+      end
+   else
+      j=1
+      for i=1,#listPoints do
+         if(listPoints[i].type == "bbox") then
+            output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} ;"
+            j=j+1
+         else
+            output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} ;"
+         end
+      end
+      for i=1,#listCircumC do
+         output = output .. "\\draw[color="..colorVoronoi.."] (" .. listCircumC[i].x ..",".. listCircumC[i].y .. ") node {$\\bullet$};"
+      end
+   end
+   return output
+end
+
+
+
+-- buildVoronoi with MP
+function buildVoronoiMPBW(chaine,mode,points,bbox,scale,tri)
+   listPoints = buildList(chaine, mode)
+   triangulation = BowyerWatson(listPoints,bbox)
+   listVoronoi = buildVoronoi(listPoints, triangulation)
+   output = traceVoronoiMP(listPoints,triangulation,listVoronoi,points,tri)
+   output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale.. ";" .. output .."endfig;\\end{mplibcode}"
+   tex.sprint(output)
+end
+
+
+-- buildVoronoi with TikZ
+function buildVoronoiTikZBW(chaine,mode,points,bbox,scale,tri,color,colorBbox,colorVoronoi)
+   listPoints = buildList(chaine, mode)
+   triangulation = BowyerWatson(listPoints,bbox)
+   listVoronoi = buildVoronoi(listPoints, triangulation)
+   output = traceVoronoiTikZ(listPoints,triangulation,listVoronoi,points,tri,color,colorBbox,colorVoronoi)
+   output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" .. output .."\\end{tikzpicture}"   tex.sprint(output)
+end
+
+
+-- buildVoronoi with MP
+function buildVoronoiMPBWinc(chaine,beginning, ending,mode,points,bbox,scale,tri)
+   listPoints = buildList(chaine, mode)
+   triangulation = BowyerWatson(listPoints,bbox)
+   listVoronoi = buildVoronoi(listPoints, triangulation)
+   output = traceVoronoiMP(listPoints,triangulation,listVoronoi,points,tri)
+   output = "\\leavevmode\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"
+   tex.sprint(output)
+end
+
+
+-- buildVoronoi with TikZ
+function buildVoronoiTikZBWinc(chaine,beginning, ending,mode,points,bbox,scale,tri,color,colorBbox,colorVoronoi)
+   listPoints = buildList(chaine, mode)
+   triangulation = BowyerWatson(listPoints,bbox)
+   listVoronoi = buildVoronoi(listPoints, triangulation)
+   output = traceVoronoiTikZ(listPoints,triangulation,listVoronoi,points,tri,color,colorBbox,colorVoronoi)
+   output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" ..beginning.. output..ending .."\\end{tikzpicture}"
+   tex.sprint(output)
+end
+
+
+
 -- trace a triangulation with TikZ
 function traceMeshTikZ(listPoints, triangulation,points,color,colorBbox)
    output = ""
+   for i=1,#listPoints do
+      output = output .. "\\coordinate (MeshPoints".. i .. ") at  (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
+   end
    for i=1,#triangulation do
       PointI = listPoints[triangulation[i][1]]
       PointJ = listPoints[triangulation[i][2]]
@@ -285,7 +530,7 @@ function traceMeshMP(listPoints, triangulation,points)
    output = "";
    output = output .. " pair MeshPoints[];"
    for i=1,#listPoints do
-      output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
+      output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ")*u;"
    end
 
    for i=1,#triangulation do
@@ -355,7 +600,7 @@ function tracePointsMP(listPoints,points)
    output = "";
    output = output .. " pair MeshPoints[];"
    for i=1,#listPoints do
-      output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
+      output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ")*u;"
    end
    if(points=="points") then
       j=1
@@ -382,6 +627,9 @@ end
 -- print points of the mesh
 function tracePointsTikZ(listPoints,points,color,colorBbox)
    output = "";
+   for i=1,#listPoints do
+      output = output .. "\\coordinate (MeshPoints".. i .. ") at  (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
+   end
    if(points=="points") then
       j=1
       for i=1,#listPoints do
@@ -423,7 +671,7 @@ function printPointsMPinc(chaine,beginning, ending, mode,points,bbox,scale)
       listPoints = buildBoundingBox(listPoints)
    end
    output = tracePointsMP(listPoints,points)
-   output = "\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"
+   output = "\\leavevmode\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"
    tex.sprint(output)
 end
 
@@ -507,6 +755,9 @@ function TeXaddOnePointTikZ(listPoints,P,step,bbox,color,colorBack, colorNew, co
    -- build the triangulation
    triangulation = BowyerWatson(listPoints,bbox)
    badTriangles = buildBadTriangles(P,triangulation)
+   for i=1,#listPoints do
+      output = output .. "\\coordinate (MeshPoints".. i .. ") at  (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
+   end
    if(step == "badT") then
       -- draw all triangle
       for i=1,#triangulation do
@@ -637,7 +888,7 @@ function TeXaddOnePointMPBW(listPoints,P,step,bbox)
    triangulation = BowyerWatson(listPoints,bbox)
    badTriangles = buildBadTriangles(P,triangulation)
    for i=1,#listPoints do
-      output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
+      output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ")*u;"
    end
    if(step == "badT") then
       -- draw all triangle