Revision: 71941
http://sourceforge.net/p/brlcad/code/71941
Author: starseeker
Date: 2018-10-03 22:13:09 +0000 (Wed, 03 Oct 2018)
Log Message:
-----------
Replace the Clarkson convex hull code with Antti Kuukka's C++11 implementation
of quickhull from https://github.com/akuukka/quickhull - a timing comparison
with the Stanford lucy model reported 2.5 seconds with the new method vs. 11
minutes 45.5 seconds with the old method. (Part of the time increase for the
old method was due to the overhead of swap space usage - memory requirements
also appear to be *much* lower for the new code.)
Modified Paths:
--------------
brlcad/trunk/CMakeLists.txt
brlcad/trunk/doc/legal/embedded/chull3d.txt
brlcad/trunk/misc/CMake/CMakeLists.txt
brlcad/trunk/src/libbg/CMakeLists.txt
Added Paths:
-----------
brlcad/trunk/misc/CMake/CheckCXX11Features/
brlcad/trunk/misc/CMake/CheckCXX11Features/cxx11-test-nullptr.cpp
brlcad/trunk/misc/CMake/CheckCXX11Features/cxx11-test-nullptr_fail_compile.cpp
brlcad/trunk/misc/CMake/CheckCXX11Features.cmake
brlcad/trunk/src/libbg/QuickHull.cxx
brlcad/trunk/src/libbg/QuickHull.hpp
brlcad/trunk/src/libbg/chull3d.cxx
Removed Paths:
-------------
brlcad/trunk/src/libbg/chull3d.cpp
Modified: brlcad/trunk/CMakeLists.txt
===================================================================
--- brlcad/trunk/CMakeLists.txt 2018-10-03 13:53:28 UTC (rev 71940)
+++ brlcad/trunk/CMakeLists.txt 2018-10-03 22:13:09 UTC (rev 71941)
@@ -2583,6 +2583,13 @@
endif (HAVE___THREAD)
endif (HAVE_THREAD_LOCAL)
+# Test for C++11 features that we need - CMake will sometimes assert that a
+# C++11 compiler flag exists, but the actual compiler support for C++11
language
+# ends up lacking in practice.
+include(CheckCXX11Features)
+cxx11_check_feature("nullptr" HAS_CXX11_NULLPTR)
+
+
# see if the compiler supports %z as a size_t print width specifier
BRLCAD_CHECK_PERCENT_Z()
@@ -2589,7 +2596,10 @@
# see if the compiler supports [static] array hints
BRLCAD_CHECK_STATIC_ARRAYS()
+# see if we have a TLS intrinsic, first check C++11 compliance
+check_cxx_source_compiles("static thread_local int i = 0; int main() { return
i; }" HAVE_THREAD_LOCAL)
+
# *******************************************************************
if(BRLCAD_PRINT_MSGS)
message("***********************************************************")
Modified: brlcad/trunk/doc/legal/embedded/chull3d.txt
===================================================================
--- brlcad/trunk/doc/legal/embedded/chull3d.txt 2018-10-03 13:53:28 UTC (rev
71940)
+++ brlcad/trunk/doc/legal/embedded/chull3d.txt 2018-10-03 22:13:09 UTC (rev
71941)
@@ -1,13 +1,10 @@
-http://www.netlib.org/voronoi/hull.html
+https://github.com/akuukka/quickhull
-Ken Clarkson wrote this. Copyright (c) 1995 by AT&T..
-Permission to use, copy, modify, and distribute this software for any
-purpose without fee is hereby granted, provided that this entire notice
-is included in all copies of any software which is or includes a copy
-or modification of this software and in all copies of the supporting
-documentation for such software.
-THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
-WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
-REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
-OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
+Implementation of the 3D QuickHull algorithm by Antti Kuukka
+Per README.md:
+
+This implementation is 100% Public Domain.
+
+Feel free to use.
+
Modified: brlcad/trunk/misc/CMake/CMakeLists.txt
===================================================================
--- brlcad/trunk/misc/CMake/CMakeLists.txt 2018-10-03 13:53:28 UTC (rev
71940)
+++ brlcad/trunk/misc/CMake/CMakeLists.txt 2018-10-03 22:13:09 UTC (rev
71941)
@@ -9,6 +9,9 @@
CMakeLists.txt
CheckCInline.cmake
CheckCSourceRuns.cmake
+ CheckCXX11Features.cmake
+ CheckCXX11Features/cxx11-test-nullptr.cpp
+ CheckCXX11Features/cxx11-test-nullptr_fail_compile.cpp
CompilerFlags.cmake
Distcheck.cmake
DocBook.cmake
Added: brlcad/trunk/misc/CMake/CheckCXX11Features/cxx11-test-nullptr.cpp
===================================================================
--- brlcad/trunk/misc/CMake/CheckCXX11Features/cxx11-test-nullptr.cpp
(rev 0)
+++ brlcad/trunk/misc/CMake/CheckCXX11Features/cxx11-test-nullptr.cpp
2018-10-03 22:13:09 UTC (rev 71941)
@@ -0,0 +1,6 @@
+int main(void)
+{
+ void *v = nullptr;
+
+ return v ? 1 : 0;
+}
Property changes on:
brlcad/trunk/misc/CMake/CheckCXX11Features/cxx11-test-nullptr.cpp
___________________________________________________________________
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
Added: svn:mime-type
## -0,0 +1 ##
+text/plain
\ No newline at end of property
Added:
brlcad/trunk/misc/CMake/CheckCXX11Features/cxx11-test-nullptr_fail_compile.cpp
===================================================================
---
brlcad/trunk/misc/CMake/CheckCXX11Features/cxx11-test-nullptr_fail_compile.cpp
(rev 0)
+++
brlcad/trunk/misc/CMake/CheckCXX11Features/cxx11-test-nullptr_fail_compile.cpp
2018-10-03 22:13:09 UTC (rev 71941)
@@ -0,0 +1,6 @@
+int main(void)
+{
+ int i = nullptr;
+
+ return 1;
+}
Property changes on:
brlcad/trunk/misc/CMake/CheckCXX11Features/cxx11-test-nullptr_fail_compile.cpp
___________________________________________________________________
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
Added: svn:mime-type
## -0,0 +1 ##
+text/plain
\ No newline at end of property
Added: brlcad/trunk/misc/CMake/CheckCXX11Features.cmake
===================================================================
--- brlcad/trunk/misc/CMake/CheckCXX11Features.cmake
(rev 0)
+++ brlcad/trunk/misc/CMake/CheckCXX11Features.cmake 2018-10-03 22:13:09 UTC
(rev 71941)
@@ -0,0 +1,110 @@
+# - Check which parts of the C++11 standard the compiler supports
+#
+# Copyright 2011,2012 Rolf Eike Beer <e...@sf-mail.de>
+# Copyright 2012 Andreas Weis
+# Copyright 2014 Ben Morgan <bmorgan.warw...@gmail.com>
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+# * Redistributions of source code must retain the above copyright
+# notice, this list of conditions and the following disclaimer.
+# * 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.
+# * Neither the name of the copyright holders nor the
+# names of its contributors may be used to endorse or promote products
+# derived from this software without specific prior written permission.
+#
+# 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 HOLDERS 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.
+
+#
+# Each feature may have up to 3 checks, every one of them in it's own file
+# FEATURE.cpp - example that must build and return 0 when run
+# FEATURE_fail.cpp - example that must build, but may not return 0
when run
+# FEATURE_fail_compile.cpp - example that must fail compilation
+#
+# The first one is mandatory, the latter 2 are optional and do not depend on
+# each other (i.e. only one may be present).
+#
+
+if(NOT CMAKE_CXX_COMPILER_LOADED)
+ message(FATAL_ERROR "CheckCXX11Features modules only works if language CXX
is enabled")
+endif()
+
+#-----------------------------------------------------------------------
+# function cxx11_check_feature(<name> <result>)
+#
+function(cxx11_check_feature FEATURE_NAME RESULT_VAR)
+ if(NOT DEFINED ${RESULT_VAR})
+ set(_bindir
"${CMAKE_CURRENT_BINARY_DIR}/CheckCXX11Features/cxx11_${FEATURE_NAME}")
+
+ set(_SRCFILE_BASE
${CMAKE_SOURCE_DIR}/misc/CMake/CheckCXX11Features/cxx11-test-${FEATURE_NAME})
+ set(_LOG_NAME "\"${FEATURE_NAME}\"")
+ message(STATUS "Checking C++11 support for ${_LOG_NAME}")
+
+ set(_SRCFILE "${_SRCFILE_BASE}.cpp")
+ set(_SRCFILE_FAIL "${_SRCFILE_BASE}_fail.cpp")
+ set(_SRCFILE_FAIL_COMPILE "${_SRCFILE_BASE}_fail_compile.cpp")
+
+ if(CMAKE_CROSSCOMPILING)
+ try_compile(${RESULT_VAR} "${_bindir}" "${_SRCFILE}"
+ COMPILE_DEFINITIONS "${CXX11_COMPILER_FLAGS}")
+ if(${RESULT_VAR} AND EXISTS ${_SRCFILE_FAIL})
+ try_compile(${RESULT_VAR} "${_bindir}_fail" "${_SRCFILE_FAIL}"
+ COMPILE_DEFINITIONS "${CXX11_COMPILER_FLAGS}")
+ endif()
+ else()
+ try_run(_RUN_RESULT_VAR _COMPILE_RESULT_VAR
+ "${_bindir}" "${_SRCFILE}"
+ #CMAKE_FLAGS "-DCMAKE_XCODE_ATTRIBUTE_CLANG_CXX_LIBRARY='libc++'"
+ COMPILE_DEFINITIONS "${CXX11_COMPILER_FLAGS}")
+ if(_COMPILE_RESULT_VAR AND NOT _RUN_RESULT_VAR)
+ set(${RESULT_VAR} TRUE)
+ else()
+ set(${RESULT_VAR} FALSE)
+ endif()
+
+ if(${RESULT_VAR} AND EXISTS ${_SRCFILE_FAIL})
+ try_run(_RUN_RESULT_VAR _COMPILE_RESULT_VAR
+ "${_bindir}_fail" "${_SRCFILE_FAIL}"
+ #CMAKE_FLAGS "-DCMAKE_XCODE_ATTRIBUTE_CLANG_CXX_LIBRARY='libc++'"
+ COMPILE_DEFINITIONS "${CXX11_COMPILER_FLAGS}")
+ if(_COMPILE_RESULT_VAR AND _RUN_RESULT_VAR)
+ set(${RESULT_VAR} TRUE)
+ else()
+ set(${RESULT_VAR} FALSE)
+ endif()
+ endif()
+ endif()
+
+ if(${RESULT_VAR} AND EXISTS ${_SRCFILE_FAIL_COMPILE})
+ try_compile(_TMP_RESULT "${_bindir}_fail_compile"
"${_SRCFILE_FAIL_COMPILE}"
+ #CMAKE_FLAGS "-DCMAKE_XCODE_ATTRIBUTE_CLANG_CXX_LIBRARY='libc++'"
+ COMPILE_DEFINITIONS "${CXX11_COMPILER_FLAGS}")
+ if(_TMP_RESULT)
+ set(${RESULT_VAR} FALSE)
+ else()
+ set(${RESULT_VAR} TRUE)
+ endif()
+ endif()
+
+ if(${RESULT_VAR})
+ message(STATUS "Checking C++11 support for ${_LOG_NAME}: works")
+ else()
+ message(FATAL_ERROR "\nChecking C++11 support for ${_LOG_NAME}: not
supported\n\n")
+ endif()
+
+ set(${RESULT_VAR} ${${RESULT_VAR}} CACHE INTERNAL "C++11 support for
${_LOG_NAME}")
+ endif()
+endfunction()
+
Property changes on: brlcad/trunk/misc/CMake/CheckCXX11Features.cmake
___________________________________________________________________
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
Added: svn:mime-type
## -0,0 +1 ##
+text/plain
\ No newline at end of property
Modified: brlcad/trunk/src/libbg/CMakeLists.txt
===================================================================
--- brlcad/trunk/src/libbg/CMakeLists.txt 2018-10-03 13:53:28 UTC (rev
71940)
+++ brlcad/trunk/src/libbg/CMakeLists.txt 2018-10-03 22:13:09 UTC (rev
71941)
@@ -9,7 +9,8 @@
set(LIBBG_SOURCES
chull.c
- chull3d.cpp
+ QuickHull.cxx
+ chull3d.cxx
obr.c
polygon.c
polygon_ear_clipping.c
@@ -29,7 +30,13 @@
add_subdirectory(tests)
-CMAKEFILES(bg_private.h pointgen.c)
+set(bg_ignore
+ QuickHull.hpp
+ bg_private.h
+ pointgen.c
+ )
+
+CMAKEFILES(${bg_ignore})
CMAKEFILES(CMakeLists.txt)
# Local Variables:
Added: brlcad/trunk/src/libbg/QuickHull.cxx
===================================================================
--- brlcad/trunk/src/libbg/QuickHull.cxx (rev 0)
+++ brlcad/trunk/src/libbg/QuickHull.cxx 2018-10-03 22:13:09 UTC (rev
71941)
@@ -0,0 +1,529 @@
+/*
+ * QuickHull.hpp
+ *
+ * QuickHull implementation from https://github.com/akuukka/quickhull by Antti
+ * Kuukka
+ *
+ * Per the README.md file:
+ *
+ * This implementation is 100% Public Domain.
+ *
+ * Feel free to use.
+ *
+ * C++11 is needed to compile it.
+ */
+
+#include "common.h"
+#include <cmath>
+#include <cassert>
+#include <iostream>
+#include <algorithm>
+#include <deque>
+#include <limits>
+#include "QuickHull.hpp"
+
+namespace quickhull {
+
+ template<>
+ const float QuickHull<float>::Epsilon = 0.0001f;
+
+ template<>
+ const double QuickHull<double>::Epsilon = 0.0000001;
+
+ /*
+ * Implementation of the algorithm
+ */
+
+ template<typename T>
+ ConvexHull<T> QuickHull<T>::getConvexHull(const
std::vector<Vector3<T>>& pointCloud, bool CCW, bool useOriginalIndices, T
epsilon) {
+ VertexDataSource<T> vertexDataSource(pointCloud);
+ return
getConvexHull(vertexDataSource,CCW,useOriginalIndices,epsilon);
+ }
+
+ template<typename T>
+ ConvexHull<T> QuickHull<T>::getConvexHull(const Vector3<T>* vertexData,
size_t vertexCount, bool CCW, bool useOriginalIndices, T epsilon) {
+ VertexDataSource<T> vertexDataSource(vertexData,vertexCount);
+ return
getConvexHull(vertexDataSource,CCW,useOriginalIndices,epsilon);
+ }
+
+ template<typename T>
+ ConvexHull<T> QuickHull<T>::getConvexHull(const T* vertexData, size_t
vertexCount, bool CCW, bool useOriginalIndices, T epsilon) {
+ VertexDataSource<T> vertexDataSource((const
vec3*)vertexData,vertexCount);
+ return
getConvexHull(vertexDataSource,CCW,useOriginalIndices,epsilon);
+ }
+
+ template<typename FloatType>
+ HalfEdgeMesh<FloatType, IndexType>
QuickHull<FloatType>::getConvexHullAsMesh(const FloatType* vertexData, size_t
vertexCount, bool CCW, FloatType epsilon) {
+ VertexDataSource<FloatType> vertexDataSource((const
vec3*)vertexData,vertexCount);
+ buildMesh(vertexDataSource, CCW, false, epsilon);
+ return HalfEdgeMesh<FloatType, IndexType>(m_mesh, m_vertexData);
+ }
+
+ template<typename T>
+ void QuickHull<T>::buildMesh(const VertexDataSource<T>& pointCloud,
bool CCW, bool useOriginalIndices, T epsilon) {
+ if (pointCloud.size()==0) {
+ m_mesh = MeshBuilder<T>();
+ return;
+ }
+ m_vertexData = pointCloud;
+
+ // Very first: find extreme values and use them to compute the
scale of the point cloud.
+ m_extremeValues = getExtremeValues();
+ m_scale = getScale(m_extremeValues);
+
+ // Epsilon we use depends on the scale
+ m_epsilon = epsilon*m_scale;
+ m_epsilonSquared = m_epsilon*m_epsilon;
+
+ // Reset diagnostics
+ m_diagnostics = DiagnosticsData();
+
+ m_planar = false; // The planar case happens when all the points
appear to lie on a two dimensional subspace of R^3.
+ createConvexHalfEdgeMesh();
+ if (m_planar) {
+ const size_t extraPointIndex = m_planarPointCloudTemp.size()-1;
+ for (auto& he : m_mesh.m_halfEdges) {
+ if (he.m_endVertex == extraPointIndex) {
+ he.m_endVertex = 0;
+ }
+ }
+ m_vertexData = pointCloud;
+ m_planarPointCloudTemp.clear();
+ }
+ }
+
+ template<typename T>
+ ConvexHull<T> QuickHull<T>::getConvexHull(const VertexDataSource<T>&
pointCloud, bool CCW, bool useOriginalIndices, T epsilon) {
+ buildMesh(pointCloud,CCW,useOriginalIndices,epsilon);
+ return ConvexHull<T>(m_mesh,m_vertexData, CCW, useOriginalIndices);
+ }
+
+ template<typename T>
+ void QuickHull<T>::createConvexHalfEdgeMesh() {
+ // Temporary variables used during iteration
+ std::vector<IndexType> visibleFaces;
+ std::vector<IndexType> horizonEdges;
+ struct FaceData {
+ IndexType m_faceIndex;
+ IndexType m_enteredFromHalfEdge; // If the face turns out not
to be visible, this half edge will be marked as horizon edge
+ FaceData(IndexType fi, IndexType he) :
m_faceIndex(fi),m_enteredFromHalfEdge(he) {}
+ };
+ std::vector<FaceData> possiblyVisibleFaces;
+
+ // Compute base tetrahedron
+ m_mesh = getInitialTetrahedron();
+ assert(m_mesh.m_faces.size()==4);
+
+ // Init face stack with those faces that have points assigned to
them
+ std::deque<IndexType> faceList;
+ for (size_t i=0;i < 4;i++) {
+ auto& f = m_mesh.m_faces[i];
+ if (f.m_pointsOnPositiveSide &&
f.m_pointsOnPositiveSide->size()>0) {
+ faceList.push_back(i);
+ f.m_inFaceStack = 1;
+ }
+ }
+
+ // Process faces until the face list is empty.
+ size_t iter = 0;
+ while (!faceList.empty()) {
+ iter++;
+ if (iter == std::numeric_limits<size_t>::max()) {
+ // Visible face traversal marks visited faces with
iteration counter (to mark that the face has been visited on this iteration)
and the max value represents unvisited faces. At this point we have to reset
iteration counter. This shouldn't be an
+ // issue on 64 bit machines.
+ iter = 0;
+ }
+
+ const IndexType topFaceIndex = faceList.front();
+ faceList.pop_front();
+
+ auto& tf = m_mesh.m_faces[topFaceIndex];
+ tf.m_inFaceStack = 0;
+
+ assert(!tf.m_pointsOnPositiveSide ||
tf.m_pointsOnPositiveSide->size()>0);
+ if (!tf.m_pointsOnPositiveSide || tf.isDisabled()) {
+ continue;
+ }
+
+ // Pick the most distant point to this triangle plane as the
point to which we extrude
+ const vec3& activePoint = m_vertexData[tf.m_mostDistantPoint];
+ const size_t activePointIndex = tf.m_mostDistantPoint;
+
+ // Find out the faces that have our active point on their
positive side (these are the "visible faces"). The face on top of the stack of
course is one of them. At the same time, we create a list of horizon edges.
+ horizonEdges.clear();
+ possiblyVisibleFaces.clear();
+ visibleFaces.clear();
+
possiblyVisibleFaces.emplace_back(topFaceIndex,std::numeric_limits<size_t>::max());
+ while (possiblyVisibleFaces.size()) {
+ const auto faceData = possiblyVisibleFaces.back();
+ possiblyVisibleFaces.pop_back();
+ auto& pvf = m_mesh.m_faces[faceData.m_faceIndex];
+ assert(!pvf.isDisabled());
+
+ if (pvf.m_visibilityCheckedOnIteration == iter) {
+ if (pvf.m_isVisibleFaceOnCurrentIteration) {
+ continue;
+ }
+ }
+ else {
+ const Plane<T>& P = pvf.m_P;
+ pvf.m_visibilityCheckedOnIteration = iter;
+ const T d = P.m_N.dotProduct(activePoint)+P.m_D;
+ if (d>0) {
+ pvf.m_isVisibleFaceOnCurrentIteration = 1;
+ pvf.m_horizonEdgesOnCurrentIteration = 0;
+ visibleFaces.push_back(faceData.m_faceIndex);
+ for (auto heIndex :
m_mesh.getHalfEdgeIndicesOfFace(pvf)) {
+ if (m_mesh.m_halfEdges[heIndex].m_opp !=
faceData.m_enteredFromHalfEdge) {
+ possiblyVisibleFaces.emplace_back(
m_mesh.m_halfEdges[m_mesh.m_halfEdges[heIndex].m_opp].m_face,heIndex );
+ }
+ }
+ continue;
+ }
+ assert(faceData.m_faceIndex != topFaceIndex);
+ }
+
+ // The face is not visible. Therefore, the halfedge we came
from is part of the horizon edge.
+ pvf.m_isVisibleFaceOnCurrentIteration = 0;
+ horizonEdges.push_back(faceData.m_enteredFromHalfEdge);
+ // Store which half edge is the horizon edge. The other
half edges of the face will not be part of the final mesh so their data slots
can by recycled.
+ const auto halfEdges =
m_mesh.getHalfEdgeIndicesOfFace(m_mesh.m_faces[m_mesh.m_halfEdges[faceData.m_enteredFromHalfEdge].m_face]);
+ const std::int8_t ind =
(halfEdges[0]==faceData.m_enteredFromHalfEdge) ? 0 :
(halfEdges[1]==faceData.m_enteredFromHalfEdge ? 1 : 2);
+
m_mesh.m_faces[m_mesh.m_halfEdges[faceData.m_enteredFromHalfEdge].m_face].m_horizonEdgesOnCurrentIteration
|= (1<<ind);
+ }
+ const size_t horizonEdgeCount = horizonEdges.size();
+
+ // Order horizon edges so that they form a loop. This may fail
due to numerical instability in which case we give up trying to solve horizon
edge for this point and accept a minor degeneration in the convex hull.
+ if (!reorderHorizonEdges(horizonEdges)) {
+ m_diagnostics.m_failedHorizonEdges++;
+ std::cerr << "Failed to solve horizon edge." << std::endl;
+ auto it =
std::find(tf.m_pointsOnPositiveSide->begin(),tf.m_pointsOnPositiveSide->end(),activePointIndex);
+ tf.m_pointsOnPositiveSide->erase(it);
+ if (tf.m_pointsOnPositiveSide->size()==0) {
+ reclaimToIndexVectorPool(tf.m_pointsOnPositiveSide);
+ }
+ continue;
+ }
+
+ // Except for the horizon edges, all half edges of the visible
faces can be marked as disabled. Their data slots will be reused.
+ // The faces will be disabled as well, but we need to remember
the points that were on the positive side of them - therefore
+ // we save pointers to them.
+ m_newFaceIndices.clear();
+ m_newHalfEdgeIndices.clear();
+ m_disabledFacePointVectors.clear();
+ size_t disableCounter = 0;
+ for (auto faceIndex : visibleFaces) {
+ auto& disabledFace = m_mesh.m_faces[faceIndex];
+ auto halfEdges =
m_mesh.getHalfEdgeIndicesOfFace(disabledFace);
+ for (size_t j=0;j<3;j++) {
+ if ((disabledFace.m_horizonEdgesOnCurrentIteration &
(1<<j)) == 0) {
+ if (disableCounter < horizonEdgeCount*2) {
+ // Use on this iteration
+ m_newHalfEdgeIndices.push_back(halfEdges[j]);
+ disableCounter++;
+ }
+ else {
+ // Mark for reusal on later iteration step
+ m_mesh.disableHalfEdge(halfEdges[j]);
+ }
+ }
+ }
+ // Disable the face, but retain pointer to the points that
were on the positive side of it. We need to assign those points
+ // to the new faces we create shortly.
+ auto t = std::move(m_mesh.disableFace(faceIndex));
+ if (t) {
+ assert(t->size()); // Because we should not assign
point vectors to faces unless needed...
+ m_disabledFacePointVectors.push_back(std::move(t));
+ }
+ }
+ if (disableCounter < horizonEdgeCount*2) {
+ const size_t newHalfEdgesNeeded =
horizonEdgeCount*2-disableCounter;
+ for (size_t i=0;i<newHalfEdgesNeeded;i++) {
+ m_newHalfEdgeIndices.push_back(m_mesh.addHalfEdge());
+ }
+ }
+
+ // Create new faces using the edgeloop
+ for (size_t i = 0; i < horizonEdgeCount; i++) {
+ const IndexType AB = horizonEdges[i];
+
+ auto horizonEdgeVertexIndices =
m_mesh.getVertexIndicesOfHalfEdge(m_mesh.m_halfEdges[AB]);
+ IndexType A,B,C;
+ A = horizonEdgeVertexIndices[0];
+ B = horizonEdgeVertexIndices[1];
+ C = activePointIndex;
+
+ const IndexType newFaceIndex = m_mesh.addFace();
+ m_newFaceIndices.push_back(newFaceIndex);
+
+ const IndexType CA = m_newHalfEdgeIndices[2*i+0];
+ const IndexType BC = m_newHalfEdgeIndices[2*i+1];
+
+ m_mesh.m_halfEdges[AB].m_next = BC;
+ m_mesh.m_halfEdges[BC].m_next = CA;
+ m_mesh.m_halfEdges[CA].m_next = AB;
+
+ m_mesh.m_halfEdges[BC].m_face = newFaceIndex;
+ m_mesh.m_halfEdges[CA].m_face = newFaceIndex;
+ m_mesh.m_halfEdges[AB].m_face = newFaceIndex;
+
+ m_mesh.m_halfEdges[CA].m_endVertex = A;
+ m_mesh.m_halfEdges[BC].m_endVertex = C;
+
+ auto& newFace = m_mesh.m_faces[newFaceIndex];
+
+ const Vector3<T> planeNormal =
mathutils::getTriangleNormal(m_vertexData[A],m_vertexData[B],activePoint);
+ newFace.m_P = Plane<T>(planeNormal,activePoint);
+ newFace.m_he = AB;
+
+ m_mesh.m_halfEdges[CA].m_opp = m_newHalfEdgeIndices[i>0 ?
i*2-1 : 2*horizonEdgeCount-1];
+ m_mesh.m_halfEdges[BC].m_opp =
m_newHalfEdgeIndices[((i+1)*2) % (horizonEdgeCount*2)];
+ }
+
+ // Assign points that were on the positive side of the disabled
faces to the new faces.
+ for (auto& disabledPoints : m_disabledFacePointVectors) {
+ assert(disabledPoints);
+ for (const auto& point : *(disabledPoints)) {
+ if (point == activePointIndex) {
+ continue;
+ }
+ for (size_t j=0;j<horizonEdgeCount;j++) {
+ if
(addPointToFace(m_mesh.m_faces[m_newFaceIndices[j]], point)) {
+ break;
+ }
+ }
+ }
+ // The points are no longer needed: we can move them to the
vector pool for reuse.
+ reclaimToIndexVectorPool(disabledPoints);
+ }
+
+ // Increase face stack size if needed
+ for (const auto newFaceIndex : m_newFaceIndices) {
+ auto& newFace = m_mesh.m_faces[newFaceIndex];
+ if (newFace.m_pointsOnPositiveSide) {
+ assert(newFace.m_pointsOnPositiveSide->size()>0);
+ if (!newFace.m_inFaceStack) {
+ faceList.push_back(newFaceIndex);
+ newFace.m_inFaceStack = 1;
+ }
+ }
+ }
+ }
+
+ // Cleanup
+ m_indexVectorPool.clear();
+ }
+
+ /*
+ * Private helper functions
+ */
+
+ template <typename T>
+ std::array<IndexType,6> QuickHull<T>::getExtremeValues() {
+ std::array<IndexType,6> outIndices{0,0,0,0,0,0};
+ T extremeVals[6] =
{m_vertexData[0].x,m_vertexData[0].x,m_vertexData[0].y,m_vertexData[0].y,m_vertexData[0].z,m_vertexData[0].z};
+ const size_t vCount = m_vertexData.size();
+ for (size_t i=1;i<vCount;i++) {
+ const Vector3<T>& pos = m_vertexData[i];
+ if (pos.x>extremeVals[0]) {
+ extremeVals[0]=pos.x;
+ outIndices[0]=(IndexType)i;
+ }
+ else if (pos.x<extremeVals[1]) {
+ extremeVals[1]=pos.x;
+ outIndices[1]=(IndexType)i;
+ }
+ if (pos.y>extremeVals[2]) {
+ extremeVals[2]=pos.y;
+ outIndices[2]=(IndexType)i;
+ }
+ else if (pos.y<extremeVals[3]) {
+ extremeVals[3]=pos.y;
+ outIndices[3]=(IndexType)i;
+ }
+ if (pos.z>extremeVals[4]) {
+ extremeVals[4]=pos.z;
+ outIndices[4]=(IndexType)i;
+ }
+ else if (pos.z<extremeVals[5]) {
+ extremeVals[5]=pos.z;
+ outIndices[5]=(IndexType)i;
+ }
+ }
+ return outIndices;
+ }
+
+ template<typename T>
+ bool QuickHull<T>::reorderHorizonEdges(std::vector<IndexType>&
horizonEdges) {
+ const size_t horizonEdgeCount = horizonEdges.size();
+ for (size_t i=0;i<horizonEdgeCount-1;i++) {
+ const IndexType endVertex = m_mesh.m_halfEdges[ horizonEdges[i]
].m_endVertex;
+ bool foundNext = false;
+ for (size_t j=i+1;j<horizonEdgeCount;j++) {
+ const IndexType beginVertex = m_mesh.m_halfEdges[
m_mesh.m_halfEdges[horizonEdges[j]].m_opp ].m_endVertex;
+ if (beginVertex == endVertex) {
+ std::swap(horizonEdges[i+1],horizonEdges[j]);
+ foundNext = true;
+ break;
+ }
+ }
+ if (!foundNext) {
+ return false;
+ }
+ }
+ assert(m_mesh.m_halfEdges[ horizonEdges[horizonEdges.size()-1]
].m_endVertex == m_mesh.m_halfEdges[ m_mesh.m_halfEdges[horizonEdges[0]].m_opp
].m_endVertex);
+ return true;
+ }
+
+ template <typename T>
+ T QuickHull<T>::getScale(const std::array<IndexType,6>& extremeValues) {
+ T s = 0;
+ for (size_t i=0;i<6;i++) {
+ const T* v = (const T*)(&m_vertexData[extremeValues[i]]);
+ v += i/2;
+ auto a = std::abs(*v);
+ if (a>s) {
+ s = a;
+ }
+ }
+ return s;
+ }
+
+ template <typename T>
+ MeshBuilder<T> QuickHull<T>::getInitialTetrahedron() {
+ const size_t vertexCount = m_vertexData.size();
+
+ // If we have at most 4 points, just return a degenerate
tetrahedron:
+ if (vertexCount <= 4) {
+ IndexType v[4] =
{0,std::min((size_t)1,vertexCount-1),std::min((size_t)2,vertexCount-1),std::min((size_t)3,vertexCount-1)};
+ const Vector3<T> N =
mathutils::getTriangleNormal(m_vertexData[v[0]],m_vertexData[v[1]],m_vertexData[v[2]]);
+ const Plane<T> trianglePlane(N,m_vertexData[v[0]]);
+ if (trianglePlane.isPointOnPositiveSide(m_vertexData[v[3]])) {
+ std::swap(v[0],v[1]);
+ }
+ return MeshBuilder<T>(v[0],v[1],v[2],v[3]);
+ }
+
+ // Find two most distant extreme points.
+ T maxD = m_epsilonSquared;
+ std::pair<IndexType,IndexType> selectedPoints;
+ for (size_t i=0;i<6;i++) {
+ for (size_t j=i+1;j<6;j++) {
+ const T d = m_vertexData[ m_extremeValues[i]
].getSquaredDistanceTo( m_vertexData[ m_extremeValues[j] ] );
+ if (d > maxD) {
+ maxD=d;
+ selectedPoints={m_extremeValues[i],m_extremeValues[j]};
+ }
+ }
+ }
+ if (EQUAL(maxD, m_epsilonSquared)) {
+ // A degenerate case: the point cloud seems to consists of a
single point
+ return
MeshBuilder<T>(0,std::min((size_t)1,vertexCount-1),std::min((size_t)2,vertexCount-1),std::min((size_t)3,vertexCount-1));
+ }
+ assert(selectedPoints.first != selectedPoints.second);
+
+ // Find the most distant point to the line between the two chosen
extreme points.
+ const Ray<T> r(m_vertexData[selectedPoints.first],
(m_vertexData[selectedPoints.second] - m_vertexData[selectedPoints.first]));
+ maxD = m_epsilonSquared;
+ size_t maxI=std::numeric_limits<size_t>::max();
+ const size_t vCount = m_vertexData.size();
+ for (size_t i=0;i<vCount;i++) {
+ const T distToRay =
mathutils::getSquaredDistanceBetweenPointAndRay(m_vertexData[i],r);
+ if (distToRay > maxD) {
+ maxD=distToRay;
+ maxI=i;
+ }
+ }
+ if (EQUAL(maxD, m_epsilonSquared)) {
+ // It appears that the point cloud belongs to a 1 dimensional
+ // subspace of R^3: convex hull has no volume => return a thin
+ // triangle Pick any point other than selectedPoints.first and
+ // selectedPoints.second as the third point of the triangle
+ auto it =
std::find_if(m_vertexData.begin(),m_vertexData.end(),[&](const vec3& ve) {
+ return ve != m_vertexData[selectedPoints.first] && ve
!= m_vertexData[selectedPoints.second];
+ });
+ const IndexType thirdPoint = (it == m_vertexData.end()) ?
selectedPoints.first : std::distance(m_vertexData.begin(),it);
+ it =
std::find_if(m_vertexData.begin(),m_vertexData.end(),[&](const vec3& ve) {
+ return ve != m_vertexData[selectedPoints.first] && ve
!= m_vertexData[selectedPoints.second] && ve != m_vertexData[thirdPoint];
+ });
+ const IndexType fourthPoint = (it == m_vertexData.end()) ?
selectedPoints.first : std::distance(m_vertexData.begin(),it);
+ return
MeshBuilder<T>(selectedPoints.first,selectedPoints.second,thirdPoint,fourthPoint);
+ }
+
+ // These three points form the base triangle for our tetrahedron.
+ assert(selectedPoints.first != maxI && selectedPoints.second !=
maxI);
+ std::array<size_t,3> baseTriangle{selectedPoints.first,
selectedPoints.second, maxI};
+ const Vector3<T> baseTriangleVertices[]={
m_vertexData[baseTriangle[0]], m_vertexData[baseTriangle[1]],
m_vertexData[baseTriangle[2]] };
+
+ // Next step is to find the 4th vertex of the tetrahedron. We
+ // naturally choose the point farthest away from the triangle
+ // plane.
+ maxD=m_epsilon;
+ maxI=0;
+ const Vector3<T> N =
mathutils::getTriangleNormal(baseTriangleVertices[0],baseTriangleVertices[1],baseTriangleVertices[2]);
+ Plane<T> trianglePlane(N,baseTriangleVertices[0]);
+ for (size_t i=0;i<vCount;i++) {
+ const T d =
std::abs(mathutils::getSignedDistanceToPlane(m_vertexData[i],trianglePlane));
+ if (d>maxD) {
+ maxD=d;
+ maxI=i;
+ }
+ }
+ if (EQUAL(maxD, m_epsilon)) {
+ // All the points seem to lie on a 2D subspace of R^3. How to
handle this? Well, let's add one extra point to the point cloud so that the
convex hull will have volume.
+ m_planar = true;
+ const vec3 N1 =
mathutils::getTriangleNormal(baseTriangleVertices[1],baseTriangleVertices[2],baseTriangleVertices[0]);
+ m_planarPointCloudTemp.clear();
+
m_planarPointCloudTemp.insert(m_planarPointCloudTemp.begin(),m_vertexData.begin(),m_vertexData.end());
+ const vec3 extraPoint = N1 + m_vertexData[0];
+ m_planarPointCloudTemp.push_back(extraPoint);
+ maxI = m_planarPointCloudTemp.size()-1;
+ m_vertexData = VertexDataSource<T>(m_planarPointCloudTemp);
+ }
+
+ // Enforce CCW orientation (if user prefers clockwise orientation,
swap two vertices in each triangle when final mesh is created)
+ const Plane<T> triPlane(N,baseTriangleVertices[0]);
+ if (triPlane.isPointOnPositiveSide(m_vertexData[maxI])) {
+ std::swap(baseTriangle[0],baseTriangle[1]);
+ }
+
+ // Create a tetrahedron half edge mesh and compute planes defined
by each triangle
+ MeshBuilder<T>
mesh(baseTriangle[0],baseTriangle[1],baseTriangle[2],maxI);
+ for (auto& f : mesh.m_faces) {
+ auto v = mesh.getVertexIndicesOfFace(f);
+ const Vector3<T>& va = m_vertexData[v[0]];
+ const Vector3<T>& vb = m_vertexData[v[1]];
+ const Vector3<T>& vc = m_vertexData[v[2]];
+ const Vector3<T> N1 = mathutils::getTriangleNormal(va, vb, vc);
+ const Plane<T> trianglePlane1(N1,va);
+ f.m_P = trianglePlane1;
+ }
+
+ // Finally we assign a face for each vertex outside the tetrahedron
(vertices inside the tetrahedron have no role anymore)
+ for (size_t i=0;i<vCount;i++) {
+ for (auto& face : mesh.m_faces) {
+ if (addPointToFace(face, i)) {
+ break;
+ }
+ }
+ }
+ return mesh;
+ }
+
+ /*
+ * Explicit template specifications for float and double
+ */
+
+ template class QuickHull<float>;
+ template class QuickHull<double>;
+}
+
+// Local Variables:
+// tab-width: 8
+// mode: C++
+// c-basic-offset: 4
+// indent-tabs-mode: t
+// c-file-style: "stroustrup"
+// End:
+// ex: shiftwidth=4 tabstop=8
+
Property changes on: brlcad/trunk/src/libbg/QuickHull.cxx
___________________________________________________________________
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
Added: svn:mime-type
## -0,0 +1 ##
+text/plain
\ No newline at end of property
Added: brlcad/trunk/src/libbg/QuickHull.hpp
===================================================================
--- brlcad/trunk/src/libbg/QuickHull.hpp (rev 0)
+++ brlcad/trunk/src/libbg/QuickHull.hpp 2018-10-03 22:13:09 UTC (rev
71941)
@@ -0,0 +1,1018 @@
+/*
+ * QuickHull.hpp
+ *
+ * Created on: May 4, 2014
+ * Author: anttiku
+ *
+ * Single header version of all header files from the QuickHull implementation
+ * at https://github.com/akuukka/quickhull by Antti Kuukka
+ *
+ * Per the README.md file:
+ *
+ * This implementation is 100% Public Domain.
+ *
+ * Feel free to use.
+ *
+ * C++11 is needed to compile it.
+ */
+
+#ifndef QUICKHULL_HPP_
+#define QUICKHULL_HPP_
+
+#include "common.h"
+
+#include <cmath>
+#include <iostream>
+
+#include <array>
+#include <cassert>
+#include <fstream>
+#include <limits>
+#include <memory>
+#include <string>
+#include <unordered_map>
+#include <vector>
+
+extern "C" {
+#include "vmath.h"
+}
+
+/*
+ * Implementation of the 3D QuickHull algorithm by Antti Kuukka
+ *
+ * No copyrights. What follows is 100% Public Domain.
+ *
+ *
+ *
+ * INPUT: a list of points in 3D space (for example, vertices of a 3D mesh)
+ *
+ * OUTPUT: a ConvexHull object which provides vertex and index buffers of the
+ * generated convex hull as a triangle mesh.
+ *
+ *
+ *
+ * The implementation is thread-safe if each thread is using its own QuickHull
+ * object.
+ *
+ *
+ * SUMMARY OF THE ALGORITHM:
+ * - Create initial simplex (tetrahedron) using extreme points. We
+ * have four faces now and they form a convex mesh M.
+ * - For each point, assign them to the first face for which they
+ * are on the positive side of (so each point is assigned to at most
+ * one face). Points inside the initial tetrahedron are left behind
+ * now and no longer affect the calculations.
+ * - Add all faces that have points assigned to them to Face Stack.
+ * - Iterate until Face Stack is empty:
+ * - Pop topmost face F from the stack
+ * - From the points assigned to F, pick the point P that is
+ * farthest away from the plane defined by F.
+ * - Find all faces of M that have P on their positive side. Let
+ * us call these the "visible faces".
+ * - Because of the way M is constructed, these faces are
+ * connected. Solve their horizon edge loop.
+ * - "Extrude to P": Create new faces by connecting P with the
+ * points belonging to the horizon edge. Add the new faces to
+ * M and remove the visible faces from M.
+ * - Each point that was assigned to visible faces is now assigned
+ * to at most one of the newly created faces.
+ * - Those new faces that have points assigned to them are added
+ * to the top of Face Stack.
+ * - M is now the convex hull.
+ *
+ * TO DO:
+ * - Implement a proper 2D QuickHull and use that to solve the degenerate 2D
+ * case (when all the points lie on the same plane in 3D space).
+ */
+
+
+namespace quickhull {
+
+ typedef size_t IndexType;
+
+ template <typename T>
+ class Vector3
+ {
+ public:
+ Vector3() = default;
+
+ Vector3(T px, T py, T pz) : x(px), y(py), z(pz) {
+
+ }
+
+ T x,y,z;
+
+ T dotProduct(const Vector3& other) const {
+ return x*other.x+y*other.y+z*other.z;
+ }
+
+ void normalize() {
+ const T len = getLength();
+ x/=len;
+ y/=len;
+ z/=len;
+ }
+
+ Vector3 getNormalized() const {
+ const T len = getLength();
+ return Vector3(x/len,y/len,z/len);
+ }
+
+ T getLength() const {
+ return std::sqrt(x*x+y*y+z*z);
+ }
+
+ Vector3 operator-(const Vector3& other) const {
+ return Vector3(x-other.x,y-other.y,z-other.z);
+ }
+
+ Vector3 operator+(const Vector3& other) const {
+ return Vector3(x+other.x,y+other.y,z+other.z);
+ }
+
+ Vector3& operator+=(const Vector3& other) {
+ x+=other.x;
+ y+=other.y;
+ z+=other.z;
+ return *this;
+ }
+ Vector3& operator-=(const Vector3& other) {
+ x-=other.x;
+ y-=other.y;
+ z-=other.z;
+ return *this;
+ }
+ Vector3& operator*=(T c) {
+ x*=c;
+ y*=c;
+ z*=c;
+ return *this;
+ }
+
+ Vector3& operator/=(T c) {
+ x/=c;
+ y/=c;
+ z/=c;
+ return *this;
+ }
+
+ Vector3 operator-() const {
+ return Vector3(-x,-y,-z);
+ }
+
+ template<typename S>
+ Vector3 operator*(S c) const {
+ return Vector3(x*c,y*c,z*c);
+ }
+
+ template<typename S>
+ Vector3 operator/(S c) const {
+ return Vector3(x/c,y/c,z/c);
+ }
+
+ T getLengthSquared() const {
+ return x*x + y*y + z*z;
+ }
+
+ bool operator!=(const Vector3& o) const {
+ return (!EQUAL(x, o.x) || !EQUAL(y, o.y) || !EQUAL(z, o.z));
+ }
+
+ // Projection onto another vector
+ Vector3 projection(const Vector3& o) const {
+ T C = dotProduct(o)/o.getLengthSquared();
+ return o*C;
+ }
+
+ Vector3 crossProduct (const Vector3& rhs ) {
+ T a = y * rhs.z - z * rhs.y ;
+ T b = z * rhs.x - x * rhs.z ;
+ T c = x * rhs.y - y * rhs.x ;
+ Vector3 product( a , b , c ) ;
+ return product ;
+ }
+
+ T getDistanceTo(const Vector3& other) const {
+ Vector3 diff = *this - other;
+ return diff.getLength();
+ }
+
+ T getSquaredDistanceTo(const Vector3& other) const {
+ const T dx = x-other.x;
+ const T dy = y-other.y;
+ const T dz = z-other.z;
+ return dx*dx+dy*dy+dz*dz;
+ }
+
+ };
+
+ // Overload also << operator for easy printing of debug data
+ template <typename T>
+ inline std::ostream& operator<<(std::ostream& os, const Vector3<T>&
vec) {
+ os << "(" << vec.x << "," << vec.y << "," << vec.z << ")";
+ return os;
+ }
+
+ template <typename T>
+ inline Vector3<T> operator*(T c, const Vector3<T>& v) {
+ return Vector3<T>(v.x*c,v.y*c,v.z*c);
+ }
+
+
+ template <typename T>
+ struct Ray {
+ const Vector3<T> m_S;
+ const Vector3<T> m_V;
+ const T m_VInvLengthSquared;
+
+ Ray(const Vector3<T>& S,const Vector3<T>& V) : m_S(S), m_V(V),
m_VInvLengthSquared(1/m_V.getLengthSquared()) {
+ }
+ };
+
+ template<typename T>
+ class Pool {
+ std::vector<std::unique_ptr<T>> m_data;
+ public:
+ void clear() {
+ m_data.clear();
+ }
+
+ void reclaim(std::unique_ptr<T>& ptr) {
+ m_data.push_back(std::move(ptr));
+ }
+
+ std::unique_ptr<T> get() {
+ if (m_data.size()==0) {
+ return std::unique_ptr<T>(new T());
+ }
+ auto it = m_data.end()-1;
+ std::unique_ptr<T> r = std::move(*it);
+ m_data.erase(it);
+ return r;
+ }
+
+ };
+
+ template<typename T>
+ class Plane {
+ public:
+ Vector3<T> m_N;
+
+ // Signed distance (if normal is of length 1) to the plane from
origin
+ T m_D;
+
+ // Normal length squared
+ T m_sqrNLength;
+
+ bool isPointOnPositiveSide(const Vector3<T>& Q) const {
+ T d = m_N.dotProduct(Q)+m_D;
+ if (d>=0) return true;
+ return false;
+ }
+
+ Plane() = default;
+
+ // Construct a plane using normal N and any point P on the plane
+ Plane(const Vector3<T>& N, const Vector3<T>& P) : m_N(N),
m_D(-N.dotProduct(P)), m_sqrNLength(m_N.x*m_N.x+m_N.y*m_N.y+m_N.z*m_N.z) {
+
+ }
+ };
+
+
+ template<typename T>
+ class VertexDataSource {
+ const Vector3<T>* m_ptr;
+ size_t m_count;
+
+ public:
+ VertexDataSource(const Vector3<T>* ptr, size_t count) : m_ptr(ptr),
m_count(count) {
+
+ }
+
+ VertexDataSource(const std::vector<Vector3<T>>& vec) :
m_ptr(&vec[0]), m_count(vec.size()) {
+
+ }
+
+ VertexDataSource() : m_ptr(nullptr), m_count(0) {
+
+ }
+
+ VertexDataSource& operator=(const VertexDataSource& other) =
default;
+
+ size_t size() const {
+ return m_count;
+ }
+
+ const Vector3<T>& operator[](size_t index) const {
+ return m_ptr[index];
+ }
+
+ const Vector3<T>* begin() const {
+ return m_ptr;
+ }
+
+ const Vector3<T>* end() const {
+ return m_ptr + m_count;
+ }
+ };
+
+ template <typename T>
+ class MeshBuilder {
+ public:
+ struct HalfEdge {
+ IndexType m_endVertex;
+ IndexType m_opp;
+ IndexType m_face;
+ IndexType m_next;
+
+ void disable() {
+ m_endVertex = std::numeric_limits<IndexType>::max();
+ }
+
+ bool isDisabled() const {
+ return m_endVertex ==
std::numeric_limits<IndexType>::max();
+ }
+ };
+
+ struct Face {
+ IndexType m_he;
+ Plane<T> m_P;
+ T m_mostDistantPointDist;
+ IndexType m_mostDistantPoint;
+ size_t m_visibilityCheckedOnIteration;
+ std::uint8_t m_isVisibleFaceOnCurrentIteration : 1;
+ std::uint8_t m_inFaceStack : 1;
+ std::uint8_t m_horizonEdgesOnCurrentIteration : 3; // Bit
for each half edge assigned to this face, each being 0 or 1 depending on
whether the edge belongs to horizon edge
+ std::unique_ptr<std::vector<IndexType>>
m_pointsOnPositiveSide;
+
+ Face() : m_he(std::numeric_limits<IndexType>::max()),
+ m_mostDistantPointDist(0),
+ m_mostDistantPoint(0),
+ m_visibilityCheckedOnIteration(0),
+ m_isVisibleFaceOnCurrentIteration(0),
+ m_inFaceStack(0),
+ m_horizonEdgesOnCurrentIteration(0)
+ {
+
+ }
+
+ void disable() {
+ m_he = std::numeric_limits<IndexType>::max();
+ }
+
+ bool isDisabled() const {
+ return m_he == std::numeric_limits<IndexType>::max();
+ }
+ };
+
+ // Mesh data
+ std::vector<Face> m_faces;
+ std::vector<HalfEdge> m_halfEdges;
+
+ // When the mesh is modified and faces and half edges are
removed from it, we do not actually remove them from the container vectors.
+ // Insted, they are marked as disabled which means that the
indices can be reused when we need to add new faces and half edges to the mesh.
+ // We store the free indices in the following vectors.
+ std::vector<IndexType> m_disabledFaces,m_disabledHalfEdges;
+
+ IndexType addFace() {
+ if (m_disabledFaces.size()) {
+ IndexType index = m_disabledFaces.back();
+ auto& f = m_faces[index];
+ assert(f.isDisabled());
+ assert(!f.m_pointsOnPositiveSide);
+ f.m_mostDistantPointDist = 0;
+ m_disabledFaces.pop_back();
+ return index;
+ }
+ m_faces.emplace_back();
+ return m_faces.size()-1;
+ }
+
+ IndexType addHalfEdge() {
+ if (m_disabledHalfEdges.size()) {
+ const IndexType index = m_disabledHalfEdges.back();
+ m_disabledHalfEdges.pop_back();
+ return index;
+ }
+ m_halfEdges.emplace_back();
+ return m_halfEdges.size()-1;
+ }
+
+ // Mark a face as disabled and return a pointer to the points
that were on the positive of it.
+ std::unique_ptr<std::vector<IndexType>> disableFace(IndexType
faceIndex) {
+ auto& f = m_faces[faceIndex];
+ f.disable();
+ m_disabledFaces.push_back(faceIndex);
+ return std::move(f.m_pointsOnPositiveSide);
+ }
+
+ void disableHalfEdge(IndexType heIndex) {
+ auto& he = m_halfEdges[heIndex];
+ he.disable();
+ m_disabledHalfEdges.push_back(heIndex);
+ }
+
+ MeshBuilder() = default;
+
+ // Create a mesh with initial tetrahedron ABCD. Dot product of
AB with the normal of triangle ABC should be negative.
+ MeshBuilder(IndexType a, IndexType b, IndexType c, IndexType d)
{
+ // Create halfedges
+ HalfEdge AB;
+ AB.m_endVertex = b;
+ AB.m_opp = 6;
+ AB.m_face = 0;
+ AB.m_next = 1;
+ m_halfEdges.push_back(AB);
+
+ HalfEdge BC;
+ BC.m_endVertex = c;
+ BC.m_opp = 9;
+ BC.m_face = 0;
+ BC.m_next = 2;
+ m_halfEdges.push_back(BC);
+
+ HalfEdge CA;
+ CA.m_endVertex = a;
+ CA.m_opp = 3;
+ CA.m_face = 0;
+ CA.m_next = 0;
+ m_halfEdges.push_back(CA);
+
+ HalfEdge AC;
+ AC.m_endVertex = c;
+ AC.m_opp = 2;
+ AC.m_face = 1;
+ AC.m_next = 4;
+ m_halfEdges.push_back(AC);
+
+ HalfEdge CD;
+ CD.m_endVertex = d;
+ CD.m_opp = 11;
+ CD.m_face = 1;
+ CD.m_next = 5;
+ m_halfEdges.push_back(CD);
+
+ HalfEdge DA;
+ DA.m_endVertex = a;
+ DA.m_opp = 7;
+ DA.m_face = 1;
+ DA.m_next = 3;
+ m_halfEdges.push_back(DA);
+
+ HalfEdge BA;
+ BA.m_endVertex = a;
+ BA.m_opp = 0;
+ BA.m_face = 2;
+ BA.m_next = 7;
+ m_halfEdges.push_back(BA);
+
+ HalfEdge AD;
+ AD.m_endVertex = d;
+ AD.m_opp = 5;
+ AD.m_face = 2;
+ AD.m_next = 8;
+ m_halfEdges.push_back(AD);
+
+ HalfEdge DB;
+ DB.m_endVertex = b;
+ DB.m_opp = 10;
+ DB.m_face = 2;
+ DB.m_next = 6;
+ m_halfEdges.push_back(DB);
+
+ HalfEdge CB;
+ CB.m_endVertex = b;
+ CB.m_opp = 1;
+ CB.m_face = 3;
+ CB.m_next = 10;
+ m_halfEdges.push_back(CB);
+
+ HalfEdge BD;
+ BD.m_endVertex = d;
+ BD.m_opp = 8;
+ BD.m_face = 3;
+ BD.m_next = 11;
+ m_halfEdges.push_back(BD);
+
+ HalfEdge DC;
+ DC.m_endVertex = c;
+ DC.m_opp = 4;
+ DC.m_face = 3;
+ DC.m_next = 9;
+ m_halfEdges.push_back(DC);
+
+ // Create faces
+ Face ABC;
+ ABC.m_he = 0;
+ m_faces.push_back(std::move(ABC));
+
+ Face ACD;
+ ACD.m_he = 3;
+ m_faces.push_back(std::move(ACD));
+
+ Face BAD;
+ BAD.m_he = 6;
+ m_faces.push_back(std::move(BAD));
+
+ Face CBD;
+ CBD.m_he = 9;
+ m_faces.push_back(std::move(CBD));
+ }
+
+ std::array<IndexType,3> getVertexIndicesOfFace(const Face& f)
const {
+ std::array<IndexType,3> v;
+ const HalfEdge* he = &m_halfEdges[f.m_he];
+ v[0] = he->m_endVertex;
+ he = &m_halfEdges[he->m_next];
+ v[1] = he->m_endVertex;
+ he = &m_halfEdges[he->m_next];
+ v[2] = he->m_endVertex;
+ return v;
+ }
+
+ std::array<IndexType,2> getVertexIndicesOfHalfEdge(const
HalfEdge& he) const {
+ return {m_halfEdges[he.m_opp].m_endVertex,he.m_endVertex};
+ }
+
+ std::array<IndexType,3> getHalfEdgeIndicesOfFace(const Face& f)
const {
+ return
{f.m_he,m_halfEdges[f.m_he].m_next,m_halfEdges[m_halfEdges[f.m_he].m_next].m_next};
+ }
+ };
+
+
+
+
+
+ template<typename FloatType, typename IndexType>
+ class HalfEdgeMesh {
+ public:
+
+ struct HalfEdge {
+ IndexType m_endVertex;
+ IndexType m_opp;
+ IndexType m_face;
+ IndexType m_next;
+ };
+
+ struct Face {
+ IndexType m_halfEdgeIndex; // Index of one of the half
edges of this face
+ };
+
+ std::vector<Vector3<FloatType>> m_vertices;
+ std::vector<Face> m_faces;
+ std::vector<HalfEdge> m_halfEdges;
+
+ HalfEdgeMesh(const MeshBuilder<FloatType>& builderObject, const
VertexDataSource<FloatType>& vertexData )
+ {
+ std::unordered_map<IndexType,IndexType> faceMapping;
+ std::unordered_map<IndexType,IndexType> halfEdgeMapping;
+ std::unordered_map<IndexType, IndexType> vertexMapping;
+
+ size_t i=0;
+ for (const auto& face : builderObject.m_faces) {
+ if (!face.isDisabled()) {
+
m_faces.push_back({static_cast<IndexType>(face.m_he)});
+ faceMapping[i] = m_faces.size()-1;
+
+ const auto heIndices =
builderObject.getHalfEdgeIndicesOfFace(face);
+ for (const auto heIndex : heIndices) {
+ const IndexType vertexIndex =
builderObject.m_halfEdges[heIndex].m_endVertex;
+ if (vertexMapping.count(vertexIndex)==0) {
+
m_vertices.push_back(vertexData[vertexIndex]);
+ vertexMapping[vertexIndex] =
m_vertices.size()-1;
+ }
+ }
+ }
+ i++;
+ }
+
+ i=0;
+ for (const auto& halfEdge : builderObject.m_halfEdges) {
+ if (!halfEdge.isDisabled()) {
+
m_halfEdges.push_back({static_cast<IndexType>(halfEdge.m_endVertex),static_cast<IndexType>(halfEdge.m_opp),static_cast<IndexType>(halfEdge.m_face),static_cast<IndexType>(halfEdge.m_next)});
+ halfEdgeMapping[i] = m_halfEdges.size()-1;
+ }
+ i++;
+ }
+
+ for (auto& face : m_faces) {
+ assert(halfEdgeMapping.count(face.m_halfEdgeIndex) ==
1);
+ face.m_halfEdgeIndex =
halfEdgeMapping[face.m_halfEdgeIndex];
+ }
+
+ for (auto& he : m_halfEdges) {
+ he.m_face = faceMapping[he.m_face];
+ he.m_opp = halfEdgeMapping[he.m_opp];
+ he.m_next = halfEdgeMapping[he.m_next];
+ he.m_endVertex = vertexMapping[he.m_endVertex];
+ }
+ }
+
+ };
+
+
+ namespace mathutils {
+
+ template <typename T>
+ inline T getSquaredDistanceBetweenPointAndRay(const Vector3<T>& p,
const Ray<T>& r) {
+ const Vector3<T> s = p-r.m_S;
+ T t = s.dotProduct(r.m_V);
+ return s.getLengthSquared() - t*t*r.m_VInvLengthSquared;
+ }
+
+ // Note that the unit of distance returned is relative to plane's
normal's length (divide by N.getNormalized() if needed to get the "real"
distance).
+ template <typename T>
+ inline T getSignedDistanceToPlane(const Vector3<T>& v, const
Plane<T>& p) {
+ return p.m_N.dotProduct(v) + p.m_D;
+ }
+
+ template <typename T>
+ inline Vector3<T> getTriangleNormal(const Vector3<T>& a,const
Vector3<T>& b,const Vector3<T>& c) {
+ // We want to get (a-c).crossProduct(b-c) without constructing
temp vectors
+ T x = a.x - c.x;
+ T y = a.y - c.y;
+ T z = a.z - c.z;
+ T rhsx = b.x - c.x;
+ T rhsy = b.y - c.y;
+ T rhsz = b.z - c.z;
+ T px = y * rhsz - z * rhsy ;
+ T py = z * rhsx - x * rhsz ;
+ T pz = x * rhsy - y * rhsx ;
+ return Vector3<T>(px,py,pz);
+ }
+
+
+ }
+
+
+ template<typename T>
+ class ConvexHull {
+ std::unique_ptr<std::vector<Vector3<T>>> m_optimizedVertexBuffer;
+ VertexDataSource<T> m_vertices;
+ std::vector<size_t> m_indices;
+ public:
+ ConvexHull() {}
+
+ // Copy constructor
+ ConvexHull(const ConvexHull& o) {
+ m_indices = o.m_indices;
+ if (o.m_optimizedVertexBuffer) {
+ m_optimizedVertexBuffer.reset(new
std::vector<Vector3<T>>(*o.m_optimizedVertexBuffer));
+ m_vertices = VertexDataSource<T>(*m_optimizedVertexBuffer);
+ }
+ else {
+ m_vertices = o.m_vertices;
+ }
+ }
+
+ ConvexHull& operator=(const ConvexHull& o) {
+ if (&o == this) {
+ return *this;
+ }
+ m_indices = o.m_indices;
+ if (o.m_optimizedVertexBuffer) {
+ m_optimizedVertexBuffer.reset(new
std::vector<Vector3<T>>(*o.m_optimizedVertexBuffer));
+ m_vertices = VertexDataSource<T>(*m_optimizedVertexBuffer);
+ }
+ else {
+ m_vertices = o.m_vertices;
+ }
+ return *this;
+ }
+
+ ConvexHull(ConvexHull&& o) {
+ m_indices = std::move(o.m_indices);
+ if (o.m_optimizedVertexBuffer) {
+ m_optimizedVertexBuffer =
std::move(o.m_optimizedVertexBuffer);
+ o.m_vertices = VertexDataSource<T>();
+ m_vertices = VertexDataSource<T>(*m_optimizedVertexBuffer);
+ }
+ else {
+ m_vertices = o.m_vertices;
+ }
+ }
+
+ ConvexHull& operator=(ConvexHull&& o) {
+ if (&o == this) {
+ return *this;
+ }
+ m_indices = std::move(o.m_indices);
+ if (o.m_optimizedVertexBuffer) {
+ m_optimizedVertexBuffer =
std::move(o.m_optimizedVertexBuffer);
+ o.m_vertices = VertexDataSource<T>();
+ m_vertices = VertexDataSource<T>(*m_optimizedVertexBuffer);
+ }
+ else {
+ m_vertices = o.m_vertices;
+ }
+ return *this;
+ }
+
+ // Construct vertex and index buffers from half edge mesh and
pointcloud
+ ConvexHull(const MeshBuilder<T>& mesh, const VertexDataSource<T>&
pointCloud, bool CCW, bool useOriginalIndices) {
+ if (!useOriginalIndices) {
+ m_optimizedVertexBuffer.reset(new
std::vector<Vector3<T>>());
+ }
+
+ std::vector<bool> faceProcessed(mesh.m_faces.size(),false);
+ std::vector<size_t> faceStack;
+ std::unordered_map<size_t,size_t> vertexIndexMapping; // Map
vertex indices from original point cloud to the new mesh vertex indices
+ for (size_t i = 0;i<mesh.m_faces.size();i++) {
+ if (!mesh.m_faces[i].isDisabled()) {
+ faceStack.push_back(i);
+ break;
+ }
+ }
+ if (faceStack.size()==0) {
+ return;
+ }
+
+ const size_t finalMeshFaceCount = mesh.m_faces.size() -
mesh.m_disabledFaces.size();
+ m_indices.reserve(finalMeshFaceCount*3);
+
+ while (faceStack.size()) {
+ auto it = faceStack.end()-1;
+ size_t top = *it;
+ assert(!mesh.m_faces[top].isDisabled());
+ faceStack.erase(it);
+ if (faceProcessed[top]) {
+ continue;
+ }
+ else {
+ faceProcessed[top]=true;
+ auto halfEdges =
mesh.getHalfEdgeIndicesOfFace(mesh.m_faces[top]);
+ size_t adjacent[] =
{mesh.m_halfEdges[mesh.m_halfEdges[halfEdges[0]].m_opp].m_face,mesh.m_halfEdges[mesh.m_halfEdges[halfEdges[1]].m_opp].m_face,mesh.m_halfEdges[mesh.m_halfEdges[halfEdges[2]].m_opp].m_face};
+ for (auto a : adjacent) {
+ if (!faceProcessed[a] &&
!mesh.m_faces[a].isDisabled()) {
+ faceStack.push_back(a);
+ }
+ }
+ auto vertices =
mesh.getVertexIndicesOfFace(mesh.m_faces[top]);
+ if (!useOriginalIndices) {
+ for (auto& v : vertices) {
+ auto lit = vertexIndexMapping.find(v);
+ if (lit == vertexIndexMapping.end()) {
+
m_optimizedVertexBuffer->push_back(pointCloud[v]);
+ vertexIndexMapping[v] =
m_optimizedVertexBuffer->size()-1;
+ v = m_optimizedVertexBuffer->size()-1;
+ }
+ else {
+ v = lit->second;
+ }
+ }
+ }
+ m_indices.push_back(vertices[0]);
+ if (CCW) {
+ m_indices.push_back(vertices[2]);
+ m_indices.push_back(vertices[1]);
+ }
+ else {
+ m_indices.push_back(vertices[1]);
+ m_indices.push_back(vertices[2]);
+ }
+ }
+ }
+
+ if (!useOriginalIndices) {
+ m_vertices = VertexDataSource<T>(*m_optimizedVertexBuffer);
+ }
+ else {
+ m_vertices = pointCloud;
+ }
+ }
+
+ std::vector<size_t>& getIndexBuffer() {
+ return m_indices;
+ }
+
+ VertexDataSource<T>& getVertexBuffer() {
+ return m_vertices;
+ }
+
+ // Export the mesh to a Waveform OBJ file
+ void writeWaveformOBJ(const std::string& filename, const
std::string& objectName = "quickhull")
+ {
+ std::ofstream objFile;
+ objFile.open (filename);
+ objFile << "o " << objectName << "\n";
+ for (const auto& v : getVertexBuffer()) {
+ objFile << "v " << v.x << " " << v.y << " " << v.z << "\n";
+ }
+ const auto& indBuf = getIndexBuffer();
+ size_t triangleCount = indBuf.size()/3;
+ for (size_t i=0;i<triangleCount;i++) {
+ objFile << "f " << indBuf[i*3]+1 << " " << indBuf[i*3+1]+1
<< " " << indBuf[i*3+2]+1 << "\n";
+ }
+ objFile.close();
+ }
+
+ };
+
+
+
+ struct DiagnosticsData {
+ size_t m_failedHorizonEdges; // How many times QuickHull failed to
solve the horizon edge. Failures lead to degenerated convex hulls.
+
+ DiagnosticsData() : m_failedHorizonEdges(0) { }
+ };
+
+ template<typename FloatType>
+ class QuickHull {
+ using vec3 = Vector3<FloatType>;
+
+ static const FloatType Epsilon;
+
+ FloatType m_epsilon, m_epsilonSquared, m_scale;
+ bool m_planar;
+ std::vector<vec3> m_planarPointCloudTemp;
+ VertexDataSource<FloatType> m_vertexData;
+ MeshBuilder<FloatType> m_mesh;
+ std::array<IndexType,6> m_extremeValues;
+ DiagnosticsData m_diagnostics;
+
+ // Temporary variables used during iteration process
+ std::vector<IndexType> m_newFaceIndices;
+ std::vector<IndexType> m_newHalfEdgeIndices;
+ std::vector< std::unique_ptr<std::vector<IndexType>> >
m_disabledFacePointVectors;
+
+ // Create a half edge mesh representing the base tetrahedron from
+ // which the QuickHull iteration proceeds. m_extremeValues must be
+ // properly set up when this is called.
+ MeshBuilder<FloatType> getInitialTetrahedron();
+
+ // Given a list of half edges, try to rearrange them so that they
+ // form a loop. Return true on success.
+ bool reorderHorizonEdges(std::vector<IndexType>& horizonEdges);
+
+ // Find indices of extreme values (max x, min x, max y, min y, max
+ // z, min z) for the given point cloud
+ std::array<IndexType,6> getExtremeValues();
+
+ // Compute scale of the vertex data.
+ FloatType getScale(const std::array<IndexType,6>& extremeValues);
+
+ // Each face contains a unique pointer to a vector of indices.
+ // However, many - often most - faces do not have any points on the
+ // positive side of them especially at the the end of the
+ // iteration. When a face is removed from the mesh, its associated
+ // point vector, if such exists, is moved to the index vector pool,
+ // and when we need to add new faces with points on the positive
+ // side to the mesh, we reuse these vectors. This reduces the
+ // amount of std::vectors we have to deal with, and impact on
+ // performance is remarkable.
+ Pool<std::vector<IndexType>> m_indexVectorPool;
+ inline std::unique_ptr<std::vector<IndexType>>
getIndexVectorFromPool();
+ inline void
reclaimToIndexVectorPool(std::unique_ptr<std::vector<IndexType>>& ptr);
+
+ // Associates a point with a face if the point resides on the
+ // positive side of the plane. Returns true if the points was on
+ // the positive side.
+ inline bool addPointToFace(typename MeshBuilder<FloatType>::Face&
f, IndexType pointIndex);
+
+ // This will update m_mesh from which we create the ConvexHull
+ // object that getConvexHull function returns
+ void createConvexHalfEdgeMesh();
+
+ // Constructs the convex hull into a MeshBuilder object which can
+ // be converted to a ConvexHull or Mesh object
+ void buildMesh(const VertexDataSource<FloatType>& pointCloud, bool
UNUSED(CCW), bool UNUSED(useOriginalIndices), FloatType eps);
+
+ // The public getConvexHull functions will setup a VertexDataSource
+ // object and call this
+ ConvexHull<FloatType> getConvexHull(const
VertexDataSource<FloatType>& pointCloud, bool CCW, bool useOriginalIndices,
FloatType eps);
+ public:
+ // Computes convex hull for a given point cloud.
+ // Params:
+ // pointCloud: a vector of of 3D points
+ //
+ // CCW: whether the output mesh triangles should have CCW
+ // orientation
+ //
+ // useOriginalIndices: should the output mesh use same vertex
+ // indices as the original point cloud. If this is false, then
+ // we generate a new vertex buffer which contains only the
+ // vertices that are part of the convex hull.
+ //
+ // eps: minimum distance to a plane to consider a point being on
+ // positive of it (for a point cloud with scale 1)
+ ConvexHull<FloatType> getConvexHull(const
std::vector<Vector3<FloatType>>& pointCloud, bool CCW, bool useOriginalIndices,
FloatType eps = Epsilon);
+
+ // Computes convex hull for a given point cloud.
+ // Params:
+ //
+ // vertexData: pointer to the first 3D point of the point cloud
+ //
+ // vertexCount: number of vertices in the point cloud
+ //
+ // CCW: whether the output mesh triangles should have CCW
+ // orientation
+ //
+ // useOriginalIndices: should the output mesh use same vertex
+ // indices as the original point cloud. If this is false, then
+ // we generate a new vertex buffer which contains only the
+ // vertices that are part of the convex hull.
+ //
+ // eps: minimum distance to a plane to consider a point being on
+ // positive of it (for a point cloud with scale 1)
+ ConvexHull<FloatType> getConvexHull(const Vector3<FloatType>*
vertexData, size_t vertexCount, bool CCW, bool useOriginalIndices, FloatType
eps = Epsilon);
+
+ // Computes convex hull for a given point cloud. This
+ // function assumes that the vertex data resides in
+ // memory in the following format:
+ // x_0,y_0,z_0,x_1,y_1,z_1,...
+ // Params:
+ // vertexData: pointer to the X component of the first point of
+ // the point cloud.
+ //
+ // vertexCount: number of vertices in the point cloud
+ //
+ // CCW: whether the output mesh triangles should have CCW
+ // orientation
+ //
+ // useOriginalIndices: should the output mesh use same vertex
+ // indices as the original point cloud. If this is false, then
+ // we generate a new vertex buffer which contains only the
+ // vertices that are part of the convex hull.
+ //
+ // eps: minimum distance to a plane to consider a point being on
+ // positive of it (for a point cloud with scale 1)
+ ConvexHull<FloatType> getConvexHull(const FloatType* vertexData,
size_t vertexCount, bool CCW, bool useOriginalIndices, FloatType eps = Epsilon);
+
+ // Computes convex hull for a given point cloud. This function
+ // assumes that the vertex data resides in memory in the following
+ // format: x_0,y_0,z_0,x_1,y_1,z_1,...
+ // Params:
+ // vertexData: pointer to the X component of the first point of
the point cloud.
+ //
+ // vertexCount: number of vertices in the point cloud
+ //
+ // CCW: whether the output mesh triangles should have CCW
orientation
+ //
+ // eps: minimum distance to a plane to consider a point being on
+ // positive of it (for a point cloud with scale 1)
+ //
+ // Returns:
+ // Convex hull of the point cloud as a mesh object with half edge
structure.
+ HalfEdgeMesh<FloatType, IndexType> getConvexHullAsMesh(const
FloatType* vertexData, size_t vertexCount, bool CCW, FloatType eps = Epsilon);
+
+ // Get diagnostics about last generated convex hull
+ const DiagnosticsData& getDiagnostics() {
+ return m_diagnostics;
+ }
+ };
+
+ /*
+ * Inline function definitions
+ */
+
+ template<typename T>
+ std::unique_ptr<std::vector<IndexType>>
QuickHull<T>::getIndexVectorFromPool() {
+ auto r = std::move(m_indexVectorPool.get());
+ r->clear();
+ return r;
+ }
+
+ template<typename T>
+ void
QuickHull<T>::reclaimToIndexVectorPool(std::unique_ptr<std::vector<IndexType>>&
ptr) {
+ const size_t oldSize = ptr->size();
+ if ((oldSize+1)*128 < ptr->capacity()) {
+ // Reduce memory usage! Huge vectors are needed at the
+ // beginning of iteration when faces have many points on their
+ // positive side. Later on, smaller vectors will suffice.
+ ptr.reset(nullptr);
+ return;
+ }
+ m_indexVectorPool.reclaim(ptr);
+ }
+
+ template<typename T>
+ bool QuickHull<T>::addPointToFace(typename MeshBuilder<T>::Face& f,
IndexType pointIndex) {
+ const T D = mathutils::getSignedDistanceToPlane(m_vertexData[
pointIndex ],f.m_P);
+ if (D>0 && D*D > m_epsilonSquared*f.m_P.m_sqrNLength) {
+ if (!f.m_pointsOnPositiveSide) {
+ f.m_pointsOnPositiveSide =
std::move(getIndexVectorFromPool());
+ }
+ f.m_pointsOnPositiveSide->push_back( pointIndex );
+ if (D > f.m_mostDistantPointDist) {
+ f.m_mostDistantPointDist = D;
+ f.m_mostDistantPoint = pointIndex;
+ }
+ return true;
+ }
+ return false;
+ }
+
+}
+
+
+#endif /* QUICKHULL_HPP_ */
+
+// Local Variables:
+// tab-width: 8
+// mode: C++
+// c-basic-offset: 4
+// indent-tabs-mode: t
+// c-file-style: "stroustrup"
+// End:
+// ex: shiftwidth=4 tabstop=8
+
Property changes on: brlcad/trunk/src/libbg/QuickHull.hpp
___________________________________________________________________
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
Added: svn:mime-type
## -0,0 +1 ##
+text/plain
\ No newline at end of property
Deleted: brlcad/trunk/src/libbg/chull3d.cpp
===================================================================
--- brlcad/trunk/src/libbg/chull3d.cpp 2018-10-03 13:53:28 UTC (rev 71940)
+++ brlcad/trunk/src/libbg/chull3d.cpp 2018-10-03 22:13:09 UTC (rev 71941)
@@ -1,1288 +0,0 @@
-/*
- * Ken Clarkson wrote this. Copyright (c) 1995 by AT&T..
- * Permission to use, copy, modify, and distribute this software for any
- * purpose without fee is hereby granted, provided that this entire notice
- * is included in all copies of any software which is or includes a copy
- * or modification of this software and in all copies of the supporting
- * documentation for such software.
- * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
- * WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
- * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
- * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
- */
-
-#include "common.h"
-
-#include <map>
-#include <cmath>
-
-#include <float.h>
-#include <locale.h>
-#include <stdarg.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <ctype.h>
-#include <time.h>
-
-/*#include "chull3d.h"*/
-#include "bu/log.h"
-#include "bu/ptbl.h"
-#include "bu/malloc.h"
-#include "bn/randmt.h"
-#include "bg/chull.h"
-
-#define BLOCKSIZE 100000
-#define MAXDIM 3
-
-/* It is faster, when dealing with large point sets, to
- * repeatedly calculate the convex hull for subsets of
- * the total set and then fine the hull of that set of
- * hulls. These parameters control that process - see
- * chull3d_intermediate_set and bg_3d_chull. In the
- * worst case this approach will slow the calculation,
- * but so long as some significant percentage of the
- * bot's vertices are not part of the hull it should
- * help */
-#define CHULL3D_COUNT_THRESHOLD 1000
-#define CHULL3D_DELTA_THRESHOLD 200
-#define CHULL3D_SUBSET_SIZE 275
-
-typedef double Coord;
-typedef Coord* point;
-typedef point site;
-typedef Coord* normalp;
-typedef Coord site_struct;
-
-typedef struct basis_s {
- struct basis_s *next; /* free list */
- int ref_count; /* storage management */
- int lscale; /* the log base 2 of total scaling of vector */
- Coord sqa, sqb; /* sums of squared norms of a part and b part */
- Coord vecs[1]; /* the actual vectors, extended by malloc'ing bigger
*/
-} basis_s;
-
-typedef struct neighbor {
- site vert; /* vertex of simplex */
- struct simplex *simp; /* neighbor sharing all vertices but vert */
- basis_s *basis; /* derived vectors */
-} neighbor;
-
-typedef struct simplex {
- struct simplex *next; /* free list */
- long visit; /* number of last site visiting this simplex */
- short mark;
- basis_s* normal; /* normal vector pointing inward */
- neighbor peak; /* if null, remaining vertices give facet */
- neighbor neigh[1]; /* neighbors of simplex */
-} simplex;
-
-typedef struct fg_node fg;
-typedef struct tree_node Tree;
-struct tree_node {
- Tree *left, *right;
- site key;
- int size; /* maintained to be the number of nodes rooted here
*/
- fg *fgs;
- Tree *next; /* freelist */
-};
-
-typedef struct fg_node {
- Tree *facets;
- double dist, vol; /* of Voronoi face dual to this */
- fg *next; /* freelist */
- short mark;
- int ref_count;
-} fg_node;
-
-
-typedef site gsitef(void *);
-typedef long site_n(void *, site);
-
-/* We use this instead of globals */
-struct chull3d_data {
- size_t simplex_size;
- simplex *simplex_list;
- size_t basis_s_size;
- basis_s *basis_s_list;
- size_t Tree_size;
- Tree *Tree_list;
- int pdim; /* point dimension */
- point *site_blocks;
- int num_blocks;
- simplex *ch_root;
- double Huge;
- int basis_vec_size;
- int exact_bits;
- float b_err_min;
- float b_err_min_sq;
- basis_s *tt_basis;
- basis_s *infinity_basis;
- Coord *hull_infinity;
- int *B;
- gsitef *get_site;
- site_n *site_num;
- fg *faces_gr_t;
- double mult_up;
- Coord *mins;
- Coord *maxs;
- site p; /* the current site */
- long pnum;
- int rdim; /* region dimension: (max) number of sites specifying
region */
- int cdim; /* number of sites currently specifying region */
- int site_size; /* size of malloc needed for a site */
- int point_size; /* size of malloc needed for a point */
- short *mi;
- short *mo;
- long *shufmat;
-
- /* these were static variables in functions */
- simplex **simplex_block_table;
- int num_simplex_blocks;
- basis_s **basis_s_block_table;
- int num_basis_s_blocks;
- Tree **Tree_block_table;
- int num_Tree_blocks;
- int lscale;
- double max_scale;
- double ldetbound;
- double Sb;
- neighbor p_neigh;
- basis_s *b;
- int fg_vn;
- long vnum;
- simplex *ns;
- simplex **st;
- long ss;
- long s_num;
- std::map<long *, int> *point_lookup;
- int free_point_lookup;
-
- /* libbn style inputs and outputs */
- const point_t *input_vert_array;
- int input_vert_cnt;
- struct bu_ptbl *output_pnts;
- int free_output_pnts;
- int next_vert;
- point_t center_pnt;
-
- /* Output containers */
- point_t *vert_array;
- int *vert_cnt;
- int *faces;
- int *face_cnt;
-};
-
-
-#define CHULL3D_VA(x) ((x)->vecs+cdata->rdim)
-#define CHULL3D_VB(x) ((x)->vecs)
-#define CHULL3D_trans(z,p,q) {int ti; for (ti=0;ti<cdata->pdim;ti++)
z[ti+cdata->rdim] = z[ti] = p[ti] - q[ti];}
-#define CHULL3D_DELIFT 0
-
-#define CHULL3D_push(x) *(cdata->st+tms++) = x;
-#define CHULL3D_pop(x) x = *(cdata->st + --tms);
-
-#define CHULL3D_NEWL(cdata, list, X, p) \
-{ \
- p = list ? list : chull3d_new_block_##X(cdata, 1); \
- /*assert(p);*/ \
- list = p->next; \
-}
-
-#define CHULL3D_NEWLRC(cdata, list, X, p) \
-{ \
- p = list ? list : chull3d_new_block_##X(cdata, 1); \
- /*assert(p);*/ \
- list = p->next; \
- p->ref_count = 1; \
-}
-
-#define CHULL3D_NULLIFY(cdata, X,v) { \
- if ((v) && --(v)->ref_count == 0) { \
- memset((v),0,cdata->X##_size); \
- (v)->next = cdata->X##_list; \
- cdata->X##_list = v; \
- } \
- v = NULL;\
-}
-
-#define CHULL3D_copy_simp(cdata, cnew, s) \
-{ \
- int mi; \
- neighbor *mrsn; \
- CHULL3D_NEWL(cdata, cdata->simplex_list, simplex, cnew); \
- memcpy(cnew,s,cdata->simplex_size); \
- for (mi=-1,mrsn=s->neigh-1;mi<(cdata->cdim);mi++,mrsn++) \
- if (mrsn->basis) mrsn->basis->ref_count++; \
-}
-
-HIDDEN simplex *
-chull3d_new_block_simplex(struct chull3d_data *cdata, int make_blocks)
-{
- int i;
- simplex *xlm, *xbt;
- if (make_blocks && cdata) {
- xbt = cdata->simplex_block_table[(cdata->num_simplex_blocks)++] =
(simplex*)malloc(cdata->input_vert_cnt * cdata->simplex_size);
- memset(xbt,0,cdata->input_vert_cnt * cdata->simplex_size);
- xlm = ((simplex*) ( (char*)xbt + (cdata->input_vert_cnt) *
cdata->simplex_size));
- for (i=0; i<cdata->input_vert_cnt; i++) {
- xlm = ((simplex*) ( (char*)xlm - cdata->simplex_size ));
- xlm->next = cdata->simplex_list;
- cdata->simplex_list = xlm;
- }
- return cdata->simplex_list;
- }
- for (i=0; i<cdata->num_simplex_blocks; i++)
free(cdata->simplex_block_table[i]);
- cdata->num_simplex_blocks = 0;
- if (cdata) cdata->simplex_list = 0;
- return 0;
-}
-
-HIDDEN void
-chull3d_free_simplex_storage(struct chull3d_data *cdata)
-{
- chull3d_new_block_simplex(cdata, 0);
-}
-
-HIDDEN basis_s *
-chull3d_new_block_basis_s(struct chull3d_data *cdata, int make_blocks)
-{
- int i;
- basis_s *xlm, *xbt;
- if (make_blocks && cdata) {
- xbt = cdata->basis_s_block_table[(cdata->num_basis_s_blocks)++] =
(basis_s*)malloc(cdata->input_vert_cnt * cdata->basis_s_size);
- memset(xbt,0,cdata->input_vert_cnt * cdata->basis_s_size);
- xlm = ((basis_s*) ( (char*)xbt + (cdata->input_vert_cnt) *
cdata->basis_s_size));
- for (i=0; i<cdata->input_vert_cnt; i++) {
- xlm = ((basis_s*) ( (char*)xlm - cdata->basis_s_size ));
- xlm->next = cdata->basis_s_list;
- cdata->basis_s_list = xlm;
- }
- return cdata->basis_s_list;
- }
- for (i=0; i<cdata->num_basis_s_blocks; i++)
free(cdata->basis_s_block_table[i]);
- cdata->num_basis_s_blocks = 0;
- if (cdata) cdata->basis_s_list = 0;
- return 0;
-}
-
-HIDDEN void
-chull3d_free_basis_s_storage(struct chull3d_data *cdata)
-{
- chull3d_new_block_basis_s(cdata, 0);
-}
-
-
-HIDDEN Tree *
-chull3d_new_block_Tree(struct chull3d_data *cdata, int make_blocks) {
- int i;
- Tree *xlm, *xbt;
- if (make_blocks && cdata) {
- xbt = cdata->Tree_block_table[cdata->num_Tree_blocks++] =
(Tree*)malloc(cdata->input_vert_cnt * cdata->Tree_size);
- memset(xbt,0,cdata->input_vert_cnt * cdata->Tree_size);
- xlm = ((Tree*) ( (char*)xbt + (cdata->input_vert_cnt) *
cdata->Tree_size));
- for (i=0; i<cdata->input_vert_cnt; i++) {
- xlm = ((Tree*) ( (char*)xlm - cdata->Tree_size ));
- xlm->next = cdata->Tree_list;
- cdata->Tree_list = xlm;
- }
- return cdata->Tree_list;
- }
- for (i=0; i<cdata->num_Tree_blocks; i++) free(cdata->Tree_block_table[i]);
- cdata->num_Tree_blocks = 0;
- if (cdata) cdata->Tree_list = 0;
- return 0;
-}
-
-HIDDEN void
-chull3d_free_Tree_storage(struct chull3d_data *cdata)
-{
- chull3d_new_block_Tree(cdata, 0);
-}
-
-/* ch.c : numerical functions for hull computation */
-
-HIDDEN Coord
-chull3d_Vec_dot(struct chull3d_data *cdata, point x, point y)
-{
- int i;
- Coord sum = 0;
- for (i=0;i<cdata->rdim;i++) sum += x[i] * y[i];
- return sum;
-}
-
-HIDDEN Coord
-chull3d_Norm2(struct chull3d_data *cdata, point x)
-{
- int i;
- Coord sum = 0;
- for (i=0;i<cdata->rdim;i++) sum += x[i] * x[i];
- return sum;
-}
-
-HIDDEN void
-chull3d_Ax_plus_y(struct chull3d_data *cdata, Coord a, point x, point y)
-{
- int i;
- for (i=0;i<cdata->rdim;i++) {
- *y++ += a * *x++;
- }
-}
-
-HIDDEN void
-chull3d_Vec_scale(struct chull3d_data *UNUSED(cdata), int n, Coord a, Coord *x)
-{
- Coord *xx = x;
- Coord *xend = xx + n;
-
- while (xx!=xend) {
- *xx *= a;
- xx++;
- }
-}
-
-
-#ifndef HAVE_CXX_LOGB
-double logb(double x) {
- if (x<=0) return -1e305;
- return log(x)/log(2.);
-}
-#endif
-
-
-/* amount by which to scale up vector, for reduce_inner */
-HIDDEN double
-chull3d_sc(struct chull3d_data *cdata, basis_s *v,simplex *s, int k, int j)
-{
- double labound;
-
- if (j<10) {
- labound = logb(v->sqa)/2;
- cdata->max_scale = cdata->exact_bits - labound - 0.66*(k-2)-1
-CHULL3D_DELIFT;
- if (cdata->max_scale<1) {
- cdata->max_scale = 1;
- }
-
- if (j==0) {
- int i;
- neighbor *sni;
- basis_s *snib;
-
- cdata->ldetbound = CHULL3D_DELIFT;
-
- cdata->Sb = 0;
- for (i=k-1,sni=s->neigh+k-1;i>0;i--,sni--) {
- snib = sni->basis;
- cdata->Sb += snib->sqb;
- cdata->ldetbound += logb(snib->sqb)/2 + 1;
- cdata->ldetbound -= snib->lscale;
- }
- }
- }
- if (cdata->ldetbound - v->lscale + logb(v->sqb)/2 + 1 < 0) {
- return 0;
- } else {
- double dlscale = floor(logb(2*cdata->Sb/(v->sqb +
v->sqa*cdata->b_err_min)))/2;
- cdata->lscale = (int)dlscale;
- if (cdata->lscale > cdata->max_scale) {
- dlscale = floor(cdata->max_scale);
- cdata->lscale = (int)dlscale;
- } else if (cdata->lscale<0) cdata->lscale = 0;
- v->lscale += cdata->lscale;
- return ( ((cdata->lscale)<20) ? 1<<(cdata->lscale) :
ldexp(1.,(cdata->lscale)) );
- }
-}
-
-
-HIDDEN int
-chull3d_reduce_inner(struct chull3d_data *cdata, basis_s *v, simplex *s, int k)
-{
- point va = CHULL3D_VA(v);
- point vb = CHULL3D_VB(v);
- int i,j;
- double dd;
- double scale;
- basis_s *snibv;
- neighbor *sni;
-
- v->sqa = v->sqb = chull3d_Norm2(cdata, vb);
- if (k<=1) {
- memcpy(vb,va,cdata->basis_vec_size);
- return 1;
- }
- for (j=0;j<250;j++) {
-
- memcpy(vb,va,cdata->basis_vec_size);
- for (i=k-1,sni=s->neigh+k-1;i>0;i--,sni--) {
- snibv = sni->basis;
- dd = -chull3d_Vec_dot(cdata, CHULL3D_VB(snibv),vb)/ snibv->sqb;
- chull3d_Ax_plus_y(cdata, dd, CHULL3D_VA(snibv), vb);
- }
- v->sqb = chull3d_Norm2(cdata, vb);
- v->sqa = chull3d_Norm2(cdata, va);
-
- if (2*v->sqb >= v->sqa) {cdata->B[j]++; return 1;}
-
- chull3d_Vec_scale(cdata, cdata->rdim,scale =
chull3d_sc(cdata,v,s,k,j),va);
-
- for (i=k-1,sni=s->neigh+k-1;i>0;i--,sni--) {
- snibv = sni->basis;
- dd = -chull3d_Vec_dot(cdata, CHULL3D_VB(snibv),va)/snibv->sqb;
- dd = floor(dd+0.5);
- chull3d_Ax_plus_y(cdata, dd, CHULL3D_VA(snibv), va);
- }
- }
- return 0;
-}
-
-HIDDEN int
-chull3d_reduce(struct chull3d_data *cdata, basis_s **v, point p, simplex *s,
int k)
-{
- point z;
- point tt = s->neigh[0].vert;
-
- if (!*v) CHULL3D_NEWLRC(cdata, cdata->basis_s_list,basis_s,(*v))
- else (*v)->lscale = 0;
- z = CHULL3D_VB(*v);
- CHULL3D_trans(z,p,tt);
- return chull3d_reduce_inner(cdata,*v,s,k);
-}
-
-HIDDEN void
-chull3d_get_basis_sede(struct chull3d_data *cdata, simplex *s)
-{
- int k=1;
- neighbor *sn = s->neigh+1,
- *sn0 = s->neigh;
- if (!sn0->basis) {
- sn0->basis = cdata->tt_basis;
- cdata->tt_basis->ref_count++;
- } else while (k < (cdata->cdim) && sn->basis) {k++;sn++;}
- while (k<(cdata->cdim)) {
- CHULL3D_NULLIFY(cdata, basis_s,sn->basis);
- chull3d_reduce(cdata, &sn->basis,sn->vert,s,k);
- k++; sn++;
- }
-}
-
-HIDDEN int
-chull3d_out_of_flat(struct chull3d_data *cdata, simplex *root, point p)
-{
- if (!cdata->p_neigh.basis) cdata->p_neigh.basis = (basis_s*)
malloc(cdata->basis_s_size);
- cdata->p_neigh.vert = p;
- (cdata->cdim)++;
- root->neigh[(cdata->cdim)-1].vert = root->peak.vert;
- CHULL3D_NULLIFY(cdata, basis_s,root->neigh[(cdata->cdim)-1].basis);
- chull3d_get_basis_sede(cdata, root);
- chull3d_reduce(cdata, &cdata->p_neigh.basis,p,root,(cdata->cdim));
- /*if (cdata->p_neigh.basis->sqa != 0) return 1;*/
- if (!ZERO(cdata->p_neigh.basis->sqa)) return 1;
- (cdata->cdim)--;
- return 0;
-}
-
-HIDDEN int chull3d_sees(struct chull3d_data *, site, simplex *);
-
-HIDDEN void
-chull3d_get_normal_sede(struct chull3d_data *cdata, simplex *s)
-{
- neighbor *rn;
- int i,j;
-
- chull3d_get_basis_sede(cdata, s);
- if (cdata->rdim==3 && (cdata->cdim)==3) {
- point c;
- point a = CHULL3D_VB(s->neigh[1].basis);
- point b = CHULL3D_VB(s->neigh[2].basis);
- CHULL3D_NEWLRC(cdata, cdata->basis_s_list, basis_s,s->normal);
- c = CHULL3D_VB(s->normal);
- c[0] = a[1]*b[2] - a[2]*b[1];
- c[1] = a[2]*b[0] - a[0]*b[2];
- c[2] = a[0]*b[1] - a[1]*b[0];
- s->normal->sqb = chull3d_Norm2(cdata, c);
- for (i=(cdata->cdim)+1,rn = cdata->ch_root->neigh+(cdata->cdim)-1; i;
i--, rn--) {
- for (j = 0; j<(cdata->cdim) && rn->vert != s->neigh[j].vert;j++) {};
- if (j<(cdata->cdim)) continue;
- if (rn->vert==cdata->hull_infinity) {
- if (c[2] > -cdata->b_err_min) continue;
- } else if (!chull3d_sees(cdata, rn->vert,s)) continue;
- c[0] = -c[0]; c[1] = -c[1]; c[2] = -c[2];
- break;
- }
- return;
- }
-
- for (i=(cdata->cdim)+1,rn = cdata->ch_root->neigh+(cdata->cdim)-1; i; i--,
rn--) {
- for (j = 0; j<(cdata->cdim) && rn->vert != s->neigh[j].vert;j++) {};
- if (j<(cdata->cdim)) continue;
- chull3d_reduce(cdata, &s->normal,rn->vert,s,(cdata->cdim));
- if (!ZERO(s->normal->sqb)) break;
- }
-}
-
-
-//HIDDEN void chull3d_get_normal(struct chull3d_data *cdata, simplex *s)
{chull3d_get_normal_sede(cdata,s); return;}
-
-HIDDEN int
-chull3d_sees(struct chull3d_data *cdata, site p, simplex *s) {
-
- point tt,zz;
- double dd,dds;
- int i;
-
- if (!cdata->b) cdata->b = (basis_s*)malloc(cdata->basis_s_size);
- else cdata->b->lscale = 0;
- zz = CHULL3D_VB(cdata->b);
- if ((cdata->cdim)==0) return 0;
- if (!s->normal) {
- chull3d_get_normal_sede(cdata, s);
- for (i=0;i<(cdata->cdim);i++) CHULL3D_NULLIFY(cdata,
basis_s,s->neigh[i].basis);
- }
- tt = s->neigh[0].vert;
- CHULL3D_trans(zz,p,tt);
- for (i=0;i<3;i++) {
- dd = chull3d_Vec_dot(cdata, zz,s->normal->vecs);
- if (ZERO(dd)) {
- return 0;
- }
- dds = dd*dd/s->normal->sqb/chull3d_Norm2(cdata, zz);
- if (dds > cdata->b_err_min_sq) return (dd<0);
- chull3d_get_basis_sede(cdata, s);
- chull3d_reduce_inner(cdata,cdata->b,s,(cdata->cdim));
- }
- return 0;
-}
-
-HIDDEN void
-chull3d_free_hull_storage(struct chull3d_data *cdata)
-{
- chull3d_free_basis_s_storage(cdata);
- chull3d_free_simplex_storage(cdata);
- chull3d_free_Tree_storage(cdata);
-}
-
-
-
-/* hull.c : "combinatorial" functions for hull computation */
-
-typedef int chull3d_test_func(struct chull3d_data *, simplex *, int, void *);
-typedef void* chull3d_visit_func(struct chull3d_data *, simplex *, void *);
-
-/* starting at s, visit simplices t such that test(s,i,0) is true,
- * and t is the i'th neighbor of s;
- * apply visit function to all visited simplices;
- * when visit returns nonNULL, exit and return its value */
-HIDDEN void *
-chull3d_visit_triang_gen(struct chull3d_data *cdata, simplex *s,
chull3d_visit_func *visit, chull3d_test_func *test)
-{
- neighbor *sn;
- void *v;
- simplex *t;
- int i;
- long tms = 0;
-
-#if 0
- /* Hey Cliff, can you check this? */
-
-/p/home/morrison/brlcad.trunk/src/libbg/chull3d.cpp(594): error #170: pointer
points outside of underlying object
- for (i=-1,sn = t->neigh-1;i<(cdata->cdim);i++,sn++)
- ^
-
-/p/home/morrison/brlcad.trunk/src/libbg/chull3d.cpp(626): error #170: pointer
points outside of underlying object
- for (i=-1,sn=s->neigh-1;i<(cdata->cdim);i++,sn++) {
- ^
-
-/p/home/morrison/brlcad.trunk/src/libbg/chull3d.cpp(634): error #170: pointer
points outside of underlying object
- for (j=-1,snn=sns->neigh-1; j<(cdata->cdim) && snn->simp!=s; j++,snn++)
{};
- ^
-
-/p/home/morrison/brlcad.trunk/src/libbg/chull3d.cpp(638): error #170: pointer
points outside of underlying object
- for (k=-1,snn=sns->neigh-1; k<(cdata->cdim); k++,snn++){
- ^
-
-/p/home/morrison/brlcad.trunk/src/libbg/chull3d.cpp(641): error #170: pointer
points outside of underlying object
- for (l=-1,sn2 = s->neigh-1;
- ^
-
-/p/home/morrison/brlcad.trunk/src/libbg/chull3d.cpp(825): error #170: pointer
points outside of underlying object
- CHULL3D_copy_simp(cdata, cdata->ns, seen);
- ^
-
-/p/home/morrison/brlcad.trunk/src/libbg/chull3d.cpp(860): error #170: pointer
points outside of underlying object
- CHULL3D_copy_simp(cdata, ns, s);
- ^
-
-/p/home/morrison/brlcad.trunk/src/libbg/chull3d.cpp(972): error #170: pointer
points outside of underlying object
- CHULL3D_copy_simp(cdata, s, root);
- ^
-#endif
-
- (cdata->vnum)--;
- if (s) CHULL3D_push(s);
- while (tms) {
- if (tms>cdata->ss) {
-
cdata->st=(simplex**)realloc(cdata->st,((cdata->ss+=cdata->ss)+MAXDIM+1)*sizeof(simplex*));
- }
- CHULL3D_pop(t);
- if (!t || t->visit == cdata->vnum) continue;
- t->visit = cdata->vnum;
- if ((v=(*visit)(cdata,t,0))) {return v;}
- for (i=-1,sn = t->neigh-1; i<(cdata->cdim); i++,sn++)
- if ((sn->simp->visit != cdata->vnum) && sn->simp &&
test(cdata,t,i,0))
- CHULL3D_push(sn->simp);
- }
- return NULL;
-}
-
-/* visit the whole triangulation */
-HIDDEN int chull3d_truet(struct chull3d_data *UNUSED(cdata), simplex
*UNUSED(s), int UNUSED(i), void *UNUSED(dum)) {return 1;}
-HIDDEN void *
-chull3d_visit_triang(struct chull3d_data *cdata, simplex *root,
chull3d_visit_func *visit)
-{
- return chull3d_visit_triang_gen(cdata, root, visit, chull3d_truet);
-}
-
-
-HIDDEN int chull3d_hullt(struct chull3d_data *UNUSED(cdata), simplex
*UNUSED(s), int i, void *UNUSED(dummy)) {return i>-1;}
-HIDDEN void *chull3d_facet_test(struct chull3d_data *UNUSED(cdata), simplex
*s, void *UNUSED(dummy)) {return (!s->peak.vert) ? s : NULL;}
-HIDDEN void *
-chull3d_visit_hull(struct chull3d_data *cdata, simplex *root,
chull3d_visit_func *visit)
-{
- return chull3d_visit_triang_gen(cdata, (struct simplex
*)chull3d_visit_triang(cdata, root, &chull3d_facet_test), visit, chull3d_hullt);
-}
-
-HIDDEN void *
-chull3d_collect_hull_pnts(struct chull3d_data *cdata, simplex *s, void
*UNUSED(p)) {
-
- point v[MAXDIM];
- int j;
- long int ip[3];
- /*point_t p1, p2, p3;*/
- const point_t *pp[3];
- int f;
-
- if (!s) return NULL;
-
- for (j=0;j<(cdata->cdim);j++) v[j] = s->neigh[j].vert;
- for (j=0;j<(cdata->cdim);j++) ip[j] = (cdata->site_num)((void *)cdata,
v[j]);
- for (j=0;j<(cdata->cdim);j++) pp[j] = &(cdata->input_vert_array[ip[j]]);
-
- /* Don't add a point if it's already added */
- for (j=0;j<(cdata->cdim);j++){
- std::map<long *, int>::iterator pli;
- pli = cdata->point_lookup->find((long *)pp[j]);
- if (pli == cdata->point_lookup->end()) {
- f = bu_ptbl_ins(cdata->output_pnts, (long *)pp[j]);
- cdata->point_lookup->insert(std::make_pair((long *)pp[j], f));
- VMOVE(cdata->vert_array[f],*pp[j]);
- (*cdata->vert_cnt)++;
- }
- }
- return NULL;
-}
-
-
-HIDDEN void *
-chull3d_collect_faces(struct chull3d_data *cdata, simplex *s, void *UNUSED(p))
{
-
- point v[MAXDIM];
- int j;
- long int ip[3];
- /*point_t p1, p2, p3;*/
- const point_t *pp[3];
- int f[3];
- vect_t a, b, normal, center_to_edge;
-
- if (!s) return NULL;
-
- /* If we're not 3d, it doesn't make sense to do this */
- if (cdata->cdim != 3) return NULL;
-
- for (j=0;j<(cdata->cdim);j++) v[j] = s->neigh[j].vert;
- for (j=0;j<(cdata->cdim);j++) ip[j] = (cdata->site_num)((void *)cdata,
v[j]);
- for (j=0;j<(cdata->cdim);j++) pp[j] = &(cdata->input_vert_array[ip[j]]);
-
- for (j=0;j<(cdata->cdim);j++){
- std::map<long *, int>::iterator pli;
- pli = cdata->point_lookup->find((long *)pp[j]);
- f[j] = pli->second;
- }
-
- /* Use cdata->center_pnt and the center pnt of the triangle to construct
- * a vector, and find the normal of the proposed face. If the face normal
- * does not point in the same direction (dot product is positive) then
- * the face needs to be reversed. Do this so the eventual bot primitive
- * won't need a bot sync. Since by definition the chull is convex, the
- * center point should be "inside" relative to all faces and this test
- * should be reliable. */
- VSUB2(a, cdata->vert_array[f[1]], cdata->vert_array[f[0]]);
- VSUB2(b, cdata->vert_array[f[2]], cdata->vert_array[f[0]]);
- VCROSS(normal, a, b);
- VUNITIZE(normal);
-
- VSUB2(center_to_edge, cdata->vert_array[f[0]], cdata->center_pnt);
- VUNITIZE(center_to_edge);
-
- cdata->faces[(*cdata->face_cnt)*3] = f[0];
- if(VDOT(normal, center_to_edge) < 0) {
- cdata->faces[(*cdata->face_cnt)*3+1] = f[2];
- cdata->faces[(*cdata->face_cnt)*3+2] = f[1];
- } else {
- cdata->faces[(*cdata->face_cnt)*3+1] = f[1];
- cdata->faces[(*cdata->face_cnt)*3+2] = f[2];
- }
- (*cdata->face_cnt)++;
- return NULL;
-}
-
-/* the neighbor entry of a containing b */
-HIDDEN neighbor *
-chull3d_op_simp(struct chull3d_data *cdata, simplex *a, simplex *b)
-{
- int i;
- neighbor *x;
- for (i=0, x = a->neigh; (x->simp != b) && (i<(cdata->cdim)) ; i++, x++) {};
- if (i<(cdata->cdim))
- return x;
- else {
- return 0;
- }
-}
-
-HIDDEN neighbor *
-chull3d_op_vert(struct chull3d_data *cdata, simplex *a, site b)
-{
- int i;
- neighbor *x;
- for (i=0, x = a->neigh; (x->vert != b) && (i<(cdata->cdim)) ; i++, x++) {};
- if (i<(cdata->cdim))
- return x;
- else {
- return 0;
- }
-}
-
-
-/* make neighbor connections between newly created simplices incident to p */
-HIDDEN void
-chull3d_connect(struct chull3d_data *cdata, simplex *s)
-{
- site xf,xb,xfi;
- simplex *sb, *sf, *seen;
- int i;
- neighbor *sn;
-
- if (!s) return;
- //assert(!s->peak.vert && s->peak.simp->peak.vert==cdata->p &&
!chull3d_op_vert(cdata,s,cdata->p)->simp->peak.vert);
- if (s->visit==cdata->pnum) return;
- s->visit = cdata->pnum;
- seen = s->peak.simp;
- xfi = chull3d_op_simp(cdata,seen,s)->vert;
- for (i=0, sn = s->neigh; i<(cdata->cdim); i++,sn++) {
- xb = sn->vert;
- if (cdata->p == xb) continue;
- sb = seen;
- sf = sn->simp;
- xf = xfi;
- if (!sf->peak.vert) { /* are we done already? */
- sf = chull3d_op_vert(cdata,seen,xb)->simp;
- if (sf->peak.vert) continue;
- } else do {
- xb = xf;
- xf = chull3d_op_simp(cdata,sf,sb)->vert;
- sb = sf;
- sf = chull3d_op_vert(cdata,sb,xb)->simp;
- } while (sf->peak.vert);
-
- sn->simp = sf;
- chull3d_op_vert(cdata,sf,xf)->simp = s;
-
- chull3d_connect(cdata, sf);
- }
-}
-
-
-/* visit simplices s with sees(p,s), and make a facet for every neighbor
- * of s not seen by p */
-HIDDEN simplex *
-chull3d_make_facets(struct chull3d_data *cdata, simplex *seen)
-{
- simplex *n;
- neighbor *bn;
- int i;
-
-
- if (!seen) return NULL;
- seen->peak.vert = cdata->p;
-
- for (i=0,bn = seen->neigh; i<(cdata->cdim); i++,bn++) {
- neighbor *ns;
- n = bn->simp;
- if (cdata->pnum != n->visit) {
- n->visit = cdata->pnum;
- if (chull3d_sees(cdata, cdata->p,n)) chull3d_make_facets(cdata, n);
- }
- if (n->peak.vert) continue;
- CHULL3D_copy_simp(cdata, cdata->ns, seen);
- cdata->ns->visit = 0;
- cdata->ns->peak.vert = 0;
- cdata->ns->normal = 0;
- cdata->ns->peak.simp = seen;
- CHULL3D_NULLIFY(cdata, basis_s,cdata->ns->neigh[i].basis);
- cdata->ns->neigh[i].vert = cdata->p;
- bn->simp = cdata->ns;
- ns = chull3d_op_simp(cdata,n,seen);
- ns->simp = cdata->ns;
- }
- return cdata->ns;
-}
-
-
-/* p lies outside flat containing previous sites;
- * make p a vertex of every current simplex, and create some new simplices */
-HIDDEN simplex *
-chull3d_extend_simplices(struct chull3d_data *cdata, simplex *s)
-{
- int i;
- int ocdim=(cdata->cdim)-1;
- simplex *ns = NULL;
- neighbor *nsn = NULL;
-
-
- if (s->visit == cdata->pnum) return s->peak.vert ? s->neigh[ocdim].simp :
s;
- s->visit = cdata->pnum;
- s->neigh[ocdim].vert = cdata->p;
- CHULL3D_NULLIFY(cdata, basis_s,s->normal);
- CHULL3D_NULLIFY(cdata, basis_s,s->neigh[0].basis);
- if (!s->peak.vert) {
- s->neigh[ocdim].simp = chull3d_extend_simplices(cdata, s->peak.simp);
- return s;
- } else {
- CHULL3D_copy_simp(cdata, ns, s);
- s->neigh[ocdim].simp = ns;
- ns->peak.vert = NULL;
- ns->peak.simp = s;
- ns->neigh[ocdim] = s->peak;
- if (s->peak.basis) s->peak.basis->ref_count++;
- for (i=0,nsn=ns->neigh;i<(cdata->cdim);i++,nsn++)
- nsn->simp = chull3d_extend_simplices(cdata, nsn->simp);
- }
- return ns;
-}
-
-
-/* return a simplex s that corresponds to a facet of the
- * current hull, and sees(p, s) */
-HIDDEN simplex *
-chull3d_search(struct chull3d_data *cdata, simplex *root)
-{
- simplex *s;
- neighbor *sn;
- int i;
- long tms = 0;
-
-
- CHULL3D_push(root->peak.simp);
- root->visit = cdata->pnum;
- if (!chull3d_sees(cdata, cdata->p,root))
- for (i=0,sn=root->neigh;i<(cdata->cdim);i++,sn++)
CHULL3D_push(sn->simp);
- while (tms) {
- if (tms>cdata->ss)
- cdata->st=(simplex**)realloc(cdata->st,
((cdata->ss+=cdata->ss)+MAXDIM+1)*sizeof(simplex*));
- CHULL3D_pop(s);
- if (s->visit == cdata->pnum) continue;
- s->visit = cdata->pnum;
- if (!chull3d_sees(cdata, cdata->p,s)) continue;
- if (!s->peak.vert) return s;
- for (i=0, sn=s->neigh; i<(cdata->cdim); i++,sn++)
CHULL3D_push(sn->simp);
- }
- return NULL;
-}
-
-
-HIDDEN point
-chull3d_get_another_site(struct chull3d_data *cdata)
-{
- point pnext;
- pnext = (*(cdata->get_site))((void *)cdata);
- if (!pnext) return NULL;
- cdata->pnum = cdata->site_num((void *)cdata, pnext)+2;
- return pnext;
-}
-
-
-
-HIDDEN void
-chull3d_buildhull(struct chull3d_data *cdata, simplex *root)
-{
- while ((cdata->cdim) < cdata->rdim) {
- cdata->p = chull3d_get_another_site(cdata);
- if (!cdata->p) return;
- if (chull3d_out_of_flat(cdata, root,cdata->p))
- chull3d_extend_simplices(cdata, root);
- else
- chull3d_connect(cdata, chull3d_make_facets(cdata,
chull3d_search(cdata, root)));
- }
- while ((cdata->p = chull3d_get_another_site(cdata)))
- chull3d_connect(cdata, chull3d_make_facets(cdata, chull3d_search(cdata,
root)));
-}
-
-
-/*
- get_s returns next site each call;
- hull construction stops when NULL returned;
- site_numm returns number of site when given site;
- dim dimension of point set;
- */
-HIDDEN simplex *
-chull3d_build_convex_hull(struct chull3d_data *cdata, gsitef *get_s, site_n
*site_numm, short dim)
-{
- double febits;
- simplex *s = NULL;
- simplex *root = NULL;
-
- if (ZERO(cdata->Huge)) cdata->Huge = DBL_MAX;
-
- (cdata->cdim) = 0;
- cdata->get_site = get_s;
- cdata->site_num = site_numm;
- cdata->pdim = dim;
-
- febits = floor(DBL_MANT_DIG*log(double(FLT_RADIX))/log(2.));
- cdata->exact_bits = (int)febits;
- cdata->b_err_min = (float)(DBL_EPSILON*MAXDIM*(1<<MAXDIM)*MAXDIM*3.01);
- cdata->b_err_min_sq = cdata->b_err_min * cdata->b_err_min;
-
- //assert(cdata->get_site!=NULL); assert(cdata->site_num!=NULL);
-
- cdata->rdim = cdata->pdim;
-
- cdata->point_size = cdata->site_size = sizeof(Coord)*cdata->pdim;
- cdata->basis_vec_size = sizeof(Coord)*cdata->rdim;
- cdata->basis_s_size = sizeof(basis_s)+ (2*cdata->rdim-1)*sizeof(Coord);
- cdata->simplex_size = sizeof(simplex) + (cdata->rdim-1)*sizeof(neighbor);
- cdata->Tree_size = sizeof(Tree);
-
- root = NULL;
- if (!(cdata->p = (*(cdata->get_site))((void *)cdata))) return 0;
-
- CHULL3D_NEWL(cdata, cdata->simplex_list, simplex,root);
-
- cdata->ch_root = root;
-
- CHULL3D_copy_simp(cdata, s, root);
- root->peak.vert = cdata->p;
- root->peak.simp = s;
- s->peak.simp = root;
-
- chull3d_buildhull(cdata, root);
- return root;
-}
-
-
-/* hullmain.c */
-
-HIDDEN long
-chull3d_site_numm(void *data, site p)
-{
- long i,j;
- struct chull3d_data *cdata = (struct chull3d_data *)data;
-
- if (!p) return -2;
- for (i=0; i<cdata->num_blocks; i++)
- if ((j=p-cdata->site_blocks[i])>=0 && j < BLOCKSIZE*cdata->pdim)
- return j/cdata->pdim+BLOCKSIZE*i;
- return -3;
-}
-
-HIDDEN site
-chull3d_new_site(struct chull3d_data *cdata, site p, long j)
-{
- //assert(cdata->num_blocks+1<cdata->input_vert_cnt);
- if (0==(j%BLOCKSIZE)) {
- //assert(cdata->num_blocks < cdata->input_vert_cnt);
-
return(cdata->site_blocks[cdata->num_blocks++]=(site)malloc(BLOCKSIZE*cdata->site_size));
- } else
- return p+cdata->pdim;
-}
-
-HIDDEN site
-chull3d_read_next_site(struct chull3d_data *cdata, long j)
-{
- int i=0;
-
- cdata->p = chull3d_new_site(cdata, cdata->p,j);
-
- if (cdata->next_vert == cdata->input_vert_cnt) return 0;
-
- cdata->p[0] = cdata->input_vert_array[cdata->next_vert][0];
- cdata->p[1] = cdata->input_vert_array[cdata->next_vert][1];
- cdata->p[2] = cdata->input_vert_array[cdata->next_vert][2];
- (cdata->next_vert)++;
- for(i = 0; i < 3; i++) {
- (cdata->p)[i] = floor(cdata->mult_up*(cdata->p)[i]+0.5);
- cdata->mins[i] = (cdata->mins[i]<(cdata->p)[i]) ? cdata->mins[i] :
(cdata->p)[i];
- cdata->maxs[i] = (cdata->maxs[i]>(cdata->p)[i]) ? cdata->maxs[i] :
(cdata->p)[i];
- }
-
- return cdata->p;
-}
-
-HIDDEN void
-chull3d_make_shuffle(struct chull3d_data *cdata)
-{
- long i, j, t;
- for (i=0;i<cdata->input_vert_cnt;i++) cdata->shufmat[i] = i;
- for (i=0;i<cdata->input_vert_cnt;i++) {
- double random_num = bn_randmt();
- t = cdata->shufmat[i];
- j = i + (cdata->input_vert_cnt - i) * (long)random_num;
- cdata->shufmat[i] = cdata->shufmat[j];
- cdata->shufmat[j] = t;
- }
-}
-
-HIDDEN site
-chull3d_get_next_site(void *data) {
- struct chull3d_data *cdata = (struct chull3d_data *)data;
- return chull3d_read_next_site(cdata, cdata->shufmat[(cdata->s_num)++]);
-}
-
-HIDDEN void
-chull3d_data_init(struct chull3d_data *data, int vert_cnt)
-{
- int i = 0;
- data->simplex_list = 0;
- data->basis_s_list = 0;
- data->Tree_list = 0;
- data->site_blocks = (point *)bu_calloc(vert_cnt * 3, sizeof(point),
"site_blocks");
- BU_GET(data->tt_basis, basis_s);
- data->tt_basis->next = 0;
- data->tt_basis->ref_count = 1;
- data->tt_basis->lscale = -1;
- data->tt_basis->sqa = 0;
- data->tt_basis->sqb = 0;
- data->tt_basis->vecs[0] = 0;
- /* point at infinity for Delaunay triang; value not used */
- data->hull_infinity = (Coord *)bu_calloc(10, sizeof(Coord),
"hull_infinity");
- data->hull_infinity[0] = 57.2;
- data->B = (int *)bu_calloc(vert_cnt, sizeof(int), "A");
- data->mins = (Coord *)bu_calloc(MAXDIM, sizeof(Coord), "mins");
- for(i = 0; i < MAXDIM; i++) data->mins[i] = DBL_MAX;
- data->maxs = (Coord *)bu_calloc(MAXDIM, sizeof(Coord), "maxs");
- for(i = 0; i < MAXDIM; i++) data->maxs[i] = -DBL_MAX;
- data->mult_up = 1000.0; /* we'll need to multiply based on a tolerance
parameter at some point... */
- data->mi = (short *)bu_calloc(vert_cnt, sizeof(short), "mi");
- data->mo = (short *)bu_calloc(vert_cnt, sizeof(short), "mo");
- data->pdim = 3;
- data->shufmat = (long*)bu_calloc(vert_cnt+1, sizeof(long), "shufmat");
-
- /* These were static variables in functions */
- data->simplex_block_table = (simplex **)bu_calloc(vert_cnt, sizeof(simplex
*), "simplex_block_table");
- data->num_simplex_blocks = 0;
- data->basis_s_block_table = (basis_s **)bu_calloc(vert_cnt, sizeof(basis_s
*), "basis_s_block_table");
- data->num_basis_s_blocks = 0;
- data->Tree_block_table = (Tree **)bu_calloc(vert_cnt, sizeof(Tree *),
"Tree_block_table");
- data->num_Tree_blocks = 0;
- data->p_neigh.vert = 0;
- data->p_neigh.simp = NULL;
- data->p_neigh.basis = NULL;
- data->b = NULL;
- data->vnum = -1;
- data->ss = vert_cnt + MAXDIM;
- data->s_num = 0;
- data->fg_vn = 0;
- data->st=(simplex**)malloc((data->ss+MAXDIM+1)*sizeof(simplex*));
- data->ns = NULL;
- data->free_point_lookup = 0;
- data->free_output_pnts = 0;
- if (!data->point_lookup) {
- data->point_lookup = new std::map<long *, int>;
- data->free_point_lookup = 1;
- }
- if (!data->output_pnts) {
- BU_GET(data->output_pnts, struct bu_ptbl);
- bu_ptbl_init(data->output_pnts, vert_cnt, "output pnts container");
- data->free_output_pnts = 1;
- }
- data->input_vert_cnt = vert_cnt;
-
- VSET(data->center_pnt,0,0,0);
-
- /* Output containers */
- if (!data->vert_array) data->vert_array = (point_t *)bu_calloc(vert_cnt,
sizeof(point_t), "vertex array");
- data->next_vert = 0;
- if (!data->faces) data->faces = (int *)bu_calloc(3*vert_cnt, sizeof(int),
"face array");
-}
-
-HIDDEN void
-chull3d_data_free(struct chull3d_data *data)
-{
- bu_free(data->site_blocks, "site_blocks");
- bu_free(data->hull_infinity, "hull_infinity");
- bu_free(data->B, "B");
- bu_free(data->mins, "mins");
- bu_free(data->maxs, "maxs");
- bu_free(data->mi, "mi");
- bu_free(data->mo, "mo");
- bu_free(data->simplex_block_table, "simplex_block_table");
- bu_free(data->basis_s_block_table, "basis_s_block_table");
- bu_free(data->Tree_block_table, "Tree_block_table");
- BU_PUT(data->tt_basis, basis_s);
- if (data->free_output_pnts) {
- bu_ptbl_free(data->output_pnts);
- BU_PUT(data->output_pnts, struct bu_ptbl);
- }
- if (data->free_point_lookup) delete data->point_lookup;
-}
-
-HIDDEN
-int
-chull3d_intermediate_set(point_t **vertices, int *num_vertices, const point_t
*input_points_3d, int num_input_pnts)
-{
- int nv = 0;
- int nf = 0;
- int curr_stp = 0;
- int pnt_stp = 0;
- int last_stp = 0;
- int pnts_per_stp = CHULL3D_SUBSET_SIZE;
- const point_t *curr_inputs;
- struct bu_ptbl *opnts;
- std::map<long *, int> *pl = new std::map<long *, int>;
- point_t *intermediate_vert_array = (point_t *)bu_calloc(num_input_pnts,
sizeof(point_t), "vertex array");
- simplex *root = NULL;
- BU_GET(opnts, struct bu_ptbl);
- bu_ptbl_init(opnts, num_input_pnts, "output pnts container");
- pnt_stp = (num_input_pnts < pnts_per_stp) ? num_input_pnts : pnts_per_stp;
- last_stp = num_input_pnts / pnts_per_stp;
- //bu_log("input point cnt: %d\n", num_input_pnts);
- while (curr_stp <= last_stp) {
- struct chull3d_data *cdata;
- //bu_log("step %d of %d\n", curr_stp, last_stp);
- curr_inputs = input_points_3d + curr_stp * pnt_stp;
- if (curr_stp == last_stp) pnt_stp = num_input_pnts - curr_stp * pnt_stp;
- BU_GET(cdata, struct chull3d_data);
- cdata->output_pnts = opnts;
- cdata->point_lookup = pl;
- cdata->vert_array = intermediate_vert_array;
- chull3d_data_init(cdata, pnt_stp);
- cdata->input_vert_array = curr_inputs;
- cdata->vert_cnt = &nv;
- cdata->face_cnt = &nf;
- chull3d_make_shuffle(cdata);
- root = chull3d_build_convex_hull(cdata, chull3d_get_next_site,
chull3d_site_numm, cdata->pdim);
- if (!root) return -1;
- chull3d_visit_hull(cdata, root, chull3d_collect_hull_pnts);
- chull3d_free_hull_storage(cdata);
- chull3d_data_free(cdata);
- BU_PUT(cdata, struct chull3d_data);
- curr_stp++;
- //bu_log("points accumulated: %d\n", nv);
- }
- //bu_log("output pnt cnt: %d\n", nv);
- bu_ptbl_free(opnts);
- BU_PUT(opnts, struct bu_ptbl);
- delete pl;
-
- (*vertices) = intermediate_vert_array;
- (*num_vertices) = nv;
- return 0;
-}
-
-int
-bg_3d_chull(int **faces, int *num_faces, point_t **vertices, int *num_vertices,
- const point_t *input_points_3d, int num_input_pnts)
-{
- int i;
- struct chull3d_data *cdata;
- simplex *root = NULL;
- int output_dim = 0;
- point_t* hull_2d;
-
- const point_t *iv_1;
- point_t *iv_2, *iva;
- int nv_1, nv_2, nva;
- iv_1 = input_points_3d;
- nv_1 = num_input_pnts;
-
- iv_2 = (point_t *) 0;
- nv_2 = 0;
-
- (void)chull3d_intermediate_set(&iv_2, &nv_2, iv_1, nv_1);
- iva = iv_2;
- nva = nv_2;
- while (nva > CHULL3D_COUNT_THRESHOLD && abs(nv_1 - nv_2) >
CHULL3D_DELTA_THRESHOLD) {
- nv_1 = nv_2;
- iv_1 = (const point_t *)iv_2;
- (void)chull3d_intermediate_set(&iv_2, &nv_2, iv_1, nv_1);
- bu_free((point_t *)iv_1, "free old set");
- iva = iv_2;
- nva = nv_2;
- }
-
- BU_GET(cdata, struct chull3d_data);
- chull3d_data_init(cdata, nva);
- cdata->input_vert_array = (const point_t *)iva;
- cdata->vert_cnt = num_vertices;
- cdata->face_cnt = num_faces;
- (*cdata->vert_cnt) = 0;
- (*cdata->face_cnt) = 0;
- chull3d_make_shuffle(cdata);
- root = chull3d_build_convex_hull(cdata, chull3d_get_next_site,
chull3d_site_numm, cdata->pdim);
@@ Diff output truncated at 100000 characters. @@
This was sent by the SourceForge.net collaborative development platform, the
world's largest Open Source development site.
_______________________________________________
BRL-CAD Source Commits mailing list
brlcad-commits@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/brlcad-commits
