if not modules then modules = { } end modules ['mlib-cnt'] = {
version = 1.001,
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
