From 94d9a3dbdd25a8338f0bfc601be0a9cbdc28e272 Mon Sep 17 00:00:00 2001
From: nicopietroni <nico.pietroni@isti.cnr.it>
Date: Mon, 5 Jan 2015 12:57:16 +0000
Subject: [PATCH] * fixed several bugs * added several functions to compute the
 quality of a polygonal meshing * added the template polygon computation as in
 "Static Aware Grid Shells" by Pietroni et Al.

---
 vcg/space/polygon3.h | 512 +++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 497 insertions(+), 15 deletions(-)

diff --git a/vcg/space/polygon3.h b/vcg/space/polygon3.h
index d79d1ed7..dd267214 100644
--- a/vcg/space/polygon3.h
+++ b/vcg/space/polygon3.h
@@ -24,33 +24,515 @@
 #ifndef POLYGON_H
 #define POLYGON_H
 
+#include <vcg/space/plane3.h>
+#include <vcg/space/fitting3.h>
+#include <vcg/space/point_matching.h>
+#include <vcg/math/matrix33.h>
+
 namespace vcg {
 
-/*
- */
+////return true if the
+//template <class CoordType>
+//bool CheckNormalizedCoords(CoordType dir)
+//{
+//    typedef typename CoordType::ScalarType ScalarType;
+//    if(isnan(dir.X()))return false;
+//    if(isnan(dir.Y()))return false;
+//    if(isnan(dir.Z()))return false;
+//    ScalarType Norm=dir.Norm();
+//    if(fabs(Norm-1.f)>0.01f)return false;
+//    return true;
+//}
 
-template<class PolygonType>
-typename PolygonType::CoordType PolygonBarycenter(PolygonType &F)
+//return per vertex Normals of a polygonal face stored as a vector of coords
+template <class CoordType>
+void GetNormals(std::vector<CoordType> &Pos,
+                std::vector<CoordType> &Norms)
 {
-  typename PolygonType::CoordType bary(0,0,0);
-  for (int i=0;i<F.VN();i++)
-    bary+=F.V(i)->P();
-
-  bary/=(typename PolygonType::ScalarType)F.VN();
-  return bary;
+    Norms.clear();
+    int size=Pos.size();
+    if (size<=2) return;
+    for (int i=0;i<size;i++)
+        Norms.push_back(Normal(Pos[i],Pos[(i+1)%size],Pos[(i+2)%size]).Normalize());
 }
 
+//return the normal of a polygonal face stored as a vector of coords
+template <class CoordType>
+CoordType Normal(std::vector<CoordType> &Pos)
+{
+    std::vector<CoordType> Norms;
+    GetNormals(Pos,Norms);
+    if (Norms.size()==0)
+        return(CoordType(1,0,0));
+
+    CoordType NSum=CoordType(0,0,0);
+    for (size_t i=0;i<Norms.size();i++)
+        NSum+=Norms[i];
+
+    NSum.Normalize();
+    return (NSum);
+}
+
+//return the area of a polygonal face stored as a vector of coords
+template <class CoordType>
+typename CoordType::ScalarType Area(const std::vector<CoordType> &Pos)
+{
+    typedef typename CoordType::ScalarType ScalarType;
+    CoordType bary=CoordType(0,0,0);
+    for (int i=0;i<Pos.size();i++)
+        bary+=Pos[i];
+
+    bary/=Pos.size();
+    ScalarType Area=0;
+    for (size_t i=0;i<Pos.size();i++)
+    {
+        CoordType p0=Pos[i];
+        CoordType p1=Pos[(i+1)% Pos.size()];
+        CoordType p2=bary;
+        vcg::Triangle3<ScalarType> T(p0,p1,p2);
+        Area+=(vcg::DoubleArea(T)/2);
+    }
+    return Area;
+}
+
+//return per vertex Normals of a polygonal face
+template<class PolygonType>
+void PolyNormals(PolygonType &F,
+                 std::vector<typename PolygonType::CoordType> &Norms)
+{
+    typedef typename PolygonType::FaceType FaceType;
+    typedef typename PolygonType::CoordType CoordType;
+    typedef typename PolygonType::ScalarType ScalarType;
+
+    Norms.clear();
+    if (F.VN()<=2) return;
+    for (int i=0;i<F.VN();i++)
+        Norms.push_back(Normal(F.P0(i),F.P1(i),F.P2(i)).Normalize());
+}
+
+//return the barycenter of a polygonal face
+template<class PolygonType>
+typename PolygonType::CoordType PolyBarycenter(PolygonType &F)
+{
+    typename PolygonType::CoordType bary(0,0,0);
+    for (int i=0;i<F.VN();i++)
+        bary+=F.V(i)->P();
+
+    bary/=(typename PolygonType::ScalarType)F.VN();
+    return bary;
+}
+
+//return the area of a polygonal face
+template<class PolygonType>
+typename PolygonType::ScalarType PolyArea(PolygonType &F)
+{
+    typedef typename PolygonType::FaceType FaceType;
+    typedef typename PolygonType::CoordType CoordType;
+    typedef typename PolygonType::ScalarType ScalarType;
+
+    CoordType bary=PolyBarycenter(F);
+    ScalarType Area=0;
+    for (size_t i=0;i<F.VN();i++)
+    {
+        CoordType p0=F.P0(i);
+        CoordType p1=F.P1(i);
+        CoordType p2=bary;
+        vcg::Triangle3<ScalarType> T(p0,p1,p2);
+        Area+=(vcg::DoubleArea(T)/2);
+    }
+    return Area;
+}
+
+//return the normal of a polygonal face
 template<class PolygonType>
 typename PolygonType::CoordType PolygonNormal(PolygonType &F)
 {
-  typename PolygonType::CoordType n(0,0,0);
+    typename PolygonType::CoordType n(0,0,0);
 
-  for (int i=0;i<F.VN();i++)
-    n+=Normal(F.P0(i),F.P(1),F.P2()).Normalize();
+    for (int i=0;i<F.VN();i++)
+        n+=Normal(F.P0(i),F.P1(i),F.P2(i)).Normalize();
 
-  return n.Normalize();
+    return n.Normalize();
 }
 
+//return the perimeter of a polygonal face
+template<class PolygonType>
+typename PolygonType::ScalarType PolyPerimeter(PolygonType &F)
+{
+    typedef typename PolygonType::FaceType FaceType;
+    typedef typename PolygonType::CoordType CoordType;
+    typedef typename PolygonType::ScalarType ScalarType;
 
- }
+    ScalarType SumL=0;
+    for (int i=0;i<F.VN();i++)
+    {
+        ScalarType L=(F.P0(i)-F.P1(i)).Norm();
+        SumL+=L;
+    }
+    return (SumL);
+}
+
+//return a Scalar value that encode the variance of the normals
+//wrt the average one (1 means hight variance, 0 no variance)
+template<class PolygonType>
+typename PolygonType::ScalarType PolyNormDeviation(PolygonType &F)
+{
+    typedef typename PolygonType::FaceType FaceType;
+    typedef typename PolygonType::CoordType CoordType;
+    typedef typename PolygonType::ScalarType ScalarType;
+
+    std::vector<typename PolygonType::CoordType> Norms;
+    PolyNormals(F,Norms);
+
+    //calculate the Avg Normal
+    CoordType AvgNorm(0,0,0);
+    for (int i=0;i<Norms.size();i++)
+        AvgNorm+=Norms[i];
+
+    AvgNorm.Normalize();
+
+    if (!CheckNormalizedCoords(AvgNorm))return 1;
+
+    ScalarType Dev=0;
+    for (int i=0;i<Norms.size();i++)
+        Dev+=pow((Norms[i]-AvgNorm).Norm()/2.0,2);
+
+    Dev/=(ScalarType)Norms.size();
+    Dev=sqrt(Dev);
+    return Dev;
+}
+
+//return a Scalar value that encode the distance wrt ideal angle for each
+//wrt the average one (1 correspond to hight variance, 0 no variance)
+template<class PolygonType>
+void PolyAngleDeviation(PolygonType &F,
+                        typename PolygonType::ScalarType &AvgDev,
+                        typename PolygonType::ScalarType &MaxDev)
+{
+    typedef typename PolygonType::FaceType FaceType;
+    typedef typename PolygonType::CoordType CoordType;
+    typedef typename PolygonType::ScalarType ScalarType;
+    assert(F.VN()>2);
+    ScalarType IdealAngle=M_PI-(2*M_PI/(ScalarType)F.VN());
+    assert(IdealAngle>0);
+
+    //then compute the angle deviation
+    MaxDev=0;
+    AvgDev=0;
+
+    for (int i=0;i<F.VN();i++)
+    {
+        CoordType dir0=F.P0(i)-F.P1(i);
+        CoordType dir1=F.P2(i)-F.P1(i);
+
+        ScalarType VAngle=vcg::Angle(dir0,dir1);
+        assert(VAngle>=0);
+        ScalarType VAngleDiff=fabs(VAngle-IdealAngle);
+
+        if (VAngleDiff>MaxDev)MaxDev=VAngleDiff;
+
+        AvgDev+=VAngleDiff;
+    }
+    AvgDev/=(ScalarType)F.VN();
+
+    AvgDev/=(M_PI/2.0);
+    MaxDev/=(M_PI/2.0);
+
+    if (AvgDev>1)AvgDev=1;
+    if (MaxDev>1)MaxDev=1;
+}
+
+//return the fitting plane of a polygonal face
+template<class PolygonType>
+vcg::Plane3<typename PolygonType::ScalarType> PolyFittingPlane(PolygonType &F)
+{
+    typedef typename PolygonType::FaceType FaceType;
+    typedef typename PolygonType::CoordType CoordType;
+    typedef typename PolygonType::ScalarType ScalarType;
+    vcg::Plane3<ScalarType> BestPL;
+    assert(F.VN()>=3);
+    std::vector<CoordType> pointVec;
+    for (int i=0;i<F.VN();i++)
+        pointVec.push_back(F.P(i));
+
+    vcg::FitPlaneToPointSet(pointVec,BestPL);
+    return BestPL;
+}
+
+//return the flatness of a polygonal face as avg distance to the best fitting plane divided by half perimeter
+template<class PolygonType>
+typename PolygonType::ScalarType PolyFlatness(PolygonType &F)
+{
+    typedef typename PolygonType::FaceType FaceType;
+    typedef typename PolygonType::CoordType CoordType;
+    typedef typename PolygonType::ScalarType ScalarType;
+
+    if (F.VN()<=3)
+        return 0;
+
+    //average lenght
+    ScalarType SumL=PolyPerimeter(F)/2.0;
+
+    //diagonal distance
+    vcg::Plane3<ScalarType> BestPL=PolyFittingPlane(F);
+
+    //then project points on the plane
+    ScalarType Flatness=0;
+    for (int i=0;i<F.VN();i++)
+    {
+        CoordType pos=F.P(i);
+        CoordType proj=BestPL.Projection(pos);
+        Flatness+=(pos-proj).Norm();
+    }
+    Flatness/=(ScalarType)F.VN();
+    return((Flatness)/SumL);
+}
+
+//evaluate the PCA directions of a polygonal face
+template<class PolygonType>
+void PolyPCA(PolygonType &F,
+             typename PolygonType::CoordType PCA[])
+{
+    typedef typename PolygonType::FaceType FaceType;
+    typedef typename PolygonType::CoordType CoordType;
+    typedef typename PolygonType::ScalarType ScalarType;
+
+    //compute the covariance matrix
+    Eigen::Matrix3d EigenCovMat;
+    //ComputeCovarianceMatrix(EigenCovMat);
+    //compute covariance matrix
+    ///compute the barycenter
+    CoordType Barycenter=PolyBarycenter(F);
+
+    // second cycle: compute the covariance matrix
+    EigenCovMat.setZero();
+    Eigen::Vector3d p;
+    for (int i=0;i<F.VN();i++)
+    {
+        (F.V(i)->P()-Barycenter).ToEigenVector(p);
+        EigenCovMat+= p*p.transpose(); // outer product
+    }
+
+    Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d > eig(EigenCovMat);
+
+    Eigen::Vector3d eval = eig.eigenvalues();
+    Eigen::Matrix3d evec = eig.eigenvectors();
+
+    eval = eval.cwiseAbs();
+    int normInd,maxInd,minInd;
+
+    ///get min and max coff ..
+    ///the minumum is the Normal
+    ///the other two the anisotropy directions
+    eval.minCoeff(&normInd);
+    eval.maxCoeff(&maxInd);
+    minInd=(maxInd+1)%3;
+
+    if (minInd==normInd)minInd=(normInd+1)%3;
+    assert((minInd!=normInd)&&(minInd!=maxInd)&&(minInd!=maxInd));
+
+    ///maximum direction of PCA
+    PCA[0][0] = evec(0,maxInd);
+    PCA[0][1] = evec(1,maxInd);
+    PCA[0][2] = evec(2,maxInd);
+    ///minimum direction of PCA
+    PCA[1][0] = evec(0,minInd);
+    PCA[1][1] = evec(1,minInd);
+    PCA[1][2] = evec(2,minInd);
+    ///Normal direction
+    PCA[2][0] = evec(0,normInd);
+    PCA[2][1] = evec(1,normInd);
+    PCA[2][2] = evec(2,normInd);
+
+    ScalarType LX=sqrt(eval[maxInd]);
+    ScalarType LY=sqrt(eval[minInd]);
+    //ScalarType LZ=sqrt(eval[normInd]);
+
+    ///scale the directions
+    PCA[0]*=LX;
+    PCA[1]*=LY;
+    //PCA[2]*=LZ;//.Normalize();
+    PCA[2].Normalize();
+}
+
+//evaluate the PCA directions of a polygonal face
+//scaled by the area of the face
+template<class PolygonType>
+void PolyScaledPCA(PolygonType &F,
+                   typename PolygonType::CoordType PCA[])
+{
+    typedef typename PolygonType::FaceType FaceType;
+    typedef typename PolygonType::CoordType CoordType;
+    typedef typename PolygonType::ScalarType ScalarType;
+
+    std::vector<CoordType> SwapPos;
+
+    ///compute the barycenter
+    //CoordType Barycenter=PolyBarycenter(F);
+
+    ///compute the Area
+    ScalarType Area=PolyArea(F);
+
+    PolyPCA(F,PCA);
+
+    ScalarType Scale=sqrt(Area/(PCA[0].Norm()*PCA[1].Norm()));
+    PCA[0]*=Scale;
+    PCA[1]*=Scale;
+
+}
+
+//return the base template polygon as
+//described by "Static Aware Grid Shells" by Pietroni et Al.
+template<class CoordType>
+void getBaseTemplatePolygon(int N,
+                            std::vector<CoordType> &TemplatePos)
+{
+    typedef typename CoordType::ScalarType ScalarType;
+    ///first find positions in the
+    ///reference frame of the passed matrix
+    ScalarType AngleInterval=2.0*M_PI/(ScalarType)N;
+    ScalarType CurrAngle=0;
+    TemplatePos.resize(N);
+    for (size_t i=0;i<TemplatePos.size();i++)
+    {
+        ///find with trigonometric functions
+        TemplatePos[i].X()=cos(CurrAngle);
+        TemplatePos[i].Y()=sin(CurrAngle);
+        TemplatePos[i].Z()=0;
+        //            TemplatePos[i].Normalize();
+        //            TemplatePos[i].X()*=Anisotropy;
+        //            TemplatePos[i].Y()*=(1-Anisotropy);
+
+        ///increment the angle
+        CurrAngle+=AngleInterval;
+    }
+}
+
+//return the rigidly aligned template polygon as
+//described by "Static Aware Grid Shells" by Pietroni et Al.
+template<class PolygonType>
+void GetPolyTemplatePos(PolygonType &F,
+                        std::vector<typename PolygonType::CoordType> &TemplatePos)
+{
+    typedef typename PolygonType::FaceType FaceType;
+    typedef typename PolygonType::CoordType CoordType;
+    typedef typename PolygonType::ScalarType ScalarType;
+    std::vector<CoordType>  UniformPos,UniformTempl;
+
+    CoordType Barycenter=PolyBarycenter(F);
+
+    getBaseTemplatePolygon(F.VN(),TemplatePos);
+
+    CoordType PCA[3];
+    PolyPCA(F,PCA);
+
+    vcg::Matrix44<ScalarType> ToPCA,ToPCAInv;
+    ToPCA.SetIdentity();
+
+    CoordType dirX=PCA[0];
+    CoordType dirY=PCA[1];
+    CoordType dirZ=PCA[2];
+
+    ///set the Rotation matrix
+    ToPCA.SetColumn(0,dirX);
+    ToPCA.SetColumn(1,dirY);
+    ToPCA.SetColumn(2,dirZ);
+    ToPCAInv=ToPCA;
+    ToPCA=vcg::Inverse(ToPCA);
+
+    ///then transform the polygon to PCA space
+    for (int i=0;i<F.VN();i++)
+    {
+        ///translate
+        CoordType Pos=F.V(i)->P()-Barycenter;
+        ///rotate
+        Pos=ToPCA*Pos;
+        //retranslate
+        UniformPos.push_back(Pos);
+    }
+
+    ///calculate the Area
+    ScalarType AreaTemplate=Area(TemplatePos);
+    ScalarType AreaUniform=Area(UniformPos);
+
+    ScalarType Scale=sqrt(AreaTemplate/AreaUniform);
+
+    for (size_t i=0;i<UniformPos.size();i++)
+        UniformPos[i]*=Scale;
+
+    ///check side
+    CoordType N0=Normal(UniformPos);
+    CoordType N1=Normal(TemplatePos);
+    if ((N0*N1)<0)std::reverse(TemplatePos.begin(),TemplatePos.end());
+
+    ///initialize
+    std::vector<CoordType> FixPoints(UniformPos.begin(),UniformPos.end());
+    std::vector<CoordType> MovPoints(TemplatePos.begin(),TemplatePos.end());
+
+    ///add displacement along Z
+    for (size_t i=0;i<FixPoints.size();i++)
+    {
+        FixPoints[i]+=CoordType(0,0,0.01);
+        MovPoints[i]+=CoordType(0,0,0.01);
+    }
+    ///add original points
+    FixPoints.insert(FixPoints.end(),UniformPos.begin(),UniformPos.end());
+    MovPoints.insert(MovPoints.end(),TemplatePos.begin(),TemplatePos.end());
+
+    ///then find the alignment
+    vcg::Matrix44<ScalarType> Rigid;
+    ///compute rigid match
+    vcg::ComputeRigidMatchMatrix<ScalarType>(FixPoints,MovPoints,Rigid);
+
+    ///then apply transformation
+    UniformTempl.resize(TemplatePos.size(),CoordType(0,0,0));
+    for (size_t i=0;i<TemplatePos.size();i++)
+        UniformTempl[i]=Rigid*TemplatePos[i];
+
+    ///then map back to 3D space
+    for (size_t i=0;i<TemplatePos.size();i++)
+    {
+        TemplatePos[i]=UniformTempl[i];
+        TemplatePos[i]*=1/Scale;
+        TemplatePos[i]=ToPCAInv*TemplatePos[i];
+    }
+
+    //        if (use_fixed_area)
+    //        {
+    //            ScalarType A0=Area(TemplatePos);
+    //            ScalarType A1=FixedArea;
+    //            //ScalarType Scale1=A1/A0;
+    //            for (size_t i=0;i<TemplatePos.size();i++)
+    //                TemplatePos[i]*=A1/A0;
+    //        }
+
+
+    for (size_t i=0;i<TemplatePos.size();i++)
+        TemplatePos[i]+=Barycenter;
+
+}
+
+//compute the aspect ratio using the rigidly aligned template polygon as
+//described by "Static Aware Grid Shells" by Pietroni et Al.
+template<class PolygonType>
+typename PolygonType::ScalarType PolyAspectRatio(PolygonType &F)
+{
+    typedef typename PolygonType::FaceType FaceType;
+    typedef typename PolygonType::CoordType CoordType;
+    typedef typename PolygonType::ScalarType ScalarType;
+    std::vector<CoordType> TemplatePos;
+
+    GetPolyTemplatePos(F,TemplatePos);
+
+    ScalarType diff=0;
+    assert((int)TemplatePos.size()==F.VN());
+
+    ScalarType AreaP=PolyArea(F);
+    for (size_t i=0;i<TemplatePos.size();i++)
+        diff+=pow((TemplatePos[i]-F.P(i)).Norm(),2)/AreaP;
+
+    return(diff);
+}
+
+}
 #endif // POLYGON_H