001/** 002 * CurveFactory -- A factory for creating NURBS curves. 003 * 004 * Copyright (C) 2009-2025, 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 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, NurbsCreator.java, by Samuel Gerber, Copyright (C) 2005. 020 */ 021package geomss.geom.nurbs; 022 023import geomss.geom.*; 024import jahuwaldt.js.param.Parameter; 025import jahuwaldt.tools.math.MathTools; 026import java.text.MessageFormat; 027import java.util.List; 028import static java.util.Objects.isNull; 029import static java.util.Objects.nonNull; 030import static java.util.Objects.requireNonNull; 031import java.util.ResourceBundle; 032import javax.measure.quantity.Angle; 033import javax.measure.quantity.Dimensionless; 034import javax.measure.quantity.Length; 035import javax.measure.unit.SI; 036import javax.measure.unit.Unit; 037import javolution.context.ArrayFactory; 038import javolution.context.ImmortalContext; 039import javolution.context.StackContext; 040import static javolution.lang.MathLib.*; 041import javolution.util.FastTable; 042import org.apache.commons.math3.linear.*; 043import org.jscience.mathematics.number.Float64; 044import org.jscience.mathematics.vector.Float64Vector; 045 046/** 047 * A collection of methods for creating NURBS curves. 048 * 049 * <p> Modified by: Joseph A. Huwaldt </p> 050 * 051 * @author Samuel Gerber, Date: May 14, 2009, Version 1.0. 052 * @version February 17, 2025 053 */ 054public final class CurveFactory { 055 056 /** 057 * The resource bundle for this package. 058 */ 059 private static final ResourceBundle RESOURCES = geomss.geom.AbstractGeomElement.RESOURCES; 060 061 /** 062 * The length threshold for a set of points being considered a degenerate curve. 063 */ 064 private static final double EPSC = 1e-10; 065 066 // The knot vector for a semi-circle. 067 private static final double[] SC_KNOTS = {0, 0, 0, 0.5, 1, 1, 1}; 068 069 // The knot vector for a full circle or ellipse. 070 private static final double[] FC_KNOTS = {0, 0, 0, 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 1, 1, 1}; 071 072 // A generic geometry space tolerance 073 private static final Parameter<Length> GTOL; 074 075 static { 076 ImmortalContext.enter(); 077 try { 078 GTOL = Parameter.valueOf(1e-8, SI.METER); 079 } finally { 080 ImmortalContext.exit(); 081 } 082 } 083 084 private CurveFactory() { } 085 086 /** 087 * Create a degenerate {@link NurbsCurve} of the specified degree that represents a 088 * point in space. 089 * 090 * @param degree The degree of the curve to create. 091 * @param p0 The location for the degenerate curve. May not be null. 092 * @return A BasicNurbsCurve that represents a point. 093 */ 094 public static BasicNurbsCurve createPoint(int degree, GeomPoint p0) { 095 requireNonNull(p0); 096 StackContext.enter(); 097 try { 098 // Create an appropriate knot vector. 099 FastTable<Float64> kvList = FastTable.newInstance(); 100 for (int i = 0; i <= degree; ++i) { 101 kvList.add(Float64.ZERO); 102 } 103 for (int i = 0; i <= degree; ++i) { 104 kvList.add(Float64.ONE); 105 } 106 107 KnotVector kv = KnotVector.newInstance(degree, Float64Vector.valueOf(kvList)); 108 109 // Create an appropriate control point list. 110 ControlPoint cp0 = ControlPoint.valueOf(p0.immutable(), 1); 111 FastTable<ControlPoint> cps = FastTable.newInstance(); 112 for (int i = 0; i < kv.length() - degree - 1; ++i) { 113 cps.add(cp0); 114 } 115 116 BasicNurbsCurve crv = BasicNurbsCurve.newInstance(cps, kv); 117 return StackContext.outerCopy(crv); 118 119 } finally { 120 StackContext.exit(); 121 } 122 } 123 124 /** 125 * Create a degenerate {@link NurbsCurve} of the specified degree that represents a 126 * point in space with the indicated number of co-located control points. 127 * 128 * @param degree The degree of the curve to create. 129 * @param numCP The number of control points to use in the degenerate curve. 130 * @param p0 The location for the degenerate curve. May not be null. 131 * @return A BasicNurbsCurve that represents a point. 132 */ 133 public static BasicNurbsCurve createPoint(int degree, int numCP, GeomPoint p0) { 134 requireNonNull(p0); 135 StackContext.enter(); 136 try { 137 // Use an even parameter spacing. 138 double[] uk = ArrayFactory.DOUBLES_FACTORY.array(numCP); 139 int numCPm1 = numCP - 1; 140 for (int i = 0; i < numCP; ++i) { 141 double s = (double)(i) / numCPm1; 142 uk[i] = s; 143 } 144 145 // Create an appropriate knot vector. 146 KnotVector uKnots = buildInterpKnotVector(uk, numCP, degree); 147 148 // Create an appropriate control point list. 149 ControlPoint cp0 = ControlPoint.valueOf(p0.immutable(), 1); 150 FastTable<ControlPoint> cps = FastTable.newInstance(); 151 for (int i = 0; i < numCP; ++i) { 152 cps.add(cp0); 153 } 154 155 BasicNurbsCurve crv = BasicNurbsCurve.newInstance(cps, uKnots); 156 return StackContext.outerCopy(crv); 157 158 } finally { 159 StackContext.exit(); 160 } 161 } 162 163 /** 164 * Create a 1-degree straight line {@link NurbsCurve} that represents the shortest 165 * distance, in Euclidean space between the two input points. 166 * 167 * @param p0 The start of the line. May not be null. 168 * @param p1 The end of the line. May not be null. 169 * @return A BasicNurbsCurve that represents a line between the input points. 170 */ 171 public static BasicNurbsCurve createLine(GeomPoint p0, GeomPoint p1) { 172 StackContext.enter(); 173 try { 174 // Make both points have the dimension of the highest dimension point. 175 int numDims = Math.max(p0.getPhyDimension(), p1.getPhyDimension()); 176 p0 = p0.toDimension(numDims); 177 p1 = p1.toDimension(numDims); 178 179 FastTable<ControlPoint> cps = FastTable.newInstance(); 180 cps.add(ControlPoint.valueOf(p0.immutable(), 1)); 181 cps.add(ControlPoint.valueOf(p1.immutable(), 1)); 182 KnotVector kv = KnotVector.newInstance(1, 0., 0., 1., 1.); 183 184 BasicNurbsCurve crv = BasicNurbsCurve.newInstance(cps, kv); 185 return StackContext.outerCopy(crv); 186 187 } finally { 188 StackContext.exit(); 189 } 190 } 191 192 /** 193 * Create a straight line {@link NurbsCurve} that represents the shortest distance, in 194 * Euclidean space between the two input points with the specified degree. 195 * 196 * @param degree The degree of the NurbsCurve to return. 197 * @param p0 The start of the line. May not be null. 198 * @param p1 The end of the line. May not be null. 199 * @return A BasicNurbsCurve of the specified degree that represents a line between 200 * the input points. 201 */ 202 public static BasicNurbsCurve createLine(int degree, GeomPoint p0, GeomPoint p1) { 203 requireNonNull(p0); 204 requireNonNull(p1); 205 return (BasicNurbsCurve)CurveUtils.elevateToDegree(createLine(p0, p1), degree); 206 } 207 208 /** 209 * Create a semi-circle {@link NurbsCurve} with the specified normal vector around the 210 * given origin with radius <code>r</code>. The curve returned has a control polygon 211 * which has 4 points and the shape of a rectangle. 212 * 213 * @param o Origin to create a semi-circle around. May not be null. 214 * @param r Radius of the semi-circle May not be null. 215 * @param n A unit plane normal vector for the plane containing the circle. If 216 * <code>null</code> is passed, a default normal vector that points along the 217 * Z axis is used. 218 * @return A BasicNurbsCurve that represents a semi-circle 219 * @throws DimensionException if the origin point or normal do not have at least 2 220 * physical dimensions. 221 */ 222 public static BasicNurbsCurve createSemiCircle(GeomPoint o, Parameter<Length> r, GeomVector<Dimensionless> n) 223 throws DimensionException { 224 225 requireNonNull(r); 226 int numDims = o.getPhyDimension(); 227 if (nonNull(n) && n.getPhyDimension() > numDims) 228 numDims = n.getPhyDimension(); 229 if (numDims < 2) 230 throw new DimensionException( 231 MessageFormat.format(RESOURCES.getString("dimensionNotAtLeast2"), 232 "input origin or normal", numDims)); 233 234 StackContext.enter(); 235 try { 236 if (nonNull(n)) { 237 o = o.toDimension(numDims); 238 n = n.toDimension(numDims); 239 } 240 241 // Make sure the normal vector exists and is a unit vector. 242 if (numDims == 2) 243 n = null; 244 else if (isNull(n)) 245 n = GeomUtil.makeZVector(numDims); 246 else 247 n = n.toUnitVector(); 248 249 // Create the yhat and xhat vectors (orthogonal Y and X in the plane of the circle). 250 Vector<Dimensionless> yhat = GeomUtil.calcYHat(n); 251 Vector<Dimensionless> xhat = GeomUtil.calcXHat(n, yhat); 252 253 // Create a list of control points. 254 FastTable<ControlPoint> cpList = FastTable.newInstance(); 255 cpList.setSize(4); 256 257 double rValue = r.getValue(SI.METER); 258 259 // Set 0 (bottom) and 3 (top) control points. 260 // p = 0,r 261 Point p = Point.valueOf(yhat); 262 p = p.times(rValue); 263 cpList.set(3, ControlPoint.valueOf(o.plus(p), 1)); 264 cpList.set(0, ControlPoint.valueOf(o.plus(p.opposite()), 1)); // p = 0,-r 265 266 // Set 1 (bottom right) and 2 (top right) control points. 267 // p = r, r 268 p = Point.valueOf(xhat.plus(yhat)); 269 p = p.times(rValue); 270 cpList.set(2, ControlPoint.valueOf(o.plus(p), 0.5)); 271 272 // p = r, -r 273 p = Point.valueOf(xhat.plus(yhat.opposite())); 274 p = p.times(rValue); 275 cpList.set(1, ControlPoint.valueOf(o.plus(p), 0.5)); 276 277 // Create the NURBS curve. 278 KnotVector kv = KnotVector.newInstance(2, SC_KNOTS); 279 return StackContext.outerCopy(BasicNurbsCurve.newInstance(cpList, kv)); 280 281 } finally { 282 StackContext.exit(); 283 } 284 } 285 286 /** 287 * Create a full-circle {@link NurbsCurve} around the given origin/center with radius 288 * <code>r</code>. The NurbsCurve returned has a control polygon which has 9 control 289 * points and the shape of a square. 290 * 291 * @param o Origin or center to create the full-circle around. May not be null. 292 * @param r Radius of the full-circle. May not be null. 293 * @param n A unit plane normal vector for the plane containing the circle. If 294 * <code>null</code> is passed, a default normal vector that points along the 295 * Z axis is used. 296 * @return A BasicNurbsCurve for a full-circle 297 * @throws DimensionException if the origin point or normal do not have at least 2 298 * physical dimensions. 299 */ 300 public static BasicNurbsCurve createCircle(GeomPoint o, Parameter<Length> r, GeomVector<Dimensionless> n) 301 throws DimensionException { 302 requireNonNull(r); 303 int numDims = o.getPhyDimension(); 304 if (nonNull(n) && n.getPhyDimension() > numDims) 305 numDims = n.getPhyDimension(); 306 if (numDims < 2) 307 throw new DimensionException( 308 MessageFormat.format(RESOURCES.getString("dimensionNotAtLeast2"), 309 "input origin or normal", numDims)); 310 311 StackContext.enter(); 312 try { 313 if (nonNull(n)) { 314 o = o.toDimension(numDims); 315 n = n.toDimension(numDims); 316 } 317 318 // Make sure the normal vector exists and is a unit vector. 319 if (numDims == 2) 320 n = null; 321 else if (isNull(n)) 322 n = GeomUtil.makeZVector(numDims); 323 else 324 n = n.toUnitVector(); 325 326 // Create the yhat and xhat vectors (orthogonal Y and X in the plane of the circle). 327 Vector<Dimensionless> yhat = GeomUtil.calcYHat(n); 328 Vector<Dimensionless> xhat = GeomUtil.calcXHat(n, yhat); 329 330 // Now just create a circle from the center, radius, and direction vectors. 331 BasicNurbsCurve crv = createEllipseCrv(o, r, r, xhat, yhat); 332 return StackContext.outerCopy(crv); 333 334 } finally { 335 StackContext.exit(); 336 } 337 } 338 339 /** 340 * Create a full-circle {@link NurbsCurve} that passes through the input (not 341 * co-linear) points. The curve will be parameterized with 0 at the starting point 342 * with the parameterization increasing toward the last point. The NurbsCurve returned 343 * has a control polygon which has 9 control points and the shape of square. 344 * 345 * @param p1 The 1st of the points that define the circle (must not be colinear with 346 * the other two). May not be null. 347 * @param p2 The 2nd of the points that define the circle (must not be colinear with 348 * the other two). May not be null. 349 * @param p3 The 3rd of the points that define the circle (must not be colinear with 350 * the other two). May not be null. 351 * @return A BasicNurbsCurve for a full-circle 352 * @throws DimensionException if one of the 3 points does not have at least 2 physical 353 * dimensions. 354 * @throws IllegalArgumentException if the input points are co-linear. 355 */ 356 public static BasicNurbsCurve createCircle(GeomPoint p1, GeomPoint p2, GeomPoint p3) throws DimensionException { 357 requireNonNull(p1); 358 requireNonNull(p2); 359 requireNonNull(p3); 360 361 // Determine the circle parameters from 3 points. 362 CircleInfo circle = GeomUtil.threePointCircle(p1, p2, p3); 363 364 // Now just create a circle from the center, radius, and direction vectors. 365 return createEllipseCrv(circle.center, circle.radius, circle.radius, circle.xhat, circle.yhat); 366 } 367 368 /** 369 * Create a full-circle {@link NurbsCurve} that is approximately centered on the 370 * specified origin point and passes through the two input points. The origin point is 371 * used to determine the plane that the circle lies in and is used to approximate the 372 * center of the circle. The true origin/center of the circle is calculated to ensure 373 * that the circle passes through the supplied edge points. The curve is parameterized 374 * such that 0 is at the first point and the parameterization increases toward the 2nd 375 * point. The NurbsCurve returned has a control polygon which has 9 control points and 376 * the shape of square. 377 * 378 * @param o Approximate origin or center to create the full-circle around. May not be 379 * null. 380 * @param p1 The 1st of the points that define the circle. May not be null. 381 * @param p2 The 2nd of the points that define the circle. May not be null. 382 * @return A BasicNurbsCurve for a full-circle 383 * @throws DimensionException if the origin or 2 points do not have at least 2 384 * physical dimensions. 385 */ 386 public static BasicNurbsCurve createCircleO(GeomPoint o, GeomPoint p1, GeomPoint p2) throws DimensionException { 387 requireNonNull(o); 388 requireNonNull(p1); 389 requireNonNull(p2); 390 391 // Determine the circle parameters from a center and 2 points. 392 CircleInfo circle = GeomUtil.centerTwoPointCircle(o, p1, p2); 393 394 // Now just create a circle from the center, radius, and direction vectors. 395 return createEllipseCrv(circle.center, circle.radius, circle.radius, circle.xhat, circle.yhat); 396 } 397 398 /** 399 * Create a full-ellipse {@link NurbsCurve} around the given origin/center with 400 * semi-major axis length <code>a</code> and semi-minor axis length <code>b</code>. 401 * The NurbsCurve returned has a control polygon which has 9 control points and the 402 * shape of a rectangle. 403 * 404 * @param o Origin or center to create the full-ellipse around. May not be null. 405 * @param a Semi-major axis length. May not be null. 406 * @param b Semi-minor axis length. May not be null. 407 * @param n A unit plane normal vector for the plane containing the circle. If 408 * <code>null</code> is passed, a default normal vector that points along the 409 * Z axis is used. 410 * @return A BasicNurbsCurve for a full-ellipse 411 * @throws DimensionException if the origin point or normal do not have at least 2 412 * physical dimensions. 413 */ 414 public static BasicNurbsCurve createEllipse(GeomPoint o, Parameter<Length> a, Parameter<Length> b, 415 GeomVector<Dimensionless> n) throws DimensionException { 416 requireNonNull(a); 417 requireNonNull(b); 418 419 int numDims = o.getPhyDimension(); 420 if (nonNull(n) && n.getPhyDimension() > numDims) 421 numDims = n.getPhyDimension(); 422 if (numDims < 2) 423 throw new DimensionException( 424 MessageFormat.format(RESOURCES.getString("dimensionNotAtLeast2"), 425 "input origin or normal", numDims)); 426 427 StackContext.enter(); 428 try { 429 if (nonNull(n)) { 430 o = o.toDimension(numDims); 431 n = n.toDimension(numDims); 432 } 433 434 // Make sure the normal vector exists and is a unit vector. 435 if (numDims == 2) 436 n = null; 437 else if (isNull(n)) 438 n = GeomUtil.makeZVector(numDims); 439 else 440 n = n.toUnitVector(); 441 442 // Create the yhat and xhat vectors (orthogonal Y and X in the plane of the circle). 443 Vector<Dimensionless> yhat = GeomUtil.calcYHat(n); 444 Vector<Dimensionless> xhat = GeomUtil.calcXHat(n, yhat); 445 446 // Now just create a circle from the center, radius, and direction vectors. 447 BasicNurbsCurve crv = createEllipseCrv(o, a, b, xhat, yhat); 448 return StackContext.outerCopy(crv); 449 450 } finally { 451 StackContext.exit(); 452 } 453 } 454 455 /** 456 * Create a full-ellipse {@link NurbsCurve} around the given origin/center with 457 * semi-major axis length <code>a</code> and semi-minor axis length <code>b</code>. 458 * The NurbsCurve returned has a control polygon which has 9 control points and the 459 * shape of a rectangle. 460 * 461 * @param o Origin or center to create the full-circle around. May not be null. 462 * @param a Semi-major axis length. May not be null. 463 * @param b Semi-minor axis length. May not be null. 464 * @param xhat A unit vector indicating the "X" direction in the plane of the ellipse. 465 * This must be orthogonal to yhat. May not be null. 466 * @param yhat A unit vector indicating the "Y" direction in the plane of the ellipse. 467 * This must be orthogonal to xhat. May not be null. 468 * @return A BasicNurbsCurve for a full-ellipse 469 */ 470 private static BasicNurbsCurve createEllipseCrv(GeomPoint o, Parameter<Length> a, Parameter<Length> b, 471 GeomVector<Dimensionless> xhat, GeomVector<Dimensionless> yhat) 472 throws DimensionException { 473 StackContext.enter(); 474 try { 475 double w = SQRT2 / 2.0; 476 double aValue = a.getValue(SI.METER); 477 double bValue = b.getValue(SI.METER); 478 479 // Create a list of control points. 480 FastTable<ControlPoint> cpList = FastTable.newInstance(); 481 cpList.setSize(9); 482 483 // Set 0 (right-start), 4 (left), and 8 (right-end) control points. 484 // p = a,0 485 Point p = Point.valueOf(xhat); 486 p = p.times(aValue); 487 cpList.set(0, ControlPoint.valueOf(o.plus(p), 1)); 488 cpList.set(4, ControlPoint.valueOf(o.plus(p.opposite()), 1)); // p = -a,0 489 cpList.set(8, ControlPoint.valueOf(o.plus(p), 1)); 490 491 // Set 1 (upper right) and 5 (lower left) control points. 492 // p = a,b 493 p = Point.valueOf(xhat.times(aValue).plus(yhat.times(bValue))); 494 cpList.set(1, ControlPoint.valueOf(o.plus(p), w)); 495 cpList.set(5, ControlPoint.valueOf(o.plus(p.opposite()), w)); // p = -a,-b 496 497 // Set 2 (top), and 6 (bottom)control points. 498 // p = 0,b 499 p = Point.valueOf(yhat); 500 p = p.times(bValue); 501 cpList.set(2, ControlPoint.valueOf(o.plus(p), 1)); 502 cpList.set(6, ControlPoint.valueOf(o.plus(p.opposite()), 1)); // p = 0,-b 503 504 // Set 3 (upper left) and 7 (lower right) control points. 505 // p = a,-b 506 p = Point.valueOf(xhat.times(aValue).plus(yhat.times(-bValue))); 507 cpList.set(7, ControlPoint.valueOf(o.plus(p), w)); 508 cpList.set(3, ControlPoint.valueOf(o.plus(p.opposite()), w)); // p = -a,b 509 510 // Create the NURBS curve. 511 KnotVector kv = KnotVector.newInstance(2, FC_KNOTS); 512 BasicNurbsCurve crv = BasicNurbsCurve.newInstance(cpList, kv); 513 return StackContext.outerCopy(crv); 514 515 } finally { 516 StackContext.exit(); 517 } 518 } 519 520 /** 521 * Create a circular arc {@link NurbsCurve} about the specified origin, with the 522 * specified radius and angular start and stop values. 523 * 524 * @param o The origin to create the arc around. May not be null. 525 * @param r The radius of the circular arc. May not be null. 526 * @param n A unit plane normal vector for the plane containing the circle. If 527 * <code>null</code> is passed, a default normal vector that points 528 * along the Z axis is used. 529 * @param thetaS The start angle of the arc (relative to the X axis). 530 * @param thetaE The end angle of the arc (relative to the X axis). If the end angle 531 * is smaller than the start angle, it is increased by increments of 532 * 2*PI until it is larger than the start angle. 533 * @return A BasicNurbsCurve for a circular arc 534 * @throws DimensionException if the origin point or normal do not have at least 2 535 * physical dimensions. 536 */ 537 public static BasicNurbsCurve createCircularArc(GeomPoint o, Parameter<Length> r, GeomVector<Dimensionless> n, 538 Parameter<Angle> thetaS, Parameter<Angle> thetaE) throws DimensionException { 539 requireNonNull(r); 540 int numDims = o.getPhyDimension(); 541 if (nonNull(n) && n.getPhyDimension() > numDims) 542 numDims = n.getPhyDimension(); 543 if (numDims < 2) 544 throw new DimensionException( 545 MessageFormat.format(RESOURCES.getString("dimensionNotAtLeast2"), 546 "input origin or normal", numDims)); 547 548 StackContext.enter(); 549 try { 550 if (nonNull(n)) { 551 o = o.toDimension(numDims); 552 n = n.toDimension(numDims); 553 } 554 555 // Make sure the normal vector exists and is a unit vector. 556 if (numDims == 2) 557 n = null; 558 else if (isNull(n)) 559 n = GeomUtil.makeZVector(numDims); 560 else 561 n = n.toUnitVector(); 562 563 // Create the yhat and xhat vectors (orthogonal Y and X in the plane of the circle). 564 Vector<Dimensionless> yhat = GeomUtil.calcYHat(n); 565 Vector<Dimensionless> xhat = GeomUtil.calcXHat(n, yhat); 566 567 BasicNurbsCurve crv = createEllipticalArcCrv(o, r, r, xhat, yhat, 568 thetaS.getValue(SI.RADIAN), thetaE.getValue(SI.RADIAN)); 569 return StackContext.outerCopy(crv); 570 571 } finally { 572 StackContext.exit(); 573 } 574 } 575 576 /** 577 * Create a circular arc {@link NurbsCurve} that passes through the input (not 578 * co-linear) points. The curve will be parameterized with 0 at the starting point and 579 * 1 at the last point. 580 * 581 * @param p1 The 1st of the points that define the arc (must not be colinear with the 582 * other two). May not be null. 583 * @param p2 The 2nd of the points that define the arc (must not be colinear with the 584 * other two). May not be null. 585 * @param p3 The 3rd of the points that define the arc (must not be colinear with the 586 * other two). May not be null. 587 * @return A BasicNurbsCurve for a circular arc through the input points 588 * @throws DimensionException if one of the 3 points does not have at least 2 physical 589 * dimensions. 590 * @throws IllegalArgumentException if the input points are co-linear. 591 */ 592 public static BasicNurbsCurve createCircularArc(GeomPoint p1, GeomPoint p2, GeomPoint p3) throws DimensionException { 593 requireNonNull(p1); 594 requireNonNull(p2); 595 requireNonNull(p3); 596 597 // Determine the circle parameters from 3 points. 598 CircleInfo circle = GeomUtil.threePointCircle(p1, p2, p3); 599 600 // Now just create a circle from the center, radius, and direction vectors. 601 return createEllipticalArcCrv(circle.center, circle.radius, circle.radius, circle.xhat, circle.yhat, 602 0, circle.angle.getValue(SI.RADIAN)); 603 } 604 605 /** 606 * Create a circular arc {@link NurbsCurve} that passes through the input points and 607 * that has the input tangents at the 1st point. The curve will be parameterized with 608 * 0 at the starting point and 1 at the last point. 609 * 610 * @param p1 The 1st of the points that define the arc. May not be null. 611 * @param t1 The tangent vector at p1 (must not point at p2). May not be null. 612 * @param p2 The 2nd of the points that define the arc. May not be null. 613 * @return A BasicNurbsCurve for a circular arc through the input points with the 614 * input tangent vector at the start. 615 * @throws DimensionException if one of the 3 points does not have at least 2 physical 616 * dimensions. 617 */ 618 public static BasicNurbsCurve createCircularArc(GeomPoint p1, GeomVector t1, GeomPoint p2) 619 throws DimensionException { 620 requireNonNull(p1); 621 requireNonNull(t1); 622 requireNonNull(p2); 623 624 // Determine the circle parameters from 3 points. 625 CircleInfo circle = GeomUtil.twoPointTangentCircle(p1, t1, p2); 626 627 // Now just create a circle from the center, radius, and direction vectors. 628 return createEllipticalArcCrv(circle.center, circle.radius, circle.radius, circle.xhat, circle.yhat, 629 0, circle.angle.getValue(SI.RADIAN)); 630 } 631 632 /** 633 * Create a circular arc {@link NurbsCurve} that is approximately centered on the 634 * specified origin point and passes through the two input points. The origin point is 635 * used to determine the plane that the circle lies in and is used to approximate the 636 * center of the circle. The true origin/center of the circle is calculated to ensure 637 * that the circle passes through the supplied edge points. The curve is parameterized 638 * such that 0 is at the first point and 1 is at the last point. The NurbsCurve 639 * returned has a control polygon which has 9 control points and the shape of square. 640 * 641 * @param o Approximate origin or center to create the arc around. May not be null. 642 * @param p1 The 1st of the points that define the arc. May not be null. 643 * @param p2 The 2nd of the points that define the arc. May not be null. 644 * @return A BasicNurbsCurve for a circular arc through the input points. 645 * @throws DimensionException if the origin or 2 points do not have at least 2 646 * physical dimensions. 647 */ 648 public static BasicNurbsCurve createCircularArcO(GeomPoint o, GeomPoint p1, GeomPoint p2) throws DimensionException { 649 requireNonNull(o); 650 requireNonNull(p1); 651 requireNonNull(p2); 652 653 // Determine the circle parameters from a center and 2 points. 654 CircleInfo circle = GeomUtil.centerTwoPointCircle(o, p1, p2); 655 656 // Now just create a circle from the center, radius, and direction vectors. 657 return createEllipticalArcCrv(circle.center, circle.radius, circle.radius, circle.xhat, circle.yhat, 658 0, circle.angle.getValue(SI.RADIAN)); 659 } 660 661 /** 662 * Create an elliptical arc {@link NurbsCurve} about the specified origin, with the 663 * specified semi-major axis length, semi-minor axis length and angular start and stop 664 * values. 665 * 666 * @param o The origin to create the arc around. May not be null. 667 * @param a The ellipse semi-major axis length (half the ellipse diameter in the 668 * xhat direction). May not be null. 669 * @param b The ellipse semi-minor axis length (half the ellipse diameter in the 670 * yhat direction). May not be null. 671 * @param xhat A unit vector indicating the "X" or semi-major axis direction in the 672 * plane of the ellipse. This must be orthogonal to yhat. May not be 673 * null. 674 * @param yhat A unit vector indicating the "Y" or semi-minor direction in the plane 675 * of the ellipse. This must be orthogonal to xhat. May not be null. 676 * @param thetaS The start angle of the arc (relative to the xhat axis). May not be 677 * null. 678 * @param thetaE The end angle of the arc (relative to the xhat axis). If the end 679 * angle is smaller than the start angle, it is increased by increments 680 * of 2*PI until it is larger than the start angle. May not be null. 681 * @return A BasicNurbsCurve for an elliptical arc 682 * @throws DimensionException if the origin point or axes do not have at least 2 683 * physical dimensions. 684 */ 685 public static BasicNurbsCurve createEllipticalArc( 686 GeomPoint o, Parameter<Length> a, Parameter<Length> b, 687 GeomVector<Dimensionless> xhat, GeomVector<Dimensionless> yhat, 688 Parameter<Angle> thetaS, Parameter<Angle> thetaE) throws DimensionException { 689 int numDims = GeomUtil.maxPhyDimension(o, xhat, yhat); 690 if (numDims < 2) 691 throw new DimensionException( 692 MessageFormat.format(RESOURCES.getString("dimensionNotAtLeast2"), 693 "input origin and axis vector", numDims)); 694 requireNonNull(a); 695 requireNonNull(b); 696 requireNonNull(thetaS); 697 requireNonNull(thetaE); 698 699 StackContext.enter(); 700 try { 701 o = o.toDimension(numDims); 702 xhat = xhat.toDimension(numDims); 703 yhat = yhat.toDimension(numDims); 704 705 // Make sure the axis vectors are unit vectors with no origin offset. 706 xhat = xhat.toUnitVector(); 707 xhat.setOrigin(Point.newInstance(numDims)); 708 yhat = yhat.toUnitVector(); 709 yhat.setOrigin(Point.newInstance(numDims)); 710 711 BasicNurbsCurve crv = createEllipticalArcCrv(o, a, b, xhat, yhat, thetaS.getValue(SI.RADIAN), thetaE.getValue(SI.RADIAN)); 712 return StackContext.outerCopy(crv); 713 714 } finally { 715 StackContext.exit(); 716 } 717 } 718 719 /** 720 * Create a canonical elliptical arc {@link NurbsCurve} about the specified origin, 721 * with the specified semi-major and semi-minor axes and angular start and stop 722 * values. 723 * 724 * @param o The origin to create the arc around. May not be null. 725 * @param a The ellipse semi-major (or minor if b is major) axis. May not be null. 726 * @param b The ellipse semi-minor (or major if a is minor) axis. May not be null. 727 * @param xhat A unit vector indicating the "X" direction in the plane of the ellipse. 728 * This must be orthogonal to yhat. May not be null. 729 * @param yhat A unit vector indicating the "Y" direction in the plane of the ellipse. 730 * This must be orthogonal to xhat. May not be null. 731 * @param ths The start angle of the arc (relative to the X axis) in radians. 732 * @param the The end angle of the arc (relative to the X axis) in radians. If the 733 * end angle is smaller than the start angle, it is increased by 734 * increments of 2*PI until it is larger than the start angle. 735 * @return A BasicNurbsCurve for an elliptical arc 736 */ 737 private static BasicNurbsCurve createEllipticalArcCrv(GeomPoint o, Parameter<Length> a, Parameter<Length> b, 738 GeomVector<Dimensionless> xhat, GeomVector<Dimensionless> yhat, 739 double ths, double the) throws DimensionException { 740 741 // Reference: Piegl, L., Tiller, W., The Nurbs Book, 2nd Edition, Springer-Verlag, Berlin, 1997. 742 // Algorithm: 7.1, pg. 308. 743 StackContext.enter(); 744 try { 745 // Get ths and the in the range 0 - 2*PI. 746 ths = GeomUtil.bound2Pi(ths); 747 the = GeomUtil.bound2Pi(the); 748 749 // Make sure the end angle is > than the start angle. 750 if (the <= ths) 751 the += TWO_PI; 752 double theta = the - ths; 753 754 // How many arc segments are required? 755 int narcs = 4; 756 if (theta <= HALF_PI) 757 narcs = 1; 758 else if (theta <= PI) 759 narcs = 2; 760 else if (theta <= PI + HALF_PI) 761 narcs = 3; 762 763 double dtheta = theta / narcs; 764 int num = 2 * narcs; // num+1 control points. 765 double w1 = cos(dtheta / 2); // dtheta/2 is the base angle 766 double aValue = a.getValue(SI.METER); 767 double bValue = b.getValue(SI.METER); 768 Unit<Length> unit = o.getUnit(); 769 770 // Create a list of control points. 771 FastTable<ControlPoint> cpList = FastTable.newInstance(); 772 cpList.setSize(num + 1); 773 774 // Create the starting control point: p0 = a*cos(ths), b*sin(ths) 775 double sinths = sin(ths), cosths = cos(ths); 776 MutablePoint p0 = MutablePoint.valueOf(xhat.times(aValue * cosths).plus(yhat.times(bValue * sinths))).to(unit); 777 cpList.set(0, ControlPoint.valueOf(o.plus(p0), 1)); 778 779 // Create starting perpendicular vector: t0 = -a*sin(ths), b*cos(ths) 780 Vector<Dimensionless> t0 = xhat.times(-aValue * sinths).plus(yhat.times(bValue * cosths)).toUnitVector(); 781 782 int index = 0; 783 double angle = ths; 784 MutablePoint p1 = MutablePoint.newInstance(o.getPhyDimension(), unit); 785 MutablePoint p2 = MutablePoint.newInstance(o.getPhyDimension(), unit); 786 for (int i = 1; i <= narcs; i++) { 787 angle += dtheta; 788 double sina = sin(angle); 789 double cosa = cos(angle); 790 791 // Create segment end control point: p2 = a*cos(angle), b*sin(angle) 792 p2.set(xhat.times(aValue * cosa).plus(yhat.times(bValue * sina))); 793 cpList.set(index + 2, ControlPoint.valueOf(o.plus(p2), 1)); 794 795 // Create end perpendicular vector: t2 = -a*sin(angle), b*cos(angle) 796 Vector<Dimensionless> t2 = xhat.times(-aValue * sina).plus(yhat.times(bValue * cosa)).toUnitVector(); 797 798 // Find the intersection of the two perpendiculars. 799 GeomUtil.lineLineIntersect(p0, t0, p2, t2, GTOL, null, p1, null); 800 801 // Create the intermediate control point. 802 cpList.set(index + 1, ControlPoint.valueOf(o.plus(p1), w1)); 803 804 index += 2; 805 if (i < narcs) { 806 p0.set(p2); 807 t0 = t2; 808 } 809 } 810 811 // Create the knot vector. 812 FastTable<Float64> uKnot = FastTable.newInstance(); 813 int j = num + 1; 814 uKnot.setSize(j + 3); 815 for (int i = 0; i < 3; i++) { 816 uKnot.set(i, Float64.ZERO); 817 uKnot.set(i + j, Float64.ONE); 818 } 819 switch (narcs) { 820 case 1: 821 break; 822 case 2: 823 uKnot.set(3, Float64.valueOf(0.5)); 824 uKnot.set(4, Float64.valueOf(0.5)); 825 break; 826 case 3: 827 uKnot.set(3, Float64.valueOf(1. / 3.)); 828 uKnot.set(4, Float64.valueOf(1. / 3.)); 829 uKnot.set(5, Float64.valueOf(2. / 3.)); 830 uKnot.set(6, Float64.valueOf(2. / 3.)); 831 break; 832 case 4: 833 uKnot.set(3, Float64.valueOf(0.25)); 834 uKnot.set(4, Float64.valueOf(0.25)); 835 uKnot.set(5, Float64.valueOf(0.5)); 836 uKnot.set(6, Float64.valueOf(0.5)); 837 uKnot.set(7, Float64.valueOf(0.75)); 838 uKnot.set(8, Float64.valueOf(0.75)); 839 break; 840 } 841 842 // Create the NURBS curve. 843 KnotVector kv = KnotVector.newInstance(2, Float64Vector.valueOf(uKnot)); 844 BasicNurbsCurve crv = BasicNurbsCurve.newInstance(cpList, kv); 845 return StackContext.outerCopy(crv); 846 847 } finally { 848 StackContext.exit(); 849 } 850 } 851 852 /** 853 * Create a planar parabolic arc {@link NurbsCurve} that joins two end points with the 854 * specified tangent vectors at each end point. If the input tangent vectors are 855 * either parallel or not coplanar, then a straight line between the two end points is 856 * returned. If the intersection of the two input tangent vectors does not fall 857 * between the two input end points, then a straight line between the two end points 858 * is returned. 859 * 860 * @param p0 The starting end of the arc. May not be null. 861 * @param t0 The tangent vector at the starting end of the arc. May not be null. 862 * @param p2 The stopping end of the arc. May not be null. 863 * @param t2 The tangent vector at the stopping end of the arc. May not be null. 864 * @return A BasicNurbsCurve for a parabolic arc between the input points. 865 * @throws DimensionException if the input points and tangents do not have at least 2 866 * physical dimensions. 867 * @see #createBlend 868 */ 869 public static BasicNurbsCurve createParabolicArc(GeomPoint p0, GeomVector<Dimensionless> t0, 870 GeomPoint p2, GeomVector<Dimensionless> t2) throws DimensionException { 871 requireNonNull(p0); 872 requireNonNull(t0); 873 requireNonNull(p2); 874 requireNonNull(t2); 875 876 // Get everything into a consistant set of dimensions. 877 int numDims = GeomUtil.maxPhyDimension(p0, t0, p2, t2); 878 if (numDims < 2) 879 throw new DimensionException( 880 MessageFormat.format(RESOURCES.getString("dimensionNotAtLeast2"), 881 "input points and tangent vectors", numDims)); 882 883 StackContext.enter(); 884 try { 885 Unit<Length> unit = p0.getUnit(); 886 Point nullPnt = Point.newInstance(numDims); 887 p0 = p0.toDimension(numDims); 888 t0 = t0.toDimension(numDims).toUnitVector(); 889 t0.setOrigin(nullPnt); 890 p2 = p2.toDimension(numDims).to(unit); 891 t2 = t2.toDimension(numDims).toUnitVector(); 892 t2.setOrigin(nullPnt); 893 894 // Find the intersection of the two tangent vector lines. 895 MutablePoint p1 = MutablePoint.newInstance(numDims, unit); 896 MutablePoint s1 = MutablePoint.newInstance(2, unit); 897 IntersectType type = GeomUtil.lineLineIntersect(p0, t0, p2, t2, GTOL, s1, p1, null); 898 if (type != IntersectType.INTERSECT) { 899 // The input tangent vectors are not co-planar or they are parellel. 900 return StackContext.outerCopy(CurveFactory.createLine(2, p0, p2)); 901 } 902 903 // Make sure the interesection falls between the two input points. 904 // Parametric positions returned from lineLineIntersect() are scaled such that 905 // 0.0 is the line origin point and any other value represents the signed distance 906 // that the intersection point is away from the origin along the direction vector. 907 double s0 = s1.getValue(0); // Distance intersection is away from 1st line origin. 908 double s2 = s1.getValue(1); // Distance intersection is away from 2nd line origin. 909 if (s0 < 0 || s2 > 0) { 910 // The intersection does not fall between the two input points. 911 return StackContext.outerCopy(CurveFactory.createLine(2, p0, p2)); 912 } 913 914 // Create the control points. 915 FastTable<ControlPoint> cps = FastTable.newInstance(); 916 cps.add(ControlPoint.valueOf(p0.immutable(), 1)); 917 cps.add(ControlPoint.valueOf(p1.immutable(), 1)); 918 cps.add(ControlPoint.valueOf(p2.immutable(), 1)); 919 920 // Create the knot vector. 921 KnotVector kv = KnotVector.newInstance(2, 0., 0., 0., 1., 1., 1.); 922 923 // Create the curve. 924 BasicNurbsCurve crv = BasicNurbsCurve.newInstance(cps, kv); 925 return StackContext.outerCopy(crv); 926 927 } finally { 928 StackContext.exit(); 929 } 930 } 931 932 /** 933 * Create the first quadrant of a common exponent super-ellipse {@link NurbsCurve} 934 * about the specified origin, with the specified semi-major and semi-minor axes and 935 * exponent. Other quadrants can be represented by reflection. A common exponent 936 * super-elliptic arc can be described using the following equation:<br> 937 * <code>(x/a)^n + (y/b)^n = 1</code><br> 938 * where a is the semi-major and b is the semi-minor axis, and n is the exponent of 939 * the super-ellipse (n > 0). Special cases include a circle (with a = b, and n = 2), 940 * an ellipse (with a ≠b, and n = 2) and a straight line (n = 1). If n is large, a 941 * box is approximated. This type of super-ellipse is represented exactly by a NURBS 942 * curve. 943 * 944 * @param o The origin to create the arc around. May not be null. 945 * @param a The super-ellipse semi-major (or minor if b is major) axis. May not be 946 * null. 947 * @param b The super-ellipse semi-minor (or major if a is minor) axis. May not be 948 * null. 949 * @param n The exponent for the super-ellipse (must be > 0). 950 * @param nhat A unit plane normal vector for the plane containing the super-ellipse. 951 * If <code>null</code> is passed, a default normal vector that points 952 * along the Z axis is used. 953 * @return A BasicNurbsCurve for the 1st quadrant of a super-ellipse. 954 */ 955 public static BasicNurbsCurve createSuperEllipticalArc( 956 GeomPoint o, Parameter<Length> a, Parameter<Length> b, double n, 957 GeomVector<Dimensionless> nhat) throws DimensionException { 958 requireNonNull(a); 959 requireNonNull(b); 960 961 int numDims = o.getPhyDimension(); 962 if (nonNull(nhat) && nhat.getPhyDimension() > numDims) 963 numDims = nhat.getPhyDimension(); 964 if (numDims < 2) 965 throw new DimensionException( 966 MessageFormat.format(RESOURCES.getString("dimensionNotAtLeast2"), 967 "input origin or normal", numDims)); 968 969 StackContext.enter(); 970 try { 971 if (nonNull(nhat)) { 972 o = o.toDimension(numDims); 973 nhat = nhat.toDimension(numDims); 974 } 975 976 // Make sure the normal vector exists and is a unit vector. 977 if (numDims == 2) 978 nhat = null; 979 else if (isNull(nhat)) 980 nhat = GeomUtil.makeZVector(numDims); 981 else 982 nhat = nhat.toUnitVector(); 983 984 // Create the yhat and xhat vectors (orthogonal Y and X in the plane of the circle). 985 Vector<Dimensionless> yhat = GeomUtil.calcYHat(nhat); 986 Vector<Dimensionless> xhat = GeomUtil.calcXHat(nhat, yhat); 987 988 BasicNurbsCurve crv = createSuperEllipticalArc(o, a, b, n, xhat, yhat); 989 return StackContext.outerCopy(crv); 990 991 } finally { 992 StackContext.exit(); 993 } 994 } 995 996 /** 997 * Create the first quadrant of a common exponent super-ellipse {@link NurbsCurve} 998 * about the specified origin, with the specified semi-major and semi-minor axes and 999 * exponent. Other quadrants can be represented by reflection. A super-elliptic arc 1000 * can be described using the following equation:<br> 1001 * <code>(x/a)^n + (y/b)^n = 1</code><br> 1002 * where a is the semi-major and b is the semi-minor axis, and n is the exponent of 1003 * the super-ellipse (n > 0). Special cases include a circle (with a = b, and n = 2), 1004 * an ellipse (with a ≠b, and n = 2) and a straight line (n = 1). If n is large, a 1005 * box is approximated. This type of super-ellipse is represented exactly by a NURBS 1006 * curve. 1007 * 1008 * @param o The origin to create the arc around. May not be null. 1009 * @param a The super-ellipse semi-major (or minor if b is major) axis. May not be 1010 * null. 1011 * @param b The super-ellipse semi-minor (or major if a is minor) axis. May not be 1012 * null. 1013 * @param n The exponent for the super-ellipse (must be > 0). 1014 * @param xhat A unit vector indicating the "X" or semi-major axis direction in the 1015 * plane of the super-ellipse. This must be orthogonal to yhat. May not be 1016 * null. 1017 * @param yhat A unit vector indicating the "Y" or semi-minor axis direction in the 1018 * plane of the super-ellipse. This must be orthogonal to xhat. May not be 1019 * null. 1020 * @return A BasicNurbsCurve for the 1st quadrant of a super-ellipse. 1021 */ 1022 public static BasicNurbsCurve createSuperEllipticalArc( 1023 GeomPoint o, Parameter<Length> a, Parameter<Length> b, double n, 1024 GeomVector<Dimensionless> xhat, GeomVector<Dimensionless> yhat) throws DimensionException { 1025 requireNonNull(o); 1026 requireNonNull(a); 1027 requireNonNull(b); 1028 requireNonNull(xhat); 1029 requireNonNull(yhat); 1030 1031 // Reference: Handbook of Grid Generation, Chapter 30, Section 30.3.5 1032 // and Joe Huwaldt's hand written derivations. 1033 if (n <= MathTools.EPS) { 1034 throw new IllegalArgumentException(MessageFormat.format(RESOURCES.getString("seExponent"), "n")); 1035 } 1036 1037 int numDims = GeomUtil.maxPhyDimension(o, xhat, yhat); 1038 if (numDims < 2) 1039 throw new DimensionException( 1040 MessageFormat.format(RESOURCES.getString("dimensionNotAtLeast2"), 1041 "input origin and axis vector", numDims)); 1042 1043 StackContext.enter(); 1044 try { 1045 o = o.toDimension(numDims); 1046 xhat = xhat.toDimension(numDims); 1047 yhat = yhat.toDimension(numDims); 1048 1049 // Make sure the axis vectors are unit vectors with no origin offset. 1050 xhat = xhat.toUnitVector(); 1051 yhat = yhat.toUnitVector(); 1052 1053 // Compute the weight for the middle control point. 1054 double aValue = a.getValue(SI.METER); 1055 double bValue = b.getValue(SI.METER); 1056 double rhoD = sqrt(aValue * aValue + bValue * bValue); 1057 double rhoM = 0.5 * rhoD; 1058 double rhoh = rhoD / pow(2, 1.0 / n); 1059 double w1 = (rhoh - rhoM) / (rhoD - rhoh); 1060 1061 // Create a list of control points. 1062 FastTable<ControlPoint> cpList = FastTable.newInstance(); 1063 1064 // Create the starting control point: p0 = a, 0 1065 Point p0 = Point.valueOf(xhat.times(aValue)); 1066 cpList.add(ControlPoint.valueOf(o.plus(p0), 1)); 1067 1068 if (w1 >= 0) { 1069 // Create the middle control point: p1 = a, b 1070 Point p1 = Point.valueOf(xhat.times(aValue).plus(yhat.times(bValue))); 1071 cpList.add(ControlPoint.valueOf(o.plus(p1), w1)); 1072 1073 } else { 1074 // Use the origin as the control point: p1 = o 1075 // And change the rhoh definition to use inverted exponents. 1076 rhoh = rhoD / pow(2, n); 1077 w1 = (rhoh - rhoM) / (rhoD - rhoh); 1078 cpList.add(ControlPoint.valueOf(o.immutable(), w1)); 1079 } 1080 1081 // Create ending control point: p2 = 0, b 1082 Point p2 = Point.valueOf(yhat.times(bValue)); 1083 cpList.add(ControlPoint.valueOf(o.plus(p2), 1)); 1084 1085 // Create the knot vector. 1086 KnotVector kv = KnotVector.newInstance(2, 0., 0., 0., 1., 1., 1.); 1087 1088 // Create the NURBS curve. 1089 BasicNurbsCurve crv = BasicNurbsCurve.newInstance(cpList, kv).to(o.getUnit()); 1090 return StackContext.outerCopy(crv); 1091 1092 } finally { 1093 StackContext.exit(); 1094 } 1095 } 1096 1097 /** 1098 * Create the first quadrant of a general super-ellipse {@link NurbsCurve} about the 1099 * specified origin, with the specified semi-major and semi-minor axes and exponents. 1100 * Other quadrants can be represented by reflection. A general super-elliptic arc can 1101 * be described using the following equation:<br> 1102 * <code>(x/a)^m + (y/b)^n = 1</code><br> 1103 * where a is the semi-major and b is the semi-minor axis, and m & n are the exponents 1104 * of the super-ellipse (m & n > 0). Special cases include a circle (with a = b, and m 1105 * = n = 2), an ellipse (with a ≠b, and m = n = 2) and a straight line (m = n = 1). 1106 * If m & n are large, a box is approximated. This type of super-ellipse is generally 1107 * only approximated by a NURBS curve to the specified tolerance. If the input 1108 * exponents are equal, then an exact representation is returned. 1109 * 1110 * @param o The origin to create the arc around. May not be null. 1111 * @param a The super-ellipse semi-major (or minor if b is major) axis. May not be 1112 * null. 1113 * @param b The super-ellipse semi-minor (or major if a is minor) axis. May not be 1114 * null. 1115 * @param m The exponent on the x/a term of the super-ellipse equation (must be > 1116 * 0). 1117 * @param n The exponent for the y/b term of the super-ellipse equation (must be > 1118 * 0). 1119 * @param xhat A unit vector indicating the "X" or semi-major axis direction in the 1120 * plane of the super-ellipse. This must be orthogonal to yhat. May not be 1121 * null. 1122 * @param yhat A unit vector indicating the "Y" or semi-minor axis direction in the 1123 * plane of the super-ellipse. This must be orthogonal to xhat. May not be 1124 * null. 1125 * @param tol A geometric tolerance on how accurately the NURBS curve must represent 1126 * the true super-ellipse. If m == n, then an exact representation is 1127 * returned and "tol" is ignored, otherwise tol may not be null. 1128 * @return A BasicNurbsCurve for the 1st quadrant of a super-ellipse. 1129 */ 1130 public static BasicNurbsCurve createSuperEllipticalArc( 1131 GeomPoint o, Parameter<Length> a, Parameter<Length> b, double m, double n, 1132 GeomVector<Dimensionless> xhat, GeomVector<Dimensionless> yhat, Parameter<Length> tol) 1133 throws DimensionException, IllegalArgumentException, SingularMatrixException, ParameterException { 1134 requireNonNull(o); 1135 requireNonNull(a); 1136 requireNonNull(b); 1137 requireNonNull(xhat); 1138 requireNonNull(yhat); 1139 1140 // Reference: Handbook of Grid Generation, Chapter 30, Section 30.3.5 1141 // and Joe Huwaldt's hand written derivations. 1142 if (m <= MathTools.EPS) { 1143 throw new IllegalArgumentException(MessageFormat.format(RESOURCES.getString("seExponent"), "m")); 1144 } 1145 if (n <= MathTools.EPS) { 1146 throw new IllegalArgumentException(MessageFormat.format(RESOURCES.getString("seExponent"), "n")); 1147 } 1148 1149 if (MathTools.isApproxEqual(m, n)) { 1150 // A common exponent super-ellipse is more accurate. 1151 return createSuperEllipticalArc(o, a, b, n, xhat, yhat); 1152 } 1153 1154 int numDims = GeomUtil.maxPhyDimension(o, xhat, yhat); 1155 if (numDims < 2) 1156 throw new DimensionException( 1157 MessageFormat.format(RESOURCES.getString("dimensionNotAtLeast2"), 1158 "input origin and axis vector", numDims)); 1159 1160 requireNonNull(tol); 1161 StackContext.enter(); 1162 try { 1163 o = o.toDimension(numDims); 1164 xhat = xhat.toDimension(numDims); 1165 yhat = yhat.toDimension(numDims); 1166 1167 // Make sure the axis vectors are unit vectors with no origin offset. 1168 xhat = xhat.toUnitVector(); 1169 yhat = yhat.toUnitVector(); 1170 1171 // Compute the weight for the middle control point. 1172 double aValue = a.getValue(SI.METER); 1173 double bValue = b.getValue(SI.METER); 1174 1175 // Create an initial, high resolution curve. 1176 int size = 51; 1177 double theta = 0.; 1178 double dTheta = 90. / (size - 1); 1179 PointString pnts = PointString.newInstance(); 1180 GeomList tangents = GeomList.newInstance(); 1181 for (int i = 0; i < size; ++i) { 1182 // Get the parametric position on the super-ellipse. 1183 double thetaR = toRadians(theta); 1184 double ctheta = cos(thetaR); 1185 double stheta = sin(thetaR); 1186 1187 // Compute the coordinates of a point on the super-ellipse. 1188 double xh = aValue * pow(abs(ctheta), 2. / m); 1189 double yh = bValue * pow(abs(stheta), 2. / n); 1190 Point p = o.plus(Point.valueOf(xhat.times(xh).plus(yhat.times(yh)))); 1191 pnts.add(p); 1192 1193 // Compute the local slope of the super-ellipse. 1194 // dz/dy = -b/x*m/n*(x/a)^m*(1 - (x/a)^m))^(1/n-1) 1195 // Unfortunately, this equation breaks down at theta=0 (xh = aValue). 1196 double xh_am = pow(xh / aValue, m); 1197 double dz_dy = -bValue / xh * m / n * xh_am * pow(1 - xh_am, 1. / n - 1.); 1198 Vector t; 1199 if (Double.isInfinite(dz_dy)) { 1200 t = Vector.valueOf(yhat.times(-1 * Math.signum(dz_dy))).toUnitVector(); 1201 1202 } else if (abs(dz_dy) < 1e-15) { 1203 if (i == 0) { 1204 // Numerically approximate it. 1205 double dthetaR = 0.017453292519943 / 2; // 1/2 degree 1206 double xi = aValue * pow(abs(cos(dthetaR)), 2. / m); 1207 double yi = bValue * pow(abs(sin(dthetaR)), 2. / n); 1208 t = Vector.valueOf(xhat.times(xi - aValue).plus(yhat.times(yi))).toUnitVector(); 1209 1210 } else 1211 t = Vector.valueOf(xhat.times(-1)).toUnitVector(); 1212 1213 } else { 1214 t = Vector.valueOf(xhat.times(-1).plus(yhat.times(-dz_dy))).toUnitVector(); 1215 } 1216 //t.setOrigin(p); 1217 tangents.add(t); 1218 1219 theta = theta + dTheta; 1220 } 1221 1222 // Create the NURBS curve. 1223 // Fit a curve to the calculated points using the calculated slopes. 1224 BasicNurbsCurve crv = fitPoints(2, pnts, tangents, true, 1.0).to(o.getUnit()); 1225 1226 // Thin down the high resolution curve to the required tolerance. 1227 crv = CurveUtils.thinKnotsToTolerance(crv, tol); 1228 1229 return StackContext.outerCopy(crv); 1230 1231 } finally { 1232 StackContext.exit(); 1233 } 1234 } 1235 1236 /** 1237 * Create the first quadrant of a general super-ellipse {@link NurbsCurve} about the 1238 * specified origin, with the specified semi-major and semi-minor axes and exponents. 1239 * Other quadrants can be represented by reflection. A general super-elliptic arc can 1240 * be described using the following equation:<br> 1241 * <code>(x/a)^m + (y/b)^n = 1</code><br> 1242 * where a is the semi-major and b is the semi-minor axis, and m & n are the exponents 1243 * of the super-ellipse. Special cases include a circle (with a = b, and m = n = 2), 1244 * an ellipse (with a ≠b, and m = n = 2) and a straight line (m = n = 1). If m & n 1245 * are large, a box is approximated. This type of super-ellipse is generally only 1246 * approximated by a NURBS curve to the specified tolerance. If the input exponents 1247 * are equal, then an exact representation is returned. 1248 * 1249 * @param o The origin to create the arc around. May not be null. 1250 * @param a The super-ellipse semi-major (or minor if b is major) axis. May not be null. 1251 * @param b The super-ellipse semi-minor (or major if a is minor) axis. May not be null. 1252 * @param m The exponent on the x/a term of the super-ellipse equation (must be > 1253 * 0). 1254 * @param n The exponent for the y/b term of the super-ellipse equation (must be > 1255 * 0). 1256 * @param nhat A unit plane normal vector for the plane containing the super-ellipse. 1257 * If <code>null</code> is passed, a default normal vector that points 1258 * along the Z axis is used. 1259 * @param tol A geometric tolerance on how accurately the NURBS curve must represent 1260 * the true super-ellipse. If m == n, then an exact representation is 1261 * returned and "tol" is ignored, otherwise "tol" may not be null.. 1262 * @return A BasicNurbsCurve for the 1st quadrant of a super-ellipse. 1263 */ 1264 public static BasicNurbsCurve createSuperEllipticalArc( 1265 GeomPoint o, Parameter<Length> a, Parameter<Length> b, double m, double n, 1266 GeomVector<Dimensionless> nhat, Parameter<Length> tol) 1267 throws DimensionException, IllegalArgumentException, SingularMatrixException, ParameterException { 1268 1269 int numDims = o.getPhyDimension(); 1270 if (nonNull(nhat) && nhat.getPhyDimension() > numDims) 1271 numDims = nhat.getPhyDimension(); 1272 if (numDims < 2) 1273 throw new DimensionException( 1274 MessageFormat.format(RESOURCES.getString("dimensionNotAtLeast2"), 1275 "input origin or normal", numDims)); 1276 1277 StackContext.enter(); 1278 try { 1279 if (nonNull(nhat)) { 1280 o = o.toDimension(numDims); 1281 nhat = nhat.toDimension(numDims); 1282 } 1283 1284 // Make sure the normal vector exists and is a unit vector. 1285 if (numDims == 2) 1286 nhat = null; 1287 else if (isNull(nhat)) 1288 nhat = GeomUtil.makeZVector(numDims); 1289 else 1290 nhat = nhat.toUnitVector(); 1291 1292 // Create the yhat and xhat vectors (orthogonal Y and X in the plane of the circle). 1293 Vector<Dimensionless> yhat = GeomUtil.calcYHat(nhat); 1294 Vector<Dimensionless> xhat = GeomUtil.calcXHat(nhat, yhat); 1295 1296 BasicNurbsCurve crv = createSuperEllipticalArc(o, a, b, m, n, xhat, yhat, tol); 1297 return StackContext.outerCopy(crv); 1298 1299 } finally { 1300 StackContext.exit(); 1301 } 1302 } 1303 1304 /** 1305 * Create a potentially non-planar cubic Bezier {@link NurbsCurve} that joins or 1306 * blends two end points with the specified tangent vectors at each end point. If unit 1307 * vectors are supplied, then they will be scaled by the distance between the end 1308 * points and the multiplicative factor provided (tanlen). 1309 * 1310 * @param p1 The starting end of the arc. May not be null. 1311 * @param t1 The tangent vector at the starting end of the arc. May not be null. 1312 * @param p2 The stopping end of the arc. May not be null. 1313 * @param t2 The tangent vector at the stopping end of the arc. May not be null. 1314 * @param unitFlag A flag indicating that the tangents are all unit vectors if 1315 * <code>true</code>. 1316 * @param tanlen A multiplicative factor for the tangent (must be greater than 0 or 1317 * the function aborts). Used only if <code>unitFlat == true</code>. 1318 * Applied to the distance between end points to scale the tangent 1319 * vectors. 1320 * @return A BasicNurbsCurve that represents a Bezier curve between the input points 1321 * with the input tangent vectors. 1322 * @throws DimensionException if the input points and tangents do not have at least 2 1323 * physical dimensions. 1324 * @see #createParabolicArc 1325 */ 1326 public static BasicNurbsCurve createBlend(GeomPoint p1, GeomVector t1, 1327 GeomPoint p2, GeomVector t2, boolean unitFlag, double tanlen) throws DimensionException { 1328 requireNonNull(p1); 1329 requireNonNull(p2); 1330 requireNonNull(t1); 1331 requireNonNull(t2); 1332 1333 int degree = 3; 1334 if (unitFlag && tanlen <= 0) 1335 throw new IllegalArgumentException( 1336 MessageFormat.format(RESOURCES.getString("posMultFactorErr"), tanlen)); 1337 1338 int numDims = GeomUtil.maxPhyDimension(p1, t1, p2, t2); 1339 if (numDims < 2) 1340 throw new DimensionException( 1341 MessageFormat.format(RESOURCES.getString("dimensionNotAtLeast2"), 1342 "input points and tangent vectors", numDims)); 1343 1344 StackContext.enter(); 1345 try { 1346 // Determine the distance between the two points. 1347 Parameter<Length> dist = p1.distance(p2); 1348 if (dist.getValue() < EPSC) 1349 return StackContext.outerCopy(createPoint(degree, p1)); 1350 1351 // Get everything into a consistant set of dimensions. 1352 p1 = p1.toDimension(numDims); 1353 p2 = p2.toDimension(numDims); 1354 t1 = t1.toDimension(numDims); 1355 t2 = t2.toDimension(numDims); 1356 if (unitFlag) { 1357 dist = dist.times(tanlen); 1358 t1 = t1.toUnitVector().times(dist); 1359 t2 = t2.toUnitVector().times(dist); 1360 } 1361 1362 t1 = t1.opposite(); 1363 1364 // Find intermediate control points. 1365 Point p3 = p1.minus(Point.valueOf(t1)); 1366 Point p4 = p2.minus(Point.valueOf(t2)); 1367 1368 // Create the control points vector. 1369 int numCP = degree + 1; 1370 FastTable<ControlPoint> cps = FastTable.newInstance(); 1371 1372 cps.add(ControlPoint.valueOf(p1.immutable(), 1)); 1373 cps.add(ControlPoint.valueOf(p3, 1)); 1374 cps.add(ControlPoint.valueOf(p4, 1)); 1375 cps.add(ControlPoint.valueOf(p2.immutable(), 1)); 1376 1377 // Create a Bezier curve knot vector of the require degree. 1378 FastTable<Float64> U = FastTable.newInstance(); 1379 int numKnots = numCP + degree + 1; 1380 int numKnotso2 = numKnots / 2; 1381 for (int i = 0; i < numKnotso2; ++i) { 1382 U.add(Float64.ZERO); 1383 } 1384 for (int i = numKnotso2; i < numKnots; ++i) { 1385 U.add(Float64.ONE); 1386 } 1387 KnotVector kv = KnotVector.newInstance(degree, Float64Vector.valueOf(U)); 1388 1389 // Create the curve. 1390 BasicNurbsCurve crv = BasicNurbsCurve.newInstance(cps, kv); 1391 return StackContext.outerCopy(crv); 1392 1393 } finally { 1394 StackContext.exit(); 1395 } 1396 } 1397 1398 /** 1399 * Create a {@link NurbsCurve} of the specified degree that is fit to, interpolates, 1400 * or passes through the specified list of geometry points in the input order. This is 1401 * a global interpolation of the points, meaning that if one of the input points is 1402 * moved, it will affect the entire curve, not just the local portion near the point. 1403 * 1404 * @param degree The degree of the NURBS curve to create (must be > 0 and ≤ the 1405 * number of data points). 1406 * @param points The string of GeomPoint objects to pass the curve through. May not be 1407 * null. 1408 * @return A NurbsCurve fit to the input list of points. 1409 * @throws IllegalArgumentException if the degree is ≤ 0 or ≥ the number of 1410 * points input. 1411 * @throws SingularMatrixException if the decomposed matrix is singular. 1412 * @see #approxPoints 1413 */ 1414 public static BasicNurbsCurve fitPoints(int degree, PointString<? extends GeomPoint> points) 1415 throws IllegalArgumentException, SingularMatrixException, ParameterException { 1416 1417 StackContext.enter(); 1418 try { 1419 if (degree <= 0) 1420 throw new IllegalArgumentException( 1421 MessageFormat.format(RESOURCES.getString("posDegreeErr"), degree)); 1422 1423 int numPoints = points.size(); 1424 if (degree >= numPoints) 1425 throw new IllegalArgumentException(RESOURCES.getString("numPointsLTDegree") + " np = " 1426 + numPoints + ", d = " + degree + "."); 1427 1428 if (points.length().getValue() < EPSC) { 1429 // We have a degenerate curve. 1430 BasicNurbsCurve curve = createPoint(degree, points.size(), points.get(0)); 1431 return StackContext.outerCopy(curve); 1432 } 1433 1434 // Turn the points into an array of parameter values along the curve. 1435 double[] uk = parameterizeString(points); 1436 parameterizationCheck(numPoints, uk); 1437 1438 // Generate the knot vector for interpolation. 1439 KnotVector uKnots = buildInterpKnotVector(uk, numPoints, degree); 1440 1441 // Now fit a curve to the points using the just calculated parameterization and knot vector. 1442 BasicNurbsCurve curve = fitPoints(degree, points, uk, uKnots); 1443 return StackContext.outerCopy(curve); 1444 1445 } finally { 1446 StackContext.exit(); 1447 } 1448 } 1449 1450 /** 1451 * Create a {@link NurbsCurve} of the specified degree that is fit to, interpolates, 1452 * or passes through the specified list of geometry points in the input order. This is 1453 * a global interpolation of the points, meaning that if one of the input points is 1454 * moved, it will affect the entire curve, not just the local portion near the point. 1455 * 1456 * @param degree The degree of the NURBS curve to create (must be > 0 and ≤ the 1457 * number of data points). 1458 * @param points The string of GeomPoint objects to pass the curve through. May not be 1459 * null. 1460 * @param uk The parameterization array to use for each of the input points. Must 1461 * be monotonically increasing and must include 0 and 1. May not be 1462 * null. 1463 * @param uKnots The knot vector to use with the output curve. May not be null. 1464 * @return A NurbsCurve fit to the input list of points. 1465 * @throws IllegalArgumentException if the degree is ≤ 0 or ≥ the number of 1466 * points input or the parameterization array is invalid. 1467 * @throws SingularMatrixException if the decomposed matrix is singular. 1468 * @see #approxPoints(int, int, geomss.geom.PointString, double[], 1469 * geomss.geom.nurbs.KnotVector) 1470 * @see #parameterizeString(geomss.geom.PointString) 1471 * @see #buildInterpKnotVector(double[], int, int) 1472 */ 1473 public static BasicNurbsCurve fitPoints(int degree, PointString<? extends GeomPoint> points, 1474 double[] uk, KnotVector uKnots) throws IllegalArgumentException, SingularMatrixException { 1475 1476 // This follows the basic process outlined here: 1477 // http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/INT-APP/CURVE-INT-global.html 1478 // Basically, you form a matrix of basis function coefficients (N). The data 1479 // points are in matrix (D) and you solve for the control point locations (P). 1480 // D = N * P --> P = N^-1 * D 1481 1482 requireNonNull(uKnots); 1483 if (degree <= 0) { 1484 throw new IllegalArgumentException( 1485 MessageFormat.format(RESOURCES.getString("posDegreeErr"), degree)); 1486 } 1487 1488 int numPoints = points.size(); 1489 if (degree >= numPoints) { 1490 throw new IllegalArgumentException(RESOURCES.getString("numPointsLTDegree") + " np = " 1491 + numPoints + ", d = " + degree + "."); 1492 } 1493 if (uk.length < numPoints) 1494 throw new IllegalArgumentException(MessageFormat.format( 1495 RESOURCES.getString("numItemMissmatchErr"), "points", numPoints, 1496 "parameter values", uk.length)); 1497 if (!MathTools.isApproxEqual(uk[0], 0.0)) 1498 throw new IllegalArgumentException(RESOURCES.getString("paramArrNotZeroErr")); 1499 if (!MathTools.isApproxEqual(uk[numPoints-1],1.0)) 1500 throw new IllegalArgumentException(RESOURCES.getString("paramArrNotOneErr")); 1501 for (int i=1; i < numPoints; ++i) { 1502 if (uk[i] <= uk[i-1]) 1503 throw new IllegalArgumentException(RESOURCES.getString("paramNotMonotonicErr")); 1504 } 1505 1506 StackContext.enter(); 1507 try { 1508 // Create the basis function coefficients matrix N. 1509 RealMatrix Nmat = createInterpBFMatrix(uKnots, uk, numPoints, degree); 1510 1511 // Create the matrix of data points (each row is a point, each column is a dimension). 1512 int numDims = points.getPhyDimension(); 1513 double[] tmpArr = new double[numDims]; 1514 RealMatrix Dmat = MatrixUtils.createRealMatrix(numPoints, numDims); 1515 Unit<Length> refUnits = points.getUnit(); 1516 for (int i = 0; i < numPoints; ++i) { 1517 GeomPoint point = points.get(i).to(refUnits); 1518 Dmat.setRow(i, point.toArray(tmpArr)); 1519 } 1520 1521 // Solve for the matrix of control points: P = N^-1 * D 1522 DecompositionSolver solver = new LUDecomposition(Nmat).getSolver(); 1523 RealMatrix Pmat = solver.solve(Dmat); 1524 1525 // Create a list of control points from the matrix of control point positions. 1526 FastTable<ControlPoint> cpList = pntMatrix2ControlPoints(Pmat, refUnits); 1527 1528 // Make sure that the 1st and the last point are exactly those from the input list. 1529 cpList.set(0, ControlPoint.valueOf(points.get(0).to(refUnits).copyToReal(), 1)); 1530 cpList.set(cpList.size() - 1, ControlPoint.valueOf(points.get(points.size() - 1).to(refUnits).copyToReal(), 1)); 1531 1532 // Create a new curve object. 1533 BasicNurbsCurve crv = BasicNurbsCurve.newInstance(cpList, uKnots); 1534 return StackContext.outerCopy(crv); 1535 1536 } finally { 1537 StackContext.exit(); 1538 } 1539 } 1540 1541 /** 1542 * Create a {@link NurbsCurve} of the specified degree that is fit to, interpolates, 1543 * or passes through the specified list of geometry points constrained to the input 1544 * list of tangent vectors at each point in the input order. This is a global 1545 * interpolation of the points, meaning that if one of the input points is moved, it 1546 * will affect the entire curve, not just the local portion near the point. 1547 * 1548 * @param degree The degree of the NURBS curve to create (must be > 0 and < the 1549 * number of data points). 1550 * @param points The string of Geom Point objects to pass the curve through. May not 1551 * be null. 1552 * @param tangents A list of tangent vectors, one for each point in "points". May not 1553 * be null. 1554 * @param unitFlag A flag indicating that the derivatives are all unit vectors if 1555 * <code>true</code>. 1556 * @param tanlen A multiplicative factor for the tangent (must be greater than 0 or 1557 * the function aborts and should be near 1.0). Used only if 1558 * <code>unitFlag == true</code>. 1559 * @return A NurbsCurve fit to the input list of points and tangent vectors. 1560 * @throws IllegalArgumentException if the degree is ≤ 0 or ≥ the number of 1561 * points input. 1562 * @throws SingularMatrixException if the decomposed matrix is singular. 1563 * @see #approxPoints 1564 */ 1565 public static BasicNurbsCurve fitPoints(int degree, PointString<? extends GeomPoint> points, 1566 List<GeomVector> tangents, boolean unitFlag, double tanlen) 1567 throws IllegalArgumentException, SingularMatrixException, ParameterException { 1568 1569 if (degree <= 0) 1570 throw new IllegalArgumentException( 1571 MessageFormat.format(RESOURCES.getString("posDegreeErr"), degree)); 1572 1573 if (tanlen <= 0) 1574 throw new IllegalArgumentException( 1575 MessageFormat.format(RESOURCES.getString("posMultFactorErr"), tanlen)); 1576 1577 int numPoints = points.size(); 1578 if (degree >= numPoints) 1579 throw new IllegalArgumentException(RESOURCES.getString("numPointsLTDegree") + " np = " 1580 + numPoints + ", d = " + degree + "."); 1581 1582 if (numPoints != tangents.size()) 1583 throw new IllegalArgumentException(MessageFormat.format( 1584 RESOURCES.getString("numItemMissmatchErr"), "points", numPoints, 1585 "tangents", tangents.size())); 1586 1587 StackContext.enter(); 1588 try { 1589 Parameter chordLength = points.length(); 1590 if (chordLength.getValue() < EPSC) { 1591 // We have a degenerate curve. 1592 return StackContext.outerCopy(createPoint(degree, points.get(0))); 1593 } 1594 if (unitFlag) 1595 chordLength = chordLength.times(tanlen); 1596 1597 int n = 2 * numPoints; 1598 1599 // Turn the points into an array of parameter values along the curve. 1600 double[] uk = parameterizeString(points); 1601 parameterizationCheck(numPoints, uk); 1602 1603 // Generate the knot vector for interpolation. 1604 KnotVector uKnots; 1605 int m = n + degree; 1606 int uLength = m + 1; 1607 double U[]; 1608 switch (degree) { 1609 case 2: 1610 U = ArrayFactory.DOUBLES_FACTORY.array(uLength); 1611 for (int i = 0; i <= degree; ++i) { 1612 U[i] = 0; 1613 U[uLength - 1 - i] = 1; 1614 } 1615 for (int i = 0; i < numPoints - 1; ++i) { 1616 U[2 * i + degree] = uk[i]; 1617 U[2 * i + degree + 1] = (uk[i] + uk[i + 1]) / 2; 1618 } 1619 1620 uKnots = array2KnotVector(degree, U, uLength); 1621 break; 1622 1623 case 3: 1624 U = ArrayFactory.DOUBLES_FACTORY.array(uLength); 1625 for (int i = 0; i <= degree; ++i) { 1626 U[i] = 0; 1627 U[uLength - 1 - i] = 1; 1628 } 1629 for (int i = 1; i < numPoints - 1; ++i) { 1630 U[degree + 2 * i] = (2 * uk[i] + uk[i + 1]) / 3; 1631 U[degree + 2 * i + 1] = (uk[i] + 2 * uk[i + 1]) / 3; 1632 } 1633 U[4] = uk[1] / 2; 1634 U[uLength - degree - 2] = (uk[numPoints - 1] + 1) / 2; 1635 1636 uKnots = array2KnotVector(degree, U, uLength); 1637 break; 1638 1639 default: 1640 int uk2Length = 2 * numPoints; 1641 double uk2[] = ArrayFactory.DOUBLES_FACTORY.array(uk2Length); 1642 for (int i = 0; i < numPoints - 1; ++i) { 1643 uk2[2 * i] = uk[i]; 1644 uk2[2 * i + 1] = (uk[i] + uk[i + 1]) / 2; 1645 } 1646 uk2[uk2Length - 2] = (uk2[uk2Length - 1] + uk2[uk2Length - 3]) / 2; 1647 uKnots = buildInterpKnotVector(uk2, uk2Length, degree); 1648 break; 1649 } 1650 1651 // Create the basis function coefficients matrix N. 1652 1653 // Create a list of lists of zeros of the right dimensions. 1654 RealMatrix Nmat = MatrixUtils.createRealMatrix(n, n); 1655 1656 // Now fill in the basis function values. 1657 for (int i = 0; i < numPoints - 1; i++) { 1658 1659 // Get the basis function values and their 1st derivatives for 1660 // this span from the knot vector. 1661 int span = uKnots.findSpan(uk[i]); 1662 double tmpD[][] = uKnots.basisFunctionDerivatives(span, uk[i], 1); 1663 1664 // Copy the basis function values into the matrix. 1665 int row1 = 2 * i; 1666 int row2 = row1 + 1; 1667 int pos = span - degree; 1668 for (int j = 0; j <= degree; ++j, ++pos) { 1669 double value = tmpD[0][j]; 1670 Nmat.setEntry(row1, pos, value); 1671 double deriv = tmpD[1][j]; 1672 Nmat.setEntry(row2, pos, deriv); 1673 } 1674 1675 } 1676 Nmat.setEntry(0, 0, 1); // N(0,0) = 1 1677 Nmat.setEntry(1, 0, -1); // N(1,0) = -1 1678 Nmat.setEntry(1, 1, 1); // N(1,1) = 1 1679 Nmat.setEntry(n - 2, n - 2, -1); // N(n-2,n-2) = -1 1680 Nmat.setEntry(n - 2, n - 1, 1); // N(n-2,n-1) = 1 1681 Nmat.setEntry(n - 1, n - 1, 1); // N(n-1,n-1) = 1 1682 1683 // Create the matrix of data points & derivatives (each row is a point/derivative, 1684 // each column is a dimension). 1685 Unit<Length> refUnits = points.getUnit(); 1686 int numDims = points.getPhyDimension(); 1687 double[] tmpArr = new double[numDims]; 1688 RealMatrix Dmat = MatrixUtils.createRealMatrix(n, numDims); 1689 int row = 0; 1690 for (int i = 0; i < numPoints; ++i, ++row) { 1691 GeomPoint qp = points.get(i).to(refUnits); 1692 GeomVector dp = tangents.get(i); 1693 Dmat.setRow(row, qp.toArray(tmpArr)); 1694 if (unitFlag) 1695 dp = dp.times(chordLength); 1696 dp = (GeomVector)dp.to(refUnits); 1697 ++row; 1698 Dmat.setRow(row, dp.toArray(tmpArr)); 1699 } 1700 1701 double d0Factor = uKnots.getValue(degree + 1) / degree; 1702 double dnFactor = (1 - uKnots.getValue(uKnots.length() - degree - 2)) / degree; 1703 GeomVector dp0 = tangents.get(0); 1704 GeomVector dpn = tangents.get(numPoints - 1); 1705 if (unitFlag) { 1706 dp0 = dp0.times(chordLength).to(refUnits); 1707 dpn = dpn.times(chordLength).to(refUnits); 1708 } 1709 GeomPoint qpn = points.get(numPoints - 1).to(refUnits); 1710 Dmat.setRow(1, dp0.times(d0Factor).toArray(tmpArr)); 1711 Dmat.setRow(n - 2, dpn.times(dnFactor).toArray(tmpArr)); 1712 Dmat.setRow(n - 1, qpn.toArray(tmpArr)); 1713 1714 // Solve for the matrix of control points: P = N^-1 * D 1715 DecompositionSolver solver = new LUDecomposition(Nmat).getSolver(); 1716 RealMatrix Pmat = solver.solve(Dmat); 1717 1718 // Create a list of control points from the matrix of control point positions. 1719 FastTable<ControlPoint> cpList = pntMatrix2ControlPoints(Pmat, refUnits); 1720 1721 // Create a new curve object. 1722 BasicNurbsCurve crv = BasicNurbsCurve.newInstance(cpList, uKnots); 1723 return StackContext.outerCopy(crv); 1724 1725 } finally { 1726 StackContext.exit(); 1727 } 1728 } 1729 1730 /** 1731 * Turn a string of points into a list of parameter values along a curve through those 1732 * points. Uses the chord length method. The returned array was created using 1733 * ArrayFactor.DOUBLES_FACTORY and can be recycled. 1734 * 1735 * @param points A string of points to be parameterized for use in a curve. May not be 1736 * null. 1737 * @return An array of the parameterizations for each point in the points string. 1738 * WARNING: the array will likely be longer than the length of the points 1739 * string and the additional values will be garbage. The returned array can be 1740 * recycled with ArrayFactor.DOUBLES_FACTORY.recycle(arr). 1741 */ 1742 public static double[] parameterizeString(PointString<? extends GeomPoint> points) { 1743 // This uses the Chord Length method which is described here: 1744 // http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/INT-APP/PARA-chord-length.html 1745 1746 Unit<Length> refUnit = points.get(0).getUnit(); 1747 int numPoints = points.size(); 1748 int n = numPoints - 1; 1749 1750 double[] uk = ArrayFactory.DOUBLES_FACTORY.array(numPoints); 1751 uk[0] = 0; 1752 uk[n] = 1; 1753 1754 StackContext.enter(); 1755 try { 1756 double tmp[] = ArrayFactory.DOUBLES_FACTORY.array(n); 1757 1758 double d = 0; 1759 for (int k = 1; k <= n; k++) { 1760 int km1 = k - 1; 1761 Parameter<Length> distance = points.get(k).distance(points.get(km1)); 1762 double distV = distance.getValue(refUnit); 1763 d += distV; 1764 tmp[km1] = distV; 1765 } 1766 if (abs(d) < EPSC) { 1767 // The points form a degenerate curve of zero length. 1768 // Just evenly space out parameters. 1769 for (int i = 1; i < n; ++i) { 1770 double u = (double)(i) / n; 1771 uk[i] = u; 1772 } 1773 1774 } else { 1775 // The points form a curve of finite length. 1776 int nm1 = n - 1; 1777 for (int i = 0, ip1 = 1; i < nm1; ++i, ++ip1) { 1778 uk[ip1] = uk[i] + tmp[i] / d; 1779 } 1780 } 1781 1782 } finally { 1783 StackContext.exit(); 1784 } 1785 1786 return uk; 1787 } 1788 1789 /** 1790 * Check the array of parameterization values and reports any problems by throwing 1791 * a ParameterException. 1792 * 1793 * @param size The number of parameters in the array. 1794 * @param uk The array of parameterization values to check. 1795 * @throws ParameterException if there is a problem with the parameterization array. 1796 */ 1797 public static void parameterizationCheck(int size, double[] uk) throws ParameterException { 1798 double sm1 = uk[0]; 1799 double s = uk[1]; 1800 if (s < sm1) { 1801 System.out.println("sm1 = " + sm1 + ", s = " + s); 1802 throw new ParameterException(RESOURCES.getString("paramNotMonotonicErr")); 1803 } 1804 1805 for (int i=2; i < size; ++i) { 1806 double sp1 = uk[i]; 1807 if (sp1 < s) { 1808 System.out.println("s = " + s + ", sp1 = " + sp1); 1809 throw new ParameterException(RESOURCES.getString("paramNotMonotonicErr")); 1810 } 1811 1812 if (sp1 - s <= MathTools.EPS || s - sm1 <= MathTools.EPS) 1813 throw new ParameterException(RESOURCES.getString("pntsToClose")); 1814 1815 // Prepare for next loop. 1816 sm1 = s; 1817 s = sp1; 1818 } 1819 } 1820 1821 /** 1822 * Create a {@link NurbsCurve} of the specified degree that is fit to, interpolates, 1823 * or passes through the specified list of control points in the input order. This is 1824 * a global interpolation of the points, meaning that if one of the input points is 1825 * moved, it will affect the entire curve, not just the local portion near the point. 1826 * 1827 * @param degree The degree of the NURBS curve to create (must be > 0 and < the 1828 * number of data points). 1829 * @param points The list of Control Point objects to pass the curve through. May not 1830 * be null. 1831 * @param uk The parameterization to use for the input points. May not be null. 1832 * @param uKnots The knot vector to use with the output curve. May not be null. 1833 * @return A NurbsCurve fit to the input list of points. 1834 * @throws IllegalArgumentException if any of the inputs are invalid. 1835 * @throws SingularMatrixException if the decomposed matrix is singular. 1836 * @see #approxPoints 1837 */ 1838 static BasicNurbsCurve fitControlPoints(int degree, List<ControlPoint> points, 1839 double[] uk, KnotVector uKnots) throws IllegalArgumentException, SingularMatrixException { 1840 1841 // This follows the basic process outlined here: 1842 // http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/INT-APP/CURVE-INT-global.html 1843 // Basically, you form a matrix of basis function coefficients (N). The data 1844 // points are in matrix (D) and you solve for the control point locations (P). 1845 // D = N * P --> P = N^-1 * D 1846 requireNonNull(uk); 1847 requireNonNull(uKnots); 1848 if (degree <= 0) 1849 throw new IllegalArgumentException( 1850 MessageFormat.format(RESOURCES.getString("posDegreeErr"), degree)); 1851 1852 int numPoints = points.size(); 1853 if (degree >= numPoints) 1854 throw new IllegalArgumentException(RESOURCES.getString("numPointsLTDegree") + " np = " 1855 + numPoints + ", d = " + degree + "."); 1856 1857 StackContext.enter(); 1858 try { 1859 // Create the basis function coefficients matrix N. 1860 RealMatrix Nmat = createInterpBFMatrix(uKnots, uk, numPoints, degree); 1861 1862 // Create the matrix of data points (each row is a point, each column is a dimension). 1863 int numDims = points.get(0).getPhyDimension(); 1864 double[] tmpArr = new double[numDims + 1]; 1865 RealMatrix Dmat = MatrixUtils.createRealMatrix(numPoints, numDims + 1); 1866 Unit<Length> refUnits = points.get(0).getUnit(); 1867 for (int i = 0; i < numPoints; ++i) { 1868 ControlPoint point = points.get(i).to(refUnits); 1869 Dmat.setRow(i, point.toArray(tmpArr)); 1870 } 1871 1872 // Solve for the matrix of control points: P = N^-1 * D 1873 DecompositionSolver solver = new LUDecomposition(Nmat).getSolver(); 1874 RealMatrix Pmat = solver.solve(Dmat); 1875 1876 // Create a list of control points from the matrix of control point coordinates. 1877 FastTable<ControlPoint> cpList = cpMatrix2ControlPoints(Pmat, refUnits); 1878 1879 // Create a new curve object. 1880 BasicNurbsCurve crv = BasicNurbsCurve.newInstance(cpList, uKnots); 1881 return StackContext.outerCopy(crv); 1882 1883 } finally { 1884 StackContext.exit(); 1885 } 1886 } 1887 1888 /** 1889 * Create a vector of knots to go with the list of parameter values for a curve that 1890 * passes through a list of points. This uses the knot averaging technique. 1891 * 1892 * @param uk The parameter values for each point in the curve. May not be null. 1893 * @param numPoints The total number of points in the curve (and parameter array). 1894 * @param p The degree of the NURBS curve. 1895 * @return A vector of knots of degree p that are parameterized for use with the input 1896 * list of point parameter values. 1897 */ 1898 public static KnotVector buildInterpKnotVector(double uk[], int numPoints, int p) { 1899 requireNonNull(uk); 1900 int m = numPoints + p; 1901 int n = numPoints - 1; 1902 int uLength = m + 1; 1903 double u[] = ArrayFactory.DOUBLES_FACTORY.array(uLength); 1904 1905 for (int i = 0; i <= p; i++) { 1906 u[i] = 0; 1907 u[uLength - 1 - i] = 1; 1908 } 1909 for (int j = 1; j <= n - p; j++) { 1910 double sum = 0; 1911 for (int i = j; i <= j + p - 1; i++) { 1912 sum += uk[i]; 1913 } 1914 u[j + p] = sum / p; 1915 } 1916 1917 KnotVector kv = array2KnotVector(p, u, uLength); 1918 ArrayFactory.DOUBLES_FACTORY.recycle(u); 1919 1920 return kv; 1921 } 1922 1923 /** 1924 * Convert an array of doubles into a KnotVector. 1925 * 1926 * @param p The degree of the KnotVector to create. 1927 * @param u The array of doubles to be converted. 1928 * @param uLength The number of elements from "u" to use in the knot vector. 1929 * @return A KnotVector of degree "p" made up from the "uLength" values of "u". 1930 */ 1931 private static KnotVector array2KnotVector(int p, double[] u, int uLength) { 1932 1933 StackContext.enter(); 1934 try { 1935 FastTable<Float64> kvValues = FastTable.newInstance(); 1936 1937 for (int i = 0; i < uLength; ++i) 1938 kvValues.add(Float64.valueOf(u[i])); 1939 KnotVector kv = KnotVector.newInstance(p, Float64Vector.valueOf(kvValues)); 1940 return StackContext.outerCopy(kv); 1941 1942 } finally { 1943 StackContext.exit(); 1944 } 1945 } 1946 1947 /** 1948 * Create a numPoints by numPoints sized matrix that contains the basis function 1949 * values for each control point for each parametric position input. 1950 * <pre> 1951 * N = [ N0,p(t0) N1,p(t0) ... Nn,p(t0)] 1952 * [ N0,p(t1) N1,p(t1) ... Nn,p(t1)] 1953 * [ . . . ] 1954 * [ N0,p(tn) N1,p(tn) ... Nn,p(tn)] 1955 * </pre> 1956 * 1957 * @param knots The knot vector for a NURBS curve. 1958 * @param tk A list of parametric parameter values. 1959 * @param numPoints The number of parametric values (the size of the data in "tk"). 1960 * @param degree The degree of the NURBS curve. 1961 */ 1962 private static RealMatrix createInterpBFMatrix(KnotVector knots, double[] tk, int numPoints, int degree) { 1963 1964 // Create an array of zeros of the right dimensions. 1965 RealMatrix Nmat = MatrixUtils.createRealMatrix(numPoints, numPoints); 1966 1967 // Now fill in the basis function values. 1968 for (int i = 0; i < numPoints; ++i) { 1969 1970 // Get the basis function values for this span from the knot vector. 1971 int span = knots.findSpan(tk[i]); 1972 double tmp[] = knots.basisFunctions(span, tk[i]); 1973 1974 // Copy the basis function values into the matrix. 1975 int pos = span - degree; 1976 for (int j = 0; j <= degree; ++j, ++pos) { 1977 Nmat.setEntry(i, pos, tmp[j]); 1978 } 1979 1980 ArrayFactory.DOUBLES_FACTORY.recycle(tmp); 1981 } 1982 1983 return Nmat; 1984 } 1985 1986 /** 1987 * Return a list of control points from the input matrix of control point positions 1988 * (points in rows, dimensions in columns). 1989 */ 1990 private static FastTable<ControlPoint> pntMatrix2ControlPoints(RealMatrix Pmat, Unit units) { 1991 // Create a list of control points from the matrix of control point positions. 1992 FastTable<ControlPoint> cpList = FastTable.newInstance(); 1993 int numPoints = Pmat.getRowDimension(); 1994 for (int i = 0; i < numPoints; ++i) { 1995 Point pnt = Point.valueOf(units, Pmat.getRow(i)); 1996 cpList.add(ControlPoint.valueOf(pnt, 1)); 1997 } 1998 return cpList; 1999 } 2000 2001 /** 2002 * Return a list of control points from the input matrix of control point coordinates 2003 * (points in rows, dimensions in columns). 2004 */ 2005 private static FastTable<ControlPoint> cpMatrix2ControlPoints(RealMatrix Pmat, Unit units) { 2006 2007 // Create a list of control points from the matrix of control point positions. 2008 FastTable<ControlPoint> cpList = FastTable.newInstance(); 2009 int numPoints = Pmat.getRowDimension(); 2010 for (int i = 0; i < numPoints; i++) { 2011 ControlPoint pnt = ControlPoint.valueOf(units, Pmat.getRow(i)); 2012 cpList.add(pnt); 2013 } 2014 2015 return cpList; 2016 } 2017 2018 /** 2019 * Create a {@link NurbsCurve} of the specified degree that approximates, in a least 2020 * squared error sense, the specified list of geometry points in the input order. This 2021 * is a global approximation of the points, meaning that if one of the input points is 2022 * moved, it will affect the entire curve, not just the local portion near the point. 2023 * 2024 * @param degree The degree of the NURBS curve to create (must be > 0 and < the 2025 * number of data points). 2026 * @param nCP The number of control points in the new curve (must be > 1 and 2027 * < the number of data points). 2028 * @param points The string of GeomPoint objects to be approximated by the curve. May 2029 * not be null. 2030 * @return A NurbsCurve approximation to the input list of points. 2031 * @throws SingularMatrixException if the decomposed matrix is singular. 2032 * @throws IllegalArgumentException if the any of the inputs is invalid. 2033 * @see #fitPoints 2034 */ 2035 public static BasicNurbsCurve approxPoints(int degree, int nCP, PointString<? extends GeomPoint> points) 2036 throws IllegalArgumentException, SingularMatrixException, ParameterException { 2037 2038 if (degree <= 0) 2039 throw new IllegalArgumentException( 2040 MessageFormat.format(RESOURCES.getString("posDegreeErr"), degree)); 2041 2042 int numPoints = points.size(); 2043 if (degree >= numPoints) 2044 throw new IllegalArgumentException(RESOURCES.getString("numPointsLTDegree") + " np = " 2045 + numPoints + ", d = " + degree + "."); 2046 2047 if (nCP <= 1) 2048 throw new IllegalArgumentException( 2049 MessageFormat.format(RESOURCES.getString("numPointsLEOneErr"), "S")); 2050 2051 if (nCP >= numPoints) 2052 throw new IllegalArgumentException(RESOURCES.getString("numPointsGTnumCPErr")); 2053 2054 // Check for a degenerate curve that collapses to a point. 2055 if (points.length().getValue() < EPSC) { 2056 return createPoint(degree, points.get(0)); 2057 } 2058 2059 // Check for the degenerate case of a degree=1 curve having only 2 points (a straight line). 2060 if (nCP - 2 <= 0) { 2061 return createLine(points.get(0), points.get(numPoints - 1)); 2062 } 2063 2064 // Turn the points into an array of parameter values along the curve. 2065 double[] ub = parameterizeString(points); 2066 parameterizationCheck(numPoints, ub); 2067 2068 // Generate the knot vector for approximation. 2069 KnotVector uKnots = buildApproxKnotVector(ub, numPoints, nCP, degree); 2070 2071 // Create the approximating curve. 2072 BasicNurbsCurve crv = approxPoints(degree, nCP, points, ub, uKnots); 2073 ArrayFactory.DOUBLES_FACTORY.recycle(ub); 2074 2075 return crv; 2076 } 2077 2078 /** 2079 * Create a {@link NurbsCurve} of the specified degree that approximates, in a least 2080 * squared error sense, the specified list of geometry points in the input order. This 2081 * is a global approximation of the points, meaning that if one of the input points is 2082 * moved, it will affect the entire curve, not just the local portion near the point. 2083 * 2084 * @param degree The degree of the NURBS curve to create (must be > 0 and < the 2085 * number of data points). 2086 * @param nCP The number of control points in the new curve (must be > 1 and < 2087 * the number of data points). 2088 * @param points The string of GeomPoint objects to be approximated by the curve. 2089 * @return A NurbsCurve approximation to the input list of points. 2090 * @throws SingularMatrixException if the decomposed matrix is singular. 2091 * @throws IllegalArgumentException if the any of the inputs is invalid. 2092 */ 2093 private static BasicNurbsCurve approxPoints(int degree, int nCP, PointString<? extends GeomPoint> points, 2094 double[] ub, KnotVector uKnots) throws IllegalArgumentException, SingularMatrixException { 2095 // Reference: Piegl, L., Tiller, W., The Nurbs Book, 2nd Edition, Springer-Verlag, Berlin, 1997. 2096 // Based on discussion in section 9.4.1 on page 410. 2097 2098 StackContext.enter(); 2099 try { 2100 // Determine the units used by the list of points (or at least by the 1st one). 2101 Unit<Length> refUnits = points.getUnit(); 2102 2103 // Create the basis function coefficients matrix N. 2104 int numPoints = points.size(); 2105 RealMatrix Nmat = createApproxBFMatrix(uKnots, ub, numPoints, nCP, degree); 2106 2107 // Create an Rk vector that represents Eqn. 9.63, pg. 411. 2108 int nm1 = nCP - 1; 2109 int mm1 = numPoints - 1; 2110 Point Q0 = points.get(0).to(refUnits).immutable(); 2111 Point Qmm1 = points.get(mm1).to(refUnits).immutable(); 2112 FastTable<Point> Rk = FastTable.newInstance(); 2113 for (int i = 0; i < numPoints; i++) { 2114 // Rk[i] = Q[i]-N(i,0)*Q[0]-N(i,n-1)*Q[m-1]; 2115 GeomPoint Qi = points.get(i).to(refUnits); 2116 double Ni0 = Nmat.getEntry(i, 0); 2117 double Ninm1 = Nmat.getEntry(i, nm1); // N(i,n-1) 2118 Point Rki = Qi.minus(Q0.times(Ni0)).minus(Qmm1.times(Ninm1)); 2119 Rk.add(Rki); 2120 } 2121 2122 // Setup the R matrix of data points (each row is a point, each column is a dimension). 2123 int numDims = Q0.getPhyDimension(); 2124 RealMatrix Rmat = MatrixUtils.createRealMatrix(nCP, numDims); 2125 double[] tmpArr = new double[numDims]; 2126 Rmat.setRow(0, Q0.toArray(tmpArr)); // R[0] = Q[0]; 2127 Rmat.setRow(nm1, Qmm1.toArray(tmpArr)); // R[n-1] = Q[m-1]; 2128 Point ZERO = Point.newInstance(numDims, refUnits); 2129 for (int i = 0; i < nCP; ++i) { 2130 Point Ri = ZERO; 2131 2132 for (int k = 1; k < mm1; ++k) { 2133 // R[i] += N(k,i)*rk[k]; 2134 double Nki = Nmat.getEntry(k, i); 2135 Ri = Ri.plus(Rk.get(k).times(Nki)); 2136 } 2137 2138 Rmat.setRow(i, Ri.toArray(tmpArr)); 2139 } 2140 2141 // Solve for the matrix of control points: N^T*N*P = R ==> P = (N^T*N)^-1 * R = M^-1 * R 2142 RealMatrix Mmat = Nmat.transpose().multiply(Nmat); // M = N^T*N 2143 DecompositionSolver solver = new LUDecomposition(Mmat).getSolver(); 2144 RealMatrix Pmat = solver.solve(Rmat); 2145 2146 // Create a list of control points from the matrix of control point positions. 2147 FastTable<ControlPoint> cpList = pntMatrix2ControlPoints(Pmat, refUnits); 2148 cpList.set(0, ControlPoint.valueOf(Q0, 1)); 2149 cpList.set(nCP - 1, ControlPoint.valueOf(Qmm1, 1)); 2150 2151 // Create a new curve object. 2152 BasicNurbsCurve crv = BasicNurbsCurve.newInstance(cpList, uKnots); 2153 return StackContext.outerCopy(crv); 2154 2155 } finally { 2156 StackContext.exit(); 2157 } 2158 } 2159 2160 /** 2161 * Create a vector of knots to go with the list of parameter values for a curve that 2162 * approximates a list of points. 2163 * 2164 * @param ub The parameter values for each point in the curve. May not be null. 2165 * @param numPoints The total number of points in the curve (and parameter array). 2166 * @param nCP The number of control points in the new curve. 2167 * @param p The degree of the NURBS curve. 2168 */ 2169 static KnotVector buildApproxKnotVector(double ub[], int numPoints, int nCP, int p) { 2170 // Similar to Eqn 9.69 in NURBS book (pg. 412), but modified slighlty. 2171 2172 requireNonNull(ub); 2173 StackContext.enter(); 2174 try { 2175 int uLength = nCP + p + 1; 2176 2177 // Initialize the ends of the knot vector. 2178 double u[] = ArrayFactory.DOUBLES_FACTORY.array(uLength); 2179 for (int i = 0; i <= p; i++) { 2180 u[i] = 0; 2181 u[uLength - 1 - i] = 1; 2182 } 2183 2184 double d = ((double)numPoints) / nCP; 2185 for (int j = 1; j < nCP - p; j++) { 2186 u[p + j] = 0; 2187 2188 for (int k = j; k < j + p; ++k) { 2189 double kd = k * d; 2190 int i = (int)kd; 2191 double a = kd - i; 2192 int i2 = (int)((k - 1) * d); 2193 u[p + j] += a * ub[i2] + (1 - a) * ub[i]; 2194 } 2195 u[p + j] /= p; 2196 } 2197 2198 FastTable<Float64> kvValues = FastTable.newInstance(); 2199 for (int i = 0; i < uLength; ++i) { 2200 kvValues.add(Float64.valueOf(u[i])); 2201 } 2202 2203 KnotVector kv = KnotVector.newInstance(p, Float64Vector.valueOf(kvValues)); 2204 return StackContext.outerCopy(kv); 2205 2206 } finally { 2207 StackContext.exit(); 2208 } 2209 } 2210 2211 /** 2212 * Create a numPoints by nCP sized matrix that contains the basis function values for 2213 * each control point for each parametric position input. 2214 * <pre> 2215 * N = [ N0,p(u0) N1,p(u0) ... Nn,p(u0)] 2216 * [ N0,p(u1) N1,p(u1) ... Nn,p(u1)] 2217 * [ . . . ] 2218 * [ N0,p(um) N1,p(um) ... Nn,p(um)] 2219 * </pre> 2220 * 2221 * @param knots The knot vector for a NURBS curve. 2222 * @param ub A list of parametric parameter values for each data point. 2223 * @param numPoints The number of points in the input data set. 2224 * @param nCP The number of control points for the output curve. 2225 * @param degree The degree of the NURBS curve. 2226 */ 2227 private static RealMatrix createApproxBFMatrix(KnotVector knots, double[] ub, int numPoints, int nCP, int degree) { 2228 // Builds up Eqn. 9.66 in the NURBS book (pg. 411). 2229 2230 RealMatrix Nmat = MatrixUtils.createRealMatrix(numPoints, nCP); 2231 Nmat.setEntry(0, 0, 1); // N(0,0) = 1.0; 2232 Nmat.setEntry(numPoints - 1, nCP - 1, 1); // N(m-1,n-1) = 1.0; 2233 2234 // Fill in the basis function values. 2235 for (int i = 0; i < numPoints; i++) { 2236 2237 // Get the basis function values for the ith span from the knot vector. 2238 int span = knots.findSpan(ub[i]); 2239 double funs[] = knots.basisFunctions(span, ub[i]); 2240 2241 // Copy the basis function values into the matrix. 2242 int pos = span - degree; 2243 for (int j = 0; j <= degree; ++j, ++pos) { 2244 Nmat.setEntry(i, pos, funs[j]); 2245 } 2246 2247 ArrayFactory.DOUBLES_FACTORY.recycle(funs); 2248 } 2249 2250 return Nmat; 2251 2252 } 2253 2254}