1 // This file is part of libigl, a simple c++ geometry processing library.
2 //
3 // Copyright (C) 2015 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 "slice_tets.h"
9 #include "LinSpaced.h"
10 #include "sort.h"
11 #include "edges.h"
12 #include "slice.h"
13 #include "cat.h"
14 #include "ismember.h"
15 #include "unique_rows.h"
16 #include <cassert>
17 #include <algorithm>
18 #include <vector>
19
20 template <
21 typename DerivedV,
22 typename DerivedT,
23 typename DerivedS,
24 typename DerivedSV,
25 typename DerivedSF,
26 typename DerivedJ,
27 typename BCType>
slice_tets(const Eigen::MatrixBase<DerivedV> & V,const Eigen::MatrixBase<DerivedT> & T,const Eigen::MatrixBase<DerivedS> & S,Eigen::PlainObjectBase<DerivedSV> & SV,Eigen::PlainObjectBase<DerivedSF> & SF,Eigen::PlainObjectBase<DerivedJ> & J,Eigen::SparseMatrix<BCType> & BC)28 IGL_INLINE void igl::slice_tets(
29 const Eigen::MatrixBase<DerivedV>& V,
30 const Eigen::MatrixBase<DerivedT>& T,
31 const Eigen::MatrixBase<DerivedS> & S,
32 Eigen::PlainObjectBase<DerivedSV>& SV,
33 Eigen::PlainObjectBase<DerivedSF>& SF,
34 Eigen::PlainObjectBase<DerivedJ>& J,
35 Eigen::SparseMatrix<BCType> & BC)
36 {
37 Eigen::MatrixXi sE;
38 Eigen::Matrix<typename DerivedSV::Scalar,Eigen::Dynamic,1> lambda;
39 igl::slice_tets(V,T,S,SV,SF,J,sE,lambda);
40 const int ns = SV.rows();
41 std::vector<Eigen::Triplet<BCType> > BCIJV(ns*2);
42 for(int i = 0;i<ns;i++)
43 {
44 BCIJV[2*i+0] = Eigen::Triplet<BCType>(i,sE(i,0), lambda(i));
45 BCIJV[2*i+1] = Eigen::Triplet<BCType>(i,sE(i,1),1.0-lambda(i));
46 }
47 BC.resize(SV.rows(),V.rows());
48 BC.setFromTriplets(BCIJV.begin(),BCIJV.end());
49 }
50
51 template <
52 typename DerivedV,
53 typename DerivedT,
54 typename DerivedS,
55 typename DerivedSV,
56 typename DerivedSF,
57 typename DerivedJ>
slice_tets(const Eigen::MatrixBase<DerivedV> & V,const Eigen::MatrixBase<DerivedT> & T,const Eigen::MatrixBase<DerivedS> & S,Eigen::PlainObjectBase<DerivedSV> & SV,Eigen::PlainObjectBase<DerivedSF> & SF,Eigen::PlainObjectBase<DerivedJ> & J)58 IGL_INLINE void igl::slice_tets(
59 const Eigen::MatrixBase<DerivedV>& V,
60 const Eigen::MatrixBase<DerivedT>& T,
61 const Eigen::MatrixBase<DerivedS> & S,
62 Eigen::PlainObjectBase<DerivedSV>& SV,
63 Eigen::PlainObjectBase<DerivedSF>& SF,
64 Eigen::PlainObjectBase<DerivedJ>& J)
65 {
66 Eigen::MatrixXi sE;
67 Eigen::Matrix<typename DerivedSV::Scalar,Eigen::Dynamic,1> lambda;
68 igl::slice_tets(V,T,S,SV,SF,J,sE,lambda);
69 }
70
71 template <
72 typename DerivedV,
73 typename DerivedT,
74 typename DerivedS,
75 typename DerivedSV,
76 typename DerivedSF,
77 typename DerivedJ,
78 typename DerivedsE,
79 typename Derivedlambda>
slice_tets(const Eigen::MatrixBase<DerivedV> & V,const Eigen::MatrixBase<DerivedT> & T,const Eigen::MatrixBase<DerivedS> & S,Eigen::PlainObjectBase<DerivedSV> & SV,Eigen::PlainObjectBase<DerivedSF> & SF,Eigen::PlainObjectBase<DerivedJ> & J,Eigen::PlainObjectBase<DerivedsE> & sE,Eigen::PlainObjectBase<Derivedlambda> & lambda)80 IGL_INLINE void igl::slice_tets(
81 const Eigen::MatrixBase<DerivedV>& V,
82 const Eigen::MatrixBase<DerivedT>& T,
83 const Eigen::MatrixBase<DerivedS> & S,
84 Eigen::PlainObjectBase<DerivedSV>& SV,
85 Eigen::PlainObjectBase<DerivedSF>& SF,
86 Eigen::PlainObjectBase<DerivedJ>& J,
87 Eigen::PlainObjectBase<DerivedsE>& sE,
88 Eigen::PlainObjectBase<Derivedlambda>& lambda)
89 {
90
91 using namespace Eigen;
92 using namespace std;
93 assert(V.cols() == 3 && "V should be #V by 3");
94 assert(T.cols() == 4 && "T should be #T by 4");
95
96 static const Eigen::Matrix<int,12,4> flipped_order =
97 (Eigen::Matrix<int,12,4>(12,4)<<
98 3,2,0,1,
99 3,1,2,0,
100 3,0,1,2,
101 2,3,1,0,
102 2,1,0,3,
103 2,0,3,1,
104 1,3,0,2,
105 1,2,3,0,
106 1,0,2,3,
107 0,3,2,1,
108 0,2,1,3,
109 0,1,3,2
110 ).finished();
111
112 // number of tets
113 const size_t m = T.rows();
114
115 typedef typename DerivedS::Scalar Scalar;
116 typedef typename DerivedT::Scalar Index;
117 typedef Matrix<Scalar,Dynamic,1> VectorXS;
118 typedef Matrix<Scalar,Dynamic,4> MatrixX4S;
119 typedef Matrix<Scalar,Dynamic,3> MatrixX3S;
120 typedef Matrix<Scalar,Dynamic,2> MatrixX2S;
121 typedef Matrix<Index,Dynamic,4> MatrixX4I;
122 typedef Matrix<Index,Dynamic,3> MatrixX3I;
123 typedef Matrix<Index,Dynamic,2> MatrixX2I;
124 typedef Matrix<Index,Dynamic,1> VectorXI;
125 typedef Array<bool,Dynamic,1> ArrayXb;
126
127 MatrixX4S IT(m,4);
128 for(size_t t = 0;t<m;t++)
129 {
130 for(size_t c = 0;c<4;c++)
131 {
132 IT(t,c) = S(T(t,c));
133 }
134 }
135
136 // Essentially, just a glorified slice(X,1)
137 //
138 // Inputs:
139 // T #T by 4 list of tet indices into V
140 // IT #IT by 4 list of isosurface values at each tet
141 // I #I list of bools whether to grab data corresponding to each tet
142 const auto & extract_rows = [](
143 const MatrixBase<DerivedT> & T,
144 const MatrixX4S & IT,
145 const ArrayXb & I,
146 MatrixX4I & TI,
147 MatrixX4S & ITI,
148 VectorXI & JI)
149 {
150 const Index num_I = std::count(I.data(),I.data()+I.size(),true);
151 TI.resize(num_I,4);
152 ITI.resize(num_I,4);
153 JI.resize(num_I,1);
154 {
155 size_t k = 0;
156 for(size_t t = 0;t<(size_t)T.rows();t++)
157 {
158 if(I(t))
159 {
160 TI.row(k) = T.row(t);
161 ITI.row(k) = IT.row(t);
162 JI(k) = t;
163 k++;
164 }
165 }
166 assert(k == num_I);
167 }
168 };
169
170 ArrayXb I13 = (IT.array()<0).rowwise().count()==1;
171 ArrayXb I31 = (IT.array()>0).rowwise().count()==1;
172 ArrayXb I22 = (IT.array()<0).rowwise().count()==2;
173 MatrixX4I T13,T31,T22;
174 MatrixX4S IT13,IT31,IT22;
175 VectorXI J13,J31,J22;
176 extract_rows(T,IT,I13,T13,IT13,J13);
177 extract_rows(T,IT,I31,T31,IT31,J31);
178 extract_rows(T,IT,I22,T22,IT22,J22);
179
180 const auto & apply_sort4 = [] (
181 const MatrixX4I & T,
182 const MatrixX4I & sJ,
183 MatrixX4I & sT)
184 {
185 sT.resize(T.rows(),4);
186 for(size_t t = 0;t<(size_t)T.rows();t++)
187 {
188 for(size_t c = 0;c<4;c++)
189 {
190 sT(t,c) = T(t,sJ(t,c));
191 }
192 }
193 };
194
195 const auto & apply_sort2 = [] (
196 const MatrixX2I & E,
197 const MatrixX2I & sJ,
198 Eigen::PlainObjectBase<DerivedsE>& sE)
199 {
200 sE.resize(E.rows(),2);
201 for(size_t t = 0;t<(size_t)E.rows();t++)
202 {
203 for(size_t c = 0;c<2;c++)
204 {
205 sE(t,c) = E(t,sJ(t,c));
206 }
207 }
208 };
209
210 const auto & one_below = [&apply_sort4](
211 const MatrixX4I & T,
212 const MatrixX4S & IT,
213 MatrixX2I & U,
214 MatrixX3I & SF)
215 {
216 // Number of tets
217 const size_t m = T.rows();
218 if(m == 0)
219 {
220 U.resize(0,2);
221 SF.resize(0,3);
222 return;
223 }
224 MatrixX4S sIT;
225 MatrixX4I sJ;
226 sort(IT,2,true,sIT,sJ);
227 MatrixX4I sT;
228 apply_sort4(T,sJ,sT);
229 U.resize(3*m,2);
230 U<<
231 sT.col(0),sT.col(1),
232 sT.col(0),sT.col(2),
233 sT.col(0),sT.col(3);
234 SF.resize(m,3);
235 for(size_t c = 0;c<3;c++)
236 {
237 SF.col(c) =
238 igl::LinSpaced<
239 Eigen::Matrix<typename DerivedSF::Scalar,Eigen::Dynamic,1> >
240 (m,0+c*m,(m-1)+c*m);
241 }
242 ArrayXb flip;
243 {
244 VectorXi _;
245 ismember_rows(sJ,flipped_order,flip,_);
246 }
247 for(int i = 0;i<m;i++)
248 {
249 if(flip(i))
250 {
251 SF.row(i) = SF.row(i).reverse().eval();
252 }
253 }
254 };
255
256 const auto & two_below = [&apply_sort4](
257 const MatrixX4I & T,
258 const MatrixX4S & IT,
259 MatrixX2I & U,
260 MatrixX3I & SF)
261 {
262 // Number of tets
263 const size_t m = T.rows();
264 if(m == 0)
265 {
266 U.resize(0,2);
267 SF.resize(0,3);
268 return;
269 }
270 MatrixX4S sIT;
271 MatrixX4I sJ;
272 sort(IT,2,true,sIT,sJ);
273 MatrixX4I sT;
274 apply_sort4(T,sJ,sT);
275 U.resize(4*m,2);
276 U<<
277 sT.col(0),sT.col(2),
278 sT.col(0),sT.col(3),
279 sT.col(1),sT.col(2),
280 sT.col(1),sT.col(3);
281 SF.resize(2*m,3);
282 SF.block(0,0,m,1) = igl::LinSpaced<VectorXI >(m,0+0*m,(m-1)+0*m);
283 SF.block(0,1,m,1) = igl::LinSpaced<VectorXI >(m,0+1*m,(m-1)+1*m);
284 SF.block(0,2,m,1) = igl::LinSpaced<VectorXI >(m,0+3*m,(m-1)+3*m);
285 SF.block(m,0,m,1) = igl::LinSpaced<VectorXI >(m,0+0*m,(m-1)+0*m);
286 SF.block(m,1,m,1) = igl::LinSpaced<VectorXI >(m,0+3*m,(m-1)+3*m);
287 SF.block(m,2,m,1) = igl::LinSpaced<VectorXI >(m,0+2*m,(m-1)+2*m);
288 ArrayXb flip;
289 {
290 VectorXi _;
291 ismember_rows(sJ,flipped_order,flip,_);
292 }
293 for(int i = 0;i<m;i++)
294 {
295 if(flip(i))
296 {
297 SF.row(i ) = SF.row(i ).reverse().eval();
298 SF.row(i+m) = SF.row(i+m).reverse().eval();
299 }
300 }
301 };
302
303 MatrixX3I SF13,SF31,SF22;
304 MatrixX2I U13,U31,U22;
305 one_below(T13, IT13,U13,SF13);
306 one_below(T31,-IT31,U31,SF31);
307 two_below(T22, IT22,U22,SF22);
308 // https://forum.kde.org/viewtopic.php?f=74&t=107974
309 const MatrixX2I U =
310 (MatrixX2I(U13.rows()+ U31.rows()+ U22.rows(),2)<<U13,U31,U22).finished();
311 MatrixX2I sU;
312 {
313 MatrixX2I _;
314 sort(U,2,true,sU,_);
315 }
316 MatrixX2I E;
317 VectorXI uI,uJ;
318 unique_rows(sU,E,uI,uJ);
319 MatrixX2S IE(E.rows(),2);
320 for(size_t t = 0;t<E.rows();t++)
321 {
322 for(size_t c = 0;c<2;c++)
323 {
324 IE(t,c) = S(E(t,c));
325 }
326 }
327 MatrixX2S sIE;
328 MatrixX2I sJ;
329 sort(IE,2,true,sIE,sJ);
330 apply_sort2(E,sJ,sE);
331 lambda = sIE.col(1).array() / (sIE.col(1)-sIE.col(0)).array();
332 SV.resize(sE.rows(),V.cols());
333 for(int e = 0;e<sE.rows();e++)
334 {
335 SV.row(e) = V.row(sE(e,0)).template cast<Scalar>()*lambda(e) +
336 V.row(sE(e,1)).template cast<Scalar>()*(1.0-lambda(e));
337 }
338 SF.resize( SF13.rows()+SF31.rows()+SF22.rows(),3);
339 SF<<
340 SF13,
341 U13.rows()+ SF31.rowwise().reverse().array(),
342 U13.rows()+U31.rows()+SF22.array();
343
344 std::for_each(
345 SF.data(),
346 SF.data()+SF.size(),
347 [&uJ](typename DerivedSF::Scalar & i){i=uJ(i);});
348
349 J.resize(SF.rows());
350 J<<J13,J31,J22,J22;
351 }
352
353 #ifdef IGL_STATIC_LIBRARY
354 // Explicit template instantiation
355 template void igl::slice_tets<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::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, double>(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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::SparseMatrix<double, 0, int>&);
356 #endif
357