File Coverage

blib/lib/Math/BSpline/Basis.pm

Criterion	Covered	Total	%
statement	106	106	100.0
branch	18	18	100.0
condition	9	9	100.0
subroutine	9	9	100.0
pod	3	3	100.0
total	145	145	100.0

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

# ABSTRACT: B-spline basis functions

use Moo 2.002005;
use List::Util 1.26 ('min');
use Ref::Util 0.010 (
    'is_ref',
    'is_plain_hashref',
    'is_blessed_hashref',
    'is_plain_arrayref',
);


around BUILDARGS => sub {
    my ($orig, $class, @args) = @_;
    my $munged_args;

    if (@args == 1) {
        if (!is_ref($args[0])) {
            # We do not understand this and dispatch to Moo (if this
            # is what $orig does, the docu is very sparse).
            return $class->$orig(@args);
        }
        elsif (
            is_plain_hashref($args[0])
            or
            is_blessed_hashref($args[0])
        ) {
            # I am trying to stay as close to Moo's default behavior
            # as I can, this is the only reason why I am supporing
            # hashrefs at all. And since Moo apparently accepts
            # blessed references, I do the same. However, I make a
            # copy, blessed or not.
            #
            # The ugly test is due to an announced change in the
            # behavior of Ref::Util. is_hashref is going to behave
            # like is_plain_hashref does now.  However, the planned
            # replacement called is_any_hashref is not there. So the
            # only future-safe implementation seems to be to use
            # both explicit functions.
            $munged_args = {%{$args[0]}};
        }
        else {
            # We do not understand this and dispatch to Moo (if this
            # is what $orig does, the docu is very sparse).
            return $class->$orig(@args);
        }
    }
    elsif (@args % 2 == 1) {
        # We do not understand this and dispatch to Moo (if this
        # is what $orig does, the docu is very sparse).
        return $class->$orig(@args);
    }
    else {
        $munged_args = {@args};
    }

    if (exists($munged_args->{knot_vector})) {
        # degree is mandatory, so we only deal with the case when it
        # is there. Otherwise we just let Moo do its job.
        if (exists($munged_args->{degree})) {
            # We do not perform any type validation etc, if the
            # attributes are there, we use them assuming that they
            # are valid.
            my $p           = $munged_args->{degree};
            my $U           = $munged_args->{knot_vector};
            my $is_modified = 0;

            # deal with empty array
            if (!defined($U) or !is_plain_arrayref($U) or @$U == 0) {
                $U = [
                    (map { 0 } (0..$p)),
                    (map { 1 } (0..$p)),
                ];
                $is_modified = 1;
            }

            # deal with unsorted
            for (my $i=1;$i<@$U;$i++) {
                if ($U->[$i] < $U->[$i-1]) {
                    $U           = [sort { $a <=> $b } @$U];
                    $is_modified = 1;
                    last;
                }
            }

            # deal with first breakpoint
            for (my $i=1;$i<=$p;$i++) {
                if ($i == @$U or $U->[$i] != $U->[$i-1]) {
                    $U = [@$U] if (!$is_modified);
                    unshift(@$U, $U->[0]);
                    $is_modified = 1;
                }
            }

            # deal with last breakpoint
            if ($U->[-1] == $U->[0]) {
                $U = [@$U] if (!$is_modified);
                push(@$U, $U->[0] + 1);
            }
            for (my $i=-2;$i>=-1-$p;$i--) {
                if ($U->[$i] != $U->[$i+1]) {
                    $U = [@$U] if (!$is_modified);
                    push(@$U, $U->[-1]);
                    $is_modified = 1;
                }
            }

            # deal with excess multiplicity
            for (my $i=$p+1;$i<@$U-1;$i++) {
                while ($i<@$U-1 and $U->[$i] == $U->[$i-$p]) {
                    $U = [@$U] if (!$is_modified);
                    splice(@$U, $i, 1);
                    $is_modified = 1;
                }
            }

            $munged_args->{knot_vector} = $U if ($is_modified);
        }
    }

    return $class->$orig($munged_args);
};



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



has 'knot_vector' => (
    is      => 'lazy',
    builder => sub {
        my ($self) = @_;
        my $p      = $self->degree;

        return [
            (map { 0 } (0..$p)),
            (map { 1 } (0..$p)),
        ]
    }
);



# I use the same variable names as in the NURBS book, although some
# of them are very generic. The use of $p, $U, $P, and $n is
# consistent throughout the relevant chapters of the book.
sub find_knot_span {
    my ($self, $u) = @_;
    my $p          = $self->degree;
    my $U          = $self->knot_vector;
    my $n          = (@$U - 1) - $p - 1;

    # We expect $u in [$U->[$p], $U->[$n+1]]. We only support
    # values outside this range for rounding errors, do not assume
    # that the result makes sense otherwise.
    return $n if ($u >= $U->[$n+1]);
    return $p if ($u <= $U->[$p]);

    # binary search
    my $low  = $p;
    my $high = $n + 1;
    my $mid  = int(($low + $high) / 2);
    while ($u < $U->[$mid] or $u >= $U->[$mid+1]) {
        if ($u < $U->[$mid]) { $high = $mid }
        else                 { $low  = $mid }
        $mid = int(($low + $high) / 2);
    }

    return $mid;
}



