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	23623	use 5.008;
	1		5
	1		42
4	1	1	6	use warnings;
	1		2
	1		27
5	1	1	5	use strict;
	1		6
	1		31
6	1	1	5	use Carp qw(carp);;
	1		1
	1		68
7	1	1	6	use Exporter();
	1		2
	1		61
8
9				our @ISA = qw(Exporter);
10				our $VERSION;
11
12				BEGIN {
13	1	1	6	use XSLoader;
	1		2
	1		46
14	1	1	2	$VERSION = '0.17';
15	1		1116	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[$_]} @{$_}],edges => [], neighbors => [], marker => undef, attributes => [], index => $tcnt++ };map {push @{$_->{elements}},$ele} @{$ele->{nodes}};$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};$isVoronoi = 1 if ($_->[0] == -1 \|\| $_->[1] == -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				return {
523				nodes => \@nodes,
524				edges => \@edges,
525				segments => \@segs,
526				elements => \@eles
527				};
528				}
529
530				sub get_point_in_polygon {
531				my $poly = shift;
532				my $point_inside;
533
534				my $bottom_left_index=0;
535				my $maxy = $poly->[$bottom_left_index]->[1];
536				for (my $i=1;$i<@{$poly};$i++) {
537				if ($poly->[$i]->[1] <= $poly->[$bottom_left_index]->[1]) {
538				if ($poly->[$i]->[1] < $poly->[$bottom_left_index]->[1] \|\|
539				$poly->[$i]->[0] < $poly->[$bottom_left_index]->[0]
540				) {
541				$bottom_left_index = $i;
542				}
543				}
544				if ($maxy < $poly->[$i]->[1]) { $maxy = $poly->[$i]->[1] }
545				}
546				my $prev_index = $bottom_left_index;
547				my $next_index = -1 * @{$poly} + $bottom_left_index;
548				--$prev_index
549				while ($poly->[$prev_index]->[0] == $poly->[$bottom_left_index]->[0] &&
550				$poly->[$prev_index]->[1] == $poly->[$bottom_left_index]->[1]
551				);
552				++$next_index
553				while ($poly->[$next_index]->[0] == $poly->[$bottom_left_index]->[0] &&
554				$poly->[$next_index]->[1] == $poly->[$bottom_left_index]->[1]
555				);
556
557				my @vec1 = ($poly->[$bottom_left_index]->[0] - $poly->[$prev_index]->[0],
558				$poly->[$bottom_left_index]->[1] - $poly->[$prev_index]->[1]);
559				my @vec2 = ($poly->[$next_index]->[0] - $poly->[$bottom_left_index]->[0],
560				$poly->[$next_index]->[1] - $poly->[$bottom_left_index]->[1]);
561				my $orient = (($vec1[0]$vec2[1] - $vec2[0]$vec1[1])>=0) ? 1:0;
562
563				my $angle1 = atan2($poly->[$prev_index]->[1] - $poly->[$bottom_left_index]->[1], $poly->[$prev_index]->[0] - $poly->[$bottom_left_index]->[0]);
564				my $angle2 = atan2($poly->[$next_index]->[1] - $poly->[$bottom_left_index]->[1], $poly->[$next_index]->[0] - $poly->[$bottom_left_index]->[0]);
565				my $angle;
566				if ($orient) {$angle = $angle2 + ($angle1 - $angle2)/2;}
567				else {$angle = $angle1 + ($angle2 - $angle1)/2;}
568				my $cosangle = cos($angle);
569				my $sinangle = sin($angle);
570				my $tanangle = $sinangle/$cosangle;
571
572				my $adequate_distance = 1.1 * ($maxy - $poly->[$bottom_left_index]->[1]) *
573				((abs($cosangle) < 0.00000001) ? 1 : 1/$sinangle);
574				my $point_wayout = [$poly->[$bottom_left_index]->[0] + $adequate_distance * $cosangle,
575				$poly->[$bottom_left_index]->[1] + $adequate_distance * $sinangle];
576
577				my @intersections = sort {$a->[2] <=> $b->[2]} ray_from_index_poly_intersections($bottom_left_index,$point_wayout,$poly,1);
578
579				if (!@intersections) {
580				print "Warning: Failed to calculate hole or region indicator point.";
581				}
582				elsif ($#intersections % 2 != 0) {
583				print "Warning: Calculated hole or region indicator point is not inside its polygon.";
584				}
585				else {
586				my $closest_intersection = $intersections[0];
587				$point_inside = [($poly->[$bottom_left_index]->[0] + $closest_intersection->[0])/2,
588				($poly->[$bottom_left_index]->[1] + $closest_intersection->[1])/2];
589				}
590				return $point_inside;
591				}
592
593				sub ray_from_index_poly_intersections {
594				my $vertind = shift;
595				my $raypt = shift;
596				my $poly = shift;
597				my $doDists = @_ ? shift : 0;
598
599				my $seg1 = [$poly->[$vertind],$raypt];
600				my $x1= $seg1->[0]->[0];
601				my $y1= $seg1->[0]->[1];
602				my $x2= $seg1->[1]->[0];
603				my $y2= $seg1->[1]->[1];
604				my @lowhix=($x2>$x1)?($x1,$x2):($x2,$x1);
605				my @lowhiy=($y2>$y1)?($y1,$y2):($y2,$y1);
606
607				my @intersections;
608
609				for (my $i = -1; $i < $#$poly; $i++) {
610
611				# skip the segs on either side of the ray base point
612				next if $i == $vertind \|\| ($i + 1) == $vertind \|\| ($vertind == $#$poly && $i == -1);
613
614				my $seg2 = [$poly->[$i],$poly->[$i+1]];
615				my @segsegret;
616
617				my $u1= $seg2->[0]->[0];
618				my $v1= $seg2->[0]->[1];
619				my $u2= $seg2->[1]->[0];
620				my $v2= $seg2->[1]->[1];
621
622				##to maybe optimize for the case where segments are
623				##expected NOT to intersect most of the time
624				#my @lowhix=($x2>$x1)?($x1,$x2):($x2,$x1);
625				#my @lowhiu=($u2>$u1)?($u1,$u2):($u2,$u1);
626				#if (
627				# $lowhix[0]>$lowhiu[1]
628				# \|\|
629				# $lowhix[1]<$lowhiu[0]
630				# ) {
631				# return;
632				# }
633				#my @lowhiy=($y2>$y1)?($y1,$y2):($y2,$y1);
634				#my @lowhiv=($v2>$v1)?($v1,$v2):($v2,$v1);
635				#if (
636				# $lowhiy[0]>$lowhiv[1]
637				# \|\|
638				# $lowhiy[1]<$lowhiv[0]
639				# ) {
640				# return;
641				# }
642
643				my $m1 = ($x2 eq $x1)?'Inf':($y2 - $y1)/($x2 - $x1);
644				my $m2 = ($u2 eq $u1)?'Inf':($v2 - $v1)/($u2 - $u1);
645
646				my $b1;
647				my $b2;
648				my $xi;
649
650				# Arranged like this to avoid m1-m2 with infinity involved, which works
651				# in other contexts, but can trigger a floating point exception here.
652				# Turns out the exception happens when we let Triangle set the FPU
653				# control word, but not when we use XPFPA.h to take care of that, but
654				# leaving it this as it is in case that doesn't fix that everywhere.
655				my $dm;
656				if ($m1 != 'Inf' && $m2 != 'Inf') {$dm = $m1 - $m2;}
657				elsif ($m1 == 'Inf' && $m2 == 'Inf') {return;}
658				else {$dm='Inf';}
659				if ($m1 == 'Inf' && $m2 != 'Inf') {$xi = $x1;$b2 = $v1 - ($m2 * $u1);}
660				elsif ($m2 == 'Inf' && $m1 != 'Inf') {$xi = $u1;$b1 = $y1 - ($m1 * $x1);}
661				elsif (abs($dm) > 0.000000000001) {
662				$b1 = $y1 - ($m1 * $x1);
663				$b2 = $v1 - ($m2 * $u1);
664				$xi=($b2-$b1)/$dm;
665				}
666				my @lowhiu=($u2>$u1)?($u1,$u2):($u2,$u1);
667				if ($m1 != 'Inf') {
668				if ($m2 == 'Inf' && ($u2<$lowhix[0] \|\| $u2>$lowhix[1]) ) {
669				next;
670				}
671				if (
672				defined $xi &&
673				($xi < $lowhix[1] \|\| $xi eq $lowhix[1]) &&
674				($xi > $lowhix[0] \|\| $xi eq $lowhix[0]) &&
675				($xi < $lowhiu[1] \|\| $xi eq $lowhiu[1]) &&
676				($xi > $lowhiu[0] \|\| $xi eq $lowhiu[0])
677				) {
678				my $y=($m1*$xi)+$b1;
679				my @lowhiv=($v2>$v1)?($v1,$v2):($v2,$v1);
680				if ($m2 == 'Inf' &&
681				($y<$lowhiv[0] \|\| $y>$lowhiv[1])
682				) {
683				next;
684				}
685				else {
686				push(@intersections,[$xi,$y]);
687				if ($y eq $v1 \|\| $y eq $v2) {
688				$i++; # avoid duplicates at endpoints
689				}
690				}
691				}
692				}
693				elsif ($m2 != 'Inf') {#so $m1 is Inf
694				if (($x1 < $lowhiu[0] \|\| $x1 > $lowhiu[1]) && ! ($x1 eq $lowhiu[0] \|\| $x1 eq $lowhiu[1])) {
695				next;
696				}
697				my @lowhiv=($v2>$v1)?($v1,$v2):($v2,$v1);
698				my $yi = ($m2*$xi)+$b2;
699				if (($yi \|\| $yi eq 0) &&
700				($yi < $lowhiy[1] \|\| $yi eq $lowhiy[1]) &&
701				($yi > $lowhiy[0] \|\| $yi eq $lowhiy[0]) &&
702				($yi < $lowhiv[1] \|\| $yi eq $lowhiv[1]) &&
703				($yi > $lowhiv[0] \|\| $yi eq $lowhiv[0])
704				) {
705				push(@intersections,[$xi,$yi]);
706				if ($xi eq $u1 \|\| $xi eq $u2) {
707				$i++; # avoid duplicates at endpoints
708				}
709				}
710				}
711				}
712				if ($doDists) {
713				foreach my $int (@intersections) {
714				push(@{$int}, sqrt( ($int->[0]-$seg1->[0]->[0])2 + ($int->[1]-$seg1->[0]->[1])2 ) );
715				}
716				}
717				return @intersections;
718				}
719
720				# Adjust location of nodes in voronoi diagram so they become
721				# centers of maximum inscribed circles, and store MIC radius in each node.
722				# This is a first step towards a better medial axis approximation.
723				# It straightens out the initial MA approx. derived from the Voronoi diagram.
724
725				sub mic_adjust {
726				my $topo = shift;
727				my $vtopo = shift;
728
729				my @new_vnode_points; # voronoi vertices moved to MIC center points
730				my @new_vnode_radii; # will be calculated and added to node data
731				my @new_vnode_tangents; # where the MIC touches the boundary PSLG
732
733				my $vnc=-1; # will use to look up triangle that corresponds to the voronoi node
734				# $vnc can probably be replaced with $vnode->{index}
735
736				foreach my $vnode (@{$vtopo->{nodes}}) {
737				$vnc++;
738
739				push(@new_vnode_points,$vnode->{point});
740				push(@new_vnode_radii,undef);
741				push(@new_vnode_tangents,[]);
742
743				my $boundary_edge;
744
745				if (@{$vnode->{edges}} == 3) {
746				my @all_edges = map {$topo->{edges}->[$_->{index}]} @{$vnode->{edges}};
747				my @boundary_edges = grep {$_->{marker}} @all_edges;
748
749				my @opposite_boundary_edges;
750				my @opposite_boundary_feet;
751
752				my $branch_node_edge_index1;
753				my $branch_node_edge_index2;
754				my $branch_node_third_tan;
755
756				# BRANCH NODE
757				# The corresponding Delaunay triangle has no edges on the PSLG boundary,
758				# but touches the boundary at its vertices. This corresponds to a
759				# node with three edges in the medial axis approximation.
760				if (@boundary_edges == 0) {
761				$vnode->{isbranchnode} = 1;
762
763				# Orient nodes in edges emanating from this branch node so that
764				# the first node is always the branch node.
765				# The notion is that direction is always outward from a branch node.
766				# This will be a problem if two branch nodes link to each other.
767				$_->{nodes} = [$vnode,$_->{nodes}->[0]!=$vnode?$_->{nodes}->[0]:$_->{nodes}->[1]] for @{$vnode->{edges}};
768
769				# Make sure corners are sorted so they are opposite the Voronoi edges in order.
770				# They may have already been that way, but couldn't determine from Triangle docs.
771				my @corners;
772				for (my $i=0;$i<@{$vnode->{edges}};$i++) {
773				push @corners, +(grep $_ != $topo->{edges}->[$vnode->{edges}->[$i]->{index}]->{nodes}->[0]
774				&& $_ != $topo->{edges}->[$vnode->{edges}->[$i]->{index}]->{nodes}->[1],
775				@{$topo->{elements}->[$vnc]->{nodes}}
776				)[0];
777				}
778
779				my @corner_boundary_edges = grep $_->{marker}, map @{$_->{edges}}, @corners;
780				my @corner_boundary_feet = map {getFoot([$_->{nodes}->[0]->{point},$_->{nodes}->[1]->{point}],$vnode->{point}->[0],$vnode->{point}->[1])} @corner_boundary_edges;
781
782				# Handle the case where the three corners are probably the three MIC tangents.
783				if (! grep $_, @corner_boundary_feet) {
784				$new_vnode_radii[-1] = dist2d($vnode->{point},$topo->{edges}->[$vnode->{edges}->[0]->{index}]->{nodes}->[0]->{point});
785				$new_vnode_tangents[-1] = [[$corners[2]->{point},
786				$corners[1]->{point}],
787				[$corners[0]->{point},
788				$corners[2]->{point}],
789				[$corners[1]->{point},
790				$corners[0]->{point}]];
791				next;
792				}
793
794				# Otherwise, set up a case similar to the non-branch node case
795				# by designating feet, a boundary edge, and an opposite edge.
796				my @boundary_feet_edges;
797				for (my $i=0;$i<@corner_boundary_feet;$i+=2) {
798				my $foot;
799				my $edge;
800				# if no foot was found on the two edges meeting at
801				# one of the triangle's vertices, make the vertex a "fake foot"
802				if (!$corner_boundary_feet[$i] && !$corner_boundary_feet[$i+1]) {
803				my $foot = ( $corner_boundary_edges[$i]->{nodes}->[0] == $corner_boundary_edges[$i+1]->{nodes}->[0]
804				\|\| $corner_boundary_edges[$i]->{nodes}->[0] == $corner_boundary_edges[$i+1]->{nodes}->[1] )
805				? $corner_boundary_edges[$i]->{nodes}->[0]->{point}
806				: $corner_boundary_edges[$i]->{nodes}->[1]->{point};
807				# edge evaluates to "false" to signal this case
808				push @boundary_feet_edges, [$foot,undef,$i/2];
809				}
810				# otherwise keep foot and edge ref
811				else {
812				if ($corner_boundary_feet[$i] ) {push @boundary_feet_edges, [$corner_boundary_feet[$i], $corner_boundary_edges[$i] ,$i/2];}
813				if ($corner_boundary_feet[$i+1]
814				&& ( !$corner_boundary_feet[$i] \|\|
815				# screen out duplicates
816				( $corner_boundary_feet[$i+1]->[0] ne $corner_boundary_feet[$i]->[0]
817				&& $corner_boundary_feet[$i+1]->[1] ne $corner_boundary_feet[$i]->[1]
818				)
819				)
820				) {push @boundary_feet_edges, [$corner_boundary_feet[$i+1],$corner_boundary_edges[$i+1],$i/2];}
821				}
822				}
823
824				@boundary_feet_edges = sort {dist2d($a->[0],$vnode->{point}) <=> dist2d($b->[0],$vnode->{point})} @boundary_feet_edges;
825
826				# closest edge treated as boundary edge, similar to non-branch node handling
827				my $first_with_edge = +(grep {$_->[1]} @boundary_feet_edges)[0];
828
829				$boundary_edge = {nodes=>[$first_with_edge->[1]->{nodes}->[0],$first_with_edge->[1]->{nodes}->[1]]};
830
831				# next closest edge treated as opposite boundary edge, similar to non-branch node handling
832				# (that is, next closest that's also not connected to the same
833				# corner as the closest)
834				my $first_not_first_with_edge = +(grep {$_ != $first_with_edge && $_->[2] != $first_with_edge->[2]} @boundary_feet_edges)[0];
835
836				$opposite_boundary_feet[0] = $first_not_first_with_edge->[0];
837				$opposite_boundary_edges[0] = $first_not_first_with_edge->[1];
838
839				# Use these later to match up tangent point pairs with Voronoi edges.
840				$branch_node_edge_index1 = $first_with_edge->[2];
841				$branch_node_edge_index2 = $first_not_first_with_edge->[2];
842
843				my @threst = grep { $_ != $first_with_edge
844				&& $_ != $first_not_first_with_edge
845				&& $_->[2] != $first_with_edge->[2]
846				&& $_->[2] != $first_not_first_with_edge->[2]
847				} @boundary_feet_edges;
848				# Not completely thought out, but results look good so far -
849				# This "tangent point" is not (generally) touched by the approximate
850				# MIC, but it might be worth trying to shift the approx MIC in
851				# a later step to come closer to this, when possible.
852				# Even without that extra shift attempt, it's useful to have
853				# this when reconstructing a polygon from the approximate
854				# medial axis enabled by mic_adjust().
855				$branch_node_third_tan = $threst[0]->[0];
856				}
857
858				# NON-BRANCH NODE
859				# The corresponding Delaunay triangle has one or two edges
860				# on the PSLG boundary. This corresponds to a node with two edges
861				# in the medial axis approximation.
862				if (@boundary_edges == 1 \|\| @boundary_edges == 2) {
863
864				$boundary_edge = $boundary_edges[0];
865
866				my $bef = getFoot([$boundary_edge->{nodes}->[0]->{point},$boundary_edge->{nodes}->[1]->{point}],$vnode->{point}->[0],$vnode->{point}->[1]);
867				my $ber = sqrt(($vnode->{point}->[0]-$bef->[0])2 + ($vnode->{point}->[1]-$bef->[1])2);
868
869				if (@boundary_edges == 2) {
870				# If the other boundary edge is closer, make that the boundary edge.
871				my $obef = getFoot([$boundary_edges[1]->{nodes}->[0]->{point},$boundary_edges[1]->{nodes}->[1]->{point}],$vnode->{point}->[0],$vnode->{point}->[1]);
872				next if (!$obef);
873				my $ober = sqrt(($vnode->{point}->[0]-$obef->[0])2 + ($vnode->{point}->[1]-$obef->[1])2);
874				if ($ober < $ber) {
875				$boundary_edge = $boundary_edges[1];
876				$bef=$obef;
877				$ber=$ober;
878				}
879				}
880
881				my @other_edges = grep $_ != $boundary_edge, map $topo->{edges}->[$_->{index}], @{$vnode->{edges}};
882				my $opposite_node = +(grep {$_ != $boundary_edge->{nodes}->[0] && $_ != $boundary_edge->{nodes}->[1]} @{$other_edges[1]->{nodes}})[0];
883				@opposite_boundary_edges = grep {$_->{marker}} @{$opposite_node->{edges}};
884				@opposite_boundary_feet = map {getFoot([$_->{nodes}->[0]->{point},$_->{nodes}->[1]->{point}],$vnode->{point}->[0],$vnode->{point}->[1])} @opposite_boundary_edges;
885
886				if (@boundary_edges == 1
887				# \|\| @boundary_edges == 2
888				) { # should apply to ==2 too, maybe, but should screen out the "opposite" that is also "adjacent" in that case
889				my @adjacent_boundary_edges = grep {$_->{marker} && $_ != $boundary_edge} map @{$_->{edges}}, @{$boundary_edge->{nodes}};
890				my @adjacent_boundary_feet = map {getFoot([$_->{nodes}->[0]->{point},$_->{nodes}->[1]->{point}],$vnode->{point}->[0],$vnode->{point}->[1])} @adjacent_boundary_edges;
891				# if any adjacent bounds closer, replace $boundary_edge with closest
892				for (my $i = 0; $i < @adjacent_boundary_feet; $i++) {
893				next if (!$adjacent_boundary_feet[$i]);
894				my $newfootdist = sqrt(($vnode->{point}->[0] - $adjacent_boundary_feet[$i]->[0])2 + ($vnode->{point}->[1]-$adjacent_boundary_feet[$i]->[1])2);
895				if ($newfootdist < $ber) {
896				$boundary_edge = $adjacent_boundary_edges[$i];
897				$bef = $adjacent_boundary_feet[$i];
898				$ber = $newfootdist;
899				}
900				}
901				}
902
903				if (grep $_, @opposite_boundary_feet) {
904				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);
905				@opposite_boundary_feet = map $opposite_boundary_feet[$_], @sortind;
906				@opposite_boundary_edges = map $opposite_boundary_edges[$_], @sortind;
907				@opposite_boundary_feet = ($opposite_boundary_feet[0]);
908				@opposite_boundary_edges = ($opposite_boundary_edges[0]);
909				my $obr = dist2d($opposite_boundary_feet[0],$vnode->{point});
910				if ($ber>$obr) {
911				my $tmp = $opposite_boundary_edges[0];
912				$opposite_boundary_edges[0] = $boundary_edge;
913				$boundary_edge = $tmp;
914				$opposite_boundary_feet[0] = $bef;
915				}
916				}
917
918				if (!grep $_, @opposite_boundary_feet) {
919				# make the foot just the opposite point
920				@opposite_boundary_feet = ($opposite_node->{point});
921				# make edge eval to false to signal this fake foot case
922				@opposite_boundary_edges = map {0} @opposite_boundary_edges;
923				}
924
925				}
926
927				# FIND THE TWO TANGENT POINTS, RADIUS, AND SHIFT POINT TO MIC CENTER
928
929				# Only one defined unique foot should be in @opposite_boundary_feet
930				# at this point, though that wasn't always the case in the past.
931				# When we're satisfied that it will be the case in the future, we'll
932				# remove the for loop.
933
934				for (my $i = 0; $i < @opposite_boundary_feet; $i++) {
935				next if !$opposite_boundary_feet[$i];
936				my $foot = $opposite_boundary_feet[$i];
937				my $opp_edge = $opposite_boundary_edges[$i];
938
939				my $a1;
940				if ($opp_edge) {
941				$a1 = atan2($opp_edge->{nodes}->[0]->{point}->[1] - $opp_edge->{nodes}->[1]->{point}->[1],
942				$opp_edge->{nodes}->[0]->{point}->[0] - $opp_edge->{nodes}->[1]->{point}->[0]);
943				}
944				else { # the "fake foot" case, where opposite edge is just a point
945				$a1 = atan2($foot->[1] - $vnode->{point}->[1],
946				$foot->[0] - $vnode->{point}->[0]);
947				$a1 -= $pi / 2;
948				}
949
950				my $a2 = atan2($boundary_edge->{nodes}->[1]->{point}->[1] - $boundary_edge->{nodes}->[0]->{point}->[1],
951				$boundary_edge->{nodes}->[1]->{point}->[0] - $boundary_edge->{nodes}->[0]->{point}->[0]);
952
953				$a1 = angle_reduce_pi($a1);
954
955				my $amid = ($a1 + $a2) / 2;
956				my $amidnorm = $amid + $pi / 2;
957
958				my $boundtanpt = line_line_intersection(
959				[ $boundary_edge->{nodes}->[0]->{point}, $boundary_edge->{nodes}->[1]->{point} ],
960				[ $foot, [$foot->[0] + 100 * cos($amidnorm), $foot->[1] + (100 * sin($amidnorm))] ],
961				);
962
963				if ($boundtanpt) {
964				my $midpt=[($foot->[0]+$boundtanpt->[0])/2,($foot->[1]+$boundtanpt->[1])/2];
965				my $nother_mid_pt=[$midpt->[0]-100cos($amid),$midpt->[1]-100sin($amid)];
966				my $center = line_line_intersection([$vnode->{point},$foot],[$midpt,$nother_mid_pt]);
967				if ($center) {
968				$new_vnode_points[-1] = $center;
969				$new_vnode_radii[-1] = dist2d($center,$foot);
970				# assign the three tangent pairs to a branch node
971				if (defined $branch_node_edge_index1) {
972				my $gets_both_found = +(grep $_ != $branch_node_edge_index1 && $_ != $branch_node_edge_index2, (0..2))[0];
973				$new_vnode_tangents[-1]->[( $gets_both_found ) % 3] = [$foot, $boundtanpt];
974				$new_vnode_tangents[-1]->[( $branch_node_edge_index1 ) % 3] = [$branch_node_third_tan, $foot];
975				$new_vnode_tangents[-1]->[( $branch_node_edge_index2 ) % 3] = [$boundtanpt, $branch_node_third_tan];
976				}
977				# assign the two tangent pairs for non-branch nodes
978				else {
979				foreach my $edge (@{$vnode->{edges}}) {
980				next if defined($edge->{vector}) && ($edge->{vector}->[0] != 0 \|\| $edge->{vector}->[1] != 0); # a ray
981				push @{$new_vnode_tangents[-1]}, [$foot, $boundtanpt];
982				}
983				}
984				}
985				}
986				}
987				}
988
989				if (!defined $new_vnode_radii[-1]) {
990				# This would be a branch node that had feet, but didn't end up
991				# getting adjusted. This is either an okay edge case or a failure.
992				# Let's watch for a while and see if we end up here.
993				$new_vnode_radii[-1] = dist2d($vnode->{point},$topo->{edges}->[$vnode->{edges}->[0]->{index}]->{nodes}->[0]->{point});
994				print "\nFailure to find radius in Math::Geometery::Delaunay::mic_adjust().\nThe developer would like to know if you come across this.\n";
995				}
996
997				}
998
999				# Can probably integrate the node point, radius, and tangents assignments
1000				# directly into the above section, and get rid of this temp array stuff.
1001				# On the other hand, if we move this toward using matched index lists (more
1002				# like Triangle's native output) might just export the lists we made, and not
1003				# update the crossref'ed hash structure.
1004				for (my $i = 0; $i < @{$vtopo->{nodes}}; $i++) {
1005				$vtopo->{nodes}->[$i]->{point} = $new_vnode_points[$i];
1006				$vtopo->{nodes}->[$i]->{radius} = $new_vnode_radii[$i];
1007				$vtopo->{nodes}->[$i]->{tangents} = $new_vnode_tangents[$i];
1008				}
1009
1010				# Fixup left-right ordering of tangent points, now that
1011				# vnodes and their edges should all be nicely centered between them.
1012
1013				for (my $i = 0; $i < @{$vtopo->{nodes}}; $i++) {
1014				my $vnode = $vtopo->{nodes}->[$i];
1015
1016				# first pass - could you do iswithin for vnode in its corresponding triangle
1017				# and if not, shift toward next/previous vnode until it is within
1018				# like to intersection of line between tri apex and bound edge
1019				# and line between point and next/prev point... or is there trajectory
1020				# based on two (usually) boundary edges where the tangents are?
1021				#
1022
1023				my @fixtangents;
1024				my @edges = grep !defined($_->{vector}) \|\| ($_->{vector}->[0] == 0 && $_->{vector}->[1] == 0), @{$vnode->{edges}};
1025
1026				for (my $j = 0; $j < @edges; $j++) {
1027				my $edge = $edges[$j];
1028
1029				my $tangents = $vnode->{tangents}->[$j]; # tangent pairs correspond to non-ray edges, in order
1030				my $other_node = ($edge->{nodes}->[0] != $vnode) ? $edge->{nodes}->[0] : $edge->{nodes}->[1];
1031				my $start_node = $vnode;
1032
1033				# Duplicate node points can happen for common reasons, so go out
1034				# to the next node if we got a duplicate. (Three-in-a-row duplicates
1035				# shouldn't happen.)
1036				if ($vnode->{point}->[0] eq $other_node->{point}->[0] && $vnode->{point}->[1] eq $other_node->{point}->[1]) {
1037				my $other_edge = +(grep $_ != $edge && !(defined $_->{vector} && ($_->{vector}->[0] != 0 \|\| $_->{vector}->[1] != 0)), @{$other_node->{edges}})[0];
1038				if ($other_edge) {
1039				$other_node = $other_edge->{nodes}->[0] != $other_node ? $other_edge->{nodes}->[0] : $other_edge->{nodes}->[1];
1040				}
1041				else {
1042				# Well, then, if this is a non-branch node, set up a test line
1043				# using the previous node, heading into this one.
1044				if (@{$vnode->{tangents}} == 2) {
1045				my $other_edge = $edges[$j == 0 ? 1 : 0];
1046				$start_node = ($other_edge->{nodes}->[0] != $vnode) ? $other_edge->{nodes}->[0] : $other_edge->{nodes}->[1];
1047				$other_node = $vnode;
1048				}
1049				#else { die "Couldn't get other edge in duplicate point case\n\n;" }
1050				}
1051				}
1052
1053				# Heading out from the start node to the other node, the first tangent associated
1054				# with this edge should be on the left (a counterclockwise turn), and the
1055				# other tangent on the right. We've got Triangle's adaptave robust
1056				# orientation code backing up counterclockwise(), for the really close cases.
1057				my $ccwl = counterclockwise($start_node->{point}, $other_node->{point}, $tangents->[0]);
1058				my $ccwr = counterclockwise($start_node->{point}, $other_node->{point}, $tangents->[1]);
1059
1060				if ($ccwl < 0) { # first tangent wasn't on the left
1061				if ($ccwr > 0) { # but the second was, so we just need to swap them
1062				@{$tangents} = reverse @{$tangents};
1063				}
1064				else { # But if both were on the right, something's wrong.
1065				# If this is a branch node, and the other two tangent pairs
1066				# don't have the same problem, we can fix up the bad one later.
1067				#push @fixtangents, $j;
1068				$vnode->{fixtangents} = [] if !defined $vnode->{fixtangents};
1069				push @{$vnode->{fixtangents}}, $j;
1070				}
1071				}
1072				elsif ($ccwr > 0) { # Both were on the left. Maybe fix up later.
1073				#push @fixtangents, $j;
1074				$vnode->{fixtangents} = [] if !defined $vnode->{fixtangents};
1075				push @{$vnode->{fixtangents}}, $j;
1076				}
1077				#elsif ($ccwl eq 0) { die "zero" } # Shouldn't happen. Leave it alone if it does.
1078				}
1079				}
1080				for (my $i = 0; $i < @{$vtopo->{nodes}}; $i++) {
1081				if (@{$vtopo->{nodes}->[$i]->{tangents}} == 3 && defined $vtopo->{nodes}->[$i]->{fixtangents} && @{$vtopo->{nodes}->[$i]->{fixtangents}} == 1) {
1082				#print "FIXUP TANGENT PAIR FOR ONE BRANCH AT BRANCH NODE\n";
1083				#$vtopo->{nodes}->[$i]->{color} = 'orange';
1084				# one bad tangent pair out of three for a branch node
1085				# infer that it just needs to be reversed if other two pairs were good
1086
1087				# debug stuff
1088				#my @edges = grep !defined($_->{vector}) \|\| ($_->{vector}->[0] == 0 && $_->{vector}->[1] == 0), @{$vtopo->{nodes}->[$i]->{edges}};
1089				#my $reconsider_edge = $edges[$vtopo->{nodes}->[$i]->{fixtangents}->[0]];
1090				#my $reconsider_node = $reconsider_edge->{nodes}->[0] != $vtopo->{nodes}->[$i] ? $reconsider_edge->{nodes}->[0]:$reconsider_edge->{nodes}->[1];
1091				#$reconsider_node->{color} = 'green';
1092
1093				$vtopo->{nodes}->[$i]->{tangents}->[$vtopo->{nodes}->[$i]->{fixtangents}->[0]] =
1094				[$vtopo->{nodes}->[$i]->{tangents}->[($vtopo->{nodes}->[$i]->{fixtangents}->[0] + 1) % 3]->[1],
1095				$vtopo->{nodes}->[$i]->{tangents}->[($vtopo->{nodes}->[$i]->{fixtangents}->[0] - 1) ]->[0]];
1096
1097				}
1098				#if (@{$vtopo->{nodes}->[$i]->{tangents}} == 3 && defined $vtopo->{nodes}->[$i]->{fixtangents} && @{$vtopo->{nodes}->[$i]->{fixtangents}} != 1) {
1099				#print " MORE NEEDED FIXUP: ",scalar(@{$vtopo->{nodes}->[$i]->{fixtangents}}),"\n";
1100				# if more than one branch had tangents mixed up...
1101				# maybe the other one or two wrong ones are really right?
1102				# ambiguous.
1103				# $vtopo->{nodes}->[$i]->{color} = 'orange';
1104				# }
1105				if ( @{$vtopo->{nodes}->[$i]->{tangents}} == 2
1106				&& $vtopo->{nodes}->[$i]->{tangents}->[0]->[0] == $vtopo->{nodes}->[$i]->{tangents}->[1]->[0]
1107				) {
1108				#print "non-branch maybe needs fixup.\n";
1109				# not sure if this a good way to approach this -
1110				# doesn't definitively capture all that can go
1111				# wrong with non-branch node orientations.
1112				$vtopo->{nodes}->[$i]->{color} = 'green';
1113				if (defined $vtopo->{nodes}->[$i]->{fixtangents} && @{$vtopo->{nodes}->[$i]->{fixtangents}} == 1) {
1114				$vtopo->{nodes}->[$i]->{tangents}->[$vtopo->{nodes}->[$i]->{fixtangents}->[0]] =
1115				[$vtopo->{nodes}->[$i]->{tangents}->[$vtopo->{nodes}->[$i]->{fixtangents}->[0] - 1]->[1],
1116				$vtopo->{nodes}->[$i]->{tangents}->[$vtopo->{nodes}->[$i]->{fixtangents}->[0] - 1]->[0]];
1117				}
1118				}
1119				}
1120				}
1121
1122				sub getFoot {
1123				my $seg = shift;
1124				my $x = shift;
1125				my $y = shift;
1126				my $foot;
1127				my $m = ($seg->[1]->[0] - $seg->[0]->[0] == 0) ? 'inf' : ($seg->[1]->[1] - $seg->[0]->[1])/($seg->[1]->[0] - $seg->[0]->[0]);
1128				my @sortx = $seg->[0]->[0] < $seg->[1]->[0] ? ($seg->[0]->[0], $seg->[1]->[0]) : ($seg->[1]->[0], $seg->[0]->[0]);
1129				my @sorty = $seg->[0]->[1] < $seg->[1]->[1] ? ($seg->[0]->[1], $seg->[1]->[1]) : ($seg->[1]->[1], $seg->[0]->[1]);
1130				if ($m == 0) {
1131				if ($x >= $sortx[0] && $x <= $sortx[1]) {$foot=[$x, $seg->[0]->[1]];}
1132				}
1133				elsif ($m == 'inf') {
1134				if ($y >= $sorty[0] && $y <= $sorty[1]) {$foot=[$seg->[0]->[0], $y];}
1135				}
1136				else {
1137				my $intersect_x = (($m$seg->[0]->[0])-($seg->[0]->[1])+((1/$m)$x)+($y))/($m+(1/$m));
1138				if ($intersect_x >= $sortx[0] && $intersect_x <= $sortx[1]) {
1139				my $intersect_y = -($seg->[0]->[0] - $intersect_x) * $m + $seg->[0]->[1];
1140				if ($intersect_y >= $sorty[0] && $intersect_y <= $sorty[1]) {
1141				$foot = [$intersect_x, $intersect_y];
1142				}
1143				}
1144				}
1145				return $foot;
1146				}
1147
1148				sub angle_reduce {
1149				my $a=shift;
1150				while($a > $pi / 2) { $a -= $pi; }
1151				while($a <= -$pi / 2) { $a += $pi; }
1152				return $a;
1153				}
1154
1155				sub angle_reduce_pi {
1156				my $a=shift;
1157				while($a > $pi) { $a -= $pi * 2; }
1158				while($a <= -$pi) { $a += $pi * 2; }
1159				return $a;
1160				}
1161
1162				sub seg_seg_intersection {
1163				my $seg1 = shift;
1164				my $seg2 = shift;
1165				my $int;
1166
1167				my $x1= $seg1->[0]->[0]; my $y1= $seg1->[0]->[1];
1168				my $x2= $seg1->[1]->[0]; my $y2= $seg1->[1]->[1];
1169				my $u1= $seg2->[0]->[0]; my $v1= $seg2->[0]->[1];
1170				my $u2= $seg2->[1]->[0]; my $v2= $seg2->[1]->[1];
1171
1172				my $m1 = ($x2 eq $x1)?'Inf':($y2 - $y1)/($x2 - $x1);
1173				my $m2 = ($u2 eq $u1)?'Inf':($v2 - $v1)/($u2 - $u1);
1174
1175				my $b1;
1176				my $b2;
1177				my $xi;
1178
1179				# Arranged like this to avoid m1-m2 with infinity involved, which works
1180				# in other contexts, but can trigger a floating point exception here.
1181				my $dm;
1182				if ($m1 != 'Inf' && $m2 != 'Inf') {$dm = $m1 - $m2;}
1183				elsif ($m1 == 'Inf' && $m2 == 'Inf') {return;}
1184				else {$dm='Inf';}
1185
1186				if ($m1 == 'Inf' && $m2 != 'Inf') {$xi = $x1;$b2 = $v1 - ($m2 * $u1);}
1187				elsif ($m2 == 'Inf' && $m1 != 'Inf') {$xi = $u1;$b1 = $y1 - ($m1 * $x1);}
1188				elsif (abs($dm) > 0.000000000001) {
1189				$b1 = $y1 - ($m1 * $x1);
1190				$b2 = $v1 - ($m2 * $u1);
1191				$xi=($b2-$b1)/$dm;
1192				}
1193				my @lowhiu=($u2>$u1)?($u1,$u2):($u2,$u1);
1194				if ($m1 != 'Inf') {
1195				my @lowhix=($x2>$x1)?($x1,$x2):($x2,$x1);
1196				if ($m2 == 'Inf' && ($u2<$lowhix[0] \|\| $u2>$lowhix[1]) ) {
1197				return;
1198				}
1199				if (
1200				($xi \|\| $xi eq 0) &&
1201				($xi < $lowhix[1] \|\| $xi eq $lowhix[1]) &&
1202				($xi > $lowhix[0] \|\| $xi eq $lowhix[0]) &&
1203				($xi < $lowhiu[1] \|\| $xi eq $lowhiu[1]) &&
1204				($xi > $lowhiu[0] \|\| $xi eq $lowhiu[0])
1205				) {
1206				my $y=($m1*$xi)+$b1;
1207				my @lowhiv=($v2>$v1)?($v1,$v2):($v2,$v1);
1208				if ($m2 == 'Inf' &&
1209				($y<$lowhiv[0] \|\| $y>$lowhiv[1])
1210				) {
1211				return;
1212				}
1213				else {
1214				$int = [$xi,$y];
1215				}
1216				}
1217				}
1218				elsif ($m2 != 'Inf') { #so $m1 is Inf
1219
1220				if ($x1 < $lowhiu[0] \|\| $x1 > $lowhiu[1] && ! ($x1 eq $lowhiu[0] \|\| $x1 eq $lowhiu[1])) {
1221				return;
1222				}
1223				my @lowhiy=($y2>$y1)?($y1,$y2):($y2,$y1);
1224				my @lowhiv=($v2>$v1)?($v1,$v2):($v2,$v1);
1225				my $yi = ($m2*$xi)+$b2;
1226				if (($yi \|\| $yi eq 0) &&
1227				($yi < $lowhiy[1] \|\| $yi eq $lowhiy[1]) &&
1228				($yi > $lowhiy[0] \|\| $yi eq $lowhiy[0]) &&
1229				($yi < $lowhiv[1] \|\| $yi eq $lowhiv[1]) &&
1230				($yi > $lowhiv[0] \|\| $yi eq $lowhiv[0])
1231				) {
1232				$int = [$xi,$yi];
1233				}
1234				}
1235				return $int;
1236				}
1237
1238				sub line_line_intersection {
1239				my $seg1 = shift;
1240				my $seg2 = shift;
1241				my $int;
1242
1243				my $x1= $seg1->[0]->[0]; my $y1= $seg1->[0]->[1];
1244				my $x2= $seg1->[1]->[0]; my $y2= $seg1->[1]->[1];
1245				my $u1= $seg2->[0]->[0]; my $v1= $seg2->[0]->[1];
1246				my $u2= $seg2->[1]->[0]; my $v2= $seg2->[1]->[1];
1247
1248				my $m1 = ($x2 eq $x1)?'Inf':($y2 - $y1)/($x2 - $x1);
1249				my $m2 = ($u2 eq $u1)?'Inf':($v2 - $v1)/($u2 - $u1);
1250
1251				my $b1;
1252				my $b2;
1253
1254				my $xi;
1255
1256				# Arranged like this to avoid m1-m2 with infinity involved, which works
1257				# in other contexts, but can trigger a floating point exception here.
1258				my $dm;
1259				if ($m1 != 'Inf' && $m2 != 'Inf') {$dm = $m1 - $m2;}
1260				elsif ($m1 == 'Inf' && $m2 == 'Inf') {return;}
1261				else {$dm='Inf';}
1262
1263				if ($m1 == 'Inf' && $m2 != 'Inf') {$xi = $x1;$b2 = $v1 - ($m2 * $u1);}
1264				elsif ($m2 == 'Inf' && $m1 != 'Inf') {$xi = $u1;$b1 = $y1 - ($m1 * $x1);}
1265				elsif (abs($dm) > 0.000000000001) {
1266				$b1 = $y1 - ($m1 * $x1);
1267				$b2 = $v1 - ($m2 * $u1);
1268				$xi=($b2-$b1)/$dm;
1269				}
1270				if ($m1 != 'Inf') {
1271				if (defined $xi) {
1272				my $y=($m1*$xi)+$b1;
1273				$int = [$xi,$y];
1274				}
1275				}
1276				elsif ($m2 != 'Inf') { # so $m1 is Inf
1277				my $yi = ($m2*$xi)+$b2;
1278				if ($yi \|\| $yi eq 0) {
1279				$int = [$xi,$yi];
1280				}
1281				}
1282				return $int;
1283				}
1284
1285				sub to_svg {
1286				my %spec = @_;
1287				my $triios = [defined($spec{topo}) ? delete $spec{topo} : undef, defined($spec{vtopo}) ? delete $spec{vtopo} : undef];
1288				my $fn = defined($spec{file}) ? delete $spec{file} : '-';
1289				my $dispsize = defined($spec{size}) ? delete $spec{size} : [800, 600];
1290				my $triio = $triios->[0];
1291				my $vorio = @{$triios}?$triios->[1]:undef;
1292				my @edges;
1293				my @segs;
1294				my @pts;
1295				my @vpts;
1296				my @vedges;
1297				my @vrays;
1298				my @circles;
1299				my @elements;
1300				my $maxx;
1301				my $maxy;
1302				my $minx;
1303				my $miny;
1304				if (!$triio) {carp "no geometry provided";return;}
1305				foreach my $key ( keys %spec ) { if (ref($spec{$key}) !~ /ARRAY/) { carp("style config for '$key' should be a reference to an array"); return; } }
1306
1307				# make copies of points, because we'll be moving and scaling them
1308
1309				if (ref($triio) =~ /HASH/ && defined $triio->{nodes}) {
1310				push @pts, map [$_,defined $spec{nodes} ? @{$spec{nodes}} : undef], map [@{$_->{point}}], @{$triio->{nodes}};
1311				}
1312				else {
1313				push @pts, map [[@{$_}],defined $spec{nodes} ? @{$spec{nodes}} : undef], ltolol(2,$triio->pointlist);
1314				}
1315				if ($vorio) {
1316				if (ref($vorio) =~ /HASH/ && defined $vorio->{nodes}) {
1317				push @vpts, map [$_,defined $spec{vnodes} ? @{$spec{vnodes}} : undef], map [@{$_->{point}}], @{$vorio->{nodes}};
1318				}
1319				else {
1320				push @vpts, map [[@{$_}],defined $spec{vnodes} ? @{$spec{vnodes}} : undef], ltolol(2,$vorio->pointlist);
1321				}
1322				}
1323
1324				$maxx = $pts[0]->[0]->[0];
1325				$minx = $pts[0]->[0]->[0];
1326				$maxy = $pts[0]->[0]->[1];
1327				$miny = $pts[0]->[0]->[1];
1328				foreach my $pt (@pts,@vpts) {
1329				if ($maxx < $pt->[0]->[0]) {$maxx = $pt->[0]->[0]}
1330				if ($maxy < $pt->[0]->[1]) {$maxy = $pt->[0]->[1]}
1331				if ($minx > $pt->[0]->[0]) {$minx = $pt->[0]->[0]}
1332				if ($miny > $pt->[0]->[1]) {$miny = $pt->[0]->[1]}
1333				}
1334
1335				# offset and scale to avoid limitations of svg renderers
1336
1337				my $dispsizex = '640';
1338				my $dispsizey = '480';
1339
1340				if (ref($dispsize) =~ /ARRAY/ && @{$dispsize} > 1) {
1341				$dispsizex = $dispsize->[0];
1342				$dispsizey = $dispsize->[1];
1343				}
1344
1345				# used to scale lines and point circle radii
1346				# so they stay visible in different viewports dimensions
1347				my $scale=(sqrt($dispsizex2+$dispsizey2))/sqrt(($maxx-$minx)2+($maxy-$miny)2);
1348
1349				foreach (@pts,@vpts) {
1350				$_->[0]->[0] -= $minx;
1351				$_->[0]->[0] *= $scale;
1352				$_->[0]->[1] -= $miny;
1353				$_->[0]->[1] *= $scale;
1354				}
1355
1356				my $scaled_maxx = ($maxx - $minx) * $scale;
1357				my $scaled_maxy = ($maxy - $miny) * $scale;
1358				my $scaled_minx = 0;
1359				my $scaled_miny = 0;
1360
1361				if (ref($triio) =~ /HASH/ && defined $triio->{nodes}) {
1362				if ($spec{edges}) {push @edges, map {[$_,@{$spec{edges}}]} map [$pts[$_->{nodes}->[0]->{index}]->[0],$pts[$_->{nodes}->[1]->{index}]->[0]], @{$triio->{edges}};}
1363				if ($spec{segments}) {push @segs, map {[$_,@{$spec{segments}}]} map [$pts[$_->{nodes}->[0]->{index}]->[0],$pts[$_->{nodes}->[1]->{index}]->[0]], @{$triio->{segments}};}
1364				#ignoring any subparametric points for elements
1365				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}};} #/
1366				}
1367				else {
1368				if ($spec{edges}) {push @edges, map {[[$pts[$_->[0]]->[0],$pts[$_->[1]]->[0]],@{$spec{edges}}]} ltolol(2,$triio->edgelist);}
1369				if ($spec{segments}) {push @segs, map {[[$pts[$_->[0]]->[0],$pts[$_->[1]]->[0]],@{$spec{segments}}]} ltolol(2,$triio->segmentlist);}
1370				#ignoring any subparametric points for elements
1371				if ($spec{elements}) {
1372				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); #/
1373				#read triangle attribute list, so at least those are available for choosing fill color in this case
1374				if (ref($spec{elements}->[0]) =~ /CODE/ && $triio->numberoftriangleattributes > 0 && $triio->numberofregions > 0) {
1375				my @eleattrs = ltolol($triio->numberoftriangleattributes,$triio->triangleattributes);
1376				for (my $i=0;$i<@elements;$i++) {
1377				# topologically-linked triangle not available here
1378				# but we'll fake it at least for the attributes list
1379				# so the color callback can still color according to region id
1380				# or whatever is in the triangle attribute list
1381				$elements[$i]->[1] = &{$spec{elements}->[0]}({attributes=>$eleattrs[$i]});
1382				}
1383				}
1384				}
1385				}
1386
1387				if ($vorio) {
1388				if (ref($vorio) =~ /HASH/ && defined $vorio->{nodes}) {
1389				# circles only available in this case
1390				if (defined $spec{circles}) {
1391				push @circles, map {
1392				[
1393				[ # this is point and radius: [x,y,r]
1394				@{$vpts[$_->{index}]->[0]},
1395				defined $_->{radius}
1396				? $_->{radius} * $scale
1397				: dist2d($vpts[$_->{index}]->[0],$edges[$_->{edges}->[0]->{index}]->[0]->[0])
1398				],
1399				@{$spec{circles}}
1400				]
1401				} @{$vorio->{nodes}};
1402				}
1403				if ($spec{vedges}) {
1404				@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}};
1405				}
1406				if (defined $spec{vrays}) {
1407				@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}};
1408				foreach my $ray (@vrays) {
1409				$ray->[0]->[1]->[0] *= $scale;
1410				$ray->[0]->[1]->[1] *= $scale;
1411				$ray->[0]->[1]->[0] += $ray->[0]->[0]->[0];
1412				$ray->[0]->[1]->[1] += $ray->[0]->[0]->[1];
1413				}
1414				}
1415				}
1416				else {
1417				if ($spec{vedges}) {
1418				my @ves =ltolol(2,$vorio->edgelist);
1419				my @vnorms=ltolol(2,$vorio->normlist);
1420				for (my $i=0;$i<@ves;$i++) {
1421				if ($ves[$i]->[0] > -1 && $ves[$i]->[1] > -1) {
1422				push @vedges, [[$vpts[$ves[$i]->[0]]->[0],$vpts[$ves[$i]->[1]]->[0]],@{$spec{vedges}}];
1423				}
1424				elsif (defined $spec{vrays}) {
1425				my $baseidx = ($ves[$i]->[0] != -1)?0:1;
1426				my $basept = $vpts[$ves[$i]->[$baseidx]]->[0];
1427				my $vec = $vnorms[$i];
1428				my $h = sqrt($vec->[0]2 + $vec->[1]2);
1429				$vec = [$vec->[0]/$h,$vec->[1]/$h];
1430				push @vrays, [
1431				[$basept,[$basept->[0]+$vec->[0]$maxx,$basept->[1]+$vec->[1]$maxx]],
1432				(defined($spec{vrays}) ? @{$spec{vrays}} : @{$spec{vedges}})];
1433				}
1434				}
1435				}
1436				if ($spec{circles}) {
1437				for (my $i=0;$i<@vpts;$i++) {
1438				push @circles, [[@{$vpts[$i]->[0]},dist2d($vpts[$i]->[0],$elements[$i]->[0])],@{$spec{circles}}];
1439				}
1440				}
1441				}
1442				}
1443
1444				my $margin_x_hi = 5;
1445				my $margin_x_lo = 5;
1446				my $margin_y_hi = 5;
1447				my $margin_y_lo = 5;
1448
1449				# extend margins to account for the radius of any circles or points
1450				my @round_things = (@circles, (map {[[$_->[0]->[0], $_->[0]->[1], $_->[2]]]} ( ($spec{nodes} ? @pts : () ), ($spec{vnodes} ? @vpts : () ))));
1451				if (scalar(@round_things) > 0) {
1452				my $cir_maxx = $round_things[0]->[0]->[0] + $round_things[0]->[0]->[2];
1453				my $cir_maxy = $round_things[0]->[0]->[1] + $round_things[0]->[0]->[2];
1454				my $cir_minx = $round_things[0]->[0]->[0] - $round_things[0]->[0]->[2];
1455				my $cir_miny = $round_things[0]->[0]->[1] - $round_things[0]->[0]->[2];
1456				foreach my $cir (@round_things) {
1457				if ($cir_maxx < $cir->[0]->[0] + $cir->[0]->[2]) {$cir_maxx = $cir->[0]->[0] + $cir->[0]->[2]}
1458				if ($cir_maxy < $cir->[0]->[1] + $cir->[0]->[2]) {$cir_maxy = $cir->[0]->[1] + $cir->[0]->[2]}
1459				if ($cir_minx > $cir->[0]->[0] - $cir->[0]->[2]) {$cir_minx = $cir->[0]->[0] - $cir->[0]->[2]}
1460				if ($cir_miny > $cir->[0]->[1] - $cir->[0]->[2]) {$cir_miny = $cir->[0]->[1] - $cir->[0]->[2]}
1461				}
1462				if ($cir_maxx-$scaled_maxx > $margin_x_hi) {$margin_x_hi = ($cir_maxx - $scaled_maxx) + 5;}
1463				if ($scaled_minx-$cir_minx > $margin_x_lo) {$margin_x_lo = ($scaled_minx - $cir_minx) + 5;}
1464				if ($cir_maxy-$scaled_maxy > $margin_y_hi) {$margin_y_hi = ($cir_maxy - $scaled_maxy) + 5;}
1465				if ($scaled_miny-$cir_miny > $margin_y_lo) {$margin_y_lo = ($scaled_miny - $cir_miny) + 5;}
1466				}
1467
1468				open(SVGO,'>'.$fn);
1469				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;
1470
1471
1472
1473
1474				EOS
1475
1476				if ($spec{elements}) {print SVGO "\n\n";}
1477				foreach my $ele (@elements) {
1478				print SVGO '',"\n"
1479				}
1480				if ($spec{edges}) {print SVGO "\n\n";}
1481				foreach my $edge (@edges) {
1482				print SVGO '',"\n"
1483				}
1484				if ($spec{segments}) {print SVGO "\n\n";}
1485				foreach my $edge (@segs) {
1486				print SVGO '',"\n"
1487				}
1488				if ($spec{vedges} && $vorio) {print SVGO "\n\n";}
1489				foreach my $edge (@vedges) {
1490				print SVGO '',"\n"
1491				}
1492				if ($spec{vrays} && $vorio) {print SVGO "\n\n";}
1493				foreach my $edge (@vrays) {
1494				print SVGO '',"\n"
1495				}
1496				if ($spec{nodes}) {
1497				print SVGO "\n\n";
1498				foreach my $pt (@pts) {
1499				print SVGO '',"\n"
1500				}
1501				}
1502				if ($spec{vnodes} && $vorio) {
1503				print SVGO "\n\n";
1504				foreach my $pt (@vpts) {
1505				print SVGO '',"\n"
1506				}
1507				}
1508				if ($spec{circles}) {print SVGO "\n\n";}
1509				foreach my $circle (@circles) {
1510				print SVGO '',"\n"
1511				}
1512
1513				if ($spec{raw}) {
1514				print SVGO "\n\n";
1515				print SVGO join("\n",map { s/ (c?x[12]?)="([0-9\.eE\-]+)"/' '.$1.'="'.(($2-$minx)*$scale).'"'/ge;
1516				s/ (c?y[12]?)="([0-9\.eE\-]+)"/' '.$1.'="'.(($2-$miny)*$scale).'"'/ge;
1517				$_;
1518				} @{$spec{raw}});
1519				}
1520				print SVGO "\n";
1521				close(SVGO);
1522				}
1523
1524				sub dist2d {sqrt(($_[0]->[0]-$_[1]->[0])2+($_[0]->[1]-$_[1]->[1])2)}
1525
1526				sub counterclockwise {
1527				my ($pa, $pb, $pc) = @_;
1528				return _counterclockwise($pa->[0],$pa->[1],$pb->[0],$pb->[1],$pc->[0],$pc->[1]);
1529				}
1530
1531				=head1 NAME
1532
1533				Math::Geometry::Delaunay - Quality Mesh Generator and Delaunay Triangulator
1534
1535				=head1 VERSION
1536
1537				Version 0.17
1538
1539				=cut
1540
1541				=head1 SYNOPSIS
1542
1543				=for html
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553				input vertices
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573				Delaunay triangulation
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592				Voronoi diagram
1593
1594
1595				use Math::Geometry::Delaunay qw(TRI_CCDT);
1596
1597				# generate Delaunay triangulation
1598				# and Voronoi diagram for a point set
1599
1600				my $point_set = [ [1,1], [7,1], [7,3],
1601				[3,3], [3,5], [1,5] ];
1602
1603				my $tri = new Math::Geometry::Delaunay();
1604				$tri->addPoints($point_set);
1605				$tri->doEdges(1);
1606				$tri->doVoronoi(1);
1607
1608				# called in void context
1609				$tri->triangulate();
1610				# populates the following lists
1611
1612				$tri->elements(); # triangles
1613				$tri->nodes(); # points
1614				$tri->edges(); # triangle edges
1615				$tri->vnodes(); # Voronoi diagram points
1616				$tri->vedges(); # Voronoi edges and rays
1617
1618
1619				=for html
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635				input PSLG
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
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				output mesh
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
1775
1776
1777
1778
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				something interesting extracted from topology
1808
1809
1810				# quality mesh of a planar straight line graph
1811				# with cross-referenced topological output
1812
1813				my $tri = new Math::Geometry::Delaunay();
1814				$tri->addPolygon($point_set);
1815				$tri->minimum_angle(23);
1816				$tri->doEdges(1);
1817
1818				# called in scalar context
1819				my $mesh_topology = $tri->triangulate(TRI_CCDT);
1820				# returns cross-referenced topology
1821
1822				# make two lists of triangles that touch boundary segments
1823
1824				my @tris_with_boundary_segment;
1825				my @tris_with_boundary_point;
1826
1827				foreach my $triangle (@{$mesh_topology->{elements}}) {
1828				my $nodes_on_boundary_count = (
1829				grep $_->{marker} == 1,
1830				@{$triangle->{nodes}}
1831				);
1832				if ($nodes_on_boundary_count == 2) {
1833				push @tris_with_boundary_segment, $triangle;
1834				}
1835				elsif ($nodes_on_boundary_count == 1) {
1836				push @tris_with_boundary_point, $triangle;
1837				}
1838				}
1839
1840
1841				=for html
1842
1843				=head1 DESCRIPTION
1844
1845				This is a Perl interface to the Jonathan Shewchuk's Triangle library.
1846
1847				"Triangle generates exact Delaunay triangulations, constrained Delaunay
1848				triangulations, conforming Delaunay triangulations, Voronoi diagrams, and
1849				high-quality triangular meshes. The latter can be generated with no small or
1850				large angles, and are thus suitable for finite element analysis."
1851				-- from L
1852
1853				=head1 EXPORTS
1854
1855				Triangle has several option switches that can be used in different combinations
1856				to choose a class of triangulation and then configure options within that class.
1857				To clarify the composition of option strings, or just to give you a head start,
1858				a few constants are supplied to configure different classes of mesh output.
1859
1860				TRI_CONSTRAINED = 'Y' for "Constrained Delaunay"
1861				TRI_CONFORMING = 'Dq0' for "Conforming Delaunay"
1862				TRI_CCDT = 'q' for "Constrained Conforming Delaunay"
1863				TRI_VORONOI = 'v' to generate the Voronoi diagram
1864
1865				For an illustration of these terms, see:
1866				L
1867
1868				=head1 CONSTRUCTOR
1869
1870				=head2 new
1871
1872				The constructor returns a Math::Geometry::Delaunay object.
1873
1874				my $tri = Math::Geometry::Delaunay->new();
1875
1876				=head1 MESH GENERATION
1877
1878				=head2 triangulate
1879
1880				Run the triangulation with specified options, and either populate the object's
1881				output lists, or return a hash reference giving access to a cross-referenced
1882				representation of the mesh topology.
1883
1884				Common options can be set prior to calling C. The full range of
1885				Triangle's options can also be passed to C as a string, or list
1886				of strings. For example:
1887
1888				my $tri = Math::Geometry::Delaunay->new('pzq0eQ');
1889
1890				my $tri = Math::Geometry::Delaunay->new(TRI_CCDT, 'q15', 'a3.5');
1891
1892				Triangle's command line switches are documented here:
1893				L
1894
1895				=head3 list output
1896
1897				After triangulate is invoked in void context, the output mesh data can be
1898				retrieved from the following methods, all of which return a reference to an
1899				array.
1900
1901				$tri->triangulate(); # void context - no return value requested
1902				# output lists now available
1903				$points = $tri->nodes(); # array of vertices
1904				$tris = $tri->elements(); # array of triangles
1905				$edges = $tri->edges(); # all the triangle edges
1906				$segs = $tri->segments(); # the PSLG segments
1907				$vpoints = $tri->vnodes(); # points in the voronoi diagram
1908				$vedges = $tri->vedges(); # edges in the voronoi diagram
1909
1910				Data may not be available for all lists, depending on which option switches were
1911				used. By default, nodes and elements are generated, while edges are not.
1912
1913				The members of the lists returned have these formats:
1914
1915				nodes: [x, y, < zero or more attributes >, < boundary marker >]
1916
1917				elements: [[x0, y0], [x1, y1], [x2, y2],
1918				< another three vertices, if "o2" switch used >,
1919				< zero or more attributes >
1920				]
1921				edges: [[x0, y0], [x1, y1], < boundary marker >]
1922
1923				segments: [[x0, y0], [x1, y1], < boundary marker >]
1924
1925				vnodes: [x, y, < zero or more attributes >]
1926
1927				vedges: [< vertex or vector >, < vertex or vector >, < ray flag >]
1928
1929				Boundary markers are 1 or 0. An edge or segment with only one end on a boundary
1930				has boundary marker 0.
1931
1932				The ray flag is 0 if the edge is not a ray, or 1 or 2, to indicate
1933				which vertex is actually a unit vector indicating the direction of the ray.
1934
1935				Import of the mesh data from the C data structures will be deferred until
1936				actually requested from the list fetching methods above. For speed and
1937				lower memory footprint, access only what you need, and consider suppressing
1938				output you don't need with option switches.
1939
1940				=head3 topological output
1941
1942				When triangulate is invoked in scalar or array context, it returns a hash ref
1943				containing the cross-referenced nodes, elements, edges, and PSLG segments of the
1944				triangulation. In array context, with the "v" switch enabled, the Voronoi
1945				topology is the second item returned.
1946
1947				my $topology = $tri->triangulate();
1948
1949				$topology now looks like this:
1950
1951				{
1952				nodes => [
1953				{ # a node
1954				point => [x0, x1],
1955				edges => [edgeref, ...],
1956				segments => [edgeref, ...], # a subset of edges
1957				elements => [elementref, ...],
1958				marker => 1 or 0 or undefined, # boundary marker
1959				attributes => [attr0, ...]
1960				},
1961				... more nodes like that
1962
1963				],
1964				elements => [
1965				{ # a triangle
1966				nodes => [noderef0, noderef1, noderef2],
1967				edges => [edgeref0, edgeref1, edgeref2],
1968				neighbors => [neighref0, neighref1, neighref2],
1969				attributes => [attrib0, ...]
1970				},
1971				... more triangles like that
1972				],
1973				edges => [
1974				{
1975				nodes => [noderef0, noderef1], # only one for a ray
1976				elements => [elemref0, elemref1], # one if on boundary
1977				vector => undefined or [x, y], # ray direction
1978				marker => 1 or 0 or undefined, # boundary marker
1979				index => # edge's index in edge list
1980				},
1981				... more edges like that
1982
1983				],
1984				segments => [
1985				{
1986				nodes => [noderef0, noderef1],
1987				elements => [elemref0, elemref1], # one if on boundary
1988				marker => 1 or 0 or undefined # boundary marker
1989				},
1990				... more segments
1991				]
1992				}
1993
1994				=head3 cross-referencing Delaunay and Voronoi
1995
1996				Corresponding Delaunay triangles and Voronoi nodes have the same index number
1997				in their respective lists.
1998
1999				In the topological output, any element in a triangulation has a record of its
2000				own index number that can by used to look up the corresponding node in the
2001				Voronoi diagram topology, or vice versa, like so:
2002
2003				($topo, $voronoi_topo) = $tri->triangulate('v');
2004
2005				# get a triangle reference where the index is not obvious
2006
2007				$element = $topo->{nodes}->[-1]->{elements}->[-1];
2008
2009				# this gets a reference to the corresponding node in the Voronoi diagram
2010
2011				$voronoi_node = $voronoi_topo->{nodes}->[$element->{index}];
2012
2013
2014				Corresponding edges in the Delaunay and Voronoi outputs have the same index
2015				number in their respective edge lists.
2016
2017				In the topological output, any edge in a triangulation has a record of its own
2018				index number that can by used to look up the corresponding edge in the Voronoi
2019				diagram topology, or vice versa, like so:
2020
2021				($topo, $voronoi_topo) = $tri->triangulate('ev');
2022
2023				# get an edge reference where it's not obvious what the edge's index is
2024
2025				$delaunay_edge = $topo->{nodes}->[-1]->{edges}->[-1];
2026
2027				# this gets a reference to the corresponding edge in the Voronoi diagram
2028
2029				$voronoi_edge = $voronoi_topo->{edges}->[$delaunay_edge->{index}];
2030
2031				=head1 METHODS TO SET SOME Triangle OPTIONS
2032
2033				=head2 area_constraint
2034
2035				Corresponds to the "a" switch.
2036
2037				With one argument, sets the maximum triangle area constraint for the
2038				triangulation. Returns the value supplied. With no argument, returns the
2039				current area constraint.
2040
2041				Passing -1 to C will disable the global area constraint.
2042
2043				=head2 minimum_angle
2044
2045				Corresponds to the "q" switch.
2046
2047				With one argument, sets the minimum angle allowed for triangles added in the
2048				triangulation. Returns the value supplied. With no argument, returns the
2049				current minimum angle constraint.
2050
2051				Passing -1 to C will cause the "q" switch to be omitted from
2052				the option string.
2053
2054				=head2 doEdges, doVoronoi, doNeighbors
2055
2056				These methods simply add or remove the corresponding letters from the
2057				option string. Pass in a true or false value to enable or disable.
2058				Invoke with no argument to read the current state.
2059
2060				=head2 quiet, verbose
2061
2062				Triangle prints a basic summary of the meshing operation to STDOUT unless
2063				the "Q" switch is present. This module includes the "Q" switch by default, but
2064				you can override this by passing a false value to C.
2065
2066				If you would like to see even more output regarding the triangulation process,
2067				there are are three levels of verbosity configurable with repeated "V"
2068				switches. Passing a number from 1 to 3 to the C method will enable
2069				the corresponding level of verbosity.
2070
2071				=head1 METHODS TO ADD VERTICES AND SEGMENTS
2072
2073				=head2 addVertices, addPoints
2074
2075				Takes a reference to an array of vertices, each vertex itself an reference to
2076				an array containing two coordinates and zero or more attributes. Attributes
2077				are floating point numbers.
2078
2079				# vertex format
2080				# [x, y, < zero or more attributes as floating point numbers >]
2081
2082				$tri->addPoints([[$x0, $y0], [$x1, $y1], ... ]);
2083
2084				Use addVertices to add vertices that are not part of a PSLG.
2085				Use addPoints to add points that are not part of a polygon or polyline.
2086				In other words, they do the same thing.
2087
2088				=head2 addSegments
2089
2090				Takes a reference to an array of segments.
2091
2092				# segment format
2093				# [[$x0, $y0], [$x1, $y1]]
2094
2095				$tri->addSegments([ $segment0, $segment1, ... ]);
2096
2097				If your segments are contiguous, it's better to use addPolyline, or addPolygon.
2098
2099				This method is provided because some point and polygon processing algorithms
2100				result in segments that represent polygons, but list the segments in a
2101				non-contiguous order, and have shared vertices repeated in each segment's record.
2102
2103				The segments added with this method will be checked for duplicate vertices, and
2104				references to these will be merged.
2105
2106				Triangle can handle duplicate vertices, but we would rather not feed them in on
2107				purpose.
2108
2109				=head2 addPolyline
2110
2111				Takes a reference to an array of vertices describing a curve.
2112				Creates PSLG segments for each pair of adjacent vertices. Adds the
2113				new segments and vertices to the triangulation input.
2114
2115				$tri->addPolyline([$vertex0, $vertex1, $vertex2, ...]);
2116
2117				=head2 addPolygon
2118
2119				Takes a reference to an array of vertices describing a polygon.
2120				Creates PSLG segments for each pair of adjacent vertices
2121				and creates and additional segment linking the last vertex to
2122				the first,to close the polygon. Adds the new segments and vertices
2123				to the triangulation input.
2124
2125				$tri->addPolygon([$vertex0, $vertex1, $vertex2, ...]);
2126
2127				=head2 addHole
2128
2129				Like addPolygon, but describing a hole or concavity - an area of the output mesh
2130				that should not be triangulated.
2131
2132				There are two ways to specify a hole. Either provide a list of vertices, like
2133				for addPolygon, or provide a single vertex that lies inside of a polygon, to
2134				identify that polygon as a hole.
2135
2136				# first way
2137				$tri->addHole([$vertex0, $vertex1, $vertex2, ...]);
2138
2139				# second way
2140				$tri->addPolygon( [ [0,0], [1,0], [1,1], [0,1] ] );
2141				$tri->addHole( [0.5,0.5] );
2142
2143				Hole marker points can also be used, in combination with the "c" option, to
2144				cause or preserve concavities in a boundary when Triangle would otherwise
2145				enclose a PSLG in a convex hull.
2146
2147				=head2 addRegion
2148
2149				Takes a polygon describing a region, and an attribute or area constraint. With
2150				both the "A" and "a" switches in effect, three arguments allow you to specify
2151				both an attribute and an optional area constraint.
2152
2153				The first argument may alternately be a single vertex that lies inside of
2154				another polygon, to identify that polygon as a region.
2155
2156				To be used in conjunction with the "A" and "a" switches.
2157
2158				# with the "A" switch
2159				$tri->addRegion(\@polygon, < attribute > );
2160
2161				# with the "a" switch
2162				$tri->addRegion(\@polygon, < area constraint > );
2163
2164				# with both "Aa"
2165				$tri->addRegion(\@polygon, < attribute >, < area constraint > );
2166
2167				If the "A" switch is used, each triangle generated within the bounds of a region
2168				will have that region's attribute added to the end of the triangle's
2169				attributes list, while each triangle not within a region will have a "0" added
2170				to the end of its attribute list.
2171
2172				If the "a" switch is used without a number following, each triangle generated
2173				within the bounds of a region will be subject to that region's area
2174				constraint.
2175
2176				If the "A" or "a" switches are not in effect, addRegion has the same effect as
2177				addPolygon.
2178
2179				=head1 METHODS TO ACCESS OUTPUT LISTS
2180
2181				The following methods retrieve the output lists after the triangulate method has
2182				been invoked in void context.
2183
2184				Triangle's output data is not imported from C to Perl until one of these methods
2185				is invoked, and then only what's needed to construct the list requested. So
2186				there may be a speed or memory advantage to accessing the output in this way -
2187				only what you need, when you need it.
2188
2189				The methods prefixed with "v" access the Voronoi diagram nodes and edges, if one
2190				was generated.
2191
2192				=head2 nodes
2193
2194				Returns a reference to a list of nodes (vertices or points).
2195
2196				my $pointlist = $tri->nodes(); # retrieve nodes/vertices/points
2197
2198				The nodes in the list have this structure:
2199
2200				[x, y, < zero or more attributes >, < boundary marker >]
2201
2202				=head2 elements
2203
2204				Returns a reference to a list of elements.
2205
2206				$triangles = $tri->elements(); # retrieve triangle list
2207
2208				The elements in the list have this structure:
2209
2210				[[x0, y0], [x1, y1], [x2, y2],
2211				< another three vertices, if "o2" switch used >
2212				< zero or more attributes >
2213				]
2214
2215				=head2 segments
2216
2217				Returns a reference to a list of segments.
2218
2219				$segs = $tri->segments(); # retrieve the PSLG segments
2220
2221				The segments in the list have this structure:
2222
2223				[[x0, y0], [x1, y1], < boundary marker >]
2224
2225				=head2 edges
2226
2227				Returns a reference to a list of edges.
2228
2229				$edges = $tri->edges(); # retrieve all the triangle edges
2230
2231				The edges in the list have this structure:
2232
2233				[[x0, y0], [x1, y1], < boundary marker >]
2234
2235				Note that the edge list is not produced by default. Request that it be generated
2236				by invoking C, or passing the 'e' switch to C.
2237
2238				=head2 vnodes
2239
2240				Returns a reference to a list of nodes in the Voronoi diagram.
2241
2242				$vpointlist = $tri->vnodes(); # retrieve Voronoi vertices
2243
2244				The Voronoi diagram nodes in the list have this structure:
2245
2246				[x, y, < zero or more attributes >]
2247
2248				=head2 vedges
2249
2250				Returns a reference to a list of edges in the Voronoi diagram. Some of these
2251				edges are actually rays.
2252
2253				$vedges = $tri->vedges(); # retrieve Voronoi diagram edges and rays
2254
2255				The Voronoi diagram edges in the list have this structure:
2256
2257				[< vertex or vector >, < vertex or vector >, < ray flag >]
2258
2259				If the edge is a true edge, the ray flag will be 0.
2260				If the edge is actually a ray, the ray flag will either be 1 or 2,
2261				to indicate whether the the first, or second vertex should be interpreted as
2262				a direction vector for the ray.
2263
2264				=head1 UTILITY FUNCTIONS
2265
2266				=head2 to_svg
2267
2268				This function is meant as a development and debugging aid, to "dump" the
2269				geometric data structures specific to this package to a graphical
2270				representation. Takes key-value pairs to specify topology hashes, output file,
2271				image dimensions, and styles for the elements in the various output lists.
2272
2273				The topology hash input for the C or C keys is just the hash
2274				returned by C. The value for the C key is a file name string.
2275				Omit C to print to STDOUT. For C, provide and array ref with width
2276				and height, in pixels. For output list styles, keys correspond to the output
2277				list names, and values consist of references to arrays containing style
2278				configurations, as demonstrated below.
2279
2280				Only geometry that has a style configuration will be displayed. The following
2281				example includes everything. To display a subset, just omit any of the style
2282				configuration key-value pairs.
2283
2284				($topo, $vtopo) = $tri->triangulate('ve');
2285
2286				to_svg( topo => $topo,
2287				vtopo => $vtopo,
2288
2289				file => "enchilada.svg", # omit for STDOUT
2290				size => [800, 600], # width, height in pixels
2291
2292				# line width or optional
2293				# svg color point radius extra CSS
2294
2295				nodes => ['black' , 0.3],
2296				edges => ['#CCCCCC', 0.7],
2297				segments => ['blue' , 0.9, 'stroke-dasharray:1 1;'],
2298				elements => ['pink'] , # string or callback; see below
2299
2300				# these require Voronoi input (vtopo)
2301
2302				vnodes => ['purple' , 0.3],
2303				vedges => ['#FF0000', 0.7],
2304				vrays => ['purple' , 0.6],
2305				circles => ['orange' , 0.6],
2306
2307				);
2308
2309				Note that for display purposes C does not include the infinite rays in
2310				the Voronoi diagram. To see the complete Voronoi diagram, including segments
2311				representing the infinite rays, you should include style configuration for the
2312				C key, as in the example above.
2313
2314				Elements (triangles) only need one style config entry, for color. (An optional
2315				second entry would be a string for additional CSS.) In this case,
2316				the first entry can also be a reference to a callback function. A reference to
2317				the triangle being processed for display will be passed to the callback
2318				function. Therefore the callback function can determine a color based on any
2319				features or relationships of that triangle.
2320
2321				Typically you might color each triangle according to the region it's in, by
2322				using Triangle's 'A' switch, and then reading the region attribute from the
2323				last item in the triangle's attribute list.
2324
2325				my $region_colors_callback = sub {
2326				my $tri_ref = shift;
2327				return ('gray','blue','green')[$tri_ref->{attributes}->[-1]];
2328				};
2329
2330				But any other data accessible through the triangle reference can be used to
2331				calculate a color. For instance, the triangle's three nodes can carry any
2332				number of attributes, which are interpolated during mesh generation. You
2333				might shade each triangle according to the average of a node attribute.
2334
2335				my $tri_nodes_average_callback = sub {
2336				my $tri_ref = shift;
2337				my $sum = 0;
2338				# calculate average of the eighth attribute in all nodes
2339				foreach my $node (@{$tri_ref->{nodes}}) {
2340				$sum += $node->{attributes}->[7];
2341				}
2342				return &attrib_val_to_grayscale_hexcode( $sum / 3 );
2343				};
2344
2345				=head2 mic_adjust
2346
2347				=for html
2348
2349
2350
2351
2352
2353
2354
2355
2356
2357
2358
2359
2360
2361
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				Voronoi edges (blue) as a poor medial axis approximation
2427
2428
2429
2430
2431
2432
2433
2434
2435
2436
2437
2438
2439
2440
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				improved approximation after calling mic_adust()
2506
2507
2508				Warning: not yet thoroughly tested; may move elsewhere
2509
2510				One use of the Voronoi diagram of a tessellated polygon is to derive an
2511				approximation of the polygon's medial axis by pruning infinite rays and perhaps
2512				trimming or refining remaining branches. The approximation improves as
2513				intervals between sample points on the polygon become shorter. But it's not
2514				always desirable to multiply the number of polygon points to achieve short
2515				intervals.
2516
2517				At any point on the true medial axis, there is a maximum inscribed circle,
2518				with it's center on the medial axis, and tangent to the polygon in at least
2519				two places.
2520
2521				The C function moves each Voronoi node so that it becomes the
2522				center of a circle that is tangent to the polygon at two points. In simple
2523				cases this is a maximum inscribed circle, and the point is on the medial axis.
2524				And when it's not, it still should be a much better approximation than the
2525				original point location. The radius to the tangent on the polygon is stored
2526				with the updated Voronoi node.
2527
2528				After calling C, the modified Voronoi topology can be used as a
2529				list of maximum inscribed circles, from which can be derive a straighter,
2530				better medial axis approximation, without having to increase the number of
2531				sample points on the polygon.
2532
2533				($topo, $voronoi_topo) = $tri->triangulate('e');
2534
2535				mic_adjust($topo, $voronoi_topo); # modifies $voronoi_topo in place
2536
2537				foreach my $node (@{$voronoi_topo->{nodes}}) {
2538				$mic_center = $node->{point};
2539				$mic_radius = $node->{radius};
2540				...
2541				}
2542
2543				Constructing a true medial axis is much more involved - a subject for a
2544				different module. Until that module appears, running topology through
2545				C and then walking and pruning the Voronoi topology might help
2546				fill the gap.
2547
2548				=head1 API STATUS
2549
2550				Currently Triangle's option strings are exposed to give more complete access to
2551				its features. More of these options, and perhaps certain common combinations of
2552				them, will likely be wrapped in method-call getter-setters. I would prefer to
2553				preserve the ability to use the option strings directly, but it may be better
2554				at some point to hide them completely behind a more descriptive interface.
2555
2556
2557				=head1 AUTHOR
2558
2559				Michael E. Sheldrake, C<< >>
2560
2561				Triangle's author is Jonathan Richard Shewchuk
2562
2563
2564				=head1 BUGS
2565
2566				Please report any bugs or feature requests to
2567
2568				C
2569
2570				or through the web interface at
2571
2572				L
2573
2574				I will be notified, and then you'll automatically be notified of progress on
2575				your bug as I make changes.
2576
2577
2578				=head1 SUPPORT
2579
2580				You can find documentation for this module with the perldoc command.
2581
2582				perldoc Math::Geometry::Delaunay
2583
2584
2585				You can also look for information at:
2586
2587				=over 4
2588
2589				=item * RT: CPAN's request tracker
2590
2591				L
2592
2593				=item * Search CPAN
2594
2595				L
2596
2597				=back
2598
2599
2600				=head1 ACKNOWLEDGEMENTS
2601
2602				Thanks go to Far Leaves Tea in Berkeley for providing oolongs and refuge,
2603				and a place for paths to intersect.
2604
2605
2606				=head1 LICENSE AND COPYRIGHT
2607
2608				Copyright 2013 Micheal E. Sheldrake.
2609
2610				This Perl binding to Triangle is free software;
2611				you can redistribute it and/or modify it under the terms of either:
2612				the GNU General Public License as published by the Free Software Foundation;
2613				or the Artistic License.
2614
2615				See http://dev.perl.org/licenses/ for more information.
2616
2617				=head2 Triangle license
2618
2619				B by Jonathan Richard Shewchuk, copyright 2005, includes the following
2620				notice in the C source code. Please refer to the C source, included in with this
2621				Perl module distribution, for the full notice.
2622
2623				This program may be freely redistributed under the condition that the
2624				copyright notices (including this entire header and the copyright
2625				notice printed when the `-h' switch is selected) are not removed, and
2626				no compensation is received. Private, research, and institutional
2627				use is free. You may distribute modified versions of this code UNDER
2628				THE CONDITION THAT THIS CODE AND ANY MODIFICATIONS MADE TO IT IN THE
2629				SAME FILE REMAIN UNDER COPYRIGHT OF THE ORIGINAL AUTHOR, BOTH SOURCE
2630				AND OBJECT CODE ARE MADE FREELY AVAILABLE WITHOUT CHARGE, AND CLEAR
2631				NOTICE IS GIVEN OF THE MODIFICATIONS. Distribution of this code as
2632				part of a commercial system is permissible ONLY BY DIRECT ARRANGEMENT
2633				WITH THE AUTHOR. (If you are not directly supplying this code to a
2634				customer, and you are instead telling them how they can obtain it for
2635				free, then you are not required to make any arrangement with me.)
2636
2637				=cut
2638
2639				1;