package CAD::Calc;
our $VERSION = '0.26';
use Math::Vec qw(
:terse
NewVec
);
use Math::Complex qw(acos);
use Math::Round::Var;
use Math::BigFloat;
# FIXME: add explicit exports to cleanup my namespace / make
# dependencies clearer
use vars qw(
$linear_precision
$angular_precision
$linr
$angr
$pi
);
$linear_precision = 1.0e-7;
$angular_precision = 1.0e-6;
$pi = atan2(1,1) * 4;
require Exporter;
@ISA='Exporter';
@EXPORT_OK = qw(
pi
distdivide
subdivide
shorten_line
dist
dist2d
line_vec
slope
segs_as_transform
chevron_to_ray
signdist
offset
shift_line
line_to_rectangle
isleft
howleft
iswithin
iswithinc
unitleft
unitright
unit_angle
angle_reduce
angle_parse
angle_quadrant
collinear
triangle_angles
intersection_data
line_intersection
seg_line_intersection
seg_seg_intersection
line_ray_intersection
seg_ray_intersection
ray_pgon_int_index
ray_pgon_closest_index
perp_through_point
foot_on_line
foot_on_segment
Determinant
pgon_as_segs
pgon_area
pgon_centroid
pgon_lengths
pgon_angles
pgon_deltas
pgon_direction
pgon_bisectors
sort_pgons_lr
pgon_start_index
re_order_pgon
order_pgon
stringify
stringify_line
pol_to_cart
cart_to_pol
print_line
point_avg
arc_2pt
);
use strict;
use Carp;
=pod
=head1 NAME
CAD::Calc - generic cad-related geometry calculations
=head1 AUTHOR
Eric L. Wilhelm <ewilhelm at cpan dot org>
http://scratchcomputing.com
=head1 COPYRIGHT
This module is copyright (C) 2004, 2005 Eric L. Wilhelm. Portions copyright
(C) 2003 by Eric L. Wilhelm and A. Zahner Co.
=head1 LICENSE
This module is distributed under the same terms as Perl. See the Perl
source package for details.
You may use this software under one of the following licenses:
(1) GNU General Public License
(found at http://www.gnu.org/copyleft/gpl.html)
(2) Artistic License
(found at http://www.perl.com/pub/language/misc/Artistic.html)
=head1 NO WARRANTY
This software is distributed with ABSOLUTELY NO WARRANTY. The author,
his former employer, and any other contributors will in no way be held
liable for any loss or damages resulting from its use.
=head1 Modifications
The source code of this module is made freely available and
distributable under the GPL or Artistic License. Modifications to and
use of this software must adhere to one of these licenses. Changes to
the code should be noted as such and this notification (as well as the
above copyright information) must remain intact on all copies of the
code.
Additionally, while the author is actively developing this code,
notification of any intended changes or extensions would be most helpful
in avoiding repeated work for all parties involved. Please contact the
author with any such development plans.
=cut
########################################################################
=head1 Configuration
Used to set package global values such as precision.
=head2 import
Not called directly. Triggered by the use() function.
import(%options, @EXPORT_TAGS);
Example:
use CAD::Calc (
-precision => 0.125,
-angular => 1.0e-6,
qw(
seg_seg_intersection
dist2d
print_line
)
);
=cut
sub import {
## print "import called with @_\n";
local @ARGV = @_; # shame that Getopt::Long isn't structured better!
use Getopt::Long ();
Getopt::Long::GetOptions( '-',
'precision=f' => \$linear_precision,
'angular=f' => \$angular_precision,
) or croak("bad import arguments");
## print "using $linear_precision for linear\n";
## print "using $angular_precision for angular\n";
$linr = Math::Round::Var->new($linear_precision);
$angr = Math::Round::Var->new($angular_precision);
## print "my linear rounding will be a ", ref($linr), "\n";
## print "my angular rounding will be a ", ref($angr), "\n";
CAD::Calc->export_to_level(1, @ARGV);
} # end subroutine import definition
########################################################################
=head1 Constants
=head2 pi
Returns the value of $CAD::Calc::pi
pi;
=cut
sub pi() {
return($pi);
} # end subroutine pi definition
########################################################################
=head1 Functions
These are all exported as options.
=head2 distdivide
Returns a list of point references resulting from dividing $line into
as many parts as possible which are at least $dist apart.
@points = distdivide(\@line, $dist);
=cut
sub distdivide {
my($line, $dist) = @_;
$dist or croak("call to distdivide would cause divide by zero");
my $ptA = NewVec(@{$line->[0]});
my $ptB = NewVec(@{$line->[1]});
my $seg = NewVec($ptB->Minus($ptA));
my $length = $seg->Length();
# optionally go for fewer points here?
my $count = $length / $dist;
$count = int($count);
return(subdivide($line, $count));
} # end subroutine distdivide definition
########################################################################
=head2 subdivide
Returns a list of point references resulting from subdividing $line
into $count parts. The list will be $count-1 items long, (does not
include $line->[0] and $line->[1]);
$line is of the form: [ [x1, y1, z1], [x2, y2, z2] ] where z1 and z2
are optional.
@points = subdivide($line, $count);
=cut
sub subdivide {
my ($line, $count) = @_;
$count || croak("cannot divide line into zero segments");
my $ptA = NewVec(@{$line->[0]});
my $ptB = NewVec(@{$line->[1]});
# print "line: @$ptA -- @$ptB\n";
my $seg = NewVec($ptB->Minus($ptA));
my @points;
for(my $st = 1; $st < $count; $st++) {
push(@points, [$ptA->Plus( [ $seg->ScalarMult($st / $count) ] ) ] );
}
return(@points);
} # end subroutine subdivide definition
########################################################################
=head2 shorten_line
Shortens the line by the distances given in $lead and $tail.
@line = shorten_line(\@line, $lead, $tail);
=cut
sub shorten_line {
my ($line, $lead, $tail) = @_;
my $ptA = NewVec(@{$line->[0]});
my $ptB = NewVec(@{$line->[1]});
# print "line: @$ptA -- @$ptB\n";
my $seg = NewVec($ptB->Minus($ptA));
my $len = $seg->Length();
($lead + $tail >= $len) && return();
# croak("CAD::Calc::shorten_line($lead, $tail)\n" .
# "\t creates inverted line from length: $len\n");
return(
[$ptA->Plus([$seg->ScalarMult($lead / $len)])],
[$ptB->Minus([$seg->ScalarMult($tail / $len)])],
);
} # end subroutine shorten_line definition
########################################################################
=head2 dist
Returns the direct distance from ptA to ptB.
dist($ptA, $ptB);
=cut
sub dist {
my($ptA, $ptB) = @_;
(ref($ptB) eq "ARRAY") || ($ptB = [0,0,0]);
my $dist = sqrt(
($ptB->[0] - $ptA->[0]) ** 2 +
($ptB->[1] - $ptA->[1]) ** 2 +
($ptB->[2] - $ptA->[2]) ** 2
);
return($dist);
} # end subroutine dist definition
########################################################################
=head2 dist2d
Purposefully ignores a z (2) coordinate.
dist2d($ptA, $ptB);
=cut
sub dist2d {
my($ptA, $ptB) = @_;
# print "ref is: ", ref($ptB), "\n";
# (ref($ptB) eq "ARRAY") || ($ptB = [0,0,0]);
# XXX why was this ^-- here?!
# print "ptB: @{$ptB}\n";
my $dist = sqrt(
($ptB->[0] - $ptA->[0]) ** 2 +
($ptB->[1] - $ptA->[1]) ** 2
);
return($dist);
} # end subroutine dist2d definition
########################################################################
=head2 line_vec
Returns a Math::Vec object representing the vector from $ptA to $ptB
(which is actually a segment.)
$vec = line_vec($ptA, $ptB);
=cut
sub line_vec {
return(NewVec(signdist(@_)));
} # end subroutine line_vec definition
########################################################################
=head2 slope
Calculates the 2D slope between points @ptA and @ptB. Slope is defined
as dy / dx (rise over run.)
If dx is 0, will return the string "inf", which Perl so kindly treats as
you would expect it to (except it doesn't like to answer the question
"what is infinity over infinity?") 5.8.? users: sorry, there seems to be
some regression here! (now we're using Math::BigFloat to return inf, so
rounding has to go through that)
$slope = slope(\@ptA, \@ptB);
=cut
sub slope {
my @line = @_;
my @delta = map({$line[1][$_] - $line[0][$_]} 0..1);
unless($delta[0]) {
if($delta[1] > 0) {
return(Math::BigFloat->binf());
}
else {
return(Math::BigFloat->binf('-'));
}
}
return($delta[1] / $delta[0]);
} # end subroutine slope definition
########################################################################
=head2 segs_as_transform
Allows two segments to specify transform data.
Returns: (\@translate, $rotate, $scale),
where:
@translate is a 2D array [$x, $y] basically describing segment @A
$rotate is the angular difference between $A[0]->$B[0] and $A[1]->$B[1]
$scale is the length of $A[1]->$B[1] divided by the length of
$A[0]->$B[0]
my ($translate, $rotate, $scale) = segs_as_transform(\@A, \@B);
=cut
sub segs_as_transform {
my ($A, $B) = @_;
my $av = line_vec(@$A);
# print_line($A);
my $sd = line_vec($A->[0], $B->[0]);
my $ed = line_vec($A->[1], $B->[1]);
my $ang = $ed->Ang() - $sd->Ang();
my $sl = $sd->Length();
$sl or croak("no length for divisor\n");
my $scale = $ed->Length() / $sl;
return([$av->[0], $av->[1]], $ang, $scale);
} # end subroutine segs_as_transform definition
########################################################################
=head2 chevron_to_ray
Converts a chevron into a directional line by finding the midpoint
between the midpoints of each edge and connecting to the middle point.
@line = chevron_to_ray(@pts);
=cut
sub chevron_to_ray {
my (@pts) = @_;
(scalar(@pts) == 3) or croak("chevron needs three points");
my @mids;
foreach my $seg (0,1) {
($mids[$seg]) = subdivide([$pts[$seg], $pts[$seg+1]], 2);
}
my ($start) = subdivide(\@mids, 2);
return($start, $pts[1]);
} # end subroutine chevron_to_ray definition
########################################################################
=head2 signdist
Returns the signed distance
signdist(\@ptA, \@ptB);
=cut
sub signdist {
my ($ptA, $ptB) = @_;
my $b = NewVec(@{$ptB});
return($b->Minus($ptA));
} # end subroutine signdist definition
########################################################################
=head2 offset
Creates a contour representing the offset of @polygon by $dist.
Positive distances are inward when @polygon is ccw.
@polygons = offset(\@polygon, $dist);
=cut
sub offset {
my ($polygon, $dist) = @_;
# this gets the OffsetPolygon routine (which still needs work)
my $helper = "Math::Geometry::Planar::Offset";
eval("require $helper;");
$@ and croak("cannot offset without $helper\n", $@);
$helper->import('OffsetPolygon');
my @pgons = OffsetPolygon($polygon, $dist);
return(@pgons);
} # end subroutine offset definition
########################################################################
=head2 intersection_data
Calculates the two numerators and the denominator which are required
for various (seg-seg, line-line, ray-ray, seg-ray, line-ray, line-seg)
intersection calculations.
($k, $l, $d) = intersection_data(\@line, \@line);
=cut
sub intersection_data {
my @l = @_;
my $n1 = Determinant(
$l[1][0][0]-$l[0][0][0],
$l[1][0][0]-$l[1][1][0],
$l[1][0][1]-$l[0][0][1],
$l[1][0][1]-$l[1][1][1],
);
my $n2 = Determinant(
$l[0][1][0]-$l[0][0][0],
$l[1][0][0]-$l[0][0][0],
$l[0][1][1]-$l[0][0][1],
$l[1][0][1]-$l[0][0][1],
);
my $d = Determinant(
$l[0][1][0]-$l[0][0][0],
$l[1][0][0]-$l[1][1][0],
$l[0][1][1]-$l[0][0][1],
$l[1][0][1]-$l[1][1][1],
);
return($n1, $n2, $d);
} # end subroutine intersection_data definition
########################################################################
=head2 line_intersection
Returns the intersection point of two lines.
@pt = line_intersection(\@line, \@line, $tolerance);
@pt or die "no intersection";
If tolerance is defined, it will be used to sprintf the parallel factor.
Beware of this, it is clunky and might change if I come up with
something better.
=cut
sub line_intersection {
my @l = (shift, shift);
my ($tol) = @_;
foreach my $should (0,1) {
# print "should have $should\n";
# print $l[$should], "\n";
(ref($l[$should]) eq "ARRAY") or warn "not good\n";
}
my ($n1, $n2, $d) = intersection_data(@l);
## print "d: $d\n";
if(defined($tol)) {
$d = sprintf("%0.${tol}f", $d);
}
if($d == 0) {
# print "parallel!\n";
return(); # parallel
}
my @pt = (
$l[0][0][0] + $n1 / $d * ($l[0][1][0] - $l[0][0][0]),
$l[0][0][1] + $n1 / $d * ($l[0][1][1] - $l[0][0][1]),
);
# print "got point: @pt\n";
return(@pt);
} # end subroutine line_intersection definition
########################################################################
=head2 seg_line_intersection
Finds the intersection of @segment and @line.
my @pt = seg_line_intersection(\@segment, \@line);
@pt or die "no intersection";
unless(defined($pt[1])) {
die "lines are parallel";
}
=cut
sub seg_line_intersection {
my (@l) = @_;
my ($n1, $n2, $d) = intersection_data(@l);
# XXX not consistent with line_intersection function
if(sprintf("%0.9f", $d) == 0) {
return(0); # lines are parallel
}
if( ! (($n1/$d <= 1) && ($n1/$d >=0)) ) {
return(); # no intersection on segment
}
my @pt = (
$l[0][0][0] + $n1 / $d * ($l[0][1][0] - $l[0][0][0]),
$l[0][0][1] + $n1 / $d * ($l[0][1][1] - $l[0][0][1]),
);
return(@pt);
} # end subroutine seg_line_intersection definition
########################################################################
=head2 seg_seg_intersection
my @pt = seg_seg_intersection(\@segmenta, \@segmentb);
=cut
sub seg_seg_intersection {
my (@l) = @_;
my ($n1, $n2, $d) = intersection_data(@l);
# print "data $n1, $n2, $d\n";
if(sprintf("%0.9f", $d) == 0) {
return(0); # lines are parallel
}
if( ! ((sprintf("%0.9f", $n1/$d) <= 1) && (sprintf("%0.9f", $n1/$d) >=0)) ) {
# warn("n1/d is ", $n1/$d);
return(); # no intersection on segment a
}
if( ! ((sprintf("%0.9f", $n2/$d) <= 1) && (sprintf("%0.9f", $n2/$d) >=0)) ) {
# warn("n2/d is ", $n2/$d);
return(); # no intersection on segment b
}
my @pt = (
$l[0][0][0] + $n1 / $d * ($l[0][1][0] - $l[0][0][0]),
$l[0][0][1] + $n1 / $d * ($l[0][1][1] - $l[0][0][1]),
);
return(@pt);
} # end subroutine seg_seg_intersection definition
########################################################################
=head2 line_ray_intersection
Intersects @line with @ray, where $ray[1] is the direction of the
infinite ray.
line_ray_intersection(\@line, \@ray);
=cut
sub line_ray_intersection {
my (@l) = @_;
my ($n1, $n2, $d) = intersection_data(@l);
# $n1 is distance along segment (must be between 0 and 1)
# $n2 is distance along ray (must be greater than 0)
if(sprintf("%0.9f", $d) == 0) {
# print "parallel\n";
return(0); # lines are parallel
}
# same as seg_ray_intersection(), but we skip the segment check
if($n2 / $d < 0) {
# print "nothing on ray\n";
# segment intersects behind ray
return();
}
my @pt = (
$l[0][0][0] + $n1 / $d * ($l[0][1][0] - $l[0][0][0]),
$l[0][0][1] + $n1 / $d * ($l[0][1][1] - $l[0][0][1]),
);
return(@pt);
} # end subroutine line_ray_intersection definition
########################################################################
=head2 seg_ray_intersection
Intersects @seg with @ray, where $ray[1] is the direction of the
infinite ray.
seg_ray_intersection(\@seg, \@ray);
=cut
sub seg_ray_intersection {
my (@l) = @_;
my ($n1, $n2, $d) = intersection_data(@l);
# $n1 is distance along segment (must be between 0 and 1)
# $n2 is distance along ray (must be greater than 0)
if(sprintf("%0.9f", $d) == 0) {
# print "parallel\n";
return(0); # lines are parallel
}
if( ! (($n1/$d <= 1) && ($n1/$d >=0)) ) {
# print "nothing on segment\n";
return(); # no intersection on segment
}
if($n2 / $d < 0) {
# print "nothing on ray\n";
# segment intersects behind ray
return();
}
my @pt = (
$l[0][0][0] + $n1 / $d * ($l[0][1][0] - $l[0][0][0]),
$l[0][0][1] + $n1 / $d * ($l[0][1][1] - $l[0][0][1]),
);
return(@pt);
} # end subroutine seg_ray_intersection definition
########################################################################
=head2 ray_pgon_int_index
Returns the first (lowest) index of @polygon which has a segment
intersected by @ray.
$index = ray_pgon_int_index(\@ray, \@polygon);
=cut
sub ray_pgon_int_index {
my ($ray, $pgon) = @_;
(scalar(@$ray) == 2) or croak("not a ray");
for(my $e = 0; $e < @$pgon; $e++) {
my $n = $e + 1;
($n > $#$pgon) && ($n -= @$pgon);
my $seg = [$pgon->[$e], $pgon->[$n]];
my @int = seg_ray_intersection($seg, $ray);
if(defined($int[1])) {
# print "intersect @int\n";
return($e);
}
}
return();
} # end subroutine ray_pgon_int_index definition
########################################################################
=head2 ray_pgon_closest_index
Returns the closest (according to dist2d) index of @polygon which has a
segment intersected by @ray.
$index = ray_pgon_closest_index(\@ray, \@polygon);
=cut
sub ray_pgon_closest_index {
my ($ray, $pgon) = @_;
(scalar(@$ray) == 2) or croak("not a ray");
my @found;
for(my $e = 0; $e < @$pgon; $e++) {
my $n = $e + 1;
($n > $#$pgon) && ($n -= @$pgon);
my $seg = [$pgon->[$e], $pgon->[$n]];
my @int = seg_ray_intersection($seg, $ray);
if(defined($int[1])) {
# print "intersect @int\n";
push(@found, [$e, dist2d($ray->[0], \@int)]);
}
}
if(@found) {
my $least = (sort({$a->[1] <=> $b->[1]} @found))[0];
return($least->[0]);
}
else {
return();
}
} # end subroutine ray_pgon_closest_index definition
########################################################################
=head2 perp_through_point
@line = perp_through_point(\@pt, \@line);
=cut
sub perp_through_point {
my ($pt, $seg) = @_;
my @nv = ( # normal vector:
$seg->[1][1] - $seg->[0][1],
- ($seg->[1][0] - $seg->[0][0]),
);
my @ep = ( # end point of ray
$pt->[0] + $nv[0],
$pt->[1] + $nv[1],
);
return($pt, \@ep);
} # end subroutine perp_through_point definition
########################################################################
=head2 foot_on_line
@pt = foot_on_line(\@pt, \@line);
=cut
sub foot_on_line {
my ($pt, $seg) = @_;
return(line_intersection($seg, [perp_through_point($pt, $seg)]));
} # end subroutine foot_on_line definition
########################################################################
=head2 foot_on_segment
Returns the perpendicular foot of @pt on @seg. See seg_ray_intersection.
@pt = foot_on_segment(\@pt, \@seg);
=cut
sub foot_on_segment {
my ($pt, $seg) = @_;
return(seg_line_intersection($seg, [perp_through_point($pt, $seg)]));
} # end subroutine foot_on_segment definition
########################################################################
=head2 Determinant
Determinant($x1, $y1, $x2, $y2);
=cut
sub Determinant {
my ($x1,$y1,$x2,$y2) = @_;
return($x1*$y2 - $x2*$y1);
} # end subroutine Determinant definition
########################################################################
=head2 pgon_as_segs
Returns a list of [[@ptA],[@ptB]] segments representing the edges of
@pgon, where segment "0" is from $pgon[0] to $pgon[1]
@segs = pgon_as_segs(@pgon);
=cut
sub pgon_as_segs {
my (@pgon) = @_;
my @segs = ([])x scalar(@pgon);
for(my $i = -1; $i < $#pgon; $i++) {
$segs[$i] = [@pgon[$i, $i+1]];
}
return(@segs);
} # end subroutine pgon_as_segs definition
########################################################################
=head2 pgon_area
Returns the area of @polygon. Returns a negative number for clockwise
polygons.
$area = pgon_area(@polygon);
=cut
sub pgon_area {
my @pgon = @_;
(@pgon > 2) or return();
my $atmp = 0;
for(my $i = -1; $i < $#pgon; $i++) {
my $j = $i+1;
my ($xi, $yi) = @{$pgon[$i]};
my ($xj, $yj) = @{$pgon[$j]};
$atmp += $xi * $yj - $xj * $yi;
}
$atmp or return(); # no area
return($atmp / 2);
} # end subroutine pgon_area definition
########################################################################
=head2 pgon_centroid
@centroid = pgon_centroid(@polygon);
=cut
sub pgon_centroid {
my @pgon = @_;
(@pgon > 2) or return();
my $atmp = 0;
my $xtmp = 0;
my $ytmp = 0;
for(my $i = -1; $i < $#pgon; $i++) {
my $j = $i+1;
my ($xi, $yi) = @{$pgon[$i]};
my ($xj, $yj) = @{$pgon[$j]};
my $ai = $xi * $yj - $xj * $yi;
$atmp += $ai;
$xtmp += ($xj + $xi) * $ai;
$ytmp += ($yj + $yi) * $ai;
}
$atmp or return(); # no area
my $area = $atmp / 2;
return($xtmp / (3 * $atmp), $ytmp / (3 * $atmp));
} # end subroutine pgon_centroid definition
########################################################################
=head2 pgon_lengths
@lengths = pgon_lengths(@pgon);
=cut
sub pgon_lengths {
my (@points) = @_;
my @lengths = (0) x scalar(@points);
for(my $i = -1; $i < $#points; $i++) {
$lengths[$i] = dist2d(@points[$i, $i+1]);
}
return(@lengths);
} # end subroutine pgon_lengths definition
########################################################################
=head2 pgon_angles
Returns the angle of each edge of polygon in xy plane. These fall
between -$pi and +$pi due to the fact that it is basically just a call
to the atan2() builtin.
Edges are numbered according to the index of the point which starts
the edge.
@angles = pgon_angles(@points);
=cut
sub pgon_angles {
my (@points) = @_;
my @angles = (0) x scalar(@points);
# print "number of angles: @angles\n";
for(my $i = -1; $i < $#points; $i++) {
my $vec = NewVec(signdist(@points[$i, $i+1]));
$angles[$i] = $vec->Ang();
}
return(@angles);
} # end subroutine pgon_angles definition
########################################################################
=head2 pgon_deltas
Returns the differences between the angles of each edge of @polygon.
These will be indexed according to the point at which they occur, and
will be positive radians for ccw angles. Summing the @deltas will yield
+/-2pi (negative for cw polygons.)
@deltas = pgon_deltas(@pgon);
=cut
sub pgon_deltas {
my (@pts) = @_;
my @angles = pgon_angles(@pts);
return(ang_deltas(@angles));
} # end subroutine pgon_deltas definition
########################################################################
=head2 ang_deltas
Returns the same thing as pgon_deltas, but saves a redundant call to
pgon_angles.
my @angs = pgon_angles(@pts);
my @dels = ang_deltas(@angs);
=cut
sub ang_deltas {
my (@angs) = @_;
my @deltas = (0)x $#angs;
for(my $i = 0; $i < @angs; $i++) {
my $ang = angle_reduce($angs[$i] - $angs[$i-1]);
$deltas[$i] = $ang;
}
return(@deltas);
} # end subroutine ang_deltas definition
########################################################################
=head2 pgon_direction
Returns 1 for counterclockwise and 0 for clockwise. Uses the sum of the
differences of angles of @polygon. If this sum is less than 0, the
polygon is clockwise.
$ang_sum = pgon_direction(@polygon);
=cut
sub pgon_direction {
my (@pgon) = @_;
my @angs = pgon_deltas(@pgon);
return(angs_direction(@angs));
} # end subroutine pgon_direction definition
########################################################################
=head2 angs_direction
Returns the same thing as pgon_direction, but saves a redundant call to
pgon_deltas.
my @angs = pgon_deltas(@pgon);
my $dir = angs_direction(@angs);
=cut
sub angs_direction {
my (@angs) = @_;
my $sum = 0;
foreach my $ang (@angs) {
$sum+= $ang;
}
return($sum > 0);
} # end subroutine angs_direction definition
########################################################################
=head2 pgon_bisectors
pgon_bisectors();
=cut
sub pgon_bisectors {
warn "unfinished";
croak "finish it";
} # end subroutine pgon_bisectors definition
########################################################################
=head2 sort_pgons_lr
Sorts polygons by their average points returning a list which reads from
left to right. (Rather odd place for this?)
@pgons = sort_pgons_lr(@pgons);
=cut
sub sort_pgons_lr {
my @pgons = @_;
# no sense calculating all for naught:
(scalar(@pgons) > 1) || return(@pgons);
my @avg;
foreach my $pgon (@pgons) {
push(@avg, [point_avg(@$pgon)]);
}
my @ord = sort({$avg[$a][0] <=> $avg[$b][0]} 0..$#avg);
return(@pgons[@ord]);
} # end subroutine sort_pgons_lr definition
########################################################################
=head2 pgon_start_index
Returns the index of pgon which is at the "lowest left".
$i = pgon_start_index(@pgon);
=cut
sub pgon_start_index {
my @pgon = @_;
my @ordered = sort({
my $c;
$c = $pgon[$a][$_] <=> $pgon[$b][$_] and return $c for 0..1;
} 0..$#pgon);
# print "index order @ordered\n";
return($ordered[0]);
} # end subroutine pgon_start_index definition
########################################################################
=head2 pgon_start_indexb
Returns the index of pgon which is at the "lowest left".
Different method (is it faster?)
$i = pgon_start_indexb(@pgon);
=cut
sub pgon_start_indexb {
my @pgon = @_;
my %sort_hash;
for(my $c = 0; $c < @pgon; $c++) {
# this is still janky, should really do something about slope?
my @pt = map({sprintf("%0.2f", $_)} @{$pgon[$c]});
$sort_hash{$pt[0]}{$pt[1]} = $c;
}
my ($least_x) = sort({$a <=> $b} keys(%sort_hash));
my ($least_y) = sort({$a <=> $b} keys(%{$sort_hash{$least_x}}));
my $index = $sort_hash{$least_x}{$least_y};
return($index);
} # end subroutine pgon_start_indexb definition
########################################################################
=head2 pgon_start_index_z
Yet another different method (is this even correct?)
pgon_start_index_z();
=cut
sub pgon_start_index_z {
my @pgon = @_;
# use the quarter-length
my @minmax = (sort({$a <=> $b} map({$_->[0]} @pgon)))[0,-1];
my $x_fourth = $minmax[0] + ($minmax[1] - $minmax[0]) / 4;
my @contend;
for(my $i = 0; $i < @pgon; $i++) {
if($pgon[$i][0] < $x_fourth) {
push(@contend, $i);
}
}
# print scalar(@contend), " contenders\n";
my @ordered = sort({$pgon[$a][1] <=> $pgon[$b][1]} @contend);
# to be even more thorough:
my @yminmax = (sort({$a <=> $b} map({$_->[1]} @pgon)))[0,-1];
# quantify with 1/4 of the yspan:
my $yspan = ($yminmax[1] - $yminmax[0]) / 4;
my $choice = shift(@ordered);
foreach my $idx (@ordered) {
# if it is below and left, then we already have it, if it above
# and left, then we might want it
if($pgon[$idx][1] < ($pgon[$choice][1] + $yspan)) {
($pgon[$idx][0] < $pgon[$choice][0]) and ($choice = $idx);
}
}
return($choice);
} # end subroutine pgon_start_index_z definition
########################################################################
=head2 re_order_pgon
Imposes counter-clockwise from "lower-left" ordering.
@pgon = re_order_pgon(@pgon);
=cut
sub re_order_pgon {
my @pgon = @_;
unless(pgon_direction(@pgon)) {
@pgon = reverse(@pgon);
}
my $index = pgon_start_index_z(@pgon);
return(order_pgon($index, \@pgon));
} # end subroutine re_order_pgon definition
########################################################################
=head2 order_pgon
Rewinds the polygon (e.g. list) to the specified $start index. This is
not restricted to polygons (just continuous (looped) lists.)
@pgon = order_pgon($start, \@pgon);
=cut
sub order_pgon {
my $index = shift;
my $pg = shift;
my @pgon = @{$pg};
($index < 0) and ($index += @pgon);
my @new;
for(my $d = 0; $d < @pgon; $d++) {
my $i = $index + $d;
($i > $#pgon) and ($i -= @pgon);
# print "using $i\n";
push(@new, $pgon[$i]);
}
return(@new);
} # end subroutine order_pgon definition
########################################################################
=head2 shift_line
Shifts line to right or left by $distance.
@line = shift_line(\@line, $distance, right|left);
=cut
sub shift_line {
my ($line, $dist, $dir) = @_;
my @line = @$line;
my $mvec;
if($dir eq "left") {
$mvec = unitleft(@line);
}
elsif($dir eq "right") {
$mvec = unitright(@line);
}
else {
croak ("direction must be \"left\" or \"right\"\n");
}
$mvec = NewVec($mvec->ScalarMult($dist));
my @newline = map({[$mvec->Plus($_)]} @line);
return(@newline);
} # end subroutine shift_line definition
########################################################################
=head2 line_to_rectangle
Creates a rectangle, centered about @line.
my @rec = line_to_rectangle(\@line, $offset, \%options);
The direction of the returned points will be counter-clockwise around
the original line, with the first point at the 'lower-left' (e.g. if
your line points up, $rec[0] will be below and to the left of
$line[0].)
Available options
ends => 1|0, # extend endpoints by $offset (default = 1)
=cut
sub line_to_rectangle {
my ($ln, $offset, $opts) = @_;
my %options = (ends => 1);
(ref($opts) eq "HASH") && (%options = %$opts);
my @line = @$ln;
($offset > 0) or
croak "offset ($offset) must be positive non-zero\n";
my $a = NewVec(@{$line[0]});
my $b = NewVec(@{$line[1]});
# unit vector of line
my $vec = NewVec(NewVec($b->Minus($a))->UnitVector());
# crossed with unit vector make unit vector left
my $perp = NewVec($vec->Cross([0,0,-1]));
my ($back, $forth);
if($options{ends}) {
$back = NewVec($a->Minus([$vec->ScalarMult($offset)]));
$forth = NewVec($b->Plus([$vec->ScalarMult($offset)]));
}
else {
$back = $a;
$forth = $b;
}
my $left = NewVec($perp->ScalarMult($offset));
my $right = NewVec($perp->ScalarMult(-$offset));
# upper and lower here only mean anything
# if line originally pointed "up"
my @ll = $back->Plus($left);
my @lr = $back->Plus($right);
my @ur = $forth->Plus($right);
my @ul = $forth->Plus($left);
return(\@ll, \@lr, \@ur, \@ul);
} # end subroutine line_to_rectangle definition
########################################################################
=head2 isleft
Returns true if @point is left of @line.
$bool = isleft(\@line, \@point);
=cut
sub isleft {
my ($line, $pt) = @_;
my $how = howleft($line, $pt);
return($how > 0);
} # end subroutine isleft definition
########################################################################
=head2 howleft
Returns positive if @point is left of @line.
$number = howleft(\@line, \@point);
=cut
sub howleft {
my ($line, $pt) = @_;
my $isleft = ($line->[1][0] - $line->[0][0]) *
($pt->[1] - $line->[0][1]) -
($line->[1][1] - $line->[0][1]) *
($pt->[0] - $line->[0][0]);
return($isleft);
} # end subroutine howleft definition
########################################################################
=head2 iswithin
Returns true if @pt is within the polygon @bound. This will be negative
for clockwise input.
$fact = iswithin(\@bound, \@pt);
=cut
sub iswithin {
my ($bnd, $pt) = @_;
my $winding = 0;
my @bound = @$bnd;
for(my $n = -1; $n < $#bound; $n ++) {
my $next = $n+1;
my @seg = ($bound[$n], $bound[$next]);
my $isleft = howleft(\@seg, $pt);
if($seg[0][1] <= $pt->[1]) {
if($seg[1][1] > $pt->[1]) {
($isleft > 0) && $winding++;
# print "winding up\n";
}
}
elsif($seg[1][1] <= $pt->[1]) {
($isleft < 0) && $winding--;
# print "winding up\n";
}
} # end for $n
# print "winding is $winding\n";
return($winding);
} # end subroutine iswithin definition
########################################################################
=head2 iswithinc
Seems to be consistently much faster than the typical winding-number
iswithin. The true return value is always positive regardless of the
polygon's direction.
$fact = iswithinc(\@bound, \@pt);
=cut
sub iswithinc {
my ($bnd, $pt) = @_;
my $c = 0;
my @bound = @$bnd;
my ($x, $y) = @$pt;
# straight from the comp.graphics.algorithms faq:
for (my $i = 0, my $j = $#bound; $i < @bound; $j = $i++) {
# print "checking from $j to $i\n";
(((($bound[$i][1]<=$y) && ($y<$bound[$j][1])) ||
(($bound[$j][1]<=$y) && ($y<$bound[$i][1]))) &&
($x < ($bound[$j][0] - $bound[$i][0]) * ($y - $bound[$i][1]) / ($bound[$j][1] - $bound[$i][1]) + $bound[$i][0]))
and ($c = !$c);
}
return($c);
} # end subroutine iswithinc definition
########################################################################
=head2 unitleft
Returns a unit vector which is perpendicular and to the left of @line.
Purposefully ignores any z-coordinates.
$vec = unitleft(@line);
=cut
sub unitleft {
my (@line) = @_;
my $ln = NewVec(
NewVec(@{$line[1]}[0,1])->Minus([@{$line[0]}[0,1]])
);
$ln = NewVec($ln->UnitVector());
my $left = NewVec($ln->Cross([0,0,-1]));
## my $isleft = isleft(\@line, [$left->Plus($line[0])]);
## print "fact said $isleft\n";
return($left);
} # end subroutine unitleft definition
########################################################################
=head2 unitright
Negative of unitleft().
$vec = unitright(@line);
=cut
sub unitright {
my $vec = unitleft(@_);
$vec = NewVec($vec->ScalarMult(-1));
return($vec);
} # end subroutine unitright definition
########################################################################
=head2 unit_angle
Returns a Math::Vec vector which has a length of one at angle $ang (in
the XY plane.) $ang is fed through angle_parse().
$vec = unit_angle($ang);
=cut
sub unit_angle {
my ($ang) = @_;
$ang = angle_parse($ang);
my $x = cos($ang);
my $y = sin($ang);
return(NewVec($x, $y));
} # end subroutine unit_angle definition
########################################################################
=head2 angle_reduce
Reduces $ang (in radians) to be between -pi and +pi.
$ang = angle_reduce($ang);
=cut
sub angle_reduce {
my $ang = shift;
while($ang > $pi) {
$ang -= 2*$pi;
}
while($ang <= -$pi) {
$ang += 2*$pi;
}
return($ang);
} # end subroutine angle_reduce definition
########################################################################
=head2 angle_parse
Parses the variable $ang and returns a variable in radians. To convert
degrees to radians: $rad = angle_parse($deg . "d")
$rad = angle_parse($ang);
=cut
sub angle_parse {
my $ang = shift;
if($ang =~ s/d$//) {
$ang *= $pi / 180;
}
return($ang);
} # end subroutine angle_parse definition
########################################################################
=head2 angle_quadrant
Returns the index of the quadrant which contains $angle. $angle is in
radians.
$q = angle_quadrant($angle);
@syms = qw(I II III IV);
print "angle is in quadrant: $syms[$q]\n";
=cut
sub angle_quadrant {
my $ang = shift;
my $x = cos($ang);
my $y = sin($ang);
my $vert = ($x < 0);
my $hori = ($y < 0);
my @list = (
[0,3],
[1,2],
);
return($list[$vert][$hori]);
} # end subroutine angle_quadrant definition
########################################################################
=head2 collinear
$fact = collinear(\@pt1, \@pt2, \@pt3);
=cut
sub collinear {
my @pts = @_;
(@pts == 3) or croak("must call with 3 points");
my ($pta, $ptb, $ptc) = @pts;
my $va = line_vec($pta, $ptb);
my $vb = line_vec($ptb, $ptc);
my $cp = NewVec($va->Cross($vb));
my $ta = $cp->Length();
# print "my vectors: @$va\n@$vb\n@$cp\n";
# print "angs: ", $va->Ang(), " and ", $vb->Ang(), "\n";
# print "ta: $ta\n";
return(abs($ta) < 0.001);
} # end subroutine collinear definition
########################################################################
=head2 triangle_angles
Calculates the angles of a triangle based on it's lengths.
@angles = triangle_angles(@lengths);
The order of the returned angle will be "the angle before the edge".
=cut
sub triangle_angles {
my @len = @_;
(@len == 3) or croak("triangle must have 3 sides");
my @angs = (
acos(
($len[2]**2 + $len[0]**2 - $len[1]**2) /
(2 * $len[2] * $len[0])
),
acos(
($len[1]**2 + $len[0]**2 - $len[2]**2) /
(2 * $len[1] * $len[0])
),
);
$angs[2] = $pi - $angs[0] - $angs[1];
print "angs: @angs\n";
} # end subroutine triangle_angles definition
########################################################################
=head2 stringify
Turns point into a string rounded according to $rnd. The optional
$count allows you to specify how many coordinates to use.
$string = stringify(\@pt, $rnd, $count);
=cut
sub stringify {
my ($pt, $rnd, $count) = @_;
# FIXME: # rounding should be able to do fancier things here:
unless(defined($count)) {
$count = scalar(@{$pt});
}
my $top = $count - 1;
my $str = join(",",
map( {sprintf("%0.${rnd}f", $_)} @{$pt}[0..$top]) );
return($str);
} # end subroutine stringify definition
########################################################################
=head2 stringify_line
Turns a line (or polyline) into a string. See stringify().
stringify_line(\@line, $char, $rnd, $count);
=cut
sub stringify_line {
my ($line, $char, $rnd, $count) = @_;
defined($char) or ($char = "\n");
defined($rnd) or ($rnd = 2);
defined($count) or ($count = 2);
return(join($char, map({stringify($_, $rnd, $count)} @$line)));
} # end subroutine stringify_line definition
########################################################################
=head2 pol_to_cart
Convert from polar to cartesian coordinates.
my ($x, $y, $z) = pol_to_cart($radius, $theta, $z);
=cut
sub pol_to_cart {
my ($r, $th, $z) = @_;
my $x = $r * cos($th);
my $y = $r * sin($th);
return($x, $y, $z);
} # end subroutine pol_to_cart definition
########################################################################
=head2 cart_to_pol
Convert from polar to cartesian coordinates.
my ($radius, $theta, $z) = cart_to_pol($x, $y, $z);
=cut
sub cart_to_pol {
my ($x, $y, $z) = @_;
my $r = sqrt($x**2 + $y**2);
my $th = atan2($y, $x);
return($r, $th, $z);
} # end subroutine cart_to_pol definition
########################################################################
=head2 print_line
print_line(\@line, $message);
=cut
sub print_line {
my ($line, $message) = @_;
unless($message) {
$message = "line:";
}
print join("\n\t", $message,
map({join(" ", @$_)} @$line)), "\n";
} # end subroutine print_line definition
########################################################################
=head2 point_avg
Averages the x and y coordinates of a list of points.
my ($x, $y) = point_avg(@points);
=cut
sub point_avg {
my(@points) = @_;
my $i;
my $num = scalar(@points);
my $x_avg = 0;
my $y_avg = 0;
# print "num is $num\n";
for($i = 0; $i < $num; $i++) {
# print "point: $points[$i][0]\n";
$x_avg += $points[$i][0];
$y_avg += $points[$i][1];
}
# print "avgs: $x_avg $y_avg\n";
$x_avg = $x_avg / $num;
$y_avg = $y_avg / $num;
return($x_avg, $y_avg);
} # end subroutine point_avg definition
=head2 arc_2pt
Given a pair of endpoints and an angle (in radians), returns an arc with
center, radius, and start/end angles.
my %arc = arc_2pt(\@pts, $angle);
=cut
sub arc_2pt {
my ($pts, $angle) = @_;
my $dir = (($angle >= 0) ? 1 : -1);
$angle = abs($angle);
my %arc;
my $chord = V(@{$pts->[1]}) - $pts->[0];
my $clen = abs($chord);
# warn "chord: $chord\n";
# warn "chord length: $clen\n";
my $eps = $angle /4;
(cos($eps) == 0) and die "ack";
my $blg = sin($eps)/cos($eps);
my $s = $clen / 2 * $blg;
my $r = (($clen/2)**2 + $s**2) / (2 * $s);
## warn "radius: $r\n";
## my $mid = $pts->[1] + $chord / 2;
my $gamma = (pi - $angle) / 2;
## warn "gamma: $gamma\n";
my $cang = $chord->Ang;
my $phi = $cang + $dir * $gamma;
## warn "phi: $phi\n";
my $conn = V(pol_to_cart($r, $phi));
my $center = $pts->[0] + $conn;
## warn "center: $center\n";
$arc{pt} = [@$center[0,1]];
$arc{rad} = $r;
$arc{angs} = [
(- $conn)->Ang,
($pts->[1] - $center)->Ang
];
($dir > 0) or ($arc{angs} = [reverse(@{$arc{angs}})]);
$arc{direction} = $dir;
return(%arc);
} # end subroutine arc_2pt definition
########################################################################
1;
syntax highlighted by Code2HTML, v. 0.9.1