From 47b39bcb2fe2307f931e6473f2d26d766b4ed6e5 Mon Sep 17 00:00:00 2001
From: Kacper Pluta <kacper.pluta@bsdlabs.com.br>
Date: Tue, 4 Aug 2020 13:28:36 +0300
Subject: [PATCH 1/3] reduce numerical errors

---
 include/hedra/polyhedral_face_normals.h | 11 ++++-------
 1 file changed, 4 insertions(+), 7 deletions(-)

diff --git a/include/hedra/polyhedral_face_normals.h b/include/hedra/polyhedral_face_normals.h
index f9d5734..431e53d 100644
--- a/include/hedra/polyhedral_face_normals.h
+++ b/include/hedra/polyhedral_face_normals.h
@@ -29,21 +29,18 @@ namespace hedra
     {
         using namespace Eigen;
         faceNormals.resize(D.rows(),3);
-        for (int i=0;i<D.rows();i++){
+        for (int i=0;i<D.rows();i++) {
             RowVector3d faceNormal; faceNormal<<0.0,0.0,0.0;
-            for (int j=0;j<D(i);j++){
+            for (int j=0;j<D(i);j++) {
                 RowVector3d vn=V.row(F(i,(j+D(i)-1)%D(i)));
                 RowVector3d v0=V.row(F(i,j));
                 RowVector3d v1=V.row(F(i,(j+1)%D(i)));
-                if (((v1-v0).cross(vn-v0)).norm()>10e-6)
-                    faceNormal=faceNormal+((v1-v0).cross(vn-v0)).normalized();
+                faceNormal += (v1-v0).cross(vn-v0);
             }
-            
             faceNormals.row(i)=faceNormal.normalized();
         }
-        
         return true;
-    }
+   }
 }
 
 

From 9083198ac48ecc5cae611b8d58e261d67ae55c53 Mon Sep 17 00:00:00 2001
From: Kacper Pluta <kacper.pluta@bsdlabs.com.br>
Date: Tue, 24 Nov 2020 19:50:59 +0200
Subject: [PATCH 2/3] add different dualizer

---
 include/hedra/dual_mesh.h | 83 +++++++++++++++++++++++++++++----------
 1 file changed, 63 insertions(+), 20 deletions(-)

diff --git a/include/hedra/dual_mesh.h b/include/hedra/dual_mesh.h
index fcc4dd8..40ed0c8 100644
--- a/include/hedra/dual_mesh.h
+++ b/include/hedra/dual_mesh.h
@@ -9,11 +9,11 @@
 #define HEDRA_DUAL_MESH_H
 #include <igl/igl_inline.h>
 #include <igl/setdiff.h>
-#include <hedra/polygonal_face_centers.h>
-#include <hedra/linear_vi_subdivision.h>
-#include <hedra/moebius_vi_subdivision.h>
-#include <hedra/subdivision_basics.h>
-#include <hedra/dcel.h>
+#include "polygonal_face_centers.h"
+#include "linear_vi_subdivision.h"
+#include "moebius_vi_subdivision.h"
+#include "subdivision_basics.h"
+#include "dcel.h"
 #include <Eigen/Core>
 #include <string>
 #include <vector>
@@ -40,7 +40,8 @@ namespace hedra
                             OneRingSubdivisionData& sd,
                             Eigen::MatrixXd& dualV,
                             Eigen::VectorXi& dualD,
-                            Eigen::MatrixXi& dualF)
+                            Eigen::MatrixXi& dualF,
+                            const bool clipBoundary)
   {
     using namespace Eigen;
     using namespace std;
@@ -76,14 +77,10 @@ namespace hedra
             currVertexinFace=k;
         }
         
-        //std::cout<<"origFaceVertices: "<<origFaceVertices<<std::endl;
         MatrixXd canonFaceVertices;
         canonFaceVertices=sd.original2Canonical(i,origFaceVertices);
-        //std::cout<<"canonFaceVertices: "<<canonFaceVertices<<std::endl;
         canonFacePoints.row(j)=canonFaceVertices.colwise().mean();
-        //std::cout<<"canonFacePoints.row(j): "<<canonFacePoints.row(j)<<std::endl;
         
