001/** 002 * SurfaceFactory -- A factory for creating NURBS surfaces. 003 * 004 * Copyright (C) 2010-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 java.text.MessageFormat; 026import java.util.List; 027import static java.util.Objects.requireNonNull; 028import java.util.ResourceBundle; 029import javax.measure.quantity.Angle; 030import javax.measure.quantity.Dimensionless; 031import javax.measure.quantity.Length; 032import javax.measure.unit.SI; 033import javax.measure.unit.Unit; 034import javolution.context.ArrayFactory; 035import javolution.context.ImmortalContext; 036import javolution.context.StackContext; 037import static javolution.lang.MathLib.*; 038import javolution.util.FastTable; 039import org.jscience.mathematics.number.Float64; 040import org.jscience.mathematics.vector.Float64Vector; 041 042/** 043 * A collection of methods for creating NURBS surfaces 044 * 045 * <p> Modified by: Joseph A. Huwaldt </p> 046 * 047 * @author Joseph A. Huwaldt, Date: June 15, 2010 048 * @version February 17, 2025 049 */ 050public final class SurfaceFactory { 051 052 /** 053 * The resource bundle for this package. 054 */ 055 private static final ResourceBundle RESOURCES = geomss.geom.AbstractGeomElement.RESOURCES; 056 057 /** 058 * The length threshold for a set of points being considered a degenerate curve. 059 */ 060 private static final double EPSC = 1e-10; 061 062 // A generic geometry space tolerance 063 private static final Parameter<Length> GTOL; 064 065 static { 066 ImmortalContext.enter(); 067 try { 068 GTOL = Parameter.valueOf(1e-8, SI.METER); 069 } finally { 070 ImmortalContext.exit(); 071 } 072 } 073 074 private SurfaceFactory() { } 075 076 /** 077 * Create a sphere {@link NurbsSurface} about the specified center point with the 078 * specified radius. 079 * 080 * @param center The center of the sphere. May not be null. 081 * @param radius The radius of the sphere. May not be null. 082 * @return A BasicNurbsSurface representing the sphere. 083 */ 084 public static BasicNurbsSurface createSphere(GeomPoint center, Parameter<Length> radius) { 085 requireNonNull(radius); 086 StackContext.enter(); 087 try { 088 // Make sure center point is at least 3D. 089 Point o = center.immutable(); 090 int numDims = o.getPhyDimension(); 091 if (numDims < 3) { 092 o = o.toDimension(3); 093 numDims = 3; 094 } 095 096 // Create a polar axis. 097 Vector<Dimensionless> polarAxis = GeomUtil.makeZVector(numDims); 098 polarAxis.setOrigin(o); 099 100 // Create a normal vector for the semi-circle. 101 Vector<Dimensionless> n = GeomUtil.makeYVector(numDims); 102 103 // Create the sphere cross-section. 104 NurbsCurve xsec = CurveFactory.createSemiCircle(o, radius, n); 105 106 // Revolve into a sphere. 107 BasicNurbsSurface srf = createRevolvedSurface(polarAxis, xsec, 0, TWO_PI); 108 return StackContext.outerCopy(srf); 109 110 } finally { 111 StackContext.exit(); 112 } 113 } 114 115 /** 116 * Create a torus {@link NurbsSurface} about the specified center point and axis. The 117 * surface will be parameterized such that the "s" direction is around the 118 * cross-section of the torus and the t-direction is around the axis. 119 * 120 * @param center The center of the torus. May not be null. 121 * @param axis The axis to place the torus about. May not be null. 122 * @param innerRadius The inner radius of the torus. May not be null. 123 * @param outerRadius The outer radius of the torus. May not be null. 124 * @return A BasicNurbsSurface representing the torus. 125 */ 126 public static BasicNurbsSurface createTorus(GeomPoint center, GeomVector<Dimensionless> axis, 127 Parameter<Length> innerRadius, Parameter<Length> outerRadius) throws DimensionException { 128 requireNonNull(innerRadius); 129 130 // Make sure the input axis is at least 3 dimensional. 131 int numDims = axis.getPhyDimension(); 132 if (numDims < 3) 133 axis = axis.toDimension(3); 134 numDims = center.getPhyDimension(); 135 if (numDims < 3) 136 center = center.toDimension(3); 137 if (center.getPhyDimension() != axis.getPhyDimension()) 138 throw new DimensionException(RESOURCES.getString("pointAxisDimErr")); 139 140 StackContext.enter(); 141 try { 142 // Create the yhat & xhat vector (Y & X axes, orthogonal to axis and each other, 143 // in the plane of the torus). 144 Vector<Dimensionless> yhat = GeomUtil.calcYHat(axis); 145 Vector<Dimensionless> xhat = GeomUtil.calcXHat(axis, yhat); 146 147 // The radius of the torus cross-section. 148 Parameter<Length> xsec_radius = outerRadius.minus(innerRadius).divide(2); 149 150 // Find the center of the x-section. 151 Point xsec_c = center.plus(Point.valueOf(xhat.times(innerRadius.plus(xsec_radius)))); 152 153 // Create a x-section circle. 154 NurbsCurve xsec = CurveFactory.createCircle(xsec_c, xsec_radius, yhat); 155 156 // Revolve into a torus. 157 BasicNurbsSurface srf = createRevolvedSurface(axis, xsec, 0, TWO_PI); 158 return StackContext.outerCopy(srf); 159 160 } finally { 161 StackContext.exit(); 162 } 163 } 164 165 /** 166 * Create a surface of revolution {@link NurbsSurface} by sweeping the specified curve 167 * about the specified axis. The surface will be parameterized such that the curve 168 * defines the "s" direction and the revolution around the axis is the "t" direction. 169 * 170 * @param axis The axis to sweep the curve about (the vector gives direction; use 171 * axis.setOrigin() to set a point on the axis). May not be null. 172 * @param curve The curve to be revolved around the axis. May not be null. 173 * @return A BasicNurbsSurface for the input curve revolved around the input axis. 174 */ 175 public static BasicNurbsSurface createRevolvedSurface(GeomVector<Dimensionless> axis, NurbsCurve curve) 176 throws DimensionException { 177 return createRevolvedSurface(axis, curve, Parameter.ZERO_ANGLE, Parameter.TWOPI_ANGLE); 178 } 179 180 /** 181 * Create a surface of revolution {@link NurbsSurface} by sweeping the specified curve 182 * about the specified axis with the specified curve and angular start and stop values 183 * (relative to the curve and axis). The surface will be parameterized such that the 184 * curve defines the "s" direction and the revolution around the axis is the "t" 185 * direction. 186 * 187 * @param axis The axis to sweep the curve about (the vector gives direction; use 188 * axis.setOrigin() to set a point on the axis). May not be null. 189 * @param curve The curve to be revolved around the axis. May not be null. 190 * @param theta The angle to revolve the curve through (starting at the curve's 191 * position). May not be null. 192 * @return A BasicNurbsSurface for the input curve revolved around the input axis. 193 */ 194 public static BasicNurbsSurface createRevolvedSurface(GeomVector<Dimensionless> axis, NurbsCurve curve, 195 Parameter<Angle> theta) throws DimensionException { 196 return createRevolvedSurface(axis, curve, Parameter.ZERO_ANGLE, theta); 197 } 198 199 /** 200 * Create a surface of revolution {@link NurbsSurface} by sweeping the specified curve 201 * about the specified axis with the specified curve and angular start and stop values 202 * (relative to the curve and axis). The surface will be parameterized such that the 203 * curve defines the "s" direction and the revolution around the axis is the "t" 204 * direction. 205 * 206 * @param axis The axis to sweep the curve about (the vector gives direction; use 207 * axis.setOrigin() to set a point on the axis). May not be null. 208 * @param curve The curve to be revolved around the axis. May not be null. 209 * @param thetaS The start angle to revolve the curve through (starting at the curve's 210 * position). May not be null. 211 * @param thetaE The end angle to revolve the curve through (starting at the curve's 212 * position). If the end angle is smaller than the start angle, it is 213 * increased by increments of 2*PI until it is larger than the start 214 * angle. May not be null. 215 * @return A BasicNurbsSurface for the input curve revolved around the input axis. 216 */ 217 public static BasicNurbsSurface createRevolvedSurface(GeomVector<Dimensionless> axis, NurbsCurve curve, 218 Parameter<Angle> thetaS, Parameter<Angle> thetaE) throws DimensionException { 219 int numDims = curve.getPhyDimension(); 220 if (numDims < 3) 221 throw new DimensionException( 222 MessageFormat.format(RESOURCES.getString("dimensionNotAtLeast3"), "curve", numDims)); 223 if (axis.getPhyDimension() != numDims) 224 throw new DimensionException(RESOURCES.getString("crvAxisDimErr")); 225 226 return createRevolvedSurface(axis, curve, thetaS.getValue(SI.RADIAN), thetaE.getValue(SI.RADIAN)); 227 } 228 229 /** 230 * Create a surface of revolution {@link NurbsSurface} by sweeping the specified curve 231 * about the specified axis with the specified curve and angular start and stop values 232 * (relative to the curve and axis). The surface will be parameterized such that the 233 * curve defines the "s" direction and the revolution around the axis is the "t" 234 * direction. 235 * 236 * @param axis The axis to sweep the curve about (the vector gives direction; use 237 * axis.setOrigin() to set a point on the axis). May not be null. 238 * @param curve The curve to be revolved around the axis. May not be null. 239 * @param ths The start angle to revolve the curve through (starting at the curve's 240 * position) in radians. 241 * @param the The end angle to revolve the curve through (starting at the curve's 242 * position) in radians. If the end angle is smaller than the start 243 * angle, it is increased by increments of 2*PI until it is larger than 244 * the start angle. 245 * @return A BasicNurbsSurface for the input curve revolved around the input axis. 246 */ 247 private static BasicNurbsSurface createRevolvedSurface(GeomVector<Dimensionless> axis, NurbsCurve curve, 248 double ths, double the) { 249 requireNonNull(axis); 250 requireNonNull(curve); 251 252 StackContext.enter(); 253 try { 254 // Get ths and the in the range 0 - 2*PI. 255 ths = GeomUtil.bound2Pi(ths); 256 the = GeomUtil.bound2Pi(the); 257 258 // Make sure the end angle is > than the start angle. 259 if (the <= ths) 260 the += TWO_PI; 261 double theta = the - ths; 262 263 // How many arc segments are required? 264 int narcs = 4; 265 if (theta <= HALF_PI) 266 narcs = 1; 267 else if (theta <= PI) 268 narcs = 2; 269 else if (theta <= PI + HALF_PI) 270 narcs = 3; 271 272 // Pre-calculate sine and cosine tables. 273 double[] cos = ArrayFactory.DOUBLES_FACTORY.array(narcs + 1); 274 double[] sin = ArrayFactory.DOUBLES_FACTORY.array(narcs + 1); 275 double angle = ths; 276 double dtheta = theta / narcs; 277 for (int i = 0; i <= narcs; ++i) { 278 sin[i] = sin(angle); 279 cos[i] = cos(angle); 280 angle += dtheta; 281 } 282 double sinths = sin(ths), cosths = cos(ths); 283 double wm = cos(dtheta / 2); 284 285 // Get the curve's control points. 286 List<ControlPoint> pj = curve.getControlPoints(); 287 int pjSize = pj.size(); 288 289 // Create a list of control points for the surface being generated. 290 FastTable<FastTable<ControlPoint>> pij = FastTable.newInstance(); 291 int num = 2 * narcs + 1; 292 for (int j = 0; j < num; ++j) { 293 FastTable<ControlPoint> tbl = FastTable.newInstance(); 294 tbl.setSize(pjSize); 295 pij.add(tbl); 296 } 297 298 // Some handy constants for use later. 299 Unit<Length> unit = curve.getUnit(); 300 int numDims = curve.getPhyDimension(); 301 GeomPoint axisOrigin = axis.getOrigin(); 302 303 // Loop over each of the curves control points. 304 MutablePoint p1 = MutablePoint.newInstance(numDims, unit); 305 for (int j = 0; j < pjSize; ++j) { 306 307 // Get the curve control point. 308 ControlPoint pjj = pj.get(j); 309 double pjjw = pjj.getWeight(); 310 Point pjjp = pjj.getPoint(); 311 312 // Determine the center of the circle for this curve control point. 313 Point center = GeomUtil.pointLineClosest(pjjp, axisOrigin, axis).to(unit); 314 315 // Determine the X-axis direction. 316 Vector<Length> xhatd = pjjp.minus(center).toGeomVector(); 317 Parameter<Length> radius = xhatd.mag(); // Circle radius 318 Vector<Dimensionless> xhat; 319 if (radius.getValue(unit) < EPSC) { 320 xhat = GeomUtil.calcYHat(axis); // X-axis in arbitrary direction perpendicular to axis. 321 radius = Parameter.ZERO_LENGTH; 322 } else { 323 xhat = xhatd.toUnitVector(); 324 } 325 326 // Determine the Y-axis direction. 327 Vector<Dimensionless> yhat = (Vector<Dimensionless>)axis.cross(xhat); 328 329 // Create the starting control point 330 Vector<Dimensionless> p0v = xhat.times(cosths).plus(yhat.times(sinths)); 331 MutablePoint p0 = MutablePoint.valueOf(center.plus(Point.valueOf(p0v.times(radius)))); 332 pij.get(0).set(j, ControlPoint.valueOf(p0.immutable(), pjjw)); 333 334 // Create starting perpendicular vector 335 Vector<Dimensionless> t0 = xhat.times(-sinths).plus(yhat.times(cosths)); 336 337 // Loop over all the required arcs around the revolution. 338 int index = 0; 339 for (int i = 1; i <= narcs; i++) { 340 341 // Create segment end control point: p2 = r*(xhat*cos[i] + yhat*sin[i]) 342 p0v = xhat.times(cos[i]).plus(yhat.times(sin[i])); 343 Point p2 = center.plus(Point.valueOf(p0v.times(radius))); 344 pij.get(index + 2).set(j, ControlPoint.valueOf(p2, pjjw)); 345 346 // Create end perpendicular vector: t2 = -xhat*r*sin[i] + yhat*r*cos[i] 347 Vector<Dimensionless> t2 = xhat.times(-sin[i]).plus(yhat.times(cos[i])); 348 349 // Find the intersection of the two perpendiculars. 350 GeomUtil.lineLineIntersect(p0, t0, p2, t2, GTOL, null, p1, null); 351 352 // Create the intermediate control point. 353 pij.get(index + 1).set(j, ControlPoint.valueOf(p1.immutable(), wm * pjjw)); 354 355 index += 2; 356 if (i < narcs) { 357 p0.set(p2); 358 t0 = t2; 359 } 360 } 361 } 362 363 // Turn the list into a network of control points for the surface. 364 ControlPointNet cpnet = ControlPointNet.valueOf(pij); 365 366 // Create the knot vector around the revolution. 367 FastTable<Float64> tKnots = FastTable.newInstance(); 368 int j = 3 + 2 * (narcs - 1); 369 tKnots.setSize(j + 3); 370 for (int i = 0; i < 3; i++) { 371 tKnots.set(i, Float64.ZERO); 372 tKnots.set(i + j, Float64.ONE); 373 } 374 switch (narcs) { 375 case 2: 376 tKnots.set(3, Float64.valueOf(0.5)); 377 tKnots.set(4, Float64.valueOf(0.5)); 378 break; 379 case 3: 380 tKnots.set(3, Float64.valueOf(1. / 3.)); 381 tKnots.set(4, Float64.valueOf(1. / 3.)); 382 tKnots.set(5, Float64.valueOf(2. / 3.)); 383 tKnots.set(6, Float64.valueOf(2. / 3.)); 384 break; 385 case 4: 386 tKnots.set(3, Float64.valueOf(0.25)); 387 tKnots.set(4, Float64.valueOf(0.25)); 388 tKnots.set(5, Float64.valueOf(0.5)); 389 tKnots.set(6, Float64.valueOf(0.5)); 390 tKnots.set(7, Float64.valueOf(0.75)); 391 tKnots.set(8, Float64.valueOf(0.75)); 392 break; 393 } 394 395 // Create the NURBS surface. 396 KnotVector kvt = KnotVector.newInstance(2, Float64Vector.valueOf(tKnots)); 397 BasicNurbsSurface srf = BasicNurbsSurface.newInstance(cpnet, curve.getKnotVector(), kvt); 398 return StackContext.outerCopy(srf); 399 400 } finally { 401 StackContext.exit(); 402 } 403 } 404 405 /** 406 * Create a surface that interpolates, is skinned or lofted, through the input list of 407 * curves. This creates a surface in NURBS form that passes through the input list of 408 * curves to within the specified tolerance. The degree in the parametric "s" 409 * direction (along the input curves) will be determined by the highest degree input 410 * curve. The degree in the "t" direction (across the input curves) is specified. If 411 * the input curves are all NURBS curves, then the "tol" parameter is ignored. 412 * Otherwise, "tol" is used to make a NURBS approximation of the input curves. It is 413 * assumed that all the input curves are oriented in an appropriate direction and 414 * order to create the desired skinned surface. 415 * 416 * @param crvs A list of curves to be skinned. May not be null. 417 * @param tDegree The degree in the t-direction of the NURBS surface to create (must 418 * be > 0 and < the number of curves in the list). 419 * @param tol The greatest possible difference between the input list of curves 420 * and the NURBS surface returned. May not be null or zero. 421 * @return A BasicNurbsSurface that interpolates the input list of curves. 422 */ 423 public static BasicNurbsSurface createSkinnedSurface(List<? extends Curve> crvs, int tDegree, Parameter<Length> tol) { 424 requireNonNull(tol); 425 int numCrvs = crvs.size(); 426 if (numCrvs < 2) 427 throw new IllegalArgumentException(RESOURCES.getString("skinnedNotEnoughCrvs")); 428 if (tDegree >= numCrvs) 429 throw new IllegalArgumentException(MessageFormat.format( 430 RESOURCES.getString("incDefiningCrvCount"), "skinned surface", tDegree + 1, numCrvs)); 431 432 //StackContext.enter(); 433 try { 434 435 // We have to create a list of compatible NURBS curves to create a NURBS surface. 436 GeomList<NurbsCurve> newCrvs = generateCompatibleCurves(crvs, tol); 437 438 // Extact the S knot vector (in common to all the compatible curves). 439 KnotVector sKnots = newCrvs.getFirst().getKnotVector(); 440 441 // Extract the control point network. 442 FastTable<FastTable<ControlPoint>> crvCParr = FastTable.newInstance(); 443 for (int i = 0; i < numCrvs; ++i) { 444 NurbsCurve crv = newCrvs.get(i); 445 crvCParr.add((FastTable<ControlPoint>)crv.getControlPoints()); 446 } 447 int numCPs = crvCParr.get(0).size(); 448 449 // Create a knot vector in the T direction by averaging the chord-type 450 // distance between the curve control points. 451 double[] d = ArrayFactory.DOUBLES_FACTORY.array(numCPs); 452 for (int i = 0; i < numCPs; ++i) { 453 d[i] = 0; 454 for (int k = 1; k < numCrvs; ++k) { 455 ControlPoint Qk = crvCParr.get(k).get(i); 456 ControlPoint Qkm1 = crvCParr.get(k - 1).get(i); 457 double distV = Qk.distance(Qkm1); 458 d[i] += distV; 459 } 460 } 461 double[] tParams = ArrayFactory.DOUBLES_FACTORY.array(numCrvs); 462 tParams[0] = 0; 463 for (int k = 1; k < numCrvs; ++k) { 464 tParams[k] = 0; 465 for (int i = 0; i < numCPs; ++i) { 466 ControlPoint Qk = crvCParr.get(k).get(i); 467 ControlPoint Qkm1 = crvCParr.get(k - 1).get(i); 468 double distV = Qk.distance(Qkm1); 469 tParams[k] += distV / d[i]; 470 } 471 tParams[k] /= numCPs; 472 tParams[k] += tParams[k - 1]; 473 } 474 tParams[numCrvs - 1] = 1; 475 476 KnotVector tKnots = CurveFactory.buildInterpKnotVector(tParams, numCrvs, tDegree); 477 478 // Convert to the surface representation. 479 FastTable<FastTable<ControlPoint>> newCPLst = FastTable.newInstance(); 480 FastTable<ControlPoint> pnts = FastTable.newInstance(); 481 for (int i = 0; i < numCPs; ++i) { 482 483 // Get the control points on each definining section. 484 pnts.clear(); 485 for (int k = 0; k < numCrvs; ++k) { 486 pnts.add(crvCParr.get(k).get(i)); 487 } 488 489 // Fit a curve to the control points in the T-direction. 490 BasicNurbsCurve iCrv = CurveFactory.fitControlPoints(tDegree, pnts, tParams, tKnots); 491 492 // Extract the control points of the i-curve. 493 FastTable<ControlPoint> newCPs = (FastTable<ControlPoint>)iCrv.getControlPoints(); 494 495 // Save off the new control points. 496 newCPLst.add(newCPs); 497 } 498 499 // Create the control point network for the surface. 500 ControlPointNet cpNet = ControlPointNet.valueOf(newCPLst).transpose(); 501 502 // Create the surface. 503 BasicNurbsSurface srf = BasicNurbsSurface.newInstance(cpNet, sKnots, tKnots); 504 return srf; //StackContext.outerCopy(srf); 505 506 } finally { 507 //StackContext.exit(); 508 } 509 } 510 511 /** 512 * Generate and return a list of NURBS surface compatible curves from the supplied 513 * list of curves. 514 */ 515 private static GeomList<NurbsCurve> generateCompatibleCurves(List<? extends Curve> crvs, Parameter<Length> tol) { 516 517 // Convert all the section curves to NURBS curves first. 518 int maxDeg = 0; 519 GeomList<NurbsCurve> lst = GeomList.newInstance(); 520 int numCrvs = crvs.size(); 521 for (int i = 0; i < numCrvs; ++i) { 522 Curve crv = crvs.get(i); 523 524 // Convert member curve to NURBS. 525 NurbsCurve nCrv = crv.toNurbs(tol); 526 lst.add(nCrv); 527 528 // Find the highest degree among the member curves. 529 if (nCrv.getDegree() > maxDeg) 530 maxDeg = nCrv.getDegree(); 531 } 532 533 // Raise all the section curves to the same degree. 534 GeomList<NurbsCurve> nCrvs = GeomList.newInstance(); 535 for (int i = 0; i < numCrvs; ++i) { 536 NurbsCurve crv = lst.get(i); 537 nCrvs.add(CurveUtils.elevateToDegree(crv, maxDeg)); 538 } 539 540 // Merge the knot vectors of all the section curves into a Union knot vector. 541 Float64Vector U = nCrvs.get(0).getKnotVector().getAll(); 542 for (int i = 1; i < numCrvs; ++i) { 543 Float64Vector Ui = nCrvs.get(i).getKnotVector().getAll(); 544 U = CurveUtils.knotVectorUnion(U, Ui); 545 } 546 547 // Merge the union knot vector with the knot vector of each section curve. 548 GeomList<NurbsCurve> compCrvs = GeomList.newInstance(); 549 for (int i = 0; i < numCrvs; ++i) { 550 NurbsCurve crv = nCrvs.get(i); 551 NurbsCurve ncrv = crv.mergeKnotVector(U); 552 compCrvs.add(ncrv); 553 } 554 555 // Finally, convert all the curves to the same units. 556 Unit<Length> refUnit = compCrvs.get(0).getUnit(); 557 GeomList<NurbsCurve> outCrvs = GeomList.newInstance(); 558 outCrvs.add(compCrvs.get(0)); 559 for (int i = 1; i < numCrvs; ++i) 560 outCrvs.add(compCrvs.get(i).to(refUnit)); 561 562 return outCrvs; 563 } 564 565 /** 566 * Create a Transfinite Interpolation (TFI) NURBS surface, or Coons patch, from the 567 * specified boundary curves. This creates a surface in NURBS form that uses the input 568 * set of curves as boundaries on all sides to within the specified tolerance. The 569 * degree in each parametric dimension will be determined by the highest degree of the 570 * input curves. If the input curves are all NURBS curves, then the "tol" parameter is 571 * ignored. Otherwise, "tol" is used to make a NURBS approximation of the input 572 * curves. It is assumed that all the input curves are oriented in an appropriate 573 * direction and order to create the desired TFI surface. The corner points are 574 * checked for coincidence. 575 * 576 * @param s0crv The curve that represents the s=0 boundary of the TFI surface. May not 577 * be null. 578 * @param t0crv The curve that represents the t=0 boundary of the TFI surface. May not 579 * be null. 580 * @param s1crv The curve that represents the s=1 boundary of the TFI surface. May not 581 * be null. 582 * @param t1crv The curve that represents the t=1 boundary of the TFI surface. May not 583 * be null. 584 * @param tol The greatest possible difference between the input list of curves and 585 * the NURBS surface returned. May not be null or zero. 586 * @return A BasicNurbsSurface that forms a TFI surface from the input curves. 587 */ 588 public static BasicNurbsSurface createTFISurface(Curve s0crv, Curve t0crv, Curve s1crv, Curve t1crv, 589 Parameter<Length> tol) throws IllegalArgumentException, ParameterException { 590 // Reference: Handbook of Grid Generation, Chapter 30, Section 30.3.8 591 592 //StackContext.enter(); 593 try { 594 // Check to make sure all the corner points are properly coincident. 595 Point p00 = s0crv.getRealPoint(0); 596 Point p1 = t0crv.getRealPoint(0); 597 if (p00.distance(p1).isGreaterThan(tol)) { 598 throw new IllegalArgumentException( 599 MessageFormat.format(RESOURCES.getString("s00t00NotCoincident"), 600 p00.distance(p1).toString())); 601 } 602 Point p10 = t0crv.getRealPoint(1); 603 p1 = s1crv.getRealPoint(0); 604 if (p10.distance(p1).isGreaterThan(tol)) { 605 throw new IllegalArgumentException( 606 MessageFormat.format(RESOURCES.getString("s10t01NotCoincident"), 607 p10.distance(p1).toString())); 608 } 609 Point p11 = s1crv.getRealPoint(1); 610 p1 = t1crv.getRealPoint(1); 611 if (p11.distance(p1).isGreaterThan(tol)) { 612 throw new IllegalArgumentException( 613 MessageFormat.format(RESOURCES.getString("s11t11NotCoincident"), 614 p11.distance(p1).toString())); 615 } 616 Point p01 = t1crv.getRealPoint(0); 617 p1 = s0crv.getRealPoint(1); 618 if (p01.distance(p1).isGreaterThan(tol)) { 619 throw new IllegalArgumentException( 620 MessageFormat.format(RESOURCES.getString("s01t10NotCoincident"), 621 p01.distance(p1).toString())); 622 } 623 624 // Create a list for the two curves in the S direction. 625 GeomList<Curve> sCrvs = GeomList.newInstance(); 626 sCrvs.add(t0crv, t1crv); 627 628 // Create a skinned intermediate surface in the S direction. 629 NurbsSurface Ps = createSkinnedSurface(sCrvs, 1, tol); 630 631 // Create a list for the two curves in the T direction. 632 GeomList<Curve> tCrvs = GeomList.newInstance(); 633 tCrvs.add(s0crv, s1crv); 634 635 // Create a skinned intermediate surface in the T direction. 636 NurbsSurface Pt = createSkinnedSurface(tCrvs, 1, tol).transpose(); 637 638 // Create a flat intermediate surface between the corner points. 639 PointString<Point> str1 = PointString.valueOf(p00, p10); 640 PointString<Point> str2 = PointString.valueOf(p01, p11); 641 PointArray<Point> arr = PointArray.valueOf(str1, str2); 642 NurbsSurface PsPt = fitPoints(1, 1, arr); 643 644 // Get the max degree in the s and t direcitons. 645 int maxDegS = Ps.getSDegree(); 646 if (PsPt.getSDegree() > maxDegS) 647 maxDegS = PsPt.getSDegree(); 648 if (Pt.getSDegree() > maxDegS) 649 maxDegS = Pt.getSDegree(); 650 int maxDegT = Ps.getTDegree(); 651 if (PsPt.getTDegree() > maxDegT) 652 maxDegT = PsPt.getTDegree(); 653 if (Pt.getTDegree() > maxDegT) 654 maxDegT = Pt.getTDegree(); 655 656 // Raise all the surfaces to the same degree in the S & T directions. 657 PsPt = SurfaceUtils.elevateToDegree(PsPt, maxDegS, maxDegT); 658 Ps = SurfaceUtils.elevateToDegree(Ps, maxDegS, maxDegT); 659 Pt = SurfaceUtils.elevateToDegree(Pt, maxDegS, maxDegT); 660 661 // Merge the knot vectors of all the surfaces into a Union knot vector in each direction. 662 Float64Vector Us = PsPt.getSKnotVector().getAll(); 663 Float64Vector Ui = Ps.getSKnotVector().getAll(); 664 Us = CurveUtils.knotVectorUnion(Us, Ui); 665 Ui = Pt.getSKnotVector().getAll(); 666 Us = CurveUtils.knotVectorUnion(Us, Ui); 667 668 Float64Vector Ut = PsPt.getTKnotVector().getAll(); 669 Ui = Ps.getTKnotVector().getAll(); 670 Ut = CurveUtils.knotVectorUnion(Ut, Ui); 671 Ui = Pt.getTKnotVector().getAll(); 672 Ut = CurveUtils.knotVectorUnion(Ut, Ui); 673 674 // Merge the union knot vectors with the knot vectors of each surface. 675 PsPt = PsPt.mergeSKnotVector(Us); 676 PsPt = PsPt.mergeTKnotVector(Ut); 677 Ps = Ps.mergeSKnotVector(Us); 678 Ps = Ps.mergeTKnotVector(Ut); 679 Pt = Pt.mergeSKnotVector(Us); 680 Pt = Pt.mergeTKnotVector(Ut); 681 682 // Get the control point networks for the intermediate surfaces. 683 ControlPointNet cpNet = Ps.getControlPoints(); 684 ControlPointNet cpNetPt = Pt.getControlPoints(); 685 ControlPointNet cpNetPsPt = PsPt.getControlPoints(); 686 687 // Construct a new control point network where: 688 // Points: P = Ps + Pt - PsPt 689 // Weights: w = Ws * Wt 690 FastTable<FastTable<ControlPoint>> cpNetArr = FastTable.newInstance(); 691 int rows = cpNet.getNumberOfRows(); 692 int cols = cpNet.getNumberOfColumns(); 693 for (int t = 0; t < cols; ++t) { 694 FastTable<ControlPoint> cpLst = FastTable.newInstance(); 695 cpNetArr.add(cpLst); 696 for (int s = 0; s < rows; ++s) { 697 ControlPoint cp = cpNet.get(s, t); 698 double wPs = cp.getWeight(); 699 ControlPoint cpPt = cpNetPt.get(s, t); 700 double wPt = cpPt.getWeight(); 701 ControlPoint cpPsPt = cpNetPsPt.get(s, t); 702 cp = cp.plus(cpPt); 703 cp = cp.minus(cpPsPt); 704 cp = cp.changeWeight(wPs * wPt); 705 cpLst.add(cp); 706 } 707 } 708 709 // Construct a new surface with the new control point network. 710 cpNet = ControlPointNet.valueOf(cpNetArr); 711 BasicNurbsSurface newSrf = BasicNurbsSurface.newInstance(cpNet, 712 Ps.getSKnotVector(), Ps.getTKnotVector()); 713 714 return newSrf; //StackContext.outerCopy(newSrf); 715 716 } finally { 717 //StackContext.exit(); 718 } 719 } 720 721 /** 722 * Create a {@link NurbsSurface} of the specified degrees that is fit to, 723 * interpolates, or passes through the specified array of geometry points in the input 724 * order. This is a global interpolation of the points, meaning that if one of the 725 * input points is moved, it will affect the entire surface, not just the local 726 * portion near the point. 727 * 728 * @param sDegree The degree in the s-direction (down the strings of points) of the 729 * NURBS surface to create (must be > 0 and < the number of 730 * points in each string). 731 * @param tDegree The degree in the t-direction (across the rows of strings) of the 732 * NURBS surface to create (must be > 0 and < the number of 733 * strings in the array). 734 * @param points The array of GeomPoint objects to pass the surface through: 735 * points[col idx][row idx] -> [t][s]. May not be null. 736 * @return A BasicNurbsSurface fit to the input array of points. 737 * @throws ParameterException if there is a problem with the parameterization 738 * resulting from fitting a curve through the supplied points. 739 * @see #approxPoints 740 */ 741 public static BasicNurbsSurface fitPoints(int sDegree, int tDegree, PointArray<? extends GeomPoint> points) 742 throws IllegalArgumentException, ParameterException { 743 int nt = points.size(); 744 if (nt < 1) 745 throw new IllegalArgumentException( 746 MessageFormat.format(RESOURCES.getString("zeroLengthListErr"), 747 "points")); 748 if (tDegree >= nt) 749 throw new IllegalArgumentException(RESOURCES.getString("numColsLEQDegErr")); 750 int ns = points.get(0).size(); 751 if (ns < 1) 752 throw new IllegalArgumentException( 753 MessageFormat.format(RESOURCES.getString("numGTZero"), 754 "number of rows", ns)); 755 if (sDegree >= ns) 756 throw new IllegalArgumentException(RESOURCES.getString("numRowsLEPDegErr")); 757 for (int i = 0; i < nt; ++i) { 758 if (points.get(i).size() != ns) 759 throw new IllegalArgumentException(RESOURCES.getString("numPointsInEachColErr")); 760 } 761 762 //StackContext.enter(); 763 try { 764 // Turn the points into an array of parameter values on the surface. 765 double[][] st = parameterizeArray(points); 766 767 // Check for bad parameterization. 768 CurveFactory.parameterizationCheck(ns, st[0]); 769 CurveFactory.parameterizationCheck(nt, st[1]); 770 771 // Generate the knot vectors. 772 KnotVector sKnots = CurveFactory.buildInterpKnotVector(st[0], ns, sDegree); 773 KnotVector tKnots = CurveFactory.buildInterpKnotVector(st[1], nt, tDegree); 774 775 // Generate the matrix of control points. 776 FastTable<FastTable<ControlPoint>> rv = FastTable.newInstance(); 777 for (int i = 0; i < ns; i++) { 778 PointString tmp = PointString.newInstance(); 779 for (int l = 0; l < nt; l++) { 780 tmp.add(points.get(l).get(i)); 781 } 782 783 NurbsCurve curve = CurveFactory.fitPoints(tDegree, tmp); 784 rv.add((FastTable<ControlPoint>)curve.getControlPoints()); 785 } 786 787 FastTable<FastTable<ControlPoint>> cp = FastTable.newInstance(); 788 for (int l = 0; l < nt; l++) { 789 PointString tmp = PointString.newInstance(); 790 for (int i = 0; i < ns; i++) { 791 tmp.add(rv.get(i).get(l).getPoint()); 792 } 793 794 NurbsCurve curve = CurveFactory.fitPoints(sDegree, tmp); 795 cp.add((FastTable<ControlPoint>)curve.getControlPoints()); 796 } 797 798 // Create the surface. 799 ControlPointNet cpNet = ControlPointNet.valueOf(cp); 800 BasicNurbsSurface srf = BasicNurbsSurface.newInstance(cpNet, sKnots, tKnots); 801 return srf; //StackContext.outerCopy(srf); 802 803 } finally { 804 //StackContext.exit(); 805 } 806 } 807 808 /** 809 * Turn an array of points into an array of parameter values on a surface through 810 * those points. The returned array was created using ArrayFactor.DOUBLES_FACTORY and 811 * can be recycled. 812 */ 813 private static double[][] parameterizeArray(PointArray<? extends GeomPoint> points) { 814 // This uses the Chord Length method which is described here: 815 // http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/INT-APP/PARA-chord-length.html 816 817 Unit<Length> refUnit = points.get(0).get(0).getUnit(); 818 int rows = points.get(0).size(); 819 int n = rows - 1; 820 int cols = points.size(); 821 int m = cols - 1; 822 double[] sk = newZeroedDoubleArray(rows); // sk = new double[rows]; 823 double[] tk = newZeroedDoubleArray(cols); // tk = new double[cols]; 824 825 StackContext.enter(); 826 try { 827 double[] cds = ArrayFactory.DOUBLES_FACTORY.array(rows * cols); // cds = new double[rows*cols]; 828 829 int num = cols; 830 sk[n] = 1; 831 for (int l = 0; l < cols; ++l) { 832 double total = 0; 833 834 for (int k = 1; k < rows; k++) { 835 int km1 = k - 1; 836 Parameter<Length> distance = points.get(l).get(k).distance(points.get(l).get(km1)); 837 double distV = distance.getValue(refUnit); 838 cds[k] = distV; 839 total += distV; 840 } 841 if (abs(total) <= EPSC) { 842 // The points form a degenerate curve of zero length. 843 844 // Just evenly space out parameters. 845 for (int k = 1; k < n; ++k) { 846 double s = (double)(k) / n; 847 sk[k] += s; 848 } 849 850 } else { 851 // The points form a curve of finite length. 852 double d = 0; 853 for (int k = 1; k < n; k++) { 854 d += cds[k]; 855 sk[k] += d / total; 856 } 857 } 858 } 859 if (num > 0) { 860 for (int k = 1; k < n; k++) { 861 sk[k] /= num; 862 } 863 } 864 865 num = rows; 866 tk[m] = 1; 867 for (int l = 0; l < rows; l++) { 868 double total = 0; 869 for (int k = 1; k < cols; k++) { 870 int km1 = k - 1; 871 Parameter<Length> distance = points.get(k).get(l).distance(points.get(km1).get(l)); 872 double distV = distance.getValue(refUnit); 873 cds[k] = distV; 874 total += distV; 875 } 876 if (abs(total) <= EPSC) { 877 // The points form a degenerate curve of zero length. 878 879 // Just evenly space out parameters. 880 for (int k = 1; k < m; ++k) { 881 double t = (double)(k) / n; 882 tk[k] += t; 883 } 884 885 } else { 886 // The points form a curve of finite length. 887 double d = 0; 888 for (int k = 1; k < m; k++) { 889 d += cds[k]; 890 tk[k] += d / total; 891 } 892 } 893 } 894 if (num > 0) { 895 for (int k = 1; k < m; k++) { 896 tk[k] /= num; 897 } 898 } 899 900 } finally { 901 StackContext.exit(); 902 } 903 904 // Create the output array of arrays. 905 double[][] res = DOUBLEDOUBLE_FACTORY.array(2); // res = new double[2][] 906 res[0] = sk; 907 res[1] = tk; 908 909 return res; 910 } 911 912 /** 913 * Create a {@link NurbsSurface} of the specified degrees that approximates, in a 914 * least squares sense, the specified array of geometry points in the input order. 915 * This is a global approximation of the points, meaning that if one of the input 916 * points is moved, it will affect the entire surface, not just the local portion near 917 * the point. 918 * 919 * @param sDegree The degree in the s-direction of the NURBS surface to create (must 920 * be > 0 and < the number of points in each string). 921 * @param tDegree The degree in the t-direction of the NURBS surface to create (must 922 * be > 0 and < the number of strings in the array). 923 * @param nS The number of surface control points to use in the parametric S 924 * direction (must be > 1 and < the number of points in each 925 * string). 926 * @param nT The number of surface control points to use in the parametric T 927 * direction (must be > 1 and < the number of strings in the 928 * array). 929 * @param points The array of GeomPoint objects approximate a surface through: 930 * points[col idx][row idx] -> [t][s]. May not be null. 931 * @return A BasicNurbsSurface approximating the input array of points. 932 * @throws ParameterException if there is a problem with the parameterization that 933 * results from approximating the supplied points. 934 * @see #fitPoints 935 */ 936 public static BasicNurbsSurface approxPoints(int sDegree, int tDegree, int nS, int nT, 937 PointArray<? extends GeomPoint> points) throws IllegalArgumentException, ParameterException { 938 int mt = points.size(); 939 if (mt < 1) 940 throw new IllegalArgumentException( 941 MessageFormat.format(RESOURCES.getString("zeroLengthListErr"), 942 "points")); 943 if (tDegree <= 0) 944 throw new IllegalArgumentException( 945 MessageFormat.format(RESOURCES.getString("numGTZero"), 946 "degree in the T-direction", tDegree)); 947 if (tDegree >= mt) 948 throw new IllegalArgumentException(RESOURCES.getString("numColsLEQDegErr")); 949 int ms = points.get(0).size(); 950 if (ms < 1) 951 throw new IllegalArgumentException( 952 MessageFormat.format(RESOURCES.getString("numGTZero"), 953 "number of rows", ms)); 954 if (sDegree <= 0) 955 throw new IllegalArgumentException( 956 MessageFormat.format(RESOURCES.getString("numGTZero"), 957 "degree in the S-direction", sDegree)); 958 if (sDegree >= ms) 959 throw new IllegalArgumentException(RESOURCES.getString("numRowsLEPDegErr")); 960 if (nS <= 1) 961 throw new IllegalArgumentException( 962 MessageFormat.format(RESOURCES.getString("numPointsLEOneErr"), "S")); 963 if (nT <= 1) 964 throw new IllegalArgumentException( 965 MessageFormat.format(RESOURCES.getString("numPointsLEOneErr"), "T")); 966 if (nS >= ms) 967 throw new IllegalArgumentException( 968 MessageFormat.format(RESOURCES.getString("numDPLEnumCPErr"), "S")); 969 if (nT >= mt) 970 throw new IllegalArgumentException( 971 MessageFormat.format(RESOURCES.getString("numDPLEnumCPErr"), "T")); 972 for (int i = 0; i < mt; ++i) { 973 if (points.get(i).size() != ms) 974 throw new IllegalArgumentException(RESOURCES.getString("numPointsInEachColErr")); 975 } 976 977 StackContext.enter(); 978 try { 979 // Turn the points into an array of parameter values on the surface. 980 double[][] st = parameterizeArray(points); 981 982 // Check the parameterizations for correctness. 983 CurveFactory.parameterizationCheck(nS, st[0]); 984 CurveFactory.parameterizationCheck(nT, st[1]); 985 986 // Generate the knot vectors. 987 KnotVector sKnots = CurveFactory.buildApproxKnotVector(st[0], ms, nS, sDegree); 988 KnotVector tKnots = CurveFactory.buildApproxKnotVector(st[1], mt, nT, tDegree); 989 990 // Generate the matrix of control points. 991 FastTable<FastTable<ControlPoint>> rv = FastTable.newInstance(); 992 for (int i = 0; i < ms; i++) { 993 PointString tmp = PointString.newInstance(); 994 for (int l = 0; l < mt; l++) { 995 tmp.add(points.get(l).get(i)); 996 } 997 998 NurbsCurve curve = CurveFactory.approxPoints(tDegree, nT, tmp); 999 rv.add((FastTable<ControlPoint>)curve.getControlPoints()); 1000 } 1001 1002 FastTable<FastTable<ControlPoint>> cp = FastTable.newInstance(); 1003 for (int l = 0; l < nT; l++) { 1004 PointString tmp = PointString.newInstance(); 1005 for (int i = 0; i < ms; i++) { 1006 List<ControlPoint> cpLst = rv.get(i); 1007 ControlPoint ctrlp = cpLst.get(l); 1008 tmp.add(ctrlp.getPoint()); 1009 } 1010 1011 NurbsCurve curve = CurveFactory.approxPoints(sDegree, nS, tmp); 1012 cp.add((FastTable<ControlPoint>)curve.getControlPoints()); 1013 } 1014 1015 // Create and return the surface. 1016 ControlPointNet cpNet = ControlPointNet.valueOf(cp); 1017 BasicNurbsSurface srf = BasicNurbsSurface.newInstance(cpNet, sKnots, tKnots); 1018 return StackContext.outerCopy(srf); 1019 1020 } finally { 1021 StackContext.exit(); 1022 } 1023 } 1024 1025 /** 1026 * Return a new double array with the values zeroed. The returned instance is created 1027 * using ArrayFactory.DOUBLES_FACTORY and may be recycled. 1028 */ 1029 private static double[] newZeroedDoubleArray(int size) { 1030 double[] arr = ArrayFactory.DOUBLES_FACTORY.array(size); 1031 for (int i = 0; i < size; ++i) 1032 arr[i] = 0; 1033 return arr; 1034 } 1035 1036 private static final ArrayFactory<double[][]> DOUBLEDOUBLE_FACTORY = new ArrayFactory<double[][]>() { 1037 @Override 1038 protected double[][] create(int size) { 1039 return new double[size][]; 1040 } 1041 }; 1042 1043}