WA cleanup

Former-commit-id: beeeb5f41c7e329c786693e12e7395af1f18c6be

mvc.h

@@ -0,0 +1,214 @@
+//
+//  mvc.h
+//
+//  Created by Olga Diamanti on 9/11/11.
+//  Copyright (c) 2011 ETH Zurich. All rights reserved.
+//
+
+#ifndef mvc_h
+#define mvc_h
+
+#include <Eigen/Core>
+
+namespace igl 
+{
+//   MVC - MEAN VALUE COORDINATES
+//  
+//   mvc(V,C,W)
+//  
+//   Inputs:
+//    V  #V x dim list of vertex positions (dim = 2 or dim = 3)
+//    C  #C x dim list of polygon vertex positions in counter-clockwise order  (dim = 2 or dim = 3)
+//  
+//   Outputs:
+//    W  weights, #V by #C matrix of weights
+//  
+    
+  void mvc(const Eigen::MatrixXd &V, const Eigen::MatrixXd &C, Eigen::MatrixXd &W);
+  
+}
+
+// Implementation
+inline void igl::mvc(const Eigen::MatrixXd &V, const Eigen::MatrixXd &C, Eigen::MatrixXd &W)
+{
+  
+  // at least three control points
+  assert(C.rows()>2);
+  
+  // dimension of points
+  assert(C.cols() == 3 || C.cols() == 2);
+  assert(V.cols() == 3 || V.cols() == 2);
+  
+  // number of polygon points
+  int num = C.rows();
+  
+  Eigen::MatrixXd V1,C1;
+  int i_prev, i_next;
+  
+  // check if either are 3D but really all z's are 0
+  bool V_flat = (V.cols() == 3) && (std::sqrt( (V.col(3)).dot(V.col(3)) ) < 1e-10);
+  bool C_flat = (C.cols() == 3) && (std::sqrt( (C.col(3)).dot(C.col(3)) ) < 1e-10);
+
+  // if both are essentially 2D then ignore z-coords
+  if((C.cols() == 2 || C_flat) && (V.cols() == 2 || V_flat))
+  {
+    // ignore z coordinate
+    V1 = V.block(0,0,V.rows(),2);
+    C1 = C.block(0,0,C.rows(),2);
+  }
+  else
+  {
+    // give dummy z coordinate to either mesh or poly
+    if(V.rows() == 2)
+    {
+      V1 = Eigen::MatrixXd(V.rows(),3);
+      V1.block(0,0,V.rows(),2) = V;
+    }
+    else
+      V1 = V;
+
+    if(C.rows() == 2)
+    {
+      C1 = Eigen::MatrixXd(C.rows(),3);
+      C1.block(0,0,C.rows(),2) = C;
+    }
+    else
+      C1 = C;
+
+    // check that C is planar
+    // average normal around poly corners
+
+    Eigen::Vector3d n = Eigen::Vector3d::Zero();
+    // take centroid as point on plane
+    Eigen::Vector3d p = Eigen::Vector3d::Zero();
+    for (int i = 0; i<num; ++i)
+    {
+      i_prev = (i>0)?(i-1):(num-1);
+      i_next = (i<num-1)?(i+1):0;
+      Eigen::Vector3d vnext = (C1.row(i_next) - C1.row(i)).transpose();
+      Eigen::Vector3d vprev = (C1.row(i_prev) - C1.row(i)).transpose();
+      n += vnext.cross(vprev);
+      p += C1.row(i);
+    }
+    p/=num;
+    n/=num;
+    // normalize n
+    n /= std::sqrt(n.dot(n));
+    
+    // check that poly is really coplanar
+    for (int i = 0; i<num; ++i)
+    {
+      double dist_to_plane_C = std::abs((C1.row(i)-p.transpose()).dot(n));
+      assert(dist_to_plane_C<1e-10);
+    }
+    
+    // check that poly is really coplanar
+    for (int i = 0; i<V1.rows(); ++i)
+    {
+      double dist_to_plane_V = std::abs((V1.row(i)-p.transpose()).dot(n));
+      if(dist_to_plane_V>1e-10)
+        std::cerr<<"Distance from V to plane of C is large..."<<std::endl;
+    }
+    
+    // change of basis
+    Eigen::Vector3d b1 = C1.row(1)-C1.row(0);
+    Eigen::Vector3d b2 = n.cross(b1);
+    // normalize basis rows
+    b1 /= std::sqrt(b1.dot(b1));
+    b2 /= std::sqrt(b2.dot(b2));
+    n /= std::sqrt(n.dot(n));
+    
+    //transpose of the basis matrix in the m-file
+    Eigen::Matrix3d basis = Eigen::Matrix3d::Zero();
+    basis.col(0) = b1;
+    basis.col(1) = b2;
+    basis.col(2) = n;
+    
+    // change basis of rows vectors by right multiplying with inverse of matrix
+    // with basis vectors as rows
+    Eigen::ColPivHouseholderQR<Eigen::Matrix3d> solver = basis.colPivHouseholderQr();
+    // Throw away coordinates in normal direction
+    V1 = solver.solve(V1.transpose()).transpose().block(0,0,V1.rows(),2);
+    C1 = solver.solve(C1.transpose()).transpose().block(0,0,C1.rows(),2);
+    
+  }
+  
+  // vectors from V to every C, where CmV(i,j,:) is the vector from domain
+  // vertex j to handle i
+  double EPSILON = 1e-10;
+  Eigen::MatrixXd WW = Eigen::MatrixXd(C1.rows(), V1.rows());
+  Eigen::MatrixXd dist_C_V (C1.rows(), V1.rows());
+  std::vector< std::pair<int,int> > on_corner(0);
+  std::vector< std::pair<int,int> > on_segment(0);
+  for (int i = 0; i<C1.rows(); ++i)
+  {
+    i_prev = (i>0)?(i-1):(num-1);
+    i_next = (i<num-1)?(i+1):0;
+    // distance from each corner in C to the next corner so that edge_length(i) 
+    // is the distance from C(i,:) to C(i+1,:) defined cyclically
+    double edge_length = std::sqrt((C1.row(i) - C1.row(i_next)).dot(C1.row(i) - C1.row(i_next)));
+    for (int j = 0; j<V1.rows(); ++j)
+    {
+      Eigen::VectorXd v = C1.row(i) - V1.row(j);
+      Eigen::VectorXd vnext = C1.row(i_next) - V1.row(j);
+      Eigen::VectorXd vprev = C1.row(i_prev) - V1.row(j);
+      // distance from V to every C, where dist_C_V(i,j) is the distance from domain
+      // vertex j to handle i
+      dist_C_V(i,j) = std::sqrt(v.dot(v));
+      double dist_C_V_next = std::sqrt(vnext.dot(vnext));
+      double a_prev = std::atan2(vprev[1],vprev[0]) - std::atan2(v[1],v[0]);
+      double a_next = std::atan2(v[1],v[0]) - std::atan2(vnext[1],vnext[0]);
+      // mean value coordinates
+      WW(i,j) = (std::tan(a_prev/2.0) + std::tan(a_next/2.0)) / dist_C_V(i,j);
+      
+      if (dist_C_V(i,j) < EPSILON)
+        on_corner.push_back(std::make_pair(j,i));
+      else
+        // only in case of no-corner (no need for checking for multiple segments afterwards --
+        // should only be on one segment (otherwise must be on a corner and we already
+        // handled that)
+        // domain vertex j is on the segment from i to i+1 if the distances from vj to
+        // pi and pi+1 are about 
+        if(abs((dist_C_V(i,j) + dist_C_V_next) / edge_length - 1) < EPSILON)
+          on_segment.push_back(std::make_pair(j,i));
+      
+    }
+  }
+  
+  // handle degenerate cases
+  // snap vertices close to corners
+  for (int i = 0; i<on_corner.size(); ++i)
+  {
+    int vi = on_corner[i].first;
+    int ci = on_corner[i].second;
+    for (int ii = 0; ii<C.rows(); ++ii)
+      WW(ii,vi) = (ii==ci)?1:0;
+  }
+  
+  // snap vertices close to segments
+  for (int i = 0; i<on_segment.size(); ++i)
+  {
+    int vi = on_segment[i].first;
+    int ci = on_segment[i].second;
+    int ci_next = (ci<num-1)?(ci+1):0;
+    for (int ii = 0; ii<C.rows(); ++ii)
+      if (ii == ci)
+        WW(ii,vi) = dist_C_V(ci_next,vi);
+      else
+      {
+        if ( ii == ci_next)
+          WW(ii,vi)  = dist_C_V(ci,vi);
+        else
+          WW(ii,vi) = 0;
+      }
+  }
+  
+  // normalize W
+  for (int i = 0; i<V.rows(); ++i)
+    WW.col(i) /= WW.col(i).sum();
+  
+  // we've made W transpose
+  W = WW.transpose();
+}
+
+#endif

