BRL-CAD/manifold
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

misc changes (#1190)

Browse Source 
* removed use of acos and fma in boolean

* add mesh debug io function

* format

* add comment

* remove another fma
This commit is contained in:
pca006132
2025-03-15 17:00:06 +08:00
committed by GitHub
parent fab833e7e1
commit beb9ca6122
7 changed files with 103 additions and 85 deletions
Download Patch File Download Diff File Expand all files Collapse all files
include/manifold/manifold.h
+4 -4
View File
@@ -437,14 +437,14 @@ class Manifold {
bool MatchesTriNormals() const;
size_t NumDegenerateTris() const;
double GetEpsilon() const;
#ifdef MANIFOLD_EXPORT
static Manifold ImportMeshGL64(std::istream& stream);
std::ostream& Dump(std::ostream& stream) const;
#endif
///@}
struct Impl;
#ifdef MANIFOLD_EXPORT
static Manifold ImportMeshGL64(std::istream& stream);
#endif
private:
Manifold(std::shared_ptr<CsgNode> pNode_);
Manifold(std::shared_ptr<Impl> pImpl_);
src/boolean3.cpp
+7 -7
View File
@@ -42,8 +42,8 @@ vec2 Interpolate(vec3 pL, vec3 pR, double x) {
if (!std::isfinite(lambda) || !std::isfinite(dLR.y) || !std::isfinite(dLR.z))
return vec2(pL.y, pL.z);
vec2 yz;
yz[0] = fma(lambda, dLR.y, useL ? pL.y : pR.y);
yz[1] = fma(lambda, dLR.z, useL ? pL.z : pR.z);
yz[0] = lambda * dLR.y + (useL ? pL.y : pR.y);
yz[1] = lambda * dLR.z + (useL ? pL.z : pR.z);
return yz;
}
@@ -57,14 +57,14 @@ vec4 Intersect(const vec3 &pL, const vec3 &pR, const vec3 &qL, const vec3 &qR) {
double lambda = (useL ? dyL : dyR) / (dyL - dyR);
if (!std::isfinite(lambda)) lambda = 0.0;
vec4 xyzz;
xyzz.x = fma(lambda, dx, useL ? pL.x : pR.x);
xyzz.x = lambda * dx + (useL ? pL.x : pR.x);
const double pDy = pR.y - pL.y;
const double qDy = qR.y - qL.y;
const bool useP = fabs(pDy) < fabs(qDy);
xyzz.y = fma(lambda, useP ? pDy : qDy,
useL ? (useP ? pL.y : qL.y) : (useP ? pR.y : qR.y));
xyzz.z = fma(lambda, pR.z - pL.z, useL ? pL.z : pR.z);
xyzz.w = fma(lambda, qR.z - qL.z, useL ? qL.z : qR.z);
xyzz.y = lambda * (useP ? pDy : qDy) +
(useL ? (useP ? pL.y : qL.y) : (useP ? pR.y : qR.y));
xyzz.z = lambda * (pR.z - pL.z) + (useL ? pL.z : pR.z);
xyzz.w = lambda * (qR.z - qL.z) + (useL ? qL.z : qR.z);
return xyzz;
}
src/csg_tree.cpp
+2
View File
@@ -119,6 +119,7 @@ std::shared_ptr<CsgLeafNode> SimpleBoolean(const Manifold::Impl &a,
<< std::endl;
std::cout << "RHS self-intersecting: " << b.IsSelfIntersecting()
<< std::endl;
#ifdef MANIFOLD_EXPORT
if (ManifoldParams().verbose) {
if (op == OpType::Add)
std::cout << "Add";
@@ -130,6 +131,7 @@ std::shared_ptr<CsgLeafNode> SimpleBoolean(const Manifold::Impl &a,
std::cout << a;
std::cout << b;
}
#endif
dump_lock.unlock();
};
try {
src/impl.cpp
+35 -15
View File
@@ -19,6 +19,7 @@
#include <map>
#include <optional>
#include "./csg_tree.h"
#include "./hashtable.h"
#include "./mesh_fixes.h"
#include "./parallel.h"
@@ -27,6 +28,7 @@
#ifdef MANIFOLD_EXPORT
#include <string.h>
#include <iomanip>
#include <iostream>
#endif
@@ -35,15 +37,20 @@ 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)) {
}
}
// Absolute error <= 6.7e-5
float acos(float x) {
float negate = float(x < 0);
x = abs(x);
float ret = -0.0187293;
ret = ret * x;
ret = ret + 0.0742610;
ret = ret * x;
ret = ret - 0.2121144;
ret = ret * x;
ret = ret + 1.5707288;
ret = ret * sqrt(1.0 - x);
ret = ret - 2 * negate * ret;
return negate * 3.14159265358979 + ret;
}
struct Transform4x3 {
@@ -561,7 +568,7 @@ void Manifold::Impl::CalculateNormals() {
// 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 : std::acos(dot));
double phi = dot >= 1 ? 0 : (dot <= -1 ? kPi : acos(dot));
normal += phi * faceNormal_[edge / 3];
});
vertNormal_[vert] = SafeNormalize(normal);
@@ -589,9 +596,9 @@ void Manifold::Impl::IncrementMeshIDs() {
UpdateMeshID({meshIDold2new.D()}));
}
#ifdef MANIFOLD_DEBUG
#ifdef MANIFOLD_EXPORT
std::ostream& operator<<(std::ostream& stream, const Manifold::Impl& impl) {
stream << std::setprecision(17); // for double precision
stream << std::setprecision(19); // for double precision
stream << "# ======= begin mesh ======" << std::endl;
stream << "# tolerance = " << impl.tolerance_ << std::endl;
stream << "# epsilon = " << impl.epsilon_ << std::endl;
@@ -610,13 +617,26 @@ std::ostream& operator<<(std::ostream& stream, const Manifold::Impl& impl) {
stream << "# ======== end mesh =======" << std::endl;
return stream;
}
#endif
#ifdef MANIFOLD_EXPORT
/**
* 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.
* Should be used with ImportMeshGL64 for reproducing issues.
*/
std::ostream& Manifold::Dump(std::ostream& stream) const {
return stream << *GetCsgLeafNode().GetImpl();
}
/**
* Import a mesh from a Wavefront OBJ file that was exported with Dump.
* This function is the counterpart to Dump and should be used with it.
* This function cannot import OBJ files not written by the Dump function.
*/
Manifold Manifold::ImportMeshGL64(std::istream& stream) {
MeshGL64 mesh;
std::optional<double> epsilon;
stream.precision(17);
stream >> std::setprecision(19);
while (true) {
char c = stream.get();
if (stream.eof()) break;
src/impl.h
+2 -2
View File
@@ -49,6 +49,8 @@ struct Manifold::Impl {
Error status_ = Error::NoError;
Vec<vec3> vertPos_;
Vec<Halfedge> halfedge_;
// Note that vertNormal_ is not precise due to the use of an approximated acos
// function
Vec<vec3> vertNormal_;
Vec<vec3> faceNormal_;
Vec<vec4> halfedgeTangent_;
@@ -376,8 +378,6 @@ struct Manifold::Impl {
void Hull(VecView<vec3> vertPos);
};
#ifdef MANIFOLD_DEBUG
extern std::mutex dump_lock;
std::ostream& operator<<(std::ostream& stream, const Manifold::Impl& impl);
#endif
} // namespace manifold
src/svd.h
+52 -56
View File
@@ -101,29 +101,26 @@ inline Givens ApproximateGivensQuaternion(Symmetric3x3& A) {
inline void JacobiConjugation(const int32_t x, const int32_t y, const int32_t z,
Symmetric3x3& S, vec4& q) {
auto g = ApproximateGivensQuaternion(S);
double scale = 1.0 / fma(g.ch, g.ch, g.sh * g.sh);
double a = fma(g.ch, g.ch, -g.sh * g.sh) * scale;
double scale = 1.0 / (g.ch * g.ch + g.sh * g.sh);
double a = (g.ch * g.ch - g.sh * g.sh) * scale;
double b = 2.0 * g.sh * g.ch * scale;
Symmetric3x3 _S = S;
// perform conjugation S = Q'*S*Q
S.m_00 =
fma(a, fma(a, _S.m_00, b * _S.m_10), b * (fma(a, _S.m_10, b * _S.m_11)));
S.m_10 = fma(a, fma(-b, _S.m_00, a * _S.m_10),
b * (fma(-b, _S.m_10, a * _S.m_11)));
S.m_11 = fma(-b, fma(-b, _S.m_00, a * _S.m_10),
a * (fma(-b, _S.m_10, a * _S.m_11)));
S.m_20 = fma(a, _S.m_20, b * _S.m_21);
S.m_21 = fma(-b, _S.m_20, a * _S.m_21);
S.m_00 = a * (a * _S.m_00 + b * _S.m_10) + b * (a * _S.m_10 + b * _S.m_11);
S.m_10 = a * (-b * _S.m_00 + a * _S.m_10) + b * (-b * _S.m_10 + a * _S.m_11);
S.m_11 = -b * (-b * _S.m_00 + a * _S.m_10) + a * (-b * _S.m_10 + a * _S.m_11);
S.m_20 = a * _S.m_20 + b * _S.m_21;
S.m_21 = -b * _S.m_20 + a * _S.m_21;
S.m_22 = _S.m_22;
// update cumulative rotation qV
vec3 tmp = g.sh * vec3(q);
g.sh *= q[3];
// (x,y,z) corresponds to ((0,1,2),(1,2,0),(2,0,1)) for (p,q) =
// ((0,1),(1,2),(0,2))
q[z] = fma(q[z], g.ch, g.sh);
q[3] = fma(q[3], g.ch, -tmp[z]); // w
q[x] = fma(q[x], g.ch, tmp[y]);
q[y] = fma(q[y], g.ch, -tmp[x]);
q[z] = q[z] * g.ch + g.sh;
q[3] = q[3] * g.ch + -tmp[z]; // w
q[x] = q[x] * g.ch + tmp[y];
q[y] = q[y] * g.ch + -tmp[x];
// re-arrange matrix for next iteration
_S.m_00 = S.m_11;
_S.m_10 = S.m_21;
@@ -148,15 +145,15 @@ inline mat3 JacobiEigenAnalysis(Symmetric3x3 S) {
JacobiConjugation(1, 2, 0, S, q);
JacobiConjugation(2, 0, 1, S, q);
}
return mat3({1.0 - 2.0 * (fma(q.y, q.y, q.z * q.z)), //
2.0 * fma(q.x, q.y, +q.w * q.z), //
2.0 * fma(q.x, q.z, -q.w * q.y)}, //
{2 * fma(q.x, q.y, -q.w * q.z), //
1 - 2 * fma(q.x, q.x, q.z * q.z), //
2 * fma(q.y, q.z, q.w * q.x)}, //
{2 * fma(q.x, q.z, q.w * q.y), //
2 * fma(q.y, q.z, -q.w * q.x), //
1 - 2 * fma(q.x, q.x, q.y * q.y)});
return mat3({1.0 - 2.0 * (q.y * q.y + q.z * q.z), //
2.0 * (q.x * q.y + +q.w * q.z), //
2.0 * (q.x * q.z + -q.w * q.y)}, //
{2 * (q.x * q.y + -q.w * q.z), //
1 - 2 * (q.x * q.x + q.z * q.z), //
2 * (q.y * q.z + q.w * q.x)}, //
{2 * (q.x * q.z + q.w * q.y), //
2 * (q.y * q.z + -q.w * q.x), //
1 - 2 * (q.x * q.x + q.y * q.y)});
}
// Implementation of Algorithm 3
inline void SortSingularValues(mat3& B, mat3& V) {
@@ -207,65 +204,64 @@ inline QR QRDecomposition(mat3& B) {
mat3 Q, R;
// first Givens rotation (ch,0,0,sh)
auto g1 = QRGivensQuaternion(B[0][0], B[0][1]);
auto a = fma(-2.0, g1.sh * g1.sh, 1.0);
auto a = -2.0 * g1.sh * g1.sh + 1.0;
auto b = 2.0 * g1.ch * g1.sh;
// apply B = Q' * B
R[0][0] = fma(a, B[0][0], b * B[0][1]);
R[1][0] = fma(a, B[1][0], b * B[1][1]);
R[2][0] = fma(a, B[2][0], b * B[2][1]);
R[0][1] = fma(-b, B[0][0], a * B[0][1]);
R[1][1] = fma(-b, B[1][0], a * B[1][1]);
R[2][1] = fma(-b, B[2][0], a * B[2][1]);
R[0][0] = a * B[0][0] + b * B[0][1];
R[1][0] = a * B[1][0] + b * B[1][1];
R[2][0] = a * B[2][0] + b * B[2][1];
R[0][1] = -b * B[0][0] + a * B[0][1];
R[1][1] = -b * B[1][0] + a * B[1][1];
R[2][1] = -b * B[2][0] + a * B[2][1];
R[0][2] = B[0][2];
R[1][2] = B[1][2];
R[2][2] = B[2][2];
// second Givens rotation (ch,0,-sh,0)
auto g2 = QRGivensQuaternion(R[0][0], R[0][2]);
a = fma(-2.0, g2.sh * g2.sh, 1.0);
a = -2.0 * g2.sh * g2.sh + 1.0;
b = 2.0 * g2.ch * g2.sh;
// apply B = Q' * B;
B[0][0] = fma(a, R[0][0], b * R[0][2]);
B[1][0] = fma(a, R[1][0], b * R[1][2]);
B[2][0] = fma(a, R[2][0], b * R[2][2]);
B[0][0] = a * R[0][0] + b * R[0][2];
B[1][0] = a * R[1][0] + b * R[1][2];
B[2][0] = a * R[2][0] + b * R[2][2];
B[0][1] = R[0][1];
B[1][1] = R[1][1];
B[2][1] = R[2][1];
B[0][2] = fma(-b, R[0][0], a * R[0][2]);
B[1][2] = fma(-b, R[1][0], a * R[1][2]);
B[2][2] = fma(-b, R[2][0], a * R[2][2]);
B[0][2] = -b * R[0][0] + a * R[0][2];
B[1][2] = -b * R[1][0] + a * R[1][2];
B[2][2] = -b * R[2][0] + a * R[2][2];
// third Givens rotation (ch,sh,0,0)
auto g3 = QRGivensQuaternion(B[1][1], B[1][2]);
a = fma(-2.0, g3.sh * g3.sh, 1.0);
a = -2.0 * g3.sh * g3.sh + 1.0;
b = 2.0 * g3.ch * g3.sh;
// R is now set to desired value
R[0][0] = B[0][0];
R[1][0] = B[1][0];
R[2][0] = B[2][0];
R[0][1] = fma(a, B[0][1], b * B[0][2]);
R[1][1] = fma(a, B[1][1], b * B[1][2]);
R[2][1] = fma(a, B[2][1], b * B[2][2]);
R[0][2] = fma(-b, B[0][1], a * B[0][2]);
R[1][2] = fma(-b, B[1][1], a * B[1][2]);
R[2][2] = fma(-b, B[2][1], a * B[2][2]);
R[0][1] = a * B[0][1] + b * B[0][2];
R[1][1] = a * B[1][1] + b * B[1][2];
R[2][1] = a * B[2][1] + b * B[2][2];
R[0][2] = -b * B[0][1] + a * B[0][2];
R[1][2] = -b * B[1][1] + a * B[1][2];
R[2][2] = -b * B[2][1] + a * B[2][2];
// construct the cumulative rotation Q=Q1 * Q2 * Q3
// the number of floating point operations for three quaternion
// multiplications is more or less comparable to the explicit form of the
// joined matrix. certainly more memory-efficient!
auto sh12 = 2.0 * fma(g1.sh, g1.sh, -0.5);
auto sh22 = 2.0 * fma(g2.sh, g2.sh, -0.5);
auto sh32 = 2.0 * fma(g3.sh, g3.sh, -0.5);
auto sh12 = 2.0 * (g1.sh * g1.sh + -0.5);
auto sh22 = 2.0 * (g2.sh * g2.sh + -0.5);
auto sh32 = 2.0 * (g3.sh * g3.sh + -0.5);
Q[0][0] = sh12 * sh22;
Q[1][0] = fma(4.0 * g2.ch * g3.ch, sh12 * g2.sh * g3.sh,
2.0 * g1.ch * g1.sh * sh32);
Q[2][0] = fma(4.0 * g1.ch * g3.ch, g1.sh * g3.sh,
-2.0 * g2.ch * sh12 * g2.sh * sh32);
Q[1][0] =
4.0 * g2.ch * g3.ch * sh12 * g2.sh * g3.sh + 2.0 * g1.ch * g1.sh * sh32;
Q[2][0] =
4.0 * g1.ch * g3.ch * g1.sh * g3.sh + -2.0 * g2.ch * sh12 * g2.sh * sh32;
Q[0][1] = -2.0 * g1.ch * g1.sh * sh22;
Q[1][1] =
fma(-8.0 * g1.ch * g2.ch * g3.ch, g1.sh * g2.sh * g3.sh, sh12 * sh32);
Q[2][1] = fma(
-2.0 * g3.ch, g3.sh,
4.0 * g1.sh * fma(g3.ch * g1.sh, g3.sh, g1.ch * g2.ch * g2.sh * sh32));
Q[1][1] = -8.0 * g1.ch * g2.ch * g3.ch * g1.sh * g2.sh * g3.sh + sh12 * sh32;
Q[2][1] =
-2.0 * g3.ch * g3.sh +
4.0 * g1.sh * (g3.ch * g1.sh * g3.sh + g1.ch * g2.ch * g2.sh * sh32);
Q[0][2] = 2.0 * g2.ch * g2.sh;
Q[1][2] = -2.0 * g3.ch * sh22 * g3.sh;
src/utils.h
+1 -1
View File
@@ -217,7 +217,7 @@ struct Negate {
inline int CCW(vec2 p0, vec2 p1, vec2 p2, double tol) {
vec2 v1 = p1 - p0;
vec2 v2 = p2 - p0;
double area = fma(v1.x, v2.y, -v1.y * v2.x);
double area = v1.x * v2.y - v1.y * v2.x;
double base2 = la::max(la::dot(v1, v1), la::dot(v2, v2));
if (area * area * 4 <= base2 * tol * tol)
return 0;