Issues when using Möller–Trumbore intersection algorithm to figure out if a vertex is within a shape

Question

I am making a program in C++ that creates a grid of vertices an then imports an stl image. I then use the Möller–Trumbore intersection algorithm on every vertex of my grid to count the number of intersections a ray cast from it makes with the triangles of the stl. If the vertex is inside the stl object, I delete it.

My problem is that some, not all, of the vertices located at the edges of my object do not get marked as inside the object.

The method I'm using is the one explained in this post with a few additions so that edge and vertex gits only mark 1 intersection.

Here is my grid when I delete all of the vertices which were detected as inside a 20x20x20 cube stl image.

Here's a top-down view

This is mesh.cpp, where I import the image count the intersections

#include "mesh.h"

using namespace std;

bool Mesh::RayIntersectsTriangle(vector3D<double> rayOrigin, 
                           vector3D<double> rayVector, 
                           Triangle* inTriangle,
                           vector3D<double>& outIntersectionPoint)
{
    const double EPSILON = 0.00000000001;
    vector3D<double> vertex0 = inTriangle->corner[0];
    vector3D<double> vertex1 = inTriangle->corner[1];  
    vector3D<double> vertex2 = inTriangle->corner[2];
    vector3D<double> edge1, edge2, h, s, q;
    double a, f, u, v;
    edge1 = vertex1 - vertex0;
    edge2 = vertex2 - vertex0;
    h = rayVector ^ edge2; //^ is cross product
    a = edge1*h; // * is dot product

    if (a > -EPSILON && a < EPSILON)
        return false;    // This ray is parallel to this triangle.

    f = 1.0 / a;
    s = rayOrigin - vertex0;
    u = f * s*h;

    if (u < 0.0 || u > 1.0)
        return false;

    q = s^edge1;
    v = f * rayVector*q;

    if (v < 0.0 || u + v > 1.0)
        return false;

    // At this stage we can compute t to find out where the intersection point is on the line.
    double t = f * edge2*q;

    if (t > EPSILON) // ray intersection
    {
        outIntersectionPoint = rayOrigin + rayVector * t;
        return true;
    }
    else // This means that there is a line intersection but not a ray intersection.
        return false;
}

//Flatten along Z axis so that I can only take into account triangles that are relatively close to my target vertex
void Mesh::check_bounding_box_collision(vector3D<double> target){
    for (const auto & trig : triangles){
        double xmax = max(max(trig.corner[0].x, trig.corner[1].x), trig.corner[2].x);
        double ymax = max(max(trig.corner[0].y, trig.corner[1].y), trig.corner[2].y);
        double xmin = min(min(trig.corner[0].x, trig.corner[1].x), trig.corner[2].x);
        double ymin = min(min(trig.corner[0].y, trig.corner[1].y), trig.corner[2].y);

        if ((target.x <= xmax && target.x >= xmin) && (target.y <= ymax && target.y >= ymin)){
            the_chosen_ones.push_back(trig);
        }

    }
}

void Mesh::displace(vector3D<double> distance){
    for (int i = 0; i < triangles.size(); ++i){
        triangles.at(i).corner[0] += distance;
        triangles.at(i).corner[1] += distance;
        triangles.at(i).corner[2] += distance;
    }
}

unsigned int Mesh::count_intersections(vector3D<double> target){
    unsigned int intersections = 0;
    int blocked = 0;
    vector3D<double> iPoint;
    set<vector3D<double>, Vector3DComparator> intersection_points;
    the_chosen_ones.clear();
    check_bounding_box_collision(target);
    for (auto& triangle : the_chosen_ones){
        if (RayIntersectsTriangle(target, {0.0, 0.0, 1.0}, &triangle, iPoint)) {//Shoot up beam along flattened axis
            if (intersection_points.find(iPoint) == intersection_points.end()) {//If this intersection point exists in the list, I've hit an edge, so I dont count the intersection for a second time.
                intersections++;
                intersection_points.insert(iPoint);
            }
        }
    }

    return intersections;
}


void Mesh::load_stl(string filename){
    stl_reader::StlMesh<double, unsigned int> mesh(filename);
    cout << "Loading " << filename << "..." << endl;

    for(size_t itri = 0; itri < mesh.num_tris(); ++itri) {
        Triangle t;
        for(size_t icorner = 0; icorner < 3; ++icorner) {
            const double* c = mesh.tri_corner_coords (itri, icorner);
            t.corner[icorner] = {c[0], c[1], c[2]};
        }
    
        const double* n = mesh.tri_normal (itri);
        t.normal = {n[0], n[1], n[2]};
        triangles.push_back(t);
    }
}

And here is mesh.h

#include <iostream>
#include <vector>
#include <set>
#include "stl_reader.h"
#include "vector3D.h"

using namespace std;

struct Triangle{
    vector3D<double> corner[3];
    vector3D<double> normal;
};

struct Vector3DComparator {
    bool operator()(const vector3D<double>& lhs, const vector3D<double>& rhs) const {
        const double TOLERANCE = 0.00000000000001;
        if (fabs(lhs.x - rhs.x) > TOLERANCE) {
            return lhs.x < rhs.x;
        } else if (fabs(lhs.y - rhs.y) > TOLERANCE) {
            return lhs.y < rhs.y;
        } else {
            return lhs.z < rhs.z;
        }
    }
};

class Mesh{
    private:
    vector<Triangle> triangles;
    vector<Triangle> the_chosen_ones;

    bool RayIntersectsTriangle(vector3D<double> rayOrigin, 
                           vector3D<double> rayVector, 
                           Triangle* inTriangle,
                           vector3D<double>& outIntersectionPoint);
    
    void check_bounding_box_collision(vector3D<double> target);

    public:
    unsigned int count_intersections(vector3D<double> target);
    void displace(vector3D<double> distance);
    void load_stl(string filename);
};

The cube im importing is just a genertic 20x20x20 cube stl