//server/stuff/Convex hull solver>
Main

Generating a convex hull from a set of planes

This piece of code shows how to extract the points making up a convex hull defined by a set of planes, using their plane equations for input, in the form ax + by + cz + d = 0. The code will test for an intersection point between every three planes of the input set, using Gauss-Jordan matrix solving by row elimination. If a single point can not be found, the Intersect function prints out the first ambigious axis, so if needed, you may print out an intersection "line" between such planes, keeping in mind that three planes may have two or three seperate lines as intersections.

If you wish to triangulate this convex hull, you will need to determine what subset of the set of points found lie on each plane, then triangulate each subset as you would for a convex polygon.
#include <iostream>

using namespace std;

#define f_eq(f1, f2) ((f1 - f2) < 0.0001f && (f2 - f1) < 0.0001f)
#define pf(x) ((!f_eq(x, 0.0f)) ? x : 0.0f)

struct Plane
{
    float a, b, c, d;
    Plane(float a, float b, float c, float d) : a(a), b(b), c(c), d(d)
    { }
    Plane() : a(0), b(0), c(0), d(0)
    { }
    
    void Scale(float s)
    { a *= s; b *= s; c *= s; d *= s; }
    
    void Add(const Plane& o, float s)
    { a += o.a * s; b += o.b * s; c += o.c * s; d += o.d * s; }
};

ostream& operator<<(ostream& o, const Plane& p)
{
    o << pf(p.a) << "x + " << pf(p.b) << "y + " << pf(p.c) << "z + " << pf(p.d) << " = 0";
    
    return o;
}

void IntersectPlanes(Plane x, Plane y, Plane z)
{
#define SwapPlane(p1, p2) Plane t = p1; p1 = p2; p2 = t;
#define Swap0(p1, p2, v) if (!f_eq(v, 0.0f)) { SwapPlane(p1, p2) }
#define ScaleOp(p, s) p.Scale(1.0f / p.s);
#define AddOp(p1, p2, s) if (!f_eq(p1.s, 0.0f)) { p1.Add(p2, -p1.s / p2.s); }
#define RemOp(p1, p2, s) if (!f_eq(p2.s, 0.0f)) { p2.Add(p1, -p2.s); }
    
    if (f_eq(x.a, 0.0f))
    {
        Swap0(x, y, y.a)
        else Swap0(x, z, z.a)
    }
    if (f_eq(y.b, 0.0f))
    {
        Swap0(y, x, x.b)
        else Swap0(y, z, z.b)
    }
    if (f_eq(x.a, 0.0f)) { cout << "Can't find x" << endl; return; }
    if (f_eq(y.b, 0.0f)) { cout << "Can't find y" << endl; return; }
    if (f_eq(z.c, 0.0f)) { cout << "Can't find z" << endl; return; }
    
    AddOp(y, x, a)
    AddOp(z, x, a)
    AddOp(z, y, b)
    
    ScaleOp(x, a)
    ScaleOp(y, b)
    ScaleOp(z, c)
    
    RemOp(y, x, b);
    RemOp(z, x, c);
    RemOp(z, y, c);
    
    cout << "X = " << -x.d << ", Y = " << -y.d << ", Z = " << -z.d << endl;

#undef SwapPlane
#undef Swap0
#undef ScaleOp
#undef AddOp
#undef RemOp
}

int main()
{
    Plane planes[256];
    int NumPlanes;
    cout << "Number of planes: ";
    cin >> NumPlanes;
    if (NumPlanes < 3)
    {
        cout << "Minimium plane count is 3" << endl;
        return 0;
    }
    for (int i = 0; i < NumPlanes; i++)
    {
        cout << "Plane " << i + 1 << ": ax + by + cz + d = 0" << endl << "Enter a, b, c and d" << endl;
        cin >> planes[i].a >> planes[i].b >> planes[i].c >> planes[i].d;
        cout << planes[i] << endl << endl;
    }
    for (int i = 0; i < NumPlanes; i++)
    {
        for (int j = i + 1; j < NumPlanes; j++)
        {
            for (int k = j + 1; k < NumPlanes; k++)
            {
                cout << i + 1 << ", " << j + 1 << ", " << k + 1 << ": ";
                IntersectPlanes(planes[i], planes[j], planes[k]);
            }
        }
    }
    
    return 0;
}