使用Möller–Trumbore算法判断顶点是否在模型内的异常问题

我正在用C++编写一个程序，该程序生成顶点网格并导入STL模型。随后对网格中的每个顶点使用Möller–Trumbore相交算法，统计从顶点射出的射线与STL模型三角形的相交次数，若顶点在STL模型内部则将其删除。

我的问题是：位于模型边缘的部分（非全部）顶点未被标记为模型内部顶点。

我采用的方法参考了Stack Overflow上的「判断点是否在3D网格内」的算法，并做了补充，确保边缘和顶点相交仅计为1次。

mesh.cpp 文件 #

#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);
    }
}

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);
};

我导入的是一个标准20x20x20的立方体STL模型。

内容的提问来源于stack exchange，提问作者dimflix 123