-- Bowyer and Watson algorithm
-- Delaunay meshing
function BowyerWatson (listPoints,bbox)
   local triangulation = {}
   local lgth = #listPoints
   -- add four points to listPoints to have a bounding box
   listPoints = buildBoundingBox(listPoints)
   -- the first triangle
   triangulation[1] = {lgth+1, lgth+2, lgth+3,type="bbox"}
   -- the second triangle
   triangulation[2] = {lgth+1, lgth+3, lgth+4,type="bbox"}
   -- add points one by one
   for i=1,lgth do
      -- find the triangles which the circumcircle contained the point to add
      badTriangles = buildBadTriangles(listPoints[i],triangulation)
      -- build the polygon of the cavity containing the point to add
      polygon = buildCavity(badTriangles, triangulation)
      -- remove the bad triangles
      for j=1,#badTriangles do
         table.remove(triangulation,badTriangles[j]-(j-1))
      end
      -- build the new triangles and add them to triangulation
      for j=1,#polygon do
         if((polygon[j][1]>lgth) or (polygon[j][2]>lgth) or (i>lgth)) then
            table.insert(triangulation,{polygon[j][1],polygon[j][2],i,type="bbox"})
         else
            table.insert(triangulation,{polygon[j][1],polygon[j][2],i,type="in"})
         end
      end
   end	-- end adding points of the listPoints
   -- remove bounding box
   if(bbox ~= "bbox") then
      triangulation = removeBoundingBox(triangulation,lgth)
      table.remove(listPoints,lgth+1)
      table.remove(listPoints,lgth+1)
      table.remove(listPoints,lgth+1)
      table.remove(listPoints,lgth+1)
   end
   return triangulation
end


function buildBoundingBox(listPoints)
   -- listPoints : list of points
   -- epsV        : parameter for the distance of the bounding box
   local xmin, xmax, ymin, ymax, eps
   xmin = 1000
   ymin = 1000
   xmax = -1000
   ymax = -1000
   for i=1,#listPoints do
      if (listPoints[i].x < xmin) then
         xmin = listPoints[i].x
      end
      if (listPoints[i].x > xmax) then
         xmax = listPoints[i].x
      end
      if (listPoints[i].y < ymin) then
         ymin = listPoints[i].y
      end
      if (listPoints[i].y > ymax) then
         ymax = listPoints[i].y
      end
   end
   eps = math.max(math.abs(xmax-xmin),math.abs(ymax-ymin))*0.15
   xmin = xmin - eps
   xmax = xmax + eps
   ymin = ymin - eps
   ymax = ymax + eps
   -- add points of the bounding box in last positions
   table.insert(listPoints,{x=xmin,y=ymin,type="bbox"})
   table.insert(listPoints,{x=xmin,y=ymax,type="bbox"})
   table.insert(listPoints,{x=xmax,y=ymax,type="bbox"})
   table.insert(listPoints,{x=xmax,y=ymin,type="bbox"})
   return listPoints
end

function removeBoundingBox(triangulation,lgth)
   -- build the four bounding box edge
   point1 = lgth+1
   point2 = lgth+2
   point3 = lgth+3
   point4 = lgth+4
   -- for all triangle
   newTriangulation = {}
   for i=1,#triangulation do
      boolE1 = pointInTriangle(point1,triangulation[i])
      boolE2 = pointInTriangle(point2,triangulation[i])
      boolE3 = pointInTriangle(point3,triangulation[i])
      boolE4 = pointInTriangle(point4,triangulation[i])
      if((not boolE1) and (not boolE2) and (not boolE3) and (not boolE4)) then
         table.insert(newTriangulation,triangulation[i])
      end
   end
   return newTriangulation
end


function buildBadTriangles(point, triangulation)
   badTriangles = {}
   for j=1,#triangulation do -- for all triangles
      A = listPoints[triangulation[j][1]]
      B = listPoints[triangulation[j][2]]
      C = listPoints[triangulation[j][3]]
      center, radius = circoncircle(A,B,C)
      CP = Vector(center,point)
      if(VectorNorm(CP)<radius) then -- the point belongs to the circoncirle
         table.insert(badTriangles,j)
      end
   end
   return badTriangles
