|
@@ -1,214 +1,33 @@
|
|
|
-//
|
|
|
-// mvc.h
|
|
|
-//
|
|
|
-// Created by Olga Diamanti on 9/11/11.
|
|
|
-// Copyright (c) 2011 ETH Zurich. All rights reserved.
|
|
|
-//
|
|
|
-
|
|
|
#ifndef IGL_MVC_H
|
|
|
#define IGL_MVC_H
|
|
|
|
|
|
-#include <Eigen/Core>
|
|
|
+#include "igl_inline.h"
|
|
|
+#include <Eigen/Dense>
|
|
|
|
|
|
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
|
|
|
-//
|
|
|
-
|
|
|
- inline void mvc(const Eigen::MatrixXd &V, const Eigen::MatrixXd &C, Eigen::MatrixXd &W);
|
|
|
+ // 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
|
|
|
+ //
|
|
|
+ // Known Bugs: implementation is listed as "Broken"
|
|
|
+ IGL_INLINE void mvc(
|
|
|
+ const Eigen::MatrixXd &V,
|
|
|
+ const Eigen::MatrixXd &C,
|
|
|
+ Eigen::MatrixXd &W);
|
|
|
|
|
|
}
|
|
|
|
|
|
-// Broken 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 (unsigned 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 (unsigned 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();
|
|
|
-}
|
|
|
+#ifdef IGL_HEADER_ONLY
|
|
|
+# include "mvc.cpp"
|
|
|
+#endif
|
|
|
|
|
|
#endif
|