|
|
@@ -7,468 +7,465 @@
|
|
|
// obtain one at http://mozilla.org/MPL/2.0/.
|
|
|
|
|
|
#include <igl/n_polyvector_general.h>
|
|
|
-#include <igl/edge_topology.h>
|
|
|
-#include <igl/local_basis.h>
|
|
|
-#include <igl/nchoosek.h>
|
|
|
-#include <igl/slice.h>
|
|
|
-#include <igl/polyroots.h>
|
|
|
-#include <Eigen/Sparse>
|
|
|
-#include <Eigen/Geometry>
|
|
|
-#include <iostream>
|
|
|
-
|
|
|
-namespace igl {
|
|
|
- template <typename DerivedV, typename DerivedF>
|
|
|
- class GeneralPolyVectorFieldFinder
|
|
|
- {
|
|
|
- private:
|
|
|
- const Eigen::PlainObjectBase<DerivedV> &V;
|
|
|
- const Eigen::PlainObjectBase<DerivedF> &F; int numF;
|
|
|
- const int n;
|
|
|
-
|
|
|
- Eigen::MatrixXi EV; int numE;
|
|
|
- Eigen::MatrixXi F2E;
|
|
|
- Eigen::MatrixXi E2F;
|
|
|
- Eigen::VectorXd K;
|
|
|
-
|
|
|
- Eigen::VectorXi isBorderEdge;
|
|
|
- int numInteriorEdges;
|
|
|
- Eigen::Matrix<int,Eigen::Dynamic,2> E2F_int;
|
|
|
- Eigen::VectorXi indInteriorToFull;
|
|
|
- Eigen::VectorXi indFullToInterior;
|
|
|
-
|
|
|
-#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
|
|
|
- Eigen::PlainObjectBase<DerivedV> B1, B2, FN;
|
|
|
-
|
|
|
- IGL_INLINE void computek();
|
|
|
- IGL_INLINE void setFieldFromGeneralCoefficients(const std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> &coeffs,
|
|
|
- std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2> > &pv);
|
|
|
- IGL_INLINE void computeCoefficientLaplacian(int n, Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &D);
|
|
|
- IGL_INLINE void getGeneralCoeffConstraints(const Eigen::VectorXi &isConstrained,
|
|
|
- const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
|
|
|
- int k,
|
|
|
- const Eigen::VectorXi &rootsIndex,
|
|
|
- Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> &Ck);
|
|
|
- IGL_INLINE void precomputeInteriorEdges();
|
|
|
-
|
|
|
-
|
|
|
- IGL_INLINE void minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &Q,
|
|
|
- const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &f,
|
|
|
- const Eigen::VectorXi isConstrained,
|
|
|
- const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &xknown,
|
|
|
- Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &x);
|
|
|
-
|
|
|
- public:
|
|
|
- IGL_INLINE GeneralPolyVectorFieldFinder(const Eigen::PlainObjectBase<DerivedV> &_V,
|
|
|
- const Eigen::PlainObjectBase<DerivedF> &_F,
|
|
|
- const int &_n);
|
|
|
- IGL_INLINE bool solve(const Eigen::VectorXi &isConstrained,
|
|
|
- const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
|
|
|
- const Eigen::VectorXi &rootsIndex,
|
|
|
- Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &output);
|
|
|
-
|
|
|
- };
|
|
|
-}
|
|
|
-
|
|
|
-template<typename DerivedV, typename DerivedF>
|
|
|
-IGL_INLINE igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
|
|
|
- GeneralPolyVectorFieldFinder(const Eigen::PlainObjectBase<DerivedV> &_V,
|
|
|
- const Eigen::PlainObjectBase<DerivedF> &_F,
|
|
|
- const int &_n):
|
|
|
-V(_V),
|
|
|
-F(_F),
|
|
|
-numF(_F.rows()),
|
|
|
-n(_n)
|
|
|
-{
|
|
|
-
|
|
|
- igl::edge_topology(V,F,EV,F2E,E2F);
|
|
|
- numE = EV.rows();
|
|
|
-
|
|
|
-
|
|
|
- precomputeInteriorEdges();
|
|
|
-
|
|
|
- igl::local_basis(V,F,B1,B2,FN);
|
|
|
-
|
|
|
- computek();
|
|
|
-
|
|
|
-};
|
|
|
-
|
|
|
-
|
|
|
-template<typename DerivedV, typename DerivedF>
|
|
|
-IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
|
|
|
-precomputeInteriorEdges()
|
|
|
-{
|
|
|
- // Flag border edges
|
|
|
- numInteriorEdges = 0;
|
|
|
- isBorderEdge.setZero(numE,1);
|
|
|
- indFullToInterior = -1*Eigen::VectorXi::Ones(numE,1);
|
|
|
-
|
|
|
- for(unsigned i=0; i<numE; ++i)
|
|
|
- {
|
|
|
- if ((E2F(i,0) == -1) || ((E2F(i,1) == -1)))
|
|
|
- isBorderEdge[i] = 1;
|
|
|
- else
|
|
|
- {
|
|
|
- indFullToInterior[i] = numInteriorEdges;
|
|
|
- numInteriorEdges++;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- E2F_int.resize(numInteriorEdges, 2);
|
|
|
- indInteriorToFull.setZero(numInteriorEdges,1);
|
|
|
- int ii = 0;
|
|
|
- for (int k=0; k<numE; ++k)
|
|
|
- {
|
|
|
- if (isBorderEdge[k])
|
|
|
- continue;
|
|
|
- E2F_int.row(ii) = E2F.row(k);
|
|
|
- indInteriorToFull[ii] = k;
|
|
|
- ii++;
|
|
|
- }
|
|
|
-
|
|
|
-}
|
|
|
-
|
|
|
-
|
|
|
-template<typename DerivedV, typename DerivedF>
|
|
|
-IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
|
|
|
-minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &Q,
|
|
|
- const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &f,
|
|
|
- const Eigen::VectorXi isConstrained,
|
|
|
- const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &xknown,
|
|
|
- Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &x)
|
|
|
-{
|
|
|
- int N = Q.rows();
|
|
|
-
|
|
|
- int nc = xknown.rows();
|
|
|
- Eigen::VectorXi known; known.setZero(nc,1);
|
|
|
- Eigen::VectorXi unknown; unknown.setZero(N-nc,1);
|
|
|
-
|
|
|
- int indk = 0, indu = 0;
|
|
|
- for (int i = 0; i<N; ++i)
|
|
|
- if (isConstrained[i])
|
|
|
- {
|
|
|
- known[indk] = i;
|
|
|
- indk++;
|
|
|
- }
|
|
|
- else
|
|
|
- {
|
|
|
- unknown[indu] = i;
|
|
|
- indu++;
|
|
|
- }
|
|
|
-
|
|
|
- Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>> Quu, Quk;
|
|
|
-
|
|
|
- igl::slice(Q,unknown, unknown, Quu);
|
|
|
- igl::slice(Q,unknown, known, Quk);
|
|
|
-
|
|
|
-
|
|
|
- std::vector<typename Eigen::Triplet<std::complex<typename DerivedV::Scalar> > > tripletList;
|
|
|
-
|
|
|
- Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > fu(N-nc,1);
|
|
|
-
|
|
|
- igl::slice(f,unknown, Eigen::VectorXi::Zero(1,1), fu);
|
|
|
-
|
|
|
- Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > rhs = (Quk*xknown).sparseView()+.5*fu;
|
|
|
-
|
|
|
- Eigen::SparseLU< Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>>> solver;
|
|
|
- solver.compute(-Quu);
|
|
|
- if(solver.info()!=Eigen::Success)
|
|
|
- {
|
|
|
- std::cerr<<"Decomposition failed!"<<std::endl;
|
|
|
- return;
|
|
|
- }
|
|
|
- Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>> b = solver.solve(rhs);
|
|
|
- if(solver.info()!=Eigen::Success)
|
|
|
- {
|
|
|
- std::cerr<<"Solving failed!"<<std::endl;
|
|
|
- return;
|
|
|
- }
|
|
|
-
|
|
|
- indk = 0, indu = 0;
|
|
|
- x.setZero(N,1);
|
|
|
- for (int i = 0; i<N; ++i)
|
|
|
- if (isConstrained[i])
|
|
|
- x[i] = xknown[indk++];
|
|
|
- else
|
|
|
- x[i] = b.coeff(indu++,0);
|
|
|
-
|
|
|
-}
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
-template<typename DerivedV, typename DerivedF>
|
|
|
-IGL_INLINE bool igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
|
|
|
- solve(const Eigen::VectorXi &isConstrained,
|
|
|
- const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
|
|
|
- const Eigen::VectorXi &rootsIndex,
|
|
|
- Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &output)
|
|
|
-{
|
|
|
-
|
|
|
- // polynomial is of the form:
|
|
|
- // z^(2n) +
|
|
|
- // -c[0]z^(2n-1) +
|
|
|
- // c[1]z^(2n-2) +
|
|
|
- // -c[2]z^(2n-3) +
|
|
|
- // ... +
|
|
|
- // (-1)^n c[n-1]
|
|
|
-
|
|
|
- std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> coeffs(n,Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>::Zero(numF, 1));
|
|
|
-
|
|
|
- for (int i =0; i<n; ++i)
|
|
|
- {
|
|
|
- int degree = i+1;
|
|
|
-
|
|
|
- Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> Ck;
|
|
|
- getGeneralCoeffConstraints(isConstrained,
|
|
|
- cfW,
|
|
|
- i,
|
|
|
- rootsIndex,
|
|
|
- Ck);
|
|
|
-
|
|
|
- Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > DD;
|
|
|
- computeCoefficientLaplacian(degree, DD);
|
|
|
- Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > f; f.resize(numF,1);
|
|
|
-
|
|
|
- if (isConstrained.sum() == numF)
|
|
|
- coeffs[i] = Ck;
|
|
|
- else
|
|
|
- minQuadWithKnownMini(DD, f, isConstrained, Ck, coeffs[i]);
|
|
|
- }
|
|
|
-
|
|
|
- std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2> > pv;
|
|
|
- setFieldFromGeneralCoefficients(coeffs, pv);
|
|
|
-
|
|
|
- output.setZero(numF,3*n);
|
|
|
- for (int fi=0; fi<numF; ++fi)
|
|
|
- {
|
|
|
- const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b1 = B1.row(fi);
|
|
|
- const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b2 = B2.row(fi);
|
|
|
- for (int i=0; i<n; ++i)
|
|
|
- output.block(fi,3*i, 1, 3) = pv[i](fi,0)*b1 + pv[i](fi,1)*b2;
|
|
|
- }
|
|
|
- return true;
|
|
|
-}
|
|
|
-
|
|
|
-template<typename DerivedV, typename DerivedF>
|
|
|
-IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::setFieldFromGeneralCoefficients(const std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> &coeffs,
|
|
|
- std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2>> &pv)
|
|
|
-{
|
|
|
- pv.assign(n, Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2>::Zero(numF, 2));
|
|
|
- for (int i = 0; i <numF; ++i)
|
|
|
- {
|
|
|
-
|
|
|
- // poly coefficients: 1, 0, -Acoeff, 0, Bcoeff
|
|
|
- // matlab code from roots (given there are no trailing zeros in the polynomial coefficients)
|
|
|
- Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> polyCoeff;
|
|
|
- polyCoeff.setZero(n+1,1);
|
|
|
- polyCoeff[0] = 1.;
|
|
|
- int sign = 1;
|
|
|
- for (int k =0; k<n; ++k)
|
|
|
- {
|
|
|
- sign = -sign;
|
|
|
- int degree = k+1;
|
|
|
- polyCoeff[degree] = (1.*sign)*coeffs[k](i);
|
|
|
- }
|
|
|
-
|
|
|
- Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> roots;
|
|
|
- igl::polyRoots<std::complex<typename DerivedV::Scalar>, typename DerivedV::Scalar >(polyCoeff,roots);
|
|
|
- for (int k=0; k<n; ++k)
|
|
|
- {
|
|
|
- pv[k](i,0) = real(roots[k]);
|
|
|
- pv[k](i,1) = imag(roots[k]);
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-}
|
|
|
-
|
|
|
-
|
|
|
-template<typename DerivedV, typename DerivedF>
|
|
|
-IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::computeCoefficientLaplacian(int n, Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &D)
|
|
|
-{
|
|
|
- std::vector<Eigen::Triplet<std::complex<typename DerivedV::Scalar> >> tripletList;
|
|
|
-
|
|
|
- // For every non-border edge
|
|
|
- for (unsigned eid=0; eid<numE; ++eid)
|
|
|
- {
|
|
|
- if (!isBorderEdge[eid])
|
|
|
- {
|
|
|
- int fid0 = E2F(eid,0);
|
|
|
- int fid1 = E2F(eid,1);
|
|
|
-
|
|
|
- tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid0,
|
|
|
- fid0,
|
|
|
- std::complex<typename DerivedV::Scalar>(1.)));
|
|
|
- tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid1,
|
|
|
- fid1,
|
|
|
- std::complex<typename DerivedV::Scalar>(1.)));
|
|
|
- tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid0,
|
|
|
- fid1,
|
|
|
- -1.*std::polar(1.,-1.*n*K[eid])));
|
|
|
- tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid1,
|
|
|
- fid0,
|
|
|
- -1.*std::polar(1.,1.*n*K[eid])));
|
|
|
-
|
|
|
- }
|
|
|
- }
|
|
|
- D.resize(numF,numF);
|
|
|
- D.setFromTriplets(tripletList.begin(), tripletList.end());
|
|
|
-
|
|
|
-
|
|
|
-}
|
|
|
-
|
|
|
-//this gives the coefficients without the (-1)^k that multiplies them
|
|
|
-template<typename DerivedV, typename DerivedF>
|
|
|
-IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::getGeneralCoeffConstraints(const Eigen::VectorXi &isConstrained,
|
|
|
- const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
|
|
|
- int k,
|
|
|
- const Eigen::VectorXi &rootsIndex,
|
|
|
- Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> &Ck)
|
|
|
-{
|
|
|
- int numConstrained = isConstrained.sum();
|
|
|
- Ck.resize(numConstrained,1);
|
|
|
- // int n = rootsIndex.cols();
|
|
|
-
|
|
|
- Eigen::MatrixXi allCombs;
|
|
|
- {
|
|
|
- Eigen::VectorXi V = Eigen::VectorXi::LinSpaced(n,0,n-1);
|
|
|
- igl::nchoosek(V,k+1,allCombs);
|
|
|
- }
|
|
|
-
|
|
|
- int ind = 0;
|
|
|
- for (int fi = 0; fi <numF; ++fi)
|
|
|
- {
|
|
|
- const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b1 = B1.row(fi);
|
|
|
- const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b2 = B2.row(fi);
|
|
|
- if(isConstrained[fi])
|
|
|
- {
|
|
|
- std::complex<typename DerivedV::Scalar> ck(0);
|
|
|
-
|
|
|
- for (int j = 0; j < allCombs.rows(); ++j)
|
|
|
- {
|
|
|
- std::complex<typename DerivedV::Scalar> tk(1.);
|
|
|
- //collect products
|
|
|
- for (int i = 0; i < allCombs.cols(); ++i)
|
|
|
- {
|
|
|
- int index = allCombs(j,i);
|
|
|
-
|
|
|
- int ri = rootsIndex[index];
|
|
|
- Eigen::Matrix<typename DerivedV::Scalar, 1, 3> w;
|
|
|
- if (ri>0)
|
|
|
- w = cfW.block(fi,3*(ri-1),1,3);
|
|
|
- else
|
|
|
- w = -cfW.block(fi,3*(-ri-1),1,3);
|
|
|
- typename DerivedV::Scalar w0 = w.dot(b1);
|
|
|
- typename DerivedV::Scalar w1 = w.dot(b2);
|
|
|
- std::complex<typename DerivedV::Scalar> u(w0,w1);
|
|
|
- tk*= u;
|
|
|
- }
|
|
|
- //collect sum
|
|
|
- ck += tk;
|
|
|
- }
|
|
|
- Ck(ind) = ck;
|
|
|
- ind ++;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
-}
|
|
|
-
|
|
|
-template<typename DerivedV, typename DerivedF>
|
|
|
-IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::computek()
|
|
|
-{
|
|
|
- K.setZero(numE);
|
|
|
- // For every non-border edge
|
|
|
- for (unsigned eid=0; eid<numE; ++eid)
|
|
|
- {
|
|
|
- if (!isBorderEdge[eid])
|
|
|
- {
|
|
|
- int fid0 = E2F(eid,0);
|
|
|
- int fid1 = E2F(eid,1);
|
|
|
-
|
|
|
- Eigen::Matrix<typename DerivedV::Scalar, 1, 3> N0 = FN.row(fid0);
|
|
|
- Eigen::Matrix<typename DerivedV::Scalar, 1, 3> N1 = FN.row(fid1);
|
|
|
-
|
|
|
- // find common edge on triangle 0 and 1
|
|
|
- int fid0_vc = -1;
|
|
|
- int fid1_vc = -1;
|
|
|
- for (unsigned i=0;i<3;++i)
|
|
|
- {
|
|
|
- if (F2E(fid0,i) == eid)
|
|
|
- fid0_vc = i;
|
|
|
- if (F2E(fid1,i) == eid)
|
|
|
- fid1_vc = i;
|
|
|
- }
|
|
|
- assert(fid0_vc != -1);
|
|
|
- assert(fid1_vc != -1);
|
|
|
-
|
|
|
- Eigen::Matrix<typename DerivedV::Scalar, 1, 3> common_edge = V.row(F(fid0,(fid0_vc+1)%3)) - V.row(F(fid0,fid0_vc));
|
|
|
- common_edge.normalize();
|
|
|
-
|
|
|
- // Map the two triangles in a new space where the common edge is the x axis and the N0 the z axis
|
|
|
- Eigen::Matrix<typename DerivedV::Scalar, 3, 3> P;
|
|
|
- Eigen::Matrix<typename DerivedV::Scalar, 1, 3> o = V.row(F(fid0,fid0_vc));
|
|
|
- Eigen::Matrix<typename DerivedV::Scalar, 1, 3> tmp = -N0.cross(common_edge);
|
|
|
- P << common_edge, tmp, N0;
|
|
|
-// P.transposeInPlace();
|
|
|
-
|
|
|
-
|
|
|
- Eigen::Matrix<typename DerivedV::Scalar, 3, 3> V0;
|
|
|
- V0.row(0) = V.row(F(fid0,0)) -o;
|
|
|
- V0.row(1) = V.row(F(fid0,1)) -o;
|
|
|
- V0.row(2) = V.row(F(fid0,2)) -o;
|
|
|
-
|
|
|
- V0 = (P*V0.transpose()).transpose();
|
|
|
-
|
|
|
-// assert(V0(0,2) < 1e-10);
|
|
|
-// assert(V0(1,2) < 1e-10);
|
|
|
-// assert(V0(2,2) < 1e-10);
|
|
|
-
|
|
|
- Eigen::Matrix<typename DerivedV::Scalar, 3, 3> V1;
|
|
|
- V1.row(0) = V.row(F(fid1,0)) -o;
|
|
|
- V1.row(1) = V.row(F(fid1,1)) -o;
|
|
|
- V1.row(2) = V.row(F(fid1,2)) -o;
|
|
|
- V1 = (P*V1.transpose()).transpose();
|
|
|
-
|
|
|
-// assert(V1(fid1_vc,2) < 10e-10);
|
|
|
-// assert(V1((fid1_vc+1)%3,2) < 10e-10);
|
|
|
-
|
|
|
- // compute rotation R such that R * N1 = N0
|
|
|
- // i.e. map both triangles to the same plane
|
|
|
- double alpha = -atan2(V1((fid1_vc+2)%3,2),V1((fid1_vc+2)%3,1));
|
|
|
-
|
|
|
- Eigen::Matrix<typename DerivedV::Scalar, 3, 3> R;
|
|
|
- R << 1, 0, 0,
|
|
|
- 0, cos(alpha), -sin(alpha) ,
|
|
|
- 0, sin(alpha), cos(alpha);
|
|
|
- V1 = (R*V1.transpose()).transpose();
|
|
|
-
|
|
|
-// assert(V1(0,2) < 1e-10);
|
|
|
-// assert(V1(1,2) < 1e-10);
|
|
|
-// assert(V1(2,2) < 1e-10);
|
|
|
-
|
|
|
- // measure the angle between the reference frames
|
|
|
- // k_ij is the angle between the triangle on the left and the one on the right
|
|
|
- Eigen::Matrix<typename DerivedV::Scalar, 1, 3> ref0 = V0.row(1) - V0.row(0);
|
|
|
- Eigen::Matrix<typename DerivedV::Scalar, 1, 3> ref1 = V1.row(1) - V1.row(0);
|
|
|
-
|
|
|
- ref0.normalize();
|
|
|
- ref1.normalize();
|
|
|
-
|
|
|
- double ktemp = atan2(ref1(1),ref1(0)) - atan2(ref0(1),ref0(0));
|
|
|
-
|
|
|
- // just to be sure, rotate ref0 using angle ktemp...
|
|
|
- Eigen::Matrix<typename DerivedV::Scalar, 2, 2> R2;
|
|
|
- R2 << cos(ktemp), -sin(ktemp), sin(ktemp), cos(ktemp);
|
|
|
-
|
|
|
- Eigen::Matrix<typename DerivedV::Scalar, 1, 2> tmp1 = R2*(ref0.head(2)).transpose();
|
|
|
-
|
|
|
-// assert(tmp1(0) - ref1(0) < 1e-10);
|
|
|
-// assert(tmp1(1) - ref1(1) < 1e-10);
|
|
|
-
|
|
|
- K[eid] = ktemp;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-}
|
|
|
+//#include <igl/edge_topology.h>
|
|
|
+//#include <igl/local_basis.h>
|
|
|
+//#include <igl/nchoosek.h>
|
|
|
+//#include <igl/slice.h>
|
|
|
+//#include <igl/polyroots.h>
|
|
|
+//#include <Eigen/Sparse>
|
|
|
+//#include <Eigen/Geometry>
|
|
|
+//#include <iostream>
|
|
|
+//#include <iostream>
|
|
|
+//
|
|
|
+//namespace igl {
|
|
|
+// template <typename DerivedV, typename DerivedF>
|
|
|
+// class GeneralPolyVectorFieldFinder
|
|
|
+// {
|
|
|
+// private:
|
|
|
+// const Eigen::PlainObjectBase<DerivedV> &V;
|
|
|
+// const Eigen::PlainObjectBase<DerivedF> &F; int numF;
|
|
|
+// const int n;
|
|
|
+//
|
|
|
+// Eigen::MatrixXi EV; int numE;
|
|
|
+// Eigen::MatrixXi F2E;
|
|
|
+// Eigen::MatrixXi E2F;
|
|
|
+// Eigen::VectorXd K;
|
|
|
+//
|
|
|
+// Eigen::VectorXi isBorderEdge;
|
|
|
+// int numInteriorEdges;
|
|
|
+// Eigen::Matrix<int,Eigen::Dynamic,2> E2F_int;
|
|
|
+// Eigen::VectorXi indInteriorToFull;
|
|
|
+// Eigen::VectorXi indFullToInterior;
|
|
|
+//
|
|
|
+////#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
|
|
|
+// Eigen::PlainObjectBase<DerivedV> B1, B2, FN;
|
|
|
+//
|
|
|
+// IGL_INLINE void computek();
|
|
|
+// IGL_INLINE void setFieldFromGeneralCoefficients(const std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> &coeffs,
|
|
|
+// std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2> > &pv);
|
|
|
+// IGL_INLINE void computeCoefficientLaplacian(int n, Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &D);
|
|
|
+// IGL_INLINE void getGeneralCoeffConstraints(const Eigen::VectorXi &isConstrained,
|
|
|
+// const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
|
|
|
+// int k,
|
|
|
+// const Eigen::VectorXi &rootsIndex,
|
|
|
+// Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> &Ck);
|
|
|
+// IGL_INLINE void precomputeInteriorEdges();
|
|
|
+//
|
|
|
+//
|
|
|
+// IGL_INLINE void minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &Q,
|
|
|
+// const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &f,
|
|
|
+// const Eigen::VectorXi isConstrained,
|
|
|
+// const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &xknown,
|
|
|
+// Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &x);
|
|
|
+//
|
|
|
+// public:
|
|
|
+// IGL_INLINE GeneralPolyVectorFieldFinder(const Eigen::PlainObjectBase<DerivedV> &_V,
|
|
|
+// const Eigen::PlainObjectBase<DerivedF> &_F,
|
|
|
+// const int &_n);
|
|
|
+// IGL_INLINE bool solve(const Eigen::VectorXi &isConstrained,
|
|
|
+// const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
|
|
|
+// const Eigen::VectorXi &rootsIndex,
|
|
|
+// Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &output);
|
|
|
+//
|
|
|
+// };
|
|
|
+//}
|
|
|
+//
|
|
|
+//template<typename DerivedV, typename DerivedF>
|
|
|
+//IGL_INLINE igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
|
|
|
+// GeneralPolyVectorFieldFinder(const Eigen::PlainObjectBase<DerivedV> &_V,
|
|
|
+// const Eigen::PlainObjectBase<DerivedF> &_F,
|
|
|
+// const int &_n):
|
|
|
+//V(_V),
|
|
|
+//F(_F),
|
|
|
+//numF(_F.rows()),
|
|
|
+//n(_n)
|
|
|
+//{
|
|
|
+//
|
|
|
+// igl::edge_topology(V,F,EV,F2E,E2F);
|
|
|
+// numE = EV.rows();
|
|
|
+//
|
|
|
+//
|
|
|
+// precomputeInteriorEdges();
|
|
|
+//
|
|
|
+// igl::local_basis(V,F,B1,B2,FN);
|
|
|
+//
|
|
|
+// computek();
|
|
|
+//
|
|
|
+//};
|
|
|
+//
|
|
|
+//
|
|
|
+//template<typename DerivedV, typename DerivedF>
|
|
|
+//IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
|
|
|
+//precomputeInteriorEdges()
|
|
|
+//{
|
|
|
+// // Flag border edges
|
|
|
+// numInteriorEdges = 0;
|
|
|
+// isBorderEdge.setZero(numE,1);
|
|
|
+// indFullToInterior = -1*Eigen::VectorXi::Ones(numE,1);
|
|
|
+//
|
|
|
+// for(unsigned i=0; i<numE; ++i)
|
|
|
+// {
|
|
|
+// if ((E2F(i,0) == -1) || ((E2F(i,1) == -1)))
|
|
|
+// isBorderEdge[i] = 1;
|
|
|
+// else
|
|
|
+// {
|
|
|
+// indFullToInterior[i] = numInteriorEdges;
|
|
|
+// numInteriorEdges++;
|
|
|
+// }
|
|
|
+// }
|
|
|
+//
|
|
|
+// E2F_int.resize(numInteriorEdges, 2);
|
|
|
+// indInteriorToFull.setZero(numInteriorEdges,1);
|
|
|
+// int ii = 0;
|
|
|
+// for (int k=0; k<numE; ++k)
|
|
|
+// {
|
|
|
+// if (isBorderEdge[k])
|
|
|
+// continue;
|
|
|
+// E2F_int.row(ii) = E2F.row(k);
|
|
|
+// indInteriorToFull[ii] = k;
|
|
|
+// ii++;
|
|
|
+// }
|
|
|
+//
|
|
|
+//}
|
|
|
+//
|
|
|
+//
|
|
|
+//template<typename DerivedV, typename DerivedF>
|
|
|
+//IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
|
|
|
+//minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &Q,
|
|
|
+// const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &f,
|
|
|
+// const Eigen::VectorXi isConstrained,
|
|
|
+// const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &xknown,
|
|
|
+// Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &x)
|
|
|
+//{
|
|
|
+// int N = Q.rows();
|
|
|
+//
|
|
|
+// int nc = xknown.rows();
|
|
|
+// Eigen::VectorXi known; known.setZero(nc,1);
|
|
|
+// Eigen::VectorXi unknown; unknown.setZero(N-nc,1);
|
|
|
+//
|
|
|
+// int indk = 0, indu = 0;
|
|
|
+// for (int i = 0; i<N; ++i)
|
|
|
+// if (isConstrained[i])
|
|
|
+// {
|
|
|
+// known[indk] = i;
|
|
|
+// indk++;
|
|
|
+// }
|
|
|
+// else
|
|
|
+// {
|
|
|
+// unknown[indu] = i;
|
|
|
+// indu++;
|
|
|
+// }
|
|
|
+//
|
|
|
+// Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>> Quu, Quk;
|
|
|
+//
|
|
|
+// igl::slice(Q,unknown, unknown, Quu);
|
|
|
+// igl::slice(Q,unknown, known, Quk);
|
|
|
+//
|
|
|
+//
|
|
|
+// std::vector<typename Eigen::Triplet<std::complex<typename DerivedV::Scalar> > > tripletList;
|
|
|
+//
|
|
|
+// Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > fu(N-nc,1);
|
|
|
+//
|
|
|
+// igl::slice(f,unknown, Eigen::VectorXi::Zero(1,1), fu);
|
|
|
+//
|
|
|
+// Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > rhs = (Quk*xknown).sparseView()+.5*fu;
|
|
|
+//
|
|
|
+// Eigen::SparseLU< Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>>> solver;
|
|
|
+// solver.compute(-Quu);
|
|
|
+// if(solver.info()!=Eigen::Success)
|
|
|
+// {
|
|
|
+// std::cerr<<"Decomposition failed!"<<std::endl;
|
|
|
+// return;
|
|
|
+// }
|
|
|
+// Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>> b = solver.solve(rhs);
|
|
|
+// if(solver.info()!=Eigen::Success)
|
|
|
+// {
|
|
|
+// std::cerr<<"Solving failed!"<<std::endl;
|
|
|
+// return;
|
|
|
+// }
|
|
|
+//
|
|
|
+// indk = 0, indu = 0;
|
|
|
+// x.setZero(N,1);
|
|
|
+// for (int i = 0; i<N; ++i)
|
|
|
+// if (isConstrained[i])
|
|
|
+// x[i] = xknown[indk++];
|
|
|
+// else
|
|
|
+// x[i] = b.coeff(indu++,0);
|
|
|
+//
|
|
|
+//}
|
|
|
+//
|
|
|
+//
|
|
|
+//
|
|
|
+//template<typename DerivedV, typename DerivedF>
|
|
|
+//IGL_INLINE bool igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
|
|
|
+// solve(const Eigen::VectorXi &isConstrained,
|
|
|
+// const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
|
|
|
+// const Eigen::VectorXi &rootsIndex,
|
|
|
+// Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &output)
|
|
|
+//{
|
|
|
+// std::cerr << "This code is broken!" << std::endl;
|
|
|
+// exit(1);
|
|
|
+// // polynomial is of the form:
|
|
|
+// // z^(2n) +
|
|
|
+// // -c[0]z^(2n-1) +
|
|
|
+// // c[1]z^(2n-2) +
|
|
|
+// // -c[2]z^(2n-3) +
|
|
|
+// // ... +
|
|
|
+// // (-1)^n c[n-1]
|
|
|
+// //std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> coeffs(n,Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>::Zero(numF, 1));
|
|
|
+// //for (int i =0; i<n; ++i)
|
|
|
+// //{
|
|
|
+// // int degree = i+1;
|
|
|
+// // Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> Ck;
|
|
|
+// // getGeneralCoeffConstraints(isConstrained,
|
|
|
+// // cfW,
|
|
|
+// // i,
|
|
|
+// // rootsIndex,
|
|
|
+// // Ck);
|
|
|
+//
|
|
|
+// // Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > DD;
|
|
|
+// // computeCoefficientLaplacian(degree, DD);
|
|
|
+// // Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > f; f.resize(numF,1);
|
|
|
+//
|
|
|
+// // if (isConstrained.sum() == numF)
|
|
|
+// // coeffs[i] = Ck;
|
|
|
+// // else
|
|
|
+// // minQuadWithKnownMini(DD, f, isConstrained, Ck, coeffs[i]);
|
|
|
+// //}
|
|
|
+// //std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2> > pv;
|
|
|
+// //setFieldFromGeneralCoefficients(coeffs, pv);
|
|
|
+// //output.setZero(numF,3*n);
|
|
|
+// //for (int fi=0; fi<numF; ++fi)
|
|
|
+// //{
|
|
|
+// // const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b1 = B1.row(fi);
|
|
|
+// // const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b2 = B2.row(fi);
|
|
|
+// // for (int i=0; i<n; ++i)
|
|
|
+// // output.block(fi,3*i, 1, 3) = pv[i](fi,0)*b1 + pv[i](fi,1)*b2;
|
|
|
+// //}
|
|
|
+// return true;
|
|
|
+//}
|
|
|
+//
|
|
|
+//template<typename DerivedV, typename DerivedF>
|
|
|
+//IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::setFieldFromGeneralCoefficients(const std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> &coeffs,
|
|
|
+// std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2>> &pv)
|
|
|
+//{
|
|
|
+// pv.assign(n, Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2>::Zero(numF, 2));
|
|
|
+// for (int i = 0; i <numF; ++i)
|
|
|
+// {
|
|
|
+//
|
|
|
+// // poly coefficients: 1, 0, -Acoeff, 0, Bcoeff
|
|
|
+// // matlab code from roots (given there are no trailing zeros in the polynomial coefficients)
|
|
|
+// Eigen::Matrix<std::complex<DerivedV::Scalar>, Eigen::Dynamic,1> polyCoeff;
|
|
|
+// polyCoeff.setZero(n+1,1);
|
|
|
+// polyCoeff[0] = 1.;
|
|
|
+// int sign = 1;
|
|
|
+// for (int k =0; k<n; ++k)
|
|
|
+// {
|
|
|
+// sign = -sign;
|
|
|
+// int degree = k+1;
|
|
|
+// polyCoeff[degree] = (1.*sign)*coeffs[k](i);
|
|
|
+// }
|
|
|
+//
|
|
|
+// Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> roots;
|
|
|
+// igl::polyRoots<std::complex<typename DerivedV::Scalar>, typename DerivedV::Scalar >(polyCoeff,roots);
|
|
|
+// for (int k=0; k<n; ++k)
|
|
|
+// {
|
|
|
+// pv[k](i,0) = real(roots[k]);
|
|
|
+// pv[k](i,1) = imag(roots[k]);
|
|
|
+// }
|
|
|
+// }
|
|
|
+//
|
|
|
+//}
|
|
|
+//
|
|
|
+//
|
|
|
+//template<typename DerivedV, typename DerivedF>
|
|
|
+//IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::computeCoefficientLaplacian(int n, Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &D)
|
|
|
+//{
|
|
|
+// std::vector<Eigen::Triplet<std::complex<typename DerivedV::Scalar> >> tripletList;
|
|
|
+//
|
|
|
+// // For every non-border edge
|
|
|
+// for (unsigned eid=0; eid<numE; ++eid)
|
|
|
+// {
|
|
|
+// if (!isBorderEdge[eid])
|
|
|
+// {
|
|
|
+// int fid0 = E2F(eid,0);
|
|
|
+// int fid1 = E2F(eid,1);
|
|
|
+//
|
|
|
+// tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid0,
|
|
|
+// fid0,
|
|
|
+// std::complex<typename DerivedV::Scalar>(1.)));
|
|
|
+// tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid1,
|
|
|
+// fid1,
|
|
|
+// std::complex<typename DerivedV::Scalar>(1.)));
|
|
|
+// tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid0,
|
|
|
+// fid1,
|
|
|
+// -1.*std::polar(1.,-1.*n*K[eid])));
|
|
|
+// tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid1,
|
|
|
+// fid0,
|
|
|
+// -1.*std::polar(1.,1.*n*K[eid])));
|
|
|
+//
|
|
|
+// }
|
|
|
+// }
|
|
|
+// D.resize(numF,numF);
|
|
|
+// D.setFromTriplets(tripletList.begin(), tripletList.end());
|
|
|
+//
|
|
|
+//
|
|
|
+//}
|
|
|
+//
|
|
|
+////this gives the coefficients without the (-1)^k that multiplies them
|
|
|
+//template<typename DerivedV, typename DerivedF>
|
|
|
+//IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::getGeneralCoeffConstraints(const Eigen::VectorXi &isConstrained,
|
|
|
+// const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
|
|
|
+// int k,
|
|
|
+// const Eigen::VectorXi &rootsIndex,
|
|
|
+// Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> &Ck)
|
|
|
+//{
|
|
|
+// int numConstrained = isConstrained.sum();
|
|
|
+// Ck.resize(numConstrained,1);
|
|
|
+// // int n = rootsIndex.cols();
|
|
|
+//
|
|
|
+// Eigen::MatrixXi allCombs;
|
|
|
+// {
|
|
|
+// Eigen::VectorXi V = Eigen::VectorXi::LinSpaced(n,0,n-1);
|
|
|
+// igl::nchoosek(V,k+1,allCombs);
|
|
|
+// }
|
|
|
+//
|
|
|
+// int ind = 0;
|
|
|
+// for (int fi = 0; fi <numF; ++fi)
|
|
|
+// {
|
|
|
+// const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b1 = B1.row(fi);
|
|
|
+// const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b2 = B2.row(fi);
|
|
|
+// if(isConstrained[fi])
|
|
|
+// {
|
|
|
+// std::complex<typename DerivedV::Scalar> ck(0);
|
|
|
+//
|
|
|
+// for (int j = 0; j < allCombs.rows(); ++j)
|
|
|
+// {
|
|
|
+// std::complex<typename DerivedV::Scalar> tk(1.);
|
|
|
+// //collect products
|
|
|
+// for (int i = 0; i < allCombs.cols(); ++i)
|
|
|
+// {
|
|
|
+// int index = allCombs(j,i);
|
|
|
+//
|
|
|
+// int ri = rootsIndex[index];
|
|
|
+// Eigen::Matrix<typename DerivedV::Scalar, 1, 3> w;
|
|
|
+// if (ri>0)
|
|
|
+// w = cfW.block(fi,3*(ri-1),1,3);
|
|
|
+// else
|
|
|
+// w = -cfW.block(fi,3*(-ri-1),1,3);
|
|
|
+// typename DerivedV::Scalar w0 = w.dot(b1);
|
|
|
+// typename DerivedV::Scalar w1 = w.dot(b2);
|
|
|
+// std::complex<typename DerivedV::Scalar> u(w0,w1);
|
|
|
+// tk*= u;
|
|
|
+// }
|
|
|
+// //collect sum
|
|
|
+// ck += tk;
|
|
|
+// }
|
|
|
+// Ck(ind) = ck;
|
|
|
+// ind ++;
|
|
|
+// }
|
|
|
+// }
|
|
|
+//
|
|
|
+//
|
|
|
+//}
|
|
|
+//
|
|
|
+//template<typename DerivedV, typename DerivedF>
|
|
|
+//IGL_INLINE void igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::computek()
|
|
|
+//{
|
|
|
+// K.setZero(numE);
|
|
|
+// // For every non-border edge
|
|
|
+// for (unsigned eid=0; eid<numE; ++eid)
|
|
|
+// {
|
|
|
+// if (!isBorderEdge[eid])
|
|
|
+// {
|
|
|
+// int fid0 = E2F(eid,0);
|
|
|
+// int fid1 = E2F(eid,1);
|
|
|
+//
|
|
|
+// Eigen::Matrix<typename DerivedV::Scalar, 1, 3> N0 = FN.row(fid0);
|
|
|
+// Eigen::Matrix<typename DerivedV::Scalar, 1, 3> N1 = FN.row(fid1);
|
|
|
+//
|
|
|
+// // find common edge on triangle 0 and 1
|
|
|
+// int fid0_vc = -1;
|
|
|
+// int fid1_vc = -1;
|
|
|
+// for (unsigned i=0;i<3;++i)
|
|
|
+// {
|
|
|
+// if (F2E(fid0,i) == eid)
|
|
|
+// fid0_vc = i;
|
|
|
+// if (F2E(fid1,i) == eid)
|
|
|
+// fid1_vc = i;
|
|
|
+// }
|
|
|
+// assert(fid0_vc != -1);
|
|
|
+// assert(fid1_vc != -1);
|
|
|
+//
|
|
|
+// Eigen::Matrix<typename DerivedV::Scalar, 1, 3> common_edge = V.row(F(fid0,(fid0_vc+1)%3)) - V.row(F(fid0,fid0_vc));
|
|
|
+// common_edge.normalize();
|
|
|
+//
|
|
|
+// // Map the two triangles in a new space where the common edge is the x axis and the N0 the z axis
|
|
|
+// Eigen::Matrix<typename DerivedV::Scalar, 3, 3> P;
|
|
|
+// Eigen::Matrix<typename DerivedV::Scalar, 1, 3> o = V.row(F(fid0,fid0_vc));
|
|
|
+// Eigen::Matrix<typename DerivedV::Scalar, 1, 3> tmp = -N0.cross(common_edge);
|
|
|
+// P << common_edge, tmp, N0;
|
|
|
+//// P.transposeInPlace();
|
|
|
+//
|
|
|
+//
|
|
|
+// Eigen::Matrix<typename DerivedV::Scalar, 3, 3> V0;
|
|
|
+// V0.row(0) = V.row(F(fid0,0)) -o;
|
|
|
+// V0.row(1) = V.row(F(fid0,1)) -o;
|
|
|
+// V0.row(2) = V.row(F(fid0,2)) -o;
|
|
|
+//
|
|
|
+// V0 = (P*V0.transpose()).transpose();
|
|
|
+//
|
|
|
+//// assert(V0(0,2) < 1e-10);
|
|
|
+//// assert(V0(1,2) < 1e-10);
|
|
|
+//// assert(V0(2,2) < 1e-10);
|
|
|
+//
|
|
|
+// Eigen::Matrix<typename DerivedV::Scalar, 3, 3> V1;
|
|
|
+// V1.row(0) = V.row(F(fid1,0)) -o;
|
|
|
+// V1.row(1) = V.row(F(fid1,1)) -o;
|
|
|
+// V1.row(2) = V.row(F(fid1,2)) -o;
|
|
|
+// V1 = (P*V1.transpose()).transpose();
|
|
|
+//
|
|
|
+//// assert(V1(fid1_vc,2) < 10e-10);
|
|
|
+//// assert(V1((fid1_vc+1)%3,2) < 10e-10);
|
|
|
+//
|
|
|
+// // compute rotation R such that R * N1 = N0
|
|
|
+// // i.e. map both triangles to the same plane
|
|
|
+// double alpha = -atan2(V1((fid1_vc+2)%3,2),V1((fid1_vc+2)%3,1));
|
|
|
+//
|
|
|
+// Eigen::Matrix<typename DerivedV::Scalar, 3, 3> R;
|
|
|
+// R << 1, 0, 0,
|
|
|
+// 0, cos(alpha), -sin(alpha) ,
|
|
|
+// 0, sin(alpha), cos(alpha);
|
|
|
+// V1 = (R*V1.transpose()).transpose();
|
|
|
+//
|
|
|
+//// assert(V1(0,2) < 1e-10);
|
|
|
+//// assert(V1(1,2) < 1e-10);
|
|
|
+//// assert(V1(2,2) < 1e-10);
|
|
|
+//
|
|
|
+// // measure the angle between the reference frames
|
|
|
+// // k_ij is the angle between the triangle on the left and the one on the right
|
|
|
+// Eigen::Matrix<typename DerivedV::Scalar, 1, 3> ref0 = V0.row(1) - V0.row(0);
|
|
|
+// Eigen::Matrix<typename DerivedV::Scalar, 1, 3> ref1 = V1.row(1) - V1.row(0);
|
|
|
+//
|
|
|
+// ref0.normalize();
|
|
|
+// ref1.normalize();
|
|
|
+//
|
|
|
+// double ktemp = atan2(ref1(1),ref1(0)) - atan2(ref0(1),ref0(0));
|
|
|
+//
|
|
|
+// // just to be sure, rotate ref0 using angle ktemp...
|
|
|
+// Eigen::Matrix<typename DerivedV::Scalar, 2, 2> R2;
|
|
|
+// R2 << cos(ktemp), -sin(ktemp), sin(ktemp), cos(ktemp);
|
|
|
+//
|
|
|
+// Eigen::Matrix<typename DerivedV::Scalar, 1, 2> tmp1 = R2*(ref0.head(2)).transpose();
|
|
|
+//
|
|
|
+//// assert(tmp1(0) - ref1(0) < 1e-10);
|
|
|
+//// assert(tmp1(1) - ref1(1) < 1e-10);
|
|
|
+//
|
|
|
+// K[eid] = ktemp;
|
|
|
+// }
|
|
|
+// }
|
|
|
+//
|
|
|
+//}
|
|
|
|
|
|
|
|
|
IGL_INLINE void igl::n_polyvector_general(const Eigen::MatrixXd &V,
|
|
|
@@ -478,17 +475,19 @@ IGL_INLINE void igl::n_polyvector_general(const Eigen::MatrixXd &V,
|
|
|
const Eigen::VectorXi &I,
|
|
|
Eigen::MatrixXd &output)
|
|
|
{
|
|
|
- Eigen::VectorXi isConstrained = Eigen::VectorXi::Constant(F.rows(),0);
|
|
|
- Eigen::MatrixXd cfW = Eigen::MatrixXd::Constant(F.rows(),bc.cols(),0);
|
|
|
-
|
|
|
- for(unsigned i=0; i<b.size();++i)
|
|
|
- {
|
|
|
- isConstrained(b(i)) = 1;
|
|
|
- cfW.row(b(i)) << bc.row(i);
|
|
|
- }
|
|
|
- int n = I.rows();
|
|
|
- igl::GeneralPolyVectorFieldFinder<Eigen::MatrixXd, Eigen::MatrixXi> pvff(V,F,n);
|
|
|
- pvff.solve(isConstrained, cfW, I, output);
|
|
|
+ // This functions is broken, please contact Olga Diamanti
|
|
|
+ assert(0);
|
|
|
+ //Eigen::VectorXi isConstrained = Eigen::VectorXi::Constant(F.rows(),0);
|
|
|
+ //Eigen::MatrixXd cfW = Eigen::MatrixXd::Constant(F.rows(),bc.cols(),0);
|
|
|
+
|
|
|
+ //for(unsigned i=0; i<b.size();++i)
|
|
|
+ //{
|
|
|
+ // isConstrained(b(i)) = 1;
|
|
|
+ // cfW.row(b(i)) << bc.row(i);
|
|
|
+ //}
|
|
|
+ //int n = I.rows();
|
|
|
+ //igl::GeneralPolyVectorFieldFinder<Eigen::MatrixXd, Eigen::MatrixXi> pvff(V,F,n);
|
|
|
+ //pvff.solve(isConstrained, cfW, I, output);
|
|
|
}
|
|
|
|
|
|
|
Alec Jacobson