CAD::Calc - generic cad-related geometry calculations

  (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)

  import(%options, @EXPORT_TAGS);

  use CAD::Calc (
  	-precision => 0.125,
	-angular   => 1.0e-6,
  	qw(
		seg_seg_intersection
		dist2d
		print_line
		)
	);

  pi;

  @points = distdivide(\@line, $dist);

  @points = subdivide($line, $count);

  @line = shorten_line(\@line, $lead, $tail);

  dist($ptA, $ptB);

  dist2d($ptA, $ptB);

  $vec = line_vec($ptA, $ptB);

  $slope = slope(\@ptA, \@ptB);

  my ($translate, $rotate, $scale) = segs_as_transform(\@A, \@B);

  @line = chevron_to_ray(@pts);

  signdist(\@ptA, \@ptB);

  @polygons = offset(\@polygon, $dist);

  ($k, $l, $d) = intersection_data(\@line, \@line);

  @pt = line_intersection(\@line, \@line, $tolerance);
  @pt or die "no intersection";

  my @pt = seg_line_intersection(\@segment, \@line);
  @pt or die "no intersection";
  unless(defined($pt[1])) {
    die "lines are parallel";
  }

  my @pt = seg_seg_intersection(\@segmenta, \@segmentb);

  line_ray_intersection(\@line, \@ray);

  seg_ray_intersection(\@seg, \@ray);

  $index = ray_pgon_int_index(\@ray, \@polygon);

  $index = ray_pgon_closest_index(\@ray, \@polygon);

  @line = perp_through_point(\@pt, \@line);

  @pt = foot_on_line(\@pt, \@line);

  @pt = foot_on_segment(\@pt, \@seg);

  Determinant($x1, $y1, $x2, $y2);

  @segs = pgon_as_segs(@pgon);

  $area = pgon_area(@polygon);

  @centroid = pgon_centroid(@polygon);

  @lengths = pgon_lengths(@pgon);

  @angles = pgon_angles(@points);

  @deltas = pgon_deltas(@pgon);

  my @angs = pgon_angles(@pts);
  my @dels = ang_deltas(@angs);

  $ang_sum = pgon_direction(@polygon);

  my @angs = pgon_deltas(@pgon);
  my $dir = angs_direction(@angs);

  pgon_bisectors();

  @pgons = sort_pgons_lr(@pgons);

  $i = pgon_start_index(@pgon);

  $i = pgon_start_indexb(@pgon);

  pgon_start_index_z();

  @pgon = re_order_pgon(@pgon);

  @pgon = order_pgon($start, \@pgon);

  @line = shift_line(\@line, $distance, right|left);

  my @rec = line_to_rectangle(\@line, $offset, \%options);

  ends => 1|0,   # extend endpoints by $offset (default = 1)

  $bool = isleft(\@line, \@point);

  $number = howleft(\@line, \@point);

  $fact = iswithin(\@bound, \@pt);

  $fact = iswithinc(\@bound, \@pt);

  $vec = unitleft(@line);

  $vec = unitright(@line);

  $vec = unit_angle($ang);

  $ang = angle_reduce($ang);

  $rad = angle_parse($ang);

  $q = angle_quadrant($angle);
  @syms = qw(I II III IV);
  print "angle is in quadrant: $syms[$q]\n";

  $fact = collinear(\@pt1, \@pt2, \@pt3);

  @angles = triangle_angles(@lengths);

  $string = stringify(\@pt, $rnd, $count);

  stringify_line(\@line, $char, $rnd, $count);

  my ($x, $y, $z) = pol_to_cart($radius, $theta, $z);

  my ($radius, $theta, $z) = cart_to_pol($x, $y, $z);

  print_line(\@line, $message);

	my ($x, $y) = point_avg(@points);

  my %arc = arc_2pt(\@pts, $angle);

package CAD::Calc;
our $VERSION = 0.27;

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;

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
########################################################################

sub pi() {
	return($pi);
} # end subroutine pi definition
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

sub line_vec {
	return(NewVec(signdist(@_)));
} # end subroutine line_vec definition
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

sub signdist {
	my ($ptA, $ptB) = @_;
	my $b = NewVec(@{$ptB});
	return($b->Minus($ptA));
} # end subroutine signdist definition
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

sub foot_on_line {
	my ($pt, $seg) = @_;
	return(line_intersection($seg, [perp_through_point($pt, $seg)]));
} # end subroutine foot_on_line definition
########################################################################

sub foot_on_segment {
	my ($pt, $seg) = @_;
	return(seg_line_intersection($seg, [perp_through_point($pt, $seg)]));
} # end subroutine foot_on_segment definition
########################################################################

sub Determinant {
	my ($x1,$y1,$x2,$y2) = @_;
	return($x1*$y2 - $x2*$y1);
} # end subroutine Determinant definition
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

sub pgon_deltas {
	my (@pts) = @_;
	my @angles = pgon_angles(@pts);
	return(ang_deltas(@angles));
} # end subroutine pgon_deltas definition
########################################################################

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
########################################################################

sub pgon_direction {
	my (@pgon) = @_;
	my @angs = pgon_deltas(@pgon);
	return(angs_direction(@angs));
} # end subroutine pgon_direction definition
########################################################################

sub angs_direction {
	my (@angs) = @_;
	my $sum = 0;
	foreach my $ang (@angs) {
		$sum+= $ang;
	}
	return($sum > 0);
} # end subroutine angs_direction definition
########################################################################

sub pgon_bisectors {
	warn "unfinished";
	croak "finish it";
} # end subroutine pgon_bisectors definition
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

sub isleft {
	my ($line, $pt) = @_;
	my $how = howleft($line, $pt);
	return($how > 0);
} # end subroutine isleft definition
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

sub unitright {
	my $vec = unitleft(@_);
	$vec = NewVec($vec->ScalarMult(-1));
	return($vec);
} # end subroutine unitright definition
########################################################################

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
########################################################################

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
########################################################################

sub angle_parse {
	my $ang = shift;
	if($ang =~ s/d$//) {
		$ang *= $pi / 180;
	}
	return($ang);
} # end subroutine angle_parse definition
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

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
########################################################################

sub print_line {
	my ($line, $message) = @_;
	unless($message) {
		$message = "line:";
	}
	print join("\n\t", $message, 
		map({join(" ", @$_)} @$line)), "\n";
} # end subroutine print_line definition
########################################################################

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

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;