orth.h

@@ -0,0 +1,55 @@
+//
+//  orth.h
+//
+//  Created by Olga Diamanti on 9/11/11.
+//  Copyright (c) 2011 ETH Zurich. All rights reserved.
+//
+
+#ifndef orth_h
+#define orth_h
+
+#include <Eigen/Core>
+
+namespace igl
+{
+//  ORTH   Orthogonalization.
+//     ORTH(A,Q) produces Q as an orthonormal basis for the range of A.
+//     That is, Q'*Q = I, the columns of Q span the same space as 
+//     the columns of A, and the number of columns of Q is the 
+//     rank of A.
+//  
+//  
+//   The algorithm  uses singular value decomposition, SVD, instead of orthogonal
+//   factorization, QR.  This doubles the computation time, but
+//   provides more reliable and consistent rank determination.
+//   Closely follows MATLAB implementation in orth.m
+//  
+  void orth(const Eigen::MatrixXd &A, Eigen::MatrixXd &Q);
+}
+
+// Implementation
+inline void igl::orth(const Eigen::MatrixXd &A, Eigen::MatrixXd &Q)
+{
+
+  //perform svd on A = U*S*V' (V is not computed and only the thin U is computed)
+  Eigen::JacobiSVD<Eigen::MatrixXd> svd(A, Eigen::ComputeThinU );
+  Eigen::MatrixXd U = svd.matrixU();
+  const Eigen::VectorXd S = svd.singularValues();
+  
+  //get rank of A
+  int m = A.rows();
+  int n = A.cols();
+  double tol = std::max(m,n) * S.maxCoeff() *  2.2204e-16;
+  int r = 0;
+  for (int i = 0; i < S.rows(); ++r,++i)
+  {
+    if (S[i] < tol)
+      break;
+  }
+  
+  //keep r first columns of U
+  Q = U.block(0,0,U.rows(),r);
+  
+}
+
+#endif