From 11716ba56b8f63264d6205ca9aa674671d40593f Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Fri, 11 Mar 2022 11:45:15 -0800 Subject: [PATCH 01/39] polymorphic device functions working --- CMakeLists.txt | 1 + sdf/CMakeLists.txt | 10 +++++++ sdf/include/sdf.cuh | 69 ++++++++++++++++++++++++++++++++++++++++++++ tools/CMakeLists.txt | 8 ++++- tools/sdfTest.cu | 16 ++++++++++ 5 files changed, 103 insertions(+), 1 deletion(-) create mode 100644 sdf/CMakeLists.txt create mode 100644 sdf/include/sdf.cuh create mode 100644 tools/sdfTest.cu diff --git a/CMakeLists.txt b/CMakeLists.txt index 6648807ef..139eae1f9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -39,6 +39,7 @@ IF(MANIFOLD_USE_CPP) ENDIF(MANIFOLD_USE_CPP) add_subdirectory(utilities) +add_subdirectory(sdf) add_subdirectory(collider) add_subdirectory(polygon) add_subdirectory(manifold) diff --git a/sdf/CMakeLists.txt b/sdf/CMakeLists.txt new file mode 100644 index 000000000..29cb6dcc0 --- /dev/null +++ b/sdf/CMakeLists.txt @@ -0,0 +1,10 @@ +project (sdf) + +add_library(${PROJECT_NAME} INTERFACE) + +target_include_directories(${PROJECT_NAME} + INTERFACE + ${PROJECT_SOURCE_DIR}/include +) + +target_link_libraries( ${PROJECT_NAME} INTERFACE utilities) diff --git a/sdf/include/sdf.cuh b/sdf/include/sdf.cuh new file mode 100644 index 000000000..d6c9d89a9 --- /dev/null +++ b/sdf/include/sdf.cuh @@ -0,0 +1,69 @@ +// Copyright 2022 Emmett Lalish +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#pragma once + +#include "structs.h" +#include "utils.cuh" +#include "vec_dh.cuh" + +namespace { +using namespace manifold; + +template +struct Compute { + const Func func; + + __host__ __device__ void operator()( + thrust::tuple inOut) const { + float& value = thrust::get<0>(inOut); + const glm::vec3 vert = thrust::get<1>(inOut); + + value = func(vert); + } +}; +} // namespace + +namespace manifold { +/** @addtogroup Core + * @{ + */ +template +class SDF { + public: + SDF(Func sdf) : sdf_{sdf} {}; + + inline __host__ __device__ float operator()(glm::vec3 point) const { + return sdf_(point); + } + + inline Mesh LevelSet(Box bounds, float edgeLength, float level = 0) const { + Mesh out; + int size = 10; + VecDH verts(size); + VecDH values(size); + thrust::fill(verts.beginD(), verts.endD(), glm::vec3(1.0f, 2.0f, 3.0f)); + thrust::for_each_n(zip(values.beginD(), verts.beginD()), size, + Compute({sdf_})); + verts.Dump(); + values.Dump(); + return out; + } + + private: + const Func sdf_; +}; + +/** @} */ +} // namespace manifold \ No newline at end of file diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index 4b502c222..156f85e6a 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -22,7 +22,13 @@ target_compile_options(perfTestCGAL PRIVATE ${MANIFOLD_FLAGS}) target_compile_features(perfTestCGAL PUBLIC cxx_std_14) # add_executable(playground playground.cpp) -# target_link_libraries(playground manifold meshIO samples) +# target_link_libraries(playground manifold meshIO samples sdf) # target_compile_options(playground PRIVATE ${MANIFOLD_FLAGS}) # target_compile_features(playground PUBLIC cxx_std_14) + +add_executable(sdfTest sdfTest.cu) +target_link_libraries(sdfTest manifold meshIO samples sdf) + +target_compile_options(sdfTest PRIVATE ${MANIFOLD_FLAGS}) +target_compile_features(sdfTest PUBLIC cxx_std_14) diff --git a/tools/sdfTest.cu b/tools/sdfTest.cu new file mode 100644 index 000000000..a7940cb1a --- /dev/null +++ b/tools/sdfTest.cu @@ -0,0 +1,16 @@ +#include "sdf.cuh" + +using namespace manifold; + +struct Test { + __host__ __device__ float operator()(glm::vec3 point) const { + return point[1]; + } +}; + +int main(int argc, char **argv) { + Test func; + SDF a(func); + Box box; + Mesh b = a.LevelSet(box, 1); +} \ No newline at end of file From e34b939a4771320edf3ac82a05790e50be037b02 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Fri, 11 Mar 2022 17:10:25 -0800 Subject: [PATCH 02/39] starting to flesh out algorithm --- sdf/include/sdf.cuh | 107 +++++++++++++++++++++++++++++++++++++------- 1 file changed, 91 insertions(+), 16 deletions(-) diff --git a/sdf/include/sdf.cuh b/sdf/include/sdf.cuh index d6c9d89a9..f0ac24c9d 100644 --- a/sdf/include/sdf.cuh +++ b/sdf/include/sdf.cuh @@ -21,18 +21,80 @@ namespace { using namespace manifold; +constexpr kOpen = std::numeric_limits::max(); + +struct GridVert { + uint64_t key = kOpen; + float distance = 0.0f / 0.0f; + int edgeIndex = -1; + glm::vec3[7] edgeIntersections; +}; + +class HashTable { + public: + HashTable(uint32_t sizeExp = 20, uint32_t step = 127) + : size_{1 << sizeExp}, step_{step}, table_{1 << sizeExp} {} + + int Entries() const { return used_.H()[0]; } + + float FilledFraction() const { + return static_cast(used_.H()[0]) / size_; + } + + __device__ __host__ bool Insert(const GridVert& vert) { + uint32_t idx = vert.key & (size_ - 1); + while (1) { + const uint64_t found = + AtomicCAS(&table_.ptrD()[idx].key, kOpen, vert.key); + if (found == kOpen) { + if (AtomicAdd(&used_.ptrD()[0], 1) * 2 > size_) { + return true; + } + table_.ptrD()[idx] = vert; + return false; + } + if (found == vert.key) return false; + idx = (idx + step_) & (size_ - 1); + } + } + + __device__ __host__ GridVert operator[](uint64_t key) const { + uint32_t idx = key & (size_ - 1); + while (1) { + const GridVert found = table_.ptrD()[idx]; + if (found.key == key) return found; + if (found.key == kOpen) return GridVert(); + idx = (idx + step_) & (size_ - 1); + } + } + + private: + const uint32_t size_; + const uint32_t step_; + VecDH table_; + VecDH used_(1, 0); +}; + template -struct Compute { +struct ComputeVerts { + HashTable gridVerts; const Func func; + const glm::vec3 origin; + const float spacing; - __host__ __device__ void operator()( - thrust::tuple inOut) const { - float& value = thrust::get<0>(inOut); - const glm::vec3 vert = thrust::get<1>(inOut); - - value = func(vert); + __host__ __device__ void operator()(uint64_t mortonCode) const { + GridVert gridVert; + const glm::ivec3 gridIndex = DecodeMorton(mortonCode); + const glm::vec3 position = origin + spacing * gridIndex; + value = func(position); } }; + +struct BuildTris { + uint64_t* halfedgeCodes; + uint64_t* index; + const HashTable gridVerts; +}; } // namespace namespace manifold { @@ -42,7 +104,7 @@ namespace manifold { template class SDF { public: - SDF(Func sdf) : sdf_{sdf} {}; + SDF(Func sdf) : sdf_{sdf} {} inline __host__ __device__ float operator()(glm::vec3 point) const { return sdf_(point); @@ -50,14 +112,27 @@ class SDF { inline Mesh LevelSet(Box bounds, float edgeLength, float level = 0) const { Mesh out; - int size = 10; - VecDH verts(size); - VecDH values(size); - thrust::fill(verts.beginD(), verts.endD(), glm::vec3(1.0f, 2.0f, 3.0f)); - thrust::for_each_n(zip(values.beginD(), verts.beginD()), size, - Compute({sdf_})); - verts.Dump(); - values.Dump(); + HashTable gridVerts(10); + glm::vec3 dim = bounds.Size(); + // Need to create a new MortonCode function that spreads when needed instead + // of always by 3, to make non-cubic domains efficient. + float maxDim = std::max(dim[0], std::max(dim[1], dim[2])); + int n = maxDim / edgeLength; + float spacing = maxDim / n; + int maxMorton = MortonCode(glm::ivec3(n)); + + thrust::for_each_n( + countAt(0), maxMorton, + ComputeVerts({gridVerts, sdf_, bounds.Min(), spacing})); + + VecDH halfedgeCodes(gridVerts.Entries() * 6); + VecDH index(1, 0); + + thrust::for_each_n( + countAt(0), maxMorton, + BuildTris({halfedgeCodes.ptrD(), index.ptrD(), gridVerts})); + + // turn halfedgeCodes into out.triVerts and gridVerts into out.vertPos return out; } From c869c50950968f9bd359a35d8f4a6344c7049e05 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Fri, 11 Mar 2022 19:03:28 -0800 Subject: [PATCH 03/39] getting closer --- sdf/include/sdf.cuh | 29 ++++++++++++++++++++--------- 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/sdf/include/sdf.cuh b/sdf/include/sdf.cuh index f0ac24c9d..e1336fa0b 100644 --- a/sdf/include/sdf.cuh +++ b/sdf/include/sdf.cuh @@ -26,8 +26,9 @@ constexpr kOpen = std::numeric_limits::max(); struct GridVert { uint64_t key = kOpen; float distance = 0.0f / 0.0f; + int vertKey = -1; int edgeIndex = -1; - glm::vec3[7] edgeIntersections; + int[7] edgeVerts; }; class HashTable { @@ -77,6 +78,8 @@ class HashTable { template struct ComputeVerts { + glm::vec3* vertPos; + int* vertIndex; HashTable gridVerts; const Func func; const glm::vec3 origin; @@ -91,9 +94,11 @@ struct ComputeVerts { }; struct BuildTris { - uint64_t* halfedgeCodes; + glm::ivec3* triVerts; uint64_t* index; const HashTable gridVerts; + + __host__ __device__ void operator()(uint64_t mortonCode) const {} }; } // namespace @@ -121,18 +126,24 @@ class SDF { float spacing = maxDim / n; int maxMorton = MortonCode(glm::ivec3(n)); - thrust::for_each_n( - countAt(0), maxMorton, - ComputeVerts({gridVerts, sdf_, bounds.Min(), spacing})); - - VecDH halfedgeCodes(gridVerts.Entries() * 6); + VecDH vertPos(gridVerts.Size() * 7); VecDH index(1, 0); thrust::for_each_n( countAt(0), maxMorton, - BuildTris({halfedgeCodes.ptrD(), index.ptrD(), gridVerts})); + ComputeVerts({vertPos.ptrD(), index.ptrD(), gridVerts, sdf_, + bounds.Min(), spacing})); + + vertPos.Resize(index.H()[0]); + VecDH triVerts(gridVerts.Entries() * 6); + index.H()[0] = 0; + + thrust::for_each_n(countAt(0), maxMorton, + BuildTris({triVerts.ptrD(), index.ptrD(), gridVerts})); + triVerts.Resize(index.H()[0]); - // turn halfedgeCodes into out.triVerts and gridVerts into out.vertPos + out.vertPos = vertPos.H(); + out.triVerts = triVerts.H(); return out; } From 677f7b87fbbc52a9362c54e6754509261dc2c45a Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Mon, 14 Mar 2022 23:32:16 -0700 Subject: [PATCH 04/39] fleshed out vert computation --- sdf/include/sdf.cuh | 35 ++++++++++++++++++++++++++++++----- 1 file changed, 30 insertions(+), 5 deletions(-) diff --git a/sdf/include/sdf.cuh b/sdf/include/sdf.cuh index e1336fa0b..febb4d1d8 100644 --- a/sdf/include/sdf.cuh +++ b/sdf/include/sdf.cuh @@ -26,7 +26,6 @@ constexpr kOpen = std::numeric_limits::max(); struct GridVert { uint64_t key = kOpen; float distance = 0.0f / 0.0f; - int vertKey = -1; int edgeIndex = -1; int[7] edgeVerts; }; @@ -81,15 +80,41 @@ struct ComputeVerts { glm::vec3* vertPos; int* vertIndex; HashTable gridVerts; - const Func func; + const Func sdf; const glm::vec3 origin; const float spacing; + const glm::vec3[7] edgeVec = {{0.5f, 0.5f, 0.5f}, {1.0f, 0.0f, 0.0f}, + {0.0f, 1.0f, 0.0f}, {0.0f, 0.0f, 1.0f}, + {-0.5f, 0.5f, 0.5f}, {0.5f, -0.5f, 0.5f}, + {0.5f, 0.5f, -0.5f}}; - __host__ __device__ void operator()(uint64_t mortonCode) const { + inline __host__ __device__ void operator()(uint64_t mortonCode) { GridVert gridVert; + gridVert.key = mortonCode; const glm::ivec3 gridIndex = DecodeMorton(mortonCode); const glm::vec3 position = origin + spacing * gridIndex; - value = func(position); + gridVert.distance = sdf(position); + float minDist2 = 0.25 * 0.25; + for (int i = 0; i < 14; ++i) { + const int j = i < 7 ? i : i - 7; + const float val = sdf(position + (i < 7 ? 1 : -1) * spacing * edgeVec[j]); + if (val * gridVert.distance > 0 || (val == 0 && gridVert.distance == 0)) + continue; + + const glm::vec3 delta = + edgeVec[j] * (1 - val / (val - gridVert.distance)); + const float dist2 = glm::dot(delta, delta); + if (dist2 < minDist2) { + gridVert.edgeIndex = i; + minDist2 = dist2; + } + if (i < 7) { + const int idx = AtomicAdd(vertIndex, 1); + vertPos[idx] = position + spacing * delta; + gridVert.edgeVerts[i] = idx; + } + } + gridVerts.Insert(gridVert); } }; @@ -98,7 +123,7 @@ struct BuildTris { uint64_t* index; const HashTable gridVerts; - __host__ __device__ void operator()(uint64_t mortonCode) const {} + __host__ __device__ void operator()(uint64_t mortonCode) {} }; } // namespace From 055035fc7925204e76143b8c08af993d2deac380 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Sat, 30 Apr 2022 13:05:38 -0700 Subject: [PATCH 05/39] more progress --- sdf/include/sdf.cuh | 102 ++++++++++++++++++++++++++++++++------------ 1 file changed, 74 insertions(+), 28 deletions(-) diff --git a/sdf/include/sdf.cuh b/sdf/include/sdf.cuh index febb4d1d8..c9ad00f41 100644 --- a/sdf/include/sdf.cuh +++ b/sdf/include/sdf.cuh @@ -25,21 +25,17 @@ constexpr kOpen = std::numeric_limits::max(); struct GridVert { uint64_t key = kOpen; - float distance = 0.0f / 0.0f; + float distance = NAN; int edgeIndex = -1; int[7] edgeVerts; }; -class HashTable { +class HashTableD { public: - HashTable(uint32_t sizeExp = 20, uint32_t step = 127) - : size_{1 << sizeExp}, step_{step}, table_{1 << sizeExp} {} + HashTableD(uint32_t step = 127, VecDH& alloc, VecDH& used) + : step_{step}, table_{alloc}, used_{used} {} - int Entries() const { return used_.H()[0]; } - - float FilledFraction() const { - return static_cast(used_.H()[0]) / size_; - } + __device__ __host__ int Size() const { return table_.size(); } __device__ __host__ bool Insert(const GridVert& vert) { uint32_t idx = vert.key & (size_ - 1); @@ -69,10 +65,30 @@ class HashTable { } private: - const uint32_t size_; const uint32_t step_; - VecDH table_; + VecD table_; + VecD used_(1, 0); +}; + +class HashTable { + public: + HashTable(uint32_t sizeExp = 20, uint32_t step = 127) + : alloc_{1 << sizeExp}, table_{step, alloc_, used_} {} + + HashTableD D() { return table_; } + + int Entries() const { return used_.H()[0]; } + + int Size() const { return table_.size(); } + + float FilledFraction() const { + return static_cast(used_.H()[0]) / table_.size(); + } + + private: + VecDH alloc_; VecDH used_(1, 0); + HashTableD table_; }; template @@ -81,26 +97,49 @@ struct ComputeVerts { int* vertIndex; HashTable gridVerts; const Func sdf; + const glm::ivec3 gridSize; const glm::vec3 origin; - const float spacing; + const glm::vec3 spacing; const glm::vec3[7] edgeVec = {{0.5f, 0.5f, 0.5f}, {1.0f, 0.0f, 0.0f}, {0.0f, 1.0f, 0.0f}, {0.0f, 0.0f, 1.0f}, {-0.5f, 0.5f, 0.5f}, {0.5f, -0.5f, 0.5f}, {0.5f, 0.5f, -0.5f}}; + inline __host__ __device__ bool AtBounds(glm::ivec3 gridIndex) const { + return gridIndex.x == 0 || gridIndex.x == gridSize.x || gridIndex.y == 0 || + gridIndex.y == gridSize.y || gridIndex.z == 0 || + gridIndex.z == gridSize.z; + } + + inline __host__ __device__ float Sdf(glm::vec3 base, int i) const { + sdf(base + (i < 7 ? 1 : -1) * spacing * edgeVec[j]) + } + + // inline __host__ __device__ float BoundedSdf(glm::vec3 base, int i) const {} + inline __host__ __device__ void operator()(uint64_t mortonCode) { GridVert gridVert; gridVert.key = mortonCode; const glm::ivec3 gridIndex = DecodeMorton(mortonCode); + + // const auto sdfFunc = + // AtBounds(gridIndex) ? &ComputeVerts::BoundedSdf : &ComputeVerts::Sdf; + const glm::vec3 position = origin + spacing * gridIndex; gridVert.distance = sdf(position); + + bool keep = false; float minDist2 = 0.25 * 0.25; for (int i = 0; i < 14; ++i) { const int j = i < 7 ? i : i - 7; - const float val = sdf(position + (i < 7 ? 1 : -1) * spacing * edgeVec[j]); + const float val = Sdf(position, i); if (val * gridVert.distance > 0 || (val == 0 && gridVert.distance == 0)) continue; + keep = true; + // Record the nearest intersection of all 14 edges, only if it is close + // enough to allow this gridVert to safely move to it without inverting + // any tetrahedra. const glm::vec3 delta = edgeVec[j] * (1 - val / (val - gridVert.distance)); const float dist2 = glm::dot(delta, delta); @@ -108,22 +147,25 @@ struct ComputeVerts { gridVert.edgeIndex = i; minDist2 = dist2; } + + // These seven edges are uniquely owned by this gridVert; any of them + // which intersect the surface create a vert. if (i < 7) { const int idx = AtomicAdd(vertIndex, 1); vertPos[idx] = position + spacing * delta; gridVert.edgeVerts[i] = idx; } } - gridVerts.Insert(gridVert); + if (keep) gridVerts.Insert(gridVert); } }; struct BuildTris { glm::ivec3* triVerts; - uint64_t* index; - const HashTable gridVerts; + int* index; + const HashTableD gridVerts; - __host__ __device__ void operator()(uint64_t mortonCode) {} + __host__ __device__ void operator()(int idx) {} }; } // namespace @@ -142,29 +184,33 @@ class SDF { inline Mesh LevelSet(Box bounds, float edgeLength, float level = 0) const { Mesh out; - HashTable gridVerts(10); + glm::vec3 dim = bounds.Size(); // Need to create a new MortonCode function that spreads when needed instead // of always by 3, to make non-cubic domains efficient. float maxDim = std::max(dim[0], std::max(dim[1], dim[2])); - int n = maxDim / edgeLength; - float spacing = maxDim / n; - int maxMorton = MortonCode(glm::ivec3(n)); + glm::ivec3 gridSize = dim / edgeLength; + glm::vec3 spacing = dim / gridSize; + int maxMorton = MortonCode(gridSize); + + HashTable gridVerts(10); // maxMorton^(2/3)? Some heuristic with ability to + // enlarge if it gets too full. VecDH vertPos(gridVerts.Size() * 7); VecDH index(1, 0); thrust::for_each_n( countAt(0), maxMorton, - ComputeVerts({vertPos.ptrD(), index.ptrD(), gridVerts, sdf_, - bounds.Min(), spacing})); - + ComputeVerts({vertPos.ptrD(), index.ptrD(), gridVerts.D(), sdf_, + gridSize, bounds.Min(), spacing})); vertPos.Resize(index.H()[0]); - VecDH triVerts(gridVerts.Entries() * 6); - index.H()[0] = 0; - thrust::for_each_n(countAt(0), maxMorton, - BuildTris({triVerts.ptrD(), index.ptrD(), gridVerts})); + VecDH triVerts(gridVerts.Entries() * 12); // worst case + + index.H()[0] = 0; + thrust::for_each_n( + countAt(0), gridVerts.Size(), + BuildTris({triVerts.ptrD(), index.ptrD(), gridVerts.D()})); triVerts.Resize(index.H()[0]); out.vertPos = vertPos.H(); From 5efb4325981a21849f1ed12fd3487232485cacca Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Mon, 2 May 2022 23:25:21 -0700 Subject: [PATCH 06/39] BuildTris --- sdf/include/sdf.cuh | 133 ++++++++++++++++++++++++++++++++---- utilities/include/structs.h | 5 +- 2 files changed, 124 insertions(+), 14 deletions(-) diff --git a/sdf/include/sdf.cuh b/sdf/include/sdf.cuh index c9ad00f41..98a414331 100644 --- a/sdf/include/sdf.cuh +++ b/sdf/include/sdf.cuh @@ -21,13 +21,57 @@ namespace { using namespace manifold; -constexpr kOpen = std::numeric_limits::max(); +constexpr uint64_t kOpen = std::numeric_limits::max(); + +constexpr glm::ivec3[16] tetTri0 = {{-1, -1, -1}, // + {0, 3, 4}, // + {0, 1, 5}, // + {1, 5, 3}, // + {1, 4, 2}, // + {1, 0, 3}, // + {2, 5, 0}, // + {5, 3, 2}, // + {2, 3, 5}, // + {0, 5, 2}, // + {3, 0, 1}, // + {2, 4, 1}, // + {3, 5, 1}, // + {5, 1, 0}, // + {4, 3, 0}, // + {-1, -1, -1}}; + +constexpr glm::ivec3[16] tetTri1 = {{-1, -1, -1}, // + {-1, -1, -1}, // + {-1, -1, -1}, // + {3, 4, 1}, // + {-1, -1, -1}, // + {3, 2, 1}, // + {0, 4, 2}, // + {-1, -1, -1}, // + {-1, -1, -1}, // + {2, 4, 0}, // + {1, 2, 3}, // + {-1, -1, -1}, // + {1, 4, 3}, // + {-1, -1, -1}, // + {-1, -1, -1}, // + {-1, -1, -1}}; + +constexpr glm::vec3[7] edgeVec = {{0.5f, 0.5f, 0.5f}, // + {1.0f, 0.0f, 0.0f}, // + {0.0f, 1.0f, 0.0f}, // + {0.0f, 0.0f, 1.0f}, // + {-0.5f, 0.5f, 0.5f}, // + {0.5f, -0.5f, 0.5f}, // + {0.5f, 0.5f, -0.5f}}; struct GridVert { uint64_t key = kOpen; float distance = NAN; int edgeIndex = -1; - int[7] edgeVerts; + int[7] edgeVerts = {-1, -1, -1, -1, -1, -1, -1}; + + int Inside() const { return edgeIndex == -1 ? (distance >= 0 ? 1 : -1) : 0 } }; class HashTableD { @@ -40,13 +84,12 @@ class HashTableD { __device__ __host__ bool Insert(const GridVert& vert) { uint32_t idx = vert.key & (size_ - 1); while (1) { - const uint64_t found = - AtomicCAS(&table_.ptrD()[idx].key, kOpen, vert.key); + const uint64_t found = AtomicCAS(&table_[idx].key, kOpen, vert.key); if (found == kOpen) { - if (AtomicAdd(&used_.ptrD()[0], 1) * 2 > size_) { + if (AtomicAdd(&used_[0], 1) * 2 > size_) { return true; } - table_.ptrD()[idx] = vert; + table_[idx] = vert; return false; } if (found == vert.key) return false; @@ -57,13 +100,15 @@ class HashTableD { __device__ __host__ GridVert operator[](uint64_t key) const { uint32_t idx = key & (size_ - 1); while (1) { - const GridVert found = table_.ptrD()[idx]; + const GridVert found = table_[idx]; if (found.key == key) return found; if (found.key == kOpen) return GridVert(); idx = (idx + step_) & (size_ - 1); } } + __device__ __host__ GridVert At(int index) const { return table_[idx]; } + private: const uint32_t step_; VecD table_; @@ -100,10 +145,6 @@ struct ComputeVerts { const glm::ivec3 gridSize; const glm::vec3 origin; const glm::vec3 spacing; - const glm::vec3[7] edgeVec = {{0.5f, 0.5f, 0.5f}, {1.0f, 0.0f, 0.0f}, - {0.0f, 1.0f, 0.0f}, {0.0f, 0.0f, 1.0f}, - {-0.5f, 0.5f, 0.5f}, {0.5f, -0.5f, 0.5f}, - {0.5f, 0.5f, -0.5f}}; inline __host__ __device__ bool AtBounds(glm::ivec3 gridIndex) const { return gridIndex.x == 0 || gridIndex.x == gridSize.x || gridIndex.y == 0 || @@ -162,10 +203,76 @@ struct ComputeVerts { struct BuildTris { glm::ivec3* triVerts; - int* index; + int* triIndex; const HashTableD gridVerts; - __host__ __device__ void operator()(int idx) {} + __host__ __device__ void CreateTris(const glm::ivec4& tet, + const int[6] edges) { + const int i = (tet[0] > 0 ? 1 : 0) + (tet[1] > 0 ? 2 : 0) + + (tet[2] > 0 ? 4 : 0) + (tet[3] > 0 ? 8 : 0); + glm::vec3 tri = tetTri0[i]; + if (tri[0] < 0) return; + int idx = AtomicAdd(triIndex, 1); + triVerts[idx] = {edges[tri[0]], edges[tri[1]], edges[tri[2]]}; + + tri = tetTri1[i]; + if (tri[0] < 0) return; + idx = AtomicAdd(triIndex, 1); + triVerts[idx] = {edges[tri[0]], edges[tri[1]], edges[tri[2]]}; + } + + __host__ __device__ void operator()(int idx) { + const GridVert& base = gridVerts.At(idx); + if (base.key == kOpen) return; + + const glm::ivec4 baseIndex = DecodeMorton(base.key); + + glm::ivec4 leadIndex = baseIndex; + if (leadIndex.w == 0) + leadIndex.w = 1; + else { + leadIndex += 1; + leadIndex.w = 0; + } + const GridVert& leadVert = gridVerts[MortonCode(leadIndex)]; + + // This GridVert is in charge of the 6 tetrahedra surrounding its edge in + // the (1,1,1) direction, attached to leadVert. + glm::ivec4 tet(leadVert.Inside(), base.Inside(), -2, -2); + glm::ivec4 thisIndex = baseIndex; + thisIndex[i] += 1; + GridVert& thisVert = gridVerts[MortonCode(thisIndex)]; + tet[2] = thisVert.Inside(); + int[6] edges = {base.edgeVerts[0], -1, -1, -1, -1, -1}; + for (const int i : {0, 1, 2}) { + edges[1] = base.edgeVerts[i + 1]; + edges[4] = thisVert.edgeVerts[i + 4]; + edges[5] = base.edgeVerts[prev3[i] + 4]; + + thisIndex = leadIndex; + thisIndex[prev3[i]] -= 1; + thisVert = gridVerts[MortonCode(thisIndex)]; + tet[3] = thisVert.Inside(); + edges[2] = thisVert.edgeVerts[next3[i] + 4]; + edges[3] = thisVert.edgeVerts[prev3[i] + 1]; + CreateTris(tet, edges); + + edges[1] = edges[5]; + edges[2] = thisVert.edgeVerts[i + 4]; + edges[4] = edges[3]; + edges[5] = base.edgeVerts[next3[i] + 1]; + + tet[2] = tet[3]; + glm::ivec4 thisIndex = baseIndex; + thisIndex[next3[i]] += 1; + thisVert = gridVerts[MortonCode(thisIndex)]; + tet[3] = thisVert.Inside(); + edges[3] = thisVert.edgeVerts[next3[i] + 4]; + CreateTris(tet, edges); + + tet[2] = tet[3]; + } + } }; } // namespace diff --git a/utilities/include/structs.h b/utilities/include/structs.h index c84d835f4..579cdaa1b 100644 --- a/utilities/include/structs.h +++ b/utilities/include/structs.h @@ -15,13 +15,13 @@ #pragma once #define GLM_FORCE_EXPLICIT_CTOR #include -#include #include #include #include #include #include #include +#include #include #include @@ -68,6 +68,9 @@ void AlwaysAssert(bool condition, const char* file, int line, #define HOST_DEVICE #endif +constexpr glm::ivec3 next3(1, 2, 0); +constexpr glm::ivec3 prev3(2, 0, 1); + inline HOST_DEVICE int Signum(float val) { return (val > 0) - (val < 0); } inline HOST_DEVICE int CCW(glm::vec2 p0, glm::vec2 p1, glm::vec2 p2, From 0dd1f5752fc0decc138a46759ec997b03de61d0b Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Tue, 3 May 2022 23:24:03 -0700 Subject: [PATCH 07/39] fixed some syntax --- sdf/include/sdf.cuh | 56 ++++++++++++++++++++++----------------------- 1 file changed, 28 insertions(+), 28 deletions(-) diff --git a/sdf/include/sdf.cuh b/sdf/include/sdf.cuh index 98a414331..81ca21d8e 100644 --- a/sdf/include/sdf.cuh +++ b/sdf/include/sdf.cuh @@ -23,7 +23,7 @@ using namespace manifold; constexpr uint64_t kOpen = std::numeric_limits::max(); -constexpr glm::ivec3[16] tetTri0 = {{-1, -1, -1}, // +constexpr glm::ivec3 tetTri0[16] = {{-1, -1, -1}, // {0, 3, 4}, // {0, 1, 5}, // {1, 5, 3}, // @@ -40,7 +40,7 @@ constexpr glm::ivec3[16] tetTri0 = {{-1, -1, -1}, // {4, 3, 0}, // {-1, -1, -1}}; -constexpr glm::ivec3[16] tetTri1 = {{-1, -1, -1}, // +constexpr glm::ivec3 tetTri1[16] = {{-1, -1, -1}, // {-1, -1, -1}, // {-1, -1, -1}, // {3, 4, 1}, // @@ -57,7 +57,7 @@ constexpr glm::ivec3[16] tetTri1 = {{-1, -1, -1}, // {-1, -1, -1}, // {-1, -1, -1}}; -constexpr glm::vec3[7] edgeVec = {{0.5f, 0.5f, 0.5f}, // +constexpr glm::vec3 edgeVec[7] = {{0.5f, 0.5f, 0.5f}, // {1.0f, 0.0f, 0.0f}, // {0.0f, 1.0f, 0.0f}, // {0.0f, 0.0f, 1.0f}, // @@ -69,70 +69,70 @@ struct GridVert { uint64_t key = kOpen; float distance = NAN; int edgeIndex = -1; - int[7] edgeVerts = {-1, -1, -1, -1, -1, -1, -1}; + int edgeVerts[7] = {-1, -1, -1, -1, -1, -1, -1}; - int Inside() const { return edgeIndex == -1 ? (distance >= 0 ? 1 : -1) : 0 } + int Inside() const { return edgeIndex == -1 ? (distance >= 0 ? 1 : -1) : 0; } }; class HashTableD { public: - HashTableD(uint32_t step = 127, VecDH& alloc, VecDH& used) + HashTableD(VecDH& alloc, VecDH& used, uint32_t step = 127) : step_{step}, table_{alloc}, used_{used} {} __device__ __host__ int Size() const { return table_.size(); } __device__ __host__ bool Insert(const GridVert& vert) { - uint32_t idx = vert.key & (size_ - 1); + uint32_t idx = vert.key & (Size() - 1); while (1) { const uint64_t found = AtomicCAS(&table_[idx].key, kOpen, vert.key); if (found == kOpen) { - if (AtomicAdd(&used_[0], 1) * 2 > size_) { + if (AtomicAdd(&used_[0], 1) * 2 > Size()) { return true; } table_[idx] = vert; return false; } if (found == vert.key) return false; - idx = (idx + step_) & (size_ - 1); + idx = (idx + step_) & (Size() - 1); } } __device__ __host__ GridVert operator[](uint64_t key) const { - uint32_t idx = key & (size_ - 1); + uint32_t idx = key & (Size() - 1); while (1) { const GridVert found = table_[idx]; if (found.key == key) return found; if (found.key == kOpen) return GridVert(); - idx = (idx + step_) & (size_ - 1); + idx = (idx + step_) & (Size() - 1); } } - __device__ __host__ GridVert At(int index) const { return table_[idx]; } + __device__ __host__ GridVert At(int idx) const { return table_[idx]; } private: const uint32_t step_; VecD table_; - VecD used_(1, 0); + VecD used_; }; class HashTable { public: HashTable(uint32_t sizeExp = 20, uint32_t step = 127) - : alloc_{1 << sizeExp}, table_{step, alloc_, used_} {} + : alloc_{1 << sizeExp}, table_{alloc_, used_, step} {} HashTableD D() { return table_; } int Entries() const { return used_.H()[0]; } - int Size() const { return table_.size(); } + int Size() const { return table_.Size(); } float FilledFraction() const { - return static_cast(used_.H()[0]) / table_.size(); + return static_cast(used_.H()[0]) / table_.Size(); } private: VecDH alloc_; - VecDH used_(1, 0); + VecDH used_ = VecDH(1, 0); HashTableD table_; }; @@ -207,10 +207,10 @@ struct BuildTris { const HashTableD gridVerts; __host__ __device__ void CreateTris(const glm::ivec4& tet, - const int[6] edges) { + const int edges[6]) { const int i = (tet[0] > 0 ? 1 : 0) + (tet[1] > 0 ? 2 : 0) + (tet[2] > 0 ? 4 : 0) + (tet[3] > 0 ? 8 : 0); - glm::vec3 tri = tetTri0[i]; + glm::ivec3 tri = tetTri0[i]; if (tri[0] < 0) return; int idx = AtomicAdd(triIndex, 1); triVerts[idx] = {edges[tri[0]], edges[tri[1]], edges[tri[2]]}; @@ -240,10 +240,10 @@ struct BuildTris { // the (1,1,1) direction, attached to leadVert. glm::ivec4 tet(leadVert.Inside(), base.Inside(), -2, -2); glm::ivec4 thisIndex = baseIndex; - thisIndex[i] += 1; + thisIndex[0] += 1; GridVert& thisVert = gridVerts[MortonCode(thisIndex)]; tet[2] = thisVert.Inside(); - int[6] edges = {base.edgeVerts[0], -1, -1, -1, -1, -1}; + int edges[6] = {base.edgeVerts[0], -1, -1, -1, -1, -1}; for (const int i : {0, 1, 2}) { edges[1] = base.edgeVerts[i + 1]; edges[4] = thisVert.edgeVerts[i + 4]; @@ -296,8 +296,8 @@ class SDF { // Need to create a new MortonCode function that spreads when needed instead // of always by 3, to make non-cubic domains efficient. float maxDim = std::max(dim[0], std::max(dim[1], dim[2])); - glm::ivec3 gridSize = dim / edgeLength; - glm::vec3 spacing = dim / gridSize; + glm::ivec3 gridSize(dim / edgeLength); + glm::vec3 spacing = dim / glm::vec3(gridSize); int maxMorton = MortonCode(gridSize); HashTable gridVerts(10); // maxMorton^(2/3)? Some heuristic with ability to @@ -309,8 +309,8 @@ class SDF { thrust::for_each_n( countAt(0), maxMorton, ComputeVerts({vertPos.ptrD(), index.ptrD(), gridVerts.D(), sdf_, - gridSize, bounds.Min(), spacing})); - vertPos.Resize(index.H()[0]); + gridSize, bounds.min, spacing})); + vertPos.resize(index.H()[0]); VecDH triVerts(gridVerts.Entries() * 12); // worst case @@ -318,10 +318,10 @@ class SDF { thrust::for_each_n( countAt(0), gridVerts.Size(), BuildTris({triVerts.ptrD(), index.ptrD(), gridVerts.D()})); - triVerts.Resize(index.H()[0]); + triVerts.resize(index.H()[0]); - out.vertPos = vertPos.H(); - out.triVerts = triVerts.H(); + out.vertPos.insert(out.vertPos.end(), vertPos.begin(), vertPos.end()); + out.triVerts.insert(out.triVerts.end(), triVerts.begin(), triVerts.end()); return out; } From 085971706c6e1793ea745eb80fba62bc9c66e580 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Tue, 10 May 2022 23:25:00 -0700 Subject: [PATCH 08/39] updating to cpp --- sdf/include/{sdf.cuh => sdf.h} | 4 ++-- tools/CMakeLists.txt | 26 ++++++++++++++++++++------ tools/{sdfTest.cu => sdfTest.cpp} | 2 +- 3 files changed, 23 insertions(+), 9 deletions(-) rename sdf/include/{sdf.cuh => sdf.h} (99%) rename tools/{sdfTest.cu => sdfTest.cpp} (92%) diff --git a/sdf/include/sdf.cuh b/sdf/include/sdf.h similarity index 99% rename from sdf/include/sdf.cuh rename to sdf/include/sdf.h index 81ca21d8e..a670d4ffd 100644 --- a/sdf/include/sdf.cuh +++ b/sdf/include/sdf.h @@ -15,8 +15,8 @@ #pragma once #include "structs.h" -#include "utils.cuh" -#include "vec_dh.cuh" +#include "utils.h" +#include "vec_dh.h" namespace { using namespace manifold; diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index 51946efa5..82092713b 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -29,12 +29,26 @@ endif() # target_compile_options(playground PRIVATE ${MANIFOLD_FLAGS}) # target_compile_features(playground PUBLIC cxx_std_14) -add_executable(sdfTest sdfTest.cu) +add_executable(sdfTest sdfTest.cpp) target_link_libraries(sdfTest manifold meshIO samples sdf) -set_property(TARGET sdfTest PROPERTY CUDA_ARCHITECTURES 61) -target_compile_options(sdfTest PRIVATE ${MANIFOLD_NVCC_FLAGS}) -target_compile_options(sdfTest - PRIVATE "$<$:${MANIFOLD_NVCC_RELEASE_FLAGS}>" "$<$:${MANIFOLD_NVCC_DEBUG_FLAGS}>" -) +if(THRUST_BACKEND STREQUAL "CUDA") + set_source_files_properties( + sdfTest.cpp + + PROPERTIES LANGUAGE CUDA + ) + set_property(TARGET sdfTest PROPERTY CUDA_ARCHITECTURES 61) + target_compile_options(sdfTest + PRIVATE ${MANIFOLD_NVCC_FLAGS} + ) + target_compile_options(sdfTest + PRIVATE "$<$:${MANIFOLD_NVCC_RELEASE_FLAGS}>" "$<$:${MANIFOLD_NVCC_DEBUG_FLAGS}>" + ) +else() + target_compile_options(sdfTest + PRIVATE ${MANIFOLD_FLAGS} + ) +endif() + target_compile_features(sdfTest PUBLIC cxx_std_14) diff --git a/tools/sdfTest.cu b/tools/sdfTest.cpp similarity index 92% rename from tools/sdfTest.cu rename to tools/sdfTest.cpp index a7940cb1a..37e0a4c5c 100644 --- a/tools/sdfTest.cu +++ b/tools/sdfTest.cpp @@ -1,4 +1,4 @@ -#include "sdf.cuh" +#include "sdf.h" using namespace manifold; From 0280a4099d58d9aea97940284bf24d35f9aaf075 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Mon, 13 Jun 2022 23:00:00 -0700 Subject: [PATCH 09/39] added VecDc --- collider/src/collider.cpp | 2 +- sdf/include/sdf.h | 25 +++++++++++++------------ utilities/include/vec_dh.h | 17 +++++++++++++++-- 3 files changed, 29 insertions(+), 15 deletions(-) diff --git a/collider/src/collider.cpp b/collider/src/collider.cpp index 1aefef44a..a0f76805f 100644 --- a/collider/src/collider.cpp +++ b/collider/src/collider.cpp @@ -69,7 +69,7 @@ __host__ __device__ int Leaf2Node(int leaf) { return leaf * 2; } struct CreateRadixTree { int* nodeParent_; thrust::pair* internalChildren_; - const VecD leafMorton_; + const VecDc leafMorton_; __host__ __device__ int PrefixLength(uint32_t a, uint32_t b) const { // count-leading-zeros is used to find the number of identical highest-order diff --git a/sdf/include/sdf.h b/sdf/include/sdf.h index a670d4ffd..8da4c9b3d 100644 --- a/sdf/include/sdf.h +++ b/sdf/include/sdf.h @@ -86,7 +86,7 @@ class HashTableD { while (1) { const uint64_t found = AtomicCAS(&table_[idx].key, kOpen, vert.key); if (found == kOpen) { - if (AtomicAdd(&used_[0], 1) * 2 > Size()) { + if (AtomicAdd(used_[0], 1) * 2 > Size()) { return true; } table_[idx] = vert; @@ -122,12 +122,12 @@ class HashTable { HashTableD D() { return table_; } - int Entries() const { return used_.H()[0]; } + int Entries() const { return used_[0]; } int Size() const { return table_.Size(); } float FilledFraction() const { - return static_cast(used_.H()[0]) / table_.Size(); + return static_cast(used_[0]) / table_.Size(); } private: @@ -140,7 +140,7 @@ template struct ComputeVerts { glm::vec3* vertPos; int* vertIndex; - HashTable gridVerts; + HashTableD gridVerts; const Func sdf; const glm::ivec3 gridSize; const glm::vec3 origin; @@ -153,7 +153,8 @@ struct ComputeVerts { } inline __host__ __device__ float Sdf(glm::vec3 base, int i) const { - sdf(base + (i < 7 ? 1 : -1) * spacing * edgeVec[j]) + return sdf(base + + (i < 7 ? 1.0f : -1.0f) * spacing * edgeVec[i < 7 ? i : i - 7]); } // inline __host__ __device__ float BoundedSdf(glm::vec3 base, int i) const {} @@ -161,7 +162,7 @@ struct ComputeVerts { inline __host__ __device__ void operator()(uint64_t mortonCode) { GridVert gridVert; gridVert.key = mortonCode; - const glm::ivec3 gridIndex = DecodeMorton(mortonCode); + const glm::vec3 gridIndex = DecodeMorton(mortonCode); // const auto sdfFunc = // AtBounds(gridIndex) ? &ComputeVerts::BoundedSdf : &ComputeVerts::Sdf; @@ -192,7 +193,7 @@ struct ComputeVerts { // These seven edges are uniquely owned by this gridVert; any of them // which intersect the surface create a vert. if (i < 7) { - const int idx = AtomicAdd(vertIndex, 1); + const int idx = AtomicAdd(*vertIndex, 1); vertPos[idx] = position + spacing * delta; gridVert.edgeVerts[i] = idx; } @@ -212,12 +213,12 @@ struct BuildTris { (tet[2] > 0 ? 4 : 0) + (tet[3] > 0 ? 8 : 0); glm::ivec3 tri = tetTri0[i]; if (tri[0] < 0) return; - int idx = AtomicAdd(triIndex, 1); + int idx = AtomicAdd(*triIndex, 1); triVerts[idx] = {edges[tri[0]], edges[tri[1]], edges[tri[2]]}; tri = tetTri1[i]; if (tri[0] < 0) return; - idx = AtomicAdd(triIndex, 1); + idx = AtomicAdd(*triIndex, 1); triVerts[idx] = {edges[tri[0]], edges[tri[1]], edges[tri[2]]}; } @@ -310,15 +311,15 @@ class SDF { countAt(0), maxMorton, ComputeVerts({vertPos.ptrD(), index.ptrD(), gridVerts.D(), sdf_, gridSize, bounds.min, spacing})); - vertPos.resize(index.H()[0]); + vertPos.resize(index[0]); VecDH triVerts(gridVerts.Entries() * 12); // worst case - index.H()[0] = 0; + index[0] = 0; thrust::for_each_n( countAt(0), gridVerts.Size(), BuildTris({triVerts.ptrD(), index.ptrD(), gridVerts.D()})); - triVerts.resize(index.H()[0]); + triVerts.resize(index[0]); out.vertPos.insert(out.vertPos.end(), vertPos.begin(), vertPos.end()); out.triVerts.insert(out.triVerts.end(), triVerts.begin(), triVerts.end()); diff --git a/utilities/include/vec_dh.h b/utilities/include/vec_dh.h index 0c3b761b9..aa515a20e 100644 --- a/utilities/include/vec_dh.h +++ b/utilities/include/vec_dh.h @@ -418,9 +418,9 @@ class VecDH { }; template -class VecD { +class VecDc { public: - VecD(const VecDH &vec) : ptr_(vec.ptrD()), size_(vec.size()) {} + VecDc(const VecDH &vec) : ptr_(vec.ptrD()), size_(vec.size()) {} __host__ __device__ const T &operator[](int i) const { return ptr_[i]; } __host__ __device__ int size() const { return size_; } @@ -429,5 +429,18 @@ class VecD { T const *const ptr_; const int size_; }; + +template +class VecD { + public: + VecD(VecDH &vec) : ptr_(vec.ptrD()), size_(vec.size()) {} + + __host__ __device__ T &operator[](int i) const { return ptr_[i]; } + __host__ __device__ int size() const { return size_; } + + private: + T *ptr_; + const int size_; +}; /** @} */ } // namespace manifold From 09cb0cf83bc5f83f4478ab5e8c23bbddbe0aa992 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Sat, 25 Jun 2022 09:20:52 -0700 Subject: [PATCH 10/39] SDF compiles --- sdf/include/sdf.h | 51 +++++++++++++++++++++++++++++++++------ utilities/include/utils.h | 12 +++++++++ 2 files changed, 55 insertions(+), 8 deletions(-) diff --git a/sdf/include/sdf.h b/sdf/include/sdf.h index 8da4c9b3d..c894b7953 100644 --- a/sdf/include/sdf.h +++ b/sdf/include/sdf.h @@ -65,6 +65,43 @@ constexpr glm::vec3 edgeVec[7] = {{0.5f, 0.5f, 0.5f}, // {0.5f, -0.5f, 0.5f}, // {0.5f, 0.5f, -0.5f}}; +__host__ __device__ uint64_t SpreadBits3(uint64_t v) { + v = v & 0x1fffff; + v = (v | v << 32) & 0x1f00000000ffff; + v = (v | v << 16) & 0x1f0000ff0000ff; + v = (v | v << 8) & 0x100f00f00f00f00f; + v = (v | v << 4) & 0x10c30c30c30c30c3; + v = (v | v << 2) & 0x1249249249249249; + return v; +} + +__host__ __device__ uint64_t SqeezeBits3(uint64_t v) { + v = v & 0x1249249249249249; + v = (v ^ (v >> 2)) & 0x10c30c30c30c30c3; + v = (v ^ (v >> 4)) & 0x100f00f00f00f00f; + v = (v ^ (v >> 8)) & 0x1f0000ff0000ff; + v = (v ^ (v >> 16)) & 0x1f00000000ffff; + v = (v ^ (v >> 32)) & 0x1fffff; + return v; +} + +// This is a modified 3D MortonCode, where the xyz code is shifted by one bit +// and the w bit is added as the least significant. This allows 21 bits per x, +// y, and z channel and 1 for w, filling the 64 bit total. +__device__ __host__ uint64_t MortonCode(const glm::ivec4& index) { + return static_cast(index.w) | (SpreadBits3(index.x) << 1) | + (SpreadBits3(index.y) << 2) | (SpreadBits3(index.z) << 3); +} + +__device__ __host__ glm::ivec4 DecodeMorton(uint64_t code) { + glm::ivec4 index; + index.x = SqeezeBits3(code >> 1); + index.y = SqeezeBits3(code >> 2); + index.z = SqeezeBits3(code >> 3); + index.w = code & 0x1u; + return index; +} + struct GridVert { uint64_t key = kOpen; float distance = NAN; @@ -84,9 +121,9 @@ class HashTableD { __device__ __host__ bool Insert(const GridVert& vert) { uint32_t idx = vert.key & (Size() - 1); while (1) { - const uint64_t found = AtomicCAS(&table_[idx].key, kOpen, vert.key); + const uint64_t found = AtomicCAS(table_[idx].key, kOpen, vert.key); if (found == kOpen) { - if (AtomicAdd(used_[0], 1) * 2 > Size()) { + if (AtomicAdd(used_[0], 0x1u) * 2 > Size()) { return true; } table_[idx] = vert; @@ -162,12 +199,12 @@ struct ComputeVerts { inline __host__ __device__ void operator()(uint64_t mortonCode) { GridVert gridVert; gridVert.key = mortonCode; - const glm::vec3 gridIndex = DecodeMorton(mortonCode); + const glm::ivec4 gridIndex = DecodeMorton(mortonCode); // const auto sdfFunc = // AtBounds(gridIndex) ? &ComputeVerts::BoundedSdf : &ComputeVerts::Sdf; - const glm::vec3 position = origin + spacing * gridIndex; + const glm::vec3 position = origin + spacing * glm::vec3(gridIndex); gridVert.distance = sdf(position); bool keep = false; @@ -242,7 +279,7 @@ struct BuildTris { glm::ivec4 tet(leadVert.Inside(), base.Inside(), -2, -2); glm::ivec4 thisIndex = baseIndex; thisIndex[0] += 1; - GridVert& thisVert = gridVerts[MortonCode(thisIndex)]; + GridVert thisVert = gridVerts[MortonCode(thisIndex)]; tet[2] = thisVert.Inside(); int edges[6] = {base.edgeVerts[0], -1, -1, -1, -1, -1}; for (const int i : {0, 1, 2}) { @@ -294,12 +331,10 @@ class SDF { Mesh out; glm::vec3 dim = bounds.Size(); - // Need to create a new MortonCode function that spreads when needed instead - // of always by 3, to make non-cubic domains efficient. float maxDim = std::max(dim[0], std::max(dim[1], dim[2])); glm::ivec3 gridSize(dim / edgeLength); glm::vec3 spacing = dim / glm::vec3(gridSize); - int maxMorton = MortonCode(gridSize); + int maxMorton = MortonCode(glm::ivec4(gridSize, 1)); HashTable gridVerts(10); // maxMorton^(2/3)? Some heuristic with ability to // enlarge if it gets too full. diff --git a/utilities/include/utils.h b/utilities/include/utils.h index 8abf3b163..a15c94a52 100644 --- a/utilities/include/utils.h +++ b/utilities/include/utils.h @@ -108,6 +108,18 @@ thrust::counting_iterator countAt(T i) { return thrust::make_counting_iterator(i); } +template +__host__ __device__ T AtomicCAS(T& target, T compare, T val) { +#ifdef __CUDA_ARCH__ + return atomicCAS(&target, compare, val); +#else + std::atomic& tar = reinterpret_cast&>(target); + T old_val = tar.load(); + tar.compare_exchange_weak(compare, val); + return old_val; +#endif +} + template __host__ __device__ T AtomicAdd(T& target, T add) { #ifdef __CUDA_ARCH__ From 38889736f7596bf51000fb1ccdef73dd2a1fffee Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Sat, 25 Jun 2022 17:02:52 -0700 Subject: [PATCH 11/39] SDF runs --- sdf/include/sdf.h | 41 +++++++++++++++++++++++++++++------------ tools/sdfTest.cpp | 8 ++++++-- 2 files changed, 35 insertions(+), 14 deletions(-) diff --git a/sdf/include/sdf.h b/sdf/include/sdf.h index c894b7953..85bf0d6b6 100644 --- a/sdf/include/sdf.h +++ b/sdf/include/sdf.h @@ -138,8 +138,7 @@ class HashTableD { uint32_t idx = key & (Size() - 1); while (1) { const GridVert found = table_[idx]; - if (found.key == key) return found; - if (found.key == kOpen) return GridVert(); + if (found.key == key || found.key == kOpen) return found; idx = (idx + step_) & (Size() - 1); } } @@ -155,7 +154,7 @@ class HashTableD { class HashTable { public: HashTable(uint32_t sizeExp = 20, uint32_t step = 127) - : alloc_{1 << sizeExp}, table_{alloc_, used_, step} {} + : alloc_{1 << sizeExp, {}}, table_{alloc_, used_, step} {} HashTableD D() { return table_; } @@ -182,6 +181,7 @@ struct ComputeVerts { const glm::ivec3 gridSize; const glm::vec3 origin; const glm::vec3 spacing; + const float level; inline __host__ __device__ bool AtBounds(glm::ivec3 gridIndex) const { return gridIndex.x == 0 || gridIndex.x == gridSize.x || gridIndex.y == 0 || @@ -191,7 +191,8 @@ struct ComputeVerts { inline __host__ __device__ float Sdf(glm::vec3 base, int i) const { return sdf(base + - (i < 7 ? 1.0f : -1.0f) * spacing * edgeVec[i < 7 ? i : i - 7]); + (i < 7 ? 1.0f : -1.0f) * spacing * edgeVec[i < 7 ? i : i - 7]) - + level; } // inline __host__ __device__ float BoundedSdf(glm::vec3 base, int i) const {} @@ -204,8 +205,9 @@ struct ComputeVerts { // const auto sdfFunc = // AtBounds(gridIndex) ? &ComputeVerts::BoundedSdf : &ComputeVerts::Sdf; - const glm::vec3 position = origin + spacing * glm::vec3(gridIndex); - gridVert.distance = sdf(position); + const glm::vec3 position = origin + spacing * glm::vec3(gridIndex) + + (gridIndex.w == 1 ? 0.5f : 0.0f) * glm::vec3(1); + gridVert.distance = sdf(position) - level; bool keep = false; float minDist2 = 0.25 * 0.25; @@ -250,11 +252,15 @@ struct BuildTris { (tet[2] > 0 ? 4 : 0) + (tet[3] > 0 ? 8 : 0); glm::ivec3 tri = tetTri0[i]; if (tri[0] < 0) return; + // for (int i : {0, 1, 2}) + // if (edges[tri[i]] < 0) return; int idx = AtomicAdd(*triIndex, 1); triVerts[idx] = {edges[tri[0]], edges[tri[1]], edges[tri[2]]}; tri = tetTri1[i]; if (tri[0] < 0) return; + // for (int i : {0, 1, 2}) + // if (edges[tri[i]] < 0) return; idx = AtomicAdd(*triIndex, 1); triVerts[idx] = {edges[tri[0]], edges[tri[1]], edges[tri[2]]}; } @@ -272,14 +278,19 @@ struct BuildTris { leadIndex += 1; leadIndex.w = 0; } + const GridVert& leadVert = gridVerts[MortonCode(leadIndex)]; + if (leadVert.key == kOpen) return; // This GridVert is in charge of the 6 tetrahedra surrounding its edge in // the (1,1,1) direction, attached to leadVert. glm::ivec4 tet(leadVert.Inside(), base.Inside(), -2, -2); glm::ivec4 thisIndex = baseIndex; - thisIndex[0] += 1; + thisIndex.x += 1; + GridVert thisVert = gridVerts[MortonCode(thisIndex)]; + bool skipTet = thisVert.key == kOpen; + tet[2] = thisVert.Inside(); int edges[6] = {base.edgeVerts[0], -1, -1, -1, -1, -1}; for (const int i : {0, 1, 2}) { @@ -290,23 +301,29 @@ struct BuildTris { thisIndex = leadIndex; thisIndex[prev3[i]] -= 1; thisVert = gridVerts[MortonCode(thisIndex)]; + tet[3] = thisVert.Inside(); edges[2] = thisVert.edgeVerts[next3[i] + 4]; edges[3] = thisVert.edgeVerts[prev3[i] + 1]; - CreateTris(tet, edges); + if (!skipTet && thisVert.key != kOpen) CreateTris(tet, edges); + skipTet = thisVert.key == kOpen; + + tet[2] = tet[3]; edges[1] = edges[5]; edges[2] = thisVert.edgeVerts[i + 4]; edges[4] = edges[3]; edges[5] = base.edgeVerts[next3[i] + 1]; - tet[2] = tet[3]; glm::ivec4 thisIndex = baseIndex; thisIndex[next3[i]] += 1; thisVert = gridVerts[MortonCode(thisIndex)]; + tet[3] = thisVert.Inside(); edges[3] = thisVert.edgeVerts[next3[i] + 4]; - CreateTris(tet, edges); + + if (!skipTet && thisVert.key != kOpen) CreateTris(tet, edges); + skipTet = thisVert.key == kOpen; tet[2] = tet[3]; } @@ -337,7 +354,7 @@ class SDF { int maxMorton = MortonCode(glm::ivec4(gridSize, 1)); HashTable gridVerts(10); // maxMorton^(2/3)? Some heuristic with ability to - // enlarge if it gets too full. + // enlarge if it gets too full. VecDH vertPos(gridVerts.Size() * 7); VecDH index(1, 0); @@ -345,7 +362,7 @@ class SDF { thrust::for_each_n( countAt(0), maxMorton, ComputeVerts({vertPos.ptrD(), index.ptrD(), gridVerts.D(), sdf_, - gridSize, bounds.min, spacing})); + gridSize, bounds.min, spacing, level})); vertPos.resize(index[0]); VecDH triVerts(gridVerts.Entries() * 12); // worst case diff --git a/tools/sdfTest.cpp b/tools/sdfTest.cpp index 37e0a4c5c..3f1c40573 100644 --- a/tools/sdfTest.cpp +++ b/tools/sdfTest.cpp @@ -1,3 +1,4 @@ +#include "meshIO.h" #include "sdf.h" using namespace manifold; @@ -11,6 +12,9 @@ struct Test { int main(int argc, char **argv) { Test func; SDF a(func); - Box box; - Mesh b = a.LevelSet(box, 1); + Box box({-1, -1, -1}, {1, 1, 1}); + Mesh b = a.LevelSet(box, 1, 0.1); + Dump(b.vertPos); + Dump(b.triVerts); + ExportMesh("sdf.gltf", b, {}); } \ No newline at end of file From 245338a8794e91401b5cf86c24dafbc3094b8d72 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Sat, 25 Jun 2022 23:45:45 -0700 Subject: [PATCH 12/39] fixed some bugs --- sdf/include/sdf.h | 39 ++++++++++++++++++++++----------------- tools/sdfTest.cpp | 6 +++--- 2 files changed, 25 insertions(+), 20 deletions(-) diff --git a/sdf/include/sdf.h b/sdf/include/sdf.h index 85bf0d6b6..b587a606d 100644 --- a/sdf/include/sdf.h +++ b/sdf/include/sdf.h @@ -191,29 +191,33 @@ struct ComputeVerts { inline __host__ __device__ float Sdf(glm::vec3 base, int i) const { return sdf(base + - (i < 7 ? 1.0f : -1.0f) * spacing * edgeVec[i < 7 ? i : i - 7]) - - level; + (i < 7 ? 1.0f : -1.0f) * spacing * edgeVec[i < 7 ? i : i - 7]); } // inline __host__ __device__ float BoundedSdf(glm::vec3 base, int i) const {} inline __host__ __device__ void operator()(uint64_t mortonCode) { + const glm::ivec4 gridIndex = DecodeMorton(mortonCode); + if (gridIndex.x > gridSize.x || gridIndex.y > gridSize.y || + gridIndex.z > gridSize.z) + return; + GridVert gridVert; gridVert.key = mortonCode; - const glm::ivec4 gridIndex = DecodeMorton(mortonCode); // const auto sdfFunc = // AtBounds(gridIndex) ? &ComputeVerts::BoundedSdf : &ComputeVerts::Sdf; - const glm::vec3 position = origin + spacing * glm::vec3(gridIndex) + - (gridIndex.w == 1 ? 0.5f : 0.0f) * glm::vec3(1); + const glm::vec3 position = + origin + spacing * (glm::vec3(gridIndex) + + (gridIndex.w == 1 ? 0.5f : 0.0f) * glm::vec3(1)); gridVert.distance = sdf(position) - level; bool keep = false; float minDist2 = 0.25 * 0.25; for (int i = 0; i < 14; ++i) { const int j = i < 7 ? i : i - 7; - const float val = Sdf(position, i); + const float val = Sdf(position, i) - level; if (val * gridVert.distance > 0 || (val == 0 && gridVert.distance == 0)) continue; keep = true; @@ -252,15 +256,11 @@ struct BuildTris { (tet[2] > 0 ? 4 : 0) + (tet[3] > 0 ? 8 : 0); glm::ivec3 tri = tetTri0[i]; if (tri[0] < 0) return; - // for (int i : {0, 1, 2}) - // if (edges[tri[i]] < 0) return; int idx = AtomicAdd(*triIndex, 1); triVerts[idx] = {edges[tri[0]], edges[tri[1]], edges[tri[2]]}; tri = tetTri1[i]; if (tri[0] < 0) return; - // for (int i : {0, 1, 2}) - // if (edges[tri[i]] < 0) return; idx = AtomicAdd(*triIndex, 1); triVerts[idx] = {edges[tri[0]], edges[tri[1]], edges[tri[2]]}; } @@ -315,7 +315,7 @@ struct BuildTris { edges[4] = edges[3]; edges[5] = base.edgeVerts[next3[i] + 1]; - glm::ivec4 thisIndex = baseIndex; + thisIndex = baseIndex; thisIndex[next3[i]] += 1; thisVert = gridVerts[MortonCode(thisIndex)]; @@ -347,11 +347,16 @@ class SDF { inline Mesh LevelSet(Box bounds, float edgeLength, float level = 0) const { Mesh out; - glm::vec3 dim = bounds.Size(); - float maxDim = std::max(dim[0], std::max(dim[1], dim[2])); - glm::ivec3 gridSize(dim / edgeLength); - glm::vec3 spacing = dim / glm::vec3(gridSize); - int maxMorton = MortonCode(glm::ivec4(gridSize, 1)); + const glm::vec3 dim = bounds.Size(); + const float maxDim = std::max(dim[0], std::max(dim[1], dim[2])); + const glm::ivec3 gridSize(dim / edgeLength); // + 1.0f); + const glm::vec3 spacing = dim / (glm::vec3(gridSize)); // - 0.5f); + + const glm::vec3 origin = bounds.min; // - 0.5f * spacing; + std::cout << spacing << ", " << origin << ", " + << origin + glm::vec3(gridSize) * spacing + 0.5f * spacing + << std::endl; + const int maxMorton = MortonCode(glm::ivec4(gridSize, 1)); HashTable gridVerts(10); // maxMorton^(2/3)? Some heuristic with ability to // enlarge if it gets too full. @@ -362,7 +367,7 @@ class SDF { thrust::for_each_n( countAt(0), maxMorton, ComputeVerts({vertPos.ptrD(), index.ptrD(), gridVerts.D(), sdf_, - gridSize, bounds.min, spacing, level})); + gridSize, origin, spacing, level})); vertPos.resize(index[0]); VecDH triVerts(gridVerts.Entries() * 12); // worst case diff --git a/tools/sdfTest.cpp b/tools/sdfTest.cpp index 3f1c40573..06bebb9c3 100644 --- a/tools/sdfTest.cpp +++ b/tools/sdfTest.cpp @@ -13,8 +13,8 @@ int main(int argc, char **argv) { Test func; SDF a(func); Box box({-1, -1, -1}, {1, 1, 1}); - Mesh b = a.LevelSet(box, 1, 0.1); - Dump(b.vertPos); - Dump(b.triVerts); + Mesh b = a.LevelSet(box, 0.5, 0.6); + // Dump(b.vertPos); + // Dump(b.triVerts); ExportMesh("sdf.gltf", b, {}); } \ No newline at end of file From 997fb73406f3520e9d775b2c8643e25af6420fc5 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Sun, 26 Jun 2022 10:42:50 -0700 Subject: [PATCH 13/39] manifold output --- sdf/include/sdf.h | 84 +++++++++++++++++++------------------ tools/sdfTest.cpp | 6 +-- utilities/include/structs.h | 8 ++++ 3 files changed, 54 insertions(+), 44 deletions(-) diff --git a/sdf/include/sdf.h b/sdf/include/sdf.h index b587a606d..f4de3e71e 100644 --- a/sdf/include/sdf.h +++ b/sdf/include/sdf.h @@ -108,7 +108,7 @@ struct GridVert { int edgeIndex = -1; int edgeVerts[7] = {-1, -1, -1, -1, -1, -1, -1}; - int Inside() const { return edgeIndex == -1 ? (distance >= 0 ? 1 : -1) : 0; } + int Inside() const { return edgeIndex == -1 ? (distance > 0 ? 1 : -1) : 0; } }; class HashTableD { @@ -178,46 +178,44 @@ struct ComputeVerts { int* vertIndex; HashTableD gridVerts; const Func sdf; - const glm::ivec3 gridSize; - const glm::vec3 origin; + const Box bounds; + const Box innerBounds; + const Box outerBounds; const glm::vec3 spacing; const float level; - inline __host__ __device__ bool AtBounds(glm::ivec3 gridIndex) const { - return gridIndex.x == 0 || gridIndex.x == gridSize.x || gridIndex.y == 0 || - gridIndex.y == gridSize.y || gridIndex.z == 0 || - gridIndex.z == gridSize.z; + inline __host__ __device__ float BoundedSDF(glm::vec3 pos) const { + const float d = sdf(pos) - level; + return innerBounds.Contains(pos) ? d : glm::min(d, 0.0f); } - inline __host__ __device__ float Sdf(glm::vec3 base, int i) const { - return sdf(base + - (i < 7 ? 1.0f : -1.0f) * spacing * edgeVec[i < 7 ? i : i - 7]); + inline __host__ __device__ float AdjacentSDF(glm::vec3 base, int i) const { + return BoundedSDF(base + (i < 7 ? 1.0f : -1.0f) * spacing * + edgeVec[i < 7 ? i : i - 7]); } - // inline __host__ __device__ float BoundedSdf(glm::vec3 base, int i) const {} - inline __host__ __device__ void operator()(uint64_t mortonCode) { const glm::ivec4 gridIndex = DecodeMorton(mortonCode); - if (gridIndex.x > gridSize.x || gridIndex.y > gridSize.y || - gridIndex.z > gridSize.z) - return; - - GridVert gridVert; - gridVert.key = mortonCode; // const auto sdfFunc = // AtBounds(gridIndex) ? &ComputeVerts::BoundedSdf : &ComputeVerts::Sdf; const glm::vec3 position = - origin + spacing * (glm::vec3(gridIndex) + - (gridIndex.w == 1 ? 0.5f : 0.0f) * glm::vec3(1)); - gridVert.distance = sdf(position) - level; + bounds.min + + spacing * (-0.5f + glm::vec3(gridIndex) + + (gridIndex.w == 1 ? 0.5f : 0.0f) * glm::vec3(1)); + + if (!outerBounds.Contains(position)) return; + + GridVert gridVert; + gridVert.key = mortonCode; + gridVert.distance = BoundedSDF(position); bool keep = false; float minDist2 = 0.25 * 0.25; for (int i = 0; i < 14; ++i) { const int j = i < 7 ? i : i - 7; - const float val = Sdf(position, i) - level; + const float val = AdjacentSDF(position, i); if (val * gridVert.distance > 0 || (val == 0 && gridVert.distance == 0)) continue; keep = true; @@ -229,7 +227,7 @@ struct ComputeVerts { edgeVec[j] * (1 - val / (val - gridVert.distance)); const float dist2 = glm::dot(delta, delta); if (dist2 < minDist2) { - gridVert.edgeIndex = i; + // gridVert.edgeIndex = i; minDist2 = dist2; } @@ -250,19 +248,19 @@ struct BuildTris { int* triIndex; const HashTableD gridVerts; - __host__ __device__ void CreateTris(const glm::ivec4& tet, - const int edges[6]) { - const int i = (tet[0] > 0 ? 1 : 0) + (tet[1] > 0 ? 2 : 0) + - (tet[2] > 0 ? 4 : 0) + (tet[3] > 0 ? 8 : 0); - glm::ivec3 tri = tetTri0[i]; + __host__ __device__ void CreateTri(const glm::ivec3& tri, + const int edges[6]) { if (tri[0] < 0) return; int idx = AtomicAdd(*triIndex, 1); triVerts[idx] = {edges[tri[0]], edges[tri[1]], edges[tri[2]]}; + } - tri = tetTri1[i]; - if (tri[0] < 0) return; - idx = AtomicAdd(*triIndex, 1); - triVerts[idx] = {edges[tri[0]], edges[tri[1]], edges[tri[2]]}; + __host__ __device__ void CreateTris(const glm::ivec4& tet, + const int edges[6]) { + const int i = (tet[0] > 0 ? 1 : 0) + (tet[1] > 0 ? 2 : 0) + + (tet[2] > 0 ? 4 : 0) + (tet[3] > 0 ? 8 : 0); + CreateTri(tetTri0[i], edges); + CreateTri(tetTri1[i], edges); } __host__ __device__ void operator()(int idx) { @@ -349,14 +347,18 @@ class SDF { const glm::vec3 dim = bounds.Size(); const float maxDim = std::max(dim[0], std::max(dim[1], dim[2])); - const glm::ivec3 gridSize(dim / edgeLength); // + 1.0f); - const glm::vec3 spacing = dim / (glm::vec3(gridSize)); // - 0.5f); + const glm::ivec3 gridSize(dim / edgeLength); + const glm::vec3 spacing = dim / (glm::vec3(gridSize)); + + const Box innerBounds(bounds.min + 0.25f * spacing, + bounds.max - 0.25f * spacing); + const Box outerBounds(bounds.min - 0.75f * spacing, + bounds.max + 0.75f * spacing); - const glm::vec3 origin = bounds.min; // - 0.5f * spacing; - std::cout << spacing << ", " << origin << ", " - << origin + glm::vec3(gridSize) * spacing + 0.5f * spacing - << std::endl; - const int maxMorton = MortonCode(glm::ivec4(gridSize, 1)); + const int maxMorton = MortonCode(glm::ivec4(gridSize + 1, 0)); + // const int maxSize = glm::max(gridSize.x, glm::max(gridSize.y, + // gridSize.z)); const int maxMorton = + // MortonCode(glm::ivec4(glm::ivec3(maxSize + 1), 0)); HashTable gridVerts(10); // maxMorton^(2/3)? Some heuristic with ability to // enlarge if it gets too full. @@ -365,9 +367,9 @@ class SDF { VecDH index(1, 0); thrust::for_each_n( - countAt(0), maxMorton, + countAt(0), maxMorton + 1, ComputeVerts({vertPos.ptrD(), index.ptrD(), gridVerts.D(), sdf_, - gridSize, origin, spacing, level})); + bounds, innerBounds, outerBounds, spacing, level})); vertPos.resize(index[0]); VecDH triVerts(gridVerts.Entries() * 12); // worst case diff --git a/tools/sdfTest.cpp b/tools/sdfTest.cpp index 06bebb9c3..9ab2561a5 100644 --- a/tools/sdfTest.cpp +++ b/tools/sdfTest.cpp @@ -5,15 +5,15 @@ using namespace manifold; struct Test { __host__ __device__ float operator()(glm::vec3 point) const { - return point[1]; + return glm::dot(point, {1, 1, 1}); } }; int main(int argc, char **argv) { Test func; SDF a(func); - Box box({-1, -1, -1}, {1, 1, 1}); - Mesh b = a.LevelSet(box, 0.5, 0.6); + Box box({-1, -1, -1}, {1, 2, 1}); + Mesh b = a.LevelSet(box, 0.5, -0.1); // Dump(b.vertPos); // Dump(b.triVerts); ExportMesh("sdf.gltf", b, {}); diff --git a/utilities/include/structs.h b/utilities/include/structs.h index e2137e697..f069f68ea 100644 --- a/utilities/include/structs.h +++ b/utilities/include/structs.h @@ -314,6 +314,14 @@ struct Box { return glm::max(absMax.x, glm::max(absMax.y, absMax.z)); } + /** + * Does this box contain (includes equal) the given point? + */ + HOST_DEVICE bool Contains(const glm::vec3& p) const { + return glm::all(glm::greaterThanEqual(p, min)) && + glm::all(glm::greaterThanEqual(max, p)); + } + /** * Does this box contain (includes equal) the given box? */ From f4eeca739ff429e0c8c2ba11a1a14b98dc44f7a1 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Mon, 27 Jun 2022 13:30:56 -0700 Subject: [PATCH 14/39] fixed out of memory bug --- sdf/include/sdf.h | 16 +++++++++------- tools/sdfTest.cpp | 2 +- 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/sdf/include/sdf.h b/sdf/include/sdf.h index f4de3e71e..037fefef9 100644 --- a/sdf/include/sdf.h +++ b/sdf/include/sdf.h @@ -239,7 +239,9 @@ struct ComputeVerts { gridVert.edgeVerts[i] = idx; } } - if (keep) gridVerts.Insert(gridVert); + + if (keep && gridVerts.Insert(gridVert)) + std::cout << "out of space!" << std::endl; } }; @@ -355,13 +357,13 @@ class SDF { const Box outerBounds(bounds.min - 0.75f * spacing, bounds.max + 0.75f * spacing); - const int maxMorton = MortonCode(glm::ivec4(gridSize + 1, 0)); - // const int maxSize = glm::max(gridSize.x, glm::max(gridSize.y, - // gridSize.z)); const int maxMorton = - // MortonCode(glm::ivec4(glm::ivec3(maxSize + 1), 0)); + const int maxMorton = MortonCode(glm::ivec4(gridSize + 1, 1)); - HashTable gridVerts(10); // maxMorton^(2/3)? Some heuristic with ability to - // enlarge if it gets too full. + const int tableSize = + glm::min(maxMorton, static_cast(100 * glm::pow(maxMorton, 0.667))); + std::cout << "maxMorton: " << maxMorton + << ", hash table size: " << tableSize << std::endl; + HashTable gridVerts(tableSize); VecDH vertPos(gridVerts.Size() * 7); VecDH index(1, 0); diff --git a/tools/sdfTest.cpp b/tools/sdfTest.cpp index 9ab2561a5..980176092 100644 --- a/tools/sdfTest.cpp +++ b/tools/sdfTest.cpp @@ -12,7 +12,7 @@ struct Test { int main(int argc, char **argv) { Test func; SDF a(func); - Box box({-1, -1, -1}, {1, 2, 1}); + Box box({-1, -1, -1}, {1, 2, 3}); Mesh b = a.LevelSet(box, 0.5, -0.1); // Dump(b.vertPos); // Dump(b.triVerts); From 941fd62bc9905698cfc6e8ae5f0062f9c4a0907f Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Wed, 29 Jun 2022 23:06:08 -0700 Subject: [PATCH 15/39] flattened sides --- sdf/include/sdf.h | 122 ++++++++++++++++++++++++++++++---------------- tools/sdfTest.cpp | 2 +- 2 files changed, 81 insertions(+), 43 deletions(-) diff --git a/sdf/include/sdf.h b/sdf/include/sdf.h index 037fefef9..6639b80dc 100644 --- a/sdf/include/sdf.h +++ b/sdf/include/sdf.h @@ -57,13 +57,20 @@ constexpr glm::ivec3 tetTri1[16] = {{-1, -1, -1}, // {-1, -1, -1}, // {-1, -1, -1}}; -constexpr glm::vec3 edgeVec[7] = {{0.5f, 0.5f, 0.5f}, // - {1.0f, 0.0f, 0.0f}, // - {0.0f, 1.0f, 0.0f}, // - {0.0f, 0.0f, 1.0f}, // - {-0.5f, 0.5f, 0.5f}, // - {0.5f, -0.5f, 0.5f}, // - {0.5f, 0.5f, -0.5f}}; +constexpr glm::ivec4 neighbors[14] = {{0, 0, 0, 1}, // + {1, 0, 0, 0}, // + {0, 1, 0, 0}, // + {0, 0, 1, 0}, // + {-1, 0, 0, 1}, // + {0, -1, 0, 1}, // + {0, 0, -1, 1}, // + {-1, -1, -1, 1}, // + {-1, 0, 0, 0}, // + {0, -1, 0, 0}, // + {0, 0, -1, 0}, // + {0, -1, -1, 1}, // + {-1, 0, -1, 1}, // + {-1, -1, 0, 1}}; __host__ __device__ uint64_t SpreadBits3(uint64_t v) { v = v & 0x1fffff; @@ -77,11 +84,11 @@ __host__ __device__ uint64_t SpreadBits3(uint64_t v) { __host__ __device__ uint64_t SqeezeBits3(uint64_t v) { v = v & 0x1249249249249249; - v = (v ^ (v >> 2)) & 0x10c30c30c30c30c3; - v = (v ^ (v >> 4)) & 0x100f00f00f00f00f; - v = (v ^ (v >> 8)) & 0x1f0000ff0000ff; - v = (v ^ (v >> 16)) & 0x1f00000000ffff; - v = (v ^ (v >> 32)) & 0x1fffff; + v = (v ^ v >> 2) & 0x10c30c30c30c30c3; + v = (v ^ v >> 4) & 0x100f00f00f00f00f; + v = (v ^ v >> 8) & 0x1f0000ff0000ff; + v = (v ^ v >> 16) & 0x1f00000000ffff; + v = (v ^ v >> 32) & 0x1fffff; return v; } @@ -153,8 +160,8 @@ class HashTableD { class HashTable { public: - HashTable(uint32_t sizeExp = 20, uint32_t step = 127) - : alloc_{1 << sizeExp, {}}, table_{alloc_, used_, step} {} + HashTable(uint32_t size, uint32_t step = 127) + : alloc_{1 << (int)ceil(log2(size)), {}}, table_{alloc_, used_, step} {} HashTableD D() { return table_; } @@ -178,44 +185,79 @@ struct ComputeVerts { int* vertIndex; HashTableD gridVerts; const Func sdf; - const Box bounds; - const Box innerBounds; - const Box outerBounds; + const glm::vec3 origin; + const glm::ivec3 gridSize; const glm::vec3 spacing; const float level; - inline __host__ __device__ float BoundedSDF(glm::vec3 pos) const { - const float d = sdf(pos) - level; - return innerBounds.Contains(pos) ? d : glm::min(d, 0.0f); + inline __host__ __device__ glm::vec3 Position(glm::ivec4 gridIndex) const { + return origin + spacing * (-0.5f + glm::vec3(gridIndex) + + (gridIndex.w == 1 ? 0.5f : 0.0f) * glm::vec3(1)); } - inline __host__ __device__ float AdjacentSDF(glm::vec3 base, int i) const { - return BoundedSDF(base + (i < 7 ? 1.0f : -1.0f) * spacing * - edgeVec[i < 7 ? i : i - 7]); + inline __host__ __device__ float BoundedSDF(glm::ivec4 gridIndex) const { + glm::vec3 pos = Position(gridIndex); + + const float d = sdf(pos) - level; + + if (gridIndex.w == 1) { + if (d > 0) { + const glm::ivec3 xyz(gridIndex); + if (glm::any(glm::equal(xyz, glm::ivec3(0))) || + glm::any(glm::equal(xyz, gridSize - 1))) + return 0.0f; + else if (glm::any(glm::equal(xyz, gridSize))) + return -1.0f; + } + } else { + bool mirror = true; + if (gridIndex.x == 0) + pos.x += spacing.x; + else if (gridIndex.y == 0) + pos.y += spacing.y; + else if (gridIndex.z == 0) + pos.z += spacing.z; + else if (gridIndex.x == gridSize.x) + pos.x -= spacing.x; + else if (gridIndex.y == gridSize.y) + pos.y -= spacing.y; + else if (gridIndex.z == gridSize.z) + pos.z -= spacing.z; + else + mirror = false; + if (mirror) { + const float dMirror = sdf(pos) - level; + return dMirror > 0 ? glm::min(d, -dMirror) : -1.0f; + } + } + return d; } inline __host__ __device__ void operator()(uint64_t mortonCode) { const glm::ivec4 gridIndex = DecodeMorton(mortonCode); + if (glm::any(glm::greaterThan(glm::ivec3(gridIndex), gridSize))) return; + // const auto sdfFunc = // AtBounds(gridIndex) ? &ComputeVerts::BoundedSdf : &ComputeVerts::Sdf; - const glm::vec3 position = - bounds.min + - spacing * (-0.5f + glm::vec3(gridIndex) + - (gridIndex.w == 1 ? 0.5f : 0.0f) * glm::vec3(1)); - - if (!outerBounds.Contains(position)) return; + const glm::vec3 position = Position(gridIndex); GridVert gridVert; gridVert.key = mortonCode; - gridVert.distance = BoundedSDF(position); + gridVert.distance = BoundedSDF(gridIndex); bool keep = false; float minDist2 = 0.25 * 0.25; for (int i = 0; i < 14; ++i) { const int j = i < 7 ? i : i - 7; - const float val = AdjacentSDF(position, i); + glm::ivec4 neighborIndex = gridIndex + neighbors[i]; + if (neighborIndex.w == 2) { + neighborIndex += 1; + neighborIndex.w = 0; + } + const glm::vec3 iPos = Position(neighborIndex); + const float val = BoundedSDF(neighborIndex); if (val * gridVert.distance > 0 || (val == 0 && gridVert.distance == 0)) continue; keep = true; @@ -223,8 +265,9 @@ struct ComputeVerts { // Record the nearest intersection of all 14 edges, only if it is close // enough to allow this gridVert to safely move to it without inverting // any tetrahedra. - const glm::vec3 delta = - edgeVec[j] * (1 - val / (val - gridVert.distance)); + const glm::vec3 interp = (val * position - gridVert.distance * iPos) / + (val - gridVert.distance); + const glm::vec3 delta = interp - position; const float dist2 = glm::dot(delta, delta); if (dist2 < minDist2) { // gridVert.edgeIndex = i; @@ -235,7 +278,7 @@ struct ComputeVerts { // which intersect the surface create a vert. if (i < 7) { const int idx = AtomicAdd(*vertIndex, 1); - vertPos[idx] = position + spacing * delta; + vertPos[idx] = interp; gridVert.edgeVerts[i] = idx; } } @@ -352,15 +395,10 @@ class SDF { const glm::ivec3 gridSize(dim / edgeLength); const glm::vec3 spacing = dim / (glm::vec3(gridSize)); - const Box innerBounds(bounds.min + 0.25f * spacing, - bounds.max - 0.25f * spacing); - const Box outerBounds(bounds.min - 0.75f * spacing, - bounds.max + 0.75f * spacing); - const int maxMorton = MortonCode(glm::ivec4(gridSize + 1, 1)); - const int tableSize = - glm::min(maxMorton, static_cast(100 * glm::pow(maxMorton, 0.667))); + const int tableSize = glm::min( + 2 * maxMorton, static_cast(100 * glm::pow(maxMorton, 0.667))); std::cout << "maxMorton: " << maxMorton << ", hash table size: " << tableSize << std::endl; HashTable gridVerts(tableSize); @@ -371,7 +409,7 @@ class SDF { thrust::for_each_n( countAt(0), maxMorton + 1, ComputeVerts({vertPos.ptrD(), index.ptrD(), gridVerts.D(), sdf_, - bounds, innerBounds, outerBounds, spacing, level})); + bounds.min, gridSize + 1, spacing, level})); vertPos.resize(index[0]); VecDH triVerts(gridVerts.Entries() * 12); // worst case diff --git a/tools/sdfTest.cpp b/tools/sdfTest.cpp index 980176092..7774a8f46 100644 --- a/tools/sdfTest.cpp +++ b/tools/sdfTest.cpp @@ -5,7 +5,7 @@ using namespace manifold; struct Test { __host__ __device__ float operator()(glm::vec3 point) const { - return glm::dot(point, {1, 1, 1}); + return glm::dot(point, {1, 1, 2}); } }; From 6dd429a436987f2c5744c065fbd411f43aa511ad Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Mon, 15 Aug 2022 22:02:37 -0700 Subject: [PATCH 16/39] back to egg-crate --- src/sdf/include/sdf.h | 53 ++++++++++++------------------------------- 1 file changed, 14 insertions(+), 39 deletions(-) diff --git a/src/sdf/include/sdf.h b/src/sdf/include/sdf.h index a41ceb443..782f8c4ea 100644 --- a/src/sdf/include/sdf.h +++ b/src/sdf/include/sdf.h @@ -115,7 +115,7 @@ struct GridVert { int edgeIndex = -1; int edgeVerts[7] = {-1, -1, -1, -1, -1, -1, -1}; - int Inside() const { return edgeIndex == -1 ? (distance > 0 ? 1 : -1) : 0; } + int Inside() const { return distance > 0 ? 1 : -1; } }; class HashTableD { @@ -196,40 +196,15 @@ struct ComputeVerts { } inline __host__ __device__ float BoundedSDF(glm::ivec4 gridIndex) const { - glm::vec3 pos = Position(gridIndex); - - const float d = sdf(pos) - level; - - if (gridIndex.w == 1) { - if (d > 0) { - const glm::ivec3 xyz(gridIndex); - if (glm::any(glm::equal(xyz, glm::ivec3(0))) || - glm::any(glm::equal(xyz, gridSize - 1))) - return 0.0f; - else if (glm::any(glm::equal(xyz, gridSize))) - return -1.0f; - } - } else { - bool mirror = true; - if (gridIndex.x == 0) - pos.x += spacing.x; - else if (gridIndex.y == 0) - pos.y += spacing.y; - else if (gridIndex.z == 0) - pos.z += spacing.z; - else if (gridIndex.x == gridSize.x) - pos.x -= spacing.x; - else if (gridIndex.y == gridSize.y) - pos.y -= spacing.y; - else if (gridIndex.z == gridSize.z) - pos.z -= spacing.z; - else - mirror = false; - if (mirror) { - const float dMirror = sdf(pos) - level; - return dMirror > 0 ? glm::min(d, -dMirror) : -1.0f; - } - } + const float d = sdf(Position(gridIndex)) - level; + + const glm::ivec3 xyz(gridIndex); + if (glm::any(glm::equal(xyz, glm::ivec3(0))) || + glm::any(glm::greaterThanEqual(xyz, gridSize)) || + (gridIndex.w == 1 && + glm::any(glm::greaterThanEqual(xyz, gridSize - 1)))) + return glm::min(d, 0.0f); + return d; } @@ -283,8 +258,8 @@ struct ComputeVerts { } } - // if (keep && gridVerts.Insert(gridVert)) - // std::cout << "out of space!" << std::endl; + if (keep && gridVerts.Insert(gridVert)) + std::cout << "out of space!" << std::endl; } }; @@ -395,7 +370,7 @@ class SDF { const glm::ivec3 gridSize(dim / edgeLength); const glm::vec3 spacing = dim / (glm::vec3(gridSize)); - const int maxMorton = MortonCode(glm::ivec4(gridSize + 1, 1)); + const int maxMorton = MortonCode(glm::ivec4(gridSize - 1, 1)); const int tableSize = glm::min( 2 * maxMorton, static_cast(100 * glm::pow(maxMorton, 0.667))); @@ -409,7 +384,7 @@ class SDF { thrust::for_each_n( countAt(0), maxMorton + 1, ComputeVerts({vertPos.ptrD(), index.ptrD(), gridVerts.D(), sdf_, - bounds.min, gridSize + 1, spacing, level})); + bounds.min, gridSize - 1, spacing, level})); vertPos.resize(index[0]); VecDH triVerts(gridVerts.Entries() * 12); // worst case From 2c9595032e7dbe9b91ebdda6fcce892b6d95c0bd Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Mon, 15 Aug 2022 23:17:21 -0700 Subject: [PATCH 17/39] works with CUDA too --- src/sdf/include/sdf.h | 186 +++++++++++++++++++------------- src/utilities/include/structs.h | 2 - src/utilities/include/utils.h | 12 --- 3 files changed, 109 insertions(+), 91 deletions(-) diff --git a/src/sdf/include/sdf.h b/src/sdf/include/sdf.h index 782f8c4ea..50cf3322c 100644 --- a/src/sdf/include/sdf.h +++ b/src/sdf/include/sdf.h @@ -21,58 +21,91 @@ namespace { using namespace manifold; -constexpr uint64_t kOpen = std::numeric_limits::max(); - -constexpr glm::ivec3 tetTri0[16] = {{-1, -1, -1}, // - {0, 3, 4}, // - {0, 1, 5}, // - {1, 5, 3}, // - {1, 4, 2}, // - {1, 0, 3}, // - {2, 5, 0}, // - {5, 3, 2}, // - {2, 3, 5}, // - {0, 5, 2}, // - {3, 0, 1}, // - {2, 4, 1}, // - {3, 5, 1}, // - {5, 1, 0}, // - {4, 3, 0}, // - {-1, -1, -1}}; - -constexpr glm::ivec3 tetTri1[16] = {{-1, -1, -1}, // - {-1, -1, -1}, // - {-1, -1, -1}, // - {3, 4, 1}, // - {-1, -1, -1}, // - {3, 2, 1}, // - {0, 4, 2}, // - {-1, -1, -1}, // - {-1, -1, -1}, // - {2, 4, 0}, // - {1, 2, 3}, // - {-1, -1, -1}, // - {1, 4, 3}, // - {-1, -1, -1}, // - {-1, -1, -1}, // - {-1, -1, -1}}; - -constexpr glm::ivec4 neighbors[14] = {{0, 0, 0, 1}, // - {1, 0, 0, 0}, // - {0, 1, 0, 0}, // - {0, 0, 1, 0}, // - {-1, 0, 0, 1}, // - {0, -1, 0, 1}, // - {0, 0, -1, 1}, // - {-1, -1, -1, 1}, // - {-1, 0, 0, 0}, // - {0, -1, 0, 0}, // - {0, 0, -1, 0}, // - {0, -1, -1, 1}, // - {-1, 0, -1, 1}, // - {-1, -1, 0, 1}}; - -__host__ __device__ uint64_t SpreadBits3(uint64_t v) { +typedef unsigned long long int Uint64; + +constexpr Uint64 kOpen = std::numeric_limits::max(); + +__host__ __device__ Uint64 AtomicCAS(Uint64& target, Uint64 compare, + Uint64 val) { +#ifdef __CUDA_ARCH__ + return atomicCAS(&target, compare, val); +#else + std::atomic& tar = reinterpret_cast&>(target); + Uint64 old_val = tar.load(); + tar.compare_exchange_weak(compare, val); + return old_val; +#endif +} + +__host__ __device__ int Next3(int i) { + constexpr glm::ivec3 next3(1, 2, 0); + return next3[i]; +} + +__host__ __device__ int Prev3(int i) { + constexpr glm::ivec3 prev3(2, 0, 1); + return prev3[i]; +} + +__host__ __device__ glm::ivec3 TetTri0(int i) { + constexpr glm::ivec3 tetTri0[16] = {{-1, -1, -1}, // + {0, 3, 4}, // + {0, 1, 5}, // + {1, 5, 3}, // + {1, 4, 2}, // + {1, 0, 3}, // + {2, 5, 0}, // + {5, 3, 2}, // + {2, 3, 5}, // + {0, 5, 2}, // + {3, 0, 1}, // + {2, 4, 1}, // + {3, 5, 1}, // + {5, 1, 0}, // + {4, 3, 0}, // + {-1, -1, -1}}; + return tetTri0[i]; +} + +__host__ __device__ glm::ivec3 TetTri1(int i) { + constexpr glm::ivec3 tetTri1[16] = {{-1, -1, -1}, // + {-1, -1, -1}, // + {-1, -1, -1}, // + {3, 4, 1}, // + {-1, -1, -1}, // + {3, 2, 1}, // + {0, 4, 2}, // + {-1, -1, -1}, // + {-1, -1, -1}, // + {2, 4, 0}, // + {1, 2, 3}, // + {-1, -1, -1}, // + {1, 4, 3}, // + {-1, -1, -1}, // + {-1, -1, -1}, // + {-1, -1, -1}}; + return tetTri1[i]; +} + +__host__ __device__ glm::ivec4 Neighbors(int i) { + constexpr glm::ivec4 neighbors[14] = {{0, 0, 0, 1}, // + {1, 0, 0, 0}, // + {0, 1, 0, 0}, // + {0, 0, 1, 0}, // + {-1, 0, 0, 1}, // + {0, -1, 0, 1}, // + {0, 0, -1, 1}, // + {-1, -1, -1, 1}, // + {-1, 0, 0, 0}, // + {0, -1, 0, 0}, // + {0, 0, -1, 0}, // + {0, -1, -1, 1}, // + {-1, 0, -1, 1}, // + {-1, -1, 0, 1}}; + return neighbors[i]; +} + +__host__ __device__ Uint64 SpreadBits3(Uint64 v) { v = v & 0x1fffff; v = (v | v << 32) & 0x1f00000000ffff; v = (v | v << 16) & 0x1f0000ff0000ff; @@ -82,7 +115,7 @@ __host__ __device__ uint64_t SpreadBits3(uint64_t v) { return v; } -__host__ __device__ uint64_t SqeezeBits3(uint64_t v) { +__host__ __device__ Uint64 SqeezeBits3(Uint64 v) { v = v & 0x1249249249249249; v = (v ^ v >> 2) & 0x10c30c30c30c30c3; v = (v ^ v >> 4) & 0x100f00f00f00f00f; @@ -95,12 +128,12 @@ __host__ __device__ uint64_t SqeezeBits3(uint64_t v) { // This is a modified 3D MortonCode, where the xyz code is shifted by one bit // and the w bit is added as the least significant. This allows 21 bits per x, // y, and z channel and 1 for w, filling the 64 bit total. -__device__ __host__ uint64_t MortonCode(const glm::ivec4& index) { - return static_cast(index.w) | (SpreadBits3(index.x) << 1) | +__host__ __device__ Uint64 MortonCode(const glm::ivec4& index) { + return static_cast(index.w) | (SpreadBits3(index.x) << 1) | (SpreadBits3(index.y) << 2) | (SpreadBits3(index.z) << 3); } -__device__ __host__ glm::ivec4 DecodeMorton(uint64_t code) { +__host__ __device__ glm::ivec4 DecodeMorton(Uint64 code) { glm::ivec4 index; index.x = SqeezeBits3(code >> 1); index.y = SqeezeBits3(code >> 2); @@ -110,12 +143,12 @@ __device__ __host__ glm::ivec4 DecodeMorton(uint64_t code) { } struct GridVert { - uint64_t key = kOpen; + Uint64 key = kOpen; float distance = NAN; int edgeIndex = -1; int edgeVerts[7] = {-1, -1, -1, -1, -1, -1, -1}; - int Inside() const { return distance > 0 ? 1 : -1; } + __host__ __device__ int Inside() const { return distance > 0 ? 1 : -1; } }; class HashTableD { @@ -123,12 +156,13 @@ class HashTableD { HashTableD(VecDH& alloc, VecDH& used, uint32_t step = 127) : step_{step}, table_{alloc}, used_{used} {} - __device__ __host__ int Size() const { return table_.size(); } + __host__ __device__ int Size() const { return table_.size(); } - __device__ __host__ bool Insert(const GridVert& vert) { + __host__ __device__ bool Insert(const GridVert& vert) { uint32_t idx = vert.key & (Size() - 1); while (1) { - const uint64_t found = AtomicCAS(table_[idx].key, kOpen, vert.key); + Uint64& key = table_[idx].key; + const Uint64 found = AtomicCAS(key, kOpen, vert.key); if (found == kOpen) { if (AtomicAdd(used_[0], 0x1u) * 2 > Size()) { return true; @@ -141,7 +175,7 @@ class HashTableD { } } - __device__ __host__ GridVert operator[](uint64_t key) const { + __host__ __device__ GridVert operator[](Uint64 key) const { uint32_t idx = key & (Size() - 1); while (1) { const GridVert found = table_[idx]; @@ -150,7 +184,7 @@ class HashTableD { } } - __device__ __host__ GridVert At(int idx) const { return table_[idx]; } + __host__ __device__ GridVert At(int idx) const { return table_[idx]; } private: const uint32_t step_; @@ -208,7 +242,7 @@ struct ComputeVerts { return d; } - inline __host__ __device__ void operator()(uint64_t mortonCode) { + inline __host__ __device__ void operator()(Uint64 mortonCode) { const glm::ivec4 gridIndex = DecodeMorton(mortonCode); if (glm::any(glm::greaterThan(glm::ivec3(gridIndex), gridSize))) return; @@ -225,8 +259,7 @@ struct ComputeVerts { bool keep = false; float minDist2 = 0.25 * 0.25; for (int i = 0; i < 14; ++i) { - const int j = i < 7 ? i : i - 7; - glm::ivec4 neighborIndex = gridIndex + neighbors[i]; + glm::ivec4 neighborIndex = gridIndex + Neighbors(i); if (neighborIndex.w == 2) { neighborIndex += 1; neighborIndex.w = 0; @@ -258,8 +291,7 @@ struct ComputeVerts { } } - if (keep && gridVerts.Insert(gridVert)) - std::cout << "out of space!" << std::endl; + if (keep && gridVerts.Insert(gridVert)) printf("out of space!\n"); } }; @@ -279,8 +311,8 @@ struct BuildTris { const int edges[6]) { const int i = (tet[0] > 0 ? 1 : 0) + (tet[1] > 0 ? 2 : 0) + (tet[2] > 0 ? 4 : 0) + (tet[3] > 0 ? 8 : 0); - CreateTri(tetTri0[i], edges); - CreateTri(tetTri1[i], edges); + CreateTri(TetTri0(i), edges); + CreateTri(TetTri1(i), edges); } __host__ __device__ void operator()(int idx) { @@ -314,15 +346,15 @@ struct BuildTris { for (const int i : {0, 1, 2}) { edges[1] = base.edgeVerts[i + 1]; edges[4] = thisVert.edgeVerts[i + 4]; - edges[5] = base.edgeVerts[prev3[i] + 4]; + edges[5] = base.edgeVerts[Prev3(i) + 4]; thisIndex = leadIndex; - thisIndex[prev3[i]] -= 1; + thisIndex[Prev3(i)] -= 1; thisVert = gridVerts[MortonCode(thisIndex)]; tet[3] = thisVert.Inside(); - edges[2] = thisVert.edgeVerts[next3[i] + 4]; - edges[3] = thisVert.edgeVerts[prev3[i] + 1]; + edges[2] = thisVert.edgeVerts[Next3(i) + 4]; + edges[3] = thisVert.edgeVerts[Prev3(i) + 1]; if (!skipTet && thisVert.key != kOpen) CreateTris(tet, edges); skipTet = thisVert.key == kOpen; @@ -331,14 +363,14 @@ struct BuildTris { edges[1] = edges[5]; edges[2] = thisVert.edgeVerts[i + 4]; edges[4] = edges[3]; - edges[5] = base.edgeVerts[next3[i] + 1]; + edges[5] = base.edgeVerts[Next3(i) + 1]; thisIndex = baseIndex; - thisIndex[next3[i]] += 1; + thisIndex[Next3(i)] += 1; thisVert = gridVerts[MortonCode(thisIndex)]; tet[3] = thisVert.Inside(); - edges[3] = thisVert.edgeVerts[next3[i] + 4]; + edges[3] = thisVert.edgeVerts[Next3(i) + 4]; if (!skipTet && thisVert.key != kOpen) CreateTris(tet, edges); skipTet = thisVert.key == kOpen; diff --git a/src/utilities/include/structs.h b/src/utilities/include/structs.h index 6f1b776a4..c6b806a06 100644 --- a/src/utilities/include/structs.h +++ b/src/utilities/include/structs.h @@ -32,8 +32,6 @@ namespace manifold { constexpr float kTolerance = 1e-5; -constexpr glm::ivec3 next3(1, 2, 0); -constexpr glm::ivec3 prev3(2, 0, 1); #ifdef __CUDACC__ #define HOST_DEVICE __host__ __device__ diff --git a/src/utilities/include/utils.h b/src/utilities/include/utils.h index 1157b1151..2ac530310 100644 --- a/src/utilities/include/utils.h +++ b/src/utilities/include/utils.h @@ -114,18 +114,6 @@ thrust::counting_iterator countAt(T i) { return thrust::make_counting_iterator(i); } -template -__host__ __device__ T AtomicCAS(T& target, T compare, T val) { -#ifdef __CUDA_ARCH__ - return atomicCAS(&target, compare, val); -#else - std::atomic& tar = reinterpret_cast&>(target); - T old_val = tar.load(); - tar.compare_exchange_weak(compare, val); - return old_val; -#endif -} - template __host__ __device__ T AtomicAdd(T& target, T add) { #ifdef __CUDA_ARCH__ From ba2348f38ed752d49eebdf14a73b5e1fd7a4827b Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Tue, 16 Aug 2022 09:21:22 -0700 Subject: [PATCH 18/39] moved SDF --- extras/CMakeLists.txt | 2 +- src/CMakeLists.txt | 1 - src/{sdf => manifold}/include/sdf.h | 2 +- src/sdf/CMakeLists.txt | 10 ---------- 4 files changed, 2 insertions(+), 13 deletions(-) rename src/{sdf => manifold}/include/sdf.h (99%) delete mode 100644 src/sdf/CMakeLists.txt diff --git a/extras/CMakeLists.txt b/extras/CMakeLists.txt index 8ee47740b..7affde9c3 100644 --- a/extras/CMakeLists.txt +++ b/extras/CMakeLists.txt @@ -33,7 +33,7 @@ if(BUILD_TEST_CGAL) endif() add_executable(sdfTest sdfTest.cpp) -target_link_libraries(sdfTest manifold meshIO samples sdf) +target_link_libraries(sdfTest manifold meshIO samples) if(MANIFOLD_USE_CUDA) set_source_files_properties(sdfTest.cpp PROPERTIES LANGUAGE CUDA) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 989db9710..0a502d5e4 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -17,4 +17,3 @@ add_subdirectory(utilities) add_subdirectory(collider) add_subdirectory(polygon) add_subdirectory(manifold) -add_subdirectory(sdf) diff --git a/src/sdf/include/sdf.h b/src/manifold/include/sdf.h similarity index 99% rename from src/sdf/include/sdf.h rename to src/manifold/include/sdf.h index 50cf3322c..a2d8698af 100644 --- a/src/sdf/include/sdf.h +++ b/src/manifold/include/sdf.h @@ -14,7 +14,7 @@ #pragma once -#include "structs.h" +#include "manifold.h" #include "utils.h" #include "vec_dh.h" diff --git a/src/sdf/CMakeLists.txt b/src/sdf/CMakeLists.txt deleted file mode 100644 index 29cb6dcc0..000000000 --- a/src/sdf/CMakeLists.txt +++ /dev/null @@ -1,10 +0,0 @@ -project (sdf) - -add_library(${PROJECT_NAME} INTERFACE) - -target_include_directories(${PROJECT_NAME} - INTERFACE - ${PROJECT_SOURCE_DIR}/include -) - -target_link_libraries( ${PROJECT_NAME} INTERFACE utilities) From 059f86d6fbe0f817004297d74e9f37b8a5b37a9a Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Tue, 16 Aug 2022 09:37:28 -0700 Subject: [PATCH 19/39] Revert "moved SDF" This reverts commit ba2348f38ed752d49eebdf14a73b5e1fd7a4827b. --- extras/CMakeLists.txt | 2 +- src/CMakeLists.txt | 1 + src/sdf/CMakeLists.txt | 10 ++++++++++ src/{manifold => sdf}/include/sdf.h | 2 +- 4 files changed, 13 insertions(+), 2 deletions(-) create mode 100644 src/sdf/CMakeLists.txt rename src/{manifold => sdf}/include/sdf.h (99%) diff --git a/extras/CMakeLists.txt b/extras/CMakeLists.txt index 7affde9c3..8ee47740b 100644 --- a/extras/CMakeLists.txt +++ b/extras/CMakeLists.txt @@ -33,7 +33,7 @@ if(BUILD_TEST_CGAL) endif() add_executable(sdfTest sdfTest.cpp) -target_link_libraries(sdfTest manifold meshIO samples) +target_link_libraries(sdfTest manifold meshIO samples sdf) if(MANIFOLD_USE_CUDA) set_source_files_properties(sdfTest.cpp PROPERTIES LANGUAGE CUDA) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 0a502d5e4..989db9710 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -17,3 +17,4 @@ add_subdirectory(utilities) add_subdirectory(collider) add_subdirectory(polygon) add_subdirectory(manifold) +add_subdirectory(sdf) diff --git a/src/sdf/CMakeLists.txt b/src/sdf/CMakeLists.txt new file mode 100644 index 000000000..29cb6dcc0 --- /dev/null +++ b/src/sdf/CMakeLists.txt @@ -0,0 +1,10 @@ +project (sdf) + +add_library(${PROJECT_NAME} INTERFACE) + +target_include_directories(${PROJECT_NAME} + INTERFACE + ${PROJECT_SOURCE_DIR}/include +) + +target_link_libraries( ${PROJECT_NAME} INTERFACE utilities) diff --git a/src/manifold/include/sdf.h b/src/sdf/include/sdf.h similarity index 99% rename from src/manifold/include/sdf.h rename to src/sdf/include/sdf.h index a2d8698af..50cf3322c 100644 --- a/src/manifold/include/sdf.h +++ b/src/sdf/include/sdf.h @@ -14,7 +14,7 @@ #pragma once -#include "manifold.h" +#include "structs.h" #include "utils.h" #include "vec_dh.h" From 8b1a17e759865bdd7b6d6b20622d856c8fc32bc3 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Tue, 16 Aug 2022 10:43:03 -0700 Subject: [PATCH 20/39] added gyroid sample --- extras/CMakeLists.txt | 2 +- samples/CMakeLists.txt | 38 +++++++++++++++++++++++++------- samples/include/samples.h | 2 ++ samples/src/gyroid_module.cpp | 41 +++++++++++++++++++++++++++++++++++ src/sdf/CMakeLists.txt | 20 ++++++++++++++--- src/sdf/include/sdf.h | 2 +- test/CMakeLists.txt | 4 ++-- test/samples_test.cpp | 8 +++++++ 8 files changed, 102 insertions(+), 15 deletions(-) create mode 100644 samples/src/gyroid_module.cpp diff --git a/extras/CMakeLists.txt b/extras/CMakeLists.txt index 8ee47740b..8ce895d92 100644 --- a/extras/CMakeLists.txt +++ b/extras/CMakeLists.txt @@ -33,7 +33,7 @@ if(BUILD_TEST_CGAL) endif() add_executable(sdfTest sdfTest.cpp) -target_link_libraries(sdfTest manifold meshIO samples sdf) +target_link_libraries(sdfTest manifold meshIO sdf) if(MANIFOLD_USE_CUDA) set_source_files_properties(sdfTest.cpp PROPERTIES LANGUAGE CUDA) diff --git a/samples/CMakeLists.txt b/samples/CMakeLists.txt index bee760c8a..fe16a0761 100644 --- a/samples/CMakeLists.txt +++ b/samples/CMakeLists.txt @@ -4,7 +4,7 @@ # you may not use this file except in compliance with the License. # You may obtain a copy of the License at # -# https://www.apache.org/licenses/LICENSE-2.0 +# https://www.apache.org/licenses/LICENSE-2.0 # # Unless required by applicable law or agreed to in writing, software # distributed under the License is distributed on an "AS IS" BASIS, @@ -12,18 +12,40 @@ # See the License for the specific language governing permissions and # limitations under the License. -project (samples) +project(samples) -file(GLOB_RECURSE SOURCE_FILES CONFIGURE_DEPENDS *.cpp) -add_library(${PROJECT_NAME} ${SOURCE_FILES}) +add_library(samples src/bracelet.cpp src/knot.cpp src/menger_sponge.cpp src/rounded_frame.cpp src/scallop.cpp src/tet_puzzle.cpp) -target_include_directories( ${PROJECT_NAME} +target_include_directories(samples PUBLIC ${PROJECT_SOURCE_DIR}/include ) -target_link_libraries( ${PROJECT_NAME} +target_link_libraries(samples PUBLIC manifold ) -target_compile_options(${PROJECT_NAME} PRIVATE ${MANIFOLD_FLAGS}) -target_compile_features(${PROJECT_NAME} PUBLIC cxx_std_14) +target_compile_options(samples PRIVATE ${MANIFOLD_FLAGS}) +target_compile_features(samples PUBLIC cxx_std_14) + +project(samplesGPU) + +add_library(samplesGPU src/gyroid_module.cpp) + +if(MANIFOLD_USE_CUDA) + set_source_files_properties(src/gyroid_module.cpp PROPERTIES LANGUAGE CUDA) + set_property(TARGET samplesGPU PROPERTY CUDA_ARCHITECTURES 61) +endif() + +target_include_directories(samplesGPU + PUBLIC ${PROJECT_SOURCE_DIR}/include +) + +target_link_libraries(samplesGPU + PUBLIC manifold sdf +) + +target_compile_options(samplesGPU PRIVATE ${MANIFOLD_FLAGS}) +target_compile_features(samplesGPU + PUBLIC cxx_std_14 + PRIVATE cxx_std_17 +) diff --git a/samples/include/samples.h b/samples/include/samples.h index 202ae4feb..291dea844 100644 --- a/samples/include/samples.h +++ b/samples/include/samples.h @@ -49,5 +49,7 @@ Manifold RoundedFrame(float edgeLength, float radius, int circularSegments = 0); Manifold TetPuzzle(float edgeLength, float gap, int nDivisions); Manifold Scallop(); + +Manifold GyroidModule(int n); /** @} */ // end of Samples } // namespace manifold \ No newline at end of file diff --git a/samples/src/gyroid_module.cpp b/samples/src/gyroid_module.cpp new file mode 100644 index 000000000..4fe4a7d95 --- /dev/null +++ b/samples/src/gyroid_module.cpp @@ -0,0 +1,41 @@ +// Copyright 2022 The Manifold Authors. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "samples.h" +#include "sdf.h" + +namespace { + +struct Gyroid { + __host__ __device__ float operator()(glm::vec3 p) const { + return cos(p.x) * sin(p.y) + cos(p.y) * sin(p.z) + cos(p.z) * sin(p.x); + } +}; +} // namespace + +namespace manifold { + +/** + * Creates a + */ +Manifold GyroidModule(int n) { + Gyroid func; + const SDF gyroidSDF(func); + + const float size = glm::two_pi(); + Manifold gyroid( + gyroidSDF.LevelSet({glm::vec3(-size), glm::vec3(size)}, size / n)); + return gyroid; +} +} \ No newline at end of file diff --git a/src/sdf/CMakeLists.txt b/src/sdf/CMakeLists.txt index 29cb6dcc0..4e0ad3a99 100644 --- a/src/sdf/CMakeLists.txt +++ b/src/sdf/CMakeLists.txt @@ -1,10 +1,24 @@ -project (sdf) +# Copyright 2022 The Manifold Authors. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +project(sdf) add_library(${PROJECT_NAME} INTERFACE) target_include_directories(${PROJECT_NAME} INTERFACE - ${PROJECT_SOURCE_DIR}/include + ${PROJECT_SOURCE_DIR}/include ) -target_link_libraries( ${PROJECT_NAME} INTERFACE utilities) +target_link_libraries(${PROJECT_NAME} INTERFACE utilities) diff --git a/src/sdf/include/sdf.h b/src/sdf/include/sdf.h index 50cf3322c..19ea4956d 100644 --- a/src/sdf/include/sdf.h +++ b/src/sdf/include/sdf.h @@ -1,4 +1,4 @@ -// Copyright 2022 Emmett Lalish +// Copyright 2022 The Manifold Authors. // // Licensed under the Apache License, Version 2.0 (the "License"); // you may not use this file except in compliance with the License. diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index cfc7fa114..eb58cb247 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -4,7 +4,7 @@ # you may not use this file except in compliance with the License. # You may obtain a copy of the License at # -# https://www.apache.org/licenses/LICENSE-2.0 +# https://www.apache.org/licenses/LICENSE-2.0 # # Unless required by applicable law or agreed to in writing, software # distributed under the License is distributed on an "AS IS" BASIS, @@ -21,7 +21,7 @@ enable_testing() set(SOURCE_FILES polygon_test.cpp mesh_test.cpp samples_test.cpp test_main.cpp) add_executable(${PROJECT_NAME} ${SOURCE_FILES}) -target_link_libraries(${PROJECT_NAME} polygon GTest::GTest manifold meshIO samples) +target_link_libraries(${PROJECT_NAME} polygon GTest::GTest manifold meshIO samples samplesGPU) target_compile_options(${PROJECT_NAME} PRIVATE ${MANIFOLD_FLAGS}) target_compile_features(${PROJECT_NAME} PUBLIC cxx_std_14) diff --git a/test/samples_test.cpp b/test/samples_test.cpp index f9ec9ec90..21953d7c7 100644 --- a/test/samples_test.cpp +++ b/test/samples_test.cpp @@ -162,6 +162,14 @@ TEST(Samples, Bracelet) { if (options.exportModels) ExportMesh("bracelet.glb", bracelet.GetMesh(), {}); } +TEST(Samples, GyroidModule) { + Manifold gyroid = GyroidModule(50); + CheckManifold(gyroid); + EXPECT_LE(gyroid.NumDegenerateTris(), 0); + EXPECT_EQ(gyroid.Genus(), 1); + if (options.exportModels) ExportMesh("gyroid.gltf", gyroid.GetMesh(), {}); +} + TEST(Samples, Sponge1) { Manifold sponge = MengerSponge(1); CheckManifold(sponge); From 73e395d7584003bbef1796d65e36f0c0938aa246 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Tue, 16 Aug 2022 11:42:26 -0700 Subject: [PATCH 21/39] added sdf_test --- samples/src/gyroid_module.cpp | 5 ++-- test/CMakeLists.txt | 18 +++++------- test/samples_test.cpp | 2 +- test/sdf_test.cpp | 53 +++++++++++++++++++++++++++++++++++ 4 files changed, 64 insertions(+), 14 deletions(-) create mode 100644 test/sdf_test.cpp diff --git a/samples/src/gyroid_module.cpp b/samples/src/gyroid_module.cpp index 4fe4a7d95..5e083581e 100644 --- a/samples/src/gyroid_module.cpp +++ b/samples/src/gyroid_module.cpp @@ -19,6 +19,7 @@ namespace { struct Gyroid { __host__ __device__ float operator()(glm::vec3 p) const { + p -= glm::pi() / 4; return cos(p.x) * sin(p.y) + cos(p.y) * sin(p.z) + cos(p.z) * sin(p.x); } }; @@ -35,7 +36,7 @@ Manifold GyroidModule(int n) { const float size = glm::two_pi(); Manifold gyroid( - gyroidSDF.LevelSet({glm::vec3(-size), glm::vec3(size)}, size / n)); + gyroidSDF.LevelSet({glm::vec3(-size), glm::vec3(size)}, size / n, 0.5)); return gyroid; } -} \ No newline at end of file +} // namespace manifold \ No newline at end of file diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index eb58cb247..8e4bd53ea 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -19,9 +19,14 @@ add_subdirectory(meshIO) enable_testing() -set(SOURCE_FILES polygon_test.cpp mesh_test.cpp samples_test.cpp test_main.cpp) +set(SOURCE_FILES polygon_test.cpp mesh_test.cpp samples_test.cpp sdf_test.cpp test_main.cpp) add_executable(${PROJECT_NAME} ${SOURCE_FILES}) -target_link_libraries(${PROJECT_NAME} polygon GTest::GTest manifold meshIO samples samplesGPU) +target_link_libraries(${PROJECT_NAME} polygon GTest::GTest manifold meshIO samples samplesGPU sdf) + +if(MANIFOLD_USE_CUDA) + set_source_files_properties(sdf_test.cpp PROPERTIES LANGUAGE CUDA) + set_property(TARGET ${PROJECT_NAME} PROPERTY CUDA_ARCHITECTURES 61) +endif() target_compile_options(${PROJECT_NAME} PRIVATE ${MANIFOLD_FLAGS}) target_compile_features(${PROJECT_NAME} PUBLIC cxx_std_14) @@ -35,12 +40,3 @@ if(EMSCRIPTEN) "-s ASSERTIONS=1 -s DEMANGLE_SUPPORT=1 --preload-file data --bind") else() endif() - -if(APPLE) - # All executables that link to CUDA need this. If you ever get - # "CUDA driver version is insufficient for CUDA runtime version", - # this is probably what's missing. - set_property(TARGET ${PROJECT_NAME} - PROPERTY - BUILD_RPATH ${CMAKE_CUDA_IMPLICIT_LINK_DIRECTORIES}) -endif() diff --git a/test/samples_test.cpp b/test/samples_test.cpp index 21953d7c7..9a0ba8cfa 100644 --- a/test/samples_test.cpp +++ b/test/samples_test.cpp @@ -163,7 +163,7 @@ TEST(Samples, Bracelet) { } TEST(Samples, GyroidModule) { - Manifold gyroid = GyroidModule(50); + Manifold gyroid = GyroidModule(20); CheckManifold(gyroid); EXPECT_LE(gyroid.NumDegenerateTris(), 0); EXPECT_EQ(gyroid.Genus(), 1); diff --git a/test/sdf_test.cpp b/test/sdf_test.cpp new file mode 100644 index 000000000..04ef3887d --- /dev/null +++ b/test/sdf_test.cpp @@ -0,0 +1,53 @@ +// Copyright 2022 The Manifold Authors. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "sdf.h" + +#include "manifold.h" +#include "meshIO.h" +#include "test.h" + +using namespace manifold; + +struct CubeVoid { + __host__ __device__ float operator()(glm::vec3 p) const { + const glm::vec3 min = glm::vec3(1) - p; + const glm::vec3 max = p - glm::vec3(1); + const float min3 = glm::max(min.x, glm::max(min.y, min.z)); + const float max3 = glm::max(max.x, glm::max(max.y, max.z)); + return glm::max(min3, max3); + } +}; + +TEST(SDF, Position) { + CubeVoid func; + const SDF voidSDF(func); + + const float size = 4; + const float edgeLength = 0.5; + + Manifold cubeVoid(voidSDF.LevelSet( + {glm::vec3(-size / 2), glm::vec3(size / 2)}, edgeLength)); + Box bounds = cubeVoid.BoundingBox(); + if (options.exportModels) ExportMesh("cubeVoid.gltf", cubeVoid.GetMesh(), {}); + + EXPECT_TRUE(cubeVoid.IsManifold()); + EXPECT_EQ(cubeVoid.Genus(), -1); + EXPECT_NEAR(bounds.min.x, -1 - edgeLength / 2, 0.001); + EXPECT_NEAR(bounds.min.y, -1 - edgeLength / 2, 0.001); + EXPECT_NEAR(bounds.min.z, -1 - edgeLength / 2, 0.001); + EXPECT_NEAR(bounds.max.x, 1 + edgeLength / 2, 0.001); + EXPECT_NEAR(bounds.max.y, 1 + edgeLength / 2, 0.001); + EXPECT_NEAR(bounds.max.z, 1 + edgeLength / 2, 0.001); +} \ No newline at end of file From 0fc8c8f012cf558577cdea221bcad919f06ad45b Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Thu, 18 Aug 2022 10:54:19 -0700 Subject: [PATCH 22/39] added SDF tests --- extras/CMakeLists.txt | 10 ---------- extras/sdfTest.cpp | 20 -------------------- src/sdf/CMakeLists.txt | 19 +++++++++++++------ src/sdf/include/sdf.h | 10 ++++++---- test/sdf_test.cpp | 34 ++++++++++++++++++++++++---------- 5 files changed, 43 insertions(+), 50 deletions(-) delete mode 100644 extras/sdfTest.cpp diff --git a/extras/CMakeLists.txt b/extras/CMakeLists.txt index 8ce895d92..e9c1bd317 100644 --- a/extras/CMakeLists.txt +++ b/extras/CMakeLists.txt @@ -31,13 +31,3 @@ if(BUILD_TEST_CGAL) target_compile_options(perfTestCGAL PRIVATE ${MANIFOLD_FLAGS}) target_compile_features(perfTestCGAL PUBLIC cxx_std_14) endif() - -add_executable(sdfTest sdfTest.cpp) -target_link_libraries(sdfTest manifold meshIO sdf) - -if(MANIFOLD_USE_CUDA) - set_source_files_properties(sdfTest.cpp PROPERTIES LANGUAGE CUDA) - set_property(TARGET sdfTest PROPERTY CUDA_ARCHITECTURES 61) -endif() - -target_compile_features(sdfTest PUBLIC cxx_std_14) diff --git a/extras/sdfTest.cpp b/extras/sdfTest.cpp deleted file mode 100644 index 7774a8f46..000000000 --- a/extras/sdfTest.cpp +++ /dev/null @@ -1,20 +0,0 @@ -#include "meshIO.h" -#include "sdf.h" - -using namespace manifold; - -struct Test { - __host__ __device__ float operator()(glm::vec3 point) const { - return glm::dot(point, {1, 1, 2}); - } -}; - -int main(int argc, char **argv) { - Test func; - SDF a(func); - Box box({-1, -1, -1}, {1, 2, 3}); - Mesh b = a.LevelSet(box, 0.5, -0.1); - // Dump(b.vertPos); - // Dump(b.triVerts); - ExportMesh("sdf.gltf", b, {}); -} \ No newline at end of file diff --git a/src/sdf/CMakeLists.txt b/src/sdf/CMakeLists.txt index 4e0ad3a99..2038246bd 100644 --- a/src/sdf/CMakeLists.txt +++ b/src/sdf/CMakeLists.txt @@ -14,11 +14,18 @@ project(sdf) -add_library(${PROJECT_NAME} INTERFACE) +file(GLOB_RECURSE SOURCE_FILES CONFIGURE_DEPENDS *.cpp) +add_library(${PROJECT_NAME} ${SOURCE_FILES}) -target_include_directories(${PROJECT_NAME} - INTERFACE - ${PROJECT_SOURCE_DIR}/include -) +if(MANIFOLD_USE_CUDA) + set_source_files_properties(${SOURCE_FILES} PROPERTIES LANGUAGE CUDA) + set_property(TARGET ${PROJECT_NAME} PROPERTY CUDA_ARCHITECTURES 61) +endif() + +target_include_directories(${PROJECT_NAME} PUBLIC ${PROJECT_SOURCE_DIR}/include) +target_link_libraries(${PROJECT_NAME} PUBLIC utilities) -target_link_libraries(${PROJECT_NAME} INTERFACE utilities) +target_compile_features(${PROJECT_NAME} + PUBLIC cxx_std_14 + PRIVATE cxx_std_17 +) diff --git a/src/sdf/include/sdf.h b/src/sdf/include/sdf.h index 19ea4956d..78bf138d5 100644 --- a/src/sdf/include/sdf.h +++ b/src/sdf/include/sdf.h @@ -247,9 +247,6 @@ struct ComputeVerts { if (glm::any(glm::greaterThan(glm::ivec3(gridIndex), gridSize))) return; - // const auto sdfFunc = - // AtBounds(gridIndex) ? &ComputeVerts::BoundedSdf : &ComputeVerts::Sdf; - const glm::vec3 position = Position(gridIndex); GridVert gridVert; @@ -382,6 +379,10 @@ struct BuildTris { } // namespace namespace manifold { + +void RemoveUnreferencedVerts(VecDH& vertPos, + VecDH& triVerts); + /** @addtogroup Core * @{ */ @@ -427,6 +428,8 @@ class SDF { BuildTris({triVerts.ptrD(), index.ptrD(), gridVerts.D()})); triVerts.resize(index[0]); + RemoveUnreferencedVerts(vertPos, triVerts); + out.vertPos.insert(out.vertPos.end(), vertPos.begin(), vertPos.end()); out.triVerts.insert(out.triVerts.end(), triVerts.begin(), triVerts.end()); return out; @@ -435,6 +438,5 @@ class SDF { private: const Func sdf_; }; - /** @} */ } // namespace manifold \ No newline at end of file diff --git a/test/sdf_test.cpp b/test/sdf_test.cpp index 04ef3887d..6853cc21d 100644 --- a/test/sdf_test.cpp +++ b/test/sdf_test.cpp @@ -23,13 +23,27 @@ using namespace manifold; struct CubeVoid { __host__ __device__ float operator()(glm::vec3 p) const { const glm::vec3 min = glm::vec3(1) - p; - const glm::vec3 max = p - glm::vec3(1); - const float min3 = glm::max(min.x, glm::max(min.y, min.z)); - const float max3 = glm::max(max.x, glm::max(max.y, max.z)); - return glm::max(min3, max3); + const glm::vec3 max = p + glm::vec3(1); + const float min3 = glm::min(min.x, glm::min(min.y, min.z)); + const float max3 = glm::min(max.x, glm::min(max.y, max.z)); + return -1.0f * glm::min(min3, max3); } }; +TEST(SDF, CubeVoid) { + CubeVoid func; + const SDF voidSDF(func); + + EXPECT_EQ(voidSDF({0, 0, 0}), -1); + EXPECT_EQ(voidSDF({0, 0, 1}), 0); + EXPECT_EQ(voidSDF({0, 1, 1}), 0); + EXPECT_EQ(voidSDF({-1, 0, 0}), 0); + EXPECT_EQ(voidSDF({1, 1, -1}), 0); + EXPECT_EQ(voidSDF({2, 0, 0}), 1); + EXPECT_EQ(voidSDF({2, -2, 0}), 1); + EXPECT_EQ(voidSDF({-2, 2, 2}), 1); +} + TEST(SDF, Position) { CubeVoid func; const SDF voidSDF(func); @@ -44,10 +58,10 @@ TEST(SDF, Position) { EXPECT_TRUE(cubeVoid.IsManifold()); EXPECT_EQ(cubeVoid.Genus(), -1); - EXPECT_NEAR(bounds.min.x, -1 - edgeLength / 2, 0.001); - EXPECT_NEAR(bounds.min.y, -1 - edgeLength / 2, 0.001); - EXPECT_NEAR(bounds.min.z, -1 - edgeLength / 2, 0.001); - EXPECT_NEAR(bounds.max.x, 1 + edgeLength / 2, 0.001); - EXPECT_NEAR(bounds.max.y, 1 + edgeLength / 2, 0.001); - EXPECT_NEAR(bounds.max.z, 1 + edgeLength / 2, 0.001); + EXPECT_NEAR(bounds.min.x, -size / 2 - edgeLength / 2, 0.001); + EXPECT_NEAR(bounds.min.y, -size / 2 - edgeLength / 2, 0.001); + EXPECT_NEAR(bounds.min.z, -size / 2 - edgeLength / 2, 0.001); + EXPECT_NEAR(bounds.max.x, size / 2 + edgeLength / 2, 0.001); + EXPECT_NEAR(bounds.max.y, size / 2 + edgeLength / 2, 0.001); + EXPECT_NEAR(bounds.max.z, size / 2 + edgeLength / 2, 0.001); } \ No newline at end of file From 01bdb70c116e393d57c94fe09a0457098f73c220 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Thu, 18 Aug 2022 11:20:29 -0700 Subject: [PATCH 23/39] tests pass --- src/sdf/include/sdf.h | 6 ++-- src/sdf/src/sdf.cpp | 68 +++++++++++++++++++++++++++++++++++++++++++ test/sdf_test.cpp | 45 +++++++++++++++++++++++----- 3 files changed, 108 insertions(+), 11 deletions(-) create mode 100644 src/sdf/src/sdf.cpp diff --git a/src/sdf/include/sdf.h b/src/sdf/include/sdf.h index 78bf138d5..5e314db5b 100644 --- a/src/sdf/include/sdf.h +++ b/src/sdf/include/sdf.h @@ -403,12 +403,10 @@ class SDF { const glm::ivec3 gridSize(dim / edgeLength); const glm::vec3 spacing = dim / (glm::vec3(gridSize)); - const int maxMorton = MortonCode(glm::ivec4(gridSize - 1, 1)); + const int maxMorton = MortonCode(glm::ivec4(gridSize + 1, 1)); const int tableSize = glm::min( 2 * maxMorton, static_cast(100 * glm::pow(maxMorton, 0.667))); - std::cout << "maxMorton: " << maxMorton - << ", hash table size: " << tableSize << std::endl; HashTable gridVerts(tableSize); VecDH vertPos(gridVerts.Size() * 7); @@ -417,7 +415,7 @@ class SDF { thrust::for_each_n( countAt(0), maxMorton + 1, ComputeVerts({vertPos.ptrD(), index.ptrD(), gridVerts.D(), sdf_, - bounds.min, gridSize - 1, spacing, level})); + bounds.min, gridSize + 1, spacing, level})); vertPos.resize(index[0]); VecDH triVerts(gridVerts.Entries() * 12); // worst case diff --git a/src/sdf/src/sdf.cpp b/src/sdf/src/sdf.cpp new file mode 100644 index 000000000..c3b7cfbcc --- /dev/null +++ b/src/sdf/src/sdf.cpp @@ -0,0 +1,68 @@ +// Copyright 2022 The Manifold Authors. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "sdf.h" + +namespace { +using namespace manifold; + +struct MarkVerts { + int* vert; + + __host__ __device__ void operator()(glm::ivec3 triVerts) { + for (int i : {0, 1, 2}) { + vert[triVerts[i]] = 1; + } + } +}; + +struct ReindexVerts { + const int* old2new; + + __host__ __device__ void operator()(glm::ivec3& triVerts) { + for (int i : {0, 1, 2}) { + triVerts[i] = old2new[triVerts[i]]; + } + } +}; + +} // namespace + +namespace manifold { +/** + * Sorts the vertices according to their Morton code. + */ +void RemoveUnreferencedVerts(VecDH& vertPos, + VecDH& triVerts) { + const int numVert = vertPos.size(); + VecDH vertOld2New(numVert + 1, 0); + auto policy = autoPolicy(numVert); + for_each(policy, triVerts.cbegin(), triVerts.cend(), + MarkVerts({vertOld2New.ptrD() + 1})); + + const VecDH oldVertPos = vertPos; + vertPos.resize(copy_if( + policy, oldVertPos.cbegin(), oldVertPos.cend(), + vertOld2New.cbegin() + 1, vertPos.begin(), + thrust::identity()) - + vertPos.begin()); + + inclusive_scan(policy, vertOld2New.begin() + 1, vertOld2New.end(), + vertOld2New.begin() + 1); + + for_each(policy, triVerts.begin(), triVerts.end(), + ReindexVerts({vertOld2New.cptrD()})); +} + +} // namespace manifold \ No newline at end of file diff --git a/test/sdf_test.cpp b/test/sdf_test.cpp index 6853cc21d..b361d33e6 100644 --- a/test/sdf_test.cpp +++ b/test/sdf_test.cpp @@ -44,7 +44,7 @@ TEST(SDF, CubeVoid) { EXPECT_EQ(voidSDF({-2, 2, 2}), 1); } -TEST(SDF, Position) { +TEST(SDF, Bounds) { CubeVoid func; const SDF voidSDF(func); @@ -54,14 +54,45 @@ TEST(SDF, Position) { Manifold cubeVoid(voidSDF.LevelSet( {glm::vec3(-size / 2), glm::vec3(size / 2)}, edgeLength)); Box bounds = cubeVoid.BoundingBox(); + const float precision = cubeVoid.Precision(); if (options.exportModels) ExportMesh("cubeVoid.gltf", cubeVoid.GetMesh(), {}); EXPECT_TRUE(cubeVoid.IsManifold()); EXPECT_EQ(cubeVoid.Genus(), -1); - EXPECT_NEAR(bounds.min.x, -size / 2 - edgeLength / 2, 0.001); - EXPECT_NEAR(bounds.min.y, -size / 2 - edgeLength / 2, 0.001); - EXPECT_NEAR(bounds.min.z, -size / 2 - edgeLength / 2, 0.001); - EXPECT_NEAR(bounds.max.x, size / 2 + edgeLength / 2, 0.001); - EXPECT_NEAR(bounds.max.y, size / 2 + edgeLength / 2, 0.001); - EXPECT_NEAR(bounds.max.z, size / 2 + edgeLength / 2, 0.001); + const float outerBound = size / 2 + edgeLength / 2; + EXPECT_NEAR(bounds.min.x, -outerBound, precision); + EXPECT_NEAR(bounds.min.y, -outerBound, precision); + EXPECT_NEAR(bounds.min.z, -outerBound, precision); + EXPECT_NEAR(bounds.max.x, outerBound, precision); + EXPECT_NEAR(bounds.max.y, outerBound, precision); + EXPECT_NEAR(bounds.max.z, outerBound, precision); +} + +TEST(SDF, Surface) { + CubeVoid func; + const SDF voidSDF(func); + + const float size = 4; + const float edgeLength = 0.5; + + Manifold cubeVoid(voidSDF.LevelSet( + {glm::vec3(-size / 2), glm::vec3(size / 2)}, edgeLength)); + + Manifold cube = Manifold::Cube(glm::vec3(size), true); + cube -= cubeVoid; + Box bounds = cube.BoundingBox(); + const float precision = cube.Precision(); + if (options.exportModels) ExportMesh("cube.gltf", cube.GetMesh(), {}); + + EXPECT_TRUE(cube.IsManifold()); + EXPECT_EQ(cube.Genus(), 0); + auto prop = cube.GetProperties(); + EXPECT_NEAR(prop.volume, 8, 0.001); + EXPECT_NEAR(prop.surfaceArea, 24, 0.001); + EXPECT_NEAR(bounds.min.x, -1, precision); + EXPECT_NEAR(bounds.min.y, -1, precision); + EXPECT_NEAR(bounds.min.z, -1, precision); + EXPECT_NEAR(bounds.max.x, 1, precision); + EXPECT_NEAR(bounds.max.y, 1, precision); + EXPECT_NEAR(bounds.max.z, 1, precision); } \ No newline at end of file From fbcb707a4ed9d836d2316498588c915a4addf365 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Thu, 18 Aug 2022 11:55:39 -0700 Subject: [PATCH 24/39] updated gyroid_module --- samples/include/samples.h | 2 +- samples/src/gyroid_module.cpp | 24 +++++++++++++++++++----- 2 files changed, 20 insertions(+), 6 deletions(-) diff --git a/samples/include/samples.h b/samples/include/samples.h index 291dea844..439586dc8 100644 --- a/samples/include/samples.h +++ b/samples/include/samples.h @@ -50,6 +50,6 @@ Manifold TetPuzzle(float edgeLength, float gap, int nDivisions); Manifold Scallop(); -Manifold GyroidModule(int n); +Manifold GyroidModule(float size = 20, int n = 20); /** @} */ // end of Samples } // namespace manifold \ No newline at end of file diff --git a/samples/src/gyroid_module.cpp b/samples/src/gyroid_module.cpp index 5e083581e..3f1806425 100644 --- a/samples/src/gyroid_module.cpp +++ b/samples/src/gyroid_module.cpp @@ -23,6 +23,14 @@ struct Gyroid { return cos(p.x) * sin(p.y) + cos(p.y) * sin(p.z) + cos(p.z) * sin(p.x); } }; + +Manifold RhombicDodecahedron(float size) { + Manifold box = + Manifold::Cube(size * glm::sqrt(2.0f) * glm::vec3(1, 1, 2), true); + Manifold result = box.Rotate(45) ^ box.Rotate(0, 45); + return result ^ box.Rotate(90, 0, 45); +} + } // namespace namespace manifold { @@ -30,13 +38,19 @@ namespace manifold { /** * Creates a */ -Manifold GyroidModule(int n) { +Manifold GyroidModule(float size, int n) { Gyroid func; const SDF gyroidSDF(func); - const float size = glm::two_pi(); - Manifold gyroid( - gyroidSDF.LevelSet({glm::vec3(-size), glm::vec3(size)}, size / n, 0.5)); - return gyroid; + auto gyroid = [&](float level) { + const float period = glm::two_pi(); + return Manifold(gyroidSDF.LevelSet({glm::vec3(-period), glm::vec3(period)}, + period / n, level)) + .Scale(glm::vec3(size / period)); + }; + + Manifold result = (RhombicDodecahedron(size) ^ gyroid(-0.5)) - gyroid(0.5); + + return result.Rotate(-45, 0, 90).Translate({0, 0, size / glm::sqrt(2.0f)}); } } // namespace manifold \ No newline at end of file From ea0c7d0c499e3e699b41a10f6baf85e26a9a7bf9 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Thu, 18 Aug 2022 12:48:42 -0700 Subject: [PATCH 25/39] fixed gyroid module --- samples/src/gyroid_module.cpp | 6 +++--- test/samples_test.cpp | 15 +++++++++++---- 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/samples/src/gyroid_module.cpp b/samples/src/gyroid_module.cpp index 3f1806425..3e6cf8f26 100644 --- a/samples/src/gyroid_module.cpp +++ b/samples/src/gyroid_module.cpp @@ -27,8 +27,8 @@ struct Gyroid { Manifold RhombicDodecahedron(float size) { Manifold box = Manifold::Cube(size * glm::sqrt(2.0f) * glm::vec3(1, 1, 2), true); - Manifold result = box.Rotate(45) ^ box.Rotate(0, 45); - return result ^ box.Rotate(90, 0, 45); + Manifold result = box.Rotate(90, 45) ^ box.Rotate(90, 45, 90); + return result ^ box.Rotate(0, 0, 45); } } // namespace @@ -49,7 +49,7 @@ Manifold GyroidModule(float size, int n) { .Scale(glm::vec3(size / period)); }; - Manifold result = (RhombicDodecahedron(size) ^ gyroid(-0.5)) - gyroid(0.5); + Manifold result = (RhombicDodecahedron(size) ^ gyroid(-0.4)) - gyroid(0.4); return result.Rotate(-45, 0, 90).Translate({0, 0, size / glm::sqrt(2.0f)}); } diff --git a/test/samples_test.cpp b/test/samples_test.cpp index 9a0ba8cfa..d7ed73260 100644 --- a/test/samples_test.cpp +++ b/test/samples_test.cpp @@ -163,11 +163,18 @@ TEST(Samples, Bracelet) { } TEST(Samples, GyroidModule) { - Manifold gyroid = GyroidModule(20); + const float size = 20; + Manifold gyroid = GyroidModule(size); CheckManifold(gyroid); - EXPECT_LE(gyroid.NumDegenerateTris(), 0); - EXPECT_EQ(gyroid.Genus(), 1); - if (options.exportModels) ExportMesh("gyroid.gltf", gyroid.GetMesh(), {}); + EXPECT_LE(gyroid.NumDegenerateTris(), 2); + EXPECT_EQ(gyroid.Genus(), 15); + + const Box bounds = gyroid.BoundingBox(); + const float precision = gyroid.Precision(); + EXPECT_NEAR(bounds.min.z, 0, precision); + EXPECT_NEAR(bounds.max.z, size * glm::sqrt(2.0f), precision); + if (options.exportModels) + ExportMesh("gyroidModule.gltf", gyroid.GetMesh(), {}); } TEST(Samples, Sponge1) { From 7a57b0f7be7258c45acf69c7f4239efc8afb2830 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Thu, 18 Aug 2022 13:09:39 -0700 Subject: [PATCH 26/39] added docs --- samples/src/gyroid_module.cpp | 6 +++++- src/sdf/include/sdf.h | 34 ++++++++++++++++++++++++++++++++++ 2 files changed, 39 insertions(+), 1 deletion(-) diff --git a/samples/src/gyroid_module.cpp b/samples/src/gyroid_module.cpp index 3e6cf8f26..1f4a46ecb 100644 --- a/samples/src/gyroid_module.cpp +++ b/samples/src/gyroid_module.cpp @@ -36,7 +36,11 @@ Manifold RhombicDodecahedron(float size) { namespace manifold { /** - * Creates a + * Creates a rhombic dodecahedral module of a gyroid manifold, which can be + * assembled together to tile space continuously. This one is designed to be + * 3D-printable, as it is oriented with minimal overhangs. This sample + * demonstrates the use of a Signed Distance Function (SDF) to create smooth, + * complex manifolds. */ Manifold GyroidModule(float size, int n) { Gyroid func; diff --git a/src/sdf/include/sdf.h b/src/sdf/include/sdf.h index 5e314db5b..7b692b806 100644 --- a/src/sdf/include/sdf.h +++ b/src/sdf/include/sdf.h @@ -389,12 +389,46 @@ void RemoveUnreferencedVerts(VecDH& vertPos, template class SDF { public: + /** + * Create a signed-distance function object, which can produce Manifold level + * sets. + * + * @param sdf The signed-distance functor, containing this function signature: + * __host__ __device__ float operator()(glm::vec3 point), which returns the + * signed distance of a given point in R^3. Positive values are inside, + * negative outside. The __host__ __device__ is only needed if you compile for + * CUDA. If you are using a large grid, the advantage of a GPU speedup is + * quite significant. + */ SDF(Func sdf) : sdf_{sdf} {} + /** + * A convenience function for calling the given SDF on a single point. + * + * @param point + */ inline __host__ __device__ float operator()(glm::vec3 point) const { return sdf_(point); } + /** + * Constructs a level-set Mesh from the SDF. This uses a form of Marching + * Tetrahedra (akin to Marching Cubes, but better for manifoldness). Instead + * of using a cubic grid, it uses a body-centered cubic grid (two shifted + * cubic grids). This means if your function's interior exceeds the given + * bounds, you will see a kind of egg-crate shape closing off the manifold, + * which is due to the underlying grid. + * + * @param bounds An axis-aligned box that defines the extent of the grid. + * @param edgeLength Approximate maximum edge length of the triangles in the + * final result. This affects grid spacing, and hence has a strong effect on + * performance. + * @param level You can inset your Mesh by using a positive value, or outset + * it with a negative value. + * @return Mesh This class does not depend on Manifold, so it just returns a + * Mesh, but it is guaranteed to be manifold and so can always be used as + * input to the Manifold contructor for further operations. + */ inline Mesh LevelSet(Box bounds, float edgeLength, float level = 0) const { Mesh out; From cd57ffee733006c42d3e83af366f22efb58e7b75 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Mon, 22 Aug 2022 23:09:29 -0700 Subject: [PATCH 27/39] addressing feedback --- src/sdf/include/sdf.h | 50 +++++++++++++++---------------------------- 1 file changed, 17 insertions(+), 33 deletions(-) diff --git a/src/sdf/include/sdf.h b/src/sdf/include/sdf.h index 7b692b806..b6b1980d7 100644 --- a/src/sdf/include/sdf.h +++ b/src/sdf/include/sdf.h @@ -31,9 +31,8 @@ __host__ __device__ Uint64 AtomicCAS(Uint64& target, Uint64 compare, return atomicCAS(&target, compare, val); #else std::atomic& tar = reinterpret_cast&>(target); - Uint64 old_val = tar.load(); - tar.compare_exchange_weak(compare, val); - return old_val; + tar.compare_exchange_strong(compare, val); + return compare; #endif } @@ -145,7 +144,6 @@ __host__ __device__ glm::ivec4 DecodeMorton(Uint64 code) { struct GridVert { Uint64 key = kOpen; float distance = NAN; - int edgeIndex = -1; int edgeVerts[7] = {-1, -1, -1, -1, -1, -1, -1}; __host__ __device__ int Inside() const { return distance > 0 ? 1 : -1; } @@ -225,19 +223,19 @@ struct ComputeVerts { const float level; inline __host__ __device__ glm::vec3 Position(glm::ivec4 gridIndex) const { - return origin + spacing * (-0.5f + glm::vec3(gridIndex) + - (gridIndex.w == 1 ? 0.5f : 0.0f) * glm::vec3(1)); + return origin + + spacing * (glm::vec3(gridIndex) + (gridIndex.w == 1 ? 0.0f : -0.5f)); } inline __host__ __device__ float BoundedSDF(glm::ivec4 gridIndex) const { const float d = sdf(Position(gridIndex)) - level; const glm::ivec3 xyz(gridIndex); - if (glm::any(glm::equal(xyz, glm::ivec3(0))) || - glm::any(glm::greaterThanEqual(xyz, gridSize)) || - (gridIndex.w == 1 && - glm::any(glm::greaterThanEqual(xyz, gridSize - 1)))) - return glm::min(d, 0.0f); + const bool onLowerBound = glm::any(glm::equal(xyz, glm::ivec3(0))); + const bool onUpperBound = glm::any(glm::greaterThanEqual(xyz, gridSize)); + const bool onHalfBound = + gridIndex.w == 1 && glm::any(glm::greaterThanEqual(xyz, gridSize - 1)); + if (onLowerBound || onUpperBound || onHalfBound) return glm::min(d, 0.0f); return d; } @@ -254,36 +252,23 @@ struct ComputeVerts { gridVert.distance = BoundedSDF(gridIndex); bool keep = false; - float minDist2 = 0.25 * 0.25; for (int i = 0; i < 14; ++i) { glm::ivec4 neighborIndex = gridIndex + Neighbors(i); if (neighborIndex.w == 2) { neighborIndex += 1; neighborIndex.w = 0; } - const glm::vec3 iPos = Position(neighborIndex); const float val = BoundedSDF(neighborIndex); - if (val * gridVert.distance > 0 || (val == 0 && gridVert.distance == 0)) - continue; + if ((val > 0) == (gridVert.distance > 0)) continue; keep = true; - // Record the nearest intersection of all 14 edges, only if it is close - // enough to allow this gridVert to safely move to it without inverting - // any tetrahedra. - const glm::vec3 interp = (val * position - gridVert.distance * iPos) / - (val - gridVert.distance); - const glm::vec3 delta = interp - position; - const float dist2 = glm::dot(delta, delta); - if (dist2 < minDist2) { - // gridVert.edgeIndex = i; - minDist2 = dist2; - } - // These seven edges are uniquely owned by this gridVert; any of them // which intersect the surface create a vert. if (i < 7) { const int idx = AtomicAdd(*vertIndex, 1); - vertPos[idx] = interp; + vertPos[idx] = + (val * position - gridVert.distance * Position(neighborIndex)) / + (val - gridVert.distance); gridVert.edgeVerts[i] = idx; } } @@ -446,8 +431,8 @@ class SDF { VecDH vertPos(gridVerts.Size() * 7); VecDH index(1, 0); - thrust::for_each_n( - countAt(0), maxMorton + 1, + for_each_n( + autoPolicy(maxMorton), countAt(0), maxMorton + 1, ComputeVerts({vertPos.ptrD(), index.ptrD(), gridVerts.D(), sdf_, bounds.min, gridSize + 1, spacing, level})); vertPos.resize(index[0]); @@ -455,9 +440,8 @@ class SDF { VecDH triVerts(gridVerts.Entries() * 12); // worst case index[0] = 0; - thrust::for_each_n( - countAt(0), gridVerts.Size(), - BuildTris({triVerts.ptrD(), index.ptrD(), gridVerts.D()})); + for_each_n(autoPolicy(gridVerts.Size()), countAt(0), gridVerts.Size(), + BuildTris({triVerts.ptrD(), index.ptrD(), gridVerts.D()})); triVerts.resize(index[0]); RemoveUnreferencedVerts(vertPos, triVerts); From c1f825fe2440a255a0d2d6932bf8f9e5bc51c898 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Mon, 22 Aug 2022 23:30:23 -0700 Subject: [PATCH 28/39] fix cuda detection --- src/manifold/src/impl.cpp | 3 --- src/sdf/include/sdf.h | 5 +++-- src/utilities/include/utils.h | 11 ----------- src/utilities/src/detect_cuda.cpp | 4 ++++ 4 files changed, 7 insertions(+), 16 deletions(-) diff --git a/src/manifold/src/impl.cpp b/src/manifold/src/impl.cpp index 655bb8ed2..f11b43ae2 100644 --- a/src/manifold/src/impl.cpp +++ b/src/manifold/src/impl.cpp @@ -283,9 +283,6 @@ Manifold::Impl::Impl(const Mesh& mesh, const std::vector& properties, const std::vector& propertyTolerance) : vertPos_(mesh.vertPos), halfedgeTangent_(mesh.halfedgeTangent) { -#ifdef MANIFOLD_DEBUG - CheckDevice(); -#endif CalculateBBox(); if (!IsFinite()) { MarkFailure(Error::NON_FINITE_VERTEX); diff --git a/src/sdf/include/sdf.h b/src/sdf/include/sdf.h index b6b1980d7..c1dadcefc 100644 --- a/src/sdf/include/sdf.h +++ b/src/sdf/include/sdf.h @@ -423,6 +423,7 @@ class SDF { const glm::vec3 spacing = dim / (glm::vec3(gridSize)); const int maxMorton = MortonCode(glm::ivec4(gridSize + 1, 1)); + const auto policy = autoPolicy(maxMorton); const int tableSize = glm::min( 2 * maxMorton, static_cast(100 * glm::pow(maxMorton, 0.667))); @@ -432,7 +433,7 @@ class SDF { VecDH index(1, 0); for_each_n( - autoPolicy(maxMorton), countAt(0), maxMorton + 1, + policy, countAt(0), maxMorton + 1, ComputeVerts({vertPos.ptrD(), index.ptrD(), gridVerts.D(), sdf_, bounds.min, gridSize + 1, spacing, level})); vertPos.resize(index[0]); @@ -440,7 +441,7 @@ class SDF { VecDH triVerts(gridVerts.Entries() * 12); // worst case index[0] = 0; - for_each_n(autoPolicy(gridVerts.Size()), countAt(0), gridVerts.Size(), + for_each_n(policy, countAt(0), gridVerts.Size(), BuildTris({triVerts.ptrD(), index.ptrD(), gridVerts.D()})); triVerts.resize(index[0]); diff --git a/src/utilities/include/utils.h b/src/utilities/include/utils.h index 2ac530310..5e813156c 100644 --- a/src/utilities/include/utils.h +++ b/src/utilities/include/utils.h @@ -45,17 +45,6 @@ inline void MemUsage() { #endif } -inline void CheckDevice() { -#if THRUST_DEVICE_SYSTEM == THRUST_DEVICE_SYSTEM_CUDA - if (CUDA_ENABLED == -1) check_cuda_available(); - if (CUDA_ENABLED) { - cudaError_t error = cudaGetLastError(); - if (error != cudaSuccess) - throw std::runtime_error(cudaGetErrorString(error)); - } -#endif -} - struct Timer { #if THRUST_DEVICE_SYSTEM == THRUST_DEVICE_SYSTEM_CUDA cudaEvent_t start, end; diff --git a/src/utilities/src/detect_cuda.cpp b/src/utilities/src/detect_cuda.cpp index 878228245..e432909ec 100644 --- a/src/utilities/src/detect_cuda.cpp +++ b/src/utilities/src/detect_cuda.cpp @@ -21,6 +21,10 @@ void check_cuda_available() { int device_count = 0; cudaError_t error = cudaGetDeviceCount(&device_count); CUDA_ENABLED = device_count != 0; + if (CUDA_ENABLED) { + cudaError_t error = cudaGetLastError(); + CUDA_ENABLED = error == cudaSuccess; + } } } // namespace manifold #else From 4d890884cdb2bbce6fea4eb5c5aa1a55d4188bfc Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Tue, 23 Aug 2022 09:32:39 -0700 Subject: [PATCH 29/39] clean up detect cuda --- src/utilities/include/par.h | 9 ++------- src/utilities/include/vec_dh.h | 7 +++---- src/utilities/src/detect_cuda.cpp | 22 ++++++++++++---------- 3 files changed, 17 insertions(+), 21 deletions(-) diff --git a/src/utilities/include/par.h b/src/utilities/include/par.h index 5fbcafa0a..23b5b79e4 100644 --- a/src/utilities/include/par.h +++ b/src/utilities/include/par.h @@ -41,12 +41,7 @@ namespace manifold { -void check_cuda_available(); -#ifdef MANIFOLD_USE_CUDA -extern int CUDA_ENABLED; -#else -constexpr int CUDA_ENABLED = 0; -#endif +bool CudaEnabled(); enum class ExecutionPolicy { ParUnseq, @@ -67,7 +62,7 @@ inline ExecutionPolicy autoPolicy(int size) { if (size <= (1 << 12)) { return Seq; } - if (size <= (1 << 16) || CUDA_ENABLED != 1) { + if (size <= (1 << 16) || !CudaEnabled()) { return Par; } return ParUnseq; diff --git a/src/utilities/include/vec_dh.h b/src/utilities/include/vec_dh.h index 89764ef0e..a1ed004cc 100644 --- a/src/utilities/include/vec_dh.h +++ b/src/utilities/include/vec_dh.h @@ -238,8 +238,7 @@ class ManagedVec { static void mallocManaged(T **ptr, size_t bytes) { #ifdef MANIFOLD_USE_CUDA - if (CUDA_ENABLED == -1) check_cuda_available(); - if (CUDA_ENABLED) + if (CudaEnabled()) cudaMallocManaged(reinterpret_cast(ptr), bytes); else #endif @@ -248,7 +247,7 @@ class ManagedVec { static void freeManaged(T *ptr) { #ifdef MANIFOLD_USE_CUDA - if (CUDA_ENABLED) + if (CudaEnabled()) cudaFree(ptr); else #endif @@ -257,7 +256,7 @@ class ManagedVec { static void prefetch(T *ptr, int bytes, bool onHost) { #ifdef MANIFOLD_USE_CUDA - if (bytes > 0 && CUDA_ENABLED) + if (bytes > 0 && CudaEnabled()) cudaMemPrefetchAsync(ptr, std::min(bytes, DEVICE_MAX_BYTES), onHost ? cudaCpuDeviceId : 0); #endif diff --git a/src/utilities/src/detect_cuda.cpp b/src/utilities/src/detect_cuda.cpp index e432909ec..8e9b284ae 100644 --- a/src/utilities/src/detect_cuda.cpp +++ b/src/utilities/src/detect_cuda.cpp @@ -15,20 +15,22 @@ #ifdef MANIFOLD_USE_CUDA #include +namespace { +int CUDA_DEVICES = -1; +} namespace manifold { -int CUDA_ENABLED = -1; -void check_cuda_available() { - int device_count = 0; - cudaError_t error = cudaGetDeviceCount(&device_count); - CUDA_ENABLED = device_count != 0; - if (CUDA_ENABLED) { - cudaError_t error = cudaGetLastError(); - CUDA_ENABLED = error == cudaSuccess; - } + +bool CudaEnabled() { + if (CUDA_DEVICES >= 0) return CUDA_DEVICES > 0; + + cudaError_t error = cudaGetDeviceCount(&CUDA_DEVICES); + if (error != cudaSuccess) CUDA_DEVICES = 0; + + return CUDA_DEVICES > 0; } } // namespace manifold #else namespace manifold { -void check_cuda_available() {} +bool CudaEnabled() { return false; } } // namespace manifold #endif From a2b282f4a3744f0b05be6b960600cd877beb4acf Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Tue, 23 Aug 2022 09:57:52 -0700 Subject: [PATCH 30/39] fixed CUDA failure --- test/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 8e4bd53ea..a18369f58 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -21,7 +21,7 @@ enable_testing() set(SOURCE_FILES polygon_test.cpp mesh_test.cpp samples_test.cpp sdf_test.cpp test_main.cpp) add_executable(${PROJECT_NAME} ${SOURCE_FILES}) -target_link_libraries(${PROJECT_NAME} polygon GTest::GTest manifold meshIO samples samplesGPU sdf) +target_link_libraries(${PROJECT_NAME} polygon GTest::GTest manifold meshIO samples samplesGPU) if(MANIFOLD_USE_CUDA) set_source_files_properties(sdf_test.cpp PROPERTIES LANGUAGE CUDA) From b91db7107cc6f196be2991d509a06e9b0c70bf8d Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Sat, 27 Aug 2022 10:13:13 -0700 Subject: [PATCH 31/39] consistent execution policy in Boolean --- src/manifold/src/boolean3.cpp | 65 ++++++++++++++++------------- src/manifold/src/boolean3.h | 1 + src/manifold/src/boolean_result.cpp | 51 ++++++++++------------ src/utilities/include/sparse.h | 9 ++-- 4 files changed, 63 insertions(+), 63 deletions(-) diff --git a/src/manifold/src/boolean3.cpp b/src/manifold/src/boolean3.cpp index 80ea9fcf8..d8eeed6a1 100644 --- a/src/manifold/src/boolean3.cpp +++ b/src/manifold/src/boolean3.cpp @@ -83,9 +83,9 @@ struct CopyFaceEdges { }; SparseIndices Filter11(const Manifold::Impl &inP, const Manifold::Impl &inQ, - const SparseIndices &p1q2, const SparseIndices &p2q1) { + const SparseIndices &p1q2, const SparseIndices &p2q1, + ExecutionPolicy policy) { SparseIndices p1q1(3 * p1q2.size() + 3 * p2q1.size()); - auto policy = autoPolicy(p1q2.size()); for_each_n(policy, zip(countAt(0), p1q2.begin(0), p1q2.begin(1)), p1q2.size(), CopyFaceEdges({p1q1.ptrDpq(), inQ.halfedge_.cptrD()})); @@ -94,7 +94,7 @@ SparseIndices Filter11(const Manifold::Impl &inP, const Manifold::Impl &inQ, p2q1.size(), CopyFaceEdges({p1q1.ptrDpq(), inP.halfedge_.cptrD()})); p1q1.SwapPQ(); - p1q1.Unique(); + p1q1.Unique(policy); return p1q1; } @@ -241,11 +241,12 @@ struct Kernel11 { std::tuple, VecDH> Shadow11(SparseIndices &p1q1, const Manifold::Impl &inP, const Manifold::Impl &inQ, - float expandP) { + float expandP, + ExecutionPolicy policy) { VecDH s11(p1q1.size()); VecDH xyzz11(p1q1.size()); - for_each_n(autoPolicy(p1q1.size()), + for_each_n(policy, zip(xyzz11.begin(), s11.begin(), p1q1.begin(0), p1q1.begin(1)), p1q1.size(), Kernel11({inP.vertPos_.cptrD(), inQ.vertPos_.cptrD(), @@ -336,14 +337,15 @@ struct Kernel02 { std::tuple, VecDH> Shadow02(const Manifold::Impl &inP, const Manifold::Impl &inQ, SparseIndices &p0q2, bool forward, - float expandP) { + float expandP, + ExecutionPolicy policy) { VecDH s02(p0q2.size()); VecDH z02(p0q2.size()); auto vertNormalP = forward ? inP.vertNormal_.cptrD() : inQ.vertNormal_.cptrD(); for_each_n( - autoPolicy(p0q2.size()), + policy, zip(s02.begin(), z02.begin(), p0q2.begin(!forward), p0q2.begin(forward)), p0q2.size(), Kernel02({inP.vertPos_.cptrD(), inQ.halfedge_.cptrD(), @@ -448,12 +450,12 @@ std::tuple, VecDH> Intersect12( const Manifold::Impl &inP, const Manifold::Impl &inQ, const VecDH &s02, const SparseIndices &p0q2, const VecDH &s11, const SparseIndices &p1q1, const VecDH &z02, const VecDH &xyzz11, - SparseIndices &p1q2, bool forward) { + SparseIndices &p1q2, bool forward, ExecutionPolicy policy) { VecDH x12(p1q2.size()); VecDH v12(p1q2.size()); for_each_n( - autoPolicy(p1q2.size()), + policy, zip(x12.begin(), v12.begin(), p1q2.begin(!forward), p1q2.begin(forward)), p1q2.size(), Kernel12({p0q2.ptrDpq(), s02.ptrD(), z02.cptrD(), p0q2.size(), @@ -467,11 +469,10 @@ std::tuple, VecDH> Intersect12( }; VecDH Winding03(const Manifold::Impl &inP, SparseIndices &p0q2, - VecDH &s02, bool reverse) { + VecDH &s02, bool reverse, ExecutionPolicy policy) { // verts that are not shadowed (not in p0q2) have winding number zero. VecDH w03(inP.NumVert(), 0); - auto policy = autoPolicy(p0q2.size()); if (!is_sorted(policy, p0q2.begin(reverse), p0q2.end(reverse))) sort_by_key(policy, p0q2.begin(reverse), p0q2.end(reverse), s02.begin()); VecDH w03val(w03.size()); @@ -481,11 +482,10 @@ VecDH Winding03(const Manifold::Impl &inP, SparseIndices &p0q2, thrust::pair>( policy, p0q2.begin(reverse), p0q2.end(reverse), s02.begin(), w03vert.begin(), w03val.begin()); - scatter(autoPolicy(endPair.second - w03val.begin()), w03val.begin(), - endPair.second, w03vert.begin(), w03.begin()); + scatter(policy, w03val.begin(), endPair.second, w03vert.begin(), w03.begin()); if (reverse) - transform(autoPolicy(w03.size()), w03.begin(), w03.end(), w03.begin(), + transform(policy, w03.begin(), w03.end(), w03.begin(), thrust::negate()); return w03; }; @@ -494,7 +494,10 @@ VecDH Winding03(const Manifold::Impl &inP, SparseIndices &p0q2, namespace manifold { Boolean3::Boolean3(const Manifold::Impl &inP, const Manifold::Impl &inQ, Manifold::OpType op) - : inP_(inP), inQ_(inQ), expandP_(op == Manifold::OpType::ADD ? 1.0 : -1.0) { + : inP_(inP), + inQ_(inQ), + expandP_(op == Manifold::OpType::ADD ? 1.0 : -1.0), + policy_(autoPolicy(glm::max(inP.NumEdge(), inQ.NumEdge()))) { // Symbolic perturbation: // Union -> expand inP // Difference, Intersection -> contract inP @@ -514,27 +517,29 @@ Boolean3::Boolean3(const Manifold::Impl &inP, const Manifold::Impl &inQ, // Level 3 // Find edge-triangle overlaps (broad phase) p1q2_ = inQ_.EdgeCollisions(inP_); - p1q2_.Sort(); + p2q1_ = inP_.EdgeCollisions(inQ_); + + policy_ = autoPolicy(glm::max(p1q2_.size(), p2q1_.size())); + p1q2_.Sort(policy_); PRINT("p1q2 size = " << p1q2_.size()); - p2q1_ = inP_.EdgeCollisions(inQ_); p2q1_.SwapPQ(); - p2q1_.Sort(); + p2q1_.Sort(policy_); PRINT("p2q1 size = " << p2q1_.size()); // Level 2 // Find vertices that overlap faces in XY-projection SparseIndices p0q2 = inQ.VertexCollisionsZ(inP.vertPos_); - p0q2.Sort(); + p0q2.Sort(policy_); PRINT("p0q2 size = " << p0q2.size()); SparseIndices p2q0 = inP.VertexCollisionsZ(inQ.vertPos_); p2q0.SwapPQ(); - p2q0.Sort(); + p2q0.Sort(policy_); PRINT("p2q0 size = " << p2q0.size()); // Find involved edge pairs from Level 3 - SparseIndices p1q1 = Filter11(inP_, inQ_, p1q2_, p2q1_); + SparseIndices p1q1 = Filter11(inP_, inQ_, p1q2_, p2q1_, policy_); PRINT("p1q1 size = " << p1q1.size()); #ifdef MANIFOLD_DEBUG @@ -548,37 +553,37 @@ Boolean3::Boolean3(const Manifold::Impl &inP, const Manifold::Impl &inQ, // each edge, keeping only those whose intersection exists. VecDH s11; VecDH xyzz11; - std::tie(s11, xyzz11) = Shadow11(p1q1, inP, inQ, expandP_); + std::tie(s11, xyzz11) = Shadow11(p1q1, inP, inQ, expandP_, policy_); PRINT("s11 size = " << s11.size()); // Build up Z-projection of vertices onto triangles, keeping only those that // fall inside the triangle. VecDH s02; VecDH z02; - std::tie(s02, z02) = Shadow02(inP, inQ, p0q2, true, expandP_); + std::tie(s02, z02) = Shadow02(inP, inQ, p0q2, true, expandP_, policy_); PRINT("s02 size = " << s02.size()); VecDH s20; VecDH z20; - std::tie(s20, z20) = Shadow02(inQ, inP, p2q0, false, expandP_); + std::tie(s20, z20) = Shadow02(inQ, inP, p2q0, false, expandP_, policy_); PRINT("s20 size = " << s20.size()); // Level 3 // Build up the intersection of the edges and triangles, keeping only those // that intersect, and record the direction the edge is passing through the // triangle. - std::tie(x12_, v12_) = - Intersect12(inP, inQ, s02, p0q2, s11, p1q1, z02, xyzz11, p1q2_, true); + std::tie(x12_, v12_) = Intersect12(inP, inQ, s02, p0q2, s11, p1q1, z02, + xyzz11, p1q2_, true, policy_); PRINT("x12 size = " << x12_.size()); - std::tie(x21_, v21_) = - Intersect12(inQ, inP, s20, p2q0, s11, p1q1, z20, xyzz11, p2q1_, false); + std::tie(x21_, v21_) = Intersect12(inQ, inP, s20, p2q0, s11, p1q1, z20, + xyzz11, p2q1_, false, policy_); PRINT("x21 size = " << x21_.size()); // Sum up the winding numbers of all vertices. - w03_ = Winding03(inP, p0q2, s02, false); + w03_ = Winding03(inP, p0q2, s02, false, policy_); - w30_ = Winding03(inQ, p2q0, s20, true); + w30_ = Winding03(inQ, p2q0, s20, true, policy_); #ifdef MANIFOLD_DEBUG intersections.Stop(); diff --git a/src/manifold/src/boolean3.h b/src/manifold/src/boolean3.h index 094197d73..5fffcbe8a 100644 --- a/src/manifold/src/boolean3.h +++ b/src/manifold/src/boolean3.h @@ -57,5 +57,6 @@ class Boolean3 { SparseIndices p1q2_, p2q1_; VecDH x12_, x21_, w03_, w30_; VecDH v12_, v21_; + ExecutionPolicy policy_; }; } // namespace manifold diff --git a/src/manifold/src/boolean_result.cpp b/src/manifold/src/boolean_result.cpp index a8b930f01..d16eb2522 100644 --- a/src/manifold/src/boolean_result.cpp +++ b/src/manifold/src/boolean_result.cpp @@ -76,14 +76,11 @@ std::tuple, VecDH> SizeOutput( Manifold::Impl &outR, const Manifold::Impl &inP, const Manifold::Impl &inQ, const VecDH &i03, const VecDH &i30, const VecDH &i12, const VecDH &i21, const SparseIndices &p1q2, const SparseIndices &p2q1, - bool invertQ) { + bool invertQ, ExecutionPolicy policy) { VecDH sidesPerFacePQ(inP.NumTri() + inQ.NumTri(), 0); auto sidesPerFaceP = sidesPerFacePQ.ptrD(); auto sidesPerFaceQ = sidesPerFacePQ.ptrD() + inP.NumTri(); - auto policy = - autoPolicy(std::max(inP.halfedge_.size(), inQ.halfedge_.size())); - for_each(policy, inP.halfedge_.begin(), inP.halfedge_.end(), CountVerts({sidesPerFaceP, i03.cptrD()})); for_each(policy, inQ.halfedge_.begin(), inQ.halfedge_.end(), @@ -403,9 +400,9 @@ struct DuplicateHalfedges { void AppendWholeEdges(Manifold::Impl &outR, VecDH &facePtrR, VecDH &halfedgeRef, const Manifold::Impl &inP, const VecDH wholeHalfedgeP, const VecDH &i03, - const VecDH &vP2R, const int *faceP2R, - bool forward) { - for_each_n(autoPolicy(inP.halfedge_.size()), + const VecDH &vP2R, const int *faceP2R, bool forward, + ExecutionPolicy policy) { + for_each_n(policy, zip(wholeHalfedgeP.begin(), inP.halfedge_.begin(), countAt(0)), inP.halfedge_.size(), DuplicateHalfedges({outR.halfedge_.ptrD(), halfedgeRef.ptrD(), @@ -483,7 +480,7 @@ struct CreateBarycentric { std::pair, VecDH> CalculateMeshRelation( Manifold::Impl &outR, const VecDH &halfedgeRef, const Manifold::Impl &inP, const Manifold::Impl &inQ, int firstNewVert, - int numFaceR, bool invertQ) { + int numFaceR, bool invertQ, ExecutionPolicy policy) { outR.meshRelation_.barycentric.resize(outR.halfedge_.size()); VecDH faceRef(numFaceR); VecDH halfedgeBary(halfedgeRef.size()); @@ -491,7 +488,7 @@ std::pair, VecDH> CalculateMeshRelation( const int offsetQ = Manifold::Impl::meshIDCounter_; VecDH idx(1, 0); for_each_n( - autoPolicy(halfedgeRef.size()), + policy, zip(halfedgeBary.begin(), halfedgeRef.begin(), outR.halfedge_.cbegin()), halfedgeRef.size(), CreateBarycentric( @@ -540,28 +537,26 @@ Manifold::Impl Boolean3::Result(Manifold::OpType op) const { VecDH i21(x21_.size()); VecDH i03(w03_.size()); VecDH i30(w30_.size()); - auto policy = autoPolicy(std::max(std::max(x12_.size(), x21_.size()), - std::max(w03_.size(), w30_.size()))); - transform(policy, x12_.begin(), x12_.end(), i12.begin(), c3 * _1); - transform(policy, x21_.begin(), x21_.end(), i21.begin(), c3 * _1); - transform(policy, w03_.begin(), w03_.end(), i03.begin(), c1 + c3 * _1); - transform(policy, w30_.begin(), w30_.end(), i30.begin(), c2 + c3 * _1); + transform(policy_, x12_.begin(), x12_.end(), i12.begin(), c3 * _1); + transform(policy_, x21_.begin(), x21_.end(), i21.begin(), c3 * _1); + transform(policy_, w03_.begin(), w03_.end(), i03.begin(), c1 + c3 * _1); + transform(policy_, w30_.begin(), w30_.end(), i30.begin(), c2 + c3 * _1); VecDH vP2R(inP_.NumVert()); - exclusive_scan(policy, i03.begin(), i03.end(), vP2R.begin(), 0, AbsSum()); + exclusive_scan(policy_, i03.begin(), i03.end(), vP2R.begin(), 0, AbsSum()); int numVertR = AbsSum()(vP2R.back(), i03.back()); const int nPv = numVertR; VecDH vQ2R(inQ_.NumVert()); - exclusive_scan(policy, i30.begin(), i30.end(), vQ2R.begin(), numVertR, + exclusive_scan(policy_, i30.begin(), i30.end(), vQ2R.begin(), numVertR, AbsSum()); numVertR = AbsSum()(vQ2R.back(), i30.back()); const int nQv = numVertR - nPv; VecDH v12R(v12_.size()); if (v12_.size() > 0) { - exclusive_scan(policy, i12.begin(), i12.end(), v12R.begin(), numVertR, + exclusive_scan(policy_, i12.begin(), i12.end(), v12R.begin(), numVertR, AbsSum()); numVertR = AbsSum()(v12R.back(), i12.back()); } @@ -569,7 +564,7 @@ Manifold::Impl Boolean3::Result(Manifold::OpType op) const { VecDH v21R(v21_.size()); if (v21_.size() > 0) { - exclusive_scan(policy, i21.begin(), i21.end(), v21R.begin(), numVertR, + exclusive_scan(policy_, i21.begin(), i21.end(), v21R.begin(), numVertR, AbsSum()); numVertR = AbsSum()(v21R.back(), i21.back()); } @@ -585,14 +580,14 @@ Manifold::Impl Boolean3::Result(Manifold::OpType op) const { outR.vertPos_.resize(numVertR); // Add vertices, duplicating for inclusion numbers not in [-1, 1]. // Retained vertices from P and Q: - for_each_n(policy, zip(i03.begin(), vP2R.begin(), inP_.vertPos_.begin()), + for_each_n(policy_, zip(i03.begin(), vP2R.begin(), inP_.vertPos_.begin()), inP_.NumVert(), DuplicateVerts({outR.vertPos_.ptrD()})); - for_each_n(policy, zip(i30.begin(), vQ2R.begin(), inQ_.vertPos_.begin()), + for_each_n(policy_, zip(i30.begin(), vQ2R.begin(), inQ_.vertPos_.begin()), inQ_.NumVert(), DuplicateVerts({outR.vertPos_.ptrD()})); // New vertices created from intersections: - for_each_n(policy, zip(i12.begin(), v12R.begin(), v12_.begin()), i12.size(), + for_each_n(policy_, zip(i12.begin(), v12R.begin(), v12_.begin()), i12.size(), DuplicateVerts({outR.vertPos_.ptrD()})); - for_each_n(policy, zip(i21.begin(), v21R.begin(), v21_.begin()), i21.size(), + for_each_n(policy_, zip(i21.begin(), v21R.begin(), v21_.begin()), i21.size(), DuplicateVerts({outR.vertPos_.ptrD()})); PRINT(nPv << " verts from inP"); @@ -617,8 +612,8 @@ Manifold::Impl Boolean3::Result(Manifold::OpType op) const { // Level 4 VecDH faceEdge; VecDH facePQ2R; - std::tie(faceEdge, facePQ2R) = - SizeOutput(outR, inP_, inQ_, i03, i30, i12, i21, p1q2_, p2q1_, invertQ); + std::tie(faceEdge, facePQ2R) = SizeOutput( + outR, inP_, inQ_, i03, i30, i12, i21, p1q2_, p2q1_, invertQ, policy_); const int numFaceR = faceEdge.size() - 1; // This gets incremented for each halfedge that's added to a face so that the @@ -640,14 +635,14 @@ Manifold::Impl Boolean3::Result(Manifold::OpType op) const { inP_.NumTri()); AppendWholeEdges(outR, facePtrR, halfedgeRef, inP_, wholeHalfedgeP, i03, vP2R, - facePQ2R.cptrD(), true); + facePQ2R.cptrD(), true, policy_); AppendWholeEdges(outR, facePtrR, halfedgeRef, inQ_, wholeHalfedgeQ, i30, vQ2R, - facePQ2R.cptrD() + inP_.NumTri(), false); + facePQ2R.cptrD() + inP_.NumTri(), false, policy_); VecDH faceRef; VecDH halfedgeBary; std::tie(faceRef, halfedgeBary) = CalculateMeshRelation( - outR, halfedgeRef, inP_, inQ_, nPv + nQv, numFaceR, invertQ); + outR, halfedgeRef, inP_, inQ_, nPv + nQv, numFaceR, invertQ, policy_); #ifdef MANIFOLD_DEBUG assemble.Stop(); diff --git a/src/utilities/include/sparse.h b/src/utilities/include/sparse.h index 386ad3f66..ef40ea882 100644 --- a/src/utilities/include/sparse.h +++ b/src/utilities/include/sparse.h @@ -62,18 +62,17 @@ class SparseIndices { int size() const { return p.size(); } void SwapPQ() { p.swap(q); } - void Sort() { sort(autoPolicy(size()), beginPQ(), endPQ()); } + void Sort(ExecutionPolicy policy) { sort(policy, beginPQ(), endPQ()); } void Resize(int size) { p.resize(size, -1); q.resize(size, -1); } - void Unique() { - Sort(); + void Unique(ExecutionPolicy policy) { + Sort(policy); int newSize = - unique(autoPolicy(size()), beginPQ(), endPQ()) - - beginPQ(); + unique(policy, beginPQ(), endPQ()) - beginPQ(); Resize(newSize); } From 94a3c697b7ba86c21060733eea5bf2cffbd04e02 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Sat, 27 Aug 2022 10:37:01 -0700 Subject: [PATCH 32/39] remove SDF wrapper class --- samples/src/gyroid_module.cpp | 7 +- src/sdf/include/sdf.h | 141 +++++++++++++++------------------- test/CMakeLists.txt | 2 +- test/sdf_test.cpp | 17 ++-- 4 files changed, 68 insertions(+), 99 deletions(-) diff --git a/samples/src/gyroid_module.cpp b/samples/src/gyroid_module.cpp index 1f4a46ecb..cf8b6ed89 100644 --- a/samples/src/gyroid_module.cpp +++ b/samples/src/gyroid_module.cpp @@ -43,13 +43,10 @@ namespace manifold { * complex manifolds. */ Manifold GyroidModule(float size, int n) { - Gyroid func; - const SDF gyroidSDF(func); - auto gyroid = [&](float level) { const float period = glm::two_pi(); - return Manifold(gyroidSDF.LevelSet({glm::vec3(-period), glm::vec3(period)}, - period / n, level)) + return Manifold(LevelSet(Gyroid(), {glm::vec3(-period), glm::vec3(period)}, + period / n, level)) .Scale(glm::vec3(size / period)); }; diff --git a/src/sdf/include/sdf.h b/src/sdf/include/sdf.h index c1dadcefc..bcb54ae04 100644 --- a/src/sdf/include/sdf.h +++ b/src/sdf/include/sdf.h @@ -371,89 +371,68 @@ void RemoveUnreferencedVerts(VecDH& vertPos, /** @addtogroup Core * @{ */ + +/** + * Constructs a level-set Mesh from the input Signed-Distance Function (SDF). + * This uses a form of Marching Tetrahedra (akin to Marching Cubes, but better + * for manifoldness). Instead of using a cubic grid, it uses a body-centered + * cubic grid (two shifted cubic grids). This means if your function's interior + * exceeds the given bounds, you will see a kind of egg-crate shape closing off + * the manifold, which is due to the underlying grid. + * + * @param sdf The signed-distance functor, containing this function signature: + * __host__ __device__ float operator()(glm::vec3 point), which returns the + * signed distance of a given point in R^3. Positive values are inside, + * negative outside. The __host__ __device__ is only needed if you compile for + * CUDA. If you are using a large grid, the advantage of a GPU speedup is + * quite significant. + * @param bounds An axis-aligned box that defines the extent of the grid. + * @param edgeLength Approximate maximum edge length of the triangles in the + * final result. This affects grid spacing, and hence has a strong effect on + * performance. + * @param level You can inset your Mesh by using a positive value, or outset + * it with a negative value. + * @return Mesh This class does not depend on Manifold, so it just returns a + * Mesh, but it is guaranteed to be manifold and so can always be used as + * input to the Manifold contructor for further operations. + */ template -class SDF { - public: - /** - * Create a signed-distance function object, which can produce Manifold level - * sets. - * - * @param sdf The signed-distance functor, containing this function signature: - * __host__ __device__ float operator()(glm::vec3 point), which returns the - * signed distance of a given point in R^3. Positive values are inside, - * negative outside. The __host__ __device__ is only needed if you compile for - * CUDA. If you are using a large grid, the advantage of a GPU speedup is - * quite significant. - */ - SDF(Func sdf) : sdf_{sdf} {} - - /** - * A convenience function for calling the given SDF on a single point. - * - * @param point - */ - inline __host__ __device__ float operator()(glm::vec3 point) const { - return sdf_(point); - } +inline Mesh LevelSet(Func sdf, Box bounds, float edgeLength, float level = 0) { + Mesh out; - /** - * Constructs a level-set Mesh from the SDF. This uses a form of Marching - * Tetrahedra (akin to Marching Cubes, but better for manifoldness). Instead - * of using a cubic grid, it uses a body-centered cubic grid (two shifted - * cubic grids). This means if your function's interior exceeds the given - * bounds, you will see a kind of egg-crate shape closing off the manifold, - * which is due to the underlying grid. - * - * @param bounds An axis-aligned box that defines the extent of the grid. - * @param edgeLength Approximate maximum edge length of the triangles in the - * final result. This affects grid spacing, and hence has a strong effect on - * performance. - * @param level You can inset your Mesh by using a positive value, or outset - * it with a negative value. - * @return Mesh This class does not depend on Manifold, so it just returns a - * Mesh, but it is guaranteed to be manifold and so can always be used as - * input to the Manifold contructor for further operations. - */ - inline Mesh LevelSet(Box bounds, float edgeLength, float level = 0) const { - Mesh out; - - const glm::vec3 dim = bounds.Size(); - const float maxDim = std::max(dim[0], std::max(dim[1], dim[2])); - const glm::ivec3 gridSize(dim / edgeLength); - const glm::vec3 spacing = dim / (glm::vec3(gridSize)); - - const int maxMorton = MortonCode(glm::ivec4(gridSize + 1, 1)); - const auto policy = autoPolicy(maxMorton); - - const int tableSize = glm::min( - 2 * maxMorton, static_cast(100 * glm::pow(maxMorton, 0.667))); - HashTable gridVerts(tableSize); - - VecDH vertPos(gridVerts.Size() * 7); - VecDH index(1, 0); - - for_each_n( - policy, countAt(0), maxMorton + 1, - ComputeVerts({vertPos.ptrD(), index.ptrD(), gridVerts.D(), sdf_, - bounds.min, gridSize + 1, spacing, level})); - vertPos.resize(index[0]); - - VecDH triVerts(gridVerts.Entries() * 12); // worst case - - index[0] = 0; - for_each_n(policy, countAt(0), gridVerts.Size(), - BuildTris({triVerts.ptrD(), index.ptrD(), gridVerts.D()})); - triVerts.resize(index[0]); - - RemoveUnreferencedVerts(vertPos, triVerts); - - out.vertPos.insert(out.vertPos.end(), vertPos.begin(), vertPos.end()); - out.triVerts.insert(out.triVerts.end(), triVerts.begin(), triVerts.end()); - return out; - } + const glm::vec3 dim = bounds.Size(); + const float maxDim = std::max(dim[0], std::max(dim[1], dim[2])); + const glm::ivec3 gridSize(dim / edgeLength); + const glm::vec3 spacing = dim / (glm::vec3(gridSize)); - private: - const Func sdf_; -}; + const int maxMorton = MortonCode(glm::ivec4(gridSize + 1, 1)); + const auto policy = autoPolicy(maxMorton); + + const int tableSize = glm::min( + 2 * maxMorton, static_cast(100 * glm::pow(maxMorton, 0.667))); + HashTable gridVerts(tableSize); + + VecDH vertPos(gridVerts.Size() * 7); + VecDH index(1, 0); + + for_each_n( + policy, countAt(0), maxMorton + 1, + ComputeVerts({vertPos.ptrD(), index.ptrD(), gridVerts.D(), sdf, + bounds.min, gridSize + 1, spacing, level})); + vertPos.resize(index[0]); + + VecDH triVerts(gridVerts.Entries() * 12); // worst case + + index[0] = 0; + for_each_n(policy, countAt(0), gridVerts.Size(), + BuildTris({triVerts.ptrD(), index.ptrD(), gridVerts.D()})); + triVerts.resize(index[0]); + + RemoveUnreferencedVerts(vertPos, triVerts); + + out.vertPos.insert(out.vertPos.end(), vertPos.begin(), vertPos.end()); + out.triVerts.insert(out.triVerts.end(), triVerts.begin(), triVerts.end()); + return out; +} /** @} */ } // namespace manifold \ No newline at end of file diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index a18369f58..f7a550b64 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -19,7 +19,7 @@ add_subdirectory(meshIO) enable_testing() -set(SOURCE_FILES polygon_test.cpp mesh_test.cpp samples_test.cpp sdf_test.cpp test_main.cpp) +set(SOURCE_FILES polygon_test.cpp mesh_test.cpp sdf_test.cpp samples_test.cpp test_main.cpp) add_executable(${PROJECT_NAME} ${SOURCE_FILES}) target_link_libraries(${PROJECT_NAME} polygon GTest::GTest manifold meshIO samples samplesGPU) diff --git a/test/sdf_test.cpp b/test/sdf_test.cpp index b361d33e6..1a3ea31dd 100644 --- a/test/sdf_test.cpp +++ b/test/sdf_test.cpp @@ -31,8 +31,7 @@ struct CubeVoid { }; TEST(SDF, CubeVoid) { - CubeVoid func; - const SDF voidSDF(func); + CubeVoid voidSDF; EXPECT_EQ(voidSDF({0, 0, 0}), -1); EXPECT_EQ(voidSDF({0, 0, 1}), 0); @@ -45,14 +44,11 @@ TEST(SDF, CubeVoid) { } TEST(SDF, Bounds) { - CubeVoid func; - const SDF voidSDF(func); - const float size = 4; const float edgeLength = 0.5; - Manifold cubeVoid(voidSDF.LevelSet( - {glm::vec3(-size / 2), glm::vec3(size / 2)}, edgeLength)); + Manifold cubeVoid(LevelSet( + CubeVoid(), {glm::vec3(-size / 2), glm::vec3(size / 2)}, edgeLength)); Box bounds = cubeVoid.BoundingBox(); const float precision = cubeVoid.Precision(); if (options.exportModels) ExportMesh("cubeVoid.gltf", cubeVoid.GetMesh(), {}); @@ -69,14 +65,11 @@ TEST(SDF, Bounds) { } TEST(SDF, Surface) { - CubeVoid func; - const SDF voidSDF(func); - const float size = 4; const float edgeLength = 0.5; - Manifold cubeVoid(voidSDF.LevelSet( - {glm::vec3(-size / 2), glm::vec3(size / 2)}, edgeLength)); + Manifold cubeVoid(LevelSet( + CubeVoid(), {glm::vec3(-size / 2), glm::vec3(size / 2)}, edgeLength)); Manifold cube = Manifold::Cube(glm::vec3(size), true); cube -= cubeVoid; From c626f4e8f1583654d7be364e0f9e982f3a5fa76d Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Sat, 27 Aug 2022 11:24:26 -0700 Subject: [PATCH 33/39] simplified tet logic --- src/sdf/include/sdf.h | 51 +++++++++++++++++++++---------------------- 1 file changed, 25 insertions(+), 26 deletions(-) diff --git a/src/sdf/include/sdf.h b/src/sdf/include/sdf.h index bcb54ae04..c0d733045 100644 --- a/src/sdf/include/sdf.h +++ b/src/sdf/include/sdf.h @@ -324,38 +324,37 @@ struct BuildTris { bool skipTet = thisVert.key == kOpen; tet[2] = thisVert.Inside(); - int edges[6] = {base.edgeVerts[0], -1, -1, -1, -1, -1}; for (const int i : {0, 1, 2}) { - edges[1] = base.edgeVerts[i + 1]; - edges[4] = thisVert.edgeVerts[i + 4]; - edges[5] = base.edgeVerts[Prev3(i) + 4]; - thisIndex = leadIndex; thisIndex[Prev3(i)] -= 1; - thisVert = gridVerts[MortonCode(thisIndex)]; - - tet[3] = thisVert.Inside(); - edges[2] = thisVert.edgeVerts[Next3(i) + 4]; - edges[3] = thisVert.edgeVerts[Prev3(i) + 1]; - - if (!skipTet && thisVert.key != kOpen) CreateTris(tet, edges); - skipTet = thisVert.key == kOpen; - - tet[2] = tet[3]; - edges[1] = edges[5]; - edges[2] = thisVert.edgeVerts[i + 4]; - edges[4] = edges[3]; - edges[5] = base.edgeVerts[Next3(i) + 1]; + GridVert nextVert = gridVerts[MortonCode(thisIndex)]; + tet[3] = nextVert.Inside(); + + const int edges1[6] = {base.edgeVerts[0], + base.edgeVerts[i + 1], + nextVert.edgeVerts[Next3(i) + 4], + nextVert.edgeVerts[Prev3(i) + 1], + thisVert.edgeVerts[i + 4], + base.edgeVerts[Prev3(i) + 4]}; + thisVert = nextVert; + if (!skipTet && nextVert.key != kOpen) CreateTris(tet, edges1); + skipTet = nextVert.key == kOpen; thisIndex = baseIndex; thisIndex[Next3(i)] += 1; - thisVert = gridVerts[MortonCode(thisIndex)]; - - tet[3] = thisVert.Inside(); - edges[3] = thisVert.edgeVerts[Next3(i) + 4]; - - if (!skipTet && thisVert.key != kOpen) CreateTris(tet, edges); - skipTet = thisVert.key == kOpen; + nextVert = gridVerts[MortonCode(thisIndex)]; + tet[2] = tet[3]; + tet[3] = nextVert.Inside(); + + const int edges2[6] = {base.edgeVerts[0], + edges1[5], + thisVert.edgeVerts[i + 4], + nextVert.edgeVerts[Next3(i) + 4], + edges1[3], + base.edgeVerts[Next3(i) + 1]}; + thisVert = nextVert; + if (!skipTet && nextVert.key != kOpen) CreateTris(tet, edges2); + skipTet = nextVert.key == kOpen; tet[2] = tet[3]; } From 55de478d737174a4bbe960a2cd19b23d5f732627 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Sat, 27 Aug 2022 12:54:35 -0700 Subject: [PATCH 34/39] use NeighborInside() --- src/sdf/include/sdf.h | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/sdf/include/sdf.h b/src/sdf/include/sdf.h index c0d733045..c0a6bdeb1 100644 --- a/src/sdf/include/sdf.h +++ b/src/sdf/include/sdf.h @@ -147,6 +147,10 @@ struct GridVert { int edgeVerts[7] = {-1, -1, -1, -1, -1, -1, -1}; __host__ __device__ int Inside() const { return distance > 0 ? 1 : -1; } + + __host__ __device__ int NeighborInside(int i) const { + return Inside() * (edgeVerts[i] < 0 ? 1 : -1); + } }; class HashTableD { @@ -311,24 +315,21 @@ struct BuildTris { leadIndex.w = 0; } - const GridVert& leadVert = gridVerts[MortonCode(leadIndex)]; - if (leadVert.key == kOpen) return; - // This GridVert is in charge of the 6 tetrahedra surrounding its edge in - // the (1,1,1) direction, attached to leadVert. - glm::ivec4 tet(leadVert.Inside(), base.Inside(), -2, -2); + // the (1,1,1) direction (edge 0). + glm::ivec4 tet(base.NeighborInside(0), base.Inside(), -2, -2); glm::ivec4 thisIndex = baseIndex; thisIndex.x += 1; GridVert thisVert = gridVerts[MortonCode(thisIndex)]; bool skipTet = thisVert.key == kOpen; - tet[2] = thisVert.Inside(); + tet[2] = base.NeighborInside(1); for (const int i : {0, 1, 2}) { thisIndex = leadIndex; thisIndex[Prev3(i)] -= 1; GridVert nextVert = gridVerts[MortonCode(thisIndex)]; - tet[3] = nextVert.Inside(); + tet[3] = base.NeighborInside(Prev3(i) + 4); const int edges1[6] = {base.edgeVerts[0], base.edgeVerts[i + 1], @@ -344,7 +345,7 @@ struct BuildTris { thisIndex[Next3(i)] += 1; nextVert = gridVerts[MortonCode(thisIndex)]; tet[2] = tet[3]; - tet[3] = nextVert.Inside(); + tet[3] = base.NeighborInside(Next3(i) + 1); const int edges2[6] = {base.edgeVerts[0], edges1[5], From 1ba085e0ff043310d64e38b325c95833baca35b4 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Tue, 30 Aug 2022 10:22:47 -0700 Subject: [PATCH 35/39] optimized SDF --- src/sdf/include/sdf.h | 60 ++++++++++++++++++------------------------- test/sdf_test.cpp | 11 ++++---- 2 files changed, 31 insertions(+), 40 deletions(-) diff --git a/src/sdf/include/sdf.h b/src/sdf/include/sdf.h index c0a6bdeb1..c25084c0f 100644 --- a/src/sdf/include/sdf.h +++ b/src/sdf/include/sdf.h @@ -87,20 +87,13 @@ __host__ __device__ glm::ivec3 TetTri1(int i) { } __host__ __device__ glm::ivec4 Neighbors(int i) { - constexpr glm::ivec4 neighbors[14] = {{0, 0, 0, 1}, // - {1, 0, 0, 0}, // - {0, 1, 0, 0}, // - {0, 0, 1, 0}, // - {-1, 0, 0, 1}, // - {0, -1, 0, 1}, // - {0, 0, -1, 1}, // - {-1, -1, -1, 1}, // - {-1, 0, 0, 0}, // - {0, -1, 0, 0}, // - {0, 0, -1, 0}, // - {0, -1, -1, 1}, // - {-1, 0, -1, 1}, // - {-1, -1, 0, 1}}; + constexpr glm::ivec4 neighbors[7] = {{0, 0, 0, 1}, // + {1, 0, 0, 0}, // + {0, 1, 0, 0}, // + {0, 0, 1, 0}, // + {-1, 0, 0, 1}, // + {0, -1, 0, 1}, // + {0, 0, -1, 1}}; return neighbors[i]; } @@ -235,7 +228,7 @@ struct ComputeVerts { const float d = sdf(Position(gridIndex)) - level; const glm::ivec3 xyz(gridIndex); - const bool onLowerBound = glm::any(glm::equal(xyz, glm::ivec3(0))); + const bool onLowerBound = glm::any(glm::lessThanEqual(xyz, glm::ivec3(0))); const bool onUpperBound = glm::any(glm::greaterThanEqual(xyz, gridSize)); const bool onHalfBound = gridIndex.w == 1 && glm::any(glm::greaterThanEqual(xyz, gridSize - 1)); @@ -256,7 +249,9 @@ struct ComputeVerts { gridVert.distance = BoundedSDF(gridIndex); bool keep = false; - for (int i = 0; i < 14; ++i) { + // These seven edges are uniquely owned by this gridVert; any of them + // which intersect the surface create a vert. + for (int i = 0; i < 7; ++i) { glm::ivec4 neighborIndex = gridIndex + Neighbors(i); if (neighborIndex.w == 2) { neighborIndex += 1; @@ -266,15 +261,11 @@ struct ComputeVerts { if ((val > 0) == (gridVert.distance > 0)) continue; keep = true; - // These seven edges are uniquely owned by this gridVert; any of them - // which intersect the surface create a vert. - if (i < 7) { - const int idx = AtomicAdd(*vertIndex, 1); - vertPos[idx] = - (val * position - gridVert.distance * Position(neighborIndex)) / - (val - gridVert.distance); - gridVert.edgeVerts[i] = idx; - } + const int idx = AtomicAdd(*vertIndex, 1); + vertPos[idx] = + (val * position - gridVert.distance * Position(neighborIndex)) / + (val - gridVert.distance); + gridVert.edgeVerts[i] = idx; } if (keep && gridVerts.Insert(gridVert)) printf("out of space!\n"); @@ -322,13 +313,16 @@ struct BuildTris { thisIndex.x += 1; GridVert thisVert = gridVerts[MortonCode(thisIndex)]; - bool skipTet = thisVert.key == kOpen; tet[2] = base.NeighborInside(1); for (const int i : {0, 1, 2}) { thisIndex = leadIndex; - thisIndex[Prev3(i)] -= 1; - GridVert nextVert = gridVerts[MortonCode(thisIndex)]; + --thisIndex[Prev3(i)]; + // MortonCodes take unsigned input, so check for negatives, given the + // decrement. + GridVert nextVert = thisIndex[Prev3(i)] < 0 + ? GridVert() + : gridVerts[MortonCode(thisIndex)]; tet[3] = base.NeighborInside(Prev3(i) + 4); const int edges1[6] = {base.edgeVerts[0], @@ -338,11 +332,10 @@ struct BuildTris { thisVert.edgeVerts[i + 4], base.edgeVerts[Prev3(i) + 4]}; thisVert = nextVert; - if (!skipTet && nextVert.key != kOpen) CreateTris(tet, edges1); - skipTet = nextVert.key == kOpen; + CreateTris(tet, edges1); thisIndex = baseIndex; - thisIndex[Next3(i)] += 1; + ++thisIndex[Next3(i)]; nextVert = gridVerts[MortonCode(thisIndex)]; tet[2] = tet[3]; tet[3] = base.NeighborInside(Next3(i) + 1); @@ -354,8 +347,7 @@ struct BuildTris { edges1[3], base.edgeVerts[Next3(i) + 1]}; thisVert = nextVert; - if (!skipTet && nextVert.key != kOpen) CreateTris(tet, edges2); - skipTet = nextVert.key == kOpen; + CreateTris(tet, edges2); tet[2] = tet[3]; } @@ -428,8 +420,6 @@ inline Mesh LevelSet(Func sdf, Box bounds, float edgeLength, float level = 0) { BuildTris({triVerts.ptrD(), index.ptrD(), gridVerts.D()})); triVerts.resize(index[0]); - RemoveUnreferencedVerts(vertPos, triVerts); - out.vertPos.insert(out.vertPos.end(), vertPos.begin(), vertPos.end()); out.triVerts.insert(out.triVerts.end(), triVerts.begin(), triVerts.end()); return out; diff --git a/test/sdf_test.cpp b/test/sdf_test.cpp index 1a3ea31dd..4ca7801f6 100644 --- a/test/sdf_test.cpp +++ b/test/sdf_test.cpp @@ -47,13 +47,14 @@ TEST(SDF, Bounds) { const float size = 4; const float edgeLength = 0.5; - Manifold cubeVoid(LevelSet( - CubeVoid(), {glm::vec3(-size / 2), glm::vec3(size / 2)}, edgeLength)); + Mesh levelSet = LevelSet( + CubeVoid(), {glm::vec3(-size / 2), glm::vec3(size / 2)}, edgeLength); + Manifold cubeVoid(levelSet); Box bounds = cubeVoid.BoundingBox(); const float precision = cubeVoid.Precision(); - if (options.exportModels) ExportMesh("cubeVoid.gltf", cubeVoid.GetMesh(), {}); + if (options.exportModels) ExportMesh("cubeVoid.gltf", levelSet, {}); - EXPECT_TRUE(cubeVoid.IsManifold()); + EXPECT_EQ(cubeVoid.Status(), Manifold::Error::NO_ERROR); EXPECT_EQ(cubeVoid.Genus(), -1); const float outerBound = size / 2 + edgeLength / 2; EXPECT_NEAR(bounds.min.x, -outerBound, precision); @@ -77,7 +78,7 @@ TEST(SDF, Surface) { const float precision = cube.Precision(); if (options.exportModels) ExportMesh("cube.gltf", cube.GetMesh(), {}); - EXPECT_TRUE(cube.IsManifold()); + EXPECT_EQ(cubeVoid.Status(), Manifold::Error::NO_ERROR); EXPECT_EQ(cube.Genus(), 0); auto prop = cube.GetProperties(); EXPECT_NEAR(prop.volume, 8, 0.001); From 0aefceb1b02ac97e4dabb804432937ca19324ccd Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Tue, 30 Aug 2022 10:49:39 -0700 Subject: [PATCH 36/39] moved RemoveUnreferencedVerts --- src/manifold/src/impl.cpp | 53 ++++++++++++++++++++++++++---- src/manifold/src/impl.h | 1 + src/sdf/CMakeLists.txt | 19 ++++------- src/sdf/include/sdf.h | 3 -- src/sdf/src/sdf.cpp | 68 --------------------------------------- 5 files changed, 54 insertions(+), 90 deletions(-) delete mode 100644 src/sdf/src/sdf.cpp diff --git a/src/manifold/src/impl.cpp b/src/manifold/src/impl.cpp index f11b43ae2..e5272be2f 100644 --- a/src/manifold/src/impl.cpp +++ b/src/manifold/src/impl.cpp @@ -137,6 +137,26 @@ struct LinkHalfedges { } }; +struct MarkVerts { + int* vert; + + __host__ __device__ void operator()(glm::ivec3 triVerts) { + for (int i : {0, 1, 2}) { + vert[triVerts[i]] = 1; + } + } +}; + +struct ReindexTriVerts { + const int* old2new; + + __host__ __device__ void operator()(glm::ivec3& triVerts) { + for (int i : {0, 1, 2}) { + triVerts[i] = old2new[triVerts[i]]; + } + } +}; + struct InitializeBaryRef { const int meshID; const Halfedge* halfedge; @@ -283,6 +303,13 @@ Manifold::Impl::Impl(const Mesh& mesh, const std::vector& properties, const std::vector& propertyTolerance) : vertPos_(mesh.vertPos), halfedgeTangent_(mesh.halfedgeTangent) { + VecDH triVerts = mesh.triVerts; + if (!IsIndexInBounds(triVerts)) { + MarkFailure(Error::VERTEX_INDEX_OUT_OF_BOUNDS); + return; + } + RemoveUnreferencedVerts(triVerts); + CalculateBBox(); if (!IsFinite()) { MarkFailure(Error::NON_FINITE_VERTEX); @@ -290,12 +317,6 @@ Manifold::Impl::Impl(const Mesh& mesh, } SetPrecision(); - VecDH triVerts = mesh.triVerts; - if (!IsIndexInBounds(triVerts)) { - MarkFailure(Error::VERTEX_INDEX_OUT_OF_BOUNDS); - return; - } - CreateHalfedges(triVerts); if (!IsManifold()) { MarkFailure(Error::NOT_MANIFOLD); @@ -359,6 +380,26 @@ Manifold::Impl::Impl(Shape shape) { InitializeNewReference(); } +void Manifold::Impl::RemoveUnreferencedVerts(VecDH& triVerts) { + VecDH vertOld2New(NumVert() + 1, 0); + auto policy = autoPolicy(NumVert()); + for_each(policy, triVerts.cbegin(), triVerts.cend(), + MarkVerts({vertOld2New.ptrD() + 1})); + + const VecDH oldVertPos = vertPos_; + vertPos_.resize(copy_if( + policy, oldVertPos.cbegin(), oldVertPos.cend(), + vertOld2New.cbegin() + 1, vertPos_.begin(), + thrust::identity()) - + vertPos_.begin()); + + inclusive_scan(policy, vertOld2New.begin() + 1, vertOld2New.end(), + vertOld2New.begin() + 1); + + for_each(policy, triVerts.begin(), triVerts.end(), + ReindexTriVerts({vertOld2New.cptrD()})); +} + void Manifold::Impl::ReinitializeReference(int meshID) { // instead of storing the meshID, we store 0 and set the mapping to // 0 -> meshID, because the meshID after boolean operation also starts from 0. diff --git a/src/manifold/src/impl.h b/src/manifold/src/impl.h index bd3f89d79..421a98ef5 100644 --- a/src/manifold/src/impl.h +++ b/src/manifold/src/impl.h @@ -59,6 +59,7 @@ struct Manifold::Impl { const std::vector& properties = std::vector(), const std::vector& propertyTolerance = std::vector()); + void RemoveUnreferencedVerts(VecDH& triVerts); void ReinitializeReference(int meshID); void CreateHalfedges(const VecDH& triVerts); void CalculateNormals(); diff --git a/src/sdf/CMakeLists.txt b/src/sdf/CMakeLists.txt index 2038246bd..b4f3b5877 100644 --- a/src/sdf/CMakeLists.txt +++ b/src/sdf/CMakeLists.txt @@ -14,18 +14,11 @@ project(sdf) -file(GLOB_RECURSE SOURCE_FILES CONFIGURE_DEPENDS *.cpp) -add_library(${PROJECT_NAME} ${SOURCE_FILES}) +add_library(${PROJECT_NAME} INTERFACE) -if(MANIFOLD_USE_CUDA) - set_source_files_properties(${SOURCE_FILES} PROPERTIES LANGUAGE CUDA) - set_property(TARGET ${PROJECT_NAME} PROPERTY CUDA_ARCHITECTURES 61) -endif() - -target_include_directories(${PROJECT_NAME} PUBLIC ${PROJECT_SOURCE_DIR}/include) -target_link_libraries(${PROJECT_NAME} PUBLIC utilities) - -target_compile_features(${PROJECT_NAME} - PUBLIC cxx_std_14 - PRIVATE cxx_std_17 +target_include_directories(${PROJECT_NAME} + INTERFACE + ${PROJECT_SOURCE_DIR}/include ) + +target_link_libraries(${PROJECT_NAME} INTERFACE utilities) \ No newline at end of file diff --git a/src/sdf/include/sdf.h b/src/sdf/include/sdf.h index c25084c0f..588fbb6d2 100644 --- a/src/sdf/include/sdf.h +++ b/src/sdf/include/sdf.h @@ -357,9 +357,6 @@ struct BuildTris { namespace manifold { -void RemoveUnreferencedVerts(VecDH& vertPos, - VecDH& triVerts); - /** @addtogroup Core * @{ */ diff --git a/src/sdf/src/sdf.cpp b/src/sdf/src/sdf.cpp deleted file mode 100644 index c3b7cfbcc..000000000 --- a/src/sdf/src/sdf.cpp +++ /dev/null @@ -1,68 +0,0 @@ -// Copyright 2022 The Manifold Authors. -// -// Licensed under the Apache License, Version 2.0 (the "License"); -// you may not use this file except in compliance with the License. -// You may obtain a copy of the License at -// -// http://www.apache.org/licenses/LICENSE-2.0 -// -// Unless required by applicable law or agreed to in writing, software -// distributed under the License is distributed on an "AS IS" BASIS, -// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -// See the License for the specific language governing permissions and -// limitations under the License. - -#include "sdf.h" - -namespace { -using namespace manifold; - -struct MarkVerts { - int* vert; - - __host__ __device__ void operator()(glm::ivec3 triVerts) { - for (int i : {0, 1, 2}) { - vert[triVerts[i]] = 1; - } - } -}; - -struct ReindexVerts { - const int* old2new; - - __host__ __device__ void operator()(glm::ivec3& triVerts) { - for (int i : {0, 1, 2}) { - triVerts[i] = old2new[triVerts[i]]; - } - } -}; - -} // namespace - -namespace manifold { -/** - * Sorts the vertices according to their Morton code. - */ -void RemoveUnreferencedVerts(VecDH& vertPos, - VecDH& triVerts) { - const int numVert = vertPos.size(); - VecDH vertOld2New(numVert + 1, 0); - auto policy = autoPolicy(numVert); - for_each(policy, triVerts.cbegin(), triVerts.cend(), - MarkVerts({vertOld2New.ptrD() + 1})); - - const VecDH oldVertPos = vertPos; - vertPos.resize(copy_if( - policy, oldVertPos.cbegin(), oldVertPos.cend(), - vertOld2New.cbegin() + 1, vertPos.begin(), - thrust::identity()) - - vertPos.begin()); - - inclusive_scan(policy, vertOld2New.begin() + 1, vertOld2New.end(), - vertOld2New.begin() + 1); - - for_each(policy, triVerts.begin(), triVerts.end(), - ReindexVerts({vertOld2New.cptrD()})); -} - -} // namespace manifold \ No newline at end of file From 3ccbba958bfbf84b3a28f908077e884b80af852d Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Tue, 30 Aug 2022 11:13:34 -0700 Subject: [PATCH 37/39] put printf statements behind MANIFOLD_DEBUG --- src/manifold/src/boolean3.cpp | 23 +++++++++++++++++++++-- src/manifold/src/edge_op.cpp | 1 - src/manifold/src/impl.cpp | 5 +++-- src/manifold/src/manifold.cpp | 6 +++--- src/manifold/src/properties.cpp | 2 ++ 5 files changed, 29 insertions(+), 8 deletions(-) diff --git a/src/manifold/src/boolean3.cpp b/src/manifold/src/boolean3.cpp index d8eeed6a1..a60a495ae 100644 --- a/src/manifold/src/boolean3.cpp +++ b/src/manifold/src/boolean3.cpp @@ -30,7 +30,9 @@ namespace { __host__ __device__ glm::vec2 Interpolate(glm::vec3 pL, glm::vec3 pR, float x) { float dxL = x - pL.x; float dxR = x - pR.x; +#ifdef MANIFOLD_DEBUG if (dxL * dxR > 0) printf("Not in domain!\n"); +#endif bool useL = fabs(dxL) < fabs(dxR); float lambda = (useL ? dxL : dxR) / (pR.x - pL.x); if (!isfinite(lambda)) return glm::vec2(pL.y, pL.z); @@ -46,7 +48,9 @@ __host__ __device__ glm::vec4 Intersect(const glm::vec3 &pL, const glm::vec3 &qR) { float dyL = qL.y - pL.y; float dyR = qR.y - pR.y; +#ifdef MANIFOLD_DEBUG if (dyL * dyR > 0) printf("No intersection!\n"); +#endif bool useL = fabs(dyL) < fabs(dyR); float dx = pR.x - pL.x; float lambda = (useL ? dyL : dyR) / (dyL - dyR); @@ -102,6 +106,17 @@ __host__ __device__ bool Shadows(float p, float q, float dir) { return p == q ? dir < 0 : p < q; } +/** + * Since this function is called from two different places, it is necessary that + * it returns identical results for identical input to keep consistency. + * Normally this is trivial as computers make floating-point errors, but are + * at least deterministic. However, in the case of CUDA, these functions can be + * compiled by two different compilers (for the CPU and GPU). We have found that + * the different compilers can cause slightly different rounding errors, so it + * is critical that the two places this function is called both use the same + * compiled function (they must agree on CPU or GPU). This is now taken care of + * by the shared policy_ member. + */ __host__ __device__ thrust::pair Shadow01( const int p0, const int q1, const glm::vec3 *vertPosP, const glm::vec3 *vertPosQ, const Halfedge *halfedgeQ, const float expandP, @@ -218,11 +233,12 @@ struct Kernel11 { if (s11 == 0) { // No intersection xyzz11 = glm::vec4(NAN); } else { +#ifdef MANIFOLD_DEBUG // Assert left and right were both found if (k != 2) { printf("k = %d\n", k); } - +#endif xyzz11 = Intersect(pRL[0], pRL[1], qRL[0], qRL[1]); const int p1s = halfedgeP[p1].startVert; @@ -316,11 +332,12 @@ struct Kernel02 { if (s02 == 0) { // No intersection z02 = NAN; } else { +#ifdef MANIFOLD_DEBUG // Assert left and right were both found if (k != 2) { printf("k = %d\n", k); } - +#endif glm::vec3 vertPos = vertPosP[p0]; z02 = Interpolate(yzzRL[0], yzzRL[1], vertPos.y)[1]; if (forward) { @@ -433,10 +450,12 @@ struct Kernel12 { if (x12 == 0) { // No intersection v12 = glm::vec3(NAN); } else { +#ifdef MANIFOLD_DEBUG // Assert left and right were both found if (k != 2) { printf("k = %d\n", k); } +#endif const glm::vec4 xzyy = Intersect(xzyLR0[0], xzyLR0[1], xzyLR1[0], xzyLR1[1]); v12.x = xzyy[0]; diff --git a/src/manifold/src/edge_op.cpp b/src/manifold/src/edge_op.cpp index 548e6d6c5..e259f04eb 100644 --- a/src/manifold/src/edge_op.cpp +++ b/src/manifold/src/edge_op.cpp @@ -96,7 +96,6 @@ struct SwappableEdge { glm::vec2 v[3]; for (int i : {0, 1, 2}) v[i] = projection * vertPos[halfedge[triedge[i]].startVert]; - // if (CCW(v[0], v[1], v[2], precision) < 0) printf("tri %d is CW!\n", tri); if (CCW(v[0], v[1], v[2], precision) > 0 || !Is01Longest(v[0], v[1], v[2])) return false; diff --git a/src/manifold/src/impl.cpp b/src/manifold/src/impl.cpp index e5272be2f..334becc3c 100644 --- a/src/manifold/src/impl.cpp +++ b/src/manifold/src/impl.cpp @@ -295,8 +295,9 @@ namespace manifold { std::atomic Manifold::Impl::meshIDCounter_(1); /** - * Create a manifold from an input triangle Mesh. Will throw if the Mesh is not - * manifold. TODO: update halfedgeTangent during SimplifyTopology. + * Create a manifold from an input triangle Mesh. Will return an empty Manifold + * and set an Error Status if the Mesh is not manifold or otherwise invalid. + * TODO: update halfedgeTangent during SimplifyTopology. */ Manifold::Impl::Impl(const Mesh& mesh, const std::vector& triProperties, diff --git a/src/manifold/src/manifold.cpp b/src/manifold/src/manifold.cpp index d143f9ee8..4402216bd 100644 --- a/src/manifold/src/manifold.cpp +++ b/src/manifold/src/manifold.cpp @@ -82,9 +82,9 @@ CsgLeafNode& Manifold::GetCsgLeafNode() const { } /** - * Convert a Mesh into a Manifold. Will throw a topologyErr exception if the - * input is not an oriented 2-manifold. Will collapse degenerate triangles and - * unnecessary vertices. + * Convert a Mesh into a Manifold. Will return an empty Manifold + * and set an Error Status if the Mesh is not an oriented 2-manifold. Will + * collapse degenerate triangles and unnecessary vertices. * * The three optional inputs should all be specified if any are. These define * any properties you may have on this mesh. These properties are not saved in diff --git a/src/manifold/src/properties.cpp b/src/manifold/src/properties.cpp index 7cefc0bcb..3560bdcb8 100644 --- a/src/manifold/src/properties.cpp +++ b/src/manifold/src/properties.cpp @@ -205,6 +205,7 @@ struct CheckCCW { int ccw = CCW(v[0], v[1], v[2], glm::abs(tol)); bool check = tol > 0 ? ccw >= 0 : ccw == 0; +#ifdef MANIFOLD_DEBUG if (tol > 0 && !check) { glm::vec2 v1 = v[1] - v[0]; glm::vec2 v2 = v[2] - v[0]; @@ -225,6 +226,7 @@ struct CheckCCW { norm.y, norm.z, halfedges[3 * face].startVert, halfedges[3 * face + 1].startVert, halfedges[3 * face + 2].startVert); } +#endif return check; } }; From c56f1bcd271f78a9b83c0be1ad3a26ac8a77a9e3 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Tue, 30 Aug 2022 12:16:51 -0700 Subject: [PATCH 38/39] enable hashtable resize --- src/sdf/include/sdf.h | 62 ++++++++++++++++++++++------------ src/utilities/include/vec_dh.h | 2 +- 2 files changed, 41 insertions(+), 23 deletions(-) diff --git a/src/sdf/include/sdf.h b/src/sdf/include/sdf.h index 588fbb6d2..7efecefd5 100644 --- a/src/sdf/include/sdf.h +++ b/src/sdf/include/sdf.h @@ -153,19 +153,19 @@ class HashTableD { __host__ __device__ int Size() const { return table_.size(); } - __host__ __device__ bool Insert(const GridVert& vert) { + __host__ __device__ bool Full() const { return used_[0] * 2 > Size(); } + + __host__ __device__ void Insert(const GridVert& vert) { uint32_t idx = vert.key & (Size() - 1); while (1) { Uint64& key = table_[idx].key; const Uint64 found = AtomicCAS(key, kOpen, vert.key); if (found == kOpen) { - if (AtomicAdd(used_[0], 0x1u) * 2 > Size()) { - return true; - } + AtomicAdd(used_[0], 0x1u); table_[idx] = vert; - return false; + return; } - if (found == vert.key) return false; + if (found == vert.key) return; idx = (idx + step_) & (Size() - 1); } } @@ -182,7 +182,7 @@ class HashTableD { __host__ __device__ GridVert At(int idx) const { return table_[idx]; } private: - const uint32_t step_; + uint32_t step_; VecD table_; VecD used_; }; @@ -198,9 +198,9 @@ class HashTable { int Size() const { return table_.Size(); } - float FilledFraction() const { - return static_cast(used_[0]) / table_.Size(); - } + bool Full() const { return used_[0] * 2 > Size(); } + + float FilledFraction() const { return static_cast(used_[0]) / Size(); } private: VecDH alloc_; @@ -238,6 +238,8 @@ struct ComputeVerts { } inline __host__ __device__ void operator()(Uint64 mortonCode) { + if (gridVerts.Full()) return; + const glm::ivec4 gridIndex = DecodeMorton(mortonCode); if (glm::any(glm::greaterThan(glm::ivec3(gridIndex), gridSize))) return; @@ -268,7 +270,7 @@ struct ComputeVerts { gridVert.edgeVerts[i] = idx; } - if (keep && gridVerts.Insert(gridVert)) printf("out of space!\n"); + if (keep) gridVerts.Insert(gridVert); } }; @@ -394,25 +396,41 @@ inline Mesh LevelSet(Func sdf, Box bounds, float edgeLength, float level = 0) { const glm::ivec3 gridSize(dim / edgeLength); const glm::vec3 spacing = dim / (glm::vec3(gridSize)); - const int maxMorton = MortonCode(glm::ivec4(gridSize + 1, 1)); + const Uint64 maxMorton = MortonCode(glm::ivec4(gridSize + 1, 1)); const auto policy = autoPolicy(maxMorton); - const int tableSize = glm::min( - 2 * maxMorton, static_cast(100 * glm::pow(maxMorton, 0.667))); + int tableSize = glm::min( + 2 * maxMorton, static_cast(100 * glm::pow(maxMorton, 0.667))); HashTable gridVerts(tableSize); - VecDH vertPos(gridVerts.Size() * 7); - VecDH index(1, 0); - for_each_n( - policy, countAt(0), maxMorton + 1, - ComputeVerts({vertPos.ptrD(), index.ptrD(), gridVerts.D(), sdf, - bounds.min, gridSize + 1, spacing, level})); - vertPos.resize(index[0]); + while (1) { + VecDH index(1, 0); + for_each_n( + policy, countAt(0), maxMorton + 1, + ComputeVerts({vertPos.ptrD(), index.ptrD(), gridVerts.D(), sdf, + bounds.min, gridSize + 1, spacing, level})); + + if (gridVerts.Full()) { // Resize HashTable + const glm::vec3 lastVert = vertPos[index[0] - 1]; + const Uint64 lastMorton = + MortonCode(glm::ivec4((lastVert - bounds.min) / spacing, 1)); + const float ratio = static_cast(maxMorton) / lastMorton; + if (ratio > 1000) // do not trust the ratio if it is too large + tableSize *= 2; + else + tableSize *= 2 * ratio; + gridVerts = HashTable(tableSize); + vertPos = VecDH(gridVerts.Size() * 7); + } else { // Success + vertPos.resize(index[0]); + break; + } + } VecDH triVerts(gridVerts.Entries() * 12); // worst case - index[0] = 0; + VecDH index(1, 0); for_each_n(policy, countAt(0), gridVerts.Size(), BuildTris({triVerts.ptrD(), index.ptrD(), gridVerts.D()})); triVerts.resize(index[0]); diff --git a/src/utilities/include/vec_dh.h b/src/utilities/include/vec_dh.h index a1ed004cc..cc181af89 100644 --- a/src/utilities/include/vec_dh.h +++ b/src/utilities/include/vec_dh.h @@ -439,7 +439,7 @@ class VecD { private: T *ptr_; - const int size_; + int size_; }; /** @} */ } // namespace manifold From 7bdaa27298256edbd101dd0cfb4a10218e3ec4dd Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Tue, 30 Aug 2022 12:59:01 -0700 Subject: [PATCH 39/39] added test for hashtable resizing --- src/sdf/include/sdf.h | 4 ++-- test/sdf_test.cpp | 16 ++++++++++++++++ 2 files changed, 18 insertions(+), 2 deletions(-) diff --git a/src/sdf/include/sdf.h b/src/sdf/include/sdf.h index 7efecefd5..fc2bd2a19 100644 --- a/src/sdf/include/sdf.h +++ b/src/sdf/include/sdf.h @@ -400,7 +400,7 @@ inline Mesh LevelSet(Func sdf, Box bounds, float edgeLength, float level = 0) { const auto policy = autoPolicy(maxMorton); int tableSize = glm::min( - 2 * maxMorton, static_cast(100 * glm::pow(maxMorton, 0.667))); + 2 * maxMorton, static_cast(10 * glm::pow(maxMorton, 0.667))); HashTable gridVerts(tableSize); VecDH vertPos(gridVerts.Size() * 7); @@ -419,7 +419,7 @@ inline Mesh LevelSet(Func sdf, Box bounds, float edgeLength, float level = 0) { if (ratio > 1000) // do not trust the ratio if it is too large tableSize *= 2; else - tableSize *= 2 * ratio; + tableSize *= ratio; gridVerts = HashTable(tableSize); vertPos = VecDH(gridVerts.Size() * 7); } else { // Success diff --git a/test/sdf_test.cpp b/test/sdf_test.cpp index 4ca7801f6..b0c289103 100644 --- a/test/sdf_test.cpp +++ b/test/sdf_test.cpp @@ -30,6 +30,13 @@ struct CubeVoid { } }; +struct Layers { + __host__ __device__ float operator()(glm::vec3 p) const { + int a = glm::mod(glm::round(2 * p.z), 4.0f); + return a == 0 ? 1 : (a == 2 ? -1 : 0); + } +}; + TEST(SDF, CubeVoid) { CubeVoid voidSDF; @@ -89,4 +96,13 @@ TEST(SDF, Surface) { EXPECT_NEAR(bounds.max.x, 1, precision); EXPECT_NEAR(bounds.max.y, 1, precision); EXPECT_NEAR(bounds.max.z, 1, precision); +} + +TEST(SDF, Resize) { + const float size = 20; + Manifold layers(LevelSet(Layers(), {glm::vec3(0), glm::vec3(size)}, 1)); + if (options.exportModels) ExportMesh("layers.gltf", layers.GetMesh(), {}); + + EXPECT_EQ(layers.Status(), Manifold::Error::NO_ERROR); + EXPECT_EQ(layers.Genus(), -8); } \ No newline at end of file