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	sub	time	code
1				package Math::Geometry::Delaunay;
2
3	1	1	14218	use 5.008;
	1		3
	1		34
4	1	1	3	use warnings;
	1		1
	1		18
5	1	1	4	use strict;
	1		3
	1		22
6	1	1	4	use Carp qw(carp);;
	1		0
	1		46
7	1	1	4	use Exporter();
	1		2
	1		43
8
9				our @ISA = qw(Exporter);
10				our $VERSION;
11
12				BEGIN {
13	1	1	4	use XSLoader;
	1		1
	1		33
14	1	1	1	$VERSION = '0.18';
15	1		602	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]-100cos($amid),$midpt->[1]-100sin($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($dispsizex2+$dispsizey2))/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.18
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;