This file is indexed.

/usr/share/octave/packages/mapping-1.2.1/shapedraw.m is in octave-mapping 1.2.1-3.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
## Copyright (C) 2014,2015 Philip Nienhuis
## 
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
## 
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
## 
## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <http://www.gnu.org/licenses/>.

## -*- texinfo -*- 
## @deftypefn  {Function File} [@var{h}] = shapedraw (@var{shp})
## @deftypefnx {Function File} [@var{h}] = shapedraw (@var{shp}, @var{clr})
## @deftypefnx {Function File} [@var{h}] = shapedraw (@var{shp}, @var{clr}, ...)
## @deftypefnx {Function File} [@var{h}] = shapedraw (@var{shp}, ...)
## Plot contents of a shapefile, a map-/geostruct or a struct made by
## shaperead.m.
##
## @var{shp} can be a shapefile (will be read and plotted), a struct made by 
## shaperead, or a mapstruct or geostruct, possibly made by some other
## function.  shapeplot.m will try to determine the type.  Points, MultiPoints,
## Polylines, Polygons and MultiPatch shape features can be plotted.
##
## The optional argument @var{clr} can be a predefined color name ('b', 'green',
## etc.) or an RGB triplet.  The default is [0.6, 0.6, 0.6] which produces a
## grey plot.  Polygons and MultiPatches can also be plotted as solid patches;
## then @var{clr} needs to have a second row indicating the fill color.  Octave
## does not support transparent fills yet.  Single-character color codes can be
## combined with linestyle indicators ":", "-", "--", "-.", ".-" and/or marker
## style indicators
## "*", ".", "+", "@@", "v", "^", ">", "<", "d", "h", "o", "p", "s"
## to modify the linestyle for polylines.
##
## Other graphics properties for drawing can be supplied either instead of,
## or after the color argument and will be conveyed as-is to the actual drawing
## routines.  Depending on shapetype, the following proqperties are accepted:
##
## @itemize
## @item All shape types: 
## Visible, LineStyle, LineWidth, Marker, MarkerEdgeColor, MarkerFaceColor,
## MarkerSize
##
## @item Point, MultiPoints, Line, Polyline: 
## Color
##
## @item Polygon, MultiPatch: 
## FaceColor, EdgeColor
## @end itemize
##
## Polygons with holes can be properly plotted provided the holes are separate
## shape feature parts comprising counterclockwise polylines; the first partial
## feature must be the clockwise outer polygon.  The Octave-Forge geometry
## package is required to assess whether multipart polygons have holes and to
## properly draw them.  shapedraw.m will search for the geometry package the
## first time it is instructed to plot filled polygons.  To initiate a new
## search later on (e.g., after the geometry package has been loaded), simply
## invoke shapedraw without any arguments.
##
## Optional output argument @var{h} is the figure handle of the plot.
##
## @seealso{geoshow, mapshow, shapeinfo, shaperead}
## @end deftypefn

