// 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 <atomic>
#include <set>
#include "disjoint_sets.h"
#include "impl.h"
#include "parallel.h"
#include "shared.h"
namespace {
using namespace manifold;
constexpr uint32_t kNoCode = 0xFFFFFFFFu;
uint32_t MortonCode(vec3 position, Box bBox) {
// Unreferenced vertices are marked NaN, and this will sort them to the end
// (the Morton code only uses the first 30 of 32 bits).
if (std::isnan(position.x)) return kNoCode;
return Collider::MortonCode(position, bBox);
}
struct ReindexFace {
VecView<Halfedge> halfedge;
VecView<vec4> halfedgeTangent;
VecView<const Halfedge> oldHalfedge;
VecView<const vec4> oldHalfedgeTangent;
VecView<const int> faceNew2Old;
VecView<const int> faceOld2New;
void operator()(int newFace) {
const int oldFace = faceNew2Old[newFace];
for (const int i : {0, 1, 2}) {
const int oldEdge = 3 * oldFace + i;
Halfedge edge = oldHalfedge[oldEdge];
const int pairedFace = edge.pairedHalfedge / 3;
const int offset = edge.pairedHalfedge - 3 * pairedFace;
edge.pairedHalfedge = 3 * faceOld2New[pairedFace] + offset;
const int newEdge = 3 * newFace + i;
halfedge[newEdge] = edge;
if (!oldHalfedgeTangent.empty()) {
halfedgeTangent[newEdge] = oldHalfedgeTangent[oldEdge];
}
}
}
};
template <typename Precision, typename I>
bool MergeMeshGLP(MeshGLP<Precision, I>& mesh) {
ZoneScoped;
std::multiset<std::pair<int, int>> openEdges;
std::vector<int> merge(mesh.NumVert());
std::iota(merge.begin(), merge.end(), 0);
for (size_t i = 0; i < mesh.mergeFromVert.size(); ++i) {
merge[mesh.mergeFromVert[i]] = mesh.mergeToVert[i];
}
const auto numVert = mesh.NumVert();
const auto numTri = mesh.NumTri();
const int next[3] = {1, 2, 0};
for (size_t tri = 0; tri < numTri; ++tri) {
for (int i : {0, 1, 2}) {
auto edge = std::make_pair(merge[mesh.triVerts[3 * tri + next[i]]],
merge[mesh.triVerts[3 * tri + i]]);
auto it = openEdges.find(edge);
if (it == openEdges.end()) {
std::swap(edge.first, edge.second);
openEdges.insert(edge);
} else {
openEdges.erase(it);
}
}
}
if (openEdges.empty()) {
return false;
}
const auto numOpenVert = openEdges.size();
Vec<int> openVerts(numOpenVert);
int i = 0;
for (const auto& edge : openEdges) {
const int vert = edge.first;
openVerts[i++] = vert;
}
Vec<Precision> vertPropD(mesh.vertProperties);
Box bBox;
for (const int i : {0, 1, 2}) {
auto iPos =
StridedRange(vertPropD.begin() + i, vertPropD.end(), mesh.numProp);
auto minMax = manifold::transform_reduce(
iPos.begin(), iPos.end(),
std::make_pair(std::numeric_limits<double>::infinity(),
-std::numeric_limits<double>::infinity()),
[](auto a, auto b) {
return std::make_pair(std::min(a.first, b.first),
std::max(a.second, b.second));
},
[](double f) { return std::make_pair(f, f); });
bBox.min[i] = minMax.first;
bBox.max[i] = minMax.second;
}
const double tolerance = std::max(static_cast<double>(mesh.tolerance),
(std::is_same<Precision, float>::value
? std::numeric_limits<float>::epsilon()
: kPrecision) *
bBox.Scale());
auto policy = autoPolicy(numOpenVert, 1e5);
Vec<Box> vertBox(numOpenVert);
Vec<uint32_t> vertMorton(numOpenVert);
for_each_n(policy, countAt(0), numOpenVert,
[&vertMorton, &vertBox, &openVerts, &bBox, &mesh,
tolerance](const int i) {
int vert = openVerts[i];
const vec3 center(mesh.vertProperties[mesh.numProp * vert],
mesh.vertProperties[mesh.numProp * vert + 1],
mesh.vertProperties[mesh.numProp * vert + 2]);
vertBox[i].min = center - tolerance / 2.0;
vertBox[i].max = center + tolerance / 2.0;
vertMorton[i] = MortonCode(center, bBox);
});
Vec<int> vertNew2Old(numOpenVert);
sequence(vertNew2Old.begin(), vertNew2Old.end());
stable_sort(vertNew2Old.begin(), vertNew2Old.end(),
[&vertMorton](const int& a, const int& b) {
return vertMorton[a] < vertMorton[b];
});
Permute(vertMorton, vertNew2Old);
Permute(vertBox, vertNew2Old);
Permute(openVerts, vertNew2Old);
Collider collider(vertBox, vertMorton);
DisjointSets uf(numVert);
auto f = [&uf, &openVerts](int a, int b) {
return uf.unite(openVerts[a], openVerts[b]);
};
auto recorder = MakeSimpleRecorder(f);
collider.Collisions<true>(vertBox.cview(), recorder, false);
for (size_t i = 0; i < mesh.mergeFromVert.size(); ++i) {
uf.unite(static_cast<int>(mesh.mergeFromVert[i]),
static_cast<int>(mesh.mergeToVert[i]));
}
mesh.mergeToVert.clear();
mesh.mergeFromVert.clear();
for (size_t v = 0; v < numVert; ++v) {
const size_t mergeTo = uf.find(v);
if (mergeTo != v) {
mesh.mergeFromVert.push_back(v);
mesh.mergeToVert.push_back(mergeTo);
}
}
return true;
}
} // namespace
namespace manifold {
/**
* Once halfedge_ has been filled in, this function can be called to create the
* rest of the internal data structures. This function also removes the verts
* and halfedges flagged for removal (NaN verts and -1 halfedges).
*/
void Manifold::Impl::Finish() {
if (halfedge_.size() == 0) return;
CalculateBBox();
SetEpsilon(epsilon_);
if (!bBox_.IsFinite()) {
// Decimated out of existence - early out.
MakeEmpty(Error::NoError);
return;
}
SortVerts();
Vec<Box> faceBox;
Vec<uint32_t> faceMorton;
GetFaceBoxMorton(faceBox, faceMorton);
SortFaces(faceBox, faceMorton);
if (halfedge_.size() == 0) return;
CompactProps();
DEBUG_ASSERT(halfedge_.size() % 6 == 0, topologyErr,
"Not an even number of faces after sorting faces!");
#ifdef MANIFOLD_DEBUG
if (ManifoldParams().intermediateChecks) {
auto MaxOrMinus = [](int a, int b) {
return std::min(a, b) < 0 ? -1 : std::max(a, b);
};
int face = 0;
Halfedge extrema = {0, 0, 0};
for (size_t i = 0; i < halfedge_.size(); i++) {
Halfedge e = halfedge_[i];
if (!e.IsForward()) std::swap(e.startVert, e.endVert);
extrema.startVert = std::min(extrema.startVert, e.startVert);
extrema.endVert = std::min(extrema.endVert, e.endVert);
extrema.pairedHalfedge =
MaxOrMinus(extrema.pairedHalfedge, e.pairedHalfedge);
face = MaxOrMinus(face, i / 3);
}
DEBUG_ASSERT(extrema.startVert >= 0, topologyErr,
"Vertex index is negative!");
DEBUG_ASSERT(extrema.endVert < static_cast<int>(NumVert()), topologyErr,
"Vertex index exceeds number of verts!");
DEBUG_ASSERT(extrema.pairedHalfedge >= 0, topologyErr,
"Halfedge index is negative!");
DEBUG_ASSERT(extrema.pairedHalfedge < 2 * static_cast<int>(NumEdge()),
topologyErr, "Halfedge index exceeds number of halfedges!");
DEBUG_ASSERT(face >= 0, topologyErr, "Face index is negative!");
DEBUG_ASSERT(face < static_cast<int>(NumTri()), topologyErr,
"Face index exceeds number of faces!");
}
#endif
DEBUG_ASSERT(meshRelation_.triRef.size() == NumTri() ||
meshRelation_.triRef.size() == 0,
logicErr, "Mesh Relation doesn't fit!");
DEBUG_ASSERT(faceNormal_.size() == NumTri() || faceNormal_.size() == 0,
logicErr,
"faceNormal size = " + std::to_string(faceNormal_.size()) +
", NumTri = " + std::to_string(NumTri()));
CalculateNormals();
collider_ = Collider(faceBox, faceMorton);
DEBUG_ASSERT(Is2Manifold(), logicErr, "mesh is not 2-manifold!");
}
/**
* Sorts the vertices according to their Morton code.
*/
void Manifold::Impl::SortVerts() {
ZoneScoped;
const auto numVert = NumVert();
Vec<uint32_t> vertMorton(numVert);
auto policy = autoPolicy(numVert, 1e5);
for_each_n(policy, countAt(0), numVert, [this, &vertMorton](const int vert) {
vertMorton[vert] = MortonCode(vertPos_[vert], bBox_);
});
Vec<int> vertNew2Old(numVert);
sequence(vertNew2Old.begin(), vertNew2Old.end());
stable_sort(vertNew2Old.begin(), vertNew2Old.end(),
[&vertMorton](const int& a, const int& b) {
return vertMorton[a] < vertMorton[b];
});
ReindexVerts(vertNew2Old, numVert);
// Verts were flagged for removal with NaNs and assigned kNoCode to sort
// them to the end, which allows them to be removed.
const auto newNumVert =
std::lower_bound(vertNew2Old.begin(), vertNew2Old.end(), kNoCode,
[&vertMorton](const int vert, const uint32_t val) {
return vertMorton[vert] < val;
}) -
vertNew2Old.begin();
vertNew2Old.resize(newNumVert);
Permute(vertPos_, vertNew2Old);
if (vertNormal_.size() == numVert) {
Permute(vertNormal_, vertNew2Old);
}
}
/**
* Updates the halfedges to point to new vert indices based on a mapping,
* vertNew2Old. This may be a subset, so the total number of original verts is
* also given.
*/
void Manifold::Impl::ReindexVerts(const Vec<int>& vertNew2Old,
size_t oldNumVert) {
ZoneScoped;
Vec<int> vertOld2New(oldNumVert);
scatter(countAt(0), countAt(static_cast<int>(NumVert())), vertNew2Old.begin(),
vertOld2New.begin());
const bool hasProp = NumProp() > 0;
for_each(autoPolicy(oldNumVert, 1e5), halfedge_.begin(), halfedge_.end(),
[&vertOld2New, hasProp](Halfedge& edge) {
if (edge.startVert < 0) return;
edge.startVert = vertOld2New[edge.startVert];
edge.endVert = vertOld2New[edge.endVert];
if (!hasProp) {
edge.propVert = edge.startVert;
}
});
}
/**
* Removes unreferenced property verts and reindexes propVerts.
*/
void Manifold::Impl::CompactProps() {
ZoneScoped;
if (numProp_ == 0) return;
const int numProp = NumProp();
const auto numVerts = properties_.size() / numProp;
Vec<int> keep(numVerts, 0);
auto policy = autoPolicy(numVerts, 1e5);
for_each(policy, halfedge_.cbegin(), halfedge_.cend(), [&keep](Halfedge h) {
reinterpret_cast<std::atomic<int>*>(&keep[h.propVert])
->store(1, std::memory_order_relaxed);
});
Vec<int> propOld2New(numVerts + 1, 0);
inclusive_scan(keep.begin(), keep.end(), propOld2New.begin() + 1);
Vec<double> oldProp = properties_;
const int numVertsNew = propOld2New[numVerts];
auto& properties = properties_;
properties.resize_nofill(numProp * numVertsNew);
for_each_n(
policy, countAt(0), numVerts,
[&properties, &oldProp, &propOld2New, &keep, &numProp](const int oldIdx) {
if (keep[oldIdx] == 0) return;
for (int p = 0; p < numProp; ++p) {
properties[propOld2New[oldIdx] * numProp + p] =
oldProp[oldIdx * numProp + p];
}
});
for_each(policy, halfedge_.begin(), halfedge_.end(),
[&propOld2New](Halfedge& edge) {
edge.propVert = propOld2New[edge.propVert];
});
}
/**
* Fills the faceBox and faceMorton input with the bounding boxes and Morton
* codes of the faces, respectively. The Morton code is based on the center of
* the bounding box.
*/
void Manifold::Impl::GetFaceBoxMorton(Vec<Box>& faceBox,
Vec<uint32_t>& faceMorton) const {
ZoneScoped;
// faceBox should be initialized
faceBox.resize(NumTri(), Box());
faceMorton.resize_nofill(NumTri());
for_each_n(autoPolicy(NumTri(), 1e5), countAt(0), NumTri(),
[this, &faceBox, &faceMorton](const int face) {
// Removed tris are marked by all halfedges having pairedHalfedge
// = -1, and this will sort them to the end (the Morton code only
// uses the first 30 of 32 bits).
if (halfedge_[3 * face].pairedHalfedge < 0) {
faceMorton[face] = kNoCode;
return;
}
vec3 center(0.0);
for (const int i : {0, 1, 2}) {
const vec3 pos = vertPos_[halfedge_[3 * face + i].startVert];
center += pos;
faceBox[face].Union(pos);
}
center /= 3;
faceMorton[face] = MortonCode(center, bBox_);
});
}
/**
* Sorts the faces of this manifold according to their input Morton code. The
* bounding box and Morton code arrays are also sorted accordingly.
*/
void Manifold::Impl::SortFaces(Vec<Box>& faceBox, Vec<uint32_t>& faceMorton) {
ZoneScoped;
Vec<int> faceNew2Old(NumTri());
sequence(faceNew2Old.begin(), faceNew2Old.end());
stable_sort(faceNew2Old.begin(), faceNew2Old.end(),
[&faceMorton](const int& a, const int& b) {
return faceMorton[a] < faceMorton[b];
});
// Tris were flagged for removal with pairedHalfedge = -1 and assigned kNoCode
// to sort them to the end, which allows them to be removed.
const int newNumTri =
std::lower_bound(faceNew2Old.begin(), faceNew2Old.end(), kNoCode,
[&faceMorton](const int face, const uint32_t val) {
return faceMorton[face] < val;
}) -
faceNew2Old.begin();
faceNew2Old.resize(newNumTri);
Permute(faceMorton, faceNew2Old);
Permute(faceBox, faceNew2Old);
GatherFaces(faceNew2Old);
}
/**
* Creates the halfedge_ vector for this manifold by copying a set of faces from
* another manifold, given by oldHalfedge. Input faceNew2Old defines the old
* faces to gather into this.
*/
void Manifold::Impl::GatherFaces(const Vec<int>& faceNew2Old) {
ZoneScoped;
const auto numTri = faceNew2Old.size();
if (meshRelation_.triRef.size() == NumTri())
Permute(meshRelation_.triRef, faceNew2Old);
if (faceNormal_.size() == NumTri()) Permute(faceNormal_, faceNew2Old);
Vec<Halfedge> oldHalfedge(std::move(halfedge_));
Vec<vec4> oldHalfedgeTangent(std::move(halfedgeTangent_));
Vec<int> faceOld2New(oldHalfedge.size() / 3);
auto policy = autoPolicy(numTri, 1e5);
scatter(countAt(0_uz), countAt(numTri), faceNew2Old.begin(),
faceOld2New.begin());
halfedge_.resize_nofill(3 * numTri);
if (oldHalfedgeTangent.size() != 0)
halfedgeTangent_.resize_nofill(3 * numTri);
for_each_n(policy, countAt(0), numTri,
ReindexFace({halfedge_, halfedgeTangent_, oldHalfedge,
oldHalfedgeTangent, faceNew2Old, faceOld2New}));
}
void Manifold::Impl::GatherFaces(const Impl& old, const Vec<int>& faceNew2Old) {
ZoneScoped;
const auto numTri = faceNew2Old.size();
meshRelation_.triRef.resize_nofill(numTri);
gather(faceNew2Old.begin(), faceNew2Old.end(),
old.meshRelation_.triRef.begin(), meshRelation_.triRef.begin());
for (const auto& pair : old.meshRelation_.meshIDtransform) {
meshRelation_.meshIDtransform[pair.first] = pair.second;
}
if (old.NumProp() > 0) {
numProp_ = old.numProp_;
properties_ = old.properties_;
}
if (old.faceNormal_.size() == old.NumTri()) {
faceNormal_.resize_nofill(numTri);
gather(faceNew2Old.begin(), faceNew2Old.end(), old.faceNormal_.begin(),
faceNormal_.begin());
}
Vec<int> faceOld2New(old.NumTri());
scatter(countAt(0_uz), countAt(numTri), faceNew2Old.begin(),
faceOld2New.begin());
halfedge_.resize_nofill(3 * numTri);
if (old.halfedgeTangent_.size() != 0)
halfedgeTangent_.resize_nofill(3 * numTri);
for_each_n(autoPolicy(numTri, 1e5), countAt(0), numTri,
ReindexFace({halfedge_, halfedgeTangent_, old.halfedge_,
old.halfedgeTangent_, faceNew2Old, faceOld2New}));
}
/**
* Updates the mergeFromVert and mergeToVert vectors in order to create a
* manifold solid. If the MeshGL is already manifold, no change will occur and
* the function will return false. Otherwise, this will merge verts along open
* edges within tolerance (the maximum of the MeshGL tolerance and the
* baseline bounding-box tolerance), keeping any from the existing merge
* vectors, and return true.
*
* There is no guarantee the result will be manifold - this is a best-effort
* helper function designed primarily to aid in the case where a manifold
* multi-material MeshGL was produced, but its merge vectors were lost due to
* a round-trip through a file format. Constructing a Manifold from the result
* will report an error status if it is not manifold.
*/
template <>
bool MeshGL::Merge() {
return MergeMeshGLP(*this);
}
template <>
bool MeshGL64::Merge() {
return MergeMeshGLP(*this);
}
} // namespace manifold