diff --git a/Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/Curvature.java b/Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/Curvature.java index cba70946b8140075e26fef83d2a1d78cb834565e..6434b2db16b4f30cb7620d13ca5803071c79b1de 100644 --- a/Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/Curvature.java +++ b/Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/Curvature.java @@ -2,7 +2,6 @@ package cz.fidentis.analyst.visitors.mesh; import cz.fidentis.analyst.mesh.MeshVisitor; import cz.fidentis.analyst.mesh.core.MeshFacet; -import cz.fidentis.analyst.mesh.core.MeshPoint; import cz.fidentis.analyst.mesh.core.MeshTriangle; import cz.fidentis.analyst.mesh.core.TriangleFan; import java.util.ArrayList; @@ -10,6 +9,7 @@ import java.util.Collections; import java.util.HashMap; import java.util.List; import java.util.Map; +import java.util.Objects; import javax.vecmath.Vector3d; /** @@ -84,9 +84,8 @@ public class Curvature extends MeshVisitor { Cache cache = new Cache(facet); - //final List<CurvTriangle> triangles = precomputeTriangles(facet); - for (int vertA = 0; vertA < facet.getNumberOfVertices(); vertA++) { - TriangleFan oneRing = facet.getOneRingNeighborhood(vertA); + for (int vertIndexA = 0; vertIndexA < facet.getNumberOfVertices(); vertIndexA++) { + TriangleFan oneRing = facet.getOneRingNeighborhood(vertIndexA); if (oneRing.isBoundary()) { this.gaussian.get(facet).add(0.0); @@ -95,29 +94,29 @@ public class Curvature extends MeshVisitor { this.maxPrincipal.get(facet).add(0.0); continue; } - + double sumAngles = 0.0; double sumArea = 0.0; Vector3d pointSum = new Vector3d(); + Vector3d vertA = facet.getVertex(vertIndexA).getPosition(); for (int i = 0; i < oneRing.getVertices().size()-1; i++) { - int vertB = oneRing.getVertices().get(i); - int vertC = oneRing.getVertices().get(i+1); - int tri = oneRing.getTriangles().get(i); + Vector3d vertB = facet.getVertex(oneRing.getVertices().get(i)).getPosition(); + Vector3d vertC = facet.getVertex(oneRing.getVertices().get(i+1)).getPosition(); // Angles: - sumAngles += cache.getAngle(vertA, vertB, vertC, tri); + sumAngles += cache.getAngle(vertA, vertB, vertC); // Area: - double alpha = cache.getAngle(vertA, vertB, vertC, tri); - double beta = cache.getAngle(vertB, vertA, vertC, tri); - double gamma = cache.getAngle(vertC, vertA, vertB, tri); + double alpha = cache.getAngle(vertA, vertB, vertC); + double beta = cache.getAngle(vertB, vertA, vertC); + double gamma = cache.getAngle(vertC, vertA, vertB); - Vector3d ab = new Vector3d(facet.getVertex(vertB).getPosition()); - ab.sub(facet.getVertex(vertA).getPosition()); + Vector3d ab = new Vector3d(vertB); + ab.sub(vertA); - Vector3d ac = new Vector3d(facet.getVertex(vertC).getPosition()); - ac.sub(facet.getVertex(vertA).getPosition()); + Vector3d ac = new Vector3d(vertC); + ac.sub(vertA); double lenSqrtAB = ab.lengthSquared(); double lenSqrtAC = ac.lengthSquared(); @@ -127,23 +126,23 @@ public class Curvature extends MeshVisitor { double area = 0.5 * lenSqrtAB * lenSqrtAC * Math.sin(alpha); sumArea += (alpha > piHalf) ? area / 2.0 : area / 4.0; } else { - double cotGamma = 1.0 / Math.tan(gamma); - double cotBeta = 1.0 / Math.tan(beta); + double cotGamma = cache.getCotan(vertC, vertA, vertB); + double cotBeta = cache.getCotan(vertB, vertA, vertC); sumArea += (lenSqrtAB * cotGamma + lenSqrtAC * cotBeta) / 8.0; } // Laplace: - int triNext = oneRing.getTriangles().get((i + 1) % oneRing.getTriangles().size()); - int vertD; + Vector3d vertD; if (i == (oneRing.getVertices().size() - 2)) { // last triangle - vertD = oneRing.getVertices().get(1); + vertD = facet.getVertex(oneRing.getVertices().get(1)).getPosition(); } else { - vertD = oneRing.getVertices().get(i + 2); + vertD = facet.getVertex(oneRing.getVertices().get(i+2)).getPosition(); } - ac.scale(cache.getCtan(vertA, vertB, vertC, tri) + cache.getCtan(vertA, vertC, vertD, triNext)); + //ac.scale(cache2.get(tri).getCotan(vertB) + cache2.get(triNext).getCotan(vertD)); + ac.scale(cache.getCotan(vertB, vertA, vertC) + cache.getCotan(vertD, vertA, vertC)); pointSum.add(ac); } - + double gaussVal = (2.0 * Math.PI - sumAngles) / sumArea; double meanVal = 0.25 * sumArea * pointSum.length(); double delta = Math.max(0, Math.pow(meanVal, 2) - gaussVal); @@ -154,55 +153,35 @@ public class Curvature extends MeshVisitor { this.maxPrincipal.get(facet).add(meanVal + Math.sqrt(delta)); } } - - /* - protected List<CurvTriangle> precomputeTriangles(MeshFacet facet) { - List<CurvTriangle> ret = new ArrayList<>(facet.getNumTriangles()); - for (MeshTriangle tri: facet) { - ret.add(new CurvTriangle(facet, tri)); - } - return ret; - } - */ - - /** - * Helper class that caches triangle characteristics used multiples times during the curvature computation. - * <ul> - * <li>A = central point of the 1-ring neighborhood</li> - * <li>B = point "on the left" (previous point in the clockwise orientation of the triangle</li> - * <li>C = point "on the right" (next point in the clockwise orientation of the triangle</li> - * </ul> - * - * @author Radek Oslejsek - */ - protected class CurvTriangle { - - } /** * Cache key. * @author Radek Oslejsek */ private class CacheKey { - - private final int vertIndex; - private final int triIndex; + + private final Vector3d vCentral; + private final Vector3d v2; + private final Vector3d v3; /** * Constructor. - * @param v vertex - * @param i triangle + * @param vCentral The vertex for which the values are stored + * @param v2 Other vertex in the triangle + * @param v3 Other vertex in the triangle */ - CacheKey(int v, int i) { - this.vertIndex = v; - this.triIndex = i; + CacheKey(Vector3d vCentral, Vector3d v2, Vector3d v3) { + this.vCentral = vCentral; + this.v2 = v2; + this.v3 = v3; } @Override public int hashCode() { int hash = 7; - hash = 61 * hash + this.vertIndex; - hash = 61 * hash + this.triIndex; + hash = 11 * hash + Objects.hashCode(this.vCentral); + hash = 11 * hash + Objects.hashCode(this.v2); + hash = 11 * hash + Objects.hashCode(this.v3); return hash; } @@ -218,76 +197,127 @@ public class Curvature extends MeshVisitor { return false; } final CacheKey other = (CacheKey) obj; - if (this.vertIndex != other.vertIndex) { + if (!Objects.equals(this.vCentral, other.vCentral)) { return false; } - if (this.triIndex != other.triIndex) { + if (!Objects.equals(this.v2, other.v2)) { + return false; + } + if (!Objects.equals(this.v3, other.v3)) { return false; } return true; } + + @Override + public String toString() { + return vCentral + " | " + v2 + " | " + v3; + } + } /** - * Cache. + * Helper class that caches triangle characteristics used multiples times during the curvature computation. * @author Radek Oslejsek */ private class Cache { - private final MeshFacet facet; - private final Map<CacheKey, Double> cache = new HashMap<>(); - + private final Map<CacheKey, List<Double>> cache = new HashMap<>(); + + /** + * Constructor. + * @param facet Mesh facet + */ Cache(MeshFacet facet) { - this.facet = facet; + for (MeshTriangle tri: facet.getTriangles()) { + computeValues(facet, tri); + } + } + + /** + * Returns cached angle in the vertex {@code vCentral} at the triangle + * containing vertices {@code v2} and {@code v3}. + * + * @param vCentral Central vertex + * @param v2 Another vertex in the triangle + * @param v3 Another vertex in the triangle + * @return Angle at vertex {@code vCentral} + */ + public double getAngle(Vector3d vCentral, Vector3d v2, Vector3d v3) { + CacheKey key = new CacheKey(vCentral, v2, v3); + if (cache.containsKey(key)) { + return cache.get(key).get(0); + } + + key = new CacheKey(vCentral, v3, v2); + if (cache.containsKey(key)) { + return cache.get(key).get(0); + } + + throw new IllegalArgumentException("[" + vCentral + ", " + v2 + ", " + v3 + "]"); } - public double getAngle(int vertA, int vertB, int vertC, int tri) { - CacheKey key = new CacheKey(vertA, tri); + /** + * Returns cached cotangent of the angle in the vertex {@code vCentral} at the triangle + * containing vertices {@code v2} and {@code v3}. + * + * @param vCentral Central vertex + * @param v2 Another vertex in the triangle + * @param v3 Another vertex in the triangle + * @return Cotangent of the angle at vertex {@code vCentral} + */ + public double getCotan(Vector3d vCentral, Vector3d v2, Vector3d v3) { + CacheKey key = new CacheKey(vCentral, v2, v3); if (cache.containsKey(key)) { - return cache.get(key); + return cache.get(key).get(1); } - //System.out.println("XXX"+vertA+"|"+vertB+"|"+vertC); - //System.out.println("YYY"+facet.getVertex(vertA).getPosition()+"|"+facet.getVertex(vertB).getPosition()+"|"+facet.getVertex(vertC).getPosition()); + key = new CacheKey(vCentral, v3, v2); + if (cache.containsKey(key)) { + return cache.get(key).get(1); + } - Vector3d ab = new Vector3d(facet.getVertex(vertB).getPosition()); - ab.sub(facet.getVertex(vertA).getPosition()); - ab.normalize(); + throw new IllegalArgumentException("[" + vCentral + ", " + v2 + ", " + v3 + "]"); + } + + private void computeValues(MeshFacet facet, MeshTriangle tri) { + Vector3d ab = new Vector3d(facet.getVertex(tri.index2).getPosition()); + ab.sub(facet.getVertex(tri.index1).getPosition()); - Vector3d ac = new Vector3d(facet.getVertex(vertC).getPosition()); - ac.sub(facet.getVertex(vertA).getPosition()); - ac.normalize(); + Vector3d ac = new Vector3d(facet.getVertex(tri.index3).getPosition()); + ac.sub(facet.getVertex(tri.index1).getPosition()); - Vector3d bc = new Vector3d(facet.getVertex(vertC).getPosition()); - bc.sub(facet.getVertex(vertB).getPosition()); - bc.normalize(); + Vector3d bc = new Vector3d(facet.getVertex(tri.index3).getPosition()); + bc.sub(facet.getVertex(tri.index2).getPosition()); - //System.out.println("AAA"+ac+"|"+ab+"|"+bc); + ab.normalize(); + ac.normalize(); + bc.normalize(); - double cosAlpha = ab.dot(ac); + double cosA = ab.dot(ac); ab.scale(-1.0); - double cosBeta = ab.dot(bc); + double cosB = ab.dot(bc); - bc.scale(-1.0); ac.scale(-1.0); - double cosGamma = bc.dot(ac); + bc.scale(-1.0); + double cosC = ac.dot(bc); - //System.out.println("BBB"+cosAlpha+"|"+cosBeta+"|"+cosGamma); + List<Double> list1 = new ArrayList<>(); + List<Double> list2 = new ArrayList<>(); + List<Double> list3 = new ArrayList<>(); - double ret = Math.acos(cosAlpha); + list1.add(Math.acos(cosA)); + list2.add(Math.acos(cosB)); + list3.add(Math.acos(cosC)); - cache.put(new CacheKey(vertA, tri), ret); - cache.put(new CacheKey(vertB, tri), Math.acos(cosBeta)); - cache.put(new CacheKey(vertC, tri), Math.acos(cosGamma)); + list1.add(1.0 / Math.tan(list1.get(0))); + list2.add(1.0 / Math.tan(list2.get(0))); + list3.add(1.0 / Math.tan(list3.get(0))); - return ret; - } - - public double getCtan(int vertA, int vertB, int vertC, int tri) { - return 1.0 / Math.tan(getAngle(vertA, vertB, vertC, tri)); + cache.put(new CacheKey(tri.getVertex1(), tri.getVertex2(), tri.getVertex3()), list1); + cache.put(new CacheKey(tri.getVertex2(), tri.getVertex1(), tri.getVertex3()), list2); + cache.put(new CacheKey(tri.getVertex3(), tri.getVertex2(), tri.getVertex1()), list3); } - } - - + } } diff --git a/GUI/src/main/java/cz/fidentis/analyst/tests/HeatMapsTestApp.java b/GUI/src/main/java/cz/fidentis/analyst/tests/HeatMapsTestApp.java index 1a95ee80eaf3dc497581bce53d67240dbe60129c..c67674657e98a463d8b513d50b8de425a00a7b2a 100644 --- a/GUI/src/main/java/cz/fidentis/analyst/tests/HeatMapsTestApp.java +++ b/GUI/src/main/java/cz/fidentis/analyst/tests/HeatMapsTestApp.java @@ -127,15 +127,12 @@ public class HeatMapsTestApp extends GeneralGLEventListener { System.err.println(); System.err.println("MODEL: " + getModel()); - Set<MeshPoint> uniqueVerts = new HashSet<>(getModel().getFacets().get(0).getVertices()); - System.err.println("UNIQUE VERTICES: " + uniqueVerts.size()); - - /* - Set<Vector3d> unique2 = new HashSet<>(); + Set<Vector3d> unique = new HashSet<>(); for (MeshPoint p: getModel().getFacets().get(0).getVertices()) { - unique2.add(p.getPosition()); + unique.add(p.getPosition()); } - System.out.println("UNIQUE VERTICES: " + unique2.size()); + System.out.println("UNIQUE VERTICES: " + unique.size()); + /* for (Vector3d v: unique2) { for (MeshPoint p: uniqueVerts) { if (v.equals(p.getPosition())) { @@ -181,8 +178,8 @@ public class HeatMapsTestApp extends GeneralGLEventListener { long duration = System.currentTimeMillis() - startTime; System.err.println(duration + "\tmsec: Mean curvature"); - minDistance = 0.0; - maxDistance = 2.0; + minDistance = -0.001; + maxDistance = 0.2; } /** diff --git a/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/CornerTable.java b/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/CornerTable.java index 93fb252704d0bec963e057c1602da2b7e95bb780..045a8666dadb95061252a558543c3226d992c4ba 100644 --- a/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/CornerTable.java +++ b/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/CornerTable.java @@ -1,6 +1,7 @@ package cz.fidentis.analyst.mesh.core; import java.util.ArrayList; +import java.util.Collections; import java.util.HashMap; import java.util.List; import java.util.Map; @@ -215,6 +216,15 @@ public class CornerTable { .map(rowIndex -> rows.get(rowIndex)) .collect(Collectors.toList()); } + + /** + * Returns cache. + * + * @return cache + */ + public Map<Integer, List<Integer>> getVertToCornerCache() { + return Collections.unmodifiableMap(vertexToRow); + } /** * returns list of indexes of faces that have a corner at specific vertex diff --git a/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshFacetImpl.java b/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshFacetImpl.java index b2d53377a7cf32dfd337ae3f0c0086304d23dfc4..af454d2e3ef578962638d7070e4cc4037c1cb2d4 100644 --- a/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshFacetImpl.java +++ b/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshFacetImpl.java @@ -131,13 +131,13 @@ public class MeshFacetImpl implements MeshFacet { List<Integer> adjacentTrianglesI = cornerTable.getTriangleIndexesByVertexIndex(vertexIndex); for (Integer triI: adjacentTrianglesI) { List<Integer> triVerticesI = cornerTable.getIndexesOfVerticesByTriangleIndex(triI); - MeshTriangle tri = new MeshTriangle( + MeshTriangle tri = new MeshTriangle( this, triVerticesI.get(0), triVerticesI.get(1), triVerticesI.get(2)); - ret.add(tri); - } + ret.add(tri); + } return ret; } diff --git a/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/TriangleFan.java b/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/TriangleFan.java index cbf29e2b9ed14320134826d2a55c04d3300c42c8..eec8bd25d17a278533061b3bf14731335732bdd5 100644 --- a/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/TriangleFan.java +++ b/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/TriangleFan.java @@ -35,6 +35,10 @@ public class TriangleFan implements Iterable<MeshTriangle> { centralVertex = vert; this.facet = facet; + // get first appearance of the vertex in the corner table: + int vertRow = facet.getCornerTable().getVertToCornerCache().get(vert).get(0); + + /* int vertRow = -1; for (int i = 0; i < facet.getCornerTable().getSize(); i++) { if (facet.getCornerTable().getRow(i).getVertexIndex() == vert) { @@ -45,6 +49,7 @@ public class TriangleFan implements Iterable<MeshTriangle> { if (vertRow == -1) { throw new IllegalArgumentException("vert"); } + */ computeOneRingData(vertRow); } @@ -101,6 +106,9 @@ public class TriangleFan implements Iterable<MeshTriangle> { break; // we reached an open end } triangles.add(ct.getIndexOfFace(ringVertRow)); + if (vertices.size() > ct.getSize()) { //emergency break + throw new RuntimeException("Error in mesh topology"); + } } ringVertRow = ct.getIndexOfPreviousCornerInFace(vertRow); @@ -114,6 +122,9 @@ public class TriangleFan implements Iterable<MeshTriangle> { break; // we reached an open end } triangles.add(ct.getIndexOfFace(ringVertRow)); + if (vertices.size() > ct.getSize()) { //emergency break + throw new RuntimeException("Error in mesh topology"); + } } }