# The variable names are inspired by the theory as laid out in the
# NURBS book. We want to calculate N_{i,p}, that inspires $N and
# $p. U is the knot vector, left and right are inspired by the
# terms in the formulas used in the theoretical derivation.
sub evaluate_basis_functions {
    my ($self, $i, $u) = @_;
    my $p              = $self->degree;
    my $U              = $self->knot_vector;
    my $n              = (@$U - 1) - $p - 1;

    if ($u < $U->[$p] or $u > $U->[$n+1]) {
        return [map { 0 } (0..$p)];
    }

    my $N     = [1];
    my $left  = [];
    my $right = [];
    for (my $j=1;$j<=$p;$j++) {
        $left->[$j]  = $u - $U->[$i+1-$j];
        $right->[$j] = $U->[$i+$j] - $u;
        my $saved = 0;
        for (my $r=0;$r<$j;$r++) {
            my $temp = $N->[$r] / ($right->[$r+1] + $left->[$j-$r]);
            $N->[$r] = $saved + $right->[$r+1] * $temp;
            $saved   = $left->[$j-$r] * $temp;
        }
        $N->[$j] = $saved;
    }

    return $N;
}



sub evaluate_basis_derivatives {
    my ($self, $i, $u, $d) = @_;
    my $p                  = $self->degree;
    my $U                  = $self->knot_vector;
    my $n                  = (@$U - 1) - $p - 1;
    my $result             = [];

    $d = min($d, $p);

    if ($u < $U->[$p] or $u > $U->[$n+1]) {
        for (my $k=0;$k<=$d;$k++) {
            push(@$result, [map { 0 } (0..$p)]);
        }
        return $result;
    }

    my $ndu   = [[1]];
    my $left  = [];
    my $right = [];
    for (my $j=1;$j<=$p;$j++) {
        $left->[$j]  = $u - $U->[$i+1-$j];
        $right->[$j] = $U->[$i+$j] - $u;
        my $saved = 0;
        for (my $r=0;$r<$j;$r++) {
            $ndu->[$j]->[$r] = $right->[$r+1] + $left->[$j-$r];
            my $temp = $ndu->[$r]->[$j-1] / $ndu->[$j]->[$r];
            $ndu->[$r]->[$j] = $saved + $right->[$r+1] * $temp;
            $saved           = $left->[$j-$r] * $temp;
        }
        $ndu->[$j]->[$j] = $saved;
    }

    # $result->[0] holds the function values (0th derivatives)
    for (my $j=0;$j<=$p;$j++) {
        $result->[0]->[$j] = $ndu->[$j]->[$p];
    }

    for (my $r=0;$r<=$p;$r++) {
        my $a         = [[1]];
        my ($l1, $l2) = (0, 1);  # alternating indices to address $a

        # compute $result->[$k] (kth derivative)
        for (my $k=1;$k<=$d;$k++) {
            my $sum = 0;
            my $rk  = $r - $k;
            my $pk  = $p - $k;
            if ($rk >= 0) {
                $a->[$l2]->[0] = $a->[$l1]->[0] / $ndu->[$pk+1]->[$rk];
                $sum = $a->[$l2]->[0] * $ndu->[$rk]->[$pk];
            }

            my $j_min = $rk >= -1      ? 1      : -$rk;
            my $j_max = $r  <= $pk + 1 ? $k - 1 : $p - $r;
            for (my $j=$j_min;$j<=$j_max;$j++) {
                $a->[$l2]->[$j] = ($a->[$l1]->[$j] - $a->[$l1]->[$j-1])
                    / $ndu->[$pk+1]->[$rk+$j];
                $sum += $a->[$l2]->[$j] * $ndu->[$rk+$j]->[$pk];
            }

            if ($r <= $pk) {
                $a->[$l2]->[$k] = -$a->[$l1]->[$k-1]
                    / $ndu->[$pk+1]->[$r];
                $sum += $a->[$l2]->[$k] * $ndu->[$r]->[$pk];
            }

            $result->[$k]->[$r] = $sum;
            ($l1, $l2) = ($l2, $l1);
        }
    }

    my $multiplicity = $p;
    for (my $k=1;$k<=$d;$k++) {
        for (my $j=0;$j<=$p;$j++) {
            $result->[$k]->[$j] *= $multiplicity;
        }
        $multiplicity *= ($p - $k);
    }

    return $result;
}


1;

__END__