|
7042
|
1 ## Copyright (C) 2007 Kai Habel |
|
|
2 ## Copyright (C) 2003 Shai Ayal |
|
|
3 ## |
|
7164
|
4 ## This file is part of Octave. |
|
|
5 ## |
|
|
6 ## Octave is free software; you can redistribute it and/or modify it |
|
7042
|
7 ## under the terms of the GNU General Public License as published by |
|
7164
|
8 ## the Free Software Foundation; either version 3 of the License, or (at |
|
|
9 ## your option) any later version. |
|
7042
|
10 ## |
|
7164
|
11 ## Octave is distributed in the hope that it will be useful, but |
|
7042
|
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
14 ## General Public License for more details. |
|
|
15 ## |
|
|
16 ## You should have received a copy of the GNU General Public License |
|
7164
|
17 ## along with Octave; see the file COPYING. If not, see |
|
|
18 ## <http://www.gnu.org/licenses/>. |
|
7042
|
19 |
|
|
20 ## -*- texinfo -*- |
|
|
21 ## @deftypefn {Function File} {[@var{c}, @var{h}] =} contourf (@var{x}, @var{y}, @var{z}, @var{lvl}) |
|
|
22 ## @deftypefnx {Function File} {[@var{c}, @var{h}] =} contourf (@var{x}, @var{y}, @var{z}, @var{n}) |
|
|
23 ## @deftypefnx {Function File} {[@var{c}, @var{h}] =} contourf (@var{x}, @var{y}, @var{z}) |
|
|
24 ## @deftypefnx {Function File} {[@var{c}, @var{h}] =} contourf (@var{z}, @var{n}) |
|
|
25 ## @deftypefnx {Function File} {[@var{c}, @var{h}] =} contourf (@var{z}, @var{lvl}) |
|
|
26 ## @deftypefnx {Function File} {[@var{c}, @var{h}] =} contourf (@var{z}) |
|
|
27 ## @deftypefnx {Function File} {[@var{c}, @var{h}] =} contourf (@var{ax}, @dots{}) |
|
|
28 ## @deftypefnx {Function File} {[@var{c}, @var{h}] =} contourf (@dots{}, @var{"property"}, @var{val}) |
|
|
29 ## Compute and plot filled contours of the matrix @var{z}. |
|
|
30 ## Parameters @var{x}, @var{y} and @var{n} or @var{lvl} are optional. |
|
|
31 ## |
|
|
32 ## The return value @var{c} is a 2xn matrix containing the contour lines |
|
|
33 ## as described in the help to the contourc function. |
|
|
34 ## |
|
|
35 ## The return value @var{h} is handle-vector to the patch objects creating |
|
|
36 ## the filled contours. |
|
|
37 ## |
|
|
38 ## If @var{x} and @var{y} are ommited they are taken as the row/column |
|
|
39 ## index of @var{z}. @var{n} is a scalar denoting the number of lines |
|
|
40 ## to compute. Alternatively @var{lvl} is a vector containing the |
|
|
41 ## contour levels. If only one value (e.g. lvl0) is wanted, set |
|
|
42 ## @var{lvl} to [lvl0, lvl0]. If both @var{n} or @var{lvl} are omitted |
|
|
43 ## a default value of 10 contour level is assumed. |
|
|
44 ## |
|
|
45 ## If provided, the filled contours are added to the axes object |
|
|
46 ## @var{ax} instead of the current axis. |
|
|
47 ## |
|
|
48 ## The following example plots filled contours of the @code{peaks} |
|
|
49 ## function. |
|
|
50 ## @example |
|
|
51 ## [x, y, z] = peaks (50); |
|
|
52 ## contourf (x, y, z, -7:9) |
|
|
53 ## @end example |
|
|
54 ## @seealso{contour, contourc, patch} |
|
|
55 ## @end deftypefn |
|
|
56 |
|
|
57 ## Author: Kai Habel <kai.habel@gmx.de> |
|
7327
|
58 ## Author: Shai Ayal <shaiay@users.sourceforge.net> |
|
7042
|
59 |
|
|
60 function varargout = contourf (varargin) |
|
|
61 |
|
7215
|
62 [ax, varargin] = __plt_get_axis_arg__ ("contourf", varargin{:}); |
|
7216
|
63 |
|
7215
|
64 [X, Y, Z, lvl, patch_props] = parse_args (varargin); |
|
7216
|
65 |
|
7042
|
66 [nr, nc] = size (Z); |
|
7216
|
67 |
|
7042
|
68 [minx, maxx] = deal (min (X(:)), max (X(:))); |
|
|
69 [miny, maxy] = deal (min (Y(:)), max (Y(:))); |
|
|
70 |
|
|
71 if (diff (lvl) < 10*eps) |
|
|
72 lvl_eps = 1e-6; |
|
|
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 = len:-1:1 |
|
|
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 if (nargout > 0) |
|
|
184 varargout{2} = h; |
|
|
185 varargout{1} = c; |
|
|
186 endif |
|
|
187 |
|
|
188 endfunction |
|
|
189 |
|
7215
|
190 function [X, Y, Z, lvl, patch_props] = parse_args (arg) |
|
7042
|
191 |
|
|
192 patch_props = {}; |
|
|
193 nolvl = false; |
|
|
194 |
|
|
195 for n = 1:numel (arg) |
|
|
196 if (ischar (arg{n})) |
|
|
197 patch_props = arg(n:end); |
|
|
198 arg(n:end) = []; |
|
|
199 break; |
|
|
200 endif |
|
|
201 endfor |
|
|
202 |
|
|
203 if (mod (numel (patch_props), 2) != 0) |
|
|
204 error ("patch: property value is missing"); |
|
|
205 endif |
|
|
206 |
|
|
207 if (numel (arg) < 3) |
|
|
208 Z = arg{1}; |
|
|
209 [X, Y] = meshgrid (1:columns (Z), 1:rows (Z)); |
|
|
210 endif |
|
|
211 |
|
|
212 if (numel (arg) == 1) |
|
|
213 nolvl = true; |
|
|
214 arg(1) = []; |
|
|
215 elseif (numel (arg) == 2) |
|
|
216 lvl = arg{2}; |
|
|
217 arg(1:2) = []; |
|
|
218 elseif (numel (arg) == 3) |
|
|
219 arg{1:3}; |
|
|
220 [X, Y, Z] = deal (arg{1:3}); |
|
|
221 arg(1:3) = []; |
|
|
222 nolvl = true; |
|
|
223 elseif (numel (arg) == 4) |
|
|
224 [X, Y, Z, lvl] = deal (arg{1:4}); |
|
|
225 arg(1:4) = []; |
|
|
226 endif |
|
|
227 |
|
|
228 if (any (size (X) != size (Y))) |
|
|
229 error ("patch: X and Y must be of same size") |
|
|
230 endif |
|
|
231 |
|
|
232 if (isvector (X) || isvector (Y)) |
|
|
233 [X, Y] = meshgrid (X, Y); |
|
|
234 endif |
|
|
235 |
|
|
236 Z_no_inf = Z(!isinf (Z)); |
|
|
237 [minz, maxz] = deal (min (Z_no_inf(:)), max (Z_no_inf(:))); |
|
|
238 Z(isnan (Z)) = -Inf; |
|
|
239 |
|
|
240 if (nolvl) |
|
|
241 lvl = linspace (minz, maxz, 12); |
|
|
242 endif |
|
|
243 |
|
|
244 if (isscalar (lvl)) |
|
|
245 lvl = linspace (minz, maxz, lvl + 2)(1:end-1); |
|
|
246 else |
|
|
247 idx1 = find(lvl < minz); |
|
|
248 idx2 = find(lvl > maxz); |
|
|
249 lvl(idx1(1:end-1)) = []; |
|
|
250 lvl(idx2) = []; |
|
|
251 if (isempty (lvl)) |
|
|
252 lvl = [minz, minz]; |
|
|
253 endif |
|
|
254 endif |
|
|
255 |
|
|
256 endfunction |
|
7245
|
257 |
|
|
258 %!demo |
|
|
259 %! [x, y, z] = peaks (50); |
|
|
260 %! contourf (x, y, z, -7:9) |
|
7327
|
261 |
|
|
262 %!demo |
|
7331
|
263 %! [theta, r] = meshgrid (linspace (0, 2*pi, 64), linspace(0,1,64)); |
|
|
264 %! [X, Y] = pol2cart (theta, r); |
|
|
265 %! Z = sin(2*theta).*(1-r); |
|
|
266 %! contourf(X, Y, abs(Z), 10) |
