1/* Copyright (c) 2011 Khaled Mamou (kmamou at gmail dot com)
2 All rights reserved.
3
4
5 Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
6
7 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
8
9 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
10
11 3. The names of the contributors may not be used to endorse or promote products derived from this software without specific prior written permission.
12
13 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
14 */
15#pragma once
16#ifndef VHACD_VOLUME_H
17#define VHACD_VOLUME_H
18#include "vhacdMesh.h"
19#include "vhacdVector.h"
20#include <assert.h>
21
22#ifdef _MSC_VER
23#pragma warning(push)
24#pragma warning(disable:4456 4701)
25#endif
26
27namespace VHACD {
28
29enum VOXEL_VALUE {
30 PRIMITIVE_UNDEFINED = 0,
31 PRIMITIVE_OUTSIDE_SURFACE = 1,
32 PRIMITIVE_INSIDE_SURFACE = 2,
33 PRIMITIVE_ON_SURFACE = 3
34};
35
36struct Voxel {
37public:
38 short m_coord[3];
39 short m_data;
40};
41
42class PrimitiveSet {
43public:
44 virtual ~PrimitiveSet(){};
45 virtual PrimitiveSet* Create() const = 0;
46 virtual const size_t GetNPrimitives() const = 0;
47 virtual const size_t GetNPrimitivesOnSurf() const = 0;
48 virtual const size_t GetNPrimitivesInsideSurf() const = 0;
49 virtual const double GetEigenValue(AXIS axis) const = 0;
50 virtual const double ComputeMaxVolumeError() const = 0;
51 virtual const double ComputeVolume() const = 0;
52 virtual void Clip(const Plane& plane, PrimitiveSet* const positivePart,
53 PrimitiveSet* const negativePart) const = 0;
54 virtual void Intersect(const Plane& plane, SArray<Vec3<double> >* const positivePts,
55 SArray<Vec3<double> >* const negativePts, const size_t sampling) const = 0;
56 virtual void ComputeExteriorPoints(const Plane& plane, const Mesh& mesh,
57 SArray<Vec3<double> >* const exteriorPts) const = 0;
58 virtual void ComputeClippedVolumes(const Plane& plane, double& positiveVolume,
59 double& negativeVolume) const = 0;
60 virtual void SelectOnSurface(PrimitiveSet* const onSurfP) const = 0;
61 virtual void ComputeConvexHull(Mesh& meshCH, const size_t sampling = 1) const = 0;
62 virtual void ComputeBB() = 0;
63 virtual void ComputePrincipalAxes() = 0;
64 virtual void AlignToPrincipalAxes() = 0;
65 virtual void RevertAlignToPrincipalAxes() = 0;
66 virtual void Convert(Mesh& mesh, const VOXEL_VALUE value) const = 0;
67 const Mesh& GetConvexHull() const { return m_convexHull; };
68 Mesh& GetConvexHull() { return m_convexHull; };
69private:
70 Mesh m_convexHull;
71};
72
73//!
74class VoxelSet : public PrimitiveSet {
75 friend class Volume;
76
77public:
78 //! Destructor.
79 ~VoxelSet(void);
80 //! Constructor.
81 VoxelSet();
82
83 const size_t GetNPrimitives() const { return m_voxels.Size(); }
84 const size_t GetNPrimitivesOnSurf() const { return m_numVoxelsOnSurface; }
85 const size_t GetNPrimitivesInsideSurf() const { return m_numVoxelsInsideSurface; }
86 const double GetEigenValue(AXIS axis) const { return m_D[axis][axis]; }
87 const double ComputeVolume() const { return m_unitVolume * m_voxels.Size(); }
88 const double ComputeMaxVolumeError() const { return m_unitVolume * m_numVoxelsOnSurface; }
89 const Vec3<short>& GetMinBBVoxels() const { return m_minBBVoxels; }
90 const Vec3<short>& GetMaxBBVoxels() const { return m_maxBBVoxels; }
91 const Vec3<double>& GetMinBB() const { return m_minBB; }
92 const double& GetScale() const { return m_scale; }
93 const double& GetUnitVolume() const { return m_unitVolume; }
94 Vec3<double> GetPoint(Vec3<short> voxel) const
95 {
96 return Vec3<double>(voxel[0] * m_scale + m_minBB[0],
97 voxel[1] * m_scale + m_minBB[1],
98 voxel[2] * m_scale + m_minBB[2]);
99 }
100 Vec3<double> GetPoint(const Voxel& voxel) const
101 {
102 return Vec3<double>(voxel.m_coord[0] * m_scale + m_minBB[0],
103 voxel.m_coord[1] * m_scale + m_minBB[1],
104 voxel.m_coord[2] * m_scale + m_minBB[2]);
105 }
106 Vec3<double> GetPoint(Vec3<double> voxel) const
107 {
108 return Vec3<double>(voxel[0] * m_scale + m_minBB[0],
109 voxel[1] * m_scale + m_minBB[1],
110 voxel[2] * m_scale + m_minBB[2]);
111 }
112 void GetPoints(const Voxel& voxel, Vec3<double>* const pts) const;
113 void ComputeConvexHull(Mesh& meshCH, const size_t sampling = 1) const;
114 void Clip(const Plane& plane, PrimitiveSet* const positivePart, PrimitiveSet* const negativePart) const;
115 void Intersect(const Plane& plane, SArray<Vec3<double> >* const positivePts,
116 SArray<Vec3<double> >* const negativePts, const size_t sampling) const;
117 void ComputeExteriorPoints(const Plane& plane, const Mesh& mesh,
118 SArray<Vec3<double> >* const exteriorPts) const;
119 void ComputeClippedVolumes(const Plane& plane, double& positiveVolume, double& negativeVolume) const;
120 void SelectOnSurface(PrimitiveSet* const onSurfP) const;
121 void ComputeBB();
122 void Convert(Mesh& mesh, const VOXEL_VALUE value) const;
123 void ComputePrincipalAxes();
124 PrimitiveSet* Create() const
125 {
126 return new VoxelSet();
127 }
128 void AlignToPrincipalAxes(){};
129 void RevertAlignToPrincipalAxes(){};
130 Voxel* const GetVoxels() { return m_voxels.Data(); }
131 const Voxel* const GetVoxels() const { return m_voxels.Data(); }
132
133private:
134 size_t m_numVoxelsOnSurface;
135 size_t m_numVoxelsInsideSurface;
136 Vec3<double> m_minBB;
137 double m_scale;
138 SArray<Voxel, 8> m_voxels;
139 double m_unitVolume;
140 Vec3<double> m_minBBPts;
141 Vec3<double> m_maxBBPts;
142 Vec3<short> m_minBBVoxels;
143 Vec3<short> m_maxBBVoxels;
144 Vec3<short> m_barycenter;
145 double m_Q[3][3];
146 double m_D[3][3];
147 Vec3<double> m_barycenterPCA;
148};
149
150struct Tetrahedron {
151public:
152 Vec3<double> m_pts[4];
153 unsigned char m_data;
154};
155
156//!
157class TetrahedronSet : public PrimitiveSet {
158 friend class Volume;
159
160public:
161 //! Destructor.
162 ~TetrahedronSet(void);
163 //! Constructor.
164 TetrahedronSet();
165
166 const size_t GetNPrimitives() const { return m_tetrahedra.Size(); }
167 const size_t GetNPrimitivesOnSurf() const { return m_numTetrahedraOnSurface; }
168 const size_t GetNPrimitivesInsideSurf() const { return m_numTetrahedraInsideSurface; }
169 const Vec3<double>& GetMinBB() const { return m_minBB; }
170 const Vec3<double>& GetMaxBB() const { return m_maxBB; }
171 const Vec3<double>& GetBarycenter() const { return m_barycenter; }
172 const double GetEigenValue(AXIS axis) const { return m_D[axis][axis]; }
173 const double GetSacle() const { return m_scale; }
174 const double ComputeVolume() const;
175 const double ComputeMaxVolumeError() const;
176 void ComputeConvexHull(Mesh& meshCH, const size_t sampling = 1) const;
177 void ComputePrincipalAxes();
178 void AlignToPrincipalAxes();
179 void RevertAlignToPrincipalAxes();
180 void Clip(const Plane& plane, PrimitiveSet* const positivePart, PrimitiveSet* const negativePart) const;
181 void Intersect(const Plane& plane, SArray<Vec3<double> >* const positivePts,
182 SArray<Vec3<double> >* const negativePts, const size_t sampling) const;
183 void ComputeExteriorPoints(const Plane& plane, const Mesh& mesh,
184 SArray<Vec3<double> >* const exteriorPts) const;
185 void ComputeClippedVolumes(const Plane& plane, double& positiveVolume, double& negativeVolume) const;
186 void SelectOnSurface(PrimitiveSet* const onSurfP) const;
187 void ComputeBB();
188 void Convert(Mesh& mesh, const VOXEL_VALUE value) const;
189 inline bool Add(Tetrahedron& tetrahedron);
190 PrimitiveSet* Create() const
191 {
192 return new TetrahedronSet();
193 }
194 static const double EPS;
195
196private:
197 void AddClippedTetrahedra(const Vec3<double> (&pts)[10], const int32_t nPts);
198
199 size_t m_numTetrahedraOnSurface;
200 size_t m_numTetrahedraInsideSurface;
201 double m_scale;
202 Vec3<double> m_minBB;
203 Vec3<double> m_maxBB;
204 Vec3<double> m_barycenter;
205 SArray<Tetrahedron, 8> m_tetrahedra;
206 double m_Q[3][3];
207 double m_D[3][3];
208};
209
210//!
211class Volume {
212public:
213 //! Destructor.
214 ~Volume(void);
215
216 //! Constructor.
217 Volume();
218
219 //! Voxelize
220 template <class T>
221 void Voxelize(const T* const points, const uint32_t stridePoints, const uint32_t nPoints,
222 const int32_t* const triangles, const uint32_t strideTriangles, const uint32_t nTriangles,
223 const size_t dim, const Vec3<double>& barycenter, const double (&rot)[3][3]);
224 unsigned char& GetVoxel(const size_t i, const size_t j, const size_t k)
225 {
226 assert(i < m_dim[0] || i >= 0);
227 assert(j < m_dim[0] || j >= 0);
228 assert(k < m_dim[0] || k >= 0);
229 return m_data[i + j * m_dim[0] + k * m_dim[0] * m_dim[1]];
230 }
231 const unsigned char& GetVoxel(const size_t i, const size_t j, const size_t k) const
232 {
233 assert(i < m_dim[0] || i >= 0);
234 assert(j < m_dim[0] || j >= 0);
235 assert(k < m_dim[0] || k >= 0);
236 return m_data[i + j * m_dim[0] + k * m_dim[0] * m_dim[1]];
237 }
238 const size_t GetNPrimitivesOnSurf() const { return m_numVoxelsOnSurface; }
239 const size_t GetNPrimitivesInsideSurf() const { return m_numVoxelsInsideSurface; }
240 void Convert(Mesh& mesh, const VOXEL_VALUE value) const;
241 void Convert(VoxelSet& vset) const;
242 void Convert(TetrahedronSet& tset) const;
243 void AlignToPrincipalAxes(double (&rot)[3][3]) const;
244
245private:
246 void FillOutsideSurface(const size_t i0, const size_t j0, const size_t k0, const size_t i1,
247 const size_t j1, const size_t k1);
248 void FillInsideSurface();
249 template <class T>
250 void ComputeBB(const T* const points, const uint32_t stridePoints, const uint32_t nPoints,
251 const Vec3<double>& barycenter, const double (&rot)[3][3]);
252 void Allocate();
253 void Free();
254
255 Vec3<double> m_minBB;
256 Vec3<double> m_maxBB;
257 double m_scale;
258 size_t m_dim[3]; //>! dim
259 size_t m_numVoxelsOnSurface;
260 size_t m_numVoxelsInsideSurface;
261 size_t m_numVoxelsOutsideSurface;
262 unsigned char* m_data;
263};
264int32_t TriBoxOverlap(const Vec3<double>& boxcenter, const Vec3<double>& boxhalfsize, const Vec3<double>& triver0,
265 const Vec3<double>& triver1, const Vec3<double>& triver2);
266template <class T>
267inline void ComputeAlignedPoint(const T* const points, const uint32_t idx, const Vec3<double>& barycenter,
268 const double (&rot)[3][3], Vec3<double>& pt){};
269template <>
270inline void ComputeAlignedPoint<float>(const float* const points, const uint32_t idx, const Vec3<double>& barycenter, const double (&rot)[3][3], Vec3<double>& pt)
271{
272 double x = points[idx + 0] - barycenter[0];
273 double y = points[idx + 1] - barycenter[1];
274 double z = points[idx + 2] - barycenter[2];
275 pt[0] = rot[0][0] * x + rot[1][0] * y + rot[2][0] * z;
276 pt[1] = rot[0][1] * x + rot[1][1] * y + rot[2][1] * z;
277 pt[2] = rot[0][2] * x + rot[1][2] * y + rot[2][2] * z;
278}
279template <>
280inline void ComputeAlignedPoint<double>(const double* const points, const uint32_t idx, const Vec3<double>& barycenter, const double (&rot)[3][3], Vec3<double>& pt)
281{
282 double x = points[idx + 0] - barycenter[0];
283 double y = points[idx + 1] - barycenter[1];
284 double z = points[idx + 2] - barycenter[2];
285 pt[0] = rot[0][0] * x + rot[1][0] * y + rot[2][0] * z;
286 pt[1] = rot[0][1] * x + rot[1][1] * y + rot[2][1] * z;
287 pt[2] = rot[0][2] * x + rot[1][2] * y + rot[2][2] * z;
288}
289template <class T>
290void Volume::ComputeBB(const T* const points, const uint32_t stridePoints, const uint32_t nPoints,
291 const Vec3<double>& barycenter, const double (&rot)[3][3])
292{
293 Vec3<double> pt;
294 ComputeAlignedPoint(points, 0, barycenter, rot, pt);
295 m_maxBB = pt;
296 m_minBB = pt;
297 for (uint32_t v = 1; v < nPoints; ++v) {
298 ComputeAlignedPoint(points, v * stridePoints, barycenter, rot, pt);
299 for (int32_t i = 0; i < 3; ++i) {
300 if (pt[i] < m_minBB[i])
301 m_minBB[i] = pt[i];
302 else if (pt[i] > m_maxBB[i])
303 m_maxBB[i] = pt[i];
304 }
305 }
306}
307template <class T>
308void Volume::Voxelize(const T* const points, const uint32_t stridePoints, const uint32_t nPoints,
309 const int32_t* const triangles, const uint32_t strideTriangles, const uint32_t nTriangles,
310 const size_t dim, const Vec3<double>& barycenter, const double (&rot)[3][3])
311{
312 if (nPoints == 0) {
313 return;
314 }
315 ComputeBB(points, stridePoints, nPoints, barycenter, rot);
316
317 double d[3] = { m_maxBB[0] - m_minBB[0], m_maxBB[1] - m_minBB[1], m_maxBB[2] - m_minBB[2] };
318 double r;
319 if (d[0] >= d[1] && d[0] >= d[2]) {
320 r = d[0];
321 m_dim[0] = dim;
322 m_dim[1] = 2 + static_cast<size_t>(dim * d[1] / d[0]);
323 m_dim[2] = 2 + static_cast<size_t>(dim * d[2] / d[0]);
324 }
325 else if (d[1] >= d[0] && d[1] >= d[2]) {
326 r = d[1];
327 m_dim[1] = dim;
328 m_dim[0] = 2 + static_cast<size_t>(dim * d[0] / d[1]);
329 m_dim[2] = 2 + static_cast<size_t>(dim * d[2] / d[1]);
330 }
331 else {
332 r = d[2];
333 m_dim[2] = dim;
334 m_dim[0] = 2 + static_cast<size_t>(dim * d[0] / d[2]);
335 m_dim[1] = 2 + static_cast<size_t>(dim * d[1] / d[2]);
336 }
337
338 m_scale = r / (dim - 1);
339 double invScale = (dim - 1) / r;
340
341 Allocate();
342 m_numVoxelsOnSurface = 0;
343 m_numVoxelsInsideSurface = 0;
344 m_numVoxelsOutsideSurface = 0;
345
346 Vec3<double> p[3];
347 size_t i, j, k;
348 size_t i0, j0, k0;
349 size_t i1, j1, k1;
350 Vec3<double> boxcenter;
351 Vec3<double> pt;
352 const Vec3<double> boxhalfsize(0.5, 0.5, 0.5);
353 for (size_t t = 0, ti = 0; t < nTriangles; ++t, ti += strideTriangles) {
354 Vec3<int32_t> tri(triangles[ti + 0],
355 triangles[ti + 1],
356 triangles[ti + 2]);
357 for (int32_t c = 0; c < 3; ++c) {
358 ComputeAlignedPoint(points, tri[c] * stridePoints, barycenter, rot, pt);
359 p[c][0] = (pt[0] - m_minBB[0]) * invScale;
360 p[c][1] = (pt[1] - m_minBB[1]) * invScale;
361 p[c][2] = (pt[2] - m_minBB[2]) * invScale;
362 i = static_cast<size_t>(p[c][0] + 0.5);
363 j = static_cast<size_t>(p[c][1] + 0.5);
364 k = static_cast<size_t>(p[c][2] + 0.5);
365 assert(i < m_dim[0] && i >= 0 && j < m_dim[1] && j >= 0 && k < m_dim[2] && k >= 0);
366
367 if (c == 0) {
368 i0 = i1 = i;
369 j0 = j1 = j;
370 k0 = k1 = k;
371 }
372 else {
373 if (i < i0)
374 i0 = i;
375 if (j < j0)
376 j0 = j;
377 if (k < k0)
378 k0 = k;
379 if (i > i1)
380 i1 = i;
381 if (j > j1)
382 j1 = j;
383 if (k > k1)
384 k1 = k;
385 }
386 }
387 if (i0 > 0)
388 --i0;
389 if (j0 > 0)
390 --j0;
391 if (k0 > 0)
392 --k0;
393 if (i1 < m_dim[0])
394 ++i1;
395 if (j1 < m_dim[1])
396 ++j1;
397 if (k1 < m_dim[2])
398 ++k1;
399 for (size_t i = i0; i < i1; ++i) {
400 boxcenter[0] = (double)i;
401 for (size_t j = j0; j < j1; ++j) {
402 boxcenter[1] = (double)j;
403 for (size_t k = k0; k < k1; ++k) {
404 boxcenter[2] = (double)k;
405 int32_t res = TriBoxOverlap(boxcenter, boxhalfsize, p[0], p[1], p[2]);
406 unsigned char& value = GetVoxel(i, j, k);
407 if (res == 1 && value == PRIMITIVE_UNDEFINED) {
408 value = PRIMITIVE_ON_SURFACE;
409 ++m_numVoxelsOnSurface;
410 }
411 }
412 }
413 }
414 }
415 FillOutsideSurface(0, 0, 0, m_dim[0], m_dim[1], 1);
416 FillOutsideSurface(0, 0, m_dim[2] - 1, m_dim[0], m_dim[1], m_dim[2]);
417 FillOutsideSurface(0, 0, 0, m_dim[0], 1, m_dim[2]);
418 FillOutsideSurface(0, m_dim[1] - 1, 0, m_dim[0], m_dim[1], m_dim[2]);
419 FillOutsideSurface(0, 0, 0, 1, m_dim[1], m_dim[2]);
420 FillOutsideSurface(m_dim[0] - 1, 0, 0, m_dim[0], m_dim[1], m_dim[2]);
421 FillInsideSurface();
422}
423}
424
425#ifdef _MSC_VER
426#pragma warning(pop)
427#endif
428
429
430#endif // VHACD_VOLUME_H
431