File Coverage

blib/lib/Math/Geometry/Planar/Offset.pm
Criterion Covered Total %
statement 246 318 77.3
branch 61 106 57.5
condition 20 30 66.6
subroutine 7 7 100.0
pod 2 2 100.0
total 336 463 72.5


line stmt bran cond sub pod time code
1             package Math::Geometry::Planar::Offset;
2              
3             require 5.006;
4              
5             our $VERSION = '1.05';
6              
7 3     3   73418 use strict;
  3         7  
  3         115  
8 3     3   16 use warnings;
  3         4  
  3         75  
9              
10 3     3   16 use Carp;
  3         8  
  3         10830  
11              
12             our $debug = 0;
13             our $precision = 7;
14              
15             require Exporter;
16             our @ISA = qw(Exporter);
17             our @EXPORT = qw();
18             our @EXPORT_OK = qw(
19             OffsetPolygon
20             $precision
21             $debug
22             );
23              
24             =pod
25              
26             =head1 NAME
27              
28             Math::Geometry::Planar::Offset - Calculate offset polygons
29              
30             =head1 SYNOPSIS
31              
32             use Math::Geometry::Planar::Offset;
33              
34             my (@results) = OffsetPolygon(\@points, $distance);
35             foreach my $polygon (@results) {
36             # do something with @$polygon
37             }
38              
39             =head1 AUTHOR
40              
41             Eric Wilhelm @
42              
43             =head1 COPYRIGHT NOTICE
44              
45             Copyright (C) 2003-2007 Eric L. Wilhelm. All rights reserved.
46              
47             =head1 NO WARRANTY
48              
49             Absolutely, positively NO WARRANTY, neither express or implied, is
50             offered with this software. You use this software at your own risk. In
51             case of loss, neither Eric Wilhelm, nor anyone else, owes you anything
52             whatseover. You have been warned.
53              
54             Note that this includes NO GUARANTEE of MATHEMATICAL CORRECTNESS. If
55             you are going to use this code in a production environment, it is YOUR
56             RESPONSIBILITY to verify that the methods return the correct values.
57              
58             =head1 LICENSE
59              
60             You may use this software under one of the following licenses:
61              
62             (1) GNU General Public License
63             (found at http://www.gnu.org/copyleft/gpl.html)
64             (2) Artistic License
65             (found at http://www.perl.com/pub/language/misc/Artistic.html)
66              
67             =head1 BUGS
68              
69             There are currently some problems with concurrent edge events on outward
70             (and maybe inward) offsets. Some significant changes need to be made.
71              
72             =cut
73              
74             =head1 METHODS
75              
76             These methods are actually defined in Math::Geometry::Planar, which uses
77             this module.
78              
79             =head2 offset_polygon
80              
81             Returns reference to an array of polygons representing the original polygon
82             offsetted by $distance
83              
84             $polygon->offset_polygon($distance);
85              
86             =cut
87              
88             my $delta = 10 ** (-$precision);
89              
90             my $offset_depth = 0;
91             my $screen_height = 600;
92             my $flag;
93             my @bisectors;
94              
95             =head1 Functions
96              
97             Only OffsetPolygon is exported.
98              
99             =head2 pi
100              
101             Returns the constant pi
102              
103             =cut
104             sub pi {
105 148     148 1 318 atan2(1,1) * 4;
106             }
107              
108             =head2 OffsetPolygon
109              
110             Make offset polygon subroutine.
111              
112             Call with offset distance and ref to array of points for original
113             polygon polygon input must be pre-wrapped so point[n]=point[0]
114              
115             Will return a list of polygons (as refs.) The number of polygons in the
116             output depends on the shape of your input polygon. It may split into
117             several pieces during the offset.
118              
119             my (@results) = OffsetPolygon(\@points, $distance);
120             foreach my $polygon (@results) {
121             # do something with @$polygon
122             }
123              
124             =cut
125             sub OffsetPolygon {
126 6     6 1 4293 my ($points, $offset, $canvas) = @_;
127 6         10 my %intersects;
128 6         9 my ($count,$n,$n2);
129 0         0 my ($time, $time_static);
130             # my $key1,$key2;
131 0         0 my @angles;
132 0         0 my @directions;
133 0         0 my @phi;
134 0         0 my @bis_dir;
135 0         0 my @bis_end;
136 0         0 my @bis_scale;
137 0         0 my @result1;
138 0         0 my @result2;
139 0         0 my @outlist;
140 0         0 my ($cut1,$cut2,$cut);
141 6         21 $#outlist = 0;
142 6         9 my @newpoints;
143             my @newpoints1;
144              
145             # set limit on number of recursions
146             # does perl have its own limit to the number of scopes generated?
147 6 50       16 if($offset_depth > 100){
148 0         0 print "reached recursion depth_limit\n";
149 0         0 return();
150             }
151             #ew: $offset_depth would have had to be reset by calling script
152             # (now is incremented and decremented on each side of recursive call)
153 6         7 my $npoints = @{$points};
  6         9  
154 6 50       15 $debug && print "points: ", $npoints, "\n";
155 6 50       18 $debug && print "number of points $npoints\n";
156             # print "points: @{$points->[0]}\n";
157             # exit;
158             # my $bis_length = 15;
159             # All polygons should be counter-clockwise !!!
160 6         19 _find_direction($points,\@angles,\@directions);
161 6 50       13 $debug and print "starting recurse: ",$offset_depth -1,"\n";
162 6 50       13 $debug and print "offset: $offset\n";
163             # ew removed junk related to null points
164             # bail out if not at least a triangle
165 6 50       19 if($npoints < 3) {
166 0         0 carp("Need at least 3 non-colinear points to offset");
167 0         0 return;
168             }
169             # Now to make the bisectors
170 6         17 for($n = 0;$n < $npoints; $n++) {
171 37         48 $phi[$n] = ((pi) - $angles[$n]) / 2;
172 37         86 $bis_scale[$n] = abs(1 / sin($phi[$n]));
173 37         54 $bis_dir[$n] = $directions[$n] + $phi[$n];
174 37         130 $bis_end[$n][0] = $points->[$n][0] + $offset * $bis_scale[$n] * cos($bis_dir[$n]); # X-coordinate
175 37         112 $bis_end[$n][1] = $points->[$n][1] + $offset * $bis_scale[$n] * sin($bis_dir[$n]); # Y-coordinate
176             }
177              
178             # Draw bisectors
179 6 50       20 if ($canvas) {
180             # first delete existing bisectors
181 0         0 for($n = 0; $n<@bisectors;$n++) {
182 0         0 $canvas ->delete($bisectors[$n]);
183             }
184 0         0 @bisectors = ();
185             # then draw the bisector points from the vertices.
186 0         0 for($n = 0; $n < $npoints; $n++) {
187 0         0 my $x1 = $points->[$n][0];
188 0         0 my $y1 = $screen_height - $points->[$n][1];
189 0         0 my $x2 = $bis_end[$n][0];
190 0         0 my $y2 = $screen_height-$bis_end[$n][1];
191 0         0 $bisectors[$n] = $canvas->create('line', $x1,$y1,$x2,$y2, '-fill' => 'blue','-arrow'=>'last');
192             }
193             }
194            
195             # Need to find whether a bisector first intersects a bisector,
196             # or a bisector first intersects a "ghost" edge. Also, must know
197             # which one and at what time:
198 6         9 my ($Ax1,$Ay1,$Ax2,$Ay2,$Bx1,$By1,$Bx2,$By2);
199 0         0 my ($delAx,$delAy,$delBx,$delBy,$Am,$Ab,$Bm,$Bb);
200 0         0 my ($x_int,$y_int);
201 0         0 my ($first_join,$first_split,$split_seg);
202 6         11 my $first_event = "";
203 6         9 my $first_time = abs($offset)+1;
204 6         10 my $time_dir = $offset / abs($offset); # +1 or -1 to give direction
205 6         6 my $new_offset;
206 6         8 my $join_time = $first_time;
207 6         7 my $split_time = $first_time;
208            
209             # first find intersections of adjacent bisectors (if any)
210 6         16 for($n = 0;$n <$npoints;$n++) {
211 37         47 $Ax1 = $points->[$n][0];
212 37         43 $Ay1 = $points->[$n][1];
213 37         38 $Ax2 = $bis_end[$n][0];
214 37         39 $Ay2 = $bis_end[$n][1];
215 37         34 $delAx = $Ax2 - $Ax1;
216 37         31 $delAy = $Ay2 - $Ay1;
217             # elw: what amateur wrote this!
218 37 50       63 if($delAx) {
219 37         34 $Am = $delAy / $delAx; $Ab = $Ay1 - $Ax1 * $Am;
  37         42  
220             }
221             else {
222 0         0 $Am = "inf";
223             }
224             # Check against the next bisector:
225 37         50 $Bx1 = $points->[$n+1-$npoints][0];
226 37         47 $By1 = $points->[$n+1-$npoints][1];
227 37         595 $Bx2 = $bis_end[$n+1-$npoints][0];
228 37         577 $By2 = $bis_end[$n+1-$npoints][1];
229 37         35 $delBx = $Bx2 - $Bx1;
230 37         35 $delBy = $By2 - $By1;
231             # elw: maybe you are getting closer to zen now:)
232 37 50       48 if($delBx) {
233 37         39 $Bm = $delBy / $delBx;$Bb = $By1 - $Bx1 * $Bm;
  37         41  
234             }
235             else{
236 0         0 $Bm = "inf";
237             }
238 37 50 33     252 if( ($Am == $Bm) || ($Am eq $Bm) ) {
239 0         0 next;
240             }
241             # Calculate determinants to find if intersection is within the bisector segment.
242 37 100       67 if (! _do_cross($Ax1,$Ay1,$Ax2,$Ay2,$Bx1,$By1,$Bx2,$By2) ){
243 35         75 next;
244             }
245             # if we have an intersection of the skeleton and only a triangle,
246             # then it has collapsed to a point:
247 2 50       6 if($npoints < 4) {
248 0 0       0 $debug && print "collapsed\n";
249 0         0 return();
250             }
251 2 50       19 if($Am eq "inf") { # Slope of first line is infinite, so use Ax1.
    50          
252 0         0 $x_int = $Ax1; # Will always be on vertical line.
253 0         0 $y_int = $Bm * $x_int + $Bb;
254             }
255             elsif($Bm eq "inf") { # Slope of second line is infinite, so use Bx1.
256 0         0 $x_int = $Bx1;$y_int = $Am * $x_int + $Ab;
  0         0  
257             }
258             else {
259 2         3 $x_int = ($Bb - $Ab) / ($Am - $Bm);$y_int = $Am *$x_int + $Ab;
  2         3  
260             }
261             # Now, if the lines are not parallel, and the intersection
262             # happens on the line segments, we are here with
263             # the x and y coordinates of the intersection of the two lines.
264             ## $debug and printf ("intersection: %6.0f,%6.0f\n",$x_int, $y_int);
265              
266             # Let's find the time of intersect.
267             # distance formula with adjustment for speed
268 2         7 $time = sqrt( ($x_int - $points->[$n][0])**2 + ($y_int - $points->[$n][1])**2 ) / $bis_scale[$n];
269 2 100       7 if( (abs($time) < $first_time) )
270             # && ( $time / abs($time)== $offset / abs($offset) ) )
271             {
272             # note that none of the times loaded here have a sign
273 1         2 $first_time = abs($time);
274 1         1 $join_time = $first_time;
275 1         2 $first_join=$n;
276 1         3 $first_event="join";
277             }
278              
279              
280             } # end for n (for each bisector vs next neighbor)
281             # Time is smallest relevant offset distance before first join.
282             # first join is address of point to be joined.
283             # first event is "join" if controlling (could be over-ridden by split)
284             # If this is the controlling case, create the joined polygon at that offset time
285             # and make a recursive call.
286             #
287             #
288             # Now to check for intersections of bisectors with edges
289             # Have to check against all segments except adjacent ones.
290             # Nothing will happen if it is a triangle.
291 6 50       14 if($npoints > 3) {
292 6         14 for($n = 0; $n < $npoints; $n++) {
293 37         91 $Ax1 = $points->[$n][0]; # Get bisector endpoints
294 37         44 $Ay1 = $points->[$n][1];
295 37         42 $Ax2 = $bis_end[$n][0];
296 37         39 $Ay2 = $bis_end[$n][1];
297 37 50       63 $debug && print "starting at $Ax1 $Ay1\n";
298 37 50       53 $debug && print "bisector toward $Ax2 $Ay2\n";
299             # I need it to be a ray, so I am just making a really long segment here.
300 37         45 $delAx = 1000* $offset* ($Ax2 - $Ax1);
301 37         46 $delAy = 1000* $offset* ($Ay2 - $Ay1);
302 37         29 $Ax2 = $Ax1+ $delAx;
303 37         37 $Ay2 = $Ay1 + $delAy;
304 37 50       51 if($delAx) {$Am = $delAy / $delAx; $Ab = $Ay1 - $Ax1 * $Am;}else{$Am = "inf";};
  37         35  
  37         45  
  0         0  
305             # this loop has to compare to all of the polygon sides except the adjacent ones
306 37         79 for($n2 = 0; $n2 < $npoints; $n2++) {
307 249         265 $time = abs($offset) +1;
308             # Can't split adjacent segments
309 249 100 100     1141 if( ($n2 == $n) || ($n2 == $n-1) || ( ($n==0)&&($n2 == $npoints-1) ) ){
      100        
      66        
310 74         156 next;
311             }
312 175 50       286 $debug && print "inner loop at n2 = $n2\n";
313 175         209 $Bx1 = $points->[$n2][0]; # Get edge endpoints
314 175         184 $By1 = $points->[$n2][1];
315 175         260 $Bx2 = $points->[$n2 + 1 - $npoints][0];
316 175         214 $By2 = $points->[$n2+1-$npoints][1];
317 175         172 $delBx = $Bx2 - $Bx1;
318 175         149 $delBy = $By2 - $By1;
319 175 100       249 if($delBx) {
320 89         87 $Bm = $delBy / $delBx;$Bb = $By1 - $Bx1 * $Bm;
  89         106  
321             }
322             else {
323 86         113 $Bm = "inf";
324             }
325 175 50 33     1136 if( ($Am == $Bm) || ($Am eq $Bm) ){
326 0         0 next;
327             }
328             # Note that I am not using the dot product here (it
329             # shows up below. I need to know where the infinite ray
330             # of the bisector crosses the infinite line of the
331             # polygon edge (but this may cause problems as well)
332 175 50       643 if($Am eq "inf") {
    100          
333             # Slope of first line is infinite, so use Ax1.
334 0         0 $x_int = $Ax1;$y_int = $Bm * $x_int + $Bb;
  0         0  
335             }
336             elsif($Bm eq "inf") {
337             # Slope of second line is infinite, so use Bx1.
338 86         90 $x_int = $Bx1;$y_int = $Am * $x_int + $Ab;
  86         102  
339             }
340             else {
341 89         118 $x_int = ($Bb - $Ab) / ($Am - $Bm);$y_int = $Am *$x_int + $Ab;
  89         108  
342             }
343             # now make sure that the intersection happened on the right side
344             # (point the ray in the offset direction)
345 175         188 my $delxto_int = $x_int - $Ax1;
346 175         167 my $delyto_int = $y_int - $Ay1;
347 175         183 my $next = $n2+1;
348 175 100       367 if ($next == $npoints) {
349 25         29 $next = 0;
350             }
351             # This method is not handling the split which happens to edges terminating
352             # at a concave point. This event is a third case because the time calculated
353             # below will be smaller for the extension of the edge before the concavity than
354             # for the edge after the concavity, which is the one to properly split.
355              
356              
357             # I have suspicions about the accuracy here:
358             # Original intent was to make sure that the ray was properly treated
359             # now handled by using a really long segment (1000*offset) for the ray.
360             # if( $bis_dir[$n] == (atan2($delyto_int,$delxto_int)))
361             {
362             # direction is good, calculate distance
363             # note that the int_scale should be naturally signed
364             # to prevent strangeness on internal offsets.
365             # dvdp:
366             # $int_scale = 1 / (sin( $directions[$n2] - $bis_dir[$n])) ;
367             # This could lead to a devision by 0 so we calculate the reverse and check first
368 175         155 my $int_scale = (sin( $directions[$n2] - $bis_dir[$n])) ;
  175         263  
369 175 50       267 $debug && print "bisector scale: $bis_scale[$n]\n";
370 175 50       252 $debug && print "intersection scale: 1 / $int_scale\n";
371 175 50       210 if ($int_scale) {
372 175         184 $int_scale = 1 / $int_scale;
373 175 100       286 if($bis_scale[$n]+$int_scale) {
374 167         420 $time_static = sqrt( ($x_int - $points->[$n][0])**2 +
375             ($y_int - $points->[$n][1])**2 ) / ($bis_scale[$n]+$int_scale);
376             }
377             else {
378             # intersection happens only at infinite offset time
379             # die "PGON offset dying at 281 with $bis_scale[$n] and $int_scale\n";
380             }
381             }
382             else {
383             # 1 / $int_scale is arbitrary large so division reduces to 0
384 0         0 $time_static = 0;
385             }
386             # once you have calculated the time, you need to see if the
387             # intersected segment is still in that place at that time to be struck by the ray
388             # Does time_dir belong here? (yes for test case A)
389 175         299 $Bx1 = $points->[$n2][0] + $time_static * $bis_scale[$n2] * cos($bis_dir[$n2]);
390 175         259 $By1 = $points->[$n2][1] + $time_static * $bis_scale[$n2] * sin($bis_dir[$n2]);
391 175         302 $Bx2 = $points->[$n2 + 1 - $npoints][0] + $time_static * $bis_scale[$next] * cos($bis_dir[$next]);
392 175         277 $By2 = $points->[$n2 + 1 - $npoints][1] + $time_static * $bis_scale[$next] * sin($bis_dir[$next]);
393 175 100       278 if(! _do_cross($Ax1,$Ay1,$Ax2,$Ay2,$Bx1,$By1,$Bx2,$By2) ) {
394 125         402 next;
395             }
396 50         57 $time = $time_static;
397             # make sure the case controls and is in the right direction. (required for case A)
398             # dvdp: replaced division by abs($time) by a multiplication to overcome
399             # potential divide by 0
400 50 100 66     196 if( (abs($time)<$first_time) and ( $time == $time_dir * abs($time) ) ) {
401             # note that none of the times loaded here have a sign
402 4         5 $first_time = abs($time);
403 4         13 $split_time = $first_time;
404 4         39 $first_split=$n;
405 4         6 $split_seg = $n2;
406 4         15 $first_event="split";
407             }
408             }
409             } # end for n2
410             } # end for n (each bisector verses each edge)
411             } # End if not triangle
412            
413             # Time is the smallest relavent split time (unless join controls)
414             # first_split is address of point to be split.
415             # first_event is "split" if controlling (could be controlled by join)(or nothing)
416             # If this is the controlling case, create the split polygons at that offset time
417             # and make a recursive call.
418            
419 6 100       18 if($first_time <= abs($offset)) {
420 2 50       4 if($debug) {
421 0 0       0 if($first_event eq "join") {
422 0         0 printf("join vertex %d at time %6.2f\n",$first_join,$first_time);
423             }
424 0 0       0 if($first_event eq "split") {
425 0         0 printf(
426             "split segment %d with vertex %d at time %6.2f\n",
427             $split_seg,$first_split,$first_time
428             );
429             }
430 0         0 print "\ttime: $first_time\n";
431             }
432             # get a list of points for the offset polygon at first_time
433              
434 2 100       8 if($first_event eq "join") {
    50          
435             # How are coincident events handled?
436             # remove the offending point and call yourself with time adjusted.
437 1         2 @newpoints = ();
438 1         4 for($n = 0; $n<$npoints; $n++) {
439 5         13 $newpoints[$n][0] = $points->[$n][0] + $first_time *
440             $time_dir* $bis_scale[$n] * cos($bis_dir[$n]);
441 5         21 $newpoints[$n][1] = $points->[$n][1] + $first_time *
442             $time_dir* $bis_scale[$n] * sin($bis_dir[$n]);
443             }
444             # keep the one with the smaller scale.
445 1 50       3 if($bis_scale[$first_join]<$bis_scale[$first_join+1]) {
446 1         2 splice(@newpoints, $first_join + 1, 1);
447             }
448             else {
449 0         0 splice(@newpoints, $first_join , 1);
450             }
451 1         2 $npoints--;
452 1         2 $new_offset = $offset - $first_time * $time_dir; # put a direction on it.
453 1         1 $offset_depth++;
454 1         26 @outlist = OffsetPolygon(\@newpoints,$new_offset) ;
455 1         3 $offset_depth--;
456 1         8 return(@outlist);
457             }
458             elsif ($first_event eq "split") {
459             # make two polygons and call yourself with time adjusted.
460 1         2 @newpoints = ();
461 1         1 my @split_check;
462             my @split_check1;
463 1         4 for($n = 0; $n<$npoints; $n++) {
464 8         24 $newpoints[$n][0] = $points->[$n][0] + $first_time *
465             $time_dir * $bis_scale[$n] * cos($bis_dir[$n]);
466 8         20 $newpoints[$n][1] = $points->[$n][1] + $first_time *
467             $time_dir * $bis_scale[$n] * sin($bis_dir[$n]);
468 8         19 $split_check[$n] = $n;
469             }
470 1         2 $new_offset = $offset - $first_time * $time_dir;
471 1         1 @newpoints1 = ();
472 1         3 @newpoints1 = @newpoints;
473 1         1 @split_check1 = ();
474 1         2 @split_check1 = @split_check;
475             # have first_split and split_seg, must find a way to wrap around.
476 1         2 my $cut_startA;
477             my $cut_startB;
478 0         0 my $cutB;
479 1 50       2 if($split_seg < $first_split) {
480 0         0 $cut1 = $#newpoints - $first_split;
481 0         0 $cut2 = $split_seg+1;
482 0         0 $cut_startA=$first_split+1;
483 0         0 $cut_startB = $split_seg+1;
484 0         0 $cutB = $first_split - $split_seg -1;
485             }
486             else {
487 1         2 $cut1 = $#newpoints - $split_seg;
488 1         2 $cut2 = $first_split;
489 1         1 $cut_startA = $split_seg +1;
490 1         1 $cut_startB = $first_split + 1;
491 1         2 $cutB = $split_seg - $first_split;
492             }
493 1         2 splice(@newpoints,$cut_startA,$cut1);
494 1         2 splice(@newpoints,0,$cut2);
495 1         11 splice(@split_check,$cut_startA,$cut1);
496 1         1 splice(@split_check,0,$cut2);
497 1         2 splice(@newpoints1,$cut_startB,$cutB);
498 1         2 splice(@split_check1,$cut_startB,$cutB);
499 1         1 my $num = @newpoints;
500 1         2 my $num1 = @newpoints1;
501              
502 1 50       3 if($debug) {
503 0         0 print "split was:\n@split_check\n@split_check1\n";
504 0         0 print "first splitted:\n";
505 0         0 for($n = 0; $n<$num;$n++) {
506 0         0 printf ("point: %d (%6.2f,%6.2f)\n",$n,@{$newpoints[$n]});
  0         0  
507             }
508 0         0 print "second splitted:\n";
509 0         0 for($n = 0; $n < $num1;$n++) {
510 0         0 printf ("point: %d (%6.2f,%6.2f)\n",$n,@{$newpoints1[$n]});
  0         0  
511             }
512             }
513 1         3 $#result1 = 0;
514 1         2 shift(@result1);
515 1         2 $#result2 = 0;
516 1         1 shift(@result2);
517 1 50       2 if($num>2) {
518 1         2 $offset_depth++;
519 1         34 @result1 = OffsetPolygon(\@newpoints,$new_offset);
520 1         2 $offset_depth--;
521             }
522             else {
523 0 0       0 $debug and print "too few on first\n";
524             }
525 1 50       5 if($num1>2) {
526 1         1 $offset_depth++;
527 1         2 @result2 = OffsetPolygon(\@newpoints1,$new_offset);
528 1         2 $offset_depth--;
529             }
530             else {
531 0 0       0 $debug and print "too few on second\n";
532             }
533 1         14 return(@result1, @result2); # This concatenates the list.
534             } # end elsif split
535             } # end if time <= offset
536             else {
537             # No splits or joins needed, so offset the thing.
538             # Bisector endpoints should be available.
539 4         7 @newpoints = ();
540 4         57 for($n = 0; $n<$npoints; $n++) {
541 24         82 $newpoints[$n][0] = $points->[$n][0] + $offset * $bis_scale[$n] * cos($bis_dir[$n]);
542 24         79 $newpoints[$n][1] = $points->[$n][1] + $offset * $bis_scale[$n] * sin($bis_dir[$n]);
543             }
544 4         49 return( \@newpoints);
545             }
546             } # End offset polygon subroutine definition
547             #########################################################################
548             # Find direction subroutine.
549             # Called with (,, ).
550             # Returns total change in angle (as radians (everything is as radians)).
551             # Second two array refs are optional.
552             # Forcing counter-clockwise currently not done here, but maybe should.
553             sub _find_direction {
554 6     6   11 my($coords,$del_theta,$thetas) = @_;
555 6         6 my $sum_theta=0;
556 6         9 my $sum_delta_theta = 0;
557 6         8 @{$del_theta} = (); # Clear arrays referenced for in-place modification.
  6         11  
558 6         7 @{$thetas} = ();
  6         8  
559 6         10 my $n = 0;
560 6         8 for(;;) {
561 43 100       40 last if ($n == @{$coords});
  43         83  
562             # dvdp: take advantage of negative indexing in Perl
563 37         35 my $n2 = $n - @{$coords} + 1;
  37         46  
564 37         56 my $n3 = $n - 1;
565 37         44 my $Ax1 = ${$coords}[$n][0]; # Line leaving point
  37         51  
566 37         107 my $Ay1 = ${$coords}[$n][1]; # this is the one belonging to the point
  37         44  
567 37         32 my $Ax2 = ${$coords}[$n2][0];
  37         47  
568 37         35 my $Ay2 = ${$coords}[$n2][1];
  37         43  
569 37         40 my $delAx = $Ax2 - $Ax1;
570 37         38 my $delAy = $Ay2 - $Ay1;
571 37         30 my $Bx1 = ${$coords}[$n3][0]; # Line coming to point n
  37         88  
572 37         32 my $By1 = ${$coords}[$n3][1];
  37         46  
573 37         35 my $Bx2 = ${$coords}[$n][0];
  37         41  
574 37         31 my $By2 = ${$coords}[$n][1];
  37         43  
575 37         44 my $delBx = $Bx2 - $Bx1;
576 37         34 my $delBy = $By2 - $By1;
577             # dvdp remove coinciding points
578 37 50 66     223 if ( ((abs($delAx) < $delta) && (abs($delAy) < $delta)) ||
      66        
      33        
579             ((abs($delBx) < $delta) && (abs($delBy) < $delta)) ) {
580 0         0 splice @{$coords},$n,1;
  0         0  
581 0         0 next;
582             }
583 37         57 my $theta_A = atan2($delAy,$delAx); # Angle of leaving line
584 37         58 my $theta_B = atan2($delBy,$delBx); # Angle of coming line
585 37         41 my $theta_AB = $theta_A - $theta_B;
586 37 100       56 if ($theta_AB < -(pi)) {
    50          
587 6         9 $theta_AB += 2 * (pi);
588             }
589             elsif ($theta_AB > (pi)) {
590 0         0 $theta_AB -= 2 * (pi);
591             }
592 37 50       63 if( (abs($theta_AB - (pi)) < $delta) ) {
593 0         0 splice @{$coords},$n,1;
  0         0  
594 0 0       0 $n and $n--; # need to recalc thetas if spike removed !
595 0         0 next;
596             }
597 37         37 ${$thetas}[$n] = $theta_A;
  37         56  
598 37         41 ${$del_theta}[$n] = $theta_AB;
  37         47  
599 37         47 $n++
600             }
601             } # End _find_direction sub
602             #########################################################################
603             # Determines if a pair of line segments have an intersection.
604             sub _do_cross {
605 212     212   314 my ($Ax1,$Ay1,$Ax2,$Ay2,$Bx1,$By1,$Bx2,$By2) = @_;
606             # Calculate four relevant minor determinants
607 212         308 my $det_123=($Ax2 - $Ax1)*($By1 - $Ay1) - ($Bx1 - $Ax1)*($Ay2 - $Ay1);
608 212         292 my $det_124=($Ax2 - $Ax1)*($By2 - $Ay1) - ($Bx2 - $Ax1)*($Ay2 - $Ay1);
609 212         262 my $det_341=($Bx1 - $Ax1)*($By2 - $Ay1) - ($Bx2 - $Ax1)*($By1 - $Ay1);
610 212         227 my $det_342 =$det_123-$det_124+$det_341;
611 212 100 100     648 if( ($det_123*$det_124 > 0) || ($det_341*$det_342 > 0) ) {
612 160         373 return(0); # segments only intersect if the above two products are both negative
613             }
614             else {
615 52         121 return(1);
616             }
617             } # End _do_cross subroutine definition
618             ########################################################################
619              
620             1;