From 83aee42b0a32e7a6b926b4575fd4a547cf6d256d Mon Sep 17 00:00:00 2001
From: nico <nico.pietroni@gmail.com>
Date: Thu, 28 Dec 2017 15:20:42 +1100
Subject: [PATCH] Several changes and one bug solved

---
 vcg/complex/algorithms/polygonal_algorithms.h | 145 ++++++++++++++----
 1 file changed, 114 insertions(+), 31 deletions(-)

diff --git a/vcg/complex/algorithms/polygonal_algorithms.h b/vcg/complex/algorithms/polygonal_algorithms.h
index 8c2f8b5b..7a7593ed 100644
--- a/vcg/complex/algorithms/polygonal_algorithms.h
+++ b/vcg/complex/algorithms/polygonal_algorithms.h
@@ -29,7 +29,7 @@
 #include <vcg/complex/algorithms/update/color.h>
 #include <vcg/complex/algorithms/closest.h>
 #include <vcg/complex/algorithms/point_sampling.h>
-
+#include <vcg/complex/algorithms/update/quality.h>
 #include <wrap/io_trimesh/export_obj.h>
 
 //define a temporary triangle mesh type
@@ -362,14 +362,14 @@ public:
         //CoordType AVN(0,0,0);
         //CoordType AVN0(0,0,0);
         CoordType Origin(0,0,0);
-//        for (int j=0;j<f.VN();j++)
-//            AVN0=AVN0+f.V(j)->N();
+        //        for (int j=0;j<f.VN();j++)
+        //            AVN0=AVN0+f.V(j)->N();
 
         CoordType AVN=vcg::PolygonNormal(f);
         //AVN0.Normalize();
-//        std::cout<<"AVN "<<AVN.X()<<","<<AVN.Y()<<","<<AVN.Z()<<std::endl;
-//        std::cout<<"AVN0 "<<AVN0.X()<<","<<AVN0.Y()<<","<<AVN0.Z()<<std::endl;
-//        std::cout<<"NormT "<<NormT.X()<<","<<NormT.Y()<<","<<NormT.Z()<<std::endl;
+        //        std::cout<<"AVN "<<AVN.X()<<","<<AVN.Y()<<","<<AVN.Z()<<std::endl;
+        //        std::cout<<"AVN0 "<<AVN0.X()<<","<<AVN0.Y()<<","<<AVN0.Z()<<std::endl;
+        //        std::cout<<"NormT "<<NormT.X()<<","<<NormT.Y()<<","<<NormT.Z()<<std::endl;
 
         for (size_t j=0;j<TemplatePos.size();j++)
             Origin+=TemplatePos[j];
@@ -397,7 +397,8 @@ public:
                           bool FixS=false,
                           bool isotropic=true,
                           ScalarType smoothTerm=0.1,
