[ggt-devel] Plane vs. Plane intersection test -> Intersection.h

[Benjamin S. <bes...@go...>]

hello,

it works, but is probably not the fastest implementation, but
the fastest way to do it general.
would like to see this in the trunk

--------------------------------------------------------------------------
   /**
    * Tests if the given Planes intersect with each other.
    *
    * @param plane1    first plane
    * @param plane2    second plane
    * @param ray      if successfull, ray will contain the line of
intersection
    *
    * @return  true if the planes intersect; false otherwise
    */
   template<class DATA_TYPE>
   bool intersect(const Plane<DATA_TYPE>& plane1, const Plane<DATA_TYPE>&
plane2, Ray<DATA_TYPE> & ray)
   {
    const Vec<DATA_TYPE, 3> & n1 = plane1.getNormal();
    const Vec<DATA_TYPE, 3> & n2 = plane2.getNormal();
    gmtlASSERT( isNormalized(n1, GMTL_EPSILON) );
    gmtlASSERT( isNormalized(n2, GMTL_EPSILON) );

    Vec<DATA_TYPE, 3> u;
    cross(u, n1, n2);

    float    ax = (u[0] >= 0 ? u[0] : -u[0]);
    float    ay = (u[1] >= 0 ? u[1] : -u[1]);
    float    az = (u[2] >= 0 ? u[2] : -u[2]);

    // test if the two planes are parallel
    if ((ax+ay+az) < GMTL_EPSILON)
    {       // Pn1 and Pn2 are near parallel
        return false;        //could check why it failes -> in which
relations the planes lie.
    }

    // first determine max abs coordinate of cross product
    // choosing the max coordinate gives best computational stability
    int      maxc;        // max coordinate
    if (ax > ay)
    {
        if (ax > az)
            maxc = 0;
        else maxc = 2;
    }
    else
    {
        if (ay > az)
            maxc = 1;
        else maxc = 2;
    }

    // next, to get a point on the intersect line
    // zero the max coord, and solve for the other two
    Point3f p;                    // intersect point
    const float & d1 = plane1.getOffset();
    const float & d2 = plane2.getOffset();

    switch (maxc)                // select max coordinate
    {
        case 0:                // intersect with x=0
            p[0] = 0;
            p[1] = (d2*n1[2] - d1*n2[2]) / u[0];
            p[2] = (d1*n2[1] - d2*n1[1]) / u[0];
            break;
        case 1:                // intersect with y=0
            p[0] = (d1*n2[2] - d2*n1[2]) / u[1];
            p[1] = 0;
            p[2] = (d2*n1[0] - d1*n2[0]) / u[1];
            break;
        case 2:                // intersect with z=0
            p[0] = (d2*n1[1] - d1*n2[1]) / u[2];
            p[1] = (d1*n2[0] - d2*n1[0]) / u[2];
            p[2] = 0;
    }

    ray.setOrigin(p);
    ray.setDir(u);

    return true;
   }
------------------------------------------