|
5683
|
1 /* Copyright (C) 2009 Carlo de Falco |
|
|
2 |
|
6958
|
3 This program is free software: you can redistribute it and/or modify |
|
5683
|
4 it under the terms of the GNU General Public License as published by |
|
11634
|
5 the Free Software Foundation, either version 3 of the License, or |
|
5683
|
6 (at your option) any later version. |
|
6958
|
7 |
|
5683
|
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. |
|
|
12 |
|
|
13 You should have received a copy of the GNU General Public License |
|
6958
|
14 along with this program. If not, see <http://www.gnu.org/licenses/>. |
|
5683
|
15 */ |
|
|
16 |
|
|
17 #include <octave/oct.h> |
|
|
18 #include "low_level_functions.h" |
|
7966
|
19 #include <omp.h> |
|
5683
|
20 |
|
|
21 static bool bspeval_bad_arguments(const octave_value_list& args); |
|
|
22 |
|
|
23 DEFUN_DLD(bspeval, args, nargout,"\ |
|
|
24 BSPEVAL: Evaluate B-Spline at parametric points\n\ |
|
|
25 \n\ |
|
|
26 \n\ |
|
|
27 Calling Sequence:\n\ |
|
|
28 \n\ |
|
|
29 p = bspeval(d,c,k,u)\n\ |
|
|
30 \n\ |
|
|
31 INPUT:\n\ |
|
|
32 \n\ |
|
|
33 d - Degree of the B-Spline.\n\ |
|
|
34 c - Control Points, matrix of size (dim,nc).\n\ |
|
|
35 k - Knot sequence, row vector of size nk.\n\ |
|
|
36 u - Parametric evaluation points, row vector of size nu.\n\ |
|
|
37 \n\ |
|
|
38 OUTPUT:\n\ |
|
|
39 \n\ |
|
|
40 p - Evaluated points, matrix of size (dim,nu)\n\ |
|
|
41 ") |
|
|
42 { |
|
|
43 |
|
|
44 octave_value_list retval; |
|
7203
|
45 if (!bspeval_bad_arguments (args)) |
|
|
46 { |
|
9760
|
47 int d = args(0).int_value(); |
|
|
48 Matrix c = args(1).matrix_value(); |
|
|
49 RowVector k = args(2).row_vector_value(); |
|
|
50 NDArray u = args(3).array_value(); |
|
7203
|
51 |
|
12672
|
52 octave_idx_type nu = u.numel (); |
|
7203
|
53 octave_idx_type mc = c.rows(), |
|
|
54 nc = c.cols(); |
|
|
55 |
|
|
56 Matrix p(mc, nu, 0.0); |
|
|
57 |
|
|
58 if (!error_state) |
|
|
59 { |
|
12672
|
60 if (nc + d == k.numel () - 1) |
|
7203
|
61 { |
|
12672
|
62 //#pragma omp parallel default (none) shared (d, c, k, u, nu, mc, nc, p) |
|
7966
|
63 { |
|
|
64 RowVector N(d+1,0.0); |
|
|
65 int s, tmp1; |
|
|
66 double tmp2; |
|
12672
|
67 //#pragma omp for |
|
7966
|
68 for (octave_idx_type col=0; col<nu; col++) |
|
|
69 { |
|
|
70 //printf ("thread %d, col %d\n", omp_get_thread_num (), col); |
|
|
71 s = findspan (nc-1, d, u(col), k); |
|
|
72 basisfun (s, u(col), d, k, N); |
|
|
73 tmp1 = s - d; |
|
|
74 for (octave_idx_type row(0); row<mc; row++) |
|
|
75 { |
|
|
76 tmp2 = 0.0; |
|
|
77 for ( octave_idx_type i(0); i<=d; i++) |
|
|
78 tmp2 += N(i)*c(row,tmp1+i); |
|
|
79 p(row,col) = tmp2; |
|
|
80 } |
|
|
81 } |
|
|
82 }// end omp |
|
|
83 } |
|
7203
|
84 else |
|
|
85 { |
|
12672
|
86 error("inconsistent bspline data, d + columns(c) != numel(k) - 1."); |
|
7203
|
87 } |
|
|
88 retval(0) = octave_value(p); |
|
|
89 } |
|
|
90 } |
|
5683
|
91 return retval; |
|
|
92 } |
|
|
93 |
|
7203
|
94 static bool bspeval_bad_arguments (const octave_value_list& args) |
|
5683
|
95 { |
|
12672
|
96 if (args.length () != 4) |
|
5683
|
97 { |
|
7203
|
98 error("bspeval: wrong number of input arguments."); |
|
5683
|
99 return true; |
|
|
100 } |
|
|
101 if (!args(0).is_real_scalar()) |
|
|
102 { |
|
7203
|
103 error("bspeval: degree should be a scalar."); |
|
5683
|
104 return true; |
|
|
105 } |
|
|
106 if (!args(1).is_real_matrix()) |
|
|
107 { |
|
7203
|
108 error("bspeval: the control net should be a matrix of doubles."); |
|
5683
|
109 return true; |
|
|
110 } |
|
|
111 if (!args(2).is_real_matrix()) |
|
|
112 { |
|
7203
|
113 error("bspeval: the knot vector should be a real vector."); |
|
5683
|
114 return true; |
|
|
115 } |
|
|
116 if (!args(3).is_real_type()) |
|
|
117 { |
|
7203
|
118 error("bspeval: the set of parametric points should be an array of doubles."); |
|
5683
|
119 return true; |
|
|
120 } |
|
|
121 return false; |
|
|
122 } |
