octave_packages/m/plot/private/__contour__.m

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