-        //std::cout<<"sd.canonical2Original(i,canonFacePoints.row(j)): "<<sd.canonical2Original(i,canonFacePoints.row(j))<<std::endl;
         
         //Lifting to candidate points
         candidateFacePoints.block(sd.ringFaces(i,j), 3*currVertexinFace,1,3)=sd.canonical2Original(i,canonFacePoints.row(j));  //not entirely optimized
@@ -91,18 +88,11 @@ namespace hedra
       
       MatrixXd canonEdgePoints = sd.canonicalEdgePoints(i, canonCenter, canonStarVertices, canonFacePoints);
       
-      /*if (isBoundaryVertex(i)){
-       canonEdgePoints.row(0)=(canonCenter+canonStarVertices.row(0))/2.0;
-       canonEdgePoints.row(vertexValences(i)-1)=(canonCenter+canonStarVertices.row(vertexValences(i)-1))/2.0;
-       }*/
-      
       //Lifting to candidate points
       MatrixXd localCandidateEdgePoints;
       localCandidateEdgePoints=sd.canonical2Original(i,canonEdgePoints);
       for (int j=0;j<sd.vertexValences(i);j++){
         int onEdge=sd.HE(sd.starHalfedges(i,j));
-        // std::cout<<"sd.starHalfedges(i,j): "<<sd.starHalfedges(i,j)<<std::endl;
-        // std::cout<<"onEdge: "<<onEdge<<std::endl;
         int inEdge = (sd.EV(onEdge,0)==i ? 0 : 1);
         candidateEdgePoints.block(onEdge,3*inEdge,1,3)=localCandidateEdgePoints.row(j);  //WRONG!!!
       }
@@ -200,6 +190,58 @@ namespace hedra
       for (int j=0;j<dualFList[i].size();j++)
         dualF(i,j)=dualFList[i][j];
     
+    //addendum: killing dual faces (all faces that are a boundary)
+    if (clipBoundary){
+      MatrixXi EV, FE, EF, EFi;
+      MatrixXd FEs;
+      VectorXi innerEdges;
+      hedra::polygonal_edge_topology(dualD, dualF, EV, FE, EF, EFi, FEs, innerEdges);
+      VectorXi isBoundaryFace=VectorXi::Zero(dualF.rows());
+      for (int i=0;i<EF.rows();i++){
+        if (EF(i,0)==-1)
+          isBoundaryFace(EF(i,1))=1;
+        if (EF(i,1)==-1)
+          isBoundaryFace(EF(i,0))=1;
+      }
+      cout<<"Clipping "<<isBoundaryFace.sum()<<" boundary faces"<<endl;
+      MatrixXi newDualF(dualF.rows()-isBoundaryFace.sum(),dualF.cols());
+      VectorXi newDualD(dualF.rows()-isBoundaryFace.sum());
+      int fCount=0;
+      for (int i=0;i<dualF.rows();i++){
+        if (!isBoundaryFace(i)){
+          newDualD(fCount)=dualD(i);
+          newDualF.row(fCount++)=dualF.row(i);
+        }
+      }
+      
+      dualF=newDualF;
+      dualD=newDualD;
+      //removing isolated vertices
+      VectorXi isUnrefVertex=VectorXi::Constant(dualV.rows(),1);
+      for (int i=0;i<dualD.rows();i++)
+        for (int j=0;j<dualD(i);j++)
+          isUnrefVertex(dualF(i,j))=0;
+      
+      int vCount=0;
+      VectorXi transVertices=VectorXi::Constant(dualV.rows(),-1);
+      for (int i=0;i<dualV.rows();i++)
+        if (!isUnrefVertex(i))
+          transVertices(i)=vCount++;
+      
+      for (int i=0;i<dualD.rows();i++)
+        for (int j=0;j<dualD(i);j++)
+          dualF(i,j)=transVertices(dualF(i,j));
+      
+      MatrixXd newDualV(vCount, 3);
+      for (int i=0;i<dualV.rows();i++)
+        if (!isUnrefVertex(i))
+          newDualV.row(transVertices(i))=dualV.row(i);
+      
+      dualV = newDualV;
+      
+    }
+    
+    
     return true;
   }
   
@@ -212,7 +254,8 @@ namespace hedra
                             const int& st,
                             Eigen::MatrixXd& fineV,
                             Eigen::VectorXi& fineD,
