diff options
| author | Howard Trickey <howard.trickey@gmail.com> | 2020-08-28 17:56:44 +0300 |
|---|---|---|
| committer | Howard Trickey <howard.trickey@gmail.com> | 2020-08-28 18:01:06 +0300 |
| commit | 9e09b5c418c0a436e3c84ccf38c065527988b0a0 (patch) | |
| tree | c0e81b8834aaf27c22a2734e364452fa2af5c6c1 /source/blender/blenlib/intern/mesh_intersect.cc | |
| parent | 4a17508c6d2a24dfb7c018ae49c80f12b4d3e610 (diff) | |
Merge newboolean branch into master.
This is for design task T67744, Boolean Redesign.
It adds a choice of solver to the Boolean modifier and the
Intersect (Boolean) and Intersect (Knife) tools.
The 'Fast' choice is the current Bmesh boolean.
The new 'Exact' choice is a more advanced algorithm that supports
overlapping geometry and uses more robust calculations, but is
slower than the Fast choice.
The default with this commit is set to 'Exact'. We can decide before
the 2.91 release whether or not this is the right choice, but this
choice now will get us more testing and feedback on the new code.
Diffstat (limited to 'source/blender/blenlib/intern/mesh_intersect.cc')
| -rw-r--r-- | source/blender/blenlib/intern/mesh_intersect.cc | 3304 |
1 files changed, 3304 insertions, 0 deletions
diff --git a/source/blender/blenlib/intern/mesh_intersect.cc b/source/blender/blenlib/intern/mesh_intersect.cc new file mode 100644 index 00000000000..d973775a797 --- /dev/null +++ b/source/blender/blenlib/intern/mesh_intersect.cc @@ -0,0 +1,3304 @@ +/* + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/** \file + * \ingroup bli + */ + +/* The #blender::meshintersect API needs GMP. */ +#ifdef WITH_GMP + +# include <algorithm> +# include <fstream> +# include <iostream> + +# include "BLI_allocator.hh" +# include "BLI_array.hh" +# include "BLI_assert.h" +# include "BLI_delaunay_2d.h" +# include "BLI_double3.hh" +# include "BLI_float3.hh" +# include "BLI_hash.hh" +# include "BLI_kdopbvh.h" +# include "BLI_map.hh" +# include "BLI_math_boolean.hh" +# include "BLI_math_mpq.hh" +# include "BLI_mpq2.hh" +# include "BLI_mpq3.hh" +# include "BLI_span.hh" +# include "BLI_task.h" +# include "BLI_threads.h" +# include "BLI_vector.hh" +# include "BLI_vector_set.hh" + +# include "PIL_time.h" + +# include "BLI_mesh_intersect.hh" + +// # define PERFDEBUG + +namespace blender::meshintersect { + +# ifdef PERFDEBUG +static void perfdata_init(void); +static void incperfcount(int countnum); +static void bumpperfcount(int countnum, int amt); +static void doperfmax(int maxnum, int val); +static void dump_perfdata(void); +# endif + +/** For debugging, can disable threading in intersect code with this static constant. */ +static constexpr bool intersect_use_threading = true; + +Vert::Vert(const mpq3 &mco, const double3 &dco, int id, int orig) + : co_exact(mco), co(dco), id(id), orig(orig) +{ +} + +bool Vert::operator==(const Vert &other) const +{ + return this->co_exact == other.co_exact; +} + +uint64_t Vert::hash() const +{ + return co_exact.hash(); +} + +std::ostream &operator<<(std::ostream &os, const Vert *v) +{ + os << "v" << v->id; + if (v->orig != NO_INDEX) { + os << "o" << v->orig; + } + os << v->co; + return os; +} + +bool Plane::operator==(const Plane &other) const +{ + return norm_exact == other.norm_exact && d_exact == other.d_exact; +} + +void Plane::make_canonical() +{ + if (norm_exact[0] != 0) { + mpq_class den = norm_exact[0]; + norm_exact = mpq3(1, norm_exact[1] / den, norm_exact[2] / den); + d_exact = d_exact / den; + } + else if (norm_exact[1] != 0) { + mpq_class den = norm_exact[1]; + norm_exact = mpq3(0, 1, norm_exact[2] / den); + d_exact = d_exact / den; + } + else { + if (norm_exact[2] != 0) { + mpq_class den = norm_exact[2]; + norm_exact = mpq3(0, 0, 1); + d_exact = d_exact / den; + } + else { + /* A degenerate plane. */ + d_exact = 0; + } + } + norm = double3(norm_exact[0].get_d(), norm_exact[1].get_d(), norm_exact[2].get_d()); + d = d_exact.get_d(); +} + +Plane::Plane(const mpq3 &norm_exact, const mpq_class &d_exact) + : norm_exact(norm_exact), d_exact(d_exact) +{ + norm = double3(norm_exact[0].get_d(), norm_exact[1].get_d(), norm_exact[2].get_d()); + d = d_exact.get_d(); +} + +Plane::Plane(const double3 &norm, const double d) : norm(norm), d(d) +{ + norm_exact = mpq3(0, 0, 0); /* Marks as "exact not yet populated". */ +} + +/** This is wrong for degenerate planes, but we don't expect to call it on those. */ +bool Plane::exact_populated() const +{ + return norm_exact[0] != 0 || norm_exact[1] != 0 || norm_exact[2] != 0; +} + +uint64_t Plane::hash() const +{ + constexpr uint64_t h1 = 33; + constexpr uint64_t h2 = 37; + constexpr uint64_t h3 = 39; + uint64_t hashx = hash_mpq_class(this->norm_exact.x); + uint64_t hashy = hash_mpq_class(this->norm_exact.y); + uint64_t hashz = hash_mpq_class(this->norm_exact.z); + uint64_t hashd = hash_mpq_class(this->d_exact); + uint64_t ans = hashx ^ (hashy * h1) ^ (hashz * h1 * h2) ^ (hashd * h1 * h2 * h3); + return ans; +} + +std::ostream &operator<<(std::ostream &os, const Plane *plane) +{ + os << "[" << plane->norm << ";" << plane->d << "]"; + return os; +} + +Face::Face( + Span<const Vert *> verts, int id, int orig, Span<int> edge_origs, Span<bool> is_intersect) + : vert(verts), edge_orig(edge_origs), is_intersect(is_intersect), id(id), orig(orig) +{ +} + +Face::Face(Span<const Vert *> verts, int id, int orig) : vert(verts), id(id), orig(orig) +{ +} + +void Face::populate_plane(bool need_exact) +{ + if (plane != nullptr) { + if (!need_exact || plane->exact_populated()) { + return; + } + } + if (need_exact) { + mpq3 normal_exact; + if (vert.size() > 3) { + Array<mpq3> co(vert.size()); + for (int i : index_range()) { + co[i] = vert[i]->co_exact; + } + normal_exact = mpq3::cross_poly(co); + } + else { + mpq3 tr02 = vert[0]->co_exact - vert[2]->co_exact; + mpq3 tr12 = vert[1]->co_exact - vert[2]->co_exact; + normal_exact = mpq3::cross(tr02, tr12); + } + mpq_class d_exact = -mpq3::dot(normal_exact, vert[0]->co_exact); + plane = new Plane(normal_exact, d_exact); + } + else { + double3 normal; + if (vert.size() > 3) { + Array<double3> co(vert.size()); + for (int i : index_range()) { + co[i] = vert[i]->co; + } + normal = double3::cross_poly(co); + } + else { + double3 tr02 = vert[0]->co - vert[2]->co; + double3 tr12 = vert[1]->co - vert[2]->co; + normal = double3::cross_high_precision(tr02, tr12); + } + double d = -double3::dot(normal, vert[0]->co); + plane = new Plane(normal, d); + } +} + +Face::~Face() +{ + if (plane != nullptr) { + delete plane; + } +} + +bool Face::operator==(const Face &other) const +{ + if (this->size() != other.size()) { + return false; + } + for (FacePos i : index_range()) { + /* Can test pointer equality since we will have + * unique vert pointers for unique co_equal's. */ + if (this->vert[i] != other.vert[i]) { + return false; + } + } + return true; +} + +bool Face::cyclic_equal(const Face &other) const +{ + if (this->size() != other.size()) { + return false; + } + int flen = this->size(); + for (FacePos start : index_range()) { + for (FacePos start_other : index_range()) { + bool ok = true; + for (int i = 0; ok && i < flen; ++i) { + FacePos p = (start + i) % flen; + FacePos p_other = (start_other + i) % flen; + if (this->vert[p] != other.vert[p_other]) { + ok = false; + } + } + if (ok) { + return true; + } + } + } + return false; +} + +std::ostream &operator<<(std::ostream &os, const Face *f) +{ + os << "f" << f->id << "o" << f->orig << "["; + for (const Vert *v : *f) { + os << "v" << v->id; + if (v->orig != NO_INDEX) { + os << "o" << v->orig; + } + if (v != f->vert[f->size() - 1]) { + os << " "; + } + } + os << "]"; + if (f->orig != NO_INDEX) { + os << "o" << f->orig; + } + os << " e_orig["; + for (int i : f->index_range()) { + os << f->edge_orig[i]; + if (f->is_intersect[i]) { + os << "#"; + } + if (i != f->size() - 1) { + os << " "; + } + } + os << "]"; + return os; +} + +/** + * Un-comment the following to try using a spin-lock instead of + * a mutex in the arena allocation routines. + * Initial tests showed that it doesn't seem to help very much, + * if at all, to use a spin-lock. + */ +// #define USE_SPINLOCK + +/** + * #IMeshArena is the owner of the Vert and Face resources used + * during a run of one of the mesh-intersect main functions. + * It also keeps has a hash table of all Verts created so that it can + * ensure that only one instance of a Vert with a given co_exact will + * exist. I.e., it de-duplicates the vertices. + */ +class IMeshArena::IMeshArenaImpl : NonCopyable, NonMovable { + + /** + * Don't use Vert itself as key since resizing may move + * pointers to the Vert around, and we need to have those pointers + * stay the same throughout the lifetime of the #IMeshArena. + */ + struct VSetKey { + Vert *vert; + + VSetKey(Vert *p) : vert(p) + { + } + + uint32_t hash() const + { + return vert->hash(); + } + + bool operator==(const VSetKey &other) const + { + return *this->vert == *other.vert; + } + }; + + VectorSet<VSetKey> vset_; /* TODO: replace with Set */ + + /** + * Ownership of the Vert memory is here, so destroying this reclaims that memory. + * + * TODO: replace these with pooled allocation, and just destroy the pools at the end. + */ + Vector<std::unique_ptr<Vert>> allocated_verts_; + Vector<std::unique_ptr<Face>> allocated_faces_; + + /* Use these to allocate ids when Verts and Faces are allocated. */ + int next_vert_id_ = 0; + int next_face_id_ = 0; + + /* Need a lock when multi-threading to protect allocation of new elements. */ +# ifdef USE_SPINLOCK + SpinLock lock_; +# else + ThreadMutex *mutex_; +# endif + + public: + IMeshArenaImpl() + { + if (intersect_use_threading) { +# ifdef USE_SPINLOCK + BLI_spin_init(&lock_); +# else + mutex_ = BLI_mutex_alloc(); +# endif + } + } + ~IMeshArenaImpl() + { + if (intersect_use_threading) { +# ifdef USE_SPINLOCK + BLI_spin_end(&lock_); +# else + BLI_mutex_free(mutex_); +# endif + } + } + + void reserve(int vert_num_hint, int face_num_hint) + { + vset_.reserve(vert_num_hint); + allocated_verts_.reserve(vert_num_hint); + allocated_faces_.reserve(face_num_hint); + } + + int tot_allocated_verts() const + { + return allocated_verts_.size(); + } + + int tot_allocated_faces() const + { + return allocated_faces_.size(); + } + + const Vert *add_or_find_vert(const mpq3 &co, int orig) + { + double3 dco(co[0].get_d(), co[1].get_d(), co[2].get_d()); + return add_or_find_vert(co, dco, orig); + } + + const Vert *add_or_find_vert(const double3 &co, int orig) + { + mpq3 mco(co[0], co[1], co[2]); + return add_or_find_vert(mco, co, orig); + } + + Face *add_face(Span<const Vert *> verts, int orig, Span<int> edge_origs, Span<bool> is_intersect) + { + Face *f = new Face(verts, next_face_id_++, orig, edge_origs, is_intersect); + if (intersect_use_threading) { +# ifdef USE_SPINLOCK + BLI_spin_lock(&lock_); +# else + BLI_mutex_lock(mutex_); +# endif + } + allocated_faces_.append(std::unique_ptr<Face>(f)); + if (intersect_use_threading) { +# ifdef USE_SPINLOCK + BLI_spin_unlock(&lock_); +# else + BLI_mutex_unlock(mutex_); +# endif + } + return f; + } + + Face *add_face(Span<const Vert *> verts, int orig, Span<int> edge_origs) + { + Array<bool> is_intersect(verts.size(), false); + return add_face(verts, orig, edge_origs, is_intersect); + } + + Face *add_face(Span<const Vert *> verts, int orig) + { + Array<int> edge_origs(verts.size(), NO_INDEX); + Array<bool> is_intersect(verts.size(), false); + return add_face(verts, orig, edge_origs, is_intersect); + } + + const Vert *find_vert(const mpq3 &co) + { + const Vert *ans; + Vert vtry(co, double3(), NO_INDEX, NO_INDEX); + VSetKey vskey(&vtry); + if (intersect_use_threading) { +# ifdef USE_SPINLOCK + BLI_spin_lock(&lock_); +# else + BLI_mutex_lock(mutex_); +# endif + } + int i = vset_.index_of_try(vskey); + if (i == -1) { + ans = nullptr; + } + else { + ans = vset_[i].vert; + } + if (intersect_use_threading) { +# ifdef USE_SPINLOCK + BLI_spin_unlock(&lock_); +# else + BLI_mutex_unlock(mutex_); +# endif + } + return ans; + } + + /** + * This is slow. Only used for unit tests right now. + * Since it is only used for that purpose, access is not lock-protected. + * The argument vs can be a cyclic shift of the actual stored Face. + */ + const Face *find_face(Span<const Vert *> vs) + { + Array<int> eorig(vs.size(), NO_INDEX); + Array<bool> is_intersect(vs.size(), false); + Face ftry(vs, NO_INDEX, NO_INDEX, eorig, is_intersect); + for (const int i : allocated_faces_.index_range()) { + if (ftry.cyclic_equal(*allocated_faces_[i])) { + return allocated_faces_[i].get(); + } + } + return nullptr; + } + + private: + const Vert *add_or_find_vert(const mpq3 &mco, const double3 &dco, int orig) + { + /* Don't allocate Vert yet, in case it is already there. */ + Vert vtry(mco, dco, NO_INDEX, NO_INDEX); + const Vert *ans; + VSetKey vskey(&vtry); + if (intersect_use_threading) { +# ifdef USE_SPINLOCK + BLI_spin_lock(&lock_); +# else + BLI_mutex_lock(mutex_); +# endif + } + int i = vset_.index_of_try(vskey); + if (i == -1) { + vskey.vert = new Vert(mco, dco, next_vert_id_++, orig); + vset_.add_new(vskey); + allocated_verts_.append(std::unique_ptr<Vert>(vskey.vert)); + ans = vskey.vert; + } + else { + /* It was a duplicate, so return the existing one. + * Note that the returned Vert may have a different orig. + * This is the intended semantics: if the Vert already + * exists then we are merging verts and using the first-seen + * one as the canonical one. */ + ans = vset_[i].vert; + } + if (intersect_use_threading) { +# ifdef USE_SPINLOCK + BLI_spin_unlock(&lock_); +# else + BLI_mutex_unlock(mutex_); +# endif + } + return ans; + }; +}; + +IMeshArena::IMeshArena() +{ + pimpl_ = std::unique_ptr<IMeshArenaImpl>(new IMeshArenaImpl()); +} + +IMeshArena::~IMeshArena() +{ +} + +void IMeshArena::reserve(int vert_num_hint, int face_num_hint) +{ + pimpl_->reserve(vert_num_hint, face_num_hint); +} + +int IMeshArena::tot_allocated_verts() const +{ + return pimpl_->tot_allocated_verts(); +} + +int IMeshArena::tot_allocated_faces() const +{ + return pimpl_->tot_allocated_faces(); +} + +const Vert *IMeshArena::add_or_find_vert(const mpq3 &co, int orig) +{ + return pimpl_->add_or_find_vert(co, orig); +} + +Face *IMeshArena::add_face(Span<const Vert *> verts, + int orig, + Span<int> edge_origs, + Span<bool> is_intersect) +{ + return pimpl_->add_face(verts, orig, edge_origs, is_intersect); +} + +Face *IMeshArena::add_face(Span<const Vert *> verts, int orig, Span<int> edge_origs) +{ + return pimpl_->add_face(verts, orig, edge_origs); +} + +Face *IMeshArena::add_face(Span<const Vert *> verts, int orig) +{ + return pimpl_->add_face(verts, orig); +} + +const Vert *IMeshArena::add_or_find_vert(const double3 &co, int orig) +{ + return pimpl_->add_or_find_vert(co, orig); +} + +const Vert *IMeshArena::find_vert(const mpq3 &co) const +{ + return pimpl_->find_vert(co); +} + +const Face *IMeshArena::find_face(Span<const Vert *> verts) const +{ + return pimpl_->find_face(verts); +} + +void IMesh::set_faces(Span<Face *> faces) +{ + face_ = faces; +} + +int IMesh::lookup_vert(const Vert *v) const +{ + BLI_assert(vert_populated_); + return vert_to_index_.lookup_default(v, NO_INDEX); +} + +void IMesh::populate_vert() +{ + /* This is likely an overestimate, since verts are shared between + * faces. It is ok if estimate is over or even under. */ + constexpr int ESTIMATE_VERTS_PER_FACE = 4; + int estimate_num_verts = ESTIMATE_VERTS_PER_FACE * face_.size(); + populate_vert(estimate_num_verts); +} + +void IMesh::populate_vert(int max_verts) +{ + if (vert_populated_) { + return; + } + vert_to_index_.reserve(max_verts); + int next_allocate_index = 0; + for (const Face *f : face_) { + for (const Vert *v : *f) { + if (v->id == 1) { + } + int index = vert_to_index_.lookup_default(v, NO_INDEX); + if (index == NO_INDEX) { + BLI_assert(next_allocate_index < UINT_MAX - 2); + vert_to_index_.add(v, next_allocate_index++); + } + } + } + int tot_v = next_allocate_index; + vert_ = Array<const Vert *>(tot_v); + for (auto item : vert_to_index_.items()) { + int index = item.value; + BLI_assert(index < tot_v); + vert_[index] = item.key; + } + /* Easier debugging (at least when there are no merged input verts) + * if output vert order is same as input, with new verts at the end. + * TODO: when all debugged, set fix_order = false. */ + const bool fix_order = true; + if (fix_order) { + std::sort(vert_.begin(), vert_.end(), [](const Vert *a, const Vert *b) { + if (a->orig != NO_INDEX && b->orig != NO_INDEX) { + return a->orig < b->orig; + } + if (a->orig != NO_INDEX) { + return true; + } + if (b->orig != NO_INDEX) { + return false; + } + return a->id < b->id; + }); + for (int i : vert_.index_range()) { + const Vert *v = vert_[i]; + vert_to_index_.add_overwrite(v, i); + } + } + vert_populated_ = true; +} + +void IMesh::erase_face_positions(int f_index, Span<bool> face_pos_erase, IMeshArena *arena) +{ + const Face *cur_f = this->face(f_index); + int cur_len = cur_f->size(); + int num_to_erase = 0; + for (int i : cur_f->index_range()) { + if (face_pos_erase[i]) { + ++num_to_erase; + } + } + if (num_to_erase == 0) { + return; + } + int new_len = cur_len - num_to_erase; + if (new_len < 3) { + /* Invalid erase. Don't do anything. */ + return; + } + Array<const Vert *> new_vert(new_len); + Array<int> new_edge_orig(new_len); + Array<bool> new_is_intersect(new_len); + int new_index = 0; + for (int i : cur_f->index_range()) { + if (!face_pos_erase[i]) { + new_vert[new_index] = (*cur_f)[i]; + new_edge_orig[new_index] = cur_f->edge_orig[i]; + new_is_intersect[new_index] = cur_f->is_intersect[i]; + ++new_index; + } + } + BLI_assert(new_index == new_len); + this->face_[f_index] = arena->add_face(new_vert, cur_f->orig, new_edge_orig, new_is_intersect); +} + +std::ostream &operator<<(std::ostream &os, const IMesh &mesh) +{ + if (mesh.has_verts()) { + os << "Verts:\n"; + int i = 0; + for (const Vert *v : mesh.vertices()) { + os << i << ": " << v << "\n"; + ++i; + } + } + os << "\nFaces:\n"; + int i = 0; + for (const Face *f : mesh.faces()) { + os << i << ": " << f << "\n"; + if (f->plane != nullptr) { + os << " plane=" << f->plane << " eorig=["; + for (Face::FacePos p = 0; p < f->size(); ++p) { + os << f->edge_orig[p] << " "; + } + os << "]\n"; + } + ++i; + } + return os; +} + +struct BoundingBox { + float3 min{FLT_MAX, FLT_MAX, FLT_MAX}; + float3 max{-FLT_MAX, -FLT_MAX, -FLT_MAX}; + + BoundingBox() = default; + BoundingBox(const float3 &min, const float3 &max) : min(min), max(max) + { + } + BoundingBox(const BoundingBox &other) : min(other.min), max(other.max) + { + } + BoundingBox(BoundingBox &&other) noexcept : min(std::move(other.min)), max(std::move(other.max)) + { + } + ~BoundingBox() = default; + BoundingBox operator=(const BoundingBox &other) + { + if (this != &other) { + min = other.min; + max = other.max; + } + return *this; + } + BoundingBox operator=(BoundingBox &&other) noexcept + { + min = std::move(other.min); + max = std::move(other.max); + return *this; + } + + void combine(const float3 &p) + { + min.x = min_ff(min.x, p.x); + min.y = min_ff(min.y, p.y); + min.z = min_ff(min.z, p.z); + max.x = max_ff(max.x, p.x); + max.y = max_ff(max.y, p.y); + max.z = max_ff(max.z, p.z); + } + + void combine(const double3 &p) + { + min.x = min_ff(min.x, static_cast<float>(p.x)); + min.y = min_ff(min.y, static_cast<float>(p.y)); + min.z = min_ff(min.z, static_cast<float>(p.z)); + max.x = max_ff(max.x, static_cast<float>(p.x)); + max.y = max_ff(max.y, static_cast<float>(p.y)); + max.z = max_ff(max.z, static_cast<float>(p.z)); + } + + void combine(const BoundingBox &bb) + { + min.x = min_ff(min.x, bb.min.x); + min.y = min_ff(min.y, bb.min.y); + min.z = min_ff(min.z, bb.min.z); + max.x = max_ff(max.x, bb.max.x); + max.y = max_ff(max.y, bb.max.y); + max.z = max_ff(max.z, bb.max.z); + } + + void expand(float pad) + { + min.x -= pad; + min.y -= pad; + min.z -= pad; + max.x += pad; + max.y += pad; + max.z += pad; + } +}; + +/** + * Assume bounding boxes have been expanded by a sufficient epsilon on all sides + * so that the comparisons against the bb bounds are sufficient to guarantee that + * if an overlap or even touching could happen, this will return true. + */ +static bool bbs_might_intersect(const BoundingBox &bb_a, const BoundingBox &bb_b) +{ + return isect_aabb_aabb_v3(bb_a.min, bb_a.max, bb_b.min, bb_b.max); +} + +/** + * Data and functions to calculate bounding boxes and pad them, in parallel. + * The bounding box calculation has the additional task of calculating the maximum + * absolute value of any coordinate in the mesh, which will be used to calculate + * the pad value. + */ +struct BBChunkData { + double max_abs_val = 0.0; +}; + +struct BBCalcData { + const IMesh &im; + Array<BoundingBox> *face_bounding_box; + + BBCalcData(const IMesh &im, Array<BoundingBox> *fbb) : im(im), face_bounding_box(fbb){}; +}; + +static void calc_face_bb_range_func(void *__restrict userdata, + const int iter, + const TaskParallelTLS *__restrict tls) +{ + BBCalcData *bbdata = static_cast<BBCalcData *>(userdata); + double max_abs = 0.0; + const Face &face = *bbdata->im.face(iter); + BoundingBox &bb = (*bbdata->face_bounding_box)[iter]; + for (const Vert *v : face) { + bb.combine(v->co); + for (int i = 0; i < 3; ++i) { + max_abs = max_dd(max_abs, fabs(v->co[i])); + } + } + BBChunkData *chunk = static_cast<BBChunkData *>(tls->userdata_chunk); + chunk->max_abs_val = max_dd(max_abs, chunk->max_abs_val); +} + +struct BBPadData { + Array<BoundingBox> *face_bounding_box; + double pad; + + BBPadData(Array<BoundingBox> *fbb, double pad) : face_bounding_box(fbb), pad(pad){}; +}; + +static void pad_face_bb_range_func(void *__restrict userdata, + const int iter, + const TaskParallelTLS *__restrict UNUSED(tls)) +{ + BBPadData *pad_data = static_cast<BBPadData *>(userdata); + (*pad_data->face_bounding_box)[iter].expand(pad_data->pad); +} + +static void calc_face_bb_reduce(const void *__restrict UNUSED(userdata), + void *__restrict chunk_join, + void *__restrict chunk) +{ + BBChunkData *bbchunk_join = static_cast<BBChunkData *>(chunk_join); + BBChunkData *bbchunk = static_cast<BBChunkData *>(chunk); + bbchunk_join->max_abs_val = max_dd(bbchunk_join->max_abs_val, bbchunk->max_abs_val); +} + +/** + * We will expand the bounding boxes by an epsilon on all sides so that + * the "less than" tests in isect_aabb_aabb_v3 are sufficient to detect + * touching or overlap. + */ +static Array<BoundingBox> calc_face_bounding_boxes(const IMesh &m) +{ + int n = m.face_size(); + Array<BoundingBox> ans(n); + TaskParallelSettings settings; + BBCalcData data(m, &ans); + BBChunkData chunk_data; + BLI_parallel_range_settings_defaults(&settings); + settings.userdata_chunk = &chunk_data; + settings.userdata_chunk_size = sizeof(chunk_data); + settings.func_reduce = calc_face_bb_reduce; + settings.min_iter_per_thread = 1000; + settings.use_threading = intersect_use_threading; + BLI_task_parallel_range(0, n, &data, calc_face_bb_range_func, &settings); + double max_abs_val = chunk_data.max_abs_val; + constexpr float pad_factor = 10.0f; + float pad = max_abs_val == 0.0f ? FLT_EPSILON : 2 * FLT_EPSILON * max_abs_val; + pad *= pad_factor; /* For extra safety. */ + TaskParallelSettings pad_settings; + BLI_parallel_range_settings_defaults(&pad_settings); + settings.min_iter_per_thread = 1000; + settings.use_threading = intersect_use_threading; + BBPadData pad_data(&ans, pad); + BLI_task_parallel_range(0, n, &pad_data, pad_face_bb_range_func, &pad_settings); + return ans; +} + +/** + * A cluster of co-planar triangles, by index. + * A pair of triangles T0 and T1 is said to "non-trivially co-planar-intersect" + * if they are co-planar, intersect, and their intersection is not just existing + * elements (verts, edges) of both triangles. + * A co-planar cluster is said to be "nontrivial" if it has more than one triangle + * and every triangle in it non-trivially co-planar-intersects with at least one other + * triangle in the cluster. + */ +class CoplanarCluster { + Vector<int> tris_; + BoundingBox bb_; + + public: + CoplanarCluster() = default; + CoplanarCluster(int t, const BoundingBox &bb) + { + this->add_tri(t, bb); + } + CoplanarCluster(const CoplanarCluster &other) : tris_(other.tris_), bb_(other.bb_) + { + } + CoplanarCluster(CoplanarCluster &&other) noexcept + : tris_(std::move(other.tris_)), bb_(std::move(other.bb_)) + { + } + ~CoplanarCluster() = default; + CoplanarCluster &operator=(const CoplanarCluster &other) + { + if (this != &other) { + tris_ = other.tris_; + bb_ = other.bb_; + } + return *this; + } + CoplanarCluster &operator=(CoplanarCluster &&other) noexcept + { + tris_ = std::move(other.tris_); + bb_ = std::move(other.bb_); + return *this; + } + + /* Assume that caller knows this will not be a duplicate. */ + void add_tri(int t, const BoundingBox &bb) + { + tris_.append(t); + bb_ = bb; + } + int tot_tri() const + { + return tris_.size(); + } + int tri(int index) const + { + return tris_[index]; + } + const int *begin() const + { + return tris_.begin(); + } + const int *end() const + { + return tris_.end(); + } + + const BoundingBox &bounding_box() const + { + return bb_; + } +}; + +/** + * Maintains indexed set of #CoplanarCluster, with the added ability + * to efficiently find the cluster index of any given triangle + * (the max triangle index needs to be given in the initializer). + * The #tri_cluster(t) function returns -1 if t is not part of any cluster. + */ +class CoplanarClusterInfo { + Vector<CoplanarCluster> clusters_; + Array<int> tri_cluster_; + + public: + CoplanarClusterInfo() = default; + explicit CoplanarClusterInfo(int numtri) : tri_cluster_(Array<int>(numtri)) + { + tri_cluster_.fill(-1); + } + + int tri_cluster(int t) const + { + BLI_assert(t < tri_cluster_.size()); + return tri_cluster_[t]; + } + + int add_cluster(CoplanarCluster cl) + { + int c_index = clusters_.append_and_get_index(cl); + for (int t : cl) { + BLI_assert(t < tri_cluster_.size()); + tri_cluster_[t] = c_index; + } + return c_index; + } + + int tot_cluster() const + { + return clusters_.size(); + } + + const CoplanarCluster *begin() + { + return clusters_.begin(); + } + + const CoplanarCluster *end() + { + return clusters_.end(); + } + + IndexRange index_range() const + { + return clusters_.index_range(); + } + + const CoplanarCluster &cluster(int index) const + { + BLI_assert(index < clusters_.size()); + return clusters_[index]; + } +}; + +static std::ostream &operator<<(std::ostream &os, const CoplanarCluster &cl); + +static std::ostream &operator<<(std::ostream &os, const CoplanarClusterInfo &clinfo); + +enum ITT_value_kind { INONE, IPOINT, ISEGMENT, ICOPLANAR }; + +struct ITT_value { + enum ITT_value_kind kind; + mpq3 p1; /* Only relevant for IPOINT and ISEGMENT kind. */ + mpq3 p2; /* Only relevant for ISEGMENT kind. */ + int t_source; /* Index of the source triangle that intersected the target one. */ + + ITT_value() : kind(INONE), t_source(-1) + { + } + ITT_value(ITT_value_kind k) : kind(k), t_source(-1) + { + } + ITT_value(ITT_value_kind k, int tsrc) : kind(k), t_source(tsrc) + { + } + ITT_value(ITT_value_kind k, const mpq3 &p1) : kind(k), p1(p1), t_source(-1) + { + } + ITT_value(ITT_value_kind k, const mpq3 &p1, const mpq3 &p2) + : kind(k), p1(p1), p2(p2), t_source(-1) + { + } + ITT_value(const ITT_value &other) + : kind(other.kind), p1(other.p1), p2(other.p2), t_source(other.t_source) + { + } + ITT_value(ITT_value &&other) noexcept + : kind(other.kind), + p1(std::move(other.p1)), + p2(std::move(other.p2)), + t_source(other.t_source) + { + } + ~ITT_value() + { + } + ITT_value &operator=(const ITT_value &other) + { + if (this != &other) { + kind = other.kind; + p1 = other.p1; + p2 = other.p2; + t_source = other.t_source; + } + return *this; + } + ITT_value &operator=(ITT_value &&other) noexcept + { + kind = other.kind; + p1 = std::move(other.p1); + p2 = std::move(other.p2); + t_source = other.t_source; + return *this; + } +}; + +static std::ostream &operator<<(std::ostream &os, const ITT_value &itt); + +/** + * Project a 3d vert to a 2d one by eliding proj_axis. This does not create + * degeneracies as long as the projection axis is one where the corresponding + * component of the originating plane normal is non-zero. + */ +static mpq2 project_3d_to_2d(const mpq3 &p3d, int proj_axis) +{ + mpq2 p2d; + switch (proj_axis) { + case (0): { + p2d[0] = p3d[1]; + p2d[1] = p3d[2]; + break; + } + case (1): { + p2d[0] = p3d[0]; + p2d[1] = p3d[2]; + break; + } + case (2): { + p2d[0] = p3d[0]; + p2d[1] = p3d[1]; + break; + } + default: + BLI_assert(false); + } + return p2d; +} + +/** + Is a point in the interior of a 2d triangle or on one of its + * edges but not either endpoint of the edge? + * orient[pi][i] is the orientation test of the point pi against + * the side of the triangle starting at index i. + * Assume the triangle is non-degenerate and CCW-oriented. + * Then answer is true if p is left of or on all three of triangle a's edges, + * and strictly left of at least on of them. + */ +static bool non_trivially_2d_point_in_tri(const int orients[3][3], int pi) +{ + int p_left_01 = orients[pi][0]; + int p_left_12 = orients[pi][1]; + int p_left_20 = orients[pi][2]; + return (p_left_01 >= 0 && p_left_12 >= 0 && p_left_20 >= 0 && + (p_left_01 + p_left_12 + p_left_20) >= 2); +} + +/** + * Given orients as defined in non_trivially_2d_intersect, do the triangles + * overlap in a "hex" pattern? That is, the overlap region is a hexagon, which + * one gets by having, each point of one triangle being strictly right-of one + * edge of the other and strictly left of the other two edges; and vice versa. + */ +static bool non_trivially_2d_hex_overlap(int orients[2][3][3]) +{ + for (int ab = 0; ab < 2; ++ab) { + for (int i = 0; i < 3; ++i) { + bool ok = orients[ab][i][0] + orients[ab][i][1] + orients[ab][i][2] == 1 && + orients[ab][i][0] != 0 && orients[ab][i][1] != 0 && orients[i][2] != 0; + if (!ok) { + return false; + } + } + } + return true; +} + +/** + * Given orients as defined in non_trivially_2d_intersect, do the triangles + * have one shared edge in a "folded-over" configuration? + * As well as a shared edge, the third vertex of one triangle needs to be + * right-of one and left-of the other two edges of the other triangle. + */ +static bool non_trivially_2d_shared_edge_overlap(int orients[2][3][3], + const mpq2 *a[3], + const mpq2 *b[3]) +{ + for (int i = 0; i < 3; ++i) { + int in = (i + 1) % 3; + int inn = (i + 2) % 3; + for (int j = 0; j < 3; ++j) { + int jn = (j + 1) % 3; + int jnn = (j + 2) % 3; + if (*a[i] == *b[j] && *a[in] == *b[jn]) { + /* Edge from a[i] is shared with edge from b[j]. */ + /* See if a[inn] is right-of or on one of the other edges of b. + * If it is on, then it has to be right-of or left-of the shared edge, + * depending on which edge it is. */ + if (orients[0][inn][jn] < 0 || orients[0][inn][jnn] < 0) { + return true; + } + if (orients[0][inn][jn] == 0 && orients[0][inn][j] == 1) { + return true; + } + if (orients[0][inn][jnn] == 0 && orients[0][inn][j] == -1) { + return true; + } + /* Similarly for `b[jnn]`. */ + if (orients[1][jnn][in] < 0 || orients[1][jnn][inn] < 0) { + return true; + } + if (orients[1][jnn][in] == 0 && orients[1][jnn][i] == 1) { + return true; + } + if (orients[1][jnn][inn] == 0 && orients[1][jnn][i] == -1) { + return true; + } + } + } + } + return false; +} + +/** + * Are the triangles the same, perhaps with some permutation of vertices? + */ +static bool same_triangles(const mpq2 *a[3], const mpq2 *b[3]) +{ + for (int i = 0; i < 3; ++i) { + if (a[0] == b[i] && a[1] == b[(i + 1) % 3] && a[2] == b[(i + 2) % 3]) { + return true; + } + } + return false; +} + +/** + * Do 2d triangles (a[0], a[1], a[2]) and (b[0], b[1], b2[2]) intersect at more than just shared + * vertices or a shared edge? This is true if any point of one triangle is non-trivially inside the + * other. NO: that isn't quite sufficient: there is also the case where the verts are all mutually + * outside the other's triangle, but there is a hexagonal overlap region where they overlap. + */ +static bool non_trivially_2d_intersect(const mpq2 *a[3], const mpq2 *b[3]) +{ + /* TODO: Could experiment with trying bounding box tests before these. + * TODO: Find a less expensive way than 18 orient tests to do this. */ + + /* `orients[0][ai][bi]` is orient of point `a[ai]` compared to segment starting at `b[bi]`. + * `orients[1][bi][ai]` is orient of point `b[bi]` compared to segment starting at `a[ai]`. */ + int orients[2][3][3]; + for (int ab = 0; ab < 2; ++ab) { + for (int ai = 0; ai < 3; ++ai) { + for (int bi = 0; bi < 3; ++bi) { + if (ab == 0) { + orients[0][ai][bi] = orient2d(*b[bi], *b[(bi + 1) % 3], *a[ai]); + } + else { + orients[1][bi][ai] = orient2d(*a[ai], *a[(ai + 1) % 3], *b[bi]); + } + } + } + } + return non_trivially_2d_point_in_tri(orients[0], 0) || + non_trivially_2d_point_in_tri(orients[0], 1) || + non_trivially_2d_point_in_tri(orients[0], 2) || + non_trivially_2d_point_in_tri(orients[1], 0) || + non_trivially_2d_point_in_tri(orients[1], 1) || + non_trivially_2d_point_in_tri(orients[1], 2) || non_trivially_2d_hex_overlap(orients) || + non_trivially_2d_shared_edge_overlap(orients, a, b) || same_triangles(a, b); + return true; +} + +/** + * Does triangle t in tm non-trivially non-co-planar intersect any triangle + * in `CoplanarCluster cl`? Assume t is known to be in the same plane as all + * the triangles in cl, and that proj_axis is a good axis to project down + * to solve this problem in 2d. + */ +static bool non_trivially_coplanar_intersects(const IMesh &tm, + int t, + const CoplanarCluster &cl, + int proj_axis) +{ + const Face &tri = *tm.face(t); + mpq2 v0 = project_3d_to_2d(tri[0]->co_exact, proj_axis); + mpq2 v1 = project_3d_to_2d(tri[1]->co_exact, proj_axis); + mpq2 v2 = project_3d_to_2d(tri[2]->co_exact, proj_axis); + if (orient2d(v0, v1, v2) != 1) { + mpq2 tmp = v1; + v1 = v2; + v2 = tmp; + } + for (const int cl_t : cl) { + const Face &cl_tri = *tm.face(cl_t); + mpq2 ctv0 = project_3d_to_2d(cl_tri[0]->co_exact, proj_axis); + mpq2 ctv1 = project_3d_to_2d(cl_tri[1]->co_exact, proj_axis); + mpq2 ctv2 = project_3d_to_2d(cl_tri[2]->co_exact, proj_axis); + if (orient2d(ctv0, ctv1, ctv2) != 1) { + mpq2 tmp = ctv1; + ctv1 = ctv2; + ctv2 = tmp; + } + const mpq2 *v[] = {&v0, &v1, &v2}; + const mpq2 *ctv[] = {&ctv0, &ctv1, &ctv2}; + if (non_trivially_2d_intersect(v, ctv)) { + return true; + } + } + return false; +} + +/* Keeping this code for a while, but for now, almost all + * trivial intersects are found before calling intersect_tri_tri now. + */ +# if 0 +/** + * Do tri1 and tri2 intersect at all, and if so, is the intersection + * something other than a common vertex or a common edge? + * The \a itt value is the result of calling intersect_tri_tri on tri1, tri2. + */ +static bool non_trivial_intersect(const ITT_value &itt, const Face * tri1, const Face * tri2) +{ + if (itt.kind == INONE) { + return false; + } + const Face * tris[2] = {tri1, tri2}; + if (itt.kind == IPOINT) { + bool has_p_as_vert[2] {false, false}; + for (int i = 0; i < 2; ++i) { + for (const Vert * v : *tris[i]) { + if (itt.p1 == v->co_exact) { + has_p_as_vert[i] = true; + break; + } + } + } + return !(has_p_as_vert[0] && has_p_as_vert[1]); + } + if (itt.kind == ISEGMENT) { + bool has_seg_as_edge[2] = {false, false}; + for (int i = 0; i < 2; ++i) { + const Face &t = *tris[i]; + for (int pos : t.index_range()) { + int nextpos = t.next_pos(pos); + if ((itt.p1 == t[pos]->co_exact && itt.p2 == t[nextpos]->co_exact) || + (itt.p2 == t[pos]->co_exact && itt.p1 == t[nextpos]->co_exact)) { + has_seg_as_edge[i] = true; + break; + } + } + } + return !(has_seg_as_edge[0] && has_seg_as_edge[1]); + } + BLI_assert(itt.kind == ICOPLANAR); + /* TODO: refactor this common code with code above. */ + int proj_axis = mpq3::dominant_axis(tri1->plane.norm_exact); + mpq2 tri_2d[2][3]; + for (int i = 0; i < 2; ++i) { + mpq2 v0 = project_3d_to_2d((*tris[i])[0]->co_exact, proj_axis); + mpq2 v1 = project_3d_to_2d((*tris[i])[1]->co_exact, proj_axis); + mpq2 v2 = project_3d_to_2d((*tris[i])[2]->co_exact, proj_axis); + if (mpq2::orient2d(v0, v1, v2) != 1) { + mpq2 tmp = v1; + v1 = v2; + v2 = tmp; + } + tri_2d[i][0] = v0; + tri_2d[i][1] = v1; + tri_2d[i][2] = v2; + } + const mpq2 *va[] = {&tri_2d[0][0], &tri_2d[0][1], &tri_2d[0][2]}; + const mpq2 *vb[] = {&tri_2d[1][0], &tri_2d[1][1], &tri_2d[1][2]}; + return non_trivially_2d_intersect(va, vb); +} +# endif + +/** + * The sup and index functions are defined in the paper: + * EXACT GEOMETRIC COMPUTATION USING CASCADING, by + * Burnikel, Funke, and Seel. They are used to find absolute + * bounds on the error due to doing a calculation in double + * instead of exactly. For calculations involving only +, -, and *, + * the supremum is the same function except using absolute values + * on inputs and using + instead of -. + * The index function follows these rules: + * index(x op y) = 1 + max(index(x), index(y)) for op + or - + * index(x * y) = 1 + index(x) + index(y) + * index(x) = 0 if input x can be represented exactly as a double + * index(x) = 1 otherwise. + * + * With these rules in place, we know an absolute error bound: + * + * |E_exact - E| <= supremum(E) * index(E) * DBL_EPSILON + * + * where E_exact is what would have been the exact value of the + * expression and E is the one calculated with doubles. + * + * So the sign of E is the same as the sign of E_exact if + * |E| > supremum(E) * index(E) * DBL_EPSILON + * + * Note: a possible speedup would be to have a simple function + * that calculates the error bound if one knows that all values + * are less than some global maximum - most of the function would + * be calculated ahead of time. The global max could be passed + * from above. + */ +static double supremum_dot_cross(const double3 &a, const double3 &b) +{ + double3 abs_a = double3::abs(a); + double3 abs_b = double3::abs(b); + double3 c; + /* This is dot(cross(a, b), cross(a,b)) but using absolute values for a and b + * and always using + when operation is + or -. */ + c[0] = abs_a[1] * abs_b[2] + abs_a[2] * abs_b[1]; + c[1] = abs_a[2] * abs_b[0] + abs_a[0] * abs_b[2]; + c[2] = abs_a[0] * abs_b[1] + abs_a[1] * abs_b[0]; + return double3::dot(c, c); +} + +/** + * Used with supremum to get error bound. See Burnikel et al paper. + * index_plane_coord is the index of a plane coordinate calculated + * for a triangle using the usual formula, assuming the input + * coordinates have index 1. + * index_cross is the index of each coordinate of the cross product. + * It is actually 2 + 2 * (max index of input coords). + * index_dot_cross is the index of the dot product of two cross products. + * It is actually 7 + 4 * (max index of input coords) + */ +constexpr int index_dot_cross = 11; + +static double supremum_dot(const double3 &a, const double3 &b) +{ + double3 abs_a = double3::abs(a); + double3 abs_b = double3::abs(b); + return double3::dot(abs_a, abs_b); +} + +/* Actually index_dot = 3 + 2 * (max index of input coordinates). */ +/* The index of dot when inputs are plane_coords with index 1 is much higher. + * Plane coords have index 6. + */ +constexpr int index_dot_plane_coords = 15; + +static double supremum_orient3d(const double3 &a, + const double3 &b, + const double3 &c, + const double3 &d) +{ + double3 abs_a = double3::abs(a); + double3 abs_b = double3::abs(b); + double3 abs_c = double3::abs(c); + double3 abs_d = double3::abs(d); + double adx = abs_a[0] + abs_d[0]; + double bdx = abs_b[0] + abs_d[0]; + double cdx = abs_c[0] + abs_d[0]; + double ady = abs_a[1] + abs_d[1]; + double bdy = abs_b[1] + abs_d[1]; + double cdy = abs_c[1] + abs_d[1]; + double adz = abs_a[2] + abs_d[2]; + double bdz = abs_b[2] + abs_d[2]; + double cdz = abs_c[2] + abs_d[2]; + + double bdxcdy = bdx * cdy; + double cdxbdy = cdx * bdy; + + double cdxady = cdx * ady; + double adxcdy = adx * cdy; + + double adxbdy = adx * bdy; + double bdxady = bdx * ady; + + double det = adz * (bdxcdy + cdxbdy) + bdz * (cdxady + adxcdy) + cdz * (adxbdy + bdxady); + return det; +} + +/** Actually index_orient3d = 10 + 4 * (max degree of input coordinates) */ +constexpr int index_orient3d = 14; + +/** + * Return the approximate orient3d of the four double3's, with + * the guarantee that if the value is -1 or 1 then the underlying + * mpq3 test would also have returned that value. + * When the return value is 0, we are not sure of the sign. + */ +static int filter_orient3d(const double3 &a, const double3 &b, const double3 &c, const double3 &d) +{ + double o3dfast = orient3d_fast(a, b, c, d); + if (o3dfast == 0.0) { + return 0; + } + double err_bound = supremum_orient3d(a, b, c, d) * index_orient3d * DBL_EPSILON; + if (fabs(o3dfast) > err_bound) { + return o3dfast > 0.0 ? 1 : -1; + } + return 0; +} + +/** + * Return the approximate orient3d of the triangle plane points and v, with + * the guarantee that if the value is -1 or 1 then the underlying + * mpq3 test would also have returned that value. + * When the return value is 0, we are not sure of the sign. + */ +static int filter_tri_plane_vert_orient3d(const Face &tri, const Vert *v) +{ + return filter_orient3d(tri[0]->co, tri[1]->co, tri[2]->co, v->co); +} + +/** + * Are vectors a and b parallel or nearly parallel? + * This routine should only return false if we are certain + * that they are not parallel, taking into account the + * possible numeric errors and input value approximation. + */ +static bool near_parallel_vecs(const double3 &a, const double3 &b) +{ + double3 cr = double3::cross_high_precision(a, b); + double cr_len_sq = cr.length_squared(); + if (cr_len_sq == 0.0) { + return true; + } + double err_bound = supremum_dot_cross(a, b) * index_dot_cross * DBL_EPSILON; + if (cr_len_sq > err_bound) { + return false; + } + return true; +} + +/** + * Return true if we are sure that dot(a,b) > 0, taking into + * account the error bounds due to numeric errors and input value + * approximation. + */ +static bool dot_must_be_positive(const double3 &a, const double3 &b) +{ + double d = double3::dot(a, b); + if (d <= 0.0) { + return false; + } + double err_bound = supremum_dot(a, b) * index_dot_plane_coords * DBL_EPSILON; + if (d > err_bound) { + return true; + } + return false; +} + +/** + * Return the approximate side of point p on a plane with normal plane_no and point plane_p. + * The answer will be 1 if p is definitely above the plane, -1 if it is definitely below. + * If the answer is 0, we are unsure about which side of the plane (or if it is on the plane). + * In exact arithmetic, the answer is just `sgn(dot(p - plane_p, plane_no))`. + * + * The plane_no input is constructed, so has a higher index. + */ +constexpr int index_plane_side = 3 + 2 * index_dot_plane_coords; + +static int filter_plane_side(const double3 &p, + const double3 &plane_p, + const double3 &plane_no, + const double3 &abs_p, + const double3 &abs_plane_p, + const double3 &abs_plane_no) +{ + double d = double3::dot(p - plane_p, plane_no); + if (d == 0.0) { + return 0; + } + double supremum = double3::dot(abs_p + abs_plane_p, abs_plane_no); + double err_bound = supremum * index_plane_side * DBL_EPSILON; + if (d > err_bound) { + return d > 0 ? 1 : -1; + } + return 0; +} + +/** + * A fast, non-exhaustive test for non_trivial intersection. + * If this returns false then we are sure that tri1 and tri2 + * do not intersect. If it returns true, they may or may not + * non-trivially intersect. + * We assume that bounding box overlap tests have already been + * done, so don't repeat those here. This routine is checking + * for the very common cases (when doing mesh self-intersect) + * where triangles share an edge or a vertex, but don't + * otherwise intersect. + */ +static bool may_non_trivially_intersect(Face *t1, Face *t2) +{ + Face &tri1 = *t1; + Face &tri2 = *t2; + BLI_assert(t1->plane_populated() && t2->plane_populated()); + Face::FacePos share1_pos[3]; + Face::FacePos share2_pos[3]; + int n_shared = 0; + for (Face::FacePos p1 = 0; p1 < 3; ++p1) { + const Vert *v1 = tri1[p1]; + for (Face::FacePos p2 = 0; p2 < 3; ++p2) { + const Vert *v2 = tri2[p2]; + if (v1 == v2) { + share1_pos[n_shared] = p1; + share2_pos[n_shared] = p2; + ++n_shared; + } + } + } + if (n_shared == 2) { + /* t1 and t2 share an entire edge. + * If their normals are not parallel, they cannot non-trivially intersect. */ + if (!near_parallel_vecs(tri1.plane->norm, tri2.plane->norm)) { + return false; + } + /* The normals are parallel or nearly parallel. + * If the normals are in the same direction and the edges have opposite + * directions in the two triangles, they cannot non-trivially intersect. */ + bool erev1 = tri1.prev_pos(share1_pos[0]) == share1_pos[1]; + bool erev2 = tri2.prev_pos(share2_pos[0]) == share2_pos[1]; + if (erev1 != erev2) { + if (dot_must_be_positive(tri1.plane->norm, tri2.plane->norm)) { + return false; + } + } + } + else if (n_shared == 1) { + /* t1 and t2 share a vertex, but not an entire edge. + * If the two non-shared verts of t2 are both on the same + * side of tri1's plane, then they cannot non-trivially intersect. + * (There are some other cases that could be caught here but + * they are more expensive to check). */ + Face::FacePos p = share2_pos[0]; + const Vert *v2a = p == 0 ? tri2[1] : tri2[0]; + const Vert *v2b = (p == 0 || p == 1) ? tri2[2] : tri2[1]; + int o1 = filter_tri_plane_vert_orient3d(tri1, v2a); + int o2 = filter_tri_plane_vert_orient3d(tri1, v2b); + if (o1 == o2 && o1 != 0) { + return false; + } + p = share1_pos[0]; + const Vert *v1a = p == 0 ? tri1[1] : tri1[0]; + const Vert *v1b = (p == 0 || p == 1) ? tri1[2] : tri1[1]; + o1 = filter_tri_plane_vert_orient3d(tri2, v1a); + o2 = filter_tri_plane_vert_orient3d(tri2, v1b); + if (o1 == o2 && o1 != 0) { + return false; + } + } + /* We weren't able to prove that any intersection is trivial. */ + return true; +} + +/* + * interesect_tri_tri and helper functions. + * This code uses the algorithm of Guigue and Devillers, as described + * in "Faster Triangle-Triangle Intersection Tests". + * Adapted from github code by Eric Haines: + * github.com/erich666/jgt-code/tree/master/Volume_08/Number_1/Guigue2003 + */ + +/** + * Return the point on ab where the plane with normal n containing point c intersects it. + * Assumes ab is not perpendicular to n. + * This works because the ratio of the projections of ab and ac onto n is the same as + * the ratio along the line ab of the intersection point to the whole of ab. + */ +static inline mpq3 tti_interp(const mpq3 &a, const mpq3 &b, const mpq3 &c, const mpq3 &n) +{ + mpq3 ab = a - b; + mpq_class den = mpq3::dot(ab, n); + BLI_assert(den != 0); + mpq_class alpha = mpq3::dot(a - c, n) / den; + return a - alpha * ab; +} + +/** + * Return +1, 0, -1 as a + ad is above, on, or below the oriented plane containing a, b, c in CCW + * order. This is the same as -oriented(a, b, c, a + ad), but uses fewer arithmetic operations. + * TODO: change arguments to `const Vert *` and use floating filters. + */ +static inline int tti_above(const mpq3 &a, const mpq3 &b, const mpq3 &c, const mpq3 &ad) +{ + mpq3 n = mpq3::cross(b - a, c - a); + return sgn(mpq3::dot(ad, n)); +} + +/** + * Given that triangles (p1, q1, r1) and (p2, q2, r2) are in canonical order, + * use the classification chart in the Guigue and Devillers paper to find out + * how the intervals [i,j] and [k,l] overlap, where [i,j] is where p1r1 and p1q1 + * intersect the plane-plane intersection line, L, and [k,l] is where p2q2 and p2r2 + * intersect L. By the canonicalization, those segments intersect L exactly once. + * Canonicalization has made it so that for p1, q1, r1, either: + * (a)) p1 is off the second triangle's plane and both q1 and r1 are either + * on the plane or on the other side of it from p1; or + * (b) p1 is on the plane both q1 and r1 are on the same side + * of the plane and at least one of q1 and r1 are off the plane. + * Similarly for p2, q2, r2 with respect to the first triangle's plane. + */ +static ITT_value itt_canon2(const mpq3 &p1, + const mpq3 &q1, + const mpq3 &r1, + const mpq3 &p2, + const mpq3 &q2, + const mpq3 &r2, + const mpq3 &n1, + const mpq3 &n2) +{ + constexpr int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "\ntri_tri_intersect_canon:\n"; + std::cout << "p1=" << p1 << " q1=" << q1 << " r1=" << r1 << "\n"; + std::cout << "p2=" << p2 << " q2=" << q2 << " r2=" << r2 << "\n"; + std::cout << "n1=" << n1 << " n2=" << n2 << "\n"; + std::cout << "approximate values:\n"; + std::cout << "p1=(" << p1[0].get_d() << "," << p1[1].get_d() << "," << p1[2].get_d() << ")\n"; + std::cout << "q1=(" << q1[0].get_d() << "," << q1[1].get_d() << "," << q1[2].get_d() << ")\n"; + std::cout << "r1=(" << r1[0].get_d() << "," << r1[1].get_d() << "," << r1[2].get_d() << ")\n"; + std::cout << "p2=(" << p2[0].get_d() << "," << p2[1].get_d() << "," << p2[2].get_d() << ")\n"; + std::cout << "q2=(" << q2[0].get_d() << "," << q2[1].get_d() << "," << q2[2].get_d() << ")\n"; + std::cout << "r2=(" << r2[0].get_d() << "," << r2[1].get_d() << "," << r2[2].get_d() << ")\n"; + std::cout << "n1=(" << n1[0].get_d() << "," << n1[1].get_d() << "," << n1[2].get_d() << ")\n"; + std::cout << "n2=(" << n2[0].get_d() << "," << n2[1].get_d() << "," << n2[2].get_d() << ")\n"; + } + mpq3 p1p2 = p2 - p1; + mpq3 intersect_1; + mpq3 intersect_2; + bool no_overlap = false; + /* Top test in classification tree. */ + if (tti_above(p1, q1, r2, p1p2) > 0) { + /* Middle right test in classification tree. */ + if (tti_above(p1, r1, r2, p1p2) <= 0) { + /* Bottom right test in classification tree. */ + if (tti_above(p1, r1, q2, p1p2) > 0) { + /* Overlap is [k [i l] j]. */ + if (dbg_level > 0) { + std::cout << "overlap [k [i l] j]\n"; + } + /* i is intersect with p1r1. l is intersect with p2r2. */ + intersect_1 = tti_interp(p1, r1, p2, n2); + intersect_2 = tti_interp(p2, r2, p1, n1); + } + else { + /* Overlap is [i [k l] j]. */ + if (dbg_level > 0) { + std::cout << "overlap [i [k l] j]\n"; + } + /* k is intersect with p2q2. l is intersect is p2r2. */ + intersect_1 = tti_interp(p2, q2, p1, n1); + intersect_2 = tti_interp(p2, r2, p1, n1); + } + } + else { + /* No overlap: [k l] [i j]. */ + if (dbg_level > 0) { + std::cout << "no overlap: [k l] [i j]\n"; + } + no_overlap = true; + } + } + else { + /* Middle left test in classification tree. */ + if (tti_above(p1, q1, q2, p1p2) < 0) { + /* No overlap: [i j] [k l]. */ + if (dbg_level > 0) { + std::cout << "no overlap: [i j] [k l]\n"; + } + no_overlap = true; + } + else { + /* Bottom left test in classification tree. */ + if (tti_above(p1, r1, q2, p1p2) >= 0) { + /* Overlap is [k [i j] l]. */ + if (dbg_level > 0) { + std::cout << "overlap [k [i j] l]\n"; + } + /* i is intersect with p1r1. j is intersect with p1q1. */ + intersect_1 = tti_interp(p1, r1, p2, n2); + intersect_2 = tti_interp(p1, q1, p2, n2); + } + else { + /* Overlap is [i [k j] l]. */ + if (dbg_level > 0) { + std::cout << "overlap [i [k j] l]\n"; + } + /* k is intersect with p2q2. j is intersect with p1q1. */ + intersect_1 = tti_interp(p2, q2, p1, n1); + intersect_2 = tti_interp(p1, q1, p2, n2); + } + } + } + if (no_overlap) { + return ITT_value(INONE); + } + if (intersect_1 == intersect_2) { + if (dbg_level > 0) { + std::cout << "single intersect: " << intersect_1 << "\n"; + } + return ITT_value(IPOINT, intersect_1); + } + if (dbg_level > 0) { + std::cout << "intersect segment: " << intersect_1 << ", " << intersect_2 << "\n"; + } + return ITT_value(ISEGMENT, intersect_1, intersect_2); +} + +/* Helper function for intersect_tri_tri. Args have been canonicalized for triangle 1. */ + +static ITT_value itt_canon1(const mpq3 &p1, + const mpq3 &q1, + const mpq3 &r1, + const mpq3 &p2, + const mpq3 &q2, + const mpq3 &r2, + const mpq3 &n1, + const mpq3 &n2, + int sp2, + int sq2, + int sr2) +{ + constexpr int dbg_level = 0; + if (sp2 > 0) { + if (sq2 > 0) { + return itt_canon2(p1, r1, q1, r2, p2, q2, n1, n2); + } + if (sr2 > 0) { + return itt_canon2(p1, r1, q1, q2, r2, p2, n1, n2); + } + return itt_canon2(p1, q1, r1, p2, q2, r2, n1, n2); + } + if (sp2 < 0) { + if (sq2 < 0) { + return itt_canon2(p1, q1, r1, r2, p2, q2, n1, n2); + } + if (sr2 < 0) { + return itt_canon2(p1, q1, r1, q2, r2, p2, n1, n2); + } + return itt_canon2(p1, r1, q1, p2, q2, r2, n1, n2); + } + if (sq2 < 0) { + if (sr2 >= 0) { + return itt_canon2(p1, r1, q1, q2, r2, p2, n1, n2); + } + return itt_canon2(p1, q1, r1, p2, q2, r2, n1, n2); + } + if (sq2 > 0) { + if (sr2 > 0) { + return itt_canon2(p1, r1, q1, p2, q2, r2, n1, n2); + } + return itt_canon2(p1, q1, r1, q2, r2, p2, n1, n2); + } + if (sr2 > 0) { + return itt_canon2(p1, q1, r1, r2, p2, q2, n1, n2); + } + if (sr2 < 0) { + return itt_canon2(p1, r1, q1, r2, p2, q2, n1, n2); + } + if (dbg_level > 0) { + std::cout << "triangles are co-planar\n"; + } + return ITT_value(ICOPLANAR); +} + +static ITT_value intersect_tri_tri(const IMesh &tm, int t1, int t2) +{ + constexpr int dbg_level = 0; +# ifdef PERFDEBUG + incperfcount(1); /* Intersect_tri_tri calls. */ +# endif + const Face &tri1 = *tm.face(t1); + const Face &tri2 = *tm.face(t2); + BLI_assert(tri1.plane_populated() && tri2.plane_populated()); + const Vert *vp1 = tri1[0]; + const Vert *vq1 = tri1[1]; + const Vert *vr1 = tri1[2]; + const Vert *vp2 = tri2[0]; + const Vert *vq2 = tri2[1]; + const Vert *vr2 = tri2[2]; + if (dbg_level > 0) { + std::cout << "\nINTERSECT_TRI_TRI t1=" << t1 << ", t2=" << t2 << "\n"; + std::cout << " p1 = " << vp1 << "\n"; + std::cout << " q1 = " << vq1 << "\n"; + std::cout << " r1 = " << vr1 << "\n"; + std::cout << " p2 = " << vp2 << "\n"; + std::cout << " q2 = " << vq2 << "\n"; + std::cout << " r2 = " << vr2 << "\n"; + } + + /* Get signs of t1's vertices' distances to plane of t2 and vice versa. */ + + /* Try first getting signs with double arithmetic, with error bounds. + * If the signs calculated in this section are not 0, they are the same + * as what they would be using exact arithmetic. */ + const double3 &d_p1 = vp1->co; + const double3 &d_q1 = vq1->co; + const double3 &d_r1 = vr1->co; + const double3 &d_p2 = vp2->co; + const double3 &d_q2 = vq2->co; + const double3 &d_r2 = vr2->co; + const double3 &d_n2 = tri2.plane->norm; + + const double3 &abs_d_p1 = double3::abs(d_p1); + const double3 &abs_d_q1 = double3::abs(d_q1); + const double3 &abs_d_r1 = double3::abs(d_r1); + const double3 &abs_d_r2 = double3::abs(d_r2); + const double3 &abs_d_n2 = double3::abs(d_n2); + + int sp1 = filter_plane_side(d_p1, d_r2, d_n2, abs_d_p1, abs_d_r2, abs_d_n2); + int sq1 = filter_plane_side(d_q1, d_r2, d_n2, abs_d_q1, abs_d_r2, abs_d_n2); + int sr1 = filter_plane_side(d_r1, d_r2, d_n2, abs_d_r1, abs_d_r2, abs_d_n2); + if ((sp1 > 0 && sq1 > 0 && sr1 > 0) || (sp1 < 0 && sq1 < 0 && sr1 < 0)) { +# ifdef PERFDEBUG + incperfcount(2); /* Tri tri intersects decided by filter plane tests. */ +# endif + if (dbg_level > 0) { + std::cout << "no intersection, all t1's verts above or below t2\n"; + } + return ITT_value(INONE); + } + + const double3 &d_n1 = tri1.plane->norm; + const double3 &abs_d_p2 = double3::abs(d_p2); + const double3 &abs_d_q2 = double3::abs(d_q2); + const double3 &abs_d_n1 = double3::abs(d_n1); + + int sp2 = filter_plane_side(d_p2, d_r1, d_n1, abs_d_p2, abs_d_r1, abs_d_n1); + int sq2 = filter_plane_side(d_q2, d_r1, d_n1, abs_d_q2, abs_d_r1, abs_d_n1); + int sr2 = filter_plane_side(d_r2, d_r1, d_n1, abs_d_r2, abs_d_r1, abs_d_n1); + if ((sp2 > 0 && sq2 > 0 && sr2 > 0) || (sp2 < 0 && sq2 < 0 && sr2 < 0)) { +# ifdef PERFDEBUG + incperfcount(2); /* Tri tri intersects decided by filter plane tests. */ +# endif + if (dbg_level > 0) { + std::cout << "no intersection, all t2's verts above or below t1\n"; + } + return ITT_value(INONE); + } + + const mpq3 &p1 = vp1->co_exact; + const mpq3 &q1 = vq1->co_exact; + const mpq3 &r1 = vr1->co_exact; + const mpq3 &p2 = vp2->co_exact; + const mpq3 &q2 = vq2->co_exact; + const mpq3 &r2 = vr2->co_exact; + + const mpq3 &n2 = tri2.plane->norm_exact; + if (sp1 == 0) { + sp1 = sgn(mpq3::dot(p1 - r2, n2)); + } + if (sq1 == 0) { + sq1 = sgn(mpq3::dot(q1 - r2, n2)); + } + if (sr1 == 0) { + sr1 = sgn(mpq3::dot(r1 - r2, n2)); + } + + if (dbg_level > 1) { + std::cout << " sp1=" << sp1 << " sq1=" << sq1 << " sr1=" << sr1 << "\n"; + } + + if ((sp1 * sq1 > 0) && (sp1 * sr1 > 0)) { + if (dbg_level > 0) { + std::cout << "no intersection, all t1's verts above or below t2 (exact)\n"; + } +# ifdef PERFDEBUG + incperfcount(3); /* Tri tri intersects decided by exact plane tests. */ +# endif + return ITT_value(INONE); + } + + /* Repeat for signs of t2's vertices with respect to plane of t1. */ + const mpq3 &n1 = tri1.plane->norm_exact; + if (sp2 == 0) { + sp2 = sgn(mpq3::dot(p2 - r1, n1)); + } + if (sq2 == 0) { + sq2 = sgn(mpq3::dot(q2 - r1, n1)); + } + if (sr2 == 0) { + sr2 = sgn(mpq3::dot(r2 - r1, n1)); + } + + if (dbg_level > 1) { + std::cout << " sp2=" << sp2 << " sq2=" << sq2 << " sr2=" << sr2 << "\n"; + } + + if ((sp2 * sq2 > 0) && (sp2 * sr2 > 0)) { + if (dbg_level > 0) { + std::cout << "no intersection, all t2's verts above or below t1 (exact)\n"; + } +# ifdef PERFDEBUG + incperfcount(3); /* Tri tri intersects decided by exact plane tests. */ +# endif + return ITT_value(INONE); + } + + /* Do rest of the work with vertices in a canonical order, where p1 is on + * positive side of plane and q1, r1 are not, or p1 is on the plane and + * q1 and r1 are off the plane on the same side. */ + ITT_value ans; + if (sp1 > 0) { + if (sq1 > 0) { + ans = itt_canon1(r1, p1, q1, p2, r2, q2, n1, n2, sp2, sr2, sq2); + } + else if (sr1 > 0) { + ans = itt_canon1(q1, r1, p1, p2, r2, q2, n1, n2, sp2, sr2, sq2); + } + else { + ans = itt_canon1(p1, q1, r1, p2, q2, r2, n1, n2, sp2, sq2, sr2); + } + } + else if (sp1 < 0) { + if (sq1 < 0) { + ans = itt_canon1(r1, p1, q1, p2, q2, r2, n1, n2, sp2, sq2, sr2); + } + else if (sr1 < 0) { + ans = itt_canon1(q1, r1, p1, p2, q2, r2, n1, n2, sp2, sq2, sr2); + } + else { + ans = itt_canon1(p1, q1, r1, p2, r2, q2, n1, n2, sp2, sr2, sq2); + } + } + else { + if (sq1 < 0) { + if (sr1 >= 0) { + ans = itt_canon1(q1, r1, p1, p2, r2, q2, n1, n2, sp2, sr2, sq2); + } + else { + ans = itt_canon1(p1, q1, r1, p2, q2, r2, n1, n2, sp2, sq2, sr2); + } + } + else if (sq1 > 0) { + if (sr1 > 0) { + ans = itt_canon1(p1, q1, r1, p2, r2, q2, n1, n2, sp2, sr2, sq2); + } + else { + ans = itt_canon1(q1, r1, p1, p2, q2, r2, n1, n2, sp2, sq2, sr2); + } + } + else { + if (sr1 > 0) { + ans = itt_canon1(r1, p1, q1, p2, q2, r2, n1, n2, sp2, sq2, sr2); + } + else if (sr1 < 0) { + ans = itt_canon1(r1, p1, q1, p2, r2, q2, n1, n2, sp2, sr2, sq2); + } + else { + if (dbg_level > 0) { + std::cout << "triangles are co-planar\n"; + } + ans = ITT_value(ICOPLANAR); + } + } + } + if (ans.kind == ICOPLANAR) { + ans.t_source = t2; + } + +# ifdef PERFDEBUG + if (ans.kind != INONE) { + incperfcount(4); + } +# endif + return ans; +} + +struct CDT_data { + const Plane *t_plane; + Vector<mpq2> vert; + Vector<std::pair<int, int>> edge; + Vector<Vector<int>> face; + /** Parallels face, gives id from input #IMesh of input face. */ + Vector<int> input_face; + /** Parallels face, says if input face orientation is opposite. */ + Vector<bool> is_reversed; + /** Result of running CDT on input with (vert, edge, face). */ + CDT_result<mpq_class> cdt_out; + int proj_axis; +}; + +/** + * We could de-duplicate verts here, but CDT routine will do that anyway. + */ +static int prepare_need_vert(CDT_data &cd, const mpq3 &p3d) +{ + mpq2 p2d = project_3d_to_2d(p3d, cd.proj_axis); + int v = cd.vert.append_and_get_index(p2d); + return v; +} + +/** + * To un-project a 2d vert that was projected along cd.proj_axis, we copy the coordinates + * from the two axes not involved in the projection, and use the plane equation of the + * originating 3d plane, cd.t_plane, to derive the coordinate of the projected axis. + * The plane equation says a point p is on the plane if dot(p, plane.n()) + plane.d() == 0. + * Assume that the projection axis is such that plane.n()[proj_axis] != 0. + */ +static mpq3 unproject_cdt_vert(const CDT_data &cd, const mpq2 &p2d) +{ + mpq3 p3d; + BLI_assert(cd.t_plane->exact_populated()); + BLI_assert(cd.t_plane->norm_exact[cd.proj_axis] != 0); + const mpq3 &n = cd.t_plane->norm_exact; + const mpq_class &d = cd.t_plane->d_exact; + switch (cd.proj_axis) { + case (0): { + mpq_class num = n[1] * p2d[0] + n[2] * p2d[1] + d; + num = -num; + p3d[0] = num / n[0]; + p3d[1] = p2d[0]; + p3d[2] = p2d[1]; + break; + } + case (1): { + p3d[0] = p2d[0]; + mpq_class num = n[0] * p2d[0] + n[2] * p2d[1] + d; + num = -num; + p3d[1] = num / n[1]; + p3d[2] = p2d[1]; + break; + } + case (2): { + p3d[0] = p2d[0]; + p3d[1] = p2d[1]; + mpq_class num = n[0] * p2d[0] + n[1] * p2d[1] + d; + num = -num; + p3d[2] = num / n[2]; + break; + } + default: + BLI_assert(false); + } + return p3d; +} + +static void prepare_need_edge(CDT_data &cd, const mpq3 &p1, const mpq3 &p2) +{ + int v1 = prepare_need_vert(cd, p1); + int v2 = prepare_need_vert(cd, p2); + cd.edge.append(std::pair<int, int>(v1, v2)); +} + +static void prepare_need_tri(CDT_data &cd, const IMesh &tm, int t) +{ + const Face &tri = *tm.face(t); + int v0 = prepare_need_vert(cd, tri[0]->co_exact); + int v1 = prepare_need_vert(cd, tri[1]->co_exact); + int v2 = prepare_need_vert(cd, tri[2]->co_exact); + bool rev; + /* How to get CCW orientation of projected triangle? Note that when look down y axis + * as opposed to x or z, the orientation of the other two axes is not right-and-up. */ + BLI_assert(cd.t_plane->exact_populated()); + if (tri.plane->norm_exact[cd.proj_axis] >= 0) { + rev = cd.proj_axis == 1; + } + else { + rev = cd.proj_axis != 1; + } + int cd_t = cd.face.append_and_get_index(Vector<int>()); + cd.face[cd_t].append(v0); + if (rev) { + cd.face[cd_t].append(v2); + cd.face[cd_t].append(v1); + } + else { + cd.face[cd_t].append(v1); + cd.face[cd_t].append(v2); + } + cd.input_face.append(t); + cd.is_reversed.append(rev); +} + +static CDT_data prepare_cdt_input(const IMesh &tm, int t, const Vector<ITT_value> itts) +{ + CDT_data ans; + BLI_assert(tm.face(t)->plane_populated()); + ans.t_plane = tm.face(t)->plane; + BLI_assert(ans.t_plane->exact_populated()); + ans.proj_axis = mpq3::dominant_axis(ans.t_plane->norm_exact); + prepare_need_tri(ans, tm, t); + for (const ITT_value &itt : itts) { + switch (itt.kind) { + case INONE: + break; + case IPOINT: { + prepare_need_vert(ans, itt.p1); + break; + } + case ISEGMENT: { + prepare_need_edge(ans, itt.p1, itt.p2); + break; + } + case ICOPLANAR: { + prepare_need_tri(ans, tm, itt.t_source); + break; + } + } + } + return ans; +} + +static CDT_data prepare_cdt_input_for_cluster(const IMesh &tm, + const CoplanarClusterInfo &clinfo, + int c, + const Vector<ITT_value> itts) +{ + CDT_data ans; + BLI_assert(c < clinfo.tot_cluster()); + const CoplanarCluster &cl = clinfo.cluster(c); + BLI_assert(cl.tot_tri() > 0); + int t0 = cl.tri(0); + BLI_assert(tm.face(t0)->plane_populated()); + ans.t_plane = tm.face(t0)->plane; + BLI_assert(ans.t_plane->exact_populated()); + ans.proj_axis = mpq3::dominant_axis(ans.t_plane->norm_exact); + for (const int t : cl) { + prepare_need_tri(ans, tm, t); + } + for (const ITT_value &itt : itts) { + switch (itt.kind) { + case IPOINT: { + prepare_need_vert(ans, itt.p1); + break; + } + case ISEGMENT: { + prepare_need_edge(ans, itt.p1, itt.p2); + break; + } + default: + break; + } + } + return ans; +} + +/** + * Fills in cd.cdt_out with result of doing the cdt calculation on (vert, edge, face). + */ +static void do_cdt(CDT_data &cd) +{ + constexpr int dbg_level = 0; + CDT_input<mpq_class> cdt_in; + cdt_in.vert = Span<mpq2>(cd.vert); + cdt_in.edge = Span<std::pair<int, int>>(cd.edge); + cdt_in.face = Span<Vector<int>>(cd.face); + if (dbg_level > 0) { + std::cout << "CDT input\nVerts:\n"; + for (int i : cdt_in.vert.index_range()) { + std::cout << "v" << i << ": " << cdt_in.vert[i] << "=(" << cdt_in.vert[i][0].get_d() << "," + << cdt_in.vert[i][1].get_d() << ")\n"; + } + std::cout << "Edges:\n"; + for (int i : cdt_in.edge.index_range()) { + std::cout << "e" << i << ": (" << cdt_in.edge[i].first << ", " << cdt_in.edge[i].second + << ")\n"; + } + std::cout << "Tris\n"; + for (int f : cdt_in.face.index_range()) { + std::cout << "f" << f << ": "; + for (int j : cdt_in.face[f].index_range()) { + std::cout << cdt_in.face[f][j] << " "; + } + std::cout << "\n"; + } + } + cdt_in.epsilon = 0; /* TODO: needs attention for non-exact T. */ + cd.cdt_out = blender::meshintersect::delaunay_2d_calc(cdt_in, CDT_INSIDE); + if (dbg_level > 0) { + std::cout << "\nCDT result\nVerts:\n"; + for (int i : cd.cdt_out.vert.index_range()) { + std::cout << "v" << i << ": " << cd.cdt_out.vert[i] << "=(" << cd.cdt_out.vert[i][0].get_d() + << "," << cd.cdt_out.vert[i][1].get_d() << "\n"; + } + std::cout << "Tris\n"; + for (int f : cd.cdt_out.face.index_range()) { + std::cout << "f" << f << ": "; + for (int j : cd.cdt_out.face[f].index_range()) { + std::cout << cd.cdt_out.face[f][j] << " "; + } + std::cout << "orig: "; + for (int j : cd.cdt_out.face_orig[f].index_range()) { + std::cout << cd.cdt_out.face_orig[f][j] << " "; + } + std::cout << "\n"; + } + std::cout << "Edges\n"; + for (int e : cd.cdt_out.edge.index_range()) { + std::cout << "e" << e << ": (" << cd.cdt_out.edge[e].first << ", " + << cd.cdt_out.edge[e].second << ") "; + std::cout << "orig: "; + for (int j : cd.cdt_out.edge_orig[e].index_range()) { + std::cout << cd.cdt_out.edge_orig[e][j] << " "; + } + std::cout << "\n"; + } + } +} + +static int get_cdt_edge_orig( + int i0, int i1, const CDT_data &cd, const IMesh &in_tm, bool *r_is_intersect) +{ + int foff = cd.cdt_out.face_edge_offset; + *r_is_intersect = false; + for (int e : cd.cdt_out.edge.index_range()) { + std::pair<int, int> edge = cd.cdt_out.edge[e]; + if ((edge.first == i0 && edge.second == i1) || (edge.first == i1 && edge.second == i0)) { + /* Pick an arbitrary orig, but not one equal to NO_INDEX, if we can help it. */ + /* TODO: if edge has origs from more than on part of the nary input, + * then want to set *r_is_intersect to true. */ + for (int orig_index : cd.cdt_out.edge_orig[e]) { + /* orig_index encodes the triangle and pos within the triangle of the input edge. */ + if (orig_index >= foff) { + int in_face_index = (orig_index / foff) - 1; + int pos = orig_index % foff; + /* We need to retrieve the edge orig field from the Face used to populate the + * in_face_index'th face of the CDT, at the pos'th position of the face. */ + int in_tm_face_index = cd.input_face[in_face_index]; + BLI_assert(in_tm_face_index < in_tm.face_size()); + const Face *facep = in_tm.face(in_tm_face_index); + BLI_assert(pos < facep->size()); + bool is_rev = cd.is_reversed[in_face_index]; + int eorig = is_rev ? facep->edge_orig[2 - pos] : facep->edge_orig[pos]; + if (eorig != NO_INDEX) { + return eorig; + } + } + else { + /* This edge came from an edge input to the CDT problem, + * so it is an intersect edge. */ + *r_is_intersect = true; + /* TODO: maybe there is an orig index: + * This happens if an input edge was formed by an input face having + * an edge that is co-planar with the cluster, while the face as a whole is not. */ + return NO_INDEX; + } + } + return NO_INDEX; + } + } + return NO_INDEX; +} + +/** + * Using the result of CDT in cd.cdt_out, extract an #IMesh representing the subdivision + * of input triangle t, which should be an element of cd.input_face. + */ +static IMesh extract_subdivided_tri(const CDT_data &cd, + const IMesh &in_tm, + int t, + IMeshArena *arena) +{ + const CDT_result<mpq_class> &cdt_out = cd.cdt_out; + int t_in_cdt = -1; + for (int i = 0; i < cd.input_face.size(); ++i) { + if (cd.input_face[i] == t) { + t_in_cdt = i; + } + } + if (t_in_cdt == -1) { + std::cout << "Could not find " << t << " in cdt input tris\n"; + BLI_assert(false); + return IMesh(); + } + int t_orig = in_tm.face(t)->orig; + constexpr int inline_buf_size = 20; + Vector<Face *, inline_buf_size> faces; + for (int f : cdt_out.face.index_range()) { + if (cdt_out.face_orig[f].contains(t_in_cdt)) { + BLI_assert(cdt_out.face[f].size() == 3); + int i0 = cdt_out.face[f][0]; + int i1 = cdt_out.face[f][1]; + int i2 = cdt_out.face[f][2]; + mpq3 v0co = unproject_cdt_vert(cd, cdt_out.vert[i0]); + mpq3 v1co = unproject_cdt_vert(cd, cdt_out.vert[i1]); + mpq3 v2co = unproject_cdt_vert(cd, cdt_out.vert[i2]); + /* No need to provide an original index: if coord matches + * an original one, then it will already be in the arena + * with the correct orig field. */ + const Vert *v0 = arena->add_or_find_vert(v0co, NO_INDEX); + const Vert *v1 = arena->add_or_find_vert(v1co, NO_INDEX); + const Vert *v2 = arena->add_or_find_vert(v2co, NO_INDEX); + Face *facep; + bool is_isect0; + bool is_isect1; + bool is_isect2; + if (cd.is_reversed[t_in_cdt]) { + int oe0 = get_cdt_edge_orig(i0, i2, cd, in_tm, &is_isect0); + int oe1 = get_cdt_edge_orig(i2, i1, cd, in_tm, &is_isect1); + int oe2 = get_cdt_edge_orig(i1, i0, cd, in_tm, &is_isect2); + facep = arena->add_face( + {v0, v2, v1}, t_orig, {oe0, oe1, oe2}, {is_isect0, is_isect1, is_isect2}); + } + else { + int oe0 = get_cdt_edge_orig(i0, i1, cd, in_tm, &is_isect0); + int oe1 = get_cdt_edge_orig(i1, i2, cd, in_tm, &is_isect1); + int oe2 = get_cdt_edge_orig(i2, i0, cd, in_tm, &is_isect2); + facep = arena->add_face( + {v0, v1, v2}, t_orig, {oe0, oe1, oe2}, {is_isect0, is_isect1, is_isect2}); + } + facep->populate_plane(false); + faces.append(facep); + } + } + return IMesh(faces); +} + +static IMesh extract_single_tri(const IMesh &tm, int t) +{ + Face *f = tm.face(t); + return IMesh({f}); +} + +static bool bvhtreeverlap_cmp(const BVHTreeOverlap &a, const BVHTreeOverlap &b) +{ + if (a.indexA < b.indexA) { + return true; + } + if ((a.indexA == b.indexA) & (a.indexB < b.indexB)) { + return true; + } + return false; +} +class TriOverlaps { + BVHTree *tree_{nullptr}; + BVHTree *tree_b_{nullptr}; + BVHTreeOverlap *overlap_{nullptr}; + uint overlap_tot_{0}; + + struct CBData { + const IMesh &tm; + std::function<int(int)> shape_fn; + int nshapes; + bool use_self; + }; + + public: + TriOverlaps(const IMesh &tm, + const Array<BoundingBox> &tri_bb, + int nshapes, + std::function<int(int)> shape_fn, + bool use_self) + { + constexpr int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "TriOverlaps construction\n"; + } + /* Tree type is 8 => octtree; axis = 6 => using XYZ axes only. */ + tree_ = BLI_bvhtree_new(tm.face_size(), FLT_EPSILON, 8, 6); + /* In the common case of a binary boolean and no self intersection in + * each shape, we will use two trees and simple bounding box overlap. */ + bool two_trees_no_self = nshapes == 2 && !use_self; + if (two_trees_no_self) { + tree_b_ = BLI_bvhtree_new(tm.face_size(), FLT_EPSILON, 8, 6); + } + float bbpts[6]; + for (int t : tm.face_index_range()) { + const BoundingBox &bb = tri_bb[t]; + copy_v3_v3(bbpts, bb.min); + copy_v3_v3(bbpts + 3, bb.max); + int shape = shape_fn(tm.face(t)->orig); + if (two_trees_no_self) { + if (shape == 0) { + BLI_bvhtree_insert(tree_, t, bbpts, 2); + } + else if (shape == 1) { + BLI_bvhtree_insert(tree_b_, t, bbpts, 2); + } + } + else { + if (shape != -1) { + BLI_bvhtree_insert(tree_, t, bbpts, 2); + } + } + } + BLI_bvhtree_balance(tree_); + if (two_trees_no_self) { + BLI_bvhtree_balance(tree_b_); + /* Don't expect a lot of trivial intersects in this case. */ + overlap_ = BLI_bvhtree_overlap(tree_, tree_b_, &overlap_tot_, NULL, NULL); + } + else { + CBData cbdata{tm, shape_fn, nshapes, use_self}; + if (nshapes == 1 && use_self) { + /* Expect a lot of trivial intersects from quads that are triangulated + * and faces that share vertices. + * Filter them out with a callback. */ + overlap_ = BLI_bvhtree_overlap( + tree_, tree_, &overlap_tot_, only_nontrivial_intersects, &cbdata); + } + else { + overlap_ = BLI_bvhtree_overlap( + tree_, tree_, &overlap_tot_, only_different_shapes, &cbdata); + } + } + /* The rest of the code is simpler and easier to parallelize if, in the two-trees case, + * we repeat the overlaps with indexA and indexB reversed. It is important that + * in the repeated part, sorting will then bring things with indexB together. */ + if (two_trees_no_self) { + overlap_ = static_cast<BVHTreeOverlap *>( + MEM_reallocN(overlap_, 2 * overlap_tot_ * sizeof(overlap_[0]))); + for (uint i = 0; i < overlap_tot_; ++i) { + overlap_[overlap_tot_ + i].indexA = overlap_[i].indexB; + overlap_[overlap_tot_ + i].indexB = overlap_[i].indexA; + } + overlap_tot_ += overlap_tot_; + } + /* Sort the overlaps to bring all the intersects with a given indexA together. */ + std::sort(overlap_, overlap_ + overlap_tot_, bvhtreeverlap_cmp); + if (dbg_level > 0) { + std::cout << overlap_tot_ << " overlaps found:\n"; + for (BVHTreeOverlap ov : overlap()) { + std::cout << "A: " << ov.indexA << ", B: " << ov.indexB << "\n"; + } + } + } + + ~TriOverlaps() + { + if (tree_) { + BLI_bvhtree_free(tree_); + } + if (tree_b_) { + BLI_bvhtree_free(tree_b_); + } + if (overlap_) { + MEM_freeN(overlap_); + } + } + + Span<BVHTreeOverlap> overlap() const + { + return Span<BVHTreeOverlap>(overlap_, overlap_tot_); + } + + private: + static bool only_nontrivial_intersects(void *userdata, + int index_a, + int index_b, + int UNUSED(thread)) + { + CBData *cbdata = static_cast<CBData *>(userdata); + return may_non_trivially_intersect(cbdata->tm.face(index_a), cbdata->tm.face(index_b)); + } + + static bool only_different_shapes(void *userdata, int index_a, int index_b, int UNUSED(thread)) + { + CBData *cbdata = static_cast<CBData *>(userdata); + return cbdata->tm.face(index_a)->orig != cbdata->tm.face(index_b)->orig; + } +}; + +/** + * Data needed for parallelization of #calc_overlap_itts. + */ +struct OverlapIttsData { + Vector<std::pair<int, int>> intersect_pairs; + Map<std::pair<int, int>, ITT_value> &itt_map; + const IMesh &tm; + IMeshArena *arena; + + OverlapIttsData(Map<std::pair<int, int>, ITT_value> &itt_map, const IMesh &tm, IMeshArena *arena) + : itt_map(itt_map), tm(tm), arena(arena) + { + } +}; + +/** + * Return a std::pair containing a and b in canonical order: + * With a <= b. + */ +static std::pair<int, int> canon_int_pair(int a, int b) +{ + if (a > b) { + std::swap(a, b); + } + return std::pair<int, int>(a, b); +} + +static void calc_overlap_itts_range_func(void *__restrict userdata, + const int iter, + const TaskParallelTLS *__restrict UNUSED(tls)) +{ + constexpr int dbg_level = 0; + OverlapIttsData *data = static_cast<OverlapIttsData *>(userdata); + std::pair<int, int> tri_pair = data->intersect_pairs[iter]; + int a = tri_pair.first; + int b = tri_pair.second; + if (dbg_level > 0) { + std::cout << "calc_overlap_itts_range_func a=" << a << ", b=" << b << "\n"; + } + ITT_value itt = intersect_tri_tri(data->tm, a, b); + if (dbg_level > 0) { + std::cout << "result of intersecting " << a << " and " << b << " = " << itt << "\n"; + } + BLI_assert(data->itt_map.contains(tri_pair)); + data->itt_map.add_overwrite(tri_pair, itt); +} + +/** + * Fill in itt_map with the vector of ITT_values that result from intersecting the triangles in ov. + * Use a canonical order for triangles: (a,b) where a < b. + * Don't bother doing this if both a and b are part of the same co-planar cluster, as the + * intersection for those will be handled by CDT, later. + */ +static void calc_overlap_itts(Map<std::pair<int, int>, ITT_value> &itt_map, + const IMesh &tm, + const CoplanarClusterInfo &clinfo, + const TriOverlaps &ov, + IMeshArena *arena) +{ + OverlapIttsData data(itt_map, tm, arena); + /* Put dummy values in `itt_map` initially, + * so map entries will exist when doing the range function. + * This means we won't have to protect the `itt_map.add_overwrite` function with a lock. */ + for (const BVHTreeOverlap &olap : ov.overlap()) { + std::pair<int, int> key = canon_int_pair(olap.indexA, olap.indexB); + if (!itt_map.contains(key)) { + int ca = clinfo.tri_cluster(key.first); + int cb = clinfo.tri_cluster(key.second); + if (ca == NO_INDEX || ca != cb) { + itt_map.add_new(key, ITT_value()); + data.intersect_pairs.append(key); + } + } + } + int tot_intersect_pairs = data.intersect_pairs.size(); + TaskParallelSettings settings; + BLI_parallel_range_settings_defaults(&settings); + settings.min_iter_per_thread = 1000; + settings.use_threading = intersect_use_threading; + BLI_task_parallel_range(0, tot_intersect_pairs, &data, calc_overlap_itts_range_func, &settings); +} + +/** + * Data needed for parallelization of calc_subdivided_tris. + */ +struct OverlapTriRange { + int tri_index; + int overlap_start; + int len; +}; +struct SubdivideTrisData { + Array<IMesh> &r_tri_subdivided; + const IMesh &tm; + const Map<std::pair<int, int>, ITT_value> &itt_map; + Span<BVHTreeOverlap> overlap; + IMeshArena *arena; + + /* This vector gives, for each triangle in tm that has an intersection + * we want to calculate: what the index of that triangle in tm is, + * where it starts in the ov structure as indexA, and how many + * overlap pairs have that same indexA (they will be continuous). */ + Vector<OverlapTriRange> overlap_tri_range; + + SubdivideTrisData(Array<IMesh> &r_tri_subdivided, + const IMesh &tm, + const Map<std::pair<int, int>, ITT_value> &itt_map, + Span<BVHTreeOverlap> overlap, + IMeshArena *arena) + : r_tri_subdivided(r_tri_subdivided), + tm(tm), + itt_map(itt_map), + overlap(overlap), + arena(arena), + overlap_tri_range{} + { + } +}; + +static void calc_subdivided_tri_range_func(void *__restrict userdata, + const int iter, + const TaskParallelTLS *__restrict UNUSED(tls)) +{ + constexpr int dbg_level = 0; + SubdivideTrisData *data = static_cast<SubdivideTrisData *>(userdata); + OverlapTriRange &otr = data->overlap_tri_range[iter]; + int t = otr.tri_index; + if (dbg_level > 0) { + std::cout << "calc_subdivided_tri_range_func\nt=" << t << " start=" << otr.overlap_start + << " len=" << otr.len << "\n"; + } + constexpr int inline_capacity = 100; + Vector<ITT_value, inline_capacity> itts(otr.len); + for (int j = otr.overlap_start; j < otr.overlap_start + otr.len; ++j) { + int t_other = data->overlap[j].indexB; + std::pair<int, int> key = canon_int_pair(t, t_other); + ITT_value itt; + if (data->itt_map.contains(key)) { + itt = data->itt_map.lookup(key); + } + if (itt.kind != INONE) { + itts.append(itt); + } + if (dbg_level > 0) { + std::cout << " tri t" << t_other << "; result = " << itt << "\n"; + } + } + if (itts.size() > 0) { + CDT_data cd_data = prepare_cdt_input(data->tm, t, itts); + do_cdt(cd_data); + data->r_tri_subdivided[t] = extract_subdivided_tri(cd_data, data->tm, t, data->arena); + if (dbg_level > 0) { + std::cout << "subdivide output\n" << data->r_tri_subdivided[t]; + } + } +} + +/** + * For each triangle in tm, fill in the corresponding slot in + * r_tri_subdivided with the result of intersecting it with + * all the other triangles in the mesh, if it intersects any others. + * But don't do this for triangles that are part of a cluster. + * Also, do nothing here if the answer is just the triangle itself. + */ +static void calc_subdivided_tris(Array<IMesh> &r_tri_subdivided, + const IMesh &tm, + const Map<std::pair<int, int>, ITT_value> &itt_map, + const CoplanarClusterInfo &clinfo, + const TriOverlaps &ov, + IMeshArena *arena) +{ + const int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "\nCALC_SUBDIVIDED_TRIS\n\n"; + } + Span<BVHTreeOverlap> overlap = ov.overlap(); + SubdivideTrisData data(r_tri_subdivided, tm, itt_map, overlap, arena); + int overlap_tot = overlap.size(); + data.overlap_tri_range = Vector<OverlapTriRange>(); + data.overlap_tri_range.reserve(overlap_tot); + int overlap_index = 0; + while (overlap_index < overlap_tot) { + int t = overlap[overlap_index].indexA; + int i = overlap_index; + while (i + 1 < overlap_tot && overlap[i + 1].indexA == t) { + ++i; + } + /* Now overlap[overlap_index] to overlap[i] have indexA == t. + * We only record ranges for triangles that are not in clusters, + * because the ones in clusters are handled separately. */ + if (clinfo.tri_cluster(t) == NO_INDEX) { + int len = i - overlap_index + 1; + if (!(len == 1 && overlap[overlap_index].indexB == t)) { + OverlapTriRange range = {t, overlap_index, len}; + data.overlap_tri_range.append(range); +# ifdef PERFDEBUG + bumpperfcount(0, len); /* Non-cluster overlaps. */ +# endif + } + } + overlap_index = i + 1; + } + int overlap_tri_range_tot = data.overlap_tri_range.size(); + TaskParallelSettings settings; + BLI_parallel_range_settings_defaults(&settings); + settings.min_iter_per_thread = 50; + settings.use_threading = intersect_use_threading; + BLI_task_parallel_range( + 0, overlap_tri_range_tot, &data, calc_subdivided_tri_range_func, &settings); +} + +/* Get first index in ov where indexA == t. Assuming sorted on indexA. */ +static int find_first_overlap_index(const TriOverlaps &ov, int t) +{ + Span<BVHTreeOverlap> span = ov.overlap(); + if (span.size() == 0) { + return -1; + } + BVHTreeOverlap bo{t, -1}; + const BVHTreeOverlap *p = std::lower_bound( + span.begin(), span.end(), bo, [](const BVHTreeOverlap &o1, const BVHTreeOverlap &o2) { + return o1.indexA < o2.indexA; + }); + if (p != span.end()) { + return p - span.begin(); + } + return -1; +} + +static CDT_data calc_cluster_subdivided(const CoplanarClusterInfo &clinfo, + int c, + const IMesh &tm, + const TriOverlaps &ov, + const Map<std::pair<int, int>, ITT_value> &itt_map, + IMeshArena *UNUSED(arena)) +{ + constexpr int dbg_level = 0; + BLI_assert(c < clinfo.tot_cluster()); + const CoplanarCluster &cl = clinfo.cluster(c); + /* Make a CDT input with triangles from C and intersects from other triangles in tm. */ + if (dbg_level > 0) { + std::cout << "CALC_CLUSTER_SUBDIVIDED for cluster " << c << " = " << cl << "\n"; + } + /* Get vector itts of all intersections of a triangle of cl with any triangle of tm not + * in cl and not co-planar with it (for that latter, if there were an intersection, + * it should already be in cluster cl). */ + Vector<ITT_value> itts; + Span<BVHTreeOverlap> ovspan = ov.overlap(); + for (int t : cl) { + if (dbg_level > 0) { + std::cout << "find intersects with triangle " << t << " of cluster\n"; + } + int first_i = find_first_overlap_index(ov, t); + if (first_i == -1) { + continue; + } + for (int i = first_i; i < ovspan.size() && ovspan[i].indexA == t; ++i) { + int t_other = ovspan[i].indexB; + if (clinfo.tri_cluster(t_other) != c) { + if (dbg_level > 0) { + std::cout << "use intersect(" << t << "," << t_other << "\n"; + } + std::pair<int, int> key = canon_int_pair(t, t_other); + if (itt_map.contains(key)) { + ITT_value itt = itt_map.lookup(key); + if (itt.kind != INONE && itt.kind != ICOPLANAR) { + itts.append(itt); + if (dbg_level > 0) { + std::cout << " itt = " << itt << "\n"; + } + } + } + } + } + } + /* Use CDT to subdivide the cluster triangles and the points and segs in itts. */ + CDT_data cd_data = prepare_cdt_input_for_cluster(tm, clinfo, c, itts); + do_cdt(cd_data); + return cd_data; +} + +static IMesh union_tri_subdivides(const blender::Array<IMesh> &tri_subdivided) +{ + int tot_tri = 0; + for (const IMesh &m : tri_subdivided) { + tot_tri += m.face_size(); + } + Array<Face *> faces(tot_tri); + int face_index = 0; + for (const IMesh &m : tri_subdivided) { + for (Face *f : m.faces()) { + faces[face_index++] = f; + } + } + return IMesh(faces); +} + +static CoplanarClusterInfo find_clusters(const IMesh &tm, const Array<BoundingBox> &tri_bb) +{ + constexpr int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "FIND_CLUSTERS\n"; + } + CoplanarClusterInfo ans(tm.face_size()); + /* There can be more than one #CoplanarCluster per plane. Accumulate them in + * a Vector. We will have to merge some elements of the Vector as we discover + * triangles that form intersection bridges between two or more clusters. */ + Map<Plane, Vector<CoplanarCluster>> plane_cls; + plane_cls.reserve(tm.face_size()); + for (int t : tm.face_index_range()) { + /* Use a canonical version of the plane for map index. + * We can't just store the canonical version in the face + * since canonicalizing loses the orientation of the normal. */ + Plane tplane = *tm.face(t)->plane; + BLI_assert(tplane.exact_populated()); + tplane.make_canonical(); + if (dbg_level > 0) { + std::cout << "plane for tri " << t << " = " << &tplane << "\n"; + } + /* Assume all planes are in canonical from (see canon_plane()). */ + if (plane_cls.contains(tplane)) { + Vector<CoplanarCluster> &curcls = plane_cls.lookup(tplane); + if (dbg_level > 0) { + std::cout << "already has " << curcls.size() << " clusters\n"; + } + int proj_axis = mpq3::dominant_axis(tplane.norm_exact); + /* Partition `curcls` into those that intersect t non-trivially, and those that don't. */ + Vector<CoplanarCluster *> int_cls; + Vector<CoplanarCluster *> no_int_cls; + for (CoplanarCluster &cl : curcls) { + if (bbs_might_intersect(tri_bb[t], cl.bounding_box()) && + non_trivially_coplanar_intersects(tm, t, cl, proj_axis)) { + if (dbg_level > 1) { + std::cout << "append to int_cls\n"; + } + int_cls.append(&cl); + } + else { + if (dbg_level > 1) { + std::cout << "append to no_int_cls\n"; + } + no_int_cls.append(&cl); + } + } + if (int_cls.size() == 0) { + /* t doesn't intersect any existing cluster in its plane, so make one just for it. */ + if (dbg_level > 1) { + std::cout << "no intersecting clusters for t, make a new one\n"; + } + curcls.append(CoplanarCluster(t, tri_bb[t])); + } + else if (int_cls.size() == 1) { + /* t intersects exactly one existing cluster, so can add t to that cluster. */ + if (dbg_level > 1) { + std::cout << "exactly one existing cluster, " << int_cls[0] << ", adding to it\n"; + } + int_cls[0]->add_tri(t, tri_bb[t]); + } + else { + /* t intersections 2 or more existing clusters: need to merge them and replace all the + * originals with the merged one in `curcls`. */ + if (dbg_level > 1) { + std::cout << "merging\n"; + } + CoplanarCluster mergecl; + mergecl.add_tri(t, tri_bb[t]); + for (CoplanarCluster *cl : int_cls) { + for (int t : *cl) { + mergecl.add_tri(t, tri_bb[t]); + } + } + Vector<CoplanarCluster> newvec; + newvec.append(mergecl); + for (CoplanarCluster *cl_no_int : no_int_cls) { + newvec.append(*cl_no_int); + } + plane_cls.add_overwrite(tplane, newvec); + } + } + else { + if (dbg_level > 0) { + std::cout << "first cluster for its plane\n"; + } + plane_cls.add_new(tplane, Vector<CoplanarCluster>{CoplanarCluster(t, tri_bb[t])}); + } + } + /* Does this give deterministic order for cluster ids? I think so, since + * hash for planes is on their values, not their addresses. */ + for (auto item : plane_cls.items()) { + for (const CoplanarCluster &cl : item.value) { + if (cl.tot_tri() > 1) { + ans.add_cluster(cl); + } + } + } + + return ans; +} + +static bool face_is_degenerate(const Face *f) +{ + const Face &face = *f; + const Vert *v0 = face[0]; + const Vert *v1 = face[1]; + const Vert *v2 = face[2]; + if (v0 == v1 || v0 == v2 || v1 == v2) { + return true; + } + double3 da = v2->co - v0->co; + double3 db = v2->co - v1->co; + double3 dab = double3::cross_high_precision(da, db); + double dab_length_squared = dab.length_squared(); + double err_bound = supremum_dot_cross(dab, dab) * index_dot_cross * DBL_EPSILON; + if (dab_length_squared > err_bound) { + return false; + } + mpq3 a = v2->co_exact - v0->co_exact; + mpq3 b = v2->co_exact - v1->co_exact; + mpq3 ab = mpq3::cross(a, b); + if (ab.x == 0 && ab.y == 0 && ab.z == 0) { + return true; + } + + return false; +} + +/* Data and functions to test triangle degeneracy in parallel. */ +struct DegenData { + const IMesh &tm; +}; + +struct DegenChunkData { + bool has_degenerate_tri = false; +}; + +static void degenerate_range_func(void *__restrict userdata, + const int iter, + const TaskParallelTLS *__restrict tls) +{ + DegenData *data = static_cast<DegenData *>(userdata); + DegenChunkData *chunk_data = static_cast<DegenChunkData *>(tls->userdata_chunk); + const Face *f = data->tm.face(iter); + bool is_degenerate = face_is_degenerate(f); + chunk_data->has_degenerate_tri |= is_degenerate; +} + +static void degenerate_reduce(const void *__restrict UNUSED(userdata), + void *__restrict chunk_join, + void *__restrict chunk) +{ + DegenChunkData *degen_chunk_join = static_cast<DegenChunkData *>(chunk_join); + DegenChunkData *degen_chunk = static_cast<DegenChunkData *>(chunk); + degen_chunk_join->has_degenerate_tri |= degen_chunk->has_degenerate_tri; +} + +/* Does triangle #IMesh tm have any triangles with zero area? */ +static bool has_degenerate_tris(const IMesh &tm) +{ + DegenData degen_data = {tm}; + DegenChunkData degen_chunk_data; + TaskParallelSettings settings; + BLI_parallel_range_settings_defaults(&settings); + settings.userdata_chunk = °en_chunk_data; + settings.userdata_chunk_size = sizeof(degen_chunk_data); + settings.func_reduce = degenerate_reduce; + settings.min_iter_per_thread = 1000; + settings.use_threading = intersect_use_threading; + BLI_task_parallel_range(0, tm.face_size(), °en_data, degenerate_range_func, &settings); + return degen_chunk_data.has_degenerate_tri; +} + +static IMesh remove_degenerate_tris(const IMesh &tm_in) +{ + IMesh ans; + Vector<Face *> new_faces; + new_faces.reserve(tm_in.face_size()); + for (Face *f : tm_in.faces()) { + if (!face_is_degenerate(f)) { + new_faces.append(f); + } + } + ans.set_faces(new_faces); + return ans; +} + +/* This is the main routine for calculating the self_intersection of a triangle mesh. */ +IMesh trimesh_self_intersect(const IMesh &tm_in, IMeshArena *arena) +{ + return trimesh_nary_intersect( + tm_in, 1, [](int UNUSED(t)) { return 0; }, true, arena); +} + +IMesh trimesh_nary_intersect(const IMesh &tm_in, + int nshapes, + std::function<int(int)> shape_fn, + bool use_self, + IMeshArena *arena) +{ + constexpr int dbg_level = 0; + if (dbg_level > 0) { + std::cout << "\nTRIMESH_NARY_INTERSECT nshapes=" << nshapes << " use_self=" << use_self + << "\n"; + for (const Face *f : tm_in.faces()) { + BLI_assert(f->is_tri()); + UNUSED_VARS_NDEBUG(f); + } + if (dbg_level > 1) { + std::cout << "input mesh:\n" << tm_in; + for (int t : tm_in.face_index_range()) { + std::cout << "shape(" << t << ") = " << shape_fn(tm_in.face(t)->orig) << "\n"; + } + } + } +# ifdef PERFDEBUG + perfdata_init(); + double start_time = PIL_check_seconds_timer(); + std::cout << "trimesh_nary_intersect start\n"; +# endif + /* Usually can use tm_in but if it has degenerate or illegal triangles, + * then need to work on a copy of it without those triangles. */ + const IMesh *tm_clean = &tm_in; + IMesh tm_cleaned; + if (has_degenerate_tris(tm_in)) { + if (dbg_level > 0) { + std::cout << "cleaning degenerate triangles\n"; + } + tm_cleaned = remove_degenerate_tris(tm_in); + tm_clean = &tm_cleaned; + if (dbg_level > 1) { + std::cout << "cleaned input mesh:\n" << tm_cleaned; + } + } + /* Temporary, while developing: populate all plane normals exactly. */ + for (Face *f : tm_clean->faces()) { + f->populate_plane(true); + } +# ifdef PERFDEBUG + double clean_time = PIL_check_seconds_timer(); + std::cout << "cleaned, time = " << clean_time - start_time << "\n"; +# endif + Array<BoundingBox> tri_bb = calc_face_bounding_boxes(*tm_clean); +# ifdef PERFDEBUG + double bb_calc_time = PIL_check_seconds_timer(); + std::cout << "bbs calculated, time = " << bb_calc_time - clean_time << "\n"; +# endif + TriOverlaps tri_ov(*tm_clean, tri_bb, nshapes, shape_fn, use_self); +# ifdef PERFDEBUG + double overlap_time = PIL_check_seconds_timer(); + std::cout << "intersect overlaps calculated, time = " << overlap_time - bb_calc_time << "\n"; + +# endif + CoplanarClusterInfo clinfo = find_clusters(*tm_clean, tri_bb); + if (dbg_level > 1) { + std::cout << clinfo; + } +# ifdef PERFDEBUG + double find_cluster_time = PIL_check_seconds_timer(); + std::cout << "clusters found, time = " << find_cluster_time - overlap_time << "\n"; + doperfmax(0, tm_in.face_size()); + doperfmax(1, clinfo.tot_cluster()); + doperfmax(2, tri_ov.overlap().size()); +# endif + /* itt_map((a,b)) will hold the intersection value resulting from intersecting + * triangles with indices a and b, where a < b. */ + Map<std::pair<int, int>, ITT_value> itt_map; + itt_map.reserve(tri_ov.overlap().size()); + calc_overlap_itts(itt_map, *tm_clean, clinfo, tri_ov, arena); +# ifdef PERFDEBUG + double itt_time = PIL_check_seconds_timer(); + std::cout << "itts found, time = " << itt_time - find_cluster_time << "\n"; +# endif + Array<IMesh> tri_subdivided(tm_clean->face_size()); + calc_subdivided_tris(tri_subdivided, *tm_clean, itt_map, clinfo, tri_ov, arena); +# ifdef PERFDEBUG + double subdivided_tris_time = PIL_check_seconds_timer(); + std::cout << "subdivided tris found, time = " << subdivided_tris_time - itt_time << "\n"; +# endif + Array<CDT_data> cluster_subdivided(clinfo.tot_cluster()); + for (int c : clinfo.index_range()) { + cluster_subdivided[c] = calc_cluster_subdivided(clinfo, c, *tm_clean, tri_ov, itt_map, arena); + } +# ifdef PERFDEBUG + double cluster_subdivide_time = PIL_check_seconds_timer(); + std::cout << "subdivided clusters found, time = " + << cluster_subdivide_time - subdivided_tris_time << "\n"; +# endif + for (int t : tm_clean->face_index_range()) { + int c = clinfo.tri_cluster(t); + if (c != NO_INDEX) { + BLI_assert(tri_subdivided[t].face_size() == 0); + tri_subdivided[t] = extract_subdivided_tri(cluster_subdivided[c], *tm_clean, t, arena); + } + else if (tri_subdivided[t].face_size() == 0) { + tri_subdivided[t] = extract_single_tri(*tm_clean, t); + } + } +# ifdef PERFDEBUG + double extract_time = PIL_check_seconds_timer(); + std::cout << "triangles extracted, time = " << extract_time - cluster_subdivide_time << "\n"; +# endif + IMesh combined = union_tri_subdivides(tri_subdivided); + if (dbg_level > 1) { + std::cout << "TRIMESH_NARY_INTERSECT answer:\n"; + std::cout << combined; + } +# ifdef PERFDEBUG + double end_time = PIL_check_seconds_timer(); + std::cout << "triangles combined, time = " << end_time - extract_time << "\n"; + std::cout << "trimesh_nary_intersect done, total time = " << end_time - start_time << "\n"; + dump_perfdata(); +# endif + return combined; +} + +static std::ostream &operator<<(std::ostream &os, const CoplanarCluster &cl) +{ + os << "cl("; + bool first = true; + for (const int t : cl) { + if (first) { + first = false; + } + else { + os << ","; + } + os << t; + } + os << ")"; + return os; +} + +static std::ostream &operator<<(std::ostream &os, const CoplanarClusterInfo &clinfo) +{ + os << "Coplanar Cluster Info:\n"; + for (int c : clinfo.index_range()) { + os << c << ": " << clinfo.cluster(c) << "\n"; + } + return os; +} + +static std::ostream &operator<<(std::ostream &os, const ITT_value &itt) +{ + switch (itt.kind) { + case INONE: + os << "none"; + break; + case IPOINT: + os << "point " << itt.p1; + break; + case ISEGMENT: + os << "segment " << itt.p1 << " " << itt.p2; + break; + case ICOPLANAR: + os << "co-planar t" << itt.t_source; + break; + } + return os; +} + +/** + * Writing the obj_mesh has the side effect of populating verts. + */ +void write_obj_mesh(IMesh &m, const std::string &objname) +{ + /* Would like to use #BKE_tempdir_base() here, but that brings in dependence on kernel library. + * This is just for developer debugging anyway, + * and should never be called in production Blender. */ +# ifdef _WIN_32 + const char *objdir = BLI_getenv("HOME"); +# else + const char *objdir = "/tmp/"; +# endif + if (m.face_size() == 0) { + return; + } + + std::string fname = std::string(objdir) + objname + std::string(".obj"); + std::ofstream f; + f.open(fname); + if (!f) { + std::cout << "Could not open file " << fname << "\n"; + return; + } + + if (!m.has_verts()) { + m.populate_vert(); + } + for (const Vert *v : m.vertices()) { + const double3 dv = v->co; + f << "v " << dv[0] << " " << dv[1] << " " << dv[2] << "\n"; + } + int i = 0; + for (const Face *face : m.faces()) { + /* OBJ files use 1-indexing for vertices. */ + f << "f "; + for (const Vert *v : *face) { + int i = m.lookup_vert(v); + BLI_assert(i != NO_INDEX); + /* OBJ files use 1-indexing for vertices. */ + f << i + 1 << " "; + } + f << "\n"; + ++i; + } + f.close(); +} + +# ifdef PERFDEBUG +struct PerfCounts { + Vector<int> count; + Vector<const char *> count_name; + Vector<int> max; + Vector<const char *> max_name; +}; + +static PerfCounts *perfdata = nullptr; + +static void perfdata_init(void) +{ + perfdata = new PerfCounts; + + /* count 0. */ + perfdata->count.append(0); + perfdata->count_name.append("Non-cluster overlaps"); + + /* count 1. */ + perfdata->count.append(0); + perfdata->count_name.append("intersect_tri_tri calls"); + + /* count 2. */ + perfdata->count.append(0); + perfdata->count_name.append("tri tri intersects decided by filter plane tests"); + + /* count 3. */ + perfdata->count.append(0); + perfdata->count_name.append("tri tri intersects decided by exact plane tests"); + + /* count 4. */ + perfdata->count.append(0); + perfdata->count_name.append("final non-NONE intersects"); + + /* max 0. */ + perfdata->max.append(0); + perfdata->max_name.append("total faces"); + + /* max 1. */ + perfdata->max.append(0); + perfdata->max_name.append("total clusters"); + + /* max 2. */ + perfdata->max.append(0); + perfdata->max_name.append("total overlaps"); +} + +static void incperfcount(int countnum) +{ + perfdata->count[countnum]++; +} + +static void bumpperfcount(int countnum, int amt) +{ + perfdata->count[countnum] += amt; +} + +static void doperfmax(int maxnum, int val) +{ + perfdata->max[maxnum] = max_ii(perfdata->max[maxnum], val); +} + +static void dump_perfdata(void) +{ + std::cout << "\nPERFDATA\n"; + for (int i : perfdata->count.index_range()) { + std::cout << perfdata->count_name[i] << " = " << perfdata->count[i] << "\n"; + } + for (int i : perfdata->max.index_range()) { + std::cout << perfdata->max_name[i] << " = " << perfdata->max[i] << "\n"; + } + delete perfdata; +} +# endif + +}; // namespace blender::meshintersect + +#endif // WITH_GMP |