DGtal 1.3.0
Loading...
Searching...
No Matches
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// ---------------------------------------------------------
35template <typename TDigitalSet, int Separation>
36inline
37double
38DGtal::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// ---------------------------------------------------------
57template <typename TDigitalSet, int Separation>
58inline
59typename DGtal::MeshVoxelizer<TDigitalSet,Separation>::TriangleOrientation
60DGtal::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// ---------------------------------------------------------
89template <typename TDigitalSet, int Separation>
90inline
91bool
92DGtal::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// ---------------------------------------------------------
105template <typename TDigitalSet, int Separation>
106inline
107void
108DGtal::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) != OUTSIDE)
160 {
161 if (outputSet.domain().isInside( v ) )
162 outputSet.insert(v);
163 }
164 }
165 }
166}
167
168// ---------------------------------------------------------
169template <typename TDigitalSet, int Separation>
170template <typename MeshPoint>
171inline
172void
173DGtal::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// ---------------------------------------------------------
214template <typename TDigitalSet, int Separation>
215template <typename MeshPoint>
216inline
217void
218DGtal::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(unsigned int i = 0; i < aMesh.nbFaces(); i++)
230 {
231 DigitalSet currentSet{rawEmpty};
232 MeshFace currentFace = aMesh.getFace(i);
233 for(unsigned int 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}