From 83f0deca4d28cf3ac59f9ec053b6571ebc7edaa5 Mon Sep 17 00:00:00 2001
From: nicopietroni <nico.pietroni@isti.cnr.it>
Date: Mon, 8 Feb 2016 13:42:29 +0000
Subject: [PATCH] added functionalities to evaluate distortion wrt a cross
 field

---
 .../algorithms/parametrization/distortion.h   | 572 ++++++++++--------
 1 file changed, 309 insertions(+), 263 deletions(-)

diff --git a/vcg/complex/algorithms/parametrization/distortion.h b/vcg/complex/algorithms/parametrization/distortion.h
index 1abfce66..f5ec2393 100644
--- a/vcg/complex/algorithms/parametrization/distortion.h
+++ b/vcg/complex/algorithms/parametrization/distortion.h
@@ -24,272 +24,318 @@
 #ifndef VCG_PARAM_DISTORTION
 #define VCG_PARAM_DISTORTION
 #include <vcg/complex/algorithms/parametrization/uv_utils.h>
+#include <vcg/complex/algorithms/parametrization/tangent_field_operators.h>
 
 namespace vcg {
-    namespace tri{
-        template <class MeshType, bool PerWedgeFlag>
-        class Distortion
-        {
-        public:
-            typedef typename MeshType::FaceType FaceType;
-            typedef typename MeshType::VertexType VertexType;
-            typedef typename MeshType::CoordType CoordType;
-            typedef typename MeshType::ScalarType ScalarType;
-            typedef typename MeshType::FaceType::TexCoordType::ScalarType TexScalarType;
+namespace tri{
+template <class MeshType, bool PerWedgeFlag>
+class Distortion
+{
+public:
+    typedef typename MeshType::FaceType FaceType;
+    typedef typename MeshType::VertexType VertexType;
+    typedef typename MeshType::CoordType CoordType;
+    typedef typename MeshType::ScalarType ScalarType;
+    typedef typename MeshType::FaceType::TexCoordType::ScalarType TexScalarType;
 
-            static ScalarType Area3D(const FaceType *f)
-            {
-                return DoubleArea(*f)*(0.5);
-            }
-
-            static ScalarType AreaUV(const FaceType *f)
-            {
-                Point2<TexScalarType> uv0,uv1,uv2;
-                if(PerWedgeFlag) {
-                  uv0=f->cWT(0).P();
-                  uv1=f->cWT(1).P();
-                  uv2=f->cWT(2).P();
-                } else {
-                  uv0=f->cV(0)->T().P();
-                  uv1=f->cV(1)->T().P();
-                  uv2=f->cV(2)->T().P();
-                }
-                ScalarType AreaUV=((uv1-uv0)^(uv2-uv0))/2.0;
-                return AreaUV;
-            }
-
-            static ScalarType EdgeLenght3D(const FaceType *f,int e)
-            {
-                assert((e>=0)&&(e<3));
-                ScalarType length=(f->cP0(e)-f->cP1(e)).Norm();
-                return (length);
-            }
-
-            static ScalarType EdgeLenghtUV(const FaceType *f,int e)
-            {
-                assert((e>=0)&&(e<3));
-                Point2<TexScalarType> uv0,uv1;
-                if(PerWedgeFlag) {
-                  uv0=f->cWT(e+0).P();
-                  uv1=f->cWT((e+1)%3).P();
-                } else {
-                  uv0=f->cV0(e)->T().P();
-                  uv1=f->cV1(e)->T().P();
-                }
-                ScalarType UVlength=Distance(uv0,uv1);
-                return UVlength;
-            }
-
-            static ScalarType AngleCos3D(const FaceType *f,int e)
-            {
-                assert((e>=0)&&(e<3));
-                CoordType p0=f->P((e+2)%3);
-                CoordType p1=f->P(e);
-                CoordType p2=f->P((e+1)%3);
-                typedef typename CoordType::ScalarType ScalarType;
-                CoordType dir0=p2-p1;
-                CoordType dir1=p0-p1;
-                dir0.Normalize();
-                dir1.Normalize();
-                ScalarType angle=dir0*dir1;
-                return angle;
-            }
-
-            static ScalarType AngleCosUV(const FaceType *f,int e)
-            {
-              Point2<ScalarType> uv0,uv1,uv2;
-              if(PerWedgeFlag) {
-                uv0=f->cWT((e+2)%3).P();
-                uv1=f->cWT((e+0)%3).P();
-                uv2=f->cWT((e+1)%3).P();
-              } else {
-                uv0=f->V2(e)->T().P();
-                uv1=f->V0(e)->T().P();
-                uv2=f->V1(e)->T().P();
-              }
-                vcg::Point2<ScalarType> dir0=uv2-uv1;
-                vcg::Point2<ScalarType> dir1=uv0-uv1;
-                dir0.Normalize();
-                dir1.Normalize();
-                ScalarType angle=dir0*dir1;
-                return angle;
-            }
-
-            static ScalarType AngleRad3D(const FaceType *f,int e)
-            {
-                assert((e>=0)&&(e<3));
-                CoordType p0=f->cP((e+2)%3);
-                CoordType p1=f->cP(e);
-                CoordType p2=f->cP((e+1)%3);
-                typedef typename CoordType::ScalarType ScalarType;
-                CoordType dir0=p2-p1;
-                CoordType dir1=p0-p1;
-                return Angle(dir0,dir1);
-            }
-
-            static ScalarType AngleRadUV(const FaceType *f,int e)
-            {
-              Point2<TexScalarType> uv0,uv1,uv2;
-              if(PerWedgeFlag) {
-                uv0=f->cWT((e+2)%3).P();
-                uv1=f->cWT((e+0)%3).P();
-                uv2=f->cWT((e+1)%3).P();
-              } else {
-                uv0=f->cV2(e)->T().P();
-                uv1=f->cV0(e)->T().P();
-                uv2=f->cV1(e)->T().P();
-              }
-                vcg::Point2<TexScalarType> dir0=uv2-uv1;
-                vcg::Point2<TexScalarType> dir1=uv0-uv1;
-                dir0.Normalize();
-                dir1.Normalize();
-                ScalarType t=dir0*dir1;
-                if(t>1) t = 1;
-                else if(t<-1) t = -1;
-                return acos(t);
-            }
-
-
-        public:
-            enum DistType{AreaDist,EdgeDist,AngleDist};
-
-            ///return the absolute difference between angle in 3D space and texture space
-            ///Actually the difference in cos space
-            static ScalarType AngleCosDistortion(const FaceType *f,int e)
-            {
-                ScalarType Angle_3D=AngleCos3D(f,e);
-                ScalarType Angle_UV=AngleCosUV(f,e);
-                ScalarType diff=fabs(Angle_3D-Angle_UV);///Angle_3D;
-                return diff;
-            }
-            ///return the absolute difference between angle in 3D space and texture space
-            ///Actually the difference in cos space
-            static ScalarType AngleRadDistortion(const FaceType *f,int e)
-            {
-                ScalarType Angle_3D=AngleRad3D(f,e);
-                ScalarType Angle_UV=AngleRadUV(f,e);
-                ScalarType diff=fabs(Angle_3D-Angle_UV);///Angle_3D;
-                return diff;
-            }
-
-            ///return the variance of angle, normalized
-            ///in absolute value
-            static ScalarType AngleDistortion(const FaceType *f)
-            {
-                return  AngleRadDistortion(f,0) +
-                        AngleRadDistortion(f,1) +
-                        AngleRadDistortion(f,2);
-            }
-
-            ///return the global scaling factors  from 3D to UV
-            static void MeshScalingFactor(const MeshType &m,
-                                          ScalarType &AreaScale,
-                                          ScalarType &EdgeScale)
-            {
-                ScalarType SumArea3D=0;
-                ScalarType SumArea2D=0;
-                ScalarType SumEdge3D=0;
-                ScalarType SumEdge2D=0;
-                for (size_t i=0;i<m.face.size();i++)
-                {
-                    SumArea3D+=Area3D(&m.face[i]);
-                    SumArea2D+=AreaUV(&m.face[i]);
-                    for (int j=0;j<3;j++)
-                    {
-                        SumEdge3D+=EdgeLenght3D(&m.face[i],j);
-                        SumEdge2D+=EdgeLenghtUV(&m.face[i],j);
-                    }
-                }
-                AreaScale=SumArea3D/SumArea2D;
-                EdgeScale=SumEdge3D/SumEdge2D;
-            }
-
-            ///return the variance of edge length, normalized in absolute value,
-            // the needed scaling factor EdgeScaleVal may be calculated
-            ///by using the ScalingFactor function
-            static ScalarType EdgeDistortion(const FaceType *f,int e,
-                ScalarType EdgeScaleVal)
-            {
-                ScalarType edgeUV=EdgeLenghtUV(f,e)*EdgeScaleVal;
-                ScalarType edge3D=EdgeLenght3D(f,e);
-                assert(edge3D > 0);
-                ScalarType diff=fabs(edge3D-edgeUV)/edge3D;
-                assert(!math::IsNAN(diff));
-                return diff;
-            }
-
-            ///return the variance of area, normalized
-            ///in absolute value, the scalar AreaScaleVal may be calculated
-            ///by using the ScalingFactor function
-            static ScalarType AreaDistortion(const FaceType *f,
-                                             ScalarType AreaScaleVal)
-            {
-              ScalarType areaUV=AreaUV(f)*AreaScaleVal;
-              ScalarType area3D=Area3D(f);
-              assert(area3D > 0);
-              ScalarType diff=fabs(areaUV-area3D)/area3D;
-              assert(!math::IsNAN(diff));
-              return diff;
-            }
-
-            ///return the number of folded faces
-            static bool Folded(const FaceType *f)
-            {
-                ScalarType areaUV=AreaUV(f);
-                /*if (areaUV<0)
-                    printf("area %5.5f \n",areaUV);*/
-                return (areaUV<0);
-            }
-
-            static int Folded(const MeshType &m)
-            {
-                int folded=0;
-                for (size_t i=0;i<m.face.size();i++)
-                {
-                    if (m.face[i].IsD())continue;
-                    if(Folded(&m.face[i]))folded++;
-                }
-                return folded;
-            }
-
-            static bool GloballyUnFolded(const MeshType &m)
-            {
-                int num=Folded(m);
-                return (num>(m.fn)/2);
-            }
-
-            static ScalarType MeshAngleDistortion(const MeshType &m)
-            {
-                ScalarType UDdist=0;
-                for (int i=0;i<m.face.size();i++)
-                {
-                    if (m.face[i].IsD())continue;
-                    const FaceType *f=&(m.face[i]);
-                    UDdist+=AngleDistortion(f)*Area3D(f);
-                }
-                return UDdist;
-            }
-
-            static void SetQasDistorsion(MeshType &m,
-                                         DistType DType=AreaDist)
-            {
-                ScalarType edge_scale,area_scale;
-                MeshScalingFactor(m,area_scale,edge_scale);
-                for (int i=0;i<m.face.size();i++)
-                {
-                    if (m.face[i].IsD())continue;
-                    if (DType==AreaDist)
-                        m.face[i].Q()=1-AreaDistortion(&m.face[i],area_scale);
-                    else
-                    if (DType==AngleDist)
-                        m.face[i].Q()=1-AngleDistortion(&m.face[i]);
-                    else
-                        m.face[i].Q()=3-EdgeDistortion(&m.face[i],0,edge_scale)-
-                                        EdgeDistortion(&m.face[i],1,edge_scale)-
-                                        EdgeDistortion(&m.face[i],2,edge_scale);
-                }
-            }
-        };
+    static ScalarType Area3D(const FaceType *f)
+    {
+        return DoubleArea(*f)*(0.5);
     }
+
+    static ScalarType AreaUV(const FaceType *f)
+    {
+        Point2<TexScalarType> uv0,uv1,uv2;
+        if(PerWedgeFlag) {
+            uv0=f->cWT(0).P();
+            uv1=f->cWT(1).P();
+            uv2=f->cWT(2).P();
+        } else {
+            uv0=f->cV(0)->T().P();
+            uv1=f->cV(1)->T().P();
+            uv2=f->cV(2)->T().P();
+        }
+        ScalarType AreaUV=((uv1-uv0)^(uv2-uv0))/2.0;
+        return AreaUV;
+    }
+
+    static ScalarType EdgeLenght3D(const FaceType *f,int e)
+    {
+        assert((e>=0)&&(e<3));
+        ScalarType length=(f->cP0(e)-f->cP1(e)).Norm();
+        return (length);
+    }
+
+    static ScalarType EdgeLenghtUV(const FaceType *f,int e)
+    {
+        assert((e>=0)&&(e<3));
+        Point2<TexScalarType> uv0,uv1;
+        if(PerWedgeFlag) {
+            uv0=f->cWT(e+0).P();
+            uv1=f->cWT((e+1)%3).P();
+        } else {
+            uv0=f->cV0(e)->T().P();
+            uv1=f->cV1(e)->T().P();
+        }
+        ScalarType UVlength=Distance(uv0,uv1);
+        return UVlength;
+    }
+
+    static ScalarType AngleCos3D(const FaceType *f,int e)
+    {
+        assert((e>=0)&&(e<3));
+        CoordType p0=f->P((e+2)%3);
+        CoordType p1=f->P(e);
+        CoordType p2=f->P((e+1)%3);
+        typedef typename CoordType::ScalarType ScalarType;
+        CoordType dir0=p2-p1;
+        CoordType dir1=p0-p1;
+        dir0.Normalize();
+        dir1.Normalize();
+        ScalarType angle=dir0*dir1;
+        return angle;
+    }
+
+    static ScalarType AngleCosUV(const FaceType *f,int e)
+    {
+        Point2<ScalarType> uv0,uv1,uv2;
+        if(PerWedgeFlag) {
+            uv0=f->cWT((e+2)%3).P();
+            uv1=f->cWT((e+0)%3).P();
+            uv2=f->cWT((e+1)%3).P();
+        } else {
+            uv0=f->V2(e)->T().P();
+            uv1=f->V0(e)->T().P();
+            uv2=f->V1(e)->T().P();
+        }
+        vcg::Point2<ScalarType> dir0=uv2-uv1;
+        vcg::Point2<ScalarType> dir1=uv0-uv1;
+        dir0.Normalize();
+        dir1.Normalize();
+        ScalarType angle=dir0*dir1;
+        return angle;
+    }
+
+    static ScalarType AngleRad3D(const FaceType *f,int e)
+    {
+        assert((e>=0)&&(e<3));
+        CoordType p0=f->cP((e+2)%3);
+        CoordType p1=f->cP(e);
+        CoordType p2=f->cP((e+1)%3);
+        typedef typename CoordType::ScalarType ScalarType;
+        CoordType dir0=p2-p1;
+        CoordType dir1=p0-p1;
+        return Angle(dir0,dir1);
+    }
+
+    static ScalarType AngleRadUV(const FaceType *f,int e)
+    {
+        Point2<TexScalarType> uv0,uv1,uv2;
+        if(PerWedgeFlag) {
+            uv0=f->cWT((e+2)%3).P();
+            uv1=f->cWT((e+0)%3).P();
+            uv2=f->cWT((e+1)%3).P();
+        } else {
+            uv0=f->cV2(e)->T().P();
+            uv1=f->cV0(e)->T().P();
+            uv2=f->cV1(e)->T().P();
+        }
+        vcg::Point2<TexScalarType> dir0=uv2-uv1;
+        vcg::Point2<TexScalarType> dir1=uv0-uv1;
+        dir0.Normalize();
+        dir1.Normalize();
+        ScalarType t=dir0*dir1;
+        if(t>1) t = 1;
+        else if(t<-1) t = -1;
+        return acos(t);
+    }
+
+
+public:
+    enum DistType{AreaDist,EdgeDist,AngleDist,CrossDist};
+
+    ///return the absolute difference between angle in 3D space and texture space
+    ///Actually the difference in cos space
+    static ScalarType AngleCosDistortion(const FaceType *f,int e)
+    {
+        ScalarType Angle_3D=AngleCos3D(f,e);
+        ScalarType Angle_UV=AngleCosUV(f,e);
+        ScalarType diff=fabs(Angle_3D-Angle_UV);///Angle_3D;
+        return diff;
+    }
+    ///return the absolute difference between angle in 3D space and texture space
+    ///Actually the difference in cos space
+    static ScalarType AngleRadDistortion(const FaceType *f,int e)
+    {
+        ScalarType Angle_3D=AngleRad3D(f,e);
+        ScalarType Angle_UV=AngleRadUV(f,e);
+        ScalarType diff=fabs(Angle_3D-Angle_UV)/Angle_3D;///Angle_3D;
+        return diff;
+    }
+
+    ///return the variance of angle, normalized
+    ///in absolute value
+    static ScalarType AngleDistortion(const FaceType *f)
+    {
+        return  (AngleRadDistortion(f,0) +
+                AngleRadDistortion(f,1) +
+                AngleRadDistortion(f,2))/3.0;
+    }
+
+    ///return the global scaling factors  from 3D to UV
+    static void MeshScalingFactor(const MeshType &m,
+                                  ScalarType &AreaScale,
+                                  ScalarType &EdgeScale)
+    {
+        ScalarType SumArea3D=0;
+        ScalarType SumArea2D=0;
+        ScalarType SumEdge3D=0;
+        ScalarType SumEdge2D=0;
+        for (size_t i=0;i<m.face.size();i++)
+        {
+            SumArea3D+=Area3D(&m.face[i]);
+            SumArea2D+=AreaUV(&m.face[i]);
+            for (int j=0;j<3;j++)
+            {
+                SumEdge3D+=EdgeLenght3D(&m.face[i],j);
+                SumEdge2D+=EdgeLenghtUV(&m.face[i],j);
+            }
+        }
+        AreaScale=SumArea3D/SumArea2D;
+        EdgeScale=SumEdge3D/SumEdge2D;
+    }
+
+    ///return the variance of edge length, normalized in absolute value,
+    // the needed scaling factor EdgeScaleVal may be calculated
+    ///by using the ScalingFactor function
+    static ScalarType EdgeDistortion(const FaceType *f,int e,
+                                     ScalarType EdgeScaleVal)
+    {
+        ScalarType edgeUV=EdgeLenghtUV(f,e)*EdgeScaleVal;
+        ScalarType edge3D=EdgeLenght3D(f,e);
+        assert(edge3D > 0);
+        ScalarType diff=fabs(edge3D-edgeUV)/edge3D;
+        assert(!math::IsNAN(diff));
+        return diff;
+    }
+
+    ///return the variance of area, normalized
+    ///in absolute value, the scalar AreaScaleVal may be calculated
+    ///by using the ScalingFactor function
+    static ScalarType AreaDistortion(const FaceType *f,
+                                     ScalarType AreaScaleVal)
+    {
+        ScalarType areaUV=AreaUV(f)*AreaScaleVal;
+        ScalarType area3D=Area3D(f);
+        assert(area3D > 0);
+        ScalarType diff=fabs(areaUV-area3D)/area3D;
+        assert(!math::IsNAN(diff));
+        return diff;
+    }
+
+    ///return the number of folded faces
+    static bool Folded(const FaceType *f)
+    {
+        ScalarType areaUV=AreaUV(f);
+        /*if (areaUV<0)
+                    printf("area %5.5f \n",areaUV);*/
+        return (areaUV<0);
+    }
+
+    static int Folded(const MeshType &m)
+    {
+        int folded=0;
+        for (size_t i=0;i<m.face.size();i++)
+        {
+            if (m.face[i].IsD())continue;
+            if(Folded(&m.face[i]))folded++;
+        }
+        return folded;
+    }
+
+    static bool GloballyUnFolded(const MeshType &m)
+    {
+        int num=Folded(m);
+        return (num>(m.fn)/2);
+    }
+
+    static ScalarType MeshAngleDistortion(const MeshType &m)
+    {
+        ScalarType UDdist=0;
+        for (int i=0;i<m.face.size();i++)
+        {
+            if (m.face[i].IsD())continue;
+            const FaceType *f=&(m.face[i]);
+            UDdist+=AngleDistortion(f)*Area3D(f);
+        }
+        return UDdist;
+    }
+
+    static void SetQasCrossDirDistortion(MeshType &m)
+    {
+        //first save the old UV dir
+        std::vector<CoordType> Dir1,Dir2;
+        for (size_t i=0;i<m.face.size();i++)
+        {
+            Dir1.push_back(m.face[i].PD1());
+            Dir2.push_back(m.face[i].PD2());
+        }
+        vcg::tri::CrossField<MeshType>::InitDirFromWEdgeUV(m);
+
+        //then compute angle deficit
+        for (size_t i=0;i<m.face.size();i++)
+        {
+            CoordType transfPD1=vcg::tri::CrossField<MeshType>::K_PI(Dir1[i],
+                                                                     m.face[i].PD1(),
+                                                                     m.face[i].N());
+            transfPD1.Normalize();
+            ScalarType AngleDeficit=vcg::Angle(transfPD1,m.face[i].PD1());
+            AngleDeficit=math::ToDeg(AngleDeficit);
+            if ((AngleDeficit>45)||(AngleDeficit<0))
+            {
+                std::cout<<"Warnign A Deficit "<<AngleDeficit<<std::endl;
+            }
+//            assert(AngleDeficit<45);
+//            assert(AngleDeficit>=0);
+            m.face[i].Q()=(AngleDeficit)/(ScalarType)45;
+        }
+
+        //finally restore the original directions
+        for (size_t i=0;i<m.face.size();i++)
+        {
+            m.face[i].PD1()=Dir1[i];
+            m.face[i].PD2()=Dir2[i];
+        }
+    }
+
+    static void SetQasDistorsion(MeshType &m,
+                                 DistType DType=AreaDist)
+    {
+        if (DType==CrossDist)
+        {
+            SetQasCrossDirDistortion(m);
+            vcg::tri::UpdateQuality<MeshType>::VertexFromFace(m,true);
+            return;
+        }
+
+        ScalarType edge_scale,area_scale;
+        MeshScalingFactor(m,area_scale,edge_scale);
+        for (int i=0;i<m.face.size();i++)
+        {
+            if (m.face[i].IsD())continue;
+            if (DType==AreaDist)
+                m.face[i].Q()=AreaDistortion(&m.face[i],area_scale);
+            else
+                if (DType==AngleDist)
+                    m.face[i].Q()=AngleDistortion(&m.face[i]);
+                else
+                    m.face[i].Q()=(EdgeDistortion(&m.face[i],0,edge_scale)+
+                            EdgeDistortion(&m.face[i],1,edge_scale)+
+                            EdgeDistortion(&m.face[i],2,edge_scale))/3.0;
+        }
+        vcg::tri::UpdateQuality<MeshType>::VertexFromFace(m,true);
+    }
+};
+}
 }
 #endif