1if not modules then modules = { } end modules ['mlib-cnt'] = {
2 version = 1.001,
3 optimize = true,
4 comment = "companion to mlib-ctx.mkiv",
5 author = "Hans Hagen, PRAGMA-ADE, Hasselt NL",
6 copyright = "PRAGMA ADE / ConTeXt Development Team",
7 license = "see context related readme files",
8}
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48local next, type, tostring = next, type, tostring
49local round, abs, min, max, floor = math.round, math.abs, math.min, math.max, math.floor
50local concat, move = table.concat, table.move
51
52local starttiming = statistics.starttiming
53local stoptiming = statistics.stoptiming
54local resettiming = statistics.resettiming
55local elapsedtiming = statistics.elapsed
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
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
194
195local nofcontours = 0
196
197
198
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 resettiming("lmt_contours")
253 starttiming("lmt_contours")
254
255 nofcontours = nofcontours + 1
256
257 local p = getparameterset()
258
259 local xmin = p.xmin
260 local xmax = p.xmax
261 local ymin = p.ymin
262 local ymax = p.ymax
263 local xstep = p.xstep
264 local ystep = p.ystep
265 local levels = p.levels
266 local colors = p.colors
267 local nx = 0
268 local ny = 0
269 local means = setmetatableindex("number")
270 local values = setmetatableindex("number")
271 local data = setmetatableindex("table")
272 local minmean = false
273 local maxmean = false
274
275 local cached, hash, name = getcache(p)
276
277 local function setcolors(preamble,levels,minz,maxz,color,nofvalues)
278 if colors then
279 local function f(k)
280 return #colors[1] == 1 and 0 or { 0, 0, 0 }
281 end
282 setmetatableindex(colors, function(t,k)
283 local v = f(k)
284 t[k] = v
285 return v
286 end)
287 local c = { }
288 local n = 1
289 for i=1,nofvalues do
290 c[i] = colors[n]
291 n = n + 1
292 end
293 return c, f
294 else
295 local tcolor = f_color(levels,minz,maxz,preamble,color)
296 local colors = { }
297 local fcolor = load(tcolor)
298 if type(fcolor) ~= "function" then
299 report("error in color function, case %i: %s",1,color)
300 fcolor = false
301 else
302 fcolor = fcolor()
303 if type(fcolor) ~= "function" then
304 report("error in color function, case %i: %s",2,color)
305 fcolor = false
306 end
307 end
308 if not fcolor then
309 local color_step = mp.lmt_color_functions.step
310 fcolor = function(l)
311 return color_step(l,levels,0.25,0.50,0.75)
312 end
313 end
314 for i=1,nofvalues do
315 colors[i] = { fcolor(i) }
316 end
317 if attributes.colors.model == "cmyk" then
318 for i=1,#colors do
319 local c = colors[i]
320 colors[i] = { 1 - c[1], 1 - c[2], 1 - c[3], 0 }
321 end
322 end
323 return colors, fcolor
324 end
325 end
326
327 if cached then
328 local result = p.result
329 local colors, color = setcolors(
330 p.preamble,
331 p.levels,
332 result.minz,
333 result.maxz,
334 p.color,
335 result.nofvalues
336 )
337 result.color = color
338 result.colors = colors
339 result.cached = true
340 return
341 end
342
343 local functioncode = p["function"]
344 local functionrange = p.range
345 local functionlist = functionrange and p.functions
346 local preamble = p.preamble
347
348 if xstep == 0 then xstep = (xmax - xmin)/100 end
349 if ystep == 0 then ystep = (ymax - ymin)/100 end
350
351 local nxmin = round(xmin/xstep)
352 local nxmax = round(xmax/xstep)
353 local nymin = round(ymin/ystep)
354 local nymax = round(ymax/ystep)
355 local nx = nxmax - nxmin + 1
356 local ny = nymax - nymin + 1
357
358 local function executed(data,code)
359 local fcode = p.check and f_function_y or f_function_n
360 fcode = fcode(preamble,code)
361 fcode = load(fcode)
362 if type(fcode) == "function" then
363 fcode = fcode()
364 end
365 if type(fcode) == "function" then
366 local er = fcode(
367 data, nx, ny,
368 nxmin, nxmax, xstep,
369 nymin, nymax, ystep,
370 p.defaultnan, p.defaultinf, report
371 )
372 if er > 0 then
373 report("%i errors in: %s",er,code)
374 end
375 return true
376 else
377 return false
378 end
379 end
380
381
382
383
384
385 if functionlist then
386
387 local datalist = { }
388 local minr = functionrange[1]
389 local maxr = functionrange[2] or minr
390 local size = #functionlist
391
392 for l=1,size do
393
394 local func = setmetatableindex("table")
395 if not executed(func,functionlist[l]) then
396 report("error in executing function %i from functionlist",l)
397 return
398 end
399
400 local bit = l
401
402 if l == 1 then
403 for i=1,nx do
404 local di = data[i]
405 local fi = func[i]
406 for j=1,ny do
407 local f = fi[j]
408 if f >= minr and f <= maxr then
409 di[j] = bit
410 else
411 di[j] = 0
412 end
413 end
414 end
415 else
416 for i=1,nx do
417 local di = data[i]
418 local fi = func[i]
419 for j=1,ny do
420 local f = fi[j]
421 if f >= minr and f <= maxr then
422 di[j] = di[j] | bit
423 end
424 end
425 end
426 end
427
428 end
429
430
431
432 elseif functioncode then
433
434 if not executed(data,functioncode) then
435 report("error in executing function")
436 return
437 end
438
439 else
440
441 report("no valid function(s)")
442 return
443
444 end
445
446 minz = data[1][1]
447 maxz = minz
448
449 for i=1,nx do
450 local d = data[i]
451 for j=1,ny do
452 local v = d[j]
453 if v < minz then
454 minz = v
455 elseif v > maxz then
456 maxz = v
457 end
458 end
459 end
460
461 if functionlist then
462
463 for i=minz,maxz do
464 values[i] = i
465 end
466
467 levels = maxz
468
469 minmean = minz
470 maxmean = maxz
471
472 else
473
474 local snap = (maxz - minz + 1) / levels
475
476 for i=1,nx do
477 local d = data[i]
478 for j=1,ny do
479 local dj = d[j]
480 local v = round((dj - minz) / snap)
481 values[v] = values[v] + 1
482 means [v] = means [v] + dj
483 d[j] = v
484 end
485 end
486
487 local list = sortedkeys(values)
488 local count = values
489 local remap = { }
490
491 values = { }
492
493 for i=1,#list do
494 local v = list[i]
495 local m = means[v] / count[v]
496 remap [v] = i
497 values[i] = m
498 if not minmean then
499 minmean = m
500 maxmean = m
501 elseif m < minmean then
502 minmean = m
503 elseif m > maxmean then
504 maxmean = m
505 end
506 end
507
508 for i=1,nx do
509 local d = data[i]
510 for j=1,ny do
511 d[j] = remap[d[j]]
512 end
513 end
514
515 end
516
517
518
519
520 local nofvalues = #values
521
522 local colors = setcolors(
523 p.preamble,levels,minz,maxz,p.color,nofvalues
524 )
525
526 p.result = {
527 version = version,
528 values = values,
529 nofvalues = nofvalues,
530 minz = minz,
531 maxz = maxz,
532 minmean = minmean,
533 maxmean = maxmean,
534 data = data,
535 color = color,
536 nx = nx,
537 ny = ny,
538 colors = colors,
539 name = name,
540 hash = hash,
541 islist = functionlist and true or false,
542 }
543
544 report("index %i, nx %i, ny %i, levels %i", nofcontours, nx, ny, nofvalues)
545end
546
547function mp.lmt_contours_stop()
548 local p = getparameterset()
549 stoptiming("lmt_contours")
550 setcache(p)
551 p.result = nil
552 local f = p["function"]
553 local l = p.functions
554 report("index %i, %0.3f seconds for: %s", nofcontours, elapsedtiming("lmt_contours"), "[ " .. concat(l or { f } ," ] [ ") .. " ]")
555end
556
557function mp.lmt_contours_bitmap_set()
558 local p = getparameterset()
559 local result = p.result
560
561 local values = result.values
562 local nofvalues = result.nofvalues
563 local rawdata = result.data
564 local nx = result.nx
565 local ny = result.ny
566 local colors = result.colors
567 local depth = #colors[1]
568
569
570
571 local bitmap = graphics.bitmaps.new(
572 nx, ny,
573 (depth == 3 and "rgb") or (depth == 4 and "cmyk") or "gray",
574 1, false, true
575 )
576
577 local palette = bitmap.index or { }
578 local data = bitmap.data
579 local p = 0
580
581 if depth == 3 or depth == 4 then
582 for i=1,nofvalues do
583 local c = colors[i]
584 local r = round((c[1] or 0) * 255)
585 local g = round((c[2] or 0) * 255)
586 local b = round((c[3] or 0) * 255)
587 local k = depth == 4 and round((c[4] or 0) * 255) or nil
588 palette[p] = {
589 (r > 255 and 255) or (r < 0 and 0) or r,
590 (g > 255 and 255) or (g < 0 and 0) or g,
591 (b > 255 and 255) or (b < 0 and 0) or b,
592 k
593 }
594 p = p + 1
595 end
596 else
597 for i=1,nofvalues do
598 local s = colors[i][1]
599 local s = round((s or 0) * 255)
600 palette[p] = (
601 (s > 255 and 255) or (s < 0 and 0) or s
602 )
603 p = p + 1
604 end
605 end
606
607
608
609
610
611 local k = 0
612 for y=ny,1,-1 do
613 k = k + 1
614 local d = data[k]
615 for x=1,nx do
616 d[x] = rawdata[x][y] - 1
617 end
618 end
619
620 result.bitmap = bitmap
621end
622
623function mp.lmt_contours_bitmap_get()
624 local p = getparameterset()
625 local result = p.result
626 local bitmap = result.bitmap
627 local box = nodes.hpack(graphics.bitmaps.flush(bitmap))
628 n_box = n_box + 1
629 nodes.boxes.savenode("contour",tostring(n_box),box)
630 return f_box(n_box,bitmap.xsize,bitmap.ysize)
631end
632
633function mp.lmt_contours_cleanup()
634 nodes.boxes.reset("contour")
635 n_box = 0
636end
637
638function mp.lmt_contours_edge_set()
639 local p = getparameterset()
640 local result = p.result
641
642 if result.cached then return end
643
644 local values = result.values
645 local nofvalues = result.nofvalues
646 local data = result.data
647 local nx = result.nx
648 local ny = result.ny
649
650 local xmin = p.xmin
651 local xmax = p.xmax
652 local ymin = p.ymin
653 local ymax = p.ymax
654 local xstep = p.xstep
655 local ystep = p.ystep
656
657 local wsp = { }
658 local edges = { }
659
660 for value=1,nofvalues do
661
662 local iwsp = 0
663 local di = data[1]
664 local dc
665 local edge = { }
666 local e = 0
667
668 for i=1,nx do
669 local di1 = data[i+1]
670 local dij = di[1]
671 local d = dij - value
672 local dij1
673 for j=1,ny do
674 if j < ny then
675 dij1 = di[j+1]
676 dc = dij1 - value
677 if (d >= 0 and dc < 0) or (d < 0 and dc >= 0) then
678 iwsp = iwsp + 1
679 local y = (d * (j+1) - dc * j) / (d - dc)
680 if i == 1 then
681 wsp[iwsp] = { i, y, 0, (i + (j-1)*nx) }
682 elseif i == nx then
683 wsp[iwsp] = { i, y, (i - 1 + (j-1)*nx), 0 }
684 else
685 local jx = (i + (j-1)*nx)
686 wsp[iwsp] = { i, y, jx - 1, jx }
687 end
688 end
689 end
690 if i < nx then
691 local dc = di1[j] - value
692 if (d >= 0 and dc < 0) or (d < 0 and dc >= 0) then
693 iwsp = iwsp + 1
694 local x = (d * (i+1) - dc * i) / (d - dc)
695 if j == 1 then
696 wsp[iwsp] = { x, j, 0, (i + (j-1)*nx) }
697 elseif j == ny then
698 wsp[iwsp] = { x, j, (i + (j-2)*nx), 0 }
699 else
700 local jx = i + (j-1)*nx
701 wsp[iwsp] = { x, j, jx - nx, jx }
702 end
703 end
704 end
705 dij = dij1
706 d = dc
707 end
708 di = di1
709 end
710
711 for i=1,iwsp do
712 local wspi = wsp[i]
713 for isq=3,4 do
714 local nsq = wspi[isq]
715 if nsq ~= 0 then
716 local px = wspi[1]
717 local py = wspi[2]
718 local p = { px, py }
719 local pn = 2
720 wspi[isq] = 0
721 while true do
722 for ii=1,iwsp do
723 local w = wsp[ii]
724 local n1 = w[3]
725 local n2 = w[4]
726 if n1 == nsq then
727 local x = w[1]
728 local y = w[2]
729 if x ~= px or y ~= py then
730 pn = pn + 1
731 p[pn] = x
732 pn = pn + 1
733 p[pn] = y
734 px = x
735 py = y
736 end
737 nsq = n2
738 w[3] = 0
739 w[4] = 0
740 if nsq == 0 then
741 if pn == 1 then
742 pn = pn + 1
743 p[pn] = w
744 end
745 goto flush
746 end
747 elseif n2 == nsq then
748 local x = w[1]
749 local y = w[2]
750 if x ~= px or y ~= py then
751 pn = pn + 1
752 p[pn] = x
753 pn = pn + 1
754 p[pn] = y
755 px = x
756 py = y
757 end
758 nsq = n1
759 w[3] = 0
760 w[4] = 0
761 if nsq == 0 then
762 goto flush
763 end
764 end
765 end
766 end
767 ::flush::
768 e = e + 1
769 edge[e] = p
770 if mpflatten then
771 mpflatten(p)
772 end
773 end
774 end
775 end
776
777
778 edges[value] = edge
779
780 end
781
782 result.edges = edges
783
784end
785
786function mp.lmt_contours_shade_set(edgetoo)
787 local p = getparameterset()
788 local result = p.result
789
790 if result.cached then return end
791
792 local values = result.values
793 local nofvalues = result.nofvalues
794 local data = result.data
795 local nx = result.nx
796 local ny = result.ny
797 local color = result.color
798
799 local edges = setmetatableindex("table")
800 local shades = setmetatableindex("table")
801
802 local sqtype = setmetatableindex("table")
803
804 local xspoly = { 0, 0, 0, 0, 0, 0 }
805 local yspoly = { 0, 0, 0, 0, 0, 0 }
806 local xrpoly = { }
807 local yrpoly = { }
808
809 local xrpoly = { }
810 local yrpoly = { }
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831 shades[1] = { { 0, 0, nx - 1, 0, nx - 1, ny - 1, 0, ny - 1 } }
832 edges [1] = { { } }
833
834
835
836 for value=1,nofvalues do
837
838
839 local edge = { }
840 local nofe = 0
841 local shade = { }
842 local nofs = 0
843
844 for i=1,nx-1 do
845 local s = sqtype[i]
846 for j=1,ny-1 do
847 s[j] = 0
848 end
849 end
850
851 local nrp = 0
852
853 local function addedge(a,b,c,d)
854 nofe = nofe + 1 edge[nofe] = a
855 nofe = nofe + 1 edge[nofe] = b
856 nofe = nofe + 1 edge[nofe] = c
857 nofe = nofe + 1 edge[nofe] = d
858 end
859 while true do
860
861 local i
862 local j
863 local d0 = data[1]
864 local d1 = data[2]
865 for ii=1,nx do
866 local s = sqtype[ii]
867 for jj=1,ny do
868 if s[jj] == 0 then
869 if d0[jj] > value then i = ii j = jj goto foundit end
870 if d1[jj] > value then i = ii j = jj goto foundit end
871 local j1 = jj + 1
872 if d1[j1] > value then i = ii j = jj goto foundit end
873 if d0[j1] > value then i = ii j = jj goto foundit end
874 end
875 end
876 d0 = d1
877 d1 = data[ii+1]
878 end
879 break
880 ::foundit::
881
882 nrp = nrp + 1
883
884 local first = true
885 local nrpoly = 0
886 local nspoly = 0
887 local nrpm = -nrp
888
889 while true do
890
891 if first then
892 first = false
893 if sqtype[i][j] == 0 then
894 goto foundit1
895 end
896 end
897 for ii=1,nx do
898 local s = sqtype[ii]
899 for jj=1,ny do
900 if s[jj] == nrpm then
901 i = ii
902 j = jj
903 goto foundit1
904 end
905 end
906 end
907 break
908 ::foundit1::
909 while true do
910
911
912
913
914
915 local i_l, i_c, i_r
916 local j_b, j_c, j_t
917
918 local i_n = i + 1
919 local j_n = j + 1
920
921 local i_p = i - 1
922 local j_p = j - 1
923
924 local d_c = data[i]
925 local d_r = data[i_n]
926
927 local sq
928
929 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
930 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
931 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
932 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
933 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
934 end end
935
936 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
937 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
938 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
939 end end
940
941 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
942 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
943 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
944 end end
945
946 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
947 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
948 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
949 end end
950
951 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
952 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
953 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
954 end end
955
956
957
958 sqtype[i][j] = nrp
959
960 break
961
962
963
964 ::foundit21::
965
966 sq[j_c] = nrpm
967
968 xspoly[1] = i_l ; yspoly[1] = j_b
969 xspoly[2] = i_c ; yspoly[2] = j_b
970 if d_r[j_c] > value then
971 xspoly[3] = i_c ; yspoly[3] = j_c
972 if d_r[j_t] > value then
973 xspoly[4] = i_l ; yspoly[4] = j_c
974 if d_c[j_t] > value then
975 nspoly = 4
976 else
977 xspoly[5] = i_l ; yspoly[5] = j_c ; nspoly = 5
978 end
979 elseif d_c[j_t] > value then
980 xspoly[4] = i_c ; yspoly[4] = j_c ;
981 xspoly[5] = i_l ; yspoly[5] = j_c ; nspoly = 5
982 else
983 xspoly[4] = i_l ; yspoly[4] = j_c ; nspoly = 4
984 if edgetoo then addedge(i_c, j_c, i_l, j_c) end
985 end
986 elseif d_r[j_t] > value then
987 xspoly[3] = i_c ; yspoly[3] = j_b
988 xspoly[4] = i_c ; yspoly[4] = j_c
989 if d_c[j_t] > value then
990 xspoly[5] = i_l ; yspoly[5] = j_c ; nspoly = 5
991 else
992 xspoly[5] = i_l ; yspoly[5] = j_c ;
993 xspoly[6] = i_l ; yspoly[6] = j_c ; nspoly = 6
994 end
995 elseif d_c[j_t] > value then
996 if edgetoo then addedge(i_c, j_b, i_c, j_c) end
997 xspoly[3] = i_c ; yspoly[3] = j_c ;
998 xspoly[4] = i_l ; yspoly[4] = j_c ; nspoly = 4
999 else
1000 if edgetoo then addedge(i_c, j_b, i_l, j_c) end
1001 xspoly[3] = i_l ; yspoly[3] = j_c ; nspoly = 3
1002 end
1003 goto done
1004
1005 ::foundit22::
1006
1007 sq[j_c] = nrpm
1008
1009 xspoly[1] = i_l ; yspoly[1] = j_c
1010 xspoly[2] = i_l ; yspoly[2] = j_b
1011 if d_c[j_c] > value then
1012 xspoly[3] = i_c ; yspoly[3] = j_b
1013 if d_r[j_c] > value then
1014 xspoly[4] = i_c ; yspoly[4] = j_c
1015 if d_r[j_t] > value then
1016 nspoly = 4
1017 else
1018 xspoly[5] = i_c ; yspoly[5] = j_c ; nspoly = 5
1019 end
1020 elseif d_r[j_t] > value then
1021 xspoly[4] = i_c ; yspoly[4] = j_b ;
1022 xspoly[5] = i_c ; yspoly[5] = j_c ; nspoly = 5
1023 else
1024 if edgetoo then addedge(i_c, j_b, i_c, j_c) end
1025 xspoly[4] = i_c ; yspoly[4] = j_c ; nspoly = 4
1026 end
1027 elseif d_r[j_c] > value then
1028 xspoly[3] = i_l ; yspoly[3] = j_b
1029 xspoly[4] = i_c ; yspoly[4] = j_b
1030 xspoly[5] = i_c ; yspoly[5] = j_c
1031 if d_r[j_t] > value then
1032 nspoly = 5
1033 else
1034 xspoly[6] = i_c ; yspoly[6] = j_c ; nspoly = 6
1035 end
1036 elseif d_r[j_t] > value then
1037 if edgetoo then addedge(i_l, j_b, i_c, j_b) end
1038 xspoly[3] = i_c ; yspoly[3] = j_b
1039 xspoly[4] = i_c ; yspoly[4] = j_c ; nspoly = 4
1040 else
1041 if edgetoo then addedge(i_l, j_b, i_c, j_c) end
1042 xspoly[3] = i_c ; yspoly[3] = j_c ; nspoly = 3
1043 end
1044 goto done
1045
1046 ::foundit23::
1047
1048 sq[j_c] = nrpm
1049
1050 xspoly[1] = i_c ; yspoly[1] = j_b
1051 xspoly[2] = i_c ; yspoly[2] = j_c
1052 if d_r[j_t] > value then
1053 xspoly[3] = i_l ; yspoly[3] = j_c
1054 if d_c[j_t] > value then
1055 xspoly[4] = i_l ; yspoly[4] = j_b
1056 if d_c[j_c] > value then
1057 nspoly = 4
1058 else
1059 xspoly[5] = i_l ; yspoly[5] = j_b ; nspoly = 5
1060 end
1061 elseif d_c[j_c] > value then
1062 xspoly[4] = i_l ; yspoly[4] = j_c
1063 xspoly[5] = i_l ; yspoly[5] = j_b ; nspoly = 5
1064 else
1065 if edgetoo then addedge(i_l, j_c, i_l, j_b) end
1066 xspoly[4] = i_l ; yspoly[4] = j_b ; nspoly = 4
1067 end
1068 elseif d_c[j_t] > value then
1069 xspoly[3] = i_c ; yspoly[3] = j_c
1070 xspoly[4] = i_l ; yspoly[4] = j_c
1071 xspoly[5] = i_l ; yspoly[5] = j_b
1072 if d_c[j_c] > value then
1073 nspoly = 5
1074 else
1075 xspoly[6] = i_l ; yspoly[6] = j_b ; nspoly = 6
1076 end
1077 elseif d_c[j_c] > value then
1078 if edgetoo then addedge(i_c, j_c, i_l, j_c) end
1079 xspoly[3] = i_l ; yspoly[3] = j_c ;
1080 xspoly[4] = i_l ; yspoly[4] = j_b ; nspoly = 4
1081 else
1082 if edgetoo then addedge(i_c, j_c, i_l, j_b) end
1083 xspoly[3] = i_l ; yspoly[3] = j_b ; nspoly = 3
1084 end
1085 goto done
1086
1087 ::foundit24::
1088
1089 sq[j_c] = nrpm
1090
1091 xspoly[1] = i_c ; yspoly[1] = j_c
1092 xspoly[2] = i_l ; yspoly[2] = j_c
1093 if d_c[j_t] > value then
1094 if d_c[j_c] > value then
1095 xspoly[3] = i_l ; yspoly[3] = j_b
1096 xspoly[4] = i_c ; yspoly[4] = j_b
1097 if d_r[j_c] > value then
1098 nspoly = 4
1099 else
1100 xspoly[5] = i_c ; yspoly[5] = j_b ; nspoly = 5
1101 end
1102 else
1103 xspoly[3] = i_l ; yspoly[3] = j_b
1104 if d_r[j_c] > value then
1105
1106 local xv34 = (dd3*i_c-dd4*i_l)/(dd3 - dd4)
1107 print("4.4 : xv34",xv34,i_c,i_l)
1108
1109
1110 xspoly[4] = xv34 ; yspoly[4] = j_b ;
1111 xspoly[5] = i_c ; yspoly[5] = j_b ; nspoly = 5
1112 else
1113 if edgetoo then addedge(i_l, j_b, i_c, j_b) end
1114 xspoly[4] = i_c ; yspoly[4] = j_b ; nspoly = 4
1115 end
1116 end
1117 elseif d_c[j_c] > value then
1118 xspoly[3] = i_l ; yspoly[3] = j_b
1119 xspoly[4] = i_l ; yspoly[4] = j_b
1120 xspoly[5] = i_c ; yspoly[5] = j_b
1121 if d_r[j_c] > value then
1122 nspoly = 5
1123 else
1124 xspoly[6] = i_c ; yspoly[6] = j_b ; nspoly = 6
1125 end
1126 elseif d_r[j_c] > value then
1127 if edgetoo then addedge(i_l, j_c, i_l, j_b) end
1128 xspoly[3] = i_l ; yspoly[3] = j_b
1129 xspoly[4] = i_c ; yspoly[4] = j_b ; nspoly = 4
1130 else
1131 if edgetoo then addedge(i_l, j_c, i_c, j_b) end
1132 xspoly[3] = i_c ; yspoly[3] = j_b ; nspoly = 3
1133 end
1134
1135
1136 ::done::
1137
1138
1139 if nrpoly == 0 then
1140
1141 for i=1,nspoly do
1142 xrpoly[i] = xspoly[i]
1143 yrpoly[i] = yspoly[i]
1144 end
1145 nrpoly = nspoly
1146 else
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198 local minx = xspoly[1]
1199 local miny = yspoly[1]
1200 local maxx = xspoly[1]
1201 local maxy = yspoly[1]
1202 for i=1,nspoly do
1203 local y = yspoly[i]
1204 if y < miny then
1205 miny = y
1206 elseif y > maxy then
1207 maxy = y
1208 end
1209 local x = xspoly[i]
1210 if x < minx then
1211 minx = y
1212 elseif x > maxx then
1213 maxx = x
1214 end
1215 end
1216
1217 local ispoly, irpoly
1218 local xr1 = xrpoly[1]
1219 local yr1 = yrpoly[1]
1220 for r0=nrpoly,1,-1 do
1221 if xr1 >= minx and xr1 <= maxx and yr1 >= miny and yr1 <= maxy then
1222 local xr0 = xrpoly[r0]
1223 local yr0 = yrpoly[r0]
1224 for s0=1,nspoly do
1225 if xr1 == xspoly[s0] and yr1 == yspoly[s0] then
1226 if s0 == nspoly then
1227 if xr0 == xspoly[1] and yr0 == yspoly[1] then
1228 ispoly = s0
1229 irpoly = r0
1230 goto foundit3
1231 end
1232 else
1233 local s1 = s0 + 1
1234 if xr0 == xspoly[s1] and yr0 == yspoly[s1] then
1235 ispoly = s0
1236 irpoly = r0
1237 goto foundit3
1238 end
1239 end
1240 end
1241 end
1242 xr1 = xr0
1243 yr1 = yr0
1244 else
1245 xr1 = xrpoly[r0]
1246 yr1 = yrpoly[r0]
1247 end
1248 end
1249
1250 goto nomatch3
1251 ::foundit3::
1252 local match1 = 0
1253 local rpoly1 = irpoly + nrpoly
1254 local spoly1 = ispoly - 1
1255 for i=2,nspoly-1 do
1256
1257 local ir = (rpoly1 - i) % nrpoly + 1
1258 local is = (spoly1 + i) % nspoly + 1
1259 if xrpoly[ir] == xspoly[is] and yrpoly[ir] == yspoly[is] then
1260 match1 = match1 + 1
1261 else
1262 break
1263 end
1264 end
1265 ::nomatch1::
1266 local match2 = 0
1267 local rpoly2 = irpoly - 1
1268 local spoly2 = ispoly + nspoly
1269 for i=2,nspoly-1 do
1270
1271 local ir = (rpoly2 + i) % nrpoly + 1
1272 local is = (spoly2 - i) % nspoly + 1
1273 if xrpoly[ir] == xspoly[is] and yrpoly[ir] == yspoly[is] then
1274 match2 = match2 + 1
1275 else
1276 break
1277 end
1278 end
1279 ::nomatch2::
1280
1281 local dnrpoly = nspoly - 2*(match1 + match2 + 1)
1282 local ispolystart = (ispoly + match1) % nspoly + 1
1283 local irpolyend = (rpoly1 - match1 - 1) % nrpoly + 1
1284 if dnrpoly ~= 0 then
1285 local irpolystart = (irpoly + match2) % nrpoly + 1
1286 if irpolystart > irpolyend then
1287
1288 if dnrpoly > 0 then
1289
1290 for i=nrpoly,irpolystart,-1 do
1291 local k = i + dnrpoly
1292 xrpoly[k] = xrpoly[i]
1293 yrpoly[k] = yrpoly[i]
1294 end
1295 else
1296
1297 for i=irpolystart,nrpoly do
1298 local k = i + dnrpoly
1299 xrpoly[k] = xrpoly[i]
1300 yrpoly[k] = yrpoly[i]
1301 end
1302 end
1303 end
1304 nrpoly = nrpoly + dnrpoly
1305 end
1306 if nrpoly < irpolyend then
1307 for i=irpolyend,nrpoly+1,-1 do
1308
1309 local k = i - nrpoly
1310 xrpoly[k] = xrpoly[i]
1311 yrpoly[k] = yrpoly[i]
1312 end
1313 end
1314 local n = nspoly - 2 - match1 - match2
1315 if n == 1 then
1316 local irpoly1 = irpolyend % nrpoly + 1
1317 local ispoly1 = ispolystart % nspoly + 1
1318 xrpoly[irpoly1] = xspoly[ispoly1]
1319 yrpoly[irpoly1] = yspoly[ispoly1]
1320 elseif n > 0 then
1321
1322 for i=1,n do
1323 local ii = i - 1
1324 local ir = (irpolyend + ii) % nrpoly + 1
1325 local is = (ispolystart + ii) % nspoly + 1
1326 xrpoly[ir] = xspoly[is]
1327 yrpoly[ir] = yspoly[is]
1328 end
1329 end
1330 ::nomatch3::
1331 end
1332 end
1333 end
1334
1335 if nrpoly > 0 then
1336 local t = { }
1337 local n = 0
1338 for i=1,nrpoly do
1339 n = n + 1 t[n] = xrpoly[i]
1340 n = n + 1 t[n] = yrpoly[i]
1341 end
1342 if mpflatten then
1343 mpflatten(t)
1344 end
1345 nofs = nofs + 1
1346 shade[nofs] = t
1347
1348 end
1349
1350 end
1351
1352 edges [value+1] = edge
1353 shades[value+1] = shade
1354
1355
1356 end
1357
1358 result.shades = shades
1359 result.shapes = edges
1360
1361end
1362
1363
1364
1365function mp.lmt_contours_nx (i) return getparameterset().result.nx end
1366function mp.lmt_contours_ny (i) return getparameterset().result.ny end
1367
1368function mp.lmt_contours_nofvalues() return getparameterset().result.nofvalues end
1369function mp.lmt_contours_value (i) return getparameterset().result.values[i] end
1370
1371function mp.lmt_contours_minz (i) return getparameterset().result.minz end
1372function mp.lmt_contours_maxz (i) return getparameterset().result.maxz end
1373
1374function mp.lmt_contours_minmean (i) return getparameterset().result.minmean end
1375function mp.lmt_contours_maxmean (i) return getparameterset().result.maxmean end
1376
1377function mp.lmt_contours_xrange () local p = getparameterset() mpstring(formatters["x = [%.3N,%.3N] ;"](p.xmin,p.xmax)) end
1378function mp.lmt_contours_yrange () local p = getparameterset() mpstring(formatters["y = [%.3N,%.3N] ;"](p.ymin,p.ymax)) end
1379
1380function mp.lmt_contours_format()
1381 local p = getparameterset()
1382 return mpstring(p.result.islist and "@i" or p.zformat or p.format)
1383end
1384
1385function mp.lmt_contours_function()
1386 local p = getparameterset()
1387 return mpstring(p.result.islist and concat(p["functions"], ", ") or p["function"])
1388end
1389
1390function mp.lmt_contours_range()
1391 local p = getparameterset()
1392 local r = p.result.islist and p.range
1393 if not r or #r == 0 then
1394 return mpstring("")
1395 elseif #r == 1 then
1396 return mpstring(r[1])
1397 else
1398 return mpstring(formatters["z = [%s,%s]"](r[1],r[2]))
1399 end
1400end
1401
1402function mp.lmt_contours_edge_paths(value)
1403 mpdraw(getparameterset().result.edges[value],true)
1404 mpflush()
1405end
1406
1407function mp.lmt_contours_shape_paths(value)
1408 mpdraw(getparameterset().result.shapes[value],false)
1409 mpflush()
1410end
1411
1412function mp.lmt_contours_shade_paths(value)
1413 mpfill(getparameterset().result.shades[value],true)
1414 mpflush()
1415end
1416
1417function mp.lmt_contours_color(value)
1418 local p = getparameterset()
1419 local color = p.result.colors[value]
1420 if color then
1421 mpcolor(color)
1422 end
1423end
1424
1425
1426
1427
1428
1429
1430
1431local d = 1/2
1432
1433local paths = {
1434 { 0, d, d, 0 },
1435 { 1, d, d, 0 },
1436 { 0, d, 1, d },
1437 { 1, d, d, 1 },
1438 { 0, d, d, 1, d, 0, 1, d },
1439 { d, 0, d, 1 },
1440 { 0, d, d, 1 },
1441 { 0, d, d, 1 },
1442 { d, 0, d, 1 },
1443 { 0, d, d, 0, 1, d, d, 1 },
1444 { 1, d, d, 1 },
1445 { 0, d, 1, d },
1446 { d, 0, 1, d },
1447 { d, 0, 0, d },
1448}
1449
1450local function whatever(data,nx,ny,threshold)
1451 local edges = { }
1452 local e = 0
1453 local d0 = data[1]
1454 for j=1,ny-1 do
1455 local d1 = data[j+1]
1456 local k = j + 1
1457 for i=1,nx-1 do
1458 local v = 0
1459 local l = i + 1
1460 local c1 = d0[i]
1461 if c1 < threshold then
1462 v = v + 8
1463 end
1464 local c2 = d0[l]
1465 if c2 < threshold then
1466 v = v + 4
1467 end
1468 local c3 = d1[l]
1469 if c3 < threshold then
1470 v = v + 2
1471 end
1472 local c4 = d1[i]
1473 if c4 < threshold then
1474 v = v + 1
1475 end
1476 if v > 0 and v < 15 then
1477 if v == 5 or v == 10 then
1478 local a = (c1 + c2 + c3 + c4) / 4
1479 if a < threshold then
1480 v = v == 5 and 10 or 5
1481 end
1482 local p = paths[v]
1483 e = e + 1 edges[e] = k - p[2]
1484 e = e + 1 edges[e] = i + p[1]
1485 e = e + 1 edges[e] = k - p[4]
1486 e = e + 1 edges[e] = i + p[3]
1487 e = e + 1 edges[e] = k - p[6]
1488 e = e + 1 edges[e] = i + p[5]
1489 e = e + 1 edges[e] = k - p[8]
1490 e = e + 1 edges[e] = i + p[7]
1491 else
1492 local p = paths[v]
1493 e = e + 1 edges[e] = k - p[2]
1494 e = e + 1 edges[e] = i + p[1]
1495 e = e + 1 edges[e] = k - p[4]
1496 e = e + 1 edges[e] = i + p[3]
1497 end
1498 end
1499 end
1500 d0 = d1
1501 end
1502 return edges
1503end
1504
1505
1506
1507function mp.lmt_contours_edge_set_by_cell()
1508 local p = getparameterset()
1509 local result = p.result
1510
1511 if result.cached then return end
1512
1513 local values = result.values
1514 local nofvalues = result.nofvalues
1515 local data = result.data
1516 local nx = result.nx
1517 local ny = result.ny
1518 local lines = { }
1519 result.lines = lines
1520 for value=1,nofvalues do
1521 lines[value] = whatever(data,ny,nx,value)
1522 end
1523end
1524
1525function mp.lmt_contours_edge_get_cell(value)
1526 mpdraw(getparameterset().result.lines[value])
1527 mpflush()
1528end
1529
1530local singles = {
1531 { d, 0, 0, 0, 0, d },
1532 { d, 0, 0, d },
1533 { 1, d, 1, 0, d, 0 },
1534 { 1, d, 1, 0, 0, 0, 0, d },
1535 { 1, d, 1, 0, d, 0, 0, d },
1536 { 1, d, d, 0 },
1537 { 1, d, d, 0, 0, 0, 0, d },
1538 { 1, d, 0, d },
1539 { d, 1, 1, 1, 1, d },
1540 false,
1541 false,
1542 { d, 1, 1, 1, 1, 0, d, 0 },
1543 { d, 1, 1, 1, 1, 0, 0, 0, 0, d },
1544 { d, 1, 1, 1, 1, 0, d, 0, 0, d },
1545 { d, 1, 1, 1, 1, d, d, 0 },
1546 { d, 1, 1, 1, 1, d, d, 0, 0, 0, 0, d },
1547 { d, 1, 1, 1, 1, d, 0, d },
1548 { d, 1, 1, d },
1549 false,
1550 false,
1551 { d, 1, 1, d, 1, 0, d, 0 },
1552 { d, 1, 1, d, 1, 0, 0, 0, 0, d },
1553 false,
1554 { d, 1, d, 0 },
1555 { d, 1, d, 0, 0, 0, 0, d },
1556 { d, 1, 0, d },
1557 { 0, 1, d, 1, 0, d },
1558 { 0, 1, d, 1, d, 0, 0, 0 },
1559 { 0, 1, d, 1, d, 0, 0, d },
1560 false,
1561 { 0, 1, d, 1, 1, d, 1, 0, 0, 0 },
1562 { 0, 1, d, 1, 1, d, 1, 0, d, 0, 0, d },
1563 false,
1564 { 0, 1, d, 1, 1, d, d, 0, 0, 0 },
1565 { 0, 1, d, 1, 1, d, 0, d },
1566 { 0, 1, 1, 1, 1, d, 0, d },
1567 { 0, 1, 1, 1, 1, d, d, 0, 0, 0 },
1568 { 0, 1, 1, 1, 1, d, d, 0, 0, d },
1569 { 0, 1, 1, 1, 1, 0, d, 0, 0, d },
1570 { 0, 1, 1, 1, 1, 0, 0, 0 },
1571 { 0, 1, 1, 1, 1, 0, d, 0, 0, d },
1572 { 0, 1, 1, 1, 1, d, d, 0, 0, d },
1573 { 0, 1, 1, 1, 1, d, d, 0, 0, 0 },
1574 { 0, 1, 1, 1, 1, d, 0, d },
1575 { 0, 1, d, 1, 1, d, 0, d },
1576 { 0, 1, d, 1, 1, d, d, 0, 0, 0 },
1577 false,
1578 { 0, 1, d, 1, 1, d, 1, 0, d, 0, 0, d },
1579 { 0, 1, d, 1, 1, d, 1, 0, 0, 0 },
1580 false,
1581 { 0, 1, d, 1, d, 0, 0, d },
1582 { 0, 1, d, 1, d, 0, 0, 0 },
1583 { 0, 1, d, 1, 0, d },
1584 { d, 1, 0, d },
1585 { d, 1, d, 0, 0, 0, 0, d },
1586 { d, 1, d, 0 },
1587 false,
1588 { d, 1, 1, d, 1, 0, 0, 0, 0, d },
1589 { d, 1, 1, d, 1, 0, d, 0 },
1590 false,
1591 false,
1592 { d, 1, 1, d },
1593 { d, 1, 1, 1, 1, d, 0, d },
1594 { d, 1, 1, 1, 1, d, d, 0, 0, 0, 0, d },
1595 { d, 1, 1, 1, 1, d, d, 0 },
1596 { d, 1, 1, 1, 1, 0, d, 0, 0, d },
1597 { d, 1, 1, 1, 1, 0, 0, 0, 0, d },
1598 { d, 1, 1, 1, 1, 0, d, 0 },
1599 false,
1600 false,
1601 { d, 1, 1, 1, 1, d },
1602 { 1, d, 0, d },
1603 { 1, d, d, 0, 0, 0, 0, d },
1604 { 1, d, d, 0 },
1605 { 1, d, 1, 0, d, 0, 0, d },
1606 { 1, d, 1, 0, 0, 0, 0, d },
1607 { 1, d, 1, 0, d, 0 },
1608 { d, 0, 0, d },
1609 { 0, d, 0, 0, d, 0 },
1610}
1611
1612local sadles = {
1613 false, false, false, false, false, false, false, false, false,
1614 { { 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 },
1615 { { d, 1, 1, 1, 1, d }, { d, 0, 0, d }, { d, 1, 1, 1, 1, d, d, 0, 0, d }, false, false, false },
1616 false, false, false, false, false, false, false,
1617 { { d, 1, 1, d }, { d, 0, 0, 0, 0, d }, { d, 1, 1, d, d, 0, 0, 0, 0, d }, false, false, false },
1618 { { 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 } },
1619 false, false,
1620 { false, false, { d, 1, 1, d, 1, 0, d, 0, 0, d }, false, { d, 1, 0,d, }, { 1, d, 1, 0,d, 0 } },
1621 false, false, false, false, false, false,
1622 { { 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 },
1623 false, false,
1624 { { 1, 0, d, 0, 0, d, }, { 1, d, d, 0 }, { 0, 1, d, 1, 1, d, d, 0, 0, d }, false, false, false },
1625 false, false, false, false, false, false, false, false, false, false, false, false, false,
1626 { false, false, { 0,1, d, 1, 1, d, d, 0, 0, d }, false, { 0,1, d, 1, 0, d }, {1, d, d, 0 } },
1627 false, false,
1628 { 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 } },
1629 false, false, false, false, false, false,
1630 { { d, 1, 0, d }, { 1, d, 1, 0, 0, d }, { d, 1, 1, d, 1, 0, d, 0, 0, d }, false, false, false },
1631 false, false,
1632 { { 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 } },
1633 { false, false, { d, 1, 1, d, d, 0, 0, 0, 0, d }, false, { d, 1, 1, d }, { d, 0, 0, 0, 0, d } },
1634 false, false, false, false, false, false, false,
1635 { false, false, { d, 1, 1, 1, 1, d, d, 0, 0, d }, false, { d, 1, 1, 1, 1, d }, { d, 0,0,d } },
1636 { 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 } },
1637}
1638
1639local function whatever(data,nx,ny,threshold,background)
1640
1641 if background then
1642
1643 local llx = 1/2
1644 local lly = llx
1645 local urx = ny + llx
1646 local ury = nx + lly
1647
1648 return { { llx, lly, urx, 0, urx, ury, 0, ury } }
1649
1650 else
1651
1652 local bands = { }
1653 local b = 0
1654
1655 local function band(s,n,x,y)
1656 if n == 6 then
1657 return {
1658 x - s[ 2], y + s[ 1], x - s[ 4], y + s[ 3], x - s[ 6], y + s[ 5],
1659 }
1660 elseif n == 8 then
1661 return {
1662 x - s[ 2], y + s[ 1], x - s[ 4], y + s[ 3], x - s[ 6], y + s[ 5],
1663 x - s[ 8], y + s[ 7],
1664 }
1665 elseif n == 10 then
1666 return {
1667 x - s[ 2], y + s[ 1], x - s[ 4], y + s[ 3], x - s[ 6], y + s[ 5],
1668 x - s[ 8], y + s[ 7], x - s[10], y + s[ 9],
1669 }
1670 elseif n == 4 then
1671 return {
1672 x - s[ 2], y + s[ 1], x - s[ 4], y + s[ 3],
1673 }
1674 else
1675 return {
1676 x - s[ 2], y + s[ 1], x - s[ 4], y + s[ 3], x - s[ 6], y + s[ 5],
1677 x - s[ 8], y + s[ 7], x - s[10], y + s[ 9], x - s[12], y + s[11],
1678 }
1679 end
1680 end
1681
1682 local pp = { }
1683
1684 local d0 = data[1]
1685 for j=1,ny-1 do
1686 local d1 = data[j+1]
1687 local k = j + 1
1688 local p = false
1689 for i=1,nx-1 do
1690 local v = 0
1691 local l = i + 1
1692 local c1 = d0[i]
1693 if c1 == threshold then
1694 v = v + 27
1695 elseif c1 > threshold then
1696 v = v + 54
1697 end
1698 local c2 = d0[l]
1699 if c2 == threshold then
1700 v = v + 9
1701 elseif c2 > threshold then
1702 v = v + 18
1703 end
1704 local c3 = d1[l]
1705 if c3 == threshold then
1706 v = v + 3
1707 elseif c3 > threshold then
1708 v = v + 6
1709 end
1710 local c4 = d1[i]
1711 if c4 == threshold then
1712 v = v + 1
1713 elseif c4 > threshold then
1714 v = v + 2
1715 end
1716 if v > 0 and v < 80 then
1717 if v == 40 then
1718
1719 if p then
1720 p[4] = l
1721 p[6] = l
1722 else
1723
1724 p = { j, i, j, l, k, l, k, i }
1725 b = b + 1 ; bands[b] = p
1726 end
1727 else
1728 local s = singles[v]
1729 if s then
1730 b = b + 1 ; bands[b] = band(s,#s,k,i)
1731 else
1732 local s = sadles[v]
1733 if s then
1734 local m = (c1 + c2 + c3 + c4) / 4
1735 if m < threshold then
1736 local s1 = s[1] if s1 then b = b + 1 ; bands[b] = band(s1,#s1,i,j) end
1737 local s2 = s[2] if s2 then b = b + 1 ; bands[b] = band(s2,#s2,i,j) end
1738 elseif m == threshold then
1739 local s3 = s[3] if s3 then b = b + 1 ; bands[b] = band(s3,#s3,i,j) end
1740 local s4 = s[4] if s4 then b = b + 1 ; bands[b] = band(s4,#s4,i,j) end
1741 else
1742 local s5 = s[5] if s5 then b = b + 1 ; bands[b] = band(s5,#s5,i,j) end
1743 local s6 = s[6] if s6 then b = b + 1 ; bands[b] = band(s6,#s6,i,j) end
1744 end
1745 end
1746 end
1747 p = false
1748 end
1749 else
1750 p = false
1751 end
1752 end
1753 d0 = d1
1754 end
1755 return bands
1756 end
1757end
1758
1759function mp.lmt_contours_edge_set_by_band(value)
1760 local p = getparameterset()
1761 local result = p.result
1762
1763 if result.cached then return end
1764
1765 local values = result.values
1766 local nofvalues = result.nofvalues
1767 local data = result.data
1768 local nx = result.nx
1769 local ny = result.ny
1770 local bands = { }
1771 result.bands = bands
1772 for value=1,nofvalues do
1773 bands[value] = whatever(data,ny,nx,value,value == 1)
1774 end
1775end
1776
1777function mp.lmt_contours_edge_get_band(value)
1778 mpfill(getparameterset().result.bands[value],true)
1779 mpflush()
1780end
1781
1782
1783
1784
1785
1786
1787
1788
1789
1790
1791
1792
1793local sqrt, sin, cos = math.sqrt, math.sin, math.cos
1794
1795local f_fill_rgb = formatters["F (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor (%.3N,%.3N,%.3N) ;"]
1796local f_draw_rgb = formatters["D (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor %.3F ;"]
1797local f_mesh_rgb = formatters["U (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor (%.3N,%.3N,%.3N) ;"]
1798local f_fill_cmy = formatters["F (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor (%.3N,%.3N,%.3N,0) ;"]
1799local f_draw_cmy = formatters["D (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor %.3F ;"]
1800local f_mesh_cmy = formatters["U (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor (%.3N,%.3N,%.3N,0) ;"]
1801
1802local f_function_n = formatters [ [[
1803 local math = math
1804 local round = math.round
1805 %s
1806 return function(x,y)
1807 return %s
1808 end
1809]] ]
1810
1811local f_function_y = formatters [ [[
1812 local math = math
1813 local round = math.round
1814 local nan = NaN
1815 local inf = math.huge
1816 local er = 0
1817 %s
1818 return function(x,y,dnan,dinf,report)
1819 local n = %s
1820 if n == nan then
1821 er = er + 1
1822 if er < 10 then
1823 report("nan at (%s,%s)",x,y)
1824 end
1825 n = dnan
1826 elseif n == inf then
1827 er = er + 1
1828 if er < 10 then
1829 report("inf at (%s,%s)",x,y)
1830 end
1831 n = dinf
1832 end
1833 dx[my] = n
1834 sy = sy + 1
1835 end
1836 return n, er
1837end
1838]] ]
1839
1840
1841
1842
1843
1844
1845
1846
1847local f_color = formatters [ [[
1848 local math = math
1849 local min = math.min
1850 local max = math.max
1851 local abs = math.abs
1852 local minz = %s
1853 local maxz = %s
1854 --
1855 local color_value = 0
1856 local color_step = mp.lmt_color_functions.step
1857 local color_shade = mp.lmt_color_functions.shade
1858
1859 local function step(...)
1860 return color_step(color_value,n,...)
1861 end
1862 local function shade(...)
1863 return color_shade(color_value,n,...)
1864 end
1865 -- local function lin(l)
1866 -- return l/n
1867 -- end
1868 %s
1869 return function(f,z)
1870 brightness_factor = f
1871 function_value = z
1872 return %s
1873 end
1874 ]] ]
1875
1876function mp.lmt_surface_do(specification)
1877
1878
1879
1880
1881
1882
1883
1884
1885 local p = getparameterset("surface")
1886
1887 local preamble = p.preamble or ""
1888 local code = p.code or "return x + y"
1889 local colorcode = p.color or "return f, f, f"
1890 local linecolor = p.linecolor or 1
1891 local xmin = p.xmin or -1
1892 local xmax = p.xmax or 1
1893 local ymin = p.ymin or -1
1894 local ymax = p.ymax or 1
1895 local xstep = p.xstep or .1
1896 local ystep = p.ystep or .1
1897 local bf = p.brightness or 100
1898 local clip = p.clip or false
1899 local lines = p.lines
1900 local ha = p.snap or 0.01
1901 local hb = 2 * ha
1902
1903 if lines == nil then lines = true end
1904
1905 if xstep == 0 then xstep = (xmax - xmin)/100 end
1906 if ystep == 0 then ystep = (ymax - ymin)/100 end
1907
1908 local nxmin = round(xmin/xstep)
1909 local nxmax = round(xmax/xstep)
1910 local nymin = round(ymin/ystep)
1911 local nymax = round(ymax/ystep)
1912 local nx = nxmax - nxmin + 1
1913 local ny = nymax - nymin + 1
1914
1915 local xvector = p.xvector or { -0.7, -0.7 }
1916 local yvector = p.yvector or { 1, 0 }
1917 local zvector = p.zvector or { 0, 1 }
1918 local light = p.light or { 3, 3, 10 }
1919
1920 local xrx, xry = xvector[1], xvector[2]
1921 local yrx, yry = yvector[1], yvector[2]
1922 local zrx, zry = zvector[1], zvector[2]
1923 local xp, yp, zp = light[1], light[2], light[3]
1924
1925 local data = setmetatableindex("table")
1926 local dx = (xmax - xmin) / nx
1927 local dy = (ymax - ymin) / ny
1928 local xt = xmin
1929
1930 local minf, maxf, minz, maxz
1931
1932
1933
1934 local fcode = load((p.check and f_function_y or f_function_n)(preamble,code))
1935 local func = type(fcode) == "function" and fcode()
1936 if type(func) ~= "function" then
1937 return false
1938 end
1939
1940 for i=0,nx do
1941 local yt = ymin
1942 for j=0,ny do
1943 local zt = func(xt,yt)
1944
1945 local x = xt * xrx + yt * yrx + zt * zrx
1946 local y = xt * xry + yt * yry + zt * zry
1947 local z = zt
1948
1949 local dfx = (func(xt+ha,yt) - func(xt-ha,yt)) / hb
1950 local dfy = (func(xt,yt+ha) - func(xt,yt-ha)) / hb
1951
1952 local ztp = zt - zp
1953 local ytp = yt - yp
1954 local xtp = xt - xp
1955 local ztp = zt - zp
1956 local ytp = yt - yp
1957 local xtp = xt - xp
1958 local ca = -ztp + dfy*ytp + dfx*xtp
1959 local cb = sqrt(1+dfx*dfx+dfy*dfy)
1960 local cc = sqrt(ztp*ztp + ytp*ytp + xtp*xtp)
1961 local fac = bf*ca/(cb*cc*cc*cc)
1962
1963 if not minf then
1964 minf = fac
1965 maxf = fac
1966 elseif fac < minf then
1967 minf = fac
1968 elseif fac > maxf then
1969 maxf = fac
1970 end
1971
1972 if not minz then
1973 minz = z
1974 maxz = z
1975 elseif z < minz then
1976 minz = z
1977 elseif z > maxz then
1978 maxz = z
1979 end
1980 data[i][j] = { x, y, fac, z }
1981
1982 yt = yt + dy
1983 end
1984 xt = xt + dx
1985 end
1986 local result = { }
1987 local r = 0
1988 local range = maxf - minf
1989 local cl = linecolor or 1
1990 local enforce = attributes.colors.model == "cmyk"
1991 local ccode = load(f_color(minz,maxz,preamble,colorcode))
1992 local color = type(ccode) == "function" and ccode()
1993 if type(color) ~= "function" then
1994 return false
1995 end
1996 for i=0,nx-1 do
1997 for j=0,ny-1 do
1998
1999 local z1 = data[i] [j]
2000 local z2 = data[i] [j+1]
2001 local z3 = data[i+1][j+1]
2002 local z4 = data[i+1][j]
2003
2004 local cf = z1[3]
2005 if clip then
2006
2007 if cf < 0 then
2008 cf = 0
2009 elseif cf > 1 then
2010 cf = 1
2011 end
2012 else
2013
2014 cf = (z1[3] - minf) / range
2015 end
2016 local z11 = z1[1]
2017 local z12 = z1[2]
2018 local z21 = z2[1]
2019 local z22 = z2[2]
2020 local z31 = z3[1]
2021 local z32 = z3[2]
2022 local z41 = z4[1]
2023 local z42 = z4[2]
2024
2025
2026
2027
2028
2029 local cr, cg, cb = color(cf,z1[4])
2030 if not cr then cr = 0 end
2031 if not cg then cg = 0 end
2032 if not cb then cb = 0 end
2033 if enforce then
2034 cr, cg, cb = 1 - cr, 1 - cg, 1 - cb
2035 r = r + 1
2036 if lines then
2037 result[r] = f_fill_cmy(z11,z12,z21,z22,z31,z32,z41,z42,cr,cg,cb)
2038 r = r + 1
2039 result[r] = f_draw_cmy(z11,z12,z21,z22,z31,z32,z41,z42,cl)
2040 else
2041 result[r] = f_mesh_cmy(z11,z12,z21,z22,z31,z32,z41,z42,cr,cg,cb)
2042 end
2043 else
2044 r = r + 1
2045 if lines then
2046 result[r] = f_fill_rgb(z11,z12,z21,z22,z31,z32,z41,z42,cr,cg,cb)
2047 r = r + 1
2048 result[r] = f_draw_rgb(z11,z12,z21,z22,z31,z32,z41,z42,cl)
2049 else
2050 result[r] = f_mesh_rgb(z11,z12,z21,z22,z31,z32,z41,z42,cr,cg,cb)
2051 end
2052 end
2053 end
2054 end
2055 mp.direct(concat(result))
2056end
2057 |