end

-- construction of the cavity composed by the bad triangles around the point to add
function buildCavity(badTriangles, triangulation)
   polygon = {}
   for j=1,#badTriangles do -- for all bad triangles
      ind = badTriangles[j]
      for k=1,3 do -- for all edges
         edge = {triangulation[ind][k],triangulation[ind][k%3+1]}
         edgeBord = false
         for l = 1,#badTriangles do -- for all badtriangles
            badInd = badTriangles[l]
            if(badInd ~= ind) then -- if not the current one
               edgeBord = edgeBord or edgeInTriangle(edge,triangulation[badInd])
            end
         end --
         -- if the edge does not belong to another bad triangle
         if(edgeBord == false) then
            -- insert the edge to the cavity
            table.insert(polygon,edge)
         end
      end --
   end --
   return polygon
end

function edgeInTriangle(e,t)
   in1 = false
   in2 = false
   for i=1,3 do
      if e[1] == t[i] then
         in1 = true
      end
      if e[2] == t[i] then
         in2 = true
      end
   end
   out = (in1 and in2)
   return out
end

function pointInTriangle(e,t)
   in1 = false
   for i=1,3 do
      if e == t[i] then
         in1 = true
      end
   end
   return in1
end


function Vector(A,B)
   local out = {x = B.x - A.x, y = B.y - A.y}
   return out
end

function VectorNorm(NP)
   return math.sqrt(NP.x*NP.x +NP.y*NP.y)
end

-- circoncircle
function circoncircle(M, N, P)
   -- Compute center and radius of the circoncircle of the triangle M N P

   -- return : (center [Point],radius [float])

   local MN = Vector(M,N)
   local NP = Vector(N,P)
   local PM = Vector(P,M)
   m = VectorNorm(NP)  -- |NP|
   n = VectorNorm(PM)  -- |PM|
   p = VectorNorm(MN)  -- |MN|

   d = (m + n + p) * (-m + n + p) * (m - n + p) * (m + n - p)
   if d > 0 then
      rad = m * n * p / math.sqrt(d)
   else
      rad = 0
   end
   d = -2 * (M.x * NP.y + N.x * PM.y + P.x * MN.y)
   O = {x=0, y=0}
   OM = Vector(O, M)
   ON = Vector(O, N)
   OP = Vector(O, P)
   om2 = math.pow(VectorNorm(OM),2)  -- |OM|**2
   on2 = math.pow(VectorNorm(ON),2)  -- |ON|**2
   op2 = math.pow(VectorNorm(OP),2)  -- |OP|**2
   x0 = -(om2 * NP.y + on2 * PM.y + op2 * MN.y) / d
   y0 = (om2 * NP.x + on2 * PM.x + op2 * MN.x) / d
   if d == 0 then
      Out = {nil, nil}
   else
      Out = {x=x0, y=y0}
   end
   return Out, rad  -- (center [Point], R [float])
end


-------------------------- TeX
-- build the list of points
function buildList(chaine, mode)
   -- if mode = int : the list is given in the chaine string (x1,y1);(x2,y2);...;(xn,yn)
   -- if mode = ext : the list is given in a file line by line with space separation
   listPoints = {}
   if mode == "int" then
      local points = string.explode(chaine, ";")
      local lgth=#points
      for i=1,lgth do
         Sx,Sy=string.match(points[i],"%((.+),(.+)%)")
         listPoints[i]={x=tonumber(Sx),y=tonumber(Sy)}
      end
   elseif mode == "ext" then
      io.input(chaine) -- open the file
      text=io.read("*all")
      lines=string.explode(text,"\n+") -- all the lines
      tablePoints={}
      for i=1,#lines do
	 xy=string.explode(lines[i]," +")
	 listPoints[i]={x=tonumber(xy[1]),y=tonumber(xy[2])}
      end
   else
      print("Non existing mode")
   end
   return listPoints
end


