tableLattice.h

Go to the documentation of this file.

#pragma once
 
#include <memory>
#include <vector>
#include <utility>
 
#include "opatIO.h"
#include "indexVector.h"
 
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>

namespace opat::lattice {

    typedef boost::numeric::ublas::matrix<double> bmat;
    typedef boost::numeric::ublas::vector<double> bvec;

    struct Simplex {
        std::size_t ID = static_cast<std::size_t>(-1); 
        std::vector<double> barycentricWeights; 
    };

    enum class InterpolationType {
        Linear,    
        Quadratic, 
        Cubic      
    };

    class TableLattice {
    public:
        explicit TableLattice(const opat::OPAT& opat);
        TableLattice(const opat::OPAT& opat, const InterpolationType& interpolationType);
 
        TableLattice(const TableLattice&) = default;
        TableLattice& operator=(const TableLattice&) = delete;
        TableLattice(TableLattice&&) noexcept = default;
        TableLattice& operator=(TableLattice&&) noexcept = delete;
        ~TableLattice() = default;

        [[nodiscard]] DataCard get(const FloatIndexVector& indexVector) const;
        [[nodiscard]] InterpolationType getInterpolationType() const;

        void setInterpolationType(InterpolationType interpolationType);

        void dumpTriangulationToAscii(const std::string &points_file, const std::string &simplices_file) const;
    private:
        const opat::OPAT &m_opat; 
        std::size_t m_indexVectorSize{}; 
        InterpolationType m_interpolationType{InterpolationType::Linear}; 
        std::vector<FloatIndexVector> m_indexVectors; 
        std::vector<std::vector<double>> m_axisValues; 
        std::size_t m_numCorners{}; 
        std::vector<std::vector<std::size_t>> m_simplices; 
        std::vector<std::vector<std::size_t>> m_simplexAdjacency; 

        void initialize();

        void buildDelaunay();
        [[nodiscard]] Simplex findContainingSimplex(const FloatIndexVector& indexVector) const;

        void validateIndexVector(const FloatIndexVector &indexVector) const;

        std::vector<double> calculateBarycentricWeights(const FloatIndexVector& queryPoint, const std::vector<FloatIndexVector>& simplexActualVertices) const;
 
        mutable Simplex m_lastFoundSimplex; 
 
    };

    bvec solveLinearSystem(bmat A, bvec b);
}