point_map.cpp

point_map.cpp

#include "libvcad/tree/leaf/point_map.h"

#include <fstream>

#include "libvcad/utils/math_utils.h"

using namespace std;
using namespace glm;
using namespace MathUtils;

PointMap::PointMap(const string& tet_mesh_path, const string& displacement_field_path,
                   vector<string> functions, vector<uint8_t> materials, bool prob_mode) :
    m_tet_mesh_path(tet_mesh_path), m_displacement_field_path(displacement_field_path),
    m_function_strings(functions), m_materials(materials), m_prob_mode(prob_mode)
{}

shared_ptr<Node> PointMap::clone()
{
    auto clone = make_shared<PointMap>();
    clone->m_tet_mesh_path = m_tet_mesh_path;
    clone->m_displacement_field_path = m_displacement_field_path;
    clone->m_function_strings = m_function_strings;
    clone->m_materials = m_materials;
    clone->m_prob_mode = m_prob_mode;
    // Shallow copy the AABB tree
    // This is okay because the AABB tree is not modified, only read after being prepared
    clone->m_aabb_tree = m_aabb_tree;

    // Build new evaluators based on the new class members
    clone->prepare_evaluators();
    return clone;
}

void PointMap::prepare()
{
    // Prepare the evaluators
    prepare_evaluators();

    // Create the AABB tree
    m_aabb_tree = make_shared<AABBTree>();

    // Load the tetrahedrons and displacement field from the files
    LoadTetrahedronsFromFile(m_tet_mesh_path, m_displacement_field_path);
}


void PointMap::prepare_evaluators()
{
    m_evaluators.clear(); // Delete the old evaluators

    // {"x", "y", "z", "rho", "phic", "r", "theta", "phis", "dx", "dy", "dz", "len"}
    const std::unordered_map<std::string, double*> variable_map =
            {
                    { "x", &m_x }, { "y", &m_y }, { "z", &m_z },
                    { "rho", &m_rho }, { "phic", &m_phic },
                    { "r", &m_r }, { "theta", &m_theta }, { "phis", &m_phis },
                    { "dx", &m_dx }, { "dy", &m_dy }, { "dz", &m_dz },
                    { "len", &m_len }
            };

    // For each function, compile the evaluator
    m_evaluators.resize(m_function_strings.size());
    for (int i = 0; i < m_function_strings.size(); ++i)
    {
        auto evaluator = make_shared<ExpressionEvaluator>();
        evaluator->compile(m_function_strings[i], variable_map);
        m_evaluators[i] = evaluator;
    }
}

double PointMap::evaluate(double x, double y, double z)
{
    SC_Point query_point(x, y, z);
    std::vector<TetrahedronPrimitive<glm::vec3>::Id> intersected_tets_ids;

    m_aabb_tree->all_intersected_primitives(query_point, std::back_inserter(intersected_tets_ids));

    for (auto& tetrahedron : intersected_tets_ids)
    {
        if(tetrahedron->is_point_inside(vec3(x, y, z)))
            return -1.0; // Point is inside the tetrahedrons
    }
    return 1.0; // Point is outside the tetrahedrons
}

vec3 PointMap::compute_displacement(double x, double y, double z)
{
    SC_Point query_point(x, y, z);
    std::vector<TetrahedronPrimitive<glm::vec3>::Id> intersected_tets_ids;

    // Determine what tetrahedrons the point *might* intersect
    m_aabb_tree->all_intersected_primitives(query_point, std::back_inserter(intersected_tets_ids));

    // Iterate over the potential intersected tetrahedrons and check if the point is inside
    for (auto& tetrahedron : intersected_tets_ids)
    {
        if(tetrahedron->is_point_inside(vec3(x, y, z)))
            return tetrahedron->sample(vec3(x, y, z)); // Interpolate the displacement at the point within the tetrahedron
    }
    return vec3(0.0f); // Point is outside the tetrahedrons
}

uint8_t PointMap::material(double x, double y, double z)
{
    auto distribution = this->distribution(x, y, z);
    return compute_new_material(distribution, m_prob_mode);
}

vector<uint8_t> PointMap::material_list()
{
    return m_materials;
}

