volume_geometry.hpp

#pragma once
 
#include "volume_complex.hpp"
 
#include <Eigen/Core>
 
#include <map>
#include <optional>
#include <vector>
 
namespace wmtk::components::tracking::volume {
 
using Positions = std::vector<Eigen::Vector3d>;
 
int orient(
    const Eigen::Vector3d& a,
    const Eigen::Vector3d& b,
    const Eigen::Vector3d& c,
    const Eigen::Vector3d& d);
 
bool inside(
    const Eigen::Vector3d& p,
    const Eigen::Vector3d& a,
    const Eigen::Vector3d& b,
    const Eigen::Vector3d& c,
    const Eigen::Vector3d& d);
Eigen::Vector4d barycentric(
    const Eigen::Vector3d& p,
    const Eigen::Vector3d& a,
    const Eigen::Vector3d& b,
    const Eigen::Vector3d& c,
    const Eigen::Vector3d& d);
Eigen::Vector4d exact_barycentric(
    const Eigen::Vector3d& p,
    const Eigen::Vector3d& a,
    const Eigen::Vector3d& b,
    const Eigen::Vector3d& c,
    const Eigen::Vector3d& d);
struct Outside
{
    int faces = 0;
    double distance = 0; // the largest
};
Outside outside_of(
    const Eigen::Vector3d& p,
    const Eigen::Vector3d& a,
    const Eigen::Vector3d& b,
    const Eigen::Vector3d& c,
    const Eigen::Vector3d& d);
struct Location
{
    size_t tet;
    Eigen::Vector4d b;
    double min_coordinate;
};
template <typename Accept>
std::optional<Location>
locate_in(const Tets& T, const Positions& X, const Eigen::Vector3d& p, const Accept& accept)
{
    if (!p.allFinite()) return std::nullopt;
    std::vector<size_t> within_rounding;
    for (size_t k = 0; k < T.size(); ++k) {
        if (!accept(k)) continue;
        const auto& t = T[k];
        // only a tet whose box (grown by a rounding) holds p can hold it or be a rounding off
        Eigen::Vector3d lo = X[t[0]], hi = lo;
        for (int c = 1; c < 4; ++c) {
            lo = lo.cwiseMin(X[t[c]]);
            hi = hi.cwiseMax(X[t[c]]);
        }
        const double size = (hi - lo).norm(), eps = 1e-9 * size;
        if (((p.array() < lo.array() - eps) || (p.array() > hi.array() + eps)).any()) continue;
        const Outside o = outside_of(p, X[t[0]], X[t[1]], X[t[2]], X[t[3]]);
        if (o.faces == 0) {
            Eigen::Vector4d b = barycentric(p, X[t[0]], X[t[1]], X[t[2]], X[t[3]]);
            if (!(b.minCoeff() >= -1e-9 && b.maxCoeff() <= 1 + 1e-9)) {
                b = exact_barycentric(p, X[t[0]], X[t[1]], X[t[2]], X[t[3]]);
            }
            b = b.cwiseMax(0.);
            return Location{k, b / b.sum(), 0.};
        }
        if (o.distance <= eps) within_rounding.push_back(k);
    }
    std::optional<Location> out;
    if (!within_rounding.empty()) {
        for (const size_t k : within_rounding) {
            const auto& t = T[k];
            const Eigen::Vector4d b = exact_barycentric(p, X[t[0]], X[t[1]], X[t[2]], X[t[3]]);
            if (!out || b.minCoeff() > out->min_coordinate) out = Location{k, b, b.minCoeff()};
        }
    }
    if (out) {
        out->b = out->b.cwiseMax(0.);
        out->b /= out->b.sum();
    }
    return out;
}
 
class StarPatch
{
public:
    static std::optional<StarPatch> make(const Tets& T, const Positions& P);
 
    const Tets& tets() const { return m_T; }
    const Positions& positions() const { return m_P; }
    const std::vector<Tri>& faces() const { return m_F; }
    double scale() const { return m_scale; }
    std::optional<Vid> boundary_apex() const;
    std::optional<Vid> interior_apex() const;
 
    bool sees_all(const Eigen::Vector3d& s) const;
    bool set_centre(const Eigen::Vector3d& s);
    const Eigen::Vector3d& centre() const { return m_s; }
    double margin() const;
 
    struct Exit
    {
        size_t face;
        Eigen::Vector3d w;
        double t;
    };
    Exit exit(const Eigen::Vector3d& x) const;
    struct Located
    {
        size_t tet;
        Eigen::Vector4d b;
    };
    Located locate(const Eigen::Vector3d& y) const;
    std::optional<size_t> index_of(const Tet& t) const;
    size_t bytes() const;
 
    std::vector<Eigen::Vector3d> normals, corners;
 
private:
    Tets m_T;
    Positions m_P;
    std::vector<Tri> m_F, m_Fo;
    std::vector<std::pair<Tet, uint32_t>> m_index; // by sorted vertices, sorted
    double m_scale = 0;
    Eigen::Vector3d m_s = Eigen::Vector3d::Zero();
};
 
std::optional<std::pair<Eigen::Vector3d, double>> chebyshev(
    const std::vector<Eigen::Vector3d>& N,
    const std::vector<double>& b,
    const Eigen::Vector3d& x0,
    double L,
    double t_max);
 
std::optional<Eigen::Vector3d> kernel(const std::vector<const StarPatch*>& patches);
std::optional<Eigen::Vector3d> apex_centre(const StarPatch& S);
std::optional<Eigen::Vector3d> centre_of(const StarPatch& S, bool* thin = nullptr);
 
std::optional<Positions> spherical(const Tets& T, const Positions& P, double fat = 0.05);
std::optional<Positions> bowl(const Tets& T, const Positions& P);
std::optional<Positions> convex(const Tets& T, const Positions& P);
 
} // namespace wmtk::components::tracking::volume