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 starttiming = statistics.starttiming local stoptiming = statistics.stoptiming local resettiming = statistics.resettiming local elapsedtiming = statistics.elapsed 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() resettiming("lmt_contours") 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] = 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() 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, elapsedtiming("lmt_contours"), "[ " .. 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 -- ]] ] local f_color = formatters [ [[ local math = math local min = math.min local max = math.max local abs = math.abs 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(f,z) brightness_factor = f function_value = z 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, minz, maxz -- -- 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 -- 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 -- if not minz then minz = z maxz = z elseif z < minz then minz = z elseif z > maxz then maxz = z end data[i][j] = { x, y, fac, z } -- 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" local ccode = load(f_color(minz,maxz,preamble,colorcode)) local color = type(ccode) == "function" and ccode() if type(color) ~= "function" then return false -- fatal error end 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,z1[4]) -- cf, zout 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