From 8228dedef99708903646cc31c415399d21f335d6 Mon Sep 17 00:00:00 2001
From: Radek Oslejsek <oslejsek@fi.muni.cz>
Date: Tue, 6 Apr 2021 20:44:52 +0200
Subject: [PATCH] Fixed curvature calculations
---
.../analyst/visitors/mesh/Curvature.java | 202 ++++++++++++++----
.../analyst/tests/EfficiencyTests.java | 36 ++--
.../fidentis/analyst/mesh/core/MeshFacet.java | 13 +-
.../analyst/mesh/core/MeshFacetImpl.java | 22 +-
.../analyst/mesh/core/TriangleFan.java | 69 +++---
5 files changed, 236 insertions(+), 106 deletions(-)
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 d83e846f..b42c398b 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
@@ -4,6 +4,7 @@ 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;
import java.util.Collections;
import java.util.HashMap;
@@ -60,7 +61,7 @@ public class Curvature extends MeshVisitor {
}
/**
- * Returns maximu principal curvatures for all inspected mesh facets. The order corresponds to
+ * Returns maximum principal curvatures for all inspected mesh facets. The order corresponds to
* the order of vertices, i.e., i-th value represents the curvature of i-th mesh vertex.
*
* @return Maximum principal curvatures of inspected mesh facets.
@@ -81,12 +82,13 @@ public class Curvature extends MeshVisitor {
maxPrincipal.put(facet, new ArrayList<>());
}
- //final List<MeshTriangle> triangles = facet.getTriangles();
+ Cache cache = new Cache(facet);
+
final List<CurvTriangle> triangles = precomputeTriangles(facet);
for (int vertA = 0; vertA < facet.getNumberOfVertices(); vertA++) {
- List<Integer> neighbouringTriangles = facet.getCornerTable().getTriangleIndexesByVertexIndex(vertA);
+ TriangleFan oneRing = facet.getOneRingNeighborhood(vertA);
- if (facet.vertexIsBoundary(vertA)) {
+ if (oneRing.isBoundary()) {
this.gaussian.get(facet).add(0.0);
this.mean.get(facet).add(0.0);
this.minPrincipal.get(facet).add(0.0);
@@ -94,30 +96,52 @@ public class Curvature extends MeshVisitor {
continue;
}
- if (neighbouringTriangles.isEmpty()) {
- this.gaussian.get(facet).add(Double.NaN);
- this.mean.get(facet).add(Double.NaN);
- this.minPrincipal.get(facet).add(Double.NaN);
- this.maxPrincipal.get(facet).add(Double.NaN);
- continue;
- }
-
- double sumArea = 0.0;
double sumAngles = 0.0;
+ double sumArea = 0.0;
Vector3d pointSum = new Vector3d();
-
- // for all surrounding triangles:
- for (int i = 0; i < neighbouringTriangles.size(); i++) {
- CurvTriangle ctri = triangles.get(neighbouringTriangles.get(i));
- CurvTriangle tNext = triangles.get(neighbouringTriangles.get((i + 1) % neighbouringTriangles.size()));
-
- sumArea += computeArea(ctri, facet.getVertex(vertA));
- sumAngles += ctri.alpha(facet.getVertex(vertA));
- Vector3d aux = new Vector3d(ctri.vertC(facet.getVertex(vertA)));
- aux.sub(ctri.vertA(facet.getVertex(vertA)));
- aux.scale(ctri.betaCotan(facet.getVertex(vertA)) + tNext.gammaCotan(facet.getVertex(vertA)));
- pointSum.add(aux);
+ 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);
+
+ // Angles:
+ sumAngles += cache.getAngle(vertA, vertB, vertC, tri);
+
+ // 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);
+
+ Vector3d ab = new Vector3d(facet.getVertex(vertB).getPosition());
+ ab.sub(facet.getVertex(vertA).getPosition());
+
+ Vector3d ac = new Vector3d(facet.getVertex(vertC).getPosition());
+ ac.sub(facet.getVertex(vertA).getPosition());
+
+ double lenSqrtAB = ab.lengthSquared();
+ double lenSqrtAC = ac.lengthSquared();
+
+ double piHalf = Math.PI / 2.0; // 90 degrees
+ if (alpha >= piHalf || beta >= piHalf || gamma >= piHalf) { // check for obtuse angle
+ 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);
+ sumArea += (lenSqrtAB * cotGamma + lenSqrtAC * cotBeta) / 8.0;
+ }
+
+ // Laplace:
+ int triNext = oneRing.getTriangles().get((i + 1) % oneRing.getTriangles().size());
+ int vertD;
+ if (i == (oneRing.getVertices().size() - 2)) { // last triangle
+ vertD = oneRing.getVertices().get(1);
+ } else {
+ vertD = oneRing.getVertices().get(i + 2);
+ }
+ ac.scale(cache.getCtan(vertA, vertB, vertC, tri) + cache.getCtan(vertA, vertC, vertD, triNext));
+ pointSum.add(ac);
}
double gaussVal = (2.0 * Math.PI - sumAngles) / sumArea;
@@ -139,23 +163,6 @@ public class Curvature extends MeshVisitor {
return ret;
}
- protected double computeArea(CurvTriangle ctri, MeshPoint vertA) {
- double alpha = ctri.alpha(vertA);
- double beta = ctri.beta(vertA);
- double gamma = ctri.gamma(vertA);
-
- double piHalf = Math.PI / 2.0; // 90 degrees
- if (alpha >= piHalf || beta >= piHalf || gamma >= piHalf) { // check for obtuse angle
- return (alpha > piHalf) ? ctri.area() / 2.0 : ctri.area() / 4.0;
- } else {
- double cotBeta = ctri.betaCotan(vertA);
- double cotGamma = ctri.gammaCotan(vertA);
- double ab = ctri.lengthSquaredAB(vertA);
- double ac = ctri.lengthSquaredAC(vertA);
- return (ab * cotGamma + ac * cotBeta) / 8.0;
- }
- }
-
/**
* Helper class that caches triangle characteristics used multiples times during the curvature computation.
* <ul>
@@ -437,4 +444,115 @@ public class Curvature extends MeshVisitor {
}
}
+ /**
+ * Cache key.
+ * @author Radek Oslejsek
+ */
+ private class CacheKey {
+
+ private final int vertIndex;
+ private final int triIndex;
+
+ /**
+ * Constructor.
+ * @param v vertex
+ * @param i triangle
+ */
+ CacheKey(int v, int i) {
+ this.vertIndex = v;
+ this.triIndex = i;
+ }
+
+ @Override
+ public int hashCode() {
+ int hash = 7;
+ hash = 61 * hash + this.vertIndex;
+ hash = 61 * hash + this.triIndex;
+ return hash;
+ }
+
+ @Override
+ public boolean equals(Object obj) {
+ if (this == obj) {
+ return true;
+ }
+ if (obj == null) {
+ return false;
+ }
+ if (getClass() != obj.getClass()) {
+ return false;
+ }
+ final CacheKey other = (CacheKey) obj;
+ if (this.vertIndex != other.vertIndex) {
+ return false;
+ }
+ if (this.triIndex != other.triIndex) {
+ return false;
+ }
+ return true;
+ }
+ }
+
+ /**
+ * Cache.
+ * @author Radek Oslejsek
+ */
+ private class Cache {
+
+ private final MeshFacet facet;
+ private final Map<CacheKey, Double> cache = new HashMap<>();
+
+ Cache(MeshFacet facet) {
+ this.facet = facet;
+ }
+
+ public double getAngle(int vertA, int vertB, int vertC, int tri) {
+ CacheKey key = new CacheKey(vertA, tri);
+ if (cache.containsKey(key)) {
+ return cache.get(key);
+ }
+
+ //System.out.println("XXX"+vertA+"|"+vertB+"|"+vertC);
+ //System.out.println("YYY"+facet.getVertex(vertA).getPosition()+"|"+facet.getVertex(vertB).getPosition()+"|"+facet.getVertex(vertC).getPosition());
+
+ Vector3d ab = new Vector3d(facet.getVertex(vertB).getPosition());
+ ab.sub(facet.getVertex(vertA).getPosition());
+ ab.normalize();
+
+ Vector3d ac = new Vector3d(facet.getVertex(vertC).getPosition());
+ ac.sub(facet.getVertex(vertA).getPosition());
+ ac.normalize();
+
+ Vector3d bc = new Vector3d(facet.getVertex(vertC).getPosition());
+ bc.sub(facet.getVertex(vertB).getPosition());
+ bc.normalize();
+
+ //System.out.println("AAA"+ac+"|"+ab+"|"+bc);
+
+ double cosAlpha = ab.dot(ac);
+
+ ab.scale(-1.0);
+ double cosBeta = ab.dot(bc);
+
+ bc.scale(-1.0);
+ ac.scale(-1.0);
+ double cosGamma = bc.dot(ac);
+
+ //System.out.println("BBB"+cosAlpha+"|"+cosBeta+"|"+cosGamma);
+
+ double ret = Math.acos(cosAlpha);
+
+ 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));
+
+ return ret;
+ }
+
+ public double getCtan(int vertA, int vertB, int vertC, int tri) {
+ return 1.0 / Math.tan(getAngle(vertA, vertB, vertC, tri));
+ }
+ }
+
+
}
diff --git a/GUI/src/main/java/cz/fidentis/analyst/tests/EfficiencyTests.java b/GUI/src/main/java/cz/fidentis/analyst/tests/EfficiencyTests.java
index 64a50054..aa40e41d 100644
--- a/GUI/src/main/java/cz/fidentis/analyst/tests/EfficiencyTests.java
+++ b/GUI/src/main/java/cz/fidentis/analyst/tests/EfficiencyTests.java
@@ -40,32 +40,32 @@ public class EfficiencyTests {
face2 = new HumanFace(boyFile);
boolean relativeDist = false;
- boolean printDetails = true;
+ boolean printDetails = false;
//face1.getMeshModel().compute(new GaussianCurvature()); // initialize everything, then measure
- //System.out.println("Symmetry plane calculation:");
- //printSymmetryPlane(face1, true, SignificantPoints.CurvatureAlg.MEAN);
- //printSymmetryPlane(face1, true, SignificantPoints.CurvatureAlg.GAUSSIAN);
- //printSymmetryPlane(face1, false, SignificantPoints.CurvatureAlg.MEAN);
- //printSymmetryPlane(face1, false, SignificantPoints.CurvatureAlg.GAUSSIAN);
+ System.out.println("Symmetry plane calculation:");
+ printSymmetryPlane(face1, true, SignificantPoints.CurvatureAlg.MEAN);
+ printSymmetryPlane(face1, true, SignificantPoints.CurvatureAlg.GAUSSIAN);
+ printSymmetryPlane(face1, false, SignificantPoints.CurvatureAlg.MEAN);
+ printSymmetryPlane(face1, false, SignificantPoints.CurvatureAlg.GAUSSIAN);
- //System.out.println();
- //System.out.println(measureKdTreeCreation(face1) + "\tmsec:\tKd-tree creation of first face");
- //System.out.println(measureKdTreeCreation(face2) + "\tmsec:\tKd-tree creation of second face");
+ System.out.println();
+ System.out.println(measureKdTreeCreation(face1) + "\tmsec:\tKd-tree creation of first face");
+ System.out.println(measureKdTreeCreation(face2) + "\tmsec:\tKd-tree creation of second face");
System.out.println();
System.out.println("Hausdorff distance sequentially:");
- testAndPrint(face2, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_POINT, relativeDist, false), printDetails);
- //testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_POINT_DISTANCE_ONLY, false, false), printDetails);
- //testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_POINT, relativeDist, false), printDetails);
- //testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_TRIANGLE_APPROXIMATE, relativeDist, false), printDetails);
+ testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_POINT, relativeDist, false), printDetails);
+ testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_POINT_DISTANCE_ONLY, false, false), printDetails);
+ testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_POINT, relativeDist, false), printDetails);
+ testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_TRIANGLE_APPROXIMATE, relativeDist, false), printDetails);
- //System.out.println();
- //System.out.println("Hausdorff distance concurrently:");
- //testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_POINT_DISTANCE_ONLY, false, true), printDetails);
- //testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_POINT, relativeDist, true), printDetails);
- //testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_TRIANGLE_APPROXIMATE, relativeDist, true), printDetails);
+ System.out.println();
+ System.out.println("Hausdorff distance concurrently:");
+ testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_POINT_DISTANCE_ONLY, false, true), printDetails);
+ testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_POINT, relativeDist, true), printDetails);
+ testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_TRIANGLE_APPROXIMATE, relativeDist, true), printDetails);
}
protected static void testAndPrint(HumanFace face, HausdorffDistance vis, boolean printDetails) {
diff --git a/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshFacet.java b/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshFacet.java
index c9dcefae..27cfd798 100644
--- a/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshFacet.java
+++ b/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshFacet.java
@@ -5,7 +5,7 @@ import cz.fidentis.analyst.mesh.MeshVisitor;
import javax.vecmath.Vector3d;
/**
- * An ancapsulated mesh plate, i.e., multiple triangles sharing vertices.
+ * An encapsulated mesh plate, i.e., multiple triangles sharing vertices.
* Mesh facet is iterable (triangle as returned) and visitable by
* {@link cz.fidentis.analyst.mesh.MeshVisitor}s.
*
@@ -114,8 +114,8 @@ public interface MeshFacet extends Iterable<MeshTriangle> {
void accept(MeshVisitor visitor);
/**
- * Computes centers of circumcircle of all triaqngles.
- * These points represent the point of Voronoi area used for Delaunay triangulation, for instenace.
+ * Computes centers of circumcircle of all triangles.
+ * These points represent the point of Voronoi area used for Delaunay triangulation, for instance.
* The list is computed only once (during the first call) and than cached.
* The order corresponds to the order of triangles, i.e., the i-th point
* is the Voronoi point of i-th triangle.
@@ -125,9 +125,10 @@ public interface MeshFacet extends Iterable<MeshTriangle> {
List<Vector3d> calculateVoronoiPoints();
/**
- * Returns {@code true} if the given vertex lie at the facet boundary.
+ * Returns 1-ring neighborhood, i.e., triangles around the given mesh point.
+ *
* @param vertexIndex Index of mesh vertex
- * @return {@code true} if the given vertex lie at the facet boundary.
+ * @return Triangles around the vertex or {@code null}
*/
- boolean vertexIsBoundary(int vertexIndex);
+ TriangleFan getOneRingNeighborhood(int vertexIndex);
}
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 b009ee71..6511eb24 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
@@ -7,10 +7,8 @@ import java.util.List;
import java.util.Map;
import javax.vecmath.Vector3d;
import cz.fidentis.analyst.mesh.MeshVisitor;
-import java.util.HashSet;
import java.util.Iterator;
import java.util.NoSuchElementException;
-import java.util.Set;
/**
* Mash facet is a compact triangular mesh without duplicated vertices.
@@ -146,18 +144,6 @@ public class MeshFacetImpl implements MeshFacet {
return ret;
}
- @Override
- public boolean vertexIsBoundary(int vertexIndex) {
- List<MeshTriangle> tris = getAdjacentTriangles(vertexIndex);
- Set<Vector3d> vset = new HashSet<>();
- for (MeshTriangle tri: tris) {
- for (MeshPoint p: tri) {
- vset.add(p.getPosition());
- }
- }
- return (tris.size()+1) != vset.size();
- }
-
@Override
public Vector3d getClosestAdjacentPoint(Vector3d point, int vertexIndex) {
double dist = Double.POSITIVE_INFINITY;
@@ -236,5 +222,13 @@ public class MeshFacetImpl implements MeshFacet {
}
return Collections.unmodifiableList(voronoiPoints);
}
+
+ @Override
+ public TriangleFan getOneRingNeighborhood(int vertexIndex) {
+ if (vertexIndex < 0 || vertexIndex >= this.getNumberOfVertices()) {
+ return null;
+ }
+ return new TriangleFan(this, vertexIndex);
+ }
}
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 c8df70d9..9536cbbb 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
@@ -16,56 +16,83 @@ public class TriangleFan {
private final List<Integer> vertices = new LinkedList<>();
private final List<Integer> triangles = new LinkedList<>();
+ /**
+ * Constructor.
+ *
+ * @param facet Mesh facet
+ * @param vert Central vertex for the triangle fan
+ * @throws IllegalArgumentException if some parameter is missing or wrong
+ */
public TriangleFan(MeshFacet facet, int vert) {
if (facet == null) {
throw new IllegalArgumentException("facet");
}
- //System.out.println(facet.getVertices());
-
- CornerTable ct = facet.getCornerTable();
-
- int vertRow = getVertexRow(ct, vert);
- //System.out.println("DDD " + vert + " " + vertRow + " " + ct.getRow(vertRow).getVertexIndex());
+ int vertRow = -1;
+ for (int i = 0; i < facet.getCornerTable().getSize(); i++) {
+ if (facet.getCornerTable().getRow(i).getVertexIndex() == vert) {
+ vertRow = i;
+ break;
+ }
+ }
if (vertRow == -1) {
throw new IllegalArgumentException("vert");
}
- computeOneRingData(ct, vertRow, vertices, triangles);
+ computeOneRingData(facet, vertRow);
}
+ /**
+ * Returns {@code true} if the triangle fan is enclosed (i.e., is not boundary).
+ *
+ * @return {@code true} if the triangle fan is enclosed (i.e., is not boundary).
+ */
public boolean isEnclosed() {
return Objects.equals(vertices.get(0), vertices.get(vertices.size()-1));
}
- public boolean isBoundery() {
+ /**
+ * Returns {@code true} if the triangle fan is boundary (i.e., is not enclosed).
+ *
+ * @return {@code true} if the triangle fan is enclosed (i.e., is not boundary).
+ */
+ public boolean isBoundary() {
return !isEnclosed();
}
+ /**
+ * Vertex indices of the 1-ring neighborhood. If the triangle fan is enclosed
+ * then the first and last vertices are the same.
+ *
+ * @return Vertex indices of the 1-ring neighborhood.
+ */
public List<Integer> getVertices() {
return Collections.unmodifiableList(vertices);
}
+ /**
+ * Triangle indices of the 1-ring neighborhood. The first triangle is related to
+ * the edge delimited by first and second vertex, etc. The number of triangles
+ * corresponds to the number of vertices minus one.
+ *
+ * @return Triangle indices of the 1-ring neighborhood.
+ */
public List<Integer> getTriangles() {
return Collections.unmodifiableList(triangles);
}
- private void computeOneRingData(CornerTable ct, int vertRow, List<Integer> vertices, List<Integer> triangles) {
- //System.out.println(ct);
- //System.out.println("AAA " + ct.getRow(vertRow).getVertexIndex());
+ private void computeOneRingData(MeshFacet facet, int vertRow) {
+ CornerTable ct = facet.getCornerTable();
int ringVertRow = ct.getIndexOfNextCornerInFace(vertRow);
while (true) {
- //System.out.println("BBB " + ct.getRow(ringVertRow).getVertexIndex());
- vertices.add(ct.getRow(ringVertRow).getVertexIndex());
+ vertices.add(facet.getCornerTable().getRow(ringVertRow).getVertexIndex());
if (vertices.size() > 1 && isEnclosed()) {
return; // the ring is closed; we are done
}
- //System.out.println("CCC " + ct.getRow(ct.getIndexOfNextCornerInFace(ringVertRow)).getVertexIndex());
ringVertRow = ct.getIndexOfOppositeCorner(ct.getIndexOfNextCornerInFace(ringVertRow));
if (ringVertRow == -1) {
break; // we reached an open end
}
- //System.out.println("DDD " + ct.getRow(ringVertRow).getVertexIndex());
triangles.add(ct.getIndexOfFace(ringVertRow));
}
@@ -81,15 +108,5 @@ public class TriangleFan {
}
triangles.add(ct.getIndexOfFace(ringVertRow));
}
- }
-
- private int getVertexRow(CornerTable cornerTable, int vert) {
- for (int i = 0; i < cornerTable.getSize(); i++) {
- if (cornerTable.getRow(i).getVertexIndex() == vert) {
- return i;
- }
- }
- return -1;
- }
-
+ }
}
--
GitLab