-                            Eigen::MatrixXi& fineF)
+                            Eigen::MatrixXi& fineF,
+                            const bool clipBoundary)
   
   {
     using namespace Eigen;
@@ -221,12 +264,12 @@ namespace hedra
     switch (st){
       case hedra::LINEAR_SUBDIVISION: {
         hedra::LinearVISubdivisionData lsd;
-        dual_mesh(V, D,F,lsd, fineV, fineD, fineF);
+        dual_mesh(V, D,F,lsd, fineV, fineD, fineF,clipBoundary);
         break;
       }
       case hedra::CANONICAL_MOEBIUS_SUBDIVISION: {
         hedra::MoebiusVISubdivisionData msd;
-        dual_mesh(V, D,F,msd, fineV, fineD, fineF);
+        dual_mesh(V, D,F,msd, fineV, fineD, fineF,clipBoundary);
         break;
       }
       default: return false;

From ede727e377a30a54cc03b2943934baf22eeb48dd Mon Sep 17 00:00:00 2001
From: Kacper Pluta <kacper.pluta@bsdlabs.com.br>
Date: Tue, 24 Nov 2020 20:46:46 +0200
Subject: [PATCH 3/3] Revert "add different dualizer"

This reverts commit 9083198ac48ecc5cae611b8d58e261d67ae55c53.
---
 include/hedra/dual_mesh.h | 83 ++++++++++-----------------------------
 1 file changed, 20 insertions(+), 63 deletions(-)

diff --git a/include/hedra/dual_mesh.h b/include/hedra/dual_mesh.h
index 40ed0c8..fcc4dd8 100644
--- a/include/hedra/dual_mesh.h
+++ b/include/hedra/dual_mesh.h
@@ -9,11 +9,11 @@
 #define HEDRA_DUAL_MESH_H
 #include <igl/igl_inline.h>
 #include <igl/setdiff.h>
-#include "polygonal_face_centers.h"
-#include "linear_vi_subdivision.h"
-#include "moebius_vi_subdivision.h"
-#include "subdivision_basics.h"
-#include "dcel.h"
+#include <hedra/polygonal_face_centers.h>
+#include <hedra/linear_vi_subdivision.h>
+#include <hedra/moebius_vi_subdivision.h>
+#include <hedra/subdivision_basics.h>
+#include <hedra/dcel.h>
 #include <Eigen/Core>
 #include <string>
 #include <vector>
@@ -40,8 +40,7 @@ namespace hedra
                             OneRingSubdivisionData& sd,
                             Eigen::MatrixXd& dualV,
                             Eigen::VectorXi& dualD,
-                            Eigen::MatrixXi& dualF,
-                            const bool clipBoundary)
+                            Eigen::MatrixXi& dualF)
   {
     using namespace Eigen;
     using namespace std;
@@ -77,10 +76,14 @@ namespace hedra
             currVertexinFace=k;
         }
         
+        //std::cout<<"origFaceVertices: "<<origFaceVertices<<std::endl;
         MatrixXd canonFaceVertices;
         canonFaceVertices=sd.original2Canonical(i,origFaceVertices);
+        //std::cout<<"canonFaceVertices: "<<canonFaceVertices<<std::endl;
         canonFacePoints.row(j)=canonFaceVertices.colwise().mean();
+        //std::cout<<"canonFacePoints.row(j): "<<canonFacePoints.row(j)<<std::endl;
         
+        //std::cout<<"sd.canonical2Original(i,canonFacePoints.row(j)): "<<sd.canonical2Original(i,canonFacePoints.row(j))<<std::endl;
         
         //Lifting to candidate points
         candidateFacePoints.block(sd.ringFaces(i,j), 3*currVertexinFace,1,3)=sd.canonical2Original(i,canonFacePoints.row(j));  //not entirely optimized
@@ -88,11 +91,18 @@ namespace hedra
       
       MatrixXd canonEdgePoints = sd.canonicalEdgePoints(i, canonCenter, canonStarVertices, canonFacePoints);
       
