001/* 002 * MeshUtil -- Static utility methods used to work with triangular meshes. 003 * 004 * Copyright (C) 2018, by Joseph A. Huwaldt 005 * All rights reserved. 006 * 007 * This library is free software; you can redistribute it and/or 008 * modify it under the terms of the GNU Lesser General Public 009 * License as published by the Free Software Foundation; either 010 * version 2.1 of the License, or (at your option) any later version. 011 * 012 * This library is distributed in the hope that it will be useful, 013 * but WITHOUT ANY WARRANTY; without even the implied warranty of 014 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 015 * Lesser General Public License for more details. 016 * 017 * You should have received a copy of the GNU Lesser General Public 018 * License along with this library; if not, write to the Free Software 019 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, 020 * MA 02110-1301 USA 021 */ 022package geomss.geom; 023 024import geomss.GeomSSScene; 025import geomss.app.MainWindow; 026import static geomss.geom.GeomUtil.SQRT_EPS; 027import geomss.ui.Colors; 028import jahuwaldt.js.param.Parameter; 029import jahuwaldt.swing.MDIApplication; 030import jahuwaldt.tools.math.MathTools; 031import jahuwaldt.util.PointRegionQuadTree; 032import java.util.*; 033import java.util.logging.Level; 034import java.util.logging.Logger; 035import javax.measure.quantity.Area; 036import javax.measure.quantity.Dimensionless; 037import javax.measure.quantity.Length; 038import javax.measure.unit.Unit; 039import javolution.context.StackContext; 040 041/** 042 * A collection of static methods used to work with triangular unstructured meshes. 043 * 044 * <p> Modified by: Joseph A. Huwaldt </p> 045 * 046 * @author Joseph A. Huwaldt, Date: February 4, 2018 047 * @version April 11, 2018 048 */ 049public final class MeshUtil { 050 051 private static final boolean DEBUG = false; 052 053 private static final double SQRT3 = Math.sqrt(3); 054 private static final double HALF_SQRT3 = 0.5*SQRT3; // Height of equilateral triangle with side length of 1. 055 private static final double IDEAL_AR = 3./SQRT3; // Aspect ratio of an equilateral triangle. 056 057 058 /** 059 * Generate an unstructured mesh/grid on a planar region defined by the input boundary 060 * points. This version of this method uses a default minimum triangle area and does 061 * not apply any smoothing. 062 * 063 * @param boundaryPnts The ordered boundary points for the region being meshed. The 064 * physical dimension of the points must be either 2 or 3. The 065 * points must be planar and must form a closed loop with the 066 * "inside" to the left of the unit vectors between each pair of 067 * points (when looking down the unit vector). 068 * @param tol The geometric tolerance on the input points. 069 * @return A list of triangles that represent the final unstructured mesh. 070 * @see #advancingFrontPlanar(geomss.geom.PointString, jahuwaldt.js.param.Parameter, jahuwaldt.js.param.Parameter, jahuwaldt.js.param.Parameter) 071 */ 072 public static TriangleList<Triangle> advancingFrontPlanar(PointString<? extends GeomPoint> boundaryPnts, Parameter<Length> tol) { 073 Parameter<Area> minArea = tol.times(tol); 074 return advancingFrontPlanar(boundaryPnts, tol, minArea, null); 075 } 076 077 /** 078 * Generate an unstructured mesh/grid on a planar region defined by the input boundary 079 * points. 080 * 081 * @param boundaryPnts The ordered boundary points for the region being meshed. The 082 * physical dimension of the points must be either 2 or 3. The 083 * points must be planar and must form a closed loop with the 084 * "inside" to the left of the unit vectors between each pair of 085 * points (when looking down the unit vector). 086 * @param tol The geometric tolerance on the input points. 087 * @param minArea The minimum area triangle to allow in the mesh. 088 * @param smoothingEps The tolerance for stopping the smoothing process. Smoothing 089 * stops when the maximum node position change in the mesh is less 090 * than this value. If this value is null or zero, then no 091 * smoothing is done. 092 * @return A list of triangles that represent the final unstructured mesh. 093 */ 094 @SuppressWarnings({"null", "SleepWhileInLoop"}) 095 public static TriangleList<Triangle> advancingFrontPlanar(PointString<? extends GeomPoint> boundaryPnts, 096 Parameter<Length> tol, Parameter<Area> minArea, Parameter<Length> smoothingEps) { 097 098 // Reference: The VGRID algorithm. 099 // Parikh,P., Pirzadeh, S., "Package for 3D Unstructured Grid Generation", NASA-CR-182090, Sept. 1990. 100 101 TriangleList<Triangle> output = TriangleList.newInstance(); 102 int numPnts = boundaryPnts.size(); 103 if (numPnts < 3) 104 return output; 105 106 int ndim = boundaryPnts.getPhyDimension(); 107 if (ndim > 3 || ndim < 2) 108 throw new IllegalArgumentException("Front physical dimension must be either 2 or 3."); 109 110 if (numPnts == 3) { 111 // Degenerate case with only 3 points input (a single triangle output). 112 Triangle tri = Triangle.valueOf(boundaryPnts); 113 output.add(tri); 114 return output; 115 } 116 117 // The string of points must form a closed loop. 118 if (!boundaryPnts.get(0).isApproxEqual(boundaryPnts.get(-1), tol)) 119 throw new IllegalArgumentException("Front points must form a closed loop with the 1st and last points equal in position."); 120 121 if (numPnts == 4) { 122 // Degenerate case with only 4 points input (1st point == last point; a single triangle output). 123 Triangle tri = Triangle.valueOf(boundaryPnts.getRange(0,2)); 124 output.add(tri); 125 return output; 126 } 127 128 Triangle[] tri_arr = null; 129 StackContext.enter(); 130 try { 131 TriangleList<Triangle> mesh = TriangleList.newInstance(); 132 133 GTransform T_ZN = GTransform.IDENTITY; 134 Parameter<Length> Zpos = null; 135 if (ndim == 3) { 136 // Figure out if the front is planar. 137 if (!boundaryPnts.isPlanar(tol)) 138 throw new IllegalArgumentException("The boundary front must be planar to within tol."); 139 140 // Compute the normal vector for the point string. 141 Vector<Dimensionless> nhat = getPointStringNormal(boundaryPnts); 142 143 if (!nhat.get(Point.Z).isApproxEqual(Parameter.ONE)) { 144 // Not already in the XY plane, so rotate it. 145 146 // Create a unit vector in the coord. system Z axis direction. 147 Vector<Dimensionless> zhat = Vector.valueOf(Dimensionless.UNIT, 0, 0, 1); 148 149 // Determine the rotation from nhat to zhat (Z wrt N). 150 T_ZN = GTransform.valueOf(nhat, zhat); 151 152 // Transform the front from an arbitrary plane to parallelt to the XY plane. 153 boundaryPnts = boundaryPnts.getTransformed(T_ZN).copyToReal(); 154 } 155 156 // Get the Z-position of the planar point string. 157 Zpos = boundaryPnts.get(0).get(Point.Z); 158 159 // Create a transformation for translating the Z position to 0. 160 Parameter<Length> ZERO = Parameter.ZERO_LENGTH; 161 T_ZN = GTransform.newTranslation(ZERO, ZERO, Zpos.opposite()).times(T_ZN); 162 163 // Convert the front to 2D (drop the Z component). 164 boundaryPnts = boundaryPnts.toDimension(2); 165 } 166 167 // Set up consistent units. 168 Unit<Length> unit = boundaryPnts.getUnit(); 169 tol = tol.to(unit); 170 minArea = minArea.to((Unit<Area>)unit.pow(2)); 171 172 // Convert front points to be immutable and in the same units. 173 // Last point is redundant since it must be equal to the 1st. 174 PointString<Point> boundaryNodes = PointString.newInstance(); 175 numPnts = numPnts - 1; 176 for (int i = 0; i < numPnts; ++i) 177 boundaryNodes.add(boundaryPnts.get(i).immutable().to(unit)); 178 179 // Get the bounds of the boundary front points. 180 Point pMin = boundaryNodes.getBoundsMin(); 181 Vector<Length> diag = Vector.valueOf(boundaryNodes.getBoundsMax().minus(pMin)); 182 183 // Convert the list of input front points into a Heap List of line segments (faces) 184 // where the length of the segments determines the order in the list. 185 // Also link all the points to the faces they are connected to. 186 PriorityQueue<LineSeg> front = createFrontHeapList(boundaryNodes); 187 188 // Determine the maximum starting front face length. 189 double Lmax = 0; 190 for (LineSeg face : front) { 191 double L = face.getArcLength(0).getValue(); 192 if (L > Lmax) 193 Lmax = L; 194 } 195 196 // Place all the existing points in a quad-tree. 197 double tolv = tol.to(unit).getValue(); 198 PointRegionQuadTree<Point> tree = new PointRegionQuadTree( 199 pMin.minus(Point.valueOf(tol, tol)), diag.getValue(Point.X) + 2 * tolv, diag.getValue(Point.Y) + 2 * tolv); 200 for (Point p : boundaryNodes) { 201 boolean r = tree.insert(p); 202 if (!r) 203 System.out.println("Error, couldn't add point " + p); 204 } 205 206 // Create a list to store all of the iterior nodes in the mesh. 207 GeomList<Point> interiorNodes = GeomList.newInstance(); 208 209 GeomSSScene scene = null; 210 if (DEBUG) { 211 // Get debugging view. 212 scene = ((MainWindow)MDIApplication.getInstance().getTopWindow()).getScene(); 213 scene.draw(boundaryNodes.getAll(), true); 214 scene.centerAndZoom(); 215 scene.setPointSize(5); 216 scene.setLineWidth(5); 217 scene.setSurfaceAlpha(0.5F); 218 GeomList<LineSeg> frontSegs = GeomList.newInstance(); 219 scene.setLineColor(Colors.BLUE); 220 frontSegs.addAll(front); 221 scene.draw(frontSegs); 222 scene.setLineColor(Colors.GREEN); 223 } 224 225 // Loop over all the front faces until they have all been removed. 226 double Rv = 0; 227 while (!front.isEmpty()) { 228 // Find the next face to be removed from the front (the shortest one) 229 // and remove it from the front. 230 LineSeg face = front.poll(); 231 232 // Determine a "best point" for the next triangle to be formed. 233 // Form stretched equilateral triangle with the face as the short side. 234 Point p0 = face.getStart(); 235 Point p1 = face.getEnd(); 236 Point pmid = p0.plus(p1).divide(2); 237 Vector<Dimensionless> face_n = GeomUtil.normal2D(face.getTangent(0)); 238 double Lface = p0.distanceValue(p1); 239 double stretching = 1.5; 240 double h = stretching * Lface * HALF_SQRT3; // Height of equilateral triangle: h = 0.5*sqrt(3)*L 241 Point pbest = pmid.plus(Point.valueOf(face_n.times(Parameter.valueOf(h, unit)))); 242 243 if (DEBUG) { 244 System.out.println("Triangle # " + (mesh.size() + 1)); 245 scene.draw(p0); 246 scene.draw(p1); 247 scene.draw(pbest); 248 try { 249 Thread.sleep(1); 250 } catch (InterruptedException ex) { 251 Logger.getLogger(MeshUtil.class.getName()).log(Level.SEVERE, null, ex); 252 } 253 } 254 255 // Find the adjacent face(s) to the face being removed. 256 GeomList<LineSeg> adj_faces = getAdjacentFaces(face, face_n); 257 258 // Find all the close points to the new point. 259 Rv = Math.max(Rv, 1.5 * Math.max(Lmax, h)); 260 double xmin = pbest.getValue(Point.X) - Rv; 261 double ymin = pbest.getValue(Point.Y) - Rv; 262 List<Point> closePoints = tree.queryRange(xmin, ymin, 2 * Rv, 2 * Rv); 263 264 // Filter out the end points of the current face. 265 closePoints.remove(p0); 266 closePoints.remove(p1); 267 268 // Find all the faces associated with all the closest points. 269 List<LineSeg> closeFaces = GeomList.newInstance(); 270 for (Point p : closePoints) 271 addUniqueFaces(closeFaces, (List)p.getUserData("faces_AF")); 272 273 // Filter the close points to remove any that can not be the "best" point. 274 filterClosePoints(closePoints, p0, face_n, pbest, adj_faces); 275 276 int size = closePoints.size(); 277 if (size > 1) 278 // Order the candidate points by apex angle (largest to smallest). 279 closePoints.sort(new ApexAngleComparator(face)); 280 281 // Starting from the top of the list of close points, check to see if the 282 // new element formed by the point would cross any existing faces. If it does 283 // not cross, then take that point, otherwise move down the list. 284 Parameter<Length> minD = face.getArcLength(0).times(0.5); 285 Point pbest_orig = pbest; 286 pbest = null; // Will cause exception to be thrown if there are no valid points. 287 for (int i = 0; i < size; ++i) { 288 Point p = closePoints.get(i); 289 290 boolean hasIntersect = false; 291 for (LineSeg cface : closeFaces) { 292 if (lineSegsIntersect(cface, p, p0) || lineSegsIntersect(cface, p, p1)) { 293 hasIntersect = true; 294 break; 295 } 296 297 // Try to avoid a newly selected point falling too close to an existing face. 298 if (size > 1 && p == pbest_orig) { 299 if (isPointNearLineSeg(cface, p, minD)) { 300 hasIntersect = true; 301 break; 302 } 303 } 304 } 305 if (!hasIntersect) { 306 if (DEBUG) 307 scene.erase(pbest_orig); 308 pbest = p; 309 break; 310 } 311 } 312 if (pbest == null) 313 throw new NullPointerException("pbest is null -- problem with algorithm"); 314 315 if (pbest == pbest_orig) 316 // A new point is being used. 317 interiorNodes.add(pbest); 318 319 // Is the best point associated with any faces already. 320 GeomList<LineSeg> faces = (GeomList<LineSeg>)pbest.getUserData("faces_AF"); 321 if (faces == null) { 322 // Got our initial guess point. 323 324 // Create two new faces. 325 LineSeg f1 = LineSeg.valueOf(p0, pbest); 326 LineSeg f2 = LineSeg.valueOf(pbest, p1); 327 front.add(f1); 328 front.add(f2); 329 330 // Update the linked-list of points-to-faces. 331 GeomList<LineSeg>p0Faces = (GeomList)p0.getUserData("faces_AF"); 332 p0Faces.remove(face); 333 p0Faces.add(f1); 334 GeomList<LineSeg>p1Faces = (GeomList)p1.getUserData("faces_AF"); 335 p1Faces.remove(face); 336 p1Faces.add(f2); 337 GeomList<LineSeg> pbFaces = GeomList.newInstance(); 338 pbFaces.add(f1); 339 pbFaces.add(f2); 340 pbest.putUserData("faces_AF", pbFaces); 341 pbest.putUserData("active_AF", Boolean.TRUE); 342 343 // Add the new point to the tree. 344 boolean r = tree.insert(pbest); 345 if (!r) 346 System.out.println("Error, couldn't add new point " + pbest); 347 348 if (DEBUG) { 349 scene.setLineColor(Colors.BLUE); 350 scene.draw(f1); 351 scene.draw(f2); 352 scene.setLineColor(Colors.GREEN); 353 } 354 355 } else { 356 // Got a point that is already on the front. 357 358 GeomList<LineSeg> p0Faces = (GeomList)p0.getUserData("faces_AF"); 359 p0Faces.remove(face); 360 GeomList<LineSeg>p1Faces = (GeomList)p1.getUserData("faces_AF"); 361 p1Faces.remove(face); 362 GeomList<LineSeg> pbFaces = (GeomList<LineSeg>)pbest.getUserData("faces_AF"); 363 LineSeg f1 = commonFace(pbFaces, p0Faces); 364 if (f1 == null) { 365 // Create a new face. 366 f1 = LineSeg.valueOf(p0, pbest); 367 368 // Add the new face to the front. 369 front.add(f1); 370 p0Faces.add(f1); 371 pbFaces.add(f1); 372 373 if (DEBUG) { 374 scene.setLineColor(Colors.BLUE); 375 scene.draw(f1); 376 scene.setLineColor(Colors.GREEN); 377 } 378 379 } else { 380 // Face already exists and will be removed from front by placing this triangle. 381 front.remove(f1); 382 p0Faces.remove(f1); 383 pbFaces.remove(f1); 384 385 if (DEBUG) 386 scene.erase(f1); 387 } 388 389 LineSeg f2 = commonFace(pbFaces, p1Faces); 390 if (f2 == null) { 391 // Create a new face. 392 f2 = LineSeg.valueOf(pbest, p1); 393 394 // Add the new face to the front. 395 front.add(f2); 396 p1Faces.add(f2); 397 pbFaces.add(f2); 398 399 if (DEBUG) { 400 scene.setLineColor(Colors.BLUE); 401 scene.draw(f2); 402 scene.setLineColor(Colors.GREEN); 403 } 404 405 } else { 406 // Face already exists and will be removed from front by placing this triangle. 407 front.remove(f2); 408 p1Faces.remove(f2); 409 pbFaces.remove(f2); 410 411 if (DEBUG) 412 scene.erase(f2); 413 } 414 415 // If a point no longer has any faces in it's faces list, then it is 416 // no longer on the active front. 417 if (p0Faces.isEmpty()) 418 p0.putUserData("active_AF", Boolean.FALSE); 419 if (p1Faces.isEmpty()) 420 p1.putUserData("active_AF", Boolean.FALSE); 421 } 422 423 // Construct a new triangle. 424 Triangle tri = Triangle.valueOf(p0, p1, pbest); 425 mesh.add(tri); 426 427 // Assign this new triangle to each of it's verticies. 428 assignTriangle2Node(p0, tri); 429 assignTriangle2Node(p1, tri); 430 assignTriangle2Node(pbest, tri); 431 432 if (DEBUG) { 433 scene.erase(face); 434 scene.draw(tri); 435 scene.erase(p0); 436 scene.erase(p1); 437 } 438 439 } // end while(!front.isEmpty()) 440 441 // Determine some geometric characteristics for each triangle. 442 calcTriangleChar(mesh); 443 444 // Size before cleaning. 445 int size = mesh.size(); 446 447 // Clean the mesh by removing any triangles that are too small in area. 448 removeSmallAreas(mesh, boundaryNodes, interiorNodes, minArea); 449 450 // Clean the mesh by removing short triangle edegs. 451 removeShortEdges(mesh, boundaryNodes, interiorNodes, 0.35); 452 453 // Smooth the mesh if requested. 454 smoothUnstructuredMesh2D(mesh, interiorNodes, smoothingEps); 455 456 // Remove all the user data from the triangle nodes. 457 for (Point p : boundaryNodes) { 458 p.removeUserData("active_AF"); 459 p.removeUserData("faces_AF"); 460 p.removeUserData("tris_AF"); 461 } 462 for (Point p : interiorNodes) { 463 p.removeUserData("active_AF"); 464 p.removeUserData("faces_AF"); 465 p.removeUserData("tris_AF"); 466 } 467 468 if (DEBUG) { 469 // Sort the input list of triangles by the shortest side to longest side ratio 470 // (largest to smallest). 471 Comparator<Triangle> comparator = new Comparator<Triangle>() { 472 @Override 473 public int compare(Triangle o1, Triangle o2) { 474 Double AR1 = (Double)o1.getUserData("short_long_AF"); 475 Double AR2 = (Double)o2.getUserData("short_long_AF"); 476 return -AR1.compareTo(AR2); 477 } 478 }; 479 mesh.sort(comparator); 480 System.out.println("Triangles removed by cleaning = " + (size - mesh.size())); 481 } 482 483 // Convert the geometry back to 3D and original orientation if necessary. 484 if (Zpos != null) 485 mesh = mesh.toDimension(3).getTransformed(T_ZN.transpose()).copyToReal(); 486 487 // Copy the mesh triangles to the outer context. 488 size = mesh.size(); 489 tri_arr = new Triangle[size]; 490 for (int i=0; i < size; ++i) { 491 Triangle t = mesh.get(i); 492 tri_arr[i] = StackContext.outerCopy(t); 493 } 494 495 if (DEBUG) { 496 scene.setPointSize(2); 497 scene.setLineWidth(1); 498 scene.setSurfaceAlpha(1); 499 } 500 501 } finally { 502 StackContext.exit(); 503 } 504 505 // Add the array of triangles to the outut list. 506 output.addAll(tri_arr); 507 508 return output; 509 } 510 511 /** 512 * Filter out any points from the input closePoints list that are not on the active 513 * front, on the wrong side of the current face (p0, and face_n) or that are on the 514 * wrong side of any adjacent faces to the current face. 515 * 516 * @param closePoints The list of close points to be filtered. 517 * @param p0 The starting point for the current face. 518 * @param face_n The normal vector for the current face. 519 * @param pbest A guess at the "best" next point. 520 * @param adj_faces The list of adjacent faces to the current face. 521 */ 522 private static void filterClosePoints(List<Point> closePoints, Point p0, 523 Vector<Dimensionless> face_n, Point pbest, GeomList<LineSeg> adj_faces) { 524 // Filter out any close points that are not on the front or that are on the 525 // wrong side of the face. 526 int size = closePoints.size(); 527 for (int i = size - 1; i >= 0; --i) { 528 Point p = closePoints.get(i); 529 530 // Is the point inactive? 531 if (!(Boolean)p.getUserData("active_AF")) 532 closePoints.remove(i); 533 534 else { 535 // Is the point on the wrong side of the face. 536 if (wrongSide(p0, face_n, p)) 537 closePoints.remove(i); 538 } 539 } 540 541 // Add the guess at the best point to the list of candidate points. 542 closePoints.add(pbest); 543 544 // Filter any points that are on the wrong side of any of the adjacent faces. 545 size = closePoints.size(); 546 for (int i = size - 1; i >= 0; --i) { 547 Point p = closePoints.get(i); 548 for (LineSeg adj : adj_faces) { 549 Vector<Dimensionless> adj_n = GeomUtil.normal2D(adj.getTangent(0)); 550 if (wrongSide(adj.getStart(), adj_n, p)) { 551 closePoints.remove(i); 552 break; 553 } 554 } 555 } 556 } 557 558 /** 559 * Create an initial front of line segments from the ordered boundary nodes. 560 * 561 * @param boundaryNodes The initial ordered boundary nodes. 562 * @return A priority queue/heap list of the faces of the initial front. 563 */ 564 private static PriorityQueue<LineSeg> createFrontHeapList(PointString<Point> boundaryNodes) { 565 // A PriorityQueue is Java's "Heap List" structure. 566 int numPnts = boundaryNodes.size(); 567 PriorityQueue<LineSeg> front = new PriorityQueue(numPnts, new ShortestLengthComparator()); 568 GeomList<LineSeg> p0Faces = GeomList.newInstance(); 569 Point p0 = boundaryNodes.get(0); 570 p0.putUserData("faces_AF", p0Faces); 571 p0.putUserData("active_AF", Boolean.TRUE); // Point is actively a part of the front. 572 for (int i = 1; i < numPnts; ++i) { 573 GeomList<LineSeg> p1Faces = GeomList.newInstance(); 574 Point p1 = boundaryNodes.get(i); 575 p1.putUserData("faces_AF", p1Faces); 576 p1.putUserData("active_AF", Boolean.TRUE); 577 578 LineSeg face = LineSeg.valueOf(p0, p1); 579 front.add(face); 580 p0Faces.add(face); 581 p1Faces.add(face); 582 p0 = p1; 583 p0Faces = p1Faces; 584 } 585 Point p1 = boundaryNodes.get(0); // Because last point must equal 1st point (closed list of points). 586 GeomList<LineSeg> p1Faces = (GeomList<LineSeg>)p1.getUserData("faces_AF"); 587 LineSeg face = LineSeg.valueOf(p0, p1); 588 front.add(face); 589 p0Faces.add(face); 590 p1Faces.add(0, face); 591 return front; 592 } 593 594 /** 595 * Smooth the unstructured mesh by shifting each interior node to the center of the 596 * surrounding polygon. 597 * 598 * @param tris The list of all the triangles in the mesh. 599 * @param nodes The list of all the interior nodes/points in the mesh. 600 * @param eps The tolerance for stopping the smoothing process. Smoothing stops when 601 * the maximum node position change in the mesh is less than this value. 602 * If this value is null or zero, then no smoothing is done. 603 */ 604 private static void smoothUnstructuredMesh2D(TriangleList<Triangle> tris, 605 GeomList<Point> nodes, Parameter<Length> eps) { 606 if (eps == null || eps.isApproxZero()) 607 return; 608 609 // Reference: Naji, H.S., "A C# Algorithm for Triangular Meshes of Highly-Irregular 610 // 2D Domains using Advancing Front Technique", JKAU, V17N2, pp 41-72, 2006. 611 612 // Loop over the smoothing process until the maximum node movement is less than eps. 613 double eps2 = eps.getValue(tris.getUnit()); 614 eps2 = eps2*eps2; 615 double ri2max = Double.MAX_VALUE; 616 int iterations = 0; 617 while (iterations < 100 && ri2max > eps2) { 618 ++iterations; 619 ri2max = 0; 620 int size = nodes.size(); 621 for (int i=size-1; i >= 0; --i) { 622 Point pi = nodes.get(i); 623 624 // Get the triangles surrounding node pi 625 GeomList<Triangle> tris_pi = (GeomList)pi.getUserData("tris_AF"); 626 627 // Get the points surrounding pi (the surrounding polygon). 628 PointString<Point> etai = PointString.newInstance(); 629 for (Triangle t : tris_pi) { 630 GeomList<Point> tpnts = (GeomList)t.getAll(); 631 addUniquePoints(etai, tpnts); 632 GeomList.recycle(tpnts); 633 } 634 etai.remove(pi); 635 636 // The new location for pi is average of surrounding points. 637 Point pi2 = etai.getAverage(); 638 639 // Distance squared between old and new interior node position. 640 double ri2 = pi2.distanceSqValue(pi); 641 ri2max = Math.max(ri2max, ri2); 642 643 // If the distance is greater than threshold, then replace the original node. 644 if (ri2 >= eps2) { 645 // Modify all the triangles associated with pi to replace pi with pi2. 646 replaceVertex(tris, pi, pi2); 647 nodes.remove(pi); 648 nodes.add(pi2); 649 } 650 } 651 652 if (DEBUG) { 653 System.out.println("smoothing iterations = " + iterations + ", rmax = " + Math.sqrt(ri2max)); 654 } 655 } 656 657 } 658 659 /** 660 * Replaces a mesh node/vertex with a new point while adjusting all the triangles that 661 * share that point. 662 * 663 * @param tris The list of all the triangles in the mesh. 664 * @param p1 The point that is being replaced. 665 * @param p2 The point that is replacing p1. 666 */ 667 private static void replaceVertex(TriangleList<Triangle> tris, Point p1, Point p2) { 668 GeomList<Triangle> tris_p1 = (GeomList)p1.getUserData("tris_AF"); 669 GeomList<Triangle> tris_p2 = GeomList.newInstance(); 670 671 // Loop over all the triangles associated with the point being replaced. 672 for (Triangle t : tris_p1) { 673 // Remove this triangle from the triangles list (it is being replaced). 674 tris.remove(t); 675 676 // Replace p1 with p2 in the list of vertices from this triangle. 677 GeomList<Point> verts = (GeomList)t.getAll(); 678 int idx = verts.indexOf(p1); 679 verts.remove(idx); 680 verts.add(idx, p2); 681 682 // Create a new triangle with p1 replaced by p2. 683 Triangle t2 = Triangle.valueOf(verts); 684 calcTriangleChar(t2); 685 686 // Save off the new triangle in the triangle list. 687 tris.add(t2); 688 tris_p2.add(t2); 689 690 // Add the new triangle to the triangle lists for the other 2 points. 691 for (int i=idx+1; i < 3; ++i) { 692 Point pi = verts.get(i); 693 GeomList<Triangle> tris_pi = (GeomList)pi.getUserData("tris_AF"); 694 tris_pi.remove(t); 695 tris_pi.add(t2); 696 } 697 for (int i=0; i < idx; ++i) { 698 Point pi = verts.get(i); 699 GeomList<Triangle> tris_pi = (GeomList)pi.getUserData("tris_AF"); 700 tris_pi.remove(t); 701 tris_pi.add(t2); 702 } 703 //GeomList.recycle(verts); 704 } 705 706 // Associate the list of triangles that use p2 with node p2. 707 p2.putUserData("tris_AF", tris_p2); 708 } 709 710 /** 711 * Remove any non-boundary triangles where the shortest-to-longest side length ratio 712 * is below a threshold and adjust neighboring triangles to fill in the gap. 713 * 714 * @param tris The list of all the triangles with pre-computed 715 * characteristics "short_long_AF" and "shortest_side_AF". 716 * @param boundaryNodes The node points that formed the original boundary of the mesh. 717 * @param interiorNodes The node points that are interior to the mesh. 718 * @param threshold The shortest-to-longest side length ratio threshold for 719 * removing triangles. 720 */ 721 private static void removeShortEdges(TriangleList<Triangle> tris, 722 PointString<Point> boundaryNodes, GeomList<Point> interiorNodes, double threshold) { 723 724 // Loop over all the triangles. If any have a shortest edge to longest edge 725 // ratio < 0.35, then collapse the shortest edge to it's start (which removes two 726 // triangles and adjust any triangles that connect to the end of that edge 727 // to instead point to the start. 728 int size = tris.size(); 729 for (int i = size - 1; i >= 0; --i) { 730 Triangle t = tris.get(i); 731 double lratio = (double)t.getUserData("short_long_AF"); 732 if (lratio < threshold) { 733 // Find the shortest edge of this triangle. 734 int sedge = (int)t.getUserData("shortest_side_AF"); 735 736 // Do not remove an edge of the original boundary. 737 if (!isBoundaryEdge(boundaryNodes, t, sedge)) { 738 // Remove the triangle by collapsing it's shortest edge. 739 int n = collapseTriangleEdge(tris, interiorNodes, t, sedge); 740 i -= n; 741 } 742 } 743 } 744 745 } 746 747 /** 748 * Remove any non-boundary triangles that are not on the boundary where the area is 749 * below a threshold and adjust neighboring triangles to fill in the gap. 750 * 751 * @param tris The list of all the triangles with pre-computed 752 * characteristics "short_long_AF" and "shortest_side_AF". 753 * @param boundaryNodes The node points that formed the original boundary of the mesh. 754 * @param interiorNodes The node points that are interior to the mesh. 755 * @param threshold The triangle area threshold for removing triangles. 756 */ 757 private static void removeSmallAreas(TriangleList<Triangle> tris, 758 PointString<Point> boundaryNodes, GeomList<Point> interiorNodes, Parameter<Area> threshold) { 759 760 if (DEBUG) { 761 Comparator<Triangle> comparator = new Comparator<Triangle>() { 762 @Override 763 public int compare(Triangle o1, Triangle o2) { 764 return -o1.getArea().compareTo(o2.getArea()); 765 } 766 }; 767 tris.sort(comparator); 768 } 769 770 // Loop over all the triangles. If any have a shortest edge to longest edge 771 // ratio < 0.35, then collapse the shortest edge to it's start (which removes two 772 // triangles and adjust any triangles that connect to the end of that edge 773 // to instead point to the start. 774 int size = tris.size(); 775 for (int i = size - 1; i >= 0; --i) { 776 Triangle t = tris.get(i); 777 Parameter<Area> area = t.getArea(); 778 if (area.isLessThan(threshold)) { 779 // Find the shortest edge of this triangle. 780 int sedge = (int)t.getUserData("shortest_side_AF"); 781 782 // Do not remove an edge of the original boundary. 783 if (!isBoundaryEdge(boundaryNodes, t, sedge)) { 784 // Remove the triangle by collapsing it's shortest edge. 785 int n = collapseTriangleEdge(tris, interiorNodes, t, sedge); 786 i -= n; 787 } 788 } 789 } 790 791 } 792 793 /** 794 * Return true if the specified edge of the specified triangle falls on the original 795 * boundary of the mesh. 796 * 797 * @param boundary The points that make up the original boundary of the mesh. 798 * @param tri The triangle being checked to see if it contains a boundary edge. 799 * @param edge The edge to check to see if it is on the outer boundary of the 800 * mesh. 801 * @return <code>true</code> if the specified edge of the triangle is on the outside 802 * boundary of the mesh. 803 */ 804 private static boolean isBoundaryEdge(PointString<Point> boundary, Triangle tri, int edge) { 805 806 // Get the start and end points on the specified edge. 807 Point p0, p1; 808 switch (edge) { 809 case 0: 810 p0 = tri.getP1(); 811 p1 = tri.getP2(); 812 break; 813 case 1: 814 p0 = tri.getP2(); 815 p1 = tri.getP3(); 816 break; 817 default: 818 p0 = tri.getP3(); 819 p1 = tri.getP1(); 820 break; 821 } 822 823 return boundary.contains(p0) && boundary.contains(p1); 824 } 825 826 /** 827 * Remove the specified triangle by collapsing the specified edge while adjusting all 828 * the adjacent triangles to fill in the gap. This will result in the deletion of the 829 * triangle that shares the specified edge and the adjustment of any triangles that 830 * connect to the end of that edge to instead connect to the start of that edge. 831 * 832 * @param tris The list of all the triangles in the mesh. 833 * @param interiorNodes The node points that are interior to the mesh. 834 * @param t1 The triangle being deleted. 835 * @param edge The edge to be collapsed in order to delete the triangle. 836 * @return The number of triangles deleted from the mesh. Will be either 0 (if t1 837 * could not be deleted) or 2. 838 */ 839 private static int collapseTriangleEdge(TriangleList<Triangle> tris, GeomList<Point> interiorNodes, 840 Triangle t1, int edge) { 841 842 // Get the start and end points on the specified edge. 843 Point p0, p1; 844 switch (edge) { 845 case 0: 846 p0 = t1.getP1(); 847 p1 = t1.getP2(); 848 break; 849 case 1: 850 p0 = t1.getP2(); 851 p1 = t1.getP3(); 852 break; 853 default: 854 p0 = t1.getP3(); 855 p1 = t1.getP1(); 856 break; 857 } 858 859 // Get all the triangles assocated with these points. 860 GeomList<Triangle> tris_p0 = (GeomList)p0.getUserData("tris_AF"); 861 GeomList<Triangle> tris_p1 = (GeomList)p1.getUserData("tris_AF"); 862 863 // Find the "other" triangle in common between these points besides "t1". 864 Triangle t2 = null; 865 outer: 866 for (Triangle t_0 : tris_p0) { 867 if (t_0 == t1) 868 continue; 869 for (Triangle t_1 : tris_p1) { 870 if (t_1 == t1) 871 continue; 872 if (t_1 == t_0) { 873 t2 = t_1; 874 break outer; 875 } 876 } 877 } 878 if (t2 == null) { 879 System.out.println("Unexpected t2 == null in MeshUtil.removeShortEdges()."); 880 return 0; 881 } 882 883 // Delete t1 and t2 from the triangle lists. 884 tris_p0.remove(t1); 885 tris_p0.remove(t2); 886 tris_p1.remove(t1); 887 tris_p1.remove(t2); 888 tris.remove(t1); 889 tris.remove(t2); 890 891 // Delete the point p1 from the list of interior nodes. 892 interiorNodes.remove(p1); 893 894 // Replace all the other triangles that touch p1 with versions that are 895 // adjusted to touch p0 instead. 896 for (Triangle t : tris_p1) { 897 // Remove this triangle from the triangles list (it is being replaced). 898 tris.remove(t); 899 900 // Replace p1 with p0 in the list of vertices from this triangle. 901 GeomList<Point> verts = (GeomList)t.getAll(); 902 int idx = verts.indexOf(p1); 903 verts.remove(idx); 904 verts.add(idx, p0); 905 906 // Create a new triangle with p1 replaced by p0. 907 Triangle t_2 = Triangle.valueOf(verts); 908 calcTriangleChar(t_2); 909 910 // Save off the new triangle in the triangle list. 911 tris.add(t_2); 912 tris_p0.add(t_2); 913 914 // Add the new triangle to the triangle lists for the other 2 points. 915 for (int i=idx+1; i < 3; ++i) { 916 Point pi = verts.get(i); 917 GeomList<Triangle> tris_pi = (GeomList)pi.getUserData("tris_AF"); 918 tris_pi.remove(t1); 919 tris_pi.remove(t2); 920 tris_pi.remove(t); 921 tris_pi.add(t_2); 922 } 923 for (int i=0; i < idx; ++i) { 924 Point pi = verts.get(i); 925 GeomList<Triangle> tris_pi = (GeomList)pi.getUserData("tris_AF"); 926 tris_pi.remove(t1); 927 tris_pi.remove(t2); 928 tris_pi.remove(t); 929 tris_pi.add(t_2); 930 } 931 GeomList.recycle(verts); 932 } 933 934 return 2; 935 } 936 937 /** 938 * Calculate the characteristics for a list of triangles. The computed values are 939 * stored in each triangle's user data with the following keys: 940 * <pre> 941 * aspect_ratio_AF = (double) Aspect ratio scaled relative to that of an equilateral triangle. 942 * short_long_AF = (double) Ratio of the shortest to longest sides of the triangle. 943 * shortest_side_AF = (int) Index of the shortest side of the triangle. 944 * </pre> 945 * 946 * @param tris The list of triangles to have the characteristics calculated for. 947 */ 948 private static void calcTriangleChar(TriangleList<Triangle> tris) { 949 for (Triangle t : tris) { 950 calcTriangleChar(t); 951 } 952 } 953 954 /** 955 * Calculate some geometric characteristics of the input triangle. The computed values 956 * are stored in the triangle's user data with the following keys: 957 * <pre> 958 * aspect_ratio_AF = (double) Aspect ratio scaled relative to that of an equilateral triangle. 959 * short_long_AF = (double) Ratio of the shortest to longest sides of the triangle. 960 * shortest_side_AF = (int) Index of the shortest side of the triangle. 961 * </pre> 962 * 963 * @param tri The triangle to have the characteristics calculated for. 964 */ 965 private static void calcTriangleChar(Triangle tri) { 966 // The diameter of the incircle of the triangle. 967 // Ref: https://en.wikipedia.org/wiki/Incircle_and_excircles_of_a_triangle#Relation_to_area_of_the_triangle 968 Point p1 = tri.getP1(); 969 Point p2 = tri.getP2(); 970 Point p3 = tri.getP3(); 971 double a = p1.distanceValue(p2); 972 double b = p2.distanceValue(p3); 973 double c = p3.distanceValue(p1); 974 double s = 0.5 * (a + b + c); // Semi-permieter of triangle. 975 double d = 2*tri.getArea().getValue()/s; // In-circle diameter. 976 977 // The triangle "scaled aspect ratio". 978 double lmax = Math.max(Math.max(a, b), c); 979 double AR = lmax/d/IDEAL_AR; // Scaled relative to an equilateral triangle. 980 tri.putUserData("aspect_ratio_AF", AR); 981 982 // The shortest side of the triangle. 983 s = Math.min(Math.min(a, b), c); 984 int shortest = 2; 985 if (s == a) 986 shortest = 0; 987 else if (s == b) 988 shortest = 1; 989 tri.putUserData("shortest_side_AF", shortest); 990 tri.putUserData("short_long_AF", s / lmax); 991 } 992 993 /** 994 * Assign the input triangle to the list of associated triangles stored in the user 995 * data for the input point. If the input point doesn't have a list of triangles, it 996 * is created and the input triangle is added to it. 997 * 998 * @param p The point to have the triangle added to the list of associated 999 * triangles. 1000 * @param tri The triangle associated with the point. 1001 */ 1002 private static void assignTriangle2Node(Point p, Triangle tri) { 1003 GeomList<Triangle> tris = (GeomList)p.getUserData("tris_AF"); 1004 if (tris == null) { 1005 tris = GeomList.newInstance(); 1006 p.putUserData("tris_AF", tris); 1007 } 1008 tris.add(tri); 1009 } 1010 1011 /** 1012 * Return a list of adjacent faces that are oriented correctly and on the "inside" of 1013 * the target face. 1014 * 1015 * @param face The face to return the adjacent faces for. 1016 * @param face_n The pre-calculated face normal for face. 1017 * @return A list of adjacent faces. May contain 1 or two faces. 1018 */ 1019 private static GeomList<LineSeg> getAdjacentFaces(LineSeg face, geomss.geom.Vector<Dimensionless> face_n) { 1020 Point p0 = face.getStart(); 1021 Point p1 = face.getEnd(); 1022 1023 GeomList<LineSeg> adj_faces = GeomList.newInstance(); 1024 adj_faces.addAll((GeomList<LineSeg>)p0.getUserData("faces_AF")); 1025 adj_faces.remove(face); 1026 adj_faces.addAll((GeomList<LineSeg>)p1.getUserData("faces_AF")); 1027 adj_faces.remove(face); 1028 1029 // Remove any adjacent faces that are oriented the wrong way or are on the wrong 1030 // side of the input face. 1031 int size = adj_faces.size(); 1032 for (int i = size - 1; i >= 0; --i) { 1033 LineSeg adj = adj_faces.get(i); 1034 Point ps = adj.getStart(); 1035 Point pe = adj.getEnd(); 1036 if (ps.isApproxEqual(p0)) 1037 // Adjacent-face starts at the start of the input face. 1038 // This adjacent-face is oriented the wrong way, remove it. 1039 adj_faces.remove(i); 1040 else if (ps.isApproxEqual(p1)) { 1041 // Adjacent-face starts at the end of the face. 1042 // The other end of the adjacent-face must be on the inward side of the input face. 1043 if (wrongSide(p0, face_n, pe)) 1044 adj_faces.remove(i); 1045 } else if (pe.isApproxEqual(p0)) { 1046 // Adjacent-face ends at the start of the face. 1047 // The other end of the adjacent-face must be on the inward side of the input face. 1048 if (wrongSide(p0, face_n, ps)) 1049 adj_faces.remove(i); 1050 } else { // else if (pe.isApproxEqual(p1)) { 1051 // Adjacent-face ends at the end of the input face. 1052 // This adjacent-face is oriented the wrong way, remove it. 1053 adj_faces.remove(i); 1054 } 1055 } 1056 1057 return adj_faces; 1058 } 1059 1060 /** 1061 * Find and return the 1st edge/face in common between the two lists of faces. 1062 * 1063 * @param faces1 The 1st list of faces. 1064 * @param faces2 The 2nd list of faces. 1065 * @return The first face found in common between the two input lists. 1066 */ 1067 private static LineSeg commonFace(GeomList<LineSeg> faces1, GeomList<LineSeg> faces2) { 1068 LineSeg common = null; 1069 for (LineSeg f : faces1) { 1070 if (faces2.contains(f)) { 1071 common = f; 1072 break; 1073 } 1074 } 1075 return common; 1076 } 1077 1078 /** 1079 * Return true if the input point is on the "wrong side" of the input face. 1080 * If the point falls on the input line segment, it is not on the wrong side. 1081 * 1082 * @param p0 Any point on the face being tested against. 1083 * @param face_n The normal vector for the face the point is being tested against. 1084 * @param p The point being tested. 1085 * @return True if the point is on the wrong side of the face. 1086 */ 1087 private static boolean wrongSide(Point p0, Vector<Dimensionless> face_n, Point p) { 1088 StackContext.enter(); 1089 try { 1090 Vector<Length> v = p.minus(p0).toGeomVector(); 1091 Parameter dotp = v.dot(face_n); 1092 return dotp.getValue() < -MathTools.EPS; 1093 } finally { 1094 StackContext.exit(); 1095 } 1096 } 1097 1098 /** 1099 * Return true if the two specified line segments intersect. 1100 * 1101 * @param seg1 The 1st line segment being tested for intersection. 1102 * @param s2p0 The start point of the 2nd line segment being tested for intersection. 1103 * @param s2p1 The end point of the 2nd line segment being tested for intersection. 1104 * @return true if the two line segments intersect. Touching at the ends is not 1105 * considered an intersection. 1106 */ 1107 private static boolean lineSegsIntersect(LineSeg seg1, Point s2p0, Point s2p1) { 1108 // Ref: https://martin-thoma.com/how-to-check-if-two-line-segments-intersect/ 1109 StackContext.enter(); 1110 try { 1111 LineSeg seg2 = LineSeg.valueOf(s2p0, s2p1); 1112 1113 // Check for intersection. 1114 return GeomUtil.boundsIntersect(seg1, seg2) && 1115 lineSegmentCrossesLine2D(seg1, seg2) && 1116 lineSegmentCrossesLine2D(seg2, seg1); 1117 1118 } finally { 1119 StackContext.exit(); 1120 } 1121 } 1122 1123 /** 1124 * Check if line segment a crosses the line that is defined by line segment b. If the 1125 * ends touch, that does not count as "crossing". 1126 * 1127 * @param a first line segment interpreted as line 1128 * @param b second line segment interpreted as line 1129 * @return <code>true</code> if line segment a crosses line b in the interior, 1130 * <code>false</code> otherwise. 1131 */ 1132 private static boolean lineSegmentCrossesLine2D(LineSeg a, LineSeg b) { 1133 return !isPointOnLine(a, b.getStart()) && !isPointOnLine(a,b.getEnd()) && 1134 isPointRightOfLine2D(a, b.getStart()) ^ isPointRightOfLine2D(a, b.getEnd()); 1135 } 1136 1137 /** 1138 * Checks if a point is to the right of a line. If the point is on the line, it is not 1139 * right of the line. 1140 * 1141 * @param a line segment interpreted as a line 1142 * @param b the point 1143 * @return <code>true</code> if the point is right of the line, <code>false</code> 1144 * otherwise 1145 */ 1146 private static boolean isPointRightOfLine2D(LineSeg a, Point b) { 1147 // Move the line segment and point, so that a.getStart() is on (0,0). 1148 Point ap0 = a.getStart(); 1149 Point ap1 = a.getEnd(); 1150 double x0 = ap0.getValue(Point.X); 1151 double y0 = ap0.getValue(Point.Y); 1152 double ax = ap1.getValue(Point.X) - x0; 1153 double ay = ap1.getValue(Point.Y) - y0; 1154 double bx = b.getValue(Point.X) - x0; 1155 double by = b.getValue(Point.Y) - y0; 1156 double r = cross2D(ax,ay, bx,by); 1157 boolean result = r < -MathTools.EPS; 1158 return result; 1159 } 1160 1161 /** 1162 * Checks if a Point is on a line. 1163 * 1164 * @param a line (interpreted as line, although given as line segment). 1165 * @param b point tested to see if it is on the line. 1166 * @return <code>true</code> if point is on line, otherwise <code>false</code> 1167 */ 1168 private static boolean isPointOnLine(LineSeg a, Point b) { 1169 // Move the line segment and point, so that a.getStart() is on (0,0). 1170 Point ap0 = a.getStart(); 1171 Point ap1 = a.getEnd(); 1172 double x0 = ap0.getValue(Point.X); 1173 double y0 = ap0.getValue(Point.Y); 1174 double ax = ap1.getValue(Point.X) - x0; 1175 double ay = ap1.getValue(Point.Y) - y0; 1176 double bx = b.getValue(Point.X) - x0; 1177 double by = b.getValue(Point.Y) - y0; 1178 double r = cross2D(ax,ay, bx,by); 1179 boolean result = MathTools.isApproxZero(r); 1180 return result; 1181 } 1182 1183 /** 1184 * Return the cross product between two 2D vectors. 1185 * 1186 * @param ax The X-component of the 1st vector. 1187 * @param ay The Y-component of the 1st vector. 1188 * @param bx The X-component of the 2nd vector. 1189 * @param by The Y-component of the 2nd vector. 1190 * @return The cross product of the two 2D vectors. 1191 */ 1192 private static double cross2D(double ax, double ay, double bx, double by) { 1193 return ax * by - bx * ay; 1194 } 1195 1196 /** 1197 * Return true if a point is near a line segment to within the specified tolerance. 1198 * 1199 * @param seg The line segment to see if a point is near. 1200 * @param p The point to test for being near the line segment. 1201 * @param tol The tolerance for being near the line segment. 1202 * @return <code>true</code> if the point is closer to the line segment than tol. 1203 */ 1204 private static boolean isPointNearLineSeg(LineSeg seg, Point p, Parameter<Length> tol) { 1205 StackContext.enter(); 1206 try { 1207 1208 Parameter<Length> d = GeomUtil.pointLineSegDistance(p, seg.getStart(), seg.getEnd()); 1209 return d.isLessThan(tol); 1210 1211 } finally { 1212 StackContext.exit(); 1213 } 1214 } 1215 1216 /** 1217 * Add any faces in list2 to list1 that are not already in list1. 1218 * 1219 * @param list1 The list to add the faces to. 1220 * @param list2 The list to pull any unique faces from. 1221 */ 1222 private static void addUniqueFaces(List<LineSeg> list1, List<LineSeg> list2) { 1223 for (LineSeg f : list2) { 1224 if (!list1.contains(f)) 1225 list1.add(f); 1226 } 1227 } 1228 1229 /** 1230 * Add any points in list2 to list1 that are not already in list1. 1231 * 1232 * @param list1 The list to add the points to. 1233 * @param list2 The list to pull any unique points from. 1234 */ 1235 private static void addUniquePoints(List<Point> list1, List<Point> list2) { 1236 for (Point f : list2) { 1237 if (!list1.contains(f)) 1238 list1.add(f); 1239 } 1240 } 1241 1242 /** 1243 * A Comparator that compares the length of a pair of line segments and identifies 1244 * the one with the shortest length. 1245 */ 1246 private static class ShortestLengthComparator implements Comparator<LineSeg> { 1247 1248 @Override 1249 public int compare(LineSeg o1, LineSeg o2) { 1250 Parameter len1 = o1.getArcLength(0); 1251 Parameter len2 = o2.getArcLength(0); 1252 return len1.compareTo(len2); 1253 } 1254 1255 } 1256 1257 /** 1258 * A Comparator that compares the apex angles of the triangle formed by the defined 1259 * face segment and a pair of points and identifies the one with the largest apex angle. 1260 */ 1261 private static class ApexAngleComparator implements Comparator<Point> { 1262 private final LineSeg face; 1263 1264 public ApexAngleComparator(LineSeg face) { 1265 this.face = face; 1266 } 1267 1268 @Override 1269 public int compare(Point o1, Point o2) { 1270 Point p0 = face.getStart(); 1271 Point p1 = face.getEnd(); 1272 1273 double cosphi1 = cosC(p0, p1, o1); 1274 double cosphi2 = cosC(p0, p1, o2); 1275 1276 // If cos(phi1) < cos(phi2) then phi1 > phi2. 1277 return Double.compare(cosphi1, cosphi2); 1278 } 1279 } 1280 1281 /** 1282 * Returns the cosine of the interior angle in corner "C" formed by a triangle 1283 * consisting of points pA, pB and pC. 1284 * 1285 * @param pA The first point of a triangle. 1286 * @param pB The second point of a triangle. 1287 * @param pC The third point of a triangle where the interior angle is to be 1288 * calculated. 1289 * @return The cosine of the angle in the corner of a triangle at point pC. 1290 */ 1291 private static double cosC(Point pA, Point pB, Point pC) { 1292 double a = pB.distanceValue(pC); 1293 double b = pA.distanceValue(pC); 1294 double c2 = pB.distanceSqValue(pA); 1295 double cosC = (a*a + b*b - c2)/(2*a*b); 1296 return cosC; 1297 } 1298 1299 /** 1300 * Returns the normal vector for a planar PointString. No check is made here for 1301 * the point string being planar. 1302 * 1303 * @param pnts The point string to return the normal vector for. 1304 * @return The normal vector for the planar PointString. 1305 */ 1306 private static Vector<Dimensionless> getPointStringNormal(PointString<? extends GeomPoint> pnts) { 1307 1308 StackContext.enter(); 1309 try { 1310 // Extract 3 points from the string of points and form a plane form them. 1311 GeomPoint p0 = pnts.get(0); 1312 int numPnts = pnts.size(); 1313 int idx1 = numPnts / 3; 1314 if (idx1 == 0) 1315 ++idx1; 1316 int idx2 = 2 * numPnts / 3; 1317 if (idx2 == idx1) 1318 ++idx2; 1319 GeomPoint p1 = pnts.get(idx1); 1320 GeomPoint p2 = pnts.get(idx2); 1321 Vector<Length> v1 = p1.minus(p0).toGeomVector(); 1322 Vector<Length> v2 = p2.minus(p0).toGeomVector(); 1323 Vector n = v1.cross(v2); 1324 Vector<Dimensionless> nhat; 1325 if (n.magValue() > SQRT_EPS) { 1326 nhat = n.toUnitVector(); 1327 nhat.setOrigin(Point.newInstance(p0.getPhyDimension())); 1328 1329 } else { 1330 // Try a different set of points on the string. 1331 idx1 = numPnts / 4; 1332 if (idx1 == 0) 1333 ++idx1; 1334 idx2 = 3 * numPnts / 4; 1335 if (idx2 == idx1) 1336 ++idx2; 1337 p1 = pnts.get(idx1); 1338 p2 = pnts.get(idx2); 1339 v1 = p1.minus(p0).toGeomVector(); 1340 v2 = p2.minus(p0).toGeomVector(); 1341 n = v1.cross(v2); 1342 if (n.magValue() > SQRT_EPS) { 1343 nhat = n.toUnitVector(); 1344 nhat.setOrigin(Point.newInstance(p0.getPhyDimension())); 1345 } else { 1346 // Points from input curve likely form a straight line. 1347 nhat = GeomUtil.calcYHat(p1.minus(p2).toGeomVector().toUnitVector()); 1348 nhat = v2.cross(nhat).toUnitVector(); 1349 } 1350 } 1351 1352 return StackContext.outerCopy(nhat); 1353 1354 } finally { 1355 StackContext.exit(); 1356 } 1357 } 1358 1359}