C++

Forum

 
Point In Triangle

Cplusc (457)
I am trying to generate points with equal spacing inside complex geometry. I have the STL file representing the triangulation of the surface of the complex geometry. I read the geometry, take the maximum and minimum for x,y and z coordinate. Using these maximums and minimums, I make a bounding box encompassing the complex geometry . Then, I generate points with equal spacing inside complex geometry. The next step is to identify if a points is inside complex geometry or not. To do so, I cast a line in +z direction from each point inside bounding box. If a line intersects with even number of triangles, it will be labeled as solid and if the line intersects with odd number of triangle it will be inside of the complex geometry and labeled as fluid. The only issue is when a line intersects edge or vertex. A vertex can be shared at most by 3 triangles and a edge can be at most shared with 2 triangle while in it should be considered one intersection in both scenario.
I have a few questions here.
First of all, I have doLinesIntersect function for both 2D and 3D and the same for isPointInsideTriangle. As a test case I considered a cylinder which I made in COMSOL and the result for 2D functions are more reasonable than 3D ones.

Second, I doubt if the functions for checking if 2 lines intersec and PointInsideTriangle are correct.

here is the code


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
 


bool Geometry::doLinesIntersect2D(const Line& line1, const Line& line2) {

    double x1 = line1.start.x;
    double y1 = line1.start.y;
    double z1 = line1.start.z;
    double x2 = line1.end.x;
    double y2 = line1.end.y;
    double z2 = line1.end.z;

    double x3 = line2.start.x;
    double y3 = line2.start.y;
    double z3 = line2.start.z;
    double x4 = line2.end.x;
    double y4 = line2.end.y;
    double z4 = line2.end.z;

    double ua = ((x4 - x3) * (y1 - y3) - (y4 - y3) * (x1 - x3)) / ((y4 - y3) * (x2 - x1) - (x4 - x3) * (y2 - y1));
    double ub = ((x2 - x1) * (y1 - y3) - (y2 - y1) * (x1 - x3)) / ((y4 - y3) * (x2 - x1) - (x4 - x3) * (y2 - y1));

    if (ua >= 0 && ua <= 1 && ub >= 0 && ub <= 1) {
        return true;
    }

    return false;
}


int Geometry::countIntersectingTriangles(Line& line) {
    const double threshold = 1e-7;
    int count = 0;
    for (const auto& triangle : triangles) {
        Line edge1, edge2, edge3;
        edge1.start = { triangle.x1, triangle.y1, triangle.z1 };
        edge1.end = { triangle.x2, triangle.y2, triangle.z2 };
        edge2.start = { triangle.x2, triangle.y2, triangle.z2 };
        edge2.end = { triangle.x3, triangle.y3, triangle.z3 };
        edge3.start = { triangle.x3, triangle.y3, triangle.z3 };
        edge3.end = { triangle.x1, triangle.y1, triangle.z1 };

        if (doLinesIntersect3D(line, edge1) || doLinesIntersect3D(line, edge2) || doLinesIntersect3D(line, edge3)) {
            count++;
            continue;
        }

        double px = line.start.x;
        double py = line.start.y;
        double pz = line.start.z;

        double nx = triangle.nx;
        double ny = triangle.ny;
        double nz = triangle.nz;

        double qx = triangle.x1;
        double qy = triangle.y1;
        double qz = triangle.z1;

        double t = ((qx - px) * nx + (qy - py) * ny + (qz - pz) * nz)/
            ((line.end.x - line.start.x) * nx + (line.end.y - line.start.y) * ny + (line.end.z - line.start.z) * nz);

        if (t >= 0 && t <= 1) {
            double ix = px + t * (line.end.x - line.start.x);
            double iy = py + t * (line.end.y - line.start.y);
            double iz = pz + t * (line.end.z - line.start.z);

           if (isPointInsideTriangle3D(ix, iy, iz, triangle) && t > threshold && t < 1 - threshold) {
                //if (isPointInsideTriangle2D(ix, iy, triangle) && t > threshold && t < 1 - threshold) {
                count++;
            }
        }
    }
    return count;
}



Last edited on
helios (17574)
complex geometry
To avoid confusion just say "geometry", or "arbitrary geometry" if you want to emphasize that it's complicated. "Complex geometry" has a specific meaning in mathematics.

vertex can be shared at most by 3 triangles
Is this a property that you're assuming of your input? Because there's plenty of shapes that have more than three triangles sharing the same vertex. Just on the platonic solids there's the octahedron and the icosahedron.

http://www.catb.org/~esr/faqs/smart-questions.html#code
Don't ask others to debug your broken code without giving a hint what sort of problem they should be searching for. Posting a few hundred lines of code, saying "it doesn't work", will get you ignored. Posting a dozen lines of code, saying "after line 7 I was expecting to see <x>, but <y> occurred instead" is much more likely to get you a response.

If you're not sure about the math just use a linear algebra library to solve the linear equations. Don't ask people to waste their time debugging linear algebra code that you're not even sure if it's incorrect.
Last edited on
jonnin (11437)
your code is huge in part because things like do lines intersect spends 3/4 of its body copying variables just so you can rename them. The compiler cleans it up but its pointless code bloat. Similar instead of saying return (bool expression) you bother with if(expression) return true. if the expression is true, and the result you want is boolean... just return the expression itself, turn 4 lines into 1. I don't usually like to nitpick but every line of code you write is a line you have to maintain, debug, step across, etc. A little extra code for readability is correct; compact gibberish is not good. But tons of extra code for 'my first CS class 101' level of of comprehension is equally bad.

That said, what exactly is the question? Are you having a problem understanding the math, or the code to express the math? It sounds a lot like a math problem (which is fine, but we need to know where the hangup is: you can't express in code what you don't understand!).
Cplusc (457)
@helios, I want to clarify that I didn't mention the code not working. What I meant was that when I consider the intersection of two lines in a plane instead of a 3D space, the results remain consistent with the 3D scenario. Additionally, the concept of complex geometry is self-evident.

About sharing a VERTEX in more than 3 triangles you are right. the vertex can be shared by more than 3 triangles.

My main concern is not related to debugging code but rather the correctness of the logic I applied for line intersection and pointInTriangle. I value your input, and if you feel that this discussion is a waste of time, please feel free to move on instead of lecturing me.
Last edited on
Cplusc (457)
@joinin Thanks for your answer. The problem is not programming. I am just not sure about the logic I used for line intersection and isPointInsideTriangle function. that's all my problem.
Last edited on
helios (17574)
A'ight.
Topic archived. No new replies allowed.