## Author: Philip Nienhuis <prnienhuis@users.sf.net>
## Created: 2014-11-15
## Updates:
## 2014-11-15 First draft
## 2014-11-17 Add support for geostructs
## 2014-11-29 Solid fill (patches)
## 2014-12-16 Solid fill (polygons) with holes
## 2014-12-17 MultiPatch plots (3D)
## 2014-12-18 3D Points/Polylines/-gons
## 2014-12-30 Transpose npr values for MultiPatch prior to drawing patches
## 2015-01-01 Fix function name in error message
## 2015-01-06 Improve plotting of polygons with holes - optimize branch cuts
## 2015-01-07 Largely rewritten; 
##            - Combined line and patch sections into one switch
##            - More rigorous checks on input parameters
## 2015-01-08 More rewrites, better input checks
##     ''     Allow marker indicators for (multi)point
## 2015-01-09 Add code to check if individual polygon parts are holes (not all)
##     ''     Improve code for single color code/stle arguments
## 2015-01-10 Renamed to shapedraw.m (shapeplot is already in OF-geometry)
##     ''     Fixed argument checking for "color" property (swapping if checks)
## 2015-01-20 Change "color" to "edgecolor" for multipatch
##     ''     Add undocumented "center" argument to cope with large coordinates
##            (OpenGL chokes there as it only works with single precision)
## 2015-01-21 Apply varargin to 3D-patches using "set" command
##     ''     Simplify input args, eliminate "clr" arg
##     ''     Support for extended mapstructs
## 2015-01-24 Restructured a bit, debugged MultiPatch drawing
## 2015-01-27 Texinfo header, check on Z-values
## 2015-01-30 Simplify ML style struct test; allow Point types (no BoundingBox)
## 2015-01-31 Fix wrong indexing in MultiPatch-triangle processing
## 2015-02-03 Swap checks for first color arg and graphics properties
## 2015-02-04 Check for Z and make it a column vector before calling plot3
##     ''     Morefixes for color argument checks
## 2015-02-11 Eliminate duplicate code, move to subfunc chk_props
## 2015-02-17 Markerstyle => Marker in point draw switch stmt
##     ''     Don't reshape back args in chk_props subfuc
##     ''     Improve checks for default color
## 2015-04-19 Make warning state changes local
## 2015-07-10 Add try-catch around varargin reshape to catch wrong input
##     ''     Simplify polygon plot code
##     ''     Fix 3D-plotting for "extended" map/geostructs
## 2015-12-27 Improve speed of drawing many polygons
## 2015-12-30 Improve multipatch drawing speed

