|
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 <octave/oct-map.h> |
|
|
19 #include <octave/parse.h> |
|
|
20 #include "low_level_functions.h" |
|
|
21 |
|
|
22 DEFUN_DLD(nrb_srf_basisfun__, args, nargout,"\ |
|
|
23 NRB_SRF_BASISFUN__: Undocumented private function\ |
|
|
24 ") |
|
|
25 { |
|
|
26 |
|
|
27 octave_value_list retval, newargs; |
|
|
28 |
|
|
29 const NDArray points = args(0).array_value(); |
|
12027
|
30 const octave_scalar_map nrb = args(1).scalar_map_value(); |
|
5683
|
31 |
|
|
32 if (!error_state) |
|
|
33 { |
|
|
34 |
|
12026
|
35 const Cell knots = nrb.contents("knots").cell_value(); |
|
|
36 const NDArray coefs = nrb.contents("coefs").array_value(); |
|
|
37 octave_idx_type m = static_cast<octave_idx_type> ((nrb.contents("number").vector_value())(0)) - 1; // m = size (nrb.coefs, 2) -1; |
|
|
38 octave_idx_type n = static_cast<octave_idx_type> ((nrb.contents("number").vector_value())(1)) - 1; // n = size (nrb.coefs, 3) -1; |
|
|
39 octave_idx_type p = static_cast<octave_idx_type> ((nrb.contents("order").vector_value())(0)) - 1; // p = nrb.order(1) -1; |
|
|
40 octave_idx_type q = static_cast<octave_idx_type> ((nrb.contents("order").vector_value())(1)) - 1; // q = nrb.order(2) -1; |
|
5683
|
41 |
|
8034
|
42 Array<idx_vector> idx(dim_vector (2, 1), idx_vector(':')); |
|
5787
|
43 idx(0) = 0; |
|
5683
|
44 const NDArray u(points.index (idx).squeeze ()); // u = points(1,:); |
|
|
45 |
|
5787
|
46 idx(0) = 1; |
|
5683
|
47 const NDArray v(points.index (idx).squeeze ()); // v = points(2,:); |
|
|
48 |
|
12672
|
49 octave_idx_type npt = u.numel (); // npt = length(u); |
|
5683
|
50 RowVector M(p+1, 0.0), N (q+1, 0.0); |
|
|
51 Matrix RIkJk(npt, (p+1)*(q+1), 0.0); |
|
|
52 Matrix indIkJk(npt, (p+1)*(q+1), 0.0); |
|
|
53 RowVector denom(npt, 0.0); |
|
|
54 |
|
|
55 const RowVector U(knots(0).row_vector_value ()); // U = nrb.knots{1}; |
|
|
56 |
|
|
57 const RowVector V(knots(1).row_vector_value ()); // V = nrb.knots{2}; |
|
|
58 |
|
8034
|
59 Array<idx_vector> idx2(dim_vector (3, 1), idx_vector(':')); idx2(0) = 3; |
|
5683
|
60 NDArray w (coefs.index (idx2).squeeze ()); // w = squeeze(nrb.coefs(4,:,:)); |
|
|
61 |
|
|
62 RowVector spu(u); |
|
|
63 for (octave_idx_type ii(0); ii < npt; ii++) |
|
|
64 { |
|
|
65 spu(ii) = findspan(m, p, u(ii), U); |
|
|
66 } // spu = findspan (m, p, u, U); |
|
|
67 |
|
|
68 newargs(3) = U; newargs(2) = p; newargs(1) = u; newargs(0) = spu; |
|
|
69 Matrix Ik = feval (std::string("numbasisfun"), newargs, 1)(0).matrix_value (); // Ik = numbasisfun (spu, u, p, U); |
|
|
70 |
|
|
71 RowVector spv(v); |
|
12672
|
72 for (octave_idx_type ii(0); ii < v.numel (); ii++) |
|
5683
|
73 { |
|
|
74 spv(ii) = findspan(n, q, v(ii), V); |
|
|
75 } // spv = findspan (n, q, v, V); |
|
|
76 |
|
|
77 newargs(3) = V; newargs(2) = q; newargs(1) = v; newargs(0) = spv; |
|
|
78 Matrix Jk = feval (std::string("numbasisfun"), newargs, 1)(0).matrix_value (); // Jk = numbasisfun (spv, v, q, V); |
|
|
79 |
|
|
80 Matrix NuIkuk(npt, p+1, 0.0); |
|
|
81 for (octave_idx_type ii(0); ii < npt; ii++) |
|
|
82 { |
|
|
83 basisfun (int(spu(ii)), u(ii), p, U, M); |
|
|
84 NuIkuk.insert (M, ii, 0); |
|
|
85 } // NuIkuk = basisfun (spu, u, p, U); |
|
|
86 |
|
12672
|
87 Matrix NvJkvk(v.numel (), q+1, 0.0); |
|
5683
|
88 for (octave_idx_type ii(0); ii < npt; ii++) |
|
|
89 { |
|
|
90 basisfun(int(spv(ii)), v(ii), q, V, N); |
|
|
91 NvJkvk.insert (N, ii, 0); |
|
|
92 } // NvJkvk = basisfun (spv, v, q, V); |
|
|
93 |
|
|
94 |
|
|
95 for (octave_idx_type k(0); k < npt; k++) |
|
|
96 for (octave_idx_type ii(0); ii < p+1; ii++) |
|
|
97 for (octave_idx_type jj(0); jj < q+1; jj++) |
|
9923
|
98 denom(k) += NuIkuk(k, ii) * NvJkvk(k, jj) * |
|
|
99 w(static_cast<octave_idx_type> (Ik(k, ii)), |
|
|
100 static_cast<octave_idx_type> (Jk(k, jj))); |
|
5683
|
101 |
|
|
102 |
|
|
103 for (octave_idx_type k(0); k < npt; k++) |
|
|
104 for (octave_idx_type ii(0); ii < p+1; ii++) |
|
|
105 for (octave_idx_type jj(0); jj < q+1; jj++) |
|
|
106 { |
|
9923
|
107 |
|
|
108 RIkJk(k, octave_idx_type(ii+(p+1)*jj)) = NuIkuk(k, ii) * NvJkvk(k, jj) * |
|
|
109 w(static_cast<octave_idx_type> (Ik(k, ii)), static_cast<octave_idx_type> (Jk(k, jj))) |
|
|
110 / denom(k); |
|
|
111 |
|
|
112 indIkJk(k, octave_idx_type(ii+(p+1)*jj))= Ik(k, ii) + (m+1) * Jk(k, jj) + 1; |
|
5683
|
113 } |
|
|
114 |
|
|
115 // for k=1:npt |
|
|
116 // [Jkb, Ika] = meshgrid(Jk(k, :), Ik(k, :)); |
|
|
117 // indIkJk(k, :) = sub2ind([m+1, n+1], Ika(:)+1, Jkb(:)+1); |
|
|
118 // wIkaJkb(1:p+1, 1:q+1) = reshape (w(indIkJk(k, :)), p+1, q+1); |
|
|
119 |
|
|
120 // NuIkukaNvJkvk(1:p+1, 1:q+1) = (NuIkuk(k, :).' * NvJkvk(k, :)); |
|
|
121 // RIkJk(k, :) = (NuIkukaNvJkvk .* wIkaJkb ./ sum(sum(NuIkukaNvJkvk .* wIkaJkb)))(:).'; |
|
|
122 // end |
|
|
123 |
|
|
124 retval(0) = RIkJk; // B = RIkJk; |
|
|
125 retval(1) = indIkJk; // N = indIkJk; |
|
|
126 |
|
|
127 } |
|
|
128 return retval; |
|
|
129 } |