--
function rectangleList(a,b,nbrA,nbrB)
   stepA = a/nbrA
   stepB = b/nbrB
   listPoints = {}
   k=1
   for i=1,(nbrA+1) do
      for j=1,(nbrB+1) do
         listPoints[k] = {x = (i-1)*stepA, y=(j-1)*stepB}
         k=k+1
      end
   end
   return listPoints
end

-- trace a triangulation with TikZ
function traceMeshTikZ(listPoints, triangulation,points,color,colorBbox)
   output = ""
   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
   if(points=="points") then
      for i=1,#listPoints do
         if(listPoints[i].type == "bbox") then
            output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
         else
            output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
         end
      end
   end
   return output
end


-- trace a triangulation with MP
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 .. ");"
   end

   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
   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
   end
   return output
end


-- buildMesh with MP
function buildMeshMPBW(chaine,mode,points,bbox,scale)
   listPoints = buildList(chaine, mode)
   triangulation = BowyerWatson(listPoints,bbox)
   output = traceMeshMP(listPoints, triangulation,points)
   output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale.. ";" .. output .."endfig;\\end{mplibcode}"
   tex.sprint(output)
end

-- buildMesh with MP include code
function buildMeshMPBWinc(chaine,beginning, ending,mode,points,bbox,scale)
   listPoints = buildList(chaine, mode)
   triangulation = BowyerWatson(listPoints,bbox)
   output = traceMeshMP(listPoints, triangulation,points)
   output = "\\leavevmode\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"
   tex.sprint(output)
end

-- buildMesh with TikZ
function buildMeshTikZBW(chaine,mode,points,bbox,scale,color,colorBbox)
   listPoints = buildList(chaine, mode)
   triangulation = BowyerWatson(listPoints,bbox)
   output = traceMeshTikZ(listPoints, triangulation,points,color,colorBbox)
   output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" .. output .."\\end{tikzpicture}"
   tex.sprint(output)
end

-- buildMesh with TikZ
function buildMeshTikZBWinc(chaine,beginning, ending,mode,points,bbox,scale,color,colorBbox)
   listPoints = buildList(chaine, mode)
   triangulation = BowyerWatson(listPoints,bbox)
   output = traceMeshTikZ(listPoints, triangulation,points,color,colorBbox)
   output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" ..beginning.. output..ending .."\\end{tikzpicture}"
   tex.sprint(output)
end


-- print points of the mesh
function tracePointsMP(listPoints,points)
   output = "";
   output = output .. " pair MeshPoints[];"
   for i=1,#listPoints do
      output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
   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
   else
      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;"
         else
            output = output .. "drawdot (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u  withcolor \\luameshmpcolor withpen pencircle scaled 3;"
         end
      end
   end
   return output
end

-- print points of the mesh
function tracePointsTikZ(listPoints,points,color,colorBbox)
   output = "";
   if(points=="points") then
      for i=1,#listPoints do
         if(listPoints[i].type == "bbox") then
            output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
         else
            output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
         end
      end
   else
      for i=1,#listPoints do
         if(listPoints[i].type == "bbox") then
            output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} ;"
         else
            output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} ;"
         end
      end
   end
   return output
end

-- print points to mesh
function printPointsMP(chaine,mode,points,bbox,scale)
   listPoints = buildList(chaine, mode)
   if(bbox == "bbox" ) then
      listPoints = buildBoundingBox(listPoints)
   end
   output = tracePointsMP(listPoints,points)
   output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale.. ";" .. output .."endfig;\\end{mplibcode}"
   tex.sprint(output)
end


-- print points to mesh
function printPointsMPinc(chaine,beginning, ending, mode,points,bbox,scale)
   listPoints = buildList(chaine, mode)
   if(bbox == "bbox" ) then
      listPoints = buildBoundingBox(listPoints)
   end
   output = tracePointsMP(listPoints,points)
   output = "\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"
   tex.sprint(output)
end

