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


line stmt bran cond sub pod time code
1             package Math::BSpline::Basis;
2             $Math::BSpline::Basis::VERSION = '0.001';
3 5     5   4019 use 5.014;
  5         18  
4 5     5   30 use warnings;
  5         9  
  5         189  
5              
6             # ABSTRACT: B-spline basis functions
7              
8 5     5   2945 use Moo 2.002005;
  5         58851  
  5         33  
9 5     5   7578 use List::Util 1.26 ('min');
  5         95  
  5         601  
10             use Ref::Util 0.010 (
11 5         8368 'is_ref',
12             'is_plain_hashref',
13             'is_blessed_hashref',
14             'is_plain_arrayref',
15 5     5   2619 );
  5         8524  
16              
17              
18             around BUILDARGS => sub {
19             my ($orig, $class, @args) = @_;
20             my $munged_args;
21              
22             if (@args == 1) {
23             if (!is_ref($args[0])) {
24             # We do not understand this and dispatch to Moo (if this
25             # is what $orig does, the docu is very sparse).
26             return $class->$orig(@args);
27             }
28             elsif (
29             is_plain_hashref($args[0])
30             or
31             is_blessed_hashref($args[0])
32             ) {
33             # I am trying to stay as close to Moo's default behavior
34             # as I can, this is the only reason why I am supporing
35             # hashrefs at all. And since Moo apparently accepts
36             # blessed references, I do the same. However, I make a
37             # copy, blessed or not.
38             #
39             # The ugly test is due to an announced change in the
40             # behavior of Ref::Util. is_hashref is going to behave
41             # like is_plain_hashref does now. However, the planned
42             # replacement called is_any_hashref is not there. So the
43             # only future-safe implementation seems to be to use
44             # both explicit functions.
45             $munged_args = {%{$args[0]}};
46             }
47             else {
48             # We do not understand this and dispatch to Moo (if this
49             # is what $orig does, the docu is very sparse).
50             return $class->$orig(@args);
51             }
52             }
53             elsif (@args % 2 == 1) {
54             # We do not understand this and dispatch to Moo (if this
55             # is what $orig does, the docu is very sparse).
56             return $class->$orig(@args);
57             }
58             else {
59             $munged_args = {@args};
60             }
61              
62             if (exists($munged_args->{knot_vector})) {
63             # degree is mandatory, so we only deal with the case when it
64             # is there. Otherwise we just let Moo do its job.
65             if (exists($munged_args->{degree})) {
66             # We do not perform any type validation etc, if the
67             # attributes are there, we use them assuming that they
68             # are valid.
69             my $p = $munged_args->{degree};
70             my $U = $munged_args->{knot_vector};
71             my $is_modified = 0;
72              
73             # deal with empty array
74             if (!defined($U) or !is_plain_arrayref($U) or @$U == 0) {
75             $U = [
76             (map { 0 } (0..$p)),
77             (map { 1 } (0..$p)),
78             ];
79             $is_modified = 1;
80             }
81              
82             # deal with unsorted
83             for (my $i=1;$i<@$U;$i++) {
84             if ($U->[$i] < $U->[$i-1]) {
85             $U = [sort { $a <=> $b } @$U];
86             $is_modified = 1;
87             last;
88             }
89             }
90              
91             # deal with first breakpoint
92             for (my $i=1;$i<=$p;$i++) {
93             if ($i == @$U or $U->[$i] != $U->[$i-1]) {
94             $U = [@$U] if (!$is_modified);
95             unshift(@$U, $U->[0]);
96             $is_modified = 1;
97             }
98             }
99              
100             # deal with last breakpoint
101             if ($U->[-1] == $U->[0]) {
102             $U = [@$U] if (!$is_modified);
103             push(@$U, $U->[0] + 1);
104             }
105             for (my $i=-2;$i>=-1-$p;$i--) {
106             if ($U->[$i] != $U->[$i+1]) {
107             $U = [@$U] if (!$is_modified);
108             push(@$U, $U->[-1]);
109             $is_modified = 1;
110             }
111             }
112              
113             # deal with excess multiplicity
114             for (my $i=$p+1;$i<@$U-1;$i++) {
115             while ($i<@$U-1 and $U->[$i] == $U->[$i-$p]) {
116             $U = [@$U] if (!$is_modified);
117             splice(@$U, $i, 1);
118             $is_modified = 1;
119             }
120             }
121              
122             $munged_args->{knot_vector} = $U if ($is_modified);
123             }
124             }
125              
126             return $class->$orig($munged_args);
127             };
128              
129              
130              
131             has 'degree' => (
132             is => 'ro',
133             required => 1,
134             );
135              
136              
137              
138             has 'knot_vector' => (
139             is => 'lazy',
140             builder => sub {
141 2     2   2666 my ($self) = @_;
142 2         7 my $p = $self->degree;
143              
144             return [
145 6         12 (map { 0 } (0..$p)),
146 2         6 (map { 1 } (0..$p)),
  6         27  
147             ]
148             }
149             );
150              
151              
152              
153             # I use the same variable names as in the NURBS book, although some
154             # of them are very generic. The use of $p, $U, $P, and $n is
155             # consistent throughout the relevant chapters of the book.
156             sub find_knot_span {
157 209     209 1 1201355 my ($self, $u) = @_;
158 209         590 my $p = $self->degree;
159 209         5461 my $U = $self->knot_vector;
160 209         1888 my $n = (@$U - 1) - $p - 1;
161              
162             # We expect $u in [$U->[$p], $U->[$n+1]]. We only support
163             # values outside this range for rounding errors, do not assume
164             # that the result makes sense otherwise.
165 209 100       941 return $n if ($u >= $U->[$n+1]);
166 194 100       600 return $p if ($u <= $U->[$p]);
167              
168             # binary search
169 178         366 my $low = $p;
170 178         269 my $high = $n + 1;
171 178         465 my $mid = int(($low + $high) / 2);
172 178   100     756 while ($u < $U->[$mid] or $u >= $U->[$mid+1]) {
173 188 100       403 if ($u < $U->[$mid]) { $high = $mid }
  114         170  
174 74         116 else { $low = $mid }
175 188         774 $mid = int(($low + $high) / 2);
176             }
177              
178 178         853 return $mid;
179             }
180              
181              
182              
183             # The variable names are inspired by the theory as laid out in the
184             # NURBS book. We want to calculate N_{i,p}, that inspires $N and
185             # $p. U is the knot vector, left and right are inspired by the
186             # terms in the formulas used in the theoretical derivation.
187             sub evaluate_basis_functions {
188 43     43 1 1885 my ($self, $i, $u) = @_;
189 43         90 my $p = $self->degree;
190 43         768 my $U = $self->knot_vector;
191 43         306 my $n = (@$U - 1) - $p - 1;
192              
193 43 100 100     173 if ($u < $U->[$p] or $u > $U->[$n+1]) {
194 2         6 return [map { 0 } (0..$p)];
  6         15  
195             }
196              
197 41         85 my $N = [1];
198 41         66 my $left = [];
199 41         62 my $right = [];
200 41         98 for (my $j=1;$j<=$p;$j++) {
201 122         227 $left->[$j] = $u - $U->[$i+1-$j];
202 122         213 $right->[$j] = $U->[$i+$j] - $u;
203 122         172 my $saved = 0;
204 122         226 for (my $r=0;$r<$j;$r++) {
205 243         417 my $temp = $N->[$r] / ($right->[$r+1] + $left->[$j-$r]);
206 243         408 $N->[$r] = $saved + $right->[$r+1] * $temp;
207 243         456 $saved = $left->[$j-$r] * $temp;
208             }
209 122         289 $N->[$j] = $saved;
210             }
211              
212 41         127 return $N;
213             }
214              
215              
216              
217             sub evaluate_basis_derivatives {
218 98     98 1 1202 my ($self, $i, $u, $d) = @_;
219 98         205 my $p = $self->degree;
220 98         1782 my $U = $self->knot_vector;
221 98         731 my $n = (@$U - 1) - $p - 1;
222 98         162 my $result = [];
223              
224 98         469 $d = min($d, $p);
225              
226 98 100 100     596 if ($u < $U->[$p] or $u > $U->[$n+1]) {
227 2         8 for (my $k=0;$k<=$d;$k++) {
228 8         18 push(@$result, [map { 0 } (0..$p)]);
  32         55  
229             }
230 2         7 return $result;
231             }
232              
233 96         250 my $ndu = [[1]];
234 96         167 my $left = [];
235 96         173 my $right = [];
236 96         282 for (my $j=1;$j<=$p;$j++) {
237 359         712 $left->[$j] = $u - $U->[$i+1-$j];
238 359         735 $right->[$j] = $U->[$i+$j] - $u;
239 359         518 my $saved = 0;
240 359         758 for (my $r=0;$r<$j;$r++) {
241 897         1771 $ndu->[$j]->[$r] = $right->[$r+1] + $left->[$j-$r];
242 897         1585 my $temp = $ndu->[$r]->[$j-1] / $ndu->[$j]->[$r];
243 897         1737 $ndu->[$r]->[$j] = $saved + $right->[$r+1] * $temp;
244 897         1753 $saved = $left->[$j-$r] * $temp;
245             }
246 359         855 $ndu->[$j]->[$j] = $saved;
247             }
248              
249             # $result->[0] holds the function values (0th derivatives)
250 96         236 for (my $j=0;$j<=$p;$j++) {
251 455         910 $result->[0]->[$j] = $ndu->[$j]->[$p];
252             }
253              
254 96         249 for (my $r=0;$r<=$p;$r++) {
255 455         871 my $a = [[1]];
256 455         749 my ($l1, $l2) = (0, 1); # alternating indices to address $a
257              
258             # compute $result->[$k] (kth derivative)
259 455         844 for (my $k=1;$k<=$d;$k++) {
260 1414         1893 my $sum = 0;
261 1414         1951 my $rk = $r - $k;
262 1414         1868 my $pk = $p - $k;
263 1414 100       2378 if ($rk >= 0) {
264 820         1626 $a->[$l2]->[0] = $a->[$l1]->[0] / $ndu->[$pk+1]->[$rk];
265 820         1245 $sum = $a->[$l2]->[0] * $ndu->[$rk]->[$pk];
266             }
267              
268 1414 100       2278 my $j_min = $rk >= -1 ? 1 : -$rk;
269 1414 100       2369 my $j_max = $r <= $pk + 1 ? $k - 1 : $p - $r;
270 1414         2613 for (my $j=$j_min;$j<=$j_max;$j++) {
271 731         1520 $a->[$l2]->[$j] = ($a->[$l1]->[$j] - $a->[$l1]->[$j-1])
272             / $ndu->[$pk+1]->[$rk+$j];
273 731         1544 $sum += $a->[$l2]->[$j] * $ndu->[$rk+$j]->[$pk];
274             }
275              
276 1414 100       2282 if ($r <= $pk) {
277 820         1759 $a->[$l2]->[$k] = -$a->[$l1]->[$k-1]
278             / $ndu->[$pk+1]->[$r];
279 820         1301 $sum += $a->[$l2]->[$k] * $ndu->[$r]->[$pk];
280             }
281              
282 1414         2241 $result->[$k]->[$r] = $sum;
283 1414         3703 ($l1, $l2) = ($l2, $l1);
284             }
285             }
286              
287 96         167 my $multiplicity = $p;
288 96         229 for (my $k=1;$k<=$d;$k++) {
289 282         514 for (my $j=0;$j<=$p;$j++) {
290 1414         2412 $result->[$k]->[$j] *= $multiplicity;
291             }
292 282         521 $multiplicity *= ($p - $k);
293             }
294              
295 96         412 return $result;
296             }
297              
298              
299             1;
300              
301             __END__