DGtal  0.9.4.1
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 ONVERTEX;
80  else if(val1 < 0. || val2 < 0. || val3 < 0.)
81  return OUTSIDE;
82  else if(val1 == 0. || val2 == 0. || val3 == 0.)
83  return ONEDGE;
84  else
85  return 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<PointR3, PointR3>& 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) != 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;
180  VectorR3 n, e1, e2;
181  PointR3 A, B, C;
182 
183  //Scaling + casting to PointR3
184  A = a*scaleFactor;
185  B = b*scaleFactor;
186  C = c*scaleFactor;
187 
188  e1 = B - A;
189  e2 = C - A;
190  n = e1.crossProduct(e2).getNormalized();
191 
192  //Boundingbox
193  bbox.first = A;
194  bbox.second = A;
195  bbox.first = bbox.first.inf( B );
196  bbox.first = bbox.first.inf( C );
197  bbox.second = bbox.second.sup( B );
198  bbox.second = bbox.second.sup( C );
199 
200  ASSERT( bbox.first <= bbox.second);
201 
202  //Rounding the bbox
203  std::transform( bbox.first.begin(), bbox.first.end(), bbox.first.begin(),
204  [](typename PointR3::Component cc) { return std::floor(cc);});
205  std::transform( bbox.second.begin(), bbox.second.end(), bbox.second.begin(),
206  [](typename PointR3::Component cc) { return std::ceil(cc);});
207 
208  // voxelize current triangle to myDigitalSet
209  voxelizeTriangle( outputSet, A, B, C, n, bbox);
210 }
211 
212 // ---------------------------------------------------------
213 template <typename TDigitalSet, int Separation>
214 template <typename MeshPoint>
215 inline
216 void
217 DGtal::MeshVoxelizer<TDigitalSet, Separation>::voxelize(DigitalSet &outputSet,
218  const Mesh<MeshPoint> &aMesh,
219  const double scaleFactor)
220 {
221  DigitalSet rawEmpty{outputSet.domain()};
222 #ifdef WITH_OPENMP
223 #pragma omp parallel for schedule(dynamic)
224 #endif
225  for(unsigned int i = 0; i < aMesh.nbFaces(); i++)
226  {
227  DigitalSet currentSet{rawEmpty};
228  MeshFace currentFace = aMesh.getFace(i);
229  for(unsigned int j=0; j + 2 < currentFace.size(); ++j)
230  {
231  voxelize(currentSet, aMesh.getVertex(currentFace[0]),
232  aMesh.getVertex(currentFace[j+1]),
233  aMesh.getVertex(currentFace[j+2]),
234  scaleFactor);
235  }
236 
237 #ifdef WITH_OPENMP
238  #pragma omp critical
239 #endif
240  {
241  outputSet += currentSet;
242  }
243  }
244 }