libigl_Martin/include/igl/n_polyvector_general.cpp

// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2014 Olga Diamanti <olga.diam@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 <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>
//#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,
                             const Eigen::MatrixXi &F,
                             const Eigen::VectorXi& b,
                             const Eigen::MatrixXd& bc,
                             const Eigen::VectorXi &I,
                             Eigen::MatrixXd &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);
}


#ifdef IGL_STATIC_LIBRARY
// Explicit template specialization
#endif