diff --git a/Comparison/src/main/java/cz/fidentis/analyst/symmetry/ApproxSymmetryPlane.java b/Comparison/src/main/java/cz/fidentis/analyst/symmetry/ApproxSymmetryPlane.java
index 53304801656eb1ac98413f70db804bebf0b34100..48228ce7c7f06b0495946cb1ce94882685034250 100644
--- a/Comparison/src/main/java/cz/fidentis/analyst/symmetry/ApproxSymmetryPlane.java
+++ b/Comparison/src/main/java/cz/fidentis/analyst/symmetry/ApproxSymmetryPlane.java
@@ -49,6 +49,7 @@ public class ApproxSymmetryPlane extends Plane implements Comparable<ApproxSymme
         MeshPoint meshPointI = facetI.getVertex(sigPointI);
         MeshPoint meshPointJ = facetJ.getVertex(sigPointJ);
                 
+        // accpet only point pairs with significantly different curvatures (WTF?)
         double minRatio = config.getMinCurvRatio();
         double maxRatio = 1.0 / minRatio;
         double ratioIJ = sigPoints.getCachedCurRatio(i, j);
@@ -62,16 +63,14 @@ public class ApproxSymmetryPlane extends Plane implements Comparable<ApproxSymme
         normal.sub(p2);
         normal.normalize(); 
         
+        // accpect only point pair with oposite normals along with the plane normal:
         double normCos = sigPoints.getCachedNormCosVec(i, j).dot(normal);
         if (Math.abs(normCos) < config.getMinNormAngleCos()) {
             throw new UnsupportedOperationException();
         }
         
-        Vector3d avrg = sigPoints.getCachedAvgPos(i, j);
-        double d = -(normal.x * avrg.x) - (normal.y * avrg.y) - (normal.z * avrg.z);
-
+        setDistance(-normal.dot(sigPoints.getCachedAvgPos(i, j))); 
         setNormal(normal);
-        setDistance(d);
     }
 
     /**
@@ -115,8 +114,8 @@ public class ApproxSymmetryPlane extends Plane implements Comparable<ApproxSymme
                     continue;
                 }
                 
-                Vector3d avrg = sigPoints.getCachedAvgPos(i, j);
-                double dist = Math.abs(normal.x * avrg.x + normal.y * avrg.y + normal.z * avrg.z + d);
+                Vector3d avg = sigPoints.getCachedAvgPos(i, j);
+                double dist = Math.abs(normal.dot(avg) + d);
                 if (dist <= maxDist) {                 
                     votes++;  
                 }   
diff --git a/Comparison/src/main/java/cz/fidentis/analyst/symmetry/Plane.java b/Comparison/src/main/java/cz/fidentis/analyst/symmetry/Plane.java
index 4516b41eec3f7160e1abc635bb502085860e4f51..80737435ff3459dbd291c509b212babfd5f0fa8f 100644
--- a/Comparison/src/main/java/cz/fidentis/analyst/symmetry/Plane.java
+++ b/Comparison/src/main/java/cz/fidentis/analyst/symmetry/Plane.java
@@ -17,28 +17,18 @@ public class Plane implements Serializable {
    
     private Vector3d normal;
     private double   distance;
+    private double   normSquare; // normal.dot(normal)
 
-    /**
-     * Creates new plane.
-     * 
-     * @param a a coordinate 
-     * @param b b coordinate
-     * @param c c coordinate
-     * @param d d coordinate
-     */
-    public Plane(double a, double b, double c, double d) {
-        normal = new Vector3d(a,b,c);
-        distance = d;
-    }
-    
     /**
      * Constructor.
-     * @param normal Normal vector of the plane
+     * 
+     * @param normal Normalized (!) normal vector of the plane
      * @param dist distance
      * @throws IllegalArgumentException if the @code{plane} argument is null
      */
     public Plane(Vector3d normal, double dist) {
-        this(normal.x, normal.y, normal.z, dist);
+        setNormal(new Vector3d(normal));
+        setDistance(dist);
     }
     
     /**
@@ -86,10 +76,11 @@ public class Plane implements Serializable {
     /**
      * Normalize the plane
      */
