768 lines
25 KiB
C++
768 lines
25 KiB
C++
// 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 <map>
|
|
#include <optional>
|
|
|
|
#include "./hashtable.h"
|
|
#include "./mesh_fixes.h"
|
|
#include "./parallel.h"
|
|
#include "./svd.h"
|
|
|
|
#ifdef MANIFOLD_EXPORT
|
|
#include <string.h>
|
|
|
|
#include <iostream>
|
|
#endif
|
|
|
|
namespace {
|
|
using namespace manifold;
|
|
|
|
constexpr uint64_t kRemove = std::numeric_limits<uint64_t>::max();
|
|
|
|
void AtomicAddVec3(vec3& target, const vec3& add) {
|
|
for (int i : {0, 1, 2}) {
|
|
std::atomic<double>& tar =
|
|
reinterpret_cast<std::atomic<double>&>(target[i]);
|
|
double old_val = tar.load(std::memory_order_relaxed);
|
|
while (!tar.compare_exchange_weak(old_val, old_val + add[i],
|
|
std::memory_order_relaxed)) {
|
|
}
|
|
}
|
|
}
|
|
|
|
struct Transform4x3 {
|
|
const mat3x4 transform;
|
|
|
|
vec3 operator()(vec3 position) { return transform * vec4(position, 1.0); }
|
|
};
|
|
|
|
template <bool calculateTriNormal>
|
|
struct AssignNormals {
|
|
VecView<vec3> faceNormal;
|
|
VecView<vec3> vertNormal;
|
|
VecView<const vec3> vertPos;
|
|
VecView<const Halfedge> halfedges;
|
|
|
|
void operator()(const int face) {
|
|
vec3& triNormal = faceNormal[face];
|
|
|
|
ivec3 triVerts;
|
|
for (int i : {0, 1, 2}) triVerts[i] = halfedges[3 * face + i].startVert;
|
|
|
|
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]]);
|
|
}
|
|
|
|
if (calculateTriNormal) {
|
|
triNormal = la::normalize(la::cross(edge[0], edge[1]));
|
|
if (std::isnan(triNormal.x)) triNormal = vec3(0, 0, 1);
|
|
}
|
|
|
|
// corner angles
|
|
vec3 phi;
|
|
double dot = -la::dot(edge[2], edge[0]);
|
|
phi[0] = dot >= 1 ? 0 : (dot <= -1 ? kPi : std::acos(dot));
|
|
dot = -la::dot(edge[0], edge[1]);
|
|
phi[1] = dot >= 1 ? 0 : (dot <= -1 ? kPi : std::acos(dot));
|
|
phi[2] = kPi - phi[0] - phi[1];
|
|
|
|
// assign weighted sum
|
|
for (int i : {0, 1, 2}) {
|
|
AtomicAddVec3(vertNormal[triVerts[i]], phi[i] * triNormal);
|
|
}
|
|
}
|
|
};
|
|
|
|
struct UpdateMeshID {
|
|
const HashTableD<uint32_t> meshIDold2new;
|
|
|
|
void operator()(TriRef& ref) { ref.meshID = meshIDold2new[ref.meshID]; }
|
|
};
|
|
|
|
struct CoplanarEdge {
|
|
VecView<std::pair<int, int>> face2face;
|
|
VecView<double> triArea;
|
|
VecView<const Halfedge> halfedge;
|
|
VecView<const vec3> vertPos;
|
|
VecView<const TriRef> triRef;
|
|
VecView<const ivec3> triProp;
|
|
const int numProp;
|
|
const double epsilon;
|
|
const double tolerance;
|
|
|
|
void operator()(const int edgeIdx) {
|
|
const Halfedge edge = halfedge[edgeIdx];
|
|
const Halfedge pair = halfedge[edge.pairedHalfedge];
|
|
const int edgeFace = edgeIdx / 3;
|
|
const int pairFace = edge.pairedHalfedge / 3;
|
|
|
|
if (triRef[edgeFace].meshID != triRef[pairFace].meshID) return;
|
|
|
|
const vec3 base = vertPos[edge.startVert];
|
|
const int baseNum = edgeIdx - 3 * edgeFace;
|
|
const int jointNum = edge.pairedHalfedge - 3 * pairFace;
|
|
|
|
if (numProp > 0) {
|
|
if (triProp[edgeFace][baseNum] != triProp[pairFace][Next3(jointNum)] ||
|
|
triProp[edgeFace][Next3(baseNum)] != triProp[pairFace][jointNum])
|
|
return;
|
|
}
|
|
|
|
if (!edge.IsForward()) return;
|
|
|
|
const int edgeNum = baseNum == 0 ? 2 : baseNum - 1;
|
|
const int pairNum = jointNum == 0 ? 2 : jointNum - 1;
|
|
const vec3 jointVec = vertPos[pair.startVert] - base;
|
|
const vec3 edgeVec =
|
|
vertPos[halfedge[3 * edgeFace + edgeNum].startVert] - base;
|
|
const vec3 pairVec =
|
|
vertPos[halfedge[3 * pairFace + pairNum].startVert] - base;
|
|
|
|
const double length = std::max(la::length(jointVec), la::length(edgeVec));
|
|
const double lengthPair =
|
|
std::max(la::length(jointVec), la::length(pairVec));
|
|
vec3 normal = la::cross(jointVec, edgeVec);
|
|
const double area = la::length(normal);
|
|
const double areaPair = la::length(la::cross(pairVec, jointVec));
|
|
|
|
// make sure we only write this once
|
|
if (edgeIdx % 3 == 0) triArea[edgeFace] = area;
|
|
// Don't link degenerate triangles
|
|
if (area < length * epsilon || areaPair < lengthPair * epsilon) return;
|
|
|
|
const double volume = std::abs(la::dot(normal, pairVec));
|
|
// Only operate on coplanar triangles
|
|
if (volume > std::max(area, areaPair) * tolerance) return;
|
|
|
|
face2face[edgeIdx] = std::make_pair(edgeFace, pairFace);
|
|
}
|
|
};
|
|
|
|
struct CheckCoplanarity {
|
|
VecView<int> comp2tri;
|
|
VecView<const Halfedge> halfedge;
|
|
VecView<const vec3> vertPos;
|
|
std::vector<int>* components;
|
|
const double tolerance;
|
|
|
|
void operator()(int tri) {
|
|
const int component = (*components)[tri];
|
|
const int referenceTri =
|
|
reinterpret_cast<std::atomic<int>*>(&comp2tri[component])
|
|
->load(std::memory_order_relaxed);
|
|
if (referenceTri < 0 || referenceTri == tri) return;
|
|
|
|
const vec3 origin = vertPos[halfedge[3 * referenceTri].startVert];
|
|
const vec3 normal = la::normalize(
|
|
la::cross(vertPos[halfedge[3 * referenceTri + 1].startVert] - origin,
|
|
vertPos[halfedge[3 * referenceTri + 2].startVert] - origin));
|
|
|
|
for (const int i : {0, 1, 2}) {
|
|
const vec3 vert = vertPos[halfedge[3 * tri + i].startVert];
|
|
// If any component vertex is not coplanar with the component's reference
|
|
// triangle, unmark the entire component so that none of its triangles are
|
|
// marked coplanar.
|
|
if (std::abs(la::dot(normal, vert - origin)) > tolerance) {
|
|
reinterpret_cast<std::atomic<int>*>(&comp2tri[component])
|
|
->store(-1, std::memory_order_relaxed);
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
};
|
|
|
|
int GetLabels(std::vector<int>& components,
|
|
const Vec<std::pair<int, int>>& edges, int numNodes) {
|
|
UnionFind<> uf(numNodes);
|
|
for (auto edge : edges) {
|
|
if (edge.first == -1 || edge.second == -1) continue;
|
|
uf.unionXY(edge.first, edge.second);
|
|
}
|
|
|
|
return uf.connectedComponents(components);
|
|
}
|
|
|
|
void DedupePropVerts(manifold::Vec<ivec3>& triProp,
|
|
const Vec<std::pair<int, int>>& vert2vert,
|
|
size_t numPropVert) {
|
|
ZoneScoped;
|
|
std::vector<int> vertLabels;
|
|
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 (auto& prop : triProp)
|
|
for (int i : {0, 1, 2}) prop[i] = label2vert[vertLabels[prop[i]]];
|
|
}
|
|
} // namespace
|
|
|
|
namespace manifold {
|
|
|
|
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();
|
|
CreateFaces();
|
|
}
|
|
|
|
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(NumTri());
|
|
for_each_n(autoPolicy(NumTri(), 1e5), countAt(0), NumTri(),
|
|
[meshID, keepFaceID, &triRef](const int tri) {
|
|
triRef[tri] = {meshID, meshID, tri,
|
|
keepFaceID ? triRef[tri].faceID : tri};
|
|
});
|
|
meshRelation_.meshIDtransform.clear();
|
|
meshRelation_.meshIDtransform[meshID] = {meshID};
|
|
}
|
|
|
|
void Manifold::Impl::CreateFaces() {
|
|
ZoneScoped;
|
|
Vec<std::pair<int, int>> face2face(halfedge_.size(), {-1, -1});
|
|
Vec<std::pair<int, int>> vert2vert(halfedge_.size(), {-1, -1});
|
|
Vec<double> triArea(NumTri());
|
|
|
|
const size_t numProp = NumProp();
|
|
if (numProp > 0) {
|
|
for_each_n(
|
|
autoPolicy(halfedge_.size(), 1e4), countAt(0), halfedge_.size(),
|
|
[&vert2vert, numProp, this](const int edgeIdx) {
|
|
const Halfedge edge = halfedge_[edgeIdx];
|
|
const Halfedge pair = halfedge_[edge.pairedHalfedge];
|
|
const int edgeFace = edgeIdx / 3;
|
|
const int pairFace = edge.pairedHalfedge / 3;
|
|
|
|
if (meshRelation_.triRef[edgeFace].meshID !=
|
|
meshRelation_.triRef[pairFace].meshID)
|
|
return;
|
|
|
|
const int baseNum = edgeIdx - 3 * edgeFace;
|
|
const int jointNum = edge.pairedHalfedge - 3 * pairFace;
|
|
|
|
const int prop0 = meshRelation_.triProperties[edgeFace][baseNum];
|
|
const int prop1 =
|
|
meshRelation_
|
|
.triProperties[pairFace][jointNum == 2 ? 0 : jointNum + 1];
|
|
if (prop0 == prop1) return;
|
|
|
|
bool propEqual = true;
|
|
for (size_t p = 0; p < numProp; ++p) {
|
|
if (meshRelation_.properties[numProp * prop0 + p] !=
|
|
meshRelation_.properties[numProp * prop1 + p]) {
|
|
propEqual = false;
|
|
break;
|
|
}
|
|
}
|
|
if (propEqual) {
|
|
vert2vert[edgeIdx] = std::make_pair(prop0, prop1);
|
|
}
|
|
});
|
|
DedupePropVerts(meshRelation_.triProperties, vert2vert, NumPropVert());
|
|
}
|
|
|
|
for_each_n(autoPolicy(halfedge_.size(), 1e4), countAt(0), halfedge_.size(),
|
|
CoplanarEdge({face2face, triArea, halfedge_, vertPos_,
|
|
meshRelation_.triRef, meshRelation_.triProperties,
|
|
meshRelation_.numProp, epsilon_, tolerance_}));
|
|
|
|
std::vector<int> components;
|
|
const int numComponent = GetLabels(components, face2face, NumTri());
|
|
|
|
Vec<int> comp2tri(numComponent, -1);
|
|
for (size_t tri = 0; tri < NumTri(); ++tri) {
|
|
const int comp = components[tri];
|
|
const int current = comp2tri[comp];
|
|
if (current < 0 || triArea[tri] > triArea[current]) {
|
|
comp2tri[comp] = tri;
|
|
triArea[comp] = triArea[tri];
|
|
}
|
|
}
|
|
|
|
for_each_n(autoPolicy(halfedge_.size(), 1e4), countAt(0), NumTri(),
|
|
CheckCoplanarity(
|
|
{comp2tri, halfedge_, vertPos_, &components, tolerance_}));
|
|
|
|
Vec<TriRef>& triRef = meshRelation_.triRef;
|
|
for (size_t tri = 0; tri < NumTri(); ++tri) {
|
|
const int referenceTri = comp2tri[components[tri]];
|
|
if (referenceTri >= 0) {
|
|
triRef[tri].faceID = referenceTri;
|
|
}
|
|
}
|
|
}
|
|
|
|
/**
|
|
* Create the halfedge_ data structure from an input triVerts array like Mesh.
|
|
*/
|
|
void Manifold::Impl::CreateHalfedges(const Vec<ivec3>& triVerts) {
|
|
ZoneScoped;
|
|
const size_t numTri = triVerts.size();
|
|
const int numHalfedge = 3 * numTri;
|
|
// drop the old value first to avoid copy
|
|
halfedge_.resize(0);
|
|
halfedge_.resize(numHalfedge);
|
|
Vec<uint64_t> edge(numHalfedge);
|
|
Vec<int> ids(numHalfedge);
|
|
auto policy = autoPolicy(numTri, 1e5);
|
|
sequence(ids.begin(), ids.end());
|
|
for_each_n(policy, countAt(0), numTri,
|
|
[this, &edge, &triVerts](const int tri) {
|
|
const ivec3& verts = triVerts[tri];
|
|
for (const int i : {0, 1, 2}) {
|
|
const int j = (i + 1) % 3;
|
|
const int e = 3 * tri + i;
|
|
halfedge_[e] = {verts[i], verts[j], -1};
|
|
// Sort the forward halfedges in front of the backward ones
|
|
// by setting the highest-order bit.
|
|
edge[e] = uint64_t(verts[i] < verts[j] ? 1 : 0) << 63 |
|
|
((uint64_t)std::min(verts[i], verts[j])) << 32 |
|
|
std::max(verts[i], verts[j]);
|
|
}
|
|
});
|
|
// Stable sort is required here so that halfedges from the same face are
|
|
// paired together (the triangles were created in face order). In some
|
|
// degenerate situations the triangulator can add the same internal edge in
|
|
// two different faces, causing this edge to not be 2-manifold. These are
|
|
// fixed by duplicating verts in SimplifyTopology.
|
|
stable_sort(ids.begin(), ids.end(), [&edge](const int& a, const int& b) {
|
|
return edge[a] < edge[b];
|
|
});
|
|
|
|
// Mark opposed triangles for removal - this may strand unreferenced verts
|
|
// which are removed later by RemoveUnreferencedVerts() and Finish().
|
|
const int numEdge = numHalfedge / 2;
|
|
for (int i = 0; i < numEdge; ++i) {
|
|
const int pair0 = ids[i];
|
|
Halfedge h0 = halfedge_[pair0];
|
|
int k = i + numEdge;
|
|
while (1) {
|
|
const int pair1 = ids[k];
|
|
Halfedge h1 = halfedge_[pair1];
|
|
if (h0.startVert != h1.endVert || h0.endVert != h1.startVert) break;
|
|
if (halfedge_[NextHalfedge(pair0)].endVert ==
|
|
halfedge_[NextHalfedge(pair1)].endVert) {
|
|
h0 = {-1, -1, -1};
|
|
h1 = {-1, -1, -1};
|
|
// Reorder so that remaining edges pair up
|
|
if (k != i + numEdge) std::swap(ids[i + numEdge], ids[k]);
|
|
break;
|
|
}
|
|
++k;
|
|
if (k >= numHalfedge) break;
|
|
}
|
|
}
|
|
|
|
// Once sorted, the first half of the range is the forward halfedges, which
|
|
// correspond to their backward pair at the same offset in the second half
|
|
// of the range.
|
|
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].startVert >= 0) {
|
|
halfedge_[pair0].pairedHalfedge = pair1;
|
|
halfedge_[pair1].pairedHalfedge = pair0;
|
|
}
|
|
});
|
|
}
|
|
|
|
/**
|
|
* 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::MarkFailure(Error status) {
|
|
bBox_ = Box();
|
|
vertPos_.resize(0);
|
|
halfedge_.resize(0);
|
|
vertNormal_.resize(0);
|
|
faceNormal_.resize(0);
|
|
halfedgeTangent_.resize(0);
|
|
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()) {
|
|
MarkFailure(Error::NonFiniteVertex);
|
|
return;
|
|
}
|
|
Update();
|
|
faceNormal_.resize(0); // force recalculation of triNormal
|
|
CalculateNormals();
|
|
SetEpsilon();
|
|
Finish();
|
|
CreateFaces();
|
|
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.MarkFailure(Error::NonFiniteVertex);
|
|
return result;
|
|
}
|
|
result.collider_ = collider_;
|
|
result.meshRelation_ = meshRelation_;
|
|
result.epsilon_ = epsilon_;
|
|
result.tolerance_ = tolerance_;
|
|
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(), 1e4);
|
|
fill(vertNormal_.begin(), vertNormal_.end(), vec3(0.0));
|
|
bool calculateTriNormal = false;
|
|
if (faceNormal_.size() != NumTri()) {
|
|
faceNormal_.resize(NumTri());
|
|
calculateTriNormal = true;
|
|
}
|
|
if (calculateTriNormal)
|
|
for_each_n(
|
|
policy, countAt(0), NumTri(),
|
|
AssignNormals<true>({faceNormal_, vertNormal_, vertPos_, halfedge_}));
|
|
else
|
|
for_each_n(
|
|
policy, countAt(0), NumTri(),
|
|
AssignNormals<false>({faceNormal_, vertNormal_, vertPos_, halfedge_}));
|
|
for_each(policy, vertNormal_.begin(), vertNormal_.end(),
|
|
[](vec3& v) { v = SafeNormalize(v); });
|
|
}
|
|
|
|
/**
|
|
* 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()}));
|
|
}
|
|
|
|
/**
|
|
* Returns a sparse array of the bounding box overlaps between the edges of
|
|
* the input manifold, Q and the faces of this manifold. Returned indices only
|
|
* point to forward halfedges.
|
|
*/
|
|
SparseIndices Manifold::Impl::EdgeCollisions(const Impl& Q,
|
|
bool inverted) const {
|
|
ZoneScoped;
|
|
Vec<TmpEdge> edges = CreateTmpEdges(Q.halfedge_);
|
|
const size_t numEdge = edges.size();
|
|
Vec<Box> QedgeBB(numEdge);
|
|
const auto& vertPos = Q.vertPos_;
|
|
auto policy = autoPolicy(numEdge, 1e5);
|
|
for_each_n(
|
|
policy, countAt(0), numEdge, [&QedgeBB, &edges, &vertPos](const int e) {
|
|
QedgeBB[e] = Box(vertPos[edges[e].first], vertPos[edges[e].second]);
|
|
});
|
|
|
|
SparseIndices q1p2(0);
|
|
if (inverted)
|
|
q1p2 = collider_.Collisions<false, true>(QedgeBB.cview());
|
|
else
|
|
q1p2 = collider_.Collisions<false, false>(QedgeBB.cview());
|
|
|
|
if (inverted)
|
|
for_each(policy, countAt(0_uz), countAt(q1p2.size()),
|
|
ReindexEdge<true>({edges, q1p2}));
|
|
else
|
|
for_each(policy, countAt(0_uz), countAt(q1p2.size()),
|
|
ReindexEdge<false>({edges, q1p2}));
|
|
return q1p2;
|
|
}
|
|
|
|
/**
|
|
* Returns a sparse array of the input vertices that project inside the XY
|
|
* bounding boxes of the faces of this manifold.
|
|
*/
|
|
SparseIndices Manifold::Impl::VertexCollisionsZ(VecView<const vec3> vertsIn,
|
|
bool inverted) const {
|
|
ZoneScoped;
|
|
if (inverted)
|
|
return collider_.Collisions<false, true>(vertsIn);
|
|
else
|
|
return collider_.Collisions<false, false>(vertsIn);
|
|
}
|
|
|
|
#ifdef MANIFOLD_DEBUG
|
|
std::ostream& operator<<(std::ostream& stream, const Manifold::Impl& impl) {
|
|
stream << std::setprecision(17); // for double precision
|
|
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;
|
|
}
|
|
#endif
|
|
|
|
#ifdef MANIFOLD_EXPORT
|
|
Manifold Manifold::ImportMeshGL64(std::istream& stream) {
|
|
MeshGL64 mesh;
|
|
std::optional<double> epsilon;
|
|
stream.precision(17);
|
|
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 i : {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 (tmp[end] != 0 && end < SIZE) 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 i : {0, 1, 2}) {
|
|
double x;
|
|
stream >> x;
|
|
mesh.vertProperties.push_back(x);
|
|
}
|
|
break;
|
|
case 'f':
|
|
for (int i : {0, 1, 2}) {
|
|
uint64_t x;
|
|
stream >> x;
|
|
mesh.triVerts.push_back(x - 1);
|
|
}
|
|
break;
|
|
case '\n':
|
|
break;
|
|
default:
|
|
DEBUG_ASSERT(false, userErr, "unexpected character in MeshGL64 import");
|
|
}
|
|
}
|
|
auto m = std::make_shared<Manifold::Impl>(mesh);
|
|
if (epsilon) m->SetEpsilon(*epsilon);
|
|
return Manifold(m);
|
|
}
|
|
#endif
|
|
|
|
} // namespace manifold
|