Search code examples
3dgeometrycomputational-geometryapproximationtriangle

3d triangle approximation with rectangular prisms


I am trying to approximate 3d triangles with rectangular prisms but I cant figure out the math

Im using lua

Ive tried flattening the triangle to a 2d plane and approximating the triangles with normal rectangles and then applying rotation again but couldnt figure it out

Ive also tried voxelization but I am trying to use as little rectangular prisms as I can

Rectangle inside triangle

I have been trying to find a big rectangle inside the triangle by drawing a line from the center of the two smallest lines as seen in the picture ^ but I am struggling how to figure out how to find the two points on the bigger line (keep in mind this is 3d)


Solution

  • you got 2 options either use aligned quads (no holes) or arbitrary quads (holes but better coverage)...

    The idea is simply take biggest or alignment side of triangle and fit biggest area quad inscribed in it. here simple C++ example for both:

    //---------------------------------------------------------------------------
    #include <math.h>
    #include "list.h"
    #include "approx.h"
    #include "GLSL_math.h"
    //---------------------------------------------------------------------------
    vec3 tri[3];
    List<vec3> quad;
    //---------------------------------------------------------------------------
    void triangle2quads(vec3 P0,vec3 P1,vec3 P2)
        {
        vec3 p0,p1,p2,p3;
        static vec3 u,v,dp,dq;
        static float u0,u1,v1,w,h;
        // find biggest side (p0,p1)
        // and its size w
        // and its unit direction u
        // p2 is remaining point
        v=P1-P0; u0=length(v);            p0=P0; p1=P1; p2=P2; u=v; w=u0;
        v=P2-P1; u0=length(v); if (w<u0){ p0=P1; p1=P2; p2=P0; u=v; w=u0; }
        v=P0-P2; u0=length(v); if (w<u0){ p0=P2; p1=P0; p2=P1; u=v; w=u0; }
        u/=w;
        // v unit vector perpendicular to u and pointing towards P2
        // h is height of triangle
        v=p2-p0;
        v-=u*dot(u,v);
        h=length(v);
        v/=h;
        // min triangle area end condition
        if (w*h<400.0) return;
        // find max area
        dp=p2-p0; dp/=length(dp);
        dq=p2-p1; dq/=length(dq);
    
        static approx V1;
        static double e;
        for (V1.init(0.05*h,0.95*h,0.1*h,6,&e);!V1.done;V1.step())
            {
            v1=V1.a;
    
            u0=dot(u,dp);               // cos(a)
            u0=v1*u0/sqrt(1.0-u0*u0);   // v1*ctg(a)
    
            u1=dot(-u,dq);              // cos(a)
            u1=v1*u1/sqrt(1.0-u1*u1);   // v1*ctg(a)
            u1=w-u1;
    
            e=(u1-u0)*v1;               // quad area
            e=1.0/(e+0.001);            // error function (is bigger on smaller area)
            }
        v1=V1.aa;
    
        u0=dot(u,dp);               // cos(a)
        u0=v1*u0/sqrt(1.0-u0*u0);   // v1*ctg(a)
    
        u1=dot(-u,dq);              // cos(a)
        u1=v1*u1/sqrt(1.0-u1*u1);   // v1*ctg(a)
        u1=w-u1;
    
        P0=p0; P1=p1; P2=p2;
        p3=p0+(u0*u)+(v1*v);
        p2=p0+(u1*u)+(v1*v);
        p1=p0+(u1*u);
        p0=p0+(u0*u);
    
        quad.add(p0);
        quad.add(p1);
        quad.add(p2);
        quad.add(p3);
    
        triangle2quads(P0,p0,p3);
        triangle2quads(P1,p1,p2);
        triangle2quads(P2,p2,p3);
        }
    //---------------------------------------------------------------------------
    void _triangle2quads_aligned(vec3 P2,vec3 P0,vec3 P1)   // align quad to P0,P1
        {
        vec3 p0,p1,p2,p3;
        static vec3 u,v,dp,dq;
        static float u0,u1,v1,w,h;
        // u is P0,P1
        // and its size w
        v=P1-P0; u0=length(v); p0=P0; p1=P1; p2=P2; u=v; w=u0;
        u/=w;
        // v unit vector perpendicular to u and pointing towards P2
        // h is height of triangle
        v=p2-p0;
        v-=u*dot(u,v);
        h=length(v);
        v/=h;
        // min triangle area end condition
        if (w*h<400.0) return;
        // find max area
        dp=p2-p0; dp/=length(dp);
        dq=p2-p1; dq/=length(dq);
    
        static approx V1;
        static double e;
        for (V1.init(0.05*h,0.95*h,0.1*h,6,&e);!V1.done;V1.step())
            {
            v1=V1.a;
    
            u0=dot(u,dp);               // cos(a)
            u0=v1*u0/sqrt(1.0-u0*u0);   // v1*ctg(a)
    
            u1=dot(-u,dq);              // cos(a)
            u1=v1*u1/sqrt(1.0-u1*u1);   // v1*ctg(a)
            u1=w-u1;
    
            e=(u1-u0)*v1;               // quad area
            e=1.0/(e+0.001);            // error function (is bigger on smaller area)
            }
        v1=V1.aa;
    
        u0=dot(u,dp);               // cos(a)
        u0=v1*u0/sqrt(1.0-u0*u0);   // v1*ctg(a)
    
        u1=dot(-u,dq);              // cos(a)
        u1=v1*u1/sqrt(1.0-u1*u1);   // v1*ctg(a)
        u1=w-u1;
    
        P0=p0; P1=p1; P2=p2;
        p3=p0+(u0*u)+(v1*v);
        p2=p0+(u1*u)+(v1*v);
        p1=p0+(u1*u);
        p0=p0+(u0*u);
    
        quad.add(p0);
        quad.add(p1);
        quad.add(p2);
        quad.add(p3);
    
        _triangle2quads_aligned(P0,p0,p3);
        _triangle2quads_aligned(P1,p1,p2);
        _triangle2quads_aligned(P2,p2,p3);
        }
    //---------------------------------------------------------------------------
    void triangle2quads_aligned(vec3 P0,vec3 P1,vec3 P2)
        {
        int i;
        float l,ll;
        ll=length(P1-P0);            i=0; l=ll;
        ll=length(P2-P1); if (l<ll){ i=1; l=ll; }
        ll=length(P0-P2); if (l<ll){ i=2; l=ll; }
        if (i==0) _triangle2quads_aligned(P2,P0,P1);
        if (i==1) _triangle2quads_aligned(P0,P1,P2);
        if (i==2) _triangle2quads_aligned(P1,P2,P0);
        }
    //---------------------------------------------------------------------------
    

    preview for triangle2quads_aligned(tri[0],tri[1],tri[2]);:

    aligned

    and preview for triangle2quads(tri[0],tri[1],tri[2]);:

    unaligned

    The number in Caption of window is number of quads...

    The quad[] is output list of points 4 per quad. The code is not heavily optimized yet and the recursion can be improved a lot...

    I used mine libs so few explanations...

    I used mine dynamic list template list.h so:


    List<double> xxx; is the same as double xxx[];
    xxx.add(5); adds 5 to end of the list
    xxx[7] access array element (safe)
    xxx.dat[7] access array element (unsafe but fast direct access)
    xxx.num is the actual used size of the array
    xxx.reset() clears the array and set xxx.num=0
    xxx.allocate(100) preallocate space for 100 items

    You can use anything you have at disposal instead like vector<> or whatever (even simple array)

    for vector math I used GLSL like math GLSL_math.h which can be found in here but you can use any other like GLM instead or encode it yourself you need just basic operations like +,*,dot.

    for the search I used mine approximation search approx.h which can be found (and its described) in here:

    You could use linear search or bisection instead (binary search is not usable).

    Now back to the problem if you look into both implementations you see they almost identical first find u,v directions and then search v1 and compute u0,u1 parameters from it and remember biggest area quad.

    After that just add the quad to output list and recursively repeat for 3 remainder/leftover triengles...

    The only difference is that aligned version search u only for first triangle and from then always align to the first quad (last two points passed to _triangle2quads_aligned).

    I avoided goniometrics and use vector math instead I used identities:

     cos(a) = dot(u,v)
     sin(a) = sqrt( 1 - cos^2(a) )
     ctg(a) = cos(a)/sin(a) = cos(a)/sqrt( 1 - cos^2(a) )
    

    where a is angle between u,v and u,v are unit vectors.

    I used coordinates in pixels so in case you have different scale you have to change the ending thresholds to value usable to your case. In code lines:

    if (w*h<400.0) return;
    

    means triangle area must be bigger or equal to 400/2 pixels^2 so just change the 400 to whatever suits you.

    Now the last thing I used QUADs as output if you want boxes (prisms) instead just shift the 4 points perpendicularly by box thickness to get the 8 points of your box ... the shift direction is w=cross(u,v)

    If you want the triangle in the middle of thickness then for each guad p0,p1,p2,p3 output box:

    sh = 0.5*thickness*cross(u,v)
    q0 = p0 - sh
    q1 = p1 - sh
    q2 = p2 - sh
    q3 = p3 - sh
    q4 = p0 + sh
    q5 = p1 + sh
    q6 = p2 + sh
    q7 = p3 + sh