reverted array to scaf

include/igl/scaf.cpp

@@ -381,7 +381,7 @@ void buildRhs(const Eigen::VectorXd &sqrt_M,
   }
 }
 
-void get_complement(const Eigen::VectorXi &bnd_ids, int v_n, Eigen::VectorXi &unknown_ids)
+void get_complement(const Eigen::VectorXi &bnd_ids, int v_n, Eigen::ArrayXi &unknown_ids)
 { // get the complement of bnd_ids.
   int assign = 0, i = 0;
   for (int get = 0; i < v_n && get < bnd_ids.size(); i++)
@@ -435,8 +435,8 @@ void build_surface_linear_system(const SCAFData &s, Eigen::SparseMatrix<double>
   {
     MatrixXd bnd_pos;
     igl::slice(s.w_uv, bnd_ids, 1, bnd_pos);
-    VectorXi known_ids(bnd_ids.size() * dim);
-    VectorXi unknown_ids((v_n - bnd_ids.rows()) * dim);
+    ArrayXi known_ids(bnd_ids.size() * dim);
+    ArrayXi unknown_ids((v_n - bnd_ids.rows()) * dim);
     get_complement(bnd_ids, v_n, unknown_ids);
     VectorXd known_pos(bnd_ids.size() * dim);
     for (int d = 0; d < dim; d++)
@@ -445,7 +445,7 @@ void build_surface_linear_system(const SCAFData &s, Eigen::SparseMatrix<double>
       known_ids.segment(d * n_b, n_b) = bnd_ids.array() + d * v_n;
       known_pos.segment(d * n_b, n_b) = bnd_pos.col(d);
       unknown_ids.block(d * (v_n - n_b), 0, v_n - n_b, unknown_ids.cols()) =
-          unknown_ids.topRows(v_n - n_b).array() + d * v_n;
+          unknown_ids.topRows(v_n - n_b) + d * v_n;
     }
 
     Eigen::SparseMatrix<double> Au, Ae;
@@ -497,8 +497,8 @@ void build_scaffold_linear_system(const SCAFData &s, Eigen::SparseMatrix<double>
   MatrixXd bnd_pos;
   igl::slice(s.w_uv, bnd_ids, 1, bnd_pos);
 
-  VectorXi known_ids(bnd_ids.size() * dim);
-  VectorXi unknown_ids((v_n - bnd_ids.rows()) * dim);
+  ArrayXi known_ids(bnd_ids.size() * dim);
+  ArrayXi unknown_ids((v_n - bnd_ids.rows()) * dim);
 
   get_complement(bnd_ids, v_n, unknown_ids);
 
@@ -509,7 +509,7 @@ void build_scaffold_linear_system(const SCAFData &s, Eigen::SparseMatrix<double>
     known_ids.segment(d * n_b, n_b) = bnd_ids.array() + d * v_n;
     known_pos.segment(d * n_b, n_b) = bnd_pos.col(d);
     unknown_ids.block(d * (v_n - n_b), 0, v_n - n_b, unknown_ids.cols()) =
-        unknown_ids.topRows(v_n - n_b).array() + d * v_n;
+        unknown_ids.topRows(v_n - n_b) + d * v_n;
   }
   Eigen::VectorXd sqrt_M = s.s_M.array().sqrt();
 
@@ -545,8 +545,8 @@ void solve_weighted_arap(SCAFData &s, Eigen::MatrixXd &uv)
   MatrixXd bnd_pos;
   igl::slice(s.w_uv, bnd_ids, 1, bnd_pos);
 
-  Eigen::VectorXi known_ids(bnd_n * dim);
-  Eigen::VectorXi unknown_ids((v_n - bnd_n) * dim);
+  ArrayXi known_ids(bnd_n * dim);
+  ArrayXi unknown_ids((v_n - bnd_n) * dim);
 
   get_complement(bnd_ids, v_n, unknown_ids);
 
@@ -557,7 +557,7 @@ void solve_weighted_arap(SCAFData &s, Eigen::MatrixXd &uv)
     known_ids.segment(d * n_b, n_b) = bnd_ids.array() + d * v_n;
     known_pos.segment(d * n_b, n_b) = bnd_pos.col(d);
     unknown_ids.block(d * (v_n - n_b), 0, v_n - n_b, unknown_ids.cols()) =
-        unknown_ids.topRows(v_n - n_b).array() + d * v_n;
+        unknown_ids.topRows(v_n - n_b) + d * v_n;
   }
 
   Eigen::SparseMatrix<double> L;
@@ -588,8 +588,8 @@ void solve_weighted_arap(SCAFData &s, Eigen::MatrixXd &uv)
 
   SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
   unknown_Uc = solver.compute(L).solve(rhs);
-  igl::slice_into(unknown_Uc, unknown_ids, 1, Uc);
-  igl::slice_into(known_pos, known_ids, 1, Uc);
+  igl::slice_into(unknown_Uc, unknown_ids.matrix(), 1, Uc);
+  igl::slice_into(known_pos, known_ids.matrix(), 1, Uc);
 
   uv = Map<Matrix<double, -1, -1, Eigen::ColMajor>>(Uc.data(), v_n, dim);
 }

include/igl/slice.cpp

@@ -139,6 +139,140 @@ IGL_INLINE void igl::slice(
 #endif
 }
 
+template <
+    typename TX,
+    typename TY,
+    typename DerivedR,
+    typename DerivedC>
+IGL_INLINE void igl::slice(
+    const Eigen::SparseMatrix<TX> &X,
+    const Eigen::ArrayBase<DerivedR> &R,
+    const Eigen::MatrixBase<DerivedC> &C,
+    Eigen::SparseMatrix<TY> &Y)
+{
+  int xm = X.rows();
+  int xn = X.cols();
+  int ym = R.size();
+  int yn = C.size();
+
+  // special case when R or C is empty
+  if (ym == 0 || yn == 0)
+  {
+    Y.resize(ym, yn);
+    return;
+  }
+
+  assert(R.minCoeff() >= 0);
+  assert(R.maxCoeff() < xm);
+  assert(C.minCoeff() >= 0);
+  assert(C.maxCoeff() < xn);
+
+  // Build reindexing maps for columns and rows, -1 means not in map
+  std::vector<std::vector<typename DerivedR::Scalar>> RI;
+  RI.resize(xm);
+  for (int i = 0; i < ym; i++)
+  {
+    RI[R(i)].push_back(i);
+  }
+  std::vector<std::vector<typename DerivedC::Scalar>> CI;
+  CI.resize(xn);
+  // initialize to -1
+  for (int i = 0; i < yn; i++)
+  {
+    CI[C(i)].push_back(i);
+  }
+  // Resize output
+  Eigen::DynamicSparseMatrix<TY, Eigen::RowMajor> dyn_Y(ym, yn);
+  // Take a guess at the number of nonzeros (this assumes uniform distribution
+  // not banded or heavily diagonal)
+  dyn_Y.reserve((X.nonZeros() / (X.rows() * X.cols())) * (ym * yn));
+  // Iterate over outside
+  for (int k = 0; k < X.outerSize(); ++k)
+  {
+    // Iterate over inside
+    for (typename Eigen::SparseMatrix<TX>::InnerIterator it(X, k); it; ++it)
+    {
+      typename std::vector<typename DerivedR::Scalar>::iterator rit;
+      typename std::vector<typename DerivedC::Scalar>::iterator cit;
+      for (rit = RI[it.row()].begin(); rit != RI[it.row()].end(); rit++)
+      {
+        for (cit = CI[it.col()].begin(); cit != CI[it.col()].end(); cit++)
+        {
+          dyn_Y.coeffRef(*rit, *cit) = it.value();
+        }
+      }
+    }
+  }
+  Y = Eigen::SparseMatrix<TY>(dyn_Y);
+}
+
+template <
+    typename TX,
+    typename TY,
+    typename DerivedR,
+    typename DerivedC>
+IGL_INLINE void igl::slice(
+    const Eigen::SparseMatrix<TX> &X,
+    const Eigen::MatrixBase<DerivedR> &R,
+    const Eigen::ArrayBase<DerivedC> &C,
+    Eigen::SparseMatrix<TY> &Y)
+{
+  int xm = X.rows();
+  int xn = X.cols();
+  int ym = R.size();
+  int yn = C.size();
+
+  // special case when R or C is empty
+  if (ym == 0 || yn == 0)
+  {
+    Y.resize(ym, yn);
+    return;
+  }
+
+  assert(R.minCoeff() >= 0);
+  assert(R.maxCoeff() < xm);
+  assert(C.minCoeff() >= 0);
+  assert(C.maxCoeff() < xn);
+
+  // Build reindexing maps for columns and rows, -1 means not in map
+  std::vector<std::vector<typename DerivedR::Scalar>> RI;
+  RI.resize(xm);
+  for (int i = 0; i < ym; i++)
+  {
+    RI[R(i)].push_back(i);
+  }
+  std::vector<std::vector<typename DerivedC::Scalar>> CI;
+  CI.resize(xn);
+  // initialize to -1
+  for (int i = 0; i < yn; i++)
+  {
+    CI[C(i)].push_back(i);
+  }
+  // Resize output
+  Eigen::DynamicSparseMatrix<TY, Eigen::RowMajor> dyn_Y(ym, yn);
+  // Take a guess at the number of nonzeros (this assumes uniform distribution
+  // not banded or heavily diagonal)
+  dyn_Y.reserve((X.nonZeros() / (X.rows() * X.cols())) * (ym * yn));
+  // Iterate over outside
+  for (int k = 0; k < X.outerSize(); ++k)
+  {
+    // Iterate over inside
+    for (typename Eigen::SparseMatrix<TX>::InnerIterator it(X, k); it; ++it)
+    {
+      typename std::vector<typename DerivedR::Scalar>::iterator rit;
+      typename std::vector<typename DerivedC::Scalar>::iterator cit;
+      for (rit = RI[it.row()].begin(); rit != RI[it.row()].end(); rit++)
+      {
+        for (cit = CI[it.col()].begin(); cit != CI[it.col()].end(); cit++)
+        {
+          dyn_Y.coeffRef(*rit, *cit) = it.value();
+        }
+      }
+    }
+  }
+  Y = Eigen::SparseMatrix<TY>(dyn_Y);
+}
+
 template <typename MatX, typename DerivedR, typename MatY>
 IGL_INLINE void igl::slice(
     const MatX &X,
@@ -173,6 +307,40 @@ IGL_INLINE void igl::slice(
   }
 }
 
+template <typename MatX, typename DerivedR, typename MatY>
+IGL_INLINE void igl::slice(
+    const MatX &X,
+    const Eigen::ArrayBase<DerivedR> &R,
+    const int dim,
+    MatY &Y)
+{
+  Eigen::Matrix<typename DerivedR::Scalar, Eigen::Dynamic, 1> C;
+  switch (dim)
+  {
+  case 1:
+    // boring base case
+    if (X.cols() == 0)
+    {
+      Y.resize(R.size(), 0);
+      return;
+    }
+    igl::colon(0, X.cols() - 1, C);
+    return slice(X, R, C, Y);
+  case 2:
+    // boring base case
+    if (X.rows() == 0)
+    {
+      Y.resize(0, R.size());
+      return;
+    }
+    igl::colon(0, X.rows() - 1, C);
+    return slice(X, C, R, Y);
+  default:
+    assert(false && "Unsupported dimension");
+    return;
+  }
+}
+
 template <
     typename DerivedX,
     typename DerivedR,
@@ -347,4 +515,5 @@ template void igl::slice<Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >
 template void igl::slice<Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>>>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const &, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const &, int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> &);
 template void igl::slice<Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::slice<Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>>>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const &, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const &, int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> &);