-- print points to mesh
function printPointsTikZ(chaine,mode,points,bbox,scale,color,colorBbox)
   listPoints = buildList(chaine, mode)
   if(bbox == "bbox" ) then
      listPoints = buildBoundingBox(listPoints)
   end
   output = tracePointsTikZ(listPoints,points,color,colorBbox)
   output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" .. output .."\\end{tikzpicture}"
   tex.sprint(output)
end


-- print points to mesh
function printPointsTikZinc(chaine,beginning, ending, mode,points,bbox,scale,color,colorBbox)
   listPoints = buildList(chaine, mode)
   if(bbox == "bbox" ) then
      listPoints = buildBoundingBox(listPoints)
   end
   output = tracePointsTikZ(listPoints,points,color,colorBbox)
   output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" ..beginning.. output..ending .."\\end{tikzpicture}"
   tex.sprint(output)
end


-- buildMesh
function buildRect(largeur,a,b,nbrA, nbrB)
   listPoints = rectangleList(a,b,nbrA,nbrB)
   triangulation = BowyerWatson(listPoints,"none")
   traceTikZ(listPoints, triangulation,largeur,"none")
end

local function shallowCopy(original)
   local copy = {}
   for key, value in pairs(original) do
      copy[key] = value
   end
   return copy
end

