File Coverage

blib/lib/Math/Vector/Real.pm

Criterion	Covered	Total	%
statement	103	444	23.2
branch	19	164	11.5
condition	4	20	20.0
subroutine	23	67	34.3
pod	33	59	55.9
total	182	754	24.1

package Math::Vector::Real;

our $VERSION = '0.18';

use strict;
use warnings;
use Carp;
use POSIX ();

use Exporter qw(import);
our @EXPORT = qw(V);

our $dont_use_XS;
unless ($dont_use_XS) {
    my $xs_version = do {
	local ($@, $!, $SIG{__DIE__});
	eval {
	    require Math::Vector::Real::XS;
	    $Math::Vector::Real::XS::VERSION;
	}
    };

    if (defined $xs_version and $xs_version < 0.07) {
	croak "Old and buggy version of Math::Vector::Real::XS detected, update it!";
    }
}


our %op = (add => '+',
	   neg => 'neg',
	   sub => '-',
	   mul => '*',
	   div => '/',
	   cross => 'x',
	   add_me => '+=',
	   sub_me => '-=',
	   mul_me => '*=',
	   div_me => '/=',
	   abs => 'abs',
	   atan2 => 'atan2',
	   equal => '==',
	   nequal => '!=',
	   clone => '=',
	   as_string => '""');

our %ol;
$ol{$op{$_}} = \&{${Math::Vector::Real::}{$_}} for keys %op;

require overload;
overload->import(%ol);

sub V { bless [@_] }

sub new {
    my $class = shift;
    bless [@_], $class
}

sub new_ref {
    my $class = shift;
    bless [@{shift()}], $class;
}

sub zero {
    my ($class, $dim) = @_;
    $dim >= 0 or croak "negative dimension";
    bless [(0) x $dim], $class
}

sub is_zero {
    $_ and return 0 for @$_[0];
    return 1
}

sub cube {
    my ($class, $dim, $size) = @_;
    bless [($size) x $dim], $class;
}

sub axis_versor {
    my ($class, $dim, $ix);
    if (ref $_[0]) {
        my ($self, $ix) = @_;
        $class = ref $self;
        $dim = @$self;
    }
    else {
        ($class, $dim, $ix) = @_;
        $dim >= 0 or croak "negative dimension";
    }
    ($ix >= 0 and $ix < $dim) or croak "axis index out of range";
    my $self = [(0) x $dim];
    $self->[$ix] = 1;
    bless $self, $class
}

sub _caller_op {
    my $level = (shift||1) + 1;
    my $sub = (caller $level)[3];
    $sub =~ s/.*:://;
    my $op = $op{$sub};
    (defined $op ? $op : $sub);
}

sub _check_dim {
    local ($@, $SIG{__DIE__});
    eval { @{$_[0]} == @{$_[1]} } and return;
    my $op = _caller_op(1);
    my $loc = ($_[2] ? 'first' : 'second');
    UNIVERSAL::isa($_[1], 'ARRAY') or croak "$loc argument to vector operator '$op' is not a vector";
    croak "vector dimensions do not match";
}

sub clone { bless [@{$_[0]}] }

sub set {
    &_check_dim;
    my ($v0, $v1) = @_;
    $v0->[$_] = $v1->[$_] for 0..$#$v1;
}

