fixes from Nico Pietroni on the MIQ pipeline and examples

compilation fixes for Apple LLVM version 6.0


Former-commit-id: c5fcec4bb5c3f8a55266d96998d3c9548b38015f

include/igl/comb_line_field.cpp

@@ -0,0 +1,132 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 Nico Pietroni <nico.pietroni@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "comb_line_field.h"
+
+#include <vector>
+#include <deque>
+#include "per_face_normals.h"
+#include "is_border_vertex.h"
+#include "rotation_matrix_from_directions.h"
+
+#include "triangle_triangle_adjacency.h"
+
+namespace igl {
+template <typename DerivedV, typename DerivedF>
+class CombLine
+{
+public:
+
+    const Eigen::PlainObjectBase<DerivedV> &V;
+    const Eigen::PlainObjectBase<DerivedF> &F;
+    const Eigen::PlainObjectBase<DerivedV> &PD1;
+    Eigen::PlainObjectBase<DerivedV> N;
+
+private:
+    // internal
+    Eigen::PlainObjectBase<DerivedF> TT;
+    Eigen::PlainObjectBase<DerivedF> TTi;
+
+
+private:
+
+
+    static inline double Sign(double a){return (double)((a>0)?+1:-1);}
+
+
+private:
+
+    // returns the 180 deg rotation of a (around n) most similar to target b
+    // a and b should be in the same plane orthogonal to N
+    static inline Eigen::Matrix<typename DerivedV::Scalar, 3, 1> K_PI_line(const Eigen::Matrix<typename DerivedV::Scalar, 3, 1>& a,
+                                                                           const Eigen::Matrix<typename DerivedV::Scalar, 3, 1>& b)
+    {
+        typename DerivedV::Scalar scorea = a.dot(b);
+        if (scorea<0)
+            return -a;
+        else
+            return a;
+    }
+
+
+
+public:
+
+    inline CombLine(const Eigen::PlainObjectBase<DerivedV> &_V,
+                    const Eigen::PlainObjectBase<DerivedF> &_F,
+                    const Eigen::PlainObjectBase<DerivedV> &_PD1):
+        V(_V),
+        F(_F),
+        PD1(_PD1)
+    {
+        igl::per_face_normals(V,F,N);
+        igl::triangle_triangle_adjacency(V,F,TT,TTi);
+    }
+
+    inline void comb(Eigen::PlainObjectBase<DerivedV> &PD1out)
+    {
+        PD1out.setZero(F.rows(),3);PD1out<<PD1;
+
+        Eigen::VectorXi mark = Eigen::VectorXi::Constant(F.rows(),false);
+
+        std::deque<int> d;
+
+        d.push_back(0);
+        mark(0) = true;
+
+        while (!d.empty())
+        {
+            int f0 = d.at(0);
+            d.pop_front();
+            for (int k=0; k<3; k++)
+            {
+                int f1 = TT(f0,k);
+                if (f1==-1) continue;
+                if (mark(f1)) continue;
+
+                Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir0  = PD1out.row(f0);
+                Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir1  = PD1out.row(f1);
+                Eigen::Matrix<typename DerivedV::Scalar, 3, 1> n0    = N.row(f0);
+                Eigen::Matrix<typename DerivedV::Scalar, 3, 1> n1    = N.row(f1);
+
+                Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir0Rot = igl::rotation_matrix_from_directions(n0, n1)*dir0;
+                dir0Rot.normalize();
+                Eigen::Matrix<typename DerivedV::Scalar, 3, 1> targD   = K_PI_line(dir1,dir0Rot);
+
+                PD1out.row(f1)  = targD;
+                //PD2out.row(f1)  = n1.cross(targD).normalized();
+
+                mark(f1) = true;
+                d.push_back(f1);
+
+            }
+        }
+
+        // everything should be marked
+        for (int i=0; i<F.rows(); i++)
+        {
+            assert(mark(i));
+        }
+    }
+
+};
+}
+
+template <typename DerivedV, typename DerivedF>
+IGL_INLINE void igl::comb_line_field(const Eigen::PlainObjectBase<DerivedV> &V,
+                                     const Eigen::PlainObjectBase<DerivedF> &F,
+                                     const Eigen::PlainObjectBase<DerivedV> &PD1,
+                                     Eigen::PlainObjectBase<DerivedV> &PD1out)
+{
+    igl::CombLine<DerivedV, DerivedF> cmb(V, F, PD1);
+    cmb.comb(PD1out);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+#endif

include/igl/comb_line_field.h

@@ -0,0 +1,39 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 Nico Pietroni <nico.pietroni@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef IGL_COMB_LINE_FIELD_H
+#define IGL_COMB_LINE_FIELD_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Computes principal matchings of the vectors of a cross field across face edges,
+  // and generates a combed cross field defined on the mesh faces
+
+  // Inputs:
+  //   V          #V by 3 eigen Matrix of mesh vertex 3D positions
+  //   F          #F by 4 eigen Matrix of face (quad) indices
+  //   PD1in      #F by 3 eigen Matrix of the first per face cross field vector
+  //   PD2in      #F by 3 eigen Matrix of the second per face cross field vector
+  // Output:
+  //   PD1out      #F by 3 eigen Matrix of the first combed cross field vector
+  //   PD2out      #F by 3 eigen Matrix of the second combed cross field vector
+  //
+
+
+  template <typename DerivedV, typename DerivedF>
+  IGL_INLINE void comb_line_field(const Eigen::PlainObjectBase<DerivedV> &V,
+                                  const Eigen::PlainObjectBase<DerivedF> &F,
+                                  const Eigen::PlainObjectBase<DerivedV> &PD1in,
+                                  Eigen::PlainObjectBase<DerivedV> &PD1out);
+}
+#ifndef IGL_STATIC_LIBRARY
+#include "comb_line_field.cpp"
+#endif
+
+#endif

