Merge branch '174-introduce-a-uniform-space-sampling' into 'master'

Resolve "Introduce a uniform space sampling"

Closes #174

See merge request grp-fidentis/analyst2!187

Comparison/src/main/java/cz/fidentis/analyst/face/HumanFace.java

@@ -13,7 +13,8 @@ import cz.fidentis.analyst.symmetry.Plane;
 import cz.fidentis.analyst.visitors.face.HumanFaceVisitor;
 import cz.fidentis.analyst.visitors.mesh.BoundingBox;
 import cz.fidentis.analyst.visitors.mesh.BoundingBox.BBox;
-import cz.fidentis.analyst.visitors.mesh.Curvature;
+import cz.fidentis.analyst.visitors.mesh.CurvatureCalculator;
+import java.awt.image.BufferedImage;
 
 import java.io.File;
 import java.io.FileInputStream;
@@ -58,8 +59,10 @@ public class HumanFace implements Serializable {
     private final transient EventBus eventBus;
     
     private final String id;
-   
-    private transient Curvature curvature;
+    
+    private transient BufferedImage preview;
+    
+    private boolean hasCurvature = false;
     
     /**
      * Fast (de)serialization handler
@@ -228,15 +231,6 @@ public class HumanFace implements Serializable {
         return symmetryPlane;
     }
     
-    /**
-     * Returns rectangular mesh facet of the symmetry plane, if exists.
-     * @return a rectangular mesh facet of the symmetry plane or {@code null}
-     */
-    //@Deprecated
-    //public MeshRectangleFacet getSymmetryPlaneFacet() {
-    //    return (symmetryPlane == null) ? null : symmetryPlane.getMesh(bbox);
-    //}
-    
     /**
      * Returns {@code true} if the face has the symmetry plane computed.
      * @return {@code true} if the face has the symmetry plane computed.
@@ -253,17 +247,25 @@ public class HumanFace implements Serializable {
         return bbox;
     }
     
+    /** 
+     * Returns {@code true} if the curvature is computed and stored.
+     * 
+     * @return {@code true} if the curvature is computed and stored.
+     */
+    public boolean hasCurvature() {
+        return this.hasCurvature;
+    }
+    
     /**
-     * Computes and returns curvature.
+     * Computes curvature
      * 
-     * @return curvature
+     * @param force If {@code} then the curvature is re-computed even if it already exists.
      */
-    public Curvature getCurvature() {
-        if (this.curvature == null) {
-            this.curvature = new Curvature();
-            getMeshModel().compute(curvature);
+    public void computeCurvature(boolean force) {
+        if (force || !hasCurvature) {
+            meshModel.compute(new CurvatureCalculator());
+            this.hasCurvature = true;
         }
-        return this.curvature;
     }
 
     /**

Comparison/src/main/java/cz/fidentis/analyst/face/HumanFaceUtils.java

@@ -11,6 +11,7 @@ import java.util.zip.DataFormatException;
 import javax.vecmath.Matrix3d;
 import javax.vecmath.Matrix4d;
 import javax.vecmath.Point3d;
+import javax.vecmath.Tuple3d;
 import javax.vecmath.Vector3d;
 
 /**
@@ -109,11 +110,13 @@ public class HumanFaceUtils {
     public static void transformFace(HumanFace face, Vector3d rotation, Vector3d translation, double scale, boolean recomputeKdTree) {
         Quaternion rot = new Quaternion(rotation.x, rotation.y, rotation.z, 1.0);
         
+        // update mesh
         face.getMeshModel().getFacets().stream()
                 .map(facet -> facet.getVertices())
                 .flatMap(meshPoint -> meshPoint.parallelStream())
                 .forEach(meshPoint -> {
                     transformPoint(meshPoint.getPosition(), rot, translation, scale);
+                    transformNormal(meshPoint.getNormal(), rot);
                 });
                     
         
@@ -165,9 +168,15 @@ public class HumanFaceUtils {
         Matrix4d trMat = statPlane.getAlignmentMatrix(tranPlane, preserveUpDir);
         
         // Transform mesh vertices:
+        Matrix3d rotMat = new Matrix3d();
         transformedFace.getMeshModel().getFacets().forEach(f -> {
             f.getVertices().stream().forEach(p -> {
                 trMat.transform(p.getPosition());
+                if (p.getNormal() != null) {
+                    trMat.getRotationScale(rotMat);
+                    rotMat.transform(p.getNormal());
+                    p.getNormal().normalize();
+                }
                 //transformPoint(p.getPosition(), rotation, translation, 1.0);
                 //rotMat.transform(p.getPosition()); // rotate around the plane's normal
             });
@@ -299,7 +308,7 @@ public class HumanFaceUtils {
      * @param translation translation
      * @param scale scale
      */
-    protected static void transformPoint(Point3d point, Quaternion rotation, Vector3d translation, double scale) {
+    protected static void transformPoint(Tuple3d point, Quaternion rotation, Vector3d translation, double scale) {
         Quaternion rotQuat = new Quaternion(point.x, point.y, point.z, 1);
 
         if (rotation != null) {
@@ -313,6 +322,12 @@ public class HumanFaceUtils {
                 rotQuat.z * scale + translation.z
         );
     }
+    
+    protected static void transformNormal(Tuple3d normal, Quaternion rotation) {
+        if (normal != null) {
+            transformPoint(normal, rotation, new Vector3d(0, 0, 0), 1.0); // rotate only
+        }
+    }
 
     /**
      * Transforms the whole plane, i.e., its normal and position.
@@ -325,7 +340,7 @@ public class HumanFaceUtils {
      */
     protected static Plane transformPlane(Plane plane, Quaternion rot, Vector3d translation, double scale) {
         Point3d point = new Point3d(plane.getNormal());
-        transformPoint(point, rot, new Vector3d(0, 0, 0), 1.0); // rotate only
+        transformNormal(point, rot);
         Plane retPlane = new Plane(point, plane.getDistance());
 
         // ... then translate and scale a point projected on the rotate plane:

Comparison/src/main/java/cz/fidentis/analyst/icp/IcpTransformation.java

@@ -94,4 +94,22 @@ public class IcpTransformation {
             );
         }
     }
+    
+    /**
+     * Rotates the normal vector accordingly.
+     *
+     * @param normal Normal vector
+     * @return transformed point or {@code null}
+     */
+    public Vector3d transformNormal(Vector3d normal) {
+        if (normal == null) {
+            return null;
+        }
+        
+        Quaternion rotQuat = new Quaternion(normal.x, normal.y, normal.z, 1);
+        Quaternion rotationCopy = Quaternion.multiply(rotQuat, getRotation().getConjugate());
+        rotQuat = Quaternion.multiply(getRotation(), rotationCopy);
+        
+        return new Vector3d(rotQuat.x, rotQuat.y, rotQuat.z);
+    }
 }

Comparison/src/main/java/cz/fidentis/analyst/icp/IcpTransformer.java

@@ -256,7 +256,7 @@ public class IcpTransformer extends MeshVisitor   {
                 HausdorffDistance.Strategy.POINT_TO_POINT, 
                 false, // relative distance
                 true,  // parallel computation
-                false  // auto cut
+                false  // auto crop
         );
         
         MeshFacet reducedFacet = new UndersampledMeshFacet(transformedFacet, samplingStrategy);
@@ -284,6 +284,9 @@ public class IcpTransformer extends MeshVisitor   {
             transformation = computeIcpTransformation(nearestPoints, distances, reducedFacet.getVertices());
             transformations.get(transformedFacet).add(transformation);
             applyTransformation(transformedFacet, transformation);
+            if (!samplingStrategy.isBackedByOrigMesh()) { // transform samples as well
+                applyTransformation(reducedFacet, transformation);
+            }
             currentIteration ++;
         }
     }
@@ -379,6 +382,9 @@ public class IcpTransformer extends MeshVisitor   {
         transformedFacet.getVertices().parallelStream()
                 .forEach(p -> {
                     p.setPosition(transformation.transformPoint(p.getPosition(), scale));
+                    if (p.getNormal() != null) {
+                        p.setNormal(transformation.transformNormal(p.getNormal()));
+                    }
                 }
         );
     }

Comparison/src/main/java/cz/fidentis/analyst/symmetry/ApproxSymmetryPlane.java

-package cz.fidentis.analyst.symmetry;
-
-import cz.fidentis.analyst.mesh.core.MeshPoint;
-import java.util.List;
-import javax.vecmath.Vector3d;
-
-/**
- * Symmetry plane with votes used for the decision about symmetry estimate of 3D models.
- *
- * @author Natalia Bebjakova
- * 
- */
-public class ApproxSymmetryPlane extends Plane implements Comparable<ApproxSymmetryPlane> {
-    
-    private int votes;
-
-    /**
-     * Constructor.
-     * 
-     * @param vertices Mesh vertices 
-     * @param cache Precomputed values form mesh vertices
-     * @param config Symmetry plane configuration
-     * @param i index of the first significant point used for the plane construction
-     * @param j index of the second significant point used for the plane construction
-     * @param maxDist Distance limit
-     * @throws UnsupportedOperationException if the symmetry plane cannot be created
-     */
-    public ApproxSymmetryPlane(List<MeshPoint> vertices, SymmetryCache cache, SymmetryConfig config, int i, int j, double maxDist) throws UnsupportedOperationException {
-        if (i == j) {
-            throw new UnsupportedOperationException();
-        }
-        setNormAndDist(vertices, cache, config, i, j);
-        computeVotes(vertices, cache, config, maxDist);
-    }
-
-    /**
-     * returns number of votes that were given to plane while computing the symmetry 
-     * 
-     * @return Number of votes 
-     */
-    public int getVotes() {
-        return votes;
-    }
-    
-    private void setNormAndDist(List<MeshPoint> vertices, SymmetryCache cache, SymmetryConfig config, int i, int j) {
-        MeshPoint meshPointI = vertices.get(i);
-        MeshPoint meshPointJ = vertices.get(j);
-                
-        // accpet only point pairs with significantly different curvatures (WTF?)
-        double minRatio = config.getMinCurvRatio();
-        double maxRatio = 1.0 / minRatio;
-        double ratioIJ = cache.getCurRatio(i, j);
-        if (Double.isFinite(ratioIJ) && (ratioIJ < minRatio || ratioIJ > maxRatio)) {
-            throw new UnsupportedOperationException();
-        }
-                
-        Vector3d p1 = new Vector3d(meshPointI.getPosition());
-        Vector3d p2 = new Vector3d(meshPointJ.getPosition());
-        Vector3d normal = new Vector3d(p1);
-        normal.sub(p2);
-        normal.normalize(); 
-        
-        // accpect only point pair with oposite normals along with the plane normal:
-        double normCos = cache.getNormCosVec(i, j).dot(normal);
-        if (Math.abs(normCos) < config.getMinNormAngleCos()) {
-            throw new UnsupportedOperationException();
-        }
-        
-        setDistance(-normal.dot(cache.getAvgPos(i, j))); 
-        setNormal(normal);
-    }
-
-    /**
-     * Computes votes for given plane 
-     *
-     * @param sigPoints Mesh vertices with the most significant curvature
-     * @param config Symmetry plane configuration
-     * @param maxDist Distance limit
-     */
-    private void computeVotes(List<MeshPoint> vertices, SymmetryCache cache, SymmetryConfig config, double maxDist) {
-        normalize();
-        Vector3d normal = getNormal();
-        double d = getDistance();
-        double maxCurvRatio = 1.0 / config.getMinCurvRatio();
-        
-        for (int i = 0; i < vertices.size(); i++) {
-            for (int j = 0; j < vertices.size(); j++) {
-                if (i == j) {
-                    continue;
-                }
-                
-                double ratioIJ = cache.getCurRatio(i, j);
-                if (Double.isFinite(ratioIJ) && (ratioIJ < config.getMinCurvRatio() || ratioIJ > maxCurvRatio)) {
-                    continue;
-                }
-                
-                double normCos = cache.getNormCosVec(i, j).dot(normal);
-                if (Math.abs(normCos) < config.getMinNormAngleCos()) {
-                    continue;
-                }
-                
-                // Caching this part doesn't improve efficiency anymore:
-                Vector3d vec = new Vector3d(vertices.get(i).getPosition());
-                vec.sub(vertices.get(j).getPosition());
-                vec.normalize();
-                double cos = vec.dot(normal);
-                if (Math.abs(cos) < config.getMinAngleCos()) {
-                    continue;
-                }
-                
-                Vector3d avg = cache.getAvgPos(i, j);
-                double dist = Math.abs(normal.dot(avg) + d);
-                if (dist <= maxDist) {                 
-                    votes++;  
-                }   
-            }
-        }
-    }
-    
-    /**
-     * Enables to compare two approximate planes due to number of votes 
-     * 
-     * @param other plane to be compared 
-     * @return number that decides which plane has more votes 
-     */
-    @Override
-    public int compareTo(ApproxSymmetryPlane other) {
-        return Integer.compare(votes, other.votes);
-    }
-    
-    @Override
-    public String toString() {
-        return this.getNormal() + " " + getDistance() + " " + votes;
-    }
-}

Comparison/src/main/java/cz/fidentis/analyst/symmetry/Plane.java

@@ -4,6 +4,7 @@ import cz.fidentis.analyst.mesh.core.MeshRectangleFacet;
 import cz.fidentis.analyst.visitors.mesh.BoundingBox.BBox;
 import java.io.Serializable;
 import java.util.List;
+import java.util.Objects;
 import javax.vecmath.Matrix4d;
 import javax.vecmath.Point3d;
 import javax.vecmath.Tuple3d;
@@ -42,7 +43,7 @@ public class Plane implements Serializable {
     }
     
     /**
-     * Constructor from a 3D point.
+     * Construction from a 3D point.
      * 
      * @param point point in space
      * @throws IllegalArgumentException if the @code{plane} argument is null
@@ -53,6 +54,24 @@ public class Plane implements Serializable {
         this.normal.normalize();
     }
     
+    /**
+     * Symmetry plane constructed from two points.
+     * 
+     * @param point1 point in space
+     * @param point2 point in space
+     * @throws IllegalArgumentException if the @code{plane} argument is null
+     */
+    public Plane(Tuple3d point1, Tuple3d point2) {
+        this.normal = new Vector3d(point1);
+        this.normal.sub(point2);
+        this.normal.normalize();
+        this.distance = this.normal.dot(new Vector3d(
+                (point1.x + point2.x) / 2.0,
+                (point1.y + point2.y) / 2.0,
+                (point1.z + point2.z) / 2.0)
+        );
+    }
+    
     /**
      * Creates average plane from existing planes.
      * 
@@ -85,13 +104,42 @@ public class Plane implements Serializable {
     
     protected Plane() {
     }
+
+    @Override
+    public int hashCode() {
+        int hash = 5;
+        hash = 97 * hash + Objects.hashCode(this.normal);
+        hash = 97 * hash + (int) (Double.doubleToLongBits(this.distance) ^ (Double.doubleToLongBits(this.distance) >>> 32));
+        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 Plane other = (Plane) obj;
+        if (Double.doubleToLongBits(this.distance) != Double.doubleToLongBits(other.distance)) {
+            return false;
+        }
+        if (!Objects.equals(this.normal, other.normal)) {
+            return false;
+        }
+        return true;
+    }
     
     /**
      * Normalize the plane
      */
     protected final void normalize() {
         double normalLength = normal.length();
-        normal.normalize();
+        normal.scale(1.0 / normalLength);
         distance /= normalLength; // Do we really want this? --ro
     }
     
@@ -105,6 +153,10 @@ public class Plane implements Serializable {
         return normal + " dist " + distance;
     }
     
+    /**
+     * Return a copy of the normal vector
+     * @return a copy of the normal vector
+     */
     public Vector3d getNormal() {
         return new Vector3d(normal);
     }
@@ -310,24 +362,30 @@ public class Plane implements Serializable {
      * @return a point on the opposite side of the plane.
      * @throws NullPointerException if the {@code point} is {@code null}
      */
-    public Point3d reflectOverPlane(Point3d point) {
+    public Point3d reflectPointOverPlane(Point3d point) {
         Point3d ret = new Point3d(normal);
         ret.scale(-2.0 * getPointDistance(point));
         ret.add(point);
         return ret;
+    }
+
+    /**
+     * Reflects the give unit vector over this plane.
+     * 
+     * @param vector A normalized 3D vector
+     * @return reflected vector
+     * @throws NullPointerException if the {@code point} is {@code null}
+     */
+    public Vector3d reflectUnitVectorOverPlane(Vector3d vector) {
+        double s = 2 * this.getNormal().dot(vector);
+        s /= this.getNormal().dot(this.getNormal());
         
-        /*
-        double shiftDist = 
-                ((normal.x * point.x) + 
-                  (normal.y * point.y) + 
-                  (normal.z * point.z) +
-                  distance) / normSquare;
-        //System.out.println("HHH "+ shiftDist);
-        Point3d ret = new Point3d(normal);
-        ret.scale(-2.0 * shiftDist);
-        ret.add(point);
-        return ret;
-        */
+        Vector3d v = new Vector3d(this.getNormal());
+        v.scale(s);
+        
+        Vector3d rv = new Vector3d(vector);
+        rv.sub(v);
+        return rv;
     }
     
     /**
@@ -368,6 +426,14 @@ public class Plane implements Serializable {
     public MeshRectangleFacet getMesh(BBox bbox) {
         return getMesh(bbox.getMidPoint(), bbox.getDiagonalLength(), bbox.getDiagonalLength());
     }
+    
+    /**
+     * Returns reference to the normal vector
+     * @return reference to the normal vector
+     */
+    protected Vector3d getNormalReference() {
+        return normal;
+    }
 
     /**
      * Changes the normal vector.

Comparison/src/main/java/cz/fidentis/analyst/symmetry/SymmetryCache.java

-package cz.fidentis.analyst.symmetry;
-
-import cz.fidentis.analyst.mesh.core.MeshPoint;
-import java.util.ArrayList;
-import java.util.List;
-import javax.vecmath.Vector3d;
-
-/**
- * Precomputed values for the symmetry plane estimation.
- * 
- * @author Radek Oslejsek
- */
-public class SymmetryCache {
-    
-    private List<List<Vector3d>> normCosVecCache;
-    private List<List<Vector3d>> avgPosCache;
-    private List<List<Double>> curRatioCache;
-    
-    /**
-     * Constructor.
-     * 
-     * @param vertices Mesh vertices 
-     * @param curvatures Curvatures in the vertices
-     */
-    public SymmetryCache(List<MeshPoint> vertices, List<Double> curvatures) {
-        if (vertices == null || vertices.isEmpty()) {
-            throw new IllegalArgumentException("points");
-        }
-        
-        normCosVecCache = new ArrayList<>(vertices.size());
-        avgPosCache = new ArrayList<>(vertices.size());
-        curRatioCache = (curvatures == null) ? null : new ArrayList<>(vertices.size());
-        
-        for (int i = 0; i < vertices.size(); i++) {
-            List<Vector3d> cosArray = new ArrayList<>(vertices.size());
-            normCosVecCache.add(cosArray);
-            
-            List<Vector3d> posArray = new ArrayList<>(vertices.size());
-            avgPosCache.add(posArray);
-            
-            List<Double> curArray = null;
-            if (curvatures != null) {
-                curArray = new ArrayList<>(vertices.size());
-                curRatioCache.add(curArray);
-            }
-            
-            for (int j = 0; j < vertices.size(); j++) {
-                MeshPoint meshPointI = vertices.get(i);
-                MeshPoint meshPointJ = vertices.get(j);
-                
-                Vector3d ni = new Vector3d(meshPointI.getNormal());
-                Vector3d nj = new Vector3d(meshPointJ.getNormal());
-                ni.normalize();
-                nj.normalize();
-                ni.sub(nj);
-                ni.normalize();
-                cosArray.add(ni);
-                
-                Vector3d avrg = new Vector3d(meshPointI.getPosition());
-                Vector3d aux = new Vector3d(meshPointJ.getPosition());
-                avrg.add(aux);
-                avrg.scale(0.5);
-                posArray.add(avrg);
-                if (curvatures != null) {
-                    curArray.add(curvatures.get(i) / curvatures.get(j));
-                }
-            }
-        }
-    }
-    
-    /**
-     * Returns cached normal vector for cos.
-     * @param i index of the i-th significant curvature point
-     * @param j index of the j-th significant curvature point
-     * @return cached value or null
-     */
-    public Vector3d getNormCosVec(int i, int j) {
-        if (normCosVecCache == null || i >= normCosVecCache.size() || i >= normCosVecCache.size() || i < 0 || j < 0) {
-            return null;
-        }
-        return normCosVecCache.get(i).get(j);
-    }
-    
-    /**
-     * Returns cached vector for the computation average normal.
-     * @param i index of the i-th significant curvature point
-     * @param j index of the j-th significant curvature point
-     * @return cached value or null
-     */
-    public Vector3d getAvgPos(int i, int j) {
-        if (avgPosCache == null || i >= avgPosCache.size() || i >= avgPosCache.size() || i < 0 || j < 0) {
-            return null;
-        }
-        return avgPosCache.get(i).get(j);
-    }
-    
-    /**
-     * Returns cached curvature ratio.
-     * @param i index of the i-th significant curvature point
-     * @param j index of the j-th significant curvature point
-     * @return cached value or {@code Double.NaN}
-     */
-    public double getCurRatio(int i, int j) {
-        if (curRatioCache == null || i >= curRatioCache.size() || i >= curRatioCache.size() || i < 0 || j < 0) {
-            return Double.NaN;
-        }
-        return curRatioCache.get(i).get(j);
-    }
-    
-}

Comparison/src/main/java/cz/fidentis/analyst/symmetry/SymmetryConfig.java

-package cz.fidentis.analyst.symmetry;
-
-import cz.fidentis.analyst.visitors.mesh.sampling.CurvatureSampling;
-
-/**
- * Representation of configuration for symmetry estimate.
- * Default numbers are given due to the best results on tested data.
- * On many different 3D models, it exists other values of config that will have 
- * better impact on results in estimate of symmetry. 
- * 
- * @author Natalia Bebjakova
- */
-public class SymmetryConfig {
-    
-    private static final double DEFAULT_MIN_CURV_RATIO = 0.5;
-    private static final double DEFAULT_MIN_ANGLE_COS = 0.985;
-    private static final double DEFAULT_MIN_NORM_ANGLE_COS = 0.985;
-    private static final double DEFAULT_MAX_REL_DISTANCE = 1.0 / 100.0;
-    private static final int DEFAULT_SIGNIFICANT_POINT_COUNT = 200;
-    private static final boolean DEFAULT_AVERAGING = true;
-    private static final CurvatureSampling.CurvatureAlg DEFAULT_CURVATURE_ALGORITHM = CurvatureSampling.CurvatureAlg.GAUSSIAN;
-    
-    private double minCurvRatio;
-    private double minAngleCos;
-    private double minNormAngleCos;
-    private double maxRelDistance;
-    private int significantPointCount;
-    private boolean averaging;
-    private CurvatureSampling.CurvatureAlg curvatureAlg;
-    
-    /**
-     * Creates configuration with default values 
-     */
-    public SymmetryConfig() {
-        minCurvRatio = DEFAULT_MIN_CURV_RATIO;
-        minAngleCos = DEFAULT_MIN_ANGLE_COS;
-        minNormAngleCos = DEFAULT_MIN_NORM_ANGLE_COS;
-        maxRelDistance = DEFAULT_MAX_REL_DISTANCE;
-        significantPointCount = DEFAULT_SIGNIFICANT_POINT_COUNT;
-        averaging = DEFAULT_AVERAGING;
-        curvatureAlg = DEFAULT_CURVATURE_ALGORITHM;
-    }
-    
-    /**
-     * Copy constructor.
-     * 
-     * @param conf Original configuration
-     */
-    public SymmetryConfig(SymmetryConfig conf) {
-        copy(conf);
-    }
-    
-    /**
-     * Copies values from another configuration object.
-     * 
-     * @param conf Another configuration
-     */
-    public void copy(SymmetryConfig conf) {
-        minCurvRatio = conf.minCurvRatio;
-        minAngleCos = conf.minAngleCos;
-        minNormAngleCos = conf.minNormAngleCos;
-        maxRelDistance = conf.maxRelDistance;
-        significantPointCount = conf.significantPointCount;
-        averaging = conf.averaging;
-        curvatureAlg = conf.curvatureAlg;
-    }
-    
-    /**
-     * Parameter which denotes how similar the Gaussian curvatures in the two vertices
-     * must be for this criteria to be satisfied.
-     * The higher the value is the more similar they must be.
-     * 
-     * @return minimal similarity of curvatures to satisfy the criteria  
-     */
-    public double getMinCurvRatio() {
-        return minCurvRatio;
-    }
-
-    /**
-     * 
-     * @param minCurvRatio new minimal similarity of curvatures to satisfy the criteria  
-     */
-    public void setMinCurvRatio(double minCurvRatio) {
-        this.minCurvRatio = minCurvRatio;
-    }
-
-    /**
-     * MinAngleCos ∈ (0,1)
-     * It is the angle between the vector (xk − xl) and the normal vector nij of the plane ρij
-     * (which is the vector (xi −xj))
-     * Returns parameter which denotes how large the angle αij can be for this criteria to be satisfied.
-     * 
-     * @return minimal angle satisfy criteria 
-     */
-    public double getMinAngleCos() {
-        return minAngleCos;
-    }
-
-    /**
-     * 
-     * @param minAngleCos new minimal angle to satisfy the criteria 
-     */
-    public void setMinAngleCos(double minAngleCos) {
-        this.minAngleCos = minAngleCos;
-    }
-
-    /**
-     * MinNormAngleCos ∈ (0,1)
-     * It is angle between vectors (xi − xj) and (ni − nj), where xi, xj are vertices and ni, nj its normals
-     * Returns parameter which denotes how large the angle αij can be for this criteria to be satisfied.
-     * 
-     * @return minimal angle to satisfy criteria 
-     */
-    public double getMinNormAngleCos() {
-        return minNormAngleCos;
-    }
-
-    /**
-     * 
-     * @param minNormAngleCos new minimal angle to satisfy the criteria 
-     */
-    public void setMinNormAngleCos(double minNormAngleCos) {
-        this.minNormAngleCos = minNormAngleCos;
-    }
-
-    /**
-     * Parameter which denotes how far (relatively to the length of the bounding box diagonal)
-     * the middle point of the two vertices can be from the plane in order to satisfy this criteria.
-     * 
-     * @return relative distance
-     */
-    public double getMaxRelDistance() {
-        return maxRelDistance;
-    }
-
-    /**
-     * 
-     * @param maxRelDistance new relative distance 
-     */
-    public void setMaxRelDistance(double maxRelDistance) {
-        this.maxRelDistance = maxRelDistance;
-    }
-
-    /**
-     * Returns number of vertices with the highest Gaussian curvature.
-     * 
-     * @return number of significant points for computing the symmetry 
-     */
-    public int getSignificantPointCount() {
-        return significantPointCount;
-    }
-
-    /**
-     * 
-     * @param significantPointCount new number of significant points for computing the symmetry 
-     */
-    public void setSignificantPointCount(int significantPointCount) {
-        this.significantPointCount = significantPointCount;
-    }
-
-    /**
-     * If there are more planes with the same highest number of votes while computing symmetry,
-     * we can average them all together. 
-     * Returns parameter that decides whether to average them or not.
-     * 
-     * @return true if planes will be averaged 
-     */
-    public boolean isAveraging() {
-        return averaging;
-    }
-
-    /**
-     * 
-     * @param averaging new averaging flag  
-     */
-    public void setAveraging(boolean averaging) {
-        this.averaging = averaging;
-    }
-    
-    /**
-     * Returns curvature algorithm.
-     * @return curvature algorithm 
-     */
-    public CurvatureSampling.CurvatureAlg getCurvatureAlg() {
-        return this.curvatureAlg;
-    }
-    
-    /**
-     * Sets  curvature algorithm.
-     * @param alg curvature algorithm
-     */
-    public void setCurvatureAlg(CurvatureSampling.CurvatureAlg alg) {
-        this.curvatureAlg = alg;
-    }
-
-    /**
-     * 
-     * @return String representation of configuration
-     */
-    @Override 
-    public String toString() {
-        String str = "PARAMETERS: ";
-        str += "\n";
-        str += "Min curvature ratio: " + minCurvRatio + "\n";
-        str += "Min angle cosine: " + minAngleCos + "\n";
-        str += "Min norm angle cosine: " + minNormAngleCos + "\n";
-        str += "Max relative distance: " + maxRelDistance + "\n";
-        str += "Significant points: " + significantPointCount + "\n";
-        str += "Averaging: " + averaging + "\n";
-        str += "Curvature: " + curvatureAlg + "\n";
-        return str;
-    }
-}
\ No newline at end of file

Comparison/src/main/java/cz/fidentis/analyst/symmetry/SymmetryEstimator.java

 package cz.fidentis.analyst.symmetry;
 
 import cz.fidentis.analyst.mesh.MeshVisitor;
-import cz.fidentis.analyst.mesh.core.MeshFacet;
-import cz.fidentis.analyst.mesh.core.MeshPoint;
-import cz.fidentis.analyst.visitors.mesh.BoundingBox;
-import cz.fidentis.analyst.visitors.mesh.BoundingBox.BBox;
-import cz.fidentis.analyst.visitors.mesh.sampling.CurvatureSampling;
-import cz.fidentis.analyst.visitors.mesh.sampling.PointSampling;
-import java.util.ArrayList;
-import java.util.List;
-import java.util.concurrent.Callable;
-import java.util.concurrent.ExecutionException;
-import java.util.concurrent.ExecutorService;
-import java.util.concurrent.Executors;
-import java.util.concurrent.Future;
-import java.util.logging.Level;
 
 /**
- * Main class for computing approximate plane of symmetry of the 3D model.
- * Default values of the configuration are given due to the best results on tested objects.
- * On many different 3D models, it exists other values of config that will have better impact on result. 
- * <p>
- * This visitor <b>is not thread-safe</b>, i.e., a single instance of the visitor 
- * cannot be used to inspect multiple meshes simultaneously  (sequential inspection is okay). 
- * It because the underlying {@link SignificantPoints} visitor is not thread-safe.
- * </p>
- * <p>
- * The main symmetry plane computation is performed by the {@link #getSymmetryPlane()} method,
- * not the {@link #visitMeshFacet(cz.fidentis.analyst.mesh.core.MeshFacet)}.
- * </p>
+ * Symmetry estimator.
  * 
- * @author Natalia Bebjakova
  * @author Radek Oslejsek
  */
-public class SymmetryEstimator extends MeshVisitor {
-    
-    private final SymmetryConfig config;
-    private final PointSampling samplingStrategy;
-    private final BoundingBox bbVisitor = new BoundingBox();
-    
-    private Plane symmetryPlane;
-    
+public abstract class SymmetryEstimator extends MeshVisitor {
+   
     /**
-     * Constructor.
-     * 
-     * @param config Algorithm options
-     * @param samplingStrategy Downsampling strategy. Must not be {@code null}
-     * @throws IllegalArgumentException if some input parameter is missing
+     * Returns a symmetry plane.
+     * @return a symmetry plane or {@code null{
      */
-    public SymmetryEstimator(SymmetryConfig config, PointSampling samplingStrategy) {
-        if (config == null) {
-            throw new IllegalArgumentException("config");
-        }
-        if (samplingStrategy == null) {
-            throw new IllegalArgumentException("samplingStrategy");
-        }
-        this.config = config;
-        this.samplingStrategy = samplingStrategy;
-    }
-    
-    @Override
-    public boolean isThreadSafe() {
-        return false;
-    }
-    
-    @Override
-    public void visitMeshFacet(MeshFacet facet) {
-        // We need vertex normals for the plane computation
-        synchronized (this) {
-            if (!facet.hasVertexNormals()) {
-                facet.calculateVertexNormals();
-            }
-        }
-        samplingStrategy.visitMeshFacet(facet);
-        bbVisitor.visitMeshFacet(facet);
-    }
-    
-    /**
-     * Computes the symmetry plane concurrently.The plane is computed only once, then the same instance is returned.
-     * 
-     * @return Symmetry plane
-     */
-    public Plane getSymmetryPlane() {
-        return getSymmetryPlane(true);
-    }
-    
-    /**
-     * Computes and returns the symmetry plane. The plane is computed only once, 
-     * then the same instance is returned.
-     * 
-     * @param concurrently If {@code true}, then parallel computation is used utilizing all CPU cores.
-     * @return the symmetry plane or {@code null}
-     */
-    public Plane getSymmetryPlane(boolean concurrently) {
-        if (symmetryPlane == null) {
-            this.calculateSymmetryPlane(concurrently);
-        }
-        return symmetryPlane;
-    }
-    
-    /**
-     * Returns the bounding box computed during the {@link SymmetryEstimator#calculateSymmetryPlane}.
-     * @return the bounding box or {@code null}
-     */
-    public BBox getBoundingBox() {
-        return bbVisitor.getBoundingBox();
-    }
-    
-    /**
-     * Calculates the symmetry plane.
-     * @param concurrently If {@code true}, then parallel computation is used utilizing all CPU cores.
-     * Otherwise, the computation is sequential.
-     */
-    protected void calculateSymmetryPlane(boolean concurrently) {
-        final List<Plane> planes = new ArrayList<>();
-        final double maxDistance = bbVisitor.getBoundingBox().getDiagonalLength() * config.getMaxRelDistance();
-        
-        List<MeshPoint> sigVertices = samplingStrategy.getSamples();
-        SymmetryCache cache = new SymmetryCache(
-                sigVertices, 
-                (samplingStrategy.getClass() == CurvatureSampling.class) 
-                        ? ((CurvatureSampling)samplingStrategy).getSampledCurvatures()
-                        : null
-        );
-        
-        /*
-         * Sequential computation
-         */
-        if (!concurrently) {
-            int maxVotes = 0;
-            for (int i = 0; i < sigVertices.size(); i++) {
-                for (int j = 0; j < sigVertices.size(); j++) {
-                    ApproxSymmetryPlane newPlane;
-                    try {
-                        newPlane = new ApproxSymmetryPlane(sigVertices, cache, config, i, j, maxDistance);
-                    } catch(UnsupportedOperationException ex) {
-                        continue;
-                    }
-                    maxVotes = checkAndUpdatePlanes(planes, newPlane, maxVotes);
-                }
-            }
-            setSymmetryPlane(planes);
-            return;
-        }
-        
-        /*
-         * Concurrent computation
-         */
-        
-        // Initiate structure for concurrent computation:
-        ExecutorService executor = Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors());
-        final List<Future<ApproxSymmetryPlane>> results = new ArrayList<>(sigVertices.size() * sigVertices.size());
-        
-        // Compute candidate planes concurrently:
-        for (int i = 0; i < sigVertices.size(); i++) {
-            for (int j = 0; j < sigVertices.size(); j++) {
-                int finalI = i;
-                int finalJ = j;
-                
-                Callable<ApproxSymmetryPlane> worker = new  Callable<ApproxSymmetryPlane>() {
-                    @Override
-                    public ApproxSymmetryPlane call() throws Exception {
-                        try {
-                            return new ApproxSymmetryPlane(sigVertices, cache, config, finalI, finalJ, maxDistance);
-                        } catch(UnsupportedOperationException ex) {
-                            return null;
-                        }
-                    }
-                };
-                
-                results.add(executor.submit(worker));
-            }
-        }
-        
-        // Wait until all symmetry planes are computed:
-        executor.shutdown();
-        while (!executor.isTerminated()){}
-        
-        // Get the best-fitting planes:
-        try {
-            int maxVotes = 0;
-            for (Future<ApproxSymmetryPlane> res: results) {
-                ApproxSymmetryPlane newPlane = res.get();
-                if (newPlane != null) {
-                    maxVotes = checkAndUpdatePlanes(planes, newPlane, maxVotes);
-                }
-            }
-        } catch (final InterruptedException | ExecutionException ex) {
-            java.util.logging.Logger.getLogger(SymmetryEstimator.class.getName()).log(Level.SEVERE, null, ex);
-        }
-        
-        setSymmetryPlane(planes);
-    }
-    
-    protected int checkAndUpdatePlanes(List<Plane> planes, ApproxSymmetryPlane newPlane, int maxVotes) {
-        if (newPlane.getVotes() > maxVotes) {
-            planes.clear();
-            planes.add(newPlane);
-            return newPlane.getVotes();
-        } else if (newPlane.getVotes() == maxVotes) {
-            planes.add(newPlane);
-        }
-        return maxVotes;
-    }
-
-    protected void setSymmetryPlane(List<Plane> planes) {
-        if (config.isAveraging() && !planes.isEmpty()) {
-            symmetryPlane = new Plane(planes);
-        } else {
-            symmetryPlane = planes.isEmpty() ? null : planes.get(0);
-        }
-    }
+    public abstract Plane getSymmetryPlane();
 }

Comparison/src/main/java/cz/fidentis/analyst/symmetry/SymmetryEstimatorMesh.java

+package cz.fidentis.analyst.symmetry;
+
+import cz.fidentis.analyst.mesh.core.Curvature;
+import cz.fidentis.analyst.mesh.core.MeshFacet;
+import cz.fidentis.analyst.mesh.core.MeshPoint;
+import cz.fidentis.analyst.visitors.mesh.BoundingBox;
+import cz.fidentis.analyst.visitors.mesh.CurvatureCalculator;
+import cz.fidentis.analyst.visitors.mesh.sampling.PointSampling;
+import java.util.ArrayList;
+import java.util.List;
+import java.util.concurrent.Callable;
+import java.util.concurrent.ExecutionException;
+import java.util.concurrent.ExecutorService;
+import java.util.concurrent.Executors;
+import java.util.concurrent.Future;
+import java.util.logging.Level;
+import javax.vecmath.Vector3d;
+
+/**
+ * A old implementation of the symmetry plane estimator taken from the old FIDENTS.
+ * Conceptually, it is similar to the {@link SymmetryEstimatorRobustMesh}.
+ * It also uses Gaussian curvature and normal vectors to prune candidate planes 
+ * and measure their quality. But some additional criteria are used as well. 
+ * Also, the calculation of "votes" differs.
+ * This implementation has the following properties:
+ * <ul>
+ * <li>Fast. A bit faster than the new robust algorithms of {@link SymmetryEstimatorRobust}
+ * and much faster that the {@link SymmetryEstimatorRobustMesh}.</li>
+ * <li>Best candidate planes are chosen based on the similarity of Gauss curvatures 
+ * on both sides of the plane, inverse direction of normal vectors, 
+ * and the distance of plane from centroid.</li>
+ * <li>Because the space of candidate planes is roughly approximated (several true-false rules), 
+ * the final symmetry plane is computed by averaging best candidates (candidates with the most votes).</li>
+ * <li>Best results are achieved with Uniform Grid sampling with 200-300 points and averaging turn on.</li>
+ * <li>Should return reasonable results also for incomplete faces (was not tested).</li>
+ * </ul>
+ * <p>
+ * This visitor <b>is not thread-safe</b>, i.e., a single instance of the visitor 
+ * cannot be used to inspect multiple meshes simultaneously  (sequential inspection is okay). 
+ * It because the underlying visitors may not be thread-safe.
+ * </p>
+ * <p>
+ * The main symmetry plane computation is performed by the {@link #getSymmetryPlane()} method,
+ * not the {@link #visitMeshFacet(cz.fidentis.analyst.mesh.core.MeshFacet)}.
+ * </p>
+ * 
+ * @author Natalia Bebjakova
+ * @author Radek Oslejsek
+ */
+public class SymmetryEstimatorMesh extends SymmetryEstimator {
+    
+    public static final double MAX_REL_DISTANCE = 0.01;
+    
+    /**
+     * If {@code true}, then the final symmetry plane is computed by averaging the best candidates.
+     * If {@code false}, then the random top candidate is selected.
+     */
+    public static final boolean AVERAGING = true;
+    
+    private final PointSampling samplingStrategy;
+    private final int samplingLimit;
+    private final BoundingBox bbVisitor = new BoundingBox();
+    
+    private Plane symmetryPlane;
+    
+    /**
+     * Constructor.
+     * 
+     * @param samplingStrategy Downsampling strategy. Must not be {@code null}
+     * @param samplingLimit Desired number of samples for udenrsampling
+     * @throws IllegalArgumentException if some input parameter is missing
+     */
+    public SymmetryEstimatorMesh(PointSampling samplingStrategy, int samplingLimit) {
+        if (samplingStrategy == null) {
+            throw new IllegalArgumentException("samplingStrategy");
+        }
+        this.samplingStrategy = samplingStrategy;
+        this.samplingLimit = samplingLimit;
+    }
+    
+    @Override
+    public boolean isThreadSafe() {
+        return false;
+    }
+    
+    @Override
+    public void visitMeshFacet(MeshFacet facet) {
+        // We need vertex normals for the plane computation
+        synchronized (this) {
+            if (!facet.hasVertexNormals()) {
+                facet.calculateVertexNormals();
+            }
+        }
+        
+        // We need Gaussian curvatures
+        synchronized (this) {
+            if (facet.getVertex(0).getCurvature() == null) {
+                facet.accept(new CurvatureCalculator());
+            }
+        }
+        
+        // We want to downsample the mesh
+        samplingStrategy.visitMeshFacet(facet);
+        
+        // We need bounding box to compute max distance
+        bbVisitor.visitMeshFacet(facet);
+    }
+    
+    /**
+     * Computes the symmetry plane concurrently.The plane is computed only once, then the same instance is returned.
+     * 
+     * @return Symmetry plane
+     */
+    public Plane getSymmetryPlane() {
+        return getSymmetryPlane(true);
+    }
+    
+    /**
+     * Computes and returns the symmetry plane. The plane is computed only once, 
+     * then the same instance is returned.
+     * 
+     * @param concurrently If {@code true}, then parallel computation is used utilizing all CPU cores.
+     * @return the symmetry plane or {@code null}
+     */
+    public Plane getSymmetryPlane(boolean concurrently) {
+        if (symmetryPlane == null) {
+            this.calculateSymmetryPlane(concurrently);
+        }
+        return symmetryPlane;
+    }
+    
+    /**
+     * Calculates the symmetry plane.
+     * @param concurrently If {@code true}, then parallel computation is used utilizing all CPU cores.
+     * Otherwise, the computation is sequential.
+     */
+    protected void calculateSymmetryPlane(boolean concurrently) {
+        final List<Plane> planes = new ArrayList<>();
+        final double maxDistance = bbVisitor.getBoundingBox().getDiagonalLength() * MAX_REL_DISTANCE;
+        
+        samplingStrategy.setRequiredSamples(samplingLimit);
+        List<MeshPoint> sigVertices = samplingStrategy.getSamples();
+        SymmetryCache cache = new SymmetryCache(sigVertices);
+        
+        /*
+         * Sequential computation
+         */
+        if (!concurrently) {
+            int maxVotes = 0;
+            for (int i = 0; i < sigVertices.size(); i++) {
+                for (int j = 0; j < sigVertices.size(); j++) {
+                    CandidatePlane newPlane;
+                    try {
+                        newPlane = new CandidatePlane(sigVertices, cache, i, j, maxDistance);
+                    } catch(UnsupportedOperationException ex) {
+                        continue;
+                    }
+                    maxVotes = checkAndUpdatePlanes(planes, newPlane, maxVotes);
+                }
+            }
+            setSymmetryPlane(planes);
+            return;
+        }
+        
+        /*
+         * Concurrent computation
+         */
+        
+        // Initiate structure for concurrent computation:
+        ExecutorService executor = Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors());
+        final List<Future<CandidatePlane>> results = new ArrayList<>(sigVertices.size() * sigVertices.size());
+        
+        // Compute candidate planes concurrently:
+        for (int i = 0; i < sigVertices.size(); i++) {
+            for (int j = 0; j < sigVertices.size(); j++) {
+                int finalI = i;
+                int finalJ = j;
+                
+                Callable<CandidatePlane> worker = new  Callable<CandidatePlane>() {
+                    @Override
+                    public CandidatePlane call() throws Exception {
+                        try {
+                            return new CandidatePlane(sigVertices, cache, finalI, finalJ, maxDistance);
+                        } catch(UnsupportedOperationException ex) {
+                            return null;
+                        }
+                    }
+                };
+                
+                results.add(executor.submit(worker));
+            }
+        }
+        
+        // Wait until all symmetry planes are computed:
+        executor.shutdown();
+        while (!executor.isTerminated()){}
+        
+        // Get the best-fitting planes:
+        try {
+            int maxVotes = 0;
+            for (Future<CandidatePlane> res: results) {
+                CandidatePlane newPlane = res.get();
+                if (newPlane != null) {
+                    maxVotes = checkAndUpdatePlanes(planes, newPlane, maxVotes);
+                }
+            }
+        } catch (final InterruptedException | ExecutionException ex) {
+            java.util.logging.Logger.getLogger(SymmetryEstimatorMesh.class.getName()).log(Level.SEVERE, null, ex);
+        }
+        
+        setSymmetryPlane(planes);
+    }
+    
+    protected int checkAndUpdatePlanes(List<Plane> planes, CandidatePlane newPlane, int maxVotes) {
+        if (newPlane.votes > maxVotes) {
+            planes.clear();
+            planes.add(newPlane);
+            return newPlane.votes;
+        } else if (newPlane.votes == maxVotes) {
+            planes.add(newPlane);
+        }
+        return maxVotes;
+    }
+
+    protected void setSymmetryPlane(List<Plane> planes) {
+        if (AVERAGING && !planes.isEmpty()) {
+            symmetryPlane = new Plane(planes);
+        } else {
+            symmetryPlane = planes.isEmpty() ? null : planes.get(0);
+        }
+    }
+    
+    /**
+     * Symmetry plane with votes used for the decision about symmetry estimate
+     * of 3D models.
+     *
+     * @author Natalia Bebjakova
+     *
+     */
+    protected class CandidatePlane extends Plane implements Comparable<CandidatePlane> {
+
+        public static final double MIN_ANGLE_COS = 0.985;
+        public static final double AVG_CURVATURE_MULTIPLICATOR = 0.01;
+
+        private int votes;
+
+        /**
+         * Constructor.
+         *
+         * @param vertices Mesh vertices
+         * @param cache Precomputed values form mesh vertices
+         * @param i index of the first significant point used for the plane
+         * construction
+         * @param j index of the second significant point used for the plane
+         * construction
+         * @param maxDist Distance limit
+         * @throws UnsupportedOperationException if the symmetry plane cannot be
+         * created
+         */
+        protected CandidatePlane(List<MeshPoint> vertices, SymmetryCache cache, int i, int j, double maxDist) throws UnsupportedOperationException {
+            if (i == j) {
+                throw new UnsupportedOperationException();
+            }
+            setNormAndDist(vertices, cache, i, j);
+            computeVotes(vertices, cache, maxDist);
+        }
+
+        /**
+         *
+         * @param vertices Vertices
+         * @param cache Cache
+         * @param i Index of the first point
+         * @param j Index of the second point
+         */
+        private void setNormAndDist(List<MeshPoint> vertices, SymmetryCache cache, int i, int j) {
+            MeshPoint meshPointI = vertices.get(i);
+            MeshPoint meshPointJ = vertices.get(j);
+
+            // accept only points with similar Gaussian curvature
+            if (!similarCurvature(meshPointI.getCurvature(), meshPointJ.getCurvature(), cache.avgGaussianCurvature * AVG_CURVATURE_MULTIPLICATOR)) {
+                throw new UnsupportedOperationException();
+            }
+
+            Vector3d planeNormal = new Vector3d(meshPointI.getPosition());
+            planeNormal.sub(meshPointJ.getPosition());
+            planeNormal.normalize();
+
+            // accpect only point pair with oposite normals along with the plane normal:
+            double normCos = cache.normCosVec[i][j].dot(planeNormal);
+            if (Math.abs(normCos) < MIN_ANGLE_COS) {
+                throw new UnsupportedOperationException();
+            }
+
+            setDistance(-planeNormal.dot(cache.avgPos[i][j]));
+            setNormal(planeNormal);
+        }
+
+        /**
+         * Computes votes for given plane
+         *
+         * @param sigPoints Mesh vertices with the most significant curvature
+         * @param maxDist Distance limit
+         */
+        private void computeVotes(List<MeshPoint> vertices, SymmetryCache cache, double maxDist) {
+            normalize();
+            Vector3d normal = getNormal();
+            double d = getDistance();
+
+            for (int i = 0; i < vertices.size(); i++) {
+                for (int j = 0; j < vertices.size(); j++) {
+                    if (i == j) {
+                        continue;
+                    }
+
+                    if (!similarCurvature(vertices.get(i).getCurvature(), vertices.get(j).getCurvature(), cache.avgGaussianCurvature * AVG_CURVATURE_MULTIPLICATOR)) {
+                        continue;
+                    }
+
+                    double normCos = cache.normCosVec[i][j].dot(normal);
+                    if (Math.abs(normCos) < MIN_ANGLE_COS) {
+                        continue;
+                    }
+
+                    Vector3d vec = new Vector3d(vertices.get(i).getPosition());
+                    vec.sub(vertices.get(j).getPosition());
+                    vec.normalize();
+                    double cos = vec.dot(normal);
+                    if (Math.abs(cos) < MIN_ANGLE_COS) {
+                        continue;
+                    }
+
+                    if (Math.abs(normal.dot(cache.avgPos[i][j]) + d) <= maxDist) {
+                        votes++;
+                    }
+                }
+            }
+        }
+
+        private boolean similarCurvature(Curvature gi, Curvature gj, double gh) {
+            if (gi == null || gj == null || gi.getGaussian() == Double.NaN || gj.getGaussian() == Double.NaN) {
+                return true; // can't decide => continue in the computation
+            }
+            return (gi.getGaussian() * gj.getGaussian() > 0)
+                    && (Math.abs(gi.getGaussian()) >= gh)
+                    && (Math.abs(gj.getGaussian()) >= gh);
+        }
+
+        /**
+         * Enables to compare two approximate planes due to number of votes
+         *
+         * @param other plane to be compared
+         * @return number that decides which plane has more votes
+         */
+        @Override
+        public int compareTo(CandidatePlane other) {
+            return Integer.compare(votes, other.votes);
+        }
+
+        @Override
+        public String toString() {
+            return this.getNormal() + " " + getDistance() + " " + votes;
+        }
+    }
+    
+    
+    /*******************************************************************
+     * Precomputed values for the symmetry plane estimation.
+     *
+     * @author Radek Oslejsek
+     */
+    protected class SymmetryCache {
+
+        private Vector3d[][] normCosVec;
+        private Vector3d[][] avgPos;
+        private double       avgGaussianCurvature;
+
+        /**
+         * Constructor.
+         *
+         * @param vertices Mesh vertices
+         */
+        protected SymmetryCache(List<MeshPoint> vertices) {
+            if (vertices == null || vertices.isEmpty()) {
+                throw new IllegalArgumentException("points");
+            }
+
+            int size = vertices.size();
+
+            normCosVec = new Vector3d[size][size];
+            avgPos = new Vector3d[size][size];
+
+            int counter = 0;
+            for (int i = 0; i < size; i++) {
+                MeshPoint meshPointI = vertices.get(i);
+
+                Curvature gi = meshPointI.getCurvature();
+                if (gi != null && gi.getGaussian() != Double.NaN) {
+                    avgGaussianCurvature += gi.getGaussian();
+                    counter++;
+                }
+
+                for (int j = 0; j < vertices.size(); j++) {
+                    MeshPoint meshPointJ = vertices.get(j);
+
+                    Vector3d ni = new Vector3d(meshPointI.getNormal());
+                    ni.sub(meshPointJ.getNormal());
+                    ni.normalize();
+                    normCosVec[i][j] = ni;
+
+                    Vector3d avrg = new Vector3d(meshPointI.getPosition());
+                    Vector3d aux = new Vector3d(meshPointJ.getPosition());
+                    avrg.add(aux);
+                    avrg.scale(0.5);
+                    avgPos[i][j] = avrg;
+                }
+            }
+            if (counter != 0) {
+                avgGaussianCurvature /= counter;
+            } else {
+                avgGaussianCurvature = Double.NaN;
+            }
+        }
+    }
+}

Comparison/src/main/java/cz/fidentis/analyst/symmetry/SymmetryEstimatorRobust.java

+package cz.fidentis.analyst.symmetry;
+
+import cz.fidentis.analyst.grid.UniformGrid3d;
+import cz.fidentis.analyst.grid.UniformGrid4d;
+import cz.fidentis.analyst.mesh.core.MeshFacet;
+import cz.fidentis.analyst.mesh.core.MeshPoint;
+import cz.fidentis.analyst.mesh.core.MeshPointImpl;
+import cz.fidentis.analyst.visitors.mesh.sampling.PointSampling;
+import java.util.ArrayList;
+import java.util.List;
+import java.util.Set;
+import java.util.stream.Collectors;
+import javax.vecmath.Point3d;
+import javax.vecmath.Vector3d;
+
+/**
+ * A robust implementation of symmetry plane estimation.
+ * The code is based based on the
+ * https://link.springer.com/article/10.1007/s00371-020-02034-w paper.
+ * This estimator <b>works with point clouds</b> (does not require manifold triangle mesh).
+ * It has the following properties:
+ * <ul>
+ * <li>The Wendland’s similarity functions is used to get best candidate planes.</li>
+ * <li>No additional weights are used.</li>
+ * <li>The computation is accelerated by using a Uniform Grid for pruning candidate planes
+ * and using two level downsampling: the radical downsampling for the generation 
+ * of candidate planes (cca. 100 points), and less radical (cca. 1000 point) 
+ * downsampling for the selection of the best candidate(s).</li>
+ * <li>Best results are achieved with Uniform Grid sampling with 100 and 1000 
+ * points (for the two phases).</li>
+ * </ul>
+ * 
+ * @author Radek Oslejsek
+ */
+public class SymmetryEstimatorRobust extends SymmetryEstimator {
+
+    private final PointSampling samplingStrategy;
+    private final int samplingLimit1;
+    private final int samplingLimit2;
+    
+    private Plane symmetryPlane;
+    
+    
+    /**
+     * Constructor.
+     * 
+     * @param samplingStrategy Downsampling strategy. Must not be {@code null}
+     *        Use {@code NoSumpling} strategy to avoid downsampling.
+     * @param samplingLimit1 Desired number of samples for finding candidate planes. 
+     * @param samplingLimit2 Desired number of samples for finding the best candidate. 
+     */
+    public SymmetryEstimatorRobust(PointSampling samplingStrategy, int samplingLimit1, int samplingLimit2) {
+        if (samplingStrategy == null) {
+            throw new IllegalArgumentException("samplingStrategy");
+        }
+        this.samplingStrategy = samplingStrategy;
+        this.samplingLimit1 = samplingLimit1;
+        this.samplingLimit2 = samplingLimit2;
+    }
+    
+    /**
+     * Computes and returns the symmetry plane. The plane is computed only once, 
+     * then the same instance is returned until the visitor is applied to another mesh.
+     * 
+     * @return the symmetry plane or {@code null}
+     */
+    public Plane getSymmetryPlane() {
+        if (symmetryPlane == null) {
+            this.calculateSymmetryPlane();
+        }
+        return symmetryPlane;
+    }
+    
+    @Override
+    public void visitMeshFacet(MeshFacet facet) {
+        samplingStrategy.visitMeshFacet(facet);
+    }
+    
+    /**
+     * Calculates the symmetry plane.
+     */
+    protected void calculateSymmetryPlane() {
+        //
+        // Phase 1: Downsample the mesh, transform it so that the centroid is 
+        //          in the space origin, and then compute candidate planes
+        //
+        samplingStrategy.setRequiredSamples(samplingLimit1);
+        List<MeshPoint> meshSamples = samplingStrategy.getSamples();
+        Point3d origCentroid = new MeshPointImpl(meshSamples).getPosition();
+        Set<SymmetryPlane> candidates = generateCandidates(meshSamples, origCentroid);
+        
+        //
+        // Phase 2: Downsample the mesh again (ususally to more samples than before), 
+        //          transform them in the same way as before, and measure their symmetry
+        //          with respect to individual candiate planes.
+        //
+        samplingStrategy.setRequiredSamples(samplingLimit2);
+        meshSamples = samplingStrategy.getSamples();
+        measureSymmetry(meshSamples, origCentroid, candidates);
+        
+        //
+        // Phase 3: "Adjust" the best 5 candidate planes so that they really 
+        //          represent the local maxima using 
+        //          the quasi-Newton optimization method L-BFGS
+        // 
+        // TO DO: candidates.stream().sorted().limit(5).forEach(optimize);
+        
+        //
+        // Phase 4: Finally, get the best plane and move it back
+        // 
+        this.symmetryPlane = candidates.stream().sorted().limit(1).findAny().orElse(null);
+        if (this.symmetryPlane != null) {
+            Point3d planePoint = this.symmetryPlane.getPlanePoint();
+            planePoint.add(origCentroid);
+            Vector3d normal = this.symmetryPlane.getNormal();
+            double dist = ((normal.x * planePoint.x) + (normal.y * planePoint.y) + (normal.z * planePoint.z))
+                    / Math.sqrt(normal.dot(normal)); // distance of tranformed surface point in the plane's mormal direction
+            this.symmetryPlane = new Plane(normal, dist);
+        }
+    }
+    
+    /**
+     * Copies mesh samples, moves them to the space origin, and then computes candidate planes.
+     * 
+     * @param meshSamples Downsampled mesh
+     * @param centroid Centroid of the downsampled mesh
+     * @return Candidate planes
+     */
+    protected Set<SymmetryPlane> generateCandidates(List<MeshPoint> meshSamples, Point3d centroid) {
+        ProcessedCloud cloud = new ProcessedCloud(meshSamples, centroid);
+        UniformGrid4d<SymmetryPlane> planesCache = new UniformGrid4d<>(SymmetryPlane.GRID_SIZE);
+        
+        for (int i = 0; i < cloud.points.size(); i++) {
+            for (int j = i; j < cloud.points.size(); j++) { // from i !!!
+                if (i == j) {
+                    continue;
+                }
+                SymmetryPlane candPlane = new SymmetryPlane(
+                        cloud.points.get(i).getPosition(),
+                        cloud.points.get(j).getPosition(),
+                        cloud.avgDistance);
+                SymmetryPlane closestPlane = candPlane.getClosestPlane(planesCache);
+                if (closestPlane == null) { // add
+                    planesCache.store(candPlane.getEstimationVector(), candPlane);
+                } else { // replace with averaged plane 
+                    SymmetryPlane avgPlane = new SymmetryPlane(closestPlane, candPlane);
+                    planesCache.remove(closestPlane.getEstimationVector(), closestPlane);
+                    planesCache.store(avgPlane.getEstimationVector(), avgPlane);
+                }
+            }
+        }
+
+        return planesCache.getAll().stream()
+                .filter(plane -> plane.getNumAverages() >= 4)
+                .collect(Collectors.toSet()); 
+    }
+    
+    /**
+     * Copies mesh samples, moves them to the space origin, and then measures the quality of 
+     * candidate planes. The results are stored in the candidate planes.
+     * 
+     * @param meshSamples Downsampled mesh
+     * @param centroid Centroid of the downsampled mesh
+     * @param candidates Candidate planes
+     */
+    protected void measureSymmetry(List<MeshPoint> meshSamples, Point3d centroid, Set<SymmetryPlane> candidates) {
+        ProcessedCloud cloud = new ProcessedCloud(meshSamples, centroid);
+        double alpha = 15.0 / cloud.avgDistance;
+        UniformGrid3d<MeshPoint> samplesGrid = new UniformGrid3d<>(2.6 / alpha, cloud.points, (MeshPoint p) -> p.getPosition());
+        candidates.parallelStream().forEach(plane -> {
+            plane.measureSymmetry(cloud.points, samplesGrid, alpha);
+        });
+    }
+    
+    
+    
+    
+    /********************************************************************
+     * A helper class that copies input mesh point and moves them 
+     * so that the given centroid is in the space origin.
+     * 
+     * @author Radek Oslejsek
+     */
+    protected class ProcessedCloud {
+
+        private List<MeshPoint> points;
+        private double avgDistance;
+
+        /**
+         * Moves orig points so that the centroid is in the space origin, copies
+         * them into a new list. Also, computes the average distance of orig
+         * points to the centroid (= average distance of shifted points into the
+         * space origin).
+         *
+         * @param centroid Centroid.
+         * @param points Original mesh point
+         */
+        ProcessedCloud(List<MeshPoint> points, Point3d centroid) {
+            this.points = new ArrayList<>(points.size());
+            for (int i = 0; i < points.size(); i++) {
+                MeshPoint mp = points.get(i);
+                Point3d p = new Point3d(mp.getPosition());
+                p.sub(centroid);
+                this.points.add(new MeshPointImpl(p, mp.getNormal(), mp.getTexCoord(), mp.getCurvature()));
+                avgDistance += Math.sqrt(p.x * p.x + p.y * p.y + p.z * p.z) / points.size();
+            }
+        }
+    }
+}

Comparison/src/main/java/cz/fidentis/analyst/symmetry/SymmetryEstimatorRobustMesh.java

+package cz.fidentis.analyst.symmetry;
+
+import cz.fidentis.analyst.grid.UniformGrid3d;
+import cz.fidentis.analyst.mesh.core.MeshFacet;
+import cz.fidentis.analyst.mesh.core.MeshPoint;
+import cz.fidentis.analyst.mesh.core.MeshPointImpl;
+import cz.fidentis.analyst.visitors.mesh.CurvatureCalculator;
+import cz.fidentis.analyst.visitors.mesh.sampling.PointSampling;
+import java.util.ArrayList;
+import java.util.HashSet;
+import java.util.List;
+import java.util.Set;
+import javax.vecmath.Point3d;
+import javax.vecmath.Vector3d;
+
+/**
+ * A robust implementation of symmetry plane estimation for manifold triangle meshes.
+ * The code is based based on the
+ * https://link.springer.com/article/10.1007/s00371-020-02034-w paper.
+ * This extension uses the similarity of Gaussian curvature values (the static weight)
+ * and the symmetry of normal vectors (the dynamic weight) in the computation
+ * along with the Wendland’s similarity function. The computation is bit slower
+ * in general (depends on parameters set, of course) that the super-class, 
+ * but the weights should help the algorithm to <b>better deal with demaged faces, 
+ * i.e., when only a small fragment of a face is available</b>.
+ * <p>
+ * Although the pruning of candidate planes and the measurement of their quality
+ * is similar to the {@link SymmetryEstimatorMesh} implementation, this code 
+ * is much slower with similar quality of results :-( Probably the price for 
+ * robustness (need to be further researched).
+ * </p>
+ * 
+ * @author Radek Oslejsek
+ */
+public class SymmetryEstimatorRobustMesh extends SymmetryEstimatorRobust {
+    
+    /**
+     * Constructor.
+     * 
+     * @param samplingStrategy Downsampling strategy. Must not be {@code null}
+     *        Use {@code NoSumpling} strategy to avoid downsampling.
+     * @param samplingLimit1 Desired number of samples for finding candidate planes. 
+     * @param samplingLimit2 Desired number of samples for finding the best candidate. 
+     */
+    public SymmetryEstimatorRobustMesh(PointSampling samplingStrategy, int samplingLimit1, int samplingLimit2) {
+        super(samplingStrategy, samplingLimit1, samplingLimit2);
+    }
+    
+    @Override
+    public void visitMeshFacet(MeshFacet facet) {
+        // We need vertex normals
+        synchronized (this) {
+            if (!facet.hasVertexNormals()) {
+                facet.calculateVertexNormals();
+            }
+        }
+        
+        // We need Gaussian curvatures
+        if (facet.getVertex(0).getCurvature() == null) {
+            facet.accept(new CurvatureCalculator());
+        }
+        
+        super.visitMeshFacet(facet);
+    }
+    
+    /**
+     * Copies mesh samples, moves them to the space origin, and then computes candidate planes.
+     * This implementation uses curvature and normal vectors to prune candidates.
+     * 
+     * @param meshSamples Downsampled mesh
+     * @param centroid Centroid of the downsampled mesh
+     * @return Candidate planes
+     */
+    protected Set<SymmetryPlane> generateCandidates(List<MeshPoint> meshSamples, Point3d centroid) {
+        ProcessedCloud cloud = new ProcessedCloud(meshSamples, centroid);
+        Set<SymmetryPlane> ret = new HashSet<>();
+        
+        for (int i = 0; i < cloud.points.size(); i++) {
+            MeshPoint pi = cloud.points.get(i);
+            for (int j = i; j < cloud.points.size(); j++) { // from i !!!
+                if (i == j) {
+                    continue;
+                }
+                MeshPoint pj = cloud.points.get(j);
+                
+                if (cloud.getCurvatureSumilarity(pi, pj) <= 0.0) {
+                    continue;
+                }
+                
+                SymmetryPlaneCur candPlane = new SymmetryPlaneCur(
+                        pi.getPosition(),
+                        pj.getPosition(),
+                        cloud.avgDistance
+                );
+                
+                Vector3d nri = candPlane.reflectUnitVectorOverPlane(pi.getNormal());
+                if (candPlane.getNormalVectorsWeight(nri, pj.getNormal()) <= 0.25) {
+                    continue;
+                }
+                
+                ret.add(candPlane);
+            }
+        }
+        
+        return ret;
+    }
+    
+    /**
+     * Copies mesh samples, moves them to the space origin, and then measures the quality of 
+     * candidate planes. The results are stored in the candidate planes.
+     * 
+     * @param meshSamples Downsampled mesh
+     * @param centroid Centroid of the downsampled mesh
+     * @param candidates Candidate planes
+     */
+    @Override
+    protected void measureSymmetry(List<MeshPoint> meshSamples, Point3d centroid, Set<SymmetryPlane> candidates) {
+        ProcessedCloud cache = new ProcessedCloud(meshSamples, centroid);
+        double alpha = 15.0 / cache.avgDistance;
+        UniformGrid3d<MeshPoint> samplesGrid = new UniformGrid3d<>(2.6 / alpha, cache.points, (MeshPoint p) -> p.getPosition());
+        candidates.parallelStream().forEach(plane -> {
+            ((SymmetryPlaneCur) plane).measureWeightedSymmetry(cache, samplesGrid, alpha);
+        });
+    }
+    
+    /********************************************************************
+     * A helper class that copies input mesh point and moves them 
+     * so that the given centroid is in the space origin.
+     * Also, it computes some additional data
+     * 
+     * @author Radek Oslejsek
+     */
+    public class ProcessedCloud {
+        
+        public static final double AVG_CUR_MAGIC_MULTIPLIER = 0.01;
+
+        private List<MeshPoint> points;
+        private double avgDistance;
+        
+        /**
+         * The average of absolute values of Gaussian curvatures in all points
+         */
+        private double avgCurvature;
+
+        /**
+         * Moves orig points so that the centroid is in the space origin, copies
+         * them into a new list. Also, computes the average distance of orig
+         * points to the centroid (= average distance of shifted points into the
+         * space origin).
+         *
+         * @param centroid Centroid.
+         * @param points Original mesh point
+         */
+        public ProcessedCloud(List<MeshPoint> points, Point3d centroid) {
+            int size = points.size();
+            this.points = new ArrayList<>(size);
+            
+            int counter = 0;
+            for (int i = 0; i < size; i++) {
+                MeshPoint mp = points.get(i);
+                Point3d p = new Point3d(mp.getPosition());
+                p.sub(centroid);
+                this.points.add(new MeshPointImpl(p, mp.getNormal(), mp.getTexCoord(), mp.getCurvature()));
+                avgDistance += Math.sqrt(p.x * p.x + p.y * p.y + p.z * p.z) / size;
+                
+                if (mp.getCurvature() != null && mp.getCurvature().getGaussian() != Double.NaN) {
+                    avgCurvature += Math.abs(mp.getCurvature().getGaussian());
+                    counter++;
+                }
+            }
+            
+            if (counter != 0) {
+                avgCurvature /= counter;
+            } else {
+                avgCurvature = Double.NaN;
+            }
+        }
+        
+        /**
+         * Returns similarity of Gaussian curvature.
+         * 
+         * @param p1 First point
+         * @param p2 Second point
+         * @return similarity of Gaussian curvature.
+         */
+        public double getCurvatureSumilarity(MeshPoint p1, MeshPoint p2) {
+            double avgh = avgCurvature * AVG_CUR_MAGIC_MULTIPLIER;
+            double g1 = p1.getCurvature().getGaussian();
+            double g2 = p2.getCurvature().getGaussian();
+            
+            if (g1 == Double.NaN || g2 == Double.NaN) {
+                return 0;
+            }
+
+            if (g1 * g1 <= 0) {
+                return 0;
+            }
+
+            if (Math.abs(g1) < avgh) {
+                return 0;
+            }
+
+            if (Math.abs(g2) < avgh) {
+                return 0;
+            }
+
+            double weight = Math.min(Math.abs(g1), Math.abs(g2)) / Math.max(Math.abs(g1), Math.abs(g2));
+
+            if (weight <= 0) {
+                return 0;
+            }
+
+            return weight;
+        }
+        
+        /**
+         * Returns i-th transformed point
+         * @param i index
+         * @return i-th transformed point
+         */
+        public MeshPoint getPoint(int i) {
+            return points.get(i);
+        }
+        
+        /**
+         * Returns the number of points.
+         * @return the number of points.
+         */
+        public int getNumPoints() {
+            return points.size();
+        }
+        
+        public double getAvgCurvature() {
+            return this.avgCurvature;
+        }
+        
+        public double getAvgDistance() {
+            return this.avgDistance;
+        }
+        
+    }        
+}

Comparison/src/main/java/cz/fidentis/analyst/symmetry/SymmetryPlane.java

+package cz.fidentis.analyst.symmetry;
+
+import cz.fidentis.analyst.grid.UniformGrid3d;
+import cz.fidentis.analyst.grid.UniformGrid4d;
+import cz.fidentis.analyst.mesh.core.MeshPoint;
+import java.util.List;
+import javax.vecmath.Point3d;
+import javax.vecmath.Tuple3d;
+import javax.vecmath.Vector3d;
+import javax.vecmath.Vector4d;
+
+/**
+ * A symmetry plane extends standard plane with functions related to the
+ * similarity of planes and the measure of their quality (symmetry precision).
+ *
+ * @author Radek Oslejsek
+ */
+public class SymmetryPlane extends Plane implements Comparable<SymmetryPlane> {
+
+    public static final double GRID_SIZE = 0.1; // similarity epsilon, should be 0.1
+
+    private final double avgDist;
+    private int numAvg; // psi
+    private double symmetry = 0;
+
+    /**
+     * New candidate symmetry plane constructed from two points.
+     *
+     * @param point1 point in space
+     * @param point2 point in space
+     * @param avgDist the average distance between vertices and their centroid
+     * @throws IllegalArgumentException if the @code{plane} argument is null
+     */
+    public SymmetryPlane(Tuple3d point1, Tuple3d point2, double avgDist) {
+        super(point1, point2);
+        this.avgDist = avgDist;
+        this.numAvg = 0;
+    }
+
+    /**
+     * Candidate symmetry plane constructed by averaging two planes.
+     *
+     * @param pu the closest existing candidate
+     * @param pv newly created candidate
+     */
+    public SymmetryPlane(SymmetryPlane pu, SymmetryPlane pv) {
+        Vector3d pun = pu.getNormal();
+        Vector3d pvn = pv.getNormal();
+        double pud = pu.getDistance();
+        double pvd = pv.getDistance();
+
+        // new created plane has numAvg set to zero and normalized normal cetor
+        if (pu.getNumAverages() >= 1) {
+            double scale = 1.0 / pu.getNormal().length();
+            pun.scale(scale);
+            pud *= scale;
+            if (pu.getNumAverages() > 1) {
+                scale = 1.0 / pu.getNumAverages();
+                pvn.scale(scale);
+                pvd *= scale;
+            }
+        }
+
+        // change pun and pud values to correspond to the average plane:
+        double dot = pun.x * pvn.x + pun.y * pvn.y + pun.z * pvn.z;
+        if (dot >= 0) {
+            pun.add(pvn);
+            pud += pvd;
+        } else {
+            pun.sub(pvn);
+            pud -= pvd;
+        }
+
+        setNormal(pun);
+        setDistance(pud);
+        this.avgDist = pu.avgDist;
+        this.numAvg = pu.getNumAverages() + 1;
+    }
+
+    /**
+     * Returns the closest plane from the cache to this plane.
+     *
+     * @param cache
+     * @return the closest plane from the cache to this plane.
+     */
+    public SymmetryPlane getClosestPlane(UniformGrid4d<SymmetryPlane> cache) {
+        Vector4d pp = this.getEstimationVector();
+        List<SymmetryPlane> closePlanes = cache.getClosest(pp);
+        pp.scale(-1.0);
+        closePlanes.addAll(cache.getClosest(pp));
+
+        SymmetryPlane closest = null;
+        double dist = Double.POSITIVE_INFINITY;
+        for (var plane : closePlanes) {
+            double d = this.distance(plane);
+            if (d < GRID_SIZE && d < dist) {
+                dist = d;
+                closest = plane;
+            }
+        }
+
+        return (dist == Double.POSITIVE_INFINITY) ? null : closest;
+    }
+    
+    /**
+     * Returns a 4D vector usable as a location for the storage in the 
+     * {@link UniformGrid4d}
+     * 
+     * @return a 4D vector usable as a location for the storage in the 
+     * {@link UniformGrid4d}
+     */
+    public Vector4d getEstimationVector() {
+        Vector3d n = getNormalReference();
+        Vector4d ret = new Vector4d(n);
+        ret.w = getDistance() / avgDist;
+        if (numAvg > 0) { // a non-averaged plane has to be normalized
+            ret.scale(1.0 / n.length());
+        }
+        return ret;
+    }
+
+    public int getNumAverages() {
+        return numAvg;
+    }
+
+    /**
+     *
+     * @param pv the second plane
+     * @return planes distance
+     */
+    protected double distance(SymmetryPlane pv) {
+        Vector4d pud = this.getEstimationVector();
+        Vector4d pvd = pv.getEstimationVector();
+        if (this.getNormalReference().dot(pv.getNormalReference()) >= 0.0) {
+            pud.sub(pvd);
+        } else {
+            pud.add(pvd);
+        }
+        return pud.length();
+    }
+
+    /**
+     * Symmetry measurement based on the Wendland’s function without additional
+     * weights
+     * <b>The plane is normalized and {@code numAvg} set to zero!</b>
+     *
+     * @param points Downsampled point cloud
+     * @param grid The same cloud stored in the uniform grid
+     * @param alpha the average distance between vertices and their centroid
+     */
+    public void measureSymmetry(List<MeshPoint> points, UniformGrid3d<MeshPoint> grid, double alpha) {
+        normalizeIfNeeded();
+        symmetry = 0.0;
+
+        for (int i = 0; i < points.size(); i++) {
+            MeshPoint p1 = points.get(i); // original point
+            Point3d rp1 = reflectPointOverPlane(points.get(i).getPosition()); // reflected point
+            List<MeshPoint> closest = grid.getClosest(rp1);
+
+            for (int j = 0; j < closest.size(); j++) {
+                MeshPoint p2 = closest.get(j);
+
+                if (p1 == p2) {
+                    continue;
+                }
+
+                double phi = similarityFunction(rp1.distance(p2.getPosition()), alpha);
+                if (phi == 0) {
+                    continue;
+                }
+
+                symmetry += phi;
+            }
+        }
+    }
+
+    /**
+     * Returns the symmetry measure computed by the {@link #getSymmetryMeasure()}.
+     * @return the symmetry measure
+     */
+    public double getSymmetryMeasure() {
+        return this.symmetry;
+    }
+
+    @Override
+    public int compareTo(SymmetryPlane o) {
+        return Double.compare(o.symmetry, this.symmetry);
+    }
+
+    @Override
+    public String toString() {
+        return "S: " + symmetry + " " + super.toString();
+    }
+
+    /**
+     * A similarity function Phi.
+     *
+     * @param length
+     * @param alpha
+     * @return similarity value
+     */
+    protected static double similarityFunction(double length, double alpha) {
+        double al = alpha * length;
+        if (al == 0 || al > 2.6) {
+            return 0;
+        }
+        double q = al / 2.6;
+        return Math.pow(1.0 - q, 5) * (8 * q * q + 5 * q + 1);
+    }
+
+    protected void setSymmetryMeasure(double symmetry) {
+        this.symmetry = symmetry;
+    }
+
+    protected void normalizeIfNeeded() {
+        if (numAvg > 0) {
+            normalize();
+            numAvg = 0;
+        }
+    }
+}

Comparison/src/main/java/cz/fidentis/analyst/symmetry/SymmetryPlaneCur.java

+package cz.fidentis.analyst.symmetry;
+
+import cz.fidentis.analyst.grid.UniformGrid3d;
+import cz.fidentis.analyst.mesh.core.MeshPoint;
+import cz.fidentis.analyst.symmetry.SymmetryEstimatorRobustMesh.ProcessedCloud;
+import java.util.List;
+import javax.vecmath.Point3d;
+import javax.vecmath.Vector3d;
+
+/**
+ * This symmetry plane changes the behavior so that similarity of Gaussian
+ * curvature values and the symmetry of normal vectors are used to check the
+ * quality of the symmetry plane.
+ *
+ * @author Radek Oslejsek
+ */
+public class SymmetryPlaneCur extends SymmetryPlane {
+    
+    /**
+     * New candidate symmetry plane constructed from two points.
+     *
+     * @param point1 point in space
+     * @param point2 point in space
+     * @param avgDist the average distance between vertices and their centroid
+     * @throws IllegalArgumentException if the @code{plane} argument is null
+     */
+    public SymmetryPlaneCur(Point3d point1, Point3d point2, double avgDist) {
+        super(point1, point2, avgDist);
+    }
+
+    /**
+     * Symmetry measurement based on the similarity of Gaussian curvatures and
+     * the symmetry of normal vectors, in addition to Wendland’s similarity function.<b>The plane is normalized and {@code numAvg} set to zero!</b>
+     *
+     * @param cache Precomputed values, including downsampled points
+     * @param grid The same cloud stored in the uniform grid
+     * @param alpha the average distance between vertices and their centroid
+     */
+    public void measureWeightedSymmetry(ProcessedCloud cache, UniformGrid3d<MeshPoint> grid, double alpha) {
+        normalizeIfNeeded();
+        setSymmetryMeasure(0.0);
+
+        for (int i = 0; i < cache.getNumPoints(); i++) {
+            MeshPoint p1 = cache.getPoint(i); // original point
+            Point3d rp1 = reflectPointOverPlane(p1.getPosition()); // reflected point
+            Vector3d rn1 = reflectUnitVectorOverPlane(p1.getNormal()); // reflected normal
+            List<MeshPoint> closest = grid.getClosest(rp1);
+
+            for (int j = 0; j < closest.size(); j++) {
+                MeshPoint p2 = closest.get(j);
+
+                if (p1 == p2) {
+                    continue;
+                }
+
+                double ws = cache.getCurvatureSumilarity(p1, p2);
+                if (ws == 0) {
+                    continue;
+                }
+
+                double wd = getNormalVectorsWeight(rn1, p2.getNormal());
+                if (wd == 0) {
+                    continue;
+                }
+
+                double phi = similarityFunction(rp1.distance(p2.getPosition()), alpha);
+                if (phi == 0) {
+                    continue;
+                }
+
+                setSymmetryMeasure(getSymmetryMeasure() + ws * wd * phi); // increment the measure
+            }
+        }
+    }
+    
+    /**
+     * Returns the symmetry of normal vectors
+     * 
+     * @param rv1 A normal vector already reflected over this plane
+     * @param v2 A second normal vector for comparison
+     * @return the symmetry of the given normal vectors
+     */
+    protected double getNormalVectorsWeight(Vector3d rv1, Vector3d v2) {
+        double cosN = rv1.dot(v2);
+        cosN = (cosN >= 0) ? Math.min(1.0, cosN) : Math.max(-1.0, cosN);
+        double acos = Math.acos(cosN);
+        return similarityFunction(acos, 4.0);
+    }
+}

Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/Curvature.java → Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/CurvatureCalculator.java

@@ -16,7 +16,8 @@ import javax.vecmath.Point3d;
 import javax.vecmath.Vector3d;
 
 /**
- * Abstract class for algorithms calculating curvatures of mesh vertices.
+ * A visitor that computes and sets curvatures of visited meshes. 
+ * It returns nothing.
  * @see https://computergraphics.stackexchange.com/questions/1718/what-is-the-simplest-way-to-compute-principal-curvature-for-a-mesh-triangle
  * @see http://rodolphe-vaillant.fr/?e=20
  * <p>
@@ -26,12 +27,7 @@ import javax.vecmath.Vector3d;
  * @author Natalia Bebjakova
  * @author Radek Oslejsek
  */
-public class Curvature extends MeshVisitor {
-    
-    private final Map<MeshFacet, List<Double>> gaussian = new HashMap<>();
-    private final Map<MeshFacet, List<Double>> mean = new HashMap<>();
-    private final Map<MeshFacet, List<Double>> minPrincipal = new HashMap<>();
-    private final Map<MeshFacet, List<Double>> maxPrincipal = new HashMap<>();
+public class CurvatureCalculator extends MeshVisitor {
     
     private final boolean parallel;
     
@@ -40,7 +36,7 @@ public class Curvature extends MeshVisitor {
      * 
      * @param parallel If {@code true}, then the algorithm runs concurrently utilizing all CPU cores
      */
-    public Curvature(boolean parallel) {
+    public CurvatureCalculator(boolean parallel) {
         this.parallel = parallel;
     }
     
@@ -48,63 +44,22 @@ public class Curvature extends MeshVisitor {
      * Constructor for parallel computation.
      * 
      */
-    public Curvature() {
+    public CurvatureCalculator() {
         this(true);
     }
     
-    /**
-     * Returns Gaussian 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 Gaussian curvatures of inspected mesh facets.
-     */
-    public Map<MeshFacet, List<Double>> getGaussianCurvatures() {
-        return Collections.unmodifiableMap(gaussian);
-    }
-    
-    /**
-     * Returns mean 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 Mean curvatures of inspected mesh facets.
-     */
-    public Map<MeshFacet, List<Double>> getMeanCurvatures() {
-        return Collections.unmodifiableMap(mean);
-    }
-    
-    /**
-     * Returns minimum 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 Minimum principal curvatures of inspected mesh facets.
-     */
-    public Map<MeshFacet, List<Double>> getMinPrincipalCurvatures() {
-        return Collections.unmodifiableMap(minPrincipal);
-    }
-    
-    /**
-     * 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.
-     */
-    public Map<MeshFacet, List<Double>> getMaxPrincipalCurvatures() {
-        return Collections.unmodifiableMap(maxPrincipal);
+    @Override
+    public boolean isThreadSafe() {
+        return false;
     }
     
     @Override
     public void visitMeshFacet(final MeshFacet facet) {
-        synchronized (this) {
-            if (gaussian.containsKey(facet)) {
-                return; // already visited facet
-            }
-            // prefill the arrays - required for concurrent computation.
-            int numVert = facet.getNumberOfVertices();
-            gaussian.put(facet, new ArrayList<Double>(Collections.nCopies(numVert, Double.NaN)));
-            mean.put(facet, new ArrayList<Double>(Collections.nCopies(numVert, Double.NaN)));
-            minPrincipal.put(facet, new ArrayList<Double>(Collections.nCopies(numVert, Double.NaN)));
-            maxPrincipal.put(facet, new ArrayList<Double>(Collections.nCopies(numVert, Double.NaN)));
-        }
+        int numVert = facet.getNumberOfVertices();
+        List<Double> gaussian = new ArrayList<>(Collections.nCopies(numVert, Double.NaN));
+        List<Double> mean = new ArrayList<>(Collections.nCopies(numVert, Double.NaN));
+        List<Double> minPrincipal = new ArrayList<>(Collections.nCopies(numVert, Double.NaN));
+        List<Double> maxPrincipal = new ArrayList<>(Collections.nCopies(numVert, Double.NaN));
         
         final Cache cache = new Cache(facet);
         
@@ -115,7 +70,7 @@ public class Curvature extends MeshVisitor {
                 Runnable worker = new Runnable() {
                     @Override
                     public void run() {
-                        computeCurvature(facet, cache, index);
+                        computeCurvature(facet, cache, index, minPrincipal, maxPrincipal, mean, gaussian);
                     }
                 };
                 executor.execute(worker);
@@ -123,21 +78,33 @@ public class Curvature extends MeshVisitor {
             executor.shutdown();
             while (!executor.isTerminated()){}
         } else {
-            for (int vertIndexA = 0; vertIndexA < facet.getNumberOfVertices(); vertIndexA++) {
-                computeCurvature(facet, cache, vertIndexA);
+            for (int i = 0; i < facet.getNumberOfVertices(); i++) {
+                computeCurvature(facet, cache, i, minPrincipal, maxPrincipal, mean, gaussian);
             }
         }
+        
+        for (int i = 0; i < numVert; i++) {
+            facet.getVertex(i).setCurvature(
+                    new cz.fidentis.analyst.mesh.core.Curvature(
+                            minPrincipal.get(i),
+                            maxPrincipal.get(i),
+                            gaussian.get(i),
+                            mean.get(i)
+                    )
+            );
+        }
     }
     
-    protected void computeCurvature(MeshFacet facet, Cache cache, int vertIndexA) {
+    protected void computeCurvature(MeshFacet facet, Cache cache, int vertIndexA,
+            List<Double> minPrincipal, List<Double> maxPrincipal, List<Double> mean, List<Double> gaussian) {
         TriangleFan oneRing = facet.getOneRingNeighborhood(vertIndexA);
 
         if (oneRing.isBoundary()) {
             synchronized (this) {
-                this.gaussian.get(facet).set(vertIndexA, 0.0);
-                this.mean.get(facet).set(vertIndexA, 0.0);
-                this.minPrincipal.get(facet).set(vertIndexA, 0.0);
-                this.maxPrincipal.get(facet).set(vertIndexA, 0.0);
+                gaussian.set(vertIndexA, 0.0);
+                mean.set(vertIndexA, 0.0);
+                minPrincipal.set(vertIndexA, 0.0);
+                maxPrincipal.set(vertIndexA, 0.0);
                 return;
             }
         }
@@ -195,10 +162,10 @@ public class Curvature extends MeshVisitor {
         double delta = Math.max(0, Math.pow(meanVal, 2) - gaussVal);
 
         synchronized (this) {
-            this.gaussian.get(facet).set(vertIndexA, gaussVal);
-            this.mean.get(facet).set(vertIndexA, meanVal);
-            this.minPrincipal.get(facet).set(vertIndexA, meanVal - Math.sqrt(delta));
-            this.maxPrincipal.get(facet).set(vertIndexA, meanVal + Math.sqrt(delta));
+            gaussian.set(vertIndexA, gaussVal);
+            mean.set(vertIndexA, meanVal);
+            minPrincipal.set(vertIndexA, meanVal - Math.sqrt(delta));
+            maxPrincipal.set(vertIndexA, meanVal + Math.sqrt(delta));
         }
     }
 
@@ -302,7 +269,8 @@ public class Curvature extends MeshVisitor {
                 return cache.get(key).get(0);
             } 
             
-            throw new IllegalArgumentException("[" + vCentral + ", " + v2 + ", " + v3 + "]");
+            //throw new IllegalArgumentException("[" + vCentral + ", " + v2 + ", " + v3 + "]");
+            return Double.NaN; // should not happen, but it appears for some models that have a triangle with two points at the same location.
         }
         
         /**
@@ -325,7 +293,8 @@ public class Curvature extends MeshVisitor {
                 return cache.get(key).get(1);
             } 
             
-            throw new IllegalArgumentException("[" + vCentral + ", " + v2 + ", " + v3 + "]");
+            //throw new IllegalArgumentException("[" + vCentral + ", " + v2 + ", " + v3 + "]");
+            return Double.NaN; // should not happen, but it appears for some models that have a triangle with two points at the same location.
         }
         
         private void computeValues(MeshFacet facet, MeshTriangle tri) {

Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/sampling/CurvatureSampling.java

 package cz.fidentis.analyst.visitors.mesh.sampling;
 
+import cz.fidentis.analyst.mesh.core.Curvature;
 import cz.fidentis.analyst.mesh.core.MeshFacet;
 import cz.fidentis.analyst.mesh.core.MeshPoint;
-import cz.fidentis.analyst.visitors.mesh.Curvature;
 import java.util.ArrayList;
 import java.util.List;
-import java.util.Map;
-import java.util.PriorityQueue;
 import java.util.stream.Collectors;
 
 /**
  * A relevance-based point sampling using highest curvature.
  * {@see https://graphics.stanford.edu/papers/zipper/zipper.pdf}
- * 
- * <p>
- * This visitor <b>is not thread-safe</b> for the efficiency reasons, i.e., 
- * a single instance of the visitor cannot be used to inspect multiple meshes 
- * simultaneously (sequential inspection is okay).
- * </p>
+ * If curvature is not set in inspected facets, then it is computed and set automatically.
+ * This visitor is tread-safe.
  * 
  * @author Radek Oslejsek
  * @author Natalia Bebjakova
@@ -36,93 +30,74 @@ public class CurvatureSampling extends PointSampling {
         MIN
     }
 
-    //private final int maxPoints;
-    
-    private final Curvature curvatureVisitor;
-    
     private final CurvatureAlg curvatureAlg;
+    private List<MeshPoint> allVertices = new ArrayList<>();
     
-    private List<VertCurvature> significantPoints;
-    
-    /**
-     * Constructor.
+    /** 
+     * Constructor for no point sampling (100 %).
      * 
      * @param cAlg Curvature strategy to be used for the selection of significant points.
-     * @param max Maximal number of significant points
-     * @throws IllegalArgumentException if the {@code curbatureVisitor} is missing
      */
-    public CurvatureSampling(CurvatureAlg cAlg, int max) {
-        this(new Curvature(), cAlg, max);
+    public CurvatureSampling(CurvatureAlg cAlg) {
+        this(cAlg, 1.0);
     }
-    
-    /**
-     * Constructor.
-     * 
-     * @param cAlg Curvature strategy to be used for the selection of significant points.
-     * @param perc Maximal number of significant points as percents of the original size
-     */
-    public CurvatureSampling(CurvatureAlg cAlg, double perc) {
-        this(new Curvature(), cAlg, perc);
-    }
-    
+
     /**
      * Constructor.
      * 
-     * @param curvatureVisitor Curvature. If the object has
-     * pre-filled curvatures of meshes, then they are also used for the computation
-     * of significant points. In this way, it is possible to reuse already computed
-     * curvatures and skip calling the {@link #visitMeshFacet} inspection method.
      * @param cAlg Curvature strategy to be used for the selection of significant points.
      * @param max Maximal number of significant points
      * @throws IllegalArgumentException if the {@code curbatureVisitor} is missing
      */
-    public CurvatureSampling(Curvature curvatureVisitor, CurvatureAlg cAlg, int max) {
+    public CurvatureSampling(CurvatureAlg cAlg, int max) {
         super(max);
-        
-        if (curvatureVisitor == null) {
-            throw new IllegalArgumentException("curvatureVisitor");
-        }
-
-        this.curvatureVisitor = curvatureVisitor;
         this.curvatureAlg = cAlg;
     }
     
     /**
      * Constructor.
      * 
-     * @param curvatureVisitor Curvature. If the object has
-     * pre-filled curvatures of meshes, then they are also used for the computation
-     * of significant points. In this way, it is possible to reuse already computed
-     * curvatures and skip calling the {@link #visitMeshFacet} inspection method.
      * @param cAlg Curvature strategy to be used for the selection of significant points.
      * @param perc Maximal number of significant points as percents of the original size
      */
-    public CurvatureSampling(Curvature curvatureVisitor, CurvatureAlg cAlg, double perc) {
+    public CurvatureSampling(CurvatureAlg cAlg, double perc) {
         super(perc);
-        
-        if (curvatureVisitor == null) {
-            throw new IllegalArgumentException("curvatureVisitor");
-        }
-
-        this.curvatureVisitor = curvatureVisitor;
         this.curvatureAlg = cAlg;
     }
     
-    @Override
-    public boolean isThreadSafe() {
-        return false;
-    }
-    
     @Override
     public void visitMeshFacet(MeshFacet facet) {
-        significantPoints = null; // clear previous results
-        curvatureVisitor.visitMeshFacet(facet); // compute curvature for new inspected facet
+        if (facet != null) {
+            synchronized(this) {
+                if (facet.getVertex(0).getCurvature() == null) {
+                    facet.accept(new cz.fidentis.analyst.visitors.mesh.CurvatureCalculator());
+                }
+                allVertices.addAll(facet.getVertices());
+            }
+        }
     }
     
     @Override
     public List<MeshPoint> getSamples() {
-        checkAndComputeSignificantPoints();
-        return significantPoints.stream().map(vc -> vc.point).collect(Collectors.toList());
+        return allVertices.stream()
+                .sorted((MeshPoint mp1, MeshPoint mp2) -> {
+                    Curvature c1 = mp1.getCurvature();
+                    Curvature c2 = mp2.getCurvature();
+                    switch (curvatureAlg) {
+                        case MEAN:
+                            return Double.compare(c1.getMean(), c2.getMean());
+                        case GAUSSIAN:
+                            return Double.compare(c1.getGaussian(), c2.getGaussian());
+                        case MAX:
+                            return Double.compare(c1.getMaxPrincipal(), c2.getMaxPrincipal());
+                        case MIN:
+                            return Double.compare(c1.getMinPrincipal(), c2.getMinPrincipal());
+                        default:
+                            throw new IllegalArgumentException("curvatureAlg");
+                    }
+                })
+                .limit(getNumDownsampledPoints(allVertices.size()))
+                .collect(Collectors.toList());
     }
     
     public CurvatureAlg getCurvatureAlg() {
@@ -134,98 +109,14 @@ public class CurvatureSampling extends PointSampling {
      * 
      * @return curvatures of selected samples
      */
-    public List<Double> getSampledCurvatures() {
-        checkAndComputeSignificantPoints();
-        return significantPoints.stream().map(vc -> vc.curvature).collect(Collectors.toList());
-    }
+    //public List<Double> getSampledCurvatures() {
+    //    checkAndComputeSignificantPoints();
+    //    return significantPoints.stream().map(vc -> vc.curvature).collect(Collectors.toList());
+    //}
     
     @Override
     public String toString() {
         return this.curvatureAlg + " curvature " + super.toString();
     }
-
-    protected void checkAndComputeSignificantPoints() {
-        if (significantPoints != null) {
-            return; // already computed
-        }
-        
-        // fill the priority queue
-        final PriorityQueue<VertCurvature> priorityQueue = new PriorityQueue<>();
-        
-        Map<MeshFacet, List<Double>> curvMap = null;
-        switch (curvatureAlg) {
-            case MEAN: 
-                curvMap = curvatureVisitor.getMeanCurvatures();
-                break;
-            case GAUSSIAN: 
-                curvMap = curvatureVisitor.getGaussianCurvatures();
-                break;
-            case MAX: 
-                curvMap = curvatureVisitor.getMaxPrincipalCurvatures();
-                break;
-            case MIN: 
-                curvMap = curvatureVisitor.getMinPrincipalCurvatures();
-                break;
-            default:
-                throw new IllegalArgumentException("curvatureAlg");
-        }
-
-        for (MeshFacet facet: curvMap.keySet()) {
-            List<Double> curvs = curvMap.get(facet);
-            for (int i = 0; i < curvs.size(); i++) { 
-                // store helper objects, replace NaN with 0.0:
-                double c = curvs.get(i);
-                priorityQueue.add(new VertCurvature(facet.getVertex(i), (Double.isNaN(c) ? 0.0 : c)));
-            }  
-        }
-        
-        // select top significant points
-        int maxPoints = getNumDownsampledPoints(priorityQueue.size());
-        significantPoints = new ArrayList<>(maxPoints);
-        for (int i = 0; i < maxPoints; i++) {
-            VertCurvature vc = priorityQueue.poll();
-            if (vc == null) {
-                break; // no more points available
-            } else {
-                significantPoints.add(vc);
-            }
-        }
-    }
-    
-    /**
-     * Helper class for sorting points with respect to their curvature.
-     * @author Radek Oslejsek
-     */
-    private class VertCurvature implements Comparable<VertCurvature> {
-        
-        private final double curvature;
-        private final MeshPoint point;
-        
-        VertCurvature(MeshPoint point, double curvature) {
-            this.curvature = curvature;
-            this.point = point;
-        }
-
-        /**
-         * Compares this object with the specified object for order. 
-         * Curvature is taken into consideration as the primary value (descendant ordering).
-         * If the curvature equals, then the objects are sorted randomly (1 is always returned).
-         * This approach preserves all points in sorted collections.
-         * 
-         * @param arg the object to be compared.
-         * @return a negative integer, zero, or a positive integer as this object 
-         * is less than, equal to, or greater than the specified object.
-         */
-        @Override
-        public int compareTo(VertCurvature arg) {
-            int comp = Double.compare(arg.curvature, this.curvature);
-            return (comp == 0) ? 1 : comp;
-        }
-        
-        @Override
-        public String toString() {
-            return ""+curvature;
-        }
-    }
     
 }

Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/sampling/NoSampling.java

@@ -19,9 +19,9 @@ public class NoSampling extends PointSampling {
     public void visitMeshFacet(MeshFacet facet) {
         if (facet != null) {
             allVertices.addAll(facet.getVertices());
-        }
     }
-
+    }
+    
     @Override
     public List<MeshPoint> getSamples() {
         return Collections.unmodifiableList(allVertices);

Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/sampling/PointSampling.java

@@ -25,8 +25,8 @@ public abstract class PointSampling extends MeshVisitor {
         MAX_VERTICES
     };
     
-    private final PointSamplingType samplingType;
-    private final double samplingLimit;
+    private PointSamplingType samplingType;
+    private double samplingLimit;
     
     /** 
      * Constructor for no point sampling.
@@ -43,11 +43,7 @@ public abstract class PointSampling extends MeshVisitor {
      * @throws IllegalArgumentException if the input parameter is wrong
      */
     public PointSampling(double perc) {
-        if (perc <= 0.0 || perc > 1) {
-            throw new IllegalArgumentException("perc");
-        }
-        this.samplingType = PointSamplingType.PERCENTAGE;
-        this.samplingLimit = perc;
+        setRequiredSamples(perc);
     }
     
     /**
@@ -57,11 +53,7 @@ public abstract class PointSampling extends MeshVisitor {
      * @throws IllegalArgumentException if the input parameter is wrong
      */
     public PointSampling(int max) {
-        if (max <= 0) {
-            throw new IllegalArgumentException("max");
-        }
-        this.samplingType = PointSamplingType.MAX_VERTICES;
-        this.samplingLimit = max;
+        setRequiredSamples(max);
     }
     
     @Override
@@ -75,6 +67,43 @@ public abstract class PointSampling extends MeshVisitor {
      */
     public abstract List<MeshPoint> getSamples();
     
+    /**
+     * If {@code true}, then the returned points samples are points
+     * from the original mesh. Therefore, any change changes the original mesh!
+     * If {@code}, then new points are returned.
+     * 
+     * @return {@code true} if the point samples include points of the original mesh
+     */
+    public boolean isBackedByOrigMesh() {
+        return true;
+    }
+    
+    /**
+     * Changes the number of required samples.
+     * 
+     * @param max Maximal number of vertices. Must be bigger than zero
+     */
+    public final void setRequiredSamples(int max) {
+        if (max <= 0) {
+            throw new IllegalArgumentException("max");
+        }
+        this.samplingType = PointSamplingType.MAX_VERTICES;
+        this.samplingLimit = max;
+    }
+    
+    /**
+     * Changes the number of required samples.
+     * 
+     * @param perc Percentage - a value in (0.0, 1.0&gt;
+     */
+    public final void setRequiredSamples(double perc) {
+        if (perc <= 0.0 || perc > 1) {
+            throw new IllegalArgumentException("perc");
+        }
+        this.samplingType = PointSamplingType.PERCENTAGE;
+        this.samplingLimit = perc;
+    }
+    
     @Override
     public String toString() {
         if (this.samplingType == PointSamplingType.PERCENTAGE) {

Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/sampling/RandomSampling.java

@@ -5,7 +5,10 @@ import cz.fidentis.analyst.mesh.core.MeshPoint;
 import java.util.ArrayList;
 import java.util.Arrays;
 import java.util.Collections;
+import java.util.HashSet;
 import java.util.List;
+import java.util.Random;
+import java.util.Set;
 import java.util.stream.Collectors;
 import java.util.stream.IntStream;
 
@@ -17,6 +20,13 @@ import java.util.stream.IntStream;
 public class RandomSampling extends PointSampling {
     
     private List<MeshPoint> allVertices = new ArrayList<>();
+    
+    /** 
+     * Constructor for no point sampling (100 %).
+     */
+    public RandomSampling() {
+        this(1.0);
+    }
 
     /**
      * Constructor for PERCENTAGE point sampling type.
@@ -43,16 +53,29 @@ public class RandomSampling extends PointSampling {
         if (facet != null) {
             allVertices.addAll(facet.getVertices());
         }
-    }
+    }        
     
     @Override
     public List<MeshPoint> getSamples() {
         int numReducedVertices = getNumDownsampledPoints(allVertices.size());
         
-        if (allVertices.size() == numReducedVertices) {
+        if (allVertices.size() <= numReducedVertices) {
             return Collections.unmodifiableList(allVertices);
         }
         
+        if (numReducedVertices/allVertices.size() <= 0.2) {
+            Random random = new Random();
+            List<MeshPoint> ret = new ArrayList<>(numReducedVertices);
+            Set<Integer> usedIndex = new HashSet<>();
+            while (ret.size() < numReducedVertices) {
+                int rand = random.nextInt(allVertices.size());
+                if (!usedIndex.add(rand)) {
+                    ret.add(allVertices.get(rand));
+                }
+            }
+            return ret;
+        }
+        
         // generate randomly ordered indexes:
         List<Integer> range = IntStream.range(0, allVertices.size()).boxed().collect(Collectors.toCollection(ArrayList::new));
         Collections.shuffle(range);
@@ -69,12 +92,12 @@ public class RandomSampling extends PointSampling {
                     i -> copy.set(i, null)
             );
             return copy.parallelStream().filter(p -> p != null).collect(Collectors.<MeshPoint>toList());
-        }        
+        } 
     }
-    
+
     @Override
     public String toString() {
         return "random " + super.toString();
     }
-    
+
 }

Comparison/src/main/java/cz/fidentis/analyst/visitors/mesh/sampling/UniformSpaceSampling.java

+package cz.fidentis.analyst.visitors.mesh.sampling;
+
+import cz.fidentis.analyst.grid.UniformGrid3d;
+import cz.fidentis.analyst.mesh.core.MeshFacet;
+import cz.fidentis.analyst.mesh.core.MeshPoint;
+import cz.fidentis.analyst.mesh.core.MeshPointImpl;
+import java.util.ArrayList;
+import java.util.Collection;
+import java.util.Collections;
+import java.util.List;
+import java.util.stream.Collectors;
+import javax.vecmath.Point3d;
+
+/**
+ * This downsampling algorithm divides the space uniformly using 3D grid.
+ * Then, it computes and returns average position (centroid) of points in grid cells.
+ * The number of samples is often slightly higher the the required number.
+ * 
+ * @author Radek Oslejsek
+ */
+public class UniformSpaceSampling extends PointSampling {
+    
+    
+    private static final double CELL_RESIZE_INIT_VALUE = 1.0;
+    private static final double CELL_RESIZE_INIT_STEP = 1.0;
+    private static final double CELL_RESIZE_STEP_CHANGE = 2.0;
+    private static final double DOWNSAMPLING_TOLERANCE = 50;
+    private static final double MAX_ITERATIONS = 100;
+    
+    private final List<MeshPoint> allVertices = new ArrayList<>();
+    
+    /** 
+     * Constructor for no point sampling (100 %).
+     */
+    public UniformSpaceSampling() {
+        this(1.0);
+    }
+    
+    /**
+     * Constructor for PERCENTAGE point sampling type.
+     * 
+     * @param perc Percentage - a value in (0.0, 1.0&gt;
+     * @throws IllegalArgumentException if the input parameter is wrong
+     */
+    public UniformSpaceSampling(double perc) {
+        super(perc);
+    }
+    
+    /**
+     * Constructor for MAX_VERTICES point sampling type.
+     * 
+     * @param max Required number of samples. Must be bigger than zero
+     * @throws IllegalArgumentException if the input parameter is wrong
+     */
+    public UniformSpaceSampling(int max) {
+        super(max);
+    }
+    
+    @Override
+    public boolean isBackedByOrigMesh() {
+        return false;
+    }
+     
+    @Override
+    public void visitMeshFacet(MeshFacet facet) {
+        if (facet != null) {
+            allVertices.addAll(facet.getVertices());
+        }
+    }
+    
+    @Override
+    public List<MeshPoint> getSamples() {
+        int numReducedVertices = getNumDownsampledPoints(allVertices.size());
+        if (allVertices.size() == numReducedVertices) {
+            return Collections.unmodifiableList(allVertices);
+        }
+        
+        MeshPoint centroid = new MeshPointImpl(allVertices);
+        double avgDist = getAvgDist(allVertices, centroid.getPosition());
+        
+        double k = CELL_RESIZE_INIT_VALUE;
+        double step = CELL_RESIZE_INIT_STEP;
+        int numCells = 0;
+        UniformGrid3d<MeshPoint> grid;
+        int counter = 0;
+        do {
+            double cellSize = avgDist / k;
+            grid =  new UniformGrid3d<>(cellSize, allVertices, (MeshPoint mp) -> mp.getPosition());
+            numCells = grid.numOccupiedCells();
+            if (step > 0 && numCells > numReducedVertices) {
+                step /= -CELL_RESIZE_STEP_CHANGE;
+            } else if (step < 0 && numCells < numReducedVertices) {
+                step /= -CELL_RESIZE_STEP_CHANGE;
+            } else if (step == 0) {
+                break;
+            }
+            if (counter++ > MAX_ITERATIONS) {
+                break;
+            }
+            k += step;
+        } while (numCells < numReducedVertices || numCells > numReducedVertices + DOWNSAMPLING_TOLERANCE);
+        
+        return grid.getNonEmptyCells().stream()
+                        .map(list -> new MeshPointImpl(list)) // compute controid of cell's points
+                        .collect(Collectors.toList());
+    }
+    
+    @Override
+    public String toString() {
+        return "uniform space " + super.toString();
+    }
+    
+    /**
+     * Computes the average distance between vertices and their centroid.
+     * 
+     * @param vertices vertices
+     * @param centroid their centroid
+     * @return the average distance between vertices and their centroid.
+     */
+    protected double getAvgDist(Collection<MeshPoint> vertices, Point3d centroid) {
+        return vertices.stream()
+                .map(p -> p.getPosition())
+                .mapToDouble(p -> centroid.distance(p))
+                .average()
+                .orElse(0);
+    }
+}
\ No newline at end of file