001/** 002 * BasicNurbsSurface -- A basic implementation of NURBS surface. 003 * 004 * Copyright (C) 2010-2015, by Joseph A. Huwaldt. All rights reserved. 005 * 006 * This library is free software; you can redistribute it and/or modify it under the terms 007 * of the GNU Lesser General Public License as published by the Free Software Foundation; 008 * either version 2.1 of the License, or (at your option) any later version. 009 * 010 * This library is distributed in the hope that it will be useful, but WITHOUT ANY 011 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A 012 * PARTICULAR PURPOSE. See the GNU Library General Public License for more details. 013 * 014 * You should have received a copy of the GNU Lesser General Public License along with 015 * this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - 016 * Suite 330, Boston, MA 02111-1307, USA. Or visit: http://www.gnu.org/licenses/lgpl.html 017 * 018 * This source file is based on, but heavily modified from, the LGPL jgeom (Geometry 019 * Library for Java) code by Samuel Gerber, Copyright (C) 2005. 020 */ 021package geomss.geom.nurbs; 022 023import geomss.geom.Point; 024import geomss.geom.Vector; 025import jahuwaldt.js.math.BinomialCoef; 026import java.text.MessageFormat; 027import java.util.List; 028import java.util.Objects; 029import static java.util.Objects.isNull; 030import javax.measure.converter.ConversionException; 031import javax.measure.quantity.Length; 032import javax.measure.unit.Unit; 033import javolution.context.ArrayFactory; 034import javolution.context.ObjectFactory; 035import javolution.context.StackContext; 036import javolution.lang.MathLib; 037import javolution.lang.ValueType; 038import javolution.util.FastTable; 039import javolution.xml.XMLFormat; 040import javolution.xml.stream.XMLStreamException; 041 042/** 043 * A basic implementation of a parametric NURBS surface. 044 * 045 * <p> Modified by: Joseph A. Huwaldt </p> 046 * 047 * @author Samuel Gerber, Date: June 14, 2010, Version 1.2. 048 * @version November 28, 2015 049 */ 050@SuppressWarnings({"serial", "CloneableImplementsClone"}) 051public final class BasicNurbsSurface extends NurbsSurface implements ValueType { 052 053 private static final double TOL = 1.0e-14; 054 055 private ControlPointNet _cnet; // The network (matrix) of control points. 056 private KnotVector _sKnots; // The knot vector in the u-direction. 057 private KnotVector _tKnots; // The knot vector in the v-direciton. 058 059 /* 060 * References: 061 * 1.) Michael E. Mortenson, "Geometric Modeling, Third Edition", ISBN:9780831132989, 2008. 062 * 2.) Piegl, L., Tiller, W., The Nurbs Book, 2nd Edition, Springer-Verlag, Berlin, 1997. 063 */ 064 065 /** 066 * Create a NURBS surface from the given control points, knots and degree. 067 * 068 * @param cps Matrix of control points: cps[u][v]. May not be null. 069 * @param p s-degree of the NURBS curve 070 * @param q t-degree of the NURBS curve 071 * @param sK Knot values in s-direction. May not be null. 072 * @param tK Knot values in t-direction. May not be null. 073 * @return A new BasicNurbsSurface from the given control points, knots and degree. 074 * @throws IllegalArgumentException if the knot vectors are not valid. 075 */ 076 public static BasicNurbsSurface newInstance(ControlPoint cps[][], int p, int q, double[] sK, double[] tK) { 077 if (cps.length < 1) 078 throw new IllegalArgumentException(RESOURCES.getString("noControlPointsErr")); 079 080 if (tK.length != q + cps.length + 1) 081 throw new IllegalArgumentException( 082 MessageFormat.format(RESOURCES.getString("srfWrongNumTKnotsErr"), 083 tK.length, q + cps.length + 1)); 084 085 if (sK.length != p + cps[0].length + 1) 086 throw new IllegalArgumentException( 087 MessageFormat.format(RESOURCES.getString("srfWrongNumSKnotsErr"), 088 sK.length, p + cps[0].length + 1)); 089 090 ControlPointNet cpNet = ControlPointNet.valueOf(cps); 091 KnotVector sKnots = KnotVector.newInstance(p, sK); 092 KnotVector tKnots = KnotVector.newInstance(q, tK); 093 094 return BasicNurbsSurface.newInstance(cpNet, sKnots, tKnots); 095 } 096 097 /** 098 * Generate a NURBS surface from the given control points and the given knot vectors. 099 * 100 * @param cpNet Matrix of control points. May not be null. 101 * @param sKnots Knot vector in the s-direction for the surface. May not be null. 102 * @param tKnots Knot vector in the t-direction for the surface. May not be null. 103 * @return A new BasicNurbsSurface from the given control points and the given knot 104 * vectors. 105 * @throws IllegalArgumentException if the knot vector is not consistent with the 106 * number of control points. 107 */ 108 public static BasicNurbsSurface newInstance(ControlPointNet cpNet, KnotVector sKnots, KnotVector tKnots) { 109 if (cpNet.size() < 1) 110 throw new IllegalArgumentException(RESOURCES.getString("noControlPointsErr")); 111 112 if (tKnots.length() != tKnots.getDegree() + cpNet.getNumberOfColumns() + 1) 113 throw new IllegalArgumentException( 114 MessageFormat.format(RESOURCES.getString("srfWrongNumTKnotsErr"), 115 tKnots.length(), tKnots.getDegree() + cpNet.getNumberOfColumns() + 1)); 116 117 if (sKnots.length() != sKnots.getDegree() + cpNet.getNumberOfRows() + 1) 118 throw new IllegalArgumentException( 119 MessageFormat.format(RESOURCES.getString("srfWrongNumSKnotsErr"), 120 sKnots.length(), sKnots.getDegree() + cpNet.getNumberOfRows() + 1)); 121 122 BasicNurbsSurface o = FACTORY.object(); 123 o._cnet = ControlPointNet.valueOf(cpNet); 124 o._sKnots = sKnots; 125 o._tKnots = tKnots; 126 127 return o; 128 } 129 130 /** 131 * Returns the number of child-elements that make up this geometry element. This 132 * implementation returns the number of control points in this NURBS surface. 133 * 134 * @return The number of child-elements that make up this geometry element. 135 */ 136 @Override 137 public int size() { 138 return _cnet.size(); 139 } 140 141 /** 142 * Returns the number of physical dimensions of the geometry element. 143 * 144 * @return The number of physical dimensions of the geometry element. 145 */ 146 @Override 147 public int getPhyDimension() { 148 return _cnet.get(0, 0).getPhyDimension(); 149 } 150 151 /** 152 * Return a matrix or network of control points for this surface. 153 * 154 * @return the ordered control points 155 */ 156 @Override 157 public ControlPointNet getControlPoints() { 158 return _cnet; 159 } 160 161 /** 162 * Return the control point matrix size in the s-direction (down a column of control 163 * points). 164 * 165 * @return The control point matrix size in the s-direction. 166 */ 167 @Override 168 public int getNumberOfRows() { 169 return _cnet.getNumberOfRows(); 170 } 171 172 /** 173 * Return the control point matrix size in the t-direction (across the columns of 174 * control points). 175 * 176 * @return The control point matrix size in the t-direction. 177 */ 178 @Override 179 public int getNumberOfColumns() { 180 return _cnet.getNumberOfColumns(); 181 } 182 183 /** 184 * Return the s-direction knot vector of this surface. 185 * 186 * @return The s-knot vector. 187 */ 188 @Override 189 public KnotVector getSKnotVector() { 190 return _sKnots; 191 } 192 193 /** 194 * Return the t-direction knot vector of this surface. 195 * 196 * @return The t-knot vector. 197 */ 198 @Override 199 public KnotVector getTKnotVector() { 200 return _tKnots; 201 } 202 203 /** 204 * Return the degree of the NURBS surface in the s-direction. 205 * 206 * @return degree of surface in s-direction 207 */ 208 @Override 209 public int getSDegree() { 210 return _sKnots.getDegree(); 211 } 212 213 /** 214 * Return the degree of the NURBS surface in the t-direction. 215 * 216 * @return degree of surface in t-direction 217 */ 218 @Override 219 public int getTDegree() { 220 return _tKnots.getDegree(); 221 } 222 223 /** 224 * Calculate a point on the surface for the given parametric position on the surface. 225 * 226 * @param s 1st parametric dimension distance to calculate a point for (0.0 to 1.0 227 * inclusive). 228 * @param t 2nd parametric dimension distance to calculate a point for (0.0 to 1.0 229 * inclusive). 230 * @return The calculated point on the surface at the specified parameter values. 231 * @throws IllegalArgumentException if there is any problem with the parameter values. 232 */ 233 @Override 234 public Point getRealPoint(double s, double t) { 235 validateParameter(s, t); 236 237 // Corner-point interpolation: S(0,0)=P(0,0); S(1,0)=P(n,0); S(0,1)=P(0,m) and S(1,1)=P(n,m). 238 if (parNearStart(s, TOL_ST)) { 239 if (parNearStart(t, TOL_ST)) 240 return _cnet.get(0, 0).getPoint(); 241 else if (parNearEnd(t, TOL_ST)) 242 return _cnet.get(0, _cnet.getNumberOfColumns() - 1).getPoint(); 243 } else if (parNearEnd(s, TOL_ST)) { 244 if (parNearStart(t, TOL_ST)) 245 return _cnet.get(_cnet.getNumberOfRows() - 1, 0).getPoint(); 246 else if (parNearEnd(t, TOL_ST)) 247 return _cnet.get(_cnet.getNumberOfRows() - 1, _cnet.getNumberOfColumns() - 1).getPoint(); 248 } 249 250 StackContext.enter(); 251 try { 252 // Algorithm: A4.3, Ref. 2. 253 /* 254 This routine calculates the following: 255 P(s,t) = sum_{i=1}^{N} sum_{j=1}^M ( Pij*Wij*Nik(s)Njk(t) ) / 256 sum_{i=1}^{N} sum_{j=1}^{M} ( Wij*Nik(s)*Njk(t) ) 257 */ 258 Unit<Length> unit = getUnit(); 259 int s_span = _sKnots.findSpan(s); 260 int p = _sKnots.getDegree(); 261 int sind = s_span - p; 262 int t_span = _tKnots.findSpan(t); 263 int q = _tKnots.getDegree(); 264 265 // Apply the basis functions to each degree of the surface in each direction. 266 // Calculate "sw" point as sum_{i=1}^N sum_{j=1}^M (Pij*Wij*Nik(s)*Njk(t)) and 267 // "sw" weight as sum_{i=1}^N sum_{j=1}^M (Wij*Nik(s)*Njk(t)). 268 269 // Get the basis function values in the s-direction (Nk(s)). 270 double[] Nks = _sKnots.basisFunctions(s_span, s); 271 272 // Get the basis function values in the t-direction (Nk(t)). 273 double[] Nkt = _tKnots.basisFunctions(t_span, t); 274 275 // Compute the composite control point. 276 ControlPoint sw = ControlPoint.newInstance(getPhyDimension(), unit); 277 for (int i = 0; i <= q; i++) { 278 ControlPoint temp = ControlPoint.newInstance(getPhyDimension(), unit); 279 int tind = t_span - q + i; 280 for (int k = 0; k <= p; ++k) { 281 ControlPoint pw = _cnet.get(sind + k, tind); 282 Point pwp = pw.getPoint(); 283 double Nksk = Nks[k]; 284 double pww = pw.getWeight() * Nksk; 285 pwp = pwp.times(pww); 286 Point tmpp = temp.getPoint(); 287 double tmpw = temp.getWeight(); 288 temp = ControlPoint.valueOf(tmpp.plus(pwp), tmpw + pww); 289 } 290 temp = temp.times(Nkt[i]); 291 sw = sw.plus(temp); 292 } 293 294 // Convert the control point to a geometric point. 295 // The "sw" weight is the sum of Wi*Ni,k calculated above. 296 Point outPoint = sw.getPoint(); 297 if (outPoint.normValue() > TOL) 298 outPoint = outPoint.divide(sw.getWeight()); 299 300 return StackContext.outerCopy(outPoint.to(unit)); 301 302 } finally { 303 StackContext.exit(); 304 } 305 } 306 307 /** 308 * Calculate all the derivatives from <code>0</code> to <code>grade</code> with 309 * respect to parametric s-position on the surface for the given parametric position 310 * on the surface, <code>d^{grade}p(s,t)/d^{grade}s</code>. 311 * <p> 312 * Example:<br> 313 * 1st derivative (grade = 1), this returns <code>[p(s,t), dp(s,t)/ds]</code>;<br> 314 * 2nd derivative (grade = 2), this returns <code>[p(s,t), dp(s,t)/ds, d^2p(s,t)/d^2s]</code>; etc. 315 * </p> 316 * 317 * @param s 1st parametric dimension distance to calculate derivative for (0.0 to 318 * 1.0 inclusive). 319 * @param t 2nd parametric dimension distance to calculate derivative for (0.0 to 320 * 1.0 inclusive). 321 * @param grade The maximum grade to calculate the u-derivatives for (1=1st 322 * derivative, 2=2nd derivative, etc) 323 * @param scaled Pass <code>true</code> for properly scaled derivatives or 324 * <code>false</code> if the magnitude of the derivative vector is not 325 * required -- this sometimes results in faster calculation times. 326 * @return A list of s-derivatives up to the specified grade of the surface at the 327 * specified parametric position. 328 * @throws IllegalArgumentException if the grade is < 0 or the parameter values are 329 * invalid. 330 */ 331 @Override 332 public List<Vector<Length>> getSDerivatives(double s, double t, int grade, boolean scaled) { 333 validateParameter(s, t); 334 335 // Handle roundoff 336 s = roundParNearEnds(s); 337 t = roundParNearEnds(t); 338 339 if (grade < 0) 340 throw new IllegalArgumentException(RESOURCES.getString("gradeLTZeroErr")); 341 342 int s_span = _sKnots.findSpan(s); 343 int p = _sKnots.getDegree(); 344 int sind = s_span - p; 345 int ds = MathLib.min(grade, p); 346 int t_span = _tKnots.findSpan(t); 347 int q = _tKnots.getDegree(); 348 int gradeP1 = grade + 1; 349 350 Vector[] cOut = Vector.allocateArray(gradeP1); // Created outside of StackContext. 351 StackContext.enter(); 352 try { 353 // Get the basis function derivative values in the s-direction (dNk(s)). 354 double[][] dNks = _sKnots.basisFunctionDerivatives(s_span, s, ds); 355 356 // Calculate z = sum_i^N sum_j^M ( Wij*dNik(s)*Njk(t) ) for each basis function derivative and 357 // pnts = sum_i^N sum_j^M ( Pij*Wij*dNik(s)*Njk(t) ) for each control point and basis function derivative. 358 int porder = p + 1; 359 int qorder = q + 1; 360 FastTable<double[]> zList = FastTable.newInstance(); 361 FastTable<Point[]> pntsList = FastTable.newInstance(); 362 for (int l = 0; l < qorder; l++) { 363 double[] z = newZeroedDoubleArray(gradeP1); 364 zList.add(z); 365 Point[] pnts = newZeroedPointArray(gradeP1); 366 pntsList.add(pnts); 367 int tind = t_span - q + l; 368 369 for (int k = 0; k <= ds; ++k) { 370 for (int i = 0; i < porder; i++) { 371 ControlPoint cpi = _cnet.get(sind + i, tind); 372 double Wi = cpi.getWeight(); 373 double WiNk = Wi * dNks[k][i]; 374 z[k] += WiNk; 375 pnts[k] = pnts[k].plus(cpi.getPoint().times(WiNk)); 376 } 377 } 378 } 379 380 // Get the basis function values in the t-direction (Nk(t)). 381 double[] Nkt = _tKnots.basisFunctions(t_span, t); 382 383 // Apply Nk(t) 384 for (int l = 0; l < qorder; l++) { 385 Point[] pnts = pntsList.get(l); 386 @SuppressWarnings("MismatchedReadAndWriteOfArray") 387 double[] z = zList.get(l); 388 389 for (int k = 0; k <= ds; ++k) { 390 double Nktl = Nkt[l]; 391 z[k] *= Nktl; 392 pnts[k] = pnts[k].times(Nktl); 393 } 394 } 395 396 // Sum up the lists of values collected above. 397 double[] z = sumDerivativeWeights(zList, grade); 398 Point[] pnts = sumPoints(pntsList, grade); 399 400 // Compute the weighted derivatives by stepping through a binomial expansion for each derivative. 401 Point[] c = binomial(pnts, z, gradeP1); 402 403 // Copy the final points out of the StackContext. 404 for (int i = gradeP1-1; i >= 0; --i) { 405 cOut[i] = StackContext.outerCopy(c[i].toGeomVector()); 406 } 407 408 } finally { 409 StackContext.exit(); 410 } 411 412 // Convert the array point points into a list of vectors. 413 FastTable<Vector<Length>> output = FastTable.newInstance(); 414 Point origin = Point.valueOf(cOut[0]); 415 for (int i = 0; i < gradeP1; ++i) { 416 Vector<Length> v = cOut[i]; 417 v.setOrigin(origin); 418 output.add(v); 419 } 420 421 Vector.recycleArray(cOut); 422 423 return output; 424 } 425 426 /** 427 * Calculate all the derivatives from <code>0</code> to <code>grade</code> with 428 * respect to parametric t-position on the surface for the given parametric position 429 * on the surface, <code>d^{grade}p(s,t)/d^{grade}t</code>. 430 * <p> 431 * Example:<br> 432 * 1st derivative (grade = 1), this returns <code>[p(s,t), dp(s,t)/dt]</code>;<br> 433 * 2nd derivative (grade = 2), this returns <code>[p(s,t), dp(s,t)/dt, d^2p(s,t)/d^2t]</code>; etc. 434 * </p> 435 * 436 * @param s 1st parametric dimension distance to calculate derivative for (0.0 to 437 * 1.0 inclusive). 438 * @param t 2nd parametric dimension distance to calculate derivative for (0.0 to 439 * 1.0 inclusive). 440 * @param grade The maximum grade to calculate the v-derivatives for (1=1st 441 * derivative, 2=2nd derivative, etc) 442 * @param scaled Pass <code>true</code> for properly scaled derivatives or 443 * <code>false</code> if the magnitude of the derivative vector is not 444 * required -- this sometimes results in faster calculation times. 445 * @return A list of t-derivatives up to the specified grade of the surface at the 446 * specified parametric position. 447 * @throws IllegalArgumentException if the grade is < 0 or the parameter values are 448 * invalid. 449 */ 450 @Override 451 public List<Vector<Length>> getTDerivatives(double s, double t, int grade, boolean scaled) { 452 validateParameter(s, t); 453 454 // Handle roundoff 455 s = roundParNearEnds(s); 456 t = roundParNearEnds(t); 457 458 if (grade < 0) 459 throw new IllegalArgumentException(RESOURCES.getString("gradeLTZeroErr")); 460 461 int s_span = _sKnots.findSpan(s); 462 int p = _sKnots.getDegree(); 463 int sind = s_span - p; 464 int t_span = _tKnots.findSpan(t); 465 int q = _tKnots.getDegree(); 466 int dt = MathLib.min(grade, q); 467 int gradeP1 = grade + 1; 468 469 Vector[] cOut = Vector.allocateArray(gradeP1); // Created outside of StackContext. 470 StackContext.enter(); 471 try { 472 // Get the basis function values for a B-Spline (unweighted derivatives). 473 double[] Nks = _sKnots.basisFunctions(s_span, s); 474 475 // Calculate z = sum_i^N sum_j^M ( Wij*Nik(s)*dNjk(t) ) for each basis function derivative and 476 // pnts = sum_i^N sum_j^M ( Pij*Wij*Nik(s)*dNjk(t) ) for each control point and basis function derivative. 477 int porder = p + 1; 478 int qorder = q + 1; 479 FastTable<double[]> zList = FastTable.newInstance(); 480 FastTable<Point[]> pntsList = FastTable.newInstance(); 481 for (int l = 0; l < qorder; l++) { 482 double[] z = newZeroedDoubleArray(gradeP1); 483 zList.add(z); 484 Point[] pnts = newZeroedPointArray(gradeP1); 485 pntsList.add(pnts); 486 int tind = t_span - q + l; 487 488 for (int k = 0; k <= dt; ++k) { 489 for (int i = 0; i < porder; i++) { 490 ControlPoint cpi = _cnet.get(sind + i, tind); 491 double Wi = cpi.getWeight(); 492 double WiNk = Wi * Nks[i]; 493 z[k] += WiNk; 494 pnts[k] = pnts[k].plus(cpi.getPoint().times(WiNk)); 495 } 496 } 497 } 498 499 // Get the basis function derivative values in the t-direction (dNk(t)). 500 double[][] dNkt = _tKnots.basisFunctionDerivatives(t_span, t, dt); 501 502 // Apply dNk(t) 503 for (int l = 0; l < qorder; l++) { 504 Point[] pnts = pntsList.get(l); 505 @SuppressWarnings("MismatchedReadAndWriteOfArray") 506 double[] z = zList.get(l); 507 508 for (int k = 0; k <= dt; ++k) { 509 double dNktl = dNkt[k][l]; 510 z[k] *= dNktl; 511 pnts[k] = pnts[k].times(dNktl); 512 } 513 } 514 515 // Sum up the lists of values collected above. 516 double[] z = sumDerivativeWeights(zList, grade); 517 Point[] pnts = sumPoints(pntsList, grade); 518 519 // Compute the weighted derivatives by stepping through a binomial expansion for each derivative. 520 Point[] c = binomial(pnts, z, gradeP1); 521 522 // Copy the final points out of the StackContext. 523 for (int i = gradeP1-1; i >= 0; --i) { 524 cOut[i] = StackContext.outerCopy(c[i].toGeomVector()); 525 } 526 527 } finally { 528 StackContext.exit(); 529 } 530 531 // Convert the array point points into a list of vectors. 532 FastTable<Vector<Length>> output = FastTable.newInstance(); 533 Point origin = Point.valueOf(cOut[0]); 534 for (int i = 0; i < gradeP1; ++i) { 535 Vector<Length> v = cOut[i]; 536 v.setOrigin(origin); 537 output.add(v); 538 } 539 540 Vector.recycleArray(cOut); 541 542 return output; 543 } 544 545 /** 546 * Calculate the twist vector (d^2P/(ds*dt) = d(dP/ds)/dt) for this surface at the 547 * specified position on this surface. 548 * 549 * @param s 1st parametric dimension distance to calculate twist vector for (0.0 to 550 * 1.0 inclusive). 551 * @param t 2nd parametric dimension distance to calculate twist vector for (0.0 to 552 * 1.0 inclusive). 553 * @return The twist vector of this surface at the specified parametric position. 554 * @throws IllegalArgumentException if the parameter values are invalid. 555 */ 556 @Override 557 public Vector<Length> getTwistVector(double s, double t) { 558 validateParameter(s, t); 559 560 // Handle roundoff 561 s = roundParNearEnds(s); 562 t = roundParNearEnds(t); 563 564 int s_span = _sKnots.findSpan(s); 565 int p = _sKnots.getDegree(); 566 int sind = s_span - p; 567 int t_span = _tKnots.findSpan(t); 568 int q = _tKnots.getDegree(); 569 570 Vector<Length> tvec; 571 Point origin; 572 StackContext.enter(); 573 try { 574 // Get the basis function derivative values in the s-direction (dNk(s)). 575 double[][] dNks = _sKnots.basisFunctionDerivatives(s_span, s, 1); 576 577 // Calculate z = sum_i^N sum_j^M ( Wij*dNik(s)*dNjk(t) ) for each basis function derivative and 578 // pnts = sum_i^N sum_j^M ( Pij*Wij*dNik(s)*dNjk(t) ) for each control point and basis function derivative. 579 int porder = p + 1; 580 int qorder = q + 1; 581 FastTable<double[]> zList = FastTable.newInstance(); 582 FastTable<Point[]> pntsList = FastTable.newInstance(); 583 for (int l = 0; l < qorder; l++) { 584 double[] z = newZeroedDoubleArray(1 + 1); 585 zList.add(z); 586 Point[] pnts = newZeroedPointArray(1 + 1); 587 pntsList.add(pnts); 588 int tind = t_span - q + l; 589 590 for (int k = 0; k <= 1; ++k) { 591 for (int i = 0; i < porder; i++) { 592 ControlPoint cpi = _cnet.get(sind + i, tind); 593 double Wi = cpi.getWeight(); 594 double WiNk = Wi * dNks[k][i]; 595 z[k] += WiNk; 596 pnts[k] = pnts[k].plus(cpi.getPoint().times(WiNk)); 597 } 598 } 599 } 600 601 // Get the basis function derivative values in the t-direction (dNk(t)). 602 double[][] dNkt = _tKnots.basisFunctionDerivatives(t_span, t, 1); 603 604 // Apply dNk(t) 605 for (int l = 0; l < qorder; l++) { 606 Point[] pnts = pntsList.get(l); 607 @SuppressWarnings("MismatchedReadAndWriteOfArray") 608 double[] z = zList.get(l); 609 610 for (int k = 0; k <= 1; ++k) { 611 double dNktl = dNkt[k][l]; 612 z[k] *= dNktl; 613 pnts[k] = pnts[k].times(dNktl); 614 } 615 } 616 617 // Sum up the lists of values collected above. 618 double[] z = sumDerivativeWeights(zList, 1); 619 Point[] pnts = sumPoints(pntsList, 1); 620 621 // Compute the weighted derivatives by stepping through a binomial expansion for each derivative. 622 Point[] c = binomial(pnts, z, 1 + 1); 623 624 // Copy out the vector origin point. 625 origin = StackContext.outerCopy(c[0]); 626 627 // Copy the required derivative point to the outer context. 628 tvec = StackContext.outerCopy(c[1].toGeomVector()); 629 630 } finally { 631 StackContext.exit(); 632 } 633 634 tvec.setOrigin(origin); 635 636 return tvec; 637 } 638 639 /** 640 * Compute the weighted derivatives by stepping through a binomial expansion for each 641 * derivative. e.g.: 642 * <pre> 643 * surface point, c(0) = p(0)/z(0) 644 * 1st derivative, c(1) = (p(1) - c(0)*z(1))/z(0), 645 * 2nd derivative, c(2) = (p(2) - c(0)*z(2) - 2*c(1)*z(1))/z(0), 646 * 3rd derivative, c(3) = (p(3) - c(0)*z(3) - 3*c(1)*z(2) - 3*c(2)*z(1))/z(0), 647 * etc 648 * </pre> 649 * 650 * @param pnts The list of weighted control points: p[0..grade+1]. 651 * @param z The derivative weights (Wi*Nik): z[0..grade+1]. 652 * @param gradeP1 The maximum grade to calculate the v-derivatives for (1=1st 653 * derivative, 2=2nd derivative, etc) + 1. 654 */ 655 private static Point[] binomial(Point[] pnts, double[] z, int gradeP1) { 656 657 BinomialCoef bin = BinomialCoef.newInstance(gradeP1); // Get the binomial coefficients. 658 659 Point[] c = Point.allocateArray(gradeP1); 660 for (int k = 0; k < gradeP1; ++k) { 661 Point v = pnts[k]; 662 for (int i = k; i > 0; --i) { 663 v = v.minus(c[k - i].times(bin.get(k, i) * z[i])); 664 } 665 if (v.normValue() > TOL) // Handles 0/0 == 0. 666 c[k] = v.divide(z[0]); 667 else 668 c[k] = v.times(0); 669 } 670 671 return c; 672 } 673 674 /** 675 * Return a new double array with the values zeroed. The returned instance is created 676 * using ArrayFactory.DOUBLES_FACTORY and may be recycled. 677 */ 678 private static double[] newZeroedDoubleArray(int size) { 679 double[] arr = ArrayFactory.DOUBLES_FACTORY.array(size); 680 for (int i = size-1; i >= 0; --i) { 681 arr[i] = 0; 682 } 683 return arr; 684 } 685 686 /** 687 * Return a new array filled with Point objects with their values zeroed. The returned 688 * instance is created using Point.allocateArray() and may be recycled. 689 * @see Point#allocateArray(int) 690 * @see Point#recycleArray(geomss.geom.Point[]) 691 */ 692 private Point[] newZeroedPointArray(int size) { 693 Point[] pnts = Point.allocateArray(size); 694 Unit unit = getUnit(); 695 int dim = getPhyDimension(); 696 for (int i = size-1; i >= 0; --i) 697 pnts[i] = Point.newInstance(dim, unit); 698 return pnts; 699 } 700 701 /** 702 * Sum up the list of point arrays into a single point array. 703 * 704 * @param pntsList A list of arrays of Point objects, each "grade+1" long. 705 * @param grade The length-1 of each Point array in pntsList. 706 */ 707 private Point[] sumPoints(List<Point[]> pntsList, int grade) { 708 int order = pntsList.size(); 709 710 Point[] pnts = newZeroedPointArray(grade + 1); 711 for (int l = 0; l < order; ++l) { 712 Point[] pntsTmp = pntsList.get(l); 713 714 for (int k = grade; k >= 0; --k) { 715 pnts[k] = pnts[k].plus(pntsTmp[k]); 716 } 717 } 718 719 return pnts; 720 } 721 722 /** 723 * Sum up the list of Wi*Nik arrays into a single array. 724 * 725 * @param arrList A list of double arrays, each "grade+1" long. 726 * @param grade The length-1 of each double array in arrList. 727 */ 728 private static double[] sumDerivativeWeights(List<double[]> arrList, int grade) { 729 int order = arrList.size(); 730 731 double[] z = newZeroedDoubleArray(grade + 1); 732 for (int l = 0; l < order; ++l) { 733 double[] ztmp = arrList.get(l); 734 735 for (int k = grade; k >= 0; --k) { 736 z[k] += ztmp[k]; 737 } 738 } 739 740 return z; 741 } 742 743 /** 744 * Return the coordinate point representing the minimum bounding box corner of this 745 * geometry element (e.g.: min X, min Y, min Z). 746 * 747 * @return The minimum bounding box coordinate for this geometry element. 748 * @throws IndexOutOfBoundsException if this list contains no geometry. 749 */ 750 @Override 751 public Point getBoundsMin() { 752 return _cnet.getBoundsMin(); 753 } 754 755 /** 756 * Return the coordinate point representing the maximum bounding box corner (e.g.: max 757 * X, max Y, max Z). 758 * 759 * @return The maximum bounding box coordinate for this geometry element. 760 * @throws IndexOutOfBoundsException if this list contains no elements. 761 */ 762 @Override 763 public Point getBoundsMax() { 764 return _cnet.getBoundsMax(); 765 } 766 767 /** 768 * Returns the unit in which the control points in this surface are stated. 769 * 770 * @return The unit in which the control points in this surface are stated. 771 */ 772 @Override 773 public Unit<Length> getUnit() { 774 return _cnet.getUnit(); 775 } 776 777 /** 778 * Returns the equivalent to this surface but stated in the specified unit. 779 * 780 * @param unit The length unit of the surface to be returned. May not be null. 781 * @return An equivalent to this surface but stated in the specified unit. 782 * @throws ConversionException if the the input unit is not a length unit. 783 */ 784 @Override 785 public BasicNurbsSurface to(Unit<Length> unit) throws ConversionException { 786 if (unit.equals(getUnit())) 787 return this; 788 789 ControlPointNet nCNet = _cnet.to(unit); 790 BasicNurbsSurface srf = BasicNurbsSurface.newInstance(nCNet, getSKnotVector(), getTKnotVector()); 791 return copyState(srf); // Copy over the super-class state for this object to the new one. 792 } 793 794 /** 795 * Return the equivalent of this surface converted to the specified number of physical 796 * dimensions. If the number of dimensions is greater than this element, then zeros 797 * are added to the additional dimensions. If the number of dimensions is less than 798 * this element, then the extra dimensions are simply dropped (truncated). If the new 799 * dimensions are the same as the dimension of this element, then this element is 800 * simply returned. 801 * 802 * @param newDim The dimension of the surface to return. 803 * @return The equivalent of this surface converted to the new dimensions. 804 */ 805 @Override 806 public BasicNurbsSurface toDimension(int newDim) { 807 if (getPhyDimension() == newDim) 808 return this; 809 810 ControlPointNet nCNet = _cnet.toDimension(newDim); 811 BasicNurbsSurface srf = BasicNurbsSurface.newInstance(nCNet, getSKnotVector(), getTKnotVector()); 812 813 return copyState(srf); // Copy over the super-class state for this object to the new one. 814 } 815 816 /** 817 * Returns a copy of this <code>BasicNurbsSurface</code> instance 818 * {@link javolution.context.AllocatorContext allocated} by the calling thread 819 * (possibly on the stack). 820 * 821 * @return an identical and independent copy of this object. 822 */ 823 @Override 824 public BasicNurbsSurface copy() { 825 return copyOf(this); 826 } 827 828 /** 829 * Return a copy of this object with any transformations or subranges removed 830 * (applied). 831 * 832 * @return A copy of this object with any transformations or subranges removed. 833 */ 834 @Override 835 public BasicNurbsSurface copyToReal() { 836 return copy(); 837 } 838 839 /** 840 * Compares this BasicNurbsSurface against the specified object for strict equality 841 * (same values and same units). 842 * 843 * @param obj the object to compare with. 844 * @return <code>true</code> if this object is identical to that object; 845 * <code>false</code> otherwise. 846 */ 847 @Override 848 public boolean equals(Object obj) { 849 if (this == obj) 850 return true; 851 if ((obj == null) || (obj.getClass() != this.getClass())) 852 return false; 853 854 BasicNurbsSurface that = (BasicNurbsSurface)obj; 855 return this._cnet.equals(that._cnet) 856 && this._sKnots.equals(that._sKnots) 857 && this._tKnots.equals(that._tKnots) 858 && super.equals(obj); 859 } 860 861 /** 862 * Returns the hash code for this parameter. 863 * 864 * @return the hash code value. 865 */ 866 @Override 867 public int hashCode() { 868 return 31*super.hashCode() + Objects.hash(_cnet, _sKnots, _tKnots); 869 } 870 871 /** 872 * Recycles a <code>BasicNurbsSurface</code> instance immediately (on the stack when 873 * executing in a StackContext). 874 * 875 * @param instance The instance to be recycled. 876 */ 877 public static void recycle(BasicNurbsSurface instance) { 878 FACTORY.recycle(instance); 879 } 880 881 /** 882 * Holds the default XML representation for this object. 883 */ 884 @SuppressWarnings("FieldNameHidesFieldInSuperclass") 885 protected static final XMLFormat<BasicNurbsSurface> XML = new XMLFormat<BasicNurbsSurface>(BasicNurbsSurface.class) { 886 @Override 887 public BasicNurbsSurface newInstance(Class<BasicNurbsSurface> cls, InputElement xml) throws XMLStreamException { 888 return FACTORY.object(); 889 } 890 891 @Override 892 public void read(InputElement xml, BasicNurbsSurface obj) throws XMLStreamException { 893 NurbsSurface.XML.read(xml, obj); // Call parent read. 894 895 ControlPointNet cpNet = xml.get("CPNet", ControlPointNet.class); 896 KnotVector sKnots = xml.get("sKnots", KnotVector.class); 897 KnotVector tKnots = xml.get("tKnots", KnotVector.class); 898 899 if (isNull(cpNet) || cpNet.size() < 1) 900 throw new XMLStreamException(RESOURCES.getString("noControlPointsErr")); 901 902 if (tKnots.length() != tKnots.getDegree() + cpNet.getNumberOfColumns() + 1) 903 throw new XMLStreamException( 904 MessageFormat.format(RESOURCES.getString("srfWrongNumTKnotsErr"), 905 tKnots.length(), tKnots.getDegree() + cpNet.getNumberOfColumns() + 1)); 906 907 if (sKnots.length() != sKnots.getDegree() + cpNet.getNumberOfRows() + 1) 908 throw new XMLStreamException( 909 MessageFormat.format(RESOURCES.getString("srfWrongNumSKnotsErr"), 910 sKnots.length(), sKnots.getDegree() + cpNet.getNumberOfRows() + 1)); 911 912 obj._cnet = cpNet; 913 obj._sKnots = sKnots; 914 obj._tKnots = tKnots; 915 } 916 917 @Override 918 public void write(BasicNurbsSurface obj, OutputElement xml) throws XMLStreamException { 919 NurbsSurface.XML.write(obj, xml); // Call parent write. 920 921 xml.add(obj._cnet, "CPNet", ControlPointNet.class); 922 xml.add(obj._sKnots, "sKnots", KnotVector.class); 923 xml.add(obj._tKnots, "tKnots", KnotVector.class); 924 925 } 926 }; 927 928 /////////////////////// 929 // Factory creation. // 930 /////////////////////// 931 private BasicNurbsSurface() { } 932 933 @SuppressWarnings("unchecked") 934 private static final ObjectFactory<BasicNurbsSurface> FACTORY = new ObjectFactory<BasicNurbsSurface>() { 935 @Override 936 protected BasicNurbsSurface create() { 937 return new BasicNurbsSurface(); 938 } 939 940 @Override 941 protected void cleanup(BasicNurbsSurface obj) { 942 obj.reset(); 943 obj._cnet = null; 944 obj._sKnots = null; 945 obj._tKnots = null; 946 } 947 }; 948 949 @SuppressWarnings("unchecked") 950 private static BasicNurbsSurface copyOf(BasicNurbsSurface original) { 951 BasicNurbsSurface o = FACTORY.object(); 952 o._cnet = original._cnet.copy(); 953 o._sKnots = original._sKnots.copy(); 954 o._tKnots = original._tKnots.copy(); 955 original.copyState(o); 956 return o; 957 } 958 959}