1if not modules then modules = { } end modules [ ' mlib-cnt ' ] = {
2 version = 1 . 001 ,
3 comment = " companion to mlib-ctx.mkiv " ,
4 author = " Hans Hagen, PRAGMA-ADE, Hasselt NL " ,
5 copyright = " PRAGMA ADE / ConTeXt Development Team " ,
6 license = " see context related readme files " ,
7}
8
9
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
47local next , type , tostring = next , type , tostring
48local round , abs , min , max , floor = math . round , math . abs , math . min , math . max , math . floor
49local concat , move = table . concat , table . move
50
51local bor = bit32 . bor
52
53local starttiming = statistics . starttiming
54local stoptiming = statistics . stoptiming
55local elapsedtime = statistics . elapsedtime
56
57local formatters = string . formatters
58local setmetatableindex = table . setmetatableindex
59local sortedkeys = table . sortedkeys
60local sequenced = table . sequenced
61
62local metapost = metapost or { }
63local mp = mp or { }
64
65local getparameterset = metapost . getparameterset
66
67local mpflush = mp . flush
68local mpcolor = mp . color
69local mpstring = mp . string
70local mpdraw = mp . draw
71local mpfill = mp . fill
72local mpflatten = mp . flatten
73
74local report = logs . reporter ( " mfun contour " )
75
76local version = 0 . 007
77
78
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 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
377 end
378 end
379
380
381
382
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
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
430
431 elseif functioncode then
432
433 if not executed ( data , functioncode ) then
434 report ( " error in executing function " )
435 return
436 end
437
438 else
439
440 report ( " no valid function(s) " )
441 return
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
517
518
519 local nofvalues = # values
520
521 local colors = setcolors (
522 p . preamble , levels , minz , maxz , p . color , nofvalues
523 )
524
525 p . result = {
526 version = version ,
527 values = values ,
528 nofvalues = nofvalues ,
529 minz = minz ,
530 maxz = maxz ,
531 minmean = minmean ,
532 maxmean = maxmean ,
533 data = data ,
534 color = color ,
535 nx = nx ,
536 ny = ny ,
537 colors = colors ,
538 name = name ,
539 hash = hash ,
540 islist = functionlist and true or false ,
541 }
542
543 report ( " index %i, nx %i, ny %i, levels %i " , nofcontours , nx , ny , nofvalues )
544end
545
546function mp . lmt_contours_stop ( )
547 local p = getparameterset ( )
548 local e = stoptiming ( " lmt_contours " )
549 setcache ( p )
550 p . result = nil
551 local f = p [ " function " ]
552 local l = p . functions
553 report ( " index %i, %0.3f seconds for: %s " ,
554 nofcontours , e , " [ " . . concat ( l or { f } , " ] [ " ) . . " ] "
555 )
556end
557
558function mp . lmt_contours_bitmap_set ( )
559 local p = getparameterset ( )
560 local result = p . result
561
562 local values = result . values
563 local nofvalues = result . nofvalues
564 local rawdata = result . data
565 local nx = result . nx
566 local ny = result . ny
567 local colors = result . colors
568 local depth = # colors [ 1 ]
569
570
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 { }
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
609
610
611
612 local k = 0
613 for y = ny , 1 , -1 do
614 k = k + 1
615 local d = data [ k ]
616 for x = 1 , nx do
617 d [ x ] = rawdata [ x ] [ y ] - 1
618 end
619 end
620
621 result . bitmap = bitmap
622end
623
624function mp . lmt_contours_bitmap_get ( )
625 local p = getparameterset ( )
626 local result = p . result
627 local bitmap = result . bitmap
628 local box = nodes . hpack ( graphics . bitmaps . flush ( bitmap ) )
629 n_box = n_box + 1
630 nodes . boxes . savenode ( " contour " , tostring ( n_box ) , box )
631 return f_box ( n_box , bitmap . xsize , bitmap . ysize )
632end
633
634function mp . lmt_contours_cleanup ( )
635 nodes . boxes . reset ( " contour " )
636 n_box = 0
637end
638
639function mp . lmt_contours_edge_set ( )
640 local p = getparameterset ( )
641 local result = p . result
642
643 if result . cached then return end
644
645 local values = result . values
646 local nofvalues = result . nofvalues
647 local data = result . data
648 local nx = result . nx
649 local ny = result . ny
650
651 local xmin = p . xmin
652 local xmax = p . xmax
653 local ymin = p . ymin
654 local ymax = p . ymax
655 local xstep = p . xstep
656 local ystep = p . ystep
657
658 local wsp = { }
659 local edges = { }
660
661 for value = 1 , nofvalues do
662
663 local iwsp = 0
664 local di = data [ 1 ]
665 local dc
666 local edge = { }
667 local e = 0
668
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
712 for i = 1 , iwsp do
713 local wspi = wsp [ i ]
714 for isq = 3 , 4 do
715 local nsq = wspi [ isq ]
716 if nsq ~ = 0 then
717 local px = wspi [ 1 ]
718 local py = wspi [ 2 ]
719 local p = { px , py }
720 local pn = 2
721 wspi [ isq ] = 0
722 while true do
723 for ii = 1 , iwsp do
724 local w = wsp [ ii ]
725 local n1 = w [ 3 ]
726 local n2 = w [ 4 ]
727 if n1 = = nsq then
728 local x = w [ 1 ]
729 local y = w [ 2 ]
730 if x ~ = px or y ~ = py then
731 pn = pn + 1
732 p [ pn ] = x
733 pn = pn + 1
734 p [ pn ] = y
735 px = x
736 py = y
737 end
738 nsq = n2
739 w [ 3 ] = 0
740 w [ 4 ] = 0
741 if nsq = = 0 then
742 if pn = = 1 then
743 pn = pn + 1
744 p [ pn ] = w
745 end
746 goto flush
747 end
748 elseif n2 = = nsq then
749 local x = w [ 1 ]
750 local y = w [ 2 ]
751 if x ~ = px or y ~ = py then
752 pn = pn + 1
753 p [ pn ] = x
754 pn = pn + 1
755 p [ pn ] = y
756 px = x
757 py = y
758 end
759 nsq = n1
760 w [ 3 ] = 0
761 w [ 4 ] = 0
762 if nsq = = 0 then
763 goto flush
764 end
765 end
766 end
767 end
768 :: flush ::
769 e = e + 1
770 edge [ e ] = p
771 if mpflatten then
772 mpflatten ( p )
773 end
774 end
775 end
776 end
777
778
779 edges [ value ] = edge
780
781 end
782
783 result . edges = edges
784
785end
786
787function mp . lmt_contours_shade_set ( edgetoo )
788 local p = getparameterset ( )
789 local result = p . result
790
791 if result . cached then return end
792
793 local values = result . values
794 local nofvalues = result . nofvalues
795 local data = result . data
796 local nx = result . nx
797 local ny = result . ny
798 local color = result . color
799
800 local edges = setmetatableindex ( " table " )
801 local shades = setmetatableindex ( " table " )
802
803 local sqtype = setmetatableindex ( " table " )
804
805 local xspoly = { 0 , 0 , 0 , 0 , 0 , 0 }
806 local yspoly = { 0 , 0 , 0 , 0 , 0 , 0 }
807 local xrpoly = { }
808 local yrpoly = { }
809
810 local xrpoly = { }
811 local yrpoly = { }
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832 shades [ 1 ] = { { 0 , 0 , nx - 1 , 0 , nx - 1 , ny - 1 , 0 , ny - 1 } }
833 edges [ 1 ] = { { } }
834
835
836
837 for value = 1 , nofvalues do
838
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
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
883 nrp = nrp + 1
884
885 local first = true
886 local nrpoly = 0
887 local nspoly = 0
888 local nrpm = - nrp
889
890 while true do
891
892 if first then
893 first = false
894 if sqtype [ i ] [ j ] = = 0 then
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
913
914
915
916 local i_l , i_c , i_r
917 local j_b , j_c , j_t
918
919 local i_n = i + 1
920 local j_n = j + 1
921
922 local i_p = i - 1
923 local j_p = j - 1
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
958
959 sqtype [ i ] [ j ] = nrp
960
961 break
962
963
964
965 :: foundit21 ::
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
972 xspoly [ 3 ] = i_c ; yspoly [ 3 ] = j_c
973 if d_r [ j_t ] > value then
974 xspoly [ 4 ] = i_l ; yspoly [ 4 ] = j_c
975 if d_c [ j_t ] > value then
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
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
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
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
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 ::
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
1013 xspoly [ 3 ] = i_c ; yspoly [ 3 ] = j_b
1014 if d_r [ j_c ] > value then
1015 xspoly [ 4 ] = i_c ; yspoly [ 4 ] = j_c
1016 if d_r [ j_t ] > value then
1017 nspoly = 4
1018 else
1019 xspoly [ 5 ] = i_c ; yspoly [ 5 ] = j_c ; nspoly = 5
1020 end
1021 elseif d_r [ j_t ] > value then
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
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
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
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 ::
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
1054 xspoly [ 3 ] = i_l ; yspoly [ 3 ] = j_c
1055 if d_c [ j_t ] > value then
1056 xspoly [ 4 ] = i_l ; yspoly [ 4 ] = j_b
1057 if d_c [ j_c ] > value then
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
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
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
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
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 ::
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
1095 if d_c [ j_c ] > value then
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
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
1106
1107 local xv34 = ( dd3 * i_c - dd4 * i_l ) / ( dd3 - dd4 )
1108 print ( " 4.4 : xv34 " , xv34 , i_c , i_l )
1109
1110
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
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
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
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
1136
1137 :: done ::
1138
1139
1140 if nrpoly = = 0 then
1141
1142 for i = 1 , nspoly do
1143 xrpoly [ i ] = xspoly [ i ]
1144 yrpoly [ i ] = yspoly [ i ]
1145 end
1146 nrpoly = nspoly
1147 else
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
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
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
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
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
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
1278 end
1279 end
1280 :: nomatch2 ::
1281
1282 local dnrpoly = nspoly - 2 * ( match1 + match2 + 1 )
1283 local ispolystart = ( ispoly + match1 ) % nspoly + 1
1284 local irpolyend = ( rpoly1 - match1 - 1 ) % nrpoly + 1
1285 if dnrpoly ~ = 0 then
1286 local irpolystart = ( irpoly + match2 ) % nrpoly + 1
1287 if irpolystart > irpolyend then
1288
1289 if dnrpoly > 0 then
1290
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
1297
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
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
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 )
1345 end
1346 nofs = nofs + 1
1347 shade [ nofs ] = t
1348
1349 end
1350
1351 end
1352
1353 edges [ value + 1 ] = edge
1354 shades [ value + 1 ] = shade
1355
1356
1357 end
1358
1359 result . shades = shades
1360 result . shapes = edges
1361
1362end
1363
1364
1365
1366function mp . lmt_contours_nx ( i ) return getparameterset ( ) . result . nx end
1367function mp . lmt_contours_ny ( i ) return getparameterset ( ) . result . ny end
1368
1369function mp . lmt_contours_nofvalues ( ) return getparameterset ( ) . result . nofvalues end
1370function mp . lmt_contours_value ( i ) return getparameterset ( ) . result . values [ i ] end
1371
1372function mp . lmt_contours_minz ( i ) return getparameterset ( ) . result . minz end
1373function mp . lmt_contours_maxz ( i ) return getparameterset ( ) . result . maxz end
1374
1375function mp . lmt_contours_minmean ( i ) return getparameterset ( ) . result . minmean end
1376function mp . lmt_contours_maxmean ( i ) return getparameterset ( ) . result . maxmean end
1377
1378function mp . lmt_contours_xrange ( ) local p = getparameterset ( ) mpstring ( formatters [ " x = [%.3N,%.3N] ; " ] ( p . xmin , p . xmax ) ) end
1379function mp . lmt_contours_yrange ( ) local p = getparameterset ( ) mpstring ( formatters [ " y = [%.3N,%.3N] ; " ] ( p . ymin , p . ymax ) ) end
1380
1381function mp . lmt_contours_format ( )
1382 local p = getparameterset ( )
1383 return mpstring ( p . result . islist and " @i " or p . zformat or p . format )
1384end
1385
1386function mp . lmt_contours_function ( )
1387 local p = getparameterset ( )
1388 return mpstring ( p . result . islist and concat ( p [ " functions " ] , " , " ) or p [ " function " ] )
1389end
1390
1391function mp . lmt_contours_range ( )
1392 local p = getparameterset ( )
1393 local r = p . result . islist and p . range
1394 if not r or # r = = 0 then
1395 return mpstring ( " " )
1396 elseif # r = = 1 then
1397 return mpstring ( r [ 1 ] )
1398 else
1399 return mpstring ( formatters [ " z = [%s,%s] " ] ( r [ 1 ] , r [ 2 ] ) )
1400 end
1401end
1402
1403function mp . lmt_contours_edge_paths ( value )
1404 mpdraw ( getparameterset ( ) . result . edges [ value ] , true )
1405 mpflush ( )
1406end
1407
1408function mp . lmt_contours_shape_paths ( value )
1409 mpdraw ( getparameterset ( ) . result . shapes [ value ] , false )
1410 mpflush ( )
1411end
1412
1413function mp . lmt_contours_shade_paths ( value )
1414 mpfill ( getparameterset ( ) . result . shades [ value ] , true )
1415 mpflush ( )
1416end
1417
1418function mp . lmt_contours_color ( value )
1419 local p = getparameterset ( )
1420 local color = p . result . colors [ value ]
1421 if color then
1422 mpcolor ( color )
1423 end
1424end
1425
1426
1427
1428
1429
1430
1431
1432local d = 1 / 2
1433
1434local paths = {
1435 { 0 , d , d , 0 } ,
1436 { 1 , d , d , 0 } ,
1437 { 0 , d , 1 , d } ,
1438 { 1 , d , d , 1 } ,
1439 { 0 , d , d , 1 , d , 0 , 1 , d } ,
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 } ,
1445 { 1 , d , d , 1 } ,
1446 { 0 , d , 1 , d } ,
1447 { d , 0 , 1 , d } ,
1448 { d , 0 , 0 , d } ,
1449}
1450
1451local function whatever ( data , nx , ny , threshold )
1452 local edges = { }
1453 local e = 0
1454 local d0 = data [ 1 ]
1455 for j = 1 , ny -1 do
1456 local d1 = data [ j + 1 ]
1457 local k = j + 1
1458 for i = 1 , nx -1 do
1459 local v = 0
1460 local l = i + 1
1461 local c1 = d0 [ i ]
1462 if c1 < threshold then
1463 v = v + 8
1464 end
1465 local c2 = d0 [ l ]
1466 if c2 < threshold then
1467 v = v + 4
1468 end
1469 local c3 = d1 [ l ]
1470 if c3 < threshold then
1471 v = v + 2
1472 end
1473 local c4 = d1 [ i ]
1474 if c4 < threshold then
1475 v = v + 1
1476 end
1477 if v > 0 and v < 15 then
1478 if v = = 5 or v = = 10 then
1479 local a = ( c1 + c2 + c3 + c4 ) / 4
1480 if a < threshold then
1481 v = v = = 5 and 10 or 5
1482 end
1483 local p = paths [ v ]
1484 e = e + 1 edges [ e ] = k - p [ 2 ]
1485 e = e + 1 edges [ e ] = i + p [ 1 ]
1486 e = e + 1 edges [ e ] = k - p [ 4 ]
1487 e = e + 1 edges [ e ] = i + p [ 3 ]
1488 e = e + 1 edges [ e ] = k - p [ 6 ]
1489 e = e + 1 edges [ e ] = i + p [ 5 ]
1490 e = e + 1 edges [ e ] = k - p [ 8 ]
1491 e = e + 1 edges [ e ] = i + p [ 7 ]
1492 else
1493 local p = paths [ v ]
1494 e = e + 1 edges [ e ] = k - p [ 2 ]
1495 e = e + 1 edges [ e ] = i + p [ 1 ]
1496 e = e + 1 edges [ e ] = k - p [ 4 ]
1497 e = e + 1 edges [ e ] = i + p [ 3 ]
1498 end
1499 end
1500 end
1501 d0 = d1
1502 end
1503 return edges
1504end
1505
1506
1507
1508function mp . lmt_contours_edge_set_by_cell ( )
1509 local p = getparameterset ( )
1510 local result = p . result
1511
1512 if result . cached then return end
1513
1514 local values = result . values
1515 local nofvalues = result . nofvalues
1516 local data = result . data
1517 local nx = result . nx
1518 local ny = result . ny
1519 local lines = { }
1520 result . lines = lines
1521 for value = 1 , nofvalues do
1522 lines [ value ] = whatever ( data , ny , nx , value )
1523 end
1524end
1525
1526function mp . lmt_contours_edge_get_cell ( value )
1527 mpdraw ( getparameterset ( ) . result . lines [ value ] )
1528 mpflush ( )
1529end
1530
1531local singles = {
1532 { d , 0 , 0 , 0 , 0 , d } ,
1533 { d , 0 , 0 , d } ,
1534 { 1 , d , 1 , 0 , d , 0 } ,
1535 { 1 , d , 1 , 0 , 0 , 0 , 0 , d } ,
1536 { 1 , d , 1 , 0 , d , 0 , 0 , d } ,
1537 { 1 , d , d , 0 } ,
1538 { 1 , d , d , 0 , 0 , 0 , 0 , d } ,
1539 { 1 , d , 0 , d } ,
1540 { d , 1 , 1 , 1 , 1 , d } ,
1541 false ,
1542 false ,
1543 { d , 1 , 1 , 1 , 1 , 0 , d , 0 } ,
1544 { d , 1 , 1 , 1 , 1 , 0 , 0 , 0 , 0 , d } ,
1545 { d , 1 , 1 , 1 , 1 , 0 , d , 0 , 0 , d } ,
1546 { d , 1 , 1 , 1 , 1 , d , d , 0 } ,
1547 { d , 1 , 1 , 1 , 1 , d , d , 0 , 0 , 0 , 0 , d } ,
1548 { d , 1 , 1 , 1 , 1 , d , 0 , d } ,
1549 { d , 1 , 1 , d } ,
1550 false ,
1551 false ,
1552 { d , 1 , 1 , d , 1 , 0 , d , 0 } ,
1553 { d , 1 , 1 , d , 1 , 0 , 0 , 0 , 0 , d } ,
1554 false ,
1555 { d , 1 , d , 0 } ,
1556 { d , 1 , d , 0 , 0 , 0 , 0 , d } ,
1557 { d , 1 , 0 , d } ,
1558 { 0 , 1 , d , 1 , 0 , d } ,
1559 { 0 , 1 , d , 1 , d , 0 , 0 , 0 } ,
1560 { 0 , 1 , d , 1 , d , 0 , 0 , d } ,
1561 false ,
1562 { 0 , 1 , d , 1 , 1 , d , 1 , 0 , 0 , 0 } ,
1563 { 0 , 1 , d , 1 , 1 , d , 1 , 0 , d , 0 , 0 , d } ,
1564 false ,
1565 { 0 , 1 , d , 1 , 1 , d , d , 0 , 0 , 0 } ,
1566 { 0 , 1 , d , 1 , 1 , d , 0 , d } ,
1567 { 0 , 1 , 1 , 1 , 1 , d , 0 , d } ,
1568 { 0 , 1 , 1 , 1 , 1 , d , d , 0 , 0 , 0 } ,
1569 { 0 , 1 , 1 , 1 , 1 , d , d , 0 , 0 , d } ,
1570 { 0 , 1 , 1 , 1 , 1 , 0 , d , 0 , 0 , d } ,
1571 { 0 , 1 , 1 , 1 , 1 , 0 , 0 , 0 } ,
1572 { 0 , 1 , 1 , 1 , 1 , 0 , d , 0 , 0 , d } ,
1573 { 0 , 1 , 1 , 1 , 1 , d , d , 0 , 0 , d } ,
1574 { 0 , 1 , 1 , 1 , 1 , d , d , 0 , 0 , 0 } ,
1575 { 0 , 1 , 1 , 1 , 1 , d , 0 , d } ,
1576 { 0 , 1 , d , 1 , 1 , d , 0 , d } ,
1577 { 0 , 1 , d , 1 , 1 , d , d , 0 , 0 , 0 } ,
1578 false ,
1579 { 0 , 1 , d , 1 , 1 , d , 1 , 0 , d , 0 , 0 , d } ,
1580 { 0 , 1 , d , 1 , 1 , d , 1 , 0 , 0 , 0 } ,
1581 false ,
1582 { 0 , 1 , d , 1 , d , 0 , 0 , d } ,
1583 { 0 , 1 , d , 1 , d , 0 , 0 , 0 } ,
1584 { 0 , 1 , d , 1 , 0 , d } ,
1585 { d , 1 , 0 , d } ,
1586 { d , 1 , d , 0 , 0 , 0 , 0 , d } ,
1587 { d , 1 , d , 0 } ,
1588 false ,
1589 { d , 1 , 1 , d , 1 , 0 , 0 , 0 , 0 , d } ,
1590 { d , 1 , 1 , d , 1 , 0 , d , 0 } ,
1591 false ,
1592 false ,
1593 { d , 1 , 1 , d } ,
1594 { d , 1 , 1 , 1 , 1 , d , 0 , d } ,
1595 { d , 1 , 1 , 1 , 1 , d , d , 0 , 0 , 0 , 0 , d } ,
1596 { d , 1 , 1 , 1 , 1 , d , d , 0 } ,
1597 { d , 1 , 1 , 1 , 1 , 0 , d , 0 , 0 , d } ,
1598 { d , 1 , 1 , 1 , 1 , 0 , 0 , 0 , 0 , d } ,
1599 { d , 1 , 1 , 1 , 1 , 0 , d , 0 } ,
1600 false ,
1601 false ,
1602 { d , 1 , 1 , 1 , 1 , d } ,
1603 { 1 , d , 0 , d } ,
1604 { 1 , d , d , 0 , 0 , 0 , 0 , d } ,
1605 { 1 , d , d , 0 } ,
1606 { 1 , d , 1 , 0 , d , 0 , 0 , d } ,
1607 { 1 , d , 1 , 0 , 0 , 0 , 0 , d } ,
1608 { 1 , d , 1 , 0 , d , 0 } ,
1609 { d , 0 , 0 , d } ,
1610 { 0 , d , 0 , 0 , d , 0 } ,
1611}
1612
1613local sadles = {
1614 false , false , false , false , false , false , false , false , false ,
1615 { { d , 1 , 1 , 1 , 1 , d } , { d , 0 , 0 , 0 , 0 , d } , { d , 1 , 1 , 1 , 1 , d , d , 0 , 0 , 0 , 0 , d } , false , false , false } ,
1616 { { d , 1 , 1 , 1 , 1 , d } , { d , 0 , 0 , d } , { d , 1 , 1 , 1 , 1 , d , d , 0 , 0 , d } , false , false , false } ,
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 } ,
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 } } ,
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 } } ,
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 } ,
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 } ,
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 } } ,
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 } } ,
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 } ,
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 } } ,
1634 { false , false , { d , 1 , 1 , d , d , 0 , 0 , 0 , 0 , d } , false , { d , 1 , 1 , d } , { d , 0 , 0 , 0 , 0 , d } } ,
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 } } ,
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 } } ,
1638}
1639
1640local function whatever ( data , nx , ny , threshold , background )
1641
1642 if background then
1643
1644 local llx = 1 / 2
1645 local lly = llx
1646 local urx = ny + llx
1647 local ury = nx + lly
1648
1649 return { { llx , lly , urx , 0 , urx , ury , 0 , ury } }
1650
1651 else
1652
1653 local bands = { }
1654 local b = 0
1655
1656 local function band ( s , n , x , y )
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
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
1720 if p then
1721 p [ 4 ] = l
1722 p [ 6 ] = l
1723 else
1724
1725 p = { j , i , j , l , k , l , k , i }
1726 b = b + 1 ; bands [ b ] = p
1727 end
1728 else
1729 local s = singles [ v ]
1730 if s then
1731 b = b + 1 ; bands [ b ] = band ( s , # s , k , i )
1732 else
1733 local s = sadles [ v ]
1734 if s then
1735 local m = ( c1 + c2 + c3 + c4 ) / 4
1736 if m < threshold then
1737 local s1 = s [ 1 ] if s1 then b = b + 1 ; bands [ b ] = band ( s1 , # s1 , i , j ) end
1738 local s2 = s [ 2 ] if s2 then b = b + 1 ; bands [ b ] = band ( s2 , # s2 , i , j ) end
1739 elseif m = = threshold then
1740 local s3 = s [ 3 ] if s3 then b = b + 1 ; bands [ b ] = band ( s3 , # s3 , i , j ) end
1741 local s4 = s [ 4 ] if s4 then b = b + 1 ; bands [ b ] = band ( s4 , # s4 , i , j ) end
1742 else
1743 local s5 = s [ 5 ] if s5 then b = b + 1 ; bands [ b ] = band ( s5 , # s5 , i , j ) end
1744 local s6 = s [ 6 ] if s6 then b = b + 1 ; bands [ b ] = band ( s6 , # s6 , i , j ) end
1745 end
1746 end
1747 end
1748 p = false
1749 end
1750 else
1751 p = false
1752 end
1753 end
1754 d0 = d1
1755 end
1756 return bands
1757 end
1758end
1759
1760function mp . lmt_contours_edge_set_by_band ( value )
1761 local p = getparameterset ( )
1762 local result = p . result
1763
1764 if result . cached then return end
1765
1766 local values = result . values
1767 local nofvalues = result . nofvalues
1768 local data = result . data
1769 local nx = result . nx
1770 local ny = result . ny
1771 local bands = { }
1772 result . bands = bands
1773 for value = 1 , nofvalues do
1774 bands [ value ] = whatever ( data , ny , nx , value , value = = 1 )
1775 end
1776end
1777
1778function mp . lmt_contours_edge_get_band ( value )
1779 mpfill ( getparameterset ( ) . result . bands [ value ] , true )
1780 mpflush ( )
1781end
1782
1783
1784
1785
1786
1787
1788
1789
1790
1791
1792
1793
1794local sqrt , sin , cos = math . sqrt , math . sin , math . cos
1795
1796local f_fill_rgb = formatters [ " F (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor (%.3N,%.3N,%.3N) ; " ]
1797local f_draw_rgb = formatters [ " D (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor %.3F ; " ]
1798local f_mesh_rgb = formatters [ " U (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor (%.3N,%.3N,%.3N) ; " ]
1799local f_fill_cmy = formatters [ " F (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor (%.3N,%.3N,%.3N,0) ; " ]
1800local f_draw_cmy = formatters [ " D (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor %.3F ; " ]
1801local f_mesh_cmy = formatters [ " U (%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--(%.6N,%.6N)--C withcolor (%.3N,%.3N,%.3N,0) ; " ]
1802
1803local f_function_n = formatters [ [[
1804 local math = math
1805 local round = math.round
1806 %s
1807 return function(x,y)
1808 return %s
1809 end
1810 ]] ]
1811
1812local f_function_y = formatters [ [[
1813 local math = math
1814 local round = math.round
1815 local nan = NaN
1816 local inf = math.huge
1817 local er = 0
1818 %s
1819 return function(x,y,dnan,dinf,report)
1820 local n = %s
1821 if n == nan then
1822 er = er + 1
1823 if er < 10 then
1824 report("nan at (%s,%s)",x,y)
1825 end
1826 n = dnan
1827 elseif n == inf then
1828 er = er + 1
1829 if er < 10 then
1830 report("inf at (%s,%s)",x,y)
1831 end
1832 n = dinf
1833 end
1834 dx[my] = n
1835 sy = sy + 1
1836 end
1837 return n, er
1838end
1839 ]] ]
1840
1841local f_color = formatters [ [[
1842 local math = math
1843 return function(f)
1844 return %s
1845 end
1846 ]] ]
1847
1848function mp . lmt_surface_do ( specification )
1849
1850
1851
1852
1853
1854
1855
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
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
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
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
1923 local x = xt * xrx + yt * yrx + zt * zrx
1924 local y = xt * xry + yt * yry + zt * zry
1925 local z = zt
1926
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
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
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
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
1969 local cf = z1 [ 3 ]
1970 if clip then
1971
1972 if cf < 0 then
1973 cf = 0
1974 elseif cf > 1 then
1975 cf = 1
1976 end
1977 else
1978
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
1990
1991
1992
1993
1994 local cr , cg , cb = color ( cf )
1995 if not cr then cr = 0 end
1996 if not cg then cg = 0 end
1997 if not cb then cb = 0 end
1998 if enforce then
1999 cr , cg , cb = 1 - cr , 1 - cg , 1 - cb
2000 r = r + 1
2001 if lines then
2002 result [ r ] = f_fill_cmy ( z11 , z12 , z21 , z22 , z31 , z32 , z41 , z42 , cr , cg , cb )
2003 r = r + 1
2004 result [ r ] = f_draw_cmy ( z11 , z12 , z21 , z22 , z31 , z32 , z41 , z42 , cl )
2005 else
2006 result [ r ] = f_mesh_cmy ( z11 , z12 , z21 , z22 , z31 , z32 , z41 , z42 , cr , cg , cb )
2007 end
2008 else
2009 r = r + 1
2010 if lines then
2011 result [ r ] = f_fill_rgb ( z11 , z12 , z21 , z22 , z31 , z32 , z41 , z42 , cr , cg , cb )
2012 r = r + 1
2013 result [ r ] = f_draw_rgb ( z11 , z12 , z21 , z22 , z31 , z32 , z41 , z42 , cl )
2014 else
2015 result [ r ] = f_mesh_rgb ( z11 , z12 , z21 , z22 , z31 , z32 , z41 , z42 , cr , cg , cb )
2016 end
2017 end
2018 end
2019 end
2020 mp . direct ( concat ( result ) )
2021end
2022 |