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
16#ifndef _CRT_SECURE_NO_WARNINGS
17#define _CRT_SECURE_NO_WARNINGS
18#endif
19
20#include <algorithm>
21#include <fstream>
22#include <iomanip>
23#include <limits>
24#include <sstream>
25#if _OPENMP
26#include <omp.h>
27#endif // _OPENMP
28
29#include "../public/VHACD.h"
30#include "btConvexHullComputer.h"
31#include "vhacdICHull.h"
32#include "vhacdMesh.h"
33#include "vhacdSArray.h"
34#include "vhacdTimer.h"
35#include "vhacdVHACD.h"
36#include "vhacdVector.h"
37#include "vhacdVolume.h"
38#include "FloatMath.h"
39
40#define MAX(a, b) (((a) > (b)) ? (a) : (b))
41#define MIN(a, b) (((a) < (b)) ? (a) : (b))
42#define ABS(a) (((a) < 0) ? -(a) : (a))
43#define ZSGN(a) (((a) < 0) ? -1 : (a) > 0 ? 1 : 0)
44#define MAX_DOUBLE (1.79769e+308)
45
46#ifdef _MSC_VER
47#pragma warning(disable:4267 4100 4244 4456)
48#endif
49
50#ifdef USE_SSE
51#include <immintrin.h>
52
53const int32_t SIMD_WIDTH = 4;
54inline int32_t FindMinimumElement(const float* const d, float* const _, const int32_t n)
55{
56 // Min within vectors
57 __m128 min_i = _mm_set1_ps(-1.0f);
58 __m128 min_v = _mm_set1_ps(std::numeric_limits<float>::max());
59 for (int32_t i = 0; i <= n - SIMD_WIDTH; i += SIMD_WIDTH) {
60 const __m128 data = _mm_load_ps(&d[i]);
61 const __m128 pred = _mm_cmplt_ps(data, min_v);
62
63 min_i = _mm_blendv_ps(min_i, _mm_set1_ps(i), pred);
64 min_v = _mm_min_ps(data, min_v);
65 }
66
67 /* Min within vector */
68 const __m128 min1 = _mm_shuffle_ps(min_v, min_v, _MM_SHUFFLE(1, 0, 3, 2));
69 const __m128 min2 = _mm_min_ps(min_v, min1);
70 const __m128 min3 = _mm_shuffle_ps(min2, min2, _MM_SHUFFLE(0, 1, 0, 1));
71 const __m128 min4 = _mm_min_ps(min2, min3);
72 float min_d = _mm_cvtss_f32(min4);
73
74 // Min index
75 const int32_t min_idx = __builtin_ctz(_mm_movemask_ps(_mm_cmpeq_ps(min_v, min4)));
76 int32_t ret = min_i[min_idx] + min_idx;
77
78 // Trailing elements
79 for (int32_t i = (n & ~(SIMD_WIDTH - 1)); i < n; ++i) {
80 if (d[i] < min_d) {
81 min_d = d[i];
82 ret = i;
83 }
84 }
85
86 *m = min_d;
87 return ret;
88}
89
90inline int32_t FindMinimumElement(const float* const d, float* const m, const int32_t begin, const int32_t end)
91{
92 // Leading elements
93 int32_t min_i = -1;
94 float min_d = std::numeric_limits<float>::max();
95 const int32_t aligned = (begin & ~(SIMD_WIDTH - 1)) + ((begin & (SIMD_WIDTH - 1)) ? SIMD_WIDTH : 0);
96 for (int32_t i = begin; i < std::min(end, aligned); ++i) {
97 if (d[i] < min_d) {
98 min_d = d[i];
99 min_i = i;
100 }
101 }
102
103 // Middle and trailing elements
104 float r_m = std::numeric_limits<float>::max();
105 const int32_t n = end - aligned;
106 const int32_t r_i = (n > 0) ? FindMinimumElement(&d[aligned], &r_m, n) : 0;
107
108 // Pick the lowest
109 if (r_m < min_d) {
110 *m = r_m;
111 return r_i + aligned;
112 }
113 else {
114 *m = min_d;
115 return min_i;
116 }
117}
118#else
119inline int32_t FindMinimumElement(const float* const d, float* const m, const int32_t begin, const int32_t end)
120{
121 int32_t idx = -1;
122 float min = (std::numeric_limits<float>::max)();
123 for (size_t i = begin; i < size_t(end); ++i) {
124 if (d[i] < min) {
125 idx = i;
126 min = d[i];
127 }
128 }
129
130 *m = min;
131 return idx;
132}
133#endif
134
135//#define OCL_SOURCE_FROM_FILE
136#ifndef OCL_SOURCE_FROM_FILE
137const char* oclProgramSource = "\
138__kernel void ComputePartialVolumes(__global short4 * voxels, \
139 const int numVoxels, \
140 const float4 plane, \
141 const float4 minBB, \
142 const float4 scale, \
143 __local uint4 * localPartialVolumes, \
144 __global uint4 * partialVolumes) \
145{ \
146 int localId = get_local_id(0); \
147 int groupSize = get_local_size(0); \
148 int i0 = get_global_id(0) << 2; \
149 float4 voxel; \
150 uint4 v; \
151 voxel = convert_float4(voxels[i0]); \
152 v.s0 = (dot(plane, mad(scale, voxel, minBB)) >= 0.0f) * (i0 < numVoxels);\
153 voxel = convert_float4(voxels[i0 + 1]); \
154 v.s1 = (dot(plane, mad(scale, voxel, minBB)) >= 0.0f) * (i0 + 1 < numVoxels);\
155 voxel = convert_float4(voxels[i0 + 2]); \
156 v.s2 = (dot(plane, mad(scale, voxel, minBB)) >= 0.0f) * (i0 + 2 < numVoxels);\
157 voxel = convert_float4(voxels[i0 + 3]); \
158 v.s3 = (dot(plane, mad(scale, voxel, minBB)) >= 0.0f) * (i0 + 3 < numVoxels);\
159 localPartialVolumes[localId] = v; \
160 barrier(CLK_LOCAL_MEM_FENCE); \
161 for (int i = groupSize >> 1; i > 0; i >>= 1) \
162 { \
163 if (localId < i) \
164 { \
165 localPartialVolumes[localId] += localPartialVolumes[localId + i]; \
166 } \
167 barrier(CLK_LOCAL_MEM_FENCE); \
168 } \
169 if (localId == 0) \
170 { \
171 partialVolumes[get_group_id(0)] = localPartialVolumes[0]; \
172 } \
173} \
174__kernel void ComputePartialSums(__global uint4 * data, \
175 const int dataSize, \
176 __local uint4 * partialSums) \
177{ \
178 int globalId = get_global_id(0); \
179 int localId = get_local_id(0); \
180 int groupSize = get_local_size(0); \
181 int i; \
182 if (globalId < dataSize) \
183 { \
184 partialSums[localId] = data[globalId]; \
185 } \
186 else \
187 { \
188 partialSums[localId] = (0, 0, 0, 0); \
189 } \
190 barrier(CLK_LOCAL_MEM_FENCE); \
191 for (i = groupSize >> 1; i > 0; i >>= 1) \
192 { \
193 if (localId < i) \
194 { \
195 partialSums[localId] += partialSums[localId + i]; \
196 } \
197 barrier(CLK_LOCAL_MEM_FENCE); \
198 } \
199 if (localId == 0) \
200 { \
201 data[get_group_id(0)] = partialSums[0]; \
202 } \
203}";
204#endif //OCL_SOURCE_FROM_FILE
205
206namespace VHACD {
207IVHACD* CreateVHACD(void)
208{
209 return new VHACD();
210}
211bool VHACD::OCLInit(void* const oclDevice, IUserLogger* const logger)
212{
213#ifdef CL_VERSION_1_1
214 m_oclDevice = (cl_device_id*)oclDevice;
215 cl_int error;
216 m_oclContext = clCreateContext(NULL, 1, m_oclDevice, NULL, NULL, &error);
217 if (error != CL_SUCCESS) {
218 if (logger) {
219 logger->Log("Couldn't create context\n");
220 }
221 return false;
222 }
223
224#ifdef OCL_SOURCE_FROM_FILE
225 std::string cl_files = OPENCL_CL_FILES;
226// read kernal from file
227#ifdef _WIN32
228 std::replace(cl_files.begin(), cl_files.end(), '/', '\\');
229#endif // _WIN32
230
231 FILE* program_handle = fopen(cl_files.c_str(), "rb");
232 fseek(program_handle, 0, SEEK_END);
233 size_t program_size = ftell(program_handle);
234 rewind(program_handle);
235 char* program_buffer = new char[program_size + 1];
236 program_buffer[program_size] = '\0';
237 fread(program_buffer, sizeof(char), program_size, program_handle);
238 fclose(program_handle);
239 // create program
240 m_oclProgram = clCreateProgramWithSource(m_oclContext, 1, (const char**)&program_buffer, &program_size, &error);
241 delete[] program_buffer;
242#else
243 size_t program_size = strlen(oclProgramSource);
244 m_oclProgram = clCreateProgramWithSource(m_oclContext, 1, (const char**)&oclProgramSource, &program_size, &error);
245#endif
246 if (error != CL_SUCCESS) {
247 if (logger) {
248 logger->Log("Couldn't create program\n");
249 }
250 return false;
251 }
252
253 /* Build program */
254 error = clBuildProgram(m_oclProgram, 1, m_oclDevice, "-cl-denorms-are-zero", NULL, NULL);
255 if (error != CL_SUCCESS) {
256 size_t log_size;
257 /* Find Size of log and print to std output */
258 clGetProgramBuildInfo(m_oclProgram, *m_oclDevice, CL_PROGRAM_BUILD_LOG, 0, NULL, &log_size);
259 char* program_log = new char[log_size + 2];
260 program_log[log_size] = '\n';
261 program_log[log_size + 1] = '\0';
262 clGetProgramBuildInfo(m_oclProgram, *m_oclDevice, CL_PROGRAM_BUILD_LOG, log_size + 1, program_log, NULL);
263 if (logger) {
264 logger->Log("Couldn't build program\n");
265 logger->Log(program_log);
266 }
267 delete[] program_log;
268 return false;
269 }
270
271 delete[] m_oclQueue;
272 delete[] m_oclKernelComputePartialVolumes;
273 delete[] m_oclKernelComputeSum;
274 m_oclQueue = new cl_command_queue[m_ompNumProcessors];
275 m_oclKernelComputePartialVolumes = new cl_kernel[m_ompNumProcessors];
276 m_oclKernelComputeSum = new cl_kernel[m_ompNumProcessors];
277
278 const char nameKernelComputePartialVolumes[] = "ComputePartialVolumes";
279 const char nameKernelComputeSum[] = "ComputePartialSums";
280 for (int32_t k = 0; k < m_ompNumProcessors; ++k) {
281 m_oclKernelComputePartialVolumes[k] = clCreateKernel(m_oclProgram, nameKernelComputePartialVolumes, &error);
282 if (error != CL_SUCCESS) {
283 if (logger) {
284 logger->Log("Couldn't create kernel\n");
285 }
286 return false;
287 }
288 m_oclKernelComputeSum[k] = clCreateKernel(m_oclProgram, nameKernelComputeSum, &error);
289 if (error != CL_SUCCESS) {
290 if (logger) {
291 logger->Log("Couldn't create kernel\n");
292 }
293 return false;
294 }
295 }
296
297 error = clGetKernelWorkGroupInfo(m_oclKernelComputePartialVolumes[0],
298 *m_oclDevice,
299 CL_KERNEL_WORK_GROUP_SIZE,
300 sizeof(size_t),
301 &m_oclWorkGroupSize,
302 NULL);
303 size_t workGroupSize = 0;
304 error = clGetKernelWorkGroupInfo(m_oclKernelComputeSum[0],
305 *m_oclDevice,
306 CL_KERNEL_WORK_GROUP_SIZE,
307 sizeof(size_t),
308 &workGroupSize,
309 NULL);
310 if (error != CL_SUCCESS) {
311 if (logger) {
312 logger->Log("Couldn't query work group info\n");
313 }
314 return false;
315 }
316
317 if (workGroupSize < m_oclWorkGroupSize) {
318 m_oclWorkGroupSize = workGroupSize;
319 }
320
321 for (int32_t k = 0; k < m_ompNumProcessors; ++k) {
322 m_oclQueue[k] = clCreateCommandQueue(m_oclContext, *m_oclDevice, 0 /*CL_QUEUE_PROFILING_ENABLE*/, &error);
323 if (error != CL_SUCCESS) {
324 if (logger) {
325 logger->Log("Couldn't create queue\n");
326 }
327 return false;
328 }
329 }
330 return true;
331#else //CL_VERSION_1_1
332 return false;
333#endif //CL_VERSION_1_1
334}
335bool VHACD::OCLRelease(IUserLogger* const logger)
336{
337#ifdef CL_VERSION_1_1
338 cl_int error;
339 if (m_oclKernelComputePartialVolumes) {
340 for (int32_t k = 0; k < m_ompNumProcessors; ++k) {
341 error = clReleaseKernel(m_oclKernelComputePartialVolumes[k]);
342 if (error != CL_SUCCESS) {
343 if (logger) {
344 logger->Log("Couldn't release kernal\n");
345 }
346 return false;
347 }
348 }
349 delete[] m_oclKernelComputePartialVolumes;
350 }
351 if (m_oclKernelComputeSum) {
352 for (int32_t k = 0; k < m_ompNumProcessors; ++k) {
353 error = clReleaseKernel(m_oclKernelComputeSum[k]);
354 if (error != CL_SUCCESS) {
355 if (logger) {
356 logger->Log("Couldn't release kernal\n");
357 }
358 return false;
359 }
360 }
361 delete[] m_oclKernelComputeSum;
362 }
363 if (m_oclQueue) {
364 for (int32_t k = 0; k < m_ompNumProcessors; ++k) {
365 error = clReleaseCommandQueue(m_oclQueue[k]);
366 if (error != CL_SUCCESS) {
367 if (logger) {
368 logger->Log("Couldn't release queue\n");
369 }
370 return false;
371 }
372 }
373 delete[] m_oclQueue;
374 }
375 error = clReleaseProgram(m_oclProgram);
376 if (error != CL_SUCCESS) {
377 if (logger) {
378 logger->Log("Couldn't release program\n");
379 }
380 return false;
381 }
382 error = clReleaseContext(m_oclContext);
383 if (error != CL_SUCCESS) {
384 if (logger) {
385 logger->Log("Couldn't release context\n");
386 }
387 return false;
388 }
389
390 return true;
391#else //CL_VERSION_1_1
392 return false;
393#endif //CL_VERSION_1_1
394}
395void VHACD::ComputePrimitiveSet(const Parameters& params)
396{
397 if (GetCancel()) {
398 return;
399 }
400 m_timer.Tic();
401
402 m_stage = "Compute primitive set";
403 m_operation = "Convert volume to pset";
404
405 std::ostringstream msg;
406 if (params.m_logger) {
407 msg << "+ " << m_stage << std::endl;
408 params.m_logger->Log(msg.str().c_str());
409 }
410
411 Update(0.0, 0.0, params);
412 if (params.m_mode == 0) {
413 VoxelSet* vset = new VoxelSet;
414 m_volume->Convert(*vset);
415 m_pset = vset;
416 }
417 else {
418 TetrahedronSet* tset = new TetrahedronSet;
419 m_volume->Convert(*tset);
420 m_pset = tset;
421 }
422
423 delete m_volume;
424 m_volume = 0;
425
426 if (params.m_logger) {
427 msg.str("");
428 msg << "\t # primitives " << m_pset->GetNPrimitives() << std::endl;
429 msg << "\t # inside surface " << m_pset->GetNPrimitivesInsideSurf() << std::endl;
430 msg << "\t # on surface " << m_pset->GetNPrimitivesOnSurf() << std::endl;
431 params.m_logger->Log(msg.str().c_str());
432 }
433
434 m_overallProgress = 15.0;
435 Update(100.0, 100.0, params);
436 m_timer.Toc();
437 if (params.m_logger) {
438 msg.str("");
439 msg << "\t time " << m_timer.GetElapsedTime() / 1000.0 << "s" << std::endl;
440 params.m_logger->Log(msg.str().c_str());
441 }
442}
443bool VHACD::Compute(const double* const points, const uint32_t nPoints,
444 const uint32_t* const triangles,const uint32_t nTriangles, const Parameters& params)
445{
446 return ComputeACD(points, nPoints, triangles, nTriangles, params);
447}
448bool VHACD::Compute(const float* const points,const uint32_t nPoints,
449 const uint32_t* const triangles,const uint32_t nTriangles, const Parameters& params)
450{
451 return ComputeACD(points, nPoints, triangles, nTriangles, params);
452}
453double ComputePreferredCuttingDirection(const PrimitiveSet* const tset, Vec3<double>& dir)
454{
455 double ex = tset->GetEigenValue(AXIS_X);
456 double ey = tset->GetEigenValue(AXIS_Y);
457 double ez = tset->GetEigenValue(AXIS_Z);
458 double vx = (ey - ez) * (ey - ez);
459 double vy = (ex - ez) * (ex - ez);
460 double vz = (ex - ey) * (ex - ey);
461 if (vx < vy && vx < vz) {
462 double e = ey * ey + ez * ez;
463 dir[0] = 1.0;
464 dir[1] = 0.0;
465 dir[2] = 0.0;
466 return (e == 0.0) ? 0.0 : 1.0 - vx / e;
467 }
468 else if (vy < vx && vy < vz) {
469 double e = ex * ex + ez * ez;
470 dir[0] = 0.0;
471 dir[1] = 1.0;
472 dir[2] = 0.0;
473 return (e == 0.0) ? 0.0 : 1.0 - vy / e;
474 }
475 else {
476 double e = ex * ex + ey * ey;
477 dir[0] = 0.0;
478 dir[1] = 0.0;
479 dir[2] = 1.0;
480 return (e == 0.0) ? 0.0 : 1.0 - vz / e;
481 }
482}
483void ComputeAxesAlignedClippingPlanes(const VoxelSet& vset, const short downsampling, SArray<Plane>& planes)
484{
485 const Vec3<short> minV = vset.GetMinBBVoxels();
486 const Vec3<short> maxV = vset.GetMaxBBVoxels();
487 Vec3<double> pt;
488 Plane plane;
489 const short i0 = minV[0];
490 const short i1 = maxV[0];
491 plane.m_a = 1.0;
492 plane.m_b = 0.0;
493 plane.m_c = 0.0;
494 plane.m_axis = AXIS_X;
495 for (short i = i0; i <= i1; i += downsampling) {
496 pt = vset.GetPoint(Vec3<double>(i + 0.5, 0.0, 0.0));
497 plane.m_d = -pt[0];
498 plane.m_index = i;
499 planes.PushBack(plane);
500 }
501 const short j0 = minV[1];
502 const short j1 = maxV[1];
503 plane.m_a = 0.0;
504 plane.m_b = 1.0;
505 plane.m_c = 0.0;
506 plane.m_axis = AXIS_Y;
507 for (short j = j0; j <= j1; j += downsampling) {
508 pt = vset.GetPoint(Vec3<double>(0.0, j + 0.5, 0.0));
509 plane.m_d = -pt[1];
510 plane.m_index = j;
511 planes.PushBack(plane);
512 }
513 const short k0 = minV[2];
514 const short k1 = maxV[2];
515 plane.m_a = 0.0;
516 plane.m_b = 0.0;
517 plane.m_c = 1.0;
518 plane.m_axis = AXIS_Z;
519 for (short k = k0; k <= k1; k += downsampling) {
520 pt = vset.GetPoint(Vec3<double>(0.0, 0.0, k + 0.5));
521 plane.m_d = -pt[2];
522 plane.m_index = k;
523 planes.PushBack(plane);
524 }
525}
526void ComputeAxesAlignedClippingPlanes(const TetrahedronSet& tset, const short downsampling, SArray<Plane>& planes)
527{
528 const Vec3<double> minV = tset.GetMinBB();
529 const Vec3<double> maxV = tset.GetMaxBB();
530 const double scale = tset.GetSacle();
531 const short i0 = 0;
532 const short j0 = 0;
533 const short k0 = 0;
534 const short i1 = static_cast<short>((maxV[0] - minV[0]) / scale + 0.5);
535 const short j1 = static_cast<short>((maxV[1] - minV[1]) / scale + 0.5);
536 const short k1 = static_cast<short>((maxV[2] - minV[2]) / scale + 0.5);
537
538 Plane plane;
539 plane.m_a = 1.0;
540 plane.m_b = 0.0;
541 plane.m_c = 0.0;
542 plane.m_axis = AXIS_X;
543 for (short i = i0; i <= i1; i += downsampling) {
544 double x = minV[0] + scale * i;
545 plane.m_d = -x;
546 plane.m_index = i;
547 planes.PushBack(plane);
548 }
549 plane.m_a = 0.0;
550 plane.m_b = 1.0;
551 plane.m_c = 0.0;
552 plane.m_axis = AXIS_Y;
553 for (short j = j0; j <= j1; j += downsampling) {
554 double y = minV[1] + scale * j;
555 plane.m_d = -y;
556 plane.m_index = j;
557 planes.PushBack(plane);
558 }
559 plane.m_a = 0.0;
560 plane.m_b = 0.0;
561 plane.m_c = 1.0;
562 plane.m_axis = AXIS_Z;
563 for (short k = k0; k <= k1; k += downsampling) {
564 double z = minV[2] + scale * k;
565 plane.m_d = -z;
566 plane.m_index = k;
567 planes.PushBack(plane);
568 }
569}
570void RefineAxesAlignedClippingPlanes(const VoxelSet& vset, const Plane& bestPlane, const short downsampling,
571 SArray<Plane>& planes)
572{
573 const Vec3<short> minV = vset.GetMinBBVoxels();
574 const Vec3<short> maxV = vset.GetMaxBBVoxels();
575 Vec3<double> pt;
576 Plane plane;
577
578 if (bestPlane.m_axis == AXIS_X) {
579 const short i0 = MAX(minV[0], bestPlane.m_index - downsampling);
580 const short i1 = MIN(maxV[0], bestPlane.m_index + downsampling);
581 plane.m_a = 1.0;
582 plane.m_b = 0.0;
583 plane.m_c = 0.0;
584 plane.m_axis = AXIS_X;
585 for (short i = i0; i <= i1; ++i) {
586 pt = vset.GetPoint(Vec3<double>(i + 0.5, 0.0, 0.0));
587 plane.m_d = -pt[0];
588 plane.m_index = i;
589 planes.PushBack(plane);
590 }
591 }
592 else if (bestPlane.m_axis == AXIS_Y) {
593 const short j0 = MAX(minV[1], bestPlane.m_index - downsampling);
594 const short j1 = MIN(maxV[1], bestPlane.m_index + downsampling);
595 plane.m_a = 0.0;
596 plane.m_b = 1.0;
597 plane.m_c = 0.0;
598 plane.m_axis = AXIS_Y;
599 for (short j = j0; j <= j1; ++j) {
600 pt = vset.GetPoint(Vec3<double>(0.0, j + 0.5, 0.0));
601 plane.m_d = -pt[1];
602 plane.m_index = j;
603 planes.PushBack(plane);
604 }
605 }
606 else {
607 const short k0 = MAX(minV[2], bestPlane.m_index - downsampling);
608 const short k1 = MIN(maxV[2], bestPlane.m_index + downsampling);
609 plane.m_a = 0.0;
610 plane.m_b = 0.0;
611 plane.m_c = 1.0;
612 plane.m_axis = AXIS_Z;
613 for (short k = k0; k <= k1; ++k) {
614 pt = vset.GetPoint(Vec3<double>(0.0, 0.0, k + 0.5));
615 plane.m_d = -pt[2];
616 plane.m_index = k;
617 planes.PushBack(plane);
618 }
619 }
620}
621void RefineAxesAlignedClippingPlanes(const TetrahedronSet& tset, const Plane& bestPlane, const short downsampling,
622 SArray<Plane>& planes)
623{
624 const Vec3<double> minV = tset.GetMinBB();
625 const Vec3<double> maxV = tset.GetMaxBB();
626 const double scale = tset.GetSacle();
627 Plane plane;
628
629 if (bestPlane.m_axis == AXIS_X) {
630 const short i0 = MAX(0, bestPlane.m_index - downsampling);
631 const short i1 = static_cast<short>(MIN((maxV[0] - minV[0]) / scale + 0.5, bestPlane.m_index + downsampling));
632 plane.m_a = 1.0;
633 plane.m_b = 0.0;
634 plane.m_c = 0.0;
635 plane.m_axis = AXIS_X;
636 for (short i = i0; i <= i1; ++i) {
637 double x = minV[0] + scale * i;
638 plane.m_d = -x;
639 plane.m_index = i;
640 planes.PushBack(plane);
641 }
642 }
643 else if (bestPlane.m_axis == AXIS_Y) {
644 const short j0 = MAX(0, bestPlane.m_index - downsampling);
645 const short j1 = static_cast<short>(MIN((maxV[1] - minV[1]) / scale + 0.5, bestPlane.m_index + downsampling));
646 plane.m_a = 0.0;
647 plane.m_b = 1.0;
648 plane.m_c = 0.0;
649 plane.m_axis = AXIS_Y;
650 for (short j = j0; j <= j1; ++j) {
651 double y = minV[1] + scale * j;
652 plane.m_d = -y;
653 plane.m_index = j;
654 planes.PushBack(plane);
655 }
656 }
657 else {
658 const short k0 = MAX(0, bestPlane.m_index - downsampling);
659 const short k1 = static_cast<short>(MIN((maxV[2] - minV[2]) / scale + 0.5, bestPlane.m_index + downsampling));
660 plane.m_a = 0.0;
661 plane.m_b = 0.0;
662 plane.m_c = 1.0;
663 plane.m_axis = AXIS_Z;
664 for (short k = k0; k <= k1; ++k) {
665 double z = minV[2] + scale * k;
666 plane.m_d = -z;
667 plane.m_index = k;
668 planes.PushBack(plane);
669 }
670 }
671}
672inline double ComputeLocalConcavity(const double volume, const double volumeCH)
673{
674 return fabs(volumeCH - volume) / volumeCH;
675}
676inline double ComputeConcavity(const double volume, const double volumeCH, const double volume0)
677{
678 return fabs(volumeCH - volume) / volume0;
679}
680
681//#define DEBUG_TEMP
682void VHACD::ComputeBestClippingPlane(const PrimitiveSet* inputPSet, const double volume, const SArray<Plane>& planes,
683 const Vec3<double>& preferredCuttingDirection, const double w, const double alpha, const double beta,
684 const int32_t convexhullDownsampling, const double progress0, const double progress1, Plane& bestPlane,
685 double& minConcavity, const Parameters& params)
686{
687 if (GetCancel()) {
688 return;
689 }
690 char msg[256];
691 size_t nPrimitives = inputPSet->GetNPrimitives();
692 bool oclAcceleration = (nPrimitives > OCL_MIN_NUM_PRIMITIVES && params.m_oclAcceleration && params.m_mode == 0) ? true : false;
693 int32_t iBest = -1;
694 int32_t nPlanes = static_cast<int32_t>(planes.Size());
695 bool cancel = false;
696 int32_t done = 0;
697 double minTotal = MAX_DOUBLE;
698 double minBalance = MAX_DOUBLE;
699 double minSymmetry = MAX_DOUBLE;
700 minConcavity = MAX_DOUBLE;
701
702 SArray<Vec3<double> >* chPts = new SArray<Vec3<double> >[2 * m_ompNumProcessors];
703 Mesh* chs = new Mesh[2 * m_ompNumProcessors];
704 PrimitiveSet* onSurfacePSet = inputPSet->Create();
705 inputPSet->SelectOnSurface(onSurfacePSet);
706
707 PrimitiveSet** psets = 0;
708 if (!params.m_convexhullApproximation) {
709 psets = new PrimitiveSet*[2 * m_ompNumProcessors];
710 for (int32_t i = 0; i < 2 * m_ompNumProcessors; ++i) {
711 psets[i] = inputPSet->Create();
712 }
713 }
714
715#ifdef CL_VERSION_1_1
716 // allocate OpenCL data structures
717 cl_mem voxels;
718 cl_mem* partialVolumes = 0;
719 size_t globalSize = 0;
720 size_t nWorkGroups = 0;
721 double unitVolume = 0.0;
722 if (oclAcceleration) {
723 VoxelSet* vset = (VoxelSet*)inputPSet;
724 const Vec3<double> minBB = vset->GetMinBB();
725 const float fMinBB[4] = { (float)minBB[0], (float)minBB[1], (float)minBB[2], 1.0f };
726 const float fSclae[4] = { (float)vset->GetScale(), (float)vset->GetScale(), (float)vset->GetScale(), 0.0f };
727 const int32_t nVoxels = (int32_t)nPrimitives;
728 unitVolume = vset->GetUnitVolume();
729 nWorkGroups = (nPrimitives + 4 * m_oclWorkGroupSize - 1) / (4 * m_oclWorkGroupSize);
730 globalSize = nWorkGroups * m_oclWorkGroupSize;
731 cl_int error;
732 voxels = clCreateBuffer(m_oclContext,
733 CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
734 sizeof(Voxel) * nPrimitives,
735 vset->GetVoxels(),
736 &error);
737 if (error != CL_SUCCESS) {
738 if (params.m_logger) {
739 params.m_logger->Log("Couldn't create buffer\n");
740 }
741 SetCancel(true);
742 }
743
744 partialVolumes = new cl_mem[m_ompNumProcessors];
745 for (int32_t i = 0; i < m_ompNumProcessors; ++i) {
746 partialVolumes[i] = clCreateBuffer(m_oclContext,
747 CL_MEM_WRITE_ONLY,
748 sizeof(uint32_t) * 4 * nWorkGroups,
749 NULL,
750 &error);
751 if (error != CL_SUCCESS) {
752 if (params.m_logger) {
753 params.m_logger->Log("Couldn't create buffer\n");
754 }
755 SetCancel(true);
756 break;
757 }
758 error = clSetKernelArg(m_oclKernelComputePartialVolumes[i], 0, sizeof(cl_mem), &voxels);
759 error |= clSetKernelArg(m_oclKernelComputePartialVolumes[i], 1, sizeof(uint32_t), &nVoxels);
760 error |= clSetKernelArg(m_oclKernelComputePartialVolumes[i], 3, sizeof(float) * 4, fMinBB);
761 error |= clSetKernelArg(m_oclKernelComputePartialVolumes[i], 4, sizeof(float) * 4, &fSclae);
762 error |= clSetKernelArg(m_oclKernelComputePartialVolumes[i], 5, sizeof(uint32_t) * 4 * m_oclWorkGroupSize, NULL);
763 error |= clSetKernelArg(m_oclKernelComputePartialVolumes[i], 6, sizeof(cl_mem), &(partialVolumes[i]));
764 error |= clSetKernelArg(m_oclKernelComputeSum[i], 0, sizeof(cl_mem), &(partialVolumes[i]));
765 error |= clSetKernelArg(m_oclKernelComputeSum[i], 2, sizeof(uint32_t) * 4 * m_oclWorkGroupSize, NULL);
766 if (error != CL_SUCCESS) {
767 if (params.m_logger) {
768 params.m_logger->Log("Couldn't kernel arguments \n");
769 }
770 SetCancel(true);
771 }
772 }
773 }
774#else // CL_VERSION_1_1
775 oclAcceleration = false;
776#endif // CL_VERSION_1_1
777
778#ifdef DEBUG_TEMP
779 Timer timerComputeCost;
780 timerComputeCost.Tic();
781#endif // DEBUG_TEMP
782
783#if USE_THREAD == 1 && _OPENMP
784#pragma omp parallel for
785#endif
786 for (int32_t x = 0; x < nPlanes; ++x) {
787 int32_t threadID = 0;
788#if USE_THREAD == 1 && _OPENMP
789 threadID = omp_get_thread_num();
790#pragma omp flush(cancel)
791#endif
792 if (!cancel) {
793 //Update progress
794 if (GetCancel()) {
795 cancel = true;
796#if USE_THREAD == 1 && _OPENMP
797#pragma omp flush(cancel)
798#endif
799 }
800 Plane plane = planes[x];
801
802 if (oclAcceleration) {
803#ifdef CL_VERSION_1_1
804 const float fPlane[4] = { (float)plane.m_a, (float)plane.m_b, (float)plane.m_c, (float)plane.m_d };
805 cl_int error = clSetKernelArg(m_oclKernelComputePartialVolumes[threadID], 2, sizeof(float) * 4, fPlane);
806 if (error != CL_SUCCESS) {
807 if (params.m_logger) {
808 params.m_logger->Log("Couldn't kernel atguments \n");
809 }
810 SetCancel(true);
811 }
812
813 error = clEnqueueNDRangeKernel(m_oclQueue[threadID], m_oclKernelComputePartialVolumes[threadID],
814 1, NULL, &globalSize, &m_oclWorkGroupSize, 0, NULL, NULL);
815 if (error != CL_SUCCESS) {
816 if (params.m_logger) {
817 params.m_logger->Log("Couldn't run kernel \n");
818 }
819 SetCancel(true);
820 }
821 int32_t nValues = (int32_t)nWorkGroups;
822 while (nValues > 1) {
823 error = clSetKernelArg(m_oclKernelComputeSum[threadID], 1, sizeof(int32_t), &nValues);
824 if (error != CL_SUCCESS) {
825 if (params.m_logger) {
826 params.m_logger->Log("Couldn't kernel atguments \n");
827 }
828 SetCancel(true);
829 }
830 size_t nWorkGroups = (nValues + m_oclWorkGroupSize - 1) / m_oclWorkGroupSize;
831 size_t globalSize = nWorkGroups * m_oclWorkGroupSize;
832 error = clEnqueueNDRangeKernel(m_oclQueue[threadID], m_oclKernelComputeSum[threadID],
833 1, NULL, &globalSize, &m_oclWorkGroupSize, 0, NULL, NULL);
834 if (error != CL_SUCCESS) {
835 if (params.m_logger) {
836 params.m_logger->Log("Couldn't run kernel \n");
837 }
838 SetCancel(true);
839 }
840 nValues = (int32_t)nWorkGroups;
841 }
842#endif // CL_VERSION_1_1
843 }
844
845 Mesh& leftCH = chs[threadID];
846 Mesh& rightCH = chs[threadID + m_ompNumProcessors];
847 rightCH.ResizePoints(0);
848 leftCH.ResizePoints(0);
849 rightCH.ResizeTriangles(0);
850 leftCH.ResizeTriangles(0);
851
852// compute convex-hulls
853#ifdef TEST_APPROX_CH
854 double volumeLeftCH1;
855 double volumeRightCH1;
856#endif //TEST_APPROX_CH
857 if (params.m_convexhullApproximation) {
858 SArray<Vec3<double> >& leftCHPts = chPts[threadID];
859 SArray<Vec3<double> >& rightCHPts = chPts[threadID + m_ompNumProcessors];
860 rightCHPts.Resize(0);
861 leftCHPts.Resize(0);
862 onSurfacePSet->Intersect(plane, &rightCHPts, &leftCHPts, convexhullDownsampling * 32);
863 inputPSet->GetConvexHull().Clip(plane, rightCHPts, leftCHPts);
864 rightCH.ComputeConvexHull((double*)rightCHPts.Data(), rightCHPts.Size());
865 leftCH.ComputeConvexHull((double*)leftCHPts.Data(), leftCHPts.Size());
866#ifdef TEST_APPROX_CH
867 Mesh leftCH1;
868 Mesh rightCH1;
869 VoxelSet right;
870 VoxelSet left;
871 onSurfacePSet->Clip(plane, &right, &left);
872 right.ComputeConvexHull(rightCH1, convexhullDownsampling);
873 left.ComputeConvexHull(leftCH1, convexhullDownsampling);
874
875 volumeLeftCH1 = leftCH1.ComputeVolume();
876 volumeRightCH1 = rightCH1.ComputeVolume();
877#endif //TEST_APPROX_CH
878 }
879 else {
880 PrimitiveSet* const right = psets[threadID];
881 PrimitiveSet* const left = psets[threadID + m_ompNumProcessors];
882 onSurfacePSet->Clip(plane, right, left);
883 right->ComputeConvexHull(rightCH, convexhullDownsampling);
884 left->ComputeConvexHull(leftCH, convexhullDownsampling);
885 }
886 double volumeLeftCH = leftCH.ComputeVolume();
887 double volumeRightCH = rightCH.ComputeVolume();
888
889 // compute clipped volumes
890 double volumeLeft = 0.0;
891 double volumeRight = 0.0;
892 if (oclAcceleration) {
893#ifdef CL_VERSION_1_1
894 uint32_t volumes[4];
895 cl_int error = clEnqueueReadBuffer(m_oclQueue[threadID], partialVolumes[threadID], CL_TRUE,
896 0, sizeof(uint32_t) * 4, volumes, 0, NULL, NULL);
897 size_t nPrimitivesRight = volumes[0] + volumes[1] + volumes[2] + volumes[3];
898 size_t nPrimitivesLeft = nPrimitives - nPrimitivesRight;
899 volumeRight = nPrimitivesRight * unitVolume;
900 volumeLeft = nPrimitivesLeft * unitVolume;
901 if (error != CL_SUCCESS) {
902 if (params.m_logger) {
903 params.m_logger->Log("Couldn't read buffer \n");
904 }
905 SetCancel(true);
906 }
907#endif // CL_VERSION_1_1
908 }
909 else {
910 inputPSet->ComputeClippedVolumes(plane, volumeRight, volumeLeft);
911 }
912 double concavityLeft = ComputeConcavity(volumeLeft, volumeLeftCH, m_volumeCH0);
913 double concavityRight = ComputeConcavity(volumeRight, volumeRightCH, m_volumeCH0);
914 double concavity = (concavityLeft + concavityRight);
915
916 // compute cost
917 double balance = alpha * fabs(volumeLeft - volumeRight) / m_volumeCH0;
918 double d = w * (preferredCuttingDirection[0] * plane.m_a + preferredCuttingDirection[1] * plane.m_b + preferredCuttingDirection[2] * plane.m_c);
919 double symmetry = beta * d;
920 double total = concavity + balance + symmetry;
921
922#if USE_THREAD == 1 && _OPENMP
923#pragma omp critical
924#endif
925 {
926 if (total < minTotal || (total == minTotal && x < iBest)) {
927 minConcavity = concavity;
928 minBalance = balance;
929 minSymmetry = symmetry;
930 bestPlane = plane;
931 minTotal = total;
932 iBest = x;
933 }
934 ++done;
935 if (!(done & 127)) // reduce update frequency
936 {
937 double progress = done * (progress1 - progress0) / nPlanes + progress0;
938 Update(m_stageProgress, progress, params);
939 }
940 }
941 }
942 }
943
944#ifdef DEBUG_TEMP
945 timerComputeCost.Toc();
946 printf_s("Cost[%i] = %f\n", nPlanes, timerComputeCost.GetElapsedTime());
947#endif // DEBUG_TEMP
948
949#ifdef CL_VERSION_1_1
950 if (oclAcceleration) {
951 clReleaseMemObject(voxels);
952 for (int32_t i = 0; i < m_ompNumProcessors; ++i) {
953 clReleaseMemObject(partialVolumes[i]);
954 }
955 delete[] partialVolumes;
956 }
957#endif // CL_VERSION_1_1
958
959 if (psets) {
960 for (int32_t i = 0; i < 2 * m_ompNumProcessors; ++i) {
961 delete psets[i];
962 }
963 delete[] psets;
964 }
965 delete onSurfacePSet;
966 delete[] chPts;
967 delete[] chs;
968 if (params.m_logger) {
969 sprintf(msg, "\n\t\t\t Best %04i T=%2.6f C=%2.6f B=%2.6f S=%2.6f (%1.1f, %1.1f, %1.1f, %3.3f)\n\n", iBest, minTotal, minConcavity, minBalance, minSymmetry, bestPlane.m_a, bestPlane.m_b, bestPlane.m_c, bestPlane.m_d);
970 params.m_logger->Log(msg);
971 }
972}
973void VHACD::ComputeACD(const Parameters& params)
974{
975 if (GetCancel()) {
976 return;
977 }
978 m_timer.Tic();
979
980 m_stage = "Approximate Convex Decomposition";
981 m_stageProgress = 0.0;
982 std::ostringstream msg;
983 if (params.m_logger) {
984 msg << "+ " << m_stage << std::endl;
985 params.m_logger->Log(msg.str().c_str());
986 }
987
988 SArray<PrimitiveSet*> parts;
989 SArray<PrimitiveSet*> inputParts;
990 SArray<PrimitiveSet*> temp;
991 inputParts.PushBack(m_pset);
992 m_pset = 0;
993 SArray<Plane> planes;
994 SArray<Plane> planesRef;
995 uint32_t sub = 0;
996 bool firstIteration = true;
997 m_volumeCH0 = 1.0;
998
999 // Compute the decomposition depth based on the number of convex hulls being requested..
1000 uint32_t hullCount = 2;
1001 uint32_t depth = 1;
1002 while (params.m_maxConvexHulls > hullCount)
1003 {
1004 depth++;
1005 hullCount *= 2;
1006 }
1007 // We must always increment the decomposition depth one higher than the maximum number of hulls requested.
1008 // The reason for this is as follows.
1009 // Say, for example, the user requests 32 convex hulls exactly. This would be a decomposition depth of 5.
1010 // However, when we do that, we do *not* necessarily get 32 hulls as a result. This is because, during
1011 // the recursive descent of the binary tree, one or more of the leaf nodes may have no concavity and
1012 // will not be split. So, in this way, even with a decomposition depth of 5, you can produce fewer than
1013 // 32 hulls. So, in this case, we would set the decomposition depth to 6 (producing up to as high as 64 convex hulls).
1014 // Then, the merge step which combines over-described hulls down to the user requested amount, we will end up
1015 // getting exactly 32 convex hulls as a result.
1016 // We could just allow the artist to directly control the decomposition depth directly, but this would be a bit
1017 // too complex and the preference is simply to let them specify how many hulls they want and derive the solution
1018 // from that.
1019 depth++;
1020
1021
1022 while (sub++ < depth && inputParts.Size() > 0 && !m_cancel) {
1023 msg.str("");
1024 msg << "Subdivision level " << sub;
1025 m_operation = msg.str();
1026
1027 if (params.m_logger) {
1028 msg.str("");
1029 msg << "\t Subdivision level " << sub << std::endl;
1030 params.m_logger->Log(msg.str().c_str());
1031 }
1032
1033 double maxConcavity = 0.0;
1034 const size_t nInputParts = inputParts.Size();
1035 Update(m_stageProgress, 0.0, params);
1036 for (size_t p = 0; p < nInputParts && !m_cancel; ++p) {
1037 const double progress0 = p * 100.0 / nInputParts;
1038 const double progress1 = (p + 0.75) * 100.0 / nInputParts;
1039 const double progress2 = (p + 1.00) * 100.0 / nInputParts;
1040
1041 Update(m_stageProgress, progress0, params);
1042
1043 PrimitiveSet* pset = inputParts[p];
1044 inputParts[p] = 0;
1045 double volume = pset->ComputeVolume();
1046 pset->ComputeBB();
1047 pset->ComputePrincipalAxes();
1048 if (params.m_pca) {
1049 pset->AlignToPrincipalAxes();
1050 }
1051
1052 pset->ComputeConvexHull(pset->GetConvexHull());
1053 double volumeCH = fabs(pset->GetConvexHull().ComputeVolume());
1054 if (firstIteration) {
1055 m_volumeCH0 = volumeCH;
1056 }
1057
1058 double concavity = ComputeConcavity(volume, volumeCH, m_volumeCH0);
1059 double error = 1.01 * pset->ComputeMaxVolumeError() / m_volumeCH0;
1060
1061 if (firstIteration) {
1062 firstIteration = false;
1063 }
1064
1065 if (params.m_logger) {
1066 msg.str("");
1067 msg << "\t -> Part[" << p
1068 << "] C = " << concavity
1069 << ", E = " << error
1070 << ", VS = " << pset->GetNPrimitivesOnSurf()
1071 << ", VI = " << pset->GetNPrimitivesInsideSurf()
1072 << std::endl;
1073 params.m_logger->Log(msg.str().c_str());
1074 }
1075
1076 if (concavity > params.m_concavity && concavity > error) {
1077 Vec3<double> preferredCuttingDirection;
1078 double w = ComputePreferredCuttingDirection(pset, preferredCuttingDirection);
1079 planes.Resize(0);
1080 if (params.m_mode == 0) {
1081 VoxelSet* vset = (VoxelSet*)pset;
1082 ComputeAxesAlignedClippingPlanes(*vset, params.m_planeDownsampling, planes);
1083 }
1084 else {
1085 TetrahedronSet* tset = (TetrahedronSet*)pset;
1086 ComputeAxesAlignedClippingPlanes(*tset, params.m_planeDownsampling, planes);
1087 }
1088
1089 if (params.m_logger) {
1090 msg.str("");
1091 msg << "\t\t [Regular sampling] Number of clipping planes " << planes.Size() << std::endl;
1092 params.m_logger->Log(msg.str().c_str());
1093 }
1094
1095 Plane bestPlane;
1096 double minConcavity = MAX_DOUBLE;
1097 ComputeBestClippingPlane(pset,
1098 volume,
1099 planes,
1100 preferredCuttingDirection,
1101 w,
1102 concavity * params.m_alpha,
1103 concavity * params.m_beta,
1104 params.m_convexhullDownsampling,
1105 progress0,
1106 progress1,
1107 bestPlane,
1108 minConcavity,
1109 params);
1110 if (!m_cancel && (params.m_planeDownsampling > 1 || params.m_convexhullDownsampling > 1)) {
1111 planesRef.Resize(0);
1112
1113 if (params.m_mode == 0) {
1114 VoxelSet* vset = (VoxelSet*)pset;
1115 RefineAxesAlignedClippingPlanes(*vset, bestPlane, params.m_planeDownsampling, planesRef);
1116 }
1117 else {
1118 TetrahedronSet* tset = (TetrahedronSet*)pset;
1119 RefineAxesAlignedClippingPlanes(*tset, bestPlane, params.m_planeDownsampling, planesRef);
1120 }
1121
1122 if (params.m_logger) {
1123 msg.str("");
1124 msg << "\t\t [Refining] Number of clipping planes " << planesRef.Size() << std::endl;
1125 params.m_logger->Log(msg.str().c_str());
1126 }
1127 ComputeBestClippingPlane(pset,
1128 volume,
1129 planesRef,
1130 preferredCuttingDirection,
1131 w,
1132 concavity * params.m_alpha,
1133 concavity * params.m_beta,
1134 1, // convexhullDownsampling = 1
1135 progress1,
1136 progress2,
1137 bestPlane,
1138 minConcavity,
1139 params);
1140 }
1141 if (GetCancel()) {
1142 delete pset; // clean up
1143 break;
1144 }
1145 else {
1146 if (maxConcavity < minConcavity) {
1147 maxConcavity = minConcavity;
1148 }
1149 PrimitiveSet* bestLeft = pset->Create();
1150 PrimitiveSet* bestRight = pset->Create();
1151 temp.PushBack(bestLeft);
1152 temp.PushBack(bestRight);
1153 pset->Clip(bestPlane, bestRight, bestLeft);
1154 if (params.m_pca) {
1155 bestRight->RevertAlignToPrincipalAxes();
1156 bestLeft->RevertAlignToPrincipalAxes();
1157 }
1158 delete pset;
1159 }
1160 }
1161 else {
1162 if (params.m_pca) {
1163 pset->RevertAlignToPrincipalAxes();
1164 }
1165 parts.PushBack(pset);
1166 }
1167 }
1168
1169 Update(95.0 * (1.0 - maxConcavity) / (1.0 - params.m_concavity), 100.0, params);
1170 if (GetCancel()) {
1171 const size_t nTempParts = temp.Size();
1172 for (size_t p = 0; p < nTempParts; ++p) {
1173 delete temp[p];
1174 }
1175 temp.Resize(0);
1176 }
1177 else {
1178 inputParts = temp;
1179 temp.Resize(0);
1180 }
1181 }
1182 const size_t nInputParts = inputParts.Size();
1183 for (size_t p = 0; p < nInputParts; ++p) {
1184 parts.PushBack(inputParts[p]);
1185 }
1186
1187 if (GetCancel()) {
1188 const size_t nParts = parts.Size();
1189 for (size_t p = 0; p < nParts; ++p) {
1190 delete parts[p];
1191 }
1192 return;
1193 }
1194
1195 m_overallProgress = 90.0;
1196 Update(m_stageProgress, 100.0, params);
1197
1198 msg.str("");
1199 msg << "Generate convex-hulls";
1200 m_operation = msg.str();
1201 size_t nConvexHulls = parts.Size();
1202 if (params.m_logger) {
1203 msg.str("");
1204 msg << "+ Generate " << nConvexHulls << " convex-hulls " << std::endl;
1205 params.m_logger->Log(msg.str().c_str());
1206 }
1207
1208 Update(m_stageProgress, 0.0, params);
1209 m_convexHulls.Resize(0);
1210 for (size_t p = 0; p < nConvexHulls && !m_cancel; ++p) {
1211 Update(m_stageProgress, p * 100.0 / nConvexHulls, params);
1212 m_convexHulls.PushBack(new Mesh);
1213 parts[p]->ComputeConvexHull(*m_convexHulls[p]);
1214 size_t nv = m_convexHulls[p]->GetNPoints();
1215 double x, y, z;
1216 for (size_t i = 0; i < nv; ++i) {
1217 Vec3<double>& pt = m_convexHulls[p]->GetPoint(i);
1218 x = pt[0];
1219 y = pt[1];
1220 z = pt[2];
1221 pt[0] = m_rot[0][0] * x + m_rot[0][1] * y + m_rot[0][2] * z + m_barycenter[0];
1222 pt[1] = m_rot[1][0] * x + m_rot[1][1] * y + m_rot[1][2] * z + m_barycenter[1];
1223 pt[2] = m_rot[2][0] * x + m_rot[2][1] * y + m_rot[2][2] * z + m_barycenter[2];
1224 }
1225 }
1226
1227 const size_t nParts = parts.Size();
1228 for (size_t p = 0; p < nParts; ++p) {
1229 delete parts[p];
1230 parts[p] = 0;
1231 }
1232 parts.Resize(0);
1233
1234 if (GetCancel()) {
1235 const size_t nConvexHulls = m_convexHulls.Size();
1236 for (size_t p = 0; p < nConvexHulls; ++p) {
1237 delete m_convexHulls[p];
1238 }
1239 m_convexHulls.Clear();
1240 return;
1241 }
1242
1243 m_overallProgress = 95.0;
1244 Update(100.0, 100.0, params);
1245 m_timer.Toc();
1246 if (params.m_logger) {
1247 msg.str("");
1248 msg << "\t time " << m_timer.GetElapsedTime() / 1000.0 << "s" << std::endl;
1249 params.m_logger->Log(msg.str().c_str());
1250 }
1251}
1252void AddPoints(const Mesh* const mesh, SArray<Vec3<double> >& pts)
1253{
1254 const int32_t n = (int32_t)mesh->GetNPoints();
1255 for (int32_t i = 0; i < n; ++i) {
1256 pts.PushBack(mesh->GetPoint(i));
1257 }
1258}
1259void ComputeConvexHull(const Mesh* const ch1, const Mesh* const ch2, SArray<Vec3<double> >& pts, Mesh* const combinedCH)
1260{
1261 pts.Resize(0);
1262 AddPoints(ch1, pts);
1263 AddPoints(ch2, pts);
1264
1265 btConvexHullComputer ch;
1266 ch.compute((double*)pts.Data(), 3 * sizeof(double), (int32_t)pts.Size(), -1.0, -1.0);
1267 combinedCH->ResizePoints(0);
1268 combinedCH->ResizeTriangles(0);
1269 for (int32_t v = 0; v < ch.vertices.size(); v++) {
1270 combinedCH->AddPoint(Vec3<double>(ch.vertices[v].getX(), ch.vertices[v].getY(), ch.vertices[v].getZ()));
1271 }
1272 const int32_t nt = ch.faces.size();
1273 for (int32_t t = 0; t < nt; ++t) {
1274 const btConvexHullComputer::Edge* sourceEdge = &(ch.edges[ch.faces[t]]);
1275 int32_t a = sourceEdge->getSourceVertex();
1276 int32_t b = sourceEdge->getTargetVertex();
1277 const btConvexHullComputer::Edge* edge = sourceEdge->getNextEdgeOfFace();
1278 int32_t c = edge->getTargetVertex();
1279 while (c != a) {
1280 combinedCH->AddTriangle(Vec3<int32_t>(a, b, c));
1281 edge = edge->getNextEdgeOfFace();
1282 b = c;
1283 c = edge->getTargetVertex();
1284 }
1285 }
1286}
1287void VHACD::MergeConvexHulls(const Parameters& params)
1288{
1289 if (GetCancel()) {
1290 return;
1291 }
1292 m_timer.Tic();
1293
1294 m_stage = "Merge Convex Hulls";
1295
1296 std::ostringstream msg;
1297 if (params.m_logger) {
1298 msg << "+ " << m_stage << std::endl;
1299 params.m_logger->Log(msg.str().c_str());
1300 }
1301
1302 // Get the current number of convex hulls
1303 size_t nConvexHulls = m_convexHulls.Size();
1304 // Iteration counter
1305 int32_t iteration = 0;
1306 // While we have more than at least one convex hull and the user has not asked us to cancel the operation
1307 if (nConvexHulls > 1 && !m_cancel)
1308 {
1309 // Get the gamma error threshold for when to exit
1310 SArray<Vec3<double> > pts;
1311 Mesh combinedCH;
1312
1313 // Populate the cost matrix
1314 size_t idx = 0;
1315 SArray<float> costMatrix;
1316 costMatrix.Resize(((nConvexHulls * nConvexHulls) - nConvexHulls) >> 1);
1317 for (size_t p1 = 1; p1 < nConvexHulls; ++p1)
1318 {
1319 const float volume1 = m_convexHulls[p1]->ComputeVolume();
1320 for (size_t p2 = 0; p2 < p1; ++p2)
1321 {
1322 ComputeConvexHull(m_convexHulls[p1], m_convexHulls[p2], pts, &combinedCH);
1323 costMatrix[idx++] = ComputeConcavity(volume1 + m_convexHulls[p2]->ComputeVolume(), combinedCH.ComputeVolume(), m_volumeCH0);
1324 }
1325 }
1326
1327 // Until we cant merge below the maximum cost
1328 size_t costSize = m_convexHulls.Size();
1329 while (!m_cancel)
1330 {
1331 msg.str("");
1332 msg << "Iteration " << iteration++;
1333 m_operation = msg.str();
1334
1335 // Search for lowest cost
1336 float bestCost = (std::numeric_limits<float>::max)();
1337 const size_t addr = FindMinimumElement(costMatrix.Data(), &bestCost, 0, costMatrix.Size());
1338 if ( (costSize-1) < params.m_maxConvexHulls)
1339 {
1340 break;
1341 }
1342 const size_t addrI = (static_cast<int32_t>(sqrt(1 + (8 * addr))) - 1) >> 1;
1343 const size_t p1 = addrI + 1;
1344 const size_t p2 = addr - ((addrI * (addrI + 1)) >> 1);
1345 assert(p1 >= 0);
1346 assert(p2 >= 0);
1347 assert(p1 < costSize);
1348 assert(p2 < costSize);
1349
1350 if (params.m_logger)
1351 {
1352 msg.str("");
1353 msg << "\t\t Merging (" << p1 << ", " << p2 << ") " << bestCost << std::endl
1354 << std::endl;
1355 params.m_logger->Log(msg.str().c_str());
1356 }
1357
1358 // Make the lowest cost row and column into a new hull
1359 Mesh* cch = new Mesh;
1360 ComputeConvexHull(m_convexHulls[p1], m_convexHulls[p2], pts, cch);
1361 delete m_convexHulls[p2];
1362 m_convexHulls[p2] = cch;
1363
1364 delete m_convexHulls[p1];
1365 std::swap(m_convexHulls[p1], m_convexHulls[m_convexHulls.Size() - 1]);
1366 m_convexHulls.PopBack();
1367
1368 costSize = costSize - 1;
1369
1370 // Calculate costs versus the new hull
1371 size_t rowIdx = ((p2 - 1) * p2) >> 1;
1372 const float volume1 = m_convexHulls[p2]->ComputeVolume();
1373 for (size_t i = 0; (i < p2) && (!m_cancel); ++i)
1374 {
1375 ComputeConvexHull(m_convexHulls[p2], m_convexHulls[i], pts, &combinedCH);
1376 costMatrix[rowIdx++] = ComputeConcavity(volume1 + m_convexHulls[i]->ComputeVolume(), combinedCH.ComputeVolume(), m_volumeCH0);
1377 }
1378
1379 rowIdx += p2;
1380 for (size_t i = p2 + 1; (i < costSize) && (!m_cancel); ++i)
1381 {
1382 ComputeConvexHull(m_convexHulls[p2], m_convexHulls[i], pts, &combinedCH);
1383 costMatrix[rowIdx] = ComputeConcavity(volume1 + m_convexHulls[i]->ComputeVolume(), combinedCH.ComputeVolume(), m_volumeCH0);
1384 rowIdx += i;
1385 assert(rowIdx >= 0);
1386 }
1387
1388 // Move the top column in to replace its space
1389 const size_t erase_idx = ((costSize - 1) * costSize) >> 1;
1390 if (p1 < costSize) {
1391 rowIdx = (addrI * p1) >> 1;
1392 size_t top_row = erase_idx;
1393 for (size_t i = 0; i < p1; ++i) {
1394 if (i != p2) {
1395 costMatrix[rowIdx] = costMatrix[top_row];
1396 }
1397 ++rowIdx;
1398 ++top_row;
1399 }
1400
1401 ++top_row;
1402 rowIdx += p1;
1403 for (size_t i = p1 + 1; i < (costSize + 1); ++i) {
1404 costMatrix[rowIdx] = costMatrix[top_row++];
1405 rowIdx += i;
1406 assert(rowIdx >= 0);
1407 }
1408 }
1409 costMatrix.Resize(erase_idx);
1410 }
1411 }
1412 m_overallProgress = 99.0;
1413 Update(100.0, 100.0, params);
1414 m_timer.Toc();
1415 if (params.m_logger) {
1416 msg.str("");
1417 msg << "\t time " << m_timer.GetElapsedTime() / 1000.0 << "s" << std::endl;
1418 params.m_logger->Log(msg.str().c_str());
1419 }
1420}
1421void VHACD::SimplifyConvexHull(Mesh* const ch, const size_t nvertices, const double minVolume)
1422{
1423 if (nvertices <= 4) {
1424 return;
1425 }
1426 ICHull icHull;
1427 if (mRaycastMesh)
1428 {
1429 // We project these points onto the original source mesh to increase precision
1430 // The voxelization process drops floating point precision so returned data points are not exactly lying on the
1431 // surface of the original source mesh.
1432 // The first step is we need to compute the bounding box of the mesh we are trying to build a convex hull for.
1433 // From this bounding box, we compute the length of the diagonal to get a relative size and center for point projection
1434 uint32_t nPoints = ch->GetNPoints();
1435 Vec3<double> *inputPoints = ch->GetPointsBuffer();
1436 Vec3<double> bmin(inputPoints[0]);
1437 Vec3<double> bmax(inputPoints[1]);
1438 for (uint32_t i = 1; i < nPoints; i++)
1439 {
1440 const Vec3<double> &p = inputPoints[i];
1441 p.UpdateMinMax(bmin, bmax);
1442 }
1443 Vec3<double> center;
1444 double diagonalLength = center.GetCenter(bmin, bmax); // Get the center of the bounding box
1445 // This is the error threshold for determining if we should use the raycast result data point vs. the voxelized result.
1446 double pointDistanceThreshold = diagonalLength * 0.05;
1447 // If a new point is within 1/100th the diagonal length of the bounding volume we do not add it. To do so would create a
1448 // thin sliver in the resulting convex hull
1449 double snapDistanceThreshold = diagonalLength * 0.01;
1450 double snapDistanceThresholdSquared = snapDistanceThreshold*snapDistanceThreshold;
1451
1452 // Allocate buffer for projected vertices
1453 Vec3<double> *outputPoints = new Vec3<double>[nPoints];
1454 uint32_t outCount = 0;
1455 for (uint32_t i = 0; i < nPoints; i++)
1456 {
1457 Vec3<double> &inputPoint = inputPoints[i];
1458 Vec3<double> &outputPoint = outputPoints[outCount];
1459 // Compute the direction vector from the center of this mesh to the vertex
1460 Vec3<double> dir = inputPoint - center;
1461 // Normalize the direction vector.
1462 dir.Normalize();
1463 // Multiply times the diagonal length of the mesh
1464 dir *= diagonalLength;
1465 // Add the center back in again to get the destination point
1466 dir += center;
1467 // By default the output point is equal to the input point
1468 outputPoint = inputPoint;
1469 double pointDistance;
1470 if (mRaycastMesh->raycast(center.GetData(), dir.GetData(), inputPoint.GetData(), outputPoint.GetData(),&pointDistance) )
1471 {
1472 // If the nearest intersection point is too far away, we keep the original source data point.
1473 // Not all points lie directly on the original mesh surface
1474 if (pointDistance > pointDistanceThreshold)
1475 {
1476 outputPoint = inputPoint;
1477 }
1478 }
1479 // Ok, before we add this point, we do not want to create points which are extremely close to each other.
1480 // This will result in tiny sliver triangles which are really bad for collision detection.
1481 bool foundNearbyPoint = false;
1482 for (uint32_t j = 0; j < outCount; j++)
1483 {
1484 // If this new point is extremely close to an existing point, we do not add it!
1485 double squaredDistance = outputPoints[j].GetDistanceSquared(outputPoint);
1486 if (squaredDistance < snapDistanceThresholdSquared )
1487 {
1488 foundNearbyPoint = true;
1489 break;
1490 }
1491 }
1492 if (!foundNearbyPoint)
1493 {
1494 outCount++;
1495 }
1496 }
1497 icHull.AddPoints(outputPoints, outCount);
1498 delete[]outputPoints;
1499 }
1500 else
1501 {
1502 icHull.AddPoints(ch->GetPointsBuffer(), ch->GetNPoints());
1503 }
1504 icHull.Process((uint32_t)nvertices, minVolume);
1505 TMMesh& mesh = icHull.GetMesh();
1506 const size_t nT = mesh.GetNTriangles();
1507 const size_t nV = mesh.GetNVertices();
1508 ch->ResizePoints(nV);
1509 ch->ResizeTriangles(nT);
1510 mesh.GetIFS(ch->GetPointsBuffer(), ch->GetTrianglesBuffer());
1511}
1512void VHACD::SimplifyConvexHulls(const Parameters& params)
1513{
1514 if (m_cancel || params.m_maxNumVerticesPerCH < 4) {
1515 return;
1516 }
1517 m_timer.Tic();
1518
1519 m_stage = "Simplify convex-hulls";
1520 m_operation = "Simplify convex-hulls";
1521
1522 std::ostringstream msg;
1523 const size_t nConvexHulls = m_convexHulls.Size();
1524 if (params.m_logger) {
1525 msg << "+ Simplify " << nConvexHulls << " convex-hulls " << std::endl;
1526 params.m_logger->Log(msg.str().c_str());
1527 }
1528
1529 Update(0.0, 0.0, params);
1530 for (size_t i = 0; i < nConvexHulls && !m_cancel; ++i) {
1531 if (params.m_logger) {
1532 msg.str("");
1533 msg << "\t\t Simplify CH[" << std::setfill('0') << std::setw(5) << i << "] " << m_convexHulls[i]->GetNPoints() << " V, " << m_convexHulls[i]->GetNTriangles() << " T" << std::endl;
1534 params.m_logger->Log(msg.str().c_str());
1535 }
1536 SimplifyConvexHull(m_convexHulls[i], params.m_maxNumVerticesPerCH, m_volumeCH0 * params.m_minVolumePerCH);
1537 }
1538
1539 m_overallProgress = 100.0;
1540 Update(100.0, 100.0, params);
1541 m_timer.Toc();
1542 if (params.m_logger) {
1543 msg.str("");
1544 msg << "\t time " << m_timer.GetElapsedTime() / 1000.0 << "s" << std::endl;
1545 params.m_logger->Log(msg.str().c_str());
1546 }
1547}
1548
1549bool VHACD::ComputeCenterOfMass(double centerOfMass[3]) const
1550{
1551 bool ret = false;
1552
1553 centerOfMass[0] = 0;
1554 centerOfMass[1] = 0;
1555 centerOfMass[2] = 0;
1556 // Get number of convex hulls in the result
1557 uint32_t hullCount = GetNConvexHulls();
1558 if (hullCount) // if we have results
1559 {
1560 ret = true;
1561 double totalVolume = 0;
1562 // Initialize the center of mass to zero
1563 centerOfMass[0] = 0;
1564 centerOfMass[1] = 0;
1565 centerOfMass[2] = 0;
1566 // Compute the total volume of all convex hulls
1567 for (uint32_t i = 0; i < hullCount; i++)
1568 {
1569 ConvexHull ch;
1570 GetConvexHull(i, ch);
1571 totalVolume += ch.m_volume;
1572 }
1573 // compute the reciprocal of the total volume
1574 double recipVolume = 1.0 / totalVolume;
1575 // Add in the weighted by volume average of the center point of each convex hull
1576 for (uint32_t i = 0; i < hullCount; i++)
1577 {
1578 ConvexHull ch;
1579 GetConvexHull(i, ch);
1580 double ratio = ch.m_volume*recipVolume;
1581 centerOfMass[0] += ch.m_center[0] * ratio;
1582 centerOfMass[1] += ch.m_center[1] * ratio;
1583 centerOfMass[2] += ch.m_center[2] * ratio;
1584 }
1585 }
1586 return ret;
1587}
1588
1589} // end of VHACD namespace
1590