1 // This file is part of libigl, a simple c++ geometry processing library.
2 //
3 // Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public License
6 // v. 2.0. If a copy of the MPL was not distributed with this file, You can
7 // obtain one at http://mozilla.org/MPL/2.0/.
8 #include "minkowski_sum.h"
9 #include "mesh_boolean.h"
10
11 #include "../../slice.h"
12 #include "../../slice_mask.h"
13 #include "../../LinSpaced.h"
14 #include "../../unique_rows.h"
15 #include "../../get_seconds.h"
16 #include "../../edges.h"
17 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
18 #include <cassert>
19 #include <vector>
20 #include <iostream>
21
22
23 template <
24 typename DerivedVA,
25 typename DerivedFA,
26 typename DerivedVB,
27 typename DerivedFB,
28 typename DerivedW,
29 typename DerivedG,
30 typename DerivedJ>
minkowski_sum(const Eigen::MatrixBase<DerivedVA> & VA,const Eigen::MatrixBase<DerivedFA> & FA,const Eigen::MatrixBase<DerivedVB> & VB,const Eigen::MatrixBase<DerivedFB> & FB,const bool resolve_overlaps,Eigen::PlainObjectBase<DerivedW> & W,Eigen::PlainObjectBase<DerivedG> & G,Eigen::PlainObjectBase<DerivedJ> & J)31 IGL_INLINE void igl::copyleft::cgal::minkowski_sum(
32 const Eigen::MatrixBase<DerivedVA> & VA,
33 const Eigen::MatrixBase<DerivedFA> & FA,
34 const Eigen::MatrixBase<DerivedVB> & VB,
35 const Eigen::MatrixBase<DerivedFB> & FB,
36 const bool resolve_overlaps,
37 Eigen::PlainObjectBase<DerivedW> & W,
38 Eigen::PlainObjectBase<DerivedG> & G,
39 Eigen::PlainObjectBase<DerivedJ> & J)
40 {
41 using namespace std;
42 using namespace Eigen;
43 assert(FA.cols() == 3 && "FA must contain a closed triangle mesh");
44 assert(FB.cols() <= FA.cols() &&
45 "FB must contain lower diemnsional simplices than FA");
46 const auto tictoc = []()->double
47 {
48 static double t_start;
49 double now = igl::get_seconds();
50 double interval = now-t_start;
51 t_start = now;
52 return interval;
53 };
54 tictoc();
55 Matrix<typename DerivedFB::Scalar,Dynamic,2> EB;
56 edges(FB,EB);
57 Matrix<typename DerivedFA::Scalar,Dynamic,2> EA(0,2);
58 if(FB.cols() == 3)
59 {
60 edges(FA,EA);
61 }
62 // number of copies of A along edges of B
63 const int n_ab = EB.rows();
64 // number of copies of B along edges of A
65 const int n_ba = EA.rows();
66
67 vector<DerivedW> vW(n_ab + n_ba);
68 vector<DerivedG> vG(n_ab + n_ba);
69 vector<DerivedJ> vJ(n_ab + n_ba);
70 vector<int> offsets(n_ab + n_ba + 1);
71 offsets[0] = 0;
72 // sweep A along edges of B
73 for(int e = 0;e<n_ab;e++)
74 {
75 Matrix<typename DerivedJ::Scalar,Dynamic,1> eJ;
76 minkowski_sum(
77 VA,
78 FA,
79 VB.row(EB(e,0)).eval(),
80 VB.row(EB(e,1)).eval(),
81 false,
82 vW[e],
83 vG[e],
84 eJ);
85 assert(vG[e].rows() == eJ.rows());
86 assert(eJ.cols() == 1);
87 vJ[e].resize(vG[e].rows(),2);
88 vJ[e].col(0) = eJ;
89 vJ[e].col(1).setConstant(e);
90 offsets[e+1] = offsets[e] + vW[e].rows();
91 }
92 // sweep B along edges of A
93 for(int e = 0;e<n_ba;e++)
94 {
95 Matrix<typename DerivedJ::Scalar,Dynamic,1> eJ;
96 const int ee = n_ab+e;
97 minkowski_sum(
98 VB,
99 FB,
100 VA.row(EA(e,0)).eval(),
101 VA.row(EA(e,1)).eval(),
102 false,
103 vW[ee],
104 vG[ee],
105 eJ);
106 vJ[ee].resize(vG[ee].rows(),2);
107 vJ[ee].col(0) = eJ.array() + (FA.rows()+1);
108 vJ[ee].col(1).setConstant(ee);
109 offsets[ee+1] = offsets[ee] + vW[ee].rows();
110 }
111 // Combine meshes
112 int n=0,m=0;
113 for_each(vW.begin(),vW.end(),[&n](const DerivedW & w){n+=w.rows();});
114 for_each(vG.begin(),vG.end(),[&m](const DerivedG & g){m+=g.rows();});
115 assert(n == offsets.back());
116
117 W.resize(n,3);
118 G.resize(m,3);
119 J.resize(m,2);
120 {
121 int m_off = 0,n_off = 0;
122 for(int i = 0;i<vG.size();i++)
123 {
124 W.block(n_off,0,vW[i].rows(),3) = vW[i];
125 G.block(m_off,0,vG[i].rows(),3) = vG[i].array()+offsets[i];
126 J.block(m_off,0,vJ[i].rows(),2) = vJ[i];
127 n_off += vW[i].rows();
128 m_off += vG[i].rows();
129 }
130 assert(n == n_off);
131 assert(m == m_off);
132 }
133 if(resolve_overlaps)
134 {
135 Eigen::Matrix<typename DerivedJ::Scalar, Eigen::Dynamic,1> SJ;
136 mesh_boolean(
137 DerivedW(W),
138 DerivedG(G),
139 Matrix<typename DerivedW::Scalar,Dynamic,Dynamic>(),
140 Matrix<typename DerivedG::Scalar,Dynamic,Dynamic>(),
141 MESH_BOOLEAN_TYPE_UNION,
142 W,
143 G,
144 SJ);
145 slice(DerivedJ(J),SJ,1,J);
146 }
147 }
148
149 template <
150 typename DerivedVA,
151 typename DerivedFA,
152 typename sType, int sCols, int sOptions,
153 typename dType, int dCols, int dOptions,
154 typename DerivedW,
155 typename DerivedG,
156 typename DerivedJ>
minkowski_sum(const Eigen::MatrixBase<DerivedVA> & VA,const Eigen::MatrixBase<DerivedFA> & FA,const Eigen::Matrix<sType,1,sCols,sOptions> & s,const Eigen::Matrix<dType,1,dCols,dOptions> & d,const bool resolve_overlaps,Eigen::PlainObjectBase<DerivedW> & W,Eigen::PlainObjectBase<DerivedG> & G,Eigen::PlainObjectBase<DerivedJ> & J)157 IGL_INLINE void igl::copyleft::cgal::minkowski_sum(
158 const Eigen::MatrixBase<DerivedVA> & VA,
159 const Eigen::MatrixBase<DerivedFA> & FA,
160 const Eigen::Matrix<sType,1,sCols,sOptions> & s,
161 const Eigen::Matrix<dType,1,dCols,dOptions> & d,
162 const bool resolve_overlaps,
163 Eigen::PlainObjectBase<DerivedW> & W,
164 Eigen::PlainObjectBase<DerivedG> & G,
165 Eigen::PlainObjectBase<DerivedJ> & J)
166 {
167 using namespace Eigen;
168 using namespace std;
169 assert(s.cols() == 3 && "s should be a 3d point");
170 assert(d.cols() == 3 && "d should be a 3d point");
171 // silly base case
172 if(FA.size() == 0)
173 {
174 W.resize(0,3);
175 G.resize(0,3);
176 return;
177 }
178 const int dim = VA.cols();
179 assert(dim == 3 && "dim must be 3D");
180 assert(s.size() == 3 && "s must be 3D point");
181 assert(d.size() == 3 && "d must be 3D point");
182 // segment vector
183 const CGAL::Vector_3<CGAL::Epeck> v(d(0)-s(0),d(1)-s(1),d(2)-s(2));
184 // number of vertices
185 const int n = VA.rows();
186 // duplicate vertices at s and d, we'll remove unreferernced later
187 W.resize(2*n,dim);
188 for(int i = 0;i<n;i++)
189 {
190 for(int j = 0;j<dim;j++)
191 {
192 W (i,j) = VA(i,j) + s(j);
193 W(i+n,j) = VA(i,j) + d(j);
194 }
195 }
196 // number of faces
197 const int m = FA.rows();
198 //// Mask whether positive dot product, or negative: because of exactly zero,
199 //// these are not necessarily complementary
200 // Nevermind, actually P = !N
201 Matrix<bool,Dynamic,1> P(m,1),N(m,1);
202 // loop over faces
203 int mp = 0,mn = 0;
204 for(int f = 0;f<m;f++)
205 {
206 const CGAL::Plane_3<CGAL::Epeck> plane(
207 CGAL::Point_3<CGAL::Epeck>(VA(FA(f,0),0),VA(FA(f,0),1),VA(FA(f,0),2)),
208 CGAL::Point_3<CGAL::Epeck>(VA(FA(f,1),0),VA(FA(f,1),1),VA(FA(f,1),2)),
209 CGAL::Point_3<CGAL::Epeck>(VA(FA(f,2),0),VA(FA(f,2),1),VA(FA(f,2),2)));
210 const auto normal = plane.orthogonal_vector();
211 const auto dt = normal * v;
212 if(dt > 0)
213 {
214 P(f) = true;
215 N(f) = false;
216 mp++;
217 }else
218 //}else if(dt < 0)
219 {
220 P(f) = false;
221 N(f) = true;
222 mn++;
223 //}else
224 //{
225 // P(f) = false;
226 // N(f) = false;
227 }
228 }
229
230 typedef Matrix<typename DerivedG::Scalar,Dynamic,Dynamic> MatrixXI;
231 typedef Matrix<typename DerivedG::Scalar,Dynamic,1> VectorXI;
232 MatrixXI GT(mp+mn,3);
233 GT<< slice_mask(FA,N,1), slice_mask((FA.array()+n).eval(),P,1);
234 // J indexes FA for parts at s and m+FA for parts at d
235 J.derived() = igl::LinSpaced<DerivedJ >(m,0,m-1);
236 DerivedJ JT(mp+mn);
237 JT << slice_mask(J,P,1), slice_mask(J,N,1);
238 JT.block(mp,0,mn,1).array()+=m;
239
240 // Original non-co-planar faces with positively oriented reversed
241 MatrixXI BA(mp+mn,3);
242 BA << slice_mask(FA,P,1).rowwise().reverse(), slice_mask(FA,N,1);
243 // Quads along **all** sides
244 MatrixXI GQ((mp+mn)*3,4);
245 GQ<<
246 BA.col(1), BA.col(0), BA.col(0).array()+n, BA.col(1).array()+n,
247 BA.col(2), BA.col(1), BA.col(1).array()+n, BA.col(2).array()+n,
248 BA.col(0), BA.col(2), BA.col(2).array()+n, BA.col(0).array()+n;
249
250 MatrixXI uGQ;
251 VectorXI S,sI,sJ;
252 // Inputs:
253 // F #F by d list of polygons
254 // Outputs:
255 // S #uF list of signed incidences for each unique face
256 // uF #uF by d list of unique faces
257 // I #uF index vector so that uF = sort(F,2)(I,:)
258 // J #F index vector so that sort(F,2) = uF(J,:)
259 [](
260 const MatrixXI & F,
261 VectorXI & S,
262 MatrixXI & uF,
263 VectorXI & I,
264 VectorXI & J)
265 {
266 const int m = F.rows();
267 const int d = F.cols();
268 MatrixXI sF = F;
269 const auto MN = sF.rowwise().minCoeff().eval();
270 // rotate until smallest index is first
271 for(int p = 0;p<d;p++)
272 {
273 for(int f = 0;f<m;f++)
274 {
275 if(sF(f,0) != MN(f))
276 {
277 for(int r = 0;r<d-1;r++)
278 {
279 std::swap(sF(f,r),sF(f,r+1));
280 }
281 }
282 }
283 }
284 // swap orienation so that last index is greater than first
285 for(int f = 0;f<m;f++)
286 {
287 if(sF(f,d-1) < sF(f,1))
288 {
289 sF.block(f,1,1,d-1) = sF.block(f,1,1,d-1).reverse().eval();
290 }
291 }
292 Matrix<bool,Dynamic,1> M = Matrix<bool,Dynamic,1>::Zero(m,1);
293 {
294 VectorXI P = igl::LinSpaced<VectorXI >(d,0,d-1);
295 for(int p = 0;p<d;p++)
296 {
297 for(int f = 0;f<m;f++)
298 {
299 bool all = true;
300 for(int r = 0;r<d;r++)
301 {
302 all = all && (sF(f,P(r)) == F(f,r));
303 }
304 M(f) = M(f) || all;
305 }
306 for(int r = 0;r<d-1;r++)
307 {
308 std::swap(P(r),P(r+1));
309 }
310 }
311 }
312 unique_rows(sF,uF,I,J);
313 S = VectorXI::Zero(uF.rows(),1);
314 assert(m == J.rows());
315 for(int f = 0;f<m;f++)
316 {
317 S(J(f)) += M(f) ? 1 : -1;
318 }
319 }(MatrixXI(GQ),S,uGQ,sI,sJ);
320 assert(S.rows() == uGQ.rows());
321 const int nq = (S.array().abs()==2).count();
322 GQ.resize(nq,4);
323 {
324 int k = 0;
325 for(int q = 0;q<uGQ.rows();q++)
326 {
327 switch(S(q))
328 {
329 case -2:
330 GQ.row(k++) = uGQ.row(q).reverse().eval();
331 break;
332 case 2:
333 GQ.row(k++) = uGQ.row(q);
334 break;
335 default:
336 // do not add
337 break;
338 }
339 }
340 assert(nq == k);
341 }
342
343 G.resize(GT.rows()+2*GQ.rows(),3);
344 G<<
345 GT,
346 GQ.col(0), GQ.col(1), GQ.col(2),
347 GQ.col(0), GQ.col(2), GQ.col(3);
348 J.resize(JT.rows()+2*GQ.rows(),1);
349 J<<JT,DerivedJ::Constant(2*GQ.rows(),1,2*m+1);
350 if(resolve_overlaps)
351 {
352 Eigen::Matrix<typename DerivedJ::Scalar, Eigen::Dynamic,1> SJ;
353 mesh_boolean(
354 DerivedW(W),DerivedG(G),
355 Matrix<typename DerivedVA::Scalar,Dynamic,Dynamic>(),MatrixXI(),
356 MESH_BOOLEAN_TYPE_UNION,
357 W,G,SJ);
358 J.derived() = slice(DerivedJ(J),SJ,1);
359 }
360 }
361
362 template <
363 typename DerivedVA,
364 typename DerivedFA,
365 typename sType, int sCols, int sOptions,
366 typename dType, int dCols, int dOptions,
367 typename DerivedW,
368 typename DerivedG,
369 typename DerivedJ>
minkowski_sum(const Eigen::MatrixBase<DerivedVA> & VA,const Eigen::MatrixBase<DerivedFA> & FA,const Eigen::Matrix<sType,1,sCols,sOptions> & s,const Eigen::Matrix<dType,1,dCols,dOptions> & d,Eigen::PlainObjectBase<DerivedW> & W,Eigen::PlainObjectBase<DerivedG> & G,Eigen::PlainObjectBase<DerivedJ> & J)370 IGL_INLINE void igl::copyleft::cgal::minkowski_sum(
371 const Eigen::MatrixBase<DerivedVA> & VA,
372 const Eigen::MatrixBase<DerivedFA> & FA,
373 const Eigen::Matrix<sType,1,sCols,sOptions> & s,
374 const Eigen::Matrix<dType,1,dCols,dOptions> & d,
375 Eigen::PlainObjectBase<DerivedW> & W,
376 Eigen::PlainObjectBase<DerivedG> & G,
377 Eigen::PlainObjectBase<DerivedJ> & J)
378 {
379 return minkowski_sum(VA,FA,s,d,true,W,G,J);
380 }
381
382 #ifdef IGL_STATIC_LIBRARY
383 // Explicit template instantiation
384 // generated by autoexplicit.sh
385 template void igl::copyleft::cgal::minkowski_sum<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, CGAL::Lazy_exact_nt<CGAL::Gmpq>, 3, 1, CGAL::Lazy_exact_nt<CGAL::Gmpq>, 3, 1, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 1, 3, 1, 1, 3> const&, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 1, 3, 1, 1, 3> const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 1, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
386 // generated by autoexplicit.sh
387 template void igl::copyleft::cgal::minkowski_sum<
388 Eigen::Matrix<float, -1, 3, 1, -1, 3>,
389 Eigen::Matrix<int, -1, 3, 1, -1, 3>,
390 double, 3, 1,
391 float, 3, 1,
392 Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 1, -1, -1>,
393 Eigen::Matrix<int, -1, -1, 0, -1, -1>,
394 Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, Eigen::Matrix<float, 1, 3, 1, 1, 3> const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 1, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
395 #endif
396