001/** 002 * AbstractSurface -- Represents the implementation in common to all surfaces in nD space. 003 * 004 * Copyright (C) 2010-2025, Joseph A. Huwaldt. All rights reserved. 005 * 006 * This library is free software; you can redistribute it and/or modify it under the terms 007 * of the GNU Lesser General Public License as published by the Free Software Foundation; 008 * either version 2.1 of the License, or (at your option) any later version. 009 * 010 * This library is distributed in the hope that it will be useful, but WITHOUT ANY 011 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A 012 * PARTICULAR PURPOSE. See the GNU Lesser 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 */ 018package geomss.geom; 019 020import geomss.geom.nurbs.BasicNurbsCurve; 021import geomss.geom.nurbs.BasicNurbsSurface; 022import geomss.geom.nurbs.CurveFactory; 023import geomss.geom.nurbs.CurveUtils; 024import jahuwaldt.js.param.Parameter; 025import jahuwaldt.tools.math.*; 026import static jahuwaldt.tools.math.MathTools.isApproxEqual; 027import java.text.MessageFormat; 028import java.util.Collections; 029import java.util.Comparator; 030import java.util.List; 031import static java.util.Objects.isNull; 032import static java.util.Objects.nonNull; 033import static java.util.Objects.requireNonNull; 034import java.util.logging.Level; 035import java.util.logging.Logger; 036import javax.measure.quantity.Area; 037import javax.measure.quantity.Dimensionless; 038import javax.measure.quantity.Length; 039import javax.measure.quantity.Volume; 040import javax.measure.unit.SI; 041import javax.measure.unit.Unit; 042import javolution.context.*; 043import static javolution.lang.MathLib.abs; 044import static javolution.lang.MathLib.round; 045import static javolution.lang.MathLib.sqrt; 046import javolution.util.FastTable; 047import org.apache.commons.math3.linear.SingularMatrixException; 048 049/** 050 * The interface and implementation in common to all surfaces. 051 * 052 * <p> Modified by: Joseph A. Huwaldt </p> 053 * 054 * @author Joseph A. Huwaldt, Date: June 14, 2010 055 * @version February 14, 2025 056 * 057 * @param <T> The sub-type of this AbstractSurface element. 058 */ 059@SuppressWarnings({"serial", "CloneableImplementsClone"}) 060public abstract class AbstractSurface<T extends AbstractSurface> extends AbstractGeomElement<T> implements Surface<T> { 061 062 /** 063 * Reference 1: Michael E. Mortenson, "Geometric Modeling, Third Edition", 064 * ISBN:9780831132989, 2008. 065 */ 066 067 /** 068 * Generic/default tolerance for use in root finders, etc. 069 */ 070 public static final double GTOL = 1e-6; 071 072 // Square-Root of EPS. 073 private static final double SQRT_EPS = Parameter.SQRT_EPS; 074 075 // A Parameter version of GTOL. 076 private static final Parameter<Length> GTOLP; 077 078 static { 079 ImmortalContext.enter(); 080 try { 081 GTOLP = Parameter.valueOf(GTOL / 10, SI.METER); 082 } finally { 083 ImmortalContext.exit(); 084 } 085 } 086 087 /** 088 * Tolerance allowed on parametric spacing being == 0 or == 1. 089 */ 090 protected static final double TOL_ST = Parameter.SQRT_EPS; 091 092 /** 093 * Returns the number of parametric dimensions of the geometry element. This 094 * implementation always returns 2. 095 */ 096 @Override 097 public int getParDimension() { 098 return 2; 099 } 100 101 /** 102 * Checks the input parameters for validity (in the range 0 to 1 for example) and 103 * throws an exception if it is not valid. 104 * 105 * @param s 1st parametric dimension distance (0.0 to 1.0 inclusive). 106 * @param t 2nd parametric dimension distance (0.0 to 1.0 inclusive). 107 * @throws IllegalArgumentException if there is any problem with the parameter values. 108 */ 109 protected void validateParameter(double s, double t) { 110 if (s < -TOL_ST || s > 1 + TOL_ST) 111 throw new IllegalArgumentException( 112 MessageFormat.format(RESOURCES.getString("srfInvalidSValue"), s)); 113 if (t < -TOL_ST || t > 1 + TOL_ST) 114 throw new IllegalArgumentException( 115 MessageFormat.format(RESOURCES.getString("srfInvalidTValue"), t)); 116 } 117 118 /** 119 * Checks the input parameters for validity (in the range 0 to 1 for example) and 120 * throws an exception if it is not valid. 121 * 122 * @param st The parametric position. Must be a 2-dimensional point with each value in 123 * the range 0 to 1.0. Units are ignored. May not be null. 124 * @throws IllegalArgumentException if there is any problem with the parameter values. 125 */ 126 protected void validateParameter(GeomPoint st) { 127 if (isNull(st) || st.getPhyDimension() != 2) 128 throw new IllegalArgumentException(MessageFormat.format( 129 RESOURCES.getString("invalidParamDim"), "Surface", 2)); 130 validateParameter(st.getValue(0), st.getValue(1)); 131 } 132 133 /** 134 * If the input parameter is within tol of 0 or 1, or outside of that bounds, then the 135 * output value will be set to exactly 0 or 1. Otherwise the input value will be 136 * returned. 137 * 138 * @param s The parameter to check for approximate end values. 139 * @return The input value or 0 and 1 if the input value is within tolerance of the 140 * ends. 141 */ 142 protected static final double roundParNearEnds(double s) { 143 return (s < TOL_ST ? 0. : (s > 1. - TOL_ST ? 1. : s)); 144 } 145 146 /** 147 * Return true if the input parameter is within tolerance of the start of the 148 * parameter range (s ≤ tol). 149 * 150 * @param s The parameter to check. 151 * @param tol The tolerance to use for the check. 152 * @return true if the parameter is within tolerance of zero. 153 */ 154 protected static final boolean parNearStart(double s, double tol) { 155 return s <= tol; 156 } 157 158 /** 159 * Return true if the input parameter is within tolerance of the end of the parameter 160 * range (s ≥ 1.0 - tol). 161 * 162 * @param s The parameter to check. 163 * @param tol The tolerance to use for the check. 164 * @return true if the parameter is within tolerance of 1.0. 165 */ 166 protected static final boolean parNearEnd(double s, double tol) { 167 return s >= 1.0 - tol; 168 } 169 170 /** 171 * Return true if the input parameter is within tolerance of either the end of the 172 * parameter range or outside that bounds (s ≤ tol || s ≥ 1.0 - tol). 173 * 174 * @param s The parameter to check. 175 * @param tol The tolerance to use for the check. 176 * @return true if the parameter is within tolerance of 0.0 or 1.0. 177 */ 178 protected static final boolean parNearEnds(double s, double tol) { 179 return s <= tol || s >= 1.0 - tol; 180 } 181 182 /** 183 * Always throws UnsupportedOperationException as this method is not supported for 184 * this object type! Subclasses that do support transpose() should override this 185 * method. 186 * 187 * @return A new surface, identical to this one, but with the transpose of the 188 * parameterization. 189 * @throws UnsupportedOperationException as this method is not supported on this 190 * object type. 191 */ 192 @Override 193 public T transpose() { 194 throw new UnsupportedOperationException(RESOURCES.getString("srfTransposeUnsupported")); 195 } 196 197 /** 198 * Split this surface at the specified parametric S-position returning a list 199 * containing two new surfaces (a lower surface with smaller S-parametric positions 200 * than "s" and an upper surface with larger S-parametric positions). 201 * 202 * @param st The S-parametric position where this surface should be split (must not be 203 * 0 or 1!). Must be a 2-dimensional point with each value in the range 0 to 204 * 1, only the 1st value is used. Units are ignored. May not be null. 205 * @return A list containing two surfaces: 0 == the lower surface, 1 == the upper 206 * surface. 207 */ 208 @Override 209 public GeomList<T> splitAtS(GeomPoint st) { 210 return splitAtS(st.getValue(0)); 211 } 212 213 /** 214 * Split this surface at the specified parametric T-position returning a list 215 * containing two new surfaces (a lower surface with smaller T-parametric positions 216 * than "t" and an upper surface with larger T-parametric positions). 217 * 218 * @param st The T-parametric position where this surface should be split (must not be 219 * 0 or 1!). Must be a 2-dimensional point with each value in the range 0 to 220 * 1, only the 2nd value is used. Units are ignored. May not be null. 221 * @return A list containing two surfaces: 0 == the lower surface, 1 == the upper 222 * surface. 223 */ 224 @Override 225 public GeomList<T> splitAtT(GeomPoint st) { 226 return splitAtT(st.getValue(1)); 227 } 228 229 /** 230 * Return a subrange point on the surface for the given parametric position on the 231 * surface. 232 * 233 * @param s 1st parametric dimension distance to calculate a point for (0.0 to 1.0 234 * inclusive). 235 * @param t 2nd parametric dimension distance to calculate a point for (0.0 to 1.0 236 * inclusive). 237 * @return The subrange point on the surface at the specified parameter values. 238 * @throws IllegalArgumentException if there is any problem with the parameter values. 239 */ 240 @Override 241 public SubrangePoint getPoint(double s, double t) { 242 validateParameter(s, t); 243 s = roundParNearEnds(s); 244 t = roundParNearEnds(t); 245 return SubrangePoint.newInstance(this, Point.valueOf(s, t)); 246 } 247 248 /** 249 * Return a subrange point on the surface for the given parametric position on the 250 * surface. 251 * 252 * @param st The parametric position to calculate a point for. Must be a 2-dimensional 253 * point with each value in the range 0 to 1.0. Units are ignored. 254 * @return The subrange point on the surface at the specified parameter values. 255 * @throws IllegalArgumentException if there is any problem with the parameter values. 256 */ 257 @Override 258 public SubrangePoint getPoint(GeomPoint st) { 259 validateParameter(st); 260 return SubrangePoint.newInstance(this, st); 261 } 262 263 /** 264 * Calculate a point on the surface for the given parametric position on the surface. 265 * 266 * @param st The parametric position to calculate a point for. Must be a 2-dimensional 267 * point with each value in the range 0 to 1.0. Units are ignored. 268 * @return The calculated point on the surface at the specified parameter values. 269 * @throws IllegalArgumentException if there is any problem with the parameter values. 270 */ 271 @Override 272 public Point getRealPoint(GeomPoint st) { 273 return getRealPoint(st.getValue(0), st.getValue(1)); 274 } 275 276 /** 277 * Return a subrange curve on the surface for the given parametric position curve. 278 * 279 * @param pcrv A curve in parametric space indicating the parametric positions on the 280 * the surface represented by the subrange curve. Must be 2D with values 281 * between 0 and 1. 282 * @return The subrange curve on the surface at the specified parametric curve 283 * positions. 284 */ 285 @Override 286 public SubrangeCurve getCurve(Curve pcrv) { 287 return SubrangeCurve.newInstance(this, requireNonNull(pcrv)); 288 } 289 290 /** 291 * Return a subrange curve at a constant parametric s-value. 292 * 293 * @param s The parametric s-position to extract a subrange curve. 294 * @return The subrange curve on the surface at the specified s-value. 295 * @throws IllegalArgumentException if there is any problem with the parameter value. 296 */ 297 @Override 298 public SubrangeCurve getSCurve(double s) { 299 validateParameter(s, 0); 300 s = roundParNearEnds(s); 301 PointString<Point> pstr = PointString.valueOf(Point.valueOf(s, 0), Point.valueOf(s, 1)); 302 Curve pcrv = CurveFactory.fitPoints(1, pstr); 303 PointString.recycle(pstr); 304 return getCurve(pcrv); 305 } 306 307 /** 308 * Return a subrange curve at a constant parametric t-value. 309 * 310 * @param t The parametric t-position to extract a subrange curve. 311 * @return The subrange curve on the surface at the specified t-value. 312 * @throws IllegalArgumentException if there is any problem with the parameter value. 313 */ 314 @Override 315 public SubrangeCurve getTCurve(double t) { 316 validateParameter(0, t); 317 t = roundParNearEnds(t); 318 PointString<Point> pstr = PointString.valueOf(Point.valueOf(0, t), Point.valueOf(1, t)); 319 Curve pcrv = CurveFactory.fitPoints(1, pstr); 320 PointString.recycle(pstr); 321 return getCurve(pcrv); 322 } 323 324 /** 325 * Calculate all the derivatives from <code>0</code> to <code>grade</code> with 326 * respect to parametric position(s) on a parametric object for the given parametric 327 * position on the object, <code>d^{grade}p(s)/d^{grade}s</code>. 328 * <p> 329 * Example:<br> 330 * 1st derivative (grade = 1), this returns <code>[p(s), dp(s)/ds]</code>;<br> 331 * 2nd derivative (grade = 2), this returns <code>[p(s), dp(s)/ds, d^2p(s)/d^2s]</code>; etc. 332 * </p> 333 * 334 * @param st The parametric position to calculate the derivatives for. Must match 335 * the parametric dimension of this parametric surface and have each 336 * value in the range 0 to 1.0. Units are ignored. 337 * @param grade The maximum grade to calculate the derivatives for (1=1st derivative, 338 * 2=2nd derivative, etc) 339 * @return A list of lists of derivatives (one list for each parametric dimension). 340 * Each list contains derivatives up to the specified grade at the specified 341 * parametric position. Example: for a surface list element 0 is the array of 342 * derivatives in the 1st parametric dimension (s), then 2nd list element is 343 * the array of derivatives in the 2nd parametric dimension (t). 344 * @throws IllegalArgumentException if the grade is < 0 or the parameter values are 345 * invalid. 346 */ 347 @Override 348 public List<List<Vector<Length>>> getDerivatives(GeomPoint st, int grade) { 349 validateParameter(st); 350 351 List<List<Vector<Length>>> list = FastTable.newInstance(); 352 list.add(getSDerivatives(st.getValue(0), st.getValue(1), grade, true)); 353 list.add(getTDerivatives(st.getValue(0), st.getValue(1), grade, true)); 354 355 return list; 356 } 357 358 /** 359 * Calculate all the derivatives from <code>0</code> to <code>grade</code> with 360 * respect to parametric s-position on the surface for the given parametric position 361 * on the surface, <code>d^{grade}p(s,t)/d^{grade}s</code>. 362 * <p> 363 * Example:<br> 364 * 1st derivative (grade = 1), this returns <code>[p(s,t), dp(s,t)/ds]</code>;<br> 365 * 2nd derivative (grade = 2), this returns <code>[p(s,t), dp(s,t)/ds, d^2p(s,t)/d^2s]</code>; etc. 366 * </p> 367 * 368 * @param st The parametric position to calculate the derivatives for. Must be a 369 * 2-dimensional point with each value in the range 0 to 1.0. Units are 370 * ignored. 371 * @param grade The maximum grade to calculate the derivatives for (1=1st derivative, 372 * 2=2nd derivative, etc) 373 * @return A list of s-derivatives up to the specified grade of the surface at the 374 * specified parametric position. 375 * @throws IllegalArgumentException if the grade is < 0 or the parameter values are 376 * invalid. 377 */ 378 @Override 379 public List<Vector<Length>> getSDerivatives(GeomPoint st, int grade) { 380 validateParameter(st); 381 return getSDerivatives(st.getValue(0), st.getValue(1), grade, true); 382 } 383 384 /** 385 * Calculate all the derivatives from <code>0</code> to <code>grade</code> with 386 * respect to parametric s-position on the surface for the given parametric position 387 * on the surface, <code>d^{grade}p(s,t)/d^{grade}s</code>. 388 * <p> 389 * Example:<br> 390 * 1st derivative (grade = 1), this returns <code>[p(s,t), dp(s,t)/ds]</code>;<br> 391 * 2nd derivative (grade = 2), this returns <code>[p(s,t), dp(s,t)/ds, d^2p(s,t)/d^2s]</code>; etc. 392 * </p> 393 * 394 * @param s 1st parametric dimension distance to calculate derivative for (0.0 to 395 * 1.0 inclusive). 396 * @param t 2nd parametric dimension distance to calculate derivative for (0.0 to 397 * 1.0 inclusive). 398 * @param grade The maximum grade to calculate the u-derivatives for (1=1st 399 * derivative, 2=2nd derivative, etc) 400 * @return A list of s-derivatives up to the specified grade of the surface at the 401 * specified parametric position. 402 * @throws IllegalArgumentException if the grade is < 0 or the parameter values are 403 * invalid. 404 */ 405 @Override 406 public List<Vector<Length>> getSDerivatives(double s, double t, int grade) { 407 return getSDerivatives(s, t, grade, true); 408 } 409 410 /** 411 * Calculate all the derivatives from <code>0</code> to <code>grade</code> with 412 * respect to parametric s-position on the surface for the given parametric position 413 * on the surface, <code>d^{grade}p(s,t)/d^{grade}s</code>. 414 * <p> 415 * Example:<br> 416 * 1st derivative (grade = 1), this returns <code>[p(s,t), dp(s,t)/ds]</code>;<br> 417 * 2nd derivative (grade = 2), this returns <code>[p(s,t), dp(s,t)/ds, d^2p(s,t)/d^2s]</code>; etc. 418 * </p> 419 * 420 * @param s 1st parametric dimension distance to calculate derivative for (0.0 to 421 * 1.0 inclusive). 422 * @param t 2nd parametric dimension distance to calculate derivative for (0.0 to 423 * 1.0 inclusive). 424 * @param grade The maximum grade to calculate the u-derivatives for (1=1st 425 * derivative, 2=2nd derivative, etc) 426 * @param scaled Pass <code>true</code> for properly scaled derivatives or 427 * <code>false</code> if the magnitude of the derivative vector is not 428 * required -- this sometimes results in faster calculation times. 429 * @return A list of s-derivatives up to the specified grade of the surface at the 430 * specified parametric position. 431 * @throws IllegalArgumentException if the grade is < 0 or the parameter values are 432 * invalid. 433 */ 434 public abstract List<Vector<Length>> getSDerivatives(double s, double t, int grade, boolean scaled); 435 436 /** 437 * Calculate all the derivatives from <code>0</code> to <code>grade</code> with 438 * respect to parametric t-position on the surface for the given parametric position 439 * on the surface, <code>d^{grade}p(s,t)/d^{grade}t</code>. 440 * <p> 441 * Example:<br> 442 * 1st derivative (grade = 1), this returns <code>[p(s,t), dp(s,t)/dt]</code>;<br> 443 * 2nd derivative (grade = 2), this returns <code>[p(s,t), dp(s,t)/dt, d^2p(s,t)/d^2t]</code>; etc. 444 * </p> 445 * 446 * @param st The parametric position to calculate the derivatives for. Must be a 447 * 2-dimensional point with each value in the range 0 to 1.0. Units are 448 * ignored. 449 * @param grade The maximum grade to calculate the derivatives for (1=1st derivative, 450 * 2=2nd derivative, etc) 451 * @return A list of t-derivatives up to the specified grade of the surface at the 452 * specified parametric position. 453 * @throws IllegalArgumentException if the grade is < 0 or the parameter values are 454 * invalid. 455 */ 456 @Override 457 public List<Vector<Length>> getTDerivatives(GeomPoint st, int grade) { 458 validateParameter(st); 459 return getTDerivatives(st.getValue(0), st.getValue(1), grade, true); 460 } 461 462 /** 463 * Calculate all the derivatives from <code>0</code> to <code>grade</code> with 464 * respect to parametric t-position on the surface for the given parametric position 465 * on the surface, <code>d^{grade}p(s,t)/d^{grade}t</code>. 466 * <p> 467 * Example:<br> 468 * 1st derivative (grade = 1), this returns <code>[p(s,t), dp(s,t)/dt]</code>;<br> 469 * 2nd derivative (grade = 2), this returns <code>[p(s,t), dp(s,t)/dt, d^2p(s,t)/d^2t]</code>; etc. 470 * </p> 471 * 472 * @param s 1st parametric dimension distance to calculate derivative for (0.0 to 473 * 1.0 inclusive). 474 * @param t 2nd parametric dimension distance to calculate derivative for (0.0 to 475 * 1.0 inclusive). 476 * @param grade The maximum grade to calculate the v-derivatives for (1=1st 477 * derivative, 2=2nd derivative, etc) 478 * @return A list of t-derivatives up to the specified grade of the surface at the 479 * specified parametric position. 480 * @throws IllegalArgumentException if the grade is < 0 or the parameter values are 481 * invalid. 482 */ 483 @Override 484 public List<Vector<Length>> getTDerivatives(double s, double t, int grade) { 485 return getTDerivatives(s, t, grade, true); 486 } 487 488 /** 489 * Calculate all the derivatives from <code>0</code> to <code>grade</code> with 490 * respect to parametric t-position on the surface for the given parametric position 491 * on the surface, <code>d^{grade}p(s,t)/d^{grade}t</code>. 492 * <p> 493 * Example:<br> 494 * 1st derivative (grade = 1), this returns <code>[p(s,t), dp(s,t)/dt]</code>;<br> 495 * 2nd derivative (grade = 2), this returns <code>[p(s,t), dp(s,t)/dt, d^2p(s,t)/d^2t]</code>; etc. 496 * </p> 497 * 498 * @param s 1st parametric dimension distance to calculate derivative for (0.0 to 499 * 1.0 inclusive). 500 * @param t 2nd parametric dimension distance to calculate derivative for (0.0 to 501 * 1.0 inclusive). 502 * @param grade The maximum grade to calculate the v-derivatives for (1=1st 503 * derivative, 2=2nd derivative, etc) 504 * @param scaled Pass <code>true</code> for properly scaled derivatives or 505 * <code>false</code> if the magnitude of the derivative vector is not 506 * required -- this sometimes results in faster calculation times. 507 * @return A list of t-derivatives up to the specified grade of the surface at the 508 * specified parametric position. 509 * @throws IllegalArgumentException if the grade is < 0 or the parameter values are 510 * invalid. 511 */ 512 public abstract List<Vector<Length>> getTDerivatives(double s, double t, int grade, boolean scaled); 513 514 /** 515 * Calculate a derivative with respect to parametric s-distance of the given grade on 516 * the surface for the given parametric position on the surface, 517 * <code>d^{grade}p(s,t)/d^{grade}s</code>. 518 * <p> 519 * Example:<br> 520 * 1st derivative (grade = 1), this returns <code>dp(s,t)/ds</code>;<br> 521 * 2nd derivative (grade = 2), this returns <code>d^2p(s,t)/d^2s</code>; etc. 522 * </p> 523 * 524 * @param st The parametric position to calculate the derivative for. Must be a 525 * 2-dimensional point with each value in the range 0 to 1.0. Units are 526 * ignored. 527 * @param grade The grade to calculate the derivative for (1=1st derivative, 2=2nd 528 * derivative, etc) 529 * @return The specified s-derivative of the surface at the specified parametric 530 * position. 531 * @throws IllegalArgumentException if the grade is < 0 or the parameter values are 532 * invalid. 533 */ 534 @Override 535 public Vector<Length> getSDerivative(GeomPoint st, int grade) { 536 validateParameter(st); 537 return getSDerivative(st.getValue(0), st.getValue(1), grade, true); 538 } 539 540 /** 541 * Calculate a derivative with respect to parametric s-distance of the given grade on 542 * the surface for the given parametric position on the surface, 543 * <code>d^{grade}p(s,t)/d^{grade}s</code>. 544 * <p> 545 * Example:<br> 546 * 1st derivative (grade = 1), this returns <code>dp(s,t)/ds</code>;<br> 547 * 2nd derivative (grade = 2), this returns <code>d^2p(s,t)/d^2s</code>; etc. 548 * </p> 549 * 550 * @param s 1st parametric dimension distance to calculate derivative for (0.0 to 551 * 1.0 inclusive). 552 * @param t 2nd parametric dimension distance to calculate derivative for (0.0 to 553 * 1.0 inclusive). 554 * @param grade The grade to calculate the derivative for (1=1st derivative, 2=2nd 555 * derivative, etc) 556 * @return The specified s-derivative of the surface at the specified parametric 557 * position. 558 * @throws IllegalArgumentException if the grade is < 0 or the parameter values are 559 * invalid. 560 */ 561 @Override 562 public Vector<Length> getSDerivative(double s, double t, int grade) { 563 return getSDerivative(s, t, grade, true); 564 } 565 566 /** 567 * Calculate a derivative with respect to parametric s-distance of the given grade on 568 * the surface for the given parametric position on the surface, 569 * <code>d^{grade}p(s,t)/d^{grade}s</code>. 570 * <p> 571 * Example:<br> 572 * 1st derivative (grade = 1), this returns <code>dp(s,t)/ds</code>;<br> 573 * 2nd derivative (grade = 2), this returns <code>d^2p(s,t)/d^2s</code>; etc. 574 * </p> 575 * 576 * @param s 1st parametric dimension distance to calculate derivative for (0.0 to 577 * 1.0 inclusive). 578 * @param t 2nd parametric dimension distance to calculate derivative for (0.0 to 579 * 1.0 inclusive). 580 * @param grade The grade to calculate the derivative for (1=1st derivative, 2=2nd 581 * derivative, etc) 582 * @param scaled Pass <code>true</code> for properly scaled derivatives or 583 * <code>false</code> if the magnitude of the derivative vector is not 584 * required -- this sometimes results in faster calculation times. 585 * @return The specified s-derivative of the surface at the specified parametric 586 * position. 587 * @throws IllegalArgumentException if the grade is < 0 or the parameter values are 588 * invalid. 589 */ 590 protected Vector<Length> getSDerivative(double s, double t, int grade, boolean scaled) { 591 List<Vector<Length>> ders = getSDerivatives(s, t, grade, scaled); 592 Vector<Length> derivative = ders.get(grade); 593 FastTable.recycle((FastTable)ders); 594 return derivative; 595 } 596 597 /** 598 * Calculate a derivative with respect to parametric t-distance of the given grade on 599 * the surface for the given parametric position on the surface, 600 * <code>d^{grade}p(s,t)/d^{grade}t</code>. 601 * <p> 602 * Example:<br> 603 * 1st derivative (grade = 1), this returns <code>dp(s,t)/dt</code>;<br> 604 * 2nd derivative (grade = 2), this returns <code>d^2p(s,t)/d^2t</code>; etc. 605 * </p> 606 * 607 * @param st The parametric position to calculate the derivative for. Must be a 608 * 2-dimensional point with each value in the range 0 to 1.0. Units are 609 * ignored. 610 * @param grade The grade to calculate the derivative for (1=1st derivative, 2=2nd 611 * derivative, etc) 612 * @return The specified t-derivative of the surface at the specified parametric 613 * position. 614 * @throws IllegalArgumentException if the grade is < 0 or the parameter values are 615 * invalid. 616 */ 617 @Override 618 public Vector<Length> getTDerivative(GeomPoint st, int grade) { 619 validateParameter(st); 620 return getTDerivative(st.getValue(0), st.getValue(1), grade, true); 621 } 622 623 /** 624 * Calculate a derivative with respect to parametric t-distance of the given grade on 625 * the surface for the given parametric position on the surface, 626 * <code>d^{grade}p(s,t)/d^{grade}t</code>. 627 * <p> 628 * Example:<br> 629 * 1st derivative (grade = 1), this returns <code>dp(s,t)/dt</code>;<br> 630 * 2nd derivative (grade = 2), this returns <code>d^2p(s,t)/d^2t</code>; etc. 631 * </p> 632 * 633 * @param s 1st parametric dimension distance to calculate derivative for (0.0 to 634 * 1.0 inclusive). 635 * @param t 2nd parametric dimension distance to calculate derivative for (0.0 to 636 * 1.0 inclusive). 637 * @param grade The grade to calculate the derivative for (1=1st derivative, 2=2nd 638 * derivative, etc) 639 * @return The specified t-derivative of the surface at the specified parametric 640 * position. 641 * @throws IllegalArgumentException if the grade is < 0 or the parameter values are 642 * invalid. 643 */ 644 @Override 645 public Vector<Length> getTDerivative(double s, double t, int grade) { 646 return getTDerivative(s, t, grade, true); 647 } 648 649 /** 650 * Calculate a derivative with respect to parametric t-distance of the given grade on 651 * the surface for the given parametric position on the surface, 652 * <code>d^{grade}p(s,t)/d^{grade}t</code>. 653 * <p> 654 * Example:<br> 655 * 1st derivative (grade = 1), this returns <code>dp(s,t)/dt</code>;<br> 656 * 2nd derivative (grade = 2), this returns <code>d^2p(s,t)/d^2t</code>; etc. 657 * </p> 658 * 659 * @param s 1st parametric dimension distance to calculate derivative for (0.0 to 660 * 1.0 inclusive). 661 * @param t 2nd parametric dimension distance to calculate derivative for (0.0 to 662 * 1.0 inclusive). 663 * @param grade The grade to calculate the derivative for (1=1st derivative, 2=2nd 664 * derivative, etc) 665 * @param scaled Pass <code>true</code> for properly scaled derivatives or 666 * <code>false</code> if the magnitude of the derivative vector is not 667 * required -- this sometimes results in faster calculation times. 668 * @return The specified t-derivative of the surface at the specified parametric 669 * position. 670 * @throws IllegalArgumentException if the grade is < 0 or the parameter values are 671 * invalid. 672 */ 673 protected Vector<Length> getTDerivative(double s, double t, int grade, boolean scaled) { 674 List<Vector<Length>> ders = getTDerivatives(s, t, grade, scaled); 675 Vector<Length> derivative = ders.get(grade); 676 FastTable.recycle((FastTable)ders); 677 return derivative; 678 } 679 680 /** 681 * Calculate the twist vector (d^2P/(ds*dt) = d(dP/ds)/dt) for this surface at the 682 * specified position on this surface. 683 * 684 * @param st The parametric position to calculate the twist vector for. Must be a 685 * 2-dimensional point with each value in the range 0 to 1.0. Units are 686 * ignored. 687 * @return The twist vector of this surface at the specified parametric position. 688 * @throws IllegalArgumentException if the parameter values are invalid. 689 */ 690 @Override 691 public Vector<Length> getTwistVector(GeomPoint st) { 692 validateParameter(st); 693 return getTwistVector(st.getValue(0), st.getValue(1)); 694 } 695 696 /** 697 * Return the at least 3D normal vector for this surface at the given parametric 698 * position (s,t) on this surface. 699 * 700 * @param s 1st parametric dimension distance to calculate normal for (0.0 to 1.0 701 * inclusive). 702 * @param t 2nd parametric dimension distance to calculate normal for (0.0 to 1.0 703 * inclusive). 704 * @return The at least 3D normal vector of the surface at the specified parametric 705 * position. 706 */ 707 @Override 708 public Vector<Dimensionless> getNormal(double s, double t) { 709 // Ref. 1, Eqn. 12.61. 710 711 StackContext.enter(); 712 try { 713 // Get the surface position and derivative in the s direction. 714 Vector pu = getSDerivative(s, t, 1, false); 715 716 // Get the surface derivative in the t direction. 717 Vector pw = getTDerivative(s, t, 1, false); 718 719 // Try to deal with collapsed edges by perturbing slightly off the edge. 720 if (pu.magValue() < SQRT_EPS) { 721 double offset = sqrt(TOL_ST); 722 if (t > 1 - offset) 723 offset *= -1; 724 pu = getSDerivative(s, t + offset, 1, false); 725 } 726 if (pw.magValue() < SQRT_EPS) { 727 double offset = sqrt(TOL_ST); 728 if (s > 1 - offset) 729 offset *= -1; 730 pw = getTDerivative(s + offset, t, 1, false); 731 } 732 733 // Make sure the vectors are at least 3D. 734 if (getPhyDimension() < 3) { 735 pu = pu.toDimension(3); 736 pw = pw.toDimension(3); 737 } 738 739 // Normal vector is the cross product of pu and pw vectors. 740 Vector n = pu.cross(pw).toUnitVector(); 741 742 return StackContext.outerCopy(n); 743 744 } finally { 745 StackContext.exit(); 746 } 747 } 748 749 /** 750 * Return the normal vector for this surface at the given parametric position (s,t) on 751 * this surface. 752 * 753 * @param st The parametric position to calculate the normal for. Must be a 754 * 2-dimensional point with each value in the range 0 to 1.0. Units are 755 * ignored. 756 * @return The normal vector of the surface at the specified parametric position. 757 */ 758 @Override 759 public Vector<Dimensionless> getNormal(GeomPoint st) { 760 validateParameter(st); 761 return getNormal(st.getValue(0), st.getValue(1)); 762 } 763 764 /** 765 * Return the tangent plane to this surface at the given parametric position (s,t) on 766 * this surface. 767 * 768 * @param s 1st parametric dimension distance to calculate tangent plane for (0.0 to 769 * 1.0 inclusive). 770 * @param t 2nd parametric dimension distance to calculate tangent plane for (0.0 to 771 * 1.0 inclusive). 772 * @return The tangent plane of the surface at the specified parametric position. 773 */ 774 @Override 775 public Plane getTangentPlane(double s, double t) { 776 // Ref. 1, Eqn. 12.63. 777 778 // Get the surface normal at the specified position. 779 Vector<Dimensionless> n = getNormal(s, t); 780 781 // Get the surface point stored in the normal. 782 GeomPoint p0 = n.getOrigin(); 783 784 // Create the plane. 785 Plane p = Plane.valueOf(n, p0.immutable()); 786 787 return p; 788 } 789 790 /** 791 * Return the tangent plane to this surface at the given parametric position (s,t) on 792 * this surface. 793 * 794 * @param st The parametric position to calculate the tangent plane for. Must be a 795 * 2-dimensional point with each value in the range 0 to 1.0. Units are 796 * ignored. 797 * @return The tangent plane of the surface at the specified parametric position. 798 */ 799 @Override 800 public Plane getTangentPlane(GeomPoint st) { 801 validateParameter(st); 802 return getTangentPlane(st.getValue(0), st.getValue(1)); 803 } 804 805 /** 806 * Returns the Gaussian Curvature for this surface at the given parametric position 807 * (s,t) on this surface. 808 * 809 * @param s 1st parametric dimension distance to calculate curvature for (0.0 to 1.0 810 * inclusive). 811 * @param t 2nd parametric dimension distance to calculate curvature for (0.0 to 1.0 812 * inclusive). 813 * @return The Gaussian curvature of the surface at the specified parametric position 814 * in units of 1/Length^2. 815 */ 816 @Override 817 public Parameter getGaussianCurvature(double s, double t) { 818 validateParameter(s, t); 819 s = roundParNearEnds(s); 820 t = roundParNearEnds(t); 821 822 // Ref. 1, Eqn. 12.71. 823 StackContext.enter(); 824 try { 825 // Calculate the fundamental form coefficients. 826 FastTable<Parameter<Length>> coefs = getFundamentalFormCoeffs(s, t); 827 Parameter<Length> E = coefs.get(0); 828 Parameter<Length> F = coefs.get(1); 829 Parameter<Length> G = coefs.get(2); 830 Parameter<Length> L = coefs.get(3); 831 Parameter<Length> M = coefs.get(4); 832 Parameter<Length> N = coefs.get(5); 833 834 // Calculate the Gaussian curvature. 835 Parameter num = L.times(N).minus(M.times(M)); 836 Parameter denom = E.times(G).minus(F.times(F)); 837 Parameter K = num.divide(denom); 838 839 return StackContext.outerCopy(K); 840 841 } finally { 842 StackContext.exit(); 843 } 844 } 845 846 /** 847 * Returns the Gaussian Curvature for this surface at the given parametric position 848 * (s,t) on this surface. 849 * 850 * @param st The parametric position to calculate the curvature for. Must be a 851 * 2-dimensional point with each value in the range 0 to 1.0. Units are 852 * ignored. 853 * @return The Gaussian curvature of the surface at the specified parametric position 854 * in units of 1/Length^2. 855 */ 856 @Override 857 public Parameter getGaussianCurvature(GeomPoint st) { 858 validateParameter(st); 859 return getGaussianCurvature(st.getValue(0), st.getValue(1)); 860 } 861 862 /** 863 * Returns the Mean Curvature for this surface at the given parametric position (s,t) 864 * on this surface. 865 * 866 * @param st The parametric position to calculate the curvature for. Must be a 867 * 2-dimensional point with each value in the range 0 to 1.0. Units are 868 * ignored. 869 * @return The Mean curvature of the surface at the specified parametric position in 870 * units of 1/Length. 871 */ 872 @Override 873 public Parameter getMeanCurvature(GeomPoint st) { 874 validateParameter(st); 875 return getMeanCurvature(st.getValue(0), st.getValue(1)); 876 } 877 878 /** 879 * Returns the Mean Curvature for this surface at the given parametric position (s,t) 880 * on this surface. 881 * 882 * @param s 1st parametric dimension distance to calculate curvature for (0.0 to 1.0 883 * inclusive). 884 * @param t 2nd parametric dimension distance to calculate curvature for (0.0 to 1.0 885 * inclusive). 886 * @return The Mean curvature of the surface at the specified parametric position in 887 * units of 1/Length. 888 */ 889 @Override 890 public Parameter getMeanCurvature(double s, double t) { 891 validateParameter(s, t); 892 s = roundParNearEnds(s); 893 t = roundParNearEnds(t); 894 895 // Ref. 1, Eqn. 12.72. 896 StackContext.enter(); 897 try { 898 // Calculate the fundamental form coefficients. 899 FastTable<Parameter<Length>> coefs = getFundamentalFormCoeffs(s, t); 900 Parameter<Length> E = coefs.get(0); 901 Parameter<Length> F = coefs.get(1); 902 Parameter<Length> G = coefs.get(2); 903 Parameter<Length> L = coefs.get(3); 904 Parameter<Length> M = coefs.get(4); 905 Parameter<Length> N = coefs.get(5); 906 907 // Calculate the Mean Curvature. 908 Parameter num = E.times(N).plus(G.times(L)).minus(F.times(M).times(2)); 909 Parameter denom = E.times(G).minus(F.times(F)).times(2); 910 Parameter H = num.divide(denom); 911 912 return StackContext.outerCopy(H); 913 914 } finally { 915 StackContext.exit(); 916 } 917 } 918 919 /** 920 * Returns the coefficients of the 1st and 2nd fundamental forms. The list contains, 921 * in this order: E, F, G, L, M, N. 922 */ 923 private FastTable<Parameter<Length>> getFundamentalFormCoeffs(double s, double t) { 924 // Ref. 1, Eqn. 12.51-53 & 12.58-60. 925 926 FastTable<Parameter<Length>> list = FastTable.newInstance(); 927 928 // Get the derivatives for the surface at this parametric position. 929 List<Vector<Length>> ders = getSDerivatives(s, t, 2, true); 930 Vector pu = ders.get(1); 931 Vector puu = ders.get(2); 932 933 Vector puw = getTwistVector(s, t); 934 935 ders = getTDerivatives(s, t, 2, true); 936 Vector pw = ders.get(1); 937 Vector pww = ders.get(2); 938 939 // E 940 list.add(pu.dot(pu)); 941 942 // F 943 list.add(pu.dot(pw)); 944 945 // G 946 list.add(pw.dot(pw)); 947 948 // Construct the normal vector. 949 Vector n = pu.cross(pw).toUnitVector(); 950 951 // L 952 list.add(puu.dot(n)); 953 954 // M 955 list.add(puw.dot(n)); 956 957 // N 958 list.add(pww.dot(n)); 959 960 return list; 961 } 962 963 /** 964 * Return the surface area of this entire surface. 965 * 966 * @param eps The desired fractional accuracy on the surface area. 967 * @return the surface area of this surface. 968 */ 969 @Override 970 public Parameter<Area> getArea(double eps) { 971 return getArea(0, 0, 1, 1, eps); 972 } 973 974 /** 975 * Return the surface area of a portion of this surface. 976 * 977 * @param st1 The starting parametric position to calculate the area for. Must be a 978 * 2-dimensional point with each value in the range 0 to 1.0. Units are 979 * ignored. 980 * @param st2 The ending parametric position to calculate the area for. Must be a 981 * 2-dimensional point with each value in the range 0 to 1.0. Units are 982 * ignored. 983 * @param eps The desired fractional accuracy on the surface area. 984 * @return the surface area of a portion of this surface. 985 */ 986 @Override 987 public Parameter<Area> getArea(GeomPoint st1, GeomPoint st2, double eps) { 988 validateParameter(st1); 989 validateParameter(st2); 990 return getArea(st1.getValue(0), st1.getValue(1), st2.getValue(0), st2.getValue(1), eps); 991 } 992 993 /** 994 * Return the surface area of a portion of this surface. 995 * 996 * @param s1 The starting 1st parametric dimension distance to calculate area for 997 * (0.0 to 1.0 inclusive). 998 * @param t1 The starting 2nd parametric dimension distance to calculate area for 999 * (0.0 to 1.0 inclusive). 1000 * @param s2 The ending 1st parametric dimension distance to calculate area for (0.0 1001 * to 1.0 inclusive). 1002 * @param t2 The ending 2nd parametric dimension distance to calculate area for (0.0 1003 * to 1.0 inclusive). 1004 * @param eps The desired fractional accuracy on the surface area. 1005 * @return the surface area of a portion of this surface. 1006 */ 1007 @Override 1008 public Parameter<Area> getArea(double s1, double t1, double s2, double t2, double eps) { 1009 validateParameter(s1, t1); 1010 validateParameter(s2, t2); 1011 1012 // Ref. 1, Eqn. 12.79. 1013 // If either s1 and s2 are equal or t1 and t2 are equal, then we know the answer already. 1014 Unit<Area> unit = getUnit().times(getUnit()).asType(Area.class); 1015 if (isApproxEqual(s1, s2, TOL_ST) || isApproxEqual(t1, t2, TOL_ST)) 1016 return Parameter.ZERO_LENGTH.to(unit); 1017 1018 // Make sure that s2 is > s1. 1019 if (s1 > s2) { 1020 double temp = s1; 1021 s1 = s2; 1022 s2 = temp; 1023 } 1024 1025 // Make sure that t2 is > t1. 1026 if (t1 > t2) { 1027 double temp = t1; 1028 t1 = t2; 1029 t2 = temp; 1030 } 1031 1032 // Use a concurrent context to evaluate four quarters of the surface at once. 1033 double sm = (s1 + s2) / 2; 1034 double tm = (t1 + t2) / 2; 1035 AreaEvaluatable areaEv1 = AreaEvaluatable.newInstance(this, s1, sm, t1, tm, eps); 1036 AreaEvaluatable areaEv2 = AreaEvaluatable.newInstance(this, s1, sm, tm, t2, eps); 1037 AreaEvaluatable areaEv3 = AreaEvaluatable.newInstance(this, sm, s2, t1, tm, eps); 1038 AreaEvaluatable areaEv4 = AreaEvaluatable.newInstance(this, sm, s2, tm, t2, eps); 1039 ConcurrentContext.enter(); 1040 double area; 1041 try { 1042 // Integrate to find the area: A = int_t1^t2{ int_s1^s2{ |pu x pw| ds } dt } 1043 ConcurrentContext.execute(areaEv1); 1044 ConcurrentContext.execute(areaEv2); 1045 ConcurrentContext.execute(areaEv3); 1046 ConcurrentContext.execute(areaEv4); 1047 1048 } finally { 1049 //System.out.println("area = " + A + ", numEvals = " + areaEvaluator.numEvaluations); 1050 ConcurrentContext.exit(); // Waits for all concurrent threads to complete. 1051 area = areaEv1.getValue() + areaEv2.getValue() + areaEv3.getValue() + areaEv4.getValue(); 1052 1053 AreaEvaluatable.recycle(areaEv1); 1054 AreaEvaluatable.recycle(areaEv2); 1055 AreaEvaluatable.recycle(areaEv3); 1056 AreaEvaluatable.recycle(areaEv4); 1057 } 1058 1059 Parameter<Area> A = Parameter.valueOf(area, unit); 1060 return A; 1061 } 1062 1063 /** 1064 * Return the enclosed volume of this entire surface. 1065 * 1066 * @param eps The desired fractional accuracy on the enclosed volume. 1067 * @return the enclosed volume of this surface. 1068 */ 1069 @Override 1070 public Parameter<Volume> getVolume(double eps) { 1071 return getVolume(0, 0, 1, 1, eps); 1072 } 1073 1074 /** 1075 * Return the enclosed volume of a portion of this surface. 1076 * 1077 * @param st1 The starting parametric position to calculate the volume for. Must be a 1078 * 2-dimensional point with each value in the range 0 to 1.0. Units are 1079 * ignored. 1080 * @param st2 The ending parametric position to calculate the volume for. Must be a 1081 * 2-dimensional point with each value in the range 0 to 1.0. Units are 1082 * ignored. 1083 * @param eps The desired fractional accuracy on the enclosed volume. 1084 * @return the enclosed volume of a portion of this surface. 1085 */ 1086 @Override 1087 public Parameter<Volume> getVolume(GeomPoint st1, GeomPoint st2, double eps) { 1088 validateParameter(st1); 1089 validateParameter(st2); 1090 return getVolume(st1.getValue(0), st1.getValue(1), st2.getValue(0), st2.getValue(1), eps); 1091 } 1092 1093 /** 1094 * Return the enclosed volume of a portion of this surface. 1095 * 1096 * @param s1 The starting 1st parametric dimension distance to calculate volume for 1097 * (0.0 to 1.0 inclusive). 1098 * @param t1 The starting 2nd parametric dimension distance to calculate volume for 1099 * (0.0 to 1.0 inclusive). 1100 * @param s2 The ending 1st parametric dimension distance to calculate volume for 1101 * (0.0 to 1.0 inclusive). 1102 * @param t2 The ending 2nd parametric dimension distance to calculate volume for 1103 * (0.0 to 1.0 inclusive). 1104 * @param eps The desired fractional accuracy on the enclosed volume. 1105 * @return the enclosed volume of a portion of this surface. 1106 */ 1107 @Override 1108 public Parameter<Volume> getVolume(double s1, double t1, double s2, double t2, double eps) { 1109 validateParameter(s1, t1); 1110 validateParameter(s2, t2); 1111 1112 // Ref. 1, Eqn. 12.82. 1113 // If either s1 and s2 are equal or t1 and t2 are equal, then we know the answer already. 1114 Unit<Volume> unit = getUnit().pow(3).asType(Volume.class); 1115 if (isApproxEqual(s1, s2, TOL_ST) || isApproxEqual(t1, t2, TOL_ST)) 1116 return Parameter.ZERO_VOLUME.to(unit); 1117 1118 // Make sure that s2 is > s1. 1119 if (s1 > s2) { 1120 double temp = s1; 1121 s1 = s2; 1122 s2 = temp; 1123 } 1124 1125 // Make sure that t2 is > t1. 1126 if (t1 > t2) { 1127 double temp = t1; 1128 t1 = t2; 1129 t2 = temp; 1130 } 1131 1132 // Use a concurrent context to evaluate four quarters of the surface at once. 1133 double sm = (s1 + s2) / 2; 1134 double tm = (t1 + t2) / 2; 1135 VolumeEvaluatable volEv1 = VolumeEvaluatable.newInstance(this, s1, sm, t1, tm, eps); 1136 VolumeEvaluatable volEv2 = VolumeEvaluatable.newInstance(this, s1, sm, tm, t2, eps); 1137 VolumeEvaluatable volEv3 = VolumeEvaluatable.newInstance(this, sm, s2, t1, tm, eps); 1138 VolumeEvaluatable volEv4 = VolumeEvaluatable.newInstance(this, sm, s2, tm, t2, eps); 1139 ConcurrentContext.enter(); 1140 double volume; 1141 try { 1142 // Integrate to find the volume: V = int_t1^t2{ int_s1^s2{ 1/3*(p(s,t) DOT n(s,t) ds } dt } 1143 ConcurrentContext.execute(volEv1); 1144 ConcurrentContext.execute(volEv2); 1145 ConcurrentContext.execute(volEv3); 1146 ConcurrentContext.execute(volEv4); 1147 1148 } finally { 1149 //System.out.println("volume = " + V + ", numEvals = " + volumeEvaluator.numEvaluations); 1150 ConcurrentContext.exit(); // Waits for all concurrent threads to complete. 1151 volume = volEv1.getValue() + volEv2.getValue() + volEv3.getValue() + volEv4.getValue(); 1152 1153 VolumeEvaluatable.recycle(volEv1); 1154 VolumeEvaluatable.recycle(volEv2); 1155 VolumeEvaluatable.recycle(volEv3); 1156 VolumeEvaluatable.recycle(volEv4); 1157 } 1158 1159 Parameter<Volume> V = Parameter.valueOf(volume, unit); 1160 return V; 1161 } 1162 1163 /** 1164 * Return an array of SubrangePoint objects that are gridded onto the surface using 1165 * the specified spacings and gridding rule. <code>gridRule</code> specifies whether 1166 * the subdivision spacing is applied with respect to arc-length in real space or 1167 * parameter space. If the spacing is arc-length, then the values are interpreted as 1168 * fractions of the arc length of the edge curves from 0 to 1. 1169 * 1170 * @param gridRule Specifies whether the subdivision spacing is applied with 1171 * respect to arc-length in real space or parameter space. 1172 * May not be null. 1173 * @param bottomSpacing A list of values used to define the subdivision spacing along 1174 * the T=0 parametric boundary of the surface. May not be null. 1175 * @param topSpacing An optional list of values used to define the subdivision 1176 * spacing along the T=1 parametric boundary of the surface. If 1177 * <code>null</code> is passed, the bottomSpacing is used for the 1178 * top. If topSpacing is provided, then it must have the same 1179 * number of values as the bottomSpacing. 1180 * @param leftSpacing A list of values used to define the subdivision spacing along 1181 * the S=0 parametric boundary of the surface. May not be null. 1182 * @param rightSpacing An optional list of values used to define the subdivision 1183 * spacing along the S=1 parametric boundary of the surface. If 1184 * <code>null</code> is passed, the leftSpacing is used for the 1185 * right. If rightSpacing is provided, then it must have the same 1186 * number of values as the leftSpacing. 1187 * @return An array of SubrangePoint objects gridded onto the surface at the specified 1188 * parametric positions. 1189 */ 1190 @Override 1191 public PointArray<SubrangePoint> extractGrid(GridRule gridRule, List<Double> bottomSpacing, List<Double> topSpacing, 1192 List<Double> leftSpacing, List<Double> rightSpacing) { 1193 requireNonNull(gridRule); 1194 requireNonNull(bottomSpacing); 1195 requireNonNull(leftSpacing); 1196 1197 // Check for valid inputs. 1198 if (isNull(topSpacing)) 1199 topSpacing = bottomSpacing; 1200 else if (topSpacing.size() != bottomSpacing.size()) 1201 throw new IllegalArgumentException(MessageFormat.format( 1202 RESOURCES.getString("lstSameSizeErr"),"topSpacing","bottomSpacing")); 1203 1204 if (isNull(rightSpacing)) 1205 rightSpacing = leftSpacing; 1206 else if (rightSpacing.size() != leftSpacing.size()) 1207 throw new IllegalArgumentException(MessageFormat.format( 1208 RESOURCES.getString("lstSameSizeErr"),"rightSpacing","leftSpacing")); 1209 1210 PointArray<SubrangePoint> array = null; 1211 1212 switch (gridRule) { 1213 case PAR: 1214 // Parametric spacing. 1215 array = parametricGrid(bottomSpacing, topSpacing, leftSpacing, rightSpacing); 1216 1217 break; 1218 1219 case ARC: 1220 // Arc-length spacing. 1221 1222 // Convert the ARC length spacing into parametric spacing along the T=0 & 1 boundary. 1223 FastTable<Double> parBottomSpacing = convertSubrangeCurveARCSpacing(bottomSpacing, getT0Curve()); 1224 FastTable<Double> parTopSpacing = parBottomSpacing; 1225 if (topSpacing != bottomSpacing) 1226 parTopSpacing = convertSubrangeCurveARCSpacing(topSpacing, getT1Curve()); 1227 1228 // Convert the ARC length spacing into parametric spacing along the S=0 & 1 boundary. 1229 FastTable<Double> parLeftSpacing = convertSubrangeCurveARCSpacing(leftSpacing, getS0Curve()); 1230 FastTable<Double> parRightSpacing = parLeftSpacing; 1231 if (rightSpacing != leftSpacing) 1232 parRightSpacing = convertSubrangeCurveARCSpacing(rightSpacing, getS1Curve()); 1233 1234 // Now grid the surface in parametric space. 1235 array = parametricGrid(parBottomSpacing, parTopSpacing, parLeftSpacing, parRightSpacing); 1236 1237 break; 1238 1239 default: 1240 throw new IllegalArgumentException( 1241 MessageFormat.format(RESOURCES.getString("srfUnsupportedGridRule"), 1242 gridRule.toString())); 1243 } 1244 1245 return array; 1246 } 1247 1248 /** 1249 * Convert a list of arc-length spacing values along a curve into a set of parametric 1250 * spacings along the curve. 1251 * 1252 * @param arcSpacing A list of arc-length spacings along the curve from 0 to 1. The 1253 * contents are replaced on exit by the parametric spacing along the 1254 * curve corresponding to the input arc-length spacing. 1255 * @param crv The curve to calculate the parametric spacing along. 1256 * @return A list of parameter spacings along the curve from 0 to 1. 1257 */ 1258 private static FastTable<Double> convertSubrangeCurveARCSpacing(List<Double> arcSpacing, Curve crv) { 1259 FastTable<Double> parSpacing = FastTable.newInstance(); 1260 1261 PointString<SubrangePoint> str = crv.extractGrid(GridRule.ARC, arcSpacing); 1262 int size = str.size(); 1263 for (int i = 0; i < size; ++i) { 1264 SubrangePoint p = str.get(i); 1265 parSpacing.add(p.getParPosition().getValue(0)); 1266 SubrangePoint.recycle(p); 1267 } 1268 PointString.recycle(str); 1269 1270 return parSpacing; 1271 } 1272 1273 /** 1274 * Extract a grid on this surface using a parametric spacing. 1275 */ 1276 private PointArray<SubrangePoint> parametricGrid(List<Double> bottomSpacing, List<Double> topSpacing, 1277 List<Double> leftSpacing, List<Double> rightSpacing) { 1278 PointArray<SubrangePoint> array = PointArray.newInstance(); 1279 1280 // Interpolate the spacing in the s and t directions using the 1281 // indexes along the edges. 1282 int tSize = leftSpacing.size(); 1283 int sSize = bottomSpacing.size(); 1284 for (int tIdx = 0; tIdx < tSize; ++tIdx) { 1285 double left = leftSpacing.get(tIdx); 1286 double right = rightSpacing.get(tIdx); 1287 1288 PointString<SubrangePoint> str = PointString.newInstance(); 1289 for (int sIdx = 0; sIdx < sSize; ++sIdx) { 1290 double s = bottomSpacing.get(sIdx); 1291 if (bottomSpacing != topSpacing) 1292 s = interpSpacing(s, topSpacing.get(sIdx), tSize, tIdx); 1293 double t = left; 1294 if (leftSpacing != rightSpacing) 1295 t = interpSpacing(left, right, sSize, sIdx); 1296 1297 str.add(getPoint(s, t)); 1298 } 1299 array.add(str); 1300 } 1301 1302 return array; 1303 } 1304 1305 /** 1306 * Linearly interpolate the spacing between one side of the surface and the other 1307 * based on the index which is incrementing from spacing1 to spacing2. 1308 */ 1309 private static double interpSpacing(double spacing1, double spacing2, int size, int idx) { 1310 double f = ((double)idx) / (size - 1); 1311 double s = (1 - f) * spacing1 + f * spacing2; 1312 return s; 1313 } 1314 1315 /** 1316 * Return an array of points that are gridded onto the surface with the number of 1317 * points and spacing chosen to result in straight line segments between the points 1318 * that have mid-points that are all within the specified tolerance of this surface. 1319 * 1320 * @param tol The maximum distance that a straight line between gridded points may 1321 * deviate from this surface. May not be null or zero. 1322 * @return An array of subrange points gridded onto the surface. 1323 */ 1324 @Override 1325 public PointArray<SubrangePoint> gridToTolerance(Parameter<Length> tol) { 1326 Parameter<Area> tol2 = tol.times(tol); 1327 1328 // Start with just the corner points. 1329 List<Double> spacing = GridSpacing.linear(2); 1330 PointArray<SubrangePoint> pntArray = extractGrid(GridRule.PAR, spacing, null, spacing, null); 1331 FastTable.recycle((FastTable)spacing); 1332 1333 // Refine in the s-direction (down the strings). 1334 subdivideSDir(tol2, pntArray); 1335 1336 // Now, refine in the t-direction (across the strings). 1337 subdivideTDir(tol2, pntArray); 1338 1339 // Now guess a set of points to use based on the surface topography. 1340 PointArray<SubrangePoint> guessedArray = guessGrid2TolPoints(requireNonNull(tol)); 1341 1342 // Check guessed grid-to-tolerance points and add any that fall out of tolerance. 1343 checkGuessedGrid2TolPoints(pntArray, guessedArray, tol2); 1344 PointArray.recycle(guessedArray); 1345 1346 // Re-refine in the s-direction (down the strings). 1347 subdivideSDir(tol2, pntArray); 1348 1349 // Finally, re-refine in the t-direction (across the strings). 1350 subdivideTDir(tol2, pntArray); 1351 1352 return pntArray; 1353 } 1354 1355 /** 1356 * Guess an initial set of points to use when gridding to tolerance. The initial 1357 * points must include each corner of the surface as well as any discontinuities in 1358 * the surface. Other than that, they should be distributed as best as possible to 1359 * indicate areas of higher curvature, but should be generated as efficiently as 1360 * possible. 1361 * <p> 1362 * The default implementation places a 10x10 grid of points equally spaced in 1363 * parameter space onto the surface. Sub-classes should override this method if they 1364 * can provide a more efficient/robust method to locate points where there may be 1365 * curvature discontinuities or areas of high curvature. 1366 * </p> 1367 * 1368 * @param tol The maximum distance that a straight line between gridded points may 1369 * deviate from this surface. May not be null or zero. 1370 * @return An array of subrange points to use as a starting point for the grid to 1371 * tolerance algorithms. 1372 */ 1373 protected PointArray<SubrangePoint> guessGrid2TolPoints(Parameter<Length> tol) { 1374 1375 // Create a square grid of subrange points on the surface. 1376 List<Double> spacing = GridSpacing.linear(10); 1377 PointArray<SubrangePoint> pntArray = extractGrid(GridRule.PAR, spacing, null, spacing, null); 1378 FastTable.recycle((FastTable)spacing); 1379 1380 return pntArray; 1381 } 1382 1383 private static final int MAX_POINTS = 10000; 1384 1385 /** 1386 * Subdivide the surface in the parametric s-direction starting with an initial guess 1387 * at the grid of points to be rendered. New points are added to each of the strings 1388 * in the pntArray object as needed. 1389 */ 1390 private PointArray<SubrangePoint> subdivideSDir(Parameter<Area> tol2, PointArray<SubrangePoint> pntArray) { 1391 int nTSeg = pntArray.size(); 1392 int nSSeg = pntArray.get(0).size(); 1393 int numPnts = nTSeg * nSSeg; 1394 if (numPnts > MAX_POINTS) { 1395 Logger.getLogger(AbstractSurface.class.getName()).log(Level.WARNING, 1396 MessageFormat.format(RESOURCES.getString("maxPointsWarningMsg"), 1397 "surface", this.getID())); 1398 return pntArray; 1399 } 1400 GeomPoint[] pntArr = new GeomPoint[2]; 1401 1402 // Loop over all the rows/strings in the array. 1403 outer: 1404 for (int i = 0; i < nTSeg; ++i) { 1405 PointString<SubrangePoint> curStr = pntArray.get(i); 1406 1407 // Loop over all the columns/points in each string. 1408 int jp1 = 1; 1409 for (int j = 0; j < curStr.size() - 1; ++j) { 1410 1411 // Subdivide the current segment until the mid-point is less than tol away from the surface. 1412 boolean distGTtol = true; 1413 while (distGTtol) { 1414 1415 // Get the current point. 1416 SubrangePoint pc = curStr.get(j); 1417 GeomPoint stc = pc.getParPosition(); 1418 1419 // Get the next point. 1420 SubrangePoint pn = curStr.get(jp1); 1421 GeomPoint stn = pn.getParPosition(); 1422 1423 // Compute the average point half-way on a line from pc to pn. 1424 pntArr[0] = pc; pntArr[1] = pn; 1425 Point pa = GeomUtil.averagePoints(pntArr); 1426 1427 // Find the middle point on the surface. 1428 pntArr[0] = stc; pntArr[1] = stn; 1429 Point stm = GeomUtil.averagePoints(pntArr); // Find the average parametric position. 1430 SubrangePoint pm = getPoint(stm); 1431 1432 // Compute the distance between the middle point and the average point 1433 // and compare with the tolerance. 1434 distGTtol = pa.distanceSq(pm).isGreaterThan(tol2); 1435 if (distGTtol) { 1436 // Insert a new column of points into each string in the array. 1437 insertInterpolatedColumn(pntArray, j, stm.getValue(0)); 1438 numPnts += nTSeg; 1439 if (numPnts > MAX_POINTS) { 1440 Logger.getLogger(AbstractSurface.class.getName()).log(Level.WARNING, 1441 MessageFormat.format(RESOURCES.getString("maxPointsWarningMsg"), 1442 "surface", this.getID())); 1443 break outer; 1444 } 1445 1446 } 1447 SubrangePoint.recycle(pm); 1448 1449 } // end while() 1450 1451 ++jp1; 1452 } // Next j 1453 } // Next i 1454 1455 return pntArray; 1456 } 1457 1458 /** 1459 * Subdivide the surface in the parametric t-direction starting with an initial guess 1460 * at the grid of points to be rendered. New columns/strings are added to the pntArray 1461 * object as needed. 1462 */ 1463 private PointArray<SubrangePoint> subdivideTDir(Parameter<Area> tol2, PointArray<SubrangePoint> pntArray) { 1464 int nTSeg = pntArray.size(); 1465 int nSSeg = pntArray.get(0).size(); 1466 int numPnts = nTSeg * nSSeg; 1467 if (numPnts > MAX_POINTS) { 1468 Logger.getLogger(AbstractSurface.class.getName()).log(Level.WARNING, 1469 MessageFormat.format(RESOURCES.getString("maxPointsWarningMsg"), 1470 "surface", this.getID())); 1471 return pntArray; 1472 } 1473 GeomPoint[] pntArr = new GeomPoint[2]; 1474 1475 // Loop over all the columns/points in each string of the array. 1476 outer: 1477 for (int j = 0; j < nSSeg; ++j) { 1478 1479 // Loop over all the rows/strings in the array. 1480 int ip1 = 1; 1481 for (int i = 0; i < nTSeg - 1; ++i) { 1482 1483 // Subdivide the current segment until the mid-point is less than tol away from the surface. 1484 boolean distGTtol = true; 1485 while (distGTtol) { 1486 1487 // Get the current point. 1488 SubrangePoint pc = pntArray.get(i).get(j); 1489 GeomPoint stc = pc.getParPosition(); 1490 1491 // Get the next point. 1492 SubrangePoint pn = pntArray.get(ip1).get(j); 1493 GeomPoint stn = pn.getParPosition(); 1494 1495 // Compute the average point half-way on a line from pc to pn. 1496 pntArr[0] = pc; pntArr[1] = pn; 1497 Point pa = GeomUtil.averagePoints(pntArr); 1498 1499 // Find the middle point on the surface. 1500 pntArr[0] = stc; pntArr[1] = stn; 1501 Point stm = GeomUtil.averagePoints(pntArr); // Find the average parametric position. 1502 SubrangePoint pm = getPoint(stm); 1503 1504 // Compute the distance between the middle point and the average point 1505 // and compare with the tolerance. 1506 distGTtol = pa.distanceSq(pm).isGreaterThan(tol2); 1507 if (distGTtol) { 1508 // Insert a new row/string of points into the array. 1509 insertInterpolatedRow(pntArray, i, stm.getValue(1)); 1510 ++nTSeg; 1511 numPnts += nSSeg; 1512 if (numPnts > MAX_POINTS) { 1513 Logger.getLogger(AbstractSurface.class.getName()).log(Level.WARNING, 1514 MessageFormat.format(RESOURCES.getString("maxPointsWarningMsg"), 1515 "surface", this.getID())); 1516 break outer; 1517 } 1518 1519 } 1520 SubrangePoint.recycle(pm); 1521 1522 } // end while() 1523 1524 ++ip1; 1525 } // Next i 1526 } // Next j 1527 1528 return pntArray; 1529 } 1530 1531 /** 1532 * Insert into the input pntArray a new row of points (new PointString) between rows 1533 * "i" and "i+1" at parametric position "t". 1534 * 1535 * @param pntArray The point array to add the new row of points (string of points) to. 1536 * @param i The index on the low side of where the points should be inserted. 1537 * @param t The parametric position where the points should be inserted. Must 1538 * be > the parametric "t" position of pntArray.get(i) and < the 1539 * "t" position of pntArray.get(i+1). 1540 */ 1541 private void insertInterpolatedRow(PointArray<SubrangePoint> pntArray, int i, double t) { 1542 1543 // Get the string of points on near where the new row will be inserted (used for "s" positions). 1544 PointString<SubrangePoint> row1 = pntArray.get(i); 1545 1546 // Create an empty string to fill in with interpolated points. 1547 PointString<SubrangePoint> stra = PointString.newInstance(); 1548 1549 // Create the interpolated points along a string (row). 1550 int numPnts = row1.size(); 1551 for (int j = 0; j < numPnts; ++j) { 1552 SubrangePoint p1 = row1.get(j); 1553 double s = p1.getParPosition().getValue(0); 1554 stra.add(getPoint(s, t)); 1555 } 1556 1557 // Insert the interpolated row into the input array of points. 1558 pntArray.add(i+1, stra); 1559 } 1560 1561 /** 1562 * Insert into the input pntArray a new column of points (across all the strings) 1563 * between columns "j" and "j+1" at parametric position "s". 1564 * 1565 * @param pntArray The point array to add the new column of points (across the 1566 * strings) to. 1567 * @param j The index on the low side of where the points should be inserted. 1568 * @param s The parametric position where the points should be inserted. Must 1569 * be > the parametric "s" position of pntArray.getColumn(j) and < 1570 * the "s" position of pntArray.getColumn(j+1). 1571 */ 1572 private void insertInterpolatedColumn(PointArray<SubrangePoint> pntArray, int j, double s) { 1573 1574 // Get the column of points on either side of where the column should be added. 1575 PointString<SubrangePoint> col1 = pntArray.getColumn(j); 1576 1577 // Create and insert the interpolated points across the strings (along a column). 1578 int numPnts = col1.size(); 1579 int jp1 = j + 1; 1580 for (int i = 0; i < numPnts; ++i) { 1581 SubrangePoint p1 = col1.get(i); 1582 double t = p1.getParPosition().getValue(1); 1583 SubrangePoint pm = getPoint(s, t); 1584 pntArray.get(i).add(jp1, pm); 1585 } 1586 1587 PointString.recycle(col1); 1588 } 1589 1590 /** 1591 * A Comparator that compares the parametric S-position of a pair of subrange points. 1592 */ 1593 private static class SComparator implements Comparator<SubrangePoint> { 1594 @Override 1595 public int compare(SubrangePoint o1, SubrangePoint o2) { 1596 double s1 = o1.getParPosition().getValue(0); 1597 double s2 = o2.getParPosition().getValue(0); 1598 return Double.compare(s1, s2); 1599 } 1600 } 1601 1602 /** 1603 * A Comparator that compares the parametric T-position of a pair of subrange points. 1604 */ 1605 private static class TComparator implements Comparator<SubrangePoint> { 1606 @Override 1607 public int compare(SubrangePoint o1, SubrangePoint o2) { 1608 double t1 = o1.getParPosition().getValue(1); 1609 double t2 = o2.getParPosition().getValue(1); 1610 return Double.compare(t1, t2); 1611 } 1612 } 1613 1614 private static final SComparator S_CMPRTR = new SComparator(); 1615 private static final TComparator T_CMPRTR = new TComparator(); 1616 1617 /** 1618 * Check guessed grid-to-tolerance points and add them to pntArray if they 1619 * are greater than the tolerance distance away from straight line segments 1620 * between the points already in pntArray. 1621 * 1622 * @param pntArray The existing array of gridded points on this surface. Any 1623 * points from guessPnts that are more than tol away from straight line 1624 * segments between these points are added to this array. 1625 * @param guessPnts The existing array of segment parameter derived guessed points. 1626 * @param tol2 The square of the maximum distance that a straight line between 1627 * gridded points may deviate from this surface. May not be null. 1628 */ 1629 private void checkGuessedGrid2TolPoints(PointArray<SubrangePoint> pntArray, 1630 PointArray<SubrangePoint> guessPnts, Parameter<Area> tol2) { 1631 1632 PointString<SubrangePoint> pntArrayT0 = pntArray.getColumn(0); 1633 int nSSegs = guessPnts.get(0).size() - 1; 1634 int nTSegs = guessPnts.size() - 1; 1635 for (int j=1; j < nTSegs; ++j) { 1636 // All the points in this row should have the same T position. 1637 PointString<SubrangePoint> guessPntsList = guessPnts.get(j); 1638 SubrangePoint pt = guessPntsList.get(0); 1639 1640 // Find where this point should fall in the gridded list of points in the T-direction. 1641 int idxj = Collections.binarySearch(pntArrayT0, pt, T_CMPRTR); 1642 if (idxj < 0) 1643 idxj = -idxj - 1; 1644 1645 for (int i=1; i < nSSegs; ++i) { 1646 SubrangePoint ps = guessPntsList.get(i); 1647 double s = ps.getParPosition().getValue(0); 1648 1649 // Find where this point should fall in the gridded list of points in the S-direction. 1650 int idx = Collections.binarySearch(pntArray.get(0), ps, S_CMPRTR); 1651 if (idx < 0) { 1652 // ps is not already in the pntList, so see if it is in tolerance or not. 1653 idx = -idx - 1; 1654 1655 // Calculate a point, pa, along the line between p(i) and p(i-1) 1656 // at the position of "s". 1657 SubrangePoint pi = pntArray.get(idxj,idx); 1658 SubrangePoint pim1 = pntArray.get(idxj,idx-1); 1659 Point pa = interpSPoint(pim1,pi,s); 1660 1661 // Compute the distance between the point at "s" and the line seg point 1662 // and compare with the tolerance. 1663 boolean distGTtol = pa.distanceSq(ps).isGreaterThan(tol2); 1664 if (distGTtol) { 1665 // Insert the new column of points into each of the rows (point strings). 1666 insertInterpolatedColumn(pntArray, idx-1, s); 1667 if (pntArray.size()*pntArray.get(0).size() >= MAX_POINTS) 1668 break; 1669 } 1670 } 1671 SubrangePoint.recycle(ps); 1672 } 1673 } 1674 PointString.recycle(pntArrayT0); 1675 1676 } 1677 1678 /** 1679 * Interpolate a point at "s" between the points p1 and p2 (that must have 1680 * parametric positions surrounding "s"). 1681 * 1682 * @param p1 Point with a parametric position less than "s". 1683 * @param p2 Point with a parametric position greater than "s". 1684 * @param s The parametric position at which a point is to be interpolated 1685 * between p1 and p2. 1686 * @return The interpolated point. 1687 */ 1688 private static Point interpSPoint(SubrangePoint p1, SubrangePoint p2, double s) { 1689 StackContext.enter(); 1690 try { 1691 double s1 = p1.getParPosition().getValue(0); 1692 double s2 = p2.getParPosition().getValue(0); 1693 Point Ds = p2.minus(p1); 1694 Point pout = p1.plus(Ds.times((s - s1) / (s2 - s1))); 1695 return StackContext.outerCopy(pout); 1696 } finally { 1697 StackContext.exit(); 1698 } 1699 } 1700 1701 /** 1702 * Returns the closest point on this surface to the specified point. 1703 * 1704 * @param point The point to find the closest point on this surface to. May not be null. 1705 * @param tol Fractional tolerance to refine the distance to. 1706 * @return The {@link SubrangePoint} on this surface that is closest to the specified 1707 * point. 1708 */ 1709 @Override 1710 public SubrangePoint getClosest(GeomPoint point, double tol) { 1711 double sOut, tOut; 1712 1713 StackContext.enter(); 1714 try { 1715 // Guess a location near enough to the closest point for the root solver 1716 // to get us the rest of the way there. 1717 GeomPoint parNear = guessClosestFarthest(requireNonNull(point), true); 1718 1719 SubrangePoint p = getClosestFarthest(point, parNear.getValue(0), parNear.getValue(1), tol, true); 1720 sOut = p.getParPosition().getValue(0); 1721 tOut = p.getParPosition().getValue(1); 1722 1723 } finally { 1724 StackContext.exit(); 1725 } 1726 return getPoint(sOut, tOut); 1727 } 1728 1729 /** 1730 * Returns the closest point on this surface to the specified point near the specified 1731 * parametric position. 1732 * 1733 * @param point The point to find the closest point on this surface to. May not be null. 1734 * @param nearS The parametric s-position where the search for the closest point 1735 * should begin. 1736 * @param nearT The parametric t-position where the search for the closest point 1737 * should begin. 1738 * @param tol Fractional tolerance to refine the distance to. 1739 * @return The {@link SubrangePoint} on this surface that is closest to the specified 1740 * point near the specified parametric position. 1741 */ 1742 @Override 1743 public SubrangePoint getClosest(GeomPoint point, double nearS, double nearT, double tol) { 1744 SubrangePoint p = getClosestFarthest(requireNonNull(point), nearS, nearT, tol, true); 1745 return p; 1746 } 1747 1748 /** 1749 * Returns the array of closest points on this surface to the specified array (list of 1750 * lists) of points. 1751 * 1752 * @param points A list of lists of points to find the closest point on this surface 1753 * to. May not be null. 1754 * @param tol Fractional tolerance to refine the distance to. 1755 * @return The {@link PointArray} of {@link SubrangePoint} on this surface that is 1756 * closest to the specified list of lists of points. 1757 */ 1758 @Override 1759 public PointArray<SubrangePoint> getClosest(List<? extends List<GeomPoint>> points, double tol) { 1760 PointArray arr = getClosestFarthest(requireNonNull(points), tol, true); 1761 return arr; 1762 } 1763 1764 /** 1765 * Returns the farthest point on this surface from the specified point. 1766 * 1767 * @param point The point to find the farthest point on this surface from. May not be 1768 * null. 1769 * @param tol Fractional tolerance to refine the distance to. 1770 * @return The {@link SubrangePoint} on this surface that is farthest from the 1771 * specified point. 1772 */ 1773 @Override 1774 public SubrangePoint getFarthest(GeomPoint point, double tol) { 1775 double sOut, tOut; 1776 1777 StackContext.enter(); 1778 try { 1779 // Guess a location near enough to the farthest point for the root solver 1780 // to get us the rest of the way there. 1781 GeomPoint parNear = guessClosestFarthest(requireNonNull(point), false); 1782 1783 SubrangePoint p = getClosestFarthest(point, parNear.getValue(0), parNear.getValue(1), tol, false); 1784 sOut = p.getParPosition().getValue(0); 1785 tOut = p.getParPosition().getValue(1); 1786 1787 } finally { 1788 StackContext.exit(); 1789 } 1790 1791 return getPoint(sOut, tOut); 1792 } 1793 1794 /** 1795 * Returns the farthest point on this surface from the specified point near the 1796 * specified parametric position on the surface. 1797 * 1798 * @param point The point to find the farthest point on this surface from. May not be 1799 * null. 1800 * @param nearS The parametric s-position where the search for the closest point 1801 * should begin. 1802 * @param nearT The parametric t-position where the search for the closest point 1803 * should begin. 1804 * @param tol Fractional tolerance to refine the distance to. 1805 * @return The {@link SubrangePoint} on this surface that is farthest from the 1806 * specified point. 1807 */ 1808 @Override 1809 public SubrangePoint getFarthest(GeomPoint point, double nearS, double nearT, double tol) { 1810 SubrangePoint p = getClosestFarthest(requireNonNull(point), nearS, nearT, tol, false); 1811 return p; 1812 } 1813 1814 /** 1815 * Returns the array of farthest points on this surface from the specified array (list 1816 * of lists) of points. 1817 * 1818 * @param points A list of lists of points to find the farthest point on this surface 1819 * from. May not be null. 1820 * @param tol Fractional tolerance to refine the distance to. 1821 * @return The {@link PointArray} of {@link SubrangePoint} on this surface that is 1822 * farthest from the specified list of lists of points. 1823 */ 1824 @Override 1825 public PointArray<SubrangePoint> getFarthest(List<? extends List<GeomPoint>> points, double tol) { 1826 PointArray<SubrangePoint> arr = getClosestFarthest(requireNonNull(points), tol, false); 1827 return arr; 1828 } 1829 1830 /** 1831 * Returns the array of closest/farthest points on this surface to the specified array 1832 * (list of lists) of points. 1833 * 1834 * @param points A list of lists of points to find the closest/farthest point on this 1835 * surface to. May not be null. 1836 * @param tol Fractional tolerance to refine the distance to. 1837 * @param closest Set to <code>true</code> to return the closest point or 1838 * <code>false</code> to return the farthest point. 1839 * @return The {@link PointArray} of {@link SubrangePoint} on this surface that is 1840 * closest/farthest to the specified list of lists of points. 1841 */ 1842 private PointArray<SubrangePoint> getClosestFarthest(List<? extends List<GeomPoint>> points, double tol, boolean closest) { 1843 requireNonNull(points); 1844 PointArray output = PointArray.newInstance(); 1845 boolean firstPass = true; 1846 SubrangePoint p_old = null; 1847 1848 // Loop over all the lists of points. 1849 int numDims = getPhyDimension(); 1850 int numCols = points.size(); 1851 int numRows = points.get(0).size(); 1852 for (int col = 0; col < numCols; ++col) { 1853 // Get this column/string of target points. 1854 List<GeomPoint> column = points.get(col); 1855 1856 // Create a string of output points on the surface. 1857 PointString<SubrangePoint> str = PointString.newInstance(); 1858 1859 // Loop over all the points in the supplied string of points and find the closest to each. 1860 for (int row = 0; row < numRows; ++row) { 1861 // Get the target point. 1862 GeomPoint q = column.get(row); 1863 1864 // Make sure the point and surface have the same dimensions. 1865 if (q.getPhyDimension() > numDims) 1866 throw new IllegalArgumentException(MessageFormat.format( 1867 RESOURCES.getString("sameOrFewerDimErr"), "points", "surface")); 1868 else if (q.getPhyDimension() < numDims) 1869 q = q.toDimension(numDims); 1870 1871 double nearS, nearT; 1872 if (firstPass) { 1873 firstPass = false; 1874 1875 // Guess a location near enough to the closest point for the root solver 1876 // to get us the rest of the way there. 1877 GeomPoint parNear = guessClosestFarthest(q, closest); 1878 nearS = parNear.getValue(0); 1879 nearT = parNear.getValue(1); 1880 1881 } else { 1882 // Assume that the new point will be near the last one. 1883 @SuppressWarnings("null") 1884 GeomPoint parNear = p_old.getParPosition(); 1885 nearS = parNear.getValue(0); 1886 nearT = parNear.getValue(1); 1887 } 1888 1889 // Get the closest point to the target point. 1890 SubrangePoint p = getClosestFarthest(q, nearS, nearT, tol, closest); 1891 1892 // Add the point to the string of points. 1893 str.add(p); 1894 1895 p_old = p; 1896 } 1897 1898 // Add the string of points to the array. 1899 output.add(str); 1900 1901 firstPass = true; 1902 } 1903 1904 return output; 1905 } 1906 1907 /** 1908 * Returns the closest or farthest point on this surface to the specified point. 1909 * 1910 * @param point The point to find the closest/farthest point on this surface to. May 1911 * not be null. 1912 * @param nearS The parametric s-position where the search for the closest/farthest 1913 * point should begin. 1914 * @param nearT The parametric t-position where the search for the closest/farthest 1915 * point should begin. 1916 * @param tol Fractional tolerance to refine the distance to. 1917 * @param closest Set to <code>true</code> to return the closest point or 1918 * <code>false</code> to return the farthest point. 1919 * @return The {@link SubrangePoint} on this surface that is closest/farthest to/from 1920 * the specified point. 1921 */ 1922 private SubrangePoint getClosestFarthest(GeomPoint point, double nearS, double nearT, double tol, boolean closest) { 1923 validateParameter(nearS, nearT); 1924 1925 double sc, tc; 1926 StackContext.enter(); 1927 try { 1928 // Make sure the point and surface have the same dimensions. 1929 int numDims = getPhyDimension(); 1930 if (point.getPhyDimension() > numDims) 1931 throw new IllegalArgumentException(MessageFormat.format( 1932 RESOURCES.getString("sameOrFewerDimErr"), "point", "surface")); 1933 else if (point.getPhyDimension() < numDims) 1934 point = point.toDimension(numDims); 1935 1936 // Create a list of potential closest points and add the edge curve closest points to it. 1937 FastTable<SubrangePoint> potentials = FastTable.newInstance(); 1938 potentials.add(SubrangePoint.newInstance(this, nearS, nearT)); 1939 if (closest) { 1940 SubrangePoint p = getT0Curve().getClosest(point, nearS, tol); 1941 potentials.add(SubrangePoint.newInstance(this, p.getParPosition().getValue(0), 0)); 1942 p = getT1Curve().getClosest(point, nearS, tol); 1943 potentials.add(SubrangePoint.newInstance(this, p.getParPosition().getValue(0), 1)); 1944 p = getS0Curve().getClosest(point, nearT, tol); 1945 potentials.add(SubrangePoint.newInstance(this, 0, p.getParPosition().getValue(0))); 1946 p = getS1Curve().getClosest(point, nearT, tol); 1947 potentials.add(SubrangePoint.newInstance(this, 1, p.getParPosition().getValue(0))); 1948 1949 } else { 1950 SubrangePoint p = getT0Curve().getFarthest(point, nearS, tol); 1951 potentials.add(SubrangePoint.newInstance(this, p.getParPosition().getValue(0), 0)); 1952 p = getT1Curve().getFarthest(point, nearS, tol); 1953 potentials.add(SubrangePoint.newInstance(this, p.getParPosition().getValue(0), 1)); 1954 p = getS0Curve().getFarthest(point, nearT, tol); 1955 potentials.add(SubrangePoint.newInstance(this, 0, p.getParPosition().getValue(0))); 1956 p = getS1Curve().getFarthest(point, nearT, tol); 1957 potentials.add(SubrangePoint.newInstance(this, 1, p.getParPosition().getValue(0))); 1958 } 1959 1960 SubrangePoint output; 1961 try { 1962 // Find the closest point interior to the surface near nearS and nearT. 1963 SubrangePoint spnt = getClosestFarthestInterior(point, nearS, nearT, tol, closest); 1964 1965 // Store the point from the minimizer. 1966 potentials.add(spnt); 1967 1968 } catch (RootException err) { 1969 Logger.getLogger(AbstractSurface.class.getName()).log(Level.WARNING, 1970 "Failed to find closest/farthest interior point.", err); 1971 1972 } finally { 1973 if (potentials.isEmpty()) { 1974 output = getPoint(0, 0); 1975 1976 } else { 1977 // Loop over the potentials and find the one that is the closest/farthest. 1978 output = potentials.get(0); 1979 double dOpt = point.distanceValue(output); 1980 int size = potentials.size(); 1981 for (int i = 1; i < size; ++i) { 1982 SubrangePoint p = potentials.get(i); 1983 double dist = point.distanceValue(p); 1984 if (closest) { 1985 if (dist < dOpt) { 1986 dOpt = dist; 1987 output = p; 1988 } 1989 } else { 1990 if (dist > dOpt) { 1991 dOpt = dist; 1992 output = p; 1993 } 1994 } 1995 } 1996 1997 } 1998 1999 // Store the closest point S & T parameters. 2000 sc = output.getParPosition().getValue(0); 2001 tc = output.getParPosition().getValue(1); 2002 } 2003 } finally { 2004 StackContext.exit(); 2005 } 2006 2007 return getPoint(sc,tc); 2008 } 2009 2010 /** 2011 * Find and return the nearest point interior to the surface. This does not look at 2012 * edge curves and it assumes all the inputs are properly set up. 2013 * 2014 * @param point The point to find the closest/farthest point on this surface to. May 2015 * not be null. 2016 * @param nearS The parametric s-position where the search for the closest/farthest 2017 * point should begin. 2018 * @param nearT The parametric t-position where the search for the closest/farthest 2019 * point should begin. 2020 * @param tol Fractional tolerance to refine the distance to. 2021 * @param closest Set to <code>true</code> to return the closest point or 2022 * <code>false</code> to return the farthest point. 2023 * @return The SubrangePoint on the interior of this surface that is closest/farthest 2024 * to/from the specified point. 2025 * @throws RootException if there is a convergence problem with the root solver. 2026 */ 2027 private SubrangePoint getClosestFarthestInterior(GeomPoint point, double nearS, double nearT, 2028 double tol, boolean closest) throws RootException { 2029 validateParameter(nearS, nearT); 2030 requireNonNull(point); 2031 2032 double s, t; 2033 StackContext.enter(); 2034 try { 2035 // Create a function that calculates the distance from point to the surface. 2036 PointSurfaceDistFunction distFunc = PointSurfaceDistFunction.newInstance(this, point, closest); 2037 2038 // Use n-dimensional minimizer to locate closest point near "pos". 2039 double[] x = ArrayFactory.DOUBLES_FACTORY.array(2); 2040 x[0] = nearS; 2041 x[1] = nearT; 2042 Minimization.findND(distFunc, x, 2, tol); 2043 //System.out.println("numEvals = " + distFunc.numEvaluations); 2044 2045 // The minimizer may go slightly past the parametric bounds. Deal with that. 2046 s = x[0]; 2047 t = x[1]; 2048 if (s < 0) 2049 s = 0; 2050 else if (s > 1) 2051 s = 1; 2052 if (t < 0) 2053 t = 0; 2054 else if (t > 1) 2055 t = 1; 2056 2057 } finally { 2058 StackContext.exit(); 2059 } 2060 2061 return getPoint(s,t); 2062 } 2063 2064 /** 2065 * Method that guesses the most likely location for the closest or farthest point on a 2066 * surface and returns that location as a 2D point containing the "s" and "t" 2067 * parameter values. This is called by getClosest() and getFarthest(). This is 2068 * required in order for the root-finding algorithm to reliably refine the closest 2069 * point to the correct location. 2070 * <p> 2071 * Subclasses may override this method to provide a method for estimating the closest 2072 * point location that is most efficient for that surface type. The default 2073 * implementation extracts a 9x9 grid of subrange points, finds the closest one and 2074 * returns it's parametric position. 2075 * </p> 2076 * 2077 * @param point The point to find the closest point on this surface to. May not be null. 2078 * @param closest Set to <code>true</code> to return the closest point or 2079 * <code>false</code> to return the farthest point. 2080 * @return The 2D parametric point on this surface that is approx. closest to the 2081 * specified point. 2082 */ 2083 protected GeomPoint guessClosestFarthest(GeomPoint point, boolean closest) { 2084 // Create a 9x9 grid of subrange points on the surface. 2085 // Return the parametric position of the closest/farthest grid point. 2086 return getClosestFarthestOnGrid(this, requireNonNull(point), closest, 9).getParPosition(); 2087 } 2088 2089 /** 2090 * Extract a regularly spaced grid of parametric points on the surface and return the 2091 * grid point that is closest/farthest to/from the target point. This ignores the 2092 * boundary points as they are handled by extracting the boundary curves. 2093 * 2094 * @param srf The surface to find the closest grid point on. 2095 * @param point The point to find the closest point on this surface to. 2096 * @param closest Set to <code>true</code> to return the closest point or 2097 * <code>false</code> to return the farthest point. 2098 * @param gridSize The number of points to space out on the surface in each parametric 2099 * direction. 2100 * @return The subrange grid point on this surface that is closest to the specified 2101 * point. 2102 */ 2103 protected static SubrangePoint getClosestFarthestOnGrid(Surface srf, GeomPoint point, boolean closest, int gridSize) { 2104 double sOut, tOut; 2105 2106 StackContext.enter(); 2107 try { 2108 // Create a square grid of subrange points on the surface. 2109 List<Double> spacing = GridSpacing.linear(gridSize); 2110 spacing.remove(0); 2111 spacing.remove(spacing.size() - 1); 2112 PointArray<SubrangePoint> potentials = srf.extractGrid(GridRule.PAR, spacing, null, spacing, null); 2113 2114 // Find the closest/farthest point in the grid. 2115 SubrangePoint pOpt = potentials.get(0, 0); 2116 double dOpt2 = point.distanceSqValue(pOpt); 2117 int cols = potentials.size(); 2118 int rows = potentials.get(0).size(); 2119 for (int i = 0; i < cols; ++i) { 2120 for (int j = 0; j < rows; ++j) { 2121 SubrangePoint p = potentials.get(i, j); 2122 double dist2 = point.distanceSqValue(p); 2123 if (closest) { 2124 if (dist2 < dOpt2) { 2125 dOpt2 = dist2; 2126 pOpt = p; 2127 } 2128 } else { 2129 if (dist2 > dOpt2) { 2130 dOpt2 = dist2; 2131 pOpt = p; 2132 } 2133 } 2134 } 2135 } 2136 sOut = pOpt.getParPosition().getValue(0); 2137 tOut = pOpt.getParPosition().getValue(1); 2138 2139 } finally { 2140 StackContext.exit(); 2141 } 2142 2143 // Return the the closest/farthest point. 2144 return srf.getPoint(sOut, tOut); 2145 } 2146 2147 /** 2148 * Returns the closest points (giving the minimum distance) between this surface and 2149 * the specified curve. If the specified curve intersects this surface, then the 1st 2150 * intersection found is returned (not all the intersections are found and there is no 2151 * guarantee about which intersection will be returned). If the specified curve is 2152 * parallel to this surface (all points are equidistant away), then any point between 2153 * this surface and the specified curve could be returned as the "closest". 2154 * 2155 * @param curve The curve to find the closest points between this surface and that 2156 * curve on. May not be null. 2157 * @param tol Fractional tolerance (in parameter space) to refine the point position 2158 * to. 2159 * @return A list containing two SubrangePoint objects that represent the closest 2160 * point on this surface (index 0) and the closest point on the specified 2161 * curve (index 1). 2162 */ 2163 @Override 2164 public PointString<SubrangePoint> getClosest(Curve curve, double tol) { 2165 2166 // Make sure that curve and this surface have the same dimensions. 2167 int numDims = getPhyDimension(); 2168 if (curve.getPhyDimension() > numDims) 2169 throw new IllegalArgumentException(MessageFormat.format( 2170 RESOURCES.getString("sameOrFewerDimErr"), "curve", "surface")); 2171 else if (curve.getPhyDimension() < numDims) 2172 curve = curve.toDimension(numDims); 2173 2174 // Use the curve's "getClosest" method to find the closest points between this surface 2175 // and the target curve. 2176 PointString<SubrangePoint> str = curve.getClosest(this, tol); 2177 2178 // Reverse the string for output to match this method's contract. 2179 PointString<SubrangePoint> output = str.reverse(); 2180 PointString.recycle(str); 2181 2182 return output; 2183 } 2184 2185 /** 2186 * Returns the closest points (giving the minimum distance) between this surface and 2187 * the specified surface. If the specified surface intersects this surface, then the 2188 * 1st intersection found is returned (not all the intersections are found and there 2189 * is no guarantee about which intersection will be returned). If the specified 2190 * surface is parallel to this surface (all points are equidistant away), then any 2191 * point between this surface and the specified surface could be returned as the 2192 * "closest". 2193 * 2194 * @param surface The surface to find the closest points between this surface and that 2195 * surface on. May not be null. 2196 * @param tol Fractional tolerance (in parameter space) to refine the point 2197 * position to. 2198 * @return A list containing two SubrangePoint objects that represent the closest 2199 * point on this surface (index 0) and the closest point on the specified 2200 * surface (index 1). 2201 */ 2202 @Override 2203 public PointString<SubrangePoint> getClosest(Surface surface, double tol) { 2204 2205 // Make sure that surface and this surface have the same dimensions. 2206 int numDims = getPhyDimension(); 2207 if (surface.getPhyDimension() > numDims) 2208 throw new IllegalArgumentException(MessageFormat.format( 2209 RESOURCES.getString("sameOrFewerDimErr"), "input surface", "surface")); 2210 else if (surface.getPhyDimension() < numDims) 2211 surface = surface.toDimension(numDims); 2212 2213 // Guess a location near enough to the closest point for the root solver 2214 // to get us the rest of the way there. 2215 GeomPoint parNear = guessClosestFarthest(surface, true); 2216 2217 // Find the closest point between the two surfaces. 2218 PointString<SubrangePoint> p = getClosest(surface, parNear.getValue(0), parNear.getValue(1), tol); 2219 2220 return p; 2221 2222 } 2223 2224 /** 2225 * Method that guesses the most likely location for the closest or farthest point on a 2226 * surface and returns that location as a 2D point containing the "s" and "t" 2227 * parameter values. This is called by getClosest() and getFarthest(). This is 2228 * required in order for the root-finding algorithm to reliably refine the closest 2229 * point to the correct location. 2230 * <p> 2231 * Subclasses may override this method to provide a method for estimating the closest 2232 * point location that is most efficient for that surface type. The default 2233 * implementation extracts a 9x9 grid of subrange points, finds the closest one and 2234 * returns it's parametric position. 2235 * </p> 2236 * 2237 * @param surface The surface to find the closest point on this surface to. May not be 2238 * null. 2239 * @param closest Set to <code>true</code> to return the closest point or 2240 * <code>false</code> to return the farthest point. 2241 * @return The 2D parametric point on this surface that is approx. closest to the 2242 * specified surface. 2243 */ 2244 protected GeomPoint guessClosestFarthest(Surface surface, boolean closest) { 2245 // Create a 9x9 grid of subrange points on the surface. 2246 // Return the parametric position of the closest/farthest grid point. 2247 return getClosestFarthestOnGrid(this, requireNonNull(surface), closest, 9).getParPosition(); 2248 } 2249 2250 /** 2251 * Extract a regularly spaced grid of parametric points on the surface and return the 2252 * grid point that is closest/farthest to/from the target surface. This ignores the 2253 * boundary points as they are handled by extracting the boundary curves. 2254 * 2255 * @param thisSrf The surface to find the closest grid point on. 2256 * @param thatSrf The surface to find the closest point on this surface to. 2257 * @param closest Set to <code>true</code> to return the closest point or 2258 * <code>false</code> to return the farthest point. 2259 * @param gridSize The number of points to space out on this surface in each 2260 * parametric direction. 2261 * @return The subrange grid point on this surface that is closest to the specified 2262 * surface. 2263 */ 2264 protected static SubrangePoint getClosestFarthestOnGrid(Surface thisSrf, Surface thatSrf, boolean closest, int gridSize) { 2265 double sOut, tOut; 2266 2267 StackContext.enter(); 2268 try { 2269 // Create a square grid of subrange points on the surface. 2270 List<Double> spacing = GridSpacing.linear(gridSize); 2271 spacing.remove(0); 2272 spacing.remove(spacing.size() - 1); 2273 PointArray<SubrangePoint> potentials = thisSrf.extractGrid(GridRule.PAR, spacing, null, spacing, null); 2274 2275 // Find the closest/farthest point in the grid. 2276 SubrangePoint pOpt = potentials.get(0, 0); 2277 SubrangePoint point = thatSrf.getClosest(pOpt, GTOL); 2278 double dOpt2 = point.distanceSqValue(pOpt); 2279 int cols = potentials.size(); 2280 int rows = potentials.get(0).size(); 2281 for (int i = 0; i < cols; ++i) { 2282 for (int j = 0; j < rows; ++j) { 2283 SubrangePoint p = potentials.get(i, j); 2284 point = thatSrf.getClosest(p, GTOL); 2285 double dist2 = point.distanceSqValue(p); 2286 if (closest) { 2287 if (dist2 < dOpt2) { 2288 dOpt2 = dist2; 2289 pOpt = p; 2290 } 2291 } else { 2292 if (dist2 > dOpt2) { 2293 dOpt2 = dist2; 2294 pOpt = p; 2295 } 2296 } 2297 } 2298 } 2299 sOut = pOpt.getParPosition().getValue(0); 2300 tOut = pOpt.getParPosition().getValue(1); 2301 2302 } finally { 2303 StackContext.exit(); 2304 } 2305 2306 // Return the the closest/farthest point. 2307 return thisSrf.getPoint(sOut, tOut); 2308 } 2309 2310 /** 2311 * Returns the closest point on this surface to the specified surface. 2312 * 2313 * @param thatSrf The surface to find the closest point on this surface to. May not be 2314 * null. 2315 * @param nearS The parametric s-position where the search for the closest point 2316 * should begin. 2317 * @param nearT The parametric t-position where the search for the closest point 2318 * should begin. 2319 * @param tol Fractional tolerance (in parameter space) to refine the point 2320 * position to. 2321 * @return A list containing SubrangePoint objects representing the closest points on 2322 * this surface (index 0) and the target surface (index 1). 2323 */ 2324 private PointString<SubrangePoint> getClosest(Surface thatSrf, double nearS, double nearT, double tol) { 2325 double thisS, thisT, thatS, thatT; 2326 2327 StackContext.enter(); 2328 try { 2329 SubrangePoint thisOutput, thatOutput; 2330 2331 // Create a function that calculates the distance from this surface to that surface. 2332 SrfSrfDistFunction distFunc = SrfSrfDistFunction.newInstance(this, requireNonNull(thatSrf), tol, true); 2333 2334 // Create a list of potential closest points and add the edge curve closest points to it. 2335 FastTable<PointString<SubrangePoint>> potentials = FastTable.newInstance(); 2336 SubrangePoint p = SubrangePoint.newInstance(this, nearS, nearT); 2337 PointString<SubrangePoint> str = PointString.valueOf(p, thatSrf.getClosest(p, tol)); 2338 potentials.add(str); 2339 2340 // T = 0 curve. 2341 PointString<SubrangePoint> pc = getT0Curve().getClosest(thatSrf, tol); 2342 p = SubrangePoint.newInstance(this, pc.get(0).getParPosition().getValue(0), 0); 2343 str = PointString.valueOf(p, pc.get(1)); 2344 potentials.add(str); 2345 2346 // T = 1 curve. 2347 pc = getT1Curve().getClosest(thatSrf, tol); 2348 p = SubrangePoint.newInstance(this, pc.get(0).getParPosition().getValue(0), 1); 2349 str = PointString.valueOf(p, pc.get(1)); 2350 potentials.add(str); 2351 2352 // S = 0 curve. 2353 pc = getS0Curve().getClosest(thatSrf, tol); 2354 p = SubrangePoint.newInstance(this, 0, pc.get(0).getParPosition().getValue(0)); 2355 str = PointString.valueOf(p, pc.get(1)); 2356 potentials.add(str); 2357 2358 // S = 1 curve. 2359 pc = getS1Curve().getClosest(thatSrf, tol); 2360 p = SubrangePoint.newInstance(this, 1, pc.get(0).getParPosition().getValue(0)); 2361 str = PointString.valueOf(p, pc.get(1)); 2362 potentials.add(str); 2363 2364 try { 2365 2366 // Use n-dimensional minimizer to locate closest point near "pos". 2367 double[] x = new double[2]; 2368 x[0] = nearS; 2369 x[1] = nearT; 2370 Minimization.findND(distFunc, x, 2, tol); 2371 2372 // The minimizer may go slightly past the parametric bounds. Deal with that. 2373 double s = x[0], t = x[1]; 2374 if (s < 0) 2375 s = 0; 2376 else if (s > 1) 2377 s = 1; 2378 if (t < 0) 2379 t = 0; 2380 else if (t > 1) 2381 t = 1; 2382 2383 // Store the point from the minimizer. 2384 p = SubrangePoint.newInstance(this, s, t); 2385 str = PointString.valueOf(p, thatSrf.getClosest(p, tol)); 2386 potentials.add(str); 2387 2388 } catch (RootException err) { 2389 Logger.getLogger(AbstractSurface.class.getName()).log(Level.WARNING, 2390 "Failed to find closest point near 'pos'.", err); 2391 2392 } finally { 2393 if (potentials.isEmpty()) { 2394 thisOutput = getPoint(0, 0); 2395 thatOutput = thatSrf.getClosest(thisOutput, tol); 2396 2397 } else { 2398 // Loop over the potentials and find the one that is the closest. 2399 thisOutput = potentials.get(0).get(0); 2400 thatOutput = potentials.get(0).get(1); 2401 double dOpt = thatOutput.distanceValue(thisOutput); 2402 int size = potentials.size(); 2403 for (int i = 1; i < size; ++i) { 2404 p = potentials.get(i).get(0); 2405 SubrangePoint q = potentials.get(i).get(1); 2406 double dist = q.distanceValue(p); 2407 if (dist < dOpt) { 2408 dOpt = dist; 2409 thisOutput = p; 2410 thatOutput = q; 2411 } 2412 } 2413 2414 } 2415 } 2416 2417 // Gather the output parametric positions. 2418 thisS = thisOutput.getParPosition().getValue(0); 2419 thisT = thisOutput.getParPosition().getValue(1); 2420 thatS = thatOutput.getParPosition().getValue(0); 2421 thatT = thatOutput.getParPosition().getValue(1); 2422 2423 } finally { 2424 StackContext.exit(); 2425 } 2426 2427 PointString<SubrangePoint> output = PointString.newInstance(); 2428 output.add(getPoint(thisS, thisT)); 2429 output.add(getPoint(thatS, thatT)); 2430 2431 return output; 2432 } 2433 2434 /** 2435 * Returns the closest point (giving the minimum distance) between this surface and 2436 * the specified plane. If the specified plane intersects this surface, then the 1st 2437 * intersection found is returned (not all the intersections are found and there is no 2438 * guarantee about which intersection will be returned). If the specified plane is 2439 * parallel to this surface (all points are equidistant away), then any point between 2440 * this surface and the specified plane could be returned as the "closest". 2441 * 2442 * @param plane The plane to find the closest points between this surface and that 2443 * plane on. May not be null. 2444 * @param tol Fractional tolerance (in parameter space) to refine the point position 2445 * to. 2446 * @return A SubrangePoint that represent the closest point on this surface to the 2447 * specified plane. 2448 */ 2449 @Override 2450 public SubrangePoint getClosest(GeomPlane plane, double tol) { 2451 double sOut, tOut; 2452 2453 StackContext.enter(); 2454 try { 2455 // Make sure that plane and this surface have the same dimensions. 2456 int numDims = getPhyDimension(); 2457 if (plane.getPhyDimension() > numDims) 2458 throw new IllegalArgumentException(MessageFormat.format( 2459 RESOURCES.getString("sameOrFewerDimErr"), "plane", "surface")); 2460 else if (plane.getPhyDimension() < numDims) 2461 plane = plane.toDimension(numDims); 2462 2463 // Guess a location near enough to the closest point for the root solver 2464 // to get us the rest of the way there. 2465 GeomPoint parNear = guessClosestFarthest(plane, true); 2466 2467 // Find the closest point between this surface and that plane. 2468 SubrangePoint p = getClosest(plane, parNear.getValue(0), parNear.getValue(1), tol); 2469 sOut = p.getParPosition().getValue(0); 2470 tOut = p.getParPosition().getValue(1); 2471 2472 } finally { 2473 StackContext.exit(); 2474 } 2475 return getPoint(sOut, tOut); 2476 } 2477 2478 /** 2479 * Method that guesses the most likely location for the closest or farthest point on a 2480 * surface and returns that location as a 2D point containing the "s" and "t" 2481 * parameter values. This is called by getClosest() and getFarthest(). This is 2482 * required in order for the root-finding algorithm to reliably refine the closest 2483 * point to the correct location. 2484 * <p> 2485 * Subclasses may override this method to provide a method for estimating the closest 2486 * point location that is most efficient for that surface type. The default 2487 * implementation extracts a 9x9 grid of subrange points, finds the closest one and 2488 * returns it's parametric position. 2489 * </p> 2490 * 2491 * @param plane The plane to find the closest point on this surface to. May not be 2492 * null. 2493 * @param closest Set to <code>true</code> to return the closest point or 2494 * <code>false</code> to return the farthest point. 2495 * @return The 2D parametric point on this surface that is approx. closest to the 2496 * specified plane. 2497 */ 2498 protected GeomPoint guessClosestFarthest(GeomPlane plane, boolean closest) { 2499 // Create a 9x9 grid of subrange points on the surface. 2500 // Return the parametric position of the closest/farthest grid point. 2501 return getClosestFarthestOnGrid(this, requireNonNull(plane), closest, 9).getParPosition(); 2502 } 2503 2504 /** 2505 * Extract a regularly spaced grid of parametric points on the surface and return the 2506 * grid point that is closest/farthest to/from the target plane. This ignores the 2507 * boundary points as they are handled by extracting the boundary curves. 2508 * 2509 * @param thisSrf The surface to find the closest grid point on. 2510 * @param thatPlane The plane to find the closest point on this surface to. 2511 * @param closest Set to <code>true</code> to return the closest point or 2512 * <code>false</code> to return the farthest point. 2513 * @param gridSize The number of points to space out on this surface in each 2514 * parametric direction. 2515 * @return The subrange grid point on this surface that is closest to the specified 2516 * plane. 2517 */ 2518 protected static SubrangePoint getClosestFarthestOnGrid(Surface thisSrf, GeomPlane thatPlane, boolean closest, int gridSize) { 2519 double sOut, tOut; 2520 2521 StackContext.enter(); 2522 try { 2523 // Create a square grid of subrange points on the surface. 2524 List<Double> spacing = GridSpacing.linear(gridSize); 2525 spacing.remove(0); 2526 spacing.remove(spacing.size() - 1); 2527 PointArray<SubrangePoint> potentials = thisSrf.extractGrid(GridRule.PAR, spacing, null, spacing, null); 2528 FastTable.recycle((FastTable) spacing); 2529 2530 // Find the closest/farthest point in the grid. 2531 SubrangePoint pOpt = potentials.get(0, 0); 2532 Point point = thatPlane.getClosest(pOpt); 2533 double dOpt2 = point.distanceSqValue(pOpt); 2534 int cols = potentials.size(); 2535 int rows = potentials.get(0).size(); 2536 for (int i = 0; i < cols; ++i) { 2537 for (int j = 0; j < rows; ++j) { 2538 SubrangePoint p = potentials.get(i, j); 2539 point = thatPlane.getClosest(p); 2540 double dist2 = point.distanceSqValue(p); 2541 if (closest) { 2542 if (dist2 < dOpt2) { 2543 dOpt2 = dist2; 2544 pOpt = p; 2545 } 2546 } else { 2547 if (dist2 > dOpt2) { 2548 dOpt2 = dist2; 2549 pOpt = p; 2550 } 2551 } 2552 } 2553 } 2554 sOut = pOpt.getParPosition().getValue(0); 2555 tOut = pOpt.getParPosition().getValue(1); 2556 2557 } finally { 2558 StackContext.exit(); 2559 } 2560 // Return the the closest/farthest point. 2561 return thisSrf.getPoint(sOut, tOut); 2562 } 2563 2564 /** 2565 * Returns the closest point on this surface to the specified plane. 2566 * 2567 * @param thatPlane The plane to find the closest point on this surface to. 2568 * @param nearS The parametric s-position where the search for the closest point 2569 * should begin. 2570 * @param nearT The parametric t-position where the search for the closest point 2571 * should begin. 2572 * @param tol Fractional tolerance (in parameter space) to refine the point 2573 * position to. 2574 * @return The closest points on this surface to the specified plane. 2575 */ 2576 private SubrangePoint getClosest(GeomPlane thatPlane, double nearS, double nearT, double tol) { 2577 double sOut, tOut; 2578 2579 StackContext.enter(); 2580 try { 2581 SubrangePoint thisOutput; 2582 2583 // Create a function that calculates the distance from this surface to that plane. 2584 SrfPlaneDistFunction distFunc = SrfPlaneDistFunction.newInstance(this, thatPlane, true); 2585 2586 // Create a list of potential closest points and add the edge curve closest points to it. 2587 FastTable<PointString<SubrangePoint>> potentials = FastTable.newInstance(); 2588 SubrangePoint p = SubrangePoint.newInstance(this, nearS, nearT); 2589 PointString str = PointString.valueOf(p, thatPlane.getClosest(p)); 2590 potentials.add(str); 2591 2592 // T = 0 curve. 2593 SubrangePoint pc = getT0Curve().getClosest(thatPlane, tol); 2594 p = SubrangePoint.newInstance(this, pc.getParPosition().getValue(0), 0); 2595 str = PointString.valueOf(p, thatPlane.getClosest(pc)); 2596 potentials.add(str); 2597 2598 // T = 1 curve. 2599 pc = getT1Curve().getClosest(thatPlane, tol); 2600 p = SubrangePoint.newInstance(this, pc.getParPosition().getValue(0), 1); 2601 str = PointString.valueOf(p, thatPlane.getClosest(pc)); 2602 potentials.add(str); 2603 2604 // S = 0 curve. 2605 pc = getS0Curve().getClosest(thatPlane, tol); 2606 p = SubrangePoint.newInstance(this, 0, pc.getParPosition().getValue(0)); 2607 str = PointString.valueOf(p, thatPlane.getClosest(pc)); 2608 potentials.add(str); 2609 2610 // S = 1 curve. 2611 pc = getS1Curve().getClosest(thatPlane, tol); 2612 p = SubrangePoint.newInstance(this, 1, pc.getParPosition().getValue(0)); 2613 str = PointString.valueOf(p, thatPlane.getClosest(pc)); 2614 potentials.add(str); 2615 2616 try { 2617 2618 // Use n-dimensional minimizer to locate closest point near "pos". 2619 double[] x = new double[2]; 2620 x[0] = nearS; 2621 x[1] = nearT; 2622 Minimization.findND(distFunc, x, 2, tol); 2623 2624 // The minimizer may go slightly past the parametric bounds. Deal with that. 2625 double s = x[0], t = x[1]; 2626 if (s < 0) 2627 s = 0; 2628 else if (s > 1) 2629 s = 1; 2630 if (t < 0) 2631 t = 0; 2632 else if (t > 1) 2633 t = 1; 2634 2635 // Store the point from the minimizer. 2636 p = SubrangePoint.newInstance(this, s, t); 2637 str = PointString.valueOf(p, thatPlane.getClosest(p)); 2638 potentials.add(str); 2639 2640 } catch (RootException err) { 2641 Logger.getLogger(AbstractSurface.class.getName()).log(Level.WARNING, 2642 "Failed to find closest point near 'pos'.", err); 2643 2644 } finally { 2645 if (potentials.isEmpty()) { 2646 thisOutput = getPoint(0, 0); 2647 2648 } else { 2649 // Loop over the potentials and find the one that is the closest. 2650 thisOutput = potentials.get(0).get(0); 2651 GeomPoint q = potentials.get(0).get(1); 2652 double dOpt2 = q.distanceSqValue(thisOutput); 2653 int size = potentials.size(); 2654 for (int i = 1; i < size; ++i) { 2655 p = potentials.get(i).get(0); 2656 q = potentials.get(i).get(1); 2657 double dist2 = q.distanceSqValue(p); 2658 if (dist2 < dOpt2) { 2659 dOpt2 = dist2; 2660 thisOutput = p; 2661 } 2662 } 2663 2664 } 2665 } 2666 sOut = thisOutput.getParPosition().getValue(0); 2667 tOut = thisOutput.getParPosition().getValue(1); 2668 2669 } finally { 2670 StackContext.exit(); 2671 } 2672 2673 return getPoint(sOut, tOut); 2674 } 2675 2676 /** 2677 * Returns the most extreme point, either minimum or maximum, in the specified 2678 * coordinate direction on this surface. This is more accurate than 2679 * <code>getBoundsMax()</code> & <code>getBoundsMin()</code>, but also typically takes 2680 * longer to compute. 2681 * 2682 * @param dim An index indicating the dimension to find the min/max point for (0=X, 2683 * 1=Y, 2=Z, etc). 2684 * @param max Set to <code>true</code> to return the maximum value, <code>false</code> 2685 * to return the minimum. 2686 * @param tol Fractional tolerance (in parameter space) to refine the min/max point 2687 * position to. 2688 * @return The point found on this surface that is the min or max in the specified 2689 * coordinate direction. 2690 * @see #getBoundsMin 2691 * @see #getBoundsMax 2692 */ 2693 @Override 2694 public SubrangePoint getLimitPoint(int dim, boolean max, double tol) { 2695 // Check the input dimension for consistancy. 2696 int thisDim = getPhyDimension(); 2697 if (dim < 0 || dim >= thisDim) 2698 throw new DimensionException(MessageFormat.format( 2699 RESOURCES.getString("inputDimOutOfRange"), "surface")); 2700 double sOut, tOut; 2701 2702 StackContext.enter(); 2703 try { 2704 // Create a vector pointing in the indicated direction. 2705 MutablePoint p = MutablePoint.newInstance(thisDim); 2706 p.set(dim, Parameter.valueOf(1, SI.METER)); 2707 Vector<Dimensionless> uv = Vector.valueOf(p).toUnitVector(); 2708 2709 // Get the bounds of the surface and bias it off 1% in the specified dimension. 2710 Point boundsMax = getBoundsMax(); 2711 Point boundsMin = getBoundsMin(); 2712 Parameter<Length> range = boundsMax.get(dim).minus(boundsMin.get(dim)); 2713 MutablePoint pp; 2714 if (max) { 2715 pp = MutablePoint.valueOf(boundsMax); 2716 pp.set(dim, pp.get(dim).plus(range.times(0.01))); 2717 } else { 2718 pp = MutablePoint.valueOf(boundsMin); 2719 pp.set(dim, pp.get(dim).minus(range.times(0.01))); 2720 } 2721 2722 // Create a plane 1% of the diagonal off the bounds. 2723 Plane plane = Plane.valueOf(uv, pp.immutable()); 2724 2725 // Return the closest point to the specified plane. 2726 SubrangePoint pnt = getClosest(plane, tol); 2727 sOut = pnt.getParPosition().getValue(0); 2728 tOut = pnt.getParPosition().getValue(1); 2729 2730 } finally { 2731 StackContext.exit(); 2732 } 2733 2734 return getPoint(sOut, tOut); 2735 } 2736 2737 /** 2738 * Return the intersections between an infinite line and this surface. 2739 * 2740 * @param L0 A point on the line (origin of the line). May not be null. 2741 * @param Ldir The direction vector for the line (does not have to be a unit vector). 2742 * May not be null. 2743 * @param tol Tolerance (in physical space) to refine the point positions to. May not 2744 * be null. 2745 * @return A PointString containing zero or more subrange points made by the 2746 * intersection of this surface with the specified infinite line. If no 2747 * intersection is found, an empty PointString is returned. 2748 */ 2749 @Override 2750 public PointString<SubrangePoint> intersect(GeomPoint L0, GeomVector Ldir, Parameter<Length> tol) { 2751 requireNonNull(L0); 2752 requireNonNull(Ldir); 2753 requireNonNull(tol); 2754 PointString<SubrangePoint> output = PointString.newInstance(); 2755 2756 // First check to see if the line even comes near the surface. 2757 if (!GeomUtil.lineBoundsIntersect(L0, Ldir, this, null)) 2758 return output; 2759 2760 // Make sure the line and curve have the same dimensions. 2761 int numDims = getPhyDimension(); 2762 if (L0.getPhyDimension() > numDims) 2763 throw new IllegalArgumentException(MessageFormat.format( 2764 RESOURCES.getString("sameOrFewerDimErr"), "line", "surface")); 2765 else if (L0.getPhyDimension() < numDims) 2766 L0 = L0.toDimension(numDims); 2767 if (Ldir.getPhyDimension() > numDims) 2768 throw new IllegalArgumentException(MessageFormat.format( 2769 RESOURCES.getString("sameOrFewerDimErr"), "line", "surface")); 2770 else if (Ldir.getPhyDimension() < numDims) 2771 Ldir = Ldir.toDimension(numDims); 2772 2773 // Convert the inputs to the units of this surface and to the dimension of this surface. 2774 Unit<Length> unit = getUnit(); 2775 L0 = L0.toDimension(numDims).to(unit); 2776 GeomVector Lu = Ldir.toDimension(numDims); 2777 2778 // Make sure that the line direction vector is NOT dimensionless. 2779 if (Dimensionless.UNIT.equals(Lu.getUnit())) { 2780 Lu = Vector.valueOf(Lu.toFloat64Vector(), unit); 2781 } else { 2782 Lu = (GeomVector)Lu.to(unit); 2783 } 2784 if (Lu == Ldir) { 2785 Lu = Ldir.copy(); 2786 } 2787 Lu.setOrigin(Point.newInstance(numDims, unit)); 2788 tol = tol.to(unit); 2789 2790 // Look for intersections using recursive subdivision combined with ND root finding. 2791 LineSrfIntersect intersector = LineSrfIntersect.newInstance(this, L0, Lu, tol, output); 2792 output = intersector.intersect(); 2793 LineSrfIntersect.recycle(intersector); 2794 2795 return output; 2796 } 2797 2798 /** 2799 * Return the intersections between a line segment and this surface. 2800 * 2801 * @param line A line segment to intersect with this surface. May not be null. 2802 * @param tol Tolerance (in physical space) to refine the point positions to. May not 2803 * be null. 2804 * @return A list containing two PointString instances each containing zero or more 2805 * subrange points, on this surface and the input line segment respectively, 2806 * made by the intersection of this surface with the specified line segment. 2807 * If no intersections are found a list of two empty PointStrings are 2808 * returned. 2809 */ 2810 @Override 2811 public GeomList<PointString<SubrangePoint>> intersect(LineSegment line, Parameter<Length> tol) { 2812 requireNonNull(line); 2813 requireNonNull(tol); 2814 PointString<SubrangePoint> thisOutput = PointString.newInstance(); 2815 PointString<SubrangePoint> thatOutput = PointString.newInstance(); 2816 GeomList<PointString<SubrangePoint>> output = GeomList.newInstance(); 2817 output.add(thisOutput); 2818 output.add(thatOutput); 2819 2820 // First check to see if the line even comes near the surface. 2821 if (!GeomUtil.lineSegBoundsIntersect(line.getStart(), line.getEnd(), this)) 2822 return output; 2823 2824 // Start by intersecting an infinite line with the surface. 2825 GeomVector Lu = line.getDerivativeVector(); 2826 PointString<SubrangePoint> potentials = intersect(line.getStart(), Lu, tol); 2827 2828 // Now throw out any points that are more than "tol" distance beyond the 2829 // ends of the line. 2830 int size = potentials.size(); 2831 for (int i = 0; i < size; ++i) { 2832 SubrangePoint pnt = potentials.get(i); 2833 Parameter<Length> distance = GeomUtil.pointLineSegDistance(pnt, line.getStart(), line.getEnd()); 2834 if (!distance.isLargerThan(tol)) { 2835 thisOutput.add(pnt); 2836 thatOutput.add(line.getClosest(pnt, GTOL)); 2837 } 2838 } 2839 2840 // Clean up before leaving. 2841 PointString.recycle(potentials); 2842 2843 return output; 2844 } 2845 2846 /** 2847 * Return the intersections between a curve and this surface. 2848 * 2849 * @param curve The curve to intersect with this surface. May not be null. 2850 * @param tol Tolerance (in physical space) to refine the point positions to. May 2851 * not be null. 2852 * @return A list containing two PointString instances each containing zero or more 2853 * subrange points, on this surface and the input curve respectively, made by 2854 * the intersection of this surface with the specified curve. If no 2855 * intersections are found a list of two empty PointStrings are returned. 2856 */ 2857 @Override 2858 public GeomList<PointString<SubrangePoint>> intersect(Curve curve, Parameter<Length> tol) { 2859 GeomList<PointString<SubrangePoint>> ints = curve.intersect(this, requireNonNull(tol)); 2860 GeomList<PointString<SubrangePoint>> output = ints.reverse(); 2861 GeomList.recycle(ints); 2862 return output; 2863 } 2864 2865 /** 2866 * Return the intersections between an infinite plane and this surface. 2867 * 2868 * @param plane The infinite plane to intersect with this surface. May not be null. 2869 * @param tol Tolerance (in physical space) to refine the point positions to. May not be null. 2870 * @return A PointString containing zero or more subrange points made by the 2871 * intersection of this surface with the specified infinite plane. If no 2872 * intersection is found, an empty PointString is returned. 2873 */ 2874 @Override 2875 public GeomList<SubrangeCurve> intersect(GeomPlane plane, Parameter<Length> tol) { 2876 requireNonNull(tol); 2877 GeomList<SubrangeCurve> output = GeomList.newInstance(); 2878 2879 // First check to see if the line even comes near the surface. 2880 if (!GeomUtil.planeBoundsIntersect(requireNonNull(plane), this)) 2881 return output; 2882 2883 // Make sure the plane has the same dimensions as this surface. 2884 int numDims = getPhyDimension(); 2885 if (plane.getPhyDimension() > numDims) 2886 throw new IllegalArgumentException(MessageFormat.format( 2887 RESOURCES.getString("sameOrFewerDimErr"), "plane", "surface")); 2888 2889 // Convert the inputs to the units of this surface and to the dimension of this surface. 2890 Unit<Length> unit = getUnit(); 2891 plane = plane.toDimension(numDims).to(unit); 2892 tol = tol.to(unit); 2893 2894 // Look for and trace intersections. 2895 PlaneSrfIntersect intersector = PlaneSrfIntersect.newInstance(this, plane, tol, output); 2896 try { 2897 output = intersector.intersect(); 2898 2899 } catch (RootException err) { 2900 Logger.getLogger(AbstractSurface.class.getName()).log(Level.WARNING, 2901 "Failed to find intersection with a plane.", err); 2902 2903 } finally { 2904 PlaneSrfIntersect.recycle(intersector); 2905 } 2906 2907 return output; 2908 } 2909 2910 /** 2911 * Return the intersections between another surface and this surface. 2912 * 2913 * @param surface A surface to intersect with this surface. May not be null. 2914 * @param tol Tolerance (in physical space) to refine the point positions to. May 2915 * not be null. 2916 * @return A list containing two lists of SubrangeCurve objects. Each list contains 2917 * zero or more subrange curves, on this surface and the input surface 2918 * respectively, made by the intersection of this surface with the specified 2919 * surface. If no intersections are found a list of two empty lists is 2920 * returned. 2921 */ 2922 @Override 2923 public GeomList<GeomList<SubrangeCurve>> intersect(Surface surface, Parameter<Length> tol) { 2924 requireNonNull(surface); 2925 requireNonNull(tol); 2926 GeomList<SubrangeCurve> thisOutput = GeomList.newInstance(); 2927 GeomList<SubrangeCurve> thatOutput = GeomList.newInstance(); 2928 GeomList<GeomList<SubrangeCurve>> output = GeomList.newInstance(); 2929 output.add(thisOutput); 2930 output.add(thatOutput); 2931 2932 // First check to see if the other surface even comes near this surface. 2933 if (!GeomUtil.boundsIntersect(this, surface)) 2934 return output; 2935 2936 // Make sure the input surface has the same dimensions as this surface. 2937 int numDims = getPhyDimension(); 2938 if (surface.getPhyDimension() > numDims) 2939 throw new IllegalArgumentException(MessageFormat.format( 2940 RESOURCES.getString("sameOrFewerDimErr"), "input surface", "surface")); 2941 2942 // Convert the inputs to the units of this surface and to the dimension of this surface. 2943 Unit<Length> unit = getUnit(); 2944 surface = surface.toDimension(numDims).to(unit); 2945 tol = tol.to(unit); 2946 2947 // Look for and trace intersections. 2948 SrfSrfIntersect intersector = SrfSrfIntersect.newInstance(this, 2949 (AbstractSurface)surface, tol, thisOutput, thatOutput); 2950 try { 2951 intersector.intersect(); 2952 2953 } catch (RootException err) { 2954 Logger.getLogger(AbstractSurface.class.getName()).log(Level.WARNING, 2955 "Failed to find intersection with a surface.", err); 2956 2957 } finally { 2958 SrfSrfIntersect.recycle(intersector); 2959 } 2960 2961 return output; 2962 } 2963 2964 /** 2965 * Return a list or lists containing parameters at the start of logical patches of 2966 * this surface in each parametric direction. The first list contains the parameters 2967 * in the "S" direction and the 2nd in the "T" direction. The first element in each 2968 * list must always be 0.0 and the last element 1.0. These should be good places to 2969 * break the surface up into pieces for analysis, but otherwise are arbitrary. 2970 * Subclasses should override this method to provide better patch starting parameters. 2971 * <p> 2972 * The default implementation always splits the patch into four pieces and returns the 2973 * following parameters [0, 0.5, 1],[0, 0.50, 1]. 2974 * </p> 2975 * 2976 * @return A list of two lists containing parameters at the start of logical patches 2977 * of this surface in each parametric direction. 2978 */ 2979 protected FastTable<FastTable<Double>> getPatchParameters() { 2980 FastTable<FastTable<Double>> tbls = FastTable.newInstance(); 2981 FastTable<Double> knots = FastTable.newInstance(); 2982 knots.add(0.0); 2983 knots.add(0.5); 2984 knots.add(1.0); 2985 tbls.add(knots); 2986 knots = FastTable.newInstance(); 2987 knots.add(0.0); 2988 knots.add(0.5); 2989 knots.add(1.0); 2990 tbls.add(knots); 2991 return tbls; 2992 } 2993 2994 /** 2995 * Return the number of points per patch in the "S" direction that should be used when 2996 * analyzing surface patches by certain methods. Subclasses should override this 2997 * method to provide a more appropriate number of points per patch in the "S" 2998 * direction. 2999 * <p> 3000 * The default implementation always returns 8. 3001 * </p> 3002 * 3003 * @return The number of points per patch in the S direction to use for analysis. 3004 */ 3005 protected int getNumPointsPerPatchS() { 3006 return 8; // 2*(3 + 1) 3007 } 3008 3009 /** 3010 * Return the number of points per patch in the "T" direction that should be used when 3011 * analyzing surface patches by certain methods. Subclasses should override this 3012 * method to provide a more appropriate number of points per patch in the "T" 3013 * direction. 3014 * <p> 3015 * The default implementation always returns 8. 3016 * </p> 3017 * 3018 * @return The number of points per patch in the T direction to use for analysis. 3019 */ 3020 protected int getNumPointsPerPatchT() { 3021 return 8; // 2*(3 + 1) 3022 } 3023 3024 /** 3025 * Return <code>true</code> if this surface is degenerate (i.e.: has area less than 3026 * the specified tolerance squared). 3027 * 3028 * @param tol The tolerance for determining if this surface is degenerate. May not be null. 3029 */ 3030 @Override 3031 public boolean isDegenerate(Parameter<Length> tol) { 3032 StackContext.enter(); 3033 try { 3034 Parameter<Area> area = getArea(GTOL * 100); 3035 return area.compareTo(tol.times(tol)) <= 0; 3036 } finally { 3037 StackContext.exit(); 3038 } 3039 } 3040 3041 /** 3042 * Return <code>true</code> if this surface is planar or <code>false</code> if it is 3043 * not. 3044 * 3045 * @param tol The geometric tolerance to use in determining if the surface is planar. 3046 * May not be null. 3047 */ 3048 @Override 3049 public boolean isPlanar(Parameter<Length> tol) { 3050 requireNonNull(tol); 3051 3052 // If the surface is less than 3D, then it must be planar. 3053 int numDims = getPhyDimension(); 3054 if (numDims <= 2) 3055 return true; 3056 3057 StackContext.enter(); 3058 try { 3059 // If the surface is degenerate (zero area; a point or curve), then check edge curves for being planar. 3060 if (isDegenerate(tol)) { 3061 Curve sCrv = getT0Curve(); 3062 Curve tCrv = getS0Curve(); 3063 if (sCrv.isPlanar(tol) && tCrv.isPlanar(tol)) { 3064 return true; 3065 } 3066 } 3067 3068 // Form a plane from the diagonals of the surface. 3069 Point p0 = getRealPoint(0, 0); 3070 Point p1 = getRealPoint(1, 1); 3071 Vector<Length> d1 = p1.minus(p0).toGeomVector(); 3072 p0 = getRealPoint(1, 0); 3073 p1 = getRealPoint(0, 1); 3074 Vector<Length> d2 = p1.minus(p0).toGeomVector(); 3075 Vector<Dimensionless> nhat = d1.cross(d2).toUnitVector(); 3076 3077 // Get the number of points per patch in the S & T directions. 3078 int nPntsS = getNumPointsPerPatchS(); 3079 int nPntsT = getNumPointsPerPatchT(); 3080 3081 // Get the locations of logical patches on this surface. 3082 FastTable<FastTable<Double>> pParams = getPatchParameters(); 3083 FastTable<Double> sParams = pParams.get(0); 3084 FastTable<Double> tParams = pParams.get(1); 3085 3086 // Enrich the patch parameters by the number of points for each patch. 3087 sParams = enrichParamList(sParams, nPntsS); 3088 tParams = enrichParamList(tParams, nPntsT); 3089 3090 // Extract a grid of points 3091 PointArray<SubrangePoint> arr = this.extractGrid(GridRule.PAR, sParams, null, tParams, null); 3092 3093 // The surface is planar if all the points in the extracted grid are planar. 3094 int numStr = arr.size(); 3095 for (int i = 0; i < numStr; ++i) { 3096 PointString<SubrangePoint> str = arr.get(i); 3097 int numPnts = str.size(); 3098 for (int j = 0; j < numPnts; ++j) { 3099 SubrangePoint sPnt = str.get(j); 3100 Point p = sPnt.copyToReal(); 3101 p1 = GeomUtil.pointPlaneClosest(p, p0, nhat); 3102 if (p.distance(p1).isGreaterThan(tol)) 3103 return false; 3104 } 3105 } 3106 3107 return true; 3108 3109 } finally { 3110 StackContext.exit(); 3111 } 3112 } 3113 3114 private static FastTable<Double> enrichParamList(List<Double> oldParams, int numPerSeg) { 3115 // A list of values evenly or linearly spaced between 0 and 1. 3116 FastTable<Double> spacing = (FastTable)GridSpacing.linear(numPerSeg); 3117 spacing.remove(spacing.size() - 1); // Remove the last point to prevent duplication below. 3118 3119 FastTable<Double> newParams = FastTable.newInstance(); 3120 int size = oldParams.size() - 1; 3121 for (int i = 0; i < size; ++i) { 3122 double pi = oldParams.get(i); 3123 double pip1 = oldParams.get(i + 1); 3124 // Interpolate in a set of values between each existing value. 3125 for (Double s : spacing) { 3126 double p = pi + (pip1 - pi) * s; 3127 newParams.add(p); 3128 } 3129 } 3130 newParams.add(1.0); 3131 3132 // Clean up. 3133 FastTable.recycle(spacing); 3134 3135 return newParams; 3136 } 3137 3138 /** 3139 * Return <code>true</code> if this surface is planar or <code>false</code> if it is 3140 * not. 3141 * 3142 * @param epsFT Flatness tolerance (cosine of the angle allowable between normal 3143 * vectors). 3144 * @param epsELT Edge linearity tolerance (cosine of the angle allowable between edge 3145 * vectors). 3146 */ 3147 @Override 3148 public boolean isPlanar(double epsFT, double epsELT) { 3149 3150 StackContext.enter(); 3151 try { 3152 // Check to see if the corner and center normals are parallel. 3153 Vector<Dimensionless> ni = getNormal(0, 0); 3154 Vector<Dimensionless> nj = getNormal(1, 0); 3155 if (ni.dot(nj).getValue() <= epsFT) 3156 return false; 3157 nj = getNormal(1, 1); 3158 if (ni.dot(nj).getValue() <= epsFT) 3159 return false; 3160 nj = getNormal(0, 1); 3161 if (ni.dot(nj).getValue() <= epsFT) 3162 return false; 3163 nj = getNormal(0.5, 0.5); 3164 if (ni.dot(nj).getValue() <= epsFT) 3165 return false; 3166 3167 // Check to see if the edge vectors at each corner are parallel. 3168 Vector<Dimensionless> Ti = getSDerivative(0, 0, 1, false).toUnitVector(); 3169 Vector<Dimensionless> Tj = getSDerivative(1, 0, 1, false).toUnitVector(); 3170 if (Ti.dot(Tj).getValue() < epsELT) 3171 return false; 3172 Ti = getSDerivative(0, 1, 1, false).toUnitVector(); 3173 Tj = getSDerivative(1, 1, 1, false).toUnitVector(); 3174 if (Ti.dot(Tj).getValue() < epsELT) 3175 return false; 3176 Ti = getTDerivative(1, 0, 1, false).toUnitVector(); 3177 Tj = getTDerivative(1, 1, 1, false).toUnitVector(); 3178 if (Ti.dot(Tj).getValue() < epsELT) 3179 return false; 3180 Ti = getTDerivative(0, 1, 1, false).toUnitVector(); 3181 Tj = getTDerivative(0, 0, 1, false).toUnitVector(); 3182 3183 return Ti.dot(Tj).getValue() >= epsELT; 3184 3185 } finally { 3186 StackContext.exit(); 3187 } 3188 } 3189 3190 /** 3191 * Converts the input parametric position on a surface patch with the specified range 3192 * of parametric positions into a parametric position on the parent surface. 3193 * 3194 * @param s The parametric position on the surface patch to be converted. 3195 * @param s0 The starting parametric position of the surface patch on the parent 3196 * surface. 3197 * @param s1 The ending parametric position of the surface patch on the parent 3198 * surface. 3199 * @return The input patch parametric position converted to the parent surface. 3200 */ 3201 protected static double segmentPos2Parent(double s, double s0, double s1) { 3202 double s_out = s0 + (s1 - s0) * s; 3203 if (s_out > 1) { // Deal with roundoff issues. 3204 s_out = 1; 3205 } else if (s_out < 0) { 3206 s_out = 0; 3207 } 3208 return s_out; 3209 } 3210 3211 /** 3212 * Method that returns true if the specified surface is smaller than the specified 3213 * tolerance. 3214 */ 3215 private static boolean isSmall(Surface srf, Parameter<Length> tol) { 3216 StackContext.enter(); 3217 try { 3218 Parameter<Length> diag = srf.getBoundsMax().distance(srf.getBoundsMin()); 3219 return tol.isGreaterThan(diag); 3220 } finally { 3221 StackContext.exit(); 3222 } 3223 } 3224 3225 /** 3226 * Remove any points from the input list that are closer than the given tolerance to 3227 * each other. 3228 */ 3229 private static PointString<SubrangePoint> removeClosePoints(PointString<SubrangePoint> input, Parameter<Length> tol) { 3230 3231 // Remove any redundant points 3232 int size = input.size(); 3233 for (int j = 0; j < size; ++j) { 3234 SubrangePoint p0 = input.get(j); 3235 for (int i = size - 1; i >= 0; --i) { 3236 if (i == j) 3237 continue; 3238 SubrangePoint p1 = input.get(i); 3239 if (p0.distance(p1).isLessThan(tol)) { 3240 // Make sure we have not jumpped a parametric boundary. 3241 Point st0 = p0.getParPosition().immutable(); 3242 double s0 = st0.getValue(0); 3243 double t0 = st0.getValue(1); 3244 Point st1 = p1.getParPosition().immutable(); 3245 double s1 = st1.getValue(0); 3246 double t1 = st1.getValue(1); 3247 if ((s0 > 0.8 && s1 < 0.2) || (s0 < 0.2 && s1 > 0.8)) 3248 continue; 3249 if ((t0 > 0.8 && t1 < 0.2) || (t0 < 0.2 && t1 > 0.8)) 3250 continue; 3251 // Go ahead and remove the point. 3252 input.remove(i); 3253 } 3254 } 3255 size = input.size(); 3256 } 3257 3258 return input; 3259 } 3260 3261 /** 3262 * An Evaulatable1D that returns the area on a surface integrated in one dimension. 3263 * This is used to find the surface area of a surface. This is <code>f(x)</code> in: 3264 * <code>A = int_t1^t2 {f(x)dt} = int_t1^t2{ int_s1^s2{ |pu x pw| ds } dt }</code>. 3265 * This is also a Runnable that integrates the total area of the surface: 3266 * <code>A = int_t1^t2{ int_s1^s2{ |pu x pw| ds } dt }</code> 3267 */ 3268 private static class AreaEvaluatable extends AbstractEvaluatable1D implements Runnable { 3269 3270 private double _value; 3271 private double _t1, _t2; 3272 3273 private double _s1, _s2; 3274 private double _eps; 3275 private AreaEvaluatable2 _subAreaEval; 3276 public int numEvaluations; 3277 3278 public static AreaEvaluatable newInstance(AbstractSurface surface, double s1, double s2, 3279 double t1, double t2, double eps) { 3280 if (surface instanceof GeomTransform) 3281 surface = (AbstractSurface)surface.copyToReal(); 3282 AreaEvaluatable o = FACTORY.object(); 3283 o._s1 = s1; 3284 o._s2 = s2; 3285 o._t1 = t1; 3286 o._t2 = t2; 3287 o._eps = eps; 3288 o._subAreaEval = AreaEvaluatable2.newInstance(surface); 3289 o.numEvaluations = 0; 3290 return o; 3291 } 3292 3293 /** 3294 * Run to find the total area of the surface: 3295 * <code>A = int_t1^t2{ int_s1^s2{ |pu x pw| ds } dt }</code> 3296 */ 3297 @Override 3298 public void run() { 3299 try { 3300 // Integrate to find the area: A = int_t1^t2{ int_s1^s2{ |pu x pw| ds } dt } 3301 _value = Quadrature.adaptLobatto(this, _t1, _t2, _eps); 3302 3303 } catch (IntegratorException | RootException e) { 3304 // Fall back on Gaussian Quadrature. 3305 Logger.getLogger(AbstractSurface.class.getName()).log(Level.WARNING, 3306 "Falling back on Gaussian Quadrature for getArea()."); 3307 try { 3308 _value = Quadrature.gaussLegendre_Wx1N40(this, _t1, _t2); 3309 3310 } catch (RootException err) { 3311 Logger.getLogger(AbstractSurface.class.getName()).log(Level.WARNING, 3312 "Failed to find area of surface patch.", err); 3313 } 3314 } 3315 3316 } 3317 3318 /** 3319 * Returns the total area of the surface in surface area units after a call to 3320 * run(). 3321 * 3322 * @return The total area of the surface in surface area units. 3323 */ 3324 public double getValue() { 3325 return _value; 3326 } 3327 3328 /** 3329 * Returns the inner integral on the surface at a given value of "t": 3330 * <code>f(t) = int_s1^s2{ |pu x pw| ds }</code> 3331 */ 3332 @Override 3333 public double function(double t) throws RootException { 3334 double value = 0; 3335 _subAreaEval.t = t; 3336 _subAreaEval.numEvaluations = 0; 3337 try { 3338 3339 // Integrate value = int_s1^s2{ |pu x pw| ds } 3340 value = Quadrature.adaptLobatto(_subAreaEval, _s1, _s2, _eps / 10); 3341 3342 } catch (IntegratorException | RootException e) { 3343 // Fall back on Gaussian Quadrature. 3344 Logger.getLogger(AbstractSurface.class.getName()).log(Level.WARNING, 3345 "Falling back on Gaussian Quadrature for AreaEvaluatable."); 3346 try { 3347 value = Quadrature.gaussLegendre_Wx1N40(_subAreaEval, _s1, _s2); 3348 3349 } catch (RootException err) { 3350 throw new RootException(err.getMessage()); 3351 } 3352 3353 } finally { 3354 numEvaluations += _subAreaEval.numEvaluations; 3355 } 3356 return value; 3357 } 3358 3359 public static void recycle(AreaEvaluatable instance) { 3360 AreaEvaluatable2.recycle(instance._subAreaEval); 3361 FACTORY.recycle(instance); 3362 } 3363 3364 private AreaEvaluatable() { } 3365 3366 @SuppressWarnings("unchecked") 3367 private static final ObjectFactory<AreaEvaluatable> FACTORY = new ObjectFactory<AreaEvaluatable>() { 3368 @Override 3369 protected AreaEvaluatable create() { 3370 return new AreaEvaluatable(); 3371 } 3372 3373 @Override 3374 protected void cleanup(AreaEvaluatable obj) { 3375 obj._subAreaEval = null; 3376 } 3377 }; 3378 } 3379 3380 /** 3381 * An Evaulatable1D that returns the area of an infinitesimally small portion of the 3382 * surface at (s,t). This is used to find the surface area of a surface. This is 3383 * <code>f(x)</code> in: 3384 * <code>A = int_s1^s2 {f(x)dt} = int_s1^s2{ |pu x pw| ds }</code>. 3385 */ 3386 private static class AreaEvaluatable2 extends AbstractEvaluatable1D { 3387 3388 private AbstractSurface _srf; 3389 private int _dim; 3390 public double t; 3391 public int numEvaluations; 3392 3393 public static AreaEvaluatable2 newInstance(AbstractSurface surface) { 3394 AreaEvaluatable2 o = FACTORY.object(); 3395 o._srf = surface; 3396 o._dim = surface.getPhyDimension(); 3397 o.numEvaluations = 0; 3398 return o; 3399 } 3400 3401 @Override 3402 public double function(double s) throws RootException { 3403 ++numEvaluations; 3404 3405 StackContext.enter(); 3406 try { 3407 // Get the surface derivative in the s direction. 3408 Vector<Length> pu = _srf.getSDerivative(s, t, 1, true); 3409 3410 // Collapsed edges provide no increment in area. 3411 if (pu.magValue() < SQRT_EPS) { 3412 return 0.0; 3413 } 3414 3415 // Gert the surface derivative in the t direction. 3416 Vector<Length> pw = _srf.getTDerivative(s, t, 1, true); 3417 3418 // Collapsed edges provide no increment in area. 3419 if (pw.magValue() < SQRT_EPS) { 3420 return 0.0; 3421 } 3422 3423 if (_dim < 3) { 3424 // 2D is a special case. 3425 double pux = pu.getValue(0), puy = pu.getValue(1); 3426 double pwx = pw.getValue(0), pwy = pw.getValue(1); 3427 double puxpw = pux * pwy - puy * pwx; 3428 return puxpw; 3429 } 3430 3431 // Normal vector is the cross product of pu and pw vectors 3432 // (it is also the area swept out by the two tangent vectors). 3433 Vector n = pu.cross(pw); 3434 3435 // Return the magnitude of the normal vector. 3436 double nmag = n.magValue(); 3437 return nmag; 3438 3439 } finally { 3440 StackContext.exit(); 3441 } 3442 } 3443 3444 public static void recycle(AreaEvaluatable2 instance) { 3445 FACTORY.recycle(instance); 3446 } 3447 3448 private AreaEvaluatable2() { } 3449 3450 @SuppressWarnings("unchecked") 3451 private static final ObjectFactory<AreaEvaluatable2> FACTORY = new ObjectFactory<AreaEvaluatable2>() { 3452 @Override 3453 protected AreaEvaluatable2 create() { 3454 return new AreaEvaluatable2(); 3455 } 3456 3457 @Override 3458 protected void cleanup(AreaEvaluatable2 obj) { 3459 obj._srf = null; 3460 } 3461 }; 3462 } 3463 3464 /** 3465 * An Evaulatable1D that returns the volume on a surface integrated in one dimension. 3466 * This is used to find the enclosed volume of a surface. This is <code>f(x)</code> 3467 * in: 3468 * <code>A = int_t1^t2 {f(x)dt} = int_t1^t2{ int_s1^s2{ 1/3*(p DOT n) ds } dt }</code>. 3469 * This is also a Runnable that integrates the total volume of the surface: 3470 * <code>V = int_t1^t2{ int_s1^s2{ 1/3*(p DOT n) ds } dt</code> 3471 */ 3472 private static class VolumeEvaluatable extends AbstractEvaluatable1D implements Runnable { 3473 3474 private double _value; 3475 private double _t1, _t2; 3476 3477 private double _s1, _s2; 3478 private double _eps; 3479 private VolumeEvaluatable2 _subVolumeEval; 3480 public int numEvaluations; 3481 3482 public static VolumeEvaluatable newInstance(AbstractSurface surface, double s1, double s2, 3483 double t1, double t2, double eps) { 3484 if (surface instanceof GeomTransform) 3485 surface = (AbstractSurface)surface.copyToReal(); 3486 VolumeEvaluatable o = FACTORY.object(); 3487 o._s1 = s1; 3488 o._s2 = s2; 3489 o._t1 = t1; 3490 o._t2 = t2; 3491 o._eps = eps; 3492 o._subVolumeEval = VolumeEvaluatable2.newInstance(surface); 3493 o.numEvaluations = 0; 3494 return o; 3495 } 3496 3497 /** 3498 * Run to find the total volume of the surface: 3499 * <code>V = int_t1^t2{ int_s1^s2{ 1/3*(p DOT n) ds } dt</code> 3500 */ 3501 @Override 3502 public void run() { 3503 try { 3504 // Integrate to find the volume: V = int_t1^t2{ int_s1^s2{ 1/3*(p DOT n) ds } dt 3505 _value = Quadrature.adaptLobatto(this, _t1, _t2, _eps); 3506 3507 } catch (IntegratorException | RootException e) { 3508 // Fall back on Gaussian Quadrature. 3509 Logger.getLogger(AbstractSurface.class.getName()).log(Level.WARNING, 3510 "Fallling back on Gaussian Quadrature for getVolume()."); 3511 try { 3512 _value = Quadrature.gaussLegendre_Wx1N40(this, _t1, _t2); 3513 3514 } catch (RootException err) { 3515 Logger.getLogger(AbstractSurface.class.getName()).log(Level.WARNING, 3516 "Failed to find volume of surface patch.", err); 3517 } 3518 } 3519 3520 } 3521 3522 /** 3523 * Returns the total volume of the surface in surface volume units after a call to 3524 * run(). 3525 */ 3526 public double getValue() { 3527 return _value; 3528 } 3529 3530 @Override 3531 public double function(double t) throws RootException { 3532 double value = 0; 3533 _subVolumeEval.t = t; 3534 _subVolumeEval.numEvaluations = 0; 3535 3536 try { 3537 3538 // Integrate value = 1/3*int_s1^s2{ (p DOT n) ds } 3539 value = Quadrature.adaptLobatto(_subVolumeEval, _s1, _s2, _eps) / 3; 3540 3541 } catch (IntegratorException | RootException e) { 3542 // Fall back on Gaussian Quadrature. 3543 Logger.getLogger(AbstractSurface.class.getName()).log(Level.WARNING, 3544 "Fallling back on Gaussian Quadrature for VolumeEvaluatable."); 3545 try { 3546 value = Quadrature.gaussLegendre_Wx1N40(_subVolumeEval, _s1, _s2) / 3; 3547 3548 } catch (RootException err) { 3549 throw new RootException(err.getMessage()); 3550 } 3551 3552 } finally { 3553 numEvaluations += _subVolumeEval.numEvaluations; 3554 } 3555 3556 return value; 3557 } 3558 3559 public static void recycle(VolumeEvaluatable instance) { 3560 VolumeEvaluatable2.recycle(instance._subVolumeEval); 3561 FACTORY.recycle(instance); 3562 } 3563 3564 private VolumeEvaluatable() { } 3565 3566 @SuppressWarnings("unchecked") 3567 private static final ObjectFactory<VolumeEvaluatable> FACTORY = new ObjectFactory<VolumeEvaluatable>() { 3568 @Override 3569 protected VolumeEvaluatable create() { 3570 return new VolumeEvaluatable(); 3571 } 3572 3573 @Override 3574 protected void cleanup(VolumeEvaluatable obj) { 3575 obj._subVolumeEval = null; 3576 } 3577 }; 3578 } 3579 3580 /** 3581 * An Evaulatable1D that returns the volume of an infinitesimally small portion of the 3582 * surface at (s,t). This is used to find the enclosed volume of a surface. This is 3583 * <code>f(x)</code> in: 3584 * <code>A = int_s1^s2 {f(x)dt} = int_s1^s2{ 1/3*(p DOT n) ds }</code>. 3585 */ 3586 private static class VolumeEvaluatable2 extends AbstractEvaluatable1D { 3587 3588 private AbstractSurface _srf; 3589 public double t; 3590 private int _dim; 3591 public int numEvaluations; 3592 3593 public static VolumeEvaluatable2 newInstance(AbstractSurface surface) { 3594 VolumeEvaluatable2 o = FACTORY.object(); 3595 o._srf = surface; 3596 o._dim = surface.getPhyDimension(); 3597 return o; 3598 } 3599 3600 @Override 3601 public double function(double s) throws RootException { 3602 ++numEvaluations; 3603 3604 StackContext.enter(); 3605 try { 3606 // Get the derivatives in the s-direction. 3607 List<Vector<Length>> ders = _srf.getSDerivatives(s, t, 1, true); 3608 3609 // Get the surface point (it is in the derivative list). 3610 Vector<Length> p = ders.get(0); 3611 if (p.norm().isApproxZero()) { 3612 return 0; 3613 } 3614 3615 // Get the surface derivative in the s direction. 3616 Vector<Length> pu = ders.get(1); 3617 3618 // Collapsed edges provide no increment in volume. 3619 if (pu.magValue() < SQRT_EPS) { 3620 return 0.0; 3621 } 3622 3623 // Get the surface derivative in the t direction. 3624 Vector<Length> pw = _srf.getTDerivative(s, t, 1, true); 3625 3626 // Collapsed edges provide no increment in volume. 3627 if (pw.magValue() < SQRT_EPS) { 3628 return 0.0; 3629 } 3630 3631 // Make sure the vectors are at least 3D. 3632 if (_dim < 3) { 3633 p = p.toDimension(3); 3634 pu = pu.toDimension(3); 3635 pw = pw.toDimension(3); 3636 } 3637 3638 // Normal vector is the cross product of pu and pw vectors 3639 // (it is also the area swept out by the two tangent vectors). 3640 Vector n = pu.cross(pw); 3641 3642 // Calculate three times the volume increment. 3643 double dV = p.dot(n).getValue(); 3644 3645 // Return the 3xVolume increment. 3646 return dV; 3647 3648 } finally { 3649 StackContext.exit(); 3650 } 3651 } 3652 3653 public static void recycle(VolumeEvaluatable2 instance) { 3654 FACTORY.recycle(instance); 3655 } 3656 3657 private VolumeEvaluatable2() { } 3658 3659 @SuppressWarnings("unchecked") 3660 private static final ObjectFactory<VolumeEvaluatable2> FACTORY = new ObjectFactory<VolumeEvaluatable2>() { 3661 @Override 3662 protected VolumeEvaluatable2 create() { 3663 return new VolumeEvaluatable2(); 3664 } 3665 3666 @Override 3667 protected void cleanup(VolumeEvaluatable2 obj) { 3668 obj._srf = null; 3669 } 3670 }; 3671 } 3672 3673 /** 3674 * A ScalarFunctionND that calculates the distance squared from a target point to 3675 * points on the specified surface at s,t. This is intended to only be used by 3676 * AbstractSurface for finding the closest and farthest points. 3677 */ 3678 private static class PointSurfaceDistFunction implements ScalarFunctionND { 3679 3680 private Surface _srf; 3681 private Point _point; 3682 private boolean _closest; 3683 public int numEvaluations; 3684 3685 public static PointSurfaceDistFunction newInstance(Surface surface, GeomPoint point, boolean closest) { 3686 if (surface instanceof GeomTransform) 3687 surface = (Surface)surface.copyToReal(); 3688 PointSurfaceDistFunction o = FACTORY.object(); 3689 o._srf = surface; 3690 o._point = point.immutable(); 3691 o._closest = closest; 3692 o.numEvaluations = 0; 3693 return o; 3694 } 3695 3696 @Override 3697 public double function(double x[]) throws RootException { 3698 double s = x[0]; 3699 double t = x[1]; 3700 3701 // Deal with the minimizer going out of bounds of the surface parameterization 3702 // by penalizing such points. 3703 double fs = 1, ft = 1; 3704 if (s < 0) { 3705 fs = 1 - s; 3706 s = 0; 3707 } 3708 if (s > 1) { 3709 fs = fs * s; 3710 s = 1; 3711 } 3712 if (t < 0) { 3713 ft = 1 - t; 3714 t = 0; 3715 } 3716 if (t > 1) { 3717 ft = ft * t; 3718 t = 1; 3719 } 3720 fs *= fs; 3721 ft *= ft; 3722 3723 StackContext.enter(); 3724 try { 3725 // Get the point on the surface at s,t. 3726 Point p = _srf.getRealPoint(s, t); 3727 3728 // Calculate the distance to the surface point from our target point. 3729 double dist = _point.distanceValue(p); 3730 3731 if (_closest) 3732 // Penalize the distance if the input point is out of bounds of the surface. 3733 dist *= fs * ft; 3734 else { 3735 // Penalize the distance if the input point is out of bounds of the surface. 3736 dist /= fs * ft; 3737 3738 // Change sign on distance to find the farthest point. 3739 dist *= -1; 3740 } 3741 3742 ++numEvaluations; 3743 return dist; 3744 3745 } finally { 3746 StackContext.exit(); 3747 } 3748 } 3749 3750 @Override 3751 public boolean derivatives(double[] x, double[] dydx) { 3752 return false; 3753 } 3754 3755 public static void recycle(PointSurfaceDistFunction instance) { 3756 FACTORY.recycle(instance); 3757 } 3758 3759 private PointSurfaceDistFunction() { } 3760 3761 @SuppressWarnings("unchecked") 3762 private static final ObjectFactory<PointSurfaceDistFunction> FACTORY = new ObjectFactory<PointSurfaceDistFunction>() { 3763 @Override 3764 protected PointSurfaceDistFunction create() { 3765 return new PointSurfaceDistFunction(); 3766 } 3767 3768 @Override 3769 protected void cleanup(PointSurfaceDistFunction obj) { 3770 obj._srf = null; 3771 obj._point = null; 3772 } 3773 }; 3774 } 3775 3776 /** 3777 * A ScalarFunctionND that calculates the distance squared from a target point to 3778 * points on the specified surface at s,t. This is intended to only be used by 3779 * AbstractSurface for finding the closest and farthest points. 3780 */ 3781 private static class SrfSrfDistFunction implements ScalarFunctionND { 3782 3783 private Surface _thisSrf; 3784 private Surface _thatSrf; 3785 private double _tol; 3786 private boolean _closest; 3787 public int numEvaluations; 3788 3789 public static SrfSrfDistFunction newInstance(Surface thisSrf, Surface thatSrf, double tol, boolean closest) { 3790 if (thisSrf instanceof GeomTransform) 3791 thisSrf = (Surface)thisSrf.copyToReal(); 3792 if (thatSrf instanceof GeomTransform) 3793 thatSrf = (Surface)thatSrf.copyToReal(); 3794 SrfSrfDistFunction o = FACTORY.object(); 3795 o._thisSrf = thisSrf; 3796 o._thatSrf = thatSrf; 3797 o._tol = tol; 3798 o._closest = closest; 3799 o.numEvaluations = 0; 3800 return o; 3801 } 3802 3803 @Override 3804 public double function(double x[]) throws RootException { 3805 double s = x[0]; 3806 double t = x[1]; 3807 3808 // Deal with the minimizer going out of bounds of the surface parameterization 3809 // by penalizing such points. 3810 double fs = 1, ft = 1; 3811 if (s < 0) { 3812 fs = 1 - s; 3813 s = 0; 3814 } 3815 if (s > 1) { 3816 fs = fs * s; 3817 s = 1; 3818 } 3819 if (t < 0) { 3820 ft = 1 - t; 3821 t = 0; 3822 } 3823 if (t > 1) { 3824 ft = ft * t; 3825 t = 1; 3826 } 3827 fs *= fs; 3828 ft *= ft; 3829 3830 StackContext.enter(); 3831 try { 3832 // Get the point on the surface at s,t. 3833 Point p = _thisSrf.getRealPoint(s, t); 3834 3835 // Get the closest point on the remote surface. 3836 SubrangePoint q = _thatSrf.getClosest(p, _tol); 3837 3838 // Calculate the distance to the surface point from our target point. 3839 double dist = q.distanceValue(p); 3840 3841 if (_closest) 3842 // Penalize the distance if the input point is out of bounds of the surface. 3843 dist *= fs * ft; 3844 else { 3845 // Penalize the distance if the input point is out of bounds of the surface. 3846 dist /= fs * ft; 3847 3848 // Change sign on distance to find the farthest point. 3849 dist *= -1; 3850 } 3851 3852 ++numEvaluations; 3853 return dist; 3854 3855 } finally { 3856 StackContext.exit(); 3857 } 3858 } 3859 3860 @Override 3861 public boolean derivatives(double[] x, double[] dydx) { 3862 return false; 3863 } 3864 3865 public static void recycle(SrfSrfDistFunction instance) { 3866 FACTORY.recycle(instance); 3867 } 3868 3869 private SrfSrfDistFunction() { } 3870 3871 @SuppressWarnings("unchecked") 3872 private static final ObjectFactory<SrfSrfDistFunction> FACTORY = new ObjectFactory<SrfSrfDistFunction>() { 3873 @Override 3874 protected SrfSrfDistFunction create() { 3875 return new SrfSrfDistFunction(); 3876 } 3877 3878 @Override 3879 protected void cleanup(SrfSrfDistFunction obj) { 3880 obj._thisSrf = null; 3881 obj._thatSrf = null; 3882 } 3883 }; 3884 } 3885 3886 /** 3887 * A ScalarFunctionND that calculates the distance squared from a target plane to 3888 * points on the specified surface at s,t. This is intended to only be used by 3889 * AbstractSurface for finding the closest and farthest points. 3890 */ 3891 private static class SrfPlaneDistFunction implements ScalarFunctionND { 3892 3893 private Surface _thisSrf; 3894 private GeomPlane _thatPlane; 3895 private boolean _closest; 3896 public int numEvaluations; 3897 3898 public static SrfPlaneDistFunction newInstance(Surface thisSrf, GeomPlane thatPlane, boolean closest) { 3899 if (thisSrf instanceof GeomTransform) 3900 thisSrf = (Surface)thisSrf.copyToReal(); 3901 if (thatPlane instanceof GeomTransform) 3902 thatPlane = thatPlane.copyToReal(); 3903 SrfPlaneDistFunction o = FACTORY.object(); 3904 o._thisSrf = thisSrf; 3905 o._thatPlane = thatPlane; 3906 o._closest = closest; 3907 o.numEvaluations = 0; 3908 return o; 3909 } 3910 3911 @Override 3912 public double function(double x[]) throws RootException { 3913 double s = x[0]; 3914 double t = x[1]; 3915 3916 // Deal with the minimizer going out of bounds of the surface parameterization 3917 // by penalizing such points. 3918 double fs = 1, ft = 1; 3919 if (s < 0) { 3920 fs = 1 - s; 3921 s = 0; 3922 } 3923 if (s > 1) { 3924 fs = fs * s; 3925 s = 1; 3926 } 3927 if (t < 0) { 3928 ft = 1 - t; 3929 t = 0; 3930 } 3931 if (t > 1) { 3932 ft = ft * t; 3933 t = 1; 3934 } 3935 fs *= fs; 3936 ft *= ft; 3937 3938 StackContext.enter(); 3939 try { 3940 // Get the point on the surface at s,t. 3941 Point p = _thisSrf.getRealPoint(s, t); 3942 3943 // Get the closest point on the remote plane. 3944 Point q = _thatPlane.getClosest(p); 3945 3946 // Calculate the distance to the surface point from our target point. 3947 double dist = q.distanceValue(p); 3948 3949 if (_closest) 3950 // Penalize the distance if the input point is out of bounds of the surface. 3951 dist *= fs * ft; 3952 else { 3953 // Penalize the distance if the input point is out of bounds of the surface. 3954 dist /= fs * ft; 3955 3956 // Change sign on distance to find the farthest point. 3957 dist *= -1; 3958 } 3959 3960 ++numEvaluations; 3961 return dist; 3962 3963 } finally { 3964 StackContext.exit(); 3965 } 3966 } 3967 3968 @Override 3969 public boolean derivatives(double[] x, double[] dydx) { 3970 return false; 3971 } 3972 3973 public static void recycle(SrfPlaneDistFunction instance) { 3974 FACTORY.recycle(instance); 3975 } 3976 3977 private SrfPlaneDistFunction() { } 3978 3979 @SuppressWarnings("unchecked") 3980 private static final ObjectFactory<SrfPlaneDistFunction> FACTORY = new ObjectFactory<SrfPlaneDistFunction>() { 3981 @Override 3982 protected SrfPlaneDistFunction create() { 3983 return new SrfPlaneDistFunction(); 3984 } 3985 3986 @Override 3987 protected void cleanup(SrfPlaneDistFunction obj) { 3988 obj._thisSrf = null; 3989 obj._thatPlane = null; 3990 } 3991 }; 3992 } 3993 3994 /** 3995 * A class used to intersect an line segment and this surface. 3996 */ 3997 private static class LineSrfIntersect { 3998 3999 private static final double EPS_FT = 0.9848; // Planar flatness angle tolerance : cos(10 deg) 4000 private static final int MAXDEPTH = 15; 4001 4002 private AbstractSurface _thisSrf; 4003 private Point _L0; 4004 private Vector<Length> _L0v; 4005 private Vector<Length> _Ldir; 4006 private Parameter<Length> _tol; 4007 private PointString<SubrangePoint> _output; 4008 private PointString<SubrangePoint> _approxPnts; 4009 private int _numDims; 4010 4011 @SuppressWarnings("null") 4012 public static LineSrfIntersect newInstance(AbstractSurface thisSrf, 4013 GeomPoint L0, GeomVector<Length> Ldir, Parameter<Length> tol, PointString<SubrangePoint> output) { 4014 4015 if (thisSrf instanceof GeomTransform) 4016 thisSrf = (AbstractSurface)thisSrf.copyToReal(); 4017 4018 // Create an instance of this class and fill in the class variables. 4019 LineSrfIntersect o = FACTORY.object(); 4020 o._thisSrf = thisSrf; 4021 o._L0 = L0.immutable(); 4022 o._L0v = L0.toGeomVector(); 4023 o._Ldir = Ldir.immutable(); 4024 o._tol = tol; 4025 o._output = output; 4026 o._approxPnts = PointString.newInstance(); 4027 o._numDims = thisSrf.getPhyDimension(); 4028 4029 return o; 4030 } 4031 4032 /** 4033 * Method that actually carries out the intersection adding the intersection 4034 * points to the list provided in the constructor. 4035 * 4036 * @return A list containing a PointString instance containing zero or more 4037 * subrange points, on this surface, made by the intersection of this 4038 * surface with the specified line. If no intersections are found an empty 4039 * PointString is returned. 4040 */ 4041 public PointString<SubrangePoint> intersect() { 4042 4043 // Use a divide and conquer approach to approximating the intersection points. 4044 // This method will add approximate intersection points to "_approxPnts". 4045 divideAndConquerLine(1, _thisSrf, 0, 0, 1, 1); 4046 4047 // Refine each approximate intersection. 4048 try { 4049 if (_numDims > 3) { 4050 // Method that works for high dimensions, but is slow. 4051 lineSrfIntersectHighDim(_approxPnts); 4052 4053 } else { 4054 // A much faster method is available for 3D. 4055 lineSrfIntersect3D(_approxPnts); 4056 } 4057 } catch (RootException err) { 4058 Logger.getLogger(AbstractSurface.class.getName()).log(Level.WARNING, 4059 "Failed to refine line-Surface intersection.", err); 4060 4061 } finally { 4062 PointString.recycle(_approxPnts); 4063 _approxPnts = null; 4064 } 4065 4066 // Remove any duplicate points. 4067 for (int i = 0; i < _output.size(); ++i) { 4068 SubrangePoint pi = _output.get(i); 4069 for (int j = _output.size() - 1; j >= 0; --j) { 4070 if (i != j) { 4071 SubrangePoint pj = _output.get(j); 4072 if (pi.isApproxEqual(pj, _tol)) { 4073 _output.remove(j); 4074 --j; 4075 } 4076 } 4077 } 4078 } 4079 4080 return _output; 4081 } 4082 4083 /** 4084 * Uses a recursive "Divide and Conquer" approach to intersecting a surface patch 4085 * with a line. On each call, the following happens: 4086 * <pre> 4087 * 1) The patch is checked to see if it is approx. a flat plane. 4088 * 1b) If it is, then a line-plane intersection is performed, the 4089 * approximate intersection point added to the _approxPnts list and the method exited. 4090 * 2) The patch is divided into four quadrant patches. 4091 * 2b) Each patch is tested to see if it's bounding box is intersected by the line. 4092 * If it is, then this method is called with that patch recursively. 4093 * Otherwise, the method is exited. 4094 * </pre> 4095 * A number of class variables are used to pass information to this recursion: 4096 * <pre> 4097 * _thisSrf The full surfaces intersections are being found for. 4098 * _L0 A point at the start of the line being intersected with this surface. 4099 * _Ldir A direction vector for the line being intersected with this surface. 4100 * _tol The tolerance to use in determining if the geometry is in tolerance. 4101 * _approxPnts A list used to collect the approximate subrange intersection points. 4102 * </pre> 4103 * 4104 * @param patch The surface patch being checked for intersection with the line. 4105 * @param s0 The "s" parametric position of the start of the patch on the 4106 * master "this" surface. 4107 * @param t0 The "t" parametric position of the start of the patch on the 4108 * master "this" surface. 4109 * @param s1 The parametric "s" position of the end of the patch on the master 4110 * "this" surface. 4111 * @param t1 The parametric "t" position of the end of the patch on the master 4112 * "this" surface. 4113 */ 4114 private void divideAndConquerLine(int depth, Surface patch, double s0, double t0, double s1, double t1) { 4115 4116 // Check to see if this patch is nearly planar. 4117 if (depth > MAXDEPTH || patch.isPlanar(EPS_FT, EPS_FT)) { // Using a fast method that unfortunately ignores "tol". 4118 4119 // Form a single quadrilateral panel from the corners of the patch. 4120 Point A = patch.getRealPoint(0, 0); 4121 Point B = patch.getRealPoint(0, 1); 4122 Point C = patch.getRealPoint(1, 1); 4123 Point D = patch.getRealPoint(1, 0); 4124 4125 // Calculate a quad normal vector. 4126 Vector<Dimensionless> nhat = GeomUtil.quadNormal(A, B, C, D); 4127 4128 // Does the input line intersect this quadrilateral panel? 4129 boolean intersect = GeomUtil.lineQuadIntersect(_L0, _Ldir, A, B, C, D, nhat, null); 4130 if (intersect) { 4131 // Store a corner of the patch as the approximate intersection point. 4132 _approxPnts.add(_thisSrf.getPoint(s0, t0)); 4133 } 4134 4135 } else { 4136 // Subdivide the patch into 4 quadrants. 4137 GeomList<Surface> leftRight = patch.splitAtS(0.5); 4138 GeomList<Surface> botTop = leftRight.get(0).splitAtT(0.5); 4139 Surface p00 = botTop.get(0); // Lower Left 4140 Surface p01 = botTop.get(1); // Top Left 4141 GeomList.recycle(botTop); 4142 botTop = leftRight.get(1).splitAtT(0.5); 4143 Surface p10 = botTop.get(0); // Lower Right 4144 Surface p11 = botTop.get(1); // Top Right 4145 GeomList.recycle(botTop); 4146 GeomList.recycle(leftRight); 4147 4148 // Check for possible intersections on the lower-left patch. 4149 if (GeomUtil.lineBoundsIntersect(_L0, _Ldir, p00, null)) { 4150 // May be an intersection. 4151 double sl = s0; 4152 double sh = 0.5 * (s0 + s1); 4153 double tl = t0; 4154 double th = 0.5 * (t0 + t1); 4155 4156 // Recurse down to a finer level of detail. 4157 divideAndConquerLine(depth + 1, p00, sl, tl, sh, th); 4158 } 4159 if (p00 instanceof BasicNurbsSurface) 4160 BasicNurbsSurface.recycle((BasicNurbsSurface)p00); 4161 4162 // Check for possible intersections on the lower-right patch. 4163 if (GeomUtil.lineBoundsIntersect(_L0, _Ldir, p10, null)) { 4164 // May be an intersection. 4165 double sl = 0.5 * (s0 + s1); 4166 double sh = s1; 4167 double tl = t0; 4168 double th = 0.5 * (t0 + t1); 4169 4170 // Recurse down to a finer level of detail. 4171 divideAndConquerLine(depth + 1, p10, sl, tl, sh, th); 4172 } 4173 if (p10 instanceof BasicNurbsSurface) 4174 BasicNurbsSurface.recycle((BasicNurbsSurface)p10); 4175 4176 // Check for possible intersections on the top-left patch. 4177 if (GeomUtil.lineBoundsIntersect(_L0, _Ldir, p01, null)) { 4178 // May be an intersection. 4179 double sl = s0; 4180 double sh = 0.5 * (s0 + s1); 4181 double tl = 0.5 * (t0 + t1); 4182 double th = t1; 4183 4184 // Recurse down to a finer level of detail. 4185 divideAndConquerLine(depth + 1, p01, sl, tl, sh, th); 4186 } 4187 if (p01 instanceof BasicNurbsSurface) 4188 BasicNurbsSurface.recycle((BasicNurbsSurface)p01); 4189 4190 // Check for possible intersections on the top-right patch. 4191 if (GeomUtil.lineBoundsIntersect(_L0, _Ldir, p11, null)) { 4192 // May be an intersection. 4193 double sl = 0.5 * (s0 + s1); 4194 double sh = s1; 4195 double tl = 0.5 * (t0 + t1); 4196 double th = t1; 4197 4198 // Recurse down to a finer level of detail. 4199 divideAndConquerLine(depth + 1, p11, sl, tl, sh, th); 4200 } 4201 if (p11 instanceof BasicNurbsSurface) 4202 BasicNurbsSurface.recycle((BasicNurbsSurface)p11); 4203 4204 } 4205 4206 } 4207 4208 /** 4209 * This high physical dimension method refines the supplied approximate 4210 * line-surface intersection to be an accurate intersection point stored in 4211 * _output. 4212 * 4213 * @param approxPnts List of approximate intersection points. 4214 */ 4215 private void lineSrfIntersectHighDim(PointString<SubrangePoint> approxPnts) throws RootException { 4216 4217 // Create a list of parametric position points. 4218 Point[] parPnts = Point.allocateArray(approxPnts.size()); 4219 int numParPnts = 0; 4220 4221 StackContext.enter(); 4222 try { 4223 int size = approxPnts.size(); 4224 for (int i = 0; i < size; ++i) { 4225 SubrangePoint aPnt = approxPnts.get(i); 4226 4227 // Get a 1st guess at "u" (parametric distance along line. 4228 Vector<Length> W = aPnt.toGeomVector().minus(_L0v); 4229 double u0 = _Ldir.dot(W).getValue(); 4230 4231 // Create a function that is used to iteratively find the intersection point. 4232 LineSrfIntersectNDFun evalFun = LineSrfIntersectNDFun.newInstance(_thisSrf, _L0, _Ldir, aPnt); 4233 4234 // Find the intersection point. 4235 double u = Roots.findRoot1D(evalFun, u0 * 0.98, u0 * 1.02, GTOL); 4236 4237 // Turn the parametric position along the line to a point on the surface. 4238 Point dQu = Point.valueOf(_Ldir.times(u)).plus(_Ldir.getOrigin()); 4239 Point Qu = _L0.plus(dQu); 4240 double nearS = aPnt.getParPosition().getValue(0); 4241 double nearT = aPnt.getParPosition().getValue(1); 4242 SubrangePoint Pst = _thisSrf.getClosest(Qu, nearS, nearT, GTOL); 4243 4244 // Store parametric position of intersection for later. 4245 Point p = Pst.getParPosition().immutable(); 4246 parPnts[numParPnts++] = StackContext.outerCopy(p); 4247 } 4248 } finally { 4249 StackContext.exit(); 4250 } 4251 4252 // Convert list of parametric positions into subrange points. 4253 for (int i = 0; i < numParPnts; ++i) { 4254 Point parPnt = parPnts[i]; 4255 SubrangePoint Pst = _thisSrf.getPoint(parPnt); 4256 _output.add(Pst); 4257 } 4258 4259 Point.recycleArray(parPnts); 4260 } 4261 4262 /** 4263 * This 3D method refines the supplied approximate line-surface intersections to 4264 * be accurate intersection points stored in _output. This method ONLY works if 4265 * the geometry is 3D. 4266 * 4267 * @param approxPnts List of approximate intersection points. 4268 */ 4269 private void lineSrfIntersect3D(PointString<SubrangePoint> approxPnts) throws RootException { 4270 4271 // Create a list of parametric position points. 4272 Point[] parPnts = Point.allocateArray(approxPnts.size()); 4273 int numParPnts = 0; 4274 4275 StackContext.enter(); 4276 try { 4277 // Re-orient the geometry so that the the intersecting line is aligned with the +Z axis. 4278 GTransform T = GTransform.newTranslation(_L0v.opposite()); // Translate so origin of line is at origin. 4279 Vector zAxis = Vector.valueOf(0, 0, 1); 4280 T = GTransform.valueOf(_Ldir.toUnitVector(), zAxis).times(T); // Align with Z axis. 4281 AbstractSurface thisSrf = (AbstractSurface)_thisSrf.getTransformed(T).copyToReal(); 4282 4283 // Create a function that is used to iteratively find the intersection points. 4284 LineSrfIntersect3DFun evalFun = LineSrfIntersect3DFun.newInstance(thisSrf); 4285 4286 // Loop over all the approximate intersections and refine them each. 4287 double[] x = ArrayFactory.DOUBLES_FACTORY.array(2); 4288 int size = approxPnts.size(); 4289 for (int i = 0; i < size; ++i) { 4290 SubrangePoint aPnt = approxPnts.get(i); 4291 4292 // Get the parametric position of a point near a solution. 4293 double s = aPnt.getParPosition().getValue(0); 4294 double t = aPnt.getParPosition().getValue(1); 4295 x[0] = s; 4296 x[1] = t; // Guess at intersection point. 4297 4298 // Find the intersection point using an N-dimensional root finder. 4299 try { 4300 boolean check = Roots.findRootsND(evalFun, x, 2); 4301 if (!check) { 4302 // Store the intersection point found. 4303 s = x[0]; 4304 t = x[1]; // Actual intersection point. 4305 s = (s > 1 ? 1 : s < 0 ? 0 : s); 4306 t = (t > 1 ? 1 : t < 0 ? 0 : t); 4307 4308 // Store parametric position of intersection for later. 4309 Point p = Point.valueOf(s, t); 4310 parPnts[numParPnts++] = StackContext.outerCopy(p); 4311 } 4312 } catch (RootException e) { 4313 /* Do nothing: Go on to the next approximate point. */ 4314 } 4315 } 4316 4317 } finally { 4318 StackContext.exit(); 4319 } 4320 4321 // Convert list of parametric positions into subrange points. 4322 for (int i = 0; i < numParPnts; ++i) { 4323 Point parPnt = parPnts[i]; 4324 SubrangePoint Pst = _thisSrf.getPoint(parPnt); 4325 _output.add(Pst); 4326 } 4327 4328 Point.recycleArray(parPnts); 4329 } 4330 4331 public static void recycle(LineSrfIntersect instance) { 4332 FACTORY.recycle(instance); 4333 } 4334 4335 private LineSrfIntersect() { } 4336 4337 private static final ObjectFactory<LineSrfIntersect> FACTORY = new ObjectFactory<LineSrfIntersect>() { 4338 @Override 4339 protected LineSrfIntersect create() { 4340 return new LineSrfIntersect(); 4341 } 4342 4343 @Override 4344 protected void cleanup(LineSrfIntersect obj) { 4345 obj._thisSrf = null; 4346 obj._L0 = null; 4347 obj._L0v = null; 4348 obj._Ldir = null; 4349 obj._output = null; 4350 obj._approxPnts = null; 4351 } 4352 }; 4353 4354 } 4355 4356 /** 4357 * An Evaluatable1D function that calculates the signed distance between a point along 4358 * an infinite line and a point on a surface. This is used by a 1D root finder to 4359 * drive that distance to zero. 4360 */ 4361 private static class LineSrfIntersectNDFun extends AbstractEvaluatable1D { 4362 4363 private Surface _thisSrf; 4364 private Point _L0; 4365 private Vector _Lu; 4366 private double _nearS, _nearT; 4367 private boolean firstPass; 4368 4369 public static LineSrfIntersectNDFun newInstance(Surface thisSrf, GeomPoint L0, GeomVector Ldir, SubrangePoint aPnt) { 4370 4371 if (thisSrf instanceof GeomTransform) 4372 thisSrf = (Surface)thisSrf.copyToReal(); 4373 4374 LineSrfIntersectNDFun o = FACTORY.object(); 4375 o._thisSrf = thisSrf; 4376 o._L0 = L0.immutable(); 4377 o._Lu = Ldir.immutable(); 4378 o._nearS = aPnt.getParPosition().getValue(0); 4379 o._nearT = aPnt.getParPosition().getValue(1); 4380 o.firstPass = true; 4381 4382 return o; 4383 } 4384 4385 /** 4386 * Function that calculates the signed distance between a point on an infinite 4387 * line [ Q(u) = L0 + Lu*t ] and the closest point on a parametric surface [ 4388 * P(s,t) ]. 4389 * 4390 * @param u Parametric distance along the infinite line from L0. 4391 * @return The signed distance between the point on the line at u and the closest 4392 * point on the surface. 4393 */ 4394 @Override 4395 public double function(double u) throws RootException { 4396 // Don't compute anything if the 1st pass flag is set. 4397 if (firstPass) { 4398 firstPass = false; 4399 return 0; 4400 } 4401 4402 StackContext.enter(); 4403 try { 4404 // Compute the point on the line at u. 4405 Point dQu = Point.valueOf(_Lu.times(u)).plus(_Lu.getOrigin()); 4406 Point Qu = _L0.plus(dQu); 4407 4408 // Find the closest point on the surface near "aPnt". 4409 SubrangePoint Pst = _thisSrf.getClosest(Qu, _nearS, _nearT, GTOL); 4410 4411 // Find the signed distance between Pst and Qu. 4412 Vector<Length> dP = Qu.toGeomVector().minus(Pst.toGeomVector()); 4413 double dPsign = dP.dot(_Lu).getValue(); // Projection of dP onto Lu. 4414 dPsign = dPsign >= 0 ? 1 : -1; 4415 4416 double d = dP.magValue(); // Unsigned distance. 4417 d = d * dPsign; // Signed distance. 4418 4419 return d; 4420 4421 } finally { 4422 StackContext.exit(); 4423 } 4424 } 4425 4426 public static void recycle(LineSrfIntersectNDFun instance) { 4427 FACTORY.recycle(instance); 4428 } 4429 4430 private LineSrfIntersectNDFun() { } 4431 4432 @SuppressWarnings("unchecked") 4433 private static final ObjectFactory<LineSrfIntersectNDFun> FACTORY = new ObjectFactory<LineSrfIntersectNDFun>() { 4434 @Override 4435 protected LineSrfIntersectNDFun create() { 4436 return new LineSrfIntersectNDFun(); 4437 } 4438 4439 @Override 4440 protected void cleanup(LineSrfIntersectNDFun obj) { 4441 obj._thisSrf = null; 4442 obj._L0 = null; 4443 obj._Lu = null; 4444 } 4445 }; 4446 } 4447 4448 /** 4449 * An Evaluatable1D function that calculates the signed X,Y distance between a point 4450 * along the Z-axis and a point on a surface. This is used by an ND root finder to 4451 * drive that distance to zero. 4452 */ 4453 private static class LineSrfIntersect3DFun implements VectorFunction { 4454 4455 private Surface _thisSrf; 4456 4457 public static LineSrfIntersect3DFun newInstance(Surface thisSrf) { 4458 4459 if (thisSrf instanceof GeomTransform) 4460 thisSrf = (Surface)thisSrf.copyToReal(); 4461 4462 LineSrfIntersect3DFun o = FACTORY.object(); 4463 o._thisSrf = thisSrf; 4464 4465 return o; 4466 } 4467 4468 /** 4469 * User supplied method that calculates the function d[X,Y] = fn(x[s,t]). This 4470 * calculates the distance along the X & Y axes from a point on the surface to the 4471 * Z-axis. 4472 * 4473 * @param n The number of variables in the x & y arrays (should always be 2 for 4474 * this problem). 4475 * @param x Independent parameters to the function (s, t), passed in as input. 4476 * @param d An existing array that is filled in with the outputs of the function 4477 */ 4478 @Override 4479 public void function(int n, double x[], double[] d) throws RootException { 4480 double s = x[0]; 4481 double t = x[1]; 4482 4483 // Keep s & t in-bounds. 4484 double ff = 1; // Fudge factor for being out of bounds. 4485 if (s > 1) { 4486 ff = s; 4487 s = 1; 4488 } else if (s < 0) { 4489 ff = 1 - s; 4490 s = 0; 4491 } 4492 if (t > 1) { 4493 ff *= t; 4494 t = 1; 4495 } else if (t < 0) { 4496 ff *= 1 - t; 4497 t = 0; 4498 } 4499 ff *= ff; 4500 4501 StackContext.enter(); 4502 try { 4503 4504 // Get the point on the surface. 4505 Point Pst = _thisSrf.getRealPoint(s, t); 4506 4507 // Extract distance to the Z-axis from the surface point. 4508 d[0] = Pst.getValue(0) * ff; 4509 d[1] = Pst.getValue(1) * ff; 4510 4511 } finally { 4512 StackContext.exit(); 4513 } 4514 } 4515 4516 /** 4517 * User supplied method that calculates the Jacobian of the function. 4518 * 4519 * @param n The number of rows and columns in the Jacobian (always 2 for this 4520 * problem). 4521 * @param jac The Jacobian array to be filled in by this function. 4522 * @return True if the Jacobian was computed by this method, false if it was not. 4523 */ 4524 @Override 4525 public boolean jacobian(int n, double[] x, double[][] jac) { 4526 return false; 4527 } 4528 4529 public static void recycle(LineSrfIntersect3DFun instance) { 4530 FACTORY.recycle(instance); 4531 } 4532 4533 private LineSrfIntersect3DFun() { } 4534 4535 @SuppressWarnings("unchecked") 4536 private static final ObjectFactory<LineSrfIntersect3DFun> FACTORY = new ObjectFactory<LineSrfIntersect3DFun>() { 4537 @Override 4538 protected LineSrfIntersect3DFun create() { 4539 return new LineSrfIntersect3DFun(); 4540 } 4541 4542 @Override 4543 protected void cleanup(LineSrfIntersect3DFun obj) { 4544 obj._thisSrf = null; 4545 } 4546 }; 4547 } 4548 4549 /** 4550 * Return a new list of subrange points that contains those points in the input list 4551 * that are NOT within the specified tolerance of the specified curve. 4552 * 4553 * @param pntLst A list of points subranged onto _thisSrf to be compared to the curve. 4554 * @param crv A subrange curve on _thisSrf that is to be compared with the points. 4555 * @param tol The tolerance for determining if the input points are close to the 4556 * input curve. 4557 * @return A new list of subrange points that are not near the specified curve. 4558 */ 4559 private static PointString<SubrangePoint> removePntsNearCurve(PointString<SubrangePoint> pntLst, SubrangeCurve crv, Parameter<Length> tol) { 4560 4561 PointString<SubrangePoint> output = PointString.newInstance(); 4562 int size = pntLst.size(); 4563 for (int i = 0; i < size; ++i) { 4564 SubrangePoint p1 = pntLst.get(i); 4565 SubrangePoint p2 = crv.getClosest(p1.immutable(), GTOL); 4566 Parameter<Length> d = p1.distance(p2); 4567 if (!d.isLessThan(tol)) { 4568 output.add(p1); 4569 } 4570 SubrangePoint.recycle(p2); 4571 } 4572 4573 return output; 4574 } 4575 4576 /** 4577 * Find the index of the point in the given list that is parametrically closest to the 4578 * specified point. All points must be subranges to the same surface. 4579 * 4580 * @param pnt The point being compared for nearness. 4581 * @param pointList The list of points subranged to the same surface as "pnt" for 4582 * which the nearest point is to be found. 4583 * @return Index in the list of the point that is parametrically nearest to "pnt". 4584 */ 4585 private static int findParNearestInList(SubrangePoint pnt, List<SubrangePoint> pointList) { 4586 int idx = -1; 4587 double minD2 = Double.MAX_VALUE; 4588 4589 GeomPoint pst = pnt.getParPosition(); 4590 int numPnts = pointList.size(); 4591 for (int i = 0; i < numPnts; ++i) { 4592 SubrangePoint p2 = pointList.get(i); 4593 double d2 = pst.distanceSqValue(p2.getParPosition()); 4594 if (d2 < minD2) { 4595 minD2 = d2; 4596 idx = i; 4597 } 4598 } 4599 4600 return idx; 4601 } 4602 4603 /** 4604 * Create a subrange curve from a string of points subranged onto a surface. The 4605 * parametric curve will be a cubic curve passing through the parametric positions of 4606 * all the input subrange points. 4607 * 4608 * @param tracedPnts A list of points all subranged onto the same surface. 4609 * @return A subrange curve on the surface that passes through the list of subrange 4610 * points. 4611 * @throws IllegalArgumentException 4612 */ 4613 private static SubrangeCurve makeSubrangeCurve(PointString<SubrangePoint> tracedPnts) throws IllegalArgumentException { 4614 if (tracedPnts.size() < 2) 4615 return null; 4616 4617 // Get the surface that the points are subranged onto. 4618 Surface srf = (Surface)tracedPnts.get(0).getChild(); 4619 4620 int degree = 3; 4621 BasicNurbsCurve pcrv; 4622 StackContext.enter(); 4623 try { 4624 // Extract the parametric positions of the points. 4625 PointString<Point> parStr = PointString.newInstance(); 4626 for (SubrangePoint sp : tracedPnts) { 4627 Point par = sp.getParPosition().immutable(); 4628 parStr.add(par); 4629 } 4630 4631 // Thin out any positions that are to close together. 4632 Parameter ptol2 = Parameter.valueOf(SQRT_EPS, SI.METER).pow(2); 4633 int size = parStr.size(); 4634 if (size > 2 && parStr.get(0).distanceSq(parStr.get(1)).isLessThan(ptol2)) { 4635 parStr.remove(1); 4636 --size; 4637 } 4638 Point lastPar = parStr.get(-1); 4639 for (int i=size-2; i > 0; --i) { 4640 Point par = parStr.get(i); 4641 if (par.distanceSq(lastPar).isLessThan(ptol2)) 4642 parStr.remove(i); 4643 lastPar = par; 4644 } 4645 if (parStr.size() <= degree) 4646 degree = parStr.size() - 1; 4647 4648 // Create a parametric curve for the traced points. 4649 try { 4650 pcrv = CurveFactory.fitPoints(degree, parStr); 4651 } catch (SingularMatrixException e) { 4652 // Fall back on just using a straight line segment. 4653 pcrv = CurveFactory.createLine(parStr.get(0), parStr.get(-1)); 4654 } 4655 4656 // Make sure this parameter curve is in-bounds. 4657 Point pcrvMin = pcrv.getBoundsMin(); 4658 Point pcrvMax = pcrv.getBoundsMax(); 4659 while (pcrvMin.getValue(0) < 0 || pcrvMin.getValue(1) < 0 4660 || pcrvMax.getValue(0) > 1 || pcrvMax.getValue(1) > 1) { 4661 // Lower the degree and try again. 4662 --degree; 4663 if (degree == 0) { 4664 // Just connect the end points with a line. 4665 pcrv = CurveFactory.createLine(parStr.get(0), parStr.get(-1)); 4666 break; 4667 } 4668 pcrv = CurveFactory.fitPoints(degree, parStr); 4669 4670 pcrvMin = pcrv.getBoundsMin(); 4671 pcrvMax = pcrv.getBoundsMax(); 4672 } 4673 4674 // Remove redundant knots to within a fine tolerance. 4675 pcrv = CurveUtils.thinKnotsToTolerance(pcrv, GTOLP); 4676 4677 pcrv = StackContext.outerCopy(pcrv); 4678 4679 } finally { 4680 StackContext.exit(); 4681 } 4682 4683 // Create the subrange curve. 4684 SubrangeCurve scrv = SubrangeCurve.newInstance(srf, pcrv); 4685 4686 return scrv; 4687 } 4688 4689 /** 4690 * A class used to intersect an infinite plane and the specified surface. 4691 */ 4692 private static class PlaneSrfIntersect { 4693 4694 private static final double GTOL100 = GTOL * 100; // A coarser version of GTOL. 4695 private static final double EPS_FT = 0.866; // Planar flatness angle tolerance : cos(30 deg) 4696 private static final double DTHETA = 0.349; // Tracing step size angle: DTheta = 20 deg = 0.349 rad 4697 private static final int MAXDEPTH = 15; 4698 4699 private AbstractSurface _thisSrf; 4700 private GeomPlane _plane; 4701 private Parameter<Length> _tol; 4702 private Parameter<Length> _patchSizeTol; 4703 private Parameter<Length> _eps1; 4704 private Parameter<Length> _eps2; 4705 private Parameter<Length> _epsCRT; // Curve refinement tolerance. 4706 private GeomList<SubrangeCurve> _output; 4707 private PointString<SubrangePoint> _appInteriorPnts; 4708 4709 public static PlaneSrfIntersect newInstance(AbstractSurface thisSrf, 4710 GeomPlane plane, Parameter<Length> tol, GeomList<SubrangeCurve> output) { 4711 4712 if (thisSrf instanceof GeomTransform) 4713 thisSrf = (AbstractSurface)thisSrf.copyToReal(); 4714 if (plane instanceof GeomTransform) 4715 plane = plane.copyToReal(); 4716 4717 @SuppressWarnings("null") 4718 Parameter<Length> thisSpan = thisSrf.getBoundsMax().distance(thisSrf.getBoundsMin()); 4719 4720 // Create an instance of this class and fill in the class variables. 4721 PlaneSrfIntersect o = FACTORY.object(); 4722 o._thisSrf = thisSrf; 4723 o._plane = plane; 4724 o._tol = tol; 4725 o._patchSizeTol = tol.times(10); 4726 o._eps1 = thisSpan.divide(5000); 4727 o._eps2 = o._eps1.times(2); 4728 o._epsCRT = thisSpan.divide(50); 4729 o._output = output; 4730 o._appInteriorPnts = PointString.newInstance(); 4731 4732 /* 4733 System.out.println("tol = " + tol); 4734 System.out.println("eps1 = " + o._eps1); 4735 System.out.println("eps2 = " + o._eps2); 4736 System.out.println("epsCRT = " + o._epsCRT); 4737 */ 4738 return o; 4739 } 4740 4741 /** 4742 * Method that actually carries out the intersection adding the intersection 4743 * points to the list provided in the constructor. 4744 * 4745 * @return A list containing zero or more subrange curves, on this surface, made 4746 * by the intersection of this surface with the specified plane. If no 4747 * intersections are found an empty GeomList is returned. 4748 */ 4749 public GeomList<SubrangeCurve> intersect() throws RootException { 4750 4751 // This algorithm is inspired by the Barnhill-Kersey surface-surface intersection 4752 // algorithm, but with some fairly major differences. 4753 // Barnhill, R.E., and Kersey, S.N., "Marching Method for Parametric Surface/Surface Intersection," CAGD, 7(1-4), June 1990, 257-280. 4754 // as discussed in: Max K. Agoston, Computer Graphics and Geometric Modelling: Implementation & Algorithms 4755 4756 // The high level outline of the algorithm as implemented is as follows: 4757 // 1 - Find edge curve intersections to get initial starting points for intersection tracing. 4758 // 2 - Trace intersections that start at edge points. 4759 // 3 - Use a divide & conquer approach to finding a spattering of interior intersection points. 4760 // 4 - Remove any interior points that are associated with the intersection curves already found. 4761 // 5 - Trace intersections starting at interior points removing any interior points that 4762 // are associated with these intersections. 4763 4764 // Get any edge intersection points. 4765 //System.out.println("Extracting edge points..."); 4766 PointString<SubrangePoint> edgePoints = getEdgePoints(_thisSrf, _plane); 4767 4768 // Begin tracing from the edge points. 4769 //System.out.println("Tracing from edge points..."); 4770 GeomList<SubrangeCurve> crvSegs = traceEdgePoints(edgePoints, _plane); 4771 _output.addAll(crvSegs); 4772 4773 // Use a divide and conquer approach to find approximate interior intersection points. 4774 //System.out.println("Finding interior points..."); 4775 divideAndConquerPlane(1, _thisSrf, 0, 0, 1, 1); 4776 4777 // Refine each approximate intersection and save as a potential start point 4778 // for an internal loop. 4779 PointString<SubrangePoint> potStartPoints = PointString.newInstance(); 4780 int size = _appInteriorPnts.size(); 4781 for (int i = 0; i < size; ++i) { 4782 SubrangePoint aPnt = _appInteriorPnts.get(i); 4783 SubrangePoint rpnt = relaxPoint(aPnt, _plane, _tol); 4784 if (nonNull(rpnt)) 4785 potStartPoints.add(rpnt); 4786 } 4787 4788 // Remove any points that are near the existing intersection curves. 4789 //System.out.println("Removing redundant interior points..."); 4790 removeClosePoints(potStartPoints, _tol); 4791 Parameter<Length> tol2 = _tol.times(100); 4792 size = _output.size(); 4793 for (int i = 0; i < size; ++i) { 4794 SubrangeCurve crv = _output.get(i); 4795 potStartPoints = removePntsNearCurve(potStartPoints, crv, tol2); 4796 } 4797 4798 // Any remaining interior points are assumed to be on interior loops. 4799 //System.out.println("Tracing from interior points..."); 4800 crvSegs = traceInteriorPoints(potStartPoints, _plane, _patchSizeTol); 4801 _output.addAll(crvSegs); 4802 4803 return _output; 4804 } 4805 4806 /** 4807 * Trace intersection curves that start/end at the edges of the surface. 4808 * 4809 * @param edgePoints The surface edge points to use as starting points for tracing 4810 * intersection curves. 4811 * @param plane The plane being used for intersecting with the surface. 4812 * @return A list of intersection curves traced from/to each edge point. 4813 * @throws RootException If there is a convergence problem with a root finder. 4814 */ 4815 private GeomList<SubrangeCurve> traceEdgePoints(PointString<SubrangePoint> edgePoints, 4816 GeomPlane plane) throws RootException { 4817 4818 GeomList<SubrangeCurve> crvSegs = GeomList.newInstance(); 4819 int numPnts = edgePoints.size(); 4820 4821 while (numPnts > 0) { 4822 SubrangePoint startPnt = edgePoints.get(0); 4823 4824 // Trace the intersection starting at "startPnt". 4825 SubrangeCurve seg = traceSegment(startPnt, plane); 4826 if (nonNull(seg)) { 4827 // Don't keep degenerate curves (points). 4828 if (!seg.isDegenerate(_tol)) 4829 crvSegs.add(seg); 4830 4831 // Remove the end point from the list. 4832 Surface srf = (Surface)startPnt.getChild(); 4833 SubrangePoint endPoint = srf.getPoint(seg.getParPosition().getRealPoint(1)); 4834 int idx = findParNearestInList(endPoint, edgePoints); 4835 if (idx > 0) { 4836 edgePoints.remove(idx); 4837 --numPnts; 4838 } else { 4839 Logger.getLogger(AbstractSurface.class.getName()).log(Level.WARNING, 4840 "End point not found."); 4841 } 4842 } 4843 4844 // Remove the start point from the list. 4845 edgePoints.remove(0); 4846 --numPnts; 4847 } 4848 4849 return crvSegs; 4850 } 4851 4852 /** 4853 * Trace intersection curves that start/end on interior points. 4854 * 4855 * @param potStartPoints A list of potential start points for intersection 4856 * tracing. 4857 * @param plane The plane used to intersect the surface. 4858 * @param tol The tolerance to use when deciding if a point belongs to 4859 * an existing intersection curve. 4860 * @return A list of intersection curves traced from/to the interior points 4861 * provided. 4862 * @throws RootException If there is a convergence problem with a root finder. 4863 */ 4864 private GeomList<SubrangeCurve> traceInteriorPoints(PointString<SubrangePoint> potStartPoints, GeomPlane plane, Parameter<Length> tol) throws RootException { 4865 4866 GeomList<SubrangeCurve> crvSegs = GeomList.newInstance(); 4867 4868 while (!potStartPoints.isEmpty()) { 4869 SubrangePoint startPnt = potStartPoints.get(0); 4870 4871 // Trace the intersection starting at "startPnt". 4872 SubrangeCurve seg = traceSegment(startPnt, plane); 4873 if (nonNull(seg)) { 4874 // Don't keep degenerate curves (points). 4875 if (!seg.isDegenerate(_tol)) 4876 crvSegs.add(seg); 4877 4878 // Remove any points that are close to this curve segment. 4879 potStartPoints.remove(0); // The start point obviously. 4880 potStartPoints = removePntsNearCurve(potStartPoints, seg, tol); 4881 4882 } else { 4883 potStartPoints.remove(0); 4884 } 4885 } 4886 4887 return crvSegs; 4888 } 4889 4890 /** 4891 * Method that finds all the exact intersections of the edges of the surface (if 4892 * any). 4893 */ 4894 private PointString<SubrangePoint> getEdgePoints(Surface srf, GeomPlane plane) { 4895 4896 PointString<SubrangePoint> output = PointString.newInstance(); 4897 4898 PointString<SubrangePoint> edge = srf.getS0Curve().intersect(plane, GTOL); 4899 int size = edge.size(); 4900 for (int i = 0; i < size; ++i) { 4901 SubrangePoint cPnt = edge.get(i); 4902 double t = cPnt.getParPosition().getValue(0); 4903 SubrangePoint sPnt = srf.getPoint(0, t); 4904 output.add(sPnt); 4905 SubrangePoint.recycle(cPnt); 4906 } 4907 PointString.recycle(edge); 4908 edge = srf.getT0Curve().intersect(plane, GTOL); 4909 size = edge.size(); 4910 for (int i = 0; i < size; ++i) { 4911 SubrangePoint cPnt = edge.get(i); 4912 double s = cPnt.getParPosition().getValue(0); 4913 SubrangePoint sPnt = srf.getPoint(s, 0); 4914 output.add(sPnt); 4915 SubrangePoint.recycle(cPnt); 4916 } 4917 PointString.recycle(edge); 4918 edge = srf.getS1Curve().intersect(plane, GTOL); 4919 size = edge.size(); 4920 for (int i = 0; i < size; ++i) { 4921 SubrangePoint cPnt = edge.get(i); 4922 double t = cPnt.getParPosition().getValue(0); 4923 SubrangePoint sPnt = srf.getPoint(1, t); 4924 output.add(sPnt); 4925 SubrangePoint.recycle(cPnt); 4926 } 4927 PointString.recycle(edge); 4928 edge = srf.getT1Curve().intersect(plane, GTOL); 4929 size = edge.size(); 4930 for (int i = 0; i < size; ++i) { 4931 SubrangePoint cPnt = edge.get(i); 4932 double s = cPnt.getParPosition().getValue(0); 4933 SubrangePoint sPnt = srf.getPoint(s, 1); 4934 output.add(sPnt); 4935 SubrangePoint.recycle(cPnt); 4936 } 4937 PointString.recycle(edge); 4938 4939 // Remove any redundant points 4940 removeClosePoints(output, _tol); 4941 4942 return output; 4943 } 4944 4945 /** 4946 * Uses a recursive "Divide and Conquer" approach to intersecting a surface patch 4947 * with a plane. On each call, the following happens: 4948 * <pre> 4949 * 1) The patch is checked to see if it is approx. a flat plane. 4950 * 1b) If it is, then a plane-plane intersection is performed, the 4951 * approximate intersection points added to the _approxPnts list and the method exited. 4952 * 2) The patch is divided into four quadrant patches. 4953 * 2b) Each patch is tested to see if it's bounding box is intersected by the plane. 4954 * If it is, then this method is called with that patch recursively. 4955 * Otherwise, the method is exited. 4956 * </pre> 4957 * A number of class variables are used to pass information to this recursion: 4958 * <pre> 4959 * _thisSrf The full surfaces intersections are being found for. 4960 * _plane The plane being intersected with this surface. 4961 * _tol The tolerance to use in determining if the geometry is in tolerance. 4962 * _appInteriorPnts A list used to collect the approximate subrange intersection points. 4963 * </pre> 4964 * 4965 * @param patch The surface patch being checked for intersection with the plane. 4966 * @param s0 The "s" parametric position of the start of the patch on the 4967 * master "this" surface. 4968 * @param t0 The "t" parametric position of the start of the patch on the 4969 * master "this" surface. 4970 * @param s1 The parametric "s" position of the end of the patch on the master 4971 * "this" surface. 4972 * @param t1 The parametric "t" position of the end of the patch on the master 4973 * "this" surface. 4974 */ 4975 private void divideAndConquerPlane(int depth, Surface patch, double s0, double t0, double s1, double t1) { 4976 4977 // Check to see if this patch is nearly planar. 4978 // Using a fast method for "isPlanar()" that unfortunately ignores "tol". 4979 if (depth > MAXDEPTH || patch.isPlanar(EPS_FT, EPS_FT) || isSmall(patch, _patchSizeTol)) { 4980 4981 // Form a single quadrilateral panel from the corners of the patch. 4982 Point A = patch.getRealPoint(0, 0); 4983 Point B = patch.getRealPoint(1, 0); 4984 Point C = patch.getRealPoint(1, 1); 4985 Point D = patch.getRealPoint(0, 1); 4986 4987 // Determine if any edge of the quadrilateral polygon intersects with the plane. 4988 GeomPoint P0 = _plane.getRefPoint(); 4989 GeomVector<Dimensionless> nhat = _plane.getNormal(); 4990 boolean intersects = GeomUtil.planeQuadIntersect(P0, nhat, A, B, C, D); 4991 4992 // Store off the approximate intersection point (panel center). 4993 if (intersects) { 4994 double s = segmentPos2Parent(0.5, s0, s1); 4995 double t = segmentPos2Parent(0.5, t0, t1); 4996 SubrangePoint centerPnt = _thisSrf.getPoint(s, t); 4997 _appInteriorPnts.add(centerPnt); 4998 } 4999 5000 } else { 5001 // Subdivide the patch into 4 quadrants. 5002 GeomList<Surface> leftRight = patch.splitAtS(0.5); 5003 GeomList<Surface> botTop = leftRight.get(0).splitAtT(0.5); 5004 Surface p00 = botTop.get(0); // Lower Left 5005 Surface p01 = botTop.get(1); // Top Left 5006 GeomList.recycle(botTop); 5007 botTop = leftRight.get(1).splitAtT(0.5); 5008 Surface p10 = botTop.get(0); // Lower Right 5009 Surface p11 = botTop.get(1); // Top Right 5010 GeomList.recycle(botTop); 5011 GeomList.recycle(leftRight); 5012 5013 // Check for possible intersections on the lower-left patch. 5014 if (GeomUtil.planeBoundsIntersect(_plane, p00)) { 5015 // May be an intersection. 5016 double sl = s0; 5017 double sh = 0.5 * (s0 + s1); 5018 double tl = t0; 5019 double th = 0.5 * (t0 + t1); 5020 5021 // Recurse down to a finer level of detail. 5022 divideAndConquerPlane(depth + 1, p00, sl, tl, sh, th); 5023 } 5024 if (p00 instanceof BasicNurbsSurface) 5025 BasicNurbsSurface.recycle((BasicNurbsSurface)p00); 5026 5027 // Check for possible intersections on the lower-right patch. 5028 if (GeomUtil.planeBoundsIntersect(_plane, p10)) { 5029 // May be an intersection. 5030 double sl = 0.5 * (s0 + s1); 5031 double sh = s1; 5032 double tl = t0; 5033 double th = 0.5 * (t0 + t1); 5034 5035 // Recurse down to a finer level of detail. 5036 divideAndConquerPlane(depth + 1, p10, sl, tl, sh, th); 5037 } 5038 if (p10 instanceof BasicNurbsSurface) 5039 BasicNurbsSurface.recycle((BasicNurbsSurface)p10); 5040 5041 // Check for possible intersections on the top-left patch. 5042 if (GeomUtil.planeBoundsIntersect(_plane, p01)) { 5043 // May be an intersection. 5044 double sl = s0; 5045 double sh = 0.5 * (s0 + s1); 5046 double tl = 0.5 * (t0 + t1); 5047 double th = t1; 5048 5049 // Recurse down to a finer level of detail. 5050 divideAndConquerPlane(depth + 1, p01, sl, tl, sh, th); 5051 } 5052 if (p01 instanceof BasicNurbsSurface) 5053 BasicNurbsSurface.recycle((BasicNurbsSurface)p01); 5054 5055 // Check for possible intersections on the top-right patch. 5056 if (GeomUtil.planeBoundsIntersect(_plane, p11)) { 5057 // May be an intersection. 5058 double sl = 0.5 * (s0 + s1); 5059 double sh = s1; 5060 double tl = 0.5 * (t0 + t1); 5061 double th = t1; 5062 5063 // Recurse down to a finer level of detail. 5064 divideAndConquerPlane(depth + 1, p11, sl, tl, sh, th); 5065 } 5066 if (p11 instanceof BasicNurbsSurface) 5067 BasicNurbsSurface.recycle((BasicNurbsSurface)p11); 5068 5069 } 5070 5071 } 5072 5073 private static final int MAXIT = 500; 5074 5075 /** 5076 * Relax the approximate intersection point toward the exact intersection until it 5077 * is within "tol" of the exact intersection. 5078 * 5079 * @param approx The approximate surface intersection point being refined or 5080 * relaxed. 5081 * @param plane The intersecting plane. 5082 * @param tol THe tolerance required of the intersection points. 5083 * @return The input point relaxed onto the intersection between the plane and the 5084 * surface to within the tolerance "tol". If the plane is locally tangent 5085 * to the surface near the intersection, this may return 5086 * <code>null</code>. 5087 */ 5088 private SubrangePoint relaxPoint(SubrangePoint approx, GeomPlane plane, 5089 Parameter<Length> tol) throws RootException { 5090 5091 AbstractSurface srf = (AbstractSurface)approx.getChild(); 5092 5093 double sOut, tOut; 5094 StackContext.enter(); 5095 try { 5096 int numDims = srf.getPhyDimension(); 5097 SubrangePoint p = approx; 5098 Parameter<Length> d = GeomUtil.pointPlaneDistance(p, plane.getRefPoint(), plane.getNormal()); 5099 if (!d.abs().isLessThan(tol)) { 5100 5101 Unit<Length> unit = srf.getUnit(); 5102 MutableVector Lu = MutableVector.newInstance(numDims, Dimensionless.UNIT); 5103 MutablePoint L0 = MutablePoint.newInstance(numDims, unit); 5104 5105 int iteration = 0; 5106 do { 5107 // Get the local tangent plane on the surface. 5108 GeomPoint st = p.getParPosition(); 5109 GeomPlane Ts = srf.getTangentPlane(st); 5110 5111 // Intersect the input plane with the surface tangent plane. 5112 IntersectType type = GeomUtil.planePlaneIntersect(plane, Ts, L0, Lu); 5113 if (type != IntersectType.INTERSECT) { 5114 return null; 5115 } 5116 5117 // Find the closest point on the plane intersection line to the previous point. 5118 Point pip1 = GeomUtil.pointLineClosest(p, L0, Lu); 5119 5120 // Re-attach this point to the surface. 5121 p = srf.getClosestFarthestInterior(pip1, st.getValue(0), st.getValue(1), GTOL100, true); 5122 5123 // Update the distance between this new point and the intersecting plane. 5124 d = GeomUtil.pointPlaneDistance(p, plane.getRefPoint(), plane.getNormal()); 5125 ++iteration; 5126 5127 } while (!d.abs().isLessThan(tol) && iteration < MAXIT); 5128 5129 if (iteration == MAXIT) 5130 Logger.getLogger(AbstractSurface.class.getName()).log(Level.WARNING, 5131 "Convergence problem in PlaneSrfIntersect.relaxPoint()."); 5132 } 5133 5134 // Copy parametric position of final point to outer context. 5135 sOut = p.getParPosition().getValue(0); 5136 tOut = p.getParPosition().getValue(1); 5137 5138 } finally { 5139 StackContext.exit(); 5140 } 5141 5142 return srf.getPoint(sOut, tOut); 5143 } 5144 5145 /** 5146 * Method that traces an intersection segment from the given start point until an 5147 * edge of the surface is encountered or the start point is re-encountered. 5148 * 5149 * @param startPnt The point on the surface to start tracing the intersection 5150 * from. 5151 * @param plane The plane that is being intersected with the surface. 5152 * @return A subrange curve representing this intersection segment. 5153 */ 5154 private SubrangeCurve traceSegment(SubrangePoint startPnt, GeomPlane plane) throws RootException { 5155 AbstractSurface srf = (AbstractSurface)startPnt.getChild(); 5156 5157 BasicNurbsCurve pCrv; 5158 StackContext.enter(); 5159 try { 5160 PointString<SubrangePoint> tracedPnts = PointString.newInstance(); 5161 5162 Unit<Length> unit = srf.getUnit(); 5163 GeomVector<Dimensionless> n2 = plane.getNormal(); 5164 SubrangePoint p1 = startPnt; 5165 tracedPnts.add(p1); 5166 GeomPoint st = p1.getParPosition(); // Parametric position on the surface of the point p1. 5167 double s0 = st.getValue(0); 5168 double t0 = st.getValue(1); 5169 5170 // Is the start point an edge point? 5171 boolean startEdge = false; 5172 if (parNearEnds(s0, GTOL100) || parNearEnds(t0, GTOL100)) { 5173 startEdge = true; 5174 } 5175 5176 boolean firstPoint = true; 5177 int collapsedEdge = -1; 5178 while (tracedPnts.size() < 1000) { 5179 // Determine the step direction. 5180 Vector<Dimensionless> pu = srf.getSDerivative(st.getValue(0), st.getValue(1), 1, false).toUnitVector(); 5181 if (!pu.isValid()) { 5182 // Step away from a collapsed edge. 5183 double t = st.getValue(1); 5184 if (t < 0.5) { 5185 st = Point.valueOf(SI.METER, 0.5, GTOL100); 5186 collapsedEdge = 2; 5187 } else { 5188 st = Point.valueOf(SI.METER, 0.5, 1. - GTOL100); 5189 collapsedEdge = 3; 5190 } 5191 if (startEdge) { 5192 p1 = srf.getPoint(st); 5193 s0 = st.getValue(0); 5194 t0 = st.getValue(1); 5195 } 5196 pu = srf.getSDerivative(st.getValue(0), st.getValue(1), 1, false).toUnitVector(); 5197 } 5198 Vector<Dimensionless> pv = srf.getTDerivative(st.getValue(0), st.getValue(1), 1, false).toUnitVector(); 5199 if (!pv.isValid()) { 5200 // Step away from a collapsed edge. 5201 double s = st.getValue(0); 5202 if (s < 0.5) { 5203 st = Point.valueOf(SI.METER, GTOL100, 0.5); 5204 collapsedEdge = 0; 5205 } else { 5206 st = Point.valueOf(SI.METER, 1. - GTOL100, 0.5); 5207 collapsedEdge = 1; 5208 } 5209 pu = pv = srf.getTDerivative(st.getValue(0), st.getValue(1), 1, false).toUnitVector(); 5210 if (startEdge) { 5211 p1 = srf.getPoint(st); 5212 s0 = st.getValue(0); 5213 t0 = st.getValue(1); 5214 } 5215 } 5216 Parameter pudotn2 = pu.dot(n2); 5217 Parameter pvdotn2 = pv.dot(n2); 5218 if (pudotn2.isApproxZero() && pvdotn2.isApproxZero()) { 5219 Logger.getLogger(AbstractSurface.class.getName()).log(Level.WARNING, 5220 "Surface is tangent to plane in traceSegment()."); 5221 break; // Surface is nearly tangent to plane. TODO: Handle this better. 5222 } 5223 // Ti = pv*(pu dot n2) - pu*(pv dot n2) 5224 Vector<Dimensionless> Ti = pv.times(pudotn2).minus((Vector)pu.times(pvdotn2)); 5225 Ti = Ti.toUnitVector(); 5226 5227 // Determine the step size "L". 5228 Point p = p1.plus(Point.valueOf(Ti.times(_eps1))); 5229 SubrangePoint x = srf.getClosestFarthestInterior(p, st.getValue(0), st.getValue(1), GTOL, true); 5230 5231 // If "Ti" initially steps us off the edge of the surface, reverse it. 5232 if (startEdge) { 5233 startEdge = false; 5234 double sx = x.getParPosition().getValue(0); 5235 double tx = x.getParPosition().getValue(1); 5236 if (sx < GTOL || sx >= (1 - GTOL) || tx < GTOL || tx > (1 - GTOL) 5237 || abs(sx - s0) > 0.3 || abs(tx - t0) > 0.3) { 5238 // Ti stepped off the edge. 5239 Ti = Ti.opposite(); 5240 p = p1.plus(Point.valueOf(Ti.times(_eps1))); 5241 x = srf.getClosestFarthestInterior(p, st.getValue(0), st.getValue(1), GTOL100, true); 5242 5243 // Reverse the intersecting plane normal so this doesn't happen again. 5244 n2 = n2.opposite(); 5245 } 5246 } 5247 5248 x = relaxPoint(x, plane, _tol); 5249 if (isNull(x)) 5250 break; 5251 p = p1.plus(Point.valueOf(Ti.times(_eps2))); 5252 SubrangePoint y = srf.getClosestFarthestInterior(p, st.getValue(0), st.getValue(1), GTOL100, true); 5253 y = relaxPoint(y, plane, _tol); 5254 if (isNull(y)) 5255 break; 5256 5257 // If the 1st point was on a collapsed edge, go back and fix things up. 5258 if (collapsedEdge >= 0) { 5259 st = x.getParPosition(); 5260 switch (collapsedEdge) { 5261 case 0: 5262 st = Point.valueOf(SI.METER, 0, st.getValue(1)); 5263 break; 5264 case 1: 5265 st = Point.valueOf(SI.METER, 1., st.getValue(1)); 5266 break; 5267 case 2: 5268 st = Point.valueOf(SI.METER, st.getValue(0), 0); 5269 break; 5270 case 3: 5271 st = Point.valueOf(SI.METER, st.getValue(0), 1.); 5272 break; 5273 } 5274 p1 = srf.getPoint(st); 5275 startPnt = p1; 5276 tracedPnts.clear(); 5277 tracedPnts.add(p1); 5278 s0 = st.getValue(0); 5279 t0 = st.getValue(1); 5280 5281 collapsedEdge = -1; 5282 } 5283 5284 Vector<Length> a = x.minus(p1).toGeomVector(); 5285 Vector<Length> b = y.minus(p1).toGeomVector(); 5286 5287 double na = a.normValue(); 5288 double nb = b.normValue(); 5289 double namb = a.minus(b).normValue(); 5290 double naxb = a.cross(b).normValue(); 5291 5292 Parameter<Length> L; 5293 if (naxb > 1e-9) { 5294 // The 3 points are not co-linear. 5295 // r = |a||b||a - b|/(2*|a x b|) 5296 double r = 0.5 * na * nb * namb / naxb; 5297 L = Parameter.valueOf(r * DTHETA, unit); 5298 //System.out.print("L = " + L); 5299 if (L.isGreaterThan(_epsCRT)) 5300 L = _epsCRT; 5301 } else { 5302 L = _epsCRT; 5303 } 5304 5305 // Calculate the new intersection point location. 5306 p = p1.plus(Point.valueOf(Ti.times(L))); 5307 SubrangePoint p2 = srf.getClosestFarthestInterior(p, st.getValue(0), st.getValue(1), GTOL100, true); 5308 p2 = relaxPoint(p2, plane, _tol); 5309 if (isNull(p2)) 5310 break; 5311 //System.out.println(", d = " + p2.distance(p)); 5312 5313 // Store the new intersection point. 5314 tracedPnts.add(p2); 5315 5316 // Is this new point an edge point? 5317 st = p2.getParPosition(); 5318 double s = st.getValue(0); 5319 double t = st.getValue(1); 5320 if (!firstPoint && (parNearEnds(s, GTOL100) || parNearEnds(t, GTOL100))) 5321 break; 5322 5323 // Is this new point near the start point? 5324 if (!firstPoint && (st.distanceValue(startPnt.getParPosition()) < 0.1 5325 && p2.distance(startPnt).isLessThan(L.times(0.75)))) { 5326 tracedPnts.add(startPnt); // Repeat the 1st point. 5327 break; 5328 } 5329 5330 p1 = p2; 5331 firstPoint = false; 5332 } 5333 5334 // Create the curve segment. 5335 if (tracedPnts.size() == 1) { 5336 st = tracedPnts.get(0).getParPosition(); 5337 pCrv = StackContext.outerCopy(CurveFactory.createPoint(1, st)); 5338 5339 } else { 5340 // Create the intersection curve from the intersection points. 5341 SubrangeCurve scrv = makeSubrangeCurve(tracedPnts); 5342 pCrv = StackContext.outerCopy((BasicNurbsCurve)scrv.getParPosition()); 5343 } 5344 5345 } finally { 5346 StackContext.exit(); 5347 } 5348 5349 return SubrangeCurve.newInstance(srf, pCrv); 5350 } 5351 5352 public static void recycle(PlaneSrfIntersect instance) { 5353 FACTORY.recycle(instance); 5354 } 5355 5356 private PlaneSrfIntersect() { } 5357 5358 private static final ObjectFactory<PlaneSrfIntersect> FACTORY = new ObjectFactory<PlaneSrfIntersect>() { 5359 @Override 5360 protected PlaneSrfIntersect create() { 5361 return new PlaneSrfIntersect(); 5362 } 5363 5364 @Override 5365 protected void cleanup(PlaneSrfIntersect obj) { 5366 obj._thisSrf = null; 5367 obj._plane = null; 5368 obj._output = null; 5369 obj._appInteriorPnts = null; 5370 obj._tol = null; 5371 obj._patchSizeTol = null; 5372 obj._eps1 = null; 5373 obj._eps2 = null; 5374 obj._epsCRT = null; 5375 } 5376 }; 5377 5378 } 5379 5380 /** 5381 * A class used to intersect two arbitrary surfaces. 5382 */ 5383 private static class SrfSrfIntersect { 5384 5385 private static final double GTOL100 = GTOL * 100; // A coarser version of GTOL. 5386 private static final double EPS_FT = 0.866; // Planar flatness angle tolerance : cos(30 deg) 5387 private static final double DTHETA = 0.1745; // Tracing step size angle: DTheta = 10 deg = 0.1745 rad 5388 private static final int MAXDEPTH = 10; 5389 5390 private AbstractSurface _thisSrf; 5391 private AbstractSurface _thatSrf; 5392 private Parameter<Length> _tol; 5393 private Parameter<Length> _patchSizeTol; 5394 private Parameter<Length> _eps1; 5395 private Parameter<Length> _eps2; 5396 private Parameter<Length> _epsCRT; // Curve refinement tolerance. 5397 private GeomList<SubrangeCurve> _thisOut; 5398 private GeomList<SubrangeCurve> _thatOut; 5399 private PointString<SubrangePoint> _appInteriorPnts; 5400 5401 @SuppressWarnings("null") 5402 public static SrfSrfIntersect newInstance(AbstractSurface thisSrf, 5403 AbstractSurface thatSrf, Parameter<Length> tol, 5404 GeomList<SubrangeCurve> thisOut, GeomList<SubrangeCurve> thatOut) { 5405 5406 if (thisSrf instanceof GeomTransform) 5407 thisSrf = (AbstractSurface)thisSrf.copyToReal(); 5408 if (thatSrf instanceof GeomTransform) 5409 thatSrf = (AbstractSurface)thatSrf.copyToReal(); 5410 5411 Parameter<Length> thisSpan = thisSrf.getBoundsMax().distance(thisSrf.getBoundsMin()); 5412 Parameter<Length> thatSpan = thatSrf.getBoundsMax().distance(thatSrf.getBoundsMin()); 5413 Parameter<Length> minSpan = thisSpan.min(thatSpan); 5414 5415 // Create an instance of this class and fill in the class variables. 5416 SrfSrfIntersect o = FACTORY.object(); 5417 o._thisSrf = thisSrf; 5418 o._thatSrf = thatSrf; 5419 o._tol = tol; 5420 o._patchSizeTol = tol.times(10); 5421 o._eps1 = minSpan.divide(5000); 5422 o._eps2 = o._eps1.times(2); 5423 o._epsCRT = minSpan.divide(100); 5424 o._thisOut = thisOut; 5425 o._thatOut = thatOut; 5426 o._appInteriorPnts = PointString.newInstance(); 5427 5428 /* 5429 System.out.println("tol = " + tol); 5430 System.out.println("eps1 = " + o._eps1); 5431 System.out.println("eps2 = " + o._eps2); 5432 System.out.println("epsCRT = " + o._epsCRT); 5433 */ 5434 return o; 5435 } 5436 5437 /** 5438 * Method that actually carries out the intersection adding the intersection 5439 * points to the lists provided in the constructor. 5440 */ 5441 public void intersect() throws RootException { 5442 5443 // This algorithm is inspired by the Barnhill-Kersey surface-surface intersection 5444 // algorithm, but with some fairly major differences. 5445 // Barnhill, R.E., and Kersey, S.N., "A Marching Method for Parametric Surface/Surface Intersection," CAGD, 7(1-4), June 1990, 257-280. 5446 // as discussed in: Max K. Agoston, Computer Graphics and Geometric Modelling: Implementation & Algorithms 5447 5448 // The high level outline of the algorithm as implemented is as follows: 5449 // 1 - Find edge curve intersections to get initial starting points for intersection tracing. 5450 // 2 - Trace intersections that start at edge points. 5451 // 3 - Use a divide & conquer approach to finding a spattering of interior intersection points. 5452 // 4 - Remove any interior points that are associated with the intersection curves already found. 5453 // 5 - Trace intersections starting at interior points removing any interior points that 5454 // are associated with these intersections. 5455 5456 // Get any edge intersection points. 5457 //System.out.println("Extracting edge points..."); 5458 PointString<SubrangePoint> thisEdgePnts = getEdgePoints(_thisSrf, _thatSrf); 5459 PointString<SubrangePoint> thatEdgePnts = getEdgePoints(_thatSrf, _thisSrf); 5460 5461 // Begin tracing from the edge points. 5462 //System.out.println("Tracing from edge points..."); 5463 //System.out.println(" This surface..."); 5464 GeomList<GeomList<SubrangeCurve>> crvSegs = traceEdgePoints(thisEdgePnts, thatEdgePnts, _thatSrf); 5465 _thisOut.addAll(crvSegs.get(0)); 5466 _thatOut.addAll(crvSegs.get(1)); 5467 5468 // Remove any points that are near the existing intersection curves. 5469 //System.out.println(" Removing redundant interior points..."); 5470 Parameter<Length> tol2 = _tol.times(100); 5471 for (SubrangeCurve crv : _thatOut) { 5472 PointString<SubrangePoint> newLst = removePntsNearCurve(thatEdgePnts, crv, tol2); 5473 thatEdgePnts = newLst; 5474 } 5475 5476 //System.out.println(" That surface..."); 5477 crvSegs = traceEdgePoints(thatEdgePnts, thisEdgePnts, _thisSrf); 5478 _thatOut.addAll(crvSegs.get(0)); 5479 _thisOut.addAll(crvSegs.get(1)); 5480 5481 // Use a divide and conquer approach to find approximate interior intersection points. 5482 //System.out.println("Finding interior points..."); 5483 divideAndConquerSrf(1, _thisSrf, 0, 0, 1, 1, _thatSrf, 0, 0, 1, 1); 5484 //System.out.println(" interior points found = " + _appInteriorPnts.size()); 5485 5486 // Refine each approximate intersection and save as a potential start point 5487 // for an internal loop. 5488 PointString<SubrangePoint> thisStartPoints = PointString.newInstance(); 5489 PointString<SubrangePoint> thatStartPoints = PointString.newInstance(); 5490 PointString<SubrangePoint> pntLst = PointString.newInstance(); 5491 for (SubrangePoint thisPnt : _appInteriorPnts) { 5492 SubrangePoint thatPnt = _thatSrf.getClosest(thisPnt, GTOL100); 5493 relaxPoint(thisPnt, _thatSrf, thatPnt.getParPosition(), pntLst, _tol); 5494 if (!pntLst.isEmpty()) { 5495 thisStartPoints.add(pntLst.getFirst()); 5496 thatStartPoints.add(pntLst.getLast()); 5497 } 5498 pntLst.clear(); 5499 SubrangePoint.recycle(thatPnt); 5500 } 5501 5502 // Remove any points that are near the existing intersection curves. 5503 //System.out.println("Removing redundant interior points..."); 5504 removeClosePoints(thisStartPoints, _tol); 5505 removeClosePoints(thatStartPoints, _tol); 5506 for (SubrangeCurve crv : _thisOut) { 5507 thisStartPoints = removePntsNearCurve(thisStartPoints, crv, tol2); 5508 thatStartPoints = removePntsNearCurve(thatStartPoints, crv, tol2); 5509 } 5510 for (SubrangeCurve crv : _thatOut) { 5511 thisStartPoints = removePntsNearCurve(thisStartPoints, crv, tol2); 5512 thatStartPoints = removePntsNearCurve(thatStartPoints, crv, tol2); 5513 } 5514 5515 // Any remaining interior points are assumed to be on interior loops. 5516 //System.out.println("Tracing from interior points..."); 5517 //System.out.println(" This surface..."); 5518 //System.out.println(" thisStartPoints.size() = " + thisStartPoints.size()); 5519 crvSegs = traceInteriorPoints(thisStartPoints, _thatSrf, _patchSizeTol); 5520 _thisOut.addAll(crvSegs.get(0)); 5521 _thatOut.addAll(crvSegs.get(1)); 5522 5523 // Remove any points near the existing intersection curves. 5524 GeomList<SubrangeCurve> segLst = crvSegs.getFirst(); 5525 for (SubrangeCurve seg : segLst) { 5526 thatStartPoints = removePntsNearCurve(thatStartPoints, seg, tol2); 5527 } 5528 segLst = crvSegs.getLast(); 5529 for (SubrangeCurve seg : segLst) { 5530 thatStartPoints = removePntsNearCurve(thatStartPoints, seg, tol2); 5531 } 5532 5533 //System.out.println(" That surface..."); 5534 //System.out.println(" thatStartPoints.size() = " + thatStartPoints.size()); 5535 crvSegs = traceInteriorPoints(thatStartPoints, _thisSrf, _patchSizeTol); 5536 _thatOut.addAll(crvSegs.get(0)); 5537 _thisOut.addAll(crvSegs.get(1)); 5538 5539 } 5540 5541 /** 5542 * Trace intersection curves that start/end on interior points. 5543 * 5544 * @param thisStartPnts A list of potential start points on "thisSrf" for 5545 * intersection tracing. 5546 * @param thatSrf "That" surface being used for intersecting with "this" 5547 * surface. 5548 * @param tol The tolerance to use when deciding if a point belongs to 5549 * an existing intersection curve. 5550 * @return A list of lists of intersection curves traced from/to each edge point 5551 * on each surface starting with the edge points on "this" surface. 5552 * @throws RootException If there is a convergence problem with a root finder. 5553 */ 5554 private GeomList<GeomList<SubrangeCurve>> traceInteriorPoints(PointString<SubrangePoint> thisStartPnts, 5555 AbstractSurface thatSrf, Parameter<Length> tol) throws RootException { 5556 5557 GeomList<GeomList<SubrangeCurve>> output = GeomList.newInstance(); 5558 GeomList<SubrangeCurve> thisSegs = GeomList.newInstance(); 5559 output.add(thisSegs); 5560 GeomList<SubrangeCurve> thatSegs = GeomList.newInstance(); 5561 output.add(thatSegs); 5562 5563 while (!thisStartPnts.isEmpty()) { 5564 SubrangePoint startPnt = thisStartPnts.get(0); 5565 5566 // Trace the intersection starting at "startPnt". 5567 GeomList<SubrangeCurve> segs = traceSegment(startPnt, thatSrf); 5568 if (segs.containsGeometry()) { 5569 SubrangeCurve thisSeg = segs.get(0); 5570 SubrangeCurve thatSeg = segs.get(1); 5571 5572 // Don't add degenerate segments to the output. 5573 if (!thisSeg.isDegenerate(_tol)) 5574 thisSegs.add(thisSeg); 5575 if (!thatSeg.isDegenerate(_tol)) 5576 thatSegs.add(thatSeg); 5577 5578 // Remove any points that are close to this curve segment. 5579 thisStartPnts.remove(0); // The start point obviously. 5580 thisStartPnts = removePntsNearCurve(thisStartPnts, thisSeg, tol); 5581 5582 } else { 5583 thisStartPnts.remove(0); 5584 } 5585 } 5586 5587 return output; 5588 } 5589 5590 /** 5591 * Method that finds all the exact intersections of the edges of the surface (if 5592 * any). 5593 * 5594 * @return A list of subrange intersection points found on each edge of the 5595 * surface. 5596 */ 5597 private PointString<SubrangePoint> getEdgePoints(Surface srf1, Surface srf2) { 5598 5599 PointString<SubrangePoint> output = PointString.newInstance(); 5600 Parameter<Length> tol = _tol.divide(10); 5601 5602 GeomList<PointString<SubrangePoint>> edge = srf1.getS0Curve().intersect(srf2, tol); 5603 PointString<SubrangePoint> edgePnts = edge.get(0); 5604 int size = edgePnts.size(); 5605 for (int i = 0; i < size; ++i) { 5606 SubrangePoint cPnt = edgePnts.get(i); 5607 double t = cPnt.getParPosition().getValue(0); 5608 SubrangePoint sPnt = srf1.getPoint(0, t); 5609 output.add(sPnt); 5610 SubrangePoint.recycle(cPnt); 5611 } 5612 PointString.recycle(edge.get(0)); 5613 PointString.recycle(edge.get(1)); 5614 GeomList.recycle(edge); 5615 5616 edge = srf1.getT0Curve().intersect(srf2, tol); 5617 edgePnts = edge.get(0); 5618 size = edgePnts.size(); 5619 for (int i = 0; i < size; ++i) { 5620 SubrangePoint cPnt = edgePnts.get(i); 5621 double s = cPnt.getParPosition().getValue(0); 5622 SubrangePoint sPnt = srf1.getPoint(s, 0); 5623 output.add(sPnt); 5624 SubrangePoint.recycle(cPnt); 5625 } 5626 PointString.recycle(edge.get(0)); 5627 PointString.recycle(edge.get(1)); 5628 GeomList.recycle(edge); 5629 5630 edge = srf1.getS1Curve().intersect(srf2, tol); 5631 edgePnts = edge.get(0); 5632 size = edgePnts.size(); 5633 for (int i = 0; i < size; ++i) { 5634 SubrangePoint cPnt = edgePnts.get(i); 5635 double t = cPnt.getParPosition().getValue(0); 5636 SubrangePoint sPnt = srf1.getPoint(1, t); 5637 output.add(sPnt); 5638 SubrangePoint.recycle(cPnt); 5639 } 5640 PointString.recycle(edge.get(0)); 5641 PointString.recycle(edge.get(1)); 5642 GeomList.recycle(edge); 5643 5644 edge = srf1.getT1Curve().intersect(srf2, tol); 5645 edgePnts = edge.get(0); 5646 size = edgePnts.size(); 5647 for (int i = 0; i < size; ++i) { 5648 SubrangePoint cPnt = edgePnts.get(i); 5649 double s = cPnt.getParPosition().getValue(0); 5650 SubrangePoint sPnt = srf1.getPoint(s, 1); 5651 output.add(sPnt); 5652 SubrangePoint.recycle(cPnt); 5653 } 5654 PointString.recycle(edge.get(0)); 5655 PointString.recycle(edge.get(1)); 5656 GeomList.recycle(edge); 5657 5658 // Remove any redundant points 5659 removeClosePoints(output, _tol); 5660 5661 return output; 5662 } 5663 5664 /** 5665 * Trace intersection curves that start/end at the edges of a surface. 5666 * 5667 * @param thisEdgePnts The edge points on "this" surface to use as starting points 5668 * for tracing intersection curves. 5669 * @param thatEdgePnts The edge points on "that" surface to use as starting points 5670 * for tracing intersection curves. 5671 * @param thatSrf "That" surface being used for intersecting with "this" 5672 * surface. 5673 * @return A list of lists of intersection curves traced from/to each edge point 5674 * on each surface starting with the edge points on "this" surface. 5675 * @throws RootException If there is a convergence problem with a root finder. 5676 */ 5677 private GeomList<GeomList<SubrangeCurve>> traceEdgePoints(PointString<SubrangePoint> thisEdgePnts, 5678 PointString<SubrangePoint> thatEdgePnts, AbstractSurface thatSrf) throws RootException { 5679 5680 Parameter<Length> tol = _tol.times(10); 5681 GeomList<GeomList<SubrangeCurve>> output = GeomList.newInstance(); 5682 GeomList<SubrangeCurve> thisSegs = GeomList.newInstance(); 5683 output.add(thisSegs); 5684 GeomList<SubrangeCurve> thatSegs = GeomList.newInstance(); 5685 output.add(thatSegs); 5686 5687 int thisNumPnts = thisEdgePnts.size(); 5688 while (thisNumPnts > 0) { 5689 SubrangePoint startPnt = thisEdgePnts.get(0); 5690 5691 // Trace the intersection starting at "startPnt". 5692 GeomList<SubrangeCurve> segs = traceSegment(startPnt, thatSrf); 5693 if (segs.containsGeometry()) { 5694 SubrangeCurve thisSeg = segs.get(0); 5695 SubrangeCurve thatSeg = segs.get(1); 5696 5697 // Don't add degenerate segments to the output. 5698 if (!thisSeg.isDegenerate(_tol)) 5699 thisSegs.add(thisSeg); 5700 if (!thatSeg.isDegenerate(_tol)) 5701 thatSegs.add(thatSeg); 5702 5703 // Remove the end point from the lists (if it is in them). 5704 Surface thisSrf = (Surface)startPnt.getChild(); 5705 SubrangePoint endPoint = thisSrf.getPoint(thisSeg.getParPosition().getRealPoint(1)); 5706 int idx = findParNearestInList(endPoint, thisEdgePnts); 5707 if (idx > 0) { 5708 SubrangePoint pidx = thisEdgePnts.get(idx); 5709 Parameter<Length> d = pidx.distance(endPoint); 5710 if (d.isLessThan(tol)) { 5711 thisEdgePnts.remove(idx); 5712 --thisNumPnts; 5713 } 5714 } 5715 if (!thatEdgePnts.isEmpty()) { 5716 endPoint = thatSrf.getPoint(thatSeg.getParPosition().getRealPoint(1)); 5717 idx = findParNearestInList(endPoint, thatEdgePnts); 5718 if (idx >= 0) { 5719 SubrangePoint pidx = thatEdgePnts.get(idx); 5720 Parameter<Length> d = pidx.distance(endPoint); 5721 if (d.isLessThan(tol)) { 5722 thatEdgePnts.remove(idx); 5723 } 5724 } 5725 } 5726 } 5727 GeomList.recycle(segs); 5728 5729 // Remove the start point from this list. 5730 thisEdgePnts.remove(0); 5731 --thisNumPnts; 5732 } // end while 5733 5734 return output; 5735 } 5736 5737 /** 5738 * Method that traces an intersection segment from the given start point until an 5739 * edge of either intersecting surface is encountered or the start point is 5740 * re-encountered. 5741 * 5742 * @param startPnt The point on "this" surface to start tracing the intersection 5743 * from. 5744 * @param thatSrf "That" surface that is being intersected with "this" surface. 5745 * @return A pair of subrange curves representing this intersection segment on 5746 * each surface (this & that). 5747 */ 5748 private GeomList<SubrangeCurve> traceSegment(SubrangePoint startPnt, AbstractSurface thatSrf) throws RootException { 5749 AbstractSurface thisSrf = (AbstractSurface)startPnt.getChild(); 5750 5751 BasicNurbsCurve pCrv1, pCrv2; 5752 StackContext.enter(); 5753 try { 5754 PointString<SubrangePoint> thisTracedPnts = PointString.newInstance(); 5755 PointString<SubrangePoint> thatTracedPnts = PointString.newInstance(); 5756 PointString<SubrangePoint> pntLst = PointString.newInstance(); 5757 5758 Unit<Length> unit = thisSrf.getUnit(); 5759 SubrangePoint p1 = startPnt; 5760 thisTracedPnts.add(p1); 5761 GeomPoint thisST = p1.getParPosition(); // Parametric position on the surface of the point p1. 5762 double s01 = thisST.getValue(0); 5763 double t01 = thisST.getValue(1); 5764 5765 SubrangePoint thatStartPnt = thatSrf.getClosest(p1.copyToReal(), GTOL); 5766 thatTracedPnts.add(thatStartPnt); 5767 GeomPoint thatST = thatStartPnt.getParPosition(); // Parametric position on the surface of the point thatP1. 5768 double s2old = thatST.getValue(0); 5769 double t2old = thatST.getValue(1); 5770 5771 // Is the start point an edge point? 5772 boolean startEdge = false; 5773 if (parNearEnds(s01, GTOL100) || parNearEnds(t01, GTOL100)) { 5774 startEdge = true; 5775 } 5776 5777 boolean firstPoint = true; 5778 int collapsedEdge = -1; 5779 boolean flipNormal = false; 5780 while (thisTracedPnts.size() < 1000) { 5781 // Determine the step direction. 5782 Vector<Dimensionless> pu = thisSrf.getSDerivative(thisST.getValue(0), 5783 thisST.getValue(1), 1, false).toUnitVector(); 5784 if (!pu.isValid()) { 5785 // Step away from a collapsed edge. 5786 double t = thisST.getValue(1); 5787 if (t < 0.5) { 5788 thisST = Point.valueOf(SI.METER, thisST.getValue(0), GTOL100); 5789 collapsedEdge = 2; 5790 } else { 5791 thisST = Point.valueOf(SI.METER, thisST.getValue(0), 1. - GTOL100); 5792 collapsedEdge = 3; 5793 } 5794 if (startEdge) { 5795 p1 = thisSrf.getPoint(thisST); 5796 s01 = thisST.getValue(0); 5797 t01 = thisST.getValue(1); 5798 } 5799 pu = thisSrf.getSDerivative(thisST.getValue(0), thisST.getValue(1), 1, false).toUnitVector(); 5800 } 5801 Vector<Dimensionless> pv = thisSrf.getTDerivative(thisST.getValue(0), 5802 thisST.getValue(1), 1, false).toUnitVector(); 5803 if (!pv.isValid()) { 5804 // Step away from a collapsed edge. 5805 double s = thisST.getValue(0); 5806 if (s < 0.5) { 5807 thisST = Point.valueOf(SI.METER, GTOL100, thisST.getValue(1)); 5808 collapsedEdge = 0; 5809 } else { 5810 thisST = Point.valueOf(SI.METER, 1. - GTOL100, thisST.getValue(1)); 5811 collapsedEdge = 1; 5812 } 5813 pu = pv = thisSrf.getTDerivative(thisST.getValue(0), thisST.getValue(1), 1, false).toUnitVector(); 5814 if (startEdge) { 5815 p1 = thisSrf.getPoint(thisST); 5816 s01 = thisST.getValue(0); 5817 t01 = thisST.getValue(1); 5818 } 5819 } 5820 Vector<Dimensionless> n2 = thatSrf.getNormal(thatST); 5821 if (flipNormal) 5822 n2 = n2.opposite(); 5823 Parameter pudotn2 = pu.dot(n2); 5824 Parameter pvdotn2 = pv.dot(n2); 5825 if (pudotn2.isApproxZero() && pvdotn2.isApproxZero()) { 5826 Logger.getLogger(AbstractSurface.class.getName()).log(Level.WARNING, 5827 "Surfaces are nearly tangent to each other in traceSegment()!"); 5828 break; // Surface is nearly tangent to other surface. TODO: Handle this better. 5829 } 5830 // Ti = pv*(pu dot n2) - pu*(pv dot n2) 5831 Vector<Dimensionless> Ti = pv.times(pudotn2).minus((Vector)pu.times(pvdotn2)); 5832 Ti = Ti.toUnitVector(); 5833 5834 // Determine the step size "L". 5835 Point p = p1.plus(Point.valueOf(Ti.times(_eps1))); 5836 SubrangePoint x = thisSrf.getClosestFarthestInterior(p, thisST.getValue(0), thisST.getValue(1), GTOL, true); 5837 5838 // If "Ti" initially steps us off the edge of the surface, reverse it. 5839 if (startEdge) { 5840 startEdge = false; 5841 double sx = x.getParPosition().getValue(0); 5842 double tx = x.getParPosition().getValue(1); 5843 if (parNearEnds(sx, GTOL) || parNearEnds(tx, GTOL) 5844 || abs(sx - s01) > 0.3 || abs(tx - t01) > 0.3) { 5845 // Ti stepped off the edge. 5846 Ti = Ti.opposite(); 5847 p = p1.plus(Point.valueOf(Ti.times(_eps1))); 5848 x = thisSrf.getClosestFarthestInterior(p, thisST.getValue(0), thisST.getValue(1), GTOL100, true); 5849 5850 // Reverse the intersecting surface normal so this doesn't happen again. 5851 flipNormal = true; 5852 } 5853 } 5854 5855 relaxPoint(x, thatSrf, thatST, pntLst, _tol); 5856 if (pntLst.isEmpty()) 5857 break; 5858 x = pntLst.get(0); 5859 pntLst.clear(); 5860 p = p1.plus(Point.valueOf(Ti.times(_eps2))); 5861 SubrangePoint y = thisSrf.getClosestFarthestInterior(p, thisST.getValue(0), thisST.getValue(1), GTOL100, true); 5862 relaxPoint(y, thatSrf, thatST, pntLst, _tol); 5863 if (pntLst.isEmpty()) 5864 break; 5865 y = pntLst.get(0); 5866 pntLst.clear(); 5867 5868 // If the start point was on a collapsed edge, go back and fix things up. 5869 if (collapsedEdge >= 0) { 5870 thisST = x.getParPosition(); 5871 switch (collapsedEdge) { 5872 case 0: 5873 thisST = Point.valueOf(SI.METER, 0, thisST.getValue(1)); 5874 break; 5875 case 1: 5876 thisST = Point.valueOf(SI.METER, 1, thisST.getValue(1)); 5877 break; 5878 case 2: 5879 thisST = Point.valueOf(SI.METER, thisST.getValue(0), 0); 5880 break; 5881 case 3: 5882 thisST = Point.valueOf(SI.METER, thisST.getValue(0), 1); 5883 break; 5884 } 5885 p1 = thisSrf.getPoint(thisST); 5886 s01 = thisST.getValue(0); 5887 t01 = thisST.getValue(1); 5888 startPnt = p1; 5889 thisTracedPnts.clear(); 5890 thisTracedPnts.add(p1); 5891 5892 collapsedEdge = -1; 5893 } 5894 5895 Vector<Length> a = x.minus(p1).toGeomVector(); 5896 Vector<Length> b = y.minus(p1).toGeomVector(); 5897 5898 double na = a.normValue(); 5899 double nb = b.normValue(); 5900 double namb = a.minus(b).normValue(); 5901 double naxb = a.cross(b).normValue(); 5902// System.out.print("naxb = " + naxb); 5903 Parameter<Length> L; 5904 if (naxb > 1e-9) { 5905 // The 3 points are not co-linear. 5906 // r = |a||b||a - b|/(2*|a x b|) 5907 double r = 0.5 * na * nb * namb / naxb; 5908 L = Parameter.valueOf(r * DTHETA, unit); 5909// System.out.print(", r*DTHETA = " + L); 5910 if (L.isGreaterThan(_epsCRT)) 5911 L = _epsCRT; 5912 } else { 5913 L = _epsCRT; 5914 } 5915// System.out.println(", L = " + L); 5916 5917 // Calculate the new intersection point locations. 5918 p = p1.plus(Point.valueOf(Ti.times(L))); 5919 SubrangePoint p2 = thisSrf.getClosestFarthestInterior(p, thisST.getValue(0), thisST.getValue(1), GTOL100, true); 5920 relaxPoint(p2, thatSrf, thatST, pntLst, _tol); 5921 if (pntLst.isEmpty()) 5922 break; 5923 p2 = pntLst.get(0); 5924 SubrangePoint thatP2 = pntLst.get(1); 5925 pntLst.clear(); 5926 5927 // See if the new point on thatSrf has jumped a boundary. 5928 thatST = thatP2.getParPosition(); 5929 double s2 = thatST.getValue(0); 5930 double t2 = thatST.getValue(1); 5931 if ((s2 > 0.8 && s2old < 0.2) || (s2 < 0.2 && s2old > 0.8)) { 5932 // We have crossed the S boundary. 5933 s2 = round(s2old); 5934 thatP2 = thatSrf.getPoint(s2, t2); 5935 5936 // Relax the corrected point and the nearest thisSrf point to the intersection. 5937 relaxPoint(thatP2, thisSrf, thisST, pntLst, _tol); 5938 if (pntLst.isEmpty()) 5939 break; 5940 thatP2 = pntLst.get(0); 5941 p2 = pntLst.get(1); 5942 pntLst.clear(); 5943 5944 thatST = thatP2.getParPosition(); 5945 5946 } else if ((t2 > 0.8 && t2old < 0.2) || (t2 < 0.2 && t2old > 0.8)) { 5947 // We have crossed the T boundary. 5948 t2 = round(t2old); 5949 thatP2 = thatSrf.getPoint(s2, t2); 5950 relaxPoint(thatP2, thisSrf, thisST, pntLst, _tol); 5951 if (pntLst.isEmpty()) 5952 break; 5953 thatP2 = pntLst.get(0); 5954 p2 = pntLst.get(1); 5955 pntLst.clear(); 5956 thatST = thatP2.getParPosition(); 5957 } 5958 5959 Parameter<Length> d = p2.distance(thatP2); 5960 if (d.isLessThan(_tol)) { 5961 // Store the new intersection points. 5962 thisTracedPnts.add(p2); 5963 thatTracedPnts.add(thatP2); 5964 } 5965 5966 // Is this new point an edge point on "this" surface? 5967 thisST = p2.getParPosition(); 5968 double s1 = thisST.getValue(0); 5969 double t1 = thisST.getValue(1); 5970 if (!firstPoint) { 5971 if (parNearEnds(s1, GTOL100) || parNearEnds(t1, GTOL100)) 5972 break; 5973 5974 // Is this new point an edge point on "that" surface? 5975 if (parNearEnds(s2, GTOL100) || parNearEnds(t2, GTOL100)) 5976 break; 5977 5978 // Is this new point near the start point? 5979 if (thisST.distanceValue(startPnt.getParPosition()) < 0.1 5980 && p2.distance(startPnt).isLessThan(L.times(0.75))) { 5981 thisTracedPnts.add(startPnt); // Repeat the 1st point. 5982 thatTracedPnts.add(thatStartPnt); 5983 break; 5984 } 5985 } 5986 5987 // Prepare for next step. 5988 p1 = p2; 5989 s2old = s2; 5990 t2old = t2; 5991 firstPoint = false; 5992 } 5993 5994 // Create the curve segment on each surface. 5995 if (thisTracedPnts.size() == 1) { 5996 thisST = thisTracedPnts.get(0).getParPosition(); 5997 pCrv1 = StackContext.outerCopy(CurveFactory.createPoint(1, thisST)); 5998 5999 thatST = thatTracedPnts.get(0).getParPosition(); 6000 pCrv2 = StackContext.outerCopy(CurveFactory.createPoint(1, thatST)); 6001 6002 } else { 6003 // Create the intersection curve from the list of intersection points. 6004 SubrangeCurve thisCrv = makeSubrangeCurve(thisTracedPnts); 6005 SubrangeCurve thatCrv = makeSubrangeCurve(thatTracedPnts); 6006 6007 // Copy the parametric curve's to the outer context. 6008 pCrv1 = StackContext.outerCopy((BasicNurbsCurve)thisCrv.getParPosition()); 6009 pCrv2 = StackContext.outerCopy((BasicNurbsCurve)thatCrv.getParPosition()); 6010 } 6011 6012 } finally { 6013 StackContext.exit(); 6014 } 6015 6016 // Create the output list and add the final subrange curves to it. 6017 GeomList<SubrangeCurve> output = GeomList.newInstance(); 6018 output.add(SubrangeCurve.newInstance(thisSrf, pCrv1)); 6019 output.add(SubrangeCurve.newInstance(thatSrf, pCrv2)); 6020 6021 return output; 6022 } 6023 6024 private static final int MAXIT = 500; 6025 6026 /** 6027 * Relax the approximate intersection point toward the exact intersection until it 6028 * is within "tol" of the exact intersection. 6029 * 6030 * @param p1 The approximate surface intersection point being refined or 6031 * relaxed. 6032 * @param thatSrf The intersecting surface. 6033 * @param thatST The parametric position on "thatSrf" near "p1". 6034 * @param tol The tolerance required of the intersection points. 6035 * @param output A list to which the input point relaxed onto the intersection 6036 * between the two surfaces to within the tolerance "tol" will be 6037 * added. The 1st point is attached to "this" surface and the 2nd 6038 * point is attached to "that" surface. If the surfaces are locally 6039 * tangent near the intersection, this may add nothing to the 6040 * output list. 6041 */ 6042 private void relaxPoint(SubrangePoint p1, AbstractSurface thatSrf, GeomPoint thatST, 6043 PointString<SubrangePoint> output, 6044 Parameter<Length> tol) throws RootException { 6045 6046 AbstractSurface thisSrf = (AbstractSurface)p1.getChild(); 6047 6048 Point parPnt1, parPnt2; 6049 StackContext.enter(); 6050 try { 6051 SubrangePoint p2 = thatSrf.getClosestFarthestInterior(p1, thatST.getValue(0), thatST.getValue(1), GTOL100, true); 6052 6053 Parameter<Length> d = p1.distance(p2); 6054 if (d.isGreaterThan(tol)) { 6055 6056 Unit<Length> unit = thisSrf.getUnit(); 6057 int numDims = thisSrf.getPhyDimension(); 6058 MutableVector Lu = MutableVector.newInstance(numDims, Dimensionless.UNIT); 6059 MutablePoint L0 = MutablePoint.newInstance(numDims, unit); 6060 6061 int iteration = 0; 6062 do { 6063 // Get the local tangent plane on this surface. 6064 GeomPoint thisST = p1.getParPosition(); 6065 GeomPlane T1 = thisSrf.getTangentPlane(thisST); 6066 6067 // Get the local tangent plane on that surface. 6068 thatST = p2.getParPosition(); 6069 GeomPlane T2 = thatSrf.getTangentPlane(thatST); 6070 6071 // Intersect the surface tangent planes. 6072 IntersectType type = GeomUtil.planePlaneIntersect(T1, T2, L0, Lu); 6073 if (type != IntersectType.INTERSECT) { 6074 return; 6075 } 6076 6077 // Find the closest point on the plane intersection line to the previous points. 6078 Point q1 = GeomUtil.pointLineClosest(p1, L0, Lu); 6079 Point q2 = GeomUtil.pointLineClosest(p2, L0, Lu); 6080 Point pip1 = q1.plus(q2).divide(2); 6081 6082 // Re-attach this point to the surfaces. 6083 p1 = thisSrf.getClosestFarthestInterior(pip1, thisST.getValue(0), thisST.getValue(1), GTOL100, true); 6084 p2 = thatSrf.getClosestFarthestInterior(pip1, thatST.getValue(0), thatST.getValue(1), GTOL100, true); 6085 6086 // Update the distance between these new points. 6087 d = p1.distance(p2); 6088 //System.out.println("it = " + iteration + ", d = " + d); 6089 6090 ++iteration; 6091 6092 } while (d.isGreaterThan(tol) && iteration < MAXIT); 6093 6094 if (iteration == MAXIT) 6095 Logger.getLogger(AbstractSurface.class.getName()).log(Level.WARNING, 6096 "Convergence problem in SrfSrfIntersect.relaxPoint()."); 6097 } 6098 6099 // Copy the parametric positions of the final points to the outer context. 6100 parPnt1 = StackContext.outerCopy(p1.getParPosition().immutable()); 6101 parPnt2 = StackContext.outerCopy(p2.getParPosition().immutable()); 6102 6103 } finally { 6104 StackContext.exit(); 6105 } 6106 6107 // Store the relaxed points. 6108 output.add(thisSrf.getPoint(parPnt1)); 6109 output.add(thatSrf.getPoint(parPnt2)); 6110 } 6111 6112 /** 6113 * Uses a recursive "Divide and Conquer" approach to intersecting a surface patch 6114 * with another surface patch. On each call, the following happens: 6115 * <pre> 6116 * 1) Each patch is checked to see if it is approx. a flat plane. 6117 * 1b) If both are, then a plane-quad intersection is performed, the 6118 * approximate intersection points added to the _appInteriorPnts list 6119 * and the method exited. 6120 * 2) Each non-planar patch is divided into four quadrant patches. 6121 * 2b) Each sub-patch is tested to see if it's bounding box is intersected by any 6122 * of the other sub-patches. 6123 * If it is, then this method is called with those patches recursively. 6124 * Otherwise, the method is exited. 6125 * </pre> A number of class variables are used to pass information to this 6126 * recursion: 6127 * <pre> 6128 * _thisSrf The full surfaces intersections are being found for. 6129 * _thatSrf The full surface being intersected with this surface. 6130 * _tol The tolerance to use in determining if the geometry is in tolerance. 6131 * _appInteriorPnts A list used to collect the approximate subrange intersection 6132 * points on _thisSrf. 6133 * </pre> 6134 * 6135 * @param patch1 The first surface patch being checked for intersection with 6136 * "that" surface. 6137 * @param s0_1 The "s" parametric position of the start of patch1 on the master 6138 * "this" surface. 6139 * @param t0_1 The "t" parametric position of the start of patch1 on the master 6140 * "this" surface. 6141 * @param s1_1 The parametric "s" position of the end of patch1 on the master 6142 * "this" surface. 6143 * @param t1_1 The parametric "t" position of the end of patch1 on the master 6144 * "this" surface. 6145 * @param patch2 The 2nd surface patch being checked for intersection with "this" 6146 * surface. 6147 * @param s0_2 The "s" parametric position of the start of patch2 on the master 6148 * "that" surface. 6149 * @param t0_2 The "t" parametric position of the start of patch2 on the master 6150 * "that" surface. 6151 * @param s1_2 The parametric "s" position of the end of patch2 on the master 6152 * "that" surface. 6153 * @param t1_2 The parametric "t" position of the end of patch2 on the master 6154 * "that" surface. 6155 */ 6156 private void divideAndConquerSrf(int depth, Surface patch1, double s0_1, double t0_1, double s1_1, double t1_1, 6157 Surface patch2, double s0_2, double t0_2, double s1_2, double t1_2) throws RootException { 6158 6159 boolean p1Planar = patch1.isPlanar(EPS_FT, EPS_FT); // Using a fast method that unfortunately ignores "tol". 6160 boolean p2Planar = patch2.isPlanar(EPS_FT, EPS_FT); 6161 6162 boolean p1IsSmall = isSmall(patch1, _patchSizeTol); 6163 boolean p2IsSmall = isSmall(patch2, _patchSizeTol); 6164 6165 // Check to see if both patches are nearly planar or to small. 6166 if (depth > MAXDEPTH || (p1Planar && p2Planar) || (p1IsSmall && p2IsSmall)) { 6167 6168 // Form a quadrilateral panel (2 triangles) from the corners of patch1. 6169 Point A1 = patch1.getRealPoint(0, 0); 6170 Point B1 = patch1.getRealPoint(1, 0); 6171 Point C1 = patch1.getRealPoint(1, 1); 6172 Point D1 = patch1.getRealPoint(0, 1); 6173 6174 // Form a quadrilateral panel (2 triangles) from the corners of patch2. 6175 Point A2 = patch2.getRealPoint(0, 0); 6176 Point B2 = patch2.getRealPoint(1, 0); 6177 Point C2 = patch2.getRealPoint(1, 1); 6178 Point D2 = patch2.getRealPoint(0, 1); 6179 6180 // Determine if any edge of patch2 intersects any triangle of patch 1. 6181 // Triangle #1 = A1, B1, C1; Triangle #2 = A1, C1, D1 6182 if (GeomUtil.lineSegTriIntersect(A2, B2, A1, B1, C1, null) == IntersectType.INTERSECT 6183 || GeomUtil.lineSegTriIntersect(A2, B2, A1, C1, D1, null) == IntersectType.INTERSECT 6184 || GeomUtil.lineSegTriIntersect(B2, C2, A1, B1, C1, null) == IntersectType.INTERSECT 6185 || GeomUtil.lineSegTriIntersect(B2, C2, A1, C1, D1, null) == IntersectType.INTERSECT 6186 || GeomUtil.lineSegTriIntersect(C2, D2, A1, B1, C1, null) == IntersectType.INTERSECT 6187 || GeomUtil.lineSegTriIntersect(C2, D2, A1, C1, D1, null) == IntersectType.INTERSECT 6188 || GeomUtil.lineSegTriIntersect(D2, A2, A1, B1, C1, null) == IntersectType.INTERSECT 6189 || GeomUtil.lineSegTriIntersect(D2, A2, A1, C1, D1, null) == IntersectType.INTERSECT) { 6190 6191 // Store off the approximate intersection point (patch1 center). 6192 double s = segmentPos2Parent(0.5, s0_1, s1_1); 6193 double t = segmentPos2Parent(0.5, t0_1, t1_1); 6194 SubrangePoint centerPnt = _thisSrf.getPoint(s, t); 6195 _appInteriorPnts.add(centerPnt); 6196 } 6197 6198 //System.out.println("depth = " + depth); 6199 } else { 6200 Surface p00_1 = patch1; 6201 Surface p00_2 = patch2; 6202 6203 Surface p01_1 = null, p10_1 = null, p11_1 = null; 6204 if (!p1Planar) { 6205 // Subdivide the 1st patch into 4 quadrants. 6206 GeomList<Surface> leftRight = patch1.splitAtS(0.5); 6207 GeomList<Surface> botTop = leftRight.get(0).splitAtT(0.5); 6208 p00_1 = botTop.get(0); // Lower Left 6209 p01_1 = botTop.get(1); // Top Left 6210 GeomList.recycle(botTop); 6211 botTop = leftRight.get(1).splitAtT(0.5); 6212 p10_1 = botTop.get(0); // Lower Right 6213 p11_1 = botTop.get(1); // Top Right 6214 GeomList.recycle(botTop); 6215 GeomList.recycle(leftRight); 6216 } 6217 6218 Surface p01_2 = null, p10_2 = null, p11_2 = null; 6219 if (!p2Planar) { 6220 // Subdivide the 2nd patch into 4 quadrants. 6221 GeomList<Surface> leftRight = patch2.splitAtS(0.5); 6222 GeomList<Surface> botTop = leftRight.get(0).splitAtT(0.5); 6223 p00_2 = botTop.get(0); // Lower Left 6224 p01_2 = botTop.get(1); // Top Left 6225 GeomList.recycle(botTop); 6226 botTop = leftRight.get(1).splitAtT(0.5); 6227 p10_2 = botTop.get(0); // Lower Right 6228 p11_2 = botTop.get(1); // Top Right 6229 GeomList.recycle(botTop); 6230 GeomList.recycle(leftRight); 6231 } 6232 6233 // Mean parametric positions on each patch. 6234 double sm_1 = 0.5 * (s0_1 + s1_1); 6235 double tm_1 = 0.5 * (t0_1 + t1_1); 6236 double sm_2 = 0.5 * (s0_2 + s1_2); 6237 double tm_2 = 0.5 * (t0_2 + t1_2); 6238 6239 // Check for possible intersections on the lower-left patch. 6240 if (GeomUtil.boundsIntersect(p00_1, p00_2)) { 6241 // May be an intersection. 6242 double sl_1 = s0_1; 6243 double sh_1 = sm_1; 6244 double tl_1 = t0_1; 6245 double th_1 = tm_1; 6246 double sl_2 = s0_2; 6247 double sh_2 = sm_2; 6248 double tl_2 = t0_2; 6249 double th_2 = tm_2; 6250 6251 // Recurse down to a finer level of detail. 6252 divideAndConquerSrf(depth + 1, p00_1, sl_1, tl_1, sh_1, th_1, 6253 p00_2, sl_2, tl_2, sh_2, th_2); 6254 } 6255 if (!p2Planar) { 6256 if (GeomUtil.boundsIntersect(p00_1, p10_2)) { 6257 // May be an intersection. 6258 double sl_1 = s0_1; 6259 double sh_1 = sm_1; 6260 double tl_1 = t0_1; 6261 double th_1 = tm_1; 6262 double sl_2 = sm_2; 6263 double sh_2 = s1_2; 6264 double tl_2 = t0_2; 6265 double th_2 = tm_2; 6266 6267 // Recurse down to a finer level of detail. 6268 divideAndConquerSrf(depth + 1, p00_1, sl_1, tl_1, sh_1, th_1, 6269 p10_2, sl_2, tl_2, sh_2, th_2); 6270 } 6271 if (GeomUtil.boundsIntersect(p00_1, p01_2)) { 6272 // May be an intersection. 6273 double sl_1 = s0_1; 6274 double sh_1 = sm_1; 6275 double tl_1 = t0_1; 6276 double th_1 = tm_1; 6277 double sl_2 = s0_2; 6278 double sh_2 = sm_2; 6279 double tl_2 = tm_2; 6280 double th_2 = t1_2; 6281 6282 // Recurse down to a finer level of detail. 6283 divideAndConquerSrf(depth + 1, p00_1, sl_1, tl_1, sh_1, th_1, 6284 p01_2, sl_2, tl_2, sh_2, th_2); 6285 } 6286 if (GeomUtil.boundsIntersect(p00_1, p11_2)) { 6287 // May be an intersection. 6288 double sl_1 = s0_1; 6289 double sh_1 = sm_1; 6290 double tl_1 = t0_1; 6291 double th_1 = tm_1; 6292 double sl_2 = sm_2; 6293 double sh_2 = s1_2; 6294 double tl_2 = tm_2; 6295 double th_2 = t1_2; 6296 6297 // Recurse down to a finer level of detail. 6298 divideAndConquerSrf(depth + 1, p00_1, sl_1, tl_1, sh_1, th_1, 6299 p11_2, sl_2, tl_2, sh_2, th_2); 6300 } 6301 } // if !p2Planar 6302 if (p00_1 != patch1 && p00_1 instanceof BasicNurbsSurface) 6303 BasicNurbsSurface.recycle((BasicNurbsSurface)p00_1); 6304 6305 if (!p1Planar) { 6306 // Check for possible intersections on the lower-right patch. 6307 if (GeomUtil.boundsIntersect(p10_1, p00_2)) { 6308 // May be an intersection. 6309 double sl_1 = sm_1; 6310 double sh_1 = s1_1; 6311 double tl_1 = t0_1; 6312 double th_1 = tm_1; 6313 double sl_2 = s0_2; 6314 double sh_2 = sm_2; 6315 double tl_2 = t0_2; 6316 double th_2 = tm_2; 6317 6318 // Recurse down to a finer level of detail. 6319 divideAndConquerSrf(depth + 1, p10_1, sl_1, tl_1, sh_1, th_1, 6320 p00_2, sl_2, tl_2, sh_2, th_2); 6321 } 6322 if (!p2Planar) { 6323 if (GeomUtil.boundsIntersect(p10_1, p10_2)) { 6324 // May be an intersection. 6325 double sl_1 = sm_1; 6326 double sh_1 = s1_1; 6327 double tl_1 = t0_1; 6328 double th_1 = tm_1; 6329 double sl_2 = sm_2; 6330 double sh_2 = s1_2; 6331 double tl_2 = t0_2; 6332 double th_2 = tm_2; 6333 6334 // Recurse down to a finer level of detail. 6335 divideAndConquerSrf(depth + 1, p10_1, sl_1, tl_1, sh_1, th_1, 6336 p10_2, sl_2, tl_2, sh_2, th_2); 6337 } 6338 if (GeomUtil.boundsIntersect(p10_1, p01_2)) { 6339 // May be an intersection. 6340 double sl_1 = sm_1; 6341 double sh_1 = s1_1; 6342 double tl_1 = t0_1; 6343 double th_1 = tm_1; 6344 double sl_2 = s0_2; 6345 double sh_2 = sm_2; 6346 double tl_2 = tm_2; 6347 double th_2 = t1_2; 6348 6349 // Recurse down to a finer level of detail. 6350 divideAndConquerSrf(depth + 1, p10_1, sl_1, tl_1, sh_1, th_1, 6351 p01_2, sl_2, tl_2, sh_2, th_2); 6352 } 6353 if (GeomUtil.boundsIntersect(p10_1, p11_2)) { 6354 // May be an intersection. 6355 double sl_1 = sm_1; 6356 double sh_1 = s1_1; 6357 double tl_1 = t0_1; 6358 double th_1 = tm_1; 6359 double sl_2 = sm_2; 6360 double sh_2 = s1_2; 6361 double tl_2 = tm_2; 6362 double th_2 = t1_2; 6363 6364 // Recurse down to a finer level of detail. 6365 divideAndConquerSrf(depth + 1, p10_1, sl_1, tl_1, sh_1, th_1, 6366 p11_2, sl_2, tl_2, sh_2, th_2); 6367 } 6368 } // if !p2Planar 6369 if (p10_1 instanceof BasicNurbsSurface) 6370 BasicNurbsSurface.recycle((BasicNurbsSurface)p10_1); 6371 6372 // Check for possible intersections on the top-left patch. 6373 if (GeomUtil.boundsIntersect(p01_1, p00_2)) { 6374 // May be an intersection. 6375 double sl_1 = s0_1; 6376 double sh_1 = sm_1; 6377 double tl_1 = tm_1; 6378 double th_1 = t1_1; 6379 double sl_2 = s0_2; 6380 double sh_2 = sm_2; 6381 double tl_2 = t0_2; 6382 double th_2 = tm_2; 6383 6384 // Recurse down to a finer level of detail. 6385 divideAndConquerSrf(depth + 1, p01_1, sl_1, tl_1, sh_1, th_1, 6386 p00_2, sl_2, tl_2, sh_2, th_2); 6387 } 6388 if (!p2Planar) { 6389 if (GeomUtil.boundsIntersect(p01_1, p10_2)) { 6390 // May be an intersection. 6391 double sl_1 = s0_1; 6392 double sh_1 = sm_1; 6393 double tl_1 = tm_1; 6394 double th_1 = t1_1; 6395 double sl_2 = sm_2; 6396 double sh_2 = s1_2; 6397 double tl_2 = t0_2; 6398 double th_2 = tm_2; 6399 6400 // Recurse down to a finer level of detail. 6401 divideAndConquerSrf(depth + 1, p01_1, sl_1, tl_1, sh_1, th_1, 6402 p10_2, sl_2, tl_2, sh_2, th_2); 6403 } 6404 if (GeomUtil.boundsIntersect(p01_1, p01_2)) { 6405 // May be an intersection. 6406 double sl_1 = s0_1; 6407 double sh_1 = sm_1; 6408 double tl_1 = tm_1; 6409 double th_1 = t1_1; 6410 double sl_2 = s0_2; 6411 double sh_2 = sm_2; 6412 double tl_2 = tm_2; 6413 double th_2 = t1_2; 6414 6415 // Recurse down to a finer level of detail. 6416 divideAndConquerSrf(depth + 1, p01_1, sl_1, tl_1, sh_1, th_1, 6417 p01_2, sl_2, tl_2, sh_2, th_2); 6418 } 6419 if (GeomUtil.boundsIntersect(p01_1, p11_2)) { 6420 // May be an intersection. 6421 double sl_1 = s0_1; 6422 double sh_1 = sm_1; 6423 double tl_1 = tm_1; 6424 double th_1 = t1_1; 6425 double sl_2 = sm_2; 6426 double sh_2 = s1_2; 6427 double tl_2 = tm_2; 6428 double th_2 = t1_2; 6429 6430 // Recurse down to a finer level of detail. 6431 divideAndConquerSrf(depth + 1, p01_1, sl_1, tl_1, sh_1, th_1, 6432 p11_2, sl_2, tl_2, sh_2, th_2); 6433 } 6434 } // if !p2Planar 6435 if (p01_1 instanceof BasicNurbsSurface) 6436 BasicNurbsSurface.recycle((BasicNurbsSurface)p01_1); 6437 6438 // Check for possible intersections on the top-right patch. 6439 if (GeomUtil.boundsIntersect(p11_1, p00_2)) { 6440 // May be an intersection. 6441 double sl_1 = sm_1; 6442 double sh_1 = s1_1; 6443 double tl_1 = tm_1; 6444 double th_1 = t1_1; 6445 double sl_2 = s0_2; 6446 double sh_2 = sm_2; 6447 double tl_2 = t0_2; 6448 double th_2 = tm_2; 6449 6450 // Recurse down to a finer level of detail. 6451 divideAndConquerSrf(depth + 1, p11_1, sl_1, tl_1, sh_1, th_1, 6452 p00_2, sl_2, tl_2, sh_2, th_2); 6453 } 6454 if (p00_2 != patch2 && p00_2 instanceof BasicNurbsSurface) 6455 BasicNurbsSurface.recycle((BasicNurbsSurface)p00_2); 6456 if (!p2Planar) { 6457 if (GeomUtil.boundsIntersect(p11_1, p10_2)) { 6458 // May be an intersection. 6459 double sl_1 = sm_1; 6460 double sh_1 = s1_1; 6461 double tl_1 = tm_1; 6462 double th_1 = t1_1; 6463 double sl_2 = sm_2; 6464 double sh_2 = s1_2; 6465 double tl_2 = t0_2; 6466 double th_2 = tm_2; 6467 6468 // Recurse down to a finer level of detail. 6469 divideAndConquerSrf(depth + 1, p11_1, sl_1, tl_1, sh_1, th_1, 6470 p10_2, sl_2, tl_2, sh_2, th_2); 6471 } 6472 if (p10_2 instanceof BasicNurbsSurface) 6473 BasicNurbsSurface.recycle((BasicNurbsSurface)p10_2); 6474 if (GeomUtil.boundsIntersect(p11_1, p01_2)) { 6475 // May be an intersection. 6476 double sl_1 = sm_1; 6477 double sh_1 = s1_1; 6478 double tl_1 = tm_1; 6479 double th_1 = t1_1; 6480 double sl_2 = s0_2; 6481 double sh_2 = sm_2; 6482 double tl_2 = tm_2; 6483 double th_2 = t1_2; 6484 6485 // Recurse down to a finer level of detail. 6486 divideAndConquerSrf(depth + 1, p11_1, sl_1, tl_1, sh_1, th_1, 6487 p01_2, sl_2, tl_2, sh_2, th_2); 6488 } 6489 if (p01_2 instanceof BasicNurbsSurface) 6490 BasicNurbsSurface.recycle((BasicNurbsSurface)p01_2); 6491 if (GeomUtil.boundsIntersect(p11_1, p11_2)) { 6492 // May be an intersection. 6493 double sl_1 = sm_1; 6494 double sh_1 = s1_1; 6495 double tl_1 = tm_1; 6496 double th_1 = t1_1; 6497 double sl_2 = sm_2; 6498 double sh_2 = s1_2; 6499 double tl_2 = tm_2; 6500 double th_2 = t1_2; 6501 6502 // Recurse down to a finer level of detail. 6503 divideAndConquerSrf(depth + 1, p11_1, sl_1, tl_1, sh_1, th_1, 6504 p11_2, sl_2, tl_2, sh_2, th_2); 6505 } 6506 if (p11_2 instanceof BasicNurbsSurface) 6507 BasicNurbsSurface.recycle((BasicNurbsSurface)p11_2); 6508 } // if !p2Planar 6509 if (p11_1 instanceof BasicNurbsSurface) 6510 BasicNurbsSurface.recycle((BasicNurbsSurface)p11_1); 6511 } // if !p1Planar 6512 } 6513 6514 } 6515 6516 public static void recycle(SrfSrfIntersect instance) { 6517 FACTORY.recycle(instance); 6518 } 6519 6520 private SrfSrfIntersect() { } 6521 6522 private static final ObjectFactory<SrfSrfIntersect> FACTORY = new ObjectFactory<SrfSrfIntersect>() { 6523 @Override 6524 protected SrfSrfIntersect create() { 6525 return new SrfSrfIntersect(); 6526 } 6527 6528 @Override 6529 protected void cleanup(SrfSrfIntersect obj) { 6530 obj._thisSrf = null; 6531 obj._thatSrf = null; 6532 obj._thisOut = null; 6533 obj._thatOut = null; 6534 obj._appInteriorPnts = null; 6535 obj._tol = null; 6536 obj._patchSizeTol = null; 6537 obj._eps1 = null; 6538 obj._eps2 = null; 6539 obj._epsCRT = null; 6540 } 6541 }; 6542 6543 } 6544 6545}