|
|
@@ -0,0 +1,208 @@
|
|
|
+#include "minkowski_sum.h"
|
|
|
+#include "mesh_boolean.h"
|
|
|
+
|
|
|
+#include "../../slice_mask.h"
|
|
|
+#include "../../unique.h"
|
|
|
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
|
|
|
+
|
|
|
+template <
|
|
|
+ typename DerivedV,
|
|
|
+ typename DerivedF,
|
|
|
+ typename Deriveds,
|
|
|
+ typename Derivedd,
|
|
|
+ typename DerivedW,
|
|
|
+ typename DerivedG,
|
|
|
+ typename DerivedJ>
|
|
|
+IGL_INLINE void igl::copyleft::boolean::minkowski_sum(
|
|
|
+ const Eigen::PlainObjectBase<DerivedV> & V,
|
|
|
+ const Eigen::PlainObjectBase<DerivedF> & F,
|
|
|
+ const Eigen::PlainObjectBase<Deriveds> & s,
|
|
|
+ const Eigen::PlainObjectBase<Derivedd> & d,
|
|
|
+ Eigen::PlainObjectBase<DerivedW> & W,
|
|
|
+ Eigen::PlainObjectBase<DerivedG> & G,
|
|
|
+ Eigen::PlainObjectBase<DerivedJ> & J)
|
|
|
+{
|
|
|
+ using namespace Eigen;
|
|
|
+ using namespace std;
|
|
|
+ // silly base case
|
|
|
+ if(F.size() == 0)
|
|
|
+ {
|
|
|
+ W.resize(0,3);
|
|
|
+ G.resize(0,3);
|
|
|
+ return;
|
|
|
+ }
|
|
|
+ const int dim = V.cols();
|
|
|
+ assert(dim == 3 && "dim must be 3D");
|
|
|
+ assert(s.size() == 3 && "s must be 3D point");
|
|
|
+ assert(d.size() == 3 && "d must be 3D point");
|
|
|
+ // segment vector
|
|
|
+ const CGAL::Vector_3<CGAL::Epeck> v(d(0)-s(0),d(1)-s(1),d(2)-s(2));
|
|
|
+ // number of vertices
|
|
|
+ const int n = V.rows();
|
|
|
+ // duplicate vertices at s and d, we'll remove unreferernced later
|
|
|
+ DerivedW WW(2*n,dim);
|
|
|
+ for(int i = 0;i<n;i++)
|
|
|
+ {
|
|
|
+ for(int j = 0;j<dim;j++)
|
|
|
+ {
|
|
|
+ WW (i,j) = V(i,j) + s(j);
|
|
|
+ WW(i+n,j) = V(i,j) + d(j);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ // number of faces
|
|
|
+ const int m = F.rows();
|
|
|
+ // Mask whether positive dot product, or negative: because of exactly zero,
|
|
|
+ // these are not necessarily complementary
|
|
|
+ Matrix<bool,Dynamic,1> P(m,1),N(m,1);
|
|
|
+ // loop over faces
|
|
|
+ int mp = 0,mn = 0;
|
|
|
+ for(int f = 0;f<m;f++)
|
|
|
+ {
|
|
|
+ const CGAL::Plane_3<CGAL::Epeck> plane(
|
|
|
+ CGAL::Point_3<CGAL::Epeck>(V(F(f,0),0),V(F(f,0),1),V(F(f,0),2)),
|
|
|
+ CGAL::Point_3<CGAL::Epeck>(V(F(f,1),0),V(F(f,1),1),V(F(f,1),2)),
|
|
|
+ CGAL::Point_3<CGAL::Epeck>(V(F(f,2),0),V(F(f,2),1),V(F(f,2),2)));
|
|
|
+ const auto normal = plane.orthogonal_vector();
|
|
|
+ const auto dt = normal * v;
|
|
|
+ if(dt > 0)
|
|
|
+ {
|
|
|
+ P(f) = true;
|
|
|
+ N(f) = false;
|
|
|
+ mp++;
|
|
|
+ }else if(dt < 0)
|
|
|
+ {
|
|
|
+ P(f) = false;
|
|
|
+ N(f) = true;
|
|
|
+ mn++;
|
|
|
+ }else
|
|
|
+ {
|
|
|
+ P(f) = false;
|
|
|
+ N(f) = false;
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ typedef Matrix<typename DerivedG::Scalar,Dynamic,Dynamic> MatrixXI;
|
|
|
+ typedef Matrix<typename DerivedG::Scalar,Dynamic,1> VectorXI;
|
|
|
+ MatrixXI GT(mp+mn,3);
|
|
|
+ GT<< slice_mask(F,N,1), slice_mask((F.array()+n).eval(),P,1);
|
|
|
+ // J indexes F for parts at s and m+F for parts at d
|
|
|
+ J = DerivedJ::LinSpaced(m,0,m-1);
|
|
|
+ DerivedJ JT(mp+mn);
|
|
|
+ JT << slice_mask(J,P,1), slice_mask(J,N,1);
|
|
|
+ JT.block(mp,0,mn,1).array()+=m;
|
|
|
+
|
|
|
+ // Original non-co-planar faces with positively oriented reversed
|
|
|
+ MatrixXI BA(mp+mn,3);
|
|
|
+ BA << slice_mask(F,P,1).rowwise().reverse(), slice_mask(F,N,1);
|
|
|
+ // Quads along **all** sides
|
|
|
+ MatrixXI GQ((mp+mn)*3,4);
|
|
|
+ GQ<<
|
|
|
+ BA.col(1), BA.col(0), BA.col(0).array()+n, BA.col(1).array()+n,
|
|
|
+ BA.col(2), BA.col(1), BA.col(1).array()+n, BA.col(2).array()+n,
|
|
|
+ BA.col(0), BA.col(2), BA.col(2).array()+n, BA.col(0).array()+n;
|
|
|
+
|
|
|
+ MatrixXI uGQ;
|
|
|
+ VectorXI S,sI,sJ;
|
|
|
+ //const auto & total_signed_distance =
|
|
|
+ [](
|
|
|
+ const MatrixXI & F,
|
|
|
+ VectorXI & S,
|
|
|
+ MatrixXI & uF,
|
|
|
+ VectorXI & I,
|
|
|
+ VectorXI & J)
|
|
|
+ {
|
|
|
+ const int m = F.rows();
|
|
|
+ const int d = F.cols();
|
|
|
+ MatrixXI sF = F;
|
|
|
+ const auto MN = sF.rowwise().minCoeff().eval();
|
|
|
+ // rotate until smallest index is first
|
|
|
+ for(int p = 0;p<d;p++)
|
|
|
+ {
|
|
|
+ for(int f = 0;f<m;f++)
|
|
|
+ {
|
|
|
+ if(sF(f,0) != MN(f))
|
|
|
+ {
|
|
|
+ for(int r = 0;r<d-1;r++)
|
|
|
+ {
|
|
|
+ std::swap(sF(f,r),sF(f,r+1));
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ // swap orienation
|
|
|
+ for(int f = 0;f<m;f++)
|
|
|
+ {
|
|
|
+ if(sF(f,d-1) < sF(f,1))
|
|
|
+ {
|
|
|
+ sF.block(f,1,1,d-1) = sF.block(f,1,1,d-1).reverse().eval();
|
|
|
+ }
|
|
|
+ }
|
|
|
+ Matrix<bool,Dynamic,1> M = Matrix<bool,Dynamic,1>::Zero(m,1);
|
|
|
+ {
|
|
|
+ VectorXI P = VectorXI::LinSpaced(d,0,d-1);
|
|
|
+ for(int p = 0;p<d;p++)
|
|
|
+ {
|
|
|
+ for(int f = 0;f<m;f++)
|
|
|
+ {
|
|
|
+ bool all = true;
|
|
|
+ for(int r = 0;r<d;r++)
|
|
|
+ {
|
|
|
+ all = all && (sF(f,P(r)) == F(f,r));
|
|
|
+ }
|
|
|
+ M(f) = M(f) || all;
|
|
|
+ }
|
|
|
+ for(int r = 0;r<d-1;r++)
|
|
|
+ {
|
|
|
+ std::swap(P(r),P(r+1));
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ unique_rows(sF,uF,I,J);
|
|
|
+ S = VectorXI::Zero(uF.rows(),1);
|
|
|
+ assert(m == J.rows());
|
|
|
+ for(int f = 0;f<m;f++)
|
|
|
+ {
|
|
|
+ S(J(f)) += M(f) ? 1 : -1;
|
|
|
+ }
|
|
|
+ }(MatrixXI(GQ),S,uGQ,sI,sJ);
|
|
|
+ assert(S.rows() == uGQ.rows());
|
|
|
+ const int nq = (S.array().abs()==2).count();
|
|
|
+ GQ.resize(nq,4);
|
|
|
+ {
|
|
|
+ int k = 0;
|
|
|
+ for(int q = 0;q<uGQ.rows();q++)
|
|
|
+ {
|
|
|
+ switch(S(q))
|
|
|
+ {
|
|
|
+#warning "If (V,F) is non-manifold, does this handle all cases?"
|
|
|
+ case -2:
|
|
|
+ GQ.row(k++) = uGQ.row(q).reverse().eval();
|
|
|
+ break;
|
|
|
+ case 2:
|
|
|
+ GQ.row(k++) = uGQ.row(q);
|
|
|
+ break;
|
|
|
+ default:
|
|
|
+ // do not add
|
|
|
+ break;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ assert(nq == k);
|
|
|
+ }
|
|
|
+
|
|
|
+ MatrixXI GG(GT.rows()+2*GQ.rows(),3);
|
|
|
+ GG<<
|
|
|
+ GT,
|
|
|
+ GQ.col(0), GQ.col(1), GQ.col(2),
|
|
|
+ GQ.col(0), GQ.col(2), GQ.col(3);
|
|
|
+ J.resize(JT.rows()+2*GQ.rows(),1);
|
|
|
+ J<<JT,DerivedJ::Constant(2*GQ.rows(),1,2*m+1);
|
|
|
+ {
|
|
|
+ DerivedJ SJ;
|
|
|
+ mesh_boolean(
|
|
|
+ WW,GG,
|
|
|
+ Matrix<typename DerivedV::Scalar,Dynamic,Dynamic>(),MatrixXI(),
|
|
|
+ MESH_BOOLEAN_TYPE_UNION,
|
|
|
+ W,G,SJ);
|
|
|
+ J = slice(DerivedJ(J),SJ,1);
|
|
|
+ }
|
|
|
+}
|
Alec Jacobson