summaryrefslogtreecommitdiff
path: root/tex/context/base/mkiv/mlib-cnt.lmt
diff options
context:
space:
mode:
Diffstat (limited to 'tex/context/base/mkiv/mlib-cnt.lmt')
-rw-r--r--tex/context/base/mkiv/mlib-cnt.lmt2022
1 files changed, 2022 insertions, 0 deletions
diff --git a/tex/context/base/mkiv/mlib-cnt.lmt b/tex/context/base/mkiv/mlib-cnt.lmt
new file mode 100644
index 000000000..667384ed7
--- /dev/null
+++ b/tex/context/base/mkiv/mlib-cnt.lmt
@@ -0,0 +1,2022 @@
+if not modules then modules = { } end modules ['mlib-cnt'] = {
+ version = 1.001,
+ optimize = true,
+ comment = "companion to mlib-ctx.mkiv",
+ author = "Hans Hagen, PRAGMA-ADE, Hasselt NL",
+ copyright = "PRAGMA ADE / ConTeXt Development Team",
+ license = "see context related readme files",
+}
+
+-- The only useful reference that I could find about this topic is in wikipedia:
+--
+-- https://en.wikipedia.org/wiki/Marching_squares
+--
+-- I derived the edge code from:
+--
+-- https://physiology.arizona.edu/people/secomb/contours
+--
+-- and also here:
+--
+-- https://github.com/secomb/GreensV4
+--
+-- which has the banner:
+--
+-- TWS, November 1989. Converted to C September 2007. Revised February 2009.
+--
+-- and has a liberal licence. I optimized the code so that it runs a bit faster in Lua and
+-- there's probably more room for optimization (I tried several variants). For instance I
+-- don't use that many tables because access is costly. We don't have a compiler that does
+-- some optimizing (even then the c code can probably be made more efficient).
+--
+-- We often have the same node so it's cheaper to reuse the wsp tables and reconstruct
+-- the point in the path then to alias the original point. We can be more clever:
+-- straighten, but it's more work so maybe I will do that later; for now I only added
+-- a test for equality. There are some experiments in this file too and not all might
+-- work. It's a playground for me.
+--
+-- The code is meant for metafun so it is not general purpose. The bitmap variant is
+-- relatively fast and bitmaps compress well. When all is settled I might make a couple
+-- of helpers in C for this purpose but not now.
+--
+-- I removed optimization code (more aggressive flattening and such because it didn't really
+-- pay off now that we use combined paths with just line segments. I also moved some other
+-- code to a experimental copy. So we now only have the bare helpers needed here.
+
+-- Todo: look into zero case (lavel 1) for shapes ... omiting the background calculation
+-- speeds up quite a bit.
+
+local next, type, tostring = next, type, tostring
+local round, abs, min, max, floor = math.round, math.abs, math.min, math.max, math.floor
+local concat, move = table.concat, table.move
+
+local bor = bit32.bor -- it's really time to ditch support for luajit
+
+local starttiming = statistics.starttiming
+local stoptiming = statistics.stoptiming
+local elapsedtime = statistics.elapsedtime
+
+local formatters = string.formatters
+local setmetatableindex = table.setmetatableindex
+local sortedkeys = table.sortedkeys
+local sequenced = table.sequenced
+
+local metapost = metapost or { }
+local mp = mp or { }
+
+local getparameterset = metapost.getparameterset
+
+local mpflush = mp.flush
+local mpcolor = mp.color
+local mpstring = mp.string
+local mpdraw = mp.draw
+local mpfill = mp.fill
+local mpflatten = mp.flatten
+
+local report = logs.reporter("mfun contour")
+
+local version = 0.007
+
+-- we iterate using integers so that we get a better behaviour at zero
+
+local f_function_n = formatters [ [[
+ local math = math
+ local round = math.round
+ %s
+ return function(data,nx,ny,nxmin,nxmax,xstep,nymin,nymax,ystep)
+ local sx = nxmin
+ for mx=1,nx do
+ local dx = data[mx]
+ local x = sx * xstep
+ local sy = nymin
+ for my=1,ny do
+ local y = sy * ystep
+ dx[my] = %s
+ sy = sy + 1
+ end
+ sx = sx + 1
+ end
+ return 0
+ end
+]] ]
+
+local f_function_y = formatters [ [[
+ local math = math
+ local round = math.round
+ local nan = NaN
+ local inf = math.huge
+ %s
+ return function(data,nx,ny,nxmin,nxmax,xstep,nymin,nymax,ystep,dnan,dinf,report)
+ local sx = nxmin
+ local er = 0
+ for mx=nxmin,nxmax do
+ local dx = data[mx]
+ local x = sx * xstep
+ local sy = nymin
+ for my=nymin,nymax do
+ local y = sy * ystep
+ local n = %s
+ if n == nan then
+ er = er + 1
+ if er < 10 then
+ report("nan at (%s,%s)",x,y)
+ end
+ n = dnan
+ elseif n == inf then
+ er = er + 1
+ if er < 10 then
+ report("inf at (%s,%s)",x,y)
+ end
+ n = dinf
+ end
+ dx[my] = n
+ sy = sy + 1
+ end
+ sx = sx + 1
+ end
+ return er
+ end
+]] ]
+
+local f_color = formatters [ [[
+ local math = math
+ local min = math.min
+ local max = math.max
+ local n = %s
+ local minz = %s
+ local maxz = %s
+
+ local color_value = 0
+ local color_step = mp.lmt_color_functions.step
+ local color_shade = mp.lmt_color_functions.shade
+
+ local function step(...)
+ return color_step(color_value,n,...)
+ end
+ local function shade(...)
+ return color_shade(color_value,n,...)
+ end
+ local function lin(l)
+ return l/n
+ end
+ %s
+ return function(l)
+ color_value = l
+ return %s
+ end
+]] ]
+
+local inbetween = attributes and attributes.colors.helpers.inbetween
+
+mp.lmt_color_functions = {
+ step = function(l,n,r,g,b,s)
+ if not s then
+ s = 1
+ end
+ local f = l / n
+ local r = s * 1.5 - 4 * abs(f-r)
+ local g = s * 1.5 - 4 * abs(f-g)
+ local b = s * 1.5 - 4 * abs(f-b)
+ return min(max(r,0),1), min(max(g,0),1), min(max(b,0),1)
+ end,
+ shade = function(l,n,one,two)
+ local f = l / n
+ local r = inbetween(one,two,1,f)
+ local g = inbetween(one,two,2,f)
+ local b = inbetween(one,two,3,f)
+ return min(max(r,0),1), min(max(g,0),1), min(max(b,0),1)
+ end,
+}
+
+local f_box = formatters [ [[draw rawtexbox("contour",%s) xysized (%s,%s) ;]] ]
+
+local n_box = 0
+
+-- todo: remove old one, so we need to store old hashes
+
+local nofcontours = 0
+
+-- We don't want cosmetics like axis and labels to trigger a calculation,
+-- especially a slow one.
+
+local hashfields = {
+ "xmin", "xmax", "xstep", "ymin", "ymax", "ystep",
+ "levels", "colors", "level", "preamble", "function", "functions", "color", "values",
+ "background", "foreground", "linewidth", "backgroundcolor", "linecolor",
+}
+
+local function prepare(p)
+ local h = { }
+ for i=1,#hashfields do
+ local f = hashfields[i]
+ h[f] = p[f]
+ end
+ local hash = md5.HEX(sequenced(h))
+ local name = formatters["%s-m-c-%03i.lua"](tex.jobname,nofcontours)
+ return name, hash
+end
+
+local function getcache(p)
+ local cache = p.cache
+ if cache then
+ local name, hash = prepare(p)
+ local data = table.load(name)
+ if data and data.hash == hash and data.version == version then
+ p.result = data
+ return true
+ else
+ return false, hash, name
+ end
+ end
+ return false, nil, nil
+end
+
+local function setcache(p)
+ local result = p.result
+ local name = result.name
+ if name then
+ if result.bitmap then
+ result.bitmap = nil
+ else
+ result.data = nil
+ end
+ result.color = nil
+ result.action = nil
+ result.cached = nil
+ io.savedata(name, table.fastserialize(result))
+ else
+ os.remove((prepare(p)))
+ end
+end
+
+function mp.lmt_contours_start()
+
+ starttiming("lmt_contours")
+
+ nofcontours = nofcontours + 1
+
+ local p = getparameterset()
+
+ local xmin = p.xmin
+ local xmax = p.xmax
+ local ymin = p.ymin
+ local ymax = p.ymax
+ local xstep = p.xstep
+ local ystep = p.ystep
+ local levels = p.levels
+ local colors = p.colors
+ local nx = 0
+ local ny = 0
+ local means = setmetatableindex("number")
+ local values = setmetatableindex("number")
+ local data = setmetatableindex("table")
+ local minmean = false
+ local maxmean = false
+
+ local cached, hash, name = getcache(p)
+
+ local function setcolors(preamble,levels,minz,maxz,color,nofvalues)
+ if colors then
+ local function f(k)
+ return #colors[1] == 1 and 0 or { 0, 0, 0 }
+ end
+ setmetatableindex(colors, function(t,k)
+ local v = f(k)
+ t[k] = v
+ return v
+ end)
+ local c = { }
+ local n = 1
+ for i=1,nofvalues do
+ c[i] = colors[n]
+ n = n + 1
+ end
+ return c, f
+ else
+ local tcolor = f_color(levels,minz,maxz,preamble,color)
+ local colors = { }
+ local fcolor = load(tcolor)
+ if type(fcolor) ~= "function" then
+ report("error in color function, case %i: %s",1,color)
+ fcolor = false
+ else
+ fcolor = fcolor()
+ if type(fcolor) ~= "function" then
+ report("error in color function, case %i: %s",2,color)
+ fcolor = false
+ end
+ end
+ if not fcolor then
+ local color_step = mp.lmt_color_functions.step
+ fcolor = function(l)
+ return color_step(l,levels,0.25,0.50,0.75)
+ end
+ end
+ for i=1,nofvalues do
+ colors[i] = { fcolor(i) }
+ end
+ if attributes.colors.model == "cmyk" then
+ for i=1,#colors do
+ local c = colors[i]
+ colors[i] = { 1 - c[1], 1 - c[2], 1 - c[3], 0 }
+ end
+ end
+ return colors, fcolor
+ end
+ end
+
+ if cached then
+ local result = p.result
+ local colors, color = setcolors(
+ p.preamble,
+ p.levels,
+ result.minz,
+ result.maxz,
+ p.color,
+ result.nofvalues
+ )
+ result.color = color
+ result.colors = colors
+ result.cached = true
+ return
+ end
+
+ local functioncode = p["function"]
+ local functionrange = p.range
+ local functionlist = functionrange and p.functions
+ local preamble = p.preamble
+
+ if xstep == 0 then xstep = (xmax - xmin)/100 end
+ if ystep == 0 then ystep = (ymax - ymin)/100 end
+
+ local nxmin = round(xmin/xstep)
+ local nxmax = round(xmax/xstep)
+ local nymin = round(ymin/ystep)
+ local nymax = round(ymax/ystep)
+ local nx = nxmax - nxmin + 1
+ local ny = nymax - nymin + 1
+
+ local function executed(data,code)
+ local fcode = p.check and f_function_y or f_function_n
+ fcode = fcode(preamble,code)
+ fcode = load(fcode)
+ if type(fcode) == "function" then
+ fcode = fcode()
+ end
+ if type(fcode) == "function" then
+ local er = fcode(
+ data, nx, ny,
+ nxmin, nxmax, xstep,
+ nymin, nymax, ystep,
+ p.defaultnan, p.defaultinf, report
+ )
+ if er > 0 then
+ report("%i errors in: %s",er,code)
+ end
+ return true
+ else
+ return false -- fatal error
+ end
+ end
+
+ -- for i=1,nx do
+ -- data[i] = lua.newtable(ny,0)
+ -- end
+
+ if functionlist then
+
+ local datalist = { }
+ local minr = functionrange[1]
+ local maxr = functionrange[2] or minr
+ local size = #functionlist
+
+ for l=1,size do
+
+ local func = setmetatableindex("table")
+ if not executed(func,functionlist[l]) then
+ report("error in executing function %i from functionlist",l)
+ return
+ end
+
+ local bit = l -- + 1
+
+ if l == 1 then
+ for i=1,nx do
+ local di = data[i]
+ local fi = func[i]
+ for j=1,ny do
+ local f = fi[j]
+ if f >= minr and f <= maxr then
+ di[j] = bit
+ else
+ di[j] = 0
+ end
+ end
+ end
+ else
+ for i=1,nx do
+ local di = data[i]
+ local fi = func[i]
+ for j=1,ny do
+ local f = fi[j]
+ if f >= minr and f <= maxr then
+ di[j] = bor(di[j],bit)
+ end
+ end
+ end
+ end
+
+ end
+
+ -- we can simplify the value mess below
+
+ elseif functioncode then
+
+ if not executed(data,functioncode) then
+ report("error in executing function")
+ return -- fatal error
+ end
+
+ else
+
+ report("no valid function(s)")
+ return -- fatal error
+
+ end
+
+ minz = data[1][1]
+ maxz = minz
+
+ for i=1,nx do
+ local d = data[i]
+ for j=1,ny do
+ local v = d[j]
+ if v < minz then
+ minz = v
+ elseif v > maxz then
+ maxz = v
+ end
+ end
+ end
+
+ if functionlist then
+
+ for i=minz,maxz do
+ values[i] = i
+ end
+
+ levels = maxz
+
+ minmean = minz
+ maxmean = maxz
+
+ else
+
+ local snap = (maxz - minz + 1) / levels
+
+ for i=1,nx do
+ local d = data[i]
+ for j=1,ny do
+ local dj = d[j]
+ local v = round((dj - minz) / snap)
+ values[v] = values[v] + 1
+ means [v] = means [v] + dj
+ d[j] = v
+ end
+ end
+
+ local list = sortedkeys(values)
+ local count = values
+ local remap = { }
+
+ values = { }
+
+ for i=1,#list do
+ local v = list[i]
+ local m = means[v] / count[v]
+ remap [v] = i
+ values[i] = m
+ if not minmean then
+ minmean = m
+ maxmean = m
+ elseif m < minmean then
+ minmean = m
+ elseif m > maxmean then
+ maxmean = m
+ end
+ end
+
+ for i=1,nx do
+ local d = data[i]
+ for j=1,ny do
+ d[j] = remap[d[j]]
+ end
+ end
+
+ end
+
+ -- due to rounding we have values + 1 so we can have a floor, ceil, round
+ -- option as well as levels -1
+
+ local nofvalues = #values
+
+ local colors = setcolors(
+ p.preamble,levels,minz,maxz,p.color,nofvalues
+ )
+
+ p.result = {
+ version = version,
+ values = values,
+ nofvalues = nofvalues,
+ minz = minz,
+ maxz = maxz,
+ minmean = minmean,
+ maxmean = maxmean,
+ data = data,
+ color = color,
+ nx = nx,
+ ny = ny,
+ colors = colors,
+ name = name,
+ hash = hash,
+ islist = functionlist and true or false,
+ }
+
+ report("index %i, nx %i, ny %i, levels %i", nofcontours, nx, ny, nofvalues)
+end
+
+function mp.lmt_contours_stop()
+ local p = getparameterset()
+ local e = stoptiming("lmt_contours")
+ setcache(p)
+ p.result = nil
+ local f = p["function"]
+ local l = p.functions
+ report("index %i, %0.3f seconds for: %s",
+ nofcontours, e, "[ " .. concat(l or { f } ," ] [ ") .. " ]"
+ )
+end
+
+function mp.lmt_contours_bitmap_set()
+ local p = getparameterset()
+ local result = p.result
+
+ local values = result.values
+ local nofvalues = result.nofvalues
+ local rawdata = result.data
+ local nx = result.nx
+ local ny = result.ny
+ local colors = result.colors
+ local depth = #colors[1] -- == 3 and "rgb" or "gray"
+
+ -- i need to figure out this offset of + 1
+
+ local bitmap = graphics.bitmaps.new(
+ nx, ny,
+ (depth == 3 and "rgb") or (depth == 4 and "cmyk") or "gray",
+ 1, false, true
+ )
+
+ local palette = bitmap.index or { } -- has to start at 0
+ local data = bitmap.data
+ local p = 0
+
+ if depth == 3 or depth == 4 then
+ for i=1,nofvalues do
+ local c = colors[i]
+ local r = round((c[1] or 0) * 255)
+ local g = round((c[2] or 0) * 255)
+ local b = round((c[3] or 0) * 255)
+ local k = depth == 4 and round((c[4] or 0) * 255) or nil
+ palette[p] = {
+ (r > 255 and 255) or (r < 0 and 0) or r,
+ (g > 255 and 255) or (g < 0 and 0) or g,
+ (b > 255 and 255) or (b < 0 and 0) or b,
+ k
+ }
+ p = p + 1
+ end
+ else
+ for i=1,nofvalues do
+ local s = colors[i][1]
+ local s = round((s or 0) * 255)
+ palette[p] = (
+ (s > 255 and 255) or (s < 0 and 0) or s
+ )
+ p = p + 1
+ end
+ end
+
+ -- As (1,1) is the left top corner so we need to flip of we start in
+ -- the left bottom (we cannot loop reverse because we want a properly
+ -- indexed table.
+
+ local k = 0
+ for y=ny,1,-1 do
+ k = k + 1
+ local d = data[k]
+ for x=1,nx do
+ d[x] = rawdata[x][y] - 1
+ end
+ end
+
+ result.bitmap = bitmap
+end
+
+function mp.lmt_contours_bitmap_get()
+ local p = getparameterset()
+ local result = p.result
+ local bitmap = result.bitmap
+ local box = nodes.hpack(graphics.bitmaps.flush(bitmap))
+ n_box = n_box + 1
+ nodes.boxes.savenode("contour",tostring(n_box),box)
+ return f_box(n_box,bitmap.xsize,bitmap.ysize)
+end
+
+function mp.lmt_contours_cleanup()
+ nodes.boxes.reset("contour")
+ n_box = 0
+end
+
+function mp.lmt_contours_edge_set()
+ local p = getparameterset()
+ local result = p.result
+
+ if result.cached then return end
+
+ local values = result.values
+ local nofvalues = result.nofvalues
+ local data = result.data
+ local nx = result.nx
+ local ny = result.ny
+
+ local xmin = p.xmin
+ local xmax = p.xmax
+ local ymin = p.ymin
+ local ymax = p.ymax
+ local xstep = p.xstep
+ local ystep = p.ystep
+
+ local wsp = { }
+ local edges = { }
+
+ for value=1,nofvalues do
+
+ local iwsp = 0
+ local di = data[1]
+ local dc
+ local edge = { }
+ local e = 0
+ -- the next loop is fast
+ for i=1,nx do
+ local di1 = data[i+1]
+ local dij = di[1]
+ local d = dij - value
+ local dij1
+ for j=1,ny do
+ if j < ny then
+ dij1 = di[j+1]
+ dc = dij1 - value
+ if (d >= 0 and dc < 0) or (d < 0 and dc >= 0) then
+ iwsp = iwsp + 1
+ local y = (d * (j+1) - dc * j) / (d - dc)
+ if i == 1 then
+ wsp[iwsp] = { i, y, 0, (i + (j-1)*nx) }
+ elseif i == nx then
+ wsp[iwsp] = { i, y, (i - 1 + (j-1)*nx), 0 }
+ else
+ local jx = (i + (j-1)*nx)
+ wsp[iwsp] = { i, y, jx - 1, jx }
+ end
+ end
+ end
+ if i < nx then
+ local dc = di1[j] - value
+ if (d >= 0 and dc < 0) or (d < 0 and dc >= 0) then
+ iwsp = iwsp + 1
+ local x = (d * (i+1) - dc * i) / (d - dc)
+ if j == 1 then
+ wsp[iwsp] = { x, j, 0, (i + (j-1)*nx) }
+ elseif j == ny then
+ wsp[iwsp] = { x, j, (i + (j-2)*nx), 0 }
+ else
+ local jx = i + (j-1)*nx
+ wsp[iwsp] = { x, j, jx - nx, jx }
+ end
+ end
+ end
+ dij = dij1
+ d = dc
+ end
+ di = di1
+ end
+ -- the next loop takes time
+ for i=1,iwsp do
+ local wspi = wsp[i]
+ for isq=3,4 do
+ local nsq = wspi[isq]
+ if nsq ~= 0 then
+ local px = wspi[1]
+ local py = wspi[2]
+ local p = { px, py }
+ local pn = 2
+ wspi[isq] = 0
+ while true do
+ for ii=1,iwsp do
+ local w = wsp[ii]
+ local n1 = w[3]
+ local n2 = w[4]
+ if n1 == nsq then
+ local x = w[1]
+ local y = w[2]
+ if x ~= px or y ~= py then
+ pn = pn + 1
+ p[pn] = x
+ pn = pn + 1
+ p[pn] = y
+ px = x
+ py = y
+ end
+ nsq = n2
+ w[3] = 0
+ w[4] = 0
+ if nsq == 0 then
+ if pn == 1 then
+ pn = pn + 1
+ p[pn] = w
+ end
+ goto flush
+ end
+ elseif n2 == nsq then
+ local x = w[1]
+ local y = w[2]
+ if x ~= px or y ~= py then
+ pn = pn + 1
+ p[pn] = x
+ pn = pn + 1
+ p[pn] = y
+ px = x
+ py = y
+ end
+ nsq = n1
+ w[3] = 0
+ w[4] = 0
+ if nsq == 0 then
+ goto flush
+ end
+ end
+ end
+ end
+ ::flush::
+ e = e + 1
+ edge[e] = p
+ if mpflatten then
+ mpflatten(p)
+ end
+ end
+ end
+ end
+
+
+ edges[value] = edge
+
+ end
+
+ result.edges = edges
+
+end
+
+function mp.lmt_contours_shade_set(edgetoo)
+ local p = getparameterset()
+ local result = p.result
+
+ if result.cached then return end
+
+ local values = result.values
+ local nofvalues = result.nofvalues
+ local data = result.data
+ local nx = result.nx
+ local ny = result.ny
+ local color = result.color
+
+ local edges = setmetatableindex("table")
+ local shades = setmetatableindex("table")
+
+ local sqtype = setmetatableindex("table")
+
+ local xspoly = { 0, 0, 0, 0, 0, 0 }
+ local yspoly = { 0, 0, 0, 0, 0, 0 }
+ local xrpoly = { }
+ local yrpoly = { }
+
+ local xrpoly = { } -- lua.newtable(2000,0)
+ local yrpoly = { } -- lua.newtable(2000,0)
+
+ -- for i=1,2000 do
+ -- xrpoly[i] = 0
+ -- yrpoly[i] = 0
+ -- end
+
+ -- Unlike a c compiler lua will not optimize loops to run in parallel so we expand
+ -- some of the loops and make sure we don't calculate when not needed. Not that nice
+ -- but not that bad either. Maybe I should just write this from scratch.
+
+-- local i = 0
+-- local j = 0
+
+ -- Analyze each rectangle separately. Overwrite lower colors
+
+ -- Unrolling the loops and copying code and using constants is faster and doesn't
+ -- produce much more code in the end, also because we then can leave out the not
+ -- seen branches. One can argue about the foundit2* blobs but by stepwise optimizing
+ -- this is the result.
+
+ shades[1] = { { 0, 0, nx - 1, 0, nx - 1, ny - 1, 0, ny - 1 } }
+ edges [1] = { { } }
+
+ -- this is way too slow ... i must have messed up some loop .. what is this with value 1
+
+ for value=1,nofvalues do
+-- for value=2,nofvalues do
+
+ local edge = { }
+ local nofe = 0
+ local shade = { }
+ local nofs = 0
+
+ for i=1,nx-1 do
+ local s = sqtype[i]
+ for j=1,ny-1 do
+ s[j] = 0
+ end
+ end
+
+ local nrp = 0
+
+ local function addedge(a,b,c,d)
+ nofe = nofe + 1 edge[nofe] = a
+ nofe = nofe + 1 edge[nofe] = b
+ nofe = nofe + 1 edge[nofe] = c
+ nofe = nofe + 1 edge[nofe] = d
+ end
+ while true do
+ -- search for a square of type 0 with >= 1 corner above contour level
+ local i
+ local j
+ local d0 = data[1]
+ local d1 = data[2]
+ for ii=1,nx do
+ local s = sqtype[ii]
+ for jj=1,ny do
+ if s[jj] == 0 then
+ if d0[jj] > value then i = ii j = jj goto foundit end
+ if d1[jj] > value then i = ii j = jj goto foundit end
+ local j1 = jj + 1
+ if d1[j1] > value then i = ii j = jj goto foundit end
+ if d0[j1] > value then i = ii j = jj goto foundit end
+ end
+ end
+ d0 = d1
+ d1 = data[ii+1]
+ end
+ break
+ ::foundit::
+ -- initialize r-polygon (nrp seems to be 1 or 2)
+ nrp = nrp + 1
+
+ local first = true
+ local nrpoly = 0
+ local nspoly = 0
+ local nrpm = -nrp
+ -- this is the main loop
+ while true do
+ -- search for a square of type -nrp
+ if first then
+ first = false
+ if sqtype[i][j] == 0 then -- true anyway
+ goto foundit1
+ end
+ end
+ for ii=1,nx do
+ local s = sqtype[ii]
+ for jj=1,ny do
+ if s[jj] == nrpm then
+ i = ii
+ j = jj
+ goto foundit1
+ end
+ end
+ end
+ break
+ ::foundit1::
+ while true do
+
+ -- search current then neighboring squares for square type 0, with a corner in common with current square above contour level
+
+ -- top/bottom ... a bit cheating here
+
+ local i_l, i_c, i_r -- i left current right
+ local j_b, j_c, j_t -- j bottom current top
+
+ local i_n = i + 1 -- i next (right)
+ local j_n = j + 1 -- j next (top)
+
+ local i_p = i - 1 -- i previous (bottom)
+ local j_p = j - 1 -- j previous (right)
+
+ local d_c = data[i]
+ local d_r = data[i_n]
+
+ local sq
+
+ i_c = i ; j_c = j ; if i_c < nx and j_c < ny then sq = sqtype[i_c] if sq[j_c] == 0 then
+ if d_c[j_c] > value then i_l = i_p ; i_r = i_n ; j_b = j_p ; j_t = j_n ; goto foundit21 end
+ if d_c[j_n] > value then i_l = i_p ; i_r = i_n ; j_b = j_p ; j_t = j_n ; goto foundit22 end
+ if d_r[j_c] > value then i_l = i_p ; i_r = i_n ; j_b = j_p ; j_t = j_n ; goto foundit23 end
+ if d_r[j_n] > value then i_l = i_p ; i_r = i_n ; j_b = j_p ; j_t = j_n ; goto foundit24 end
+ end end
+
+ i_c = i_n ; j_c = j ; if i_c < nx and j_c < ny then sq = sqtype[i_c] if sq[j_c] == 0 then
+ if d_r[j_c] > value then i_l = i ; i_r = i_n + 1 ; j_b = j_p ; j_t = j_n ; d_c = d_r ; d_r = data[i_r] ; goto foundit21 end
+ if d_r[j_n] > value then i_l = i ; i_r = i_n + 1 ; j_b = j_p ; j_t = j_n ; d_c = d_r ; d_r = data[i_r] ; goto foundit22 end
+ end end
+
+ i_c = i ; j_c = j_n ; if i_c < nx and j_c < ny then sq = sqtype[i_c] if sq[j_c] == 0 then
+ if d_c[j_n] > value then i_l = i_p ; i_r = i_n ; j_b = j ; j_t = j_n + 1 ; goto foundit21 end
+ if d_r[j_n] > value then i_l = i_p ; i_r = i_n ; j_b = j ; j_t = j_n + 1 ; goto foundit23 end
+ end end
+
+ i_c = i_p ; j_c = j ; if i_c > 0 and j_c < ny then sq = sqtype[i_c] if sq[j_c] == 0 then
+ if d_c[j_c] > value then i_l = i_p - 1 ; i_r = i ; j_b = j_p ; j_t = j_n ; d_r = d_c ; d_c = data[i_p] ; goto foundit23 end
+ if d_c[j_n] > value then i_l = i_p - 1 ; i_r = i ; j_b = j_p ; j_t = j_n ; d_r = d_c ; d_c = data[i_p] ; goto foundit24 end
+ end end
+
+ i_c = i ; j_c = j_p ; if i < nx and j_c > 0 then sq = sqtype[i_c] if sq[j_c] == 0 then
+ if d_c[j] > value then i_l = i_p ; i_r = i_n ; j_b = j_p - 1 ; j_t = j ; goto foundit22 end
+ if d_r[j] > value then i_l = i_p ; i_r = i_n ; j_b = j_p - 1 ; j_t = j ; goto foundit24 end
+ end end
+
+ -- not found
+
+ sqtype[i][j] = nrp
+
+ break
+
+ -- define s-polygon for found square (i_c,j_c) - may have up to 6 sides
+
+ ::foundit21:: -- 1 2 3 4
+
+ sq[j_c] = nrpm
+
+ xspoly[1] = i_l ; yspoly[1] = j_b
+ xspoly[2] = i_c ; yspoly[2] = j_b
+ if d_r[j_c] > value then -- dd2
+ xspoly[3] = i_c ; yspoly[3] = j_c
+ if d_r[j_t] > value then -- dd3
+ xspoly[4] = i_l ; yspoly[4] = j_c
+ if d_c[j_t] > value then -- dd4
+ nspoly = 4
+ else
+ xspoly[5] = i_l ; yspoly[5] = j_c ; nspoly = 5
+ end
+ elseif d_c[j_t] > value then -- dd4
+ xspoly[4] = i_c ; yspoly[4] = j_c ;
+ xspoly[5] = i_l ; yspoly[5] = j_c ; nspoly = 5
+ else
+ xspoly[4] = i_l ; yspoly[4] = j_c ; nspoly = 4
+ if edgetoo then addedge(i_c, j_c, i_l, j_c) end
+ end
+ elseif d_r[j_t] > value then -- dd3
+ xspoly[3] = i_c ; yspoly[3] = j_b
+ xspoly[4] = i_c ; yspoly[4] = j_c
+ if d_c[j_t] > value then -- dd4
+ xspoly[5] = i_l ; yspoly[5] = j_c ; nspoly = 5
+ else
+ xspoly[5] = i_l ; yspoly[5] = j_c ;
+ xspoly[6] = i_l ; yspoly[6] = j_c ; nspoly = 6
+ end
+ elseif d_c[j_t] > value then -- dd4
+ if edgetoo then addedge(i_c, j_b, i_c, j_c) end
+ xspoly[3] = i_c ; yspoly[3] = j_c ;
+ xspoly[4] = i_l ; yspoly[4] = j_c ; nspoly = 4
+ else
+ if edgetoo then addedge(i_c, j_b, i_l, j_c) end
+ xspoly[3] = i_l ; yspoly[3] = j_c ; nspoly = 3
+ end
+ goto done
+
+ ::foundit22:: -- 4 1 2 3
+
+ sq[j_c] = nrpm
+
+ xspoly[1] = i_l ; yspoly[1] = j_c
+ xspoly[2] = i_l ; yspoly[2] = j_b
+ if d_c[j_c] > value then -- dd2
+ xspoly[3] = i_c ; yspoly[3] = j_b
+ if d_r[j_c] > value then -- dd3
+ xspoly[4] = i_c ; yspoly[4] = j_c
+ if d_r[j_t] > value then -- dd4
+ nspoly = 4
+ else
+ xspoly[5] = i_c ; yspoly[5] = j_c ; nspoly = 5 -- suspicious, the same
+ end
+ elseif d_r[j_t] > value then -- dd4
+ xspoly[4] = i_c ; yspoly[4] = j_b ;
+ xspoly[5] = i_c ; yspoly[5] = j_c ; nspoly = 5
+ else
+ if edgetoo then addedge(i_c, j_b, i_c, j_c) end
+ xspoly[4] = i_c ; yspoly[4] = j_c ; nspoly = 4
+ end
+ elseif d_r[j_c] > value then -- dd3
+ xspoly[3] = i_l ; yspoly[3] = j_b
+ xspoly[4] = i_c ; yspoly[4] = j_b
+ xspoly[5] = i_c ; yspoly[5] = j_c
+ if d_r[j_t] > value then -- dd4
+ nspoly = 5
+ else
+ xspoly[6] = i_c ; yspoly[6] = j_c ; nspoly = 6
+ end
+ elseif d_r[j_t] > value then -- dd4
+ if edgetoo then addedge(i_l, j_b, i_c, j_b) end
+ xspoly[3] = i_c ; yspoly[3] = j_b
+ xspoly[4] = i_c ; yspoly[4] = j_c ; nspoly = 4
+ else
+ if edgetoo then addedge(i_l, j_b, i_c, j_c) end
+ xspoly[3] = i_c ; yspoly[3] = j_c ; nspoly = 3
+ end
+ goto done
+
+ ::foundit23:: -- 2 3 4 1
+
+ sq[j_c] = nrpm
+
+ xspoly[1] = i_c ; yspoly[1] = j_b
+ xspoly[2] = i_c ; yspoly[2] = j_c
+ if d_r[j_t] > value then -- dd2
+ xspoly[3] = i_l ; yspoly[3] = j_c
+ if d_c[j_t] > value then -- dd3
+ xspoly[4] = i_l ; yspoly[4] = j_b
+ if d_c[j_c] > value then -- dd4
+ nspoly = 4
+ else
+ xspoly[5] = i_l ; yspoly[5] = j_b ; nspoly = 5
+ end
+ elseif d_c[j_c] > value then -- dd4
+ xspoly[4] = i_l ; yspoly[4] = j_c
+ xspoly[5] = i_l ; yspoly[5] = j_b ; nspoly = 5
+ else
+ if edgetoo then addedge(i_l, j_c, i_l, j_b) end
+ xspoly[4] = i_l ; yspoly[4] = j_b ; nspoly = 4
+ end
+ elseif d_c[j_t] > value then -- dd3
+ xspoly[3] = i_c ; yspoly[3] = j_c
+ xspoly[4] = i_l ; yspoly[4] = j_c
+ xspoly[5] = i_l ; yspoly[5] = j_b
+ if d_c[j_c] > value then -- dd4
+ nspoly = 5
+ else
+ xspoly[6] = i_l ; yspoly[6] = j_b ; nspoly = 6
+ end
+ elseif d_c[j_c] > value then -- dd4
+ if edgetoo then addedge(i_c, j_c, i_l, j_c) end
+ xspoly[3] = i_l ; yspoly[3] = j_c ;
+ xspoly[4] = i_l ; yspoly[4] = j_b ; nspoly = 4
+ else
+ if edgetoo then addedge(i_c, j_c, i_l, j_b) end
+ xspoly[3] = i_l ; yspoly[3] = j_b ; nspoly = 3
+ end
+ goto done
+
+ ::foundit24:: -- 3 4 1 2
+
+ sq[j_c] = nrpm
+
+ xspoly[1] = i_c ; yspoly[1] = j_c
+ xspoly[2] = i_l ; yspoly[2] = j_c
+ if d_c[j_t] > value then -- dd2
+ if d_c[j_c] > value then -- dd3
+ xspoly[3] = i_l ; yspoly[3] = j_b
+ xspoly[4] = i_c ; yspoly[4] = j_b
+ if d_r[j_c] > value then -- dd4
+ nspoly = 4
+ else
+ xspoly[5] = i_c ; yspoly[5] = j_b ; nspoly = 5
+ end
+ else
+ xspoly[3] = i_l ; yspoly[3] = j_b
+ if d_r[j_c] > value then -- dd4
+
+ local xv34 = (dd3*i_c-dd4*i_l)/(dd3 - dd4) -- probably i_l
+ print("4.4 : xv34",xv34,i_c,i_l)
+
+ -- if edgetoo then addedge(i_l, j_b, xv34, j_b) end
+ xspoly[4] = xv34 ; yspoly[4] = j_b ;
+ xspoly[5] = i_c ; yspoly[5] = j_b ; nspoly = 5
+ else
+ if edgetoo then addedge(i_l, j_b, i_c, j_b) end
+ xspoly[4] = i_c ; yspoly[4] = j_b ; nspoly = 4
+ end
+ end
+ elseif d_c[j_c] > value then -- dd3
+ xspoly[3] = i_l ; yspoly[3] = j_b
+ xspoly[4] = i_l ; yspoly[4] = j_b
+ xspoly[5] = i_c ; yspoly[5] = j_b
+ if d_r[j_c] > value then -- dd4
+ nspoly = 5
+ else
+ xspoly[6] = i_c ; yspoly[6] = j_b ; nspoly = 6
+ end
+ elseif d_r[j_c] > value then -- dd4
+ if edgetoo then addedge(i_l, j_c, i_l, j_b) end
+ xspoly[3] = i_l ; yspoly[3] = j_b
+ xspoly[4] = i_c ; yspoly[4] = j_b ; nspoly = 4
+ else
+ if edgetoo then addedge(i_l, j_c, i_c, j_b) end
+ xspoly[3] = i_c ; yspoly[3] = j_b ; nspoly = 3
+ end
+ -- goto done
+
+ ::done::
+ -- combine s-polygon with existing r-polygon, eliminating redundant segments
+
+ if nrpoly == 0 then
+ -- initiate r-polygon
+ for i=1,nspoly do
+ xrpoly[i] = xspoly[i]
+ yrpoly[i] = yspoly[i]
+ end
+ nrpoly = nspoly
+ else
+ -- search r-polygon and s-polygon for one side that matches
+ --
+ -- this is a bottleneck ... we keep this variant here but next go for a faster
+ -- alternative
+ --
+ -- local ispoly, irpoly
+ -- for r=nrpoly,1,-1 do
+ -- local r1
+ -- for s=1,nspoly do
+ -- local s1 = s % nspoly + 1
+ -- if xrpoly[r] == xspoly[s1] and yrpoly[r] == yspoly[s1] then
+ -- if not r1 then
+ -- r1 = r % nrpoly + 1
+ -- end
+ -- if xrpoly[r1] == xspoly[s] and yrpoly[r1] == yspoly[s] then
+ -- ispoly = s
+ -- irpoly = r
+ -- goto foundit3
+ -- end
+ -- end
+ -- end
+ -- end
+ --
+ -- local ispoly, irpoly
+ -- local xr1 = xrpoly[1]
+ -- local yr1 = yrpoly[1]
+ -- for r0=nrpoly,1,-1 do
+ -- for s0=1,nspoly do
+ -- if xr1 == xspoly[s0] and yr1 == yspoly[s0] then
+ -- if s0 == nspoly then
+ -- if xr0 == xspoly[1] and yr0 == yspoly[1] then
+ -- ispoly = s0
+ -- irpoly = r0
+ -- goto foundit3
+ -- end
+ -- else
+ -- local s1 = s0 + 1
+ -- if xr0 == xspoly[s1] and yr0 == yspoly[s1] then
+ -- ispoly = s0
+ -- irpoly = r0
+ -- goto foundit3
+ -- end
+ -- end
+ -- end
+ -- end
+ -- xr1 = xrpoly[r0]
+ -- yr1 = yrpoly[r0]
+ -- end
+ --
+ -- but ...
+ --
+ local minx = xspoly[1]
+ local miny = yspoly[1]
+ local maxx = xspoly[1]
+ local maxy = yspoly[1]
+ for i=1,nspoly do
+ local y = yspoly[i]
+ if y < miny then
+ miny = y
+ elseif y > maxy then
+ maxy = y
+ end
+ local x = xspoly[i]
+ if x < minx then
+ minx = y
+ elseif x > maxx then
+ maxx = x
+ end
+ end
+ -- we can delay accessing y ...
+ local ispoly, irpoly
+ local xr1 = xrpoly[1]
+ local yr1 = yrpoly[1]
+ for r0=nrpoly,1,-1 do
+ if xr1 >= minx and xr1 <= maxx and yr1 >= miny and yr1 <= maxy then
+ local xr0 = xrpoly[r0]
+ local yr0 = yrpoly[r0]
+ for s0=1,nspoly do
+ if xr1 == xspoly[s0] and yr1 == yspoly[s0] then
+ if s0 == nspoly then
+ if xr0 == xspoly[1] and yr0 == yspoly[1] then
+ ispoly = s0
+ irpoly = r0
+ goto foundit3
+ end
+ else
+ local s1 = s0 + 1
+ if xr0 == xspoly[s1] and yr0 == yspoly[s1] then
+ ispoly = s0
+ irpoly = r0
+ goto foundit3
+ end
+ end
+ end
+ end
+ xr1 = xr0
+ yr1 = yr0
+ else
+ xr1 = xrpoly[r0]
+ yr1 = yrpoly[r0]
+ end
+ end
+ --
+ goto nomatch3
+ ::foundit3::
+ local match1 = 0
+ local rpoly1 = irpoly + nrpoly
+ local spoly1 = ispoly - 1
+ for i=2,nspoly-1 do
+ -- search for further matches nearby
+ local ir = (rpoly1 - i) % nrpoly + 1
+ local is = (spoly1 + i) % nspoly + 1
+ if xrpoly[ir] == xspoly[is] and yrpoly[ir] == yspoly[is] then
+ match1 = match1 + 1
+ else
+ break -- goto nomatch1
+ end
+ end
+ ::nomatch1::
+ local match2 = 0
+ local rpoly2 = irpoly - 1
+ local spoly2 = ispoly + nspoly
+ for i=2,nspoly-1 do
+ -- search other way for further matches nearby
+ local ir = (rpoly2 + i) % nrpoly + 1
+ local is = (spoly2 - i) % nspoly + 1
+ if xrpoly[ir] == xspoly[is] and yrpoly[ir] == yspoly[is] then
+ match2 = match2 + 1
+ else
+ break -- goto nomatch2
+ end
+ end
+ ::nomatch2::
+ -- local dnrpoly = nspoly - 2 - 2*match1 - 2*match2
+ local dnrpoly = nspoly - 2*(match1 + match2 + 1)
+ local ispolystart = (ispoly + match1) % nspoly + 1 -- first node of s-polygon to include
+ local irpolyend = (rpoly1 - match1 - 1) % nrpoly + 1 -- last node of s-polygon to include
+ if dnrpoly ~= 0 then
+ local irpolystart = (irpoly + match2) % nrpoly + 1 -- first node of s-polygon to include
+ if irpolystart > irpolyend then
+ -- local ispolyend = (spoly1 - match2 + nspoly)%nspoly + 1 -- last node of s-polygon to include
+ if dnrpoly > 0 then
+ -- expand the arrays xrpoly and yrpoly
+ for i=nrpoly,irpolystart,-1 do
+ local k = i + dnrpoly
+ xrpoly[k] = xrpoly[i]
+ yrpoly[k] = yrpoly[i]
+ end
+ else -- if dnrpoly < 0 then
+ -- contract the arrays xrpoly and yrpoly
+ for i=irpolystart,nrpoly do
+ local k = i + dnrpoly
+ xrpoly[k] = xrpoly[i]
+ yrpoly[k] = yrpoly[i]
+ end
+ end
+ end
+ nrpoly = nrpoly + dnrpoly
+ end
+ if nrpoly < irpolyend then
+ for i=irpolyend,nrpoly+1,-1 do
+ -- otherwise these values get lost!
+ local k = i - nrpoly
+ xrpoly[k] = xrpoly[i]
+ yrpoly[k] = yrpoly[i]
+ end
+ end
+ local n = nspoly - 2 - match1 - match2
+ if n == 1 then
+ local irpoly1 = irpolyend % nrpoly + 1
+ local ispoly1 = ispolystart % nspoly + 1
+ xrpoly[irpoly1] = xspoly[ispoly1]
+ yrpoly[irpoly1] = yspoly[ispoly1]
+ elseif n > 0 then
+ -- often 2
+ for i=1,n do
+ local ii = i - 1
+ local ir = (irpolyend + ii) % nrpoly + 1
+ local is = (ispolystart + ii) % nspoly + 1
+ xrpoly[ir] = xspoly[is]
+ yrpoly[ir] = yspoly[is]
+ end
+ end
+ ::nomatch3::
+ end
+ end
+ end
+
+ if nrpoly > 0 then
+ local t = { }
+ local n = 0
+ for i=1,nrpoly do
+ n = n + 1 t[n] = xrpoly[i]
+ n = n + 1 t[n] = yrpoly[i]
+ end
+ if mpflatten then
+ mpflatten(t) -- maybe integrate
+ end
+ nofs = nofs + 1
+ shade[nofs] = t
+ -- print(value,nrpoly,#t,#t-nrpoly*2)
+ end
+
+ end
+
+ edges [value+1] = edge
+ shades[value+1] = shade
+-- edges [value] = edge
+-- shades[value] = shade
+ end
+
+ result.shades = shades
+ result.shapes = edges
+
+end
+
+-- accessors
+
+function mp.lmt_contours_nx (i) return getparameterset().result.nx end
+function mp.lmt_contours_ny (i) return getparameterset().result.ny end
+
+function mp.lmt_contours_nofvalues() return getparameterset().result.nofvalues end
+function mp.lmt_contours_value (i) return getparameterset().result.values[i] end
+
+function mp.lmt_contours_minz (i) return getparameterset().result.minz end
+function mp.lmt_contours_maxz (i) return getparameterset().result.maxz end
+
+function mp.lmt_contours_minmean (i) return getparameterset().result.minmean end
+function mp.lmt_contours_maxmean (i) return getparameterset().result.maxmean end
+
+function mp.lmt_contours_xrange () local p = getparameterset() mpstring(formatters["x = [%.3N,%.3N] ;"](p.xmin,p.xmax)) end
+function mp.lmt_contours_yrange () local p = getparameterset() mpstring(formatters["y = [%.3N,%.3N] ;"](p.ymin,p.ymax)) end
+
+function mp.lmt_contours_format()
+ local p = getparameterset()
+ return mpstring(p.result.islist and "@i" or p.zformat or p.format)
+end
+
+function mp.lmt_contours_function()
+ local p = getparameterset()
+ return mpstring(p.result.islist and concat(p["functions"], ", ") or p["function"])
+end
+
+function mp.lmt_contours_range()
+ local p = getparameterset()
+ local r = p.result.islist and p.range
+ if not r or #r == 0 then
+ return mpstring("")
+ elseif #r == 1 then
+ return mpstring(r[1])
+ else
+ return mpstring(formatters["z = [%s,%s]"](r[1],r[2]))
+ end
+end
+
+function mp.lmt_contours_edge_paths(value)
+ mpdraw(getparameterset().result.edges[value],true)
+ mpflush()
+end
+
+function mp.lmt_contours_shape_paths(value)
+ mpdraw(getparameterset().result.shapes[value],false)
+ mpflush()
+end
+
+function mp.lmt_contours_shade_paths(value)
+ mpfill(getparameterset().result.shades[value],true)
+ mpflush()
+end
+
+function mp.lmt_contours_color(value)
+ local p = getparameterset()
+ local color = p.result.colors[value]
+ if color then
+ mpcolor(color)
+ end
+end
+
+-- The next code is based on the wikipedia page. It was a bit tedius job to define the
+-- coordinates but hupefully I made no errors. I rendered all shapes independently and
+-- tripple checked bit one never knows ...
+
+-- maybe some day write from scatch, like this (axis are swapped):
+
+local d = 1/2
+
+local paths = {
+ { 0, d, d, 0 },
+ { 1, d, d, 0 },
+ { 0, d, 1, d },
+ { 1, d, d, 1 },
+ { 0, d, d, 1, d, 0, 1, d }, -- saddle
+ { d, 0, d, 1 },
+ { 0, d, d, 1 },
+ { 0, d, d, 1 },
+ { d, 0, d, 1 },
+ { 0, d, d, 0, 1, d, d, 1 }, -- saddle
+ { 1, d, d, 1 },
+ { 0, d, 1, d },
+ { d, 0, 1, d },
+ { d, 0, 0, d },
+}
+
+local function whatever(data,nx,ny,threshold)
+ local edges = { }
+ local e = 0
+ local d0 = data[1]
+ for j=1,ny-1 do
+ local d1 = data[j+1]
+ local k = j + 1
+ for i=1,nx-1 do
+ local v = 0
+ local l = i + 1
+ local c1 = d0[i]
+ if c1 < threshold then
+ v = v + 8
+ end
+ local c2 = d0[l]
+ if c2 < threshold then
+ v = v + 4
+ end
+ local c3 = d1[l]
+ if c3 < threshold then
+ v = v + 2
+ end
+ local c4 = d1[i]
+ if c4 < threshold then
+ v = v + 1
+ end
+ if v > 0 and v < 15 then
+ if v == 5 or v == 10 then
+ local a = (c1 + c2 + c3 + c4) / 4
+ if a < threshold then
+ v = v == 5 and 10 or 5
+ end
+ local p = paths[v]
+ e = e + 1 edges[e] = k - p[2]
+ e = e + 1 edges[e] = i + p[1]
+ e = e + 1 edges[e] = k - p[4]
+ e = e + 1 edges[e] = i + p[3]
+ e = e + 1 edges[e] = k - p[6]
+ e = e + 1 edges[e] = i + p[5]
+ e = e + 1 edges[e] = k - p[8]
+ e = e + 1 edges[e] = i + p[7]
+ else
+ local p = paths[v]
+ e = e + 1 edges[e] = k - p[2]
+ e = e + 1 edges[e] = i + p[1]
+ e = e + 1 edges[e] = k - p[4]
+ e = e + 1 edges[e] = i + p[3]
+ end
+ end
+ end
+ d0 = d1
+ end
+ return edges
+end
+
+-- todo: just fetch when needed, no need to cache
+
+function mp.lmt_contours_edge_set_by_cell()
+ local p = getparameterset()
+ local result = p.result
+
+ if result.cached then return end
+
+ local values = result.values
+ local nofvalues = result.nofvalues
+ local data = result.data
+ local nx = result.nx
+ local ny = result.ny
+ local lines = { }
+ result.lines = lines
+ for value=1,nofvalues do
+ lines[value] = whatever(data,ny,nx,value)
+ end
+end
+
+function mp.lmt_contours_edge_get_cell(value)
+ mpdraw(getparameterset().result.lines[value])
+ mpflush()
+end
+
+local singles = {
+ { d, 0, 0, 0, 0, d }, -- 1 0001
+ { d, 0, 0, d }, -- 2 0002
+ { 1, d, 1, 0, d, 0 }, -- 3 0010
+ { 1, d, 1, 0, 0, 0, 0, d }, -- 4 0011
+ { 1, d, 1, 0, d, 0, 0, d }, -- 5 0012
+ { 1, d, d, 0 }, -- 6 0020
+ { 1, d, d, 0, 0, 0, 0, d }, -- 7 0021
+ { 1, d, 0, d }, -- 8 0022
+ { d, 1, 1, 1, 1, d }, -- 9 0100
+ false, -- 10 0101
+ false, -- 11 0102
+ { d, 1, 1, 1, 1, 0, d, 0 }, -- 12 0110
+ { d, 1, 1, 1, 1, 0, 0, 0, 0, d }, -- 13 0111
+ { d, 1, 1, 1, 1, 0, d, 0, 0, d }, -- 14 0112
+ { d, 1, 1, 1, 1, d, d, 0 }, -- 15 0120
+ { d, 1, 1, 1, 1, d, d, 0, 0, 0, 0, d }, -- 16 0121
+ { d, 1, 1, 1, 1, d, 0, d }, -- 17 0122
+ { d, 1, 1, d }, -- 18 0200
+ false, -- 19 0201
+ false, -- 20 0202
+ { d, 1, 1, d, 1, 0, d, 0 }, -- 21 0210
+ { d, 1, 1, d, 1, 0, 0, 0, 0, d }, -- 22 0211
+ false, -- 23 0212
+ { d, 1, d, 0 }, -- 24 0220
+ { d, 1, d, 0, 0, 0, 0, d }, -- 25 0221
+ { d, 1, 0, d }, -- 26 0222
+ { 0, 1, d, 1, 0, d }, -- 27 1000
+ { 0, 1, d, 1, d, 0, 0, 0 }, -- 28 1001
+ { 0, 1, d, 1, d, 0, 0, d }, -- 29 1002
+ false, -- 30 1010
+ { 0, 1, d, 1, 1, d, 1, 0, 0, 0 }, -- 31 1011
+ { 0, 1, d, 1, 1, d, 1, 0, d, 0, 0, d }, -- 32 1012
+ false, -- 33 1020
+ { 0, 1, d, 1, 1, d, d, 0, 0, 0 }, -- 34 1021
+ { 0, 1, d, 1, 1, d, 0, d }, -- 35 1022
+ { 0, 1, 1, 1, 1, d, 0, d }, -- 36 1100
+ { 0, 1, 1, 1, 1, d, d, 0, 0, 0 }, -- 37 1101
+ { 0, 1, 1, 1, 1, d, d, 0, 0, d }, -- 38 1102
+ { 0, 1, 1, 1, 1, 0, d, 0, 0, d }, -- 39 1110
+ { 0, 1, 1, 1, 1, 0, 0, 0 }, -- 40 1111
+ { 0, 1, 1, 1, 1, 0, d, 0, 0, d }, -- 41 1112
+ { 0, 1, 1, 1, 1, d, d, 0, 0, d }, -- 42 1120
+ { 0, 1, 1, 1, 1, d, d, 0, 0, 0 }, -- 43 1121
+ { 0, 1, 1, 1, 1, d, 0, d }, -- 44 1122
+ { 0, 1, d, 1, 1, d, 0, d }, -- 45 1200
+ { 0, 1, d, 1, 1, d, d, 0, 0, 0 }, -- 46 1201
+ false, -- 47 1202
+ { 0, 1, d, 1, 1, d, 1, 0, d, 0, 0, d }, -- 48 1210
+ { 0, 1, d, 1, 1, d, 1, 0, 0, 0 }, -- 49 1211
+ false, -- 50 1212
+ { 0, 1, d, 1, d, 0, 0, d }, -- 51 1220
+ { 0, 1, d, 1, d, 0, 0, 0 }, -- 52 1221
+ { 0, 1, d, 1, 0, d }, -- 53 1222
+ { d, 1, 0, d }, -- 54 2000
+ { d, 1, d, 0, 0, 0, 0, d }, -- 55 2001
+ { d, 1, d, 0 }, -- 56 2002
+ false, -- 57 2010
+ { d, 1, 1, d, 1, 0, 0, 0, 0, d }, -- 58 2011
+ { d, 1, 1, d, 1, 0, d, 0 }, -- 59 2012
+ false, -- 60 2020
+ false, -- 61 2021
+ { d, 1, 1, d }, -- 62 2022
+ { d, 1, 1, 1, 1, d, 0, d }, -- 63 2100
+ { d, 1, 1, 1, 1, d, d, 0, 0, 0, 0, d }, -- 64 2101
+ { d, 1, 1, 1, 1, d, d, 0 }, -- 65 2102
+ { d, 1, 1, 1, 1, 0, d, 0, 0, d }, -- 66 2110
+ { d, 1, 1, 1, 1, 0, 0, 0, 0, d }, -- 67 2111
+ { d, 1, 1, 1, 1, 0, d, 0 }, -- 68 2112
+ false, -- 69 2120
+ false, -- 70 2121
+ { d, 1, 1, 1, 1, d }, -- 71 2122
+ { 1, d, 0, d }, -- 72 2200
+ { 1, d, d, 0, 0, 0, 0, d }, -- 73 2201
+ { 1, d, d, 0 }, -- 74 2202
+ { 1, d, 1, 0, d, 0, 0, d }, -- 75 2210
+ { 1, d, 1, 0, 0, 0, 0, d }, -- 76 2211
+ { 1, d, 1, 0, d, 0 }, -- 77 2212
+ { d, 0, 0, d }, -- 78 2220
+ { 0, d, 0, 0, d, 0 }, -- 79 2221
+}
+
+local sadles = {
+ false, false, false, false, false, false, false, false, false,
+ { { d, 1, 1, 1, 1, d }, { d, 0, 0, 0, 0, d }, { d, 1, 1, 1, 1, d, d, 0, 0, 0, 0, d }, false, false, false }, -- 10 0101
+ { { d, 1, 1, 1, 1, d }, { d, 0, 0, d }, { d, 1, 1, 1, 1, d, d, 0, 0, d }, false, false, false }, -- 11 0102
+ false, false, false, false, false, false, false,
+ { { d, 1, 1, d }, { d, 0, 0, 0, 0, d }, { d, 1, 1, d, d, 0, 0, 0, 0, d }, false, false, false }, -- 19 0201
+ { { d, 1, 1, d }, { d, 0, 0, d }, { d, 1, 1, d, d, 0, 0, d }, false, { d, 1, 0, d }, { 1, d, d, 0 } }, -- 20 0202
+ false, false,
+ { false, false, { d, 1, 1, d, 1, 0, d, 0, 0, d }, false, { d, 1, 0,d, }, { 1, d, 1, 0,d, 0 } }, -- 23 0212
+ false, false, false, false, false, false,
+ { { 0, 1, d, 1, 0, d }, { 1, d, 1, 0, d, 0 }, { 0, 1, d, 1, 1, d, 1, 0, d, 0, 0, d }, false, false, false }, -- 30 1010
+ false, false,
+ { { 1, 0, d, 0, 0, d, }, { 1, d, d, 0 }, { 0, 1, d, 1, 1, d, d, 0, 0, d }, false, false, false }, -- 33 1020
+ false, false, false, false, false, false, false, false, false, false, false, false, false,
+ { false, false, { 0,1, d, 1, 1, d, d, 0, 0, d }, false, { 0,1, d, 1, 0, d }, {1, d, d, 0 } }, -- 47 1202
+ false, false,
+ { false, false, { 0, 1, d, 1, 1, d, 1, 0, d, 0, 0, d }, false, { 0, 1, d, 1, 0, d }, { 1, d, 1, 0, d, 0 } }, -- 50 1212
+ false, false, false, false, false, false,
+ { { d, 1, 0, d }, { 1, d, 1, 0, 0, d }, { d, 1, 1, d, 1, 0, d, 0, 0, d }, false, false, false }, -- 57 2010
+ false, false,
+ { { d, 1, 0,d }, { 1, d, d, 0 }, { d, 1, 1, d, d, 0, 0, d }, false, { d, 1, 1, d }, { d, 0, 0, d } }, -- 60 2020
+ { false, false, { d, 1, 1, d, d, 0, 0, 0, 0, d }, false, { d, 1, 1, d }, { d, 0, 0, 0, 0, d } }, -- 61 2021
+ false, false, false, false, false, false, false,
+ { false, false, { d, 1, 1, 1, 1, d, d, 0, 0, d }, false, { d, 1, 1, 1, 1, d }, { d, 0,0,d } }, -- 69 2120
+ { false, false, { d, 1, 1, 1, 1, d, d, 0, 0, 0, 0, d }, false, { d, 1, 1, 1, 1, d }, { d, 0, 0, 0, 0, d } }, -- 70 2121
+}
+
+local function whatever(data,nx,ny,threshold,background)
+
+ if background then
+
+ local llx = 1/2
+ local lly = llx
+ local urx = ny + llx
+ local ury = nx + lly
+
+ return { { llx, lly, urx, 0, urx, ury, 0, ury } }
+
+ else
+
+ local bands = { }
+ local b = 0
+
+ local function band(s,n,x,y) -- simple. no closure so fast
+ if n == 6 then
+ return {
+ x - s[ 2], y + s[ 1], x - s[ 4], y + s[ 3], x - s[ 6], y + s[ 5],
+ }
+ elseif n == 8 then
+ return {
+ x - s[ 2], y + s[ 1], x - s[ 4], y + s[ 3], x - s[ 6], y + s[ 5],
+ x - s[ 8], y + s[ 7],
+ }
+ elseif n == 10 then
+ return {
+ x - s[ 2], y + s[ 1], x - s[ 4], y + s[ 3], x - s[ 6], y + s[ 5],
+ x - s[ 8], y + s[ 7], x - s[10], y + s[ 9],
+ }
+ elseif n == 4 then
+ return {
+ x - s[ 2], y + s[ 1], x - s[ 4], y + s[ 3],
+ }
+ else -- 12
+ return {
+ x - s[ 2], y + s[ 1], x - s[ 4], y + s[ 3], x - s[ 6], y + s[ 5],
+ x - s[ 8], y + s[ 7], x - s[10], y + s[ 9], x - s[12], y + s[11],
+ }
+ end
+ end
+
+ local pp = { }
+
+ local d0 = data[1]
+ for j=1,ny-1 do
+ local d1 = data[j+1]
+ local k = j + 1
+ local p = false
+ for i=1,nx-1 do
+ local v = 0
+ local l = i + 1
+ local c1 = d0[i]
+ if c1 == threshold then
+ v = v + 27
+ elseif c1 > threshold then
+ v = v + 54
+ end
+ local c2 = d0[l]
+ if c2 == threshold then
+ v = v + 9
+ elseif c2 > threshold then
+ v = v + 18
+ end
+ local c3 = d1[l]
+ if c3 == threshold then
+ v = v + 3
+ elseif c3 > threshold then
+ v = v + 6
+ end
+ local c4 = d1[i]
+ if c4 == threshold then
+ v = v + 1
+ elseif c4 > threshold then
+ v = v + 2
+ end
+ if v > 0 and v < 80 then
+ if v == 40 then
+ -- a little optimization: full areas appended horizontally
+ if p then
+ p[4] = l -- i + 1
+ p[6] = l -- i + 1
+ else
+ -- x-0 y+1 x-1 y+1 x-1 y+0 x-0 y+0
+ p = { j, i, j, l, k, l, k, i }
+ b = b + 1 ; bands[b] = p
+ end
+ else
+ local s = singles[v]
+ if s then
+ b = b + 1 ; bands[b] = band(s,#s,k,i)
+ else
+ local s = sadles[v]
+ if s then
+ local m = (c1 + c2 + c3 + c4) / 4
+ if m < threshold then
+ local s1 = s[1] if s1 then b = b + 1 ; bands[b] = band(s1,#s1,i,j) end
+ local s2 = s[2] if s2 then b = b + 1 ; bands[b] = band(s2,#s2,i,j) end
+ elseif m == threshold then
+ local s3 = s[3] if s3 then b = b + 1 ; bands[b] = band(s3,#s3,i,j) end
+ local s4 = s[4] if s4 then b = b + 1 ; bands[b] = band(s4,#s4,i,j) end
+ else
+ local s5 = s[5] if s5 then b = b + 1 ; bands[b] = band(s5,#s5,i,j) end
+ local s6 = s[6] if s6 then b = b + 1 ; bands[b] = band(s6,#s6,i,j) end
+ end
+ end
+ end
+ p = false
+ end
+ else
+ p = false
+ end
+ end
+ d0 = d1
+ end
+ return bands
+ end
+end
+
+function mp.lmt_contours_edge_set_by_band(value)
+ local p = getparameterset()
+ local result = p.result
+
+ if result.cached then return end
+
+ local values = result.values
+ local nofvalues = result.nofvalues
+ local data = result.data
+ local nx = result.nx
+ local ny = result.ny
+ local bands = { }
+ result.bands = bands
+ for value=1,nofvalues do
+ bands[value] = whatever(data,ny,nx,value,value == 1)
+ end
+end
+
+function mp.lmt_contours_edge_get_band(value)
+ mpfill(getparameterset().result.bands[value],true)
+ mpflush()
+end
+
+-- Because we share some code surface plots also end up here. When working on the
+-- contour macros by concidence I ran into a 3D plot in
+--
+-- https://staff.science.uva.nl/a.j.p.heck/Courses/mptut.pdf
+--
+-- The code is pure MetaPost and works quite well. With a bit of optimization
+-- performance is also ok, but in the end a Lua solution is twice as fast and also
+-- permits some more tweaking at no cost. So, below is an adaptation of an example
+-- in the mentioned link. It's one of these cases where access to pseudo arrays
+-- is slowing down MP.
+
+local sqrt, sin, cos = math.sqrt, math.sin, math.cos
+
+local f_fill_rgb = formatters["F (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor (%.3N,%.3N,%.3N) ;"]
+local f_draw_rgb = formatters["D (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor %.3F ;"]
+local f_mesh_rgb = formatters["U (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor (%.3N,%.3N,%.3N) ;"]
+local f_fill_cmy = formatters["F (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor (%.3N,%.3N,%.3N,0) ;"]
+local f_draw_cmy = formatters["D (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor %.3F ;"]
+local f_mesh_cmy = formatters["U (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor (%.3N,%.3N,%.3N,0) ;"]
+
+local f_function_n = formatters [ [[
+ local math = math
+ local round = math.round
+ %s
+ return function(x,y)
+ return %s
+ end
+]] ]
+
+local f_function_y = formatters [ [[
+ local math = math
+ local round = math.round
+ local nan = NaN
+ local inf = math.huge
+ local er = 0
+ %s
+ return function(x,y,dnan,dinf,report)
+ local n = %s
+ if n == nan then
+ er = er + 1
+ if er < 10 then
+ report("nan at (%s,%s)",x,y)
+ end
+ n = dnan
+ elseif n == inf then
+ er = er + 1
+ if er < 10 then
+ report("inf at (%s,%s)",x,y)
+ end
+ n = dinf
+ end
+ dx[my] = n
+ sy = sy + 1
+ end
+ return n, er
+end
+]] ]
+
+local f_color = formatters [ [[
+ local math = math
+ return function(f)
+ return %s
+ end
+]] ]
+
+function mp.lmt_surface_do(specification)
+ --
+ -- The projection and color brightness calculation have been inlined. We also store
+ -- differently.
+ --
+ -- todo: ignore weird paths
+ --
+ -- The prototype is now converted to use lmt parameter sets.
+ --
+ local p = getparameterset("surface")
+ --
+ local preamble = p.preamble or ""
+ local code = p.code or "return x + y"
+ local colorcode = p.color or "return f, f, f"
+ local linecolor = p.linecolor or 1
+ local xmin = p.xmin or -1
+ local xmax = p.xmax or 1
+ local ymin = p.ymin or -1
+ local ymax = p.ymax or 1
+ local xstep = p.xstep or .1
+ local ystep = p.ystep or .1
+ local bf = p.brightness or 100
+ local clip = p.clip or false
+ local lines = p.lines
+ local ha = p.snap or 0.01
+ local hb = 2 * ha
+ --
+ if lines == nil then lines = true end
+ --
+ if xstep == 0 then xstep = (xmax - xmin)/100 end
+ if ystep == 0 then ystep = (ymax - ymin)/100 end
+
+ local nxmin = round(xmin/xstep)
+ local nxmax = round(xmax/xstep)
+ local nymin = round(ymin/ystep)
+ local nymax = round(ymax/ystep)
+ local nx = nxmax - nxmin + 1
+ local ny = nymax - nymin + 1
+ --
+ local xvector = p.xvector or { -0.7, -0.7 }
+ local yvector = p.yvector or { 1, 0 }
+ local zvector = p.zvector or { 0, 1 }
+ local light = p.light or { 3, 3, 10 }
+ --
+ local xrx, xry = xvector[1], xvector[2]
+ local yrx, yry = yvector[1], yvector[2]
+ local zrx, zry = zvector[1], zvector[2]
+ local xp, yp, zp = light[1], light[2], light[3]
+ --
+ local data = setmetatableindex("table")
+ local dx = (xmax - xmin) / nx
+ local dy = (ymax - ymin) / ny
+ local xt = xmin
+ --
+ local minf, maxf
+ --
+ -- similar as contours but no data loop here
+ --
+ local fcode = load((p.check and f_function_y or f_function_n)(preamble,code))
+ local func = type(fcode) == "function" and fcode()
+ if type(func) ~= "function" then
+ return false -- fatal error
+ end
+ --
+ local ccode = load(f_color(colorcode))
+ local color = type(ccode) == "function" and ccode()
+ if type(color) ~= "function" then
+ return false -- fatal error
+ end
+ --
+ for i=0,nx do
+ local yt = ymin
+ for j=0,ny do
+ local zt = func(xt,yt)
+ -- projection from 3D to 2D coordinates
+ local x = xt * xrx + yt * yrx + zt * zrx
+ local y = xt * xry + yt * yry + zt * zry
+ local z = zt
+ -- numerical derivatives by central differences
+ local dfx = (func(xt+ha,yt) - func(xt-ha,yt)) / hb
+ local dfy = (func(xt,yt+ha) - func(xt,yt-ha)) / hb
+ -- compute brightness factor at a point
+ local ztp = zt - zp
+ local ytp = yt - yp
+ local xtp = xt - xp
+ local ztp = zt - zp
+ local ytp = yt - yp
+ local xtp = xt - xp
+ local ca = -ztp + dfy*ytp + dfx*xtp
+ local cb = sqrt(1+dfx*dfx+dfy*dfy)
+ local cc = sqrt(ztp*ztp + ytp*ytp + xtp*xtp)
+ local fac = bf*ca/(cb*cc*cc*cc)
+ -- addition: check range
+ if not minf then
+ minf = fac
+ maxf = fac
+ elseif fac < minf then
+ minf = fac
+ elseif fac > maxf then
+ maxf = fac
+ end
+ --
+ data[i][j] = { x, y, fac }
+ --
+ yt = yt + dy
+ end
+ xt = xt + dx
+ end
+ local result = { }
+ local r = 0
+ local range = maxf - minf
+ local cl = linecolor or 1
+ local enforce = attributes.colors.model == "cmyk"
+ for i=0,nx-1 do
+ for j=0,ny-1 do
+ -- points
+ local z1 = data[i] [j]
+ local z2 = data[i] [j+1]
+ local z3 = data[i+1][j+1]
+ local z4 = data[i+1][j]
+ -- color
+ local cf = z1[3]
+ if clip then
+ -- best clip here if needed
+ if cf < 0 then
+ cf = 0
+ elseif cf > 1 then
+ cf = 1
+ end
+ else
+ -- or remap when we want to
+ cf = (z1[3] - minf) / range
+ end
+ local z11 = z1[1]
+ local z12 = z1[2]
+ local z21 = z2[1]
+ local z22 = z2[2]
+ local z31 = z3[1]
+ local z32 = z3[2]
+ local z41 = z4[1]
+ local z42 = z4[2]
+ -- if lines then
+ -- -- fill first and draw then, previous shapes can be covered
+ -- else
+ -- -- fill and draw in one go to prevent artifacts
+ -- end
+ local cr, cg, cb = color(cf)
+ if not cr then cr = 0 end
+ if not cg then cg = 0 end
+ if not cb then cb = 0 end
+ if enforce then
+ cr, cg, cb = 1 - cr, 1 - cg, 1 - cb
+ r = r + 1
+ if lines then
+ result[r] = f_fill_cmy(z11,z12,z21,z22,z31,z32,z41,z42,cr,cg,cb)
+ r = r + 1
+ result[r] = f_draw_cmy(z11,z12,z21,z22,z31,z32,z41,z42,cl)
+ else
+ result[r] = f_mesh_cmy(z11,z12,z21,z22,z31,z32,z41,z42,cr,cg,cb)
+ end
+ else
+ r = r + 1
+ if lines then
+ result[r] = f_fill_rgb(z11,z12,z21,z22,z31,z32,z41,z42,cr,cg,cb)
+ r = r + 1
+ result[r] = f_draw_rgb(z11,z12,z21,z22,z31,z32,z41,z42,cl)
+ else
+ result[r] = f_mesh_rgb(z11,z12,z21,z22,z31,z32,z41,z42,cr,cg,cb)
+ end
+ end
+ end
+ end
+ mp.direct(concat(result))
+end
generated by cgit v1.2.3 (git 2.25.1) at 2024-04-29 10:06:34 +0000