include/igl/comiso/miq.cpp

@@ -974,7 +974,7 @@ IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::SolvePoisson(Eigen::Vect
 
   clearUserConstraint();
   // copy the user constraints number
-  for (int i = 0; i < hardFeatures.size(); ++i)
+  for (size_t i = 0; i < hardFeatures.size(); ++i)
   {
     addSharpEdgeConstraint(hardFeatures[i][0],hardFeatures[i][1]);
   }
@@ -1445,7 +1445,7 @@ IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::AddSingularityRound()
 template <typename DerivedV, typename DerivedF>
 IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::AddToRoundVertices(std::vector<int> ids)
 {
-  for (int i = 0; i < ids.size(); ++i)
+  for (size_t i = 0; i < ids.size(); ++i)
   {
     if (ids[i] < 0 || ids[i] >= V.rows())
       std::cerr << "WARNING: Ignored round vertex constraint, vertex " << ids[i] << " does not exist in the mesh." << std::endl;
@@ -1789,7 +1789,7 @@ IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::addSharpEdgeConstraint(i
   // prepare constraint
   std::vector<int> c(Handle_SystemInfo.num_vert_variables*2 + 1);
 
-  for (int i = 0; i < c.size(); ++i)
+  for (size_t i = 0; i < c.size(); ++i)
   {
     c[i] = 0;
   }
@@ -2271,7 +2271,9 @@ IGL_INLINE void igl::miq(const Eigen::PlainObjectBase<DerivedV> &V,
            DirectRound,
            iter,
            localIter,
-           DoRound);
+           DoRound,
+           roundVertices,
+           hardFeatures);
 
 }
 

include/igl/doublearea.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "doublearea.h"
 #include "edge_lengths.h"
@@ -12,9 +12,9 @@
 #include <iostream>
 
 template <typename DerivedV, typename DerivedF, typename DeriveddblA>
-IGL_INLINE void igl::doublearea( 
-  const Eigen::PlainObjectBase<DerivedV> & V, 
-  const Eigen::PlainObjectBase<DerivedF> & F, 
+IGL_INLINE void igl::doublearea(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
   Eigen::PlainObjectBase<DeriveddblA> & dblA)
 {
   const int dim = V.cols();
@@ -27,7 +27,7 @@ IGL_INLINE void igl::doublearea(
   // http://www.cs.berkeley.edu/~jrs/meshpapers/robnotes.pdf
 
   // Projected area helper
-  const auto & proj_doublearea = 
+  const auto & proj_doublearea =
     [&V,&F](const int x, const int y, const int f)->double
   {
     auto rx = V(F(f,0),x)-V(F(f,2),x);
@@ -47,7 +47,7 @@ IGL_INLINE void igl::doublearea(
         for(int d = 0;d<3;d++)
         {
           double dblAd = proj_doublearea(d,(d+1)%3,f);
-          dblA[f] += dblAd*dblAd;
+          dblA(f) += dblAd*dblAd;
         }
       }
       dblA = dblA.array().sqrt().eval();
@@ -75,9 +75,9 @@ template <
   typename DerivedB,
   typename DerivedC,
   typename DerivedD>
-IGL_INLINE void doublearea( 
-  const Eigen::PlainObjectBase<DerivedA> & A, 
-  const Eigen::PlainObjectBase<DerivedB> & B, 
+IGL_INLINE void doublearea(
+  const Eigen::PlainObjectBase<DerivedA> & A,
+  const Eigen::PlainObjectBase<DerivedB> & B,
   const Eigen::PlainObjectBase<DerivedC> & C,
   Eigen::PlainObjectBase<DerivedD> & D)
 {
@@ -96,8 +96,8 @@ template <
   typename DerivedB,
   typename DerivedC>
 IGL_INLINE typename DerivedA::Scalar igl::doublearea_single(
-  const Eigen::PlainObjectBase<DerivedA> & A, 
-  const Eigen::PlainObjectBase<DerivedB> & B, 
+  const Eigen::PlainObjectBase<DerivedA> & A,
+  const Eigen::PlainObjectBase<DerivedB> & B,
   const Eigen::PlainObjectBase<DerivedC> & C)
 {
   auto r = A-C;
@@ -106,7 +106,7 @@ IGL_INLINE typename DerivedA::Scalar igl::doublearea_single(
 }
 
 template <typename Derivedl, typename DeriveddblA>
-IGL_INLINE void igl::doublearea( 
+IGL_INLINE void igl::doublearea(
   const Eigen::PlainObjectBase<Derivedl> & ul,
   Eigen::PlainObjectBase<DeriveddblA> & dblA)
 {
@@ -127,12 +127,12 @@ IGL_INLINE void igl::doublearea(
   for(int i = 0;i<m;i++)
   {
     //// Heron's formula for area
-    //const typename Derivedl::Scalar arg = 
+    //const typename Derivedl::Scalar arg =
     //  s(i)*(s(i)-l(i,0))*(s(i)-l(i,1))*(s(i)-l(i,2));
     //assert(arg>=0);
     //dblA(i) = 2.0*sqrt(arg);
     // Kahan's Heron's formula
-    const typename Derivedl::Scalar arg = 
+    const typename Derivedl::Scalar arg =
       (l(i,0)+(l(i,1)+l(i,2)))*
       (l(i,2)-(l(i,0)-l(i,1)))*
       (l(i,2)+(l(i,0)-l(i,1)))*

include/igl/find_cross_field_singularities.cpp

@@ -36,7 +36,7 @@ IGL_INLINE void igl::find_cross_field_singularities(const Eigen::PlainObjectBase
   {
     ///check that is on border..
     if (V_border[vid])
-      return;
+      continue;
 
     int missmatch=0;
     for (unsigned int i=0;i<VF[vid].size();i++)

include/igl/line_field_missmatch.cpp

@@ -0,0 +1,143 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 Nico Pietroni <nico.pietroni@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "line_field_missmatch.h"
+
+#include <vector>
+#include <deque>
+#include <igl/comb_line_field.h>
+#include <igl/rotate_vectors.h>
+#include <igl/comb_cross_field.h>
+#include <igl/comb_line_field.h>
+#include <igl/per_face_normals.h>
+#include <igl/is_border_vertex.h>
+#include <igl/vertex_triangle_adjacency.h>
+#include <igl/triangle_triangle_adjacency.h>
+#include <igl/rotation_matrix_from_directions.h>
+#include <igl/local_basis.h>
+
+namespace igl {
+template <typename DerivedV, typename DerivedF, typename DerivedO>
+class MissMatchCalculatorLine
+{
+public:
+
+    const Eigen::PlainObjectBase<DerivedV> &V;
+    const Eigen::PlainObjectBase<DerivedF> &F;
+    const Eigen::PlainObjectBase<DerivedV> &PD1;
+    const Eigen::PlainObjectBase<DerivedV> &PD2;
+    Eigen::PlainObjectBase<DerivedV> N;
+
+private:
+    // internal
+    std::vector<bool> V_border; // bool
+    std::vector<std::vector<int> > VF;
+    std::vector<std::vector<int> > VFi;
+    Eigen::PlainObjectBase<DerivedF> TT;
+    Eigen::PlainObjectBase<DerivedF> TTi;
+
+
+private:
+
+    //compute the mismatch between 2 faces
+    inline int MissMatchByLine(const int f0,
+                               const int f1)
+    {
+        Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir0 = PD1.row(f0);
+        Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir1 = PD1.row(f1);
+        Eigen::Matrix<typename DerivedV::Scalar, 3, 1> n0 = N.row(f0);
+        Eigen::Matrix<typename DerivedV::Scalar, 3, 1> n1 = N.row(f1);
+
+        Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir1Rot = igl::rotation_matrix_from_directions(n1,n0)*dir1;
+        dir1Rot.normalize();
+
+        // TODO: this should be equivalent to the other code below, to check!
+        // Compute the angle between the two vectors
+        //    double a0 = atan2(dir0.dot(B2.row(f0)),dir0.dot(B1.row(f0)));
+        //    double a1 = atan2(dir1Rot.dot(B2.row(f0)),dir1Rot.dot(B1.row(f0)));
+        //
+        //    double angle_diff = a1-a0;   //VectToAngle(f0,dir1Rot);
+
+        double angle_diff = atan2(dir1Rot.dot(PD2.row(f0)),dir1Rot.dot(PD1.row(f0)));
+
+        double step=M_PI;
+        int i=(int)floor((angle_diff/step)+0.5);
+        assert((i>=-2)&&(i<=2));
+        int k=0;
+        if (i>=0)
+            k=i%2;
+        else
+            k=(2+i)%2;
+
+        assert((k==0)||(k==1));
+        return (k*2);
+    }
+
+public:
+
+    inline MissMatchCalculatorLine(const Eigen::PlainObjectBase<DerivedV> &_V,
+                               const Eigen::PlainObjectBase<DerivedF> &_F,
+                               const Eigen::PlainObjectBase<DerivedV> &_PD1,
+                               const Eigen::PlainObjectBase<DerivedV> &_PD2
+                               ):
+        V(_V),
+        F(_F),
+        PD1(_PD1),
+        PD2(_PD2)
+    {
+        igl::per_face_normals(V,F,N);
+        V_border = igl::is_border_vertex(V,F);
+        igl::vertex_triangle_adjacency(V,F,VF,VFi);
+        igl::triangle_triangle_adjacency(V,F,TT,TTi);
+    }
+
+    inline void calculateMissmatchLine(Eigen::PlainObjectBase<DerivedO> &Handle_MMatch)
+    {
+        Handle_MMatch.setConstant(F.rows(),3,-1);
+        for (unsigned int i=0;i<F.rows();i++)
+        {
+            for (int j=0;j<3;j++)
+            {
+                if (i==TT(i,j) || TT(i,j) == -1)
+                    Handle_MMatch(i,j)=0;
+                else
+                    Handle_MMatch(i,j) = MissMatchByLine(i,TT(i,j));
+            }
+        }
+    }
+
+};
+}
+
+
+template <typename DerivedV, typename DerivedF, typename DerivedO>
+IGL_INLINE void igl::line_field_missmatch(const Eigen::PlainObjectBase<DerivedV> &V,
+                                const Eigen::PlainObjectBase<DerivedF> &F,
+                                const Eigen::PlainObjectBase<DerivedV> &PD1,
+                                const bool isCombed,
+                                Eigen::PlainObjectBase<DerivedO> &missmatch)
+{
+    Eigen::PlainObjectBase<DerivedV> PD1_combed;
+    Eigen::PlainObjectBase<DerivedV> PD2_combed;
+
+    if (!isCombed)
+        igl::comb_line_field(V,F,PD1,PD1_combed);
+    else
+    {
+        PD1_combed = PD1;
+    }
+    Eigen::MatrixXd B1,B2,B3;
+    igl::local_basis(V,F,B1,B2,B3);
+    PD2_combed = igl::rotate_vectors(PD1_combed, Eigen::VectorXd::Constant(1,M_PI/2), B1, B2);
+    igl::MissMatchCalculatorLine<DerivedV, DerivedF, DerivedO> sf(V, F, PD1_combed, PD2_combed);
+    sf.calculateMissmatchLine(missmatch);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+#endif

include/igl/line_field_missmatch.h

@@ -0,0 +1,42 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 Nico Pietroni <nico.pietroni@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef IGL_LINE_FIELD_MISSMATCH_H
+#define IGL_LINE_FIELD_MISSMATCH_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Calculates the missmatch (integer), at each face edge, of a cross field defined on the mesh faces.
+  // The integer missmatch is a multiple of pi/2 that transforms the cross on one side of the edge to
+  // the cross on the other side. It represents the deviation from a Lie connection across the edge.
+
+  // Inputs:
+  //   V         #V by 3 eigen Matrix of mesh vertex 3D positions
+  //   F         #F by 3 eigen Matrix of face (quad) indices
+  //   PD1       #F by 3 eigen Matrix of the first per face cross field vector
+  //   PD2       #F by 3 eigen Matrix of the second per face cross field vector
+  //   isCombed  boolean, specifying whether the field is combed (i.e. matching has been precomputed.
+  //             If not, the field is combed first.
+  // Output:
+  //   Handle_MMatch    #F by 3 eigen Matrix containing the integer missmatch of the cross field
+  //                    across all face edges
+  //
+
+    template <typename DerivedV, typename DerivedF, typename DerivedO>
+    IGL_INLINE void line_field_missmatch(const Eigen::PlainObjectBase<DerivedV> &V,
+                                         const Eigen::PlainObjectBase<DerivedF> &F,
+                                         const Eigen::PlainObjectBase<DerivedV> &PD1,
+                                         const bool isCombed,
+                                         Eigen::PlainObjectBase<DerivedO> &missmatch);
+}
+#ifndef IGL_STATIC_LIBRARY
+#include "line_field_missmatch.cpp"
+#endif
+
+#endif

tutorial/506_FrameField/main.cpp

@@ -201,7 +201,7 @@ int main(int argc, char *argv[])
   igl::barycenter(V_deformed, F, B_deformed);
 
   // Find the closest crossfield to the deformed frame field
-  igl::frame_to_cross_field(V,F,FF1_deformed,FF2_deformed,X1_deformed);
+igl::frame_to_cross_field(V_deformed,F,FF1_deformed,FF2_deformed,X1_deformed);
 
   // Find a smooth crossfield that interpolates the deformed constraints
   MatrixXd bc_x(b.size(),3);