// Copyright 2021 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 "impl.h"
#include <algorithm>
#include <atomic>
#include <cstring>
#include <map>
#include <optional>
#include "csg_tree.h"
#include "disjoint_sets.h"
#include "hashtable.h"
#include "manifold/optional_assert.h"
#include "mesh_fixes.h"
#include "parallel.h"
#include "shared.h"
#include "svd.h"
namespace {
using namespace manifold;
/**
* Returns arc cosine of 𝑥.
*
* @return value in range [0,M_PI]
* @return NAN if 𝑥 ∈ {NAN,+INFINITY,-INFINITY}
* @return NAN if 𝑥 ∉ [-1,1]
*/
double sun_acos(double x) {
/*
* Origin of acos function: FreeBSD /usr/src/lib/msun/src/e_acos.c
* Changed the use of union to memcpy to avoid undefined behavior.
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
constexpr double pio2_hi =
1.57079632679489655800e+00; /* 0x3FF921FB, 0x54442D18 */
constexpr double pio2_lo =
6.12323399573676603587e-17; /* 0x3C91A626, 0x33145C07 */
constexpr double pS0 =
1.66666666666666657415e-01; /* 0x3FC55555, 0x55555555 */
constexpr double pS1 =
-3.25565818622400915405e-01; /* 0xBFD4D612, 0x03EB6F7D */
constexpr double pS2 =
2.01212532134862925881e-01; /* 0x3FC9C155, 0x0E884455 */
constexpr double pS3 =
-4.00555345006794114027e-02; /* 0xBFA48228, 0xB5688F3B */
constexpr double pS4 =
7.91534994289814532176e-04; /* 0x3F49EFE0, 0x7501B288 */
constexpr double pS5 =
3.47933107596021167570e-05; /* 0x3F023DE1, 0x0DFDF709 */
constexpr double qS1 =
-2.40339491173441421878e+00; /* 0xC0033A27, 0x1C8A2D4B */
constexpr double qS2 =
2.02094576023350569471e+00; /* 0x40002AE5, 0x9C598AC8 */
constexpr double qS3 =
-6.88283971605453293030e-01; /* 0xBFE6066C, 0x1B8D0159 */
constexpr double qS4 =
7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
auto R = [=](double z) {
double p, q;
p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5)))));
q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4)));
return p / q;
};
double z, w, s, c, df;
uint64_t xx;
uint32_t hx, lx, ix;
memcpy(&xx, &x, sizeof(xx));
hx = xx >> 32;
ix = hx & 0x7fffffff;
/* |x| >= 1 or nan */
if (ix >= 0x3ff00000) {
lx = xx;
if (((ix - 0x3ff00000) | lx) == 0) {
/* acos(1)=0, acos(-1)=pi */
if (hx >> 31) return 2 * pio2_hi + 0x1p-120f;
return 0;
}
return 0 / (x - x);
}
/* |x| < 0.5 */
if (ix < 0x3fe00000) {
if (ix <= 0x3c600000) /* |x| < 2**-57 */
return pio2_hi + 0x1p-120f;
return pio2_hi - (x - (pio2_lo - x * R(x * x)));
}
/* x < -0.5 */
if (hx >> 31) {
z = (1.0 + x) * 0.5;
s = sqrt(z);
w = R(z) * s - pio2_lo;
return 2 * (pio2_hi - (s + w));
}
/* x > 0.5 */
z = (1.0 - x) * 0.5;
s = sqrt(z);
memcpy(&xx, &s, sizeof(xx));
xx &= 0xffffffff00000000;
memcpy(&df, &xx, sizeof(xx));
c = (z - df * df) / (s + df);
w = R(z) * s + c;
return 2 * (df + w);
}
struct Transform4x3 {
const mat3x4 transform;
vec3 operator()(vec3 position) { return transform * vec4(position, 1.0); }
};
struct UpdateMeshID {
const HashTableD<uint32_t> meshIDold2new;
void operator()(TriRef& ref) { ref.meshID = meshIDold2new[ref.meshID]; }
};
int GetLabels(std::vector<int>& components,
const Vec<std::pair<int, int>>& edges, int numNodes) {
DisjointSets uf(numNodes);
for (auto edge : edges) {
if (edge.first == -1 || edge.second == -1) continue;
uf.unite(edge.first, edge.second);
}
return uf.connectedComponents(components);
}
} // namespace
namespace manifold {
#if (MANIFOLD_PAR == 1)
#if (TBB_VERSION_MAJOR < 2021)
tbb::task_arena gc_arena(1, 1);
#else
tbb::task_arena gc_arena(1, 1, tbb::task_arena::priority::low);
#endif
#endif
std::atomic<uint32_t> Manifold::Impl::meshIDCounter_(1);
uint32_t Manifold::Impl::ReserveIDs(uint32_t n) {
return Manifold::Impl::meshIDCounter_.fetch_add(n, std::memory_order_relaxed);
}
/**
* Create either a unit tetrahedron, cube or octahedron. The cube is in the
* first octant, while the others are symmetric about the origin.
*/
Manifold::Impl::Impl(Shape shape, const mat3x4 m) {
std::vector<vec3> vertPos;
std::vector<ivec3> triVerts;
switch (shape) {
case Shape::Tetrahedron:
vertPos = {{-1.0, -1.0, 1.0},
{-1.0, 1.0, -1.0},
{1.0, -1.0, -1.0},
{1.0, 1.0, 1.0}};
triVerts = {{2, 0, 1}, {0, 3, 1}, {2, 3, 0}, {3, 2, 1}};
break;
case Shape::Cube:
vertPos = {{0.0, 0.0, 0.0}, //
{0.0, 0.0, 1.0}, //
{0.0, 1.0, 0.0}, //
{0.0, 1.0, 1.0}, //
{1.0, 0.0, 0.0}, //
{1.0, 0.0, 1.0}, //
{1.0, 1.0, 0.0}, //
{1.0, 1.0, 1.0}};
triVerts = {{1, 0, 4}, {2, 4, 0}, //
{1, 3, 0}, {3, 1, 5}, //
{3, 2, 0}, {3, 7, 2}, //
{5, 4, 6}, {5, 1, 4}, //
{6, 4, 2}, {7, 6, 2}, //
{7, 3, 5}, {7, 5, 6}};
break;
case Shape::Octahedron:
vertPos = {{1.0, 0.0, 0.0}, //
{-1.0, 0.0, 0.0}, //
{0.0, 1.0, 0.0}, //
{0.0, -1.0, 0.0}, //
{0.0, 0.0, 1.0}, //
{0.0, 0.0, -1.0}};
triVerts = {{0, 2, 4}, {1, 5, 3}, //
{2, 1, 4}, {3, 5, 0}, //
{1, 3, 4}, {0, 5, 2}, //
{3, 0, 4}, {2, 5, 1}};
break;
}
vertPos_ = vertPos;
for (auto& v : vertPos_) v = m * vec4(v, 1.0);
CreateHalfedges(triVerts);
Finish();
InitializeOriginal();
MarkCoplanar();
}
void Manifold::Impl::RemoveUnreferencedVerts() {
ZoneScoped;
const int numVert = NumVert();
Vec<int> keep(numVert, 0);
auto policy = autoPolicy(numVert, 1e5);
for_each(policy, halfedge_.cbegin(), halfedge_.cend(), [&keep](Halfedge h) {
if (h.startVert >= 0) {
reinterpret_cast<std::atomic<int>*>(&keep[h.startVert])
->store(1, std::memory_order_relaxed);
}
});
for_each_n(policy, countAt(0), numVert, [&keep, this](int v) {
if (keep[v] == 0) {
vertPos_[v] = vec3(NAN);
}
});
}
void Manifold::Impl::InitializeOriginal(bool keepFaceID) {
const int meshID = ReserveIDs(1);
meshRelation_.originalID = meshID;
auto& triRef = meshRelation_.triRef;
triRef.resize_nofill(NumTri());
for_each_n(autoPolicy(NumTri(), 1e5), countAt(0), NumTri(),
[meshID, keepFaceID, &triRef](const int tri) {
triRef[tri] = {meshID, meshID, -1,
keepFaceID ? triRef[tri].coplanarID : tri};
});
meshRelation_.meshIDtransform.clear();
meshRelation_.meshIDtransform[meshID] = {meshID};
}
void Manifold::Impl::MarkCoplanar() {
ZoneScoped;
const int numTri = NumTri();
struct TriPriority {
double area2;
int tri;
};
Vec<TriPriority> triPriority(numTri);
for_each_n(autoPolicy(numTri), countAt(0), numTri,
[&triPriority, this](int tri) {
meshRelation_.triRef[tri].coplanarID = -1;
if (halfedge_[3 * tri].startVert < 0) {
triPriority[tri] = {0, tri};
return;
}
const vec3 v = vertPos_[halfedge_[3 * tri].startVert];
triPriority[tri] = {
length2(cross(vertPos_[halfedge_[3 * tri].endVert] - v,
vertPos_[halfedge_[3 * tri + 1].endVert] - v)),
tri};
});
stable_sort(triPriority.begin(), triPriority.end(),
[](auto a, auto b) { return a.area2 > b.area2; });
Vec<int> interiorHalfedges;
for (const auto tp : triPriority) {
if (meshRelation_.triRef[tp.tri].coplanarID >= 0) continue;
meshRelation_.triRef[tp.tri].coplanarID = tp.tri;
if (halfedge_[3 * tp.tri].startVert < 0) continue;
const vec3 base = vertPos_[halfedge_[3 * tp.tri].startVert];
const vec3 normal = faceNormal_[tp.tri];
interiorHalfedges.resize(3);
interiorHalfedges[0] = 3 * tp.tri;
interiorHalfedges[1] = 3 * tp.tri + 1;
interiorHalfedges[2] = 3 * tp.tri + 2;
while (!interiorHalfedges.empty()) {
const int h =
NextHalfedge(halfedge_[interiorHalfedges.back()].pairedHalfedge);
interiorHalfedges.pop_back();
if (meshRelation_.triRef[h / 3].coplanarID >= 0) continue;
const vec3 v = vertPos_[halfedge_[h].endVert];
if (std::abs(dot(v - base, normal)) < tolerance_) {
meshRelation_.triRef[h / 3].coplanarID = tp.tri;
if (interiorHalfedges.empty() ||
h != halfedge_[interiorHalfedges.back()].pairedHalfedge) {
interiorHalfedges.push_back(h);
} else {
interiorHalfedges.pop_back();
}
const int hNext = NextHalfedge(h);
interiorHalfedges.push_back(hNext);
}
}
}
}
/**
* Dereference duplicate property vertices if they are exactly floating-point
* equal. These unreferenced properties are then removed by CompactProps.
*/
void Manifold::Impl::DedupePropVerts() {
ZoneScoped;
const size_t numProp = NumProp();
if (numProp == 0) return;
Vec<std::pair<int, int>> vert2vert(halfedge_.size(), {-1, -1});
for_each_n(autoPolicy(halfedge_.size(), 1e4), countAt(0), halfedge_.size(),
[&vert2vert, numProp, this](const int edgeIdx) {
const Halfedge edge = halfedge_[edgeIdx];
if (edge.pairedHalfedge < 0) return;
const int edgeFace = edgeIdx / 3;
const int pairFace = edge.pairedHalfedge / 3;
if (meshRelation_.triRef[edgeFace].meshID !=
meshRelation_.triRef[pairFace].meshID)
return;
const int prop0 = halfedge_[edgeIdx].propVert;
const int prop1 =
halfedge_[NextHalfedge(edge.pairedHalfedge)].propVert;
bool propEqual = true;
for (size_t p = 0; p < numProp; ++p) {
if (properties_[numProp * prop0 + p] !=
properties_[numProp * prop1 + p]) {
propEqual = false;
break;
}
}
if (propEqual) {
vert2vert[edgeIdx] = std::make_pair(prop0, prop1);
}
});
std::vector<int> vertLabels;
const size_t numPropVert = NumPropVert();
const int numLabels = GetLabels(vertLabels, vert2vert, numPropVert);
std::vector<int> label2vert(numLabels);
for (size_t v = 0; v < numPropVert; ++v) label2vert[vertLabels[v]] = v;
for (Halfedge& edge : halfedge_)
edge.propVert = label2vert[vertLabels[edge.propVert]];
}
constexpr int kRemovedHalfedge = -2;
struct HalfedgePairData {
int largeVert;
int tri;
int edgeIndex;
bool operator<(const HalfedgePairData& other) const {
return largeVert < other.largeVert ||
(largeVert == other.largeVert && tri < other.tri);
}
};
template <bool useProp, typename F>
struct PrepHalfedges {
VecView<Halfedge> halfedges;
const VecView<ivec3> triProp;
const VecView<ivec3> triVert;
F& f;
void operator()(const int tri) {
const ivec3& props = triProp[tri];
for (const int i : {0, 1, 2}) {
const int j = Next3(i);
const int k = Next3(j);
const int e = 3 * tri + i;
const int v0 = useProp ? props[i] : triVert[tri][i];
const int v1 = useProp ? props[j] : triVert[tri][j];
DEBUG_ASSERT(v0 != v1, logicErr, "topological degeneracy");
halfedges[e] = {v0, v1, -1, props[i]};
f(e, v0, v1);
}
}
};
/**
* Create the halfedge_ data structure from a list of triangles. If the optional
* prop2vert array is missing, it's assumed these triangles are are pointing to
* both vert and propVert indices. If prop2vert is present, the triangles are
* assumed to be pointing to propVert indices only. The prop2vert array is used
* to map the propVert indices to vert indices.
*/
void Manifold::Impl::CreateHalfedges(const Vec<ivec3>& triProp,
const Vec<ivec3>& triVert) {
ZoneScoped;
const size_t numTri = triProp.size();
const int numHalfedge = 3 * numTri;
// drop the old value first to avoid copy
halfedge_.clear(true);
halfedge_.resize_nofill(numHalfedge);
auto policy = autoPolicy(numTri, 1e5);
int vertCount = static_cast<int>(vertPos_.size());
Vec<int> ids(numHalfedge);
{
ZoneScopedN("PrepHalfedges");
if (vertCount < (1 << 18)) {
// For small vertex count, it is faster to just do sorting
Vec<uint64_t> edge(numHalfedge);
auto setEdge = [&edge](int e, int v0, int v1) {
edge[e] = static_cast<uint64_t>(v0 < v1 ? 1 : 0) << 63 |
(static_cast<uint64_t>(std::min(v0, v1))) << 32 |
static_cast<uint64_t>(std::max(v0, v1));
};
if (triVert.empty()) {
for_each_n(policy, countAt(0), numTri,
PrepHalfedges<true, decltype(setEdge)>{halfedge_, triProp,
triVert, setEdge});
} else {
for_each_n(policy, countAt(0), numTri,
PrepHalfedges<false, decltype(setEdge)>{halfedge_, triProp,
triVert, setEdge});
}
sequence(ids.begin(), ids.end());
stable_sort(ids.begin(), ids.end(), [&edge](const int& a, const int& b) {
return edge[a] < edge[b];
});
} else {
// For larger vertex count, we separate the ids into slices for halfedges
// with the same smaller vertex.
// We first copy them there (as HalfedgePairData), and then do sorting
// locally for each slice.
// This helps with memory locality, and is faster for larger meshes.
Vec<HalfedgePairData> entries(numHalfedge);
Vec<int> offsets(vertCount * 2, 0);
auto setOffset = [&offsets, vertCount](int _e, int v0, int v1) {
const int offset = v0 > v1 ? 0 : vertCount;
AtomicAdd(offsets[std::min(v0, v1) + offset], 1);
};
if (triVert.empty()) {
for_each_n(policy, countAt(0), numTri,
PrepHalfedges<true, decltype(setOffset)>{
halfedge_, triProp, triVert, setOffset});
} else {
for_each_n(policy, countAt(0), numTri,
PrepHalfedges<false, decltype(setOffset)>{
halfedge_, triProp, triVert, setOffset});
}
exclusive_scan(offsets.begin(), offsets.end(), offsets.begin());
for_each_n(policy, countAt(0), numTri,
[this, &offsets, &entries, vertCount](const int tri) {
for (const int i : {0, 1, 2}) {
const int e = 3 * tri + i;
const int v0 = halfedge_[e].startVert;
const int v1 = halfedge_[e].endVert;
const int offset = v0 > v1 ? 0 : vertCount;
const int start = std::min(v0, v1);
const int index = AtomicAdd(offsets[start + offset], 1);
entries[index] = {std::max(v0, v1), tri, e};
}
});
for_each_n(policy, countAt(0), offsets.size(), [&](const int v) {
int start = v == 0 ? 0 : offsets[v - 1];
int end = offsets[v];
for (int i = start; i < end; ++i) ids[i] = i;
std::sort(ids.begin() + start, ids.begin() + end,
[&entries](int a, int b) { return entries[a] < entries[b]; });
for (int i = start; i < end; ++i) ids[i] = entries[ids[i]].edgeIndex;
});
}
}
// Mark opposed triangles for removal - this may strand unreferenced verts
// which are removed later by RemoveUnreferencedVerts() and Finish().
const int numEdge = numHalfedge / 2;
const auto body = [&](int i, int consecutiveStart, int segmentEnd) {
const int pair0 = ids[i];
Halfedge& h0 = halfedge_[pair0];
int k = consecutiveStart + numEdge;
while (1) {
const int pair1 = ids[k];
Halfedge& h1 = halfedge_[pair1];
if (h0.startVert != h1.endVert || h0.endVert != h1.startVert) break;
if (h1.pairedHalfedge != kRemovedHalfedge &&
halfedge_[NextHalfedge(pair0)].endVert ==
halfedge_[NextHalfedge(pair1)].endVert) {
h0.pairedHalfedge = h1.pairedHalfedge = kRemovedHalfedge;
// Reorder so that remaining edges pair up
if (k != i + numEdge) std::swap(ids[i + numEdge], ids[k]);
break;
}
++k;
if (k >= segmentEnd + numEdge) break;
}
if (i + 1 == segmentEnd) return consecutiveStart;
Halfedge& h1 = halfedge_[ids[i + 1]];
if (h1.startVert == h0.startVert && h1.endVert == h0.endVert)
return consecutiveStart;
return i + 1;
};
#if MANIFOLD_PAR == 1
Vec<std::pair<int, int>> ranges;
const int increment = std::min(
std::max(numEdge / tbb::this_task_arena::max_concurrency() / 2, 1024),
numEdge);
const auto duplicated = [&](int a, int b) {
const Halfedge& h0 = halfedge_[ids[a]];
const Halfedge& h1 = halfedge_[ids[b]];
return h0.startVert == h1.startVert && h0.endVert == h1.endVert;
};
int end = 0;
while (end < numEdge) {
const int start = end;
end = std::min(end + increment, numEdge);
// make sure duplicated halfedges are in the same partition
while (end < numEdge && duplicated(end - 1, end)) end++;
ranges.push_back(std::make_pair(start, end));
}
for_each(ExecutionPolicy::Par, ranges.begin(), ranges.end(),
[&](const std::pair<int, int>& range) {
const auto [start, end] = range;
int consecutiveStart = start;
for (int i = start; i < end; ++i)
consecutiveStart = body(i, consecutiveStart, end);
});
#else
int consecutiveStart = 0;
for (int i = 0; i < numEdge; ++i)
consecutiveStart = body(i, consecutiveStart, numEdge);
#endif
for_each_n(policy, countAt(0), numEdge, [this, &ids, numEdge](int i) {
const int pair0 = ids[i];
const int pair1 = ids[i + numEdge];
if (halfedge_[pair0].pairedHalfedge != kRemovedHalfedge) {
halfedge_[pair0].pairedHalfedge = pair1;
halfedge_[pair1].pairedHalfedge = pair0;
} else {
halfedge_[pair0] = halfedge_[pair1] = {-1, -1, -1};
}
});
}
/**
* Does a full recalculation of the face bounding boxes, including updating
* the collider, but does not resort the faces.
*/
void Manifold::Impl::Update() {
CalculateBBox();
Vec<Box> faceBox;
Vec<uint32_t> faceMorton;
GetFaceBoxMorton(faceBox, faceMorton);
collider_.UpdateBoxes(faceBox);
}
void Manifold::Impl::MakeEmpty(Error status) {
bBox_ = Box();
vertPos_.clear();
halfedge_.clear();
vertNormal_.clear();
faceNormal_.clear();
halfedgeTangent_.clear();
meshRelation_ = MeshRelationD();
status_ = status;
}
void Manifold::Impl::Warp(std::function<void(vec3&)> warpFunc) {
WarpBatch([&warpFunc](VecView<vec3> vecs) {
for_each(ExecutionPolicy::Seq, vecs.begin(), vecs.end(), warpFunc);
});
}
void Manifold::Impl::WarpBatch(std::function<void(VecView<vec3>)> warpFunc) {
warpFunc(vertPos_.view());
CalculateBBox();
if (!IsFinite()) {
MakeEmpty(Error::NonFiniteVertex);
return;
}
Update();
faceNormal_.clear(); // force recalculation of triNormal
SetEpsilon();
Finish();
MarkCoplanar();
meshRelation_.originalID = -1;
}
Manifold::Impl Manifold::Impl::Transform(const mat3x4& transform_) const {
ZoneScoped;
if (transform_ == mat3x4(la::identity)) return *this;
auto policy = autoPolicy(NumVert());
Impl result;
if (status_ != Manifold::Error::NoError) {
result.status_ = status_;
return result;
}
if (!all(la::isfinite(transform_))) {
result.MakeEmpty(Error::NonFiniteVertex);
return result;
}
result.collider_ = collider_;
result.meshRelation_ = meshRelation_;
result.epsilon_ = epsilon_;
result.tolerance_ = tolerance_;
result.numProp_ = numProp_;
result.properties_ = properties_;
result.bBox_ = bBox_;
result.halfedge_ = halfedge_;
result.halfedgeTangent_.resize(halfedgeTangent_.size());
result.meshRelation_.originalID = -1;
for (auto& m : result.meshRelation_.meshIDtransform) {
m.second.transform = transform_ * Mat4(m.second.transform);
}
result.vertPos_.resize(NumVert());
result.faceNormal_.resize(faceNormal_.size());
result.vertNormal_.resize(vertNormal_.size());
transform(vertPos_.begin(), vertPos_.end(), result.vertPos_.begin(),
Transform4x3({transform_}));
mat3 normalTransform = NormalTransform(transform_);
transform(faceNormal_.begin(), faceNormal_.end(), result.faceNormal_.begin(),
TransformNormals({normalTransform}));
transform(vertNormal_.begin(), vertNormal_.end(), result.vertNormal_.begin(),
TransformNormals({normalTransform}));
const bool invert = la::determinant(mat3(transform_)) < 0;
if (halfedgeTangent_.size() > 0) {
for_each_n(policy, countAt(0), halfedgeTangent_.size(),
TransformTangents({result.halfedgeTangent_, 0, mat3(transform_),
invert, halfedgeTangent_, halfedge_}));
}
if (invert) {
for_each_n(policy, countAt(0), result.NumTri(),
FlipTris({result.halfedge_}));
}
// This optimization does a cheap collider update if the transform is
// axis-aligned.
if (!result.collider_.Transform(transform_)) result.Update();
result.CalculateBBox();
// Scale epsilon by the norm of the 3x3 portion of the transform.
result.epsilon_ *= SpectralNorm(mat3(transform_));
// Maximum of inherited epsilon loss and translational epsilon loss.
result.SetEpsilon(result.epsilon_);
return result;
}
/**
* Sets epsilon based on the bounding box, and limits its minimum value
* by the optional input.
*/
void Manifold::Impl::SetEpsilon(double minEpsilon, bool useSingle) {
epsilon_ = MaxEpsilon(minEpsilon, bBox_);
double minTol = epsilon_;
if (useSingle)
minTol =
std::max(minTol, std::numeric_limits<float>::epsilon() * bBox_.Scale());
tolerance_ = std::max(tolerance_, minTol);
}
/**
* If face normals are already present, this function uses them to compute
* vertex normals (angle-weighted pseudo-normals); otherwise it also computes
* the face normals. Face normals are only calculated when needed because
* nearly degenerate faces will accrue rounding error, while the Boolean can
* retain their original normal, which is more accurate and can help with
* merging coplanar faces.
*
* If the face normals have been invalidated by an operation like Warp(),
* ensure you do faceNormal_.resize(0) before calling this function to force
* recalculation.
*/
void Manifold::Impl::CalculateNormals() {
ZoneScoped;
vertNormal_.resize(NumVert());
auto policy = autoPolicy(NumTri());
std::vector<std::atomic<int>> vertHalfedgeMap(NumVert());
for_each_n(policy, countAt(0), NumVert(), [&](const size_t vert) {
vertHalfedgeMap[vert] = std::numeric_limits<int>::max();
});
auto atomicMin = [&vertHalfedgeMap](int value, int vert) {
if (vert < 0) return;
int old = std::numeric_limits<int>::max();
while (!vertHalfedgeMap[vert].compare_exchange_strong(old, value))
if (old < value) break;
};
if (faceNormal_.size() != NumTri()) {
faceNormal_.resize(NumTri());
for_each_n(policy, countAt(0), NumTri(), [&](const int face) {
vec3& triNormal = faceNormal_[face];
if (halfedge_[3 * face].startVert < 0) {
triNormal = vec3(0, 0, 1);
return;
}
ivec3 triVerts;
for (int i : {0, 1, 2}) {
int v = halfedge_[3 * face + i].startVert;
triVerts[i] = v;
atomicMin(3 * face + i, v);
}
vec3 edge[3];
for (int i : {0, 1, 2}) {
const int j = (i + 1) % 3;
edge[i] = la::normalize(vertPos_[triVerts[j]] - vertPos_[triVerts[i]]);
}
triNormal = la::normalize(la::cross(edge[0], edge[1]));
if (std::isnan(triNormal.x)) triNormal = vec3(0, 0, 1);
});
} else {
for_each_n(policy, countAt(0), halfedge_.size(),
[&](const int i) { atomicMin(i, halfedge_[i].startVert); });
}
for_each_n(policy, countAt(0), NumVert(), [&](const size_t vert) {
int firstEdge = vertHalfedgeMap[vert].load();
// not referenced
if (firstEdge == std::numeric_limits<int>::max()) {
vertNormal_[vert] = vec3(0.0);
return;
}
vec3 normal = vec3(0.0);
ForVert(firstEdge, [&](int edge) {
ivec3 triVerts = {halfedge_[edge].startVert, halfedge_[edge].endVert,
halfedge_[NextHalfedge(edge)].endVert};
vec3 currEdge =
la::normalize(vertPos_[triVerts[1]] - vertPos_[triVerts[0]]);
vec3 prevEdge =
la::normalize(vertPos_[triVerts[0]] - vertPos_[triVerts[2]]);
// if it is not finite, this means that the triangle is degenerate, and we
// should just exclude it from the normal calculation...
if (!la::isfinite(currEdge[0]) || !la::isfinite(prevEdge[0])) return;
double dot = -la::dot(prevEdge, currEdge);
double phi = dot >= 1 ? 0 : (dot <= -1 ? kPi : sun_acos(dot));
normal += phi * faceNormal_[edge / 3];
});
vertNormal_[vert] = SafeNormalize(normal);
});
}
/**
* Remaps all the contained meshIDs to new unique values to represent new
* instances of these meshes.
*/
void Manifold::Impl::IncrementMeshIDs() {
HashTable<uint32_t> meshIDold2new(meshRelation_.meshIDtransform.size() * 2);
// Update keys of the transform map
std::map<int, Relation> oldTransforms;
std::swap(meshRelation_.meshIDtransform, oldTransforms);
const int numMeshIDs = oldTransforms.size();
int nextMeshID = ReserveIDs(numMeshIDs);
for (const auto& pair : oldTransforms) {
meshIDold2new.D().Insert(pair.first, nextMeshID);
meshRelation_.meshIDtransform[nextMeshID++] = pair.second;
}
const size_t numTri = NumTri();
for_each_n(autoPolicy(numTri, 1e5), meshRelation_.triRef.begin(), numTri,
UpdateMeshID({meshIDold2new.D()}));
}
#ifdef MANIFOLD_DEBUG
/**
* Debugging output using high precision OBJ files with specialized comments
*/
std::ostream& operator<<(std::ostream& stream, const Manifold::Impl& impl) {
stream << std::setprecision(19); // for double precision
stream << std::fixed; // for uniformity in output numbers
stream << "# ======= begin mesh ======" << std::endl;
stream << "# tolerance = " << impl.tolerance_ << std::endl;
stream << "# epsilon = " << impl.epsilon_ << std::endl;
// TODO: Mesh relation, vertex normal and face normal
for (const vec3& v : impl.vertPos_)
stream << "v " << v.x << " " << v.y << " " << v.z << std::endl;
std::vector<ivec3> triangles;
triangles.reserve(impl.halfedge_.size() / 3);
for (size_t i = 0; i < impl.halfedge_.size(); i += 3)
triangles.emplace_back(impl.halfedge_[i].startVert + 1,
impl.halfedge_[i + 1].startVert + 1,
impl.halfedge_[i + 2].startVert + 1);
sort(triangles.begin(), triangles.end());
for (const auto& tri : triangles)
stream << "f " << tri.x << " " << tri.y << " " << tri.z << std::endl;
stream << "# ======== end mesh =======" << std::endl;
return stream;
}
/**
* Import a mesh from a Wavefront OBJ file that was exported with Write. This
* function is the counterpart to Write and should be used with it. This
* function is not guaranteed to be able to import OBJ files not written by the
* Write function.
*/
Manifold Manifold::ReadOBJ(std::istream& stream) {
if (!stream.good()) return Invalid();
MeshGL64 mesh;
std::optional<double> epsilon;
stream >> std::setprecision(19);
while (true) {
char c = stream.get();
if (stream.eof()) break;
switch (c) {
case '#': {
char c = stream.get();
if (c == ' ') {
constexpr int SIZE = 10;
std::array<char, SIZE> tmp;
stream.get(tmp.data(), SIZE, '\n');
if (strncmp(tmp.data(), "tolerance", SIZE) == 0) {
// skip 3 letters
for (int _ : {0, 1, 2}) stream.get();
stream >> mesh.tolerance;
} else if (strncmp(tmp.data(), "epsilon =", SIZE) == 0) {
double tmp;
stream >> tmp;
epsilon = {tmp};
} else {
// add it back because it is not what we want
int end = 0;
while (end < SIZE && tmp[end] != 0) end++;
while (--end > -1) stream.putback(tmp[end]);
}
c = stream.get();
}
// just skip the remaining comment
while (c != '\n' && !stream.eof()) {
c = stream.get();
}
break;
}
case 'v':
for (int _ : {0, 1, 2}) {
double x;
stream >> x;
mesh.vertProperties.push_back(x);
}
break;
case 'f':
for (int _ : {0, 1, 2}) {
uint64_t x;
stream >> x;
mesh.triVerts.push_back(x - 1);
}
break;
case '\r':
case '\n':
break;
default:
DEBUG_ASSERT(false, userErr, "unexpected character in Manifold import");
}
}
auto m = std::make_shared<Manifold::Impl>(mesh);
if (epsilon) m->SetEpsilon(*epsilon);
return Manifold(m);
}
/**
* Export the mesh to a Wavefront OBJ file in a way that preserves the full
* 64-bit precision of the vertex positions, as well as storing metadata such as
* the tolerance and epsilon. Useful for debugging and testing. Files written
* by WriteOBJ should be read back in with ReadOBJ.
*/
bool Manifold::WriteOBJ(std::ostream& stream) const {
if (!stream.good()) return false;
stream << *this->GetCsgLeafNode().GetImpl();
return true;
}
#endif
} // namespace manifold