Finding Polygon of Plane-AABB Intersection

Warning! Some information on this page is older than 5 years now. I keep it for reference, but it probably doesn't reflect my current knowledge and beliefs.

# Finding Polygon of Plane-AABB Intersection

Tue
15
Feb 2011

While doing volumetric rendering according to the article Volume Rendering Techniques (Milan Ikits, Joe Kniss, Aaron Lefohn, Charles Hansen, Chapter 39 of "GPU Gems" book), I came across a problem of calculating intersection between a plane and an AABB (axis-aligned bounding box), which forms a 3- to 6-vertex polygon. The solution is actually described in the article, but here I'd like to share my implementation (using D3DX).

The solution is based on finding points of intersection between the plane and every edge of the box. So we first need a ray-to-plane intersection code. My function for this is based on: Ray-Plane Intersection.

// OutVD > 0 means ray is back-facing the plane
// returns false if there is no intersection because ray is perpedicular to plane
bool ray_to_plane(const D3DXVECTOR3 &RayOrig, const D3DXVECTOR3 &RayDir, const D3DXPLANE &Plane, float *OutT, float *OutVD)
{
    *OutVD = Plane.a * RayDir.x + Plane.b * RayDir.y + Plane.c * RayDir.z;
    if (*OutVD == 0.0f)
        return false;
    *OutT = - (Plane.a * RayOrig.x + Plane.b * RayOrig.y + Plane.c * RayOrig.z + Plane.d) / *OutVD;
    return true;
}

Now having a segment A-B between two points, we can calculate ray_to_plane intersection where RayOrig=A, RayDir=B-A and we get parameter T. The segment is a subset of the infinite ray. The intersection point lies on this segment when 0 <= T <= 1. The position of this point is A+T*B. Now it's time for plane-box intersection code:

// Maximum out_point_count == 6, so out_points must point to 6-element array.
// out_point_count == 0 mean no intersection.
// out_points are not sorted.
void calc_plane_aabb_intersection_points(const D3DXPLANE &plane,
    const D3DXVECTOR3 &aabb_min, const D3DXVECTOR3 &aabb_max,
    D3DXVECTOR3 *out_points, unsigned &out_point_count)
{
    out_point_count = 0;
    float vd, t;

    // Test edges along X axis, pointing right.
    D3DXVECTOR3 dir = D3DXVECTOR3(aabb_max.x - aabb_min.x, 0.f, 0.f);
    D3DXVECTOR3 orig = aabb_min;
    if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
        out_points[out_point_count++] = orig + dir * t;
    orig = D3DXVECTOR3(aabb_min.x, aabb_max.y, aabb_min.z);
    if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
        out_points[out_point_count++] = orig + dir * t;
    orig = D3DXVECTOR3(aabb_min.x, aabb_min.y, aabb_max.z);
    if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
        out_points[out_point_count++] = orig + dir * t;
    orig = D3DXVECTOR3(aabb_min.x, aabb_max.y, aabb_max.z);
    if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
        out_points[out_point_count++] = orig + dir * t;

    // Test edges along Y axis, pointing up.
    dir = D3DXVECTOR3(0.f, aabb_max.y - aabb_min.y, 0.f);
    orig = D3DXVECTOR3(aabb_min.x, aabb_min.y, aabb_min.z);
    if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
        out_points[out_point_count++] = orig + dir * t;
    orig = D3DXVECTOR3(aabb_max.x, aabb_min.y, aabb_min.z);
    if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
        out_points[out_point_count++] = orig + dir * t;
    orig = D3DXVECTOR3(aabb_min.x, aabb_min.y, aabb_max.z);
    if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
        out_points[out_point_count++] = orig + dir * t;
    orig = D3DXVECTOR3(aabb_max.x, aabb_min.y, aabb_max.z);
    if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
        out_points[out_point_count++] = orig + dir * t;

    // Test edges along Z axis, pointing forward.
    dir = D3DXVECTOR3(0.f, 0.f, aabb_max.z - aabb_min.z);
    orig = D3DXVECTOR3(aabb_min.x, aabb_min.y, aabb_min.z);
    if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
        out_points[out_point_count++] = orig + dir * t;
    orig = D3DXVECTOR3(aabb_max.x, aabb_min.y, aabb_min.z);
    if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
        out_points[out_point_count++] = orig + dir * t;
    orig = D3DXVECTOR3(aabb_min.x, aabb_max.y, aabb_min.z);
    if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
        out_points[out_point_count++] = orig + dir * t;
    orig = D3DXVECTOR3(aabb_max.x, aabb_max.y, aabb_min.z);
    if (ray_to_plane(orig, dir, plane, &t, &vd) && t >= 0.f && t <= 1.f)
        out_points[out_point_count++] = orig + dir * t;
}

But that's not all. Returned points are not sorted! If we don't sort them, we may in some cases end up with a mess like this:

Here comes another clever algorithm - the one for sorting vertices (counter?)clockwise. I came up with it on my own and it's a bit similar to the algorithm for finding a convex hull. We are given a set of points that all lie in a single plane and already form a convex polygon. We want to sort them clockwise "around some center point". Beginner could probably think here about calculating such center point, then for each vertex calculating and comparing some angle using atan2, of course in degrees ;)

But that's not the solution. What we need here is a way to tell whether one vector is more "to the right" then the other, relative to some reference point. Whether one vector lies on the right-hand side of the other can be easily told in 2D using the sign of a scalar returned by 2D cross-product: a.x*b.y - b.x*a.y. In 3D it's not that easy. My idea is to use 3D cross-product. It returns 3D vector, but has similar property - it is anticommutative, which means: cross(a,b) == -cross(b,a). Now if only I could tell whether the direction of the resulting vector is "right" or "opposite"... That's what the sign of dot product means! I'm just going to use dot product between the cross product result and the normal vector of the plane. Eureka! Now I have all necessary parts. I choose first vertex as a reference. I also use lambda expression from C++0x (supported in Visual C++ 2010) to more conveniently pass custom comparator to std::sort STL algorithm (otherwise I would have to write a functor). Here is the final solution:

#include <algorithm>

void sort_points(D3DXVECTOR3 *points, unsigned point_count, const D3DXPLANE &plane)
{
    if (point_count == 0) return;

    const D3DXVECTOR3 plane_normal = D3DXVECTOR3(plane.a, plane.b, plane.c);
    const D3DXVECTOR3 origin = points[0];

    std::sort(points, points + point_count, [&](const D3DXVECTOR3 &lhs, const D3DXVECTOR3 &rhs) -> bool {
        D3DXVECTOR3 v;
        D3DXVec3Cross(&v, &(lhs - origin), &(rhs - origin));
        return D3DXVec3Dot(&v, &plane_normal) < 0;
    } );
}

Comments | #directx #math Share

Comments

STAT NO AD
[Stat] [STAT NO AD] [Download] [Dropbox] [pub] [Mirror] [Privacy policy]
Copyright © 2004-2019