Generate block mesh using binary stl file

#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <queue>
using namespace std;

typedef std::vector<double>          VecDbl_t;
typedef std::vector<int>             VecInt_t;
typedef std::vector<VecDbl_t>        VecVecDbl_t;
typedef std::vector<VecInt_t>        VecVecInt_t;

struct Vertex {
    double x, y, z;
};


struct Triangle {
    double x1, y1, z1;
    double x2, y2, z2;
    double x3, y3, z3;
};

struct BlockMesh {
    std::vector<Vertex> vertices;
    VecVecInt_t blocks;
};


double distanceToTriangle(double x, double y, double z, const Triangle& triangle) {
    double ux = triangle.x2 - triangle.x1;
    double uy = triangle.y2 - triangle.y1;
    double uz = triangle.z2 - triangle.z1;
    double vx = triangle.x3 - triangle.x1;
    double vy = triangle.y3 - triangle.y1;
    double vz = triangle.z3 - triangle.z1;
    double wx = x - triangle.x1;
    double wy = y - triangle.y1;
    double wz = z - triangle.z1;

    double uu = ux * ux + uy * uy + uz * uz;
    double uv = ux * vx + uy * vy + uz * vz;
    double vv = vx * vx + vy * vy + vz * vz;
    double uw = ux * wx + uy * wy + uz * wz;
    double vw = vx * wx + vy * wy + vz * wz;

    double denom = uu * vv - uv * uv;
    double s = (uv * vw - vv * uw) / denom;
    double t = (uv * uw - uu * vw) / denom;

    double dist = std::sqrt(uu * vv - uv * uv) / std::sqrt(denom);

    if (s >= 0 && t >= 0 && s + t <= 1) {
        return dist;
    }
    else {
        double dist1 = std::sqrt((x - triangle.x1) * (x - triangle.x1) + (y - triangle.y1) * (y - triangle.y1) + (z - triangle.z1) * (z - triangle.z1));
        double dist2 = std::sqrt((x - triangle.x2) * (x - triangle.x2) + (y - triangle.y2) * (y - triangle.y2) + (z - triangle.z2) * (z - triangle.z2));
        double dist3 = std::sqrt((x - triangle.x3) * (x - triangle.x3) + (y - triangle.y3) * (y - triangle.y3) + (z - triangle.z3) * (z - triangle.z3));
        return std::min({ dist1, dist2, dist3 });
    }
}

std::vector<Triangle> readSTL(const std::string& stlFile) {
    std::vector<Triangle> triangles;

    std::ifstream file(stlFile, std::ios::in | std::ios::binary);

    if (!file) {
        std::cerr << "Error opening file " << stlFile << std::endl;
        exit(EXIT_FAILURE);
    }

    file.seekg(80, std::ios::beg);
    uint32_t numTriangles;
    file.read((char*)&numTriangles, sizeof(numTriangles));
    for (uint32_t i = 0; i < numTriangles; i++) {

        file.seekg(3 * sizeof(float), std::ios::cur);
        float x1, y1, z1, x2, y2, z2, x3, y3, z3;
        file.read((char*)&x1, sizeof(float));
        file.read((char*)&y1, sizeof(float));
        file.read((char*)&z1, sizeof(float));
        file.read((char*)&x2, sizeof(float));
        file.read((char*)&y2, sizeof(float));
        file.read((char*)&z2, sizeof(float));
        file.read((char*)&x3, sizeof(float));
        file.read((char*)&y3, sizeof(float));
        file.read((char*)&z3, sizeof(float));
        Triangle tri = { x1, y1, z1, x2, y2, z2, x3, y3, z3 };
        triangles.push_back(tri);
        file.seekg(sizeof(uint16_t), std::ios::cur);
    }
    file.close();

    std::cout << "Number of triangles: " << triangles.size() << std::endl;
    for (uint32_t i = 0; i < triangles.size(); i++) {
        std::cout << "Triangle " << i + 1 << ":" << std::endl;
        std::cout << "Vertex 1: " << triangles[i].x1 << ", " << triangles[i].y1 << ", " << triangles[i].z1 << std::endl;
        std::cout << "Vertex 2: " << triangles[i].x2 << ", " << triangles[i].y2 << ", " << triangles[i].z2 << std::endl;
        std::cout << "Vertex 3: " << triangles[i].x3 << ", " << triangles[i].y3 << ", " << triangles[i].z3 << std::endl;
    }

    return triangles;
}

BlockMesh createBlockMesh(const std::string& stlFile, double dx, double maxDist) {

    std::vector<Triangle> triangles = readSTL(stlFile);

    double minX = std::numeric_limits<double>::max();
    double minY = std::numeric_limits<double>::max();
    double minZ = std::numeric_limits<double>::max();
    double maxX = std::numeric_limits<double>::lowest();
    double maxY = std::numeric_limits<double>::lowest();
    double maxZ = std::numeric_limits<double>::lowest();
    for (const auto& triangle : triangles) {
        minX = std::min(minX, std::min(triangle.x1, std::min(triangle.x2, triangle.x3)));
        minY = std::min(minY, std::min(triangle.y1, std::min(triangle.y2, triangle.y3)));
        minZ = std::min(minZ, std::min(triangle.z1, std::min(triangle.z2, triangle.z3)));
        maxX = std::max(maxX, std::max(triangle.x1, std::max(triangle.x2, triangle.x3)));
        maxY = std::max(maxY, std::max(triangle.y1, std::max(triangle.y2, triangle.y3)));
        maxZ = std::max(maxZ, std::max(triangle.z1, std::max(triangle.z2, triangle.z3)));
    }

    double blockX = std::ceil((maxX - minX) / dx);
    double blockY = std::ceil((maxY - minY) / dx);
    double blockZ = std::ceil((maxZ - minZ) / dx);

    std::vector<Vertex> vertices;
    VecVecInt_t blocks(blockX * blockY * blockZ);
    for (int i = 0; i < blockX; i++) {
        for (int j = 0; j < blockY; j++) {
            for (int k = 0; k < blockZ; k++) {
                double x = minX + i * dx;
                double y = minY + j * dx;
                double z = minZ + k * dx;
                bool inside = false;
                for (const auto& triangle : triangles) {
                    double dist = distanceToTriangle(x, y, z, triangle);
                    if (dist < maxDist) {
                        inside = !inside;
                    }
                }
                if (inside) {
                    int idx = k + j * blockZ + i * blockZ * blockY;
                    blocks[idx].push_back(vertices.size());
                    vertices.push_back({ x, y, z });
                }
            }
        }
    }

    BlockMesh blockMesh = { vertices, blocks };
    return blockMesh;
}


int main() {

    BlockMesh mesh = createBlockMesh("cylinder.stl", 0.1,0.001);
    return 0;
}