File Coverage

blib/lib/Math/BSpline/Curve.pm

Criterion	Covered	Total	%
statement	79	91	86.8
branch	9	18	50.0
condition			n/a
subroutine	11	12	91.6
pod	3	3	100.0
total	102	124	82.2

package Math::BSpline::Curve;
$Math::BSpline::Curve::VERSION = '0.001';
use 5.014;
use warnings;

# ABSTRACT: B-spline curves

use Moo 2.002005;
use List::Util 1.26 ('min');
use Ref::Util 0.010 (
    'is_plain_arrayref',
);
use Math::BSpline::Basis 0.001;
use Math::Matrix::Banded 0.004;


has '_degree' => (
    is       => 'ro',
    required => 1,
    init_arg => 'degree',
);



has '_knot_vector' => (
    is        => 'ro',
    init_arg  => 'knot_vector',
    predicate => 1,
);



has 'control_points' => (
    is      => 'lazy',
    builder => sub { return [] },
);



has 'basis' => (
    is      => 'lazy',
    handles => [
        'degree',
        'knot_vector',
    ],
    builder => sub {
        my ($self) = @_;

        return Math::BSpline::Basis->new(
            degree      => $self->_degree,
            (
                $self->_has_knot_vector
                    ? (knot_vector => $self->_knot_vector)
                    : (),
            ),
        )
    }
);



sub evaluate {
    my ($self, $u) = @_;
    my $basis      = $self->basis;

    my $p   = $self->degree;
    my $P   = $self->control_points;
    my $s   = $basis->find_knot_span($u);
    my $Nip = $basis->evaluate_basis_functions($s, $u);

    return undef if (!@$P);
    my $value;
    if (is_plain_arrayref($P->[0])) {
        # The control points are plain arrayrefs, hence we have no
        # overloaded scalar multiplication at our disposal and have
        # to manipulate the components individually.
        my $dim = scalar(@{$P->[0]});
        $value = [map { 0 } (1..$dim)];
        for (my $i=0;$i<=$p;$i++) {
            my $c      = $Nip->[$i];
            my $this_P = $P->[$s-$p+$i];
            for (my $j=0;$j<$dim;$j++) {
                $value->[$j] += $c * $this_P->[$j];
            }
        }
    }
    else {
        # We use the first control point to initialize the value in
        # order to support all objects that overload addition and
        # scalar multiplication.
        $value = 0 * $P->[0];
        for (my $i=0;$i<=$p;$i++) {
            $value += $Nip->[$i] * $P->[$s-$p+$i];
        }
    }

    return $value;
}



sub evaluate_derivatives {
    my ($self, $u, $d) = @_;
    my $basis          = $self->basis;

    my $p = $self->degree;
    my $P = $self->control_points;
    my $s = $basis->find_knot_span($u);
    my $D = $basis->evaluate_basis_derivatives($s, $u, min($d, $p));

    return undef if (!@$P);
    my $value = [];
    if (is_plain_arrayref($P->[0])) {
        # The control points are plain arrayrefs, hence we have no
        # overloaded scalar multiplication at our disposal and have
        # to manipulate the components individually.
        my $dim = scalar(@{$P->[0]});
        for (my $k=0;$k<=$d;$k++) {
            $value->[$k] = [map { 0 } (1..$dim)];

            if ($k <= $p) {
                for (my $i=0;$i<=$p;$i++) {
                    my $c      = $D->[$k]->[$i];
                    my $this_P = $P->[$s-$p+$i];
                    for (my $j=0;$j<$dim;$j++) {
                        $value->[$k]->[$j] += $c * $this_P->[$j];
                    }
                }
            }
        }
    }
    else {
        # We use the first control point to initialize the value in
        # order to support all objects that overload addition and
        # scalar multiplication.
        for (my $k=0;$k<=$d;$k++) {
            $value->[$k] = 0 * $P->[0];

            if ($k <= $p) {
                for (my $i=0;$i<=$p;$i++) {
                    $value->[$k] += $D->[$k]->[$i] * $P->[$s-$p+$i];
                }
            }
        }
    }

    return $value;
}


sub derivative {
    my ($self) = @_;
    my $p      = $self->degree;
    my $P      = $self->control_points;
    my $U      = $self->knot_vector;

    return undef if (!@$P);

    my $q = $p - 1;
    my $V = [@$U[1..($#$U-1)]];
    my $Q = [];
    if (is_plain_arrayref($P->[0])) {
        # The control points are plain arrayrefs, hence we have no
        # overloaded scalar multiplication at our disposal and have
        # to manipulate the components individually.
        my $dim = scalar(@{$P->[0]});
        for (my $i=0;$i<@$P-1;$i++) {
            my $c = $p / ($U->[$i+$p+1] - $U->[$i+1]);
            $Q->[$i] = [];
            for (my $j=0;$j<$dim;$j++) {
                $Q->[$i]->[$j] = $c * ($P->[$i+1]->[$j] - $P->[$i]->[$j]);
            }
        }
    }
    else {
        for (my $i=0;$i<@$P-1;$i++) {
            my $c = $p / ($U->[$i+$p+1] - $U->[$i+1]);
            $Q->[$i] = $c * ($P->[$i+1] - $P->[$i]);
        }
    }

    return Math::BSpline::Curve->new(
        degree         => $q,
        knot_vector    => $V,
        control_points => $Q,
    );
}


1;

__END__