Skip to content
Snippets Groups Projects
Commit 8228dede authored by Radek Ošlejšek's avatar Radek Ošlejšek
Browse files

Fixed curvature calculations

parent ddc66d6b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 236 additions and 106 deletions
Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/Curvature.java
+ 160
− 42
View file @ 8228dede
......@@ -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));
}
}
}
GUI/src/main/java/cz/fidentis/analyst/tests/EfficiencyTests.java
+ 18
− 18
View file @ 8228dede
......@@ -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) {
......
MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshFacet.java
+ 7
− 6
View file @ 8228dede
......@@ -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);
}
MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshFacetImpl.java
+ 8
− 14
View file @ 8228dede
......@@ -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);
}
}
MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/TriangleFan.java
+ 43
− 26
View file @ 8228dede
......@@ -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;
}
}
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment