// 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 <limits>
#if MANIFOLD_PAR == 1
#include <tbb/combinable.h>
#endif
#include "impl.h"
#include "parallel.h"
#include "tri_dist.h"
namespace {
using namespace manifold;
struct CurvatureAngles {
VecView<double> meanCurvature;
VecView<double> gaussianCurvature;
VecView<double> area;
VecView<double> degree;
VecView<const Halfedge> halfedge;
VecView<const vec3> vertPos;
VecView<const vec3> triNormal;
void operator()(size_t tri) {
vec3 edge[3];
vec3 edgeLength(0.0);
for (int i : {0, 1, 2}) {
const int startVert = halfedge[3 * tri + i].startVert;
const int endVert = halfedge[3 * tri + i].endVert;
edge[i] = vertPos[endVert] - vertPos[startVert];
edgeLength[i] = la::length(edge[i]);
edge[i] /= edgeLength[i];
const int neighborTri = halfedge[3 * tri + i].pairedHalfedge / 3;
const double dihedral =
0.25 * edgeLength[i] *
std::asin(la::dot(la::cross(triNormal[tri], triNormal[neighborTri]),
edge[i]));
AtomicAdd(meanCurvature[startVert], dihedral);
AtomicAdd(meanCurvature[endVert], dihedral);
AtomicAdd(degree[startVert], 1.0);
}
vec3 phi;
phi[0] = std::acos(-la::dot(edge[2], edge[0]));
phi[1] = std::acos(-la::dot(edge[0], edge[1]));
phi[2] = kPi - phi[0] - phi[1];
const double area3 = edgeLength[0] * edgeLength[1] *
la::length(la::cross(edge[0], edge[1])) / 6;
for (int i : {0, 1, 2}) {
const int vert = halfedge[3 * tri + i].startVert;
AtomicAdd(gaussianCurvature[vert], -phi[i]);
AtomicAdd(area[vert], area3);
}
}
};
struct CheckHalfedges {
VecView<const Halfedge> halfedges;
bool operator()(size_t edge) const {
const Halfedge halfedge = halfedges[edge];
if (halfedge.startVert == -1 && halfedge.endVert == -1 &&
halfedge.pairedHalfedge == -1)
return true;
if (halfedges[NextHalfedge(edge)].startVert == -1 ||
halfedges[NextHalfedge(NextHalfedge(edge))].startVert == -1) {
return false;
}
if (halfedge.pairedHalfedge == -1) return false;
const Halfedge paired = halfedges[halfedge.pairedHalfedge];
bool good = true;
good &= paired.pairedHalfedge == static_cast<int>(edge);
good &= halfedge.startVert != halfedge.endVert;
good &= halfedge.startVert == paired.endVert;
good &= halfedge.endVert == paired.startVert;
return good;
}
};
struct CheckCCW {
VecView<const Halfedge> halfedges;
VecView<const vec3> vertPos;
VecView<const vec3> triNormal;
const double tol;
bool operator()(size_t face) const {
if (halfedges[3 * face].pairedHalfedge < 0) return true;
const mat2x3 projection = GetAxisAlignedProjection(triNormal[face]);
vec2 v[3];
for (int i : {0, 1, 2})
v[i] = projection * vertPos[halfedges[3 * face + i].startVert];
const int ccw = CCW(v[0], v[1], v[2], std::abs(tol));
return tol > 0 ? ccw >= 0 : ccw == 0;
}
};
} // namespace
namespace manifold {
/**
* Returns true if this manifold is in fact an oriented even manifold and all of
* the data structures are consistent.
*/
bool Manifold::Impl::IsManifold() const {
if (halfedge_.size() == 0) return true;
if (halfedge_.size() % 3 != 0) return false;
return all_of(countAt(0_uz), countAt(halfedge_.size()),
CheckHalfedges({halfedge_}));
}
/**
* Returns true if this manifold is in fact an oriented 2-manifold and all of
* the data structures are consistent.
*/
bool Manifold::Impl::Is2Manifold() const {
if (halfedge_.size() == 0) return true;
if (!IsManifold()) return false;
Vec<Halfedge> halfedge(halfedge_);
stable_sort(halfedge.begin(), halfedge.end());
return all_of(
countAt(0_uz), countAt(2 * NumEdge() - 1), [&halfedge](size_t edge) {
const Halfedge h = halfedge[edge];
if (h.startVert == -1 && h.endVert == -1 && h.pairedHalfedge == -1)
return true;
return h.startVert != halfedge[edge + 1].startVert ||
h.endVert != halfedge[edge + 1].endVert;
});
}
#ifdef MANIFOLD_DEBUG
std::mutex dump_lock;
#endif
/**
* Returns true if this manifold is self-intersecting.
* Note that this is not checking for epsilon-validity.
*/
bool Manifold::Impl::IsSelfIntersecting() const {
const double ep = 2 * epsilon_;
const double epsilonSq = ep * ep;
Vec<Box> faceBox;
Vec<uint32_t> faceMorton;
GetFaceBoxMorton(faceBox, faceMorton);
std::atomic<bool> intersecting(false);
auto f = [&](int tri0, int tri1) {
std::array<vec3, 3> triVerts0, triVerts1;
for (int i : {0, 1, 2}) {
triVerts0[i] = vertPos_[halfedge_[3 * tri0 + i].startVert];
triVerts1[i] = vertPos_[halfedge_[3 * tri1 + i].startVert];
}
// if triangles tri0 and tri1 share a vertex, return true to skip the
// check. we relax the sharing criteria a bit to allow for at most
// distance epsilon squared
for (int i : {0, 1, 2})
for (int j : {0, 1, 2})
if (distance2(triVerts0[i], triVerts1[j]) <= epsilonSq) return;
if (DistanceTriangleTriangleSquared(triVerts0, triVerts1) == 0.0) {
// try to move the triangles around the normal of the other face
std::array<vec3, 3> tmp0, tmp1;
for (int i : {0, 1, 2}) tmp0[i] = triVerts0[i] + ep * faceNormal_[tri1];
if (DistanceTriangleTriangleSquared(tmp0, triVerts1) > 0.0) return;
for (int i : {0, 1, 2}) tmp0[i] = triVerts0[i] - ep * faceNormal_[tri1];
if (DistanceTriangleTriangleSquared(tmp0, triVerts1) > 0.0) return;
for (int i : {0, 1, 2}) tmp1[i] = triVerts1[i] + ep * faceNormal_[tri0];
if (DistanceTriangleTriangleSquared(triVerts0, tmp1) > 0.0) return;
for (int i : {0, 1, 2}) tmp1[i] = triVerts1[i] - ep * faceNormal_[tri0];
if (DistanceTriangleTriangleSquared(triVerts0, tmp1) > 0.0) return;
#ifdef MANIFOLD_DEBUG
if (ManifoldParams().verbose > 0) {
dump_lock.lock();
std::cout << "intersecting:" << std::endl;
for (int i : {0, 1, 2}) std::cout << triVerts0[i] << " ";
std::cout << std::endl;
for (int i : {0, 1, 2}) std::cout << triVerts1[i] << " ";
std::cout << std::endl;
dump_lock.unlock();
}
#endif
intersecting.store(true);
}
};
auto recorder = MakeSimpleRecorder(f);
collider_.Collisions<true>(faceBox.cview(), recorder);
return intersecting.load();
}
/**
* Returns true if all triangles are CCW relative to their triNormals_.
*/
bool Manifold::Impl::MatchesTriNormals() const {
if (halfedge_.size() == 0 || faceNormal_.size() != NumTri()) return true;
return all_of(countAt(0_uz), countAt(NumTri()),
CheckCCW({halfedge_, vertPos_, faceNormal_, 2 * epsilon_}));
}
/**
* Returns the number of triangles that are colinear within epsilon_.
*/
int Manifold::Impl::NumDegenerateTris() const {
if (halfedge_.size() == 0 || faceNormal_.size() != NumTri()) return true;
return count_if(
countAt(0_uz), countAt(NumTri()),
CheckCCW({halfedge_, vertPos_, faceNormal_, -1 * epsilon_ / 2}));
}
double Manifold::Impl::GetProperty(Property prop) const {
ZoneScoped;
if (IsEmpty()) return 0;
auto Volume = [this](size_t tri) {
const vec3 v = vertPos_[halfedge_[3 * tri].startVert];
vec3 crossP = la::cross(vertPos_[halfedge_[3 * tri + 1].startVert] - v,
vertPos_[halfedge_[3 * tri + 2].startVert] - v);
return la::dot(crossP, v) / 6.0;
};
auto Area = [this](size_t tri) {
const vec3 v = vertPos_[halfedge_[3 * tri].startVert];
return la::length(
la::cross(vertPos_[halfedge_[3 * tri + 1].startVert] - v,
vertPos_[halfedge_[3 * tri + 2].startVert] - v)) /
2.0;
};
// Kahan summation
double value = 0;
double valueCompensation = 0;
for (size_t i = 0; i < NumTri(); ++i) {
const double value1 = prop == Property::SurfaceArea ? Area(i) : Volume(i);
const double t = value + value1;
valueCompensation += (value - t) + value1;
value = t;
}
value += valueCompensation;
return value;
}
void Manifold::Impl::CalculateCurvature(int gaussianIdx, int meanIdx) {
ZoneScoped;
if (IsEmpty()) return;
if (gaussianIdx < 0 && meanIdx < 0) return;
Vec<double> vertMeanCurvature(NumVert(), 0);
Vec<double> vertGaussianCurvature(NumVert(), kTwoPi);
Vec<double> vertArea(NumVert(), 0);
Vec<double> degree(NumVert(), 0);
auto policy = autoPolicy(NumTri(), 1e4);
for_each(policy, countAt(0_uz), countAt(NumTri()),
CurvatureAngles({vertMeanCurvature, vertGaussianCurvature, vertArea,
degree, halfedge_, vertPos_, faceNormal_}));
for_each_n(policy, countAt(0), NumVert(),
[&vertMeanCurvature, &vertGaussianCurvature, &vertArea,
°ree](const int vert) {
const double factor = degree[vert] / (6 * vertArea[vert]);
vertMeanCurvature[vert] *= factor;
vertGaussianCurvature[vert] *= factor;
});
const int oldNumProp = NumProp();
const int numProp = std::max(oldNumProp, std::max(gaussianIdx, meanIdx) + 1);
const Vec<double> oldProperties = properties_;
properties_ = Vec<double>(numProp * NumPropVert(), 0);
numProp_ = numProp;
Vec<uint8_t> counters(NumPropVert(), 0);
for_each_n(policy, countAt(0_uz), NumTri(), [&](const size_t tri) {
for (const int i : {0, 1, 2}) {
const Halfedge& edge = halfedge_[3 * tri + i];
const int vert = edge.startVert;
const int propVert = edge.propVert;
auto old = std::atomic_exchange(
reinterpret_cast<std::atomic<uint8_t>*>(&counters[propVert]),
static_cast<uint8_t>(1));
if (old == 1) continue;
for (int p = 0; p < oldNumProp; ++p) {
properties_[numProp * propVert + p] =
oldProperties[oldNumProp * propVert + p];
}
if (gaussianIdx >= 0) {
properties_[numProp * propVert + gaussianIdx] =
vertGaussianCurvature[vert];
}
if (meanIdx >= 0) {
properties_[numProp * propVert + meanIdx] = vertMeanCurvature[vert];
}
}
});
}
/**
* Calculates the bounding box of the entire manifold, which is stored
* internally to short-cut Boolean operations. Ignores NaNs.
*/
void Manifold::Impl::CalculateBBox() {
bBox_.min =
reduce(vertPos_.begin(), vertPos_.end(),
vec3(std::numeric_limits<double>::infinity()), [](auto a, auto b) {
if (std::isnan(a.x)) return b;
if (std::isnan(b.x)) return a;
return la::min(a, b);
});
bBox_.max = reduce(vertPos_.begin(), vertPos_.end(),
vec3(-std::numeric_limits<double>::infinity()),
[](auto a, auto b) {
if (std::isnan(a.x)) return b;
if (std::isnan(b.x)) return a;
return la::max(a, b);
});
}
/**
* Determines if all verts are finite. Checking just the bounding box dimensions
* is insufficient as it ignores NaNs.
*/
bool Manifold::Impl::IsFinite() const {
return transform_reduce(
vertPos_.begin(), vertPos_.end(), true,
[](bool a, bool b) { return a && b; },
[](auto v) { return la::all(la::isfinite(v)); });
}
/**
* Checks that the input triVerts array has all indices inside bounds of the
* vertPos_ array.
*/
bool Manifold::Impl::IsIndexInBounds(VecView<const ivec3> triVerts) const {
ivec2 minmax = transform_reduce(
triVerts.begin(), triVerts.end(),
ivec2(std::numeric_limits<int>::max(), std::numeric_limits<int>::min()),
[](auto a, auto b) {
a[0] = std::min(a[0], b[0]);
a[1] = std::max(a[1], b[1]);
return a;
},
[](auto tri) {
return ivec2(std::min(tri[0], std::min(tri[1], tri[2])),
std::max(tri[0], std::max(tri[1], tri[2])));
});
return minmax[0] >= 0 && minmax[1] < static_cast<int>(NumVert());
}
struct MinDistanceRecorder {
using Local = double;
const Manifold::Impl &self, &other;
#if MANIFOLD_PAR == 1
tbb::combinable<double> store = tbb::combinable<double>(
[]() { return std::numeric_limits<double>::infinity(); });
Local& local() { return store.local(); }
double get() {
double result = std::numeric_limits<double>::infinity();
store.combine_each([&](double& val) { result = std::min(result, val); });
return result;
}
#else
double result = std::numeric_limits<double>::infinity();
Local& local() { return result; }
double get() { return result; }
#endif
void record(int triOther, int tri, double& minDistance) {
std::array<vec3, 3> p;
std::array<vec3, 3> q;
for (const int j : {0, 1, 2}) {
p[j] = self.vertPos_[self.halfedge_[3 * tri + j].startVert];
q[j] = other.vertPos_[other.halfedge_[3 * triOther + j].startVert];
}
minDistance = std::min(minDistance, DistanceTriangleTriangleSquared(p, q));
}
};
/*
* Returns the minimum gap between two manifolds. Returns a double between
* 0 and searchLength.
*/
double Manifold::Impl::MinGap(const Manifold::Impl& other,
double searchLength) const {
ZoneScoped;
Vec<Box> faceBoxOther;
Vec<uint32_t> faceMortonOther;
other.GetFaceBoxMorton(faceBoxOther, faceMortonOther);
transform(faceBoxOther.begin(), faceBoxOther.end(), faceBoxOther.begin(),
[searchLength](const Box& box) {
return Box(box.min - vec3(searchLength),
box.max + vec3(searchLength));
});
MinDistanceRecorder recorder{*this, other};
collider_.Collisions<false, Box, MinDistanceRecorder>(faceBoxOther.cview(),
recorder, false);
double minDistanceSquared =
std::min(recorder.get(), searchLength * searchLength);
return sqrt(minDistanceSquared);
};
} // namespace manifold