DGtal  1.5.beta
MeshVoxelizer.ih
1 /**
2  * This program is free software: you can redistribute it and/or modify
3  * it under the terms of the GNU Lesser General Public License as
4  * published by the Free Software Foundation, either version 3 of the
5  * License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  *
15  **/
16 
17 /**
18  * @file MeshVoxelizer.ih
19  *
20  * @date 2016/01/24
21  *
22  * Implementation of inline methods defined in MeshVoxelizer.h
23  *
24  * This file is part of the DGtal library.
25  */
26 
27 /////////////////////////////////////////////////////////////////////////////
28 // IMPLEMENTATION of inline methods.
29 /////////////////////////////////////////////////////////////////////////////
30 #include <algorithm>
31 /////////////////////////////////////////////////////////////////////////////
32 // ----------------------- Standard services --------------------------------
33 
34 // ---------------------------------------------------------
35 template <typename TDigitalSet, int Separation>
36 inline
37 double
38 DGtal::MeshVoxelizer<TDigitalSet,Separation>::distance(const PointR3& M,
39  const VectorR3& n,
40  const PointZ3& voxel)
41 {
42  ASSERT( n.norm()!=0 );
43  double distance = 0.;
44  double d = n.dot(M);
45 
46  distance = n.dot(voxel) - d;
47 
48  if(distance < 0)
49  distance *= -1;
50 
51  distance /= n.norm();
52 
53  return distance;
54 }
55 
56 // ---------------------------------------------------------
57 template <typename TDigitalSet, int Separation>
58 inline
59 typename DGtal::MeshVoxelizer<TDigitalSet,Separation>::TriangleOrientation
60 DGtal::MeshVoxelizer<TDigitalSet, Separation>::pointIsInside2DTriangle(const PointR2& A,
61  const PointR2& B,
62  const PointR2& C,
63  const PointR2& v)
64 {
65  // AC
66  double val1 = (C[1] - A[1])*(v[0] - A[0]) + (C[0] - A[0])*(A[1] - v[1]);
67  // CB
68  double val2 = (B[1] - C[1])*(v[0] - C[0]) + (B[0] - C[0])*(C[1] - v[1]);
69  // BA
70  double val3 = (A[1] - B[1])*(v[0] - B[0]) + (A[0] - B[0])*(B[1] - v[1]);
71 
72  // 0 : outside
73  // 1 : inside
74  // 2 : on edge
75  // 3 : on vertex
76  if( ( val1 == 0. && val2 == 0. ) ||
77  ( val1 == 0. && val3 == 0. ) ||
78  ( val2 == 0. && val3 == 0. ) )
79  return TRIANGLE_ONVERTEX;
80  else if(val1 < 0. || val2 < 0. || val3 < 0.)
81  return TRIANGLE_OUTSIDE;
82  else if(val1 == 0. || val2 == 0. || val3 == 0.)
83  return TRIANGLE_ONEDGE;
84  else
85  return TRIANGLE_INSIDE;
86 }
87 
88 // ---------------------------------------------------------
89 template <typename TDigitalSet, int Separation>
90 inline
91 bool
92 DGtal::MeshVoxelizer<TDigitalSet, Separation>::pointIsInsideVoxel(const PointR3& P,
93  const PointZ3& v)
94 {
95  bool isInside = true;
96 
97  PointR3 PP = P - v;
98  for(int i(0); i < 3; i++)
99  isInside = isInside && ( -0.5 <= PP[i] && PP[i] <= 0.5 );
100 
101  return isInside;
102 }
103 
104 // ---------------------------------------------------------
105 template <typename TDigitalSet, int Separation>
106 inline
107 void
108 DGtal::MeshVoxelizer<TDigitalSet, Separation>::voxelizeTriangle(DigitalSet &outputSet,
109  const PointR3& A,
110  const PointR3& B,
111  const PointR3& C,
112  const VectorR3& n,
113  const std::pair<PointZ3, PointZ3>& bbox)
114 {
115  OrientationFunctor orientationFunctor;
116 
117  //geometric predicate
118  PredicateFromOrientationFunctor2<OrientationFunctor> pointPredicate( orientationFunctor );
119 
120  // foreach intersection target
121  for(unsigned int i(0); i < myIntersectionTarget().size(); i++)
122  {
123  // 2D projection of A ; B ; C
124  PointR2 AA = myIntersectionTarget.project(i, A);
125  PointR2 BB = myIntersectionTarget.project(i, B);
126  PointR2 CC = myIntersectionTarget.project(i, C);
127 
128  // check orientation
129  if(! pointPredicate(AA, BB, CC))
130  std::swap(AA, CC);
131 
132  // traverse all voxel of current face bounding box
133  PointZ3 v = bbox.first;
134  for(; v[1] <= bbox.second[1]; v[1]++)
135  for(v[0] = bbox.first[0]; v[0] <= bbox.second[0]; v[0]++)
136  for(v[2] = bbox.first[2]; v[2] <= bbox.second[2]; v[2]++)
137  {
138  auto target = myIntersectionTarget(i);
139  // check if points are on different side
140  target.myFirst += v;
141  target.mySecond += v;
142 
143  VectorR3 a2myFirst = A - target.myFirst;
144  VectorR3 a2mySecond = A - target.mySecond;
145 
146  VectorR3 w = target.mySecond - target.myFirst;
147  double den = n.dot(w);
148 
149  bool isSameSide = den == 0 // target on plane
150  || ( a2myFirst.dot(n) * a2mySecond.dot(n) > 0 ); // target on one side
151 
152  // if target is on the same side -> skip current iteration
153  if( isSameSide )
154  continue;
155 
156  PointR2 pp = myIntersectionTarget.project(i, v);
157 
158  // check if current voxel projection is inside ABC projection
159  if(pointIsInside2DTriangle(AA, BB, CC, pp) != TRIANGLE_OUTSIDE)
160  {
161  if (outputSet.domain().isInside( v ) )
162  outputSet.insert(v);
163  }
164  }
165  }
166 }
167 
168 // ---------------------------------------------------------
169 template <typename TDigitalSet, int Separation>
170 template <typename MeshPoint>
171 inline
172 void
173 DGtal::MeshVoxelizer<TDigitalSet,Separation>::voxelize(DigitalSet &outputSet,
174  const MeshPoint &a,
175  const MeshPoint &b,
176  const MeshPoint &c,
177  const double scaleFactor)
178 {
179  std::pair<PointR3, PointR3> bbox_r3;
180  std::pair<PointZ3, PointZ3> bbox_z3;
181  VectorR3 n, e1, e2;
182  PointR3 A, B, C;
183 
184  //Scaling + casting to PointR3
185  A = a*scaleFactor;
186  B = b*scaleFactor;
187  C = c*scaleFactor;
188 
189  e1 = B - A;
190  e2 = C - A;
191  n = e1.crossProduct(e2).getNormalized();
192 
193  //Boundingbox
194  bbox_r3.first = A;
195  bbox_r3.second = A;
196  bbox_r3.first = bbox_r3.first.inf( B );
197  bbox_r3.first = bbox_r3.first.inf( C );
198  bbox_r3.second = bbox_r3.second.sup( B );
199  bbox_r3.second = bbox_r3.second.sup( C );
200 
201  ASSERT( bbox_r3.first <= bbox_r3.second);
202 
203  //Rounding the r3 bbox into the z3 bbox
204  std::transform( bbox_r3.first.begin(), bbox_r3.first.end(), bbox_z3.first.begin(),
205  [](typename PointR3::Component cc) { return std::floor(cc);});
206  std::transform( bbox_r3.second.begin(), bbox_r3.second.end(), bbox_z3.second.begin(),
207  [](typename PointR3::Component cc) { return std::ceil(cc);});
208 
209  // voxelize current triangle to myDigitalSet
210  voxelizeTriangle( outputSet, A, B, C, n, bbox_z3);
211 }
212 
213 // ---------------------------------------------------------
214 template <typename TDigitalSet, int Separation>
215 template <typename MeshPoint>
216 inline
217 void
218 DGtal::MeshVoxelizer<TDigitalSet, Separation>::voxelize(DigitalSet &outputSet,
219  const Mesh<MeshPoint> &aMesh,
220  const double scaleFactor)
221 {
222  DigitalSet rawEmpty{outputSet.domain()};
223  typedef typename Mesh<MeshPoint>::Index Index;
224  typedef std::vector<Index> MeshFace;
225 
226 #ifdef WITH_OPENMP
227 #pragma omp parallel for schedule(dynamic)
228 #endif
229  for(int i = 0; i < (int)aMesh.nbFaces(); i++)
230  {
231  DigitalSet currentSet{rawEmpty};
232  MeshFace currentFace = aMesh.getFace(i);
233  for(size_t j=0; j + 2 < currentFace.size(); ++j)
234  {
235  voxelize(currentSet, aMesh.getVertex(currentFace[0]),
236  aMesh.getVertex(currentFace[j+1]),
237  aMesh.getVertex(currentFace[j+2]),
238  scaleFactor);
239  }
240 
241 #ifdef WITH_OPENMP
242  #pragma omp critical
243 #endif
244  {
245  outputSet += currentSet;
246  }
247  }
248 }