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.
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.
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/>.
18 * @file MeshVoxelizer.ih
22 * Implementation of inline methods defined in MeshVoxelizer.h
24 * This file is part of the DGtal library.
27/////////////////////////////////////////////////////////////////////////////
28// IMPLEMENTATION of inline methods.
29/////////////////////////////////////////////////////////////////////////////
31/////////////////////////////////////////////////////////////////////////////
32// ----------------------- Standard services --------------------------------
34// ---------------------------------------------------------
35template <typename TDigitalSet, int Separation>
38DGtal::MeshVoxelizer<TDigitalSet,Separation>::distance(const PointR3& M,
42 ASSERT( n.norm()!=0 );
46 distance = n.dot(voxel) - d;
56// ---------------------------------------------------------
57template <typename TDigitalSet, int Separation>
59typename DGtal::MeshVoxelizer<TDigitalSet,Separation>::TriangleOrientation
60DGtal::MeshVoxelizer<TDigitalSet, Separation>::pointIsInside2DTriangle(const PointR2& A,
66 double val1 = (C[1] - A[1])*(v[0] - A[0]) + (C[0] - A[0])*(A[1] - v[1]);
68 double val2 = (B[1] - C[1])*(v[0] - C[0]) + (B[0] - C[0])*(C[1] - v[1]);
70 double val3 = (A[1] - B[1])*(v[0] - B[0]) + (A[0] - B[0])*(B[1] - v[1]);
76 if( ( val1 == 0. && val2 == 0. ) ||
77 ( val1 == 0. && val3 == 0. ) ||
78 ( val2 == 0. && val3 == 0. ) )
80 else if(val1 < 0. || val2 < 0. || val3 < 0.)
82 else if(val1 == 0. || val2 == 0. || val3 == 0.)
88// ---------------------------------------------------------
89template <typename TDigitalSet, int Separation>
92DGtal::MeshVoxelizer<TDigitalSet, Separation>::pointIsInsideVoxel(const PointR3& P,
98 for(int i(0); i < 3; i++)
99 isInside = isInside && ( -0.5 <= PP[i] && PP[i] <= 0.5 );
104// ---------------------------------------------------------
105template <typename TDigitalSet, int Separation>
108DGtal::MeshVoxelizer<TDigitalSet, Separation>::voxelizeTriangle(DigitalSet &outputSet,
113 const std::pair<PointZ3, PointZ3>& bbox)
115 OrientationFunctor orientationFunctor;
117 //geometric predicate
118 PredicateFromOrientationFunctor2<OrientationFunctor> pointPredicate( orientationFunctor );
120 // foreach intersection target
121 for(unsigned int i(0); i < myIntersectionTarget().size(); i++)
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);
129 if(! pointPredicate(AA, BB, CC))
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]++)
138 auto target = myIntersectionTarget(i);
139 // check if points are on different side
141 target.mySecond += v;
143 VectorR3 a2myFirst = A - target.myFirst;
144 VectorR3 a2mySecond = A - target.mySecond;
146 VectorR3 w = target.mySecond - target.myFirst;
147 double den = n.dot(w);
149 bool isSameSide = den == 0 // target on plane
150 || ( a2myFirst.dot(n) * a2mySecond.dot(n) > 0 ); // target on one side
152 // if target is on the same side -> skip current iteration
156 PointR2 pp = myIntersectionTarget.project(i, v);
158 // check if current voxel projection is inside ABC projection
159 if(pointIsInside2DTriangle(AA, BB, CC, pp) != OUTSIDE)
161 if (outputSet.domain().isInside( v ) )
168// ---------------------------------------------------------
169template <typename TDigitalSet, int Separation>
170template <typename MeshPoint>
173DGtal::MeshVoxelizer<TDigitalSet,Separation>::voxelize(DigitalSet &outputSet,
177 const double scaleFactor)
179 std::pair<PointR3, PointR3> bbox_r3;
180 std::pair<PointZ3, PointZ3> bbox_z3;
184 //Scaling + casting to PointR3
191 n = e1.crossProduct(e2).getNormalized();
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 );
201 ASSERT( bbox_r3.first <= bbox_r3.second);
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);});
209 // voxelize current triangle to myDigitalSet
210 voxelizeTriangle( outputSet, A, B, C, n, bbox_z3);
213// ---------------------------------------------------------
214template <typename TDigitalSet, int Separation>
215template <typename MeshPoint>
218DGtal::MeshVoxelizer<TDigitalSet, Separation>::voxelize(DigitalSet &outputSet,
219 const Mesh<MeshPoint> &aMesh,
220 const double scaleFactor)
222 DigitalSet rawEmpty{outputSet.domain()};
223 typedef typename Mesh<MeshPoint>::Index Index;
224 typedef std::vector<Index> MeshFace;
227#pragma omp parallel for schedule(dynamic)
229 for(unsigned int i = 0; i < aMesh.nbFaces(); i++)
231 DigitalSet currentSet{rawEmpty};
232 MeshFace currentFace = aMesh.getFace(i);
233 for(unsigned int j=0; j + 2 < currentFace.size(); ++j)
235 voxelize(currentSet, aMesh.getVertex(currentFace[0]),
236 aMesh.getVertex(currentFace[j+1]),
237 aMesh.getVertex(currentFace[j+2]),
245 outputSet += currentSet;