|
6958
|
1 /* Copyright (C) 2009 Carlo de Falco |
|
|
2 |
|
|
3 This program is free software: you can redistribute it and/or modify |
|
|
4 it under the terms of the GNU General Public License as published by |
|
|
5 the Free Software Foundation, either version 2 of the License, or |
|
|
6 (at your option) any later version. |
|
6052
|
7 |
|
6958
|
8 This program is distributed in the hope that it will be useful, |
|
|
9 but WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
|
11 GNU General Public License for more details. |
|
6052
|
12 |
|
6958
|
13 You should have received a copy of the GNU General Public License |
|
|
14 along with this program. If not, see <http://www.gnu.org/licenses/>. |
|
6052
|
15 */ |
|
|
16 |
|
|
17 #include <octave/oct.h> |
|
|
18 #include <iostream> |
|
|
19 |
|
|
20 double onebasisfun__ (double u, octave_idx_type p, RowVector U) |
|
|
21 { |
|
|
22 |
|
|
23 //std::cout << "u=" << u << " " << "p=" << p << " \n" << "U=" << U; |
|
|
24 |
|
|
25 double N = 0.0; |
|
|
26 if ((u < U.min ()) || ( u > U.max ())) |
|
|
27 return (N); |
|
|
28 else if (p == 0) |
|
|
29 return (1.0); |
|
|
30 |
|
|
31 double ln = u - U(0); |
|
|
32 double ld = U(U.length () - 2) - U(0); |
|
|
33 if (ld != 0) |
|
|
34 N += ln * onebasisfun__ (u, p-1, U.extract (0, U.length () - 2))/ ld; |
|
|
35 |
|
|
36 double dn = U(U.length () - 1) - u; |
|
|
37 double dd = U(U.length () - 1) - U(1); |
|
|
38 if (dd != 0) |
|
|
39 N += dn * onebasisfun__ (u, p-1, U.extract (1, U.length () - 1))/ dd; |
|
|
40 |
|
|
41 return (N); |
|
|
42 } |
|
|
43 |
|
|
44 |
|
|
45 DEFUN_DLD(tbasisfun, args, nargout,"\ |
|
6464
|
46 TBASISFUN: Compute a B- or T-Spline basis function from its local knot vector.\n\ |
|
|
47 \n\ |
|
6052
|
48 usage:\n\ |
|
|
49 \n\ |
|
|
50 N = tbasisfun (u, p, U)\n\ |
|
|
51 N = tbasisfun ([u; v], [p q], {U, V})\n\ |
|
|
52 \n\ |
|
|
53 INPUT:\n\ |
|
|
54 u or [u; v] : points in parameter space where the basis function is to be\n\ |
|
|
55 evaluated \n\ |
|
|
56 \n\ |
|
|
57 U or {U, V} : local knot vector\n\ |
|
|
58 \n\ |
|
6958
|
59 p or [p q] : polynomial order of the basis function\n\ |
|
6052
|
60 \n\ |
|
|
61 OUTPUT:\n\ |
|
|
62 N : basis function evaluated at the given parametric points\n") |
|
|
63 |
|
|
64 { |
|
|
65 |
|
|
66 Matrix u = args(0).matrix_value (); |
|
|
67 |
|
|
68 if (! args(2).is_cell ()) |
|
|
69 { |
|
|
70 |
|
|
71 double p = args(1).idx_type_value (); |
|
|
72 RowVector U = args(2).row_vector_value (true, true); |
|
|
73 assert (U.numel () == p+2); |
|
|
74 |
|
|
75 RowVector N(u.cols ()); |
|
|
76 for (octave_idx_type ii=0; ii<u.numel (); ii++) |
|
|
77 N(ii) = onebasisfun__ (u(ii), p, U); |
|
|
78 |
|
|
79 } else { |
|
|
80 |
|
|
81 RowVector p = args(1).row_vector_value (); |
|
|
82 Cell C = args(2).cell_value (); |
|
|
83 RowVector U = C(0).row_vector_value (true, true); |
|
|
84 RowVector V = C(1).row_vector_value (true, true); |
|
|
85 |
|
|
86 |
|
|
87 RowVector N(u.cols ()); |
|
|
88 for (octave_idx_type ii=0; ii<u.cols (); ii++) |
|
|
89 { |
|
|
90 N(ii) = onebasisfun__ (u(0, ii), octave_idx_type(p(0)), U) * |
|
|
91 onebasisfun__ (u(1, ii), octave_idx_type(p(1)), V); |
|
|
92 //std::cout << "N=" << N(ii) << "\n\n\n"; |
|
|
93 } |
|
|
94 return octave_value (N); |
|
|
95 } |
|
|
96 } |
|
|
97 |
|
|
98 |
|
|
99 /* |
|
|
100 %!demo |
|
|
101 %! U = {[0 0 1/2 1 1], [0 0 0 1 1]}; |
|
|
102 %! p = [3, 3]; |
|
|
103 %! [X, Y] = meshgrid (linspace(0, 1, 30)); |
|
|
104 %! u = [X(:), Y(:)]'; |
|
|
105 %! N = tbasisfun (u, p, U); |
|
|
106 %! surf (X, Y, reshape (N, size(X))) |
|
|
107 */ |