-- function give a real polygon without repeting points
function cleanPoly(polygon)
   polyNew = {}
   polyCopy = shallowCopy(polygon)
   e1 = polyCopy[1][1]
   e2 = polyCopy[1][2]
   table.insert(polyNew, e1)
   table.insert(polyNew, e2)
   table.remove(polyCopy,1)
   j = 2
   while #polyCopy>1 do
      i=1
      find = false
      while (i<=#polyCopy and find==false) do
         bool1 = (polyCopy[i][1] == polyNew[j])
         bool2 = (polyCopy[i][2] == polyNew[j])
         if(bool1 or bool2) then -- the edge has a common point with polyNew[j]
            if(not bool1) then
               table.insert(polyNew, polyCopy[i][1])
               find = true
               table.remove(polyCopy,i)
               j = j+1
            elseif(not bool2) then
               table.insert(polyNew, polyCopy[i][2])
               find = true
               table.remove(polyCopy,i)
               j = j+1
            end
         end
         i=i+1
      end
   end
   return polyNew
end

--
function TeXaddOnePointTikZ(chaine,point,step,color,colorBack, colorNew, colorCircle,colorBbox)
   output = ""
   -- build the triangulation
   triangulation = BowyerWatson(listPoints,"none")
   badTriangles = buildBadTriangles(P,triangulation)
   if(step == "badT") then
      -- draw all triangle
      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
      -- draw and fill the bad triangle
      for i=1,#badTriangles do
         PointI = listPoints[triangulation[badTriangles[i]][1]]
         PointJ = listPoints[triangulation[badTriangles[i]][2]]
         PointK = listPoints[triangulation[badTriangles[i]][3]]
         output = output .. "\\draw[fill="..colorBack.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
      end
      -- draw the circoncircle
      for i=1,#badTriangles do
         PointI = listPoints[triangulation[badTriangles[i]][1]]
         PointJ = listPoints[triangulation[badTriangles[i]][2]]
         PointK = listPoints[triangulation[badTriangles[i]][3]]
         center, radius = circoncircle(PointI, PointJ, PointK)
         output = output .. "\\draw[dashed, color="..colorCircle.."] ("..center.x .. "," .. center.y .. ") circle ("..radius ..");"
      end
      -- mark the points
      for i=1,#listPoints do
         if(listPoints[i].type == "bbox") then
            output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
         else
            output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
         end
      end
      -- mark the point to add
      output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
   elseif(step == "cavity") then
      polygon = buildCavity(badTriangles, triangulation)
      polyNew = cleanPoly(polygon)
      -- remove the bad triangles
      for j=1,#badTriangles do
         table.remove(triangulation,badTriangles[j]-(j-1))
      end
      -- draw the triangles
      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
      -- fill and draw the cavity
      path = ""
      for i=1,#polyNew do
         PointI = listPoints[polyNew[i]]
         path = path .. "(".. PointI.x ..",".. PointI.y ..")--"
      end
      output = output .. "\\draw[color="..colorNew..",fill ="..colorBack..", thick] " .. path .. "cycle;"
      -- mark the points of the mesh
      for i=1,#listPoints do
         if(listPoints[i].type == "bbox") then
            output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
         else
            output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
         end
      end
      -- mark the adding point
      output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
   elseif(step == "newT") then
      polygon = buildCavity(badTriangles, triangulation)
      polyNew = cleanPoly(polygon)
      -- remove the bad triangles
      for j=1,#badTriangles do
         table.remove(triangulation,badTriangles[j]-(j-1))
      end
      -- draw the triangle of the triangulation
      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
      -- fill and draw the cavity
      path = ""
      for i=1,#polyNew do
         PointI = listPoints[polyNew[i]]
         path = path .. "(".. PointI.x ..",".. PointI.y ..")--"
      end
      output = output .. "\\draw[color="..colorNew..",fill ="..colorBack..", thick] " .. path .. "cycle;"
      -- draw the new triangles composed by the edges of the polygon and the added point
      for i=1,#polygon do
         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 ..");"
         output = output .. "\\draw[color="..colorNew..", thick]".."(".. listPoints[polygon[i][1]].x .. "," .. listPoints[polygon[i][1]].y .. ") -- (" .. P.x .. "," .. P.y ..");"
         output = output .. "\\draw[color="..colorNew..", thick]".."(".. listPoints[polygon[i][2]].x .. "," .. listPoints[polygon[i][2]].y .. ") -- (" .. P.x .. "," .. P.y ..");"
      end
      -- mark points
      for i=1,#listPoints do
         if(listPoints[i].type == "bbox") then
            output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
         else
            output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
         end
      end
      -- mark the added point
      output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
   end
   return output
end

function TeXaddOnePointMPBW(listPoints,P,step,bbox)
   output = "";
   output = output .. "pair MeshPoints[];"
   -- build the triangulation
   triangulation = BowyerWatson(listPoints,bbox)
   badTriangles = buildBadTriangles(P,triangulation)
   for i=1,#listPoints do
      output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
   end
   if(step == "badT") then
      -- draw all triangle
      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
      -- draw and fill the bad triangle
      for i=1,#badTriangles do
         PointI = listPoints[triangulation[badTriangles[i]][1]]
         PointJ = listPoints[triangulation[badTriangles[i]][2]]
         PointK = listPoints[triangulation[badTriangles[i]][3]]
         output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
         output = output .. "fill (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBack;"
      end
      -- draw the circoncircle
      for i=1,#badTriangles do
         PointI = listPoints[triangulation[badTriangles[i]][1]]
         PointJ = listPoints[triangulation[badTriangles[i]][2]]
         PointK = listPoints[triangulation[badTriangles[i]][3]]
         center, radius = circoncircle(PointI, PointJ, PointK)
         output = output .. "draw fullcircle scaled ("..radius .."*2u) shifted ("..center.x .. "*u," .. center.y .. "*u) dashed evenly withcolor \\luameshmpcolorCircle;"
      end
      -- mark the points
      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
      -- mark the point to add
      output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\luameshmpcolorNew;"
   elseif(step == "cavity") then
      polygon = buildCavity(badTriangles, triangulation)
      polyNew = cleanPoly(polygon)
      -- remove the bad triangles
      for j=1,#badTriangles do
         table.remove(triangulation,badTriangles[j]-(j-1))
      end
      -- draw the triangles
      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
      -- fill and draw the cavity
      path = ""
      for i=1,#polyNew do
         PointI = listPoints[polyNew[i]]
         path = path .. "(".. PointI.x ..",".. PointI.y ..")*u--"
      end
      output = output .. "fill " .. path .. "cycle withcolor \\luameshmpcolorBack;"
      output = output .. "draw " .. path .. "cycle withcolor \\luameshmpcolorNew  withpen pencircle scaled 1pt;"
      -- mark the points of the mesh
      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
      -- mark the adding point
      output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\luameshmpcolorNew ;"
   elseif(step == "newT") then
      polygon = buildCavity(badTriangles, triangulation)
      polyNew = cleanPoly(polygon)
      -- remove the bad triangles
      for j=1,#badTriangles do
         table.remove(triangulation,badTriangles[j]-(j-1))
      end
      -- draw the triangle of the triangulation
      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
      -- fill  the cavity
      path = ""
      for i=1,#polyNew do
         PointI = listPoints[polyNew[i]]
         path = path .. "(".. PointI.x ..",".. PointI.y ..")*u--"
      end
      output = output .. "fill " .. path .. "cycle withcolor \\luameshmpcolorBack;"
      -- draw the new triangles composed by the edges of the polygon and the added point
      for i=1,#polygon do
         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;"
         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;"
         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;"
      end
      -- mark points
      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
      -- mark the added point
      output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\luameshmpcolorNew ;"
   end
   return output
end


-- build the list of points extern and stop at nbr
function buildListExt(chaine, stop)
   listPoints = {}
   io.input(chaine) -- open the file
   text=io.read("*all")
   lines=string.explode(text,"\n+") -- all the lines
   for i=1,tonumber(stop) do
      xy=string.explode(lines[i]," +")
      table.insert(listPoints,{x=tonumber(xy[1]),y=tonumber(xy[2])})
   end
   xy=string.explode(lines[stop+1]," +")
   point={x=tonumber(xy[1]),y=tonumber(xy[2])}
   return point, listPoints
end


function TeXOnePointTikZBW(chaine,point,step,scale,mode,bbox,color,colorBack,colorNew,colorCircle,colorBbox)
   if(mode=="int") then
      Sx,Sy=string.match(point,"%((.+),(.+)%)")
      P = {x=Sx, y=Sy}
      listPoints = buildList(chaine, mode)
   else
      -- point is a number
      P, listPoints = buildListExt(chaine,tonumber(point))
   end
   output = TeXaddOnePointTikZ(chaine,point,step,color,colorBack,colorNew,colorCircle,colorBbox)
   output = "\\noindent\\begin{tikzpicture}[x="..scale..",y="..scale.."]".. output .. "\\end{tikzpicture}"
   tex.sprint(output)
end

function TeXOnePointTikZBWinc(chaine,point,beginning, ending,step,scale,mode,bbox,color,colorBack,colorNew,colorCircle,colorBbox)
   if(mode=="int") then
      Sx,Sy=string.match(point,"%((.+),(.+)%)")
      P = {x=Sx, y=Sy}
      listPoints = buildList(chaine, mode)
   else
      -- point is a number
      P, listPoints = buildListExt(chaine,tonumber(point))
   end
   output = TeXaddOnePointTikZ(chaine,point,step,color,colorBack,colorNew,colorCircle,colorBbox)
   output = "\\noindent\\begin{tikzpicture}[x="..scale..",y="..scale.."]".. beginning..output ..ending.. "\\end{tikzpicture}"
   tex.sprint(output)
end

function TeXOnePointMPBW(chaine,point,step,scale,mode,bbox)
   if(mode=="int") then
      Sx,Sy=string.match(point,"%((.+),(.+)%)")
      P = {x=Sx, y=Sy}
      listPoints = buildList(chaine, mode)
   else
      -- point is a number
      P, listPoints = buildListExt(chaine,tonumber(point))
   end
   output = TeXaddOnePointMPBW(listPoints,P,step,bbox)
   output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale..";".. output .. "endfig;\\end{mplibcode}"
   tex.sprint(output)
end

function TeXOnePointMPBWinc(chaine,point,beginning,ending,step,scale,mode,bbox)
   if(mode=="int") then
      Sx,Sy=string.match(point,"%((.+),(.+)%)")
      P = {x=Sx, y=Sy}
      listPoints = buildList(chaine, mode)
   else
      -- point is a number
      P, listPoints = buildListExt(chaine,tonumber(point))
   end
   output = TeXaddOnePointMPBW(listPoints,P,step,bbox)
   output = "\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"
   tex.sprint(output)
end