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