-                          bool fixB=true)
+                          bool fixB=true,
+                          bool WeightByQuality=false)
     {
         (void)isotropic;
         typedef typename PolyMeshType::FaceType PolygonType;
@@ -413,6 +414,9 @@ public:
 
         PolyMeshType TestM;
 
+        if (WeightByQuality)
+            UpdateQuality(poly_m,QTemplate);
+
         for (size_t s=0;s<(size_t)relax_step;s++)
         {
             //initialize the accumulation vector
@@ -428,7 +432,10 @@ public:
                 ScalarType val=vcg::PolyArea(poly_m.face[i]);
                 if (val<(AvgArea*0.00001))
                     val=(AvgArea*0.00001);
+
                 ScalarType W=1.0/val;
+                if (WeightByQuality)
+                    W=poly_m.face[i].Q()+0.00001;
 
                 for (size_t j=0;j<TemplatePos.size();j++)
                 {
@@ -627,13 +634,63 @@ public:
                                    TriMeshType &tri_mesh,
                                    int relaxStep=100,
                                    bool fixIrr=false,
-                                   ScalarType Damp=0.5)
+                                   ScalarType Damp=0.5,
+                                   ScalarType SharpDeg=0,
+                                   bool WeightByQuality=false)
     {
+        vcg::tri::UpdateFlags<PolyMeshType>::VertexClearS(poly_m);
+
         vcg::tri::UpdateTopology<PolyMeshType>::FaceFace(poly_m);
 
         //UpdateBorderVertexFromPFFAdj(poly_m);
         vcg::tri::UpdateFlags<PolyMeshType>::VertexBorderFromFaceAdj(poly_m);
 
+        std::vector<std::vector<vcg::Line3<ScalarType> > > SharpEdge(poly_m.vert.size());
+        //first select sharp features
+        if (SharpDeg>0)
+        {
+            for (size_t i=0;i<poly_m.face.size();i++)
+                for (size_t j=0;j<poly_m.face[i].VN();j++)
+                {
+                    //check only one side
+                    if ((&poly_m.face[i])>=poly_m.face[i].FFp(j))continue;
+
+                    CoordType N0=poly_m.face[i].N();
+                    CoordType N1=poly_m.face[i].FFp(j)->N();
+
+                    ScalarType Angle=vcg::Angle(N0,N1);
+                    if (fabs(Angle)>(SharpDeg* (M_PI / 180.0)))
+                    {
+                        CoordType Pos0=poly_m.face[i].V0(j)->P();
+                        CoordType Pos1=poly_m.face[i].V1(j)->P();
+                        CoordType Ori=Pos0;
+                        CoordType Dir=Pos1-Pos0;
+                        Dir.Normalize();
+                        vcg::Line3<ScalarType> L(Ori,Dir);
+                        int Index0=vcg::tri::Index(poly_m,poly_m.face[i].V0(j));
+                        int Index1=vcg::tri::Index(poly_m,poly_m.face[i].V1(j));
+                        SharpEdge[Index0].push_back(L);
+                        SharpEdge[Index1].push_back(L);
+                    }
+                }
+            for (size_t i=0;i<poly_m.vert.size();i++)
+            {
+                if (SharpEdge[i].size()==0)continue;
+                if (SharpEdge[i].size()>2)poly_m.vert[i].SetS();
+            }
+        }
+        if (fixIrr)
+        {
+            vcg::tri::UpdateQuality<PolyMeshType>::VertexValence(poly_m);
+            for (size_t i=0;i<poly_m.vert.size();i++)
+            {
+                if (poly_m.vert[i].IsB())continue;
+                if (poly_m.vert[i].Q()==4)continue;
+                poly_m.vert[i].SetS();
+            }
+        }
+
+
         typedef typename TriMeshType::FaceType FaceType;
         typedef vcg::GridStaticPtr<FaceType, typename TriMeshType::ScalarType> TriMeshGrid;
         TriMeshGrid grid;
@@ -656,13 +713,13 @@ public:
             typename TriMeshType::FaceType *f=NULL;
             typename TriMeshType::CoordType norm,ip;
             f=vcg::tri::GetClosestFaceBase(tri_mesh,grid,testPos,MaxD,minDist,closestPt,norm,ip);
-            poly_m.vert[i].N().Import(norm);
+            //poly_m.vert[i].N().Import(norm);
         }
 
         for(int k=0;k<relaxStep;k++)
         {
             //smooth PCA step
-            SmoothPCA(poly_m,1,Damp,fixIrr);
+            SmoothPCA(poly_m,1,Damp,true,true,0.1,true,WeightByQuality);
             //reprojection step
             //laplacian smooth step
             //Laplacian(poly_m,Damp,1);
@@ -673,13 +730,48 @@ public:
                 testPos.Import(poly_m.vert[i].P());
                 typename TriMeshType::CoordType closestPt;
                 typename TriMeshType::ScalarType minDist;
-                FaceType *f=NULL;
-                typename TriMeshType::CoordType norm,ip;
-                f=vcg::tri::GetClosestFaceBase(tri_mesh,grid,testPos,MaxD,minDist,closestPt,norm,ip);
-                poly_m.vert[i].P().Import(testPos*Damp+closestPt*(1-Damp));
-                poly_m.vert[i].N().Import(norm);
+                if (SharpEdge[i].size()==0)//reproject onto original mesh
+                {
+                    FaceType *f=NULL;
+                    typename TriMeshType::CoordType norm,ip;
+                    f=vcg::tri::GetClosestFaceBase(tri_mesh,grid,testPos,MaxD,minDist,closestPt,norm,ip);
+                    poly_m.vert[i].P().Import(testPos*Damp+closestPt*(1-Damp));
+                    //poly_m.vert[i].N().Import(norm);
+                }
+                else //reproject onto segments
+                {
+                    CoordType av_closest(0,0,0);
+                    size_t sum=0;
+                    for (size_t j=0;j<SharpEdge[i].size();j++)
+                    {
+                        CoordType closest;
+                        ScalarType dist;
+                        vcg::LinePointDistance(SharpEdge[i][j],testPos,closest,dist);
+                        av_closest+=closest;
+                        sum++;
+                    }
+                    assert(sum>0);
+                    poly_m.vert[i].P()=av_closest/sum;
+                }
             }
+            UpdateFaceNormals(poly_m);
+            vcg::tri::UpdateNormal<PolyMeshType>::PerVertexFromCurrentFaceNormal(poly_m);
         }
