diff options
Diffstat (limited to 'Master/texmf-dist/tex/lualatex/bezierplot')
| -rwxr-xr-x | Master/texmf-dist/tex/lualatex/bezierplot/bezierplot.lua | 756 | ||||
| -rw-r--r-- | Master/texmf-dist/tex/lualatex/bezierplot/bezierplot.sty | 10 |
2 files changed, 430 insertions, 336 deletions
diff --git a/Master/texmf-dist/tex/lualatex/bezierplot/bezierplot.lua b/Master/texmf-dist/tex/lualatex/bezierplot/bezierplot.lua index 8cc414f278d..12574c6eef0 100755 --- a/Master/texmf-dist/tex/lualatex/bezierplot/bezierplot.lua +++ b/Master/texmf-dist/tex/lualatex/bezierplot/bezierplot.lua @@ -1,6 +1,6 @@ #!/usr/bin/env lua -- Linus Romer, published 2018 under LPPL Version 1.3c --- version 1.1 2018-06-10 +-- version 1.2 2018-07-23 abs = math.abs acos = math.acos asin = math.asin @@ -13,6 +13,7 @@ pi = math.pi sin = math.sin sqrt = math.sqrt tan = math.tan +huge = math.huge -- cube root defined for all real numbers x function cbrt(x) @@ -42,15 +43,30 @@ local function round(num, decimals) end end +-- check if a point (x,y) satisfies xmin <= x <= xmax and ymin <= < <= ymax +local function is_in_window(x,y,xmin,xmax,ymin,ymax) + if x >= xmin and x <= xmax and y >= ymin and y <= ymax then + return true + else + return false + end +end + +local function evaluate(s) + local tempfunc = assert(load("return " .. s)) + return tempfunc() +end + + -- 5-stencil method -- return from a graph from f in the form {{x,y},...} --- the derivatives in form {{x,y,dy/dx,ddy/ddx},...} +-- the derivatives in form {{x,y,dy/dx,ddy/ddx,extrema,inflection},...} local function diffgraph(func,graph,h) local dgraph = {} local l = #graph if l < 4 then -- this is not worth the pain... for i = 1, l do - table.insert(dgraph,{graph[i][1],graph[i][2],0,0}) + table.insert(dgraph,{graph[i][1],graph[i][2],0,0,0,0}) end else local yh = func(graph[1][1]-h) @@ -149,6 +165,68 @@ local function diffgraph(func,graph,h) return dgraph end +-- simplified diffgraph function, the function is derived only once +-- return from a graph from f in the form {{x,y},...} +-- the derivatives in form {{x,y,dy/dx},...} +-- we start with index 1 and then always jump indexjump to the next +-- index +local function diffgraphsimple(func,graph,h,indexjump) + local dgraph = {} + local l = math.floor((#graph-1)/indexjump)*indexjump+1 + if l < 2 then -- this is not worth the pain... + dgraph = {graph[1][1],graph[1][2],0} + else + local yh = func(graph[1][1]-h) + local yhh = func(graph[1][1]-2*h) + if yhh > -math.huge and yhh < math.huge -- if defined at all + and yh > -math.huge and yh < math.huge then + dgraph[1] = {graph[1][1],graph[1][2], + (yhh-8*yh+8*graph[2][2]-graph[3][2])/(12*h)} + else -- take neighbour values + dgraph[1] = {graph[1][1],graph[1][2], + (graph[1][2]-8*graph[2][2]+8*graph[4][2]-graph[5][2]) + /(12*h)} + end + for i = 1+indexjump, l-1, indexjump do + table.insert(dgraph,{graph[i][1],graph[i][2], + (graph[i-2][2]-8*graph[i-1][2]+8*graph[i+1][2]-graph[i+2][2]) + /(12*h)}) + end + yh = func(graph[l][1]+h) + yhh = func(graph[l][1]+2*h) + if yhh > -math.huge and yhh < math.huge -- if defined at all + and yh > -math.huge and yh < math.huge then + + table.insert(dgraph,{graph[l][1],graph[l][2], + (graph[l-2][2]-8*graph[l-1][2]+8*yh-yhh)/(12*h)}) + else + -- take neighbour values + table.insert(dgraph,{graph[l][1],graph[l][2], + (graph[l-4][2]-8*graph[l-3][2]+8*graph[l-1][2]-graph[l][2]) + /(12*h)}) + end + end + return dgraph +end + +-- diffgraph for cubic function ax^3+bx^2+cx+d +-- returns the derivatives in form {{x,y,dy/dx,ddy/ddx},...} +-- if isinverse = true then the coordinates will be inversed +local function diffgraphcubic(graph,a,b,c,d,isinverse) + local dgraph = {} + local l = #graph + for i = 1, l do + if isinverse then + dgraph[#dgraph+1] = {graph[i][2],graph[i][1],c + +graph[i][2]*(2*b+3*a*graph[i][2]),6*a*graph[i][2]+2*b} + else + dgraph[#dgraph+1] = {graph[i][1],graph[i][2],c + +graph[i][1]*(2*b+3*a*graph[i][1]),6*a*graph[i][1]+2*b} + end + end + return dgraph +end + -- checks for 100 x, if the function given by funcstring -- fits the graph g (up to maxerror) after filling in -- the parameters a, b, c, d @@ -283,94 +361,6 @@ local function parameters_affine(xp,yp,xq,yq) return a, b end --- returns true iff the function is of type f(x)=a*x^3+b*x^2+c*x+d --- a, b, c, d being real numbers -local function is_cubic(graph,maxerror) - local l = #graph - if l < 2 then - return false - else - local a, b, c, d = parameters_cubic(graph[1][1],graph[1][2], - graph[math.floor(l/3)][1],graph[math.floor(l/3)][2], - graph[math.floor(2*l/3)][1],graph[math.floor(2*l/3)][2], - graph[l][1],graph[l][2]) - return do_parameters_fit(a,b,c,d,"a*x^3+b*x^2+c*x+d",graph, - maxerror,false) - end -end - --- returns true iff the function is of type f(x)=a*x^3+b*x^2+c*x+d --- a, b, c, d being real numbers --- this takes several graph parts --- the idea is to have a possibility to avoid tan(x) -local function are_cubic(graphs,maxerror) - if is_cubic(graphs[1],maxerror) then - if #graphs < 2 then - return true - else -- check for the next part - local a, b, c, d = parameters_cubic(graphs[1][1][1], - graphs[1][1][2],graphs[1][math.floor(l/3)][1], - graphs[1][math.floor(l/3)][2], - graphs[1][math.floor(2*l/3)][1], - graphs[1][math.floor(2*l/3)][2], - graphs[1][l][1],graphs[1][l][2]) - return do_parameters_fit(a,b,c,d,"a*x^3+b*x^2+c*x+d", - graphs[2],maxerror,false) - end - else - return false - end -end - --- returns true iff the inverse function is of type --- f(x)=a*x^3+b*x^2+c*x+d --- a, b, c, d being real numbers -local function is_cuberoot(graph,maxerror) - local l = #graph - if l < 2 then - return false - else - local a, b, c, d = parameters_cubic(graph[1][2],graph[1][1], - graph[math.floor(l/3)][2],graph[math.floor(l/3)][1], - graph[math.floor(2*l/3)][2],graph[math.floor(2*l/3)][1], - graph[l][2],graph[l][1]) - return do_parameters_fit(a,b,c,d,"a*x^3+b*x^2+c*x+d",graph, - maxerror,true) - end -end - --- returns true iff the function is of type f(x)=a*x^3+b*x^2+c*x+d --- a, b, c, d being real numbers --- this takes several graph parts --- the idea is to have a possibility to avoid tan(x) -local function are_cuberoot(graphs,maxerror) - if is_cuberoot(graphs[1],maxerror) then - if #graphs < 2 then - return true - else -- check for the next part - local a, b, c, d = parameters_cubic(graphs[1][1][2], - graphs[1][1][1],graphs[1][math.floor(l/3)][2], - graphs[1][math.floor(l/3)][1], - graphs[1][math.floor(2*l/3)][2], - graphs[1][math.floor(2*l/3)][1], - graphs[1][l][2],graphs[1][l][1]) - return do_parameters_fit(a,b,c,d,"a*x^3+b*x^2+c*x+d", - graphs[2],maxerror,true) - end - else - return false - end -end - --- returns true iff function is of type f(x)=a*x+b --- a, b being real numbers -local function is_affine(graph,maxerror) - l = #graph - local a, b = parameters_affine(graph[1][1],graph[1][2], - graph[l][1],graph[l][2]) - return do_parameters_fit(a,b,0,0,"a*x+b",graph,maxerror,false) -end - -- what is the sum of the squared error -- when comparing the bezier path -- p.. control q and r .. s @@ -378,7 +368,7 @@ end -- (looking at the points at roughly t=.33 and t=.67) local function squareerror(f,g,starti,endi,qx,qy,rx,ry) local result = 0 - for t = .33, .7, .34 do + for t = .1, .9, .1 do x = (1-t)^3*g[starti][1]+3*t*(1-t)^2*qx+3*t^2*(1-t)*rx+t^3*g[endi][1] y = (1-t)^3*g[starti][2]+3*t*(1-t)^2*qy+3*t^2*(1-t)*ry+t^3*g[endi][2] result = result + (y-f(x))^2 @@ -389,7 +379,7 @@ end -- converts a table with bezier point information -- to a string with rounded values -- the path is reversed, if rev is true --- e.g. if b = {{0,1},{2,3,4,5,6,7},{8,9,10,11,12,13}} +-- e.g. if bezierpoints = {{0,1},{2,3,4,5,6,7},{8,9,10,11,12,13}} -- then -- (0,1) .. controls (2,3) and (4,5) .. (6,7) .. controls -- (8,9) and (10,11) .. (12,13) @@ -399,50 +389,68 @@ end -- NO: 0 1 \\ 6 7 \\ 2 3 \\ 4 5 \\ \\ 6 7 \\ 12 13 \\ 8 9 \\ 10 11 \\ -- As pgfplots does not connect the bezier segments -- reverse paths are not implemented -local function beziertabletostring(b,rndx,rndy,rev,notation) +local function beziertabletostring(bezierpoints,rndx,rndy,rev,notation) local bezierstring = "" + local b = {{round(bezierpoints[1][1],rndx),round(bezierpoints[1][2],rndy)}} -- rounded and then + -- reduced points (if identical after rounding) + -- rounding + for i = 2, #bezierpoints do + -- check if x--coordinates are identical + if round(bezierpoints[i][#bezierpoints[i]-1],rndx) ~= b[#b][#b[#b]-1] then + b[#b+1] = {} + for j = 1, #bezierpoints[i] do + if j % 2 == 0 then -- x coordinate + b[#b][j] = round(bezierpoints[i][j],rndx) + else + b[#b][j] = round(bezierpoints[i][j],rndy) + end + end + end + end if #b > 1 then -- if not empty or single point - if #b == 2 and #b[2] == 2 then -- straight line + -- check if bezierstring contains only straight lines + local onlystraightlines = true + for i = 1, #b do + if #b[i] > 2 then + onlystraightlines = false + break + end + end + if onlystraightlines then if rev then - bezierstring = "(" .. round(b[2][1],rndx) .. "," - .. round(b[2][2],rndy) ..")" - .. " -- (" .. round(b[1][1],rndx) .. "," - .. round(b[1][2],rndy) ..")" + bezierstring = "(" .. b[#b][1] .. "," .. b[#b][2] ..")" + for i = #b-1, 1, -1 do + bezierstring = bezierstring .. + " -- (" .. b[i][1] .. "," .. b[i][2] ..")" + end else if notation == "pgfplots" then bezierstring = "\\addplot coordinates {(" - .. round(b[1][1],rndx) .. "," - .. round(b[1][2],rndy) .. ") (" - .. round(b[2][1],rndx) .. "," - .. round(b[2][2],rndy) .. ") (" - .. round(b[1][1],rndx) .. "," - .. round(b[1][2],rndy) .. ") (" - .. round(b[2][1],rndx) .. "," - .. round(b[2][2],rndy) .. ") }" + .. b[1][1] .. "," .. b[1][2] .. ") (" + .. b[2][1] .. "," .. b[2][2] .. ") (" + .. b[1][1] .. "," .. b[1][2] .. ") (" + .. b[2][1] .. "," .. b[2][2] .. ") }" else -- notation = tikz - bezierstring = "(" .. round(b[1][1],rndx) .. "," - .. round(b[1][2],rndy) ..")" - .. " -- (" .. round(b[2][1],rndx) .. "," - .. round(b[2][2],rndy) ..")" + bezierstring = "(" .. b[1][1] .. "," .. b[1][2] ..")" + for i = 2, #b do + bezierstring = bezierstring .. + " -- (" .. b[i][1] .. "," .. b[i][2] ..")" + end end end else if rev then - bezierstring = "(" .. round(b[#b][#b[#b]-1],rndx) .. "," - .. round(b[#b][#b[#b]],rndy) ..")" -- initial point + bezierstring = "(" .. b[#b][#b[#b]-1] .. "," + .. b[#b][#b[#b]] ..")" -- initial point for i = #b, 2, -1 do if #b[i] >= 6 then -- cubic bezier spline bezierstring = bezierstring .. " .. controls (" - .. round(b[i][3],rndx) .. "," - .. round(b[i][4],rndy) ..") and (" - .. round(b[i][1],rndx) .. "," - .. round(b[i][2],rndy) .. ") .. (" - .. round(b[i-1][#b[i-1]-1],rndx) .. "," - .. round(b[i-1][#b[i-1]],rndy)..")" + .. b[i][3] .. "," .. b[i][4] ..") and (" + .. b[i][1] .. "," .. b[i][2] .. ") .. (" + .. b[i-1][#b[i-1]-1] .. "," .. b[i-1][#b[i-1]]..")" else bezierstring = bezierstring .. " (" - .. round(b[i-1][#b[i-1]-1],rndx) .. "," - .. round(b[i-1][#b[i-1]],rndy) ..")" + .. b[i-1][#b[i-1]-1] .. "," .. b[i-1][#b[i-1]] ..")" end end else @@ -451,33 +459,25 @@ local function beziertabletostring(b,rndx,rndy,rev,notation) for i = 1, #b-1 do if #b[i+1] >= 6 then -- cubic bezier spline bezierstring = bezierstring .. "(" - .. round(b[i][#b[i]-1],rndx) .. "," - .. round(b[i][#b[i]],rndy) .. ") (" - .. round(b[i+1][5],rndx) .. "," - .. round(b[i+1][6],rndy) .. ") (" - .. round(b[i+1][1],rndx) .. "," - .. round(b[i+1][2],rndy) .. ") (" - .. round(b[i+1][3],rndx) .. "," - .. round(b[i+1][4],rndy) .. ") " + .. b[i][#b[i]-1] .. "," .. b[i][#b[i]] .. ") (" + .. b[i+1][5] .. "," .. b[i+1][6] .. ") (" + .. b[i+1][1] .. "," .. b[i+1][2] .. ") (" + .. b[i+1][3] .. "," .. b[i+1][4] .. ") " end end bezierstring = bezierstring .. "}" else -- notation = tikz - bezierstring = "(" .. round(b[1][1],rndx) .. "," - .. round(b[1][2],rndy) ..")" -- initial point + bezierstring = "(" .. b[1][1] .. "," + .. b[1][2] ..")" -- initial point for i = 2, #b do if #b[i] >= 6 then -- cubic bezier spline bezierstring = bezierstring .. " .. controls (" - .. round(b[i][1],rndx) .. "," - .. round(b[i][2],rndy) ..") and (" - .. round(b[i][3],rndx) .. "," - .. round(b[i][4],rndy) .. ") .. (" - .. round(b[i][5],rndx) .. "," - .. round(b[i][6],rndy)..")" + .. b[i][1] .. "," .. b[i][2] ..") and (" + .. b[i][3] .. "," .. b[i][4] .. ") .. (" + .. b[i][5] .. "," .. b[i][6]..")" else bezierstring = bezierstring .. " (" - .. round(b[i][1],rndx) .. "," - .. round(b[i][2],rndy) ..")" + .. b[i][1] .. "," .. b[i][2] ..")" end end end @@ -492,6 +492,7 @@ end -- without extrema or inflection points inbetween -- and try to approximate it with a cubic bezier curve -- (round to rndx and rndy when printing) +-- if maxerror <= 0, the function will not be recursive anymore local function graphtobezierapprox(f,g,starti,endi,maxerror) local px = g[starti][1] local py = g[starti][2] @@ -504,59 +505,64 @@ local function graphtobezierapprox(f,g,starti,endi,maxerror) local cy = (dp * ((ds * px) - (ds * sx) - py + sy) / (dp - ds)) + py -- now we slide q between p and c & r between s and c -- and search for the best qx and best rx - local qx = px+.05*(cx-px) - local qy = py+.05*(cy-py) - local rx = sx+.05*(cx-sx) - local ry = sy+.05*(cy-sy) + local qx = px+.01*(cx-px) + local qy = py+.01*(cy-py) + local rx = sx+.01*(cx-sx) + local ry = sy+.01*(cy-sy) local err = squareerror(f,g,starti,endi,qx,qy,rx,ry) - for i = 2, 19 do - for j = 2, 19 do - xa = px+i*.05*(cx-px) - ya = py+i*.05*(cy-py) - xb = sx+j*.05*(cx-sx) - yb = sy+j*.05*(cy-sy) + for i = 2, 99 do + for j = 2, 99 do + xa = px+i*.01*(cx-px) + ya = py+i*.01*(cy-py) + xb = sx+j*.01*(cx-sx) + yb = sy+j*.01*(cy-sy) -- now check, if xa and xb fit better - -- at roughly t=0.33 and t=0.66 for f(x) -- than the last qx and rx did -- (sum of squares must be smaller) - if squareerror(f,g,starti,endi,xa,ya,xb,yb) < err then + local newerror = squareerror(f,g,starti,endi,xa,ya,xb,yb) + if newerror < err then qx = xa qy = ya rx = xb ry = yb - err = squareerror(f,g,starti,endi,qx,qy,rx,ry) + err = newerror end end end - -- check if it is close enough: (recycling err, xa, ya) - err = 0 - for t = .1, .9, .1 do - xa = (1-t)^3*g[starti][1]+3*t*(1-t)^2*qx+3*t^2*(1-t)*rx+t^3*g[endi][1] - ya = (1-t)^3*g[starti][2]+3*t*(1-t)^2*qy+3*t^2*(1-t)*ry+t^3*g[endi][2] - if abs(ya-f(xa)) > err then - err = abs(ya-f(xa)) - end - end - if err <= maxerror then - return {qx,qy,rx,ry,sx,sy} - else - -- search for an intermediate point where the graph has the same - -- slope as the line from the start point to the end point: - local interindex = math.floor(.5*starti+.5*endi) -- will change - for i = starti + 1, endi - 1 do - if abs(g[i][3]-(g[endi][2]-g[starti][2]) - /(g[endi][1]-g[starti][1])) - < abs(g[interindex][3]-(g[endi][2]-g[starti][2]) - /(g[endi][1]-g[starti][1])) then - interindex = i + if maxerror > 0 then + -- check if it is close enough: (recycling err, xa, ya) + err = 0 + for t = .1, .9, .1 do + xa = (1-t)^3*g[starti][1]+3*t*(1-t)^2*qx+3*t^2*(1-t)*rx+t^3*g[endi][1] + ya = (1-t)^3*g[starti][2]+3*t*(1-t)^2*qy+3*t^2*(1-t)*ry+t^3*g[endi][2] + if abs(ya-f(xa)) > err then + err = abs(ya-f(xa)) + err = abs(ya-f(xa)) end end - local left = graphtobezierapprox(f,g,starti,interindex,maxerror) - local right = graphtobezierapprox(f,g,interindex,endi,maxerror) - for i=1, #right do --now append the right to the left: - left[#left+1] = right[i] + if err <= maxerror then + return {qx,qy,rx,ry,sx,sy} + else + -- search for an intermediate point where the graph has the same + -- slope as the line from the start point to the end point: + local interindex = math.floor(.5*starti+.5*endi) -- will change + for i = starti + 1, endi - 1 do + if abs(g[i][3]-(g[endi][2]-g[starti][2]) + /(g[endi][1]-g[starti][1])) + < abs(g[interindex][3]-(g[endi][2]-g[starti][2]) + /(g[endi][1]-g[starti][1])) then + interindex = i + end + end + local left = graphtobezierapprox(f,g,starti,interindex,maxerror) + local right = graphtobezierapprox(f,g,interindex,endi,maxerror) + for i=1, #right do --now append the right to the left: + left[#left+1] = right[i] + end + return left end - return left + else + return {qx,qy,rx,ry,sx,sy} end end @@ -591,31 +597,58 @@ local function printtable(t) end -- main function -function bezierplot(functionstring,xmin,xmax,ymin,ymax,notation) +function bezierplot(functionstring,xminstring,xmaxstring,yminstring,ymaxstring,samplesstring,notation) local fstringreplaced = string.gsub(functionstring, "%*%*", "^") local f = assert(load("local x = ...; return " .. fstringreplaced)) + local xmin = evaluate(xminstring) + local xmax = evaluate(xmaxstring) + local ymin = evaluate(yminstring) + local ymax = evaluate(ymaxstring) + local samples = evaluate(samplesstring) local isreverse = false if xmin > xmax then isreverse = true + elseif xmin == xmax then + xmax = xmin + 10 end xmin, xmax = math.min(xmin,xmax), math.max(xmin,xmax) - local xstep = (xmax-xmin)/20000 - -- the output of the x coordinates will be rounded to rndx digits - local rndx = math.max(0,math.floor(4.5-log(xmax-xmin)/log(10))) - local xerror = abs(xmax-xmin)/(100*10^rndx) + if ymin == ymax then + ymax = ymin + 10 + end ymin, ymax = math.min(ymin,ymax), math.max(ymin,ymax) - -- the output of the x coordinates will be rounded to rndy digits - local rndy = math.max(0,math.floor(4.5-log(ymax-ymin)/log(10))) - local yerror = (ymax-ymin)/(100*10^rndy) + local xsteps = 50000 + -- if samples < 2 the samples will be chosen as wisely as possible + local arbitrary_samples = true + if samples >= 2 then + arbitrary_samples = false + xsteps = (samples-1)*math.max(2,math.floor(xsteps/(samples-1))) + end + local xstep = (xmax-xmin)/xsteps + -- the output of the x coordinates will be rounded to rndx digits + local rndx = math.max(0,math.floor(5.5-log(xmax-xmin)/log(10))) + local xerror = abs(xmax-xmin)/(10^rndx) + -- the output of the y coordinates will be rounded to rndy digits + local rndy = math.max(0,math.floor(5.5-log(ymax-ymin)/log(10))) + local yerror = (ymax-ymin)/(10^rndy) -- determine parts of the graph that are inside window - local graphs = {} + local graphs = {} -- graph split to the connected parts + local graph = {} -- graphs concatenated (needed for function type) local outside = true -- value is outside window local i = 0 local j = 0 - for n = 0, 20000 do - local x = xmin + n/20000*(xmax-xmin) + local yminreal -- determine the real minimimum of the y coord. + local ymaxreal -- just decring + local yminrealfound = false + local ymaxrealfound = false + for n = 0, xsteps do + local x = xmin + n/xsteps*(xmax-xmin) + if n == xsteps then + x = xmax + end local y = f(x) - if y >= ymin-yerror and y <= ymax+yerror then -- inside + if (y >= ymin-yerror and ymin ~= -huge or y > ymin and ymin == -huge) + and (y <= ymax+yerror and ymax ~= huge or y < ymax and ymax == huge) + then -- inside if outside then -- if it was outside before outside = false j = 0 @@ -624,203 +657,264 @@ function bezierplot(functionstring,xmin,xmax,ymin,ymax,notation) end j = j + 1 graphs[i][j] = {x,y} + graph[#graph+1] = {x,y} + if not yminrealfound or yminrealfound and y < yminreal then + yminreal = y + yminrealfound = true + end + if not ymaxrealfound or ymaxrealfound and y > ymaxreal then + ymaxreal = y + ymaxrealfound = true + end else outside = true end end - + + -- some redefinitions + if #graph ~= 0 and yminreal ~= ymaxreal then + ymin = yminreal + ymax = ymaxreal + rndy = math.max(0,math.floor(5.5-log(ymax-ymin)/log(10))) + yerror = (ymax-ymin)/(10^rndy) + end + + -- check for the function type (for this, we need the concatenated + -- parts of the graph) + -- go through the connected parts local functiontype = "unknown" + local a, b, c, d -- possible function parameter + -- check for affine functions: + local l = #graph + a, b = parameters_affine(graph[1][1],graph[1][2], + graph[l][1],graph[l][2]) + if do_parameters_fit(a,b,0,0,"a*x+b",graph,yerror,false) then + functiontype = "affine" + else -- check for cubic functions (includes quadratic functions) + a, b, c, d = parameters_cubic(graph[1][1],graph[1][2], + graph[math.floor(l/3)][1],graph[math.floor(l/3)][2], + graph[math.floor(2*l/3)][1],graph[math.floor(2*l/3)][2], + graph[l][1],graph[l][2]) + if do_parameters_fit(a,b,c,d,"a*x^3+b*x^2+c*x+d",graph, + yerror,false) then + functiontype = "cubic" + else -- check for cuberoot functions (includes squareroots) + a, b, c, d = parameters_cubic(graph[1][2],graph[1][1], + graph[math.floor(l/3)][2],graph[math.floor(l/3)][1], + graph[math.floor(2*l/3)][2],graph[math.floor(2*l/3)][1], + graph[l][2],graph[l][1]) + if do_parameters_fit(a,b,c,d,"a*x^3+b*x^2+c*x+d",graph, + xerror,true) then + functiontype = "cuberoot" + end + end + end + local bezierpoints = {} -- the bezier path (0,1) .. controls -- (2,3) and (4,5) .. (6,7) .. controls -- (8,9) and (10,11) .. (12,13) -- will be stored as -- bezierpoints={{0,1},{2,3,4,5,6,7},{8,9,10,11,12,13}} - - -- go through the connected parts - for part = 1, #graphs do - local d = diffgraph(f,graphs[part],xstep) - --for i = 1, #d do -- just for debugging - -- print(d[i][1],d[i][2],d[i][3],d[i][4],d[i][5],d[i][6]) - --end - -- check for type of function (only for the first part) - if part == 1 then - if is_affine(d,yerror) then - functiontype = "affine" - elseif are_cubic(graphs,yerror) then - functiontype = "cubic" - elseif are_cuberoot(graphs,xerror) then - functiontype = "cuberoot" + + if functiontype == "affine" then + if arbitrary_samples then + bezierpoints = {{graph[1][1],graph[1][2]},{graph[#graph][1], + graph[#graph][2]}} + else -- we can here savely assume that graphs has only one part, + -- therefore graphs[1]=graph + for i = 1, #graph, math.floor(xsteps/(samples-1)) do + bezierpoints[#bezierpoints+1] = {graph[i][1],graph[i][2]} end end - if functiontype ~= "cuberoot" then -- start with initial point - bezierpoints[#bezierpoints+1] = {round(d[1][1],rndx), - round(d[1][2],rndy)} - end - if functiontype == "affine" then - bezierpoints[#bezierpoints+1] = {round(d[#d][1],rndx), - round(d[#d][2],rndy)} - elseif functiontype == "cubic" then - local startindex = 1 - local extremainbetween = false - for k = 2, #d do - if d[k][5] then -- extrema - extremainbetween = true - bezierpoints[#bezierpoints+1] = graphtobezier(d, - startindex,k,false) - startindex = k + elseif functiontype == "cubic" then + local extrema_inflections = {} -- store the extrema and + -- inflection points for arbitrary samples + if arbitrary_samples then + if math.abs(a) < yerror*1e-10 then -- quadratic case (one extremum) + if is_in_window(-c/(2*b),(-c^2+4*b*d)/(4*b),xmin,xmax, + ymin,ymax) then + extrema_inflections = {{-c/(2*b),(-c^2+4*b*d)/(4*b)}} + end + else -- cubic case (two extrema and one inflection point) + -- we order the points with the help of sgn + -- check for first extrema + if is_in_window((-sgn(a)*sqrt(-3*a*c+b^2)-b)/(3*a), + (2*b^3+27*a^2*d-9*a*b*c+sqrt(b^2-3*a*c)*sgn(a)* + (2*b^2-6*a*c))/(27*a^2),xmin,xmax,ymin,ymax) then + extrema_inflections[#extrema_inflections+1] = + {(-sgn(a)*sqrt(-3*a*c+b^2)-b)/(3*a),(2*b^3+27*a^2*d- + 9*a*b*c+sqrt(b^2-3*a*c)*sgn(a)*(2*b^2-6*a*c))/(27*a^2)} + end + -- check for inflection point (has to be inbetween) + if is_in_window(-b/(3*a),(2*b^3+27*a^2*d-9*a*b*c) + /(27*a^2),xmin,xmax,ymin,ymax) then + extrema_inflections[#extrema_inflections+1]={-b/(3*a), + (2*b^3+27*a^2*d-9*a*b*c)/(27*a^2)} + end + -- check for second extrema + if is_in_window((sgn(a)*sqrt(-3*a*c+b^2)-b)/(3*a), + (2*b^3+27*a^2*d-9*a*b*c+sqrt(b^2-3*a*c)*sgn(a)* + (-2*b^2+6*a*c))/(27*a^2),xmin,xmax,ymin,ymax) then + extrema_inflections[#extrema_inflections+1] = + {(sgn(a)*sqrt(-3*a*c+b^2)-b)/(3*a),(2*b^3+27*a^2*d- + 9*a*b*c+sqrt(b^2-3*a*c)*sgn(a)*(-2*b^2+6*a*c))/(27*a^2)} end end - if not extremainbetween then - for k = 2, #d do - if d[k][6] then -- inflection point - -- check, if the controlpoints are outside - -- of the bounding box defined by the vertices - -- (d[1][1],d[1][2]) and (d[#d][1],d[#d][2]) - local qx = d[1][1]+(d[#d][1]-d[1][1])/3 - local rx = d[1][1]+2*(d[#d][1]-d[1][1])/3 - local qy = d[1][2]+(qx-d[1][1])*d[1][3] - local ry = d[#d][2]+(rx-d[#d][1])*d[#d][3] - if math.max(qy,ry) > ymax - or math.min(qy,ry) < ymin then - bezierpoints[#bezierpoints+1] = graphtobezier( - d,startindex,k,false) - startindex = k - end + end + for part = 1, #graphs do + bezierpoints[#bezierpoints+1] = {graphs[part][1][1], + graphs[part][1][2]} -- initial points + local graphsamples = {}-- will be the graph reduced to the + -- samples (or the most important points) + local dg -- will be the differentiated graph + if arbitrary_samples then -- add extrema and inflection + -- points to the border points + graphsamples = {{graphs[part][1][1], + graphs[part][1][2]}} + for j = 1, #extrema_inflections do + if extrema_inflections[j][1] > math.min( + graphs[part][1][1] ,graphs[part][#graphs[part]][1]) + and extrema_inflections[j][1] < math.max( + graphs[part][1][1] ,graphs[part][#graphs[part]][1]) + then + graphsamples[#graphsamples+1] = + {extrema_inflections[j][1], + extrema_inflections[j][2]} end end + graphsamples[#graphsamples+1] = + {graphs[part][#graphs[part]][1], + graphs[part][#graphs[part]][2]} + else + for i = 1, #graphs[part], xsteps/(samples-1) do + graphsamples[#graphsamples+1] = + {graphs[part][i][1],graphs[part][i][2]} + end end - if startindex ~= #d then -- if no special points inbetween - bezierpoints[#bezierpoints+1] = graphtobezier(d, - startindex,#d,false) + dg = diffgraphcubic(graphsamples,a,b,c,d,false) + for i = 2, #dg do + bezierpoints[#bezierpoints+1] = graphtobezier(dg,i-1,i,false) end - elseif functiontype == "cuberoot" then - -- we determine a, b, c, d and then - -- get x' = 3ay^2+2by+c - local a, b, c, dd = parameters_cubic( - d[math.floor(.2*l)][2], d[math.floor(.2*l)][1], - d[math.floor(.4*l)][2], d[math.floor(.4*l)][1], - d[math.floor(.6*l)][2], d[math.floor(.6*l)][1], - d[math.floor(.8*l)][2], d[math.floor(.8*l)][1]) - -- now recalculate the graph with the inverse function: - -- we can increase the accuracy - xstep = (ymax-ymin)/100000 -- inverse redefinition - local finverse = assert(load("local x = ...; return " - ..a.."*x^3+"..b.."*x^2+"..c.."*x+"..dd)) - local graphinverse = {} - local i = 1 - for y = ymin, ymax, xstep do - local x = finverse(y) - if x > xmin and x < xmax -- inside - and abs(y-f(x)) < (ymax-ymin)/(100*10^rndy) then - graphinverse[i] = {y,x} - i = i + 1 - end + end + elseif functiontype == "cuberoot" then + local inflection = {} -- store the inflection point + if arbitrary_samples and math.abs(a) ~= 0 + and is_in_window((2*b^3+27*a^2*d-9*a*b*c)/(27*a^2),-b/(3*a), + xmin,xmax,ymin,ymax) then + inflection = {(2*b^3+27*a^2*d-9*a*b*c)/(27*a^2),-b/(3*a)} + end + -- (there cannot be more than one part) + bezierpoints[#bezierpoints+1] = {graphs[1][1][1], + graphs[1][1][2]} -- initial points + local graphsamples = {}-- will be the graph reduced to the + -- samples (or the most important points) + local dg -- will be the differentiated graph + if arbitrary_samples then -- add inflection point (if exis.) + graphsamples = {{graphs[1][1][1], + graphs[1][1][2]}} + if #inflection > 0 and inflection[1] > math.min( + graphs[1][1][1],graphs[1][#graphs[1]][1]) + and inflection[1] < math.max( + graphs[1][1][1],graphs[1][#graphs[1]][1]) + then + graphsamples[#graphsamples+1] = + {inflection[1],inflection[2]} end - d = diffgraph(finverse,graphinverse,xstep) - bezierpoints[#bezierpoints+1] = {d[1][2],d[1][1]} -- initial point - local startindex = 1 - for k = 2, #d do - if d[k][6] then -- inflection point - -- check, if the controlpoints are outside - -- of the bounding box defined by the vertices - -- (d[1][1],d[1][2]) and (d[#d][1],d[#d][2]) - local qx = d[1][1]+(d[#d][1]-d[1][1])/3 - local rx = d[1][1]+2*(d[#d][1]-d[1][1])/3 - local qy = d[1][2]+(qx-d[1][1])*d[1][3] - local ry = d[#d][2]+(rx-d[#d][1])*d[#d][3] - if math.max(qy,ry) > xmax - or math.min(qy,ry) < xmin then - bezierpoints[#bezierpoints+1] = graphtobezier( - d,startindex,k,true) + graphsamples[#graphsamples+1] = + {graphs[1][#graphs[1]][1], + graphs[1][#graphs[1]][2]} + else + for i = 1, #graphs[1], xsteps/(samples-1) do + graphsamples[#graphsamples+1] = + {graphs[1][i][1],graphs[1][i][2]} + end + end + dg = diffgraphcubic(graphsamples,a,b,c,d,true) + for i = 2, #dg do + bezierpoints[#bezierpoints+1] = graphtobezier(dg,i-1,i,true) + end + else + ---------- generic case (no special function) ---------------- + if arbitrary_samples then + -- go through the connected parts + for part = 1, #graphs do + local dg = diffgraph(f,graphs[part],xstep) + bezierpoints[#bezierpoints+1] = {dg[1][1],dg[1][2]} + local startindex = 1 + for k = 2, #dg do + if dg[k][5] or dg[k][6] then -- extrema and inflection points + local tobeadded = graphtobezierapprox( + f,dg,startindex,k,10*yerror) + -- tobeadded may contain a multiple of 6 entries + -- e.g. {1,2,3,4,5,6,7,8,9,10,11,12} + for i = 1, math.floor(#tobeadded/6) do + bezierpoints[#bezierpoints+1] = {} + for j = 1, 6 do + bezierpoints[#bezierpoints][j] = tobeadded[(i-1)*6+j] + end + end startindex = k end end - end - if startindex ~= #d then -- if no special points inbetween - bezierpoints[#bezierpoints+1] = graphtobezier(d, - startindex,#d,true) - end - else - -- standard case (nothing special) - local startindex = 1 - for k = 2, #d do - if d[k][5] or d[k][6] then -- extrema and inflection points - local tobeadded = graphtobezierapprox( - f,d,startindex,k,(ymax-ymin)/(0.5*10^rndy)) + if startindex ~= #dg then -- if no special points inbetween + local tobeadded = graphtobezierapprox(f,dg, + startindex,#dg,10*yerror) -- tobeadded may contain a multiple of 6 entries - -- e.g. {1,2,3,4,5,6,7,8,9,10,11,12} + -- e.g. {1,2,3,4,5,6,7,8,9,10,11,12} for i = 1, math.floor(#tobeadded/6) do bezierpoints[#bezierpoints+1] = {} for j = 1, 6 do bezierpoints[#bezierpoints][j] = tobeadded[(i-1)*6+j] end end - startindex = k end end - if startindex ~= #d then -- if no special points inbetween - local tobeadded = graphtobezierapprox(f,d, - startindex,#d,(ymax-ymin)/(0.5*10^rndy)) - -- tobeadded may contain a multiple of 6 entries - -- e.g. {1,2,3,4,5,6,7,8,9,10,11,12} - for i = 1, math.floor(#tobeadded/6) do - bezierpoints[#bezierpoints+1] = {} - for j = 1, 6 do - bezierpoints[#bezierpoints][j] = tobeadded[(i-1)*6+j] - end + else -- fixed samples in the generic case + -- go through the connected parts + for part = 1, #graphs do + local dg = diffgraphsimple(f,graphs[part],xstep, + math.floor(0.5+xsteps/(samples-1))) + bezierpoints[#bezierpoints+1] = {dg[1][1],dg[1][2]} -- initial points + for i = 2, #dg do + bezierpoints[#bezierpoints+1] = graphtobezier(dg,i-1,i,false) end end end end - -- only for debugging: - --for i = 1, #bezierpoints do - -- for j = 1, #bezierpoints[i] do - -- print(bezierpoints[i][j]) - -- end - ---print("\n") - --end - return beziertabletostring(bezierpoints,rndx,rndy,isreverse,notation) + return beziertabletostring(bezierpoints,rndx,rndy,isreverse,notation) end -- main program -- if not pcall(debug.getlocal, 4, 1) then +--if debug.getinfo(3) == nil then if #arg >= 1 then local xmin = -5 local xmax = 5 if #arg >= 2 then - local tempfunc = assert(load("return " .. arg[2])) - xmin = tempfunc() + xmin = arg[2] end if #arg >= 3 then - if arg[3] == arg[2] then - xmax = xmin + 10 - else - local tempfunc = assert(load("return " .. arg[3])) - xmax = tempfunc() - end + xmax = arg[3] end local ymin = -5 local ymax = 5 if #arg >= 4 then - local tempfunc = assert(load("return " .. arg[4])) - ymin = tempfunc() + ymin = arg[4] end if #arg >= 5 then - if arg[5] == arg[4] then - ymax = ymin + 10 - else - local tempfunc = assert(load("return " .. arg[5])) - ymax = tempfunc() - end + ymax = arg[5] end + local samples = 0 if #arg >= 6 then - notation = arg[6] - else - notation = "tikz" + samples = arg[6] + end + local notation = "tikz" + if #arg >= 7 then + notation = arg[7] end - print(bezierplot(arg[1],xmin,xmax,ymin,ymax,notation)) + print(bezierplot(arg[1],xmin,xmax,ymin,ymax,samples,notation)) end end diff --git a/Master/texmf-dist/tex/lualatex/bezierplot/bezierplot.sty b/Master/texmf-dist/tex/lualatex/bezierplot/bezierplot.sty index 26e76f0e25b..be4c89550e5 100644 --- a/Master/texmf-dist/tex/lualatex/bezierplot/bezierplot.sty +++ b/Master/texmf-dist/tex/lualatex/bezierplot/bezierplot.sty @@ -1,16 +1,16 @@ \NeedsTeXFormat{LaTeX2e} -\ProvidesPackage{bezierplot}[2018/06/10 bezierplot] +\ProvidesPackage{bezierplot}[2018/07/23 bezierplot] \RequirePackage{xparse} \RequirePackage{iftex} \ifLuaTeX \directlua{require("bezierplot")} - \DeclareExpandableDocumentCommand{\xbezierplot}{O{-5} O{5} O{-5} O{5} O{tikz} m}{% - \directlua{tex.sprint(bezierplot("#6",#1,#2,#3,#4,"#5"))} + \DeclareExpandableDocumentCommand{\xbezierplot}{O{-5} O{5} O{-5} O{5} O{0} O{tikz} m}{% + \directlua{tex.sprint(bezierplot("#7","#1","#2","#3","#4","#5","#6"))} } \else \let\xpandblinpt\@@input - \DeclareExpandableDocumentCommand{\xbezierplot}{O{-5} O{5} O{-5} O{5} O{tikz} m}{% - \xpandblinpt|"bezierplot '#6' #1 #2 #3 #4 '#5'" + \DeclareExpandableDocumentCommand{\xbezierplot}{O{-5} O{5} O{-5} O{5} O{0} O{tikz} m}{% + \xpandblinpt|"bezierplot '#7' '#1' '#2' '#3' '#4' '#5' '#6'" } \fi \providecommand\bezierplot{\romannumeral`\^^@\xbezierplot} |