unordered_map<uint8_t, float> PointMap::distribution(double x, double y, double z)
{
    if(evaluate(x, y, z) > 0.0)
        return {};

    // Compute the interpolated displacement of the point
    auto displacement = compute_displacement(x, y, z);
    auto length = glm::length(displacement);

    // Cartesian (x, y, z))
    m_x = x;
    m_y = y;
    m_z = z;

    // Cylindrical (rho, phi, z)
    m_rho = sqrt(x * x + y * y);
    m_phic = atan2(y, x);

    // Spherical (r, theta, phi)
    m_r = sqrt(x * x + y * y + z * z);
    m_theta = atan2(y, x);
    m_phis = acos(z / m_r);

    // Displacement
    m_dx = displacement.x;
    m_dy = displacement.y;
    m_dz = displacement.z;

    // Length
    m_len = length;

    // Evaluate the functions
    unordered_map<uint8_t, float> probabilities;
    for (int i = 0; i < m_function_strings.size(); ++i)
    {
        auto value = m_evaluators[i]->evaluate();
        probabilities[m_materials[i]] = value;
    }

    // Normalize the probabilities
    float sum = 0.0;
    for(auto& [material, probability] : probabilities) sum += probability;
    for(auto& [material, probability] : probabilities) probability /= sum;

    return probabilities;
}

std::vector<std::string> PointMap::Split(const string &s, char delimiter)
{
    std::vector<std::string> tokens;
    std::string token;
    std::istringstream tokenStream(s);
    while (std::getline(tokenStream, token, delimiter))
        tokens.push_back(token);
    return tokens;
}

void PointMap::LoadTetrahedronsFromFile(const std::string& inp_filename, const std::string& point_map_filename)
{
    std::vector<ValuedPoint> points;

    // Open the inp file
    ifstream inp_file(inp_filename);
    if(!inp_file.is_open()) { throw runtime_error("Could not open inp file: " + inp_filename); }
    string inp_file_line;

    // Open the point map file
    ifstream point_map_file(point_map_filename);
    if(!point_map_file.is_open()) { throw runtime_error("Could not open point map file: " + point_map_filename); }
    string point_map_file_line;

    // Skip lines in the inp file until we reach the *NODE section
    while(getline(inp_file, inp_file_line)) { if(inp_file_line.find("*NODE") != string::npos) { break; }}

    // Read the nodes until we "**"
    while(getline(inp_file, inp_file_line))
    {
        if(inp_file_line.find("**") != string::npos) { break; }

        // Split the inp_file_line into tokens
        vector<string> tokens = Split(inp_file_line, ',');

        // Parse the tokens
        if(tokens.size() != 4) { throw runtime_error("Invalid line in inp file: " + inp_filename); }

        // Convert the tokens to doubles
        double x = stod(tokens[1]);
        double y = stod(tokens[2]);
        double z = stod(tokens[3]);

        // Add the point to the list of points
        points.push_back(ValuedPoint{vec3(x, y, z), vec3(0.0f)});
    }

    // Now we read the point map file and populate the displacement values
    size_t point_index = 0;
    while(getline(point_map_file, point_map_file_line))
    {
        // Split the line into tokens
        vector<string> tokens = Split(point_map_file_line, ',');

        // Parse the tokens
        if(tokens.size() != 6) { throw runtime_error("Invalid line in displacement file: " + point_map_filename); }

        // Convert the tokens to doubles (0-2 are XYZ)
        double dx = stod(tokens[3]);
        double dy = stod(tokens[4]);
        double dz = stod(tokens[5]);

        // Update the displacement field
        points[point_index].displacement = vec3(dx, dy, dz);
        ++point_index;
    }

    // Going back the INP file, skip lines in the inp_file until we reach the *ELEMENT section
    while(getline(inp_file, inp_file_line)) { if(inp_file_line.find("*ELEMENT") != string::npos) { break; } }

    // Read the tetrahedrons until we reach the end of the inp_file
    while(getline(inp_file, inp_file_line))
    {
        // Split the inp_file_line into tokens
        vector<string> tokens = Split(inp_file_line, ',');

        // Parse the tokens
        if(tokens.size() != 5) { throw runtime_error("Invalid line in inp file: " + inp_filename); }

        // Convert the tokens to integers
        size_t a_idx = stoul(tokens[1]) - 1;
        size_t b_idx = stoul(tokens[2]) - 1;
        size_t c_idx = stoul(tokens[3]) - 1;
        size_t d_idx = stoul(tokens[4]) - 1;

        // Get the vertices of the tetrahedron
        const vec3& a = points[a_idx].coords;
        const vec3& b = points[b_idx].coords;
        const vec3& c = points[c_idx].coords;
        const vec3& d = points[d_idx].coords;

        const vec3& a_disp = points[a_idx].displacement;
        const vec3& b_disp = points[b_idx].displacement;
        const vec3& c_disp = points[c_idx].displacement;
        const vec3& d_disp = points[d_idx].displacement;

        // Add to the AABB tree
        m_aabb_tree->insert(new Tetrahedron<glm::vec3>(a, b, c, d, a_disp, b_disp, c_disp, d_disp));
    }
}