diff --git a/Comparison/src/main/java/cz/fidentis/analyst/EfficiencyTests.java b/Comparison/src/main/java/cz/fidentis/analyst/EfficiencyTests.java
index 8f39e0925a8482bbfec40b97521d51689a30aedf..d05318295732334fe3658e012b3b37ddc342db71 100644
--- a/Comparison/src/main/java/cz/fidentis/analyst/EfficiencyTests.java
+++ b/Comparison/src/main/java/cz/fidentis/analyst/EfficiencyTests.java
@@ -6,6 +6,7 @@ import cz.fidentis.analyst.mesh.core.MeshFacet;
import cz.fidentis.analyst.symmetry.Config;
import cz.fidentis.analyst.symmetry.Plane;
import cz.fidentis.analyst.symmetry.SymmetryEstimator;
+import cz.fidentis.analyst.visitors.mesh.BoundingBoxVisitor;
import cz.fidentis.analyst.visitors.mesh.HausdorffDistance;
import cz.fidentis.analyst.visitors.mesh.HausdorffDistance.Strategy;
import java.io.File;
@@ -45,17 +46,17 @@ public class EfficiencyTests {
boolean relativeDist = true;
boolean printDetails = false;
- System.out.println();
- System.out.println("Hausdorff distance sequentially:");
- testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_POINT_ABSOLUTE, false, false, false), printDetails);
- testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_POINT, relativeDist, false, false), printDetails);
- testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_TRIANGLE_APPROXIMATE, relativeDist, false, false), printDetails);
+ //System.out.println();
+ //System.out.println("Hausdorff distance sequentially:");
+ //testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_POINT_ABSOLUTE, false, false, false), printDetails);
+ //testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_POINT, relativeDist, false, false), printDetails);
+ //testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_TRIANGLE_APPROXIMATE, relativeDist, false, false), printDetails);
- System.out.println();
- System.out.println("Hausdorff distance consurrently:");
- testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_POINT_ABSOLUTE, false, false, true), printDetails);
- testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_POINT, relativeDist, false, true), printDetails);
- testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_TRIANGLE_APPROXIMATE, relativeDist, false, true), printDetails);
+ //System.out.println();
+ //System.out.println("Hausdorff distance consurrently:");
+ //testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_POINT_ABSOLUTE, false, false, true), printDetails);
+ //testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_POINT, relativeDist, false, true), printDetails);
+ //testAndPrint(face1, new HausdorffDistance(face2.getMeshModel(), Strategy.POINT_TO_TRIANGLE_APPROXIMATE, relativeDist, false, true), printDetails);
}
protected static void testAndPrint(HumanFace face, HausdorffDistance vis, boolean printDetails) {
@@ -87,7 +88,7 @@ public class EfficiencyTests {
private static long measureSymmetryPlane(HumanFace face, boolean printDetails) {
long startTime = System.currentTimeMillis();
- SymmetryEstimator est = new SymmetryEstimator(face.getMeshModel().getFacets().get(0), new Config());
+ SymmetryEstimator est = new SymmetryEstimator(face.getMeshModel().getFacets().get(0), Config.getDefault());
Plane plane = est.getApproxSymmetryPlane(null);
long retTime = System.currentTimeMillis() - startTime;
diff --git a/Comparison/src/main/java/cz/fidentis/analyst/symmetry/SymmetryEstimator.java b/Comparison/src/main/java/cz/fidentis/analyst/symmetry/SymmetryEstimator.java
index 61a773781b87b698bc3790d0c6f721f7248209e3..d4025e8e235b6be48cf81eb99cf656f4977ff27b 100644
--- a/Comparison/src/main/java/cz/fidentis/analyst/symmetry/SymmetryEstimator.java
+++ b/Comparison/src/main/java/cz/fidentis/analyst/symmetry/SymmetryEstimator.java
@@ -8,6 +8,9 @@ import cz.fidentis.analyst.mesh.core.MeshFacetImpl;
import cz.fidentis.analyst.mesh.core.MeshPointImpl;
import cz.fidentis.analyst.mesh.core.MeshTriangle;
import cz.fidentis.analyst.visitors.mesh.BoundingBoxVisitor;
+import cz.fidentis.analyst.visitors.mesh.CurvatureVisitor;
+import cz.fidentis.analyst.visitors.mesh.GaussianCurvatureVisitor;
+import cz.fidentis.analyst.visitors.mesh.MaxCurvatureVisitor;
import cz.fidentis.analyst.visitors.mesh.TriangleListVisitor;
import java.util.ArrayList;
import java.util.Collections;
@@ -31,7 +34,7 @@ public class SymmetryEstimator {
private final MeshFacet facet; // Facet of the model on which symmetry is computed
- private TriangleVertexAreas[] areas; // Representation for areas of Voronoi region of triangles
+ private List<Vector3d> areas; // Representation for areas of Voronoi region of triangles
private BoundingBox boundingBox; // Represent min-max box. It is automatically maintained by given point array.
@@ -49,10 +52,10 @@ public class SymmetryEstimator {
this.facet = f;
this.config = config;
- this.areas = new TriangleVertexAreas[f.getNumTriangles()];
+ this.areas = new ArrayList<>(f.getNumTriangles());
int i = 0;
for (MeshTriangle tri: f) {
- areas[i++] = computeTriangleVertexAreas(tri);
+ areas.add(tri.getVoronoiPoint());
}
TriangleListVisitor vis = new TriangleListVisitor(false);
@@ -94,7 +97,7 @@ public class SymmetryEstimator {
*/
public BoundingBox getBoundingBox() {
if (boundingBox == null) {
- BoundingBoxVisitor visitor = new BoundingBoxVisitor(true);
+ BoundingBoxVisitor visitor = new BoundingBoxVisitor(false);
facet.accept(visitor);
boundingBox = visitor.getBoundingBox();
}
@@ -103,36 +106,33 @@ public class SymmetryEstimator {
/**
* Computes the approximate plane of symmetry.
+ * TO DO: ZRUSIT VSTUPNI PARAMETR!!!
+ *
* @param panel parrent window for the progress window (can be null)
* @return approximate plane of symmtetry
*/
public Plane getApproxSymmetryPlane(JPanel panel) {
- List<ApproxSymmetryPlane> planes = new ArrayList<>();
- //List<Vector3d> normals = calculateNormals();
-
if (!facet.hasVertexNormals()) {
facet.calculateVertexNormals();
}
- double[] curvatures = initCurvatures();
- List<Double> sortedCurvatures = new ArrayList<>();
- for (int i = 0; i < curvatures.length; i++) {
- sortedCurvatures.add(curvatures[i]);
- }
+ List<ApproxSymmetryPlane> planes = new ArrayList<>();
+
+ List<Double> curvatures = initCurvatures(facet);
+ List<Double> sortedCurvatures = new ArrayList<>(curvatures);
Collections.sort(sortedCurvatures);
- if(config.getSignificantPointCount() > facet.getNumberOfVertices()) {
+ if (config.getSignificantPointCount() > facet.getNumberOfVertices()) {
config.setSignificantPointCount((facet.getNumberOfVertices()) - 1);
}
double bottomCurvature = sortedCurvatures.get(sortedCurvatures.size() - 1 - config.getSignificantPointCount());
List<Integer> significantPoints = new ArrayList<>();
List<Double> significantCurvatures = new ArrayList<>();
-
for (int i = 0; i < facet.getNumberOfVertices(); i++) {
- if (curvatures[i] >= bottomCurvature) {
- significantCurvatures.add(curvatures[i]);
+ if (curvatures.get(i) >= bottomCurvature) {
+ significantCurvatures.add(curvatures.get(i));
significantPoints.add(i);
}
}
@@ -140,74 +140,71 @@ public class SymmetryEstimator {
Plane plane = new Plane(0, 0, 0, 0);
int lastVotes = 0;
- ProgressMonitor progressMonitor = null;
- if (panel != null) {
- UIManager.put("ProgressMonitor.progressText", "Counting symmetry...");
- progressMonitor = new ProgressMonitor(panel, "Counting...", "Steps", 0, significantPoints.size());
- }
-
- double onePercent = significantCurvatures.size() / 100.0;
- double percentsPerStep = 1 / onePercent;
for (int i = 0; i < significantPoints.size(); i++) {
for (int j = 0; j < significantPoints.size(); j++) {
- if (i != j) {
- double minRatio = config.getMinCurvRatio();
- double maxRatio = 1.0 / minRatio;
- if (significantCurvatures.get(i) / significantCurvatures.get(j) >= minRatio &&
- significantCurvatures.get(i) / significantCurvatures.get(j) <= maxRatio) {
-
- Vector3d p1 = new Vector3d(facet.getVertex(significantPoints.get(i)).getPosition());
- Vector3d p2 = new Vector3d(facet.getVertex(significantPoints.get(j)).getPosition());
+ if (i == j) {
+ continue;
+ }
+
+ double minRatio = config.getMinCurvRatio();
+ double maxRatio = 1.0 / minRatio;
+ if (significantCurvatures.get(i) / significantCurvatures.get(j) < minRatio ||
+ significantCurvatures.get(i) / significantCurvatures.get(j) > maxRatio) {
+ continue;
+ }
+
+ MeshPoint meshPointI = facet.getVertex(significantPoints.get(i));
+ MeshPoint meshPointJ = facet.getVertex(significantPoints.get(j));
+
+ Vector3d p1 = new Vector3d(meshPointI.getPosition());
+ Vector3d p2 = new Vector3d(meshPointJ.getPosition());
- Vector3d avrg = new Vector3d(p1);
- avrg.add(p2);
- avrg.scale(0.5);
+ Vector3d avrg = new Vector3d(p1);
+ avrg.add(p2);
+ avrg.scale(0.5);
- Vector3d normal = new Vector3d(p1);
- normal.sub(p2);
- normal.normalize();
+ Vector3d normal = new Vector3d(p1);
+ normal.sub(p2);
+ normal.normalize();
- double d = -(normal.x * avrg.x) - (normal.y * avrg.y) - (normal.z * avrg.z);
+ double d = -(normal.x * avrg.x) - (normal.y * avrg.y) - (normal.z * avrg.z);
- Vector3d ni = new Vector3d(facet.getVertex(significantPoints.get(i)).getNormal());
- Vector3d nj = new Vector3d(facet.getVertex(significantPoints.get(j)).getNormal());
- ni.normalize();
- nj.normalize();
+ Vector3d ni = new Vector3d(meshPointI.getNormal());
+ Vector3d nj = new Vector3d(meshPointI.getNormal());
+ ni.normalize();
+ nj.normalize();
- Vector3d normVec = ni;
- normVec.sub(nj);
- normVec.normalize();
- double normCos = normVec.dot(normal);
+ Vector3d normVec = new Vector3d(ni);
+ normVec.sub(nj);
+ normVec.normalize();
+ double normCos = normVec.dot(normal);
- if (Math.abs(normCos) >= config.getMinNormAngleCos()) {
- Plane newPlane = new Plane(normal.x, normal.y, normal.z, d);
- int currentVotes = getVotes(newPlane,
- significantCurvatures,
- significantPoints,
- config.getMinCurvRatio(),
- config.getMinAngleCos(),
- config.getMinNormAngleCos(),
- getBoundingBox().getMaxDiag() * config.getMaxRelDistance());
+ if (Math.abs(normCos) < config.getMinNormAngleCos()) {
+ continue;
+ }
+
+ Plane newPlane = new Plane(normal.x, normal.y, normal.z, d);
+ int currentVotes = getVotes(
+ newPlane,
+ significantCurvatures,
+ significantPoints,
+ config.getMinCurvRatio(),
+ config.getMinAngleCos(),
+ config.getMinNormAngleCos(),
+ getBoundingBox().getMaxDiag() * config.getMaxRelDistance()
+ );
- planes.add(new ApproxSymmetryPlane(newPlane, currentVotes));
-
- if (currentVotes > lastVotes) {
- lastVotes = currentVotes;
- plane = newPlane;
- }
- }
-
- }
+ planes.add(new ApproxSymmetryPlane(newPlane, currentVotes));
+
+ if (currentVotes > lastVotes) {
+ lastVotes = currentVotes;
+ plane = newPlane;
}
}
- if (panel != null) {
- progressMonitor.setNote("Task step: " + (int) ((i + 1) * percentsPerStep));
- progressMonitor.setProgress(i);
- }
}
Collections.sort(planes);
- ArrayList<ApproxSymmetryPlane> finalPlanes = new ArrayList<>();
+ List<ApproxSymmetryPlane> finalPlanes = new ArrayList<>();
for (int i = 0; i < planes.size(); i++) {
if (planes.get(i).getVotes() == lastVotes){
finalPlanes.add(planes.get(i));
@@ -217,27 +214,47 @@ public class SymmetryEstimator {
if (config.isAveraging()){
plane = finalPlane;
}
- if (panel != null) {
- JOptionPane.showMessageDialog(panel, "Symmetry estimate done.", "Done", 0,
- new ImageIcon(getClass().getResource("/cz/fidentis/analyst/gui/resources/exportedModel.png")));
- progressMonitor.close();
- }
+
return plane;
}
- private double[] initCurvatures() {
- double[] curvatures = new double[facet.getNumberOfVertices()];
+ private List<Double> initCurvatures(MeshFacet facet) {
+ CurvatureVisitor vis;
+ if (facet.getNumberOfVertices() < 2500) {
+ // searching for Maximum curvature in vertex using Gaussian curvature
+ // and Mean curvature leads to longer calculation but better results for real face models
+ //curvatures.add(getMaxCurvature(i));
+ vis = new MaxCurvatureVisitor(false);
+ } else {
+ //curvatures.add(getGaussianCurvature(i));
+ vis = new GaussianCurvatureVisitor(false);
+ }
+ facet.accept(vis);
+ List<Double> curvatures = new ArrayList<>(vis.getCurvatures().get(facet));
+ for (int i = 0; i < curvatures.size(); i++) {
+ if (Double.isNaN(curvatures.get(i))) {
+ curvatures.set(i, Double.MIN_VALUE); // !!!!
+ }
+ }
+ return curvatures;
+
+ /*
for (int i = 0; i < facet.getNumberOfVertices(); i++) {
- if (facet.getNumberOfVertices() == 2500) {
- curvatures[i] = this.getMaxCurvature(i);
+ if (facet.getNumberOfVertices() < 2500) {
+ // searching for Maximum curvature in vertex using Gaussian curvature
+ // and Mean curvature leads to longer calculation but better results for real face models
+ //curvatures.add(getMaxCurvature(i));
+ vis = new MaxCurvatureVisitor(false);
} else {
- curvatures[i] = this.getGaussianCurvature(i);
+ //curvatures.add(getGaussianCurvature(i));
+ vis = new GaussianCurvatureVisitor(false);
}
- if (Double.isNaN(curvatures[i])){
- curvatures[i] = Double.MIN_VALUE;
+ if (Double.isNaN(curvatures.get(i))){
+ curvatures.set(i, Double.MIN_VALUE);
}
}
return curvatures;
+ */
}
private Plane computeNewPlane(List<ApproxSymmetryPlane> finalPlanes) {
@@ -364,302 +381,6 @@ public class SymmetryEstimator {
return facet;
}
- /**
- * Representation of areas of Voronoi region of triangles
- * used for computing Gaussian curvature
- *
- * @author Natália Bebjaková
- */
- private final class TriangleVertexAreas {
- public final double v1Area;
- public final double v2Area;
- public final double v3Area;
-
- private TriangleVertexAreas(double v1, double v2, double v3) {
- v1Area = v1;
- v2Area = v2;
- v3Area = v3;
- }
- }
-
- /**
- *
- * @param pointIndex index of point for which we are searching traingles in its neighborhood
- * @return triangle neighbours of given index of vertex
- */
- protected List<Integer> getNeighbours(int pointIndex) {
- return facet.getCornerTable().getTriangleIndexesByVertexIndex(pointIndex);
- }
-
- /**
- * Return area of Voronoi region of triangle
- *
- * @param t triangle of which area is computed
- * @return area of Voronoi region
- */
- private TriangleVertexAreas computeTriangleVertexAreas(MeshTriangle t) {
- double a = (t.vertex2.subtractPosition(t.vertex3)).abs();
- double b = (t.vertex3.subtractPosition(t.vertex1)).abs();
- double c = (t.vertex2.subtractPosition(t.vertex1)).abs();
-
- double d1 = a * a * (b * b + c * c - a * a);
- double d2 = b * b * (c * c + a * a - b * b);
- double d3 = c * c * (a * a + b * b - c * c);
- double dSum = d1 + d2 + d3;
-
- d1 /= dSum;
- d2 /= dSum;
- d3 /= dSum;
-
- MeshPoint v1Half = (t.vertex2.addPosition(t.vertex3)).dividePosition(2);
- MeshPoint v2Half = (t.vertex1.addPosition(t.vertex3)).dividePosition(2);
- MeshPoint v3Half = (t.vertex2.addPosition(t.vertex1)).dividePosition(2);
-
-
- //TriangleVertexAreas area = new TriangleVertexAreas();
-
- if (d1 < 0) {
- double v3Area = ((v2Half.subtractPosition(t.vertex3)).crossProduct(v1Half.subtractPosition(t.vertex3))).abs() / 2.0;
- double v2Area = ((v3Half.subtractPosition(t.vertex2)).crossProduct(v1Half.subtractPosition(t.vertex2))).abs() / 2.0;
- double v1Area = (((v1Half.subtractPosition(t.vertex1)).crossProduct(v3Half.subtractPosition(t.vertex1))).abs() / 2.0) +
- (((v1Half.subtractPosition(t.vertex1)).crossProduct(v2Half.subtractPosition(t.vertex1))).abs() / 2.0);
- return new TriangleVertexAreas(v1Area, v2Area, v3Area);
- }
- if (d2 < 0) {
- double v1Area = ((v3Half.subtractPosition(t.vertex1)).crossProduct(v2Half.subtractPosition(t.vertex1))).abs() / 2.0;
- double v3Area = ((v1Half.subtractPosition(t.vertex3)).crossProduct(v2Half.subtractPosition(t.vertex3))).abs() / 2.0;
- double v2Area = (((v2Half.subtractPosition(t.vertex2)).crossProduct(v1Half.subtractPosition(t.vertex2))).abs() / 2.0) +
- (((v2Half.subtractPosition(t.vertex2)).crossProduct(v3Half.subtractPosition(t.vertex2))).abs() / 2.0);
- return new TriangleVertexAreas(v1Area, v2Area, v3Area);
- }
- if (d3 < 0) {
- double v2Area = ((v1Half.subtractPosition(t.vertex2)).crossProduct(v3Half.subtractPosition(t.vertex2))).abs() / 2.0;
- double v1Area = ((v2Half.subtractPosition(t.vertex1)).crossProduct(v3Half.subtractPosition(t.vertex1))).abs() / 2.0;
- double v3Area = (((v3Half.subtractPosition(t.vertex3)).crossProduct(v2Half.subtractPosition(t.vertex3))).abs() / 2.0) +
- (((v3Half.subtractPosition(t.vertex3)).crossProduct(v1Half.subtractPosition(t.vertex3))).abs() / 2.0);
- return new TriangleVertexAreas(v1Area, v2Area, v3Area);
- }
-
- MeshPoint circumcenter = t.vertex1.multiplyPosition(d1).addPosition(t.vertex2.multiplyPosition(d2).addPosition(t.vertex3.multiplyPosition(d3)));
-
- double v1Area = (((v2Half.subtractPosition(t.vertex1)).crossProduct(circumcenter.subtractPosition(t.vertex1))).abs() / 2.0) +
- (((v3Half.subtractPosition(t.vertex1)).crossProduct(circumcenter.subtractPosition(t.vertex1))).abs() / 2.0);
-
- double v2Area = (((v3Half.subtractPosition(t.vertex2)).crossProduct(circumcenter.subtractPosition(t.vertex2))).abs() / 2.0) +
- (((v1Half.subtractPosition(t.vertex2)).crossProduct(circumcenter.subtractPosition(t.vertex2))).abs() / 2.0);
-
- double v3Area = (((v1Half.subtractPosition(t.vertex3)).crossProduct(circumcenter.subtractPosition(t.vertex3))).abs() / 2.0) +
- (((v2Half.subtractPosition(t.vertex3)).crossProduct(circumcenter.subtractPosition(t.vertex3))).abs() / 2.0);
- return new TriangleVertexAreas(v1Area, v2Area, v3Area);
- }
-
- /**
- * Calculates angle of the triangle according to center point
- *
- * @param centerPoint center point
- * @param t triangle
- * @return angle
- */
- private double calculateTriangleAlfa(MeshPoint centerPoint, MeshTriangle t) {
- double alfaCos = 0;
-
- if (t.vertex1 == centerPoint) {
- MeshPoint e1 = t.vertex1.subtractPosition(t.vertex3);
- e1 = e1.dividePosition(e1.abs());
-
- MeshPoint e2 = t.vertex2.subtractPosition(t.vertex3);
- e2 = e2.dividePosition(e2.abs());
-
- alfaCos = e1.dotProduct(e2);
- } else if (t.vertex2 == centerPoint) {
- MeshPoint e1 = t.vertex2.subtractPosition(t.vertex1);
- e1 = e1.dividePosition(e1.abs());
-
- MeshPoint e2 = t.vertex3.subtractPosition(t.vertex1);
- e2 = e2.dividePosition(e2.abs());
-
- alfaCos = e1.dotProduct(e2);
- } else if (t.vertex3 == centerPoint){
- MeshPoint e1 = t.vertex3.subtractPosition(t.vertex2);
- e1 = e1.dividePosition(e1.abs());
-
- MeshPoint e2 = t.vertex1.subtractPosition(t.vertex2);
- e2 = e2.dividePosition(e2.abs());
-
- alfaCos = e1.dotProduct(e2);
- }
- double alfa = Math.acos(alfaCos);
- return 1 / Math.tan(alfa);
- }
-
- /**
- * Calculates angle of the triangle according to center point
- *
- * @param centerPoint center point
- * @param t triangle
- * @return angle
- */
- private double calculateTriangleBeta(MeshPoint centerPoint, MeshTriangle t) {
- double betaCos = 0;
-
- if (t.vertex1 == centerPoint) {
- MeshPoint e1 = t.vertex1.subtractPosition(t.vertex2);
- e1 = e1.dividePosition(e1.abs());
-
- MeshPoint e2 = t.vertex3.subtractPosition(t.vertex2);
- e2 = e2.dividePosition(e2.abs());
-
- betaCos = e1.dotProduct(e2);
- } else if (t.vertex2 == centerPoint) {
- MeshPoint e1 = t.vertex2.subtractPosition(t.vertex3);
- e1 = e1.dividePosition(e1.abs());
-
- MeshPoint e2 = t.vertex1.subtractPosition(t.vertex3);
- e2 = e2.dividePosition(e2.abs());
-
- betaCos = e1.dotProduct(e2);
- } else if (t.vertex3 == centerPoint){
- MeshPoint e1 = t.vertex3.subtractPosition(t.vertex1);
- e1 = e1.dividePosition(e1.abs());
-
- MeshPoint e2 = t.vertex2.subtractPosition(t.vertex1);
- e2 = e2.dividePosition(e2.abs());
-
- betaCos = e1.dotProduct(e2);
- }
- double beta = Math.acos(betaCos);
- return 1 / Math.tan(beta);
- }
-
- /**
- * Calculates Laplacian
- *
- * @param centerIndex centerIndex
- * @return laplacian
- */
- private Vector3d calculateLaplacian(int centerIndex) {
- List<Integer> trianglesNeighbours = getNeighbours(centerIndex);
-
- if (trianglesNeighbours.isEmpty()) {
- return new Vector3d();
- }
-
- double areaSum = 0;
- Vector3d pointSum = new Vector3d();
- for (int i = 0; i < trianglesNeighbours.size(); i++) {
- MeshTriangle t = triangles.get(trianglesNeighbours.get(i));
- MeshTriangle tNext = triangles.get(trianglesNeighbours.get((i + 1) % trianglesNeighbours.size()));
- double area = 0;
-
- //MeshPoint tVertex1 = facet.getVertices().get(t.vertex1);
- //MeshPoint tVertex2 = facet.getVertices().get(t.vertex2);
- //MeshPoint tVertex3 = facet.getVertices().get(t.vertex3);
- MeshPoint centerPoint = facet.getVertices().get(centerIndex);
-
- Vector3d aux = new Vector3d();
- if (t.vertex1 == facet.getVertex(centerIndex)) {
- area = areas[trianglesNeighbours.get(i)].v1Area;
- aux = new Vector3d(t.vertex2.getPosition());
- } else if (t.vertex2 == facet.getVertex(centerIndex)) {
- area = areas[trianglesNeighbours.get(i)].v2Area;
- aux = new Vector3d(t.vertex3.getPosition());
- } else if (t.vertex3 == facet.getVertex(centerIndex)) {
- area = areas[trianglesNeighbours.get(i)].v3Area;
- aux = new Vector3d(t.vertex1.getPosition());
- }
- aux.sub(centerPoint.getPosition());
- aux.scale((calculateTriangleAlfa(centerPoint, t) + calculateTriangleBeta(centerPoint, tNext)) / 2.0);
- pointSum.add(aux);
- areaSum += area;
- }
- pointSum.scale(1.0/areaSum);
- return pointSum;
- }
-
- /**
- * Gaussian curvature in a vertex of a triangle mesh.
- * It can only be estimated because triangle mesh is not a continuous
- * but discrete representation of surface.
- *
- * @param centerIndex index of vertex in which gaussian curvature is computed
- * @return Gaussian curvature in given vertex
- */
- private double getGaussianCurvature(int centerIndex) {
-
- List<Integer> triangleNeighbours = this.getNeighbours(centerIndex);
-
- if (triangleNeighbours.isEmpty()) {
- return Double.NaN;
- }
- double sum = 2 * Math.PI;
- double areaSum = 0;
-
- for (int i = 0; i < triangleNeighbours.size(); i++) {
- Vector3d v1 = new Vector3d();
- Vector3d v2 = new Vector3d();
- MeshTriangle t = triangles.get(triangleNeighbours.get(i));
-
- double area = 0;
- if (t.vertex1 == facet.getVertex(centerIndex)) {
- v1 = new Vector3d(t.vertex2.getPosition());
- v2 = new Vector3d(t.vertex3.getPosition());
- v1.sub(t.vertex1.getPosition());
- v2.sub(t.vertex1.getPosition());
-
- area = areas[triangleNeighbours.get(i)].v1Area;
- } else if (t.vertex2 == facet.getVertex(centerIndex)) {
- v1 = new Vector3d(t.vertex1.getPosition());
- v2 = new Vector3d(t.vertex3.getPosition());
- v1.sub(t.vertex2.getPosition());
- v2.sub(t.vertex2.getPosition());
-
- area = areas[triangleNeighbours.get(i)].v2Area;
- } else if (t.vertex3 == facet.getVertex(centerIndex)) {
- v1 = new Vector3d(t.vertex2.getPosition());
- v2 = new Vector3d(t.vertex1.getPosition());
- v1.sub(t.vertex3.getPosition());
- v2.sub(t.vertex3.getPosition());
-
- area = areas[triangleNeighbours.get(i)].v3Area;
- }
-
- areaSum += area;
- v1.normalize();
- v2.normalize();
-
- sum -= Math.acos(v1.dot(v2));
- }
- double value = sum * (1.0 / areaSum);
- return value;
- }
-
- /**
- *
- * @param centerIndex center index
- * @return mean curvature
- */
- private double getMeanCurvature(int centerIndex) {
- return calculateLaplacian(centerIndex).length();
- }
-
- /**
- *
- * @param centerIndex center index
- * @return max curvature
- */
- private double getMaxCurvature(int centerIndex) {
- double gaussianCurvature = getGaussianCurvature(centerIndex);
- if (Double.isNaN(gaussianCurvature)) {
- return Double.NaN;
- }
- double meanCurvature = getMeanCurvature(centerIndex);
- double meanCurvatureSqr = meanCurvature * meanCurvature;
- if (meanCurvatureSqr <= gaussianCurvature) {
- return meanCurvature;
- }
- return meanCurvature + Math.sqrt(meanCurvatureSqr - gaussianCurvature);
- }
-
/**
* Computes votes for given plane
*
@@ -679,6 +400,7 @@ public class SymmetryEstimator {
double minAngleCos,
double minNormAngleCos,
double maxDist) {
+
plane.normalize();
Vector3d normal = plane.getNormal();
@@ -687,9 +409,9 @@ public class SymmetryEstimator {
int votes = 0;
//List<Vector3d> normals = calculateNormals();
- if (!facet.hasVertexNormals()) {
- facet.calculateVertexNormals();
- }
+ //if (!facet.hasVertexNormals()) {
+ // facet.calculateVertexNormals();
+ //}
for (int i = 0; i < curvatures.size(); i++) {
for (int j = 0; j < curvatures.size(); j++) {
diff --git a/Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/CurvatureVisitor.java b/Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/CurvatureVisitor.java
new file mode 100644
index 0000000000000000000000000000000000000000..a6f2c06ddb4649a6e893561dc90bcad708a3e6c4
--- /dev/null
+++ b/Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/CurvatureVisitor.java
@@ -0,0 +1,57 @@
+/*
+ * To change this license header, choose License Headers in Project Properties.
+ * To change this template file, choose Tools | Templates
+ * and open the template in the editor.
+ */
+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.symmetry.SymmetryEstimator;
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import javax.vecmath.Vector3d;
+
+/**
+ *
+ * @author oslejsek
+ */
+public abstract class CurvatureVisitor extends MeshVisitor {
+
+ private Map<MeshFacet, List<Double>> curvatures = new HashMap<>();
+
+ /**
+ * Constructor.
+ *
+ * @param asynchronous If {@code true}, then asynchronous visitor is created.
+ */
+ public CurvatureVisitor(boolean asynchronous) {
+ super(asynchronous);
+ }
+
+ @Override
+ protected synchronized void visitMeshFacet(MeshFacet facet) {
+ if (curvatures.containsKey(facet)) {
+ return; // already visited facet
+ }
+
+ curvatures.put(facet, new ArrayList<>());
+
+ List<MeshTriangle> triangles = facet.getTriangles();
+ for (int i = 0; i < facet.getNumberOfVertices(); i++) {
+ curvatures.get(facet).add(calculateCurvature(facet, triangles, i));
+ }
+ }
+
+ public Map<MeshFacet, List<Double>> getCurvatures() {
+ return Collections.unmodifiableMap(curvatures);
+ }
+
+ protected abstract double calculateCurvature(MeshFacet facet, List<MeshTriangle> triangles, int centerIndex);
+
+}
diff --git a/Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/GaussianCurvatureVisitor.java b/Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/GaussianCurvatureVisitor.java
new file mode 100644
index 0000000000000000000000000000000000000000..f6dd505146f0f9b029587301f55c0d0902a7f7b3
--- /dev/null
+++ b/Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/GaussianCurvatureVisitor.java
@@ -0,0 +1,77 @@
+/*
+ * To change this license header, choose License Headers in Project Properties.
+ * To change this template file, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package cz.fidentis.analyst.visitors.mesh;
+
+import cz.fidentis.analyst.mesh.core.MeshFacet;
+import cz.fidentis.analyst.mesh.core.MeshTriangle;
+import java.util.List;
+import javax.vecmath.Vector3d;
+
+/**
+ *
+ * @author oslejsek
+ */
+public class GaussianCurvatureVisitor extends CurvatureVisitor {
+
+ /**
+ * Constructor.
+ *
+ * @param asynchronous If {@code true}, then asynchronous visitor is created.
+ */
+ public GaussianCurvatureVisitor(boolean asynchronous) {
+ super(asynchronous);
+ }
+
+ protected double calculateCurvature(MeshFacet facet, List<MeshTriangle> triangles, int centerIndex) {
+ List<Integer> neighbouringTriangles = facet.getCornerTable().getTriangleIndexesByVertexIndex(centerIndex);
+
+ if (neighbouringTriangles.isEmpty()) {
+ return Double.NaN;
+ }
+
+ double sum = 2 * Math.PI;
+ double areaSum = 0;
+
+ // fro all surrounding triangles:
+ List<Vector3d> voronoiPoints = facet.calculateVoronoiPoints();
+ for (int i = 0; i < neighbouringTriangles.size(); i++) {
+ Vector3d v1 = new Vector3d();
+ Vector3d v2 = new Vector3d();
+ MeshTriangle tri = triangles.get(neighbouringTriangles.get(i));
+ Vector3d voronoiPoint = voronoiPoints.get(neighbouringTriangles.get(i));
+
+ double area = 0;
+ if (tri.vertex1 == facet.getVertex(centerIndex)) {
+ v1 = new Vector3d(tri.vertex2.getPosition());
+ v2 = new Vector3d(tri.vertex3.getPosition());
+ v1.sub(tri.vertex1.getPosition());
+ v2.sub(tri.vertex1.getPosition());
+ area = voronoiPoint.x;
+ } else if (tri.vertex2 == facet.getVertex(centerIndex)) {
+ v1 = new Vector3d(tri.vertex1.getPosition());
+ v2 = new Vector3d(tri.vertex3.getPosition());
+ v1.sub(tri.vertex2.getPosition());
+ v2.sub(tri.vertex2.getPosition());
+ area = voronoiPoint.y;
+ } else if (tri.vertex3 == facet.getVertex(centerIndex)) {
+ v1 = new Vector3d(tri.vertex2.getPosition());
+ v2 = new Vector3d(tri.vertex1.getPosition());
+ v1.sub(tri.vertex3.getPosition());
+ v2.sub(tri.vertex3.getPosition());
+ area = voronoiPoint.z;
+ }
+
+ areaSum += area;
+ v1.normalize();
+ v2.normalize();
+
+ sum -= Math.acos(v1.dot(v2));
+ }
+ double value = sum * (1.0 / areaSum);
+
+ return value;
+ }
+}
diff --git a/Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/MaxCurvatureVisitor.java b/Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/MaxCurvatureVisitor.java
new file mode 100644
index 0000000000000000000000000000000000000000..d5018e0837fa3e36c4b2a30747bdea7d87c735c7
--- /dev/null
+++ b/Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/MaxCurvatureVisitor.java
@@ -0,0 +1,168 @@
+package cz.fidentis.analyst.visitors.mesh;
+
+import cz.fidentis.analyst.mesh.core.MeshFacet;
+import cz.fidentis.analyst.mesh.core.MeshPoint;
+import cz.fidentis.analyst.mesh.core.MeshTriangle;
+import java.util.List;
+import javax.vecmath.Vector3d;
+
+/**
+ *
+ * @author oslejsek
+ */
+public class MaxCurvatureVisitor extends GaussianCurvatureVisitor{
+
+ /**
+ * Constructor.
+ *
+ * @param asynchronous If {@code true}, then asynchronous visitor is created.
+ */
+ public MaxCurvatureVisitor(boolean asynchronous) {
+ super(asynchronous);
+ }
+
+ protected double calculateCurvature(MeshFacet facet, List<MeshTriangle> triangles, int centerIndex) {
+ double gaussianCurvature = super.calculateCurvature(facet, triangles, centerIndex);
+ if (Double.isNaN(gaussianCurvature)) {
+ return Double.NaN;
+ }
+ double meanCurvature = getMeanCurvature(facet, triangles, centerIndex);
+ double meanCurvatureSqr = meanCurvature * meanCurvature;
+ if (meanCurvatureSqr <= gaussianCurvature) {
+ return meanCurvature;
+ }
+ return meanCurvature + Math.sqrt(meanCurvatureSqr - gaussianCurvature);
+ }
+
+ /**
+ *
+ * @param centerIndex center index
+ * @return mean curvature
+ */
+ private double getMeanCurvature(MeshFacet facet, List<MeshTriangle> triangles, int centerIndex) {
+ return calculateLaplacian(facet, triangles, centerIndex).length();
+ }
+
+ /**
+ * Calculates Laplacian
+ *
+ * @param centerIndex centerIndex
+ * @return laplacian
+ */
+ private Vector3d calculateLaplacian(MeshFacet facet, List<MeshTriangle> triangles, int centerIndex) {
+ List<Integer> neighbouringTriangles = facet.getCornerTable().getTriangleIndexesByVertexIndex(centerIndex);
+
+ if (neighbouringTriangles.isEmpty()) {
+ return new Vector3d();
+ }
+
+ double areaSum = 0;
+ Vector3d pointSum = new Vector3d();
+ List<Vector3d> voronoiPoints = facet.calculateVoronoiPoints();
+
+ for (int i = 0; i < neighbouringTriangles.size(); i++) {
+ MeshTriangle t = triangles.get(neighbouringTriangles.get(i));
+ MeshTriangle tNext = triangles.get(neighbouringTriangles.get((i + 1) % neighbouringTriangles.size()));
+ double area = 0;
+
+ MeshPoint centerPoint = facet.getVertices().get(centerIndex);
+ Vector3d voronoiPoint = voronoiPoints.get(neighbouringTriangles.get(i));
+
+ Vector3d aux = new Vector3d();
+ if (t.vertex1 == facet.getVertex(centerIndex)) {
+ area = voronoiPoint.x;
+ aux = new Vector3d(t.vertex2.getPosition());
+ } else if (t.vertex2 == facet.getVertex(centerIndex)) {
+ area = voronoiPoint.y;
+ aux = new Vector3d(t.vertex3.getPosition());
+ } else if (t.vertex3 == facet.getVertex(centerIndex)) {
+ area = voronoiPoint.z;
+ aux = new Vector3d(t.vertex1.getPosition());
+ }
+ aux.sub(centerPoint.getPosition());
+ aux.scale((calculateTriangleAlfa(centerPoint, t) + calculateTriangleBeta(centerPoint, tNext)) / 2.0);
+ pointSum.add(aux);
+ areaSum += area;
+ }
+ pointSum.scale(1.0/areaSum);
+ return pointSum;
+ }
+
+ /**
+ * Calculates angle of the triangle according to center point
+ *
+ * @param centerPoint center point
+ * @param t triangle
+ * @return angle
+ */
+ private double calculateTriangleAlfa(MeshPoint centerPoint, MeshTriangle t) {
+ double alfaCos = 0;
+
+ if (t.vertex1 == centerPoint) {
+ MeshPoint e1 = t.vertex1.subtractPosition(t.vertex3);
+ e1 = e1.dividePosition(e1.abs());
+
+ MeshPoint e2 = t.vertex2.subtractPosition(t.vertex3);
+ e2 = e2.dividePosition(e2.abs());
+
+ alfaCos = e1.dotProduct(e2);
+ } else if (t.vertex2 == centerPoint) {
+ MeshPoint e1 = t.vertex2.subtractPosition(t.vertex1);
+ e1 = e1.dividePosition(e1.abs());
+
+ MeshPoint e2 = t.vertex3.subtractPosition(t.vertex1);
+ e2 = e2.dividePosition(e2.abs());
+
+ alfaCos = e1.dotProduct(e2);
+ } else if (t.vertex3 == centerPoint){
+ MeshPoint e1 = t.vertex3.subtractPosition(t.vertex2);
+ e1 = e1.dividePosition(e1.abs());
+
+ MeshPoint e2 = t.vertex1.subtractPosition(t.vertex2);
+ e2 = e2.dividePosition(e2.abs());
+
+ alfaCos = e1.dotProduct(e2);
+ }
+ double alfa = Math.acos(alfaCos);
+ return 1 / Math.tan(alfa);
+ }
+
+ /**
+ * Calculates angle of the triangle according to center point
+ *
+ * @param centerPoint center point
+ * @param t triangle
+ * @return angle
+ */
+ private double calculateTriangleBeta(MeshPoint centerPoint, MeshTriangle t) {
+ double betaCos = 0;
+
+ if (t.vertex1 == centerPoint) {
+ MeshPoint e1 = t.vertex1.subtractPosition(t.vertex2);
+ e1 = e1.dividePosition(e1.abs());
+
+ MeshPoint e2 = t.vertex3.subtractPosition(t.vertex2);
+ e2 = e2.dividePosition(e2.abs());
+
+ betaCos = e1.dotProduct(e2);
+ } else if (t.vertex2 == centerPoint) {
+ MeshPoint e1 = t.vertex2.subtractPosition(t.vertex3);
+ e1 = e1.dividePosition(e1.abs());
+
+ MeshPoint e2 = t.vertex1.subtractPosition(t.vertex3);
+ e2 = e2.dividePosition(e2.abs());
+
+ betaCos = e1.dotProduct(e2);
+ } else if (t.vertex3 == centerPoint){
+ MeshPoint e1 = t.vertex3.subtractPosition(t.vertex1);
+ e1 = e1.dividePosition(e1.abs());
+
+ MeshPoint e2 = t.vertex2.subtractPosition(t.vertex1);
+ e2 = e2.dividePosition(e2.abs());
+
+ betaCos = e1.dotProduct(e2);
+ }
+ double beta = Math.acos(betaCos);
+ return 1 / Math.tan(beta);
+ }
+}
diff --git a/Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/TriangleListVisitor.java b/Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/TriangleListVisitor.java
index bd136b7656ecff18eb8a7672ef1a923aa31d81e5..cc929f0cfc3a2e9d01408079ce1066df16f28591 100644
--- a/Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/TriangleListVisitor.java
+++ b/Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/TriangleListVisitor.java
@@ -4,6 +4,7 @@ import cz.fidentis.analyst.mesh.MeshVisitor;
import cz.fidentis.analyst.mesh.core.MeshFacet;
import cz.fidentis.analyst.mesh.core.MeshTriangle;
import java.util.ArrayList;
+import java.util.Collections;
import java.util.List;
/**
@@ -30,13 +31,13 @@ public class TriangleListVisitor extends MeshVisitor {
}
@Override
- protected synchronized void visitMeshFacet(MeshFacet facet) {
- for (MeshTriangle tri : facet) {
- triangles.add(tri);
+ protected void visitMeshFacet(MeshFacet facet) {
+ synchronized (this) {
+ triangles.addAll(facet.getTriangles());
}
}
public List<MeshTriangle> getTriangles() {
- return triangles;
+ return Collections.unmodifiableList(triangles);
}
}
diff --git a/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshTriangle.java b/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshTriangle.java
index 5df388a191fcb8850905f26ae81fc0933fae9f2b..bcd476b777886e70c26514a84a515c20fa034439 100644
--- a/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshTriangle.java
+++ b/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshTriangle.java
@@ -27,6 +27,8 @@ public class MeshTriangle implements Iterable<MeshPoint> {
public final int index2;
public final int index3;
+ private Vector3d voronoiPoint;
+
/**
* Creates new triangle
*
@@ -118,6 +120,10 @@ public class MeshTriangle implements Iterable<MeshPoint> {
* @return the center of circumcircle
*/
public Vector3d getVoronoiPoint() {
+ if (voronoiPoint != null) {
+ return voronoiPoint;
+ }
+
double a = (vertex2.subtractPosition(vertex3)).abs();
double b = (vertex3.subtractPosition(vertex1)).abs();
double c = (vertex2.subtractPosition(vertex1)).abs();
@@ -141,21 +147,24 @@ public class MeshTriangle implements Iterable<MeshPoint> {
double v2Area = ((v3Half.subtractPosition(vertex2)).crossProduct(v1Half.subtractPosition(vertex2))).abs() / 2.0;
double v1Area = (((v1Half.subtractPosition(vertex1)).crossProduct(v3Half.subtractPosition(vertex1))).abs() / 2.0) +
(((v1Half.subtractPosition(vertex1)).crossProduct(v2Half.subtractPosition(vertex1))).abs() / 2.0);
- return new Vector3d(v1Area, v2Area, v3Area);
+ voronoiPoint = new Vector3d(v1Area, v2Area, v3Area);
+ return voronoiPoint;
}
if (d2 < 0) {
double v1Area = ((v3Half.subtractPosition(vertex1)).crossProduct(v2Half.subtractPosition(vertex1))).abs() / 2.0;
double v3Area = ((v1Half.subtractPosition(vertex3)).crossProduct(v2Half.subtractPosition(vertex3))).abs() / 2.0;
double v2Area = (((v2Half.subtractPosition(vertex2)).crossProduct(v1Half.subtractPosition(vertex2))).abs() / 2.0) +
(((v2Half.subtractPosition(vertex2)).crossProduct(v3Half.subtractPosition(vertex2))).abs() / 2.0);
- return new Vector3d(v1Area, v2Area, v3Area);
+ voronoiPoint = new Vector3d(v1Area, v2Area, v3Area);
+ return voronoiPoint;
}
if (d3 < 0) {
double v2Area = ((v1Half.subtractPosition(vertex2)).crossProduct(v3Half.subtractPosition(vertex2))).abs() / 2.0;
double v1Area = ((v2Half.subtractPosition(vertex1)).crossProduct(v3Half.subtractPosition(vertex1))).abs() / 2.0;
double v3Area = (((v3Half.subtractPosition(vertex3)).crossProduct(v2Half.subtractPosition(vertex3))).abs() / 2.0) +
(((v3Half.subtractPosition(vertex3)).crossProduct(v1Half.subtractPosition(vertex3))).abs() / 2.0);
- return new Vector3d(v1Area, v2Area, v3Area);
+ voronoiPoint = new Vector3d(v1Area, v2Area, v3Area);
+ return voronoiPoint;
}
MeshPoint circumcenter = vertex1.multiplyPosition(d1).addPosition(vertex2.multiplyPosition(d2).addPosition(vertex3.multiplyPosition(d3)));
@@ -168,7 +177,9 @@ public class MeshTriangle implements Iterable<MeshPoint> {
double v3Area = (((v1Half.subtractPosition(vertex3)).crossProduct(circumcenter.subtractPosition(vertex3))).abs() / 2.0) +
(((v2Half.subtractPosition(vertex3)).crossProduct(circumcenter.subtractPosition(vertex3))).abs() / 2.0);
- return new Vector3d(v1Area, v2Area, v3Area);
+ voronoiPoint = new Vector3d(v1Area, v2Area, v3Area);
+
+ return voronoiPoint;
}
@Override