SourceBrowser

1if
 not
 modules
 then
 modules
 =
 {
 }
 end
 modules
 [
'
mlib-cnt
'
]
 =
 {
2    version
   =
 1
.
001
,
3    comment
   =
 "
companion to mlib-ctx.mkiv
"
,
4    author
    =
 "
Hans Hagen, PRAGMA-ADE, Hasselt NL
"
,
5    copyright
 =
 "
PRAGMA ADE / ConTeXt Development Team
"
,
6    license
   =
 "
see context related readme files
"
,
7}
8
9-- The only useful reference that I could find about this topic is in wikipedia:
10--
11--     https://en.wikipedia.org/wiki/Marching_squares
12--
13-- I derived the edge code from:
14--
15--     https://physiology.arizona.edu/people/secomb/contours
16--
17-- and also here:
18--
19--     https://github.com/secomb/GreensV4
20--
21-- which has the banner:
22--
23--     TWS, November 1989. Converted to C September 2007. Revised February 2009.
24--
25-- and has a liberal licence. I optimized the code so that it runs a bit faster in Lua and
26-- there's probably more room for optimization (I tried several variants). For instance I
27-- don't use that many tables because access is costly. We don't have a compiler that does
28-- some optimizing (even then the c code can probably be made more efficient).
29--
30-- We often have the same node so it's cheaper to reuse the wsp tables and reconstruct
31-- the point in the path then to alias the original point. We can be more clever:
32-- straighten, but it's more work so maybe I will do that later; for now I only added
33-- a test for equality. There are some experiments in this file too and not all might
34-- work. It's a playground for me.
35--
36-- The code is meant for metafun so it is not general purpose. The bitmap variant is
37-- relatively fast and bitmaps compress well. When all is settled I might make a couple
38-- of helpers in C for this purpose but not now.
39--
40-- I removed optimization code (more aggressive flattening and such because it didn't really
41-- pay off now that we use combined paths with just line segments. I also moved some other
42-- code to a experimental copy. So we now only have the bare helpers needed here.
43
44-- Todo: look into zero case (lavel 1) for shapes ... omiting the background calculation
45-- speeds up quite a bit.
46
47local
 next
,
 type
,
 tostring
 =
 next
,
 type
,
 tostring
48local
 round
,
 abs
,
 min
,
 max
,
 floor
 =
 math
.
round
,
 math
.
abs
,
 math
.
min
,
 math
.
max
,
 math
.
floor
49local
 concat
,
 move
 =
 table
.
concat
,
 table
.
move
50
51local
 bor
 =
 bit32
.
bor
 -- it's really time to ditch support for luajit
52
53local
 starttiming
       =
 statistics
.
starttiming
54local
 stoptiming
        =
 statistics
.
stoptiming
55local
 elapsedtime
       =
 statistics
.
elapsedtime
56
57local
 formatters
        =
 string
.
formatters
58local
 setmetatableindex
 =
 table
.
setmetatableindex
59local
 sortedkeys
        =
 table
.
sortedkeys
60local
 sequenced
         =
 table
.
sequenced
61
62local
 metapost
          =
 metapost
 or
 {
 }
63local
 mp
                =
 mp
 or
 {
 }
64
65local
 getparameterset
   =
 metapost
.
getparameterset
66
67local
 mpflush
           =
 mp
.
flush
68local
 mpcolor
           =
 mp
.
color
69local
 mpstring
          =
 mp
.
string
70local
 mpdraw
            =
 mp
.
draw
71local
 mpfill
            =
 mp
.
fill
72local
 mpflatten
         =
 mp
.
flatten
73
74local
 report
            =
 logs
.
reporter
(
"
mfun contour
"
)
75
76local
 version
           =
 0
.
007
77
78-- we iterate using integers so that we get a better behaviour at zero
79
80local
 f_function_n
 =
 formatters
 [
 [[

81    local math  = math
82    local round = math.round
83    %s
84    return function(data,nx,ny,nxmin,nxmax,xstep,nymin,nymax,ystep)
85        local sx = nxmin
86        for mx=1,nx do
87            local dx = data[mx]
88            local x  = sx * xstep
89            local sy = nymin
90            for my=1,ny do
91                local y  = sy * ystep
92                dx[my] = %s
93                sy = sy + 1
94            end
95            sx = sx + 1
96        end
97        return 0
98    end
99
]]
 ]
100
101local
 f_function_y
 =
 formatters
 [
 [[

102    local math  = math
103    local round = math.round
104    local nan   = NaN
105    local inf   = math.huge
106    %s
107    return function(data,nx,ny,nxmin,nxmax,xstep,nymin,nymax,ystep,dnan,dinf,report)
108        local sx = nxmin
109        local er = 0
110        for mx=nxmin,nxmax do
111            local dx = data[mx]
112            local x  = sx * xstep
113            local sy = nymin
114            for my=nymin,nymax do
115                local y = sy * ystep
116                local n = %s
117                if n == nan then
118                    er = er + 1
119                    if er < 10 then
120                        report("nan at (%s,%s)",x,y)
121                    end
122                    n = dnan
123                elseif n == inf then
124                    er = er + 1
125                    if er < 10 then
126                        report("inf at (%s,%s)",x,y)
127                    end
128                    n = dinf
129                end
130                dx[my] = n
131                sy = sy + 1
132            end
133            sx = sx + 1
134        end
135        return er
136    end
137
]]
 ]
138
139local
 f_color
 =
 formatters
 [
 [[

140    local math = math
141    local min  = math.min
142    local max  = math.max
143    local n    = %s
144    local minz = %s
145    local maxz = %s
146
147    local color_value = 0
148    local color_step  = mp.lmt_color_functions.step
149    local color_shade = mp.lmt_color_functions.shade
150
151    local function step(...)
152        return color_step(color_value,n,...)
153    end
154    local function shade(...)
155        return color_shade(color_value,n,...)
156    end
157    local function lin(l)
158        return l/n
159    end
160    %s
161    return function(l)
162        color_value = l
163        return %s
164    end
165
]]
 ]
166
167local
 inbetween
 =
 attributes
 and
 attributes
.
colors
.
helpers
.
inbetween
168
169mp
.
lmt_color_functions
 =
 {
170    step
 =
 function
(
l
,
n
,
r
,
g
,
b
,
s
)
171        if
 not
 s
 then
172            s
 =
 1
173        end
174        local
 f
 =
 l
 /
 n
175        local
 r
 =
 s
 *
 1
.
5
 -
 4
 *
 abs
(
f
-
r
)
176        local
 g
 =
 s
 *
 1
.
5
 -
 4
 *
 abs
(
f
-
g
)
177        local
 b
 =
 s
 *
 1
.
5
 -
 4
 *
 abs
(
f
-
b
)
178        return
 min
(
max
(
r
,
0
)
,
1
)
,
 min
(
max
(
g
,
0
)
,
1
)
,
 min
(
max
(
b
,
0
)
,
1
)
179    end
,
180    shade
 =
 function
(
l
,
n
,
one
,
two
)
181        local
 f
 =
 l
 /
 n
182        local
 r
 =
 inbetween
(
one
,
two
,
1
,
f
)
183        local
 g
 =
 inbetween
(
one
,
two
,
2
,
f
)
184        local
 b
 =
 inbetween
(
one
,
two
,
3
,
f
)
185        return
 min
(
max
(
r
,
0
)
,
1
)
,
 min
(
max
(
g
,
0
)
,
1
)
,
 min
(
max
(
b
,
0
)
,
1
)
186    end
,
187}
188
189local
 f_box
 =
 formatters
 [
 [[
draw rawtexbox("contour",%s) xysized (%s,%s) ;
]]
 ]
190
191local
 n_box
 =
 0
192
193-- todo: remove old one, so we need to store old hashes
194
195local
 nofcontours
 =
 0
196
197-- We don't want cosmetics like axis and labels to trigger a  calculation,
198-- especially a slow one.
199
200local
 hashfields
 =
 {
201    "
xmin
"
,
 "
xmax
"
,
 "
xstep
"
,
 "
ymin
"
,
 "
ymax
"
,
 "
ystep
"
,
202    "
levels
"
,
 "
colors
"
,
 "
level
"
,
 "
preamble
"
,
 "
function
"
,
 "
functions
"
,
 "
color
"
,
 "
values
"
,
203    "
background
"
,
 "
foreground
"
,
 "
linewidth
"
,
 "
backgroundcolor
"
,
 "
linecolor
"
,
204}
205
206local
 function
 prepare
(
p
)
207    local
 h
 =
 {
 }
208    for
 i
=
1
,
#
hashfields
 do
209        local
 f
 =
 hashfields
[
i
]
210        h
[
f
]
 =
 p
[
f
]
211    end
212    local
 hash
 =
 md5
.
HEX
(
sequenced
(
h
)
)
213    local
 name
 =
 formatters
[
"
%s-m-c-%03i.lua
"
]
(
tex
.
jobname
,
nofcontours
)
214    return
 name
,
 hash
215end
216
217local
 function
 getcache
(
p
)
218    local
 cache
 =
 p
.
cache
219    if
 cache
 then
220        local
 name
,
 hash
 =
 prepare
(
p
)
221        local
 data
 =
 table
.
load
(
name
)
222        if
 data
 and
 data
.
hash
 =
=
 hash
 and
 data
.
version
 =
=
 version
 then
223            p
.
result
 =
 data
224            return
 true
225        else
226            return
 false
,
 hash
,
 name
227        end
228    end
229    return
 false
,
 nil
,
 nil
230end
231
232local
 function
 setcache
(
p
)
233    local
 result
 =
 p
.
result
234    local
 name
   =
 result
.
name
235    if
 name
 then
236        if
 result
.
bitmap
 then
237            result
.
bitmap
 =
 nil
238        else
239            result
.
data
 =
 nil
240        end
241        result
.
color
  =
 nil
242        result
.
action
 =
 nil
243        result
.
cached
 =
 nil
244        io
.
savedata
(
name
,
 table
.
fastserialize
(
result
)
)
245    else
246        os
.
remove
(
(
prepare
(
p
)
)
)
247    end
248end
249
250function
 mp
.
lmt_contours_start
(
)
251
252    starttiming
(
"
lmt_contours
"
)
253
254    nofcontours
 =
 nofcontours
 +
 1
255
256    local
 p
 =
 getparameterset
(
)
257
258    local
 xmin
    =
 p
.
xmin
259    local
 xmax
    =
 p
.
xmax
260    local
 ymin
    =
 p
.
ymin
261    local
 ymax
    =
 p
.
ymax
262    local
 xstep
   =
 p
.
xstep
263    local
 ystep
   =
 p
.
ystep
264    local
 levels
  =
 p
.
levels
265    local
 colors
  =
 p
.
colors
266    local
 nx
      =
 0
267    local
 ny
      =
 0
268    local
 means
   =
 setmetatableindex
(
"
number
"
)
269    local
 values
  =
 setmetatableindex
(
"
number
"
)
270    local
 data
    =
 setmetatableindex
(
"
table
"
)
271    local
 minmean
 =
 false
272    local
 maxmean
 =
 false
273
274    local
 cached
,
 hash
,
 name
 =
 getcache
(
p
)
275
276    local
 function
 setcolors
(
preamble
,
levels
,
minz
,
maxz
,
color
,
nofvalues
)
277        if
 colors
 then
278            local
 function
 f
(
k
)
279                return
 #
colors
[
1
]
 =
=
 1
 and
 0
 or
 {
 0
,
 0
,
 0
 }
280            end
281            setmetatableindex
(
colors
,
 function
(
t
,
k
)
282                local
 v
 =
 f
(
k
)
283                t
[
k
]
 =
 v
284                return
 v
285            end
)
286            local
 c
 =
 {
 }
287            local
 n
 =
 1
288            for
 i
=
1
,
nofvalues
 do
289                c
[
i
]
 =
 colors
[
n
]
290                n
 =
 n
 +
 1
291            end
292            return
 c
,
 f
293        else
294            local
 tcolor
 =
 f_color
(
levels
,
minz
,
maxz
,
preamble
,
color
)
295            local
 colors
 =
 {
 }
296            local
 fcolor
 =
 load
(
tcolor
)
297            if
 type
(
fcolor
)
 ~
=
 "
function
"
 then
298                report
(
"
error in color function, case %i: %s
"
,
1
,
color
)
299                fcolor
 =
 false
300            else
301                fcolor
 =
 fcolor
(
)
302                if
 type
(
fcolor
)
 ~
=
 "
function
"
 then
303                    report
(
"
error in color function, case %i: %s
"
,
2
,
color
)
304                    fcolor
 =
 false
305                end
306            end
307            if
 not
 fcolor
 then
308                local
 color_step
  =
 mp
.
lmt_color_functions
.
step
309                fcolor
 =
 function
(
l
)
310                    return
 color_step
(
l
,
levels
,
0
.
25
,
0
.
50
,
0
.
75
)
311                end
312            end
313            for
 i
=
1
,
nofvalues
 do
314                colors
[
i
]
 =
 {
 fcolor
(
i
)
 }
315            end
316            if
 attributes
.
colors
.
model
 =
=
 "
cmyk
"
 then
317                for
 i
=
1
,
#
colors
 do
318                    local
 c
 =
 colors
[
i
]
319                    colors
[
i
]
 =
 {
 1
 -
 c
[
1
]
,
 1
 -
 c
[
2
]
,
 1
 -
 c
[
3
]
,
 0
 }
320                end
321            end
322            return
 colors
,
 fcolor
323        end
324    end
325
326    if
 cached
 then
327        local
 result
 =
 p
.
result
328        local
 colors
,
 color
 =
 setcolors
(
329            p
.
preamble
,
330            p
.
levels
,
331            result
.
minz
,
332            result
.
maxz
,
333            p
.
color
,
334            result
.
nofvalues
335        )
336        result
.
color
  =
 color
337        result
.
colors
 =
 colors
338        result
.
cached
 =
 true
339        return
340    end
341
342    local
 functioncode
  =
 p
[
"
function
"
]
343    local
 functionrange
 =
 p
.
range
344    local
 functionlist
  =
 functionrange
 and
 p
.
functions
345    local
 preamble
      =
 p
.
preamble
346
347    if
 xstep
 =
=
 0
 then
 xstep
 =
 (
xmax
 -
 xmin
)
/
100
 end
348    if
 ystep
 =
=
 0
 then
 ystep
 =
 (
ymax
 -
 ymin
)
/
100
 end
349
350    local
 nxmin
 =
 round
(
xmin
/
xstep
)
351    local
 nxmax
 =
 round
(
xmax
/
xstep
)
352    local
 nymin
 =
 round
(
ymin
/
ystep
)
353    local
 nymax
 =
 round
(
ymax
/
ystep
)
354    local
 nx
    =
 nxmax
 -
 nxmin
 +
 1
355    local
 ny
    =
 nymax
 -
 nymin
 +
 1
356
357    local
 function
 executed
(
data
,
code
)
358        local
 fcode
 =
 p
.
check
 and
 f_function_y
 or
 f_function_n
359        fcode
 =
 fcode
(
preamble
,
code
)
360        fcode
 =
 load
(
fcode
)
361        if
 type
(
fcode
)
 =
=
 "
function
"
 then
362            fcode
 =
 fcode
(
)
363        end
364        if
 type
(
fcode
)
 =
=
 "
function
"
 then
365            local
 er
 =
 fcode
(
366                data
,
 nx
,
 ny
,
367                nxmin
,
 nxmax
,
 xstep
,
368                nymin
,
 nymax
,
 ystep
,
369                p
.
defaultnan
,
 p
.
defaultinf
,
 report
370            )
371            if
 er
 >
 0
 then
372                report
(
"
%i errors in: %s
"
,
er
,
code
)
373            end
374            return
 true
375        else
376            return
 false
 -- fatal error
377        end
378    end
379
380 -- for i=1,nx do
381 --     data[i] = lua.newtable(ny,0)
382 -- end
383
384    if
 functionlist
 then
385
386        local
 datalist
 =
 {
 }
387        local
 minr
     =
 functionrange
[
1
]
388        local
 maxr
     =
 functionrange
[
2
]
 or
 minr
389        local
 size
     =
 #
functionlist
390
391        for
 l
=
1
,
size
 do
392
393            local
 func
 =
 setmetatableindex
(
"
table
"
)
394            if
 not
 executed
(
func
,
functionlist
[
l
]
)
 then
395                report
(
"
error in executing function %i from functionlist
"
,
l
)
396                return
397            end
398
399            local
 bit
 =
 l
 -- + 1
400
401            if
 l
 =
=
 1
 then
402                for
 i
=
1
,
nx
 do
403                    local
 di
 =
 data
[
i
]
404                    local
 fi
 =
 func
[
i
]
405                    for
 j
=
1
,
ny
 do
406                        local
 f
 =
 fi
[
j
]
407                        if
 f
 >
=
 minr
 and
 f
 <
=
 maxr
 then
408                            di
[
j
]
 =
 bit
409                        else
410                            di
[
j
]
 =
 0
411                        end
412                    end
413                end
414            else
415                for
 i
=
1
,
nx
 do
416                    local
 di
 =
 data
[
i
]
417                    local
 fi
 =
 func
[
i
]
418                    for
 j
=
1
,
ny
 do
419                        local
 f
 =
 fi
[
j
]
420                        if
 f
 >
=
 minr
 and
 f
 <
=
 maxr
 then
421                            di
[
j
]
 =
 bor
(
di
[
j
]
,
bit
)
422                        end
423                    end
424                end
425            end
426
427        end
428
429        -- we can simplify the value mess below
430
431    elseif
 functioncode
 then
432
433        if
 not
 executed
(
data
,
functioncode
)
 then
434            report
(
"
error in executing function
"
)
435            return
 -- fatal error
436        end
437
438    else
439
440        report
(
"
no valid function(s)
"
)
441        return
 -- fatal error
442
443    end
444
445    minz
 =
 data
[
1
]
[
1
]
446    maxz
 =
 minz
447
448    for
 i
=
1
,
nx
 do
449        local
 d
 =
 data
[
i
]
450        for
 j
=
1
,
ny
 do
451            local
 v
 =
 d
[
j
]
452            if
 v
 <
 minz
 then
453                minz
 =
 v
454            elseif
 v
 >
 maxz
 then
455                maxz
 =
 v
456            end
457        end
458    end
459
460    if
 functionlist
 then
461
462        for
 i
=
minz
,
maxz
 do
463            values
[
i
]
 =
 i
464        end
465
466        levels
 =
 maxz
467
468        minmean
 =
 minz
469        maxmean
 =
 maxz
470
471    else
472
473        local
 snap
 =
 (
maxz
 -
 minz
 +
 1
)
 /
 levels
474
475        for
 i
=
1
,
nx
 do
476            local
 d
 =
 data
[
i
]
477            for
 j
=
1
,
ny
 do
478                local
 dj
 =
 d
[
j
]
479                local
 v
  =
 round
(
(
dj
 -
 minz
)
 /
 snap
)
480                values
[
v
]
 =
 values
[
v
]
 +
 1
481                means
 [
v
]
 =
 means
 [
v
]
 +
 dj
482                d
[
j
]
 =
 v
483            end
484        end
485
486        local
 list
  =
 sortedkeys
(
values
)
487        local
 count
 =
 values
488        local
 remap
 =
 {
 }
489
490        values
 =
 {
 }
491
492        for
 i
=
1
,
#
list
 do
493            local
 v
 =
 list
[
i
]
494            local
 m
 =
 means
[
v
]
 /
 count
[
v
]
495            remap
 [
v
]
 =
 i
496            values
[
i
]
 =
 m
497            if
 not
 minmean
 then
498                minmean
 =
 m
499                maxmean
 =
 m
500            elseif
 m
 <
 minmean
 then
501                minmean
 =
 m
502            elseif
 m
 >
 maxmean
 then
503                maxmean
 =
 m
504            end
505        end
506
507        for
 i
=
1
,
nx
 do
508            local
 d
 =
 data
[
i
]
509            for
 j
=
1
,
ny
 do
510                d
[
j
]
 =
 remap
[
d
[
j
]
]
511            end
512        end
513
514    end
515
516    -- due to rounding we have values + 1 so we can have a floor, ceil, round
517    -- option as well as levels -1
518
519    local
 nofvalues
 =
 #
values
520
521    local
 colors
 =
 setcolors
(
522        p
.
preamble
,
levels
,
minz
,
maxz
,
p
.
color
,
nofvalues
523    )
524
525    p
.
result
 =
 {
526        version
   =
 version
,
527        values
    =
 values
,
528        nofvalues
 =
 nofvalues
,
529        minz
      =
 minz
,
530        maxz
      =
 maxz
,
531        minmean
   =
 minmean
,
532        maxmean
   =
 maxmean
,
533        data
      =
 data
,
534        color
     =
 color
,
535        nx
        =
 nx
,
536        ny
        =
 ny
,
537        colors
    =
 colors
,
538        name
      =
 name
,
539        hash
      =
 hash
,
540        islist
    =
 functionlist
 and
 true
 or
 false
,
541    }
542
543    report
(
"
index %i, nx %i, ny %i, levels %i
"
,
 nofcontours
,
 nx
,
 ny
,
 nofvalues
)
544end
545
546function
 mp
.
lmt_contours_stop
(
)
547    local
 p
 =
 getparameterset
(
)
548    local
 e
 =
 stoptiming
(
"
lmt_contours
"
)
549    setcache
(
p
)
550    p
.
result
 =
 nil
551    local
 f
 =
 p
[
"
function
"
]
552    local
 l
 =
 p
.
functions
553    report
(
"
index %i, %0.3f seconds for: %s
"
,
554        nofcontours
,
 e
,
 "
[ 
"
 .
.
 concat
(
l
 or
 {
 f
 }
 ,
"
 ] [ 
"
)
 .
.
 "
 ]
"
555    )
556end
557
558function
 mp
.
lmt_contours_bitmap_set
(
)
559    local
 p
          =
 getparameterset
(
)
560    local
 result
     =
 p
.
result
561
562    local
 values
     =
 result
.
values
563    local
 nofvalues
  =
 result
.
nofvalues
564    local
 rawdata
    =
 result
.
data
565    local
 nx
         =
 result
.
nx
566    local
 ny
         =
 result
.
ny
567    local
 colors
     =
 result
.
colors
568    local
 depth
      =
 #
colors
[
1
]
 -- == 3 and "rgb" or "gray"
569
570    -- i need to figure out this offset of + 1
571
572    local
 bitmap
    =
 graphics
.
bitmaps
.
new
(
573        nx
,
 ny
,
574        (
depth
 =
=
 3
 and
 "
rgb
"
)
 or
 (
depth
 =
=
 4
 and
 "
cmyk
"
)
 or
 "
gray
"
,
575        1
,
 false
,
 true
576    )
577
578    local
 palette
   =
 bitmap
.
index
 or
 {
 }
 -- has to start at 0
579    local
 data
      =
 bitmap
.
data
580    local
 p
         =
 0
581
582    if
 depth
 =
=
 3
 or
 depth
 =
=
 4
 then
583        for
 i
=
1
,
nofvalues
 do
584            local
 c
 =
 colors
[
i
]
585            local
 r
 =
 round
(
(
c
[
1
]
 or
 0
)
 *
 255
)
586            local
 g
 =
 round
(
(
c
[
2
]
 or
 0
)
 *
 255
)
587            local
 b
 =
 round
(
(
c
[
3
]
 or
 0
)
 *
 255
)
588            local
 k
 =
 depth
 =
=
 4
 and
 round
(
(
c
[
4
]
 or
 0
)
 *
 255
)
 or
 nil
589            palette
[
p
]
 =
 {
590                (
r
 >
 255
 and
 255
)
 or
 (
r
 <
 0
 and
 0
)
 or
 r
,
591                (
g
 >
 255
 and
 255
)
 or
 (
g
 <
 0
 and
 0
)
 or
 g
,
592                (
b
 >
 255
 and
 255
)
 or
 (
b
 <
 0
 and
 0
)
 or
 b
,
593                k
594            }
595            p
 =
 p
 +
 1
596        end
597    else
598        for
 i
=
1
,
nofvalues
 do
599            local
 s
 =
 colors
[
i
]
[
1
]
600            local
 s
 =
 round
(
(
s
 or
 0
)
 *
 255
)
601            palette
[
p
]
 =
 (
602                (
s
 >
 255
 and
 255
)
 or
 (
s
 <
 0
 and
 0
)
 or
 s
603            )
604            p
 =
 p
 +
 1
605        end
606    end
607
608    -- As (1,1) is the left top corner so we need to flip of we start in
609    -- the left bottom (we cannot loop reverse because we want a properly
610    -- indexed table.
611
612    local
 k
 =
 0
613    for
 y
=
ny
,
1
,
-1
 do
614        k
 =
 k
 +
 1
615        local
 d
 =
 data
[
k
]
616        for
 x
=
1
,
nx
 do
617            d
[
x
]
 =
 rawdata
[
x
]
[
y
]
 -
 1
618        end
619    end
620
621    result
.
bitmap
 =
 bitmap
622end
623
624function
 mp
.
lmt_contours_bitmap_get
(
)
625    local
 p
      =
 getparameterset
(
)
626    local
 result
 =
 p
.
result
627    local
 bitmap
 =
 result
.
bitmap
628    local
 box
    =
 nodes
.
hpack
(
graphics
.
bitmaps
.
flush
(
bitmap
)
)
629    n_box
 =
 n_box
 +
 1
630    nodes
.
boxes
.
savenode
(
"
contour
"
,
tostring
(
n_box
)
,
box
)
631    return
 f_box
(
n_box
,
bitmap
.
xsize
,
bitmap
.
ysize
)
632end
633
634function
 mp
.
lmt_contours_cleanup
(
)
635    nodes
.
boxes
.
reset
(
"
contour
"
)
636    n_box
 =
 0
637end
638
639function
 mp
.
lmt_contours_edge_set
(
)
640    local
 p
         =
 getparameterset
(
)
641    local
 result
    =
 p
.
result
642
643    if
 result
.
cached
 then
 return
 end
644
645    local
 values
    =
 result
.
values
646    local
 nofvalues
 =
 result
.
nofvalues
647    local
 data
      =
 result
.
data
648    local
 nx
        =
 result
.
nx
649    local
 ny
        =
 result
.
ny
650
651    local
 xmin
      =
 p
.
xmin
652    local
 xmax
      =
 p
.
xmax
653    local
 ymin
      =
 p
.
ymin
654    local
 ymax
      =
 p
.
ymax
655    local
 xstep
     =
 p
.
xstep
656    local
 ystep
     =
 p
.
ystep
657
658    local
 wsp
       =
 {
 }
659    local
 edges
     =
 {
 }
660
661    for
 value
=
1
,
nofvalues
 do
662
663        local
 iwsp
 =
 0
664        local
 di
   =
 data
[
1
]
665        local
 dc
666        local
 edge
 =
 {
 }
667        local
 e
    =
 0
668        -- the next loop is fast
669        for
 i
=
1
,
nx
 do
670            local
 di1
 =
 data
[
i
+
1
]
671            local
 dij
 =
 di
[
1
]
672            local
 d
   =
 dij
 -
 value
673            local
 dij1
674            for
 j
=
1
,
ny
 do
675                if
 j
 <
 ny
 then
676                    dij1
 =
 di
[
j
+
1
]
677                    dc
 =
 dij1
 -
 value
678                    if
 (
d
 >
=
 0
 and
 dc
 <
 0
)
 or
 (
d
 <
 0
 and
 dc
 >
=
 0
)
 then
679                        iwsp
 =
 iwsp
 +
 1
680                        local
 y
 =
 (
d
 *
 (
j
+
1
)
 -
 dc
 *
 j
)
 /
 (
d
 -
 dc
)
681                        if
 i
 =
=
 1
 then
682                            wsp
[
iwsp
]
 =
 {
 i
,
 y
,
 0
,
 (
i
 +
 (
j
-1
)
*
nx
)
 }
683                        elseif
 i
 =
=
 nx
 then
684                            wsp
[
iwsp
]
 =
 {
 i
,
 y
,
 (
i
 -
 1
 +
 (
j
-1
)
*
nx
)
,
 0
 }
685                        else
686                            local
 jx
 =
 (
i
 +
 (
j
-1
)
*
nx
)
687                            wsp
[
iwsp
]
 =
 {
 i
,
 y
,
 jx
 -
 1
,
 jx
 }
688                        end
689                    end
690                end
691                if
 i
 <
 nx
 then
692                    local
 dc
 =
 di1
[
j
]
 -
 value
693                    if
 (
d
 >
=
 0
 and
 dc
 <
 0
)
 or
 (
d
 <
 0
 and
 dc
 >
=
 0
)
 then
694                        iwsp
 =
 iwsp
 +
 1
695                        local
 x
 =
 (
d
 *
 (
i
+
1
)
 -
 dc
 *
 i
)
 /
 (
d
 -
 dc
)
696                        if
 j
 =
=
 1
 then
697                            wsp
[
iwsp
]
 =
 {
 x
,
 j
,
 0
,
 (
i
 +
 (
j
-1
)
*
nx
)
 }
698                        elseif
 j
 =
=
 ny
 then
699                            wsp
[
iwsp
]
 =
 {
 x
,
 j
,
 (
i
 +
 (
j
-2
)
*
nx
)
,
 0
 }
700                        else
701                            local
 jx
 =
 i
 +
 (
j
-1
)
*
nx
702                            wsp
[
iwsp
]
 =
 {
 x
,
 j
,
 jx
 -
 nx
,
 jx
 }
703                        end
704                    end
705                end
706                dij
 =
 dij1
707                d
   =
 dc
708            end
709            di
 =
 di1
710        end
711        -- the next loop takes time
712        for
 i
=
1
,
iwsp
 do
713            local
 wspi
 =
 wsp
[
i
]
714            for
 isq
=
3
,
4
 do
715                local
 nsq
 =
 wspi
[
isq
]
716                if
 nsq
 ~
=
 0
 then
717                    local
 px
 =
 wspi
[
1
]
718                    local
 py
 =
 wspi
[
2
]
719                    local
 p
  =
 {
 px
,
 py
 }
720                    local
 pn
 =
 2
721                    wspi
[
isq
]
 =
 0
722                    while
 true
 do
723                        for
 ii
=
1
,
iwsp
 do
724                            local
 w
 =
 wsp
[
ii
]
725                            local
 n1
 =
 w
[
3
]
726                            local
 n2
 =
 w
[
4
]
727                            if
 n1
 =
=
 nsq
 then
728                                local
 x
 =
 w
[
1
]
729                                local
 y
 =
 w
[
2
]
730                                if
 x
 ~
=
 px
 or
 y
 ~
=
 py
 then
731                                    pn
    =
 pn
 +
 1
732                                    p
[
pn
]
 =
 x
733                                    pn
    =
 pn
 +
 1
734                                    p
[
pn
]
 =
 y
735                                    px
    =
 x
736                                    py
    =
 y
737                                end
738                                nsq
   =
 n2
739                                w
[
3
]
  =
 0
740                                w
[
4
]
  =
 0
741                                if
 nsq
 =
=
 0
 then
742                                    if
 pn
 =
=
 1
 then
743                                        pn
    =
 pn
 +
 1
744                                        p
[
pn
]
 =
 w
745                                    end
746                                    goto
 flush
747                                end
748                            elseif
 n2
 =
=
 nsq
 then
749                                local
 x
 =
 w
[
1
]
750                                local
 y
 =
 w
[
2
]
751                                if
 x
 ~
=
 px
 or
 y
 ~
=
 py
 then
752                                    pn
    =
 pn
 +
 1
753                                    p
[
pn
]
 =
 x
754                                    pn
    =
 pn
 +
 1
755                                    p
[
pn
]
 =
 y
756                                    px
    =
 x
757                                    py
    =
 y
758                                end
759                                nsq
   =
 n1
760                                w
[
3
]
  =
 0
761                                w
[
4
]
  =
 0
762                                if
 nsq
 =
=
 0
 then
763                                    goto
 flush
764                                end
765                            end
766                        end
767                    end
768                ::
flush
::
769                    e
 =
 e
 +
 1
770                    edge
[
e
]
 =
 p
771                    if
 mpflatten
 then
772                        mpflatten
(
p
)
773                    end
774                end
775            end
776        end
777
778
779        edges
[
value
]
 =
 edge
780
781    end
782
783    result
.
edges
 =
 edges
784
785end
786
787function
 mp
.
lmt_contours_shade_set
(
edgetoo
)
788    local
 p
        =
 getparameterset
(
)
789    local
 result
   =
 p
.
result
790
791    if
 result
.
cached
 then
 return
 end
792
793    local
 values
    =
 result
.
values
794    local
 nofvalues
 =
 result
.
nofvalues
795    local
 data
      =
 result
.
data
796    local
 nx
        =
 result
.
nx
797    local
 ny
        =
 result
.
ny
798    local
 color
     =
 result
.
color
799
800    local
 edges
     =
 setmetatableindex
(
"
table
"
)
801    local
 shades
    =
 setmetatableindex
(
"
table
"
)
802
803    local
 sqtype
    =
 setmetatableindex
(
"
table
"
)
804
805    local
 xspoly
    =
 {
 0
,
 0
,
 0
,
 0
,
 0
,
 0
 }
806    local
 yspoly
    =
 {
 0
,
 0
,
 0
,
 0
,
 0
,
 0
 }
807    local
 xrpoly
    =
 {
 }
808    local
 yrpoly
    =
 {
 }
809
810    local
 xrpoly
    =
 {
 }
 -- lua.newtable(2000,0)
811    local
 yrpoly
    =
 {
 }
 -- lua.newtable(2000,0)
812
813 -- for i=1,2000 do
814 --     xrpoly[i] = 0
815 --     yrpoly[i] = 0
816 -- end
817
818    -- Unlike a c compiler lua will not optimize loops to run in parallel so we expand
819    -- some of the loops and make sure we don't calculate when not needed. Not that nice
820    -- but not that bad either. Maybe I should just write this from scratch.
821
822--     local i = 0
823--     local j = 0
824
825    -- Analyze each rectangle separately. Overwrite lower colors
826
827    -- Unrolling the loops and copying code and using constants is faster and doesn't
828    -- produce much more code in the end, also because we then can leave out the not
829    -- seen branches. One can argue about the foundit2* blobs but by stepwise optimizing
830    -- this is the result.
831
832    shades
[
1
]
 =
 {
 {
 0
,
 0
,
 nx
 -
 1
,
 0
,
 nx
 -
 1
,
 ny
 -
 1
,
 0
,
 ny
 -
 1
 }
 }
833    edges
 [
1
]
 =
 {
 {
 }
 }
834
835    -- this is way too slow ... i must have messed up some loop .. what is this with value 1
836
837    for
 value
=
1
,
nofvalues
 do
838--     for value=2,nofvalues do
839
840        local
 edge
  =
 {
 }
841        local
 nofe
  =
 0
842        local
 shade
 =
 {
 }
843        local
 nofs
  =
 0
844
845        for
 i
=
1
,
nx
-1
 do
846            local
 s
 =
 sqtype
[
i
]
847            for
 j
=
1
,
ny
-1
 do
848                s
[
j
]
 =
 0
849            end
850        end
851
852        local
 nrp
 =
 0
853
854        local
 function
 addedge
(
a
,
b
,
c
,
d
)
855            nofe
 =
 nofe
 +
 1
 edge
[
nofe
]
 =
 a
856            nofe
 =
 nofe
 +
 1
 edge
[
nofe
]
 =
 b
857            nofe
 =
 nofe
 +
 1
 edge
[
nofe
]
 =
 c
858            nofe
 =
 nofe
 +
 1
 edge
[
nofe
]
 =
 d
859        end
860        while
 true
 do
861            -- search for a square of type 0 with >= 1 corner above contour level
862            local
 i
863            local
 j
864            local
 d0
 =
 data
[
1
]
865            local
 d1
 =
 data
[
2
]
866            for
 ii
=
1
,
nx
 do
867                local
 s
 =
 sqtype
[
ii
]
868                for
 jj
=
1
,
ny
 do
869                    if
 s
[
jj
]
 =
=
 0
 then
870                        if
 d0
[
jj
]
 >
 value
 then
 i
 =
 ii
 j
 =
 jj
 goto
 foundit
 end
871                        if
 d1
[
jj
]
 >
 value
 then
 i
 =
 ii
 j
 =
 jj
 goto
 foundit
 end
872                        local
 j1
 =
 jj
 +
 1
873                        if
 d1
[
j1
]
 >
 value
 then
 i
 =
 ii
 j
 =
 jj
 goto
 foundit
 end
874                        if
 d0
[
j1
]
 >
 value
 then
 i
 =
 ii
 j
 =
 jj
 goto
 foundit
 end
875                    end
876                end
877                d0
 =
 d1
878                d1
 =
 data
[
ii
+
1
]
879            end
880            break
881        ::
foundit
::
882            -- initialize r-polygon (nrp seems to be 1 or 2)
883            nrp
 =
 nrp
 +
 1
884
885            local
 first
  =
 true
886            local
 nrpoly
 =
 0
887            local
 nspoly
 =
 0
888            local
 nrpm
   =
 -
nrp
889            -- this is the main loop
890            while
 true
 do
891                -- search for a square of type -nrp
892                if
 first
 then
893                    first
 =
 false
894                    if
 sqtype
[
i
]
[
j
]
 =
=
 0
 then
 -- true anyway
895                        goto
 foundit1
896                    end
897                end
898                for
 ii
=
1
,
nx
 do
899                    local
 s
 =
 sqtype
[
ii
]
900                    for
 jj
=
1
,
ny
 do
901                        if
 s
[
jj
]
 =
=
 nrpm
 then
902                            i
 =
 ii
903                            j
 =
 jj
904                            goto
 foundit1
905                        end
906                    end
907                end
908                break
909            ::
foundit1
::
910                while
 true
 do
911
912                     -- search current then neighboring squares for square type 0, with a corner in common with current square above contour level
913
914                    -- top/bottom ... a bit cheating here
915
916                    local
 i_l
,
 i_c
,
 i_r
    -- i left   current right
917                    local
 j_b
,
 j_c
,
 j_t
    -- j bottom current top
918
919                    local
 i_n
 =
 i
 +
 1
      -- i next (right)
920                    local
 j_n
 =
 j
 +
 1
      -- j next (top)
921
922                    local
 i_p
 =
 i
 -
 1
      -- i previous (bottom)
923                    local
 j_p
 =
 j
 -
 1
      -- j previous (right)
924
925                    local
 d_c
 =
 data
[
i
]
926                    local
 d_r
 =
 data
[
i_n
]
927
928                    local
 sq
929
930                    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
931                        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
932                        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
933                        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
934                        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
935                    end
  end
936
937                    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
938                        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
939                        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
940                    end
 end
941
942                    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
943                        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
944                        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
945                    end
 end
946
947                    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
948                        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
949                        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
950                    end
 end
951
952                    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
953                        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
954                        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
955                    end
 end
956
957                    -- not found
958
959                    sqtype
[
i
]
[
j
]
 =
 nrp
960
961                    break
962
963                    -- define s-polygon for found square (i_c,j_c) - may have up to 6 sides
964
965                ::
foundit21
::
 -- 1 2 3 4
966
967                    sq
[
j_c
]
 =
 nrpm
968
969                    xspoly
[
1
]
 =
 i_l
 ;
 yspoly
[
1
]
 =
 j_b
970                    xspoly
[
2
]
 =
 i_c
 ;
 yspoly
[
2
]
 =
 j_b
971                    if
 d_r
[
j_c
]
 >
 value
 then
 -- dd2
972                        xspoly
[
3
]
 =
 i_c
 ;
 yspoly
[
3
]
 =
 j_c
973                        if
 d_r
[
j_t
]
 >
 value
 then
 -- dd3
974                            xspoly
[
4
]
 =
 i_l
 ;
 yspoly
[
4
]
 =
 j_c
975                            if
 d_c
[
j_t
]
 >
 value
 then
 -- dd4
976                                nspoly
 =
 4
977                            else
978                                xspoly
[
5
]
 =
 i_l
 ;
 yspoly
[
5
]
 =
 j_c
 ;
 nspoly
 =
 5
979                            end
980                        elseif
 d_c
[
j_t
]
 >
 value
 then
 -- dd4
981                            xspoly
[
4
]
 =
 i_c
 ;
 yspoly
[
4
]
 =
 j_c
 ;
982                            xspoly
[
5
]
 =
 i_l
 ;
 yspoly
[
5
]
 =
 j_c
 ;
 nspoly
 =
 5
983                        else
984                            xspoly
[
4
]
 =
 i_l
 ;
 yspoly
[
4
]
 =
 j_c
 ;
 nspoly
 =
 4
985                            if
 edgetoo
 then
 addedge
(
i_c
,
 j_c
,
 i_l
,
 j_c
)
 end
986                        end
987                    elseif
 d_r
[
j_t
]
 >
 value
 then
 -- dd3
988                        xspoly
[
3
]
 =
 i_c
 ;
 yspoly
[
3
]
 =
 j_b
989                        xspoly
[
4
]
 =
 i_c
 ;
 yspoly
[
4
]
 =
 j_c
990                        if
 d_c
[
j_t
]
 >
 value
 then
 -- dd4
991                            xspoly
[
5
]
 =
 i_l
 ;
 yspoly
[
5
]
 =
 j_c
 ;
 nspoly
 =
 5
992                        else
993                            xspoly
[
5
]
 =
 i_l
 ;
 yspoly
[
5
]
 =
 j_c
 ;
994                            xspoly
[
6
]
 =
 i_l
 ;
 yspoly
[
6
]
 =
 j_c
 ;
 nspoly
 =
 6
995                        end
996                    elseif
 d_c
[
j_t
]
 >
 value
 then
 -- dd4
997                        if
 edgetoo
 then
 addedge
(
i_c
,
 j_b
,
 i_c
,
 j_c
)
 end
998                        xspoly
[
3
]
 =
 i_c
 ;
 yspoly
[
3
]
 =
 j_c
 ;
999                        xspoly
[
4
]
 =
 i_l
 ;
 yspoly
[
4
]
 =
 j_c
 ;
 nspoly
 =
 4
1000                    else
1001                        if
 edgetoo
 then
 addedge
(
i_c
,
 j_b
,
 i_l
,
 j_c
)
 end
1002                        xspoly
[
3
]
 =
 i_l
 ;
 yspoly
[
3
]
 =
 j_c
 ;
 nspoly
 =
 3
1003                    end
1004                    goto
 done
1005
1006                ::
foundit22
::
 -- 4 1 2 3
1007
1008                    sq
[
j_c
]
 =
 nrpm
1009
1010                    xspoly
[
1
]
 =
 i_l
 ;
 yspoly
[
1
]
 =
 j_c
1011                    xspoly
[
2
]
 =
 i_l
 ;
 yspoly
[
2
]
 =
 j_b
1012                    if
 d_c
[
j_c
]
 >
 value
 then
 -- dd2
1013                        xspoly
[
3
]
 =
 i_c
 ;
 yspoly
[
3
]
 =
 j_b
1014                        if
 d_r
[
j_c
]
 >
 value
 then
 -- dd3
1015                            xspoly
[
4
]
 =
 i_c
 ;
 yspoly
[
4
]
 =
 j_c
1016                            if
 d_r
[
j_t
]
 >
 value
 then
 -- dd4
1017                                nspoly
 =
 4
1018                            else
1019                                xspoly
[
5
]
 =
 i_c
 ;
 yspoly
[
5
]
 =
 j_c
 ;
 nspoly
 =
 5
 -- suspicious, the same
1020                            end
1021                        elseif
 d_r
[
j_t
]
 >
 value
 then
 -- dd4
1022                            xspoly
[
4
]
 =
 i_c
 ;
 yspoly
[
4
]
 =
 j_b
 ;
1023                            xspoly
[
5
]
 =
 i_c
 ;
 yspoly
[
5
]
 =
 j_c
 ;
 nspoly
 =
 5
1024                        else
1025                            if
 edgetoo
 then
 addedge
(
i_c
,
 j_b
,
 i_c
,
 j_c
)
 end
1026                            xspoly
[
4
]
 =
 i_c
 ;
 yspoly
[
4
]
 =
 j_c
 ;
  nspoly
 =
 4
1027                        end
1028                    elseif
 d_r
[
j_c
]
 >
 value
 then
 -- dd3
1029                        xspoly
[
3
]
 =
 i_l
 ;
 yspoly
[
3
]
 =
 j_b
1030                        xspoly
[
4
]
 =
 i_c
 ;
 yspoly
[
4
]
 =
 j_b
1031                        xspoly
[
5
]
 =
 i_c
 ;
 yspoly
[
5
]
 =
 j_c
1032                        if
 d_r
[
j_t
]
 >
 value
 then
 -- dd4
1033                            nspoly
 =
 5
1034                        else
1035                            xspoly
[
6
]
 =
 i_c
 ;
 yspoly
[
6
]
 =
 j_c
 ;
 nspoly
 =
 6
1036                        end
1037                    elseif
 d_r
[
j_t
]
 >
 value
 then
 -- dd4
1038                        if
 edgetoo
 then
 addedge
(
i_l
,
 j_b
,
 i_c
,
 j_b
)
 end
1039                        xspoly
[
3
]
 =
 i_c
 ;
 yspoly
[
3
]
 =
 j_b
1040                        xspoly
[
4
]
 =
 i_c
 ;
 yspoly
[
4
]
 =
 j_c
 ;
 nspoly
 =
 4
1041                    else
1042                        if
 edgetoo
 then
 addedge
(
i_l
,
 j_b
,
 i_c
,
 j_c
)
 end
1043                        xspoly
[
3
]
 =
 i_c
 ;
 yspoly
[
3
]
 =
 j_c
 ;
 nspoly
 =
 3
1044                    end
1045                    goto
 done
1046
1047                ::
foundit23
::
 -- 2 3 4 1
1048
1049                    sq
[
j_c
]
 =
 nrpm
1050
1051                    xspoly
[
1
]
 =
 i_c
 ;
 yspoly
[
1
]
 =
 j_b
1052                    xspoly
[
2
]
 =
 i_c
 ;
 yspoly
[
2
]
 =
 j_c
1053                    if
 d_r
[
j_t
]
 >
 value
 then
 -- dd2
1054                        xspoly
[
3
]
 =
 i_l
 ;
 yspoly
[
3
]
 =
 j_c
1055                        if
 d_c
[
j_t
]
 >
 value
 then
 -- dd3
1056                            xspoly
[
4
]
 =
 i_l
 ;
 yspoly
[
4
]
 =
 j_b
1057                            if
 d_c
[
j_c
]
 >
 value
 then
 -- dd4
1058                                nspoly
 =
 4
1059                            else
1060                                xspoly
[
5
]
 =
 i_l
 ;
 yspoly
[
5
]
 =
 j_b
 ;
 nspoly
 =
 5
1061                            end
1062                        elseif
 d_c
[
j_c
]
 >
 value
 then
 -- dd4
1063                            xspoly
[
4
]
 =
 i_l
 ;
 yspoly
[
4
]
 =
 j_c
1064                            xspoly
[
5
]
 =
 i_l
 ;
 yspoly
[
5
]
 =
 j_b
 ;
 nspoly
 =
 5
1065                        else
1066                            if
 edgetoo
 then
 addedge
(
i_l
,
 j_c
,
 i_l
,
 j_b
)
 end
1067                            xspoly
[
4
]
 =
 i_l
 ;
 yspoly
[
4
]
 =
 j_b
 ;
 nspoly
 =
 4
1068                        end
1069                    elseif
 d_c
[
j_t
]
 >
 value
 then
 -- dd3
1070                        xspoly
[
3
]
 =
 i_c
 ;
 yspoly
[
3
]
 =
 j_c
1071                        xspoly
[
4
]
 =
 i_l
 ;
 yspoly
[
4
]
 =
 j_c
1072                        xspoly
[
5
]
 =
 i_l
 ;
 yspoly
[
5
]
 =
 j_b
1073                        if
 d_c
[
j_c
]
 >
 value
 then
 -- dd4
1074                            nspoly
 =
 5
1075                        else
1076                            xspoly
[
6
]
 =
 i_l
 ;
 yspoly
[
6
]
 =
 j_b
 ;
 nspoly
 =
 6
1077                        end
1078                    elseif
 d_c
[
j_c
]
 >
 value
 then
 -- dd4
1079                        if
 edgetoo
 then
 addedge
(
i_c
,
 j_c
,
 i_l
,
 j_c
)
 end
1080                        xspoly
[
3
]
 =
 i_l
 ;
 yspoly
[
3
]
 =
 j_c
 ;
1081                        xspoly
[
4
]
 =
 i_l
 ;
 yspoly
[
4
]
 =
 j_b
 ;
 nspoly
 =
 4
1082                    else
1083                        if
 edgetoo
 then
 addedge
(
i_c
,
 j_c
,
 i_l
,
 j_b
)
 end
1084                        xspoly
[
3
]
 =
 i_l
 ;
 yspoly
[
3
]
 =
 j_b
 ;
 nspoly
 =
 3
1085                    end
1086                    goto
 done
1087
1088                ::
foundit24
::
 -- 3 4 1 2
1089
1090                    sq
[
j_c
]
 =
 nrpm
1091
1092                    xspoly
[
1
]
 =
 i_c
 ;
 yspoly
[
1
]
 =
 j_c
1093                    xspoly
[
2
]
 =
 i_l
 ;
 yspoly
[
2
]
 =
 j_c
1094                    if
 d_c
[
j_t
]
 >
 value
 then
 -- dd2
1095                        if
 d_c
[
j_c
]
 >
 value
 then
 -- dd3
1096                            xspoly
[
3
]
 =
 i_l
 ;
 yspoly
[
3
]
 =
 j_b
1097                            xspoly
[
4
]
 =
 i_c
 ;
 yspoly
[
4
]
 =
 j_b
1098                            if
 d_r
[
j_c
]
 >
 value
 then
 -- dd4
1099                                nspoly
 =
 4
1100                            else
1101                                xspoly
[
5
]
 =
 i_c
 ;
 yspoly
[
5
]
 =
 j_b
 ;
 nspoly
 =
 5
1102                            end
1103                        else
1104                            xspoly
[
3
]
 =
 i_l
 ;
 yspoly
[
3
]
 =
 j_b
1105                            if
 d_r
[
j_c
]
 >
 value
 then
 -- dd4
1106
1107                                local
 xv34
 =
 (
dd3
*
i_c
-
dd4
*
i_l
)
/
(
dd3
 -
 dd4
)
 -- probably i_l
1108                                print
(
"
4.4 : xv34
"
,
xv34
,
i_c
,
i_l
)
1109
1110                             -- if edgetoo then addedge(i_l, j_b, xv34, j_b) end
1111                                xspoly
[
4
]
 =
 xv34
 ;
 yspoly
[
4
]
 =
 j_b
 ;
1112                                xspoly
[
5
]
 =
 i_c
 ;
 yspoly
[
5
]
 =
 j_b
 ;
 nspoly
 =
 5
1113                            else
1114                                if
 edgetoo
 then
 addedge
(
i_l
,
 j_b
,
 i_c
,
 j_b
)
 end
1115                                xspoly
[
4
]
 =
 i_c
 ;
 yspoly
[
4
]
 =
 j_b
 ;
 nspoly
 =
 4
1116                            end
1117                        end
1118                    elseif
 d_c
[
j_c
]
 >
 value
 then
 -- dd3
1119                        xspoly
[
3
]
 =
 i_l
 ;
 yspoly
[
3
]
 =
 j_b
1120                        xspoly
[
4
]
 =
 i_l
 ;
 yspoly
[
4
]
 =
 j_b
1121                        xspoly
[
5
]
 =
 i_c
 ;
 yspoly
[
5
]
 =
 j_b
1122                        if
 d_r
[
j_c
]
 >
 value
 then
 -- dd4
1123                            nspoly
 =
 5
1124                        else
1125                            xspoly
[
6
]
 =
 i_c
 ;
 yspoly
[
6
]
 =
 j_b
 ;
 nspoly
 =
 6
1126                        end
1127                    elseif
 d_r
[
j_c
]
 >
 value
 then
 -- dd4
1128                        if
 edgetoo
 then
 addedge
(
i_l
,
 j_c
,
 i_l
,
 j_b
)
 end
1129                        xspoly
[
3
]
 =
 i_l
 ;
 yspoly
[
3
]
 =
 j_b
1130                        xspoly
[
4
]
 =
 i_c
 ;
 yspoly
[
4
]
 =
 j_b
 ;
 nspoly
 =
 4
1131                    else
1132                        if
 edgetoo
 then
 addedge
(
i_l
,
 j_c
,
 i_c
,
 j_b
)
 end
1133                        xspoly
[
3
]
 =
 i_c
 ;
 yspoly
[
3
]
 =
 j_b
 ;
 nspoly
 =
 3
1134                    end
1135                 -- goto done
1136
1137                ::
done
::
1138                    -- combine s-polygon with existing r-polygon, eliminating redundant segments
1139
1140                    if
 nrpoly
 =
=
 0
 then
1141                        -- initiate r-polygon
1142                        for
 i
=
1
,
nspoly
 do
1143                            xrpoly
[
i
]
 =
 xspoly
[
i
]
1144                            yrpoly
[
i
]
 =
 yspoly
[
i
]
1145                        end
1146                        nrpoly
 =
 nspoly
1147                    else
1148                        -- search r-polygon and s-polygon for one side that matches
1149                        --
1150                        -- this is a bottleneck ... we keep this variant here but next go for a faster
1151                        -- alternative
1152                        --
1153                     -- local ispoly, irpoly
1154                     -- for r=nrpoly,1,-1 do
1155                     --     local r1
1156                     --     for s=1,nspoly do
1157                     --         local s1 = s % nspoly + 1
1158                     --         if xrpoly[r] == xspoly[s1] and yrpoly[r] == yspoly[s1] then
1159                     --             if not r1 then
1160                     --                 r1 = r % nrpoly + 1
1161                     --             end
1162                     --             if xrpoly[r1] == xspoly[s] and yrpoly[r1] == yspoly[s] then
1163                     --                 ispoly = s
1164                     --                 irpoly = r
1165                     --                 goto foundit3
1166                     --             end
1167                     --         end
1168                     --     end
1169                     -- end
1170                        --
1171                     -- local ispoly, irpoly
1172                     -- local xr1 = xrpoly[1]
1173                     -- local yr1 = yrpoly[1]
1174                     -- for r0=nrpoly,1,-1 do
1175                     --     for s0=1,nspoly do
1176                     --         if xr1 == xspoly[s0] and yr1 == yspoly[s0] then
1177                     --             if s0 == nspoly then
1178                     --                 if xr0 == xspoly[1] and yr0 == yspoly[1] then
1179                     --                     ispoly = s0
1180                     --                     irpoly = r0
1181                     --                     goto foundit3
1182                     --                 end
1183                     --             else
1184                     --                 local s1 = s0 + 1
1185                     --                 if xr0 == xspoly[s1] and yr0 == yspoly[s1] then
1186                     --                     ispoly = s0
1187                     --                     irpoly = r0
1188                     --                     goto foundit3
1189                     --                 end
1190                     --             end
1191                     --         end
1192                     --     end
1193                     --     xr1 = xrpoly[r0]
1194                     --     yr1 = yrpoly[r0]
1195                     -- end
1196                        --
1197                        -- but ...
1198                        --
1199                        local
 minx
 =
 xspoly
[
1
]
1200                        local
 miny
 =
 yspoly
[
1
]
1201                        local
 maxx
 =
 xspoly
[
1
]
1202                        local
 maxy
 =
 yspoly
[
1
]
1203                        for
 i
=
1
,
nspoly
 do
1204                            local
 y
 =
 yspoly
[
i
]
1205                            if
 y
 <
 miny
 then
1206                                miny
 =
 y
1207                            elseif
 y
 >
 maxy
 then
1208                                maxy
 =
 y
1209                            end
1210                            local
 x
 =
 xspoly
[
i
]
1211                            if
 x
 <
 minx
 then
1212                                minx
 =
 y
1213                            elseif
 x
 >
 maxx
 then
1214                                maxx
 =
 x
1215                            end
1216                        end
1217                        -- we can delay accessing y ...
1218                        local
 ispoly
,
 irpoly
1219                        local
 xr1
 =
 xrpoly
[
1
]
1220                        local
 yr1
 =
 yrpoly
[
1
]
1221                        for
 r0
=
nrpoly
,
1
,
-1
 do
1222                            if
 xr1
 >
=
 minx
 and
 xr1
 <
=
 maxx
 and
 yr1
 >
=
 miny
 and
 yr1
 <
=
 maxy
 then
1223                                local
 xr0
 =
 xrpoly
[
r0
]
1224                                local
 yr0
 =
 yrpoly
[
r0
]
1225                                for
 s0
=
1
,
nspoly
 do
1226                                    if
 xr1
 =
=
 xspoly
[
s0
]
 and
 yr1
 =
=
 yspoly
[
s0
]
 then
1227                                        if
 s0
 =
=
 nspoly
 then
1228                                            if
 xr0
 =
=
 xspoly
[
1
]
 and
 yr0
 =
=
 yspoly
[
1
]
 then
1229                                                ispoly
 =
 s0
1230                                                irpoly
 =
 r0
1231                                                goto
 foundit3
1232                                            end
1233                                        else
1234                                            local
 s1
 =
 s0
 +
 1
1235                                            if
 xr0
 =
=
 xspoly
[
s1
]
 and
 yr0
 =
=
 yspoly
[
s1
]
 then
1236                                                ispoly
 =
 s0
1237                                                irpoly
 =
 r0
1238                                                goto
 foundit3
1239                                            end
1240                                        end
1241                                    end
1242                                end
1243                                xr1
 =
 xr0
1244                                yr1
 =
 yr0
1245                            else
1246                                xr1
 =
 xrpoly
[
r0
]
1247                                yr1
 =
 yrpoly
[
r0
]
1248                            end
1249                        end
1250                        --
1251                        goto
 nomatch3
1252                    ::
foundit3
::
1253                        local
 match1
 =
 0
1254                        local
 rpoly1
 =
 irpoly
 +
 nrpoly
1255                        local
 spoly1
 =
 ispoly
 -
 1
1256                        for
 i
=
2
,
nspoly
-1
 do
1257                            -- search for further matches nearby
1258                            local
 ir
 =
 (
rpoly1
 -
 i
)
 %
 nrpoly
 +
 1
1259                            local
 is
 =
 (
spoly1
 +
 i
)
 %
 nspoly
 +
 1
1260                            if
 xrpoly
[
ir
]
 =
=
 xspoly
[
is
]
 and
 yrpoly
[
ir
]
 =
=
 yspoly
[
is
]
 then
1261                                match1
 =
 match1
 +
 1
1262                            else
1263                                break
 -- goto nomatch1
1264                            end
1265                        end
1266                    ::
nomatch1
::
1267                        local
 match2
 =
 0
1268                        local
 rpoly2
 =
 irpoly
 -
 1
1269                        local
 spoly2
 =
 ispoly
 +
 nspoly
1270                        for
 i
=
2
,
nspoly
-1
 do
1271                            -- search other way for further matches nearby
1272                            local
 ir
 =
 (
rpoly2
 +
 i
)
 %
 nrpoly
 +
 1
1273                            local
 is
 =
 (
spoly2
 -
 i
)
 %
 nspoly
 +
 1
1274                            if
 xrpoly
[
ir
]
 =
=
 xspoly
[
is
]
 and
 yrpoly
[
ir
]
 =
=
 yspoly
[
is
]
 then
1275                                match2
 =
 match2
 +
 1
1276                            else
1277                                break
 -- goto nomatch2
1278                            end
1279                        end
1280                    ::
nomatch2
::
1281                     -- local dnrpoly     = nspoly - 2 - 2*match1 - 2*match2
1282                        local
 dnrpoly
     =
 nspoly
 -
 2
*
(
match1
 +
 match2
 +
 1
)
1283                        local
 ispolystart
 =
 (
ispoly
 +
 match1
)
 %
 nspoly
 +
 1
              -- first node of s-polygon to include
1284                        local
 irpolyend
   =
 (
rpoly1
 -
 match1
 -
 1
)
 %
 nrpoly
 +
 1
          -- last node of s-polygon to include
1285                        if
 dnrpoly
 ~
=
 0
 then
1286                            local
 irpolystart
 =
 (
irpoly
 +
 match2
)
 %
 nrpoly
 +
 1
          -- first node of s-polygon to include
1287                            if
 irpolystart
 >
 irpolyend
 then
1288                             -- local ispolyend = (spoly1 - match2 + nspoly)%nspoly + 1 -- last node of s-polygon to include
1289                                if
 dnrpoly
 >
 0
 then
1290                                    -- expand the arrays xrpoly and yrpoly
1291                                    for
 i
=
nrpoly
,
irpolystart
,
-1
 do
1292                                        local
 k
 =
 i
 +
 dnrpoly
1293                                        xrpoly
[
k
]
 =
 xrpoly
[
i
]
1294                                        yrpoly
[
k
]
 =
 yrpoly
[
i
]
1295                                    end
1296                                else
 -- if dnrpoly < 0 then
1297                                    -- contract the arrays xrpoly and yrpoly
1298                                    for
 i
=
irpolystart
,
nrpoly
 do
1299                                        local
 k
 =
 i
 +
 dnrpoly
1300                                        xrpoly
[
k
]
 =
 xrpoly
[
i
]
1301                                        yrpoly
[
k
]
 =
 yrpoly
[
i
]
1302                                    end
1303                                end
1304                            end
1305                            nrpoly
 =
 nrpoly
 +
 dnrpoly
1306                        end
1307                        if
 nrpoly
 <
 irpolyend
 then
1308                            for
 i
=
irpolyend
,
nrpoly
+
1
,
-1
 do
1309                                -- otherwise these values get lost!
1310                                local
 k
 =
 i
 -
 nrpoly
1311                                xrpoly
[
k
]
 =
 xrpoly
[
i
]
1312                                yrpoly
[
k
]
 =
 yrpoly
[
i
]
1313                            end
1314                        end
1315                        local
 n
 =
 nspoly
 -
 2
 -
 match1
 -
 match2
1316                        if
 n
 =
=
 1
 then
1317                            local
 irpoly1
 =
 irpolyend
   %
 nrpoly
 +
 1
1318                            local
 ispoly1
 =
 ispolystart
 %
 nspoly
 +
 1
1319                            xrpoly
[
irpoly1
]
 =
 xspoly
[
ispoly1
]
1320                            yrpoly
[
irpoly1
]
 =
 yspoly
[
ispoly1
]
1321                        elseif
 n
 >
 0
 then
1322                            -- often 2
1323                            for
 i
=
1
,
n
 do
1324                                local
 ii
 =
 i
 -
 1
1325                                local
 ir
 =
 (
irpolyend
   +
 ii
)
 %
 nrpoly
 +
 1
1326                                local
 is
 =
 (
ispolystart
 +
 ii
)
 %
 nspoly
 +
 1
1327                                xrpoly
[
ir
]
 =
 xspoly
[
is
]
1328                                yrpoly
[
ir
]
 =
 yspoly
[
is
]
1329                            end
1330                        end
1331                ::
nomatch3
::
1332                    end
1333                end
1334            end
1335
1336            if
 nrpoly
 >
 0
 then
1337                local
 t
 =
 {
 }
1338                local
 n
 =
 0
1339                for
 i
=
1
,
nrpoly
 do
1340                    n
 =
 n
 +
 1
 t
[
n
]
 =
 xrpoly
[
i
]
1341                    n
 =
 n
 +
 1
 t
[
n
]
 =
 yrpoly
[
i
]
1342                end
1343                if
 mpflatten
 then
1344                    mpflatten
(
t
)
 -- maybe integrate
1345                end
1346                nofs
 =
 nofs
 +
 1
1347                shade
[
nofs
]
 =
 t
1348             -- print(value,nrpoly,#t,#t-nrpoly*2)
1349            end
1350
1351        end
1352
1353        edges
 [
value
+
1
]
 =
 edge
1354        shades
[
value
+
1
]
 =
 shade
1355--         edges [value] = edge
1356--         shades[value] = shade
1357    end
1358
1359    result
.
shades
 =
 shades
1360    result
.
shapes
 =
 edges
1361
1362end
1363
1364-- accessors
1365
1366function
 mp
.
lmt_contours_nx
       (
i
)
 return
 getparameterset
(
)
.
result
.
nx
 end
1367function
 mp
.
lmt_contours_ny
       (
i
)
 return
 getparameterset
(
)
.
result
.
ny
 end
1368
1369function
 mp
.
lmt_contours_nofvalues
(
)
  return
 getparameterset
(
)
.
result
.
nofvalues
 end
1370function
 mp
.
lmt_contours_value
    (
i
)
 return
 getparameterset
(
)
.
result
.
values
[
i
]
 end
1371
1372function
 mp
.
lmt_contours_minz
     (
i
)
 return
 getparameterset
(
)
.
result
.
minz
 end
1373function
 mp
.
lmt_contours_maxz
     (
i
)
 return
 getparameterset
(
)
.
result
.
maxz
 end
1374
1375function
 mp
.
lmt_contours_minmean
  (
i
)
 return
 getparameterset
(
)
.
result
.
minmean
 end
1376function
 mp
.
lmt_contours_maxmean
  (
i
)
 return
 getparameterset
(
)
.
result
.
maxmean
 end
1377
1378function
 mp
.
lmt_contours_xrange
   (
)
  local
 p
 =
 getparameterset
(
)
 mpstring
(
formatters
[
"
x = [%.3N,%.3N] ;
"
]
(
p
.
xmin
,
p
.
xmax
)
)
 end
1379function
 mp
.
lmt_contours_yrange
   (
)
  local
 p
 =
 getparameterset
(
)
 mpstring
(
formatters
[
"
y = [%.3N,%.3N] ;
"
]
(
p
.
ymin
,
p
.
ymax
)
)
 end
1380
1381function
 mp
.
lmt_contours_format
(
)
1382    local
 p
 =
 getparameterset
(
)
1383    return
 mpstring
(
p
.
result
.
islist
 and
 "
@i
"
 or
 p
.
zformat
 or
 p
.
format
)
1384end
1385
1386function
 mp
.
lmt_contours_function
(
)
1387    local
 p
 =
 getparameterset
(
)
1388    return
 mpstring
(
p
.
result
.
islist
 and
 concat
(
p
[
"
functions
"
]
,
 "
, 
"
)
 or
 p
[
"
function
"
]
)
1389end
1390
1391function
 mp
.
lmt_contours_range
(
)
1392    local
 p
 =
 getparameterset
(
)
1393    local
 r
 =
 p
.
result
.
islist
 and
 p
.
range
1394    if
 not
 r
 or
 #
r
 =
=
 0
 then
1395        return
 mpstring
(
"
"
)
1396    elseif
 #
r
 =
=
 1
 then
1397        return
 mpstring
(
r
[
1
]
)
1398    else
1399        return
 mpstring
(
formatters
[
"
z = [%s,%s]
"
]
(
r
[
1
]
,
r
[
2
]
)
)
1400    end
1401end
1402
1403function
 mp
.
lmt_contours_edge_paths
(
value
)
1404    mpdraw
(
getparameterset
(
)
.
result
.
edges
[
value
]
,
true
)
1405    mpflush
(
)
1406end
1407
1408function
 mp
.
lmt_contours_shape_paths
(
value
)
1409    mpdraw
(
getparameterset
(
)
.
result
.
shapes
[
value
]
,
false
)
1410    mpflush
(
)
1411end
1412
1413function
 mp
.
lmt_contours_shade_paths
(
value
)
1414    mpfill
(
getparameterset
(
)
.
result
.
shades
[
value
]
,
true
)
1415    mpflush
(
)
1416end
1417
1418function
 mp
.
lmt_contours_color
(
value
)
1419    local
 p
     =
 getparameterset
(
)
1420    local
 color
 =
 p
.
result
.
colors
[
value
]
1421    if
 color
 then
1422        mpcolor
(
color
)
1423    end
1424end
1425
1426-- The next code is based on the wikipedia page. It was a bit tedius job to define the
1427-- coordinates but hupefully I made no errors. I rendered all shapes independently and
1428-- tripple checked bit one never knows ...
1429
1430-- maybe some day write from scatch, like this (axis are swapped):
1431
1432local
 d
 =
 1
/
2
1433
1434local
 paths
 =
 {
1435    {
 0
,
 d
,
 d
,
 0
 }
,
1436    {
 1
,
 d
,
 d
,
 0
 }
,
1437    {
 0
,
 d
,
 1
,
 d
 }
,
1438    {
 1
,
 d
,
 d
,
 1
 }
,
1439    {
 0
,
 d
,
 d
,
 1
,
 d
,
 0
,
 1
,
 d
 }
,
 -- saddle
1440    {
 d
,
 0
,
 d
,
 1
 }
,
1441    {
 0
,
 d
,
 d
,
 1
 }
,
1442    {
 0
,
 d
,
 d
,
 1
 }
,
1443    {
 d
,
 0
,
 d
,
 1
 }
,
1444    {
 0
,
 d
,
 d
,
 0
,
 1
,
 d
,
 d
,
 1
 }
,
 -- saddle
1445    {
 1
,
 d
,
 d
,
 1
 }
,
1446    {
 0
,
 d
,
 1
,
 d
 }
,
1447    {
 d
,
 0
,
 1
,
 d
 }
,
1448    {
 d
,
 0
,
 0
,
 d
 }
,
1449}
1450
1451local
 function
 whatever
(
data
,
nx
,
ny
,
threshold
)
1452    local
 edges
 =
 {
 }
1453    local
 e
     =
 0
1454    local
 d0
 =
 data
[
1
]
1455    for
 j
=
1
,
ny
-1
 do
1456        local
 d1
 =
 data
[
j
+
1
]
1457        local
 k
 =
 j
 +
 1
1458        for
 i
=
1
,
nx
-1
 do
1459            local
 v
 =
 0
1460            local
 l
 =
 i
 +
 1
1461            local
 c1
 =
 d0
[
i
]
1462            if
 c1
 <
 threshold
 then
1463                v
 =
 v
 +
 8
1464            end
1465            local
 c2
 =
 d0
[
l
]
1466            if
 c2
 <
 threshold
 then
1467                v
 =
 v
 +
 4
1468            end
1469            local
 c3
 =
 d1
[
l
]
1470            if
 c3
 <
 threshold
 then
1471                v
 =
 v
 +
 2
1472            end
1473            local
 c4
 =
 d1
[
i
]
1474            if
 c4
 <
 threshold
 then
1475                v
 =
 v
 +
 1
1476            end
1477            if
 v
 >
 0
 and
 v
 <
 15
 then
1478                if
 v
 =
=
 5
 or
 v
 =
=
 10
 then
1479                    local
 a
 =
 (
c1
 +
 c2
 +
 c3
 +
 c4
)
 /
 4
1480                    if
 a
 <
 threshold
 then
1481                        v
 =
 v
 =
=
 5
 and
 10
 or
 5
1482                    end
1483                    local
 p
 =
 paths
[
v
]
1484                    e
 =
 e
 +
 1
 edges
[
e
]
 =
 k
 -
 p
[
2
]
1485                    e
 =
 e
 +
 1
 edges
[
e
]
 =
 i
 +
 p
[
1
]
1486                    e
 =
 e
 +
 1
 edges
[
e
]
 =
 k
 -
 p
[
4
]
1487                    e
 =
 e
 +
 1
 edges
[
e
]
 =
 i
 +
 p
[
3
]
1488                    e
 =
 e
 +
 1
 edges
[
e
]
 =
 k
 -
 p
[
6
]
1489                    e
 =
 e
 +
 1
 edges
[
e
]
 =
 i
 +
 p
[
5
]
1490                    e
 =
 e
 +
 1
 edges
[
e
]
 =
 k
 -
 p
[
8
]
1491                    e
 =
 e
 +
 1
 edges
[
e
]
 =
 i
 +
 p
[
7
]
1492                else
1493                    local
 p
 =
 paths
[
v
]
1494                    e
 =
 e
 +
 1
 edges
[
e
]
 =
 k
 -
 p
[
2
]
1495                    e
 =
 e
 +
 1
 edges
[
e
]
 =
 i
 +
 p
[
1
]
1496                    e
 =
 e
 +
 1
 edges
[
e
]
 =
 k
 -
 p
[
4
]
1497                    e
 =
 e
 +
 1
 edges
[
e
]
 =
 i
 +
 p
[
3
]
1498                end
1499            end
1500        end
1501        d0
 =
 d1
1502    end
1503    return
 edges
1504end
1505
1506-- todo: just fetch when needed, no need to cache
1507
1508function
 mp
.
lmt_contours_edge_set_by_cell
(
)
1509    local
 p
         =
 getparameterset
(
)
1510    local
 result
    =
 p
.
result
1511
1512    if
 result
.
cached
 then
 return
 end
1513
1514    local
 values
    =
 result
.
values
1515    local
 nofvalues
 =
 result
.
nofvalues
1516    local
 data
      =
 result
.
data
1517    local
 nx
        =
 result
.
nx
1518    local
 ny
        =
 result
.
ny
1519    local
 lines
     =
 {
 }
1520    result
.
lines
    =
 lines
1521    for
 value
=
1
,
nofvalues
 do
1522        lines
[
value
]
 =
 whatever
(
data
,
ny
,
nx
,
value
)
1523    end
1524end
1525
1526function
 mp
.
lmt_contours_edge_get_cell
(
value
)
1527    mpdraw
(
getparameterset
(
)
.
result
.
lines
[
value
]
)
1528    mpflush
(
)
1529end
1530
1531local
 singles
 =
 {
1532    {
 d
,
 0
,
 0
,
 0
,
 0
,
 d
 }
,
                   --  1  0001
1533    {
 d
,
 0
,
 0
,
 d
 }
,
                         --  2  0002
1534    {
 1
,
 d
,
 1
,
 0
,
 d
,
 0
 }
,
                   --  3  0010
1535    {
 1
,
 d
,
 1
,
 0
,
 0
,
 0
,
 0
,
 d
 }
,
             --  4  0011
1536    {
 1
,
 d
,
 1
,
 0
,
 d
,
 0
,
 0
,
 d
 }
,
             --  5  0012
1537    {
 1
,
 d
,
 d
,
 0
 }
,
                         --  6  0020
1538    {
 1
,
 d
,
 d
,
 0
,
 0
,
 0
,
 0
,
 d
 }
,
             --  7  0021
1539    {
 1
,
 d
,
 0
,
 d
 }
,
                         --  8  0022
1540    {
 d
,
 1
,
 1
,
 1
,
 1
,
 d
 }
,
                   --  9  0100
1541    false
,
                                  -- 10  0101
1542    false
,
                                  -- 11  0102
1543    {
 d
,
 1
,
 1
,
 1
,
 1
,
 0
,
 d
,
 0
 }
,
             -- 12  0110
1544    {
 d
,
 1
,
 1
,
 1
,
 1
,
 0
,
 0
,
 0
,
 0
,
 d
 }
,
       -- 13  0111
1545    {
 d
,
 1
,
 1
,
 1
,
 1
,
 0
,
 d
,
 0
,
 0
,
 d
 }
,
       -- 14  0112
1546    {
 d
,
 1
,
 1
,
 1
,
 1
,
 d
,
 d
,
 0
 }
,
             -- 15  0120
1547    {
 d
,
 1
,
 1
,
 1
,
 1
,
 d
,
 d
,
 0
,
 0
,
 0
,
 0
,
 d
 }
,
 -- 16  0121
1548    {
 d
,
 1
,
 1
,
 1
,
 1
,
 d
,
 0
,
 d
 }
,
             -- 17  0122
1549    {
 d
,
 1
,
 1
,
 d
 }
,
                         -- 18  0200
1550    false
,
                                  -- 19  0201
1551    false
,
                                  -- 20  0202
1552    {
 d
,
 1
,
 1
,
 d
,
 1
,
 0
,
 d
,
 0
 }
,
             -- 21  0210
1553    {
 d
,
 1
,
 1
,
 d
,
 1
,
 0
,
 0
,
 0
,
 0
,
 d
 }
,
       -- 22  0211
1554    false
,
                                  -- 23  0212
1555    {
 d
,
 1
,
 d
,
 0
 }
,
                         -- 24  0220
1556    {
 d
,
 1
,
 d
,
 0
,
 0
,
 0
,
 0
,
 d
 }
,
             -- 25  0221
1557    {
 d
,
 1
,
 0
,
 d
 }
,
                         -- 26  0222
1558    {
 0
,
 1
,
 d
,
 1
,
 0
,
 d
 }
,
                   -- 27  1000
1559    {
 0
,
 1
,
 d
,
 1
,
 d
,
 0
,
 0
,
 0
 }
,
             -- 28  1001
1560    {
 0
,
 1
,
 d
,
 1
,
 d
,
 0
,
 0
,
 d
 }
,
             -- 29  1002
1561    false
,
                                  -- 30  1010
1562    {
 0
,
 1
,
 d
,
 1
,
 1
,
 d
,
 1
,
 0
,
 0
,
 0
 }
,
       -- 31  1011
1563    {
 0
,
 1
,
 d
,
 1
,
 1
,
 d
,
 1
,
 0
,
 d
,
 0
,
 0
,
 d
 }
,
 -- 32  1012
1564    false
,
                                  -- 33  1020
1565    {
 0
,
 1
,
 d
,
 1
,
 1
,
 d
,
 d
,
 0
,
 0
,
 0
 }
,
       -- 34  1021
1566    {
 0
,
 1
,
 d
,
 1
,
 1
,
 d
,
 0
,
 d
 }
,
             -- 35  1022
1567    {
 0
,
 1
,
 1
,
 1
,
 1
,
 d
,
 0
,
 d
 }
,
             -- 36  1100
1568    {
 0
,
 1
,
 1
,
 1
,
 1
,
 d
,
 d
,
 0
,
 0
,
 0
 }
,
       -- 37  1101
1569    {
 0
,
 1
,
 1
,
 1
,
 1
,
 d
,
 d
,
 0
,
 0
,
 d
 }
,
       -- 38  1102
1570    {
 0
,
 1
,
 1
,
 1
,
 1
,
 0
,
 d
,
 0
,
 0
,
 d
 }
,
       -- 39  1110
1571    {
 0
,
 1
,
 1
,
 1
,
 1
,
 0
,
 0
,
 0
 }
,
             -- 40  1111
1572    {
 0
,
 1
,
 1
,
 1
,
 1
,
 0
,
 d
,
 0
,
 0
,
 d
 }
,
       -- 41  1112
1573    {
 0
,
 1
,
 1
,
 1
,
 1
,
 d
,
 d
,
 0
,
 0
,
 d
 }
,
       -- 42  1120
1574    {
 0
,
 1
,
 1
,
 1
,
 1
,
 d
,
 d
,
 0
,
 0
,
 0
 }
,
       -- 43  1121
1575    {
 0
,
 1
,
 1
,
 1
,
 1
,
 d
,
 0
,
 d
 }
,
             -- 44  1122
1576    {
 0
,
 1
,
 d
,
 1
,
 1
,
 d
,
 0
,
 d
 }
,
             -- 45  1200
1577    {
 0
,
 1
,
 d
,
 1
,
 1
,
 d
,
 d
,
 0
,
 0
,
 0
 }
,
       -- 46  1201
1578    false
,
                                  -- 47  1202
1579    {
 0
,
 1
,
 d
,
 1
,
 1
,
 d
,
 1
,
 0
,
 d
,
 0
,
 0
,
 d
 }
,
 -- 48  1210
1580    {
 0
,
 1
,
 d
,
 1
,
 1
,
 d
,
 1
,
 0
,
 0
,
 0
 }
,
       -- 49  1211
1581    false
,
                                  -- 50  1212
1582    {
 0
,
 1
,
 d
,
 1
,
 d
,
 0
,
 0
,
 d
 }
,
             -- 51  1220
1583    {
 0
,
 1
,
 d
,
 1
,
 d
,
 0
,
 0
,
 0
 }
,
             -- 52  1221
1584    {
 0
,
 1
,
 d
,
 1
,
 0
,
 d
 }
,
                   -- 53  1222
1585    {
 d
,
 1
,
 0
,
 d
 }
,
                         -- 54  2000
1586    {
 d
,
 1
,
 d
,
 0
,
 0
,
 0
,
 0
,
 d
 }
,
             -- 55  2001
1587    {
 d
,
 1
,
 d
,
 0
 }
,
                         -- 56  2002
1588    false
,
                                  -- 57  2010
1589    {
 d
,
 1
,
 1
,
 d
,
 1
,
 0
,
 0
,
 0
,
 0
,
 d
 }
,
       -- 58  2011
1590    {
 d
,
 1
,
 1
,
 d
,
 1
,
 0
,
 d
,
 0
 }
,
             -- 59  2012
1591    false
,
                                  -- 60  2020
1592    false
,
                                  -- 61  2021
1593    {
 d
,
 1
,
 1
,
 d
 }
,
                         -- 62  2022
1594    {
 d
,
 1
,
 1
,
 1
,
 1
,
 d
,
 0
,
 d
 }
,
             -- 63  2100
1595    {
 d
,
 1
,
 1
,
 1
,
 1
,
 d
,
 d
,
 0
,
 0
,
 0
,
 0
,
 d
 }
,
 -- 64  2101
1596    {
 d
,
 1
,
 1
,
 1
,
 1
,
 d
,
 d
,
 0
 }
,
             -- 65  2102
1597    {
 d
,
 1
,
 1
,
 1
,
 1
,
 0
,
 d
,
 0
,
 0
,
 d
 }
,
       -- 66  2110
1598    {
 d
,
 1
,
 1
,
 1
,
 1
,
 0
,
 0
,
 0
,
 0
,
 d
 }
,
       -- 67  2111
1599    {
 d
,
 1
,
 1
,
 1
,
 1
,
 0
,
 d
,
 0
 }
,
             -- 68  2112
1600    false
,
                                  -- 69  2120
1601    false
,
                                  -- 70  2121
1602    {
 d
,
 1
,
 1
,
 1
,
 1
,
 d
 }
,
                   -- 71  2122
1603    {
 1
,
 d
,
 0
,
 d
 }
,
                         -- 72  2200
1604    {
 1
,
 d
,
 d
,
 0
,
 0
,
 0
,
 0
,
 d
 }
,
             -- 73  2201
1605    {
 1
,
 d
,
 d
,
 0
 }
,
                         -- 74  2202
1606    {
 1
,
 d
,
 1
,
 0
,
 d
,
 0
,
 0
,
 d
 }
,
             -- 75  2210
1607    {
 1
,
 d
,
 1
,
 0
,
 0
,
 0
,
 0
,
 d
 }
,
             -- 76  2211
1608    {
 1
,
 d
,
 1
,
 0
,
 d
,
 0
 }
,
                   -- 77  2212
1609    {
 d
,
 0
,
 0
,
 d
 }
,
                         -- 78  2220
1610    {
 0
,
 d
,
 0
,
 0
,
 d
,
 0
 }
,
                   -- 79  2221
1611}
1612
1613local
 sadles
 =
 {
1614    false
,
 false
,
 false
,
 false
,
 false
,
 false
,
 false
,
 false
,
 false
,
1615    {
 {
 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
1616    {
 {
 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
1617    false
,
 false
,
 false
,
 false
,
 false
,
 false
,
 false
,
1618    {
 {
 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
1619    {
 {
 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
1620    false
,
 false
,
1621    {
 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
1622    false
,
 false
,
 false
,
 false
,
 false
,
 false
,
1623    {
 {
 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
1624    false
,
 false
,
1625    {
 {
 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
1626    false
,
 false
,
 false
,
 false
,
 false
,
 false
,
 false
,
 false
,
 false
,
 false
,
 false
,
 false
,
 false
,
1627    {
 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
1628    false
,
 false
,
1629    {
 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
1630    false
,
 false
,
 false
,
 false
,
 false
,
 false
,
1631    {
 {
 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
1632    false
,
 false
,
1633    {
 {
 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
1634    {
 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
1635    false
,
 false
,
 false
,
 false
,
 false
,
 false
,
 false
,
1636    {
 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
1637    {
 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
1638}
1639
1640local
 function
 whatever
(
data
,
nx
,
ny
,
threshold
,
background
)
1641
1642    if
 background
 then
1643
1644        local
 llx
 =
 1
/
2
1645        local
 lly
 =
 llx
1646        local
 urx
 =
 ny
 +
 llx
1647        local
 ury
 =
 nx
 +
 lly
1648
1649        return
 {
 {
 llx
,
 lly
,
 urx
,
 0
,
 urx
,
 ury
,
 0
,
 ury
 }
 }
1650
1651    else
1652
1653        local
 bands
 =
 {
 }
1654        local
 b
     =
 0
1655
1656        local
 function
 band
(
s
,
n
,
x
,
y
)
 -- simple. no closure so fast
1657            if
 n
 =
=
 6
 then
1658                return
 {
1659                    x
 -
 s
[
 2
]
,
 y
 +
 s
[
 1
]
,
 x
 -
 s
[
 4
]
,
 y
 +
 s
[
 3
]
,
 x
 -
 s
[
 6
]
,
 y
 +
 s
[
 5
]
,
1660                }
1661            elseif
 n
 =
=
 8
 then
1662                return
 {
1663                    x
 -
 s
[
 2
]
,
 y
 +
 s
[
 1
]
,
 x
 -
 s
[
 4
]
,
 y
 +
 s
[
 3
]
,
 x
 -
 s
[
 6
]
,
 y
 +
 s
[
 5
]
,
1664                    x
 -
 s
[
 8
]
,
 y
 +
 s
[
 7
]
,
1665                }
1666            elseif
 n
 =
=
 10
 then
1667                return
 {
1668                    x
 -
 s
[
 2
]
,
 y
 +
 s
[
 1
]
,
 x
 -
 s
[
 4
]
,
 y
 +
 s
[
 3
]
,
 x
 -
 s
[
 6
]
,
 y
 +
 s
[
 5
]
,
1669                    x
 -
 s
[
 8
]
,
 y
 +
 s
[
 7
]
,
 x
 -
 s
[
10
]
,
 y
 +
 s
[
 9
]
,
1670                }
1671            elseif
 n
 =
=
 4
 then
1672                return
 {
1673                    x
 -
 s
[
 2
]
,
 y
 +
 s
[
 1
]
,
 x
 -
 s
[
 4
]
,
 y
 +
 s
[
 3
]
,
1674                }
1675            else
 -- 12
1676                return
 {
1677                    x
 -
 s
[
 2
]
,
 y
 +
 s
[
 1
]
,
 x
 -
 s
[
 4
]
,
 y
 +
 s
[
 3
]
,
 x
 -
 s
[
 6
]
,
 y
 +
 s
[
 5
]
,
1678                    x
 -
 s
[
 8
]
,
 y
 +
 s
[
 7
]
,
 x
 -
 s
[
10
]
,
 y
 +
 s
[
 9
]
,
 x
 -
 s
[
12
]
,
 y
 +
 s
[
11
]
,
1679                }
1680            end
1681        end
1682
1683        local
 pp
 =
 {
 }
1684
1685        local
 d0
 =
 data
[
1
]
1686        for
 j
=
1
,
ny
-1
 do
1687            local
 d1
 =
 data
[
j
+
1
]
1688            local
 k
 =
 j
 +
 1
1689            local
 p
 =
 false
1690            for
 i
=
1
,
nx
-1
 do
1691                local
 v
 =
 0
1692                local
 l
 =
 i
 +
 1
1693                local
 c1
 =
 d0
[
i
]
1694                if
 c1
 =
=
 threshold
 then
1695                    v
 =
 v
 +
 27
1696                elseif
 c1
 >
 threshold
 then
1697                    v
 =
 v
 +
 54
1698                end
1699                local
 c2
 =
 d0
[
l
]
1700                if
 c2
 =
=
 threshold
 then
1701                    v
 =
 v
 +
 9
1702                elseif
 c2
 >
 threshold
 then
1703                    v
 =
 v
 +
 18
1704                end
1705                local
 c3
 =
 d1
[
l
]
1706                if
 c3
 =
=
 threshold
 then
1707                    v
 =
 v
 +
 3
1708                elseif
 c3
 >
 threshold
 then
1709                    v
 =
 v
 +
 6
1710                end
1711                local
 c4
 =
 d1
[
i
]
1712                if
 c4
 =
=
 threshold
 then
1713                    v
 =
 v
 +
 1
1714                elseif
 c4
 >
 threshold
 then
1715                    v
 =
 v
 +
 2
1716                end
1717                if
 v
 >
 0
 and
 v
 <
 80
 then
1718                    if
 v
 =
=
 40
 then
1719                        -- a little optimization: full areas appended horizontally
1720                        if
 p
 then
1721                            p
[
4
]
 =
 l
 -- i + 1
1722                            p
[
6
]
 =
 l
 -- i + 1
1723                        else
1724                            -- x-0 y+1 x-1 y+1 x-1 y+0 x-0 y+0
1725                            p
 =
 {
 j
,
 i
,
 j
,
 l
,
 k
,
 l
,
 k
,
 i
 }
1726                            b
 =
 b
 +
 1
 ;
 bands
[
b
]
 =
 p
1727                        end
1728                    else
1729                        local
 s
 =
 singles
[
v
]
1730                        if
 s
 then
1731                            b
 =
 b
 +
 1
 ;
 bands
[
b
]
 =
 band
(
s
,
#
s
,
k
,
i
)
1732                        else
1733                            local
 s
 =
 sadles
[
v
]
1734                            if
 s
 then
1735                                local
 m
 =
 (
c1
 +
 c2
 +
 c3
 +
 c4
)
 /
 4
1736                                if
 m
 <
 threshold
 then
1737                                    local
 s1
 =
 s
[
1
]
 if
 s1
 then
 b
 =
 b
 +
 1
 ;
 bands
[
b
]
 =
 band
(
s1
,
#
s1
,
i
,
j
)
 end
1738                                    local
 s2
 =
 s
[
2
]
 if
 s2
 then
 b
 =
 b
 +
 1
 ;
 bands
[
b
]
 =
 band
(
s2
,
#
s2
,
i
,
j
)
 end
1739                                elseif
 m
 =
=
 threshold
 then
1740                                    local
 s3
 =
 s
[
3
]
 if
 s3
 then
 b
 =
 b
 +
 1
 ;
 bands
[
b
]
 =
 band
(
s3
,
#
s3
,
i
,
j
)
 end
1741                                    local
 s4
 =
 s
[
4
]
 if
 s4
 then
 b
 =
 b
 +
 1
 ;
 bands
[
b
]
 =
 band
(
s4
,
#
s4
,
i
,
j
)
 end
1742                                else
1743                                    local
 s5
 =
 s
[
5
]
 if
 s5
 then
 b
 =
 b
 +
 1
 ;
 bands
[
b
]
 =
 band
(
s5
,
#
s5
,
i
,
j
)
 end
1744                                    local
 s6
 =
 s
[
6
]
 if
 s6
 then
 b
 =
 b
 +
 1
 ;
 bands
[
b
]
 =
 band
(
s6
,
#
s6
,
i
,
j
)
 end
1745                                end
1746                            end
1747                        end
1748                        p
 =
 false
1749                    end
1750                else
1751                    p
 =
 false
1752                end
1753            end
1754            d0
 =
 d1
1755        end
1756        return
 bands
1757    end
1758end
1759
1760function
 mp
.
lmt_contours_edge_set_by_band
(
value
)
1761    local
 p
         =
 getparameterset
(
)
1762    local
 result
    =
 p
.
result
1763
1764    if
 result
.
cached
 then
 return
 end
1765
1766    local
 values
    =
 result
.
values
1767    local
 nofvalues
 =
 result
.
nofvalues
1768    local
 data
      =
 result
.
data
1769    local
 nx
        =
 result
.
nx
1770    local
 ny
        =
 result
.
ny
1771    local
 bands
     =
 {
 }
1772    result
.
bands
    =
 bands
1773    for
 value
=
1
,
nofvalues
 do
1774        bands
[
value
]
 =
 whatever
(
data
,
ny
,
nx
,
value
,
value
 =
=
 1
)
1775    end
1776end
1777
1778function
 mp
.
lmt_contours_edge_get_band
(
value
)
1779    mpfill
(
getparameterset
(
)
.
result
.
bands
[
value
]
,
true
)
1780    mpflush
(
)
1781end
1782
1783-- Because we share some code surface plots also end up here. When working on the
1784-- contour macros by concidence I ran into a 3D plot in
1785--
1786-- https://staff.science.uva.nl/a.j.p.heck/Courses/mptut.pdf
1787--
1788-- The code is pure MetaPost and works quite well. With a bit of optimization
1789-- performance is also ok, but in the end a Lua solution is twice as fast and also
1790-- permits some more tweaking at no cost. So, below is an adaptation of an example
1791-- in the mentioned link. It's one of these cases where access to pseudo arrays
1792-- is slowing down MP.
1793
1794local
 sqrt
,
 sin
,
 cos
 =
 math
.
sqrt
,
 math
.
sin
,
 math
.
cos
1795
1796local
 f_fill_rgb
 =
 formatters
[
"
F (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor (%.3N,%.3N,%.3N) ;
"
]
1797local
 f_draw_rgb
 =
 formatters
[
"
D (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor %.3F ;
"
]
1798local
 f_mesh_rgb
 =
 formatters
[
"
U (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor (%.3N,%.3N,%.3N) ;
"
]
1799local
 f_fill_cmy
 =
 formatters
[
"
F (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor (%.3N,%.3N,%.3N,0) ;
"
]
1800local
 f_draw_cmy
 =
 formatters
[
"
D (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor %.3F ;
"
]
1801local
 f_mesh_cmy
 =
 formatters
[
"
U (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor (%.3N,%.3N,%.3N,0) ;
"
]
1802
1803local
 f_function_n
 =
 formatters
 [
 [[

1804    local math  = math
1805    local round = math.round
1806    %s
1807    return function(x,y)
1808        return %s
1809    end
1810
]]
 ]
1811
1812local
 f_function_y
 =
 formatters
 [
 [[

1813    local math  = math
1814    local round = math.round
1815    local nan   = NaN
1816    local inf   = math.huge
1817    local er    = 0
1818    %s
1819    return function(x,y,dnan,dinf,report)
1820        local n = %s
1821        if n == nan then
1822            er = er + 1
1823            if er < 10 then
1824                report("nan at (%s,%s)",x,y)
1825            end
1826            n = dnan
1827        elseif n == inf then
1828            er = er + 1
1829            if er < 10 then
1830                report("inf at (%s,%s)",x,y)
1831            end
1832            n = dinf
1833        end
1834        dx[my] = n
1835        sy = sy + 1
1836    end
1837    return n, er
1838end
1839
]]
 ]
1840
1841local
 f_color
 =
 formatters
 [
 [[

1842    local math = math
1843    return function(f)
1844        return %s
1845    end
1846
]]
 ]
1847
1848function
 mp
.
lmt_surface_do
(
specification
)
1849    --
1850    -- The projection and color brightness calculation have been inlined. We also store
1851    -- differently.
1852    --
1853    -- todo: ignore weird paths
1854    --
1855    -- The prototype is now converted to use lmt parameter sets.
1856    --
1857    local
 p
 =
 getparameterset
(
"
surface
"
)
1858    --
1859    local
 preamble
  =
 p
.
preamble
   or
 "
"
1860    local
 code
      =
 p
.
code
       or
 "
return x + y
"
1861    local
 colorcode
 =
 p
.
color
      or
 "
return f, f, f
"
1862    local
 linecolor
 =
 p
.
linecolor
  or
 1
1863    local
 xmin
      =
 p
.
xmin
       or
 -1
1864    local
 xmax
      =
 p
.
xmax
       or
  1
1865    local
 ymin
      =
 p
.
ymin
       or
 -1
1866    local
 ymax
      =
 p
.
ymax
       or
  1
1867    local
 xstep
     =
 p
.
xstep
      or
 .
1
1868    local
 ystep
     =
 p
.
ystep
      or
 .
1
1869    local
 bf
        =
 p
.
brightness
 or
 100
1870    local
 clip
      =
 p
.
clip
       or
 false
1871    local
 lines
     =
 p
.
lines
1872    local
 ha
        =
 p
.
snap
       or
 0
.
01
1873    local
 hb
        =
 2
 *
 ha
1874    --
1875    if
 lines
 =
=
 nil
 then
 lines
 =
 true
 end
1876    --
1877    if
 xstep
 =
=
 0
 then
 xstep
 =
 (
xmax
 -
 xmin
)
/
100
 end
1878    if
 ystep
 =
=
 0
 then
 ystep
 =
 (
ymax
 -
 ymin
)
/
100
 end
1879
1880    local
 nxmin
 =
 round
(
xmin
/
xstep
)
1881    local
 nxmax
 =
 round
(
xmax
/
xstep
)
1882    local
 nymin
 =
 round
(
ymin
/
ystep
)
1883    local
 nymax
 =
 round
(
ymax
/
ystep
)
1884    local
 nx
    =
 nxmax
 -
 nxmin
 +
 1
1885    local
 ny
    =
 nymax
 -
 nymin
 +
 1
1886    --
1887    local
 xvector
 =
 p
.
xvector
    or
 {
 -0
.
7
,
 -0
.
7
 }
1888    local
 yvector
 =
 p
.
yvector
    or
 {
 1
,
 0
 }
1889    local
 zvector
 =
 p
.
zvector
    or
 {
 0
,
 1
 }
1890    local
 light
   =
 p
.
light
      or
 {
 3
,
 3
,
 10
 }
1891    --
1892    local
 xrx
,
 xry
   =
 xvector
[
1
]
,
 xvector
[
2
]
1893    local
 yrx
,
 yry
   =
 yvector
[
1
]
,
 yvector
[
2
]
1894    local
 zrx
,
 zry
   =
 zvector
[
1
]
,
 zvector
[
2
]
1895    local
 xp
,
 yp
,
 zp
 =
 light
[
1
]
,
 light
[
2
]
,
 light
[
3
]
1896    --
1897    local
 data
 =
 setmetatableindex
(
"
table
"
)
1898    local
 dx
   =
 (
xmax
 -
 xmin
)
 /
 nx
1899    local
 dy
   =
 (
ymax
 -
 ymin
)
 /
 ny
1900    local
 xt
   =
 xmin
1901    --
1902    local
 minf
,
 maxf
1903    --
1904    -- similar as contours but no data loop here
1905    --
1906    local
 fcode
 =
 load
(
(
p
.
check
 and
 f_function_y
 or
 f_function_n
)
(
preamble
,
code
)
)
1907    local
 func
  =
 type
(
fcode
)
 =
=
 "
function
"
 and
 fcode
(
)
1908    if
 type
(
func
)
 ~
=
 "
function
"
 then
1909        return
 false
 -- fatal error
1910    end
1911    --
1912    local
 ccode
 =
 load
(
f_color
(
colorcode
)
)
1913    local
 color
 =
 type
(
ccode
)
 =
=
 "
function
"
 and
 ccode
(
)
1914    if
 type
(
color
)
 ~
=
 "
function
"
 then
1915        return
 false
 -- fatal error
1916    end
1917    --
1918    for
 i
=
0
,
nx
 do
1919        local
 yt
 =
 ymin
1920        for
 j
=
0
,
ny
 do
1921            local
 zt
 =
 func
(
xt
,
yt
)
1922            -- projection from 3D to 2D coordinates
1923            local
 x
 =
 xt
 *
 xrx
 +
 yt
 *
 yrx
 +
 zt
 *
 zrx
1924            local
 y
 =
 xt
 *
 xry
 +
 yt
 *
 yry
 +
 zt
 *
 zry
1925            local
 z
 =
 zt
1926            -- numerical derivatives by central differences
1927            local
 dfx
 =
 (
func
(
xt
+
ha
,
yt
)
 -
 func
(
xt
-
ha
,
yt
)
)
 /
 hb
1928            local
 dfy
 =
 (
func
(
xt
,
yt
+
ha
)
 -
 func
(
xt
,
yt
-
ha
)
)
 /
 hb
1929            -- compute brightness factor at a point
1930            local
 ztp
 =
 zt
 -
 zp
1931            local
 ytp
 =
 yt
 -
 yp
1932            local
 xtp
 =
 xt
 -
 xp
1933            local
 ztp
 =
 zt
 -
 zp
1934            local
 ytp
 =
 yt
 -
 yp
1935            local
 xtp
 =
 xt
 -
 xp
1936            local
 ca
  =
 -
ztp
 +
 dfy
*
ytp
 +
 dfx
*
xtp
1937            local
 cb
  =
 sqrt
(
1
+
dfx
*
dfx
+
dfy
*
dfy
)
1938            local
 cc
  =
 sqrt
(
ztp
*
ztp
 +
 ytp
*
ytp
 +
 xtp
*
xtp
)
1939            local
 fac
 =
 bf
*
ca
/
(
cb
*
cc
*
cc
*
cc
)
1940            -- addition: check range
1941            if
 not
 minf
 then
1942                minf
 =
 fac
1943                maxf
 =
 fac
1944            elseif
 fac
 <
 minf
 then
1945                minf
 =
 fac
1946            elseif
 fac
 >
 maxf
 then
1947                maxf
 =
 fac
1948            end
1949            --
1950            data
[
i
]
[
j
]
 =
 {
 x
,
 y
,
 fac
 }
1951            --
1952            yt
 =
 yt
 +
 dy
1953        end
1954        xt
 =
 xt
 +
 dx
1955    end
1956    local
 result
  =
 {
 }
1957    local
 r
       =
 0
1958    local
 range
   =
 maxf
 -
 minf
1959    local
 cl
      =
 linecolor
 or
 1
1960    local
 enforce
 =
 attributes
.
colors
.
model
 =
=
 "
cmyk
"
1961    for
 i
=
0
,
nx
-1
 do
1962        for
 j
=
0
,
ny
-1
 do
1963            -- points
1964            local
 z1
 =
 data
[
i
]
  [
j
]
1965            local
 z2
 =
 data
[
i
]
  [
j
+
1
]
1966            local
 z3
 =
 data
[
i
+
1
]
[
j
+
1
]
1967            local
 z4
 =
 data
[
i
+
1
]
[
j
]
1968            -- color
1969            local
 cf
 =
 z1
[
3
]
1970            if
 clip
 then
1971                -- best clip here if needed
1972                if
 cf
 <
 0
 then
1973                    cf
 =
 0
1974                elseif
 cf
 >
 1
 then
1975                    cf
 =
 1
1976                end
1977            else
1978                -- or remap when we want to
1979                cf
 =
 (
z1
[
3
]
 -
 minf
)
 /
 range
1980            end
1981            local
 z11
 =
 z1
[
1
]
1982            local
 z12
 =
 z1
[
2
]
1983            local
 z21
 =
 z2
[
1
]
1984            local
 z22
 =
 z2
[
2
]
1985            local
 z31
 =
 z3
[
1
]
1986            local
 z32
 =
 z3
[
2
]
1987            local
 z41
 =
 z4
[
1
]
1988            local
 z42
 =
 z4
[
2
]
1989         -- if lines then
1990         --     -- fill first and draw then, previous shapes can be covered
1991         -- else
1992         --     -- fill and draw in one go to prevent artifacts
1993         -- end
1994            local
 cr
,
 cg
,
 cb
 =
 color
(
cf
)
1995            if
 not
 cr
 then
 cr
 =
 0
 end
1996            if
 not
 cg
 then
 cg
 =
 0
 end
1997            if
 not
 cb
 then
 cb
 =
 0
 end
1998            if
 enforce
 then
1999                cr
,
 cg
,
 cb
 =
 1
 -
 cr
,
 1
 -
 cg
,
 1
 -
 cb
2000                r
 =
 r
 +
 1
2001                if
 lines
 then
2002                    result
[
r
]
 =
 f_fill_cmy
(
z11
,
z12
,
z21
,
z22
,
z31
,
z32
,
z41
,
z42
,
cr
,
cg
,
cb
)
2003                    r
 =
 r
 +
 1
2004                    result
[
r
]
 =
 f_draw_cmy
(
z11
,
z12
,
z21
,
z22
,
z31
,
z32
,
z41
,
z42
,
cl
)
2005                else
2006                    result
[
r
]
 =
 f_mesh_cmy
(
z11
,
z12
,
z21
,
z22
,
z31
,
z32
,
z41
,
z42
,
cr
,
cg
,
cb
)
2007                end
2008            else
2009                r
 =
 r
 +
 1
2010                if
 lines
 then
2011                    result
[
r
]
 =
 f_fill_rgb
(
z11
,
z12
,
z21
,
z22
,
z31
,
z32
,
z41
,
z42
,
cr
,
cg
,
cb
)
2012                    r
 =
 r
 +
 1
2013                    result
[
r
]
 =
 f_draw_rgb
(
z11
,
z12
,
z21
,
z22
,
z31
,
z32
,
z41
,
z42
,
cl
)
2014                else
2015                    result
[
r
]
 =
 f_mesh_rgb
(
z11
,
z12
,
z21
,
z22
,
z31
,
z32
,
z41
,
z42
,
cr
,
cg
,
cb
)
2016                end
2017            end
2018        end
2019    end
2020    mp
.
direct
(
concat
(
result
)
)
2021end
2022
Source Browser ?