-    public void normalize() {
+    protected final void normalize() {
         double normalLength = normal.length();
         normal.normalize();
         distance /= normalLength; // Do we really want this? --ro
+        this.normSquare = normal.dot(normal);
     }
     
     /**
@@ -123,50 +114,74 @@ public class Plane implements Serializable {
         this.distance += value;
     }
     
+    /**
+     * Returns a point laying on the plane.
+     * 
+     * @param point A 3D point, which is projected to the plane
+     * @return a point laying on the plane.
+     * @throws NullPointerException if the {@code point} is {@code null}
+     */
+    public Point3d projectToPlane(Point3d point) {
+        double shiftDist = 
+                ((normal.x * point.x) + 
+                  (normal.y * point.y) + 
+                  (normal.z * point.z) +
+                  distance) / normSquare;
+        
+        Point3d ret = new Point3d(normal);
+        ret.scale(-shiftDist);
+        ret.add(point);
+        return ret;
+    }
+    
+    /**
+     * Returns a point laying on the opposite side of the plane ("mirrors" the point).
+     * This method implements equation (1) of 
+     * <a href="https://link.springer.com/content/pdf/10.1007/s00371-020-02034-w.pdf">Hruda et al: Robust, fast and flexible symmetry plane detection based
+on differentiable symmetry measure</a>
+     * 
+     * @param point A 3D point to be reflected
+     * @return a point on the opposite side of the plane.
+     * @throws NullPointerException if the {@code point} is {@code null}
+     */
+    public Point3d reflectOverPlane(Point3d point) {
+        double shiftDist = 
+                ((normal.x * point.x) + 
+                  (normal.y * point.y) + 
+                  (normal.z * point.z) +
+                  distance) / normSquare;
+        
+        Point3d ret = new Point3d(normal);
+        ret.scale(-2.0 * shiftDist);
+        ret.add(point);
+        return ret;
+    }
+    
     /**
      * Returns rectangular mesh facet of this symmetry plane.
      * 
-     * @param midPoint A 3D point, which is projected to the plane and used as a center of the rectangular facet
+     * @param point A 3D point, which is projected to the plane and used as a center of the rectangular facet
      * @param width Width
      * @param height Height
      * @return a rectangular mesh facet of this symmetry plane
-     * @throws NullPointerException if the {@code midPoint} is {@code null}
+     * @throws NullPointerException if the {@code point} is {@code null}
      * @throws IllegalArgumentException if {@code width} or {@code height} are &lt;= 0
      */
-    public MeshRectangleFacet getMesh(Point3d midPoint, double width, double height) {
-        if (width <= 0) {
-            throw new IllegalArgumentException("width");
-        }
-        
-        if (height <= 0) {
-            throw new IllegalArgumentException("height");
-        }
-        
-        double alpha = 
-                -((normal.x * midPoint.x) + 
-                  (normal.y * midPoint.y) + 
-                  (normal.z * midPoint.z) +
-                  distance) / (normal.dot(normal));
-        
-        Point3d midPointOnPlane = new Point3d(midPoint);
-        Vector3d nn = new Vector3d(normal);
-        nn.scale(alpha);
-        midPointOnPlane.add(nn);
-        
-        return new MeshRectangleFacet(midPointOnPlane, normal, width, height);
+    public MeshRectangleFacet getMesh(Point3d point, double width, double height) {
+        return new MeshRectangleFacet(projectToPlane(point), normal, width, height);
     }
     
     /**
      * Returns rectangular mesh facet of this symmetry plane.
      * 
-     * @param midPoint A 3D point, which is projected to the plane and used as a center of the rectangular facet
+     * @param point A 3D point, which is projected to the plane and used as a center of the rectangular facet
      * @param size Edge length
      * @return a rectangular mesh facet of this symmetry plane
-     * @throws NullPointerException if the {@code midPoint} is {@code null}
+     * @throws NullPointerException if the {@code point} is {@code null}
      * @throws IllegalArgumentException if {@code size} is &lt;= 0
      */
-    public MeshRectangleFacet getMesh(Point3d midPoint, double size) {
-        return Plane.this.getMesh(midPoint, size, size);
+    public MeshRectangleFacet getMesh(Point3d point, double size) {
+        return Plane.this.getMesh(point, size, size);
     }
 
     /**
@@ -181,11 +196,12 @@ public class Plane implements Serializable {
         return Plane.this.getMesh(bbox.getMidPoint(), bbox.getDiagonalLength(), bbox.getDiagonalLength());
     }
 
-    protected void setNormal(Vector3d normal) {
+    protected final void setNormal(Vector3d normal) {
         this.normal = normal;
+        this.normSquare = normal.dot(normal);
     }
     
-    protected void setDistance(double dist) {
+    protected final void setDistance(double dist) {
         this.distance = dist;
     }
     
@@ -197,7 +213,7 @@ public class Plane implements Serializable {
      */
     public double getPointDistance(Point3d point) {
         return ((normal.x * point.x) + (normal.y * point.y) + (normal.z * point.z) + distance)
-                / Math.sqrt(normal.x * normal.x + normal.y * normal.y + normal.z * normal.z);
+                / Math.sqrt(normSquare);
     }
 
     /**
diff --git a/Comparison/src/main/java/cz/fidentis/analyst/symmetry/SymmetryEstimator.java b/Comparison/src/main/java/cz/fidentis/analyst/symmetry/SymmetryEstimator.java
index 993c2996ec49e4f4f7bfee9da9b63bd7a1dd05e4..7431560852e8ee3f643b0df08a16a4484cefc801 100644
--- a/Comparison/src/main/java/cz/fidentis/analyst/symmetry/SymmetryEstimator.java
+++ b/Comparison/src/main/java/cz/fidentis/analyst/symmetry/SymmetryEstimator.java
@@ -12,7 +12,6 @@ import java.util.concurrent.ExecutorService;
 import java.util.concurrent.Executors;
 import java.util.concurrent.Future;
 import java.util.logging.Level;
-import java.util.logging.Logger;
 
 /**
  * Main class for computing approximate plane of symmetry of the 3D model.
@@ -23,6 +22,10 @@ import java.util.logging.Logger;
  * 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>
  * 
  * @author Natalia Bebjakova
  * @author Radek Oslejsek
@@ -169,7 +172,7 @@ public class SymmetryEstimator extends MeshVisitor {
                 }
             }
         } catch (final InterruptedException | ExecutionException ex) {
-            Logger.getLogger(SymmetryEstimator.class.getName()).log(Level.SEVERE, null, ex);
+            java.util.logging.Logger.getLogger(SymmetryEstimator.class.getName()).log(Level.SEVERE, null, ex);
         }
         
         setSymmetryPlane(planes);
diff --git a/Comparison/src/main/java/cz/fidentis/analyst/visitors/face/HausdorffDistancePrioritized.java b/Comparison/src/main/java/cz/fidentis/analyst/visitors/face/HausdorffDistancePrioritized.java
index c15d2e62d3eb8a94676bb52e3d8995056726bf6c..98850ff527713a1112961f2c3b45664fc8934fd3 100644
--- a/Comparison/src/main/java/cz/fidentis/analyst/visitors/face/HausdorffDistancePrioritized.java
+++ b/Comparison/src/main/java/cz/fidentis/analyst/visitors/face/HausdorffDistancePrioritized.java
@@ -1,6 +1,5 @@
 package cz.fidentis.analyst.visitors.face;
 
-import cz.fidentis.analyst.Logger;
 import cz.fidentis.analyst.face.HumanFace;
 import cz.fidentis.analyst.feature.FeaturePoint;
 import cz.fidentis.analyst.feature.FeaturePointType;
diff --git a/GUI/src/main/java/cz/fidentis/analyst/symmetry/SymmetryAction.java b/GUI/src/main/java/cz/fidentis/analyst/symmetry/SymmetryAction.java
index 6de1ad7c66eb7d6fb0ad9e6181ff73900606ff59..81599d7b530d280517d00b9738667d2a1651909f 100644
--- a/GUI/src/main/java/cz/fidentis/analyst/symmetry/SymmetryAction.java
+++ b/GUI/src/main/java/cz/fidentis/analyst/symmetry/SymmetryAction.java
@@ -17,7 +17,9 @@ import cz.fidentis.analyst.visitors.mesh.HausdorffDistance;
 import cz.fidentis.analyst.visitors.mesh.HausdorffDistance.Strategy;
 
 import java.awt.event.ActionEvent;
+import java.util.ArrayList;
 import java.util.HashMap;
+import java.util.List;
 import java.util.Map;
 import javax.swing.JTabbedPane;
 import javax.vecmath.Point3d;
@@ -126,62 +128,41 @@ public class SymmetryAction extends ControlPanelAction implements HumanFaceListe
             return;
         }
         
-        //FeaturePointTypeProvider.getInstance();
         Map<String, Point3d> processFP = new HashMap<>();
-        Vector3d dir = new Vector3d();
-        int numDirs = 0;
-        Point3d centroid = new Point3d();
+        List<Plane> planes = new ArrayList<>();
         
         for (FeaturePoint fp: face.getFeaturePoints()) {
-            Point3d auxCentroid = new Point3d(fp.getPosition());
-                auxCentroid.scale(1.0 / face.getFeaturePoints().size());
-            if (numDirs == 0) {
-                centroid.set(auxCentroid);
-            } else {
-                
-                centroid.add(auxCentroid);
-            }
-
             String code = fp.getFeaturePointType().getCode();
+            String pairCode;
             if (code.endsWith("_R")) {
-                String pairCode = code.replaceAll("_R", "_L");
-                if (processFP.containsKey(pairCode)) {
-                    computeDir(processFP.remove(pairCode), fp.getPosition(), dir, numDirs++);
-                } else {
-                    processFP.put(code, fp.getPosition());
-                }
+                pairCode = code.replaceAll("_R", "_L");
             } else if (code.endsWith("_L")) {
-                String pairCode = code.replaceAll("_L", "_R");
-                if (processFP.containsKey(pairCode)) {
-                    computeDir(fp.getPosition(), processFP.remove(pairCode), dir, numDirs++);
-                } else {
-                    processFP.put(code, fp.getPosition());
-                }                
+                pairCode = code.replaceAll("_L", "_R");
+            } else {
+                continue;
             }
+            
+            if (!processFP.containsKey(pairCode)) {
+                processFP.put(code, fp.getPosition());
+                continue;
+            }
+            
+            Vector3d normal = new Vector3d(processFP.remove(pairCode));
+            Vector3d center = new Vector3d(normal);
+                
+            normal.sub(fp.getPosition());
+            normal.normalize();
+                    
+            center.add(fp.getPosition());
+            center.scale(0.5);
+                    
+            planes.add(new Plane(normal, normal.dot(center)));
         }
         
-        if (numDirs == 0) { // no pair found
-            return;
-        }
-        
-        dir.scale(1.0 / numDirs);
-        dir.normalize();
-        
-        face.setSymmetryPlane(new Plane(dir, dir.dot(new Vector3d(centroid))));
+        face.setSymmetryPlane(new Plane(planes)); // creates an average plane
         setDrawablePlane(face, index);
     }
     
-    protected void computeDir(Point3d left, Point3d right, Vector3d dir, int numDirs) {
-        Vector3d auxDir = new Vector3d(right);
-        auxDir.sub(left);
-        auxDir.normalize();
-        if (numDirs == 0) {
-            dir.set(auxDir);
-        } else {
-            dir.add(auxDir);
-        }
-    }
-
     protected void setDrawablePlane(HumanFace face, int index) {
         DrawablePlane symmetryPlane = new DrawablePlane(
                 face.getSymmetryPlaneFacet(),
@@ -199,13 +180,9 @@ public class SymmetryAction extends ControlPanelAction implements HumanFaceListe
         
         MeshModel clone = new MeshModel(face.getMeshModel());
         
-        final Vector3d normal = face.getSymmetryPlane().getNormal();
         for (MeshFacet facet: clone.getFacets()) {
             facet.getVertices().parallelStream().forEach(v -> {
-                double d = normal.dot(new Vector3d(v.getPosition()));
-                Vector3d shiftDir = new Vector3d(normal);
-                shiftDir.scale(-2.0 * d);
-                v.getPosition().add(shiftDir);
+                v.getPosition().set(face.getSymmetryPlane().reflectOverPlane(v.getPosition()));
             });
         }
         
diff --git a/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshModel.java b/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshModel.java
index ae5462997554ea1a954b513f1150529efe57e9dc..9d4ee4fe39f73445ec80a8b6a0b1137e4d1cd171 100644
--- a/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshModel.java
+++ b/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshModel.java
@@ -153,6 +153,9 @@ public class MeshModel implements Serializable {
      * @return Number of vertices
      */
     public long getNumVertices() {
-        return facets.stream().map(f -> f.getNumberOfVertices()).reduce(0, Integer::sum);
-    }
+        return facets
+                .stream()
+                .mapToInt(f -> f.getNumberOfVertices())
+                .sum();
+    }    
 }
diff --git a/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshRectangleFacet.java b/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshRectangleFacet.java
index cd5fc6d7ea06642d1443891868c27f301eb1d210..def63ef6de80f95b7bfb1d00a7faa737a107ddea 100644
--- a/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshRectangleFacet.java
+++ b/MeshModel/src/main/java/cz/fidentis/analyst/mesh/core/MeshRectangleFacet.java
@@ -15,10 +15,19 @@ public class MeshRectangleFacet extends MeshFacetImpl {
      * 
      * @param center Central point of the rectangle
      * @param normal Normalized normal vector 
-     * @param w Width
-     * @param h Height
+     * @param width Width
+     * @param height Height
+     * @throws IllegalArgumentException if {@code width} or {@code height} are &lt;= 0
      */
-    public MeshRectangleFacet(Point3d center, Vector3d normal, double w, double h) {
+    public MeshRectangleFacet(Point3d center, Vector3d normal, double width, double height) {
+        if (width <= 0) {
+            throw new IllegalArgumentException("width");
+        }
+        
+        if (height <= 0) {
+            throw new IllegalArgumentException("height");
+        }
+        
         Vector3d hDir = new Vector3d();
         if (Math.abs(normal.dot(new Vector3d(0.0, 1.0, 0.0))) > Math.abs(normal.dot(new Vector3d (1.0, 0.0, 0.0)))) {
             hDir.cross(normal, new Vector3d(1.0, 0.0, 0.0));
@@ -31,7 +40,7 @@ public class MeshRectangleFacet extends MeshFacetImpl {
         vDir.cross(normal, hDir);
         vDir.normalize();
         
-        initRectangle(center, hDir, vDir, normal, w, h);
+        initRectangle(center, hDir, vDir, normal, width, height);
     }
     
     /**