Merge pull request #1345 from libigl/alecjacobson-patch-2

fix https://github.com/libigl/libigl/issues/1343

include/igl/random_points_on_mesh.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-//
+// 
 // Copyright (C) 2014 Alec Jacobson <[email protected]>
-//
-// 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
+// 
+// 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 "random_points_on_mesh.h"
 #include "doublearea.h"
@@ -26,7 +26,6 @@ IGL_INLINE void igl::random_points_on_mesh(
   typedef Matrix<Scalar,Dynamic,1> VectorXs;
   VectorXs A;
   doublearea(V,F,A);
-  A /= A.array().sum();
   // Should be traingle mesh. Although Turk's method 1 generalizes...
   assert(F.cols() == 3);
   VectorXs C;
@@ -35,6 +34,10 @@ IGL_INLINE void igl::random_points_on_mesh(
   A0.bottomRightCorner(A.size(),1) = A;
   // Even faster would be to use the "Alias Table Method"
   cumsum(A0,1,C);
+  const Scalar Cmax = C(C.size()-1);
+  assert(Cmax > 0 && "Total surface area should be positive");
+  // Why is this more accurate than `C /= C(C.size()-1)` ?
+  for(int i = 0;i<C.size();i++) { C(i) = C(i)/Cmax; }
   const VectorXs R = (VectorXs::Random(n,1).array() + 1.)/2.;
   assert(R.minCoeff() >= 0);
   assert(R.maxCoeff() <= 1);
@@ -47,6 +50,31 @@ IGL_INLINE void igl::random_points_on_mesh(
   B.col(2) = S.array() * T.array().sqrt();
 }
 
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename DerivedB, 
+  typename DerivedFI,
+  typename DerivedX>
+IGL_INLINE void igl::random_points_on_mesh(
+  const int n,
+  const Eigen::MatrixBase<DerivedV > & V,
+  const Eigen::MatrixBase<DerivedF > & F,
+  Eigen::PlainObjectBase<DerivedB > & B,
+  Eigen::PlainObjectBase<DerivedFI > & FI,
+  Eigen::PlainObjectBase<DerivedX> & X)
+{
+  random_points_on_mesh(n,V,F,B,FI);
+  X = DerivedX::Zero(B.rows(),V.cols());
+  for(int x = 0;x<B.rows();x++)
+  {
+    for(int b = 0;b<B.cols();b++)
+    {
+      X.row(x) += B(x,b)*V.row(F(FI(x),b));
+    }
+  }
+}
+
 template <typename DerivedV, typename DerivedF, typename ScalarB, typename DerivedFI>
 IGL_INLINE void igl::random_points_on_mesh(
   const int n,
@@ -78,9 +106,8 @@ IGL_INLINE void igl::random_points_on_mesh(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-
+// generated by autoexplicit.sh
+template void igl::random_points_on_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -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> >(int, Eigen::MatrixBase<Eigen::Matrix<double, -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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::random_points_on_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::random_points_on_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::random_points_on_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 1, 0, 1, 1> >(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 1, 0, 1, 1> >&);
-
 #endif

include/igl/random_points_on_mesh.h

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-//
+// 
 // Copyright (C) 2013 Alec Jacobson <[email protected]>
-//
-// 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
+// 
+// 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/.
 #ifndef IGL_RANDOM_POINTS_ON_MESH_H
 #define IGL_RANDOM_POINTS_ON_MESH_H
@@ -23,7 +23,7 @@ namespace igl
   // Outputs:
   //   B  n by 3 list of barycentric coordinates, ith row are coordinates of
   //     ith sampled point in face FI(i)
-  //   FI  n list of indices into F
+  //   FI  n list of indices into F 
   //
   template <typename DerivedV, typename DerivedF, typename DerivedB, typename DerivedFI>
   IGL_INLINE void random_points_on_mesh(
@@ -33,6 +33,21 @@ namespace igl
     Eigen::PlainObjectBase<DerivedB > & B,
     Eigen::PlainObjectBase<DerivedFI > & FI);
   // Outputs:
+  //   X  n by dim list of sample positions.
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename DerivedB, 
+    typename DerivedFI,
+    typename DerivedX>
+  IGL_INLINE void random_points_on_mesh(
+    const int n,
+    const Eigen::MatrixBase<DerivedV > & V,
+    const Eigen::MatrixBase<DerivedF > & F,
+    Eigen::PlainObjectBase<DerivedB > & B,
+    Eigen::PlainObjectBase<DerivedFI > & FI,
+    Eigen::PlainObjectBase<DerivedX> & X);
+  // Outputs:
   //   B n by #V sparse matrix so that  B*V produces a list of sample points
   template <typename DerivedV, typename DerivedF, typename ScalarB, typename DerivedFI>
   IGL_INLINE void random_points_on_mesh(
@@ -48,5 +63,3 @@ namespace igl
 #endif
 
 #endif
-
-