+
+    }
+
+
+    template <class TriMeshType>
+    static void TriangulateToTriMesh(PolyMeshType &poly_m,TriMeshType &triangle_mesh)
+    {
+        triangle_mesh.Clear();
+
+        PolyMeshType PolySwap;
+        vcg::tri::Append<PolyMeshType,PolyMeshType>::Mesh(PolySwap,poly_m);
+        Triangulate(PolySwap);
+
+        //then copy onto the triangle mesh
+        vcg::tri::Append<TriMeshType,PolyMeshType>::Mesh(triangle_mesh,PolySwap);
     }
 
     /*! \brief This function performs the polygon regularization as in "Statics Aware Grid Shells"
@@ -688,18 +780,21 @@ public:
     static void SmoothReprojectPCA(PolyMeshType &poly_m,
                                    int relaxStep=100,
                                    bool fixIrr=false,
-                                   ScalarType Damp=0.5)
+                                   ScalarType Damp=0.5,
+                                   ScalarType SharpDeg=0,
+                                   bool WeightByQuality=false)
     {
         //transform into triangular
         TempMesh GuideSurf;
-        vcg::tri::PolygonSupport<TempMesh,PolyMeshType>::ImportFromPolyMesh(GuideSurf,poly_m);
+        //vcg::tri::PolygonSupport<TempMesh,PolyMeshType>:(GuideSurf,poly_m);
+        TriangulateToTriMesh<TempMesh>(poly_m,GuideSurf);
         vcg::tri::UpdateBounding<TempMesh>::Box(GuideSurf);
         vcg::tri::UpdateNormal<TempMesh>::PerVertexNormalizedPerFace(GuideSurf);
         vcg::tri::UpdateTopology<TempMesh>::FaceFace(GuideSurf);
         vcg::tri::UpdateFlags<TempMesh>::FaceBorderFromFF(GuideSurf);
 
         //optimize it
-        vcg::PolygonalAlgorithm<PolyMeshType>::SmoothReprojectPCA<TempMesh>(poly_m,GuideSurf,relaxStep,fixIrr,Damp);
+        vcg::PolygonalAlgorithm<PolyMeshType>::SmoothReprojectPCA<TempMesh>(poly_m,GuideSurf,relaxStep,fixIrr,Damp,SharpDeg,WeightByQuality);
     }
 
     /*! \brief This function return average edge size
@@ -922,7 +1017,7 @@ public:
 
         for (size_t i=0;i<poly_m.face.size();i++)
         {
-//            ScalarType AreaF=vcg::PolyArea(poly_m.face[i]);
+            //            ScalarType AreaF=vcg::PolyArea(poly_m.face[i]);
             size_t sizeV=poly_m.face[i].VN()-1;
             CoordType baryF=vcg::PolyBarycenter(poly_m.face[i]);
             for (size_t j=0;j<poly_m.face[i].VN();j++)
@@ -1027,18 +1122,6 @@ public:
             Triangulate(poly_m,i);
     }
 
-    template <class TriMeshType>
-    static void TriangulateToTriMesh(PolyMeshType &poly_m,TriMeshType &triangle_mesh)
-    {
-        triangle_mesh.Clear();
-
-        PolyMeshType PolySwap;
-        vcg::tri::Append<PolyMeshType,PolyMeshType>::Mesh(PolySwap,poly_m);
-        Triangulate(PolySwap);
-
-        //then copy onto the triangle mesh
-        vcg::tri::Append<TriMeshType,PolyMeshType>::Mesh(triangle_mesh,PolySwap);
-    }
 };
 
 }//end namespace vcg