File Coverage

blib/lib/Math/Geometry/Delaunay.pm
Criterion Covered Total %
statement 20 21 95.2
branch n/a
condition n/a
subroutine 7 7 100.0
pod n/a
total 27 28 96.4


line stmt bran cond sub pod time code
1             package Math::Geometry::Delaunay;
2              
3 1     1   24558 use 5.008;
  1         6  
  1         51  
4 1     1   7 use warnings;
  1         3  
  1         39  
5 1     1   5 use strict;
  1         8  
  1         51  
6 1     1   6 use Carp qw(carp);;
  1         2  
  1         90  
7 1     1   6 use Exporter();
  1         2  
  1         74  
8              
9             our @ISA = qw(Exporter);
10             our $VERSION;
11              
12             BEGIN {
13 1     1   7 use XSLoader;
  1         2  
  1         61  
14 1     1   3 $VERSION = '0.19';
15 1         1608 XSLoader::load('Math::Geometry::Delaunay');
16 0           exactinit();
17             }
18              
19             use constant {
20             TRI_CONSTRAINED => 'Y',
21             TRI_CONFORMING => 'Dq0',
22             TRI_CCDT => 'q',
23             TRI_VORONOI => 'v',
24             };
25              
26             our @EXPORT_OK = qw(TRI_CONSTRAINED TRI_CONFORMING TRI_CCDT TRI_VORONOI);
27             our @EXPORT = qw();
28              
29             my $pi = atan2(1,1)*4;
30              
31             sub new {
32             my $class = shift;
33             my $self = {};
34             $self->{in} = Math::Geometry::Delaunay::Triangulateio->new();
35             $self->{out} = undef;
36             $self->{vorout} = undef;
37             $self->{poly} = {
38             regions => [],
39             holes => [],
40             polylines => [],
41             points => [],
42             segments => [],
43             outnodes => [], #for cache, first time C output arrays are imported
44             voutnodes => [], #for cache
45             segptrefs => {}, #used to avoid dups of points added with addSegments
46             };
47              
48             $self->{optstr} = '';
49             # Triangle switches
50             # prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh
51             # where __ is an optional number
52             $self->{a} = -1; # max tri area
53             $self->{q} = -1; # quality min angle
54             $self->{e} = 0; # produce edges switch
55             $self->{v} = 0; # voronoi switch
56             $self->{n} = 0; # neighbors switch
57             $self->{N} = 0; # suppress node output
58             $self->{E} = 0; # suppress element output
59             $self->{O} = 0; # suppress holes - ignore input holes
60             $self->{o2}= 0; # subparametric switch (for 6 pts/tri)
61             $self->{Q} = 1; # quiet - switch for Triangle's printed output
62             $self->{V} = 0; # verbose - from 0 to 3 for increasing verbosity
63              
64             bless $self, $class;
65             return $self;
66             }
67              
68             sub reset {
69             my $self = shift;
70             # clear input
71             $self->{poly}->{$_} = [] for qw(regions holes polylines points segments);
72             $self->{poly}->{segptrefs} = {};
73             # clear any previous output
74             $self->{poly}->{$_} = [] for qw(outnodes voutnodes);
75             }
76              
77             # triangulatio interfaces
78             sub in {return $_[0]->{in};}
79             sub out {return $_[0]->{out};}
80             sub vorout {return $_[0]->{vorout};}
81              
82             # getter/setters for the triangulate switches that take numbers
83              
84             sub area_constraint { # -1 for disabled
85             if (@_>1) {$_[0]->{a}=$_[1];}
86             return $_[0]->{a};
87             }
88             sub minimum_angle { # "q" switch, in degrees, -1 for disabled
89             if (@_>1) {$_[0]->{q}=$_[1];}
90             return $_[0]->{q};
91             }
92             sub subparametric {
93             if (@_>1) {$_[0]->{o2}=$_[1]?1:0;}
94             return $_[0]->{o2};
95             }
96             sub doEdges {
97             if (@_>1) {$_[0]->{e}=$_[1]?1:0;}
98             return $_[0]->{e};
99             }
100             sub doVoronoi {
101             if (@_>1) {$_[0]->{v}=$_[1]?1:0;}
102             return $_[0]->{v};
103             }
104             sub doNeighbors {
105             if (@_>1) {$_[0]->{n}=$_[1]?1:0;}
106             return $_[0]->{n};
107             }
108             sub quiet {
109             if (@_>1) {$_[0]->{Q}=$_[1]?1:0;}
110             return $_[0]->{Q};
111             }
112             sub verbose { # 0 to 3
113             if (@_>1) {$_[0]->{V}=$_[1]?1:0;}
114             return $_[0]->{V};
115             }
116              
117             # everything to add input geometry
118              
119             sub addRegion {
120             my $self = shift;
121             my $poly = shift;
122             my $attribute = @_ ? shift : undef;
123             my $area = @_ ? shift:undef;
124             my $point_inside = @_ ? shift : undef; # not expected, but we'll use it
125              
126             if (@{$poly}==1) {
127             carp "first arg to addRegion should be a polygon, or point";
128             return;
129             }
130             elsif (@{$poly}==2 && !ref($poly->[0])) { # a region identifying point
131             $point_inside = $poly;
132             }
133             else {
134             $self->addPolygon($poly);
135             }
136              
137             my $ray; # return ray used for $point_inside calc for debugging, for now
138              
139             if (!$point_inside) {
140             ($point_inside, $ray) = get_point_in_polygon($poly);
141             }
142             if (defined $point_inside) {
143             push @{$self->{poly}->{regions}}, [ $point_inside, $attribute, ($area && $area > 0) ? $area : -1 ];
144             }
145             return $point_inside, $ray;
146             }
147              
148             sub addHole {
149             my $self = shift;
150             my $poly = shift;
151             my $point_inside = @_ ? shift : undef; # not expected, but we'll use it if available
152              
153             if (@{$poly}==1) {
154             carp "first arg to addHole should be a polygon, or point";
155             return;
156             }
157             elsif (@{$poly}==2 && !ref($poly->[0])) { # it's really the hole identifying point
158             $point_inside = $poly;
159             }
160             else {
161             $self->addPolygon($poly);
162             }
163              
164             my $ray; # return ray used for $point_inside calc for debugging, for now
165              
166             if (!$point_inside) {
167             ($point_inside, $ray) = get_point_in_polygon($poly);
168             }
169             if (defined $point_inside) {
170             push @{$self->{poly}->{holes}}, $point_inside;
171             }
172             return $point_inside, $ray;
173             }
174              
175             sub addPolygon {
176             my $self = shift;
177             my $poly = shift;
178             if (@{$poly} == 1 ) {return $self->addPoints([$poly->[0]]);}
179             push @{$self->{poly}->{polylines}}, ['polygon',$poly];
180             return;
181             }
182              
183             sub addPolyline {
184             my $self = shift;
185             my $poly = shift;
186             if (@{$poly} == 1 ) {return $self->addPoints([$poly->[0]]);}
187             push @{$self->{poly}->{polylines}}, ['polyline',$poly];
188             return;
189             }
190              
191             sub addSegments {
192             my $self = shift;
193             my $segments = shift;
194             push @{$self->{poly}->{segments}}, @{$segments};
195             return;
196             }
197              
198             sub addPoints { # points unaffiliated with PLSG segments
199             my $self = shift;
200             my $points = shift;
201             push @{$self->{poly}->{points}}, @{$points};
202             return;
203             }
204              
205             # compile all the input geometry in to Triangle-format lists
206             # set up option strings
207             # and initialize output lists
208              
209             sub prepPoly {
210             my $self = shift;
211             my $optstr = shift;
212             $self->{optstr} = '';
213             # option string options:
214             # prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh
215             # where __ is an optional number
216             $self->{optstr} .= ''.
217             $optstr.
218             ($self->{q} > -1?'q'.$self->{q}:'').
219             ($self->{a} > -1?'a'.$self->{a}:'').
220             ($self->{e} ?'e':'').
221             ($self->{v} ?'v':'').
222             ($self->{n} ?'n':'').
223             'z'. # always number everything starting with zero
224             ($self->{o2} ?'o2':'').
225             ($self->{Q} ?'Q':'').
226             ($self->{V} > 0?(($self->{V} > 2) ? 'VVV' : ($self->{V} x 'V')) : '')
227             ;
228             my @allpts;
229             my @allsegs;
230              
231             if (@{$self->{poly}->{segments}}) {
232              
233             # The list of segments is the most likely to have duplicate points.
234             # Some algorithms in this space result in lists of segments,
235             # perhaps listing subsegments of intersecting segments,
236             # or representing a boundary or polygon with out-of-order,
237             # non-contiguous segment lists, where shared vertices are
238             # repeated in each segment's record.
239              
240             # addSegments() is meant for that kind of data
241             # and this is where we boil the duplicate points down,
242             # so Triangle doesn't have to.
243              
244             # We look both for duplicate point references and duplicate coordinates.
245             # The coordinate check could collapse points that are topologically
246             # unique, or it could fail to merge points that should be considered
247             # duplicates - but we hope most of the time it does the best thing.
248              
249             foreach my $seg (@{$self->{poly}->{segments}}) {
250             if ( !defined($self->{segptrefs}->{$seg->[0]})
251             && !defined($self->{segptrefs}->{$seg->[0]->[0].','.$seg->[0]->[1]})
252             ) {
253             push @allpts, $seg->[0];
254             $self->{segptrefs}->{$seg->[0]} = $#allpts;
255             $self->{segptrefs}->{$seg->[0]->[0].','.$seg->[0]->[1]} = $#allpts;
256             }
257             push @allsegs, defined($self->{segptrefs}->{$seg->[0]})
258             ? $self->{segptrefs}->{$seg->[0]}
259             : $self->{segptrefs}->{$seg->[0]->[0].','.$seg->[0]->[1]};
260             if ( !defined($self->{segptrefs}->{$seg->[1]})
261             && !defined($self->{segptrefs}->{$seg->[1]->[0].','.$seg->[1]->[1]})
262             ) {
263             push @allpts, $seg->[1];
264             $self->{segptrefs}->{$seg->[1]} = $#allpts;
265             $self->{segptrefs}->{$seg->[1]->[0].','.$seg->[1]->[1]} = $#allpts;
266             }
267             push @allsegs, defined($self->{segptrefs}->{$seg->[1]})
268             ? $self->{segptrefs}->{$seg->[1]}
269             : $self->{segptrefs}->{$seg->[1]->[0].','.$seg->[1]->[1]};
270             }
271             }
272             $self->{segptrefs} = {};
273              
274             if (@{$self->{poly}->{polylines}} || @allsegs) {
275             # doing PSLG - add poly flag to options
276             $self->{optstr} = 'p'.$self->{optstr};
277             #set up points and segments lists for each polygon and polyline
278             foreach my $polycont (@{$self->{poly}->{polylines}}) {
279             my $poly = $polycont->[1];
280             push @allpts, $poly->[0];
281             my $startind=$#allpts;
282             foreach my $thispt (@{$poly}[1..@{$poly}-1]) {
283             push(@allsegs, $#allpts, $#allpts + 1);
284             push(@allpts, $thispt);
285             }
286             if ($polycont->[0] eq 'polygon') { # add segment to close it
287             push(@allsegs, $#allpts, $startind);
288             }
289              
290             }
291             # add segments to C struct
292             my $segs_added_count = $self->in->segmentlist(@allsegs);
293              
294             # Add region mark points and any attributes and area constraints to C struct
295             if (@{$self->{poly}->{regions}}) {
296             my $regions_added_count = $self->in->regionlist(map {grep defined, @{$_->[0]},$_->[1],$_->[2]} @{$self->{poly}->{regions}});
297             }
298             # Add hole mark points to C struct
299             if (@{$self->{poly}->{holes}}) {
300             my $holes_added_count = $self->in->holelist(map {@{$_}} @{$self->{poly}->{holes}});
301             }
302             }
303              
304             # add all points from PSLG, (possibly none)
305             # and then any other points (not associated with segments)
306             # into the C struct
307             my $points_added_count = $self->in->pointlist(map {$_->[0],$_->[1]} (@allpts, @{$self->{poly}->{points}}));
308              
309             # set up attribute array if any points have more than 2 items (the coordinates) in list
310             my $coords_plus_attrs = 2; # 2 for the coords - we'll skip over them when it's time
311             foreach my $point (@allpts, @{$self->{poly}->{points}}) {
312             if ($coords_plus_attrs < @{$point}) {$coords_plus_attrs = @{$point}}
313             }
314             if ($coords_plus_attrs > 2) {
315             # Extend / fill in the attribute lists for any points
316             # that don't have the full set of attributes defined.
317             # Set missing/undefined attributes to zero.
318             foreach my $point (@allpts, @{$self->{poly}->{points}}) {
319             if (@{$point} < $coords_plus_attrs) {
320             foreach (2 .. $coords_plus_attrs - 1) {
321             if (!defined($point->[$_])) {$point->[$_]=0;}
322             }
323             }
324             }
325             # put attributes into C struct
326             $self->in->numberofpointattributes($coords_plus_attrs - 2);
327             my $attributes_added_count = $self->in->pointattributelist(
328             map {@{$_}[2 .. $coords_plus_attrs - 1]} (@allpts, @{$self->{poly}->{points}}));
329             }
330              
331             # discard intermediate data now that it's been loaded into C arrays
332             $self->reset();
333              
334             # set up new triangulateio C structs to receive output
335             $self->{out} = Math::Geometry::Delaunay::Triangulateio->new();
336             $self->{vorout} = Math::Geometry::Delaunay::Triangulateio->new();
337              
338             return;
339             }
340              
341              
342             sub triangulate() {
343             my $self = shift;
344             my $dotopo = defined wantarray && !wantarray; # scalar or array return context
345             my $optstr = @_ ? join('',@_):'';
346             $self->prepPoly($optstr);
347             Math::Geometry::Delaunay::_triangulate($self->{optstr},$self->in->to_ptr,$self->out->to_ptr,$self->vorout->to_ptr);
348             if (defined wantarray) { # else don't do expensive topology stuff if undef/void context
349             if (wantarray && index($self->{optstr},'v') != -1) {
350             return topology($self->out), topology($self->vorout);
351             }
352             return topology($self->out);
353             }
354             return;
355             }
356              
357             # This probably performs well, but it does show up high enough in profiler
358             # reports that it's worth looking into speed-up. Consider something with unpack maybe.
359             sub ltolol {($#_<$_[0])?():map [@_[$_*$_[0]+1..$_*$_[0]+$_[0]]],0..$#_/$_[0]-1}#perl.
360              
361             sub nodes {
362             my $self = shift;
363             my $fromVouty = @_ ? shift : 0;
364             my $triio = $fromVouty ? $self->vorout : $self->out;
365             my $cachetarg = $fromVouty ? 'voutnodes' : 'outnodes';
366             if (@{$self->{poly}->{$cachetarg}} == 0) {
367             my @nodeattributes;
368             if ($triio->numberofpointattributes) {
369             @nodeattributes = ltolol($triio->numberofpointattributes,$triio->pointattributelist);
370             }
371             @{$self->{poly}->{$cachetarg}} = ltolol(2,$triio->pointlist);
372             for (my $i=0;$i<@nodeattributes;$i++) {
373             push @{$self->{poly}->{$cachetarg}->[$i]}, @{$nodeattributes[$i]};
374             }
375             if (!$fromVouty) {
376             my @nodemarkers = $triio->pointmarkerlist;
377             for (my $i=0;$i<@nodemarkers;$i++) {
378             push @{$self->{poly}->{$cachetarg}->[$i]}, $nodemarkers[$i];
379             }
380             }
381             }
382             return $self->{poly}->{$cachetarg};
383             }
384              
385             sub elements {
386             my $self = shift;
387             my $triio = $self->out;
388             my $nodes = $self->nodes;
389             my @outelements;
390             my @triangleattributes;
391             if ($triio->numberoftriangleattributes) {
392             @triangleattributes = ltolol($triio->numberoftriangleattributes,$triio->triangleattributelist);
393             }
394             @outelements = map {[map {$nodes->[$_]} @{$_}]} ltolol($triio->numberofcorners,$triio->trianglelist);
395             for (my $i=0;$i<@triangleattributes;$i++) {
396             push @{$outelements[$i]}, @{$triangleattributes[$i]};
397             }
398             return \@outelements;
399             }
400              
401             sub segments {
402             my $self = shift;
403             my $triio = $self->out;
404             my $nodes = $self->nodes;
405             my @outsegments;
406             my @segmentmarkers = $triio->segmentmarkerlist;
407             @outsegments = map {[$nodes->[$_->[0]],$nodes->[$_->[1]]]} ltolol(2,$triio->segmentlist);
408             for (my $i=0;$i<@segmentmarkers;$i++) {
409             push @{$outsegments[$i]}, $segmentmarkers[$i];
410             }
411             return \@outsegments;
412             }
413            
414             sub edges {
415             my $self = shift;
416             my $fromVouty = @_ ? shift : 0;
417             my $triio = $fromVouty ? $self->vorout : $self->out;
418             my $nodes = $self->nodes($fromVouty);
419             my @outedges;
420             @outedges = map {[map { $_==-1?0:$nodes->[$_]} @{$_}]} ltolol(2,$triio->edgelist);
421             if (!$fromVouty) {
422             my @edgemarkers = $triio->edgemarkerlist;
423             for (my $i=0;$i<@edgemarkers;$i++) {
424             push @{$outedges[$i]}, $edgemarkers[$i];
425             }
426             }
427             return \@outedges;
428             }
429              
430             sub vnodes {return $_[0]->nodes(1);}
431              
432             sub vedges {
433             my $self = shift;
434             my $vedges = $self->edges(1);
435             my $triio = $self->vorout;
436             my @outrays;
437             @outrays = ltolol(2,$triio->normlist);
438             for (my $i=0;$i<@{$vedges};$i++) {
439             # if one end was a ray (missing node ref)
440             # look up the direction vector and use that as missing point
441             # and set third element in edge array to true
442             # as a flag to identify this edge as a ray
443             if (!$vedges->[$i]->[0]) {
444             $vedges->[$i]->[0] = $outrays[$i];
445             $vedges->[$i]->[2] = 1;
446             }
447             elsif (!$vedges->[$i]->[1]) {
448             $vedges->[$i]->[1] = $outrays[$i];
449             $vedges->[$i]->[2] = 2;
450             }
451             else {
452             $vedges->[$i]->[2] = 0;
453             }
454             }
455             return $vedges;
456             }
457              
458             sub topology {
459             my $triio = shift;
460              
461             my $isVoronoi = 0; # we'll detect this when reading edges
462              
463             my $pcnt = 0; # In Voronoi diagram node index corresponds to dual Delaunay element.
464             my @nodes = map {{ point => $_, attributes => [], marker => undef, elements => [], edges => [] , segments => [], index => $pcnt++}} ltolol(2,$triio->pointlist);
465             my $tcnt = 0; # In Delaunay triangulation element index corresponds to dual Voronoi node.
466             my @eles = map {my $ele={ nodes=>[map {$nodes[$_]} @{$_}], marker => undef, edges => [], neighbors => [], attributes => [], index => $tcnt++ }; foreach (@{$ele->{nodes}}) {push(@{$_->{elements}},$ele)};$ele} ltolol($triio->numberofcorners,$triio->trianglelist);
467             my $ecnt = 0; # Corresponding edges in the Delaunay and Voronoi topologies will have the same index.
468             my @edges = map {my $edg={ nodes=>[map {$nodes[$_]} grep {$_>-1} @{$_}], marker => undef, elements => [], vector => undef, index => $ecnt++}; foreach (@{$edg->{nodes}}) {push @{$_->{edges} },$edg};if (!$isVoronoi && ($_->[0] == -1 || $_->[1] == -1)) {$isVoronoi = 1};$edg} ltolol(2,$triio->edgelist);
469             my @segs = map {my $edg={ nodes=>[map {$nodes[$_]} @{$_}], marker => undef, elements => [] }; foreach (@{$edg->{nodes}}) {push @{$_->{segments}},$edg}; $edg} ltolol(2,$triio->segmentlist);
470              
471             my @elementattributes;
472             if ($triio->numberoftriangleattributes) {
473             @elementattributes = ltolol($triio->numberoftriangleattributes,$triio->triangleattributelist);
474             }
475             for (my $i=0;$i<@elementattributes;$i++) {
476             $eles[$i]->{attributes} = $elementattributes[$i];
477             }
478             my @nodeattributes;
479             if ($triio->numberofpointattributes) {
480             @nodeattributes = ltolol($triio->numberofpointattributes,$triio->pointattributelist);
481             }
482             my @nodemarkers = $triio->pointmarkerlist; # always there for pslg, unlike attributes
483             for (my $i=0;$i<@nodemarkers;$i++) {
484             if ($triio->numberofpointattributes) {
485             $nodes[$i]->{attributes} = $nodeattributes[$i];
486             }
487             $nodes[$i]->{marker} = $nodemarkers[$i];
488             }
489             my @edgemarkers = $triio->edgemarkerlist;
490             for (my $i=0;$i<@edgemarkers;$i++) {
491             $edges[$i]->{marker} = $edgemarkers[$i];
492             }
493             my @segmentmarkers = $triio->segmentmarkerlist; # because some can be internal to boundaries
494             for (my $i=0;$i<@segmentmarkers;$i++) {
495             $segs[$i]->{marker} = $segmentmarkers[$i];
496             }
497             my @neighs = ltolol(3,$triio->neighborlist);
498             for (my $i=0;$i<@neighs;$i++) {
499             $eles[$i]->{neighbors} = [map {$eles[$_]} grep {$_ != -1} @{$neighs[$i]}];
500             }
501             if ($isVoronoi) {
502             my @edgevectors = ltolol(2,$triio->normlist); # voronoi ray vectors
503             for (my $i=0;$i<@edgevectors;$i++) {
504             $edges[$i]->{vector} = $edgevectors[$i];
505             }
506             }
507              
508             # cross reference elements and edges
509             if (@edges) { # but only if edges were generated
510             foreach my $ele (@eles) {
511             for (my $i=-1;$i<@{$ele->{nodes}}-1;$i++) {
512             foreach my $edge (@{$ele->{nodes}->[$i]->{edges}}) {
513             if ($ele->{nodes}->[$i+1] == $edge->{nodes}->[0] || $ele->{nodes}->[$i+1] == $edge->{nodes}->[1]) {
514             push @{$ele->{edges}}, $edge;
515             push @{$edge->{elements}}, $ele;
516             last;
517             }
518             }
519             }
520             }
521             }
522              
523             my $ret = {
524             nodes => \@nodes,
525             edges => \@edges,
526             segments => \@segs,
527             elements => \@eles
528             };
529             bless $ret, 'mgd_topo'; # gives the hash a DESTROY method that helps with garbage collection
530             return $ret;
531             }
532              
533             sub get_point_in_polygon {
534             my $poly = shift;
535             my $point_inside;
536              
537             my $bottom_left_index=0;
538             my $maxy = $poly->[$bottom_left_index]->[1];
539             for (my $i=1;$i<@{$poly};$i++) {
540             if ($poly->[$i]->[1] <= $poly->[$bottom_left_index]->[1]) {
541             if ($poly->[$i]->[1] < $poly->[$bottom_left_index]->[1] ||
542             $poly->[$i]->[0] < $poly->[$bottom_left_index]->[0]
543             ) {
544             $bottom_left_index = $i;
545             }
546             }
547             if ($maxy < $poly->[$i]->[1]) { $maxy = $poly->[$i]->[1] }
548             }
549             my $prev_index = $bottom_left_index;
550             my $next_index = -1 * @{$poly} + $bottom_left_index;
551             --$prev_index
552             while ($poly->[$prev_index]->[0] == $poly->[$bottom_left_index]->[0] &&
553             $poly->[$prev_index]->[1] == $poly->[$bottom_left_index]->[1]
554             );
555             ++$next_index
556             while ($poly->[$next_index]->[0] == $poly->[$bottom_left_index]->[0] &&
557             $poly->[$next_index]->[1] == $poly->[$bottom_left_index]->[1]
558             );
559              
560             my @vec1 = ($poly->[$bottom_left_index]->[0] - $poly->[$prev_index]->[0],
561             $poly->[$bottom_left_index]->[1] - $poly->[$prev_index]->[1]);
562             my @vec2 = ($poly->[$next_index]->[0] - $poly->[$bottom_left_index]->[0],
563             $poly->[$next_index]->[1] - $poly->[$bottom_left_index]->[1]);
564             my $orient = (($vec1[0]*$vec2[1] - $vec2[0]*$vec1[1])>=0) ? 1:0;
565            
566             my $angle1 = atan2($poly->[$prev_index]->[1] - $poly->[$bottom_left_index]->[1], $poly->[$prev_index]->[0] - $poly->[$bottom_left_index]->[0]);
567             my $angle2 = atan2($poly->[$next_index]->[1] - $poly->[$bottom_left_index]->[1], $poly->[$next_index]->[0] - $poly->[$bottom_left_index]->[0]);
568             my $angle;
569             if ($orient) {$angle = $angle2 + ($angle1 - $angle2)/2;}
570             else {$angle = $angle1 + ($angle2 - $angle1)/2;}
571             my $cosangle = cos($angle);
572             my $sinangle = sin($angle);
573             my $tanangle = $sinangle/$cosangle;
574              
575             my $adequate_distance = 1.1 * ($maxy - $poly->[$bottom_left_index]->[1]) *
576             ((abs($cosangle) < 0.00000001) ? 1 : 1/$sinangle);
577             my $point_wayout = [$poly->[$bottom_left_index]->[0] + $adequate_distance * $cosangle,
578             $poly->[$bottom_left_index]->[1] + $adequate_distance * $sinangle];
579              
580             my @intersections = sort {$a->[2] <=> $b->[2]} ray_from_index_poly_intersections($bottom_left_index,$point_wayout,$poly,1);
581              
582             if (!@intersections) {
583             print "Warning: Failed to calculate hole or region indicator point.";
584             }
585             elsif ($#intersections % 2 != 0) {
586             print "Warning: Calculated hole or region indicator point is not inside its polygon.";
587             }
588             else {
589             my $closest_intersection = $intersections[0];
590             $point_inside = [($poly->[$bottom_left_index]->[0] + $closest_intersection->[0])/2,
591             ($poly->[$bottom_left_index]->[1] + $closest_intersection->[1])/2];
592             }
593             return $point_inside;
594             }
595              
596             sub ray_from_index_poly_intersections {
597             my $vertind = shift;
598             my $raypt = shift;
599             my $poly = shift;
600             my $doDists = @_ ? shift : 0;
601              
602             my $seg1 = [$poly->[$vertind],$raypt];
603             my $x1= $seg1->[0]->[0];
604             my $y1= $seg1->[0]->[1];
605             my $x2= $seg1->[1]->[0];
606             my $y2= $seg1->[1]->[1];
607             my @lowhix=($x2>$x1)?($x1,$x2):($x2,$x1);
608             my @lowhiy=($y2>$y1)?($y1,$y2):($y2,$y1);
609              
610             my @intersections;
611              
612             for (my $i = -1; $i < $#$poly; $i++) {
613              
614             # skip the segs on either side of the ray base point
615             next if $i == $vertind || ($i + 1) == $vertind || ($vertind == $#$poly && $i == -1);
616              
617             my $seg2 = [$poly->[$i],$poly->[$i+1]];
618             my @segsegret;
619              
620             my $u1= $seg2->[0]->[0];
621             my $v1= $seg2->[0]->[1];
622             my $u2= $seg2->[1]->[0];
623             my $v2= $seg2->[1]->[1];
624              
625             ##to maybe optimize for the case where segments are
626             ##expected NOT to intersect most of the time
627             #my @lowhix=($x2>$x1)?($x1,$x2):($x2,$x1);
628             #my @lowhiu=($u2>$u1)?($u1,$u2):($u2,$u1);
629             #if (
630             # $lowhix[0]>$lowhiu[1]
631             # ||
632             # $lowhix[1]<$lowhiu[0]
633             # ) {
634             # return;
635             # }
636             #my @lowhiy=($y2>$y1)?($y1,$y2):($y2,$y1);
637             #my @lowhiv=($v2>$v1)?($v1,$v2):($v2,$v1);
638             #if (
639             # $lowhiy[0]>$lowhiv[1]
640             # ||
641             # $lowhiy[1]<$lowhiv[0]
642             # ) {
643             # return;
644             # }
645              
646             my $m1 = ($x2 eq $x1)?'Inf':($y2 - $y1)/($x2 - $x1);
647             my $m2 = ($u2 eq $u1)?'Inf':($v2 - $v1)/($u2 - $u1);
648              
649             my $b1;
650             my $b2;
651             my $xi;
652              
653             # Arranged like this to avoid m1-m2 with infinity involved, which works
654             # in other contexts, but can trigger a floating point exception here.
655             # Turns out the exception happens when we let Triangle set the FPU
656             # control word, but not when we use XPFPA.h to take care of that, but
657             # leaving it this as it is in case that doesn't fix that everywhere.
658             my $dm;
659             if ($m1 != 'Inf' && $m2 != 'Inf') {$dm = $m1 - $m2;}
660             elsif ($m1 == 'Inf' && $m2 == 'Inf') {return;}
661             else {$dm='Inf';}
662             if ($m1 == 'Inf' && $m2 != 'Inf') {$xi = $x1;$b2 = $v1 - ($m2 * $u1);}
663             elsif ($m2 == 'Inf' && $m1 != 'Inf') {$xi = $u1;$b1 = $y1 - ($m1 * $x1);}
664             elsif (abs($dm) > 0.000000000001) {
665             $b1 = $y1 - ($m1 * $x1);
666             $b2 = $v1 - ($m2 * $u1);
667             $xi=($b2-$b1)/$dm;
668             }
669             my @lowhiu=($u2>$u1)?($u1,$u2):($u2,$u1);
670             if ($m1 != 'Inf') {
671             if ($m2 == 'Inf' && ($u2<$lowhix[0] || $u2>$lowhix[1]) ) {
672             next;
673             }
674             if (
675             defined $xi &&
676             ($xi < $lowhix[1] || $xi eq $lowhix[1]) &&
677             ($xi > $lowhix[0] || $xi eq $lowhix[0]) &&
678             ($xi < $lowhiu[1] || $xi eq $lowhiu[1]) &&
679             ($xi > $lowhiu[0] || $xi eq $lowhiu[0])
680             ) {
681             my $y=($m1*$xi)+$b1;
682             my @lowhiv=($v2>$v1)?($v1,$v2):($v2,$v1);
683             if ($m2 == 'Inf' &&
684             ($y<$lowhiv[0] || $y>$lowhiv[1])
685             ) {
686             next;
687             }
688             else {
689             push(@intersections,[$xi,$y]);
690             if ($y eq $v1 || $y eq $v2) {
691             $i++; # avoid duplicates at endpoints
692             }
693             }
694             }
695             }
696             elsif ($m2 != 'Inf') {#so $m1 is Inf
697             if (($x1 < $lowhiu[0] || $x1 > $lowhiu[1]) && ! ($x1 eq $lowhiu[0] || $x1 eq $lowhiu[1])) {
698             next;
699             }
700             my @lowhiv=($v2>$v1)?($v1,$v2):($v2,$v1);
701             my $yi = ($m2*$xi)+$b2;
702             if (($yi || $yi eq 0) &&
703             ($yi < $lowhiy[1] || $yi eq $lowhiy[1]) &&
704             ($yi > $lowhiy[0] || $yi eq $lowhiy[0]) &&
705             ($yi < $lowhiv[1] || $yi eq $lowhiv[1]) &&
706             ($yi > $lowhiv[0] || $yi eq $lowhiv[0])
707             ) {
708             push(@intersections,[$xi,$yi]);
709             if ($xi eq $u1 || $xi eq $u2) {
710             $i++; # avoid duplicates at endpoints
711             }
712             }
713             }
714             }
715             if ($doDists) {
716             foreach my $int (@intersections) {
717             push(@{$int}, sqrt( ($int->[0]-$seg1->[0]->[0])**2 + ($int->[1]-$seg1->[0]->[1])**2 ) );
718             }
719             }
720             return @intersections;
721             }
722              
723             # Adjust location of nodes in voronoi diagram so they become
724             # centers of maximum inscribed circles, and store MIC radius in each node.
725             # This is a first step towards a better medial axis approximation.
726             # It straightens out the initial MA approx. derived from the Voronoi diagram.
727              
728             sub mic_adjust {
729             my $topo = shift;
730             my $vtopo = shift;
731              
732             my @new_vnode_points; # voronoi vertices moved to MIC center points
733             my @new_vnode_radii; # will be calculated and added to node data
734             my @new_vnode_tangents; # where the MIC touches the boundary PSLG
735              
736             my $vnc=-1; # will use to look up triangle that corresponds to the voronoi node
737             # $vnc can probably be replaced with $vnode->{index}
738              
739             foreach my $vnode (@{$vtopo->{nodes}}) {
740             $vnc++;
741              
742             push(@new_vnode_points,$vnode->{point});
743             push(@new_vnode_radii,undef);
744             push(@new_vnode_tangents,[]);
745              
746             my $boundary_edge;
747              
748             if (@{$vnode->{edges}} == 3) {
749             my @all_edges = map {$topo->{edges}->[$_->{index}]} @{$vnode->{edges}};
750             my @boundary_edges = grep {$_->{marker}} @all_edges;
751              
752             my @opposite_boundary_edges;
753             my @opposite_boundary_feet;
754              
755             my $branch_node_edge_index1;
756             my $branch_node_edge_index2;
757             my $branch_node_third_tan;
758              
759             # BRANCH NODE
760             # The corresponding Delaunay triangle has no edges on the PSLG boundary,
761             # but touches the boundary at its vertices. This corresponds to a
762             # node with three edges in the medial axis approximation.
763             if (@boundary_edges == 0) {
764             $vnode->{isbranchnode} = 1;
765              
766             # Orient nodes in edges emanating from this branch node so that
767             # the first node is always the branch node.
768             # The notion is that direction is always outward from a branch node.
769             # This will be a problem if two branch nodes link to each other.
770             $_->{nodes} = [$vnode,$_->{nodes}->[0]!=$vnode?$_->{nodes}->[0]:$_->{nodes}->[1]] for @{$vnode->{edges}};
771              
772             # Make sure corners are sorted so they are opposite the Voronoi edges in order.
773             # They may have already been that way, but couldn't determine from Triangle docs.
774             my @corners;
775             for (my $i=0;$i<@{$vnode->{edges}};$i++) {
776             push @corners, +(grep $_ != $topo->{edges}->[$vnode->{edges}->[$i]->{index}]->{nodes}->[0]
777             && $_ != $topo->{edges}->[$vnode->{edges}->[$i]->{index}]->{nodes}->[1],
778             @{$topo->{elements}->[$vnc]->{nodes}}
779             )[0];
780             }
781              
782             my @corner_boundary_edges = grep $_->{marker}, map @{$_->{edges}}, @corners;
783             my @corner_boundary_feet = map {getFoot([$_->{nodes}->[0]->{point},$_->{nodes}->[1]->{point}],$vnode->{point}->[0],$vnode->{point}->[1])} @corner_boundary_edges;
784              
785             # Handle the case where the three corners are probably the three MIC tangents.
786             if (! grep $_, @corner_boundary_feet) {
787             $new_vnode_radii[-1] = dist2d($vnode->{point},$topo->{edges}->[$vnode->{edges}->[0]->{index}]->{nodes}->[0]->{point});
788             $new_vnode_tangents[-1] = [[$corners[2]->{point},
789             $corners[1]->{point}],
790             [$corners[0]->{point},
791             $corners[2]->{point}],
792             [$corners[1]->{point},
793             $corners[0]->{point}]];
794             next;
795             }
796              
797             # Otherwise, set up a case similar to the non-branch node case
798             # by designating feet, a boundary edge, and an opposite edge.
799             my @boundary_feet_edges;
800             for (my $i=0;$i<@corner_boundary_feet;$i+=2) {
801             my $foot;
802             my $edge;
803             # if no foot was found on the two edges meeting at
804             # one of the triangle's vertices, make the vertex a "fake foot"
805             if (!$corner_boundary_feet[$i] && !$corner_boundary_feet[$i+1]) {
806             my $foot = ( $corner_boundary_edges[$i]->{nodes}->[0] == $corner_boundary_edges[$i+1]->{nodes}->[0]
807             || $corner_boundary_edges[$i]->{nodes}->[0] == $corner_boundary_edges[$i+1]->{nodes}->[1] )
808             ? $corner_boundary_edges[$i]->{nodes}->[0]->{point}
809             : $corner_boundary_edges[$i]->{nodes}->[1]->{point};
810             # edge evaluates to "false" to signal this case
811             push @boundary_feet_edges, [$foot,undef,$i/2];
812             }
813             # otherwise keep foot and edge ref
814             else {
815             if ($corner_boundary_feet[$i] ) {push @boundary_feet_edges, [$corner_boundary_feet[$i], $corner_boundary_edges[$i] ,$i/2];}
816             if ($corner_boundary_feet[$i+1]
817             && ( !$corner_boundary_feet[$i] ||
818             # screen out duplicates
819             ( $corner_boundary_feet[$i+1]->[0] ne $corner_boundary_feet[$i]->[0]
820             && $corner_boundary_feet[$i+1]->[1] ne $corner_boundary_feet[$i]->[1]
821             )
822             )
823             ) {push @boundary_feet_edges, [$corner_boundary_feet[$i+1],$corner_boundary_edges[$i+1],$i/2];}
824             }
825             }
826              
827             @boundary_feet_edges = sort {dist2d($a->[0],$vnode->{point}) <=> dist2d($b->[0],$vnode->{point})} @boundary_feet_edges;
828              
829             # closest edge treated as boundary edge, similar to non-branch node handling
830             my $first_with_edge = +(grep {$_->[1]} @boundary_feet_edges)[0];
831              
832             $boundary_edge = {nodes=>[$first_with_edge->[1]->{nodes}->[0],$first_with_edge->[1]->{nodes}->[1]]};
833              
834             # next closest edge treated as opposite boundary edge, similar to non-branch node handling
835             # (that is, next closest that's also not connected to the same
836             # corner as the closest)
837             my $first_not_first_with_edge = +(grep {$_ != $first_with_edge && $_->[2] != $first_with_edge->[2]} @boundary_feet_edges)[0];
838              
839             $opposite_boundary_feet[0] = $first_not_first_with_edge->[0];
840             $opposite_boundary_edges[0] = $first_not_first_with_edge->[1];
841              
842             # Use these later to match up tangent point pairs with Voronoi edges.
843             $branch_node_edge_index1 = $first_with_edge->[2];
844             $branch_node_edge_index2 = $first_not_first_with_edge->[2];
845              
846             my @threst = grep { $_ != $first_with_edge
847             && $_ != $first_not_first_with_edge
848             && $_->[2] != $first_with_edge->[2]
849             && $_->[2] != $first_not_first_with_edge->[2]
850             } @boundary_feet_edges;
851             # Not completely thought out, but results look good so far -
852             # This "tangent point" is not (generally) touched by the approximate
853             # MIC, but it might be worth trying to shift the approx MIC in
854             # a later step to come closer to this, when possible.
855             # Even without that extra shift attempt, it's useful to have
856             # this when reconstructing a polygon from the approximate
857             # medial axis enabled by mic_adjust().
858             $branch_node_third_tan = $threst[0]->[0];
859             }
860              
861             # NON-BRANCH NODE
862             # The corresponding Delaunay triangle has one or two edges
863             # on the PSLG boundary. This corresponds to a node with two edges
864             # in the medial axis approximation.
865             if (@boundary_edges == 1 || @boundary_edges == 2) {
866              
867             $boundary_edge = $boundary_edges[0];
868              
869             my $bef = getFoot([$boundary_edge->{nodes}->[0]->{point},$boundary_edge->{nodes}->[1]->{point}],$vnode->{point}->[0],$vnode->{point}->[1]);
870             my $ber = sqrt(($vnode->{point}->[0]-$bef->[0])**2 + ($vnode->{point}->[1]-$bef->[1])**2);
871              
872             if (@boundary_edges == 2) {
873             # If the other boundary edge is closer, make that the boundary edge.
874             my $obef = getFoot([$boundary_edges[1]->{nodes}->[0]->{point},$boundary_edges[1]->{nodes}->[1]->{point}],$vnode->{point}->[0],$vnode->{point}->[1]);
875             next if (!$obef);
876             my $ober = sqrt(($vnode->{point}->[0]-$obef->[0])**2 + ($vnode->{point}->[1]-$obef->[1])**2);
877             if ($ober < $ber) {
878             $boundary_edge = $boundary_edges[1];
879             $bef=$obef;
880             $ber=$ober;
881             }
882             }
883              
884             my @other_edges = grep $_ != $boundary_edge, map $topo->{edges}->[$_->{index}], @{$vnode->{edges}};
885             my $opposite_node = +(grep {$_ != $boundary_edge->{nodes}->[0] && $_ != $boundary_edge->{nodes}->[1]} @{$other_edges[1]->{nodes}})[0];
886             @opposite_boundary_edges = grep {$_->{marker}} @{$opposite_node->{edges}};
887             @opposite_boundary_feet = map {getFoot([$_->{nodes}->[0]->{point},$_->{nodes}->[1]->{point}],$vnode->{point}->[0],$vnode->{point}->[1])} @opposite_boundary_edges;
888              
889             if (@boundary_edges == 1
890             # || @boundary_edges == 2
891             ) { # should apply to ==2 too, maybe, but should screen out the "opposite" that is also "adjacent" in that case
892             my @adjacent_boundary_edges = grep {$_->{marker} && $_ != $boundary_edge} map @{$_->{edges}}, @{$boundary_edge->{nodes}};
893             my @adjacent_boundary_feet = map {getFoot([$_->{nodes}->[0]->{point},$_->{nodes}->[1]->{point}],$vnode->{point}->[0],$vnode->{point}->[1])} @adjacent_boundary_edges;
894             # if any adjacent bounds closer, replace $boundary_edge with closest
895             for (my $i = 0; $i < @adjacent_boundary_feet; $i++) {
896             next if (!$adjacent_boundary_feet[$i]);
897             my $newfootdist = sqrt(($vnode->{point}->[0] - $adjacent_boundary_feet[$i]->[0])**2 + ($vnode->{point}->[1]-$adjacent_boundary_feet[$i]->[1])**2);
898             if ($newfootdist < $ber) {
899             $boundary_edge = $adjacent_boundary_edges[$i];
900             $bef = $adjacent_boundary_feet[$i];
901             $ber = $newfootdist;
902             }
903             }
904             }
905              
906             if (grep $_, @opposite_boundary_feet) {
907             my @sortind = sort {dist2d($vnode->{point},$opposite_boundary_feet[$a]) <=> dist2d($vnode->{point},$opposite_boundary_feet[$b])} grep {$opposite_boundary_feet[$_]} (0..$#opposite_boundary_feet);
908             @opposite_boundary_feet = map $opposite_boundary_feet[$_], @sortind;
909             @opposite_boundary_edges = map $opposite_boundary_edges[$_], @sortind;
910             @opposite_boundary_feet = ($opposite_boundary_feet[0]);
911             @opposite_boundary_edges = ($opposite_boundary_edges[0]);
912             my $obr = dist2d($opposite_boundary_feet[0],$vnode->{point});
913             if ($ber>$obr) {
914             my $tmp = $opposite_boundary_edges[0];
915             $opposite_boundary_edges[0] = $boundary_edge;
916             $boundary_edge = $tmp;
917             $opposite_boundary_feet[0] = $bef;
918             }
919             }
920              
921             if (!grep $_, @opposite_boundary_feet) {
922             # make the foot just the opposite point
923             @opposite_boundary_feet = ($opposite_node->{point});
924             # make edge eval to false to signal this fake foot case
925             @opposite_boundary_edges = map {0} @opposite_boundary_edges;
926             }
927              
928             }
929              
930             # FIND THE TWO TANGENT POINTS, RADIUS, AND SHIFT POINT TO MIC CENTER
931            
932             # Only one defined unique foot should be in @opposite_boundary_feet
933             # at this point, though that wasn't always the case in the past.
934             # When we're satisfied that it will be the case in the future, we'll
935             # remove the for loop.
936              
937             for (my $i = 0; $i < @opposite_boundary_feet; $i++) {
938             next if !$opposite_boundary_feet[$i];
939             my $foot = $opposite_boundary_feet[$i];
940             my $opp_edge = $opposite_boundary_edges[$i];
941              
942             my $a1;
943             if ($opp_edge) {
944             $a1 = atan2($opp_edge->{nodes}->[0]->{point}->[1] - $opp_edge->{nodes}->[1]->{point}->[1],
945             $opp_edge->{nodes}->[0]->{point}->[0] - $opp_edge->{nodes}->[1]->{point}->[0]);
946             }
947             else { # the "fake foot" case, where opposite edge is just a point
948             $a1 = atan2($foot->[1] - $vnode->{point}->[1],
949             $foot->[0] - $vnode->{point}->[0]);
950             $a1 -= $pi / 2;
951             }
952              
953             my $a2 = atan2($boundary_edge->{nodes}->[1]->{point}->[1] - $boundary_edge->{nodes}->[0]->{point}->[1],
954             $boundary_edge->{nodes}->[1]->{point}->[0] - $boundary_edge->{nodes}->[0]->{point}->[0]);
955              
956             $a1 = angle_reduce_pi($a1);
957              
958             my $amid = ($a1 + $a2) / 2;
959             my $amidnorm = $amid + $pi / 2;
960              
961             my $boundtanpt = line_line_intersection(
962             [ $boundary_edge->{nodes}->[0]->{point}, $boundary_edge->{nodes}->[1]->{point} ],
963             [ $foot, [$foot->[0] + 100 * cos($amidnorm), $foot->[1] + (100 * sin($amidnorm))] ],
964             );
965              
966             if ($boundtanpt) {
967             my $midpt=[($foot->[0]+$boundtanpt->[0])/2,($foot->[1]+$boundtanpt->[1])/2];
968             my $nother_mid_pt=[$midpt->[0]-100*cos($amid),$midpt->[1]-100*sin($amid)];
969             my $center = line_line_intersection([$vnode->{point},$foot],[$midpt,$nother_mid_pt]);
970             if ($center) {
971             $new_vnode_points[-1] = $center;
972             $new_vnode_radii[-1] = dist2d($center,$foot);
973             # assign the three tangent pairs to a branch node
974             if (defined $branch_node_edge_index1) {
975             my $gets_both_found = +(grep $_ != $branch_node_edge_index1 && $_ != $branch_node_edge_index2, (0..2))[0];
976             $new_vnode_tangents[-1]->[( $gets_both_found ) % 3] = [$foot, $boundtanpt];
977             $new_vnode_tangents[-1]->[( $branch_node_edge_index1 ) % 3] = [$branch_node_third_tan, $foot];
978             $new_vnode_tangents[-1]->[( $branch_node_edge_index2 ) % 3] = [$boundtanpt, $branch_node_third_tan];
979             }
980             # assign the two tangent pairs for non-branch nodes
981             else {
982             foreach my $edge (@{$vnode->{edges}}) {
983             next if defined($edge->{vector}) && ($edge->{vector}->[0] != 0 || $edge->{vector}->[1] != 0); # a ray
984             push @{$new_vnode_tangents[-1]}, [$foot, $boundtanpt];
985             }
986             }
987             }
988             }
989             }
990             }
991              
992             if (!defined $new_vnode_radii[-1]) {
993             # This would be a branch node that had feet, but didn't end up
994             # getting adjusted. This is either an okay edge case or a failure.
995             # Let's watch for a while and see if we end up here.
996             $new_vnode_radii[-1] = dist2d($vnode->{point},$topo->{edges}->[$vnode->{edges}->[0]->{index}]->{nodes}->[0]->{point});
997             print "\nFailure to find radius in Math::Geometery::Delaunay::mic_adjust().\nThe developer would like to know if you come across this.\n";
998             }
999              
1000             }
1001              
1002             # Can probably integrate the node point, radius, and tangents assignments
1003             # directly into the above section, and get rid of this temp array stuff.
1004             # On the other hand, if we move this toward using matched index lists (more
1005             # like Triangle's native output) might just export the lists we made, and not
1006             # update the crossref'ed hash structure.
1007             for (my $i = 0; $i < @{$vtopo->{nodes}}; $i++) {
1008             $vtopo->{nodes}->[$i]->{point} = $new_vnode_points[$i];
1009             $vtopo->{nodes}->[$i]->{radius} = $new_vnode_radii[$i];
1010             $vtopo->{nodes}->[$i]->{tangents} = $new_vnode_tangents[$i];
1011             }
1012              
1013             # Fixup left-right ordering of tangent points, now that
1014             # vnodes and their edges should all be nicely centered between them.
1015              
1016             for (my $i = 0; $i < @{$vtopo->{nodes}}; $i++) {
1017             my $vnode = $vtopo->{nodes}->[$i];
1018              
1019             # first pass - could you do iswithin for vnode in its corresponding triangle
1020             # and if not, shift toward next/previous vnode until it is within
1021             # like to intersection of line between tri apex and bound edge
1022             # and line between point and next/prev point... or is there trajectory
1023             # based on two (usually) boundary edges where the tangents are?
1024             #
1025              
1026             my @fixtangents;
1027             my @edges = grep !defined($_->{vector}) || ($_->{vector}->[0] == 0 && $_->{vector}->[1] == 0), @{$vnode->{edges}};
1028              
1029             for (my $j = 0; $j < @edges; $j++) {
1030             my $edge = $edges[$j];
1031              
1032             my $tangents = $vnode->{tangents}->[$j]; # tangent pairs correspond to non-ray edges, in order
1033             my $other_node = ($edge->{nodes}->[0] != $vnode) ? $edge->{nodes}->[0] : $edge->{nodes}->[1];
1034             my $start_node = $vnode;
1035              
1036             # Duplicate node points can happen for common reasons, so go out
1037             # to the next node if we got a duplicate. (Three-in-a-row duplicates
1038             # shouldn't happen.)
1039             if ($vnode->{point}->[0] eq $other_node->{point}->[0] && $vnode->{point}->[1] eq $other_node->{point}->[1]) {
1040             my $other_edge = +(grep $_ != $edge && !(defined $_->{vector} && ($_->{vector}->[0] != 0 || $_->{vector}->[1] != 0)), @{$other_node->{edges}})[0];
1041             if ($other_edge) {
1042             $other_node = $other_edge->{nodes}->[0] != $other_node ? $other_edge->{nodes}->[0] : $other_edge->{nodes}->[1];
1043             }
1044             else {
1045             # Well, then, if this is a non-branch node, set up a test line
1046             # using the previous node, heading into this one.
1047             if (@{$vnode->{tangents}} == 2) {
1048             my $other_edge = $edges[$j == 0 ? 1 : 0];
1049             $start_node = ($other_edge->{nodes}->[0] != $vnode) ? $other_edge->{nodes}->[0] : $other_edge->{nodes}->[1];
1050             $other_node = $vnode;
1051             }
1052             #else { die "Couldn't get other edge in duplicate point case\n\n;" }
1053             }
1054             }
1055              
1056             # Heading out from the start node to the other node, the first tangent associated
1057             # with this edge should be on the left (a counterclockwise turn), and the
1058             # other tangent on the right. We've got Triangle's adaptave robust
1059             # orientation code backing up counterclockwise(), for the really close cases.
1060             my $ccwl = counterclockwise($start_node->{point}, $other_node->{point}, $tangents->[0]);
1061             my $ccwr = counterclockwise($start_node->{point}, $other_node->{point}, $tangents->[1]);
1062              
1063             if ($ccwl < 0) { # first tangent wasn't on the left
1064             if ($ccwr > 0) { # but the second was, so we just need to swap them
1065             @{$tangents} = reverse @{$tangents};
1066             }
1067             else { # But if both were on the right, something's wrong.
1068             # If this is a branch node, and the other two tangent pairs
1069             # don't have the same problem, we can fix up the bad one later.
1070             #push @fixtangents, $j;
1071             $vnode->{fixtangents} = [] if !defined $vnode->{fixtangents};
1072             push @{$vnode->{fixtangents}}, $j;
1073             }
1074             }
1075             elsif ($ccwr > 0) { # Both were on the left. Maybe fix up later.
1076             #push @fixtangents, $j;
1077             $vnode->{fixtangents} = [] if !defined $vnode->{fixtangents};
1078             push @{$vnode->{fixtangents}}, $j;
1079             }
1080             #elsif ($ccwl eq 0) { die "zero" } # Shouldn't happen. Leave it alone if it does.
1081             }
1082             }
1083             for (my $i = 0; $i < @{$vtopo->{nodes}}; $i++) {
1084             if (@{$vtopo->{nodes}->[$i]->{tangents}} == 3 && defined $vtopo->{nodes}->[$i]->{fixtangents} && @{$vtopo->{nodes}->[$i]->{fixtangents}} == 1) {
1085             #print "FIXUP TANGENT PAIR FOR ONE BRANCH AT BRANCH NODE\n";
1086             #$vtopo->{nodes}->[$i]->{color} = 'orange';
1087             # one bad tangent pair out of three for a branch node
1088             # infer that it just needs to be reversed if other two pairs were good
1089              
1090             # debug stuff
1091             #my @edges = grep !defined($_->{vector}) || ($_->{vector}->[0] == 0 && $_->{vector}->[1] == 0), @{$vtopo->{nodes}->[$i]->{edges}};
1092             #my $reconsider_edge = $edges[$vtopo->{nodes}->[$i]->{fixtangents}->[0]];
1093             #my $reconsider_node = $reconsider_edge->{nodes}->[0] != $vtopo->{nodes}->[$i] ? $reconsider_edge->{nodes}->[0]:$reconsider_edge->{nodes}->[1];
1094             #$reconsider_node->{color} = 'green';
1095              
1096             $vtopo->{nodes}->[$i]->{tangents}->[$vtopo->{nodes}->[$i]->{fixtangents}->[0]] =
1097             [$vtopo->{nodes}->[$i]->{tangents}->[($vtopo->{nodes}->[$i]->{fixtangents}->[0] + 1) % 3]->[1],
1098             $vtopo->{nodes}->[$i]->{tangents}->[($vtopo->{nodes}->[$i]->{fixtangents}->[0] - 1) ]->[0]];
1099              
1100             }
1101             #if (@{$vtopo->{nodes}->[$i]->{tangents}} == 3 && defined $vtopo->{nodes}->[$i]->{fixtangents} && @{$vtopo->{nodes}->[$i]->{fixtangents}} != 1) {
1102             #print " MORE NEEDED FIXUP: ",scalar(@{$vtopo->{nodes}->[$i]->{fixtangents}}),"\n";
1103             # if more than one branch had tangents mixed up...
1104             # maybe the other one or two wrong ones are really right?
1105             # ambiguous.
1106             # $vtopo->{nodes}->[$i]->{color} = 'orange';
1107             # }
1108             if ( @{$vtopo->{nodes}->[$i]->{tangents}} == 2
1109             && $vtopo->{nodes}->[$i]->{tangents}->[0]->[0] == $vtopo->{nodes}->[$i]->{tangents}->[1]->[0]
1110             ) {
1111             #print "non-branch maybe needs fixup.\n";
1112             # not sure if this a good way to approach this -
1113             # doesn't definitively capture all that can go
1114             # wrong with non-branch node orientations.
1115             $vtopo->{nodes}->[$i]->{color} = 'green';
1116             if (defined $vtopo->{nodes}->[$i]->{fixtangents} && @{$vtopo->{nodes}->[$i]->{fixtangents}} == 1) {
1117             $vtopo->{nodes}->[$i]->{tangents}->[$vtopo->{nodes}->[$i]->{fixtangents}->[0]] =
1118             [$vtopo->{nodes}->[$i]->{tangents}->[$vtopo->{nodes}->[$i]->{fixtangents}->[0] - 1]->[1],
1119             $vtopo->{nodes}->[$i]->{tangents}->[$vtopo->{nodes}->[$i]->{fixtangents}->[0] - 1]->[0]];
1120             }
1121             }
1122             }
1123             }
1124              
1125             sub getFoot {
1126             my $seg = shift;
1127             my $x = shift;
1128             my $y = shift;
1129             my $foot;
1130             my $m = ($seg->[1]->[0] - $seg->[0]->[0] == 0) ? 'inf' : ($seg->[1]->[1] - $seg->[0]->[1])/($seg->[1]->[0] - $seg->[0]->[0]);
1131             my @sortx = $seg->[0]->[0] < $seg->[1]->[0] ? ($seg->[0]->[0], $seg->[1]->[0]) : ($seg->[1]->[0], $seg->[0]->[0]);
1132             my @sorty = $seg->[0]->[1] < $seg->[1]->[1] ? ($seg->[0]->[1], $seg->[1]->[1]) : ($seg->[1]->[1], $seg->[0]->[1]);
1133             if ($m == 0) {
1134             if ($x >= $sortx[0] && $x <= $sortx[1]) {$foot=[$x, $seg->[0]->[1]];}
1135             }
1136             elsif ($m == 'inf') {
1137             if ($y >= $sorty[0] && $y <= $sorty[1]) {$foot=[$seg->[0]->[0], $y];}
1138             }
1139             else {
1140             my $intersect_x = (($m*$seg->[0]->[0])-($seg->[0]->[1])+((1/$m)*$x)+($y))/($m+(1/$m));
1141             if ($intersect_x >= $sortx[0] && $intersect_x <= $sortx[1]) {
1142             my $intersect_y = -($seg->[0]->[0] - $intersect_x) * $m + $seg->[0]->[1];
1143             if ($intersect_y >= $sorty[0] && $intersect_y <= $sorty[1]) {
1144             $foot = [$intersect_x, $intersect_y];
1145             }
1146             }
1147             }
1148             return $foot;
1149             }
1150              
1151             sub angle_reduce {
1152             my $a=shift;
1153             while($a > $pi / 2) { $a -= $pi; }
1154             while($a <= -$pi / 2) { $a += $pi; }
1155             return $a;
1156             }
1157              
1158             sub angle_reduce_pi {
1159             my $a=shift;
1160             while($a > $pi) { $a -= $pi * 2; }
1161             while($a <= -$pi) { $a += $pi * 2; }
1162             return $a;
1163             }
1164              
1165             sub seg_seg_intersection {
1166             my $seg1 = shift;
1167             my $seg2 = shift;
1168             my $int;
1169              
1170             my $x1= $seg1->[0]->[0]; my $y1= $seg1->[0]->[1];
1171             my $x2= $seg1->[1]->[0]; my $y2= $seg1->[1]->[1];
1172             my $u1= $seg2->[0]->[0]; my $v1= $seg2->[0]->[1];
1173             my $u2= $seg2->[1]->[0]; my $v2= $seg2->[1]->[1];
1174              
1175             my $m1 = ($x2 eq $x1)?'Inf':($y2 - $y1)/($x2 - $x1);
1176             my $m2 = ($u2 eq $u1)?'Inf':($v2 - $v1)/($u2 - $u1);
1177              
1178             my $b1;
1179             my $b2;
1180             my $xi;
1181              
1182             # Arranged like this to avoid m1-m2 with infinity involved, which works
1183             # in other contexts, but can trigger a floating point exception here.
1184             my $dm;
1185             if ($m1 != 'Inf' && $m2 != 'Inf') {$dm = $m1 - $m2;}
1186             elsif ($m1 == 'Inf' && $m2 == 'Inf') {return;}
1187             else {$dm='Inf';}
1188              
1189             if ($m1 == 'Inf' && $m2 != 'Inf') {$xi = $x1;$b2 = $v1 - ($m2 * $u1);}
1190             elsif ($m2 == 'Inf' && $m1 != 'Inf') {$xi = $u1;$b1 = $y1 - ($m1 * $x1);}
1191             elsif (abs($dm) > 0.000000000001) {
1192             $b1 = $y1 - ($m1 * $x1);
1193             $b2 = $v1 - ($m2 * $u1);
1194             $xi=($b2-$b1)/$dm;
1195             }
1196             my @lowhiu=($u2>$u1)?($u1,$u2):($u2,$u1);
1197             if ($m1 != 'Inf') {
1198             my @lowhix=($x2>$x1)?($x1,$x2):($x2,$x1);
1199             if ($m2 == 'Inf' && ($u2<$lowhix[0] || $u2>$lowhix[1]) ) {
1200             return;
1201             }
1202             if (
1203             ($xi || $xi eq 0) &&
1204             ($xi < $lowhix[1] || $xi eq $lowhix[1]) &&
1205             ($xi > $lowhix[0] || $xi eq $lowhix[0]) &&
1206             ($xi < $lowhiu[1] || $xi eq $lowhiu[1]) &&
1207             ($xi > $lowhiu[0] || $xi eq $lowhiu[0])
1208             ) {
1209             my $y=($m1*$xi)+$b1;
1210             my @lowhiv=($v2>$v1)?($v1,$v2):($v2,$v1);
1211             if ($m2 == 'Inf' &&
1212             ($y<$lowhiv[0] || $y>$lowhiv[1])
1213             ) {
1214             return;
1215             }
1216             else {
1217             $int = [$xi,$y];
1218             }
1219             }
1220             }
1221             elsif ($m2 != 'Inf') { #so $m1 is Inf
1222              
1223             if ($x1 < $lowhiu[0] || $x1 > $lowhiu[1] && ! ($x1 eq $lowhiu[0] || $x1 eq $lowhiu[1])) {
1224             return;
1225             }
1226             my @lowhiy=($y2>$y1)?($y1,$y2):($y2,$y1);
1227             my @lowhiv=($v2>$v1)?($v1,$v2):($v2,$v1);
1228             my $yi = ($m2*$xi)+$b2;
1229             if (($yi || $yi eq 0) &&
1230             ($yi < $lowhiy[1] || $yi eq $lowhiy[1]) &&
1231             ($yi > $lowhiy[0] || $yi eq $lowhiy[0]) &&
1232             ($yi < $lowhiv[1] || $yi eq $lowhiv[1]) &&
1233             ($yi > $lowhiv[0] || $yi eq $lowhiv[0])
1234             ) {
1235             $int = [$xi,$yi];
1236             }
1237             }
1238             return $int;
1239             }
1240              
1241             sub line_line_intersection {
1242             my $seg1 = shift;
1243             my $seg2 = shift;
1244             my $int;
1245              
1246             my $x1= $seg1->[0]->[0]; my $y1= $seg1->[0]->[1];
1247             my $x2= $seg1->[1]->[0]; my $y2= $seg1->[1]->[1];
1248             my $u1= $seg2->[0]->[0]; my $v1= $seg2->[0]->[1];
1249             my $u2= $seg2->[1]->[0]; my $v2= $seg2->[1]->[1];
1250              
1251             my $m1 = ($x2 eq $x1)?'Inf':($y2 - $y1)/($x2 - $x1);
1252             my $m2 = ($u2 eq $u1)?'Inf':($v2 - $v1)/($u2 - $u1);
1253              
1254             my $b1;
1255             my $b2;
1256              
1257             my $xi;
1258              
1259             # Arranged like this to avoid m1-m2 with infinity involved, which works
1260             # in other contexts, but can trigger a floating point exception here.
1261             my $dm;
1262             if ($m1 != 'Inf' && $m2 != 'Inf') {$dm = $m1 - $m2;}
1263             elsif ($m1 == 'Inf' && $m2 == 'Inf') {return;}
1264             else {$dm='Inf';}
1265              
1266             if ($m1 == 'Inf' && $m2 != 'Inf') {$xi = $x1;$b2 = $v1 - ($m2 * $u1);}
1267             elsif ($m2 == 'Inf' && $m1 != 'Inf') {$xi = $u1;$b1 = $y1 - ($m1 * $x1);}
1268             elsif (abs($dm) > 0.000000000001) {
1269             $b1 = $y1 - ($m1 * $x1);
1270             $b2 = $v1 - ($m2 * $u1);
1271             $xi=($b2-$b1)/$dm;
1272             }
1273             if ($m1 != 'Inf') {
1274             if (defined $xi) {
1275             my $y=($m1*$xi)+$b1;
1276             $int = [$xi,$y];
1277             }
1278             }
1279             elsif ($m2 != 'Inf') { # so $m1 is Inf
1280             my $yi = ($m2*$xi)+$b2;
1281             if ($yi || $yi eq 0) {
1282             $int = [$xi,$yi];
1283             }
1284             }
1285             return $int;
1286             }
1287              
1288             sub to_svg {
1289             my %spec = @_;
1290             my $triios = [defined($spec{topo}) ? delete $spec{topo} : undef, defined($spec{vtopo}) ? delete $spec{vtopo} : undef];
1291             my $fn = defined($spec{file}) ? delete $spec{file} : '-';
1292             my $dispsize = defined($spec{size}) ? delete $spec{size} : [800, 600];
1293             my $triio = $triios->[0];
1294             my $vorio = @{$triios}?$triios->[1]:undef;
1295             my @edges;
1296             my @segs;
1297             my @pts;
1298             my @vpts;
1299             my @vedges;
1300             my @vrays;
1301             my @circles;
1302             my @elements;
1303             my $maxx;
1304             my $maxy;
1305             my $minx;
1306             my $miny;
1307             if (!$triio) {carp "no geometry provided";return;}
1308             foreach my $key ( keys %spec ) { if (ref($spec{$key}) !~ /ARRAY/) { carp("style config for '$key' should be a reference to an array"); return; } }
1309              
1310             # make copies of points, because we'll be moving and scaling them
1311            
1312             if (ref($triio) =~ /HASH/ && defined $triio->{nodes}) {
1313             push @pts, map [$_,defined $spec{nodes} ? @{$spec{nodes}} : undef], map [@{$_->{point}}], @{$triio->{nodes}};
1314             }
1315             else {
1316             push @pts, map [[@{$_}],defined $spec{nodes} ? @{$spec{nodes}} : undef], ltolol(2,$triio->pointlist);
1317             }
1318             if ($vorio) {
1319             if (ref($vorio) =~ /HASH/ && defined $vorio->{nodes}) {
1320             push @vpts, map [$_,defined $spec{vnodes} ? @{$spec{vnodes}} : undef], map [@{$_->{point}}], @{$vorio->{nodes}};
1321             }
1322             else {
1323             push @vpts, map [[@{$_}],defined $spec{vnodes} ? @{$spec{vnodes}} : undef], ltolol(2,$vorio->pointlist);
1324             }
1325             }
1326              
1327             $maxx = $pts[0]->[0]->[0];
1328             $minx = $pts[0]->[0]->[0];
1329             $maxy = $pts[0]->[0]->[1];
1330             $miny = $pts[0]->[0]->[1];
1331             foreach my $pt (@pts,@vpts) {
1332             if ($maxx < $pt->[0]->[0]) {$maxx = $pt->[0]->[0]}
1333             if ($maxy < $pt->[0]->[1]) {$maxy = $pt->[0]->[1]}
1334             if ($minx > $pt->[0]->[0]) {$minx = $pt->[0]->[0]}
1335             if ($miny > $pt->[0]->[1]) {$miny = $pt->[0]->[1]}
1336             }
1337              
1338             # offset and scale to avoid limitations of svg renderers
1339              
1340             my $dispsizex = '640';
1341             my $dispsizey = '480';
1342            
1343             if (ref($dispsize) =~ /ARRAY/ && @{$dispsize} > 1) {
1344             $dispsizex = $dispsize->[0];
1345             $dispsizey = $dispsize->[1];
1346             }
1347            
1348             # used to scale lines and point circle radii
1349             # so they stay visible in different viewports dimensions
1350             my $scale=(sqrt($dispsizex**2+$dispsizey**2))/sqrt(($maxx-$minx)**2+($maxy-$miny)**2);
1351              
1352             foreach (@pts,@vpts) {
1353             $_->[0]->[0] -= $minx;
1354             $_->[0]->[0] *= $scale;
1355             $_->[0]->[1] -= $miny;
1356             $_->[0]->[1] *= $scale;
1357             }
1358              
1359             my $scaled_maxx = ($maxx - $minx) * $scale;
1360             my $scaled_maxy = ($maxy - $miny) * $scale;
1361             my $scaled_minx = 0;
1362             my $scaled_miny = 0;
1363              
1364             if (ref($triio) =~ /HASH/ && defined $triio->{nodes}) {
1365             if ($spec{edges}) {push @edges, map {[$_,@{$spec{edges}}]} map [$pts[$_->{nodes}->[0]->{index}]->[0],$pts[$_->{nodes}->[1]->{index}]->[0]], @{$triio->{edges}};}
1366             if ($spec{segments}) {push @segs, map {[$_,@{$spec{segments}}]} map [$pts[$_->{nodes}->[0]->{index}]->[0],$pts[$_->{nodes}->[1]->{index}]->[0]], @{$triio->{segments}};}
1367             #ignoring any subparametric points for elements
1368             if ($spec{elements}) {push @elements, map [[map $pts[$_->{index}]->[0], @{$_->{nodes}}[0..2]], (ref($spec{elements}->[0]) =~ /CODE/ ? &{$spec{elements}->[0]}($_) : $spec{elements}->[0]), defined($spec{elements}->[1]) ? $spec{elements}->[1] : ''], @{$triio->{elements}};} #/
1369             }
1370             else {
1371             if ($spec{edges}) {push @edges, map {[[$pts[$_->[0]]->[0],$pts[$_->[1]]->[0]],@{$spec{edges}}]} ltolol(2,$triio->edgelist);}
1372             if ($spec{segments}) {push @segs, map {[[$pts[$_->[0]]->[0],$pts[$_->[1]]->[0]],@{$spec{segments}}]} ltolol(2,$triio->segmentlist);}
1373             #ignoring any subparametric points for elements
1374             if ($spec{elements}) {
1375             push @elements, map {[[$pts[$_->[0]]->[0],$pts[$_->[1]]->[0],$pts[$_->[2]]->[0]],$spec{elements}->[0], defined($spec{elements}->[1]) ? $spec{elements}->[1] : '']} ltolol($triio->numberofcorners,$triio->trianglelist); #/
1376             #read triangle attribute list, so at least those are available for choosing fill color in this case
1377             if (ref($spec{elements}->[0]) =~ /CODE/ && $triio->numberoftriangleattributes > 0 && $triio->numberofregions > 0) {
1378             my @eleattrs = ltolol($triio->numberoftriangleattributes,$triio->triangleattributes);
1379             for (my $i=0;$i<@elements;$i++) {
1380             # topologically-linked triangle not available here
1381             # but we'll fake it at least for the attributes list
1382             # so the color callback can still color according to region id
1383             # or whatever is in the triangle attribute list
1384             $elements[$i]->[1] = &{$spec{elements}->[0]}({attributes=>$eleattrs[$i]});
1385             }
1386             }
1387             }
1388             }
1389              
1390             if ($vorio) {
1391             if (ref($vorio) =~ /HASH/ && defined $vorio->{nodes}) {
1392             # circles only available in this case
1393             if (defined $spec{circles}) {
1394             push @circles, map {
1395             [
1396             [ # this is point and radius: [x,y,r]
1397             @{$vpts[$_->{index}]->[0]},
1398             defined $_->{radius}
1399             ? $_->{radius} * $scale
1400             : dist2d($vpts[$_->{index}]->[0],$edges[$_->{edges}->[0]->{index}]->[0]->[0])
1401             ],
1402             @{$spec{circles}}
1403             ]
1404             } @{$vorio->{nodes}};
1405             }
1406             if ($spec{vedges}) {
1407             @vedges = map [[$vpts[$_->{nodes}->[0]->{index}]->[0],$vpts[$_->{nodes}->[1]->{index}]->[0]],@{$spec{vedges}}], grep $_->{vector}->[0] eq 0 && $_->{vector}->[1] eq 0, @{$vorio->{edges}};
1408             }
1409             if (defined $spec{vrays}) {
1410             @vrays = map [[$vpts[$_->{nodes}->[0]->{index}]->[0],[@{$_->{vector}}]],(defined($spec{vrays}) ? @{$spec{vrays}} : @{$spec{vedges}})], grep $_->{vector}->[0] ne 0 || $_->{vector}->[1] ne 0, @{$vorio->{edges}};
1411             foreach my $ray (@vrays) {
1412             $ray->[0]->[1]->[0] *= $scale;
1413             $ray->[0]->[1]->[1] *= $scale;
1414             $ray->[0]->[1]->[0] += $ray->[0]->[0]->[0];
1415             $ray->[0]->[1]->[1] += $ray->[0]->[0]->[1];
1416             }
1417             }
1418             }
1419             else {
1420             if ($spec{vedges}) {
1421             my @ves =ltolol(2,$vorio->edgelist);
1422             my @vnorms=ltolol(2,$vorio->normlist);
1423             for (my $i=0;$i<@ves;$i++) {
1424             if ($ves[$i]->[0] > -1 && $ves[$i]->[1] > -1) {
1425             push @vedges, [[$vpts[$ves[$i]->[0]]->[0],$vpts[$ves[$i]->[1]]->[0]],@{$spec{vedges}}];
1426             }
1427             elsif (defined $spec{vrays}) {
1428             my $baseidx = ($ves[$i]->[0] != -1)?0:1;
1429             my $basept = $vpts[$ves[$i]->[$baseidx]]->[0];
1430             my $vec = $vnorms[$i];
1431             my $h = sqrt($vec->[0]**2 + $vec->[1]**2);
1432             $vec = [$vec->[0]/$h,$vec->[1]/$h];
1433             push @vrays, [
1434             [$basept,[$basept->[0]+$vec->[0]*$maxx,$basept->[1]+$vec->[1]*$maxx]],
1435             (defined($spec{vrays}) ? @{$spec{vrays}} : @{$spec{vedges}})];
1436             }
1437             }
1438             }
1439             if ($spec{circles}) {
1440             for (my $i=0;$i<@vpts;$i++) {
1441             push @circles, [[@{$vpts[$i]->[0]},dist2d($vpts[$i]->[0],$elements[$i]->[0])],@{$spec{circles}}];
1442             }
1443             }
1444             }
1445             }
1446              
1447             my $margin_x_hi = 5;
1448             my $margin_x_lo = 5;
1449             my $margin_y_hi = 5;
1450             my $margin_y_lo = 5;
1451              
1452             # extend margins to account for the radius of any circles or points
1453             my @round_things = (@circles, (map {[[$_->[0]->[0], $_->[0]->[1], $_->[2]]]} ( ($spec{nodes} ? @pts : () ), ($spec{vnodes} ? @vpts : () ))));
1454             if (scalar(@round_things) > 0) {
1455             my $cir_maxx = $round_things[0]->[0]->[0] + $round_things[0]->[0]->[2];
1456             my $cir_maxy = $round_things[0]->[0]->[1] + $round_things[0]->[0]->[2];
1457             my $cir_minx = $round_things[0]->[0]->[0] - $round_things[0]->[0]->[2];
1458             my $cir_miny = $round_things[0]->[0]->[1] - $round_things[0]->[0]->[2];
1459             foreach my $cir (@round_things) {
1460             if ($cir_maxx < $cir->[0]->[0] + $cir->[0]->[2]) {$cir_maxx = $cir->[0]->[0] + $cir->[0]->[2]}
1461             if ($cir_maxy < $cir->[0]->[1] + $cir->[0]->[2]) {$cir_maxy = $cir->[0]->[1] + $cir->[0]->[2]}
1462             if ($cir_minx > $cir->[0]->[0] - $cir->[0]->[2]) {$cir_minx = $cir->[0]->[0] - $cir->[0]->[2]}
1463             if ($cir_miny > $cir->[0]->[1] - $cir->[0]->[2]) {$cir_miny = $cir->[0]->[1] - $cir->[0]->[2]}
1464             }
1465             if ($cir_maxx-$scaled_maxx > $margin_x_hi) {$margin_x_hi = ($cir_maxx - $scaled_maxx) + 5;}
1466             if ($scaled_minx-$cir_minx > $margin_x_lo) {$margin_x_lo = ($scaled_minx - $cir_minx) + 5;}
1467             if ($cir_maxy-$scaled_maxy > $margin_y_hi) {$margin_y_hi = ($cir_maxy - $scaled_maxy) + 5;}
1468             if ($scaled_miny-$cir_miny > $margin_y_lo) {$margin_y_lo = ($scaled_miny - $cir_miny) + 5;}
1469             }
1470              
1471             open(SVGO,'>'.$fn);
1472             print SVGO sprintf <<"EOS", $dispsizex, $dispsizey, -$margin_x_lo, -$margin_y_hi, $scaled_maxx + ($margin_x_lo + $margin_x_hi), $scaled_maxy + ($margin_y_lo + $margin_y_hi), $scaled_maxy;
1473            
1474            
1475            
1476            
1477             EOS
1478              
1479             if ($spec{elements}) {print SVGO "\n\n";}
1480             foreach my $ele (@elements) {
1481             print SVGO '',"\n"
1482             }
1483             if ($spec{edges}) {print SVGO "\n\n";}
1484             foreach my $edge (@edges) {
1485             print SVGO '',"\n"
1486             }
1487             if ($spec{segments}) {print SVGO "\n\n";}
1488             foreach my $edge (@segs) {
1489             print SVGO '',"\n"
1490             }
1491             if ($spec{vedges} && $vorio) {print SVGO "\n\n";}
1492             foreach my $edge (@vedges) {
1493             print SVGO '',"\n"
1494             }
1495             if ($spec{vrays} && $vorio) {print SVGO "\n\n";}
1496             foreach my $edge (@vrays) {
1497             print SVGO '',"\n"
1498             }
1499             if ($spec{nodes}) {
1500             print SVGO "\n\n";
1501             foreach my $pt (@pts) {
1502             print SVGO '',"\n"
1503             }
1504             }
1505             if ($spec{vnodes} && $vorio) {
1506             print SVGO "\n\n";
1507             foreach my $pt (@vpts) {
1508             print SVGO '',"\n"
1509             }
1510             }
1511             if ($spec{circles}) {print SVGO "\n\n";}
1512             foreach my $circle (@circles) {
1513             print SVGO '',"\n"
1514             }
1515              
1516             if ($spec{raw}) {
1517             print SVGO "\n\n";
1518             print SVGO join("\n",map { s/ (c?x[12]?)="([0-9\.eE\-]+)"/' '.$1.'="'.(($2-$minx)*$scale).'"'/ge;
1519             s/ (c?y[12]?)="([0-9\.eE\-]+)"/' '.$1.'="'.(($2-$miny)*$scale).'"'/ge;
1520             $_;
1521             } @{$spec{raw}});
1522             }
1523             print SVGO "\n";
1524             close(SVGO);
1525             }
1526              
1527             sub dist2d {sqrt(($_[0]->[0]-$_[1]->[0])**2+($_[0]->[1]-$_[1]->[1])**2)}
1528              
1529             sub counterclockwise {
1530             my ($pa, $pb, $pc) = @_;
1531             return _counterclockwise($pa->[0],$pa->[1],$pb->[0],$pb->[1],$pc->[0],$pc->[1]);
1532             }
1533              
1534             package mgd_topo;
1535              
1536             sub DESTROY {
1537             # circular references in the topology data seem to thwart garbage collection
1538             # so we'll undo some of the cross references to help with that
1539             my $self = shift;
1540             if (exists $self->{elements}) { undef $_->{nodes} for @{$self->{elements}};}
1541             if (exists $self->{edges}) { undef $_->{nodes} for @{$self->{edges}}; }
1542             if (exists $self->{segments}) { undef $_->{nodes} for @{$self->{segments}};}
1543             }
1544              
1545             =head1 NAME
1546              
1547             Math::Geometry::Delaunay - Quality Mesh Generator and Delaunay Triangulator
1548              
1549             =head1 VERSION
1550              
1551             Version 0.19
1552              
1553             =cut
1554              
1555             =head1 SYNOPSIS
1556              
1557             =for html
1558            
1559            
1560            
1561            
1562            
1563            
1564            
1565            
1566            
1567            
input vertices
1568            
1569            
1570            
1571            
1572            
1573            
1574            
1575            
1576            
1577            
1578            
1579            
1580            
1581            
1582            
1583            
1584            
1585            
1586            
1587            
Delaunay triangulation
1588            
1589            
1590            
1591            
1592            
1593            
1594            
1595            
1596            
1597            
1598            
1599            
1600            
1601            
1602            
1603            
1604            
1605            
1606            
Voronoi diagram
1607            
1608              
1609             use Math::Geometry::Delaunay qw(TRI_CCDT);
1610              
1611             # generate Delaunay triangulation
1612             # and Voronoi diagram for a point set
1613              
1614             my $point_set = [ [1,1], [7,1], [7,3],
1615             [3,3], [3,5], [1,5] ];
1616              
1617             my $tri = new Math::Geometry::Delaunay();
1618             $tri->addPoints($point_set);
1619             $tri->doEdges(1);
1620             $tri->doVoronoi(1);
1621            
1622             # called in void context
1623             $tri->triangulate();
1624             # populates the following lists
1625              
1626             $tri->elements(); # triangles
1627             $tri->nodes(); # points
1628             $tri->edges(); # triangle edges
1629             $tri->vnodes(); # Voronoi diagram points
1630             $tri->vedges(); # Voronoi edges and rays
1631              
1632              
1633             =for html
1634            
1635            
1636            
1637            
1638            
1639            
1640            
1641            
1642            
1643            
1644            
1645            
1646            
1647            
1648            
1649            
input PSLG
1650            
1651            
1652            
1653            
1654            
1655            
1656            
1657            
1658            
1659            
1660            
1661            
1662            
1663            
1664            
1665            
1666            
1667            
1668            
1669            
1670            
1671            
1672            
1673            
1674            
1675            
1676            
1677            
1678            
1679            
1680            
1681            
1682            
1683            
1684            
1685            
1686            
1687            
1688            
1689            
1690            
1691            
1692            
1693            
1694            
1695            
1696            
1697            
1698            
1699            
1700            
1701            
1702            
1703            
1704            
1705            
1706            
1707            
1708            
1709            
1710            
1711            
1712            
1713            
1714            
1715            
1716            
1717            
1718            
1719            
1720            
1721            
1722            
1723            
1724            
1725            
1726            
1727            
1728            
1729            
1730            
1731            
1732            
1733            
1734            
1735            
1736            
1737            
1738            
1739            
1740            
1741            
1742            
1743            
1744            
1745            
1746            
1747            
1748            
1749            
1750            
1751            
1752            
1753            
1754            
1755            
1756            
1757            
1758            
1759            
1760            
1761            
1762            
1763            
1764            
1765            
1766            
1767            
1768            
1769            
1770            
1771            
1772            
1773            
1774            
1775            
1776            
1777            
1778            
output mesh
1779            
1780            
1781            
1782            
1783            
1784            
1785            
1786            
1787            
1788            
1789            
1790            
1791            
1792            
1793            
1794            
1795            
1796            
1797            
1798            
1799            
1800            
1801            
1802            
1803            
1804            
1805            
1806            
1807            
1808            
1809            
1810            
1811            
1812            
1813            
1814            
1815            
1816            
1817            
1818            
1819            
1820            
1821            
something interesting
extracted from topology
1822            
1823              
1824             # quality mesh of a planar straight line graph
1825             # with cross-referenced topological output
1826              
1827             my $tri = new Math::Geometry::Delaunay();
1828             $tri->addPolygon($point_set);
1829             $tri->minimum_angle(23);
1830             $tri->doEdges(1);
1831              
1832             # called in scalar context
1833             my $mesh_topology = $tri->triangulate(TRI_CCDT);
1834             # returns cross-referenced topology
1835              
1836             # make two lists of triangles that touch boundary segments
1837              
1838             my @tris_with_boundary_segment;
1839             my @tris_with_boundary_point;
1840              
1841             foreach my $triangle (@{$mesh_topology->{elements}}) {
1842             my $nodes_on_boundary_count = (
1843             grep $_->{marker} == 1,
1844             @{$triangle->{nodes}}
1845             );
1846             if ($nodes_on_boundary_count == 2) {
1847             push @tris_with_boundary_segment, $triangle;
1848             }
1849             elsif ($nodes_on_boundary_count == 1) {
1850             push @tris_with_boundary_point, $triangle;
1851             }
1852             }
1853            
1854              
1855             =for html
1856              
1857             =head1 DESCRIPTION
1858              
1859             This is a Perl interface to the Jonathan Shewchuk's Triangle library.
1860              
1861             "Triangle generates exact Delaunay triangulations, constrained Delaunay
1862             triangulations, conforming Delaunay triangulations, Voronoi diagrams, and
1863             high-quality triangular meshes. The latter can be generated with no small or
1864             large angles, and are thus suitable for finite element analysis."
1865             -- from L
1866              
1867             =head1 EXPORTS
1868              
1869             Triangle has several option switches that can be used in different combinations
1870             to choose a class of triangulation and then configure options within that class.
1871             To clarify the composition of option strings, or just to give you a head start,
1872             a few constants are supplied to configure different classes of mesh output.
1873              
1874             TRI_CONSTRAINED = 'Y' for "Constrained Delaunay"
1875             TRI_CONFORMING = 'Dq0' for "Conforming Delaunay"
1876             TRI_CCDT = 'q' for "Constrained Conforming Delaunay"
1877             TRI_VORONOI = 'v' to generate the Voronoi diagram
1878              
1879             For an illustration of these terms, see:
1880             L
1881              
1882             =head1 CONSTRUCTOR
1883              
1884             =head2 new
1885              
1886             The constructor returns a Math::Geometry::Delaunay object.
1887              
1888             my $tri = Math::Geometry::Delaunay->new();
1889              
1890             =head1 MESH GENERATION
1891              
1892             =head2 triangulate
1893              
1894             Run the triangulation with specified options, and either populate the object's
1895             output lists, or return a hash reference giving access to a cross-referenced
1896             representation of the mesh topology.
1897              
1898             Common options can be set prior to calling C. The full range of
1899             Triangle's options can also be passed to C as a string, or list
1900             of strings. For example:
1901              
1902             my $tri = Math::Geometry::Delaunay->new('pzq0eQ');
1903              
1904             my $tri = Math::Geometry::Delaunay->new(TRI_CCDT, 'q15', 'a3.5');
1905              
1906             Triangle's command line switches are documented here:
1907             L
1908              
1909             =head3 list output
1910              
1911             After triangulate is invoked in void context, the output mesh data can be
1912             retrieved from the following methods, all of which return a reference to an
1913             array.
1914              
1915             $tri->triangulate(); # void context - no return value requested
1916             # output lists now available
1917             $points = $tri->nodes(); # array of vertices
1918             $tris = $tri->elements(); # array of triangles
1919             $edges = $tri->edges(); # all the triangle edges
1920             $segs = $tri->segments(); # the PSLG segments
1921             $vpoints = $tri->vnodes(); # points in the voronoi diagram
1922             $vedges = $tri->vedges(); # edges in the voronoi diagram
1923              
1924             Data may not be available for all lists, depending on which option switches were
1925             used. By default, nodes and elements are generated, while edges are not.
1926              
1927             The members of the lists returned have these formats:
1928              
1929             nodes: [x, y, < zero or more attributes >, < boundary marker >]
1930              
1931             elements: [[x0, y0], [x1, y1], [x2, y2],
1932             < another three vertices, if "o2" switch used >,
1933             < zero or more attributes >
1934             ]
1935             edges: [[x0, y0], [x1, y1], < boundary marker >]
1936              
1937             segments: [[x0, y0], [x1, y1], < boundary marker >]
1938              
1939             vnodes: [x, y, < zero or more attributes >]
1940              
1941             vedges: [< vertex or vector >, < vertex or vector >, < ray flag >]
1942              
1943             Boundary markers are 1 or 0. An edge or segment with only one end on a boundary
1944             has boundary marker 0.
1945              
1946             The ray flag is 0 if the edge is not a ray, or 1 or 2, to indicate
1947             which vertex is actually a unit vector indicating the direction of the ray.
1948              
1949             Import of the mesh data from the C data structures will be deferred until
1950             actually requested from the list fetching methods above. For speed and
1951             lower memory footprint, access only what you need, and consider suppressing
1952             output you don't need with option switches.
1953              
1954             =head3 topological output
1955              
1956             When triangulate is invoked in scalar or array context, it returns a hash ref
1957             containing the cross-referenced nodes, elements, edges, and PSLG segments of the
1958             triangulation. In array context, with the "v" switch enabled, the Voronoi
1959             topology is the second item returned.
1960              
1961             my $topology = $tri->triangulate();
1962              
1963             $topology now looks like this:
1964            
1965             {
1966             nodes => [
1967             { # a node
1968             point => [x0, x1],
1969             edges => [edgeref, ...],
1970             segments => [edgeref, ...], # a subset of edges
1971             elements => [elementref, ...],
1972             marker => 1 or 0 or undefined, # boundary marker
1973             attributes => [attr0, ...]
1974             },
1975             ... more nodes like that
1976            
1977             ],
1978             elements => [
1979             { # a triangle
1980             nodes => [noderef0, noderef1, noderef2],
1981             edges => [edgeref0, edgeref1, edgeref2],
1982             neighbors => [neighref0, neighref1, neighref2],
1983             attributes => [attrib0, ...]
1984             },
1985             ... more triangles like that
1986             ],
1987             edges => [
1988             {
1989             nodes => [noderef0, noderef1], # only one for a ray
1990             elements => [elemref0, elemref1], # one if on boundary
1991             vector => undefined or [x, y], # ray direction
1992             marker => 1 or 0 or undefined, # boundary marker
1993             index => # edge's index in edge list
1994             },
1995             ... more edges like that
1996            
1997             ],
1998             segments => [
1999             {
2000             nodes => [noderef0, noderef1],
2001             elements => [elemref0, elemref1], # one if on boundary
2002             marker => 1 or 0 or undefined # boundary marker
2003             },
2004             ... more segments
2005             ]
2006             }
2007              
2008             =head3 cross-referencing Delaunay and Voronoi
2009              
2010             Corresponding Delaunay triangles and Voronoi nodes have the same index number
2011             in their respective lists.
2012              
2013             In the topological output, any element in a triangulation has a record of its
2014             own index number that can by used to look up the corresponding node in the
2015             Voronoi diagram topology, or vice versa, like so:
2016              
2017             ($topo, $voronoi_topo) = $tri->triangulate('v');
2018              
2019             # get a triangle reference where the index is not obvious
2020            
2021             $element = $topo->{nodes}->[-1]->{elements}->[-1];
2022            
2023             # this gets a reference to the corresponding node in the Voronoi diagram
2024            
2025             $voronoi_node = $voronoi_topo->{nodes}->[$element->{index}];
2026              
2027              
2028             Corresponding edges in the Delaunay and Voronoi outputs have the same index
2029             number in their respective edge lists.
2030              
2031             In the topological output, any edge in a triangulation has a record of its own
2032             index number that can by used to look up the corresponding edge in the Voronoi
2033             diagram topology, or vice versa, like so:
2034              
2035             ($topo, $voronoi_topo) = $tri->triangulate('ev');
2036            
2037             # get an edge reference where it's not obvious what the edge's index is
2038            
2039             $delaunay_edge = $topo->{nodes}->[-1]->{edges}->[-1];
2040            
2041             # this gets a reference to the corresponding edge in the Voronoi diagram
2042            
2043             $voronoi_edge = $voronoi_topo->{edges}->[$delaunay_edge->{index}];
2044              
2045             =head1 METHODS TO SET SOME Triangle OPTIONS
2046              
2047             =head2 area_constraint
2048              
2049             Corresponds to the "a" switch.
2050              
2051             With one argument, sets the maximum triangle area constraint for the
2052             triangulation. Returns the value supplied. With no argument, returns the
2053             current area constraint.
2054              
2055             Passing -1 to C will disable the global area constraint.
2056              
2057             =head2 minimum_angle
2058              
2059             Corresponds to the "q" switch.
2060              
2061             With one argument, sets the minimum angle allowed for triangles added in the
2062             triangulation. Returns the value supplied. With no argument, returns the
2063             current minimum angle constraint.
2064              
2065             Passing -1 to C will cause the "q" switch to be omitted from
2066             the option string.
2067              
2068             =head2 doEdges, doVoronoi, doNeighbors
2069              
2070             These methods simply add or remove the corresponding letters from the
2071             option string. Pass in a true or false value to enable or disable.
2072             Invoke with no argument to read the current state.
2073              
2074             =head2 quiet, verbose
2075              
2076             Triangle prints a basic summary of the meshing operation to STDOUT unless
2077             the "Q" switch is present. This module includes the "Q" switch by default, but
2078             you can override this by passing a false value to C.
2079              
2080             If you would like to see even more output regarding the triangulation process,
2081             there are are three levels of verbosity configurable with repeated "V"
2082             switches. Passing a number from 1 to 3 to the C method will enable
2083             the corresponding level of verbosity.
2084              
2085             =head1 METHODS TO ADD VERTICES AND SEGMENTS
2086              
2087             =head2 addVertices, addPoints
2088              
2089             Takes a reference to an array of vertices, each vertex itself an reference to
2090             an array containing two coordinates and zero or more attributes. Attributes
2091             are floating point numbers.
2092            
2093             # vertex format
2094             # [x, y, < zero or more attributes as floating point numbers >]
2095              
2096             $tri->addPoints([[$x0, $y0], [$x1, $y1], ... ]);
2097              
2098             Use addVertices to add vertices that are not part of a PSLG.
2099             Use addPoints to add points that are not part of a polygon or polyline.
2100             In other words, they do the same thing.
2101              
2102             =head2 addSegments
2103              
2104             Takes a reference to an array of segments.
2105              
2106             # segment format
2107             # [[$x0, $y0], [$x1, $y1]]
2108              
2109             $tri->addSegments([ $segment0, $segment1, ... ]);
2110              
2111             If your segments are contiguous, it's better to use addPolyline, or addPolygon.
2112              
2113             This method is provided because some point and polygon processing algorithms
2114             result in segments that represent polygons, but list the segments in a
2115             non-contiguous order, and have shared vertices repeated in each segment's record.
2116              
2117             The segments added with this method will be checked for duplicate vertices, and
2118             references to these will be merged.
2119              
2120             Triangle can handle duplicate vertices, but we would rather not feed them in on
2121             purpose.
2122              
2123             =head2 addPolyline
2124              
2125             Takes a reference to an array of vertices describing a curve.
2126             Creates PSLG segments for each pair of adjacent vertices. Adds the
2127             new segments and vertices to the triangulation input.
2128              
2129             $tri->addPolyline([$vertex0, $vertex1, $vertex2, ...]);
2130              
2131             =head2 addPolygon
2132              
2133             Takes a reference to an array of vertices describing a polygon.
2134             Creates PSLG segments for each pair of adjacent vertices
2135             and creates and additional segment linking the last vertex to
2136             the first,to close the polygon. Adds the new segments and vertices
2137             to the triangulation input.
2138              
2139             $tri->addPolygon([$vertex0, $vertex1, $vertex2, ...]);
2140              
2141             =head2 addHole
2142              
2143             Like addPolygon, but describing a hole or concavity - an area of the output mesh
2144             that should not be triangulated.
2145              
2146             There are two ways to specify a hole. Either provide a list of vertices, like
2147             for addPolygon, or provide a single vertex that lies inside of a polygon, to
2148             identify that polygon as a hole.
2149              
2150             # first way
2151             $tri->addHole([$vertex0, $vertex1, $vertex2, ...]);
2152              
2153             # second way
2154             $tri->addPolygon( [ [0,0], [1,0], [1,1], [0,1] ] );
2155             $tri->addHole( [0.5,0.5] );
2156              
2157             Hole marker points can also be used, in combination with the "c" option, to
2158             cause or preserve concavities in a boundary when Triangle would otherwise
2159             enclose a PSLG in a convex hull.
2160              
2161             =head2 addRegion
2162              
2163             Takes a polygon describing a region, and an attribute or area constraint. With
2164             both the "A" and "a" switches in effect, three arguments allow you to specify
2165             both an attribute and an optional area constraint.
2166              
2167             The first argument may alternately be a single vertex that lies inside of
2168             another polygon, to identify that polygon as a region.
2169              
2170             To be used in conjunction with the "A" and "a" switches.
2171              
2172             # with the "A" switch
2173             $tri->addRegion(\@polygon, < attribute > );
2174            
2175             # with the "a" switch
2176             $tri->addRegion(\@polygon, < area constraint > );
2177              
2178             # with both "Aa"
2179             $tri->addRegion(\@polygon, < attribute >, < area constraint > );
2180              
2181             If the "A" switch is used, each triangle generated within the bounds of a region
2182             will have that region's attribute added to the end of the triangle's
2183             attributes list, while each triangle not within a region will have a "0" added
2184             to the end of its attribute list.
2185              
2186             If the "a" switch is used without a number following, each triangle generated
2187             within the bounds of a region will be subject to that region's area
2188             constraint.
2189              
2190             If the "A" or "a" switches are not in effect, addRegion has the same effect as
2191             addPolygon.
2192              
2193             =head1 METHODS TO ACCESS OUTPUT LISTS
2194              
2195             The following methods retrieve the output lists after the triangulate method has
2196             been invoked in void context.
2197              
2198             Triangle's output data is not imported from C to Perl until one of these methods
2199             is invoked, and then only what's needed to construct the list requested. So
2200             there may be a speed or memory advantage to accessing the output in this way -
2201             only what you need, when you need it.
2202              
2203             The methods prefixed with "v" access the Voronoi diagram nodes and edges, if one
2204             was generated.
2205              
2206             =head2 nodes
2207              
2208             Returns a reference to a list of nodes (vertices or points).
2209              
2210             my $pointlist = $tri->nodes(); # retrieve nodes/vertices/points
2211            
2212             The nodes in the list have this structure:
2213              
2214             [x, y, < zero or more attributes >, < boundary marker >]
2215              
2216             =head2 elements
2217              
2218             Returns a reference to a list of elements.
2219              
2220             $triangles = $tri->elements(); # retrieve triangle list
2221              
2222             The elements in the list have this structure:
2223              
2224             [[x0, y0], [x1, y1], [x2, y2],
2225             < another three vertices, if "o2" switch used >
2226             < zero or more attributes >
2227             ]
2228              
2229             =head2 segments
2230              
2231             Returns a reference to a list of segments.
2232              
2233             $segs = $tri->segments(); # retrieve the PSLG segments
2234              
2235             The segments in the list have this structure:
2236              
2237             [[x0, y0], [x1, y1], < boundary marker >]
2238              
2239             =head2 edges
2240              
2241             Returns a reference to a list of edges.
2242              
2243             $edges = $tri->edges(); # retrieve all the triangle edges
2244              
2245             The edges in the list have this structure:
2246              
2247             [[x0, y0], [x1, y1], < boundary marker >]
2248              
2249             Note that the edge list is not produced by default. Request that it be generated
2250             by invoking C, or passing the 'e' switch to C.
2251              
2252             =head2 vnodes
2253              
2254             Returns a reference to a list of nodes in the Voronoi diagram.
2255              
2256             $vpointlist = $tri->vnodes(); # retrieve Voronoi vertices
2257              
2258             The Voronoi diagram nodes in the list have this structure:
2259              
2260             [x, y, < zero or more attributes >]
2261              
2262             =head2 vedges
2263              
2264             Returns a reference to a list of edges in the Voronoi diagram. Some of these
2265             edges are actually rays.
2266              
2267             $vedges = $tri->vedges(); # retrieve Voronoi diagram edges and rays
2268              
2269             The Voronoi diagram edges in the list have this structure:
2270              
2271             [< vertex or vector >, < vertex or vector >, < ray flag >]
2272              
2273             If the edge is a true edge, the ray flag will be 0.
2274             If the edge is actually a ray, the ray flag will either be 1 or 2,
2275             to indicate whether the the first, or second vertex should be interpreted as
2276             a direction vector for the ray.
2277              
2278             =head1 UTILITY FUNCTIONS
2279              
2280             =head2 to_svg
2281              
2282             This function is meant as a development and debugging aid, to "dump" the
2283             geometric data structures specific to this package to a graphical
2284             representation. Takes key-value pairs to specify topology hashes, output file,
2285             image dimensions, and styles for the elements in the various output lists.
2286              
2287             The topology hash input for the C or C keys is just the hash
2288             returned by C. The value for the C key is a file name string.
2289             Omit C to print to STDOUT. For C, provide and array ref with width
2290             and height, in pixels. For output list styles, keys correspond to the output
2291             list names, and values consist of references to arrays containing style
2292             configurations, as demonstrated below.
2293              
2294             Only geometry that has a style configuration will be displayed. The following
2295             example includes everything. To display a subset, just omit any of the style
2296             configuration key-value pairs.
2297              
2298             ($topo, $vtopo) = $tri->triangulate('ve');
2299              
2300             to_svg( topo => $topo,
2301             vtopo => $vtopo,
2302            
2303             file => "enchilada.svg", # omit for STDOUT
2304             size => [800, 600], # width, height in pixels
2305            
2306             # line width or optional
2307             # svg color point radius extra CSS
2308            
2309             nodes => ['black' , 0.3],
2310             edges => ['#CCCCCC', 0.7],
2311             segments => ['blue' , 0.9, 'stroke-dasharray:1 1;'],
2312             elements => ['pink'] , # string or callback; see below
2313              
2314             # these require Voronoi input (vtopo)
2315              
2316             vnodes => ['purple' , 0.3],
2317             vedges => ['#FF0000', 0.7],
2318             vrays => ['purple' , 0.6],
2319             circles => ['orange' , 0.6],
2320            
2321             );
2322              
2323             Note that for display purposes C does not include the infinite rays in
2324             the Voronoi diagram. To see the complete Voronoi diagram, including segments
2325             representing the infinite rays, you should include style configuration for the
2326             C key, as in the example above.
2327              
2328             Elements (triangles) only need one style config entry, for color. (An optional
2329             second entry would be a string for additional CSS.) In this case,
2330             the first entry can also be a reference to a callback function. A reference to
2331             the triangle being processed for display will be passed to the callback
2332             function. Therefore the callback function can determine a color based on any
2333             features or relationships of that triangle.
2334              
2335             Typically you might color each triangle according to the region it's in, by
2336             using Triangle's 'A' switch, and then reading the region attribute from the
2337             last item in the triangle's attribute list.
2338              
2339             my $region_colors_callback = sub {
2340             my $tri_ref = shift;
2341             return ('gray','blue','green')[$tri_ref->{attributes}->[-1]];
2342             };
2343              
2344             But any other data accessible through the triangle reference can be used to
2345             calculate a color. For instance, the triangle's three nodes can carry any
2346             number of attributes, which are interpolated during mesh generation. You
2347             might shade each triangle according to the average of a node attribute.
2348              
2349             my $tri_nodes_average_callback = sub {
2350             my $tri_ref = shift;
2351             my $sum = 0;
2352             # calculate average of the eighth attribute in all nodes
2353             foreach my $node (@{$tri_ref->{nodes}}) {
2354             $sum += $node->{attributes}->[7];
2355             }
2356             return &attrib_val_to_grayscale_hexcode( $sum / 3 );
2357             };
2358              
2359             =head2 mic_adjust
2360              
2361             =for html
2362            
2363            
2364            
2365            
2366            
2367            
2368            
2369            
2370            
2371            
2372            
2373            
2374            
2375            
2376            
2377            
2378            
2379            
2380            
2381            
2382            
2383            
2384            
2385            
2386            
2387            
2388            
2389            
2390            
2391            
2392            
2393            
2394            
2395            
2396            
2397            
2398            
2399            
2400            
2401            
2402            
2403            
2404            
2405            
2406            
2407            
2408            
2409            
2410            
2411            
2412            
2413            
2414            
2415            
2416            
2417            
2418            
2419            
2420            
2421            
2422            
2423            
2424            
2425            
2426            
2427            
2428            
2429            
2430            
2431            
2432            
2433            
2434            
2435            
2436            
2437            
2438            
2439            
2440            
Voronoi edges (blue)
as a poor medial
axis approximation

2441            
2442            
2443            
2444            
2445            
2446            
2447            
2448            
2449            
2450            
2451            
2452            
2453            
2454            
2455            
2456            
2457            
2458            
2459            
2460            
2461            
2462            
2463            
2464            
2465            
2466            
2467            
2468            
2469            
2470            
2471            
2472            
2473            
2474            
2475            
2476            
2477            
2478            
2479            
2480            
2481            
2482            
2483            
2484            
2485            
2486            
2487            
2488            
2489            
2490            
2491            
2492            
2493            
2494            
2495            
2496            
2497            
2498            
2499            
2500            
2501            
2502            
2503            
2504            
2505            
2506            
2507            
2508            
2509            
2510            
2511            
2512            
2513            
2514            
2515            
2516            
2517            
2518            
2519            
improved approximation
after calling mic_adust()

2520            
2521              
2522             Warning: not yet thoroughly tested; may move elsewhere
2523              
2524             One use of the Voronoi diagram of a tessellated polygon is to derive an
2525             approximation of the polygon's medial axis by pruning infinite rays and perhaps
2526             trimming or refining remaining branches. The approximation improves as
2527             intervals between sample points on the polygon become shorter. But it's not
2528             always desirable to multiply the number of polygon points to achieve short
2529             intervals.
2530              
2531             At any point on the true medial axis, there is a maximum inscribed circle,
2532             with it's center on the medial axis, and tangent to the polygon in at least
2533             two places.
2534              
2535             The C function moves each Voronoi node so that it becomes the
2536             center of a circle that is tangent to the polygon at two points. In simple
2537             cases this is a maximum inscribed circle, and the point is on the medial axis.
2538             And when it's not, it still should be a much better approximation than the
2539             original point location. The radius to the tangent on the polygon is stored
2540             with the updated Voronoi node.
2541              
2542             After calling C, the modified Voronoi topology can be used as a
2543             list of maximum inscribed circles, from which can be derive a straighter,
2544             better medial axis approximation, without having to increase the number of
2545             sample points on the polygon.
2546              
2547             ($topo, $voronoi_topo) = $tri->triangulate('e');
2548              
2549             mic_adjust($topo, $voronoi_topo); # modifies $voronoi_topo in place
2550            
2551             foreach my $node (@{$voronoi_topo->{nodes}}) {
2552             $mic_center = $node->{point};
2553             $mic_radius = $node->{radius};
2554             ...
2555             }
2556              
2557             Constructing a true medial axis is much more involved - a subject for a
2558             different module. Until that module appears, running topology through
2559             C and then walking and pruning the Voronoi topology might help
2560             fill the gap.
2561              
2562             =head1 API STATUS
2563              
2564             Currently Triangle's option strings are exposed to give more complete access to
2565             its features. More of these options, and perhaps certain common combinations of
2566             them, will likely be wrapped in method-call getter-setters. I would prefer to
2567             preserve the ability to use the option strings directly, but it may be better
2568             at some point to hide them completely behind a more descriptive interface.
2569              
2570              
2571             =head1 AUTHOR
2572              
2573             Michael E. Sheldrake, C<< >>
2574              
2575             Triangle's author is Jonathan Richard Shewchuk
2576              
2577              
2578             =head1 BUGS
2579              
2580             Please report any bugs or feature requests to
2581              
2582             C
2583              
2584             or through the web interface at
2585              
2586             L
2587              
2588             I will be notified, and then you'll automatically be notified of progress on
2589             your bug as I make changes.
2590              
2591              
2592             =head1 SUPPORT
2593              
2594             You can find documentation for this module with the perldoc command.
2595              
2596             perldoc Math::Geometry::Delaunay
2597              
2598              
2599             You can also look for information at:
2600              
2601             =over 4
2602              
2603             =item * RT: CPAN's request tracker
2604              
2605             L
2606              
2607             =item * Search CPAN
2608              
2609             L
2610              
2611             =back
2612              
2613              
2614             =head1 ACKNOWLEDGEMENTS
2615              
2616             Thanks go to Far Leaves Tea in Berkeley for providing oolongs and refuge,
2617             and a place for paths to intersect.
2618              
2619              
2620             =head1 LICENSE AND COPYRIGHT
2621              
2622             Copyright 2013 Micheal E. Sheldrake.
2623              
2624             This Perl binding to Triangle is free software;
2625             you can redistribute it and/or modify it under the terms of either:
2626             the GNU General Public License as published by the Free Software Foundation;
2627             or the Artistic License.
2628              
2629             See http://dev.perl.org/licenses/ for more information.
2630              
2631             =head2 Triangle license
2632              
2633             B by Jonathan Richard Shewchuk, copyright 2005, includes the following
2634             notice in the C source code. Please refer to the C source, included in with this
2635             Perl module distribution, for the full notice.
2636              
2637             This program may be freely redistributed under the condition that the
2638             copyright notices (including this entire header and the copyright
2639             notice printed when the `-h' switch is selected) are not removed, and
2640             no compensation is received. Private, research, and institutional
2641             use is free. You may distribute modified versions of this code UNDER
2642             THE CONDITION THAT THIS CODE AND ANY MODIFICATIONS MADE TO IT IN THE
2643             SAME FILE REMAIN UNDER COPYRIGHT OF THE ORIGINAL AUTHOR, BOTH SOURCE
2644             AND OBJECT CODE ARE MADE FREELY AVAILABLE WITHOUT CHARGE, AND CLEAR
2645             NOTICE IS GIVEN OF THE MODIFICATIONS. Distribution of this code as
2646             part of a commercial system is permissible ONLY BY DIRECT ARRANGEMENT
2647             WITH THE AUTHOR. (If you are not directly supplying this code to a
2648             customer, and you are instead telling them how they can obtain it for
2649             free, then you are not required to make any arrangement with me.)
2650              
2651             =cut
2652              
2653             1;