mlib-cnt.lua /size: 72 Kb    last modification: 2020-07-01 14:35
1
if
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 47
local
next
,
type
,
tostring
=
next
,
type
,
tostring
48
local
round
,
abs
,
min
,
max
,
floor
=
math
.
round
,
math
.
abs
,
math
.
min
,
math
.
max
,
math
.
floor
49
local
concat
,
move
=
table
.
concat
,
table
.
move
50 51
local
bor
=
bit32
.
bor
-- it's really time to ditch support for luajit
52 53
local
starttiming
=
statistics
.
starttiming
54
local
stoptiming
=
statistics
.
stoptiming
55
local
elapsedtime
=
statistics
.
elapsedtime
56 57
local
formatters
=
string
.
formatters
58
local
setmetatableindex
=
table
.
setmetatableindex
59
local
sortedkeys
=
table
.
sortedkeys
60
local
sequenced
=
table
.
sequenced
61 62
local
metapost
=
metapost
or
{
}
63
local
mp
=
mp
or
{
}
64 65
local
getparameterset
=
metapost
.
getparameterset
66 67
local
mpflush
=
mp
.
flush
68
local
mpcolor
=
mp
.
color
69
local
mpstring
=
mp
.
string
70
local
mpdraw
=
mp
.
draw
71
local
mpfill
=
mp
.
fill
72
local
mpflatten
=
mp
.
flatten
73 74
local
report
=
logs
.
reporter
(
"
mfun contour
"
)
75 76
local
version
=
0
.
007
77 78
-- we iterate using integers so that we get a better behaviour at zero
79 80
local
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 101
local
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 139
local
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 167
local
inbetween
=
attributes
and
attributes
.
colors
.
helpers
.
inbetween
168 169
mp
.
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 189
local
f_box
=
formatters
[
[[
draw rawtexbox("contour",%s) xysized (%s,%s) ;
]]
]
190 191
local
n_box
=
0
192 193
-- todo: remove old one, so we need to store old hashes
194 195
local
nofcontours
=
0
196 197
-- We don't want cosmetics like axis and labels to trigger a calculation,
198
-- especially a slow one.
199 200
local
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 206
local
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
215
end
216 217
local
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
230
end
231 232
local
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
248
end
249 250
function
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
)
544
end
545 546
function
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
)
556
end
557 558
function
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
622
end
623 624
function
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
)
632
end
633 634
function
mp
.
lmt_contours_cleanup
(
)
635
nodes
.
boxes
.
reset
(
"
contour
"
)
636
n_box
=
0
637
end
638 639
function
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 785
end
786 787
function
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 1362
end
1363 1364
-- accessors
1365 1366
function
mp
.
lmt_contours_nx
(
i
)
return
getparameterset
(
)
.
result
.
nx
end
1367
function
mp
.
lmt_contours_ny
(
i
)
return
getparameterset
(
)
.
result
.
ny
end
1368 1369
function
mp
.
lmt_contours_nofvalues
(
)
return
getparameterset
(
)
.
result
.
nofvalues
end
1370
function
mp
.
lmt_contours_value
(
i
)
return
getparameterset
(
)
.
result
.
values
[
i
]
end
1371 1372
function
mp
.
lmt_contours_minz
(
i
)
return
getparameterset
(
)
.
result
.
minz
end
1373
function
mp
.
lmt_contours_maxz
(
i
)
return
getparameterset
(
)
.
result
.
maxz
end
1374 1375
function
mp
.
lmt_contours_minmean
(
i
)
return
getparameterset
(
)
.
result
.
minmean
end
1376
function
mp
.
lmt_contours_maxmean
(
i
)
return
getparameterset
(
)
.
result
.
maxmean
end
1377 1378
function
mp
.
lmt_contours_xrange
(
)
local
p
=
getparameterset
(
)
mpstring
(
formatters
[
"
x = [%.3N,%.3N] ;
"
]
(
p
.
xmin
,
p
.
xmax
)
)
end
1379
function
mp
.
lmt_contours_yrange
(
)
local
p
=
getparameterset
(
)
mpstring
(
formatters
[
"
y = [%.3N,%.3N] ;
"
]
(
p
.
ymin
,
p
.
ymax
)
)
end
1380 1381
function
mp
.
lmt_contours_format
(
)
1382
local
p
=
getparameterset
(
)
1383
return
mpstring
(
p
.
result
.
islist
and
"
@i
"
or
p
.
zformat
or
p
.
format
)
1384
end
1385 1386
function
mp
.
lmt_contours_function
(
)
1387
local
p
=
getparameterset
(
)
1388
return
mpstring
(
p
.
result
.
islist
and
concat
(
p
[
"
functions
"
]
,
"
,
"
)
or
p
[
"
function
"
]
)
1389
end
1390 1391
function
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
1401
end
1402 1403
function
mp
.
lmt_contours_edge_paths
(
value
)
1404
mpdraw
(
getparameterset
(
)
.
result
.
edges
[
value
]
,
true
)
1405
mpflush
(
)
1406
end
1407 1408
function
mp
.
lmt_contours_shape_paths
(
value
)
1409
mpdraw
(
getparameterset
(
)
.
result
.
shapes
[
value
]
,
false
)
1410
mpflush
(
)
1411
end
1412 1413
function
mp
.
lmt_contours_shade_paths
(
value
)
1414
mpfill
(
getparameterset
(
)
.
result
.
shades
[
value
]
,
true
)
1415
mpflush
(
)
1416
end
1417 1418
function
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
1424
end
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 1432
local
d
=
1
/
2
1433 1434
local
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 1451
local
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
1504
end
1505 1506
-- todo: just fetch when needed, no need to cache
1507 1508
function
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
1524
end
1525 1526
function
mp
.
lmt_contours_edge_get_cell
(
value
)
1527
mpdraw
(
getparameterset
(
)
.
result
.
lines
[
value
]
)
1528
mpflush
(
)
1529
end
1530 1531
local
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 1613
local
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 1640
local
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
1758
end
1759 1760
function
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
1776
end
1777 1778
function
mp
.
lmt_contours_edge_get_band
(
value
)
1779
mpfill
(
getparameterset
(
)
.
result
.
bands
[
value
]
,
true
)
1780
mpflush
(
)
1781
end
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 1794
local
sqrt
,
sin
,
cos
=
math
.
sqrt
,
math
.
sin
,
math
.
cos
1795 1796
local
f_fill_rgb
=
formatters
[
"
F (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor (%.3N,%.3N,%.3N) ;
"
]
1797
local
f_draw_rgb
=
formatters
[
"
D (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor %.3F ;
"
]
1798
local
f_mesh_rgb
=
formatters
[
"
U (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor (%.3N,%.3N,%.3N) ;
"
]
1799
local
f_fill_cmy
=
formatters
[
"
F (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor (%.3N,%.3N,%.3N,0) ;
"
]
1800
local
f_draw_cmy
=
formatters
[
"
D (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor %.3F ;
"
]
1801
local
f_mesh_cmy
=
formatters
[
"
U (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor (%.3N,%.3N,%.3N,0) ;
"
]
1802 1803
local
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 1812
local
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 1841
local
f_color
=
formatters
[
[[
1842 local math = math 1843 return function(f) 1844 return %s 1845 end 1846
]]
]
1847 1848
function
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
)
)
2021
end
2022