Kód Průnik paprsku s trojúhelníkem v CPP

Z CHWiki

autori:

l-BigBoss and Tutchek
ladik-BigBoss(at)seznam.cz
michal.tulacek(at)igroup.cz

zdroje:

http://www.cs.princeton.edu/courses/archive/fall00/cs426/lectures/raycast/sld017.htm
http://www.cs.princeton.edu/courses/archive/fall00/cs426/lectures/raycast/sld018.htm
http://mathworld.wolfram.com/

ukazkovy program - pouzito pro lightmapping:

http://www.hrebcin-dobrovsky.com/bigboss/lightmapping.html

kód:

class Vector
{
    double x;
    double y;
    double z;
    // a potrebne pretizeni operatoru, metody... napiste si sami :)
}

double dot(const Vector& u, const Vector& v)	// skalarni soucin 3D vektoru
{
    return (u.x*v.x + u.y*v.y + u.z*v.z);
}

Vector cross(const Vector & u, const Vector& v)	// vektorovy soucin 3D vektoru
{ 
    return Vector(u.y*v.z - u.z*v.y, u.z*v.x - u.x*v.z, u.x*v.y - u.y*v.x);
}

inline int sgn(double x)	// zjisteni znamenka cisla, sgn(0) = 0
{
     return (x == 0 ? 0 : (x > 0 ? 1 : -1 ) );
}

// A,B,C - vrcholy trojuhelniku; P0 - pocatek secky, V - vektor ke konci usecky, P - prusecik
bool intersect(Vector A, Vector B, Vector C, Vector P0, Vector V, Vector &P)	
{
// PRUSECIK S ROVINOU
    // ABC - trojuhelnik zadavajici rovinu
    // u=B-A, v=C-A, N=u×v - N je normala roviny
    Vector N=cross(B-A, C-A);
    // ax+by+cz+d=0	// x,y,z libovolny bod v rovine <=> A
    // N=(a,b,c)
    double d=-(N.x*A.x + N.y*A.y + N.z*A.z);
    double t, dot_vn;

    if ((dot_vn=dot(V, N)) != 0)		// zamezeni /0
    {
        t=-(dot(P0, N)+d)/dot_vn;
        if(t>0.0 && t<1.0)		// lezi na usecce?
            return false;
    }
    else
        return false;

    P=P0+V*t; // vypocet pruseciku

//TEST JESTLI JE PRUSECIK V TROJUHELNIKU ABC
  // P0 ustavime vrcholem jehlanu P0,A,B,C

    Vector V1=A-P0;
    Vector V2=B-P0;
    Vector V3=C-P0;

    P0=(P0*(-1.0f)); // vynasobeni -1; jinak d=dot(-P0, Nx)

    // test jestli je P ve stejnem poloprostoru delenem rovinou steny A, B, P0 jako C
    Vector Nx=cross(V1, V2);
    d=dot(P0, Nx);	// -P0
    if(sgn(dot(P, Nx)+d) != sgn(dot(C, Nx)+d))   
        return false;
	
    // test jestli je P ve stejnem poloprostoru delenem rovinou steny A, C, P0 jako B
    Nx=cross(V1, V3);
    d=dot(P0, Nx);
    if(sgn(dot(P, Nx)+d) != sgn(dot(B, Nx)+d))   
        return false;

    // test jestli je P ve stejnem poloprostoru delenem rovinou steny B, C, P0 jako A
    Nx=cross(V2, V3);
    d=dot(P0, Nx);
    if(sgn(dot(P, Nx)+d) != sgn(dot(A, Nx)+d))   
        return false;

    return true; // trojuhelnik je protinan v bode vracenem na referenci P
}