// 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 <algorithm>
#include <map>
#include <numeric>
#include "boolean3.h"
#include "csg_tree.h"
#include "impl.h"
#include "parallel.h"
#include "shared.h"
namespace {
using namespace manifold;
ExecutionParams manifoldParams;
Manifold Halfspace(Box bBox, vec3 normal, double originOffset) {
normal = la::normalize(normal);
Manifold cutter = Manifold::Cube(vec3(2.0), true).Translate({1.0, 0.0, 0.0});
double size = la::length(bBox.Center() - normal * originOffset) +
0.5 * la::length(bBox.Size());
cutter = cutter.Scale(vec3(size)).Translate({originOffset, 0.0, 0.0});
double yDeg = degrees(-std::asin(normal.z));
double zDeg = degrees(std::atan2(normal.y, normal.x));
return cutter.Rotate(0.0, yDeg, zDeg);
}
template <typename Precision, typename I>
MeshGLP<Precision, I> GetMeshGLImpl(const manifold::Manifold::Impl& impl,
int normalIdx) {
ZoneScoped;
const int numProp = impl.NumProp();
const int numVert = impl.NumPropVert();
const int numTri = impl.NumTri();
const bool isOriginal = impl.meshRelation_.originalID >= 0;
const bool updateNormals = !isOriginal && normalIdx >= 0;
MeshGLP<Precision, I> out;
out.numProp = 3 + numProp;
out.tolerance = impl.tolerance_;
if (std::is_same<Precision, float>::value)
out.tolerance =
std::max(out.tolerance,
static_cast<Precision>(std::numeric_limits<float>::epsilon() *
impl.bBox_.Scale()));
out.triVerts.resize(3 * numTri);
const int numHalfedge = impl.halfedgeTangent_.size();
out.halfedgeTangent.resize(4 * numHalfedge);
for (int i = 0; i < numHalfedge; ++i) {
const vec4 t = impl.halfedgeTangent_[i];
out.halfedgeTangent[4 * i] = t.x;
out.halfedgeTangent[4 * i + 1] = t.y;
out.halfedgeTangent[4 * i + 2] = t.z;
out.halfedgeTangent[4 * i + 3] = t.w;
}
// Sort the triangles into runs
out.faceID.resize(numTri);
std::vector<int> triNew2Old(numTri);
std::iota(triNew2Old.begin(), triNew2Old.end(), 0);
VecView<const TriRef> triRef = impl.meshRelation_.triRef;
// Don't sort originals - keep them in order
if (!isOriginal) {
std::stable_sort(triNew2Old.begin(), triNew2Old.end(),
[triRef](int a, int b) {
return triRef[a].originalID == triRef[b].originalID
? triRef[a].meshID < triRef[b].meshID
: triRef[a].originalID < triRef[b].originalID;
});
}
std::vector<mat3> runNormalTransform;
auto addRun = [updateNormals, isOriginal](
MeshGLP<Precision, I>& out,
std::vector<mat3>& runNormalTransform, int tri,
const manifold::Manifold::Impl::Relation& rel) {
out.runIndex.push_back(3 * tri);
out.runOriginalID.push_back(rel.originalID);
if (updateNormals) {
runNormalTransform.push_back(NormalTransform(rel.transform) *
(rel.backSide ? -1.0 : 1.0));
}
if (!isOriginal) {
for (const int col : {0, 1, 2, 3}) {
for (const int row : {0, 1, 2}) {
out.runTransform.push_back(rel.transform[col][row]);
}
}
}
};
auto meshIDtransform = impl.meshRelation_.meshIDtransform;
int lastID = -1;
for (int tri = 0; tri < numTri; ++tri) {
const int oldTri = triNew2Old[tri];
const auto ref = triRef[oldTri];
const int meshID = ref.meshID;
out.faceID[tri] = ref.faceID >= 0 ? ref.faceID : ref.coplanarID;
for (const int i : {0, 1, 2})
out.triVerts[3 * tri + i] = impl.halfedge_[3 * oldTri + i].startVert;
if (meshID != lastID) {
manifold::Manifold::Impl::Relation rel;
auto it = meshIDtransform.find(meshID);
if (it != meshIDtransform.end()) rel = it->second;
addRun(out, runNormalTransform, tri, rel);
meshIDtransform.erase(meshID);
lastID = meshID;
}
}
// Add runs for originals that did not contribute any faces to the output
for (const auto& pair : meshIDtransform) {
addRun(out, runNormalTransform, numTri, pair.second);
}
out.runIndex.push_back(3 * numTri);
// Early return for no props
if (numProp == 0) {
out.vertProperties.resize(3 * numVert);
for (int i = 0; i < numVert; ++i) {
const vec3 v = impl.vertPos_[i];
out.vertProperties[3 * i] = v.x;
out.vertProperties[3 * i + 1] = v.y;
out.vertProperties[3 * i + 2] = v.z;
}
return out;
}
// Duplicate verts with different props
std::vector<int> vert2idx(impl.NumVert(), -1);
std::vector<std::vector<ivec2>> vertPropPair(impl.NumVert());
out.vertProperties.reserve(numVert * static_cast<size_t>(out.numProp));
for (size_t run = 0; run < out.runOriginalID.size(); ++run) {
for (size_t tri = out.runIndex[run] / 3; tri < out.runIndex[run + 1] / 3;
++tri) {
for (const int i : {0, 1, 2}) {
const int prop = impl.halfedge_[3 * triNew2Old[tri] + i].propVert;
const int vert = out.triVerts[3 * tri + i];
auto& bin = vertPropPair[vert];
bool bFound = false;
for (const auto& b : bin) {
if (b.x == prop) {
bFound = true;
out.triVerts[3 * tri + i] = b.y;
break;
}
}
if (bFound) continue;
const int idx = out.vertProperties.size() / out.numProp;
out.triVerts[3 * tri + i] = idx;
bin.push_back({prop, idx});
for (int p : {0, 1, 2}) {
out.vertProperties.push_back(impl.vertPos_[vert][p]);
}
for (int p = 0; p < numProp; ++p) {
out.vertProperties.push_back(impl.properties_[prop * numProp + p]);
}
if (updateNormals) {
vec3 normal;
const int start = out.vertProperties.size() - out.numProp;
for (int i : {0, 1, 2}) {
normal[i] = out.vertProperties[start + 3 + normalIdx + i];
}
normal = la::normalize(runNormalTransform[run] * normal);
for (int i : {0, 1, 2}) {
out.vertProperties[start + 3 + normalIdx + i] = normal[i];
}
}
if (vert2idx[vert] == -1) {
vert2idx[vert] = idx;
} else {
out.mergeFromVert.push_back(idx);
out.mergeToVert.push_back(vert2idx[vert]);
}
}
}
}
return out;
}
} // namespace
namespace manifold {
static int circularSegments_ = DEFAULT_SEGMENTS;
static double circularAngle_ = DEFAULT_ANGLE;
static double circularEdgeLength_ = DEFAULT_LENGTH;
/**
* Sets an angle constraint the default number of circular segments for the
* CrossSection::Circle(), Manifold::Cylinder(), Manifold::Sphere(), and
* Manifold::Revolve() constructors. The number of segments will be rounded up
* to the nearest factor of four.
*
* @param angle The minimum angle in degrees between consecutive segments. The
* angle will increase if the the segments hit the minimum edge length.
* Default is 10 degrees.
*/
void Quality::SetMinCircularAngle(double angle) {
if (angle <= 0) return;
circularAngle_ = angle;
}
/**
* Sets a length constraint the default number of circular segments for the
* CrossSection::Circle(), Manifold::Cylinder(), Manifold::Sphere(), and
* Manifold::Revolve() constructors. The number of segments will be rounded up
* to the nearest factor of four.
*
* @param length The minimum length of segments. The length will
* increase if the the segments hit the minimum angle. Default is 1.0.
*/
void Quality::SetMinCircularEdgeLength(double length) {
if (length <= 0) return;
circularEdgeLength_ = length;
}
/**
* Sets the default number of circular segments for the
* CrossSection::Circle(), Manifold::Cylinder(), Manifold::Sphere(), and
* Manifold::Revolve() constructors. Overrides the edge length and angle
* constraints and sets the number of segments to exactly this value.
*
* @param number Number of circular segments. Default is 0, meaning no
* constraint is applied.
*/
void Quality::SetCircularSegments(int number) {
if (number < 3 && number != 0) return;
circularSegments_ = number;
}
/**
* Determine the result of the SetMinCircularAngle(),
* SetMinCircularEdgeLength(), and SetCircularSegments() defaults.
*
* @param radius For a given radius of circle, determine how many default
* segments there will be.
*/
int Quality::GetCircularSegments(double radius) {
if (circularSegments_ > 0) return circularSegments_;
int nSegA = 360.0 / circularAngle_;
int nSegL = 2.0 * radius * kPi / circularEdgeLength_;
int nSeg = fmin(nSegA, nSegL) + 3;
nSeg -= nSeg % 4;
return std::max(nSeg, 4);
}
/**
* Resets the circular construction parameters to their defaults if
* SetMinCircularAngle, SetMinCircularEdgeLength, or SetCircularSegments have
* been called.
*/
void Quality::ResetToDefaults() {
circularSegments_ = DEFAULT_SEGMENTS;
circularAngle_ = DEFAULT_ANGLE;
circularEdgeLength_ = DEFAULT_LENGTH;
}
/**
* Construct an empty Manifold.
*
*/
Manifold::Manifold() : pNode_{std::make_shared<CsgLeafNode>()} {}
Manifold::~Manifold() = default;
Manifold::Manifold(Manifold&&) noexcept = default;
Manifold& Manifold::operator=(Manifold&&) noexcept = default;
Manifold::Manifold(const Manifold& other) : pNode_(other.pNode_) {}
Manifold::Manifold(std::shared_ptr<CsgNode> pNode) : pNode_(pNode) {}
Manifold::Manifold(std::shared_ptr<Impl> pImpl_)
: pNode_(std::make_shared<CsgLeafNode>(pImpl_)) {}
Manifold Manifold::Invalid() {
auto pImpl_ = std::make_shared<Impl>();
pImpl_->status_ = Error::InvalidConstruction;
return Manifold(pImpl_);
}
Manifold& Manifold::operator=(const Manifold& other) {
if (this != &other) {
pNode_ = other.pNode_;
}
return *this;
}
CsgLeafNode& Manifold::GetCsgLeafNode() const {
if (pNode_->GetNodeType() != CsgNodeType::Leaf) {
pNode_ = pNode_->ToLeafNode();
}
return *std::static_pointer_cast<CsgLeafNode>(pNode_);
}
/**
* Convert a MeshGL into a Manifold, retaining its properties and merging only
* the positions according to the merge vectors. Will return an empty Manifold
* and set an Error Status if the result is not an oriented 2-manifold. Will
* collapse degenerate triangles and unnecessary vertices.
*
* All fields are read, making this structure suitable for a lossless round-trip
* of data from GetMeshGL. For multi-material input, use ReserveIDs to set a
* unique originalID for each material, and sort the materials into triangle
* runs.
*
* @param meshGL The input MeshGL.
*/
Manifold::Manifold(const MeshGL& meshGL)
: pNode_(std::make_shared<CsgLeafNode>(std::make_shared<Impl>(meshGL))) {}
/**
* Convert a MeshGL into a Manifold, retaining its properties and merging only
* the positions according to the merge vectors. Will return an empty Manifold
* and set an Error Status if the result is not an oriented 2-manifold. Will
* collapse degenerate triangles and unnecessary vertices.
*
* All fields are read, making this structure suitable for a lossless round-trip
* of data from GetMeshGL. For multi-material input, use ReserveIDs to set a
* unique originalID for each material, and sort the materials into triangle
* runs.
*
* @param meshGL64 The input MeshGL64.
*/
Manifold::Manifold(const MeshGL64& meshGL64)
: pNode_(std::make_shared<CsgLeafNode>(std::make_shared<Impl>(meshGL64))) {}
/**
* The most complete output of this library, returning a MeshGL that is designed
* to easily push into a renderer, including all interleaved vertex properties
* that may have been input. It also includes relations to all the input meshes
* that form a part of this result and the transforms applied to each.
*
* @param normalIdx If the original MeshGL inputs that formed this manifold had
* properties corresponding to normal vectors, you can specify the first of the
* three consecutive property channels forming the (x, y, z) normals, which will
* cause this output MeshGL to automatically update these normals according to
* the applied transforms and front/back side. normalIdx + 3 must be <=
* numProp, and all original MeshGLs must use the same channels for their
* normals.
*/
MeshGL Manifold::GetMeshGL(int normalIdx) const {
const Impl& impl = *GetCsgLeafNode().GetImpl();
return GetMeshGLImpl<float, uint32_t>(impl, normalIdx);
}
/**
* The most complete output of this library, returning a MeshGL that is designed
* to easily push into a renderer, including all interleaved vertex properties
* that may have been input. It also includes relations to all the input meshes
* that form a part of this result and the transforms applied to each.
*
* @param normalIdx If the original MeshGL inputs that formed this manifold had
* properties corresponding to normal vectors, you can specify the first of the
* three consecutive property channels forming the (x, y, z) normals, which will
* cause this output MeshGL to automatically update these normals according to
* the applied transforms and front/back side. normalIdx + 3 must be <=
* numProp, and all original MeshGLs must use the same channels for their
* normals.
*/
MeshGL64 Manifold::GetMeshGL64(int normalIdx) const {
const Impl& impl = *GetCsgLeafNode().GetImpl();
return GetMeshGLImpl<double, uint64_t>(impl, normalIdx);
}
/**
* Does the Manifold have any triangles?
*/
bool Manifold::IsEmpty() const { return GetCsgLeafNode().GetImpl()->IsEmpty(); }
/**
* Returns the reason for an input Mesh producing an empty Manifold. This Status
* will carry on through operations like NaN propogation, ensuring an errored
* mesh doesn't get mysteriously lost. Empty meshes may still show
* NoError, for instance the intersection of non-overlapping meshes.
*/
Manifold::Error Manifold::Status() const {
return GetCsgLeafNode().GetImpl()->status_;
}
/**
* The number of vertices in the Manifold.
*/
size_t Manifold::NumVert() const {
return GetCsgLeafNode().GetImpl()->NumVert();
}
/**
* The number of edges in the Manifold.
*/
size_t Manifold::NumEdge() const {
return GetCsgLeafNode().GetImpl()->NumEdge();
}
/**
* The number of triangles in the Manifold.
*/
size_t Manifold::NumTri() const { return GetCsgLeafNode().GetImpl()->NumTri(); }
/**
* The number of properties per vertex in the Manifold.
*/
size_t Manifold::NumProp() const {
return GetCsgLeafNode().GetImpl()->NumProp();
}
/**
* The number of property vertices in the Manifold. This will always be >=
* NumVert, as some physical vertices may be duplicated to account for different
* properties on different neighboring triangles.
*/
size_t Manifold::NumPropVert() const {
return GetCsgLeafNode().GetImpl()->NumPropVert();
}
/**
* Returns the axis-aligned bounding box of all the Manifold's vertices.
*/
Box Manifold::BoundingBox() const { return GetCsgLeafNode().GetImpl()->bBox_; }
/**
* Returns the epsilon value of this Manifold's vertices, which tracks the
* approximate rounding error over all the transforms and operations that have
* led to this state. This is the value of ε defining
* [ε-valid](https://github.com/elalish/manifold/wiki/Manifold-Library#definition-of-%CE%B5-valid).
*/
double Manifold::GetEpsilon() const {
return GetCsgLeafNode().GetImpl()->epsilon_;
}
/**
* Returns the tolerance value of this Manifold. Triangles that are coplanar
* within tolerance tend to be merged and edges shorter than tolerance tend to
* be collapsed.
*/
double Manifold::GetTolerance() const {
return GetCsgLeafNode().GetImpl()->tolerance_;
}
/**
* Return a copy of the manifold with the set tolerance value.
* This performs mesh simplification when the tolerance value is increased.
*/
Manifold Manifold::SetTolerance(double tolerance) const {
auto impl = std::make_shared<Impl>(*GetCsgLeafNode().GetImpl());
if (tolerance > impl->tolerance_) {
impl->tolerance_ = tolerance;
impl->MarkCoplanar();
impl->SimplifyTopology();
impl->Finish();
} else {
// for reducing tolerance, we need to make sure it is still at least
// equal to epsilon.
impl->tolerance_ = std::max(impl->epsilon_, tolerance);
}
return Manifold(impl);
}
/**
* Return a copy of the manifold simplified to the given tolerance, but with its
* actual tolerance value unchanged. If the tolerance is not given or is less
* than the current tolerance, the current tolerance is used for simplification.
* The result will contain a subset of the original verts and all surfaces will
* have moved by less than tolerance.
*/
Manifold Manifold::Simplify(double tolerance) const {
auto impl = std::make_shared<Impl>(*GetCsgLeafNode().GetImpl());
const double oldTolerance = impl->tolerance_;
if (tolerance == 0) tolerance = oldTolerance;
if (tolerance > oldTolerance) {
impl->tolerance_ = tolerance;
impl->MarkCoplanar();
}
impl->SimplifyTopology();
impl->Finish();
impl->tolerance_ = oldTolerance;
return Manifold(impl);
}
/**
* The genus is a topological property of the manifold, representing the number
* of "handles". A sphere is 0, torus 1, etc. It is only meaningful for a single
* mesh, so it is best to call Decompose() first.
*/
int Manifold::Genus() const {
int chi = NumVert() - NumEdge() + NumTri();
return 1 - chi / 2;
}
/**
* Returns the surface area of the manifold.
*/
double Manifold::SurfaceArea() const {
return GetCsgLeafNode().GetImpl()->GetProperty(Impl::Property::SurfaceArea);
}
/**
* Returns the volume of the manifold.
*/
double Manifold::Volume() const {
return GetCsgLeafNode().GetImpl()->GetProperty(Impl::Property::Volume);
}
/**
* If this mesh is an original, this returns its meshID that can be referenced
* by product manifolds' MeshRelation. If this manifold is a product, this
* returns -1.
*/
int Manifold::OriginalID() const {
return GetCsgLeafNode().GetImpl()->meshRelation_.originalID;
}
/**
* This removes all relations (originalID, faceID, transform) to ancestor meshes
* and this new Manifold is marked an original. It also recreates faces
* - these don't get joined at boundaries where originalID changes, so the
* reset may allow triangles of flat faces to be further collapsed with
* Simplify().
*/
Manifold Manifold::AsOriginal() const {
auto oldImpl = GetCsgLeafNode().GetImpl();
if (oldImpl->status_ != Error::NoError) {
auto newImpl = std::make_shared<Impl>();
newImpl->status_ = oldImpl->status_;
return Manifold(std::make_shared<CsgLeafNode>(newImpl));
}
auto newImpl = std::make_shared<Impl>(*oldImpl);
newImpl->InitializeOriginal();
newImpl->MarkCoplanar();
newImpl->InitializeOriginal(true);
return Manifold(std::make_shared<CsgLeafNode>(newImpl));
}
/**
* Returns the first of n sequential new unique mesh IDs for marking sets of
* triangles that can be looked up after further operations. Assign to
* MeshGL.runOriginalID vector.
*/
uint32_t Manifold::ReserveIDs(uint32_t n) {
return Manifold::Impl::ReserveIDs(n);
}
/**
* The triangle normal vectors are saved over the course of operations rather
* than recalculated to avoid rounding error. This checks that triangles still
* match their normal vectors within Precision().
*/
bool Manifold::MatchesTriNormals() const {
return GetCsgLeafNode().GetImpl()->MatchesTriNormals();
}
/**
* The number of triangles that are colinear within Precision(). This library
* attempts to remove all of these, but it cannot always remove all of them
* without changing the mesh by too much.
*/
size_t Manifold::NumDegenerateTris() const {
return GetCsgLeafNode().GetImpl()->NumDegenerateTris();
}
/**
* Move this Manifold in space. This operation can be chained. Transforms are
* combined and applied lazily.
*
* @param v The vector to add to every vertex.
*/
Manifold Manifold::Translate(vec3 v) const {
return Manifold(pNode_->Translate(v));
}
/**
* Scale this Manifold in space. This operation can be chained. Transforms are
* combined and applied lazily.
*
* @param v The vector to multiply every vertex by per component.
*/
Manifold Manifold::Scale(vec3 v) const { return Manifold(pNode_->Scale(v)); }
/**
* Applies an Euler angle rotation to the manifold, first about the X axis, then
* Y, then Z, in degrees. We use degrees so that we can minimize rounding error,
* and eliminate it completely for any multiples of 90 degrees. Additionally,
* more efficient code paths are used to update the manifold when the transforms
* only rotate by multiples of 90 degrees. This operation can be chained.
* Transforms are combined and applied lazily.
*
* @param xDegrees First rotation, degrees about the X-axis.
* @param yDegrees Second rotation, degrees about the Y-axis.
* @param zDegrees Third rotation, degrees about the Z-axis.
*/
Manifold Manifold::Rotate(double xDegrees, double yDegrees,
double zDegrees) const {
return Manifold(pNode_->Rotate(xDegrees, yDegrees, zDegrees));
}
/**
* Transform this Manifold in space. The first three columns form a 3x3 matrix
* transform and the last is a translation vector. This operation can be
* chained. Transforms are combined and applied lazily.
*
* @param m The affine transform matrix to apply to all the vertices.
*/
Manifold Manifold::Transform(const mat3x4& m) const {
return Manifold(pNode_->Transform(m));
}
/**
* Mirror this Manifold over the plane described by the unit form of the given
* normal vector. If the length of the normal is zero, an empty Manifold is
* returned. This operation can be chained. Transforms are combined and applied
* lazily.
*
* @param normal The normal vector of the plane to be mirrored over
*/
Manifold Manifold::Mirror(vec3 normal) const {
if (la::length(normal) == 0.) {
return Manifold();
}
auto n = la::normalize(normal);
auto m = mat3x4(mat3(la::identity) - 2.0 * la::outerprod(n, n), vec3());
return Manifold(pNode_->Transform(m));
}
/**
* This function does not change the topology, but allows the vertices to be
* moved according to any arbitrary input function. It is easy to create a
* function that warps a geometrically valid object into one which overlaps, but
* that is not checked here, so it is up to the user to choose their function
* with discretion.
*
* @param warpFunc A function that modifies a given vertex position.
*/
Manifold Manifold::Warp(std::function<void(vec3&)> warpFunc) const {
auto oldImpl = GetCsgLeafNode().GetImpl();
if (oldImpl->status_ != Error::NoError) {
auto pImpl = std::make_shared<Impl>();
pImpl->status_ = oldImpl->status_;
return Manifold(std::make_shared<CsgLeafNode>(pImpl));
}
auto pImpl = std::make_shared<Impl>(*oldImpl);
pImpl->Warp(warpFunc);
return Manifold(std::make_shared<CsgLeafNode>(pImpl));
}
/**
* Same as Manifold::Warp but calls warpFunc with with
* a VecView which is roughly equivalent to std::span
* pointing to all vec3 elements to be modified in-place
*
* @param warpFunc A function that modifies multiple vertex positions.
*/
Manifold Manifold::WarpBatch(
std::function<void(VecView<vec3>)> warpFunc) const {
auto oldImpl = GetCsgLeafNode().GetImpl();
if (oldImpl->status_ != Error::NoError) {
auto pImpl = std::make_shared<Impl>();
pImpl->status_ = oldImpl->status_;
return Manifold(std::make_shared<CsgLeafNode>(pImpl));
}
auto pImpl = std::make_shared<Impl>(*oldImpl);
pImpl->WarpBatch(warpFunc);
return Manifold(std::make_shared<CsgLeafNode>(pImpl));
}
/**
* Create a new copy of this manifold with updated vertex properties by
* supplying a function that takes the existing position and properties as
* input. You may specify any number of output properties, allowing creation and
* removal of channels. Note: undefined behavior will result if you read past
* the number of input properties or write past the number of output properties.
*
* If propFunc is a nullptr, this function will just set the channel to zeroes.
*
* @param numProp The new number of properties per vertex.
* @param propFunc A function that modifies the properties of a given vertex.
*/
Manifold Manifold::SetProperties(
int numProp,
std::function<void(double* newProp, vec3 position, const double* oldProp)>
propFunc) const {
auto pImpl = std::make_shared<Impl>(*GetCsgLeafNode().GetImpl());
const int oldNumProp = NumProp();
const Vec<double> oldProperties = pImpl->properties_;
if (numProp == 0) {
pImpl->properties_.clear();
} else {
pImpl->properties_ = Vec<double>(numProp * NumPropVert(), 0);
for_each_n(
propFunc == nullptr ? ExecutionPolicy::Par : ExecutionPolicy::Seq,
countAt(0), NumTri(), [&](int tri) {
for (int i : {0, 1, 2}) {
const Halfedge& edge = pImpl->halfedge_[3 * tri + i];
const int vert = edge.startVert;
const int propVert = edge.propVert;
if (propFunc == nullptr) {
for (int p = 0; p < numProp; ++p) {
pImpl->properties_[numProp * propVert + p] = 0;
}
} else {
propFunc(&pImpl->properties_[numProp * propVert],
pImpl->vertPos_[vert],
oldProperties.data() + oldNumProp * propVert);
}
}
});
}
pImpl->numProp_ = numProp;
return Manifold(std::make_shared<CsgLeafNode>(pImpl));
}
/**
* Curvature is the inverse of the radius of curvature, and signed such that
* positive is convex and negative is concave. There are two orthogonal
* principal curvatures at any point on a manifold, with one maximum and the
* other minimum. Gaussian curvature is their product, while mean
* curvature is their sum. This approximates them for every vertex and assigns
* them as vertex properties on the given channels.
*
* @param gaussianIdx The property channel index in which to store the Gaussian
* curvature. An index < 0 will be ignored (stores nothing). The property set
* will be automatically expanded to include the channel index specified.
*
* @param meanIdx The property channel index in which to store the mean
* curvature. An index < 0 will be ignored (stores nothing). The property set
* will be automatically expanded to include the channel index specified.
*/
Manifold Manifold::CalculateCurvature(int gaussianIdx, int meanIdx) const {
auto pImpl = std::make_shared<Impl>(*GetCsgLeafNode().GetImpl());
pImpl->CalculateCurvature(gaussianIdx, meanIdx);
return Manifold(std::make_shared<CsgLeafNode>(pImpl));
}
/**
* Fills in vertex properties for normal vectors, calculated from the mesh
* geometry. Flat faces composed of three or more triangles will remain flat.
*
* @param normalIdx The property channel in which to store the X
* values of the normals. The X, Y, and Z channels will be sequential. The
* property set will be automatically expanded such that NumProp will be at
* least normalIdx + 3.
*
* @param minSharpAngle Any edges with angles greater than this value will
* remain sharp, getting different normal vector properties on each side of the
* edge. By default, no edges are sharp and all normals are shared. With a value
* of zero, the model is faceted and all normals match their triangle normals,
* but in this case it would be better not to calculate normals at all.
*/
Manifold Manifold::CalculateNormals(int normalIdx, double minSharpAngle) const {
auto pImpl = std::make_shared<Impl>(*GetCsgLeafNode().GetImpl());
pImpl->SetNormals(normalIdx, minSharpAngle);
return Manifold(std::make_shared<CsgLeafNode>(pImpl));
}
/**
* Smooths out the Manifold by filling in the halfedgeTangent vectors. The
* geometry will remain unchanged until Refine or RefineToLength is called to
* interpolate the surface. This version uses the supplied vertex normal
* properties to define the tangent vectors. Faces of two coplanar triangles
* will be marked as quads, while faces with three or more will be flat.
*
* @param normalIdx The first property channel of the normals. NumProp must be
* at least normalIdx + 3. Any vertex where multiple normals exist and don't
* agree will result in a sharp edge.
*/
Manifold Manifold::SmoothByNormals(int normalIdx) const {
auto pImpl = std::make_shared<Impl>(*GetCsgLeafNode().GetImpl());
if (!IsEmpty()) {
pImpl->CreateTangents(normalIdx);
}
return Manifold(std::make_shared<CsgLeafNode>(pImpl));
}
/**
* Smooths out the Manifold by filling in the halfedgeTangent vectors. The
* geometry will remain unchanged until Refine or RefineToLength is called to
* interpolate the surface. This version uses the geometry of the triangles and
* pseudo-normals to define the tangent vectors. Faces of two coplanar triangles
* will be marked as quads.
*
* @param minSharpAngle degrees, default 60. Any edges with angles greater than
* this value will remain sharp. The rest will be smoothed to G1 continuity,
* with the caveat that flat faces of three or more triangles will always remain
* flat. With a value of zero, the model is faceted, but in this case there is
* no point in smoothing.
*
* @param minSmoothness range: 0 - 1, default 0. The smoothness applied to sharp
* angles. The default gives a hard edge, while values > 0 will give a small
* fillet on these sharp edges. A value of 1 is equivalent to a minSharpAngle of
* 180 - all edges will be smooth.
*/
Manifold Manifold::SmoothOut(double minSharpAngle, double minSmoothness) const {
auto pImpl = std::make_shared<Impl>(*GetCsgLeafNode().GetImpl());
if (!IsEmpty()) {
if (minSmoothness == 0) {
const int numProp = pImpl->numProp_;
Vec<double> properties = pImpl->properties_;
Vec<Halfedge> halfedge = pImpl->halfedge_;
pImpl->SetNormals(0, minSharpAngle);
pImpl->CreateTangents(0);
// Reset the properties to the original values, removing temporary normals
pImpl->numProp_ = numProp;
pImpl->properties_.swap(properties);
pImpl->halfedge_.swap(halfedge);
} else {
pImpl->CreateTangents(pImpl->SharpenEdges(minSharpAngle, minSmoothness));
}
}
return Manifold(std::make_shared<CsgLeafNode>(pImpl));
}
/**
* Increase the density of the mesh by splitting every edge into n pieces. For
* instance, with n = 2, each triangle will be split into 4 triangles. Quads
* will ignore their interior triangle bisector. These will all be coplanar (and
* will not be immediately collapsed) unless the Mesh/Manifold has
* halfedgeTangents specified (e.g. from the Smooth() constructor), in which
* case the new vertices will be moved to the interpolated surface according to
* their barycentric coordinates.
*
* @param n The number of pieces to split every edge into. Must be > 1.
*/
Manifold Manifold::Refine(int n) const {
auto pImpl = std::make_shared<Impl>(*GetCsgLeafNode().GetImpl());
if (n > 1) {
pImpl->Refine([n](vec3, vec4, vec4) { return n - 1; });
}
return Manifold(std::make_shared<CsgLeafNode>(pImpl));
}
/**
* Increase the density of the mesh by splitting each edge into pieces of
* roughly the input length. Interior verts are added to keep the rest of the
* triangulation edges also of roughly the same length. If halfedgeTangents are
* present (e.g. from the Smooth() constructor), the new vertices will be moved
* to the interpolated surface according to their barycentric coordinates. Quads
* will ignore their interior triangle bisector.
*
* @param length The length that edges will be broken down to.
*/
Manifold Manifold::RefineToLength(double length) const {
length = std::abs(length);
auto pImpl = std::make_shared<Impl>(*GetCsgLeafNode().GetImpl());
pImpl->Refine([length](vec3 edge, vec4, vec4) {
return static_cast<int>(la::length(edge) / length);
});
return Manifold(std::make_shared<CsgLeafNode>(pImpl));
}
/**
* Increase the density of the mesh by splitting each edge into pieces such that
* any point on the resulting triangles is roughly within tolerance of the
* smoothly curved surface defined by the tangent vectors. This means tightly
* curving regions will be divided more finely than smoother regions. If
* halfedgeTangents are not present, the result will simply be a copy of the
* original. Quads will ignore their interior triangle bisector.
*
* @param tolerance The desired maximum distance between the faceted mesh
* produced and the exact smoothly curving surface. All vertices are exactly on
* the surface, within rounding error.
*/
Manifold Manifold::RefineToTolerance(double tolerance) const {
tolerance = std::abs(tolerance);
auto pImpl = std::make_shared<Impl>(*GetCsgLeafNode().GetImpl());
if (!pImpl->halfedgeTangent_.empty()) {
pImpl->Refine(
[tolerance](vec3 edge, vec4 tangentStart, vec4 tangentEnd) {
const vec3 edgeNorm = la::normalize(edge);
// Weight heuristic
const vec3 tStart = vec3(tangentStart);
const vec3 tEnd = vec3(tangentEnd);
// Perpendicular to edge
const vec3 start = tStart - edgeNorm * la::dot(edgeNorm, tStart);
const vec3 end = tEnd - edgeNorm * la::dot(edgeNorm, tEnd);
// Circular arc result plus heuristic term for non-circular curves
const double d = 0.5 * (la::length(start) + la::length(end)) +
la::length(start - end);
return static_cast<int>(std::sqrt(3 * d / (4 * tolerance)));
},
true);
}
return Manifold(std::make_shared<CsgLeafNode>(pImpl));
}
/**
* The central operation of this library: the Boolean combines two manifolds
* into another by calculating their intersections and removing the unused
* portions.
* [ε-valid](https://github.com/elalish/manifold/wiki/Manifold-Library#definition-of-%CE%B5-valid)
* inputs will produce ε-valid output. ε-invalid input may fail
* triangulation.
*
* These operations are optimized to produce nearly-instant results if either
* input is empty or their bounding boxes do not overlap.
*
* @param second The other Manifold.
* @param op The type of operation to perform.
*/
Manifold Manifold::Boolean(const Manifold& second, OpType op) const {
return Manifold(pNode_->Boolean(second.pNode_, op));
}
/**
* Perform the given boolean operation on a list of Manifolds. In case of
* Subtract, all Manifolds in the tail are differenced from the head.
*/
Manifold Manifold::BatchBoolean(const std::vector<Manifold>& manifolds,
OpType op) {
if (manifolds.size() == 0)
return Manifold();
else if (manifolds.size() == 1)
return manifolds[0];
std::vector<std::shared_ptr<CsgNode>> children;
children.reserve(manifolds.size());
for (const auto& m : manifolds) children.push_back(m.pNode_);
return Manifold(std::make_shared<CsgOpNode>(children, op));
}
/**
* Shorthand for Boolean Union.
*/
Manifold Manifold::operator+(const Manifold& Q) const {
return Boolean(Q, OpType::Add);
}
/**
* Shorthand for Boolean Union assignment.
*/
Manifold& Manifold::operator+=(const Manifold& Q) {
*this = *this + Q;
return *this;
}
/**
* Shorthand for Boolean Difference.
*/
Manifold Manifold::operator-(const Manifold& Q) const {
return Boolean(Q, OpType::Subtract);
}
/**
* Shorthand for Boolean Difference assignment.
*/
Manifold& Manifold::operator-=(const Manifold& Q) {
*this = *this - Q;
return *this;
}
/**
* Shorthand for Boolean Intersection.
*/
Manifold Manifold::operator^(const Manifold& Q) const {
return Boolean(Q, OpType::Intersect);
}
/**
* Shorthand for Boolean Intersection assignment.
*/
Manifold& Manifold::operator^=(const Manifold& Q) {
*this = *this ^ Q;
return *this;
}
/**
* Split cuts this manifold in two using the cutter manifold. The first result
* is the intersection, second is the difference. This is more efficient than
* doing them separately.
*
* @param cutter
*/
std::pair<Manifold, Manifold> Manifold::Split(const Manifold& cutter) const {
auto impl1 = GetCsgLeafNode().GetImpl();
auto impl2 = cutter.GetCsgLeafNode().GetImpl();
Boolean3 boolean(*impl1, *impl2, OpType::Subtract);
auto result1 = std::make_shared<CsgLeafNode>(
std::make_unique<Impl>(boolean.Result(OpType::Intersect)));
auto result2 = std::make_shared<CsgLeafNode>(
std::make_unique<Impl>(boolean.Result(OpType::Subtract)));
return std::make_pair(Manifold(result1), Manifold(result2));
}
/**
* Convenient version of Split() for a half-space.
*
* @param normal This vector is normal to the cutting plane and its length does
* not matter. The first result is in the direction of this vector, the second
* result is on the opposite side.
* @param originOffset The distance of the plane from the origin in the
* direction of the normal vector.
*/
std::pair<Manifold, Manifold> Manifold::SplitByPlane(
vec3 normal, double originOffset) const {
return Split(Halfspace(BoundingBox(), normal, originOffset));
}
/**
* Identical to SplitByPlane(), but calculating and returning only the first
* result.
*
* @param normal This vector is normal to the cutting plane and its length does
* not matter. The result is in the direction of this vector from the plane.
* @param originOffset The distance of the plane from the origin in the
* direction of the normal vector.
*/
Manifold Manifold::TrimByPlane(vec3 normal, double originOffset) const {
return *this ^ Halfspace(BoundingBox(), normal, originOffset);
}
/**
* Returns the cross section of this object parallel to the X-Y plane at the
* specified Z height, defaulting to zero. Using a height equal to the bottom of
* the bounding box will return the bottom faces, while using a height equal to
* the top of the bounding box will return empty.
*/
Polygons Manifold::Slice(double height) const {
return GetCsgLeafNode().GetImpl()->Slice(height);
}
/**
* Returns polygons representing the projected outline of this object
* onto the X-Y plane. These polygons will often self-intersect, so it is
* recommended to run them through the positive fill rule of CrossSection to get
* a sensible result before using them.
*/
Polygons Manifold::Project() const {
return GetCsgLeafNode().GetImpl()->Project();
}
ExecutionParams& ManifoldParams() { return manifoldParams; }
/**
* Compute the convex hull of a set of points. If the given points are fewer
* than 4, or they are all coplanar, an empty Manifold will be returned.
*
* @param pts A vector of 3-dimensional points over which to compute a convex
* hull.
*/
Manifold Manifold::Hull(const std::vector<vec3>& pts) {
std::shared_ptr<Impl> impl = std::make_shared<Impl>();
impl->Hull(Vec<vec3>(pts));
return Manifold(std::make_shared<CsgLeafNode>(impl));
}
/**
* Compute the convex hull of this manifold.
*/
Manifold Manifold::Hull() const {
std::shared_ptr<Impl> impl = std::make_shared<Impl>();
impl->Hull(GetCsgLeafNode().GetImpl()->vertPos_);
return Manifold(std::make_shared<CsgLeafNode>(impl));
}
/**
* Compute the convex hull enveloping a set of manifolds.
*
* @param manifolds A vector of manifolds over which to compute a convex hull.
*/
Manifold Manifold::Hull(const std::vector<Manifold>& manifolds) {
return Compose(manifolds).Hull();
}
/**
* Returns the minimum gap between two manifolds. Returns a double between
* 0 and searchLength.
*
* @param other The other manifold to compute the minimum gap to.
* @param searchLength The maximum distance to search for a minimum gap.
*/
double Manifold::MinGap(const Manifold& other, double searchLength) const {
auto intersect = *this ^ other;
if (!intersect.IsEmpty()) return 0.0;
return GetCsgLeafNode().GetImpl()->MinGap(*other.GetCsgLeafNode().GetImpl(),
searchLength);
}
} // namespace manifold