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