fwn for soups; clean up; templates; tutorial entry
This commit is contained in:
Alec Jacobson
File diff suppressed because it is too large
Load Diff
@@ -15,8 +15,7 @@ namespace igl
|
||||
// box
|
||||
//
|
||||
// Inputs:
|
||||
// V #V by 3 list of vertex positions
|
||||
// F #F by 3 list of triangle indices into V
|
||||
// V #V by 3 list of vertex/point positions
|
||||
// Returns length of bounding box diagonal
|
||||
IGL_INLINE double bounding_box_diagonal( const Eigen::MatrixXd & V);
|
||||
}
|
||||
|
||||
@@ -179,3 +179,8 @@ namespace igl {
|
||||
|
||||
|
||||
|
||||
#ifdef IGL_STATIC_LIBRARY
|
||||
// Explicit template instantiation
|
||||
// generated by autoexplicit.sh
|
||||
template void igl::copyleft::cgal::point_areas<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<double, -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<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
|
||||
#endif
|
||||
|
||||
@@ -38,6 +38,8 @@ namespace igl
|
||||
// N #P by 3 list of point normals
|
||||
// Outputs:
|
||||
// A #P list of estimated areas
|
||||
//
|
||||
// See also: igl::knn
|
||||
template <typename DerivedP, typename DerivedI, typename DerivedN,
|
||||
typename DerivedA>
|
||||
IGL_INLINE void point_areas(
|
||||
|
||||
+396
-302
@@ -4,321 +4,415 @@
|
||||
#include "parallel_for.h"
|
||||
#include "PI.h"
|
||||
#include <vector>
|
||||
#include <cassert>
|
||||
|
||||
namespace igl {
|
||||
template <typename DerivedP, typename DerivedA, typename DerivedN,
|
||||
typename Index, typename DerivedCH, typename DerivedCM, typename DerivedR,
|
||||
template <
|
||||
typename DerivedP,
|
||||
typename DerivedA,
|
||||
typename DerivedN,
|
||||
typename Index,
|
||||
typename DerivedCH,
|
||||
typename DerivedCM,
|
||||
typename DerivedR,
|
||||
typename DerivedEC>
|
||||
IGL_INLINE void fast_winding_number(const Eigen::MatrixBase<DerivedP>& P,
|
||||
const Eigen::MatrixBase<DerivedN>& N,
|
||||
const Eigen::MatrixBase<DerivedA>& A,
|
||||
const std::vector<std::vector<Index> > & point_indices,
|
||||
const Eigen::MatrixBase<DerivedCH>& CH,
|
||||
const int expansion_order,
|
||||
Eigen::PlainObjectBase<DerivedCM>& CM,
|
||||
Eigen::PlainObjectBase<DerivedR>& R,
|
||||
Eigen::PlainObjectBase<DerivedEC>& EC)
|
||||
IGL_INLINE void igl::fast_winding_number(
|
||||
const Eigen::MatrixBase<DerivedP>& P,
|
||||
const Eigen::MatrixBase<DerivedN>& N,
|
||||
const Eigen::MatrixBase<DerivedA>& A,
|
||||
const std::vector<std::vector<Index> > & point_indices,
|
||||
const Eigen::MatrixBase<DerivedCH>& CH,
|
||||
const int expansion_order,
|
||||
Eigen::PlainObjectBase<DerivedCM>& CM,
|
||||
Eigen::PlainObjectBase<DerivedR>& R,
|
||||
Eigen::PlainObjectBase<DerivedEC>& EC)
|
||||
{
|
||||
typedef typename DerivedP::Scalar real_p;
|
||||
typedef typename DerivedN::Scalar real_n;
|
||||
typedef typename DerivedA::Scalar real_a;
|
||||
typedef typename DerivedCM::Scalar real_cm;
|
||||
typedef typename DerivedR::Scalar real_r;
|
||||
typedef typename DerivedEC::Scalar real_ec;
|
||||
|
||||
typedef Eigen::Matrix<real_p,1,3> RowVec3p;
|
||||
|
||||
int m = CH.size();
|
||||
int num_terms;
|
||||
|
||||
assert(expansion_order < 3 && expansion_order >= 0 && "m must be less than n");
|
||||
if(expansion_order == 0){
|
||||
num_terms = 3;
|
||||
} else if(expansion_order ==1){
|
||||
num_terms = 3 + 9;
|
||||
} else if(expansion_order == 2){
|
||||
num_terms = 3 + 9 + 27;
|
||||
}
|
||||
|
||||
R.resize(m);
|
||||
CM.resize(m,3);
|
||||
EC.resize(m,num_terms);
|
||||
EC.setZero(m,num_terms);
|
||||
std::function< void(const int) > helper;
|
||||
helper = [&helper,
|
||||
&P,&N,&A,&expansion_order,&point_indices,&CH,&EC,&R,&CM]
|
||||
(const int index)-> void
|
||||
{
|
||||
typedef typename DerivedP::Scalar real_p;
|
||||
typedef typename DerivedN::Scalar real_n;
|
||||
typedef typename DerivedA::Scalar real_a;
|
||||
typedef typename DerivedCM::Scalar real_cm;
|
||||
typedef typename DerivedR::Scalar real_r;
|
||||
typedef typename DerivedEC::Scalar real_ec;
|
||||
|
||||
typedef Eigen::Matrix<real_p,1,3> RowVec3p;
|
||||
|
||||
int m = CH.size();
|
||||
int num_terms;
|
||||
|
||||
assert(expansion_order < 3 && expansion_order >= 0 && "m must be less than n");
|
||||
if(expansion_order == 0){
|
||||
num_terms = 3;
|
||||
} else if(expansion_order ==1){
|
||||
num_terms = 3 + 9;
|
||||
} else if(expansion_order == 2){
|
||||
num_terms = 3 + 9 + 27;
|
||||
}
|
||||
|
||||
R.resize(m);
|
||||
CM.resize(m,3);
|
||||
EC.resize(m,num_terms);
|
||||
EC.setZero(m,num_terms);
|
||||
std::function< void(const int) > helper;
|
||||
helper = [&helper,
|
||||
&P,&N,&A,&expansion_order,&point_indices,&CH,&EC,&R,&CM]
|
||||
(const int index)-> void
|
||||
{
|
||||
Eigen::Matrix<real_cm,1,3> masscenter;
|
||||
masscenter << 0,0,0;
|
||||
Eigen::Matrix<real_ec,1,3> zeroth_expansion;
|
||||
zeroth_expansion << 0,0,0;
|
||||
real_p areatotal = 0.0;
|
||||
for(int j = 0; j < point_indices.at(index).size(); j++){
|
||||
int curr_point_index = point_indices.at(index).at(j);
|
||||
|
||||
areatotal += A(curr_point_index);
|
||||
masscenter += A(curr_point_index)*P.row(curr_point_index);
|
||||
zeroth_expansion += A(curr_point_index)*N.row(curr_point_index);
|
||||
}
|
||||
|
||||
masscenter = masscenter/areatotal;
|
||||
CM.row(index) = masscenter;
|
||||
EC.block(index,0,1,3) = zeroth_expansion;
|
||||
|
||||
real_r max_norm = 0;
|
||||
real_r curr_norm;
|
||||
|
||||
for(int i = 0; i < point_indices.at(index).size(); i++){
|
||||
//Get max distance from center of mass:
|
||||
int curr_point_index = point_indices.at(index).at(i);
|
||||
Eigen::Matrix<real_r,1,3> point =
|
||||
P.row(curr_point_index)-masscenter;
|
||||
curr_norm = point.norm();
|
||||
if(curr_norm > max_norm){
|
||||
max_norm = curr_norm;
|
||||
}
|
||||
|
||||
//Calculate higher order terms if necessary
|
||||
Eigen::Matrix<real_ec,3,3> TempCoeffs;
|
||||
if(EC.cols() >= (3+9)){
|
||||
TempCoeffs = A(curr_point_index)*point.transpose()*
|
||||
N.row(curr_point_index);
|
||||
EC.block(index,3,1,9) +=
|
||||
Eigen::Map<Eigen::Matrix<real_ec,1,9> >(TempCoeffs.data(),
|
||||
TempCoeffs.size());
|
||||
}
|
||||
|
||||
if(EC.cols() == (3+9+27)){
|
||||
for(int k = 0; k < 3; k++){
|
||||
TempCoeffs = 0.5 * point(k) * (A(curr_point_index)*
|
||||
point.transpose()*N.row(curr_point_index));
|
||||
EC.block(index,12+9*k,1,9) += Eigen::Map<
|
||||
Eigen::Matrix<real_ec,1,9> >(TempCoeffs.data(),
|
||||
TempCoeffs.size());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
R(index) = max_norm;
|
||||
if(CH(index,0) != -1)
|
||||
{
|
||||
for(int i = 0; i < 8; i++){
|
||||
int child = CH(index,i);
|
||||
helper(child);
|
||||
}
|
||||
}
|
||||
};
|
||||
helper(0);
|
||||
}
|
||||
|
||||
template <typename DerivedP, typename DerivedA, typename DerivedN,
|
||||
typename Index, typename DerivedCH, typename DerivedCM, typename DerivedR,
|
||||
typename DerivedEC, typename DerivedQ, typename BetaType,
|
||||
typename DerivedWN>
|
||||
IGL_INLINE void fast_winding_number(const Eigen::MatrixBase<DerivedP>& P,
|
||||
const Eigen::MatrixBase<DerivedN>& N,
|
||||
const Eigen::MatrixBase<DerivedA>& A,
|
||||
const std::vector<std::vector<Index> > & point_indices,
|
||||
const Eigen::MatrixBase<DerivedCH>& CH,
|
||||
const Eigen::MatrixBase<DerivedCM>& CM,
|
||||
const Eigen::MatrixBase<DerivedR>& R,
|
||||
const Eigen::MatrixBase<DerivedEC>& EC,
|
||||
const Eigen::MatrixBase<DerivedQ>& Q,
|
||||
const BetaType beta,
|
||||
Eigen::PlainObjectBase<DerivedWN>& WN){
|
||||
|
||||
typedef typename DerivedP::Scalar real_p;
|
||||
typedef typename DerivedN::Scalar real_n;
|
||||
typedef typename DerivedA::Scalar real_a;
|
||||
typedef typename DerivedCM::Scalar real_cm;
|
||||
typedef typename DerivedR::Scalar real_r;
|
||||
typedef typename DerivedEC::Scalar real_ec;
|
||||
typedef typename DerivedQ::Scalar real_q;
|
||||
typedef typename DerivedWN::Scalar real_wn;
|
||||
|
||||
typedef Eigen::Matrix<real_q,1,3> RowVec;
|
||||
typedef Eigen::Matrix<real_ec,3,3> EC_3by3;
|
||||
|
||||
auto direct_eval = [](const RowVec & loc,
|
||||
const Eigen::Matrix<real_ec,1,3> & anorm){
|
||||
real_wn wn = (loc(0)*anorm(0)+loc(1)*anorm(1)+loc(2)*anorm(2))
|
||||
/(4.0*igl::PI*std::pow(loc.norm(),3));
|
||||
if(std::isnan(wn)){
|
||||
return 0.5;
|
||||
}else{
|
||||
return wn;
|
||||
}
|
||||
};
|
||||
|
||||
auto expansion_eval = [&direct_eval](const RowVec & loc,
|
||||
const Eigen::RowVectorXd & EC){
|
||||
real_wn wn = direct_eval(loc,EC.head<3>());
|
||||
double r = loc.norm();
|
||||
if(EC.size()>3){
|
||||
Eigen::Matrix<real_ec,3,3> SecondDerivative =
|
||||
Eigen::Matrix<real_ec,3,3>::Identity()/(4.0*igl::PI*std::pow(r,3));
|
||||
SecondDerivative += -3.0*loc.transpose()*loc/(4.0*igl::PI*std::pow(r,5));
|
||||
Eigen::Matrix<real_ec,1,9> derivative_vector =
|
||||
Eigen::Map<Eigen::Matrix<real_ec,1,9> >(SecondDerivative.data(),
|
||||
SecondDerivative.size());
|
||||
wn += derivative_vector.cwiseProduct(EC.segment<9>(3)).sum();
|
||||
Eigen::Matrix<real_cm,1,3> masscenter;
|
||||
masscenter << 0,0,0;
|
||||
Eigen::Matrix<real_ec,1,3> zeroth_expansion;
|
||||
zeroth_expansion << 0,0,0;
|
||||
real_p areatotal = 0.0;
|
||||
for(int j = 0; j < point_indices.at(index).size(); j++){
|
||||
int curr_point_index = point_indices.at(index).at(j);
|
||||
|
||||
areatotal += A(curr_point_index);
|
||||
masscenter += A(curr_point_index)*P.row(curr_point_index);
|
||||
zeroth_expansion += A(curr_point_index)*N.row(curr_point_index);
|
||||
}
|
||||
if(EC.size()>3+9){
|
||||
Eigen::Matrix<real_ec,3,3> ThirdDerivative;
|
||||
for(int i = 0; i < 3; i++){
|
||||
ThirdDerivative =
|
||||
15.0*loc(i)*loc.transpose()*loc/(4.0*igl::PI*std::pow(r,7));
|
||||
Eigen::Matrix<real_ec,3,3> Diagonal;
|
||||
Diagonal << loc(i), 0, 0,
|
||||
0, loc(i), 0,
|
||||
0, 0, loc(i);
|
||||
Eigen::Matrix<real_ec,3,3> RowCol;
|
||||
RowCol.setZero(3,3);
|
||||
RowCol.row(i) = loc;
|
||||
Eigen::Matrix<real_ec,3,3> RowColT = RowCol.transpose();
|
||||
RowCol = RowCol + RowColT;
|
||||
ThirdDerivative +=
|
||||
-3.0/(4.0*igl::PI*std::pow(r,5))*(RowCol+Diagonal);
|
||||
Eigen::Matrix<real_ec,1,9> derivative_vector =
|
||||
Eigen::Map<Eigen::Matrix<real_ec,1,9> >(ThirdDerivative.data(),
|
||||
ThirdDerivative.size());
|
||||
wn += derivative_vector.cwiseProduct(
|
||||
EC.segment<9>(12 + i*9)).sum();
|
||||
}
|
||||
}
|
||||
return wn;
|
||||
};
|
||||
|
||||
int m = Q.rows();
|
||||
WN.resize(m,1);
|
||||
|
||||
std::function< real_wn(const RowVec, const std::vector<int>) > helper;
|
||||
helper = [&helper,
|
||||
&P,&N,&A,
|
||||
&point_indices,&CH,
|
||||
&CM,&R,&EC,&beta,
|
||||
&direct_eval,&expansion_eval]
|
||||
(const RowVec query, const std::vector<int> near_indices)-> real_wn
|
||||
{
|
||||
std::vector<int> new_near_indices;
|
||||
real_wn wn = 0;
|
||||
for(int i = 0; i < near_indices.size(); i++){
|
||||
int index = near_indices.at(i);
|
||||
//Leaf Case, Brute force
|
||||
if(CH(index,0) == -1){
|
||||
for(int j = 0; j < point_indices.at(index).size(); j++){
|
||||
int curr_row = point_indices.at(index).at(j);
|
||||
wn += direct_eval(P.row(curr_row)-query,
|
||||
N.row(curr_row)*A(curr_row));
|
||||
}
|
||||
}
|
||||
//Non-Leaf Case
|
||||
else {
|
||||
for(int child = 0; child < 8; child++){
|
||||
int child_index = CH(index,child);
|
||||
if(point_indices.at(child_index).size() > 0){
|
||||
if((CM.row(child_index)-query).norm() > beta*R(child_index)){
|
||||
if(CH(child_index,0) == -1){
|
||||
for(int j=0;j<point_indices.at(child_index).size();j++){
|
||||
int curr_row = point_indices.at(child_index).at(j);
|
||||
wn += direct_eval(P.row(curr_row)-query,
|
||||
N.row(curr_row)*A(curr_row));
|
||||
}
|
||||
}else{
|
||||
wn += expansion_eval(CM.row(child_index)-query,
|
||||
EC.row(child_index));
|
||||
}
|
||||
}else {
|
||||
new_near_indices.emplace_back(child_index);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
if(new_near_indices.size() > 0){
|
||||
wn += helper(query,new_near_indices);
|
||||
}
|
||||
return wn;
|
||||
};
|
||||
|
||||
|
||||
if(beta > 0){
|
||||
std::vector<int> near_indices_start = {0};
|
||||
igl::parallel_for(m,[&](int iter){
|
||||
WN(iter) = helper(Q.row(iter),near_indices_start);
|
||||
},1000);
|
||||
} else {
|
||||
igl::parallel_for(m,[&](int iter){
|
||||
double wn = 0;
|
||||
for(int j = 0; j <P.rows(); j++){
|
||||
wn += direct_eval(P.row(j)-Q.row(iter),N.row(j)*A(j));
|
||||
}
|
||||
WN(iter) = wn;
|
||||
},1000);
|
||||
}
|
||||
}
|
||||
|
||||
template <typename DerivedP, typename DerivedA, typename DerivedN,
|
||||
typename DerivedQ, typename BetaType, typename DerivedWN>
|
||||
IGL_INLINE void fast_winding_number(const Eigen::MatrixBase<DerivedP>& P,
|
||||
const Eigen::MatrixBase<DerivedN>& N,
|
||||
const Eigen::MatrixBase<DerivedA>& A,
|
||||
const Eigen::MatrixBase<DerivedQ>& Q,
|
||||
const int expansion_order,
|
||||
const BetaType beta,
|
||||
Eigen::PlainObjectBase<DerivedWN>& WN
|
||||
){
|
||||
typedef typename DerivedWN::Scalar real;
|
||||
|
||||
std::vector<std::vector<int> > point_indices;
|
||||
Eigen::Matrix<int,Eigen::Dynamic,8> CH;
|
||||
Eigen::Matrix<real,Eigen::Dynamic,3> CN;
|
||||
Eigen::Matrix<real,Eigen::Dynamic,1> W;
|
||||
|
||||
octree(P,point_indices,CH,CN,W);
|
||||
|
||||
Eigen::Matrix<real,Eigen::Dynamic,Eigen::Dynamic> EC;
|
||||
Eigen::Matrix<real,Eigen::Dynamic,3> CM;
|
||||
Eigen::Matrix<real,Eigen::Dynamic,1> R;
|
||||
|
||||
fast_winding_number(P,N,A,point_indices,CH,expansion_order,CM,R,EC);
|
||||
fast_winding_number(P,N,A,point_indices,CH,CM,R,EC,Q,beta,WN);
|
||||
}
|
||||
|
||||
template <typename DerivedP, typename DerivedA, typename DerivedN,
|
||||
typename DerivedQ, typename DerivedWN>
|
||||
IGL_INLINE void fast_winding_number(const Eigen::MatrixBase<DerivedP>& P,
|
||||
const Eigen::MatrixBase<DerivedN>& N,
|
||||
const Eigen::MatrixBase<DerivedA>& A,
|
||||
const Eigen::MatrixBase<DerivedQ>& Q,
|
||||
Eigen::PlainObjectBase<DerivedWN>& WN
|
||||
){
|
||||
fast_winding_number(P,N,A,Q,2,2.0,WN);
|
||||
masscenter = masscenter/areatotal;
|
||||
CM.row(index) = masscenter;
|
||||
EC.block(index,0,1,3) = zeroth_expansion;
|
||||
|
||||
real_r max_norm = 0;
|
||||
real_r curr_norm;
|
||||
|
||||
for(int i = 0; i < point_indices.at(index).size(); i++){
|
||||
//Get max distance from center of mass:
|
||||
int curr_point_index = point_indices.at(index).at(i);
|
||||
Eigen::Matrix<real_r,1,3> point =
|
||||
P.row(curr_point_index)-masscenter;
|
||||
curr_norm = point.norm();
|
||||
if(curr_norm > max_norm){
|
||||
max_norm = curr_norm;
|
||||
}
|
||||
|
||||
//Calculate higher order terms if necessary
|
||||
Eigen::Matrix<real_ec,3,3> TempCoeffs;
|
||||
if(EC.cols() >= (3+9)){
|
||||
TempCoeffs = A(curr_point_index)*point.transpose()*
|
||||
N.row(curr_point_index);
|
||||
EC.block(index,3,1,9) +=
|
||||
Eigen::Map<Eigen::Matrix<real_ec,1,9> >(TempCoeffs.data(),
|
||||
TempCoeffs.size());
|
||||
}
|
||||
|
||||
if(EC.cols() == (3+9+27)){
|
||||
for(int k = 0; k < 3; k++){
|
||||
TempCoeffs = 0.5 * point(k) * (A(curr_point_index)*
|
||||
point.transpose()*N.row(curr_point_index));
|
||||
EC.block(index,12+9*k,1,9) += Eigen::Map<
|
||||
Eigen::Matrix<real_ec,1,9> >(TempCoeffs.data(),
|
||||
TempCoeffs.size());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
R(index) = max_norm;
|
||||
if(CH(index,0) != -1)
|
||||
{
|
||||
for(int i = 0; i < 8; i++){
|
||||
int child = CH(index,i);
|
||||
helper(child);
|
||||
}
|
||||
}
|
||||
};
|
||||
helper(0);
|
||||
}
|
||||
|
||||
template <
|
||||
typename DerivedP,
|
||||
typename DerivedA,
|
||||
typename DerivedN,
|
||||
typename Index,
|
||||
typename DerivedCH,
|
||||
typename DerivedCM,
|
||||
typename DerivedR,
|
||||
typename DerivedEC,
|
||||
typename DerivedQ,
|
||||
typename BetaType,
|
||||
typename DerivedWN>
|
||||
IGL_INLINE void igl::fast_winding_number(
|
||||
const Eigen::MatrixBase<DerivedP>& P,
|
||||
const Eigen::MatrixBase<DerivedN>& N,
|
||||
const Eigen::MatrixBase<DerivedA>& A,
|
||||
const std::vector<std::vector<Index> > & point_indices,
|
||||
const Eigen::MatrixBase<DerivedCH>& CH,
|
||||
const Eigen::MatrixBase<DerivedCM>& CM,
|
||||
const Eigen::MatrixBase<DerivedR>& R,
|
||||
const Eigen::MatrixBase<DerivedEC>& EC,
|
||||
const Eigen::MatrixBase<DerivedQ>& Q,
|
||||
const BetaType beta,
|
||||
Eigen::PlainObjectBase<DerivedWN>& WN)
|
||||
{
|
||||
|
||||
typedef typename DerivedP::Scalar real_p;
|
||||
typedef typename DerivedN::Scalar real_n;
|
||||
typedef typename DerivedA::Scalar real_a;
|
||||
typedef typename DerivedCM::Scalar real_cm;
|
||||
typedef typename DerivedR::Scalar real_r;
|
||||
typedef typename DerivedEC::Scalar real_ec;
|
||||
typedef typename DerivedQ::Scalar real_q;
|
||||
typedef typename DerivedWN::Scalar real_wn;
|
||||
|
||||
typedef Eigen::Matrix<real_q,1,3> RowVec;
|
||||
typedef Eigen::Matrix<real_ec,3,3> EC_3by3;
|
||||
|
||||
auto direct_eval = [](const RowVec & loc,
|
||||
const Eigen::Matrix<real_ec,1,3> & anorm)
|
||||
{
|
||||
real_wn wn = (loc(0)*anorm(0)+loc(1)*anorm(1)+loc(2)*anorm(2))
|
||||
/(4.0*igl::PI*std::pow(loc.norm(),3));
|
||||
if(std::isnan(wn)){
|
||||
return 0.5;
|
||||
}else{
|
||||
return wn;
|
||||
}
|
||||
};
|
||||
|
||||
auto expansion_eval = [&direct_eval](const RowVec & loc,
|
||||
const Eigen::RowVectorXd & EC){
|
||||
real_wn wn = direct_eval(loc,EC.head<3>());
|
||||
double r = loc.norm();
|
||||
if(EC.size()>3){
|
||||
Eigen::Matrix<real_ec,3,3> SecondDerivative =
|
||||
Eigen::Matrix<real_ec,3,3>::Identity()/(4.0*igl::PI*std::pow(r,3));
|
||||
SecondDerivative += -3.0*loc.transpose()*loc/(4.0*igl::PI*std::pow(r,5));
|
||||
Eigen::Matrix<real_ec,1,9> derivative_vector =
|
||||
Eigen::Map<Eigen::Matrix<real_ec,1,9> >(SecondDerivative.data(),
|
||||
SecondDerivative.size());
|
||||
wn += derivative_vector.cwiseProduct(EC.segment<9>(3)).sum();
|
||||
}
|
||||
if(EC.size()>3+9){
|
||||
Eigen::Matrix<real_ec,3,3> ThirdDerivative;
|
||||
for(int i = 0; i < 3; i++){
|
||||
ThirdDerivative =
|
||||
15.0*loc(i)*loc.transpose()*loc/(4.0*igl::PI*std::pow(r,7));
|
||||
Eigen::Matrix<real_ec,3,3> Diagonal;
|
||||
Diagonal << loc(i), 0, 0,
|
||||
0, loc(i), 0,
|
||||
0, 0, loc(i);
|
||||
Eigen::Matrix<real_ec,3,3> RowCol;
|
||||
RowCol.setZero(3,3);
|
||||
RowCol.row(i) = loc;
|
||||
Eigen::Matrix<real_ec,3,3> RowColT = RowCol.transpose();
|
||||
RowCol = RowCol + RowColT;
|
||||
ThirdDerivative +=
|
||||
-3.0/(4.0*igl::PI*std::pow(r,5))*(RowCol+Diagonal);
|
||||
Eigen::Matrix<real_ec,1,9> derivative_vector =
|
||||
Eigen::Map<Eigen::Matrix<real_ec,1,9> >(ThirdDerivative.data(),
|
||||
ThirdDerivative.size());
|
||||
wn += derivative_vector.cwiseProduct(
|
||||
EC.segment<9>(12 + i*9)).sum();
|
||||
}
|
||||
}
|
||||
return wn;
|
||||
};
|
||||
|
||||
int m = Q.rows();
|
||||
WN.resize(m,1);
|
||||
|
||||
std::function< real_wn(const RowVec, const std::vector<int>) > helper;
|
||||
helper = [&helper,
|
||||
&P,&N,&A,
|
||||
&point_indices,&CH,
|
||||
&CM,&R,&EC,&beta,
|
||||
&direct_eval,&expansion_eval]
|
||||
(const RowVec query, const std::vector<int> near_indices)-> real_wn
|
||||
{
|
||||
std::vector<int> new_near_indices;
|
||||
real_wn wn = 0;
|
||||
for(int i = 0; i < near_indices.size(); i++){
|
||||
int index = near_indices.at(i);
|
||||
//Leaf Case, Brute force
|
||||
if(CH(index,0) == -1){
|
||||
for(int j = 0; j < point_indices.at(index).size(); j++){
|
||||
int curr_row = point_indices.at(index).at(j);
|
||||
wn += direct_eval(P.row(curr_row)-query,
|
||||
N.row(curr_row)*A(curr_row));
|
||||
}
|
||||
}
|
||||
//Non-Leaf Case
|
||||
else {
|
||||
for(int child = 0; child < 8; child++){
|
||||
int child_index = CH(index,child);
|
||||
if(point_indices.at(child_index).size() > 0){
|
||||
if((CM.row(child_index)-query).norm() > beta*R(child_index)){
|
||||
if(CH(child_index,0) == -1){
|
||||
for(int j=0;j<point_indices.at(child_index).size();j++){
|
||||
int curr_row = point_indices.at(child_index).at(j);
|
||||
wn += direct_eval(P.row(curr_row)-query,
|
||||
N.row(curr_row)*A(curr_row));
|
||||
}
|
||||
}else{
|
||||
wn += expansion_eval(CM.row(child_index)-query,
|
||||
EC.row(child_index));
|
||||
}
|
||||
}else {
|
||||
new_near_indices.emplace_back(child_index);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
if(new_near_indices.size() > 0){
|
||||
wn += helper(query,new_near_indices);
|
||||
}
|
||||
return wn;
|
||||
};
|
||||
|
||||
|
||||
if(beta > 0){
|
||||
std::vector<int> near_indices_start = {0};
|
||||
igl::parallel_for(m,[&](int iter){
|
||||
WN(iter) = helper(Q.row(iter),near_indices_start);
|
||||
},1000);
|
||||
} else {
|
||||
igl::parallel_for(m,[&](int iter){
|
||||
double wn = 0;
|
||||
for(int j = 0; j <P.rows(); j++){
|
||||
wn += direct_eval(P.row(j)-Q.row(iter),N.row(j)*A(j));
|
||||
}
|
||||
WN(iter) = wn;
|
||||
},1000);
|
||||
}
|
||||
}
|
||||
|
||||
template <
|
||||
typename DerivedP,
|
||||
typename DerivedA,
|
||||
typename DerivedN,
|
||||
typename DerivedQ,
|
||||
typename BetaType,
|
||||
typename DerivedWN>
|
||||
IGL_INLINE void igl::fast_winding_number(
|
||||
const Eigen::MatrixBase<DerivedP>& P,
|
||||
const Eigen::MatrixBase<DerivedN>& N,
|
||||
const Eigen::MatrixBase<DerivedA>& A,
|
||||
const Eigen::MatrixBase<DerivedQ>& Q,
|
||||
const int expansion_order,
|
||||
const BetaType beta,
|
||||
Eigen::PlainObjectBase<DerivedWN>& WN)
|
||||
{
|
||||
typedef typename DerivedWN::Scalar real;
|
||||
|
||||
std::vector<std::vector<int> > point_indices;
|
||||
Eigen::Matrix<int,Eigen::Dynamic,8> CH;
|
||||
Eigen::Matrix<real,Eigen::Dynamic,3> CN;
|
||||
Eigen::Matrix<real,Eigen::Dynamic,1> W;
|
||||
|
||||
octree(P,point_indices,CH,CN,W);
|
||||
|
||||
Eigen::Matrix<real,Eigen::Dynamic,Eigen::Dynamic> EC;
|
||||
Eigen::Matrix<real,Eigen::Dynamic,3> CM;
|
||||
Eigen::Matrix<real,Eigen::Dynamic,1> R;
|
||||
|
||||
fast_winding_number(P,N,A,point_indices,CH,expansion_order,CM,R,EC);
|
||||
fast_winding_number(P,N,A,point_indices,CH,CM,R,EC,Q,beta,WN);
|
||||
}
|
||||
|
||||
template <
|
||||
typename DerivedP,
|
||||
typename DerivedA,
|
||||
typename DerivedN,
|
||||
typename DerivedQ,
|
||||
typename DerivedWN>
|
||||
IGL_INLINE void igl::fast_winding_number(
|
||||
const Eigen::MatrixBase<DerivedP>& P,
|
||||
const Eigen::MatrixBase<DerivedN>& N,
|
||||
const Eigen::MatrixBase<DerivedA>& A,
|
||||
const Eigen::MatrixBase<DerivedQ>& Q,
|
||||
Eigen::PlainObjectBase<DerivedWN>& WN)
|
||||
{
|
||||
fast_winding_number(P,N,A,Q,2,2.0,WN);
|
||||
}
|
||||
|
||||
template <
|
||||
typename DerivedV,
|
||||
typename DerivedF,
|
||||
typename DerivedQ,
|
||||
typename DerivedW>
|
||||
IGL_INLINE void igl::fast_winding_number(
|
||||
const Eigen::MatrixBase<DerivedV> & V,
|
||||
const Eigen::MatrixBase<DerivedF> & F,
|
||||
const Eigen::MatrixBase<DerivedQ> & Q,
|
||||
Eigen::PlainObjectBase<DerivedW> & W)
|
||||
{
|
||||
igl::FastWindingNumberBVH fwn_bvh;
|
||||
int order = 2;
|
||||
igl::fast_winding_number(V,F,order,fwn_bvh);
|
||||
float accuracy_scale = 2;
|
||||
igl::fast_winding_number(fwn_bvh,accuracy_scale,Q,W);
|
||||
}
|
||||
|
||||
template <
|
||||
typename DerivedV,
|
||||
typename DerivedF>
|
||||
IGL_INLINE void igl::fast_winding_number(
|
||||
const Eigen::MatrixBase<DerivedV> & V,
|
||||
const Eigen::MatrixBase<DerivedF> & F,
|
||||
const int order,
|
||||
FastWindingNumberBVH & fwn_bvh)
|
||||
{
|
||||
assert(V.cols() == 3 && "V should be 3D");
|
||||
assert(F.cols() == 3 && "F should contain triangles");
|
||||
// Extra copies. Usuually this won't be the bottleneck.
|
||||
fwn_bvh.U.resize(V.rows());
|
||||
for(int i = 0;i<V.rows();i++)
|
||||
{
|
||||
for(int j = 0;j<3;j++)
|
||||
{
|
||||
fwn_bvh.U[i][j] = V(i,j);
|
||||
}
|
||||
}
|
||||
// Wouldn't need to copy if F is **RowMajor**
|
||||
fwn_bvh.F.resize(F.size());
|
||||
for(int f = 0;f<F.rows();f++)
|
||||
{
|
||||
for(int c = 0;c<F.cols();c++)
|
||||
{
|
||||
fwn_bvh.F[c+f*F.cols()] = F(f,c);
|
||||
}
|
||||
}
|
||||
fwn_bvh.ut_solid_angle.clear();
|
||||
fwn_bvh.ut_solid_angle.init(
|
||||
fwn_bvh.F.size()/3,
|
||||
&fwn_bvh.F[0],
|
||||
fwn_bvh.U.size(),
|
||||
&fwn_bvh.U[0],
|
||||
order);
|
||||
}
|
||||
|
||||
template <
|
||||
typename DerivedQ,
|
||||
typename DerivedW>
|
||||
IGL_INLINE void igl::fast_winding_number(
|
||||
const FastWindingNumberBVH & fwn_bvh,
|
||||
const float _accuracy_scale,
|
||||
const Eigen::MatrixBase<DerivedQ> & Q,
|
||||
Eigen::PlainObjectBase<DerivedW> & W)
|
||||
{
|
||||
const double accuracy_scale = 2;
|
||||
assert(Q.cols() == 3 && "Q should be 3D");
|
||||
W.resize(Q.rows(),1);
|
||||
igl::parallel_for(Q.rows(),[&](int p)
|
||||
{
|
||||
FastWindingNumber::HDK_Sample::UT_Vector3T<float>Qp;
|
||||
Qp[0] = Q(p,0);
|
||||
Qp[1] = Q(p,1);
|
||||
Qp[2] = Q(p,2);
|
||||
W(p) = fwn_bvh.ut_solid_angle.computeSolidAngle(
|
||||
Qp,
|
||||
accuracy_scale)
|
||||
/ (4.0*M_PI);
|
||||
},1000);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
#ifdef IGL_STATIC_LIBRARY
|
||||
// Explicit template instantiation
|
||||
// generated by autoexplicit.sh
|
||||
template void igl::fast_winding_number<Eigen::CwiseUnaryOp<Eigen::internal::scalar_cast_op<double, float>, Eigen::Matrix<double, -1, -1, 0, -1, -1> const>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_cast_op<double, float>, Eigen::Matrix<double, -1, -1, 0, -1, -1> const> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, igl::FastWindingNumberBVH&);
|
||||
// generated by autoexplicit.sh
|
||||
template void igl::fast_winding_number<Eigen::CwiseUnaryOp<Eigen::internal::scalar_cast_op<double, float>, Eigen::Matrix<double, -1, -1, 0, -1, -1> const>, Eigen::Matrix<float, -1, 1, 0, -1, 1> >(igl::FastWindingNumberBVH const&, float, Eigen::MatrixBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_cast_op<double, float>, Eigen::Matrix<double, -1, -1, 0, -1, -1> const> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> >&);
|
||||
template void igl::fast_winding_number<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, int, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, int, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
|
||||
template void igl::fast_winding_number<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, int, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
|
||||
#endif
|
||||
|
||||
@@ -1,6 +1,7 @@
|
||||
#ifndef IGL_FAST_WINDING_NUMBER
|
||||
#define IGL_FAST_WINDING_NUMBER
|
||||
#include "igl_inline.h"
|
||||
#include "FastWindingNumberForSoups.h"
|
||||
#include <Eigen/Core>
|
||||
#include <vector>
|
||||
namespace igl
|
||||
@@ -12,8 +13,8 @@ namespace igl
|
||||
// data, and an expansion order, we define a taylor series expansion at each
|
||||
// octree cell.
|
||||
//
|
||||
// The octree data is designed to come from igl::octree, and the areas
|
||||
// (if not obtained at scan time), may be calculated using
|
||||
// The octree data is designed to come from igl::octree, and the areas (if not
|
||||
// obtained at scan time), may be calculated using
|
||||
// igl::copyleft::cgal::point_areas.
|
||||
//
|
||||
// Inputs:
|
||||
@@ -32,19 +33,26 @@ namespace igl
|
||||
// EC #OctreeCells by #TaylorCoefficients list of expansion coefficients.
|
||||
// (Note that #TaylorCoefficients = ∑_{i=1}^{expansion_order} 3^i)
|
||||
//
|
||||
template <typename DerivedP, typename DerivedA, typename DerivedN,
|
||||
typename Index, typename DerivedCH, typename DerivedCM, typename DerivedR,
|
||||
// See also: igl::copyleft::cgal::point_areas, igl::knn
|
||||
template <
|
||||
typename DerivedP,
|
||||
typename DerivedA,
|
||||
typename DerivedN,
|
||||
typename Index,
|
||||
typename DerivedCH,
|
||||
typename DerivedCM,
|
||||
typename DerivedR,
|
||||
typename DerivedEC>
|
||||
IGL_INLINE void fast_winding_number(const Eigen::MatrixBase<DerivedP>& P,
|
||||
const Eigen::MatrixBase<DerivedN>& N,
|
||||
const Eigen::MatrixBase<DerivedA>& A,
|
||||
const std::vector<std::vector<Index> > & point_indices,
|
||||
const Eigen::MatrixBase<DerivedCH>& CH,
|
||||
const int expansion_order,
|
||||
Eigen::PlainObjectBase<DerivedCM>& CM,
|
||||
Eigen::PlainObjectBase<DerivedR>& R,
|
||||
Eigen::PlainObjectBase<DerivedEC>& EC);
|
||||
|
||||
IGL_INLINE void fast_winding_number(
|
||||
const Eigen::MatrixBase<DerivedP>& P,
|
||||
const Eigen::MatrixBase<DerivedN>& N,
|
||||
const Eigen::MatrixBase<DerivedA>& A,
|
||||
const std::vector<std::vector<Index> > & point_indices,
|
||||
const Eigen::MatrixBase<DerivedCH>& CH,
|
||||
const int expansion_order,
|
||||
Eigen::PlainObjectBase<DerivedCM>& CM,
|
||||
Eigen::PlainObjectBase<DerivedR>& R,
|
||||
Eigen::PlainObjectBase<DerivedEC>& EC);
|
||||
// Evaluate the fast winding number for point data, having already done the
|
||||
// the precomputation
|
||||
//
|
||||
@@ -70,22 +78,45 @@ namespace igl
|
||||
// Outputs:
|
||||
// WN #Q by 1 list of windinng number values at each query point
|
||||
//
|
||||
template <typename DerivedP, typename DerivedA, typename DerivedN,
|
||||
typename Index, typename DerivedCH, typename DerivedCM, typename DerivedR,
|
||||
typename DerivedEC, typename DerivedQ, typename BetaType,
|
||||
template <
|
||||
typename DerivedP,
|
||||
typename DerivedA,
|
||||
typename DerivedN,
|
||||
typename Index,
|
||||
typename DerivedCH,
|
||||
typename DerivedCM,
|
||||
typename DerivedR,
|
||||
typename DerivedEC,
|
||||
typename DerivedQ,
|
||||
typename BetaType,
|
||||
typename DerivedWN>
|
||||
IGL_INLINE void fast_winding_number(const Eigen::MatrixBase<DerivedP>& P,
|
||||
const Eigen::MatrixBase<DerivedN>& N,
|
||||
const Eigen::MatrixBase<DerivedA>& A,
|
||||
const std::vector<std::vector<Index> > & point_indices,
|
||||
const Eigen::MatrixBase<DerivedCH>& CH,
|
||||
const Eigen::MatrixBase<DerivedCM>& CM,
|
||||
const Eigen::MatrixBase<DerivedR>& R,
|
||||
const Eigen::MatrixBase<DerivedEC>& EC,
|
||||
const Eigen::MatrixBase<DerivedQ>& Q,
|
||||
const BetaType beta,
|
||||
Eigen::PlainObjectBase<DerivedWN>& WN);
|
||||
|
||||
IGL_INLINE void fast_winding_number(
|
||||
const Eigen::MatrixBase<DerivedP>& P,
|
||||
const Eigen::MatrixBase<DerivedN>& N,
|
||||
const Eigen::MatrixBase<DerivedA>& A,
|
||||
const std::vector<std::vector<Index> > & point_indices,
|
||||
const Eigen::MatrixBase<DerivedCH>& CH,
|
||||
const Eigen::MatrixBase<DerivedCM>& CM,
|
||||
const Eigen::MatrixBase<DerivedR>& R,
|
||||
const Eigen::MatrixBase<DerivedEC>& EC,
|
||||
const Eigen::MatrixBase<DerivedQ>& Q,
|
||||
const BetaType beta,
|
||||
Eigen::PlainObjectBase<DerivedWN>& WN);
|
||||
template <
|
||||
typename DerivedP,
|
||||
typename DerivedA,
|
||||
typename DerivedN,
|
||||
typename DerivedQ,
|
||||
typename BetaType,
|
||||
typename DerivedWN>
|
||||
IGL_INLINE void fast_winding_number(
|
||||
const Eigen::MatrixBase<DerivedP>& P,
|
||||
const Eigen::MatrixBase<DerivedN>& N,
|
||||
const Eigen::MatrixBase<DerivedA>& A,
|
||||
const Eigen::MatrixBase<DerivedQ>& Q,
|
||||
const int expansion_order,
|
||||
const BetaType beta,
|
||||
Eigen::PlainObjectBase<DerivedWN>& WN);
|
||||
// Evaluate the fast winding number for point data, with default expansion
|
||||
// order and beta (both are set to 2).
|
||||
//
|
||||
@@ -101,14 +132,80 @@ namespace igl
|
||||
// Outputs:
|
||||
// WN #Q by 1 list of windinng number values at each query point
|
||||
//
|
||||
template <typename DerivedP, typename DerivedA, typename DerivedN,
|
||||
typename DerivedQ, typename DerivedWN>
|
||||
IGL_INLINE void fast_winding_number(const Eigen::MatrixBase<DerivedP>& P,
|
||||
const Eigen::MatrixBase<DerivedN>& N,
|
||||
const Eigen::MatrixBase<DerivedA>& A,
|
||||
const Eigen::MatrixBase<DerivedQ>& Q,
|
||||
Eigen::PlainObjectBase<DerivedWN>& WN
|
||||
);
|
||||
template <
|
||||
typename DerivedP,
|
||||
typename DerivedA,
|
||||
typename DerivedN,
|
||||
typename DerivedQ,
|
||||
typename DerivedWN>
|
||||
IGL_INLINE void fast_winding_number(
|
||||
const Eigen::MatrixBase<DerivedP>& P,
|
||||
const Eigen::MatrixBase<DerivedN>& N,
|
||||
const Eigen::MatrixBase<DerivedA>& A,
|
||||
const Eigen::MatrixBase<DerivedQ>& Q,
|
||||
Eigen::PlainObjectBase<DerivedWN>& WN);
|
||||
// Class declaration
|
||||
namespace FastWindingNumber { namespace HDK_Sample{ template <typename T1, typename T2> class UT_SolidAngle;} }
|
||||
struct FastWindingNumberBVH {
|
||||
FastWindingNumber::HDK_Sample::UT_SolidAngle<float,float> ut_solid_angle;
|
||||
// Need copies of these so they stay alive between calls.
|
||||
std::vector<FastWindingNumber::HDK_Sample::UT_Vector3T<float> > U;
|
||||
std::vector<int> F;
|
||||
};
|
||||
// Compute approximate winding number of a triangle soup mesh according to
|
||||
// "Fast Winding Numbers for Soups and Clouds" [Barill et al. 2018].
|
||||
//
|
||||
// Inputs:
|
||||
// V #V by 3 list of mesh vertex positions
|
||||
// F #F by 3 list of triangle mesh indices into rows of V
|
||||
// Q #Q by 3 list of query positions
|
||||
// Outputs:
|
||||
// W #Q list of winding number values
|
||||
template <
|
||||
typename DerivedV,
|
||||
typename DerivedF,
|
||||
typename DerivedQ,
|
||||
typename DerivedW>
|
||||
IGL_INLINE void fast_winding_number(
|
||||
const Eigen::MatrixBase<DerivedV> & V,
|
||||
const Eigen::MatrixBase<DerivedF> & F,
|
||||
const Eigen::MatrixBase<DerivedQ> & Q,
|
||||
Eigen::PlainObjectBase<DerivedW> & W);
|
||||
// Precomputation for computing approximate winding numbers of a triangle
|
||||
// soup.
|
||||
//
|
||||
// Inputs:
|
||||
// V #V by 3 list of mesh vertex positions
|
||||
// F #F by 3 list of triangle mesh indices into rows of V
|
||||
// order Taylor series expansion order to use (e.g., 2)
|
||||
// Outputs:
|
||||
// fwn_bvh Precomputed bounding volume hierarchy
|
||||
//
|
||||
template <
|
||||
typename DerivedV,
|
||||
typename DerivedF>
|
||||
IGL_INLINE void fast_winding_number(
|
||||
const Eigen::MatrixBase<DerivedV> & V,
|
||||
const Eigen::MatrixBase<DerivedF> & F,
|
||||
const int order,
|
||||
FastWindingNumberBVH & fwn_bvh);
|
||||
// After precomputation, compute winding number at a each of many points in a
|
||||
// list.
|
||||
//
|
||||
// Inputs:
|
||||
// fwn_bvh Precomputed bounding volume hierarchy
|
||||
// accuracy_scale parameter controlling accuracy (e.g., 2)
|
||||
// Q #Q by 3 list of query positions
|
||||
// Outputs:
|
||||
// W #Q list of winding number values
|
||||
template <
|
||||
typename DerivedQ,
|
||||
typename DerivedW>
|
||||
IGL_INLINE void fast_winding_number(
|
||||
const FastWindingNumberBVH & fwn_bvh,
|
||||
const float accuracy_scale,
|
||||
const Eigen::MatrixBase<DerivedQ> & Q,
|
||||
Eigen::PlainObjectBase<DerivedW> & W);
|
||||
}
|
||||
#ifndef IGL_STATIC_LIBRARY
|
||||
# include "fast_winding_number.cpp"
|
||||
|
||||
+3
-1
@@ -103,6 +103,8 @@ namespace igl {
|
||||
|
||||
|
||||
#ifdef IGL_STATIC_LIBRARY
|
||||
// Explicit template instantiation
|
||||
// generated by autoexplicit.sh
|
||||
template void igl::knn<Eigen::Matrix<double, -1, -1, 0, -1, -1>, int, int, Eigen::Matrix<int, -1, -1, 0, -1, -1>, 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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
|
||||
template void igl::knn<Eigen::Matrix<double, -1, -1, 0, -1, -1>, int, int, Eigen::Matrix<int, -1, 8, 0, -1, 8>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 8, 0, -1, 8> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
|
||||
#endif
|
||||
|
||||
|
||||
@@ -176,5 +176,7 @@ namespace igl {
|
||||
|
||||
#ifdef IGL_STATIC_LIBRARY
|
||||
// Explicit template instantiation
|
||||
// generated by autoexplicit.sh
|
||||
template void igl::octree<Eigen::Matrix<double, -1, -1, 0, -1, -1>, int, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
|
||||
template void igl::octree<Eigen::Matrix<double, -1, -1, 0, -1, -1>, int, Eigen::Matrix<int, -1, 8, 0, -1, 8>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 8, 0, -1, 8> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
|
||||
#endif
|
||||
|
||||
@@ -247,6 +247,8 @@ IGL_INLINE DerivedX igl::slice(
|
||||
|
||||
#ifdef IGL_STATIC_LIBRARY
|
||||
// Explicit template instantiation
|
||||
// generated by autoexplicit.sh
|
||||
template void igl::slice<Eigen::MatrixBase<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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
|
||||
#if EIGEN_VERSION_AT_LEAST(3,3,0)
|
||||
#else
|
||||
// generated by autoexplicit.sh
|
||||
|
||||
@@ -0,0 +1,5 @@
|
||||
get_filename_component(PROJECT_NAME ${CMAKE_CURRENT_SOURCE_DIR} NAME)
|
||||
project(${PROJECT_NAME})
|
||||
|
||||
add_executable(${PROJECT_NAME}_bin main.cpp)
|
||||
target_link_libraries(${PROJECT_NAME}_bin igl::core igl::opengl igl::opengl_glfw igl::cgal tutorials)
|
||||
Executable
+200
@@ -0,0 +1,200 @@
|
||||
#include "tutorial_shared_path.h"
|
||||
#include <igl/fast_winding_number.h>
|
||||
#include <igl/read_triangle_mesh.h>
|
||||
#include <igl/slice_mask.h>
|
||||
#include <Eigen/Geometry>
|
||||
#include <igl/octree.h>
|
||||
#include <igl/barycenter.h>
|
||||
#include <igl/knn.h>
|
||||
#include <igl/random_points_on_mesh.h>
|
||||
#include <igl/bounding_box_diagonal.h>
|
||||
#include <igl/per_face_normals.h>
|
||||
#include <igl/copyleft/cgal/point_areas.h>
|
||||
#include <igl/opengl/glfw/Viewer.h>
|
||||
#include <igl/get_seconds.h>
|
||||
#include <iostream>
|
||||
#include <cstdlib>
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
const auto time = [](std::function<void(void)> func)->double
|
||||
{
|
||||
const double t_before = igl::get_seconds();
|
||||
func();
|
||||
const double t_after = igl::get_seconds();
|
||||
return t_after-t_before;
|
||||
};
|
||||
|
||||
Eigen::MatrixXd V;
|
||||
Eigen::MatrixXi F;
|
||||
igl::read_triangle_mesh(argc>1?argv[1]:TUTORIAL_SHARED_PATH "/bunny.off",V,F);
|
||||
// Sample mesh for point cloud
|
||||
Eigen::MatrixXd P,N;
|
||||
{
|
||||
Eigen::VectorXi I;
|
||||
Eigen::SparseMatrix<double> B;
|
||||
igl::random_points_on_mesh(10000,V,F,B,I);
|
||||
P = B*V;
|
||||
Eigen::MatrixXd FN;
|
||||
igl::per_face_normals(V,F,FN);
|
||||
N.resize(P.rows(),3);
|
||||
for(int p = 0;p<I.rows();p++)
|
||||
{
|
||||
N.row(p) = FN.row(I(p));
|
||||
}
|
||||
}
|
||||
// Build octree
|
||||
std::vector<std::vector<int > > O_PI;
|
||||
Eigen::MatrixXi O_CH;
|
||||
Eigen::MatrixXd O_CN;
|
||||
Eigen::VectorXd O_W;
|
||||
igl::octree(P,O_PI,O_CH,O_CN,O_W);
|
||||
Eigen::VectorXd A;
|
||||
{
|
||||
Eigen::MatrixXi I;
|
||||
igl::knn(P,20,O_PI,O_CH,O_CN,O_W,I);
|
||||
// CGAL is only used to help get point areas
|
||||
igl::copyleft::cgal::point_areas(P,I,N,A);
|
||||
}
|
||||
|
||||
if(argc<=1)
|
||||
{
|
||||
// corrupt mesh
|
||||
Eigen::MatrixXd BC;
|
||||
igl::barycenter(V,F,BC);
|
||||
Eigen::MatrixXd OV = V;
|
||||
V.resize(F.rows()*3,3);
|
||||
for(int f = 0;f<F.rows();f++)
|
||||
{
|
||||
for(int c = 0;c<3;c++)
|
||||
{
|
||||
int v = f+c*F.rows();
|
||||
// random rotation about barycenter
|
||||
Eigen::AngleAxisd R(
|
||||
0.5*static_cast <double> (rand()) / static_cast <double> (RAND_MAX),
|
||||
Eigen::Vector3d::Random(3,1));
|
||||
V.row(v) = (OV.row(F(f,c))-BC.row(f))*R.matrix()+BC.row(f);
|
||||
F(f,c) = v;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Generate a list of random query points in the bounding box
|
||||
Eigen::MatrixXd Q = Eigen::MatrixXd::Random(1000000,3);
|
||||
const Eigen::RowVector3d Vmin = V.colwise().minCoeff();
|
||||
const Eigen::RowVector3d Vmax = V.colwise().maxCoeff();
|
||||
const Eigen::RowVector3d Vdiag = Vmax-Vmin;
|
||||
for(int q = 0;q<Q.rows();q++)
|
||||
{
|
||||
Q.row(q) = (Q.row(q).array()*0.5+0.5)*Vdiag.array() + Vmin.array();
|
||||
}
|
||||
|
||||
// Positions of points inside of point cloud P
|
||||
Eigen::MatrixXd QiP;
|
||||
{
|
||||
Eigen::MatrixXd O_CM;
|
||||
Eigen::VectorXd O_R;
|
||||
Eigen::MatrixXd O_EC;
|
||||
printf(" point cloud precomputation: %g secs\n",
|
||||
time([&](){igl::fast_winding_number(P,N,A,O_PI,O_CH,2,O_CM,O_R,O_EC);}));
|
||||
Eigen::VectorXd WiP;
|
||||
printf(" point cloud evaluation: %g secs\n",
|
||||
time([&](){igl::fast_winding_number(P,N,A,O_PI,O_CH,O_CM,O_R,O_EC,Q,2,WiP);}));
|
||||
igl::slice_mask(Q,WiP.array()>0.5,1,QiP);
|
||||
}
|
||||
|
||||
// Positions of points inside of triangle soup (V,F)
|
||||
Eigen::MatrixXd QiV;
|
||||
{
|
||||
igl::FastWindingNumberBVH fwn_bvh;
|
||||
printf("triangle soup precomputation: %g secs\n",
|
||||
time([&](){igl::fast_winding_number(V.cast<float>(),F,2,fwn_bvh);}));
|
||||
Eigen::VectorXf WiV;
|
||||
printf(" triangle soup evaluation: %g secs\n",
|
||||
time([&](){igl::fast_winding_number(fwn_bvh,2,Q.cast<float>(),WiV);}));
|
||||
igl::slice_mask(Q,WiV.array()>0.5,1,QiV);
|
||||
}
|
||||
|
||||
|
||||
// Visualization
|
||||
igl::opengl::glfw::Viewer viewer;
|
||||
// For dislpaying normals as little line segments
|
||||
Eigen::MatrixXd PN(2*P.rows(),3);
|
||||
Eigen::MatrixXi E(P.rows(),2);
|
||||
const double bbd = igl::bounding_box_diagonal(V);
|
||||
for(int p = 0;p<P.rows();p++)
|
||||
{
|
||||
E(p,0) = 2*p;
|
||||
E(p,1) = 2*p+1;
|
||||
PN.row(E(p,0)) = P.row(p);
|
||||
PN.row(E(p,1)) = P.row(p)+bbd*0.01*N.row(p);
|
||||
}
|
||||
|
||||
bool show_P = false;
|
||||
int show_Q = 0;
|
||||
|
||||
int query_data = 0;
|
||||
viewer.data_list[query_data].set_mesh(V,F);
|
||||
viewer.data_list[query_data].clear();
|
||||
viewer.data_list[query_data].point_size = 2;
|
||||
viewer.append_mesh();
|
||||
int object_data = 1;
|
||||
viewer.data_list[object_data].set_mesh(V,F);
|
||||
viewer.data_list[object_data].point_size = 5;
|
||||
|
||||
const auto update = [&]()
|
||||
{
|
||||
viewer.data_list[query_data].clear();
|
||||
switch(show_Q)
|
||||
{
|
||||
case 1:
|
||||
// show all Q
|
||||
viewer.data_list[query_data].set_points(Q,Eigen::RowVector3d(0.996078,0.760784,0.760784));
|
||||
break;
|
||||
case 2:
|
||||
// show all Q inside
|
||||
if(show_P)
|
||||
{
|
||||
viewer.data_list[query_data].set_points(QiP,Eigen::RowVector3d(0.564706,0.847059,0.768627));
|
||||
}else
|
||||
{
|
||||
viewer.data_list[query_data].set_points(QiV,Eigen::RowVector3d(0.564706,0.847059,0.768627));
|
||||
}
|
||||
break;
|
||||
}
|
||||
|
||||
viewer.data_list[object_data].clear();
|
||||
if(show_P)
|
||||
{
|
||||
viewer.data_list[object_data].set_points(P,Eigen::RowVector3d(1,1,1));
|
||||
viewer.data_list[object_data].set_edges(PN,E,Eigen::RowVector3d(0.8,0.8,0.8));
|
||||
}else
|
||||
{
|
||||
viewer.data_list[object_data].set_mesh(V,F);
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
|
||||
viewer.callback_key_pressed =
|
||||
[&](igl::opengl::glfw::Viewer &, unsigned int key, int mod)
|
||||
{
|
||||
switch(key)
|
||||
{
|
||||
default:
|
||||
return false;
|
||||
case '1':
|
||||
show_P = !show_P;
|
||||
break;
|
||||
case '2':
|
||||
show_Q = (show_Q+1) % 3;
|
||||
break;
|
||||
}
|
||||
update();
|
||||
return true;
|
||||
};
|
||||
|
||||
update();
|
||||
viewer.launch();
|
||||
|
||||
}
|
||||
@@ -155,4 +155,7 @@ if(TUTORIALS_CHAPTER7)
|
||||
endif()
|
||||
add_subdirectory("715_MeshImplicitFunction")
|
||||
add_subdirectory("716_HeatGeodesics")
|
||||
if(LIBIGL_WITH_CGAL)
|
||||
add_subdirectory("717_FastWindingNumber")
|
||||
endif()
|
||||
endif()
|
||||
|
||||
Reference in New Issue
Block a user