sub add {
    &_check_dim;
    my ($v0, $v1) = @_;
    bless [map $v0->[$_] + $v1->[$_], 0..$#$v0]
}

sub add_me {
    &_check_dim;
    my ($v0, $v1) = @_;
    $v0->[$_] += $v1->[$_] for 0..$#$v0;
    $v0;
}

sub neg { bless [map -$_, @{$_[0]}] }

sub sub {
    &_check_dim;
    my ($v0, $v1) = ($_[2] ? @_[1, 0] : @_);
    bless [map $v0->[$_] - $v1->[$_], 0..$#$v0]
}

sub sub_me {
    &_check_dim;
    my ($v0, $v1) = @_;
    $v0->[$_] -= $v1->[$_] for 0..$#$v0;
    $v0;
}

sub mul {
    if (ref $_[1]) {
	&_check_dim;
	my ($v0, $v1) = @_;
	my $acu = 0;
	$acu += $v0->[$_] * $v1->[$_] for 0..$#$v0;
	$acu;
    }
    else {
	my ($v, $s) = @_;
	bless [map $s * $_, @$v];
    }
}

sub mul_me {
    ref $_[1] and croak "can not multiply by a vector in place as the result is not a vector";
    my ($v, $s) = @_;
    $_ *= $s for @$v;
    $v
}

sub div {
    croak "can't use vector as dividend"
	if ($_[2] or ref $_[1]);
    my ($v, $div) = @_;
    $div == 0 and croak "illegal division by zero";
    my $i = 1 / $div;
    bless [map $i * $_, @$v]
}

sub div_me {
    croak "can't use vector as dividend" if ref $_[1];
    my $v = shift;
    my $i = 1.0 / shift;
    $_ *= $i for @$v;
    $v;
}

sub equal {
    &_check_dim;
    my ($v0, $v1) = @_;
    $v0->[$_] == $v1->[$_] || return 0 for 0..$#$v0;
    1;
}

sub nequal {
    &_check_dim;
    my ($v0, $v1) = @_;
    $v0->[$_] == $v1->[$_] || return 1 for 0..$#$v0;
    0;
}

sub cross {
    &_check_dim;
    my ($v0, $v1) = ($_[2] ? @_[1, 0] : @_);
    my $dim = @$v0;
    if ($dim == 3) {
	return bless [$v0->[1] * $v1->[2] - $v0->[2] * $v1->[1],
		      $v0->[2] * $v1->[0] - $v0->[0] * $v1->[2],
		      $v0->[0] * $v1->[1] - $v0->[1] * $v1->[0]]
    }
    if ($dim == 7) {
	croak "cross product for dimension 7 not implemented yet, patches welcome!";
    }
    else {
	croak "cross product not defined for dimension $dim"
    }
}

sub as_string { "{" . join(", ", @{$_[0]}). "}" }

sub abs {
    my $acu = 0;
    $acu += $_ * $_ for @{$_[0]};
    sqrt $acu;
}

sub abs2 {
    my $acu = 0;
    $acu += $_ * $_ for @{$_[0]};
    $acu;
}

sub dist {
    &_check_dim;
    my ($v0, $v1) = @_;
    my $d2 = 0;
    for (0..$#$v0) {
	my $d = $v0->[$_] - $v1->[$_];
	$d2 += $d * $d;
    }
    sqrt($d2);
}

sub dist2 {
    &_check_dim;
    my ($v0, $v1) = @_;
    my $d2 = 0;
    for (0..$#$v0) {
	my $d = $v0->[$_] - $v1->[$_];
	$d2 += $d * $d;
    }
    $d2;
}

sub max_component {
    my $max = 0;
    for (@{shift()}) {
	my $abs = CORE::abs($_);
	$abs > $max and $max = $abs;
    }
    $max
}

sub min_component {
    my $self = shift;
    my $min = CORE::abs($self->[0]);
    for (@$self) {
	my $abs = CORE::abs($_);
	$abs < $min and $min = $abs;
    }
    $min
}

sub manhattan_norm {
    my $n = 0;
    $n += CORE::abs($_) for @{$_[0]};
    return $n;
}

sub manhattan_dist {
    &_check_dim;
    my ($v0, $v1) = @_;
    my $d = 0;
    $d += CORE::abs($v0->[$_] - $v1->[$_]) for 0..$#$v0;
    return $d;
}

sub chebyshev_dist {
    &_check_dim;
    my ($v0, $v1) = @_;
    my $max = 0;
    for (0..$#$v0) {
        my $d = CORE::abs($v0->[$_] - $v1->[$_]);
        $max = $d if $d > $max;
    }
    $max;
}

sub _upgrade {
    my $dim;
    map {
	my $d = eval { @{$_} };
	defined $d or croak "argument is not a vector or array";
	if (defined $dim) {
	    $d == $dim or croak "dimensions do not match";
	}
	else {
	    $dim = $d;
	}
	UNIVERSAL::isa($_, __PACKAGE__) ? $_ : clone($_);
    } @_;
}

sub atan2 {
    my ($v0, $v1) = @_;
    if (@$v0 == 2) {
        my $dot = $v0->[0] * $v1->[0] + $v0->[1] * $v1->[1];
        my $cross = $v0->[0] * $v1->[1] - $v0->[1] * $v1->[0];
        return CORE::atan2($cross, $dot);
    }
    else {
        my $a0 = &abs($v0);
        return 0 unless $a0;
        my $u0 = $v0 / $a0;
        my $p = $v1 * $u0;
        CORE::atan2(&abs($v1 - $p * $u0), $p);
    }
}

sub versor {
    my $self = shift;
    my $f = 0;
    $f += $_ * $_ for @$self;
    $f == 0 and croak "Illegal division by zero";
    $f = 1/sqrt $f;
    bless [map $f * $_, @$self]
}

sub wrap {
    my ($self, $v) = @_;
    &_check_dim;

    bless [map  { my $s = $self->[$_];
		  my $c = $v->[$_];
		  $c - $s * POSIX::floor($c/$s) } (0..$#$self)];
}

sub first_orthant_reflection {
    my $self = shift;
    bless [map CORE::abs, @$self];
}

sub sum {
    ref $_[0] or shift; # works both as a class and as an instance method
    my $sum;
    if (@_) {
        $sum = V(@{shift()});
        $sum += $_ for @_;
    }
    return $sum;
}

sub box {
    shift;
    return unless @_;
    my $min = clone(shift);
    my $max = clone($min);
    my $dim = $#$min;
    for (@_) {
        for my $ix (0..$dim) {
            my $c = $_->[$ix];
            if ($max->[$ix] < $c) {
                $max->[$ix] = $c;
            }
            elsif ($min->[$ix] > $c) {
                $min->[$ix] = $c
            }
        }
    }
    wantarray ? ($min, $max) : $max - $min;
}

sub nearest_in_box {
    my $p = shift->clone;
    my ($min, $max) = Math::Vector::Real->box(@_);
    for (0..$#$p) {
        if ($p->[$_] < $min->[$_]) {
            $p->[$_] = $min->[$_];
        }
        elsif ($p->[$_] > $max->[$_]) {
            $p->[$_] = $max->[$_];
        }
    }
    $p
}

sub dist2_to_box {
    @_ > 1 or croak 'Usage: $v->dist2_to_box($w0, ...)';
    my $p = shift;
    my $d2 = 0;
    my ($min, $max) = Math::Vector::Real->box(@_);
    for (0..$#$p) {
        if ($p->[$_] < $min->[$_]) {
            my $d = $p->[$_] - $min->[$_];
            $d2 += $d * $d;
        }
        elsif ($p->[$_] > $max->[$_]) {
            my $d = $p->[$_] - $max->[$_];
            $d2 += $d * $d;
        }
    }
    $d2;
}

sub chebyshev_dist_to_box {
    @_ > 1 or croak 'Usage $v->chebyshev_dist_to_box($w0, ...)';
    my $p = shift;
    my $d = 0;
    my ($min, $max) = Math::Vector::Real->box(@_);
    for (0..$#$p) {
        if ($p->[$_] < $min->[$_]) {
            my $delta = CORE::abs($p->[$_] - $min->[$_]);
            $d = $delta if $delta > $d;
        }
        elsif ($p->[$_] > $max->[$_]) {
            my $delta = CORE::abs($p->[$_] - $min->[$_]);
            $d = $delta if $delta > $d;
        }
    }
    $d;
}

sub chebyshev_cut_box {
    @_ > 2 or croak 'Usage $v->chebyshev_cut_box($cd, $w0, ...)';
    my $p = shift;
    my $cd = shift;
    my ($min, $max) = Math::Vector::Real->box(@_);
    for (0..$#$p) {
        my $a = $p->[$_];
        my $a_min = $a - $cd;
        my $a_max = $a + $cd;
        my $b_min = $min->[$_];
        my $b_max = $max->[$_];
        return if $b_min > $a_max or $b_max < $a_min;
        $min->[$_] = $a_min if $b_min < $a_min;
        $max->[$_] = $a_max if $b_min > $a_max;
    }
    ($min, $max);
}

sub nearest_in_box_border {
    # TODO: this method can be optimized
    my $p = shift->clone;
    my ($b0, $b1) = Math::Vector::Real->box(@_);
    my $in = 0;
    for (0..$#$p) {
        if ($p->[$_] < $b0->[$_]) {
            $p->[$_] = $b0->[$_];
        }
        elsif ($p->[$_] > $b1->[$_]) {
            $p->[$_] = $b1->[$_];
        }
        else {
            $in++;
        }
    }
    if ($in == @$p) {
        # vector was inside the box
        my $min_d = 'inf';
        my ($comp, $comp_ix);
        for my $q ($b0, $b1) {
            for (0..$#$p) {
                my $d = CORE::abs($p->[$_] - $q->[$_]);
                if ($min_d > $d) {
                    $min_d = $d;
                    $comp = $q->[$_];
                    $comp_ix = $_;
                }
            }
        }
        $p->[$comp_ix] = $comp;
    }
    $p;
}

sub max_dist2_to_box {
    @_ > 1 or croak 'Usage: $v->max_dist2_to_box($w0, ...)';
    my $p = shift;
    my ($c0, $c1) = Math::Vector::Real->box(@_);
    my $d2 = 0;
    for (0..$#$p) {
        my $d0 = CORE::abs($c0->[$_] - $p->[$_]);
        my $d1 = CORE::abs($c1->[$_] - $p->[$_]);
        $d2 += ($d0 >= $d1 ? $d0 * $d0 : $d1 * $d1);
    }
    return $d2;
}

sub dist2_between_boxes {
    my ($class, $a0, $a1, $b0, $b1) = @_;
    my ($c0, $c1) = $class->box($a0, $a1);
    my ($d0, $d1) = $class->box($b0, $b1);
    my $d2 = 0;
    for (0..$#$c0) {
        my $e0 = $d0->[$_] - $c1->[$_];
        if ($e0 >= 0) {
            $d2 += $e0 * $e0;
        }
        else {
            my $e1 = $c0->[$_] - $d1->[$_];
            if ($e1 > 0) {
                $d2 += $e1 * $e1;
            }
        }
    }
    $d2;
}

*min_dist2_between_boxes = \&dist2_between_boxes;

sub max_dist2_between_boxes {
    my ($class, $a0, $a1, $b0, $b1) = @_;
    my ($c0, $c1) = $class->box($a0, $a1);
    my ($d0, $d1) = $class->box($b0, $b1);
    my $d2 = 0;
    for (0..$#$c0) {
        my $e0 = $d1->[$_] - $c0->[$_];
        my $e1 = $d0->[$_] - $c1->[$_];
        $e0 *= $e0;
        $e1 *= $e1;
        $d2 += ($e0 > $e1 ? $e0 : $e1);
    }
    $d2;
}

sub max_component_index {
    my $self = shift;
    return unless @$self;
    my $max = 0;
    my $max_ix = 0;
    for my $ix (0..$#$self) {
        my $c = CORE::abs($self->[$ix]);
        if ($c > $max) {
            $max = $c;
            $max_ix = $ix;
        }
    }
    $max_ix;
}

sub min_component_index {
    my $self = shift;
    return unless @$self;
    my $min = CORE::abs($self->[0]);
    my $min_ix = 0;
    for my $ix (1..$#$self) {
        my $c = CORE::abs($self->[$ix]);
        if ($c < $min) {
            $min = $c;
            $min_ix = $ix
        }
    }
    $min_ix;
}

sub decompose {
    my ($u, $v) = @_;
    my $p = $u * ($u * $v)/abs2($u);
    my $n = $v - $p;
    wantarray ? ($p, $n) : $n;
}

sub canonical_base {
    my ($class, $dim) = @_;
    my @base = map { bless [(0) x $dim], $class } 1..$dim;
    $base[$_][$_] = 1 for 0..$#base;
    return @base;
}

sub rotation_base_3d {
    my $v = shift;
    @$v == 3 or croak "rotation_base_3d requires a vector with three dimensions";
    $v = $v->versor;
    my $n = [0, 0, 0];
    for (0..2) {
        if (CORE::abs($v->[$_]) > 0.57) {
            $n->[($_ + 1) % 3] = 1;
            $n = $v->decompose($n)->versor;
            return ($v, $n, $v x $n);
        }
    }
    die "internal error, all the components where smaller than 0.57!";
}

sub rotate_3d {
    my $v = shift;
    my $angle = shift;
    my $c = cos($angle); my $s = sin($angle);
    my ($i, $j, $k) = $v->rotation_base_3d;
    my $rj = $c * $j + $s * $k;
    my $rk = $c * $k - $s * $j;
    if (wantarray) {
        return map { ($_ * $i) * $i + ($_ * $j) * $rj + ($_ * $k) * $rk } @_;
    }
    else {
        my $a = shift;
        return (($a * $i) * $i + ($a * $j) * $rj + ($a * $k) * $rk);
    }
}

sub normal_base { __PACKAGE__->complementary_base(@_) }

sub complementary_base {
    shift;
    @_ or croak "complementaty_base requires at least one argument in order to determine the dimension";
    my $dim = @{$_[0]};
    if ($dim == 2 and @_ == 1) {
        my $u = versor($_[0]);
        @$u = ($u->[1], -$u->[0]);
        return $u;
    }

    my @v = map clone($_), @_;
    my @base = Math::Vector::Real->canonical_base($dim);
    for my $i (0..$#v) {
        my $u = versor($v[$i]);
        $_ = decompose($u, $_) for @v[$i+1 .. $#v];
        $_ = decompose($u, $_) for @base;
    }

    my $last = $#base - @v;
    return if $last < 0;
    for my $i (0 .. $last) {
        my $max = abs2($base[$i]);
        if ($max < 0.3) {
            for my $j ($i+1 .. $#base) {
                my $d2 = abs2($base[$j]);
                if ($d2 > $max) {
                    @base[$i, $j] = @base[$j, $i];
                    last unless $d2 < 0.3;
                    $max = $d2;
                }
            }
        }
        my $versor = $base[$i] = versor($base[$i]);
        $_ = decompose($versor, $_) for @base[$i+1..$#base];
    }
    wantarray ? @base[0..$last] : $base[0];
}

sub select_in_ball {
    my $v = shift;
    my $r = shift;
    my $r2 = $r * $r;
    grep $v->dist2($_) <= $r2, @_;
}

sub select_in_ball_ref2bitmap {
    my $v = shift;
    my $r = shift;
    my $p = shift;
    my $r2 = $r * $r;
    my $bm = "\0" x int((@$p + 7) / 8);
    for my $ix (0..$#$p) {
        vec($bm, $ix, 1) = 1 if $v->dist2($p->[$ix]) <= $r2;
    }
    return $bm;
}

sub dist2_to_segment {
    my ($p, $a, $b) = @_;
    my $ab = $a - $b;
    my $ap = $a - $p;
    my $ap_ab = $ap * $ab;
    return norm2($ap) if $ap_ab <= 0;
    my $x = $ap * $ab / ($ab * $ab);
    return dist2($ap, $ab) if $x >= 1;
    return dist2($ap, $x * $ab);
}

sub dist_to_segment { sqrt(&dist_to_segment) }

sub dist2_between_segments {
    my ($class, $a, $b, $c, $d) = @_;

    my $ab = $a - $b;
    my $cd = $c - $d;
    my $bd = $b - $d;

    if (@$a > 2) {
        my $ab_ab = $ab * $ab;
        my $ab_cd = $ab * $cd;
        my $cd_cd = $cd * $cd;

        if (CORE::abs(1.0 - ($ab_cd * $ab_cd) / ($ab_ab * $cd_cd)) > 1e-10) {
            # This method works for non-parallel segments
            my $ab_bd = $ab * $bd;
            my $bd_cd = $bd * $cd;

            my $D01 = $ab_cd * $ab_cd - $ab_ab * $cd_cd;
            my $D21 = $cd_cd * $ab_bd - $bd_cd * $ab_cd;
            my $x = $D21 / $D01;
            return dist2_to_segment($b, $c, $d) if $x < 0;
            return dist2_to_segment($a, $c, $d) if $x > 1;

            my $D02 = $ab_cd * $ab_bd - $bd_cd * $ab_ab;
            my $y = $D02 / $D01;
            return dist2_to_segment($d, $a, $b) if $y < 0;
            return dist2_to_segment($c, $a, $b) if $y > 1;

            my $p = $b + $ab * $x;
            my $q = $d + $cd * $y;

            return $p->dist2($q);
        }
    }

    # We are in 2D or lines are parallel, we consider the distance
    # between one segment to the vertices of the other one and
    # viceverse and return the minimum.
    my $min_d2 = dist2_to_segment($a, $c, $d);
    my $d2 = dist2_to_segment($b, $c, $d);
    $d2 = dist2_to_segment($c, $a, $b);
    $min_d2 = $d2 if $d2 < $min_d2;
    $d2 = dist2_to_segment($d, $a, $b);
    $min_d2 = $d2 if $d2 < $min_d2;
    return $min_d2;
}

sub dist_between_segments { sqrt(&dist2_between_segments) }

# This is run *after* Math::Vector::Real::XS is loaded!
*norm = \&abs;
*norm2 = \&abs2;
*max = \&max_component;
*min = \&min_component;
*chebyshev_norm = \&max_component;

1;
__END__