+      /*if (isBoundaryVertex(i)){
+       canonEdgePoints.row(0)=(canonCenter+canonStarVertices.row(0))/2.0;
+       canonEdgePoints.row(vertexValences(i)-1)=(canonCenter+canonStarVertices.row(vertexValences(i)-1))/2.0;
+       }*/
+      
       //Lifting to candidate points
       MatrixXd localCandidateEdgePoints;
       localCandidateEdgePoints=sd.canonical2Original(i,canonEdgePoints);
       for (int j=0;j<sd.vertexValences(i);j++){
         int onEdge=sd.HE(sd.starHalfedges(i,j));
+        // std::cout<<"sd.starHalfedges(i,j): "<<sd.starHalfedges(i,j)<<std::endl;
+        // std::cout<<"onEdge: "<<onEdge<<std::endl;
         int inEdge = (sd.EV(onEdge,0)==i ? 0 : 1);
         candidateEdgePoints.block(onEdge,3*inEdge,1,3)=localCandidateEdgePoints.row(j);  //WRONG!!!
       }
@@ -190,58 +200,6 @@ namespace hedra
       for (int j=0;j<dualFList[i].size();j++)
         dualF(i,j)=dualFList[i][j];
     
-    //addendum: killing dual faces (all faces that are a boundary)
-    if (clipBoundary){
-      MatrixXi EV, FE, EF, EFi;
-      MatrixXd FEs;
-      VectorXi innerEdges;
-      hedra::polygonal_edge_topology(dualD, dualF, EV, FE, EF, EFi, FEs, innerEdges);
-      VectorXi isBoundaryFace=VectorXi::Zero(dualF.rows());
-      for (int i=0;i<EF.rows();i++){
-        if (EF(i,0)==-1)
-          isBoundaryFace(EF(i,1))=1;
-        if (EF(i,1)==-1)
-          isBoundaryFace(EF(i,0))=1;
-      }
-      cout<<"Clipping "<<isBoundaryFace.sum()<<" boundary faces"<<endl;
-      MatrixXi newDualF(dualF.rows()-isBoundaryFace.sum(),dualF.cols());
-      VectorXi newDualD(dualF.rows()-isBoundaryFace.sum());
-      int fCount=0;
-      for (int i=0;i<dualF.rows();i++){
-        if (!isBoundaryFace(i)){
-          newDualD(fCount)=dualD(i);
-          newDualF.row(fCount++)=dualF.row(i);
-        }
-      }
-      
-      dualF=newDualF;
-      dualD=newDualD;
-      //removing isolated vertices
-      VectorXi isUnrefVertex=VectorXi::Constant(dualV.rows(),1);
-      for (int i=0;i<dualD.rows();i++)
-        for (int j=0;j<dualD(i);j++)
-          isUnrefVertex(dualF(i,j))=0;
-      
-      int vCount=0;
-      VectorXi transVertices=VectorXi::Constant(dualV.rows(),-1);
-      for (int i=0;i<dualV.rows();i++)
-        if (!isUnrefVertex(i))
-          transVertices(i)=vCount++;
-      
-      for (int i=0;i<dualD.rows();i++)
-        for (int j=0;j<dualD(i);j++)
-          dualF(i,j)=transVertices(dualF(i,j));
-      
-      MatrixXd newDualV(vCount, 3);
-      for (int i=0;i<dualV.rows();i++)
-        if (!isUnrefVertex(i))
-          newDualV.row(transVertices(i))=dualV.row(i);
-      
-      dualV = newDualV;
-      
-    }
-    
-    
     return true;
   }
   
@@ -254,8 +212,7 @@ namespace hedra
                             const int& st,
                             Eigen::MatrixXd& fineV,
                             Eigen::VectorXi& fineD,
-                            Eigen::MatrixXi& fineF,
-                            const bool clipBoundary)
+                            Eigen::MatrixXi& fineF)
   
   {
     using namespace Eigen;
@@ -264,12 +221,12 @@ namespace hedra
     switch (st){
       case hedra::LINEAR_SUBDIVISION: {
         hedra::LinearVISubdivisionData lsd;
-        dual_mesh(V, D,F,lsd, fineV, fineD, fineF,clipBoundary);
+        dual_mesh(V, D,F,lsd, fineV, fineD, fineF);
         break;
       }
       case hedra::CANONICAL_MOEBIUS_SUBDIVISION: {
         hedra::MoebiusVISubdivisionData msd;
-        dual_mesh(V, D,F,msd, fineV, fineD, fineF,clipBoundary);
+        dual_mesh(V, D,F,msd, fineV, fineD, fineF);
         break;
       }
       default: return false;