From f1129626606384a7a55a21a83531f51f8b5dee25 Mon Sep 17 00:00:00 2001 From: Hans Hagen Date: Tue, 14 Jul 2020 00:25:53 +0200 Subject: 2020-07-13 23:52:00 --- tex/context/base/mkiv/mlib-cnt.lmt | 2022 ++++++++++++++++++++++++++++++++++++ 1 file changed, 2022 insertions(+) create mode 100644 tex/context/base/mkiv/mlib-cnt.lmt (limited to 'tex/context/base/mkiv/mlib-cnt.lmt') 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 -- cgit v1.2.3