function [h] = shapedraw (shpf, varargin)

  ## Keep track of OF geometry package presence
  persistent ofgeom = [];

  ## What this function does in what order:
  ## 1. Input checks & analysis.
  ##    Get a struct & assess its type; or read data from .shp file set.
  ## 2. Morph XY[Z} in a suitable form for faster plotting (avoid for loops).
  ##    On the way, explore multipart polygons and if found, get separators.
  ## 3. Process input args; first of all marker/line & fill colors & style.
  ##    A bit involved as there are many interfering options.
  ## 4. Find out if plot is on top of existing figure or must be a new figure.
  ## 5. Depending on geometry type, draw the shape features.
  ##    For Polygons, explore and process "holes".
  ##    FIXME Holes must be implemented for Multipatch as well - pending

  ## 1. Input checks & analyses

  if (nargin < 1)
    if (! isempty (ofgeom))
      ## Reset ofgeom
      ofgeom = [];
      return;
    endif
    print_usage ()
  endif
  
  if (isempty (shpf))
    ## Nothing to plot
    return;

  if (isempty (varargin))
    ## Supply default color
    varargin = {[0.6, 0.6, 0.6]};
  endif

  elseif (isstruct (shpf))
    ## Yep. Find out what type
    fldn = fieldnames (shpf);
    if (ismember ("vals", fldn) && ismember ("shpbox", fldn))
      ## Assume it is an Octave-style struct read by shaperead
      type = 0;
    elseif (all (ismember ({"Geometry", "X", "Y"}, fldn)))
      ## Assume it is a Matlab-style mapstruct
      type = 1;
    elseif (all (ismember ({"Geometry", "Lat", "Lon"}, fldn)))
      ## Assume it is a Matlab-style geostruct
      type = 2;
    else
      ## Not a supported struct type
      error ("shapedraw: unsupported struct type.\n")
    endif

  elseif (ischar (shpf))
    ## Filename? let's see:
    [pth, fnm, ext] = fileparts (shpf);
    if (isempty (ext) || ! strcmpi (ext, ".shp"))
      ## Just add a .shp suffix
      shpf = [shpf ".shp"];
    endif
    fid = fopen (shpf);
    if (fid < 0)
      error ("File %s not found", shpf);
    else
      fclose (fid);
    endif
    ## Read shapefile in ML style into mapstruct
    shpf = shaperead (shpf, 1);
    type = 1;
  else

    error ("struct name of file name expected");
  endif

  ## 2. Morph XY[Z} data in a suitable form for fastest plotting

  ## Prepare XY plot. Get vertices & prepare some geometry data
  mp = [];
  if (type == 0)
    ## Octave style with intermediate NaN rows. Easy except MultiPatch
    geom = shpf.vals(1, 6);
    X = shpf.vals(:, 1);
    Y = shpf.vals(:, 2);
    Z = shpf.vals(:, 3);
    ## Find multipart shape features to be able to plot polygons with holes
    if (ismember (shpf.vals(1, 6), [5 15 25 31]))
      ## Yep we have polygons
      jdx = [ 0 find(isnan (shpf.vals(:, 1))') ] + 1;
      idx = shpf.idx';
      if (numel (jdx) > numel (idx))
        ## Multipart features present. Find them
        mp = find (diff (find (ismember (jdx, idx))) - 1);
        ## The last feature is a bit special...
        if (jdx(end) > idx(end))
          mp = [ mp numel(idx) ];
        endif
      endif
      npr = shpf.npr;
    endif
        
  elseif (type == 1 || type == 2)
    ## Matlab struct type
    switch lower (shpf(1).Geometry)
      case "point"
        geom = 1;
      case "multipoint"
        geom = 8;
      case {"polyline", "line"}
        geom = 3;
      case "polygon"
        geom = 5;
      case "multipatch"
        geom = 31;
      otherwise
    endswitch
    tmp = cell (1, 2*numel (shpf));
    if (type == 1)
    ## mapstruct. Needs some preparation: insert a NaN after each subshape
      [tmp(1:2:2*numel(shpf))] = deal ({shpf.X});
      [tmp(2:2:2*numel(shpf)-1)] = NaN;
      X = [tmp{:}];
      [tmp(1:2:2*numel(shpf))] = deal ({shpf.Y});
      [tmp(2:2:2*numel(shpf)-1)] = NaN;
      Y = [tmp{:}];
    else
      ## geostruct. Needs some preparation: insert a NaN after each subshape
      [tmp(1:2:2*numel(shpf))] = deal ({shpf.Lon});
      [tmp(2:2:2*numel(shpf)-1)] = NaN;
      X = [tmp{:}];
      [tmp(1:2:2*numel(shpf))] = deal ({shpf.Lat});
      [tmp(2:2:2*numel(shpf)-1)] = NaN;
      Y = [tmp{:}];
    endif
    ## For "extended" (not strictly ML-compatible) mapstructs:
    if (isfield (shpf, "Z"))
      ## Give it a try. Could be one Z per shape, rather than one per vertex
      ## if Z values came in from xBase (dbf) rather than .shp
      [tmp(1:2:2*numel(shpf))] = deal ({shpf.Z});
      [tmp(2:2:2*numel(shpf)-1)] = NaN;
      Z = [tmp{:}];
      ## Check nr. of Z values
      if (numel (Z) == numel (X))
        ## Z values apparently originate from .shp
        if (isfield (shpf, "Parts"))
          ## Multipatch. Get Parts (contains multipatch subtypes)
          geom = 31;
          npr = {shpf.Parts};
        else
          ## Bump shape type to Z-type equivalent. Could've been 1/3/5/8 or 23...
          geom = mod (geom, 10) + 10;
        endif
      else
        Z = [];
      endif
    else
      Z = [];
    endif
    ## Find multipart shape features to be able to plot polygons with holes
    if (mod (geom, 10) == 5 || geom == 31)
      ## Polygons
      jdx = [ 0 find(isnan (X)) ] + 1;
      ## A little trick....
      idx = [ 0 (cumsum (cellfun (@(x) numel(x), {shpf(1:end-1).X})) + ...
              cumsum (ones (1, numel (shpf) - 1))) ] + 1;
      if (numel (jdx) > numel (idx))
        ## Multipart features present. Find them
        mp = find (diff (find (ismember (jdx, idx))) - 1);
        ## The last feature is a bit special...
        if (jdx(end) > idx(end))
          mp = [ mp numel(idx) ];
        endif
      endif
    endif
  endif

  ## Check for "verbose" arg
  verbose = find (strncmpi (varargin, "ver", 3));
  varargin (verbose) = [];
  verbose = ! isempty (verbose);

  ## If requested, reduce coordinates (subtract means). May help to get better
  ## plots for large coordinates as OpenGL (w. single accuracy) fouls those up
  center = find (strncmpi (varargin, "cen", 3));
  varargin (center) = [];
  if (! isempty (center))
    X -= mean (X(isfinite (X)));
    Y -= mean (Y(isfinite (Y)));
    if (ismember (geom, [13, 15, 18, 31]))
      Z -= mean (Z(isfinite (Z)));
    endif
  endif

  ## 3. Process input args; mainly marker/line & fill colors

  ## Predefine some properties to be able to check later
  pt_properties = {"Marker", "Color", "LineStyle", "MarkerEdgeColor", ...
                   "MarkerFaceColor", "MarkerSize", "Visible"};
  pl_properties = {"Color", "LineStyle", "LineWidth", "Visible"};
  pg_properties = {"FaceColor", "FaceAlpha", "LineStyle", "LineWidth", ...
                   "EdgeColor", "EdgeAlpha", "Visible"};
  color_codes = { ...
  "b", "k", "c", "g", "m", "y", "w", "r", "blue", "black", "cyan", "green", ...
  "magenta", "yellow", "white", "red", "blue"};
  clrptn = [ '(^[kbcgmprwy](--|\.-|-\.|[:\-osdhv^<>*\.\+@])' '|' ...
             '^([:\-osdhv^<>*\.\+@]|--|\.-|-\.)[bcgkmprwy])' ];

  ## Check color arg
  ## FIXME support for colormaps would be great...
  if (isempty (varargin))
    ## Restore default color
    if (ismember (geom, [5, 15, 25, 31]))
      ## polygons and multipatches
      varargin = {"facecolor" [0.6 0.7 0.9] "edgecolor" [0.5, 0.5, 0.5]};
    else
      varargin = { "color" [0.6, 0.6, 0.6] };
    endif

  elseif (ischar (varargin{1}))
    ## Check if arg#1 is a color code
    if (ismember (varargin{1}(1, :), color_codes) ||
        ! isempty (cell2mat (regexp (varargin{1}(1, :), clrptn, "match"))))
      ## Color OK
      clrtyp = 1;
      if (size (varargin{1}, 1) == 1 && ...
          (size (varargin{1}, 2) == 2 || size (varargin{1}, 2) == 3))
        ## Line- or markerstyle color. Draw polygons as polylines ([3/13/23]
        if (ismember (geom, [5, 15, 25]))
          geom -= 2;
        endif
      elseif (ismember (geom, [5, 15, 25, 31]))
        colr = varargin{1};
        varargin(1) = [];
        varargin = {"edgecolor" colr(1, :) varargin{:}};
        if (size (colr, 1) > 1)
          varargin = {"facecolor" colr(2, :) varargin{:}};
        endif
      else
        varargin = {"color" varargin{1}(1, :) varargin{2:end}};
      endif
    ## Check if it's a graphics property
    elseif (ismember (lower (varargin{1}(1, :)), lower (unique ([pt_properties, ...
                                      pl_properties, pg_properties]))))
      ## Yes. Do some checks
      if (numel (varargin) > 1)
        props = lower ({reshape(varargin, 2, []){1, :}});
      endif
      gr_props_ok = ismember (props, ...
                    unique (lower ([pt_properties pl_properties pg_properties])));
      if (! all (gr_props_ok))
        warning ("shapedraw: unknown graphics properties\n")
      endif
      ## Restore default color if required
      if (ismember (geom, [5, 15, 25, 31]))
        if (! ismember ("facecolor", props))
          ## Introduce polygon edge/face colors
          varargin = { varargin{:} "facecolor" [0.5 0.7 0.9]};
        endif
        if (! ismember ("edgecolor", props))
          ## Introduce polygon edge/face colors
          varargin = { varargin{:} "edgecolor" [0.6, 0.6, 0.6]};
        endif
      elseif (! any (ismember ({"color", "markeredgecolor"}, props)))
        varargin = { "color", [0.6, 0.6, 0.6], varargin{:} };
      endif
    else
      error ("shapedraw: incomprehensible arg#1: %s\n", varargin{1});
    endif

  elseif (isnumeric (varargin{1}))
    ## Assume it's a color RGB triplet, or array of triplets
    if (size (varargin{1}, 2) != 3)
      ## Too few or too many
      error ("shapeplot: incomprehensible argument #2\n");
    else
      if (ismember (geom, [5, 15, 25, 31]))
        if (size (varargin{1}, 1) < 2)
          varargin{1}(2, :) = [0.5 0.7 0.9];
        endif
        varargin = {"edgecolor", varargin{1}(1, :) "facecolor" varargin{1}(2, :)};
      else
        varargin = {"color" clr(1, :) varargin{:}};
      endif
    endif
  endif

  ## 4. Find out if plot is on top of existing figure or must be a new figure

  hh = ishold ();
  ## ishold automatically creates an empty figure if none found
  hold ("on");
  if (nargout > 0)
    h = get (gcf);
  endif

  ## 5. Depending on geometry type, draw the shape features

  ## For each shape type check drawing options
  if (numel (varargin) > 1)
    try
      varargin = reshape (varargin, 2, []);
    catch
      varargin
      error ("shapedraw.m: ^^^^ looks like invalid drawing properties were entered");
    end_try_catch
  endif

  switch geom
    case {1, 21, 8, 28}                                     ## (Multi-)Points[M]
      varargin = chk_props (varargin, pt_properties, color_codes, clrptn, "(Multi)Point");
      ## Check for some marker style or marker color code (latter always arg #1)
      if (! ismember ("marker", lower (varargin(1, :))) && ...
          ! ismember (varargin(1), color_codes))
        varargin = { "marker", "." varargin{:} };
      endif
      plot (X, Y, varargin{:});

    case {11, 18}                                           ## (Multi-)PointsZ
      varargin = chk_props (varargin, pt_properties, color_codes, clrptn, "(Multi)PointZ");
      ## Check for some marker style or marker color code (latter always arg #1)
      if (! ismember ("marker", lower (varargin(1, :))) && ...
          ! ismember (varargin(1), color_codes))
        varargin = {"marker", "." varargin{:} };
      endif
      plot3 (X, Y, Z, varargin{:});

    case {3, 13, 23}                                        ## Polylines
      varargin = chk_props (varargin, pl_properties, color_codes, clrptn, "(Poly)Lines");
      ## Make sure we have an Nx2 matrix
      X = [X(:)]; 
      Y = [Y(:)];
      if (geom == 13 && ! isempty (Z))
        Z = [Z(:)];
        plot3 (X, Y, Z, varargin{:});
      else
        Z = [];
        plot (X, Y, varargin{:});
      endif

    case {5, 15, 25}                                        ## Polygons
      varargin = chk_props (varargin, pg_properties, color_codes, clrptn, "Polygons");
      ## Make sure we have an Nx2 matrix
      X = [X(:)]; 
      Y = [Y(:)];
      if (! isempty (Z))
        Z = [Z(:)];
      endif
      idx = [ idx (numel(X)+2) ];
      ## Check for holes in multipart poygon. Outer polygon must be clockwise
      has_holes = 0;
      ## Check if we have the OF geometry pacakage loaded
      if (isempty (ofgeom))
        ofgeom = ! isempty (which ("polygonArea"));
        if (! ofgeom)
          warning ("shapedraw: function 'polygonArea' not found.\n");
          printf  ("           (OF geometry package installed and loaded?)\n");
          printf  ("           => Holes in polygons will be filled\n");
        endif
      endif
      ipt = [ 0; (find (isnan (X))); (numel (X) + 1) ];
      if (ofgeom && ! isempty (mp))
        hdx = cell (mp, 1);
        ## Search multipart polygons for holes, backwards to save mp order
        for ii=numel(mp):-1:1
          ## Depending on struct type, get shape feature
          XT = X(idx(mp):idx(mp+1)-2);
          YT = Y(idx(mp):idx(mp+1)-2);
          jpt = [ 0; (find (isnan (XT))); (numel (XT) + 1) ];
          ## Check for hole
          h_idx = [];
          for jj=2:numel (jpt) - 1
            ## A clockwise polygon has area < 1, counterclockwise > 1
            hole = (polygonArea (XT(jpt(jj)+1:jpt(jj+1)-1), ...
                                 YT(jpt(jj)+1:jpt(jj+1)-1)) > 1);
            has_holes = has_holes || hole;
            if (hole)
              h_idx = [ h_idx jj ];
            endif
          endfor
          if (isempty (h_idx))
            mp(ii) = [];
          else
            hdx(mp) = h_idx;
          endif
        endfor
      endif

      ## Draw polygon parts one by one
      if (! ofgeom || isempty (mp))
        ## If no holes can be left open, draw each feature part separately
        idx = [ jdx (numel(X)+2) ];
      endif
      ## Prepare 'faces' argument for patch
      faces = NaN (numel(idx)-1, max (diff (find (isnan ([NaN; X; NaN])))));
      for ii=1:numel (idx) - 1
        XT = X(idx(ii):idx(ii+1)-2);
        YT = Y(idx(ii):idx(ii+1)-2);
        if (! isempty (Z))
          ZT = Z(idx(ii):idx(ii+1)-2);
        endif
        if (has_holes && ofgeom && ismember (ii, mp))
          ## To be able to not fill holes in polygons we need a trick. That
          ## comprises connecting the holes to the outline through "branch
          ## cuts" that do not run across other holes or even the outline.
          ## To obtain the latter we need to find optimal locations for the
          ## vertices on both ends of the branch cut.
          [XT, YT, ZT] = optimize_branch_cuts (XT, YT, hdx{ii});
          X(idx(ii):idx(ii+1)-2) = XT;
          Y(idx(ii):idx(ii+1)-2) = YT;
          if (! isempty (Z))
            Z(idx(ii):idx(ii+1)-2) = ZT;
          endif
          faces (ii, 1:numel(XT)) = [1:numel(XT)] + idx(ii) - 1;
        else
          faces (ii, 1:numel(XT)) = [idx(ii):idx(ii+1)-2];
        endif
      endfor
      if (geom == 5)
        p_h = patch ("vertices", [X Y],   "faces", faces, ...
                     "facecolor", "none", "edgecolor", "none");
      else
        p_h = patch ("vertices", [X Y Z], "faces", faces, ...
                     "facecolor", "none", "edgecolor", "none");
      endif
      set (p_h, varargin{:});

    case {31}                                               ## MultiPatch
      ## Multipatch. Process individual subfeatures one by one
      ## First check drawing options
      varargin = chk_props (varargin, pg_properties, color_codes, clrptn, ...
                            "MultiPatch");
      ## Add a NaN row index so below loops will work up to max loop counter
      idx = [ idx numel(X)+2 ];
      if (verbose)
        printf ("\n");
      endif
      ## Make sure we have column vectors
      X = [X(:)];
      Y = [Y(:)];
      Z = [Z(:)];
      ## Init patch 'faces' ptrs. First find out nr. of triangle strips/fans
      prts = [shpf.Parts];
      it = find (prts(2, :) <= 1);
      ## Preallocate 'faces' argument for subfeature patches.
      faces = NaN (size (prts, 2), max (diff (prts(1, :))-1));
      clear prts;                                           ## can be large
      ipf = 0;                                              ## ptr into faces
      for jj=1:numel (npr)
        ## For each feature
        sep = [ (npr{jj}(1, :)) (diff (idx)(jj)) ];
        ibase = idx(jj) - 1;
        lastyp = [];
        for ii=1:size (npr{jj}, 2)
          ## For each subfeature
          if (verbose)
            printf ("Shape part %d, subpart %d\r", jj, ii);
          endif
          switch npr{jj}(2, ii)
            case 0                                            ## Triangle strip
              faces(++ipf, 1:3) = [sep(ii)+1:sep(ii)+3] + ibase;
              ## Each following point forms next triangle with previous two.
              for kk=sep(ii)+4:sep(ii+1)-1
                faces(++ipf, 1:3) = [kk-2:kk] + ibase;
              endfor
              lastyp = 0;
            case 1                                            ## Triangle fan
              faces(++ipf, 1:3) = [sep(ii)+1:sep(ii)+3] + ibase;
              ## Each following point forms next triangle with first and last
              ## points
              for kk=sep(ii)+4:sep(ii+1)-1
                faces(++ipf, 1:3) = [1 kk-1:kk] + ibase;
              endfor
              lastyp = 1;
            case 2                                            ## Outer ring
              ## FIXME Plotting holes hasn't been implemented yet for MultiPatch
              #        (although it seems to work with files found in the wild)
              faces(++ipf, 1:(sep(ii+1)-sep(ii)-1)) = ...
               [sep(ii)+1:sep(ii+1)-1] + ibase;
              lastyp = 2;
            case 3                                            ## Inner ring
              ## FIXME Plotting holes hasn't been implemented yet for MultiPatch
              ##       (although it seems to work with files found in the wild)
              ## FIXME untested code
              faces(++ipf, 1:(sep(ii+1)-sep(ii)-1)) = ...
               [sep(ii)+1:sep(ii+1)-1] + ibase;
              lastyp = 3;
            case 4                                            ## First ring
              ## FIXME untested code
              faces(++ipf, 1:(sep(ii+1)-sep(ii)-1)) = ...
               [sep(ii)+1:sep(ii+1)-1] + ibase;
            case 5                                            ## Ring
              ## FIXME untested code
              faces(++ipf, 1:(sep(ii+1)-sep(ii)-1)) = ...
               [sep(ii)+1:sep(ii+1)-1] + ibase;
              lastyp = 4;
            otherwise
          endswitch
        endfor
        vals = [];
      endfor
      ## Just to be sure, weed out Nan rows out of faces
      faces(isnan (faces(:, 1)), :) = [];
      if (! all (isnan (faces(:))))
        p_h = patch ("vertices", [X Y Z], "faces", faces, ...
                     "facecolor", "none", "edgecolor", "none");
        set (p_h, varargin{:});
      endif
      if (verbose)
        printf ("\n");
      endif

    otherwise
      ## Other shape feature type?
  endswitch

  axis equal;

  if (! hh)
    hold ("off");
  endif

endfunction


##-----------------------------------------------------------------------------
## Copyright (C) 2015 Philip Nienhuis
##
## Optimize branch cuts for inner polygons (holes)

function [XX, YY, ZZ, kdx, mdx] = optimize_branch_cuts (XX, YY, varargin);

  ## Check args #3 and #4
  ZZ = [];
  hdx = {};
  if (! isempty (varargin))
    if (nargin == 3)
      ## Assume it's a hdx file indicating which subshapes are holes
      hdx = varargin{1};
    elseif (nargin == 4)
      ## Assume they're Z coordinates
      ZZ = varargin{1};
      hdx = varargin{2};
    endif
  endif

  ## Also keep track of Z. Z isn't (yet) in the branch cut optimization (but
  ## that could be done easily)
  if (isempty (ZZ))
    ## At least provide pointers where Z coordinates have gone in output arrays
    ZZ = [1:numel(XX)]';
  endif
  XY = [ XX YY ZZ];

  ## Find NAN separators
  ipt = [0; find(isnan (XX)); numel(XX)+1];
  if (numel (ipt) < 3)
    ## No NaN separators => no subfeatures. Return
    return
  endif
  for ii=1:numel (ipt) - 1
    XY(ipt(ii)+1:ipt(ii+1)-1, 4) = [ ipt(ii)+1:ipt(ii+1)-1 ]';
  endfor
  if (isempty (hdx))
    ## Assume all except first subfeatures are holes
    hdx = [ 2:numel(ipt)-1 ];
  endif

  ## Silence broadcasting warning
  warning ("off",  "Octave:broadcast", "local");
  
  ## 1. Locate locate optimal branch cut vertices in outer and inner polygons
  ## FIXME Below we only relate holes to the outer polygon. That bears risk
  ##       for inner holes shielded from the outer boundary by other holes.
  ##       It's better to find the hole closest to the outer boundary, start
  ##       there, and then successively process each hole in turn and find the
  ##       smallest distance to either the outer polygon or already processed
  ##       holes. Such a strategy is neither robust nor easily devised, however.
  for ii=2:numel (ipt) - 1
    b1 = distancePoints (XY(ipt(1)+1:ipt(2)-1, :), XY(ipt(ii)+1:ipt(ii+1)-1, :));
    [~, col(ii-1)] = min (min (b1));
    [~, row] = min (b1(:, col(ii-1)));
    ## Copy vertices of outer polygon down to make room for marker
    XY(row+2:end+2, :) = XY(row:end, :);
    ## Add marker & pointer to inner polygon
    XY(row+1, 1) = Inf;
    XY(row+1, 2:3) = [ii ii];
    ipt(2:end) += 2;
  endfor

  ## 2. ---- Insert inner polygons into outer ----
  ## Update ipt
  ipt = [find(isnan (XY(:, 1))); size(XY, 1)+1];
  jpt = [find(isinf (XY(:, 1))) ; ipt(1)];
  XYXY = XY(1:jpt(1)-1, :);
  for ii=1:numel(jpt) - 1
    ## First rotate inner polygon vertices. Closed polygon => no cirshift
    tmp = XY(ipt(ii)+1:ipt(ii+1)-1, :);
    tmp = [ tmp(col(ii):end-1, :) ; tmp(1:col(ii), :) ];
    XYXY = [ XYXY ; tmp ];
    XYXY = [ XYXY ; XY(jpt(ii)+1:jpt(ii+1)-1, :) ];
  endfor

  XX = XYXY(:, 1);
  YY = XYXY(:, 2);
  if (nargout == 3)
    ZZ = XYXY(:, 3);
  endif

endfunction


##-----------------------------------------------------------------------------
## Copyright (C) 2015 Philip Nienhuis
##
## Check graphics properties for various geometries

function args = chk_props (args, properties, color_codes, clrptn, geometry)

  gr_props = ! ismember (lower (args(1, :)), [lower(properties) color_codes]);
  if (any (gr_props) && isempty (cell2mat (regexp (args(1, :), ...
                                                   clrptn, "match"))))
    printf ("shapedraw: illegal drawing options for %s ignored\n", geometry);
    printf ("(%s)\n", strjoin (args(1, :)(gr_props)));
  else
    gr_props = [];
  endif
  args(:, gr_props) = [];

endfunction