Root/nanonote-files/example-files/data/Examples/lua-plplot-examples/x21.lua

1--[[ $Id: x21.lua 9533 2009-02-16 22:18:37Z smekal $
2    Grid data demo
3
4   Copyright (C) 200 Werner Smekal
5
6   This file is part of PLplot.
7
8   PLplot is free software you can redistribute it and/or modify
9   it under the terms of the GNU General Library Public License as published
10   by the Free Software Foundation either version 2 of the License, or
11   (at your option) any later version.
12
13   PLplot is distributed in the hope that it will be useful,
14   but WITHOUT ANY WARRANTY without even the implied warranty of
15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16   GNU Library General Public License for more details.
17
18   You should have received a copy of the GNU Library General Public License
19   along with PLplot if not, write to the Free Software
20   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21--]]
22
23
24-- initialise Lua bindings for PLplot examples.
25dofile("plplot_examples.lua")
26
27-- bitwise or operator from http://lua-users.org/wiki/BaseSixtyFour
28-- (c) 2006-2008 by Alex Kloss
29-- licensed under the terms of the LGPL2
30
31-- return single bit (for OR)
32function bit(x,b)
33    return (math.mod(x, 2^b) - math.mod(x,2^(b-1)) > 0)
34end
35
36-- logic OR for number values
37function lor(x,y)
38    result = 0
39    for p=1,8 do result = result + (((bit(x,p) or bit(y,p)) == true) and 2^(p-1) or 0) end
40    return result
41end
42
43-- Options data structure definition.
44pts = 500
45xp = 25
46yp = 20
47nl = 16
48knn_order = 20
49threshold = 1.001
50wmin = -1e3
51randn = 0
52rosen = 0
53
54
55function cmap1_init()
56  i = { 0, 1 } -- left and right boundary
57
58  h = { 240, 0 } -- blue -> green -> yellow -> red
59  l = { 0.6, 0.6 }
60  s = { 0.8, 0.8 }
61
62  pl.scmap1n(256)
63  pl.scmap1l(0, i, h, l, s)
64end
65
66
67function create_grid(px, py)
68  local x = {}
69  local y = {}
70
71  for i = 1, px do
72    x[i] = xm + (xM-xm)*(i-1)/(px-1)
73  end
74
75  for i = 1, py do
76    y[i] = ym + (yM-ym)*(i-1)/(py-1)
77  end
78  
79  return x, y
80end
81
82
83function create_data(pts)
84  local x = {}
85  local y = {}
86  local z = {}
87
88  for i = 1, pts do
89    xt = (xM-xm)*pl.randd()
90    yt = (yM-ym)*pl.randd()
91    if randn==0 then
92      x[i] = xt + xm
93      y[i] = yt + ym
94    else -- std=1, meaning that many points are outside the plot range
95      x[i] = math.sqrt(-2*math.log(xt)) * math.cos(2*math.pi*yt) + xm
96      y[i] = math.sqrt(-2*math.log(xt)) * math.sin(2*math.pi*yt) + ym
97    end
98    if rosen==0 then
99      r = math.sqrt(x[i]^2 + y[i]^2)
100      z[i] = math.exp(-r^2) * math.cos(2*math.pi*r)
101    else
102      z[i] = math.log((1-x[i])^2 + 100*(y[i] - x[i]^2)^2)
103    end
104  end
105  
106  return x, y, z
107end
108
109
110title = { "Cubic Spline Approximation",
111          "Delaunay Linear Interpolation",
112          "Natural Neighbors Interpolation",
113          "KNN Inv. Distance Weighted",
114          "3NN Linear Interpolation",
115          "4NN Around Inv. Dist. Weighted" }
116
117
118
119xm = -0.2
120ym = -0.2
121xM = 0.6
122yM = 0.6
123
124pl.parseopts(arg, pl.PL_PARSE_FULL)
125
126opt = { 0, 0, wmin, knn_order, threshold, 0 }
127
128-- Initialize plplot
129pl.init()
130
131-- Initialise random number generator
132pl.seed(5489)
133
134x, y, z = create_data(pts) -- the sampled data
135zmin = z[1]
136zmax = z[1]
137for i=2, pts do
138  if z[i]>zmax then zmax = z[i] end
139  if z[i]<zmin then zmin = z[i] end
140end
141
142xg, yg = create_grid(xp, yp) -- grid the data at
143clev = {}
144
145pl.col0(1)
146pl.env(xm, xM, ym, yM, 2, 0)
147pl.col0(15)
148pl.lab("X", "Y", "The original data sampling")
149pl.col0(2)
150pl.poin(x, y, 5)
151pl.adv(0)
152
153pl.ssub(3, 2)
154
155for k = 1, 2 do
156  pl.adv(0)
157  for alg=1, 6 do
158    zg = pl.griddata(x, y, z, xg, yg, alg, opt[alg])
159
160    --[[
161       - CSA can generate NaNs (only interpolates?!).
162       - DTLI and NNI can generate NaNs for points outside the convex hull
163         of the data points.
164       - NNLI can generate NaNs if a sufficiently thick triangle is not found
165     
166       PLplot should be NaN/Inf aware, but changing it now is quite a job...
167       so, instead of not plotting the NaN regions, a weighted average over
168       the neighbors is done. --]]
169     
170
171    if alg==pl.GRID_CSA or alg==pl.GRID_DTLI or alg==pl.GRID_NNLI or alg==pl.GRID_NNI then
172      for i = 1, xp do
173        for j = 1, yp do
174          if zg[i][j]~=zg[i][j] then -- average (IDW) over the 8 neighbors
175            zg[i][j] = 0
176            dist = 0
177
178            for ii=i-1, i+1 do
179              if ii<=xp then
180                for jj=j-1, j+1 do
181                  if jj<=yp then
182                    if ii>=1 and jj>=1 and zg[ii][jj]==zg[ii][jj] then
183                      if (math.abs(ii-i) + math.abs(jj-j)) == 1 then
184                        d = 1
185                      else
186                        d = 1.4142
187                      end
188                      zg[i][j] = zg[i][j] + zg[ii][jj]/(d^2)
189                      dist = dist + d
190                    end
191                  end
192                end
193              end
194            end
195            if dist~=0 then
196              zg[i][j] = zg[i][j]/dist
197            else
198              zg[i][j] = zmin
199            end
200          end
201        end
202      end
203    end
204
205    lzM, lzm = pl.MinMax2dGrid(zg)
206
207    if lzm~=lzm then lzm=zmin else lzm = math.min(lzm, zmin) end
208    if lzM~=lzM then lzM=zmax else lzM = math.max(lzM, zmax) end
209    
210    -- Increase limits slightly to prevent spurious contours
211    -- due to rounding errors
212    lzm = lzm-0.01
213    lzM = lzM+0.01
214
215    pl.col0(1)
216
217    pl.adv(alg)
218
219    if k==1 then
220      for i = 1, nl do
221        clev[i] = lzm + (lzM-lzm)/(nl-1)*(i-1)
222      end
223
224      pl.env0(xm, xM, ym, yM, 2, 0)
225      pl.col0(15)
226      pl.lab("X", "Y", title[alg])
227      pl.shades(zg, xm, xM, ym, yM, clev, 1, 0, 1, 1)
228      pl.col0(2)
229    else
230      for i = 1, nl do
231        clev[i] = lzm + (lzM-lzm)/(nl-1)*(i-1)
232      end
233
234      cmap1_init()
235      pl.vpor(0, 1, 0, 0.9)
236      pl.wind(-1.1, 0.75, -0.65, 1.20)
237      
238      -- For the comparison to be fair, all plots should have the
239      -- same z values, but to get the max/min of the data generated
240      -- by all algorithms would imply two passes. Keep it simple.
241      --
242      -- pl.w3d(1, 1, 1, xm, xM, ym, yM, zmin, zmax, 30, -60)
243       
244
245      pl.w3d(1, 1, 1, xm, xM, ym, yM, lzm, lzM, 30, -40)
246      pl.box3("bntu", "X", 0, 0,
247              "bntu", "Y", 0, 0,
248              "bcdfntu", "Z", 0.5, 0)
249      pl.col0(15)
250      pl.lab("", "", title[alg])
251      pl.plot3dc(xg, yg, zg, lor(lor(pl.DRAW_LINEXY, pl.MAG_COLOR), pl.BASE_CONT), clev)
252    end
253  end
254end
255
256pl.plend()
257

Archive Download this file

Branches:
backfire
master
release_2010-11-17
release_2010-12-14
release_2011-02-23
release_2011-05-28
release_2011-08-27
release_2011-11-13
release_2012-03-18
release_2012-04-09
release_2012-07-11
release_2012-10-11



interactive