1######################################################################## 2## 3## Copyright (C) 2007-2021 The Octave Project Developers 4## 5## See the file COPYRIGHT.md in the top-level directory of this 6## distribution or <https://octave.org/copyright/>. 7## 8## This file is part of Octave. 9## 10## Octave is free software: you can redistribute it and/or modify it 11## under the terms of the GNU General Public License as published by 12## the Free Software Foundation, either version 3 of the License, or 13## (at your option) any later version. 14## 15## Octave is distributed in the hope that it will be useful, but 16## WITHOUT ANY WARRANTY; without even the implied warranty of 17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18## GNU General Public License for more details. 19## 20## You should have received a copy of the GNU General Public License 21## along with Octave; see the file COPYING. If not, see 22## <https://www.gnu.org/licenses/>. 23## 24######################################################################## 25 26## -*- texinfo -*- 27## @deftypefn {} {[@var{c}, @var{hg}] =} __contour__ (@dots{}) 28## Undocumented internal function. 29## @end deftypefn 30 31function [c, hg] = __contour__ (varargin) 32 33 ax = varargin{1}; 34 zlevel = varargin{2}; 35 filled = "off"; 36 37 linespec.color = "auto"; 38 linespec.linestyle = "-"; 39 opts = {}; 40 i = 3; 41 while (i <= length (varargin)) 42 if (ischar (varargin{i}) || iscellstr (varargin{i})) 43 arg = varargin{i}; 44 if (i < length (varargin)) 45 if (strcmpi (arg, "fill")) 46 filled = varargin{i+1}; 47 varargin(i:i+1) = []; 48 continue; 49 elseif (strcmpi (arg, "linecolor")) 50 linespec.color = varargin{i+1}; 51 varargin(i:i+1) = []; 52 continue; 53 endif 54 endif 55 56 [lspec, valid] = __pltopt__ ("__contour__", arg, false); 57 if (valid) 58 varargin(i) = []; 59 if (! isempty (lspec.color)) 60 linespec.color = lspec.color; 61 endif 62 if (! isempty (lspec.linestyle)) 63 linespec.linestyle = lspec.linestyle; 64 endif 65 else # unrecognized option, pass unmodified in opts cell array 66 if (i < length (varargin)) 67 opts(end+(1:2)) = varargin(i:i+1); 68 varargin(i:i+1) = []; 69 else 70 error ("__contour__: Uneven number of PROP/VAL pairs"); 71 endif 72 endif 73 74 else # skip numeric arguments 75 i += 1; 76 endif 77 endwhile 78 79 if (length (varargin) < 5) 80 z1 = varargin{3}; 81 x1 = 1 : columns (z1); 82 y1 = 1 : rows (z1); 83 if (! ismatrix (z1)) 84 error ("__contour__: Z must be a 2-D matrix"); 85 endif 86 else 87 x1 = varargin{3}; 88 y1 = varargin{4}; 89 z1 = varargin{5}; 90 if (! ismatrix (z1) || ! ismatrix (x1) || ! ismatrix (y1)) 91 error ("__contour__: X, Y, and Z must be matrices"); 92 endif 93 if (isvector (x1)) 94 if (columns (z1) != length (x1)) 95 error ("__contour__: size of X must match number of columns in Z"); 96 endif 97 else 98 if (! size_equal (x1, z1)) 99 error ("__contour__: size of X and Z must match"); 100 endif 101 endif 102 if (isvector (y1)) 103 if (rows (z1) != length (y1)) 104 error ("__contour__: size of Y must match number of rows in Z"); 105 endif 106 else 107 if (! size_equal (y1, z1)) 108 error ("__contour__: size of Y and Z must match"); 109 endif 110 endif 111 endif 112 113 if (length (varargin) == 4 || length (varargin) == 6) 114 vn = varargin{end}; 115 vnauto = false; 116 else 117 vn = 10; 118 vnauto = true; 119 endif 120 121 if (isscalar (vn)) 122 ## FIXME: The levels should be determined similarly to {x,y,z}ticks 123 ## so that they aren't set at extremely odd values. 124 lvl = linspace (min (z1(! isinf (z1))), max (z1(! isinf (z1))), vn + 2); 125 ## Strip off max outlier, min must stay for contourf hole algorithm. 126 lvl = lvl(1:end-1); 127 else 128 lvl = sort (vn); 129 endif 130 131 if (strcmpi (filled, "on")) 132 if (isvector (x1) || isvector (y1)) 133 [x1, y1] = meshgrid (x1, y1); 134 endif 135 [nr, nc] = size (z1); 136 x0 = prepad (x1, nc+1, 2 * x1(1, 1) - x1(1, 2), 2); 137 x0 = postpad (x0, nc+2, 2 * x1(1, nc) - x1(1, nc - 1), 2); 138 x0 = [x0(1, :); x0; x0(1, :)]; 139 y0 = prepad (y1, nr+1, 2 * y1(1, 1) - y1(2, 1), 1); 140 y0 = postpad (y0, nr+2, 2 * y1(nr, 1) - y1(nr - 1, 1)); 141 y0 = [y0(:, 1), y0, y0(:, 1)]; 142 z0 = -Inf (nr+2, nc+2); 143 z0(2:nr+1, 2:nc+1) = z1; 144 [c, lev] = contourc (x0, y0, z0, lvl); 145 else 146 [c, lev] = contourc (x1, y1, z1, lvl); 147 endif 148 149 hg = hggroup ("__appdata__", struct ("__creator__", "__contour__")); 150 opts = __add_datasource__ ("__countour__", hg, {"x", "y", "z"}, opts{:}); 151 152 addproperty ("xdata", hg, "data", x1); 153 addproperty ("ydata", hg, "data", y1); 154 addproperty ("zdata", hg, "data", z1); 155 addproperty ("contourmatrix", hg, "data", c); 156 157 addlistener (hg, "xdata", @update_data); 158 addlistener (hg, "ydata", @update_data); 159 addlistener (hg, "zdata", @update_data); 160 addlistener (hg, "contourmatrix", @update_data); 161 162 addproperty ("fill", hg, "radio", "on|{off}", filled); 163 164 ## The properties zlevel and zlevelmode don't exist in matlab, but allow the 165 ## use of contourgroups with the contour3, meshc, and surfc functions. 166 if (isnumeric (zlevel)) 167 addproperty ("zlevelmode", hg, "radio", "{none}|auto|manual", "manual"); 168 addproperty ("zlevel", hg, "data", zlevel); 169 else 170 addproperty ("zlevelmode", hg, "radio", "{none}|auto|manual", zlevel); 171 addproperty ("zlevel", hg, "data", 0.); 172 endif 173 174 lvlstep = sum (diff (lvl)) / (length (lvl) - 1); 175 176 addproperty ("levellist", hg, "data", lev); 177 addproperty ("levelstep", hg, "double", lvlstep); 178 if (vnauto) 179 addproperty ("levellistmode", hg, "radio", "{auto}|manual", "auto"); 180 addproperty ("levelstepmode", hg, "radio", "{auto}|manual", "auto"); 181 elseif (isscalar (vn)) 182 addproperty ("levellistmode", hg, "radio", "{auto}|manual", "auto"); 183 addproperty ("levelstepmode", hg, "radio", "{auto}|manual", "manual"); 184 else 185 addproperty ("levellistmode", hg, "radio", "{auto}|manual", "manual"); 186 addproperty ("levelstepmode", hg, "radio", "{auto}|manual", "auto"); 187 endif 188 189 addproperty ("labelspacing", hg, "double", 144); 190 addproperty ("textlist", hg, "data", lev); 191 addproperty ("textlistmode", hg, "radio", "{auto}|manual", "auto"); 192 addproperty ("textstep", hg, "double", lvlstep); 193 addproperty ("textstepmode", hg, "radio", "{auto}|manual", "auto"); 194 addproperty ("showtext", hg, "radio", "on|{off}", "off"); 195 196 addproperty ("linecolor", hg, "color", linespec.color, "{auto}|none"); 197 addproperty ("linestyle", hg, "linelinestyle", linespec.linestyle); 198 addproperty ("linewidth", hg, "linelinewidth", 0.5); 199 200 ## Matlab property, although Octave does not implement it. 201 addproperty ("hittestarea", hg, "radio", "on|{off}", "off"); 202 203 addlistener (hg, "fill", {@update_data, "fill"}); 204 205 addlistener (hg, "zlevelmode", @update_zlevel); 206 addlistener (hg, "zlevel", @update_zlevel); 207 208 addlistener (hg, "levellist", {@update_data, "levellist"}); 209 addlistener (hg, "levelstep", {@update_data, "levelstep"}); 210 addlistener (hg, "levellistmode", @update_data); 211 addlistener (hg, "levelstepmode", @update_data); 212 213 addlistener (hg, "labelspacing", @update_text); 214 addlistener (hg, "textlist", {@update_text, "textlist"}); 215 addlistener (hg, "textlistmode", @update_text); 216 addlistener (hg, "textstep", {@update_text, "textstep"}); 217 addlistener (hg, "textstepmode", @update_text); 218 addlistener (hg, "showtext", @update_text); 219 220 addlistener (hg, "linecolor", @update_line); 221 addlistener (hg, "linestyle", @update_line); 222 addlistener (hg, "linewidth", @update_line); 223 224 ## Set axis before adding patches so that each new patch does not trigger 225 ## new axis calculation. No need if mode is already "manual". 226 if (all (strcmp (get (gca (), {"xlimmode", "ylimmode"}), "auto"))) 227 axis ([min(x1(:)) max(x1(:)) min(y1(:)) max(y1(:))]); 228 endif 229 230 add_patch_children (hg); 231 232 if (! isempty (opts)) 233 set (hg, opts{:}); 234 endif 235 236endfunction 237 238function add_patch_children (hg) 239 240 c = get (hg, "contourmatrix"); 241 lev = get (hg, "levellist"); 242 fill = get (hg, "fill"); 243 zlev = get (hg, "zlevel"); 244 zmode = get (hg, "zlevelmode"); 245 lc = get (hg, "linecolor"); 246 lw = get (hg, "linewidth"); 247 ls = get (hg, "linestyle"); 248 filled = get (hg, "fill"); 249 ca = gca (); 250 251 ## Turn off automatic updating of clim while adding patches 252 climmode = get (ca, "climmode"); 253 set (ca, "climmode", "manual"); 254 255 if (strcmp (lc, "auto")) 256 lc = "flat"; 257 endif 258 259 if (strcmp (filled, "on")) 260 261 lvl_eps = get_lvl_eps (lev); 262 263 ## Decode contourc output format. 264 i = 1; 265 ncont = 0; 266 while (i < columns (c)) 267 ncont += 1; 268 cont_lev(ncont) = c(1, i); 269 cont_len(ncont) = c(2, i); 270 cont_idx(ncont) = i+1; 271 ii = i + (1:cont_len(ncont)); 272 cont_area(ncont) = polyarea (c(1, ii), c(2, ii)); 273 i += cont_len(ncont) + 1; 274 endwhile 275 276 ## Handle for each level the case where we have (a) hole(s) in a patch. 277 ## Those are to be filled with the color of level below or with the 278 ## background color. 279 for k = 1:numel (lev) 280 lvl_idx = find (abs (cont_lev - lev(k)) < lvl_eps); 281 len = numel (lvl_idx); 282 if (len > 1) 283 mark = false (size (lvl_idx)); 284 a = 1; 285 while (a < len) 286 ## take 1st patch 287 pa_idx = lvl_idx(a); 288 ## get pointer to contour start, and contour length 289 curr_ct_idx = cont_idx(pa_idx); 290 curr_ct_len = cont_len(pa_idx); 291 ## get contour 292 curr_ct = c(:, curr_ct_idx:curr_ct_idx+curr_ct_len-1); 293 b_vec = (a+1):len; 294 next_ct_pt_vec = c(:, cont_idx(lvl_idx(b_vec))); 295 in = inpolygon (next_ct_pt_vec(1,:), next_ct_pt_vec(2,:), 296 curr_ct(1, :), curr_ct(2, :)); 297 mark(b_vec(in)) = ! mark(b_vec(in)); 298 a += 1; 299 endwhile 300 if (numel (mark) > 0) 301 ## All marked contours describe a hole in a larger contour of 302 ## the same level and must be filled with color of level below. 303 ma_idx = lvl_idx(mark); 304 if (k > 1) 305 ## Find color of level below. 306 tmp = find (abs (cont_lev - lev(k - 1)) < lvl_eps); 307 lvl_bel_idx = tmp(1); 308 ## Set color of patches found. 309 cont_lev(ma_idx) = cont_lev(lvl_bel_idx); 310 else 311 ## Set lowest level contour to NaN. 312 cont_lev(ma_idx) = NaN; 313 endif 314 endif 315 endif 316 endfor 317 318 ## The algorithm can create patches with the size of the plotting 319 ## area, we would like to draw only the patch with the highest level. 320 del_idx = []; 321 max_idx = find (cont_area == max (cont_area)); 322 if (numel (max_idx) > 1) 323 ## delete double entries 324 del_idx = max_idx(1:end-1); 325 cont_area(del_idx) = cont_lev(del_idx) = []; 326 cont_len(del_idx) = cont_idx(del_idx) = []; 327 endif 328 329 ## Now we have everything together and can start plotting the patches 330 ## beginning with largest area. 331 [~, svec] = sort (cont_area); 332 len = ncont - numel (del_idx); 333 h = []; 334 for n = len:-1:1 335 idx = svec(n); 336 ctmp = c(:, cont_idx(idx):cont_idx(idx) + cont_len(idx) - 1); 337 if (all (ctmp(:,1) == ctmp(:,end))) 338 ## patch() doesn't need/want closed contour. It will do it itself. 339 ctmp(:,end) = []; 340 endif 341 if (isnan (cont_lev(idx))) 342 fc = get (ca, "color"); 343 if (strcmp (fc, "none")) 344 fc = get (ancestor (ca, "figure"), "color"); 345 endif 346 else 347 fc = "flat"; 348 endif 349 h = [h; __go_patch__(ca, "xdata", ctmp(1, :)(:), "ydata", ctmp(2, :)(:), 350 "vertices", ctmp.', "faces", 1:(cont_len(idx)-1), 351 "facevertexcdata", cont_lev(idx), "facecolor", fc, 352 "cdata", cont_lev(idx), "edgecolor", lc, 353 "linestyle", ls, "linewidth", lw, 354 "parent", hg)]; 355 endfor 356 357 if (min (lev) == max (lev)) 358 set (ca, "clim", [min(lev)-1, max(lev)+1], "layer", "top"); 359 else 360 set (ca, "clim", [min(lev), max(lev)], "layer", "top"); 361 endif 362 else 363 ## Decode contourc output format. 364 h = []; 365 i = 1; 366 while (i < length (c)) 367 clev = c(1,i); 368 clen = c(2,i); 369 370 if (all (c(:,i+1) == c(:,i+clen))) 371 p = c(:, i+1:i+clen-1).'; 372 else 373 p = [c(:, i+1:i+clen), NaN(2, 1)].'; 374 endif 375 376 switch (zmode) 377 case "none" 378 h = [h; __go_patch__(ca, "xdata", p(:,1), "ydata", p(:,2), 379 "zdata", [], 380 "vertices", p, "faces", 1:rows (p), 381 "facevertexcdata", clev, "facecolor", "none", 382 "cdata", clev, "edgecolor", lc, 383 "linestyle", ls, "linewidth", lw, 384 "parent", hg)]; 385 case "auto" 386 h = [h; __go_patch__(ca, "xdata", p(:,1), "ydata", p(:,2), 387 "zdata", clev * ones (rows (p),1), 388 "vertices", [p, clev * ones(rows(p),1)], 389 "faces", 1:rows(p), 390 "facevertexcdata", clev, "facecolor", "none", 391 "cdata", clev, "edgecolor", lc, 392 "linestyle", ls, "linewidth", lw, 393 "parent", hg)]; 394 otherwise 395 h = [h; __go_patch__(ca, "xdata", p(:,1), "ydata", p(:,2), 396 "zdata", zlev * ones (rows (p), 1), 397 "vertices", [p, zlev * ones(rows(p),1)], 398 "faces", 1:rows (p), 399 "facevertexcdata", clev, "facecolor", "none", 400 "cdata", clev, "edgecolor", lc, 401 "linestyle", ls, "linewidth", lw, 402 "parent", hg)]; 403 endswitch 404 i += clen + 1; 405 endwhile 406 endif 407 408 set (ca, "climmode", climmode); 409 410endfunction 411 412function update_zlevel (h, ~) 413 414 z = get (h, "zlevel"); 415 zmode = get (h, "zlevelmode"); 416 kids = get (h, "children"); 417 418 switch (zmode) 419 case "none" 420 set (kids, "zdata", []); 421 case "auto" 422 for i = 1 : length (kids) 423 set (kids(i), "zdata", get (kids(i), "cdata") .* 424 ones (size (get (kids(i), "xdata")))); 425 endfor 426 otherwise 427 for i = 1 : length (kids) 428 set (kids(i), "zdata", z .* ones (size (get (kids(i), "xdata")))); 429 endfor 430 endswitch 431 432endfunction 433 434function update_line (h, ~) 435 436 lc = get (h, "linecolor"); 437 if (strcmp (lc, "auto")) 438 lc = "flat"; 439 endif 440 set (findobj (h, "type", "patch"), "edgecolor", lc, 441 "linewidth", get (h, "linewidth"), "linestyle", get (h, "linestyle")); 442 443endfunction 444 445function update_data (h, ~, prop = "") 446 persistent recursive = false; 447 448 if (! recursive) 449 recursive = true; 450 451 delete (get (h, "children")); 452 453 switch (prop) 454 case "levellist" 455 set (h, "levellistmode", "manual"); 456 case "levelstep" 457 set (h, "levelstepmode", "manual"); 458 case "fill" 459 ## Switching from filled ('k' linespec) to unfilled, reset linecolor 460 if (strcmp (get (h, "fill"), "off")) 461 set (h, "linecolor", "auto"); 462 else 463 set (h, "linecolor", "black"); 464 endif 465 endswitch 466 467 if (strcmp (get (h, "levellistmode"), "manual") 468 && ! strcmp (prop, "levelstep")) 469 lvl = get (h, "levellist"); 470 elseif (strcmp (get (h, "levelstepmode"), "manual")) 471 z = get (h, "zdata"); 472 lvs = get (h, "levelstep"); 473 lvl(1) = ceil (min (z(:)) / lvs) * lvs; 474 lvl(2) = floor (max (z(:)) / lvs) * lvs; 475 if (lvl(1) >= lvl(2)) 476 lvl = median (z(:)); 477 else 478 lvl = lvl(1) : lvs : lvl(2); 479 endif 480 set (h, "levellist", lvl, "levellistmode", "auto"); 481 else 482 z = get (h, "zdata"); 483 ## FIXME: The levels should be determined similarly to {x,y,z}ticks 484 ## so that they aren't set at extremely odd values. 485 lvl = linspace (min (z(! isinf (z))), max (z(! isinf (z))), 10 + 2); 486 ## Strip off max outlier, min must stay for contourf hole algorithm. 487 lvl = lvl(1:end-1); 488 endif 489 490 if (strcmp (get (h, "fill"), "on")) 491 X = get (h, "xdata"); 492 Y = get (h, "ydata"); 493 Z = get (h, "zdata"); 494 if (isvector (X) || isvector (Y)) 495 [X, Y] = meshgrid (X, Y); 496 endif 497 [nr, nc] = size (Z); 498 X0 = prepad (X, nc+1, 2 * X(1, 1) - X(1, 2), 2); 499 X0 = postpad (X0, nc+2, 2 * X(1, nc) - X(1, nc - 1), 2); 500 X0 = [X0(1, :); X0; X0(1, :)]; 501 Y0 = prepad (Y, nr+1, 2 * Y(1, 1) - Y(2, 1), 1); 502 Y0 = postpad (Y0, nr+2, 2 * Y(nr, 1) - Y(nr - 1, 1)); 503 Y0 = [Y0(:, 1), Y0, Y0(:, 1)]; 504 Z0 = -Inf (nr+2, nc+2); 505 Z0(2:nr+1, 2:nc+1) = Z; 506 [c, lev] = contourc (X0, Y0, Z0, lvl); 507 else 508 [c, lev] = contourc (get (h, "xdata"), get (h, "ydata"), 509 get (h, "zdata"), lvl); 510 endif 511 set (h, "contourmatrix", c); 512 513 if (strcmp (get (h, "levellistmode"), "manual")) 514 ## Do nothing 515 elseif (strcmp (get (h, "levelstepmode"), "manual")) 516 set (h, "levellist", lev); 517 else 518 set (h, "levellist", lev); 519 lvlstep = sum (diff (lvl)) / (length (lvl) - 1); 520 set (h, "levelstep", lvlstep); 521 endif 522 523 add_patch_children (h); 524 update_text (h); 525 recursive = false; 526 endif 527 528endfunction 529 530function update_text (h, ~, prop = "") 531 persistent recursive = false; 532 533 if (! recursive) 534 recursive = true; 535 536 delete (findobj (h, "type", "text")); 537 538 switch (prop) 539 case "textlist" 540 set (h, "textlistmode", "manual"); 541 case "textstep" 542 set (h, "textstepmode", "manual"); 543 endswitch 544 545 if (strcmp (get (h, "textlistmode"), "manual")) 546 lvl = get (h, "textlist"); 547 elseif (strcmp (get (h, "textstepmode"), "manual")) 548 lev = get (h, "levellist"); 549 550 lvl_eps = get_lvl_eps (lev); 551 552 stp = get (h, "textstep"); 553 t = [0, floor(cumsum(diff (lev)) / (abs(stp) - lvl_eps))]; 554 lvl = lev([true, t(1:end-1) != t(2:end)]); 555 set (h, "textlist", lvl); 556 else 557 lvl = get (h, "levellist"); 558 set (h, "textlist", lvl, "textstep", get (h, "levelstep")); 559 endif 560 561 if (strcmp (get (h, "showtext"), "on")) 562 switch (get (h, "zlevelmode")) 563 case "manual" 564 __clabel__ (get (h, "contourmatrix"), lvl, h, 565 get (h, "labelspacing"), get (h, "zlevel")); 566 case "auto" 567 __clabel__ (get (h, "contourmatrix"), lvl, h, 568 get (h, "labelspacing"), "auto"); 569 otherwise 570 __clabel__ (get (h, "contourmatrix"), lvl, h, 571 get (h, "labelspacing"), []); 572 endswitch 573 endif 574 575 recursive = false; 576 endif 577 578endfunction 579 580function lvl_eps = get_lvl_eps (lev) 581 582 ## FIXME: is this the right thing to do for this tolerance? Should 583 ## it be an absolute or relative tolerance, or switch from one to the 584 ## other depending on the value of lev? 585 if (isscalar (lev)) 586 lvl_eps = abs (lev) * sqrt (eps) + sqrt (eps); 587 else 588 tmp = min (abs (diff (lev))); 589 if (tmp < 10*eps) 590 lvl_eps = sqrt (eps); 591 else 592 lvl_eps = tmp / 1000.0; 593 endif 594 endif 595 596endfunction 597