+template void igl::slice<Eigen::SparseMatrix<double, 0, int>, Eigen::Array<int, -1, 1, 0, -1, 1>, Eigen::SparseMatrix<double, 0, int> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::ArrayBase<Eigen::Array<int, -1, 1, 0, -1, 1> > const&, int, Eigen::SparseMatrix<double, 0, int>&);
 #endif

include/igl/slice.h

@@ -33,6 +33,26 @@ namespace igl
     const Eigen::MatrixBase<DerivedR> & R,
     const Eigen::MatrixBase<DerivedC> & C,
     Eigen::SparseMatrix<TY>& Y);
+  template <
+    typename TX,
+    typename TY,
+    typename DerivedR,
+    typename DerivedC>
+  IGL_INLINE void slice(
+    const Eigen::SparseMatrix<TX>& X,
+    const Eigen::ArrayBase<DerivedR> & R,
+    const Eigen::MatrixBase<DerivedC> & C,
+    Eigen::SparseMatrix<TY>& Y);
+  template <
+    typename TX,
+    typename TY,
+    typename DerivedR,
+    typename DerivedC>
+  IGL_INLINE void slice(
+    const Eigen::SparseMatrix<TX>& X,
+    const Eigen::MatrixBase<DerivedR> & R,
+    const Eigen::ArrayBase<DerivedC> & C,
+    Eigen::SparseMatrix<TY>& Y);
   // Wrapper to only slice in one direction
   //
   // Inputs:
@@ -48,6 +68,15 @@ namespace igl
     const Eigen::MatrixBase<DerivedR> & R,
     const int dim,
     MatY& Y);
+  template <
+    typename MatX,
+    typename DerivedR,
+    typename MatY>
+  IGL_INLINE void slice(
+    const MatX& X,
+    const Eigen::ArrayBase<DerivedR> & R,
+    const int dim,
+    MatY& Y);
   template <
     typename DerivedX,
     typename DerivedR,

include/igl/slice_into.cpp

@@ -137,4 +137,5 @@ template void igl::slice_into<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::M
 template void igl::slice_into<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::Matrix<double, -1, 1, 0, -1, 1>&);
 template void igl::slice_into<Eigen::SparseMatrix<double, 0, int>, Eigen::SparseMatrix<double, 0, int>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::SparseMatrix<double, 0, int>&);
 template void igl::slice_into<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::slice_into<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::MatrixWrapper<Eigen::Array<int, -1, 1, 0, -1, 1> > >(Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, Eigen::MatrixBase<Eigen::MatrixWrapper<Eigen::Array<int, -1, 1, 0, -1, 1> > > const&, int, Eigen::Matrix<double, -1, 1, 0, -1, 1>&);
 #endif