Mercurial > hg > octave-nkf

comparison scripts/plot/contourf.m @ 8289:ac7f334d9652

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add contour group objects and the clabel function
author David Bateman <dbateman@free.fr>
date Thu, 30 Oct 2008 00:13:05 +0100
parents 586b02ac671e
children eb63fbe60fab
comparison
equal deleted inserted replaced
8288:2368aa769ab9 8289:ac7f334d9652
55 ## @end deftypefn 55 ## @end deftypefn
56 56
57 ## Author: Kai Habel <kai.habel@gmx.de> 57 ## Author: Kai Habel <kai.habel@gmx.de>
58 ## Author: Shai Ayal <shaiay@users.sourceforge.net> 58 ## Author: Shai Ayal <shaiay@users.sourceforge.net>
59 59
60 function varargout = contourf (varargin) 60 function [c, h] = contourf (varargin)
61 61
62 [ax, varargin] = __plt_get_axis_arg__ ("contourf", varargin{:}); 62 [xh, varargin] = __plt_get_axis_arg__ ("contour", varargin{:});
63 63
64 [X, Y, Z, lvl, patch_props] = parse_args (varargin); 64 oldh = gca ();
65 65 unwind_protect
66 [nr, nc] = size (Z); 66 axes (xh);
67 67 newplot ();
68 [minx, maxx] = deal (min (X(:)), max (X(:))); 68 [ctmp, htmp] = __contour__ (xh, "none", "fill", "on",
69 [miny, maxy] = deal (min (Y(:)), max (Y(:))); 69 "linecolor", "black", varargin{:});
70 70 unwind_protect_cleanup
71 if (diff (lvl) < 10*eps) 71 axes (oldh);
72 lvl_eps = 1e-6; 72 end_unwind_protect
73 else
74 lvl_eps = min (diff (lvl)) / 1000.0;
75 endif
76
77 X0 = prepad(X, nc+1, 2 * X(1, 1) - X(1, 2), 2);
78 X0 = postpad(X0, nc+2, 2 * X(1, nc) - X(1, nc - 1), 2);
79 X0 = [X0(1, :); X0; X0(1, :)];
80 Y0 = prepad(Y, nr+1, 2 * Y(1, 1) - Y(2, 1), 1);
81 Y0 = postpad(Y0, nr+2, 2 * Y(nr, 1) - Y(nr - 1, 1));
82 Y0 = [Y0(:, 1), Y0, Y0(:, 1)];
83
84 Z0 = -Inf(nr+2, nc+2);
85 Z0(2:nr+1, 2:nc+1) = Z;
86 [c, lev] = contourc (X0, Y0, Z0, lvl);
87 cmap = colormap ();
88
89 levx = linspace (min (lev), max (lev), size (cmap, 1));
90
91 newplot ();
92
93 ## Decode contourc output format.
94 i1 = 1;
95 ncont = 0;
96 while (i1 < columns (c))
97 ncont++;
98 cont_lev(ncont) = c(1, i1);
99 cont_len(ncont) = c(2, i1);
100 cont_idx(ncont) = i1+1;
101
102 ii = i1+1:i1+cont_len(ncont);
103 cur_cont = c(:, ii);
104 c(:, ii);
105 startidx = ii(1);
106 stopidx = ii(end);
107 cont_area(ncont) = polyarea (c(1, ii), c(2, ii));
108 i1 += c(2, i1) + 1;
109 endwhile
110
111 ## Handle for each level the case where we have (a) hole(s) in a patch.
112 ## Those are to be filled with the color of level below or with the
113 ## background colour.
114 for k = 1:numel (lev)
115 lvl_idx = find (abs (cont_lev - lev(k)) < lvl_eps);
116 len = numel (lvl_idx);
117 if (len > 1)
118 ## mark = logical(zeros(size(lvl_idx)));
119 mark = false (size (lvl_idx));
120 a = 1;
121 while (a < len)
122 # take 1st patch
123 b = a + 1;
124 pa_idx = lvl_idx(a);
125 # get pointer to contour start, and contour length
126 curr_ct_idx = cont_idx(pa_idx);
127 curr_ct_len = cont_len(pa_idx);
128 # get contour
129 curr_ct = c(:, curr_ct_idx:curr_ct_idx+curr_ct_len-1);
130 b_vec = (a+1):len;
131 next_ct_pt_vec = c(:, cont_idx(lvl_idx(b_vec)));
132 in = inpolygon (next_ct_pt_vec(1,:), next_ct_pt_vec(2,:),
133 curr_ct(1, :), curr_ct(2, :));
134 mark(b_vec(in)) = !mark(b_vec(in));
135 a++;
136 endwhile
137 if (numel (mark) > 0)
138 ## All marked contours describe a hole in a larger contour of
139 ## the same level and must be filled with colour of level below.
140 ma_idx = lvl_idx(mark);
141 if (k > 1)
142 ## Find color of level below.
143 tmp = find(abs(cont_lev - lev(k - 1)) < lvl_eps);
144 lvl_bel_idx = tmp(1);
145 ## Set color of patches found.
146 cont_lev(ma_idx) = cont_lev(lvl_bel_idx);
147 else
148 ## Set lowest level contour to NaN.
149 cont_lev(ma_idx) = NaN;
150 endif
151 endif
152 endif
153 endfor
154
155 ## The algorithm can create patches with the size of the plotting
156 ## area, we would like to draw only the patch with the highest level.
157 del_idx = [];
158 max_idx = find (cont_area == max (cont_area));
159 if (numel (max_idx) > 1)
160 # delete double entries
161 del_idx = max_idx(1:end-1);
162 cont_area(del_idx) = cont_lev(del_idx) = [];
163 cont_len(del_idx) = cont_idx(del_idx) = [];
164 endif
165
166 ## Now we have everything together and can start plotting the patches
167 ## beginning with largest area.
168 [tmp, svec] = sort (cont_area);
169 len = ncont - numel (del_idx);
170 h = zeros (1, len);
171 for n = 1:len
172 idx = svec(n);
173 ii = cont_idx(idx):cont_idx(idx) + cont_len(idx) - 2;
174 h(n) = patch (c(1, ii), c(2, ii), cont_lev(idx), patch_props{:});
175 endfor
176
177 if (min (lev) == max (lev))
178 set (gca (), "clim", [min(lev)-1, max(lev)+1]);
179 else
180 set (gca(), "clim", [min(lev), max(lev)]);
181 endif
182
183 set (gca (), "layer", "top");
184 73
185 if (nargout > 0) 74 if (nargout > 0)
186 varargout{2} = h; 75 c = ctmp;
187 varargout{1} = c; 76 h = htmp;
188 endif 77 endif
189
190 endfunction
191
192 function [X, Y, Z, lvl, patch_props] = parse_args (arg)
193
194 patch_props = {};
195 nolvl = false;
196
197 for n = 1:numel (arg)
198 if (ischar (arg{n}))
199 patch_props = arg(n:end);
200 arg(n:end) = [];
201 break;
202 endif
203 endfor
204
205 if (mod (numel (patch_props), 2) != 0)
206 error ("patch: property value is missing");
207 endif
208
209 if (numel (arg) < 3)
210 Z = arg{1};
211 [X, Y] = meshgrid (1:columns (Z), 1:rows (Z));
212 endif
213
214 if (numel (arg) == 1)
215 nolvl = true;
216 arg(1) = [];
217 elseif (numel (arg) == 2)
218 lvl = arg{2};
219 arg(1:2) = [];
220 elseif (numel (arg) == 3)
221 arg{1:3};
222 [X, Y, Z] = deal (arg{1:3});
223 arg(1:3) = [];
224 nolvl = true;
225 elseif (numel (arg) == 4)
226 [X, Y, Z, lvl] = deal (arg{1:4});
227 arg(1:4) = [];
228 endif
229
230 if (! ((isvector (X) && isvector (Y)) || size_equal (X, Y)))
231 error ("patch: X and Y must be of same size")
232 endif
233
234 if (isvector (X) || isvector (Y))
235 [X, Y] = meshgrid (X, Y);
236 endif
237
238 Z_no_inf = Z(!isinf (Z));
239 [minz, maxz] = deal (min (Z_no_inf(:)), max (Z_no_inf(:)));
240 Z(isnan (Z)) = -Inf;
241
242 if (nolvl)
243 lvl = linspace (minz, maxz, 12);
244 endif
245
246 if (isscalar (lvl))
247 lvl = linspace (minz, maxz, lvl + 2)(1:end-1);
248 else
249 idx1 = find(lvl < minz);
250 idx2 = find(lvl > maxz);
251 lvl(idx1(1:end-1)) = [];
252 lvl(idx2) = [];
253 if (isempty (lvl))
254 lvl = [minz, minz];
255 endif
256 endif
257
258 endfunction 78 endfunction
259 79
260 %!demo 80 %!demo
261 %! [x, y, z] = peaks (50); 81 %! [x, y, z] = peaks (50);
262 %! contourf (x, y, z, -7:9) 82 %! contourf (x, y, z, -7:9)