Mercurial > forge
annotate extra/nurbs/src/tbasisfun.cc @ 12672:59e8aae64812 octave-forge
prepare for release
| author | cdf |
|---|---|
| date | Mon, 17 Aug 2015 10:23:44 +0000 |
| parents | 37d08939bb7b |
| children |
| rev | line source |
|---|---|
| 6958 | 1 /* Copyright (C) 2009 Carlo de Falco |
|
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
2 Copyright (C) 2012 Rafael Vazquez |
| 6958 | 3 |
| 4 This program is free software: you can redistribute it and/or modify | |
| 5 it under the terms of the GNU General Public License as published by | |
| 11634 | 6 the Free Software Foundation, either version 3 of the License, or |
| 6958 | 7 (at your option) any later version. |
| 6052 | 8 |
| 6958 | 9 This program is distributed in the hope that it will be useful, |
| 10 but WITHOUT ANY WARRANTY; without even the implied warranty of | |
| 11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
| 12 GNU General Public License for more details. | |
| 6052 | 13 |
| 6958 | 14 You should have received a copy of the GNU General Public License |
| 15 along with this program. If not, see <http://www.gnu.org/licenses/>. | |
| 6052 | 16 */ |
| 17 | |
| 18 #include <octave/oct.h> | |
| 19 #include <iostream> | |
| 20 | |
| 9610 | 21 void onebasisfun__ (double u, octave_idx_type p, RowVector U, double *N) |
| 6052 | 22 { |
| 9610 | 23 *N = 0.0; |
|
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
24 if ((u <= U.min ()) || ( u > U.max ())) |
| 9610 | 25 return; |
| 6052 | 26 else if (p == 0) |
| 9610 | 27 { |
| 28 *N = 1.0; | |
| 29 return; | |
| 9484 | 30 } |
| 9610 | 31 else if (p == 1) |
| 32 { | |
| 33 if (u < U(1)) | |
| 34 { | |
| 35 *N = (u - U(0)) / (U(1) - U(0)); | |
| 36 return; | |
| 37 } | |
| 38 else | |
| 39 { | |
| 40 *N = (U(2) - u) / (U(2) - U(1)); | |
| 41 return; | |
| 42 } | |
| 9484 | 43 } |
| 9610 | 44 else if (p == 2) |
| 45 { | |
| 46 double ln = u - U(0); | |
| 47 double dn = U(3) - u; | |
| 48 double ld = U(2) - U(0); | |
| 49 double dd = U(3) - U(1); | |
| 50 if (u < U(1)) | |
| 51 { | |
| 52 *N = ln*ln / (ld * (U(1) - U(0))); | |
| 53 return; | |
| 54 } | |
| 55 else if (u > U(2)) | |
| 56 { | |
| 57 *N = dn*dn / (dd * (U(3) - U(2))); | |
| 58 return; | |
| 59 } | |
| 60 else | |
| 61 { | |
| 62 if (ld != 0) | |
| 63 *N += ln * (U(2) - u) / ((U(2) - U(1)) * ld); | |
| 64 | |
| 65 if (dd != 0) | |
| 66 *N += dn * (u - U(1)) / ((U(2) - U(1)) * dd); | |
| 67 return; | |
| 68 } | |
| 9484 | 69 } |
| 6052 | 70 |
| 71 double ln = u - U(0); | |
| 12672 | 72 double ld = U(U.numel () - 2) - U(0); |
| 6052 | 73 if (ld != 0) |
| 9610 | 74 { |
| 75 double tmp; | |
| 12672 | 76 onebasisfun__ (u, p-1, U.extract (0, U.numel () - 2), &tmp); |
| 9610 | 77 *N += ln * tmp / ld; |
| 78 } | |
| 12672 | 79 double dn = U(U.numel () - 1) - u; |
| 80 double dd = U(U.numel () - 1) - U(1); | |
| 6052 | 81 if (dd != 0) |
| 9610 | 82 { |
| 83 double tmp; | |
| 12672 | 84 onebasisfun__ (u, p-1, U.extract (1, U.numel () - 1), &tmp); |
| 9610 | 85 *N += dn * tmp / dd; |
| 86 } | |
| 87 return; | |
| 6052 | 88 } |
| 89 | |
| 10766 | 90 void onebasisfun__ (double u, double p, RowVector U, double *N) |
| 91 { onebasisfun__ (u, static_cast<octave_idx_type> (p), U, N); } | |
|
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
92 |
| 9610 | 93 void onebasisfunder__ (double u, octave_idx_type p, RowVector U, double *N, double *Nder) |
|
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
94 { |
| 9610 | 95 double aux; |
| 10777 | 96 *N = 0.0; *Nder = 0.0; |
| 97 if ((u <= U.min ()) || ( u > U.max ())) | |
| 9610 | 98 return; |
| 10777 | 99 else if (p == 0) |
| 100 { | |
| 101 *N = 1.0; | |
| 102 *Nder = 0.0; | |
| 9610 | 103 return; |
| 10777 | 104 } |
| 105 else { | |
| 106 double ln = u - U(0); | |
| 12672 | 107 double ld = U(U.numel () - 2) - U(0); |
| 10777 | 108 |
| 109 if (ld != 0) | |
| 9610 | 110 { |
| 12672 | 111 onebasisfun__ (u, p-1, U.extract (0, U.numel () - 2), &aux); |
| 10777 | 112 aux = aux / ld; |
| 113 *N += ln * aux; | |
| 114 *Nder += p * aux; | |
| 115 } | |
|
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
116 |
| 12672 | 117 double dn = U(U.numel () - 1) - u; |
| 118 double dd = U(U.numel () - 1) - U(1); | |
| 10777 | 119 |
| 120 if (dd != 0) | |
| 9610 | 121 { |
| 12672 | 122 onebasisfun__ (u, p-1, U.extract (1, U.numel () - 1), &aux); |
| 9610 | 123 aux = aux / dd; |
| 124 *N += dn *aux; | |
| 125 *Nder -= p * aux; | |
| 126 } | |
|
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
127 } |
|
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
128 } |
| 6052 | 129 |
| 9405 | 130 DEFUN_DLD(tbasisfun, args, nargout,"\ |
|
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
131 TBASISFUN: Compute a B- or T-Spline basis function, and its derivatives, from its local knot vector.\n\ |
| 6464 | 132 \n\ |
| 6052 | 133 usage:\n\ |
| 134 \n\ | |
|
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
135 [N, Nder] = tbasisfun (u, p, U)\n\ |
|
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
136 [N, Nder] = tbasisfun ([u; v], [p q], {U, V})\n\ |
|
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
137 [N, Nder] = tbasisfun ([u; v; w], [p q r], {U, V, W})\n\ |
| 6052 | 138 \n\ |
| 139 INPUT:\n\ | |
| 140 u or [u; v] : points in parameter space where the basis function is to be\n\ | |
| 141 evaluated \n\ | |
| 142 \n\ | |
| 143 U or {U, V} : local knot vector\n\ | |
| 144 \n\ | |
|
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
145 p or [p q] : polynomial order of the basis function\n\ |
| 6052 | 146 \n\ |
| 147 OUTPUT:\n\ | |
|
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
148 N : basis function evaluated at the given parametric points\n\ |
|
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
149 Nder : gradient of the basis function evaluated at the given points\n") |
| 6052 | 150 |
| 151 { | |
| 152 | |
| 7204 | 153 octave_value_list retval; |
| 6052 | 154 Matrix u = args(0).matrix_value (); |
| 155 | |
|
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
156 RowVector N(u.cols ()); |
| 9610 | 157 double *Nptr = N.fortran_vec (); |
| 158 | |
| 6052 | 159 if (! args(2).is_cell ()) |
| 160 { | |
| 161 | |
| 162 double p = args(1).idx_type_value (); | |
| 163 RowVector U = args(2).row_vector_value (true, true); | |
| 164 assert (U.numel () == p+2); | |
| 165 | |
| 9610 | 166 if (nargout == 1) |
| 167 for (octave_idx_type ii = 0; ii < u.numel (); ii++) | |
| 168 onebasisfun__ (u(ii), p, U, &(Nptr[ii])); | |
| 6052 | 169 |
| 9610 | 170 if (nargout == 2) |
| 171 { | |
| 172 RowVector Nder(u.cols ()); | |
| 173 double *Nderptr = Nder.fortran_vec (); | |
|
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
174 |
| 9610 | 175 for (octave_idx_type ii=0; ii<u.numel (); ii++) |
| 176 onebasisfunder__ (u(ii), p, U, &(Nptr[ii]), &(Nderptr[ii])); | |
|
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
177 |
| 9610 | 178 retval(1) = Nder; |
| 179 } | |
| 180 } | |
| 181 else | |
| 182 { | |
| 183 RowVector p = args(1).row_vector_value (); | |
| 12672 | 184 if (p.numel () == 2) |
| 9610 | 185 { |
| 186 Cell C = args(2).cell_value (); | |
| 187 RowVector U = C(0).row_vector_value (true, true); | |
| 188 RowVector V = C(1).row_vector_value (true, true); | |
| 189 | |
| 190 if (nargout == 1) | |
| 191 { | |
| 192 for (octave_idx_type ii=0; ii<u.cols (); ii++) | |
| 193 { | |
| 194 double Nu, Nv; | |
| 195 onebasisfun__ (u(0, ii), octave_idx_type(p(0)), U, &Nu); | |
| 196 onebasisfun__ (u(1, ii), octave_idx_type(p(1)), V, &Nv); | |
| 197 Nptr[ii] = Nu * Nv; | |
| 198 } | |
| 199 } | |
| 200 else if (nargout == 2) | |
| 201 { | |
| 202 double Nu, Nv, Ndu, Ndv; | |
| 203 Matrix Nder (2, u.cols()); | |
| 204 double *Nderptr = Nder.fortran_vec (); | |
| 205 for (octave_idx_type ii = 0; ii < u.cols (); ii++) | |
| 206 { | |
| 207 onebasisfunder__ (u(0, ii), octave_idx_type(p(0)), U, &Nu, &Ndu); | |
| 208 onebasisfunder__ (u(1, ii), octave_idx_type(p(1)), V, &Nv, &Ndv); | |
| 209 Nptr[ii] = Nu * Nv; | |
| 210 | |
| 211 Nderptr[0 + (2 * ii)] = Ndu * Nv; | |
| 212 Nderptr[1 + (2 * ii)] = Ndv * Nu; | |
| 213 } | |
| 214 retval(1) = Nder; | |
| 215 } | |
| 216 | |
| 217 } | |
| 12672 | 218 else if (p.numel () == 3) |
| 9610 | 219 { |
| 220 Cell C = args(2).cell_value (); | |
| 221 RowVector U = C(0).row_vector_value (true, true); | |
| 222 RowVector V = C(1).row_vector_value (true, true); | |
| 223 RowVector W = C(2).row_vector_value (true, true); | |
| 224 | |
| 225 if (nargout == 1) | |
| 226 { | |
| 227 for (octave_idx_type ii = 0; ii < u.cols (); ii++) | |
| 228 { | |
| 229 double Nu, Nv, Nw; | |
| 230 onebasisfun__ (u(0, ii), octave_idx_type(p(0)), U, &Nu); | |
| 231 onebasisfun__ (u(1, ii), octave_idx_type(p(1)), V, &Nv); | |
| 232 onebasisfun__ (u(2, ii), octave_idx_type(p(2)), W, &Nw); | |
| 233 Nptr[ii] = Nu * Nv * Nw; | |
| 234 } | |
| 235 } | |
| 236 else if (nargout == 2) | |
| 237 { | |
| 238 double Nu, Nv, Nw, Ndu, Ndv, Ndw; | |
| 239 Matrix Nder (3, u.cols()); | |
| 240 double *Nderptr = Nder.fortran_vec (); | |
| 241 for (octave_idx_type ii=0; ii<u.cols (); ii++) | |
| 242 { | |
| 243 onebasisfunder__ (u(0, ii), octave_idx_type(p(0)), U, &Nu, &Ndu); | |
| 244 onebasisfunder__ (u(1, ii), octave_idx_type(p(1)), V, &Nv, &Ndv); | |
| 245 onebasisfunder__ (u(2, ii), octave_idx_type(p(2)), W, &Nw, &Ndw); | |
| 246 Nptr[ii] = Nu * Nv * Nw; | |
| 247 Nderptr[0 + (3 * ii)] = Ndu * Nv * Nw; | |
| 248 Nderptr[1 + (3 * ii)] = Ndv * Nu * Nw; | |
| 249 Nderptr[2 + (3 * ii)] = Ndw * Nu * Nv; | |
| 250 } | |
| 251 retval(1) = Nder; | |
|
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
252 } |
| 9610 | 253 |
| 9605 | 254 } |
|
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
255 } |
|
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
256 retval(0) = octave_value (N); |
| 7204 | 257 return retval; |
| 6052 | 258 } |
| 9406 | 259 |
| 260 /* | |
| 261 | |
| 262 %!demo | |
| 263 %! U = {[0 0 1/2 1 1], [0 0 0 1 1]}; | |
| 264 %! p = [3, 3]; | |
| 265 %! [X, Y] = meshgrid (linspace(0, 1, 30)); | |
| 266 %! u = [X(:), Y(:)]'; | |
| 267 %! N = tbasisfun (u, p, U); | |
| 268 %! surf (X, Y, reshape (N, size(X))) | |
| 269 %! title('Basis function associated to a local knot vector') | |
| 270 %! hold off | |
| 271 | |
| 272 %!test | |
| 273 %! U = [0 1/2 1]; | |
| 274 %! p = 1; | |
| 275 %! u = [0.3 0.4 0.6 0.7]; | |
| 276 %! [N, Nder] = tbasisfun (u, p, U); | |
| 277 %! assert (N, [0.6 0.8 0.8 0.6], 1e-12); | |
| 278 %! assert (Nder, [2 2 -2 -2], 1e-12); | |
| 279 | |
| 280 %!test | |
| 281 %! U = {[0 1/2 1] [0 1/2 1]}; | |
| 282 %! p = [1 1]; | |
| 283 %! u = [0.3 0.4 0.6 0.7; 0.3 0.4 0.6 0.7]; | |
| 284 %! [N, Nder] = tbasisfun (u, p, U); | |
| 285 %! assert (N, [0.36 0.64 0.64 0.36], 1e-12); | |
| 286 %! assert (Nder, [1.2 1.6 -1.6 -1.2; 1.2 1.6 -1.6 -1.2], 1e-12); | |
| 287 | |
| 288 %!test | |
| 289 %! U = {[0 1/2 1] [0 1/2 1] [0 1/2 1]}; | |
| 290 %! p = [1 1 1]; | |
| 291 %! u = [0.4 0.4 0.6 0.6; 0.4 0.4 0.6 0.6; 0.4 0.6 0.4 0.6]; | |
| 292 %! [N, Nder] = tbasisfun (u, p, U); | |
| 293 %! assert (N, [0.512 0.512 0.512 0.512], 1e-12); | |
| 294 %! assert (Nder, [1.28 1.28 -1.28 -1.28; 1.28 1.28 -1.28 -1.28; 1.28 -1.28 1.28 -1.28], 1e-12); | |
| 295 | |
| 296 */ |