biVector forum

Frustum Culling in GA?

Hi. I am new to Geometric Algebra. I have watched a few of the GAME2020 talks, so I know a bit of the terminology. But I don’t have any hands-on experience yet.

I would like to know if it is possible to define a non-infinite volume and check if this volume intersects with another non-infinite volume.

In particular, I would like to check if an axis-aligned bounding box intersects with an OpenGL view frustum. For example, here’s an implementation that uses the separating axis theorem, but that involves quite a bit of work with lots of opportunities for small mistakes. Instead, would it be possible to define this intersection in GA with a meet operation or similar?

Looking at vectors, bivectors and trivectors: Those define oriented lines, areas, volumes. But if I understand it correctly, then those are defined with respect to the origin? I wouldn’t know how to solve my problem this way.

Looking at the 3D PGA cheatsheet: This shows how to define lines and planes (and I assume volumes in 4D), but those are of infinite size? This also doesn’t help me.

I would appreciate some pointers into the right direction.

Thank you
Martin

@Martin
This is an interesting question.
I think that PGA can be of use here.
First, obtain a representation of the planes p_i of the frustum satisfy p_i^2 = 1 (normalized) and p_i \vee O < 0 where O is the “mid-point” of the frustum (intersection of the space diagonals of the frustum, any interior point of the frustum will serve the same purpose). Hence a normalized point P satisfying p_i \vee P > 0 for any i lies outside the frustum. (Note: using the regressive operator \vee here produces a scalar directly, thereby avoiding having to extract the magnitude of the pseudoscalar p_i \wedge P. This trick is based on the identity (p \vee P)\mathbf{I} = p \wedge P for a plane p and a point P.)

I think the answer to your question can be obtained by a sequence of two tests; if either yields a positive result there is an intersection of volumes; if both yield a negative result, there is no intersection.

The first test is based on the observation that two convex polyedra A and B intersect whenever a vertex of A lies in the interior of B, or vice-versa. This can be easily tested using in PGA. Let A be the frustum and B the object, and further let \{P_j\} be the set of normalized vertices of B. Let d_{ij} := p_i \vee P_j be the signed distances of p_i from P_j. Then if d_{ij} < 0 for all i and some j, the vertex lies in the frustum and hence the object intersects the frustrum. If any vertex satisfies this condition, the object intersects the frustum and we are done.

Otherwise proceed to a second and final test.
Here the idea is that two convex polyedra A and B intersect whenever 1) the vertices of an edge e of A lie on opposite sides of a face plane of B, and 2) the line determined by e intersects the interior of this face.
First, pre-compute 24 edge lines of the frustum f_{mn} for m\in \{1,2,3,4,5,6\} and n\in \{1,2,3,4\}, so that \{f_{m1},f_{m2},f_{m3},f_{m4}\} are the lines obtained by a cyclic traversal of the vertices of the m^{th} quadrilateral face of the frustum. That is, if \{F_{m_1}, F_{m_2}, f_{m_3}, f_{m_4}\} are the vertices of the m^{th} quadrilateral in cyclic order, then f_{m1} := F_{m_1} \vee F_{m_2}, …, f_{m4} := F_{m_4} \vee F_{m_1}. Observation (Exercise) : a line x intersects the interior of the face f_m exactly when all the scalars f_{mi} \vee x have the same (non-zero) sign.

To carry out this second test, run through the edges e_{jk} of B. The endpoints of e_{jk} are P_j and P_k. If d_{ij} * d_{ik} >= 0, the two vertices do not lie on opposite sides of the face plane p_i, and we proceed to the next edge. If however d_{ij} * d_{ik} < 0, then use the above observation to check whether the edge line x_{jk} := P_j \vee P_k intersects the face f_i. If it does, then the object intersects the frustum. If no edge satisfies this test, then the object does not intersect the frustum.

To avoid numerical problems on the boundary I recommend using a slightly extended frustum obtained by scaling the frustum around the point O. This can be implemented by replacing the frustum plane p_i = d_i e_0 + a_i e_1 + b_i e_2 + c_i e_3 with (1+\epsilon)d_i e_0 + a_i e_1 + b_i e_2 + c_i e_3. This illustrates the general principle that adding k e_0 to a normalized plane translates the plane a (signed) distance k in its normal direction.

Thank you very much for your answer.

I was thinking about whether the first test (vertex inside frustum) is even necessary? Does the second test not catch this, or are there exceptions I have not considered?

I believe you are right, the second test suffices.
Note however that both tests depend on the signed distances d_{ij} from the vertices of the object to the planes of the frustum; once you have them test 1 is trivial to perform.
Furthermore, I can’t prove it but it seems to me that most intersections will be caught by test 1.

I just found a case, where the 1st test is necessary (at least for one vertex): If the AABB is completely inside the frustum, then the 2nd test doesn’t hit.

Thanks again. Your algorithm is really helpful.

One more thing, if you plan to implement the above algorithm:

It uses the neat identity (K \vee M)\mathbf{I} = K \wedge M, for a k-vector K and m-vector M satisfying k+m=n where n is the dimension of the underlying vector space. As I mentioned, using the \vee produces a signed scalar magnitude directly without having to extract it from the pseudoscalar produced by the \wedge. It is traditional to extract \alpha, the signed magnitude of a pseudoscalar p = \alpha\mathbf{I}, via p\mathbf{I}^{-1}, but for a degenerate metric this is not defined.

For P(\mathbb{R}^*_{3,0,1}), the algebra of interest here, n = 4. This gives rise to the two cases k=1,m=3 (plane and point) and k=2,m=2 (two lines), both used in the algorithm.

While this formula is true for the implementation described on the SIGGRAPH cheat sheets, the current implementation of ganja.js introduces a minus sign for k=1,m=3: (K \vee M)\mathbf{I} = -K \wedge M. (As far as I can tell, this is due to ganja reversing the order of the arguments to the regressive product \vee vis-a-vis the form given on the cheat sheet.)

I’ve written a short test case in ganja illustrating this behavior that can be found here.

BTW, if I’m not mistaken, all the algorithms presented above also work with infinite convex volumes. One only has to make sure that the infinite face(s) is (are) represented as a polygon (polygons) in the ideal plane and that the vertices occur in cyclic order as prescribed when computing the edges.

Note that when the ideal plane is a boundary plane of the polyhedron, then every euclidean vertex will lie on the same side of this plane. In particular any vertex lies on the same side as the midpoint O does, so every vertex will lie with O on the “inside” of the ideal plane. Which is consistent with our experience of it: it encloses euclidean space.

I think the linked webpage “SIGGRAPH cheat sheets” is missing in your answer. I tried clicking on it, but it doesn’t open.

I used the code generator from BiVector.net to get a PGA3D c++ implementation. I created an OpenGL desktop application to test out the culling, but in the end, I will need to port this code to run on the iPhone. I will probably use the klein library or port the algorithm to use apple’s simd library.

Here is the source code to my desktop test application. I am running it on macOS 11.5.2, using glfw and glm from homebrew. The main file might have to be modified a bit on other platforms to get OpenGL running (I’m using version 2.1)

Note that the code is not optimized at all. I just wanted to get the algorithm implemented (hopefully) correctly.

frustum_culling.cpp

#include <GLFW/glfw3.h>

#include <glm/glm.hpp>
#include <glm/ext.hpp>
#include <glm/gtx/quaternion.hpp>

#include <cmath>
#include <vector>
#include <iostream>

struct AABB
{
    glm::vec3 min;
    glm::vec3 max;
};

struct Viewport
{
    float x, y, width, height;
};

// https://discourse.bivector.net/t/frustum-culling-in-ga/430
std::vector<bool> intersectFrustum(const glm::mat4x4& view_projection, const std::vector<AABB>& aabbs)
{
    // http://donw.io/post/frustum-point-extraction
    PGA3D frustum_points[8];
    {
        glm::mat4x4 clip_to_world = glm::inverse(view_projection);

        // near == [0 1 2 3]
        // far  == [4 5 6 7]
        // origin at [0 0 0]
        //
        //   7----6      y
        //  /|   /|      ^
        // 3----2 |      |
        // | 4--|-5      o--> x
        // |/   |/      /
        // 0----1      z
        glm::vec4 frustum[8] =
        {
            {-1, -1, -1, 1}, // 0
            { 1, -1, -1, 1}, // 1
            { 1,  1, -1, 1}, // 2
            {-1,  1, -1, 1}, // 3
            {-1, -1,  1, 1}, // 4
            { 1, -1,  1, 1}, // 5
            { 1,  1,  1, 1}, // 6
            {-1,  1,  1, 1}  // 7
        };

        for (int i = 0; i < 8; i++)
        {
            frustum[i] = clip_to_world * frustum[i];
            frustum[i] /= frustum[i].w;

            frustum_points[i] = point(frustum[i].x, frustum[i].y, frustum[i].z);
        }
    }

    // I constructed the planes such that frustum_planes[i] v mid_point < 0
    //
    // Another option to define the planes, but I need the frustum_points anyways
    // https://fgiesen.wordpress.com/2012/08/31/frustum-planes-from-the-projection-matrix
    PGA3D frustum_planes[6];
    {
        frustum_planes[0] = (frustum_points[0] & frustum_points[3] & frustum_points[7]).normalized(); // left
        frustum_planes[1] = (frustum_points[1] & frustum_points[5] & frustum_points[6]).normalized(); // right
        frustum_planes[2] = (frustum_points[0] & frustum_points[4] & frustum_points[5]).normalized(); // bottom
        frustum_planes[3] = (frustum_points[2] & frustum_points[6] & frustum_points[7]).normalized(); // top
        frustum_planes[4] = (frustum_points[0] & frustum_points[1] & frustum_points[2]).normalized(); // near
        frustum_planes[5] = (frustum_points[4] & frustum_points[7] & frustum_points[6]).normalized(); // far

        // TODO: This should be done both for the planes and lines
        //// Extend frustum around mid_point to avoid numerical problems on the boundary
        //for (int i = 0; i < 6; ++i)
        //{
        //    frustum_planes[i][A0] += 0.01;
        //}
    }

    // Edge lines of the frustum
    PGA3D f_mn[6][4];
    {
        // left
        f_mn[0][0] = frustum_points[0] & frustum_points[3];
        f_mn[0][1] = frustum_points[3] & frustum_points[7];
        f_mn[0][2] = frustum_points[7] & frustum_points[4];
        f_mn[0][3] = frustum_points[4] & frustum_points[0];

        // right
        f_mn[1][0] = frustum_points[1] & frustum_points[5];
        f_mn[1][1] = frustum_points[5] & frustum_points[6];
        f_mn[1][2] = frustum_points[6] & frustum_points[2];
        f_mn[1][3] = frustum_points[2] & frustum_points[1];

        // bottom
        f_mn[2][0] = frustum_points[0] & frustum_points[4];
        f_mn[2][1] = frustum_points[4] & frustum_points[5];
        f_mn[2][2] = frustum_points[5] & frustum_points[1];
        f_mn[2][3] = frustum_points[1] & frustum_points[0];

        // top
        f_mn[3][0] = frustum_points[2] & frustum_points[6];
        f_mn[3][1] = frustum_points[6] & frustum_points[7];
        f_mn[3][2] = frustum_points[7] & frustum_points[3];
        f_mn[3][3] = frustum_points[3] & frustum_points[2];

        // near
        f_mn[4][0] = frustum_points[0] & frustum_points[1];
        f_mn[4][1] = frustum_points[1] & frustum_points[2];
        f_mn[4][2] = frustum_points[2] & frustum_points[3];
        f_mn[4][3] = frustum_points[3] & frustum_points[0];

        // far
        f_mn[5][0] = frustum_points[4] & frustum_points[7];
        f_mn[5][1] = frustum_points[7] & frustum_points[6];
        f_mn[5][2] = frustum_points[6] & frustum_points[5];
        f_mn[5][3] = frustum_points[5] & frustum_points[4];
    }

    // Edge indices of the AABB
    glm::ivec2 e_jk[12];
    {
        e_jk[0]  = glm::ivec2(0, 1);
        e_jk[1]  = glm::ivec2(1, 2);
        e_jk[2]  = glm::ivec2(2, 3);
        e_jk[3]  = glm::ivec2(3, 0);
        e_jk[4]  = glm::ivec2(4, 5);
        e_jk[5]  = glm::ivec2(5, 6);
        e_jk[6]  = glm::ivec2(6, 7);
        e_jk[7]  = glm::ivec2(7, 4);
        e_jk[8]  = glm::ivec2(0, 4);
        e_jk[9]  = glm::ivec2(1, 5);
        e_jk[10] = glm::ivec2(2, 6);
        e_jk[11] = glm::ivec2(3, 7);
    }

    float d_ij[6][8];
    PGA3D aabb_points[8];
    const auto intersect = [&frustum_planes, &f_mn, &e_jk, &d_ij, &aabb_points](const AABB& aabb)
    {
        aabb_points[0] = point(aabb.min.x, aabb.min.y, aabb.min.z).normalized();
        aabb_points[1] = point(aabb.max.x, aabb.min.y, aabb.min.z).normalized();
        aabb_points[2] = point(aabb.max.x, aabb.max.y, aabb.min.z).normalized();
        aabb_points[3] = point(aabb.min.x, aabb.max.y, aabb.min.z).normalized();
        aabb_points[4] = point(aabb.min.x, aabb.min.y, aabb.max.z).normalized();
        aabb_points[5] = point(aabb.max.x, aabb.min.y, aabb.max.z).normalized();
        aabb_points[6] = point(aabb.max.x, aabb.max.y, aabb.max.z).normalized();
        aabb_points[7] = point(aabb.min.x, aabb.max.y, aabb.max.z).normalized();

        // Check if any aabb point is inside the frustum
        for (int j = 0; j < 8; ++j)
        {
            const PGA3D& P_j = aabb_points[j];

            d_ij[0][j] = (frustum_planes[0] & P_j)[AO];
            d_ij[1][j] = (frustum_planes[1] & P_j)[AO];
            d_ij[2][j] = (frustum_planes[2] & P_j)[AO];
            d_ij[3][j] = (frustum_planes[3] & P_j)[AO];
            d_ij[4][j] = (frustum_planes[4] & P_j)[AO];
            d_ij[5][j] = (frustum_planes[5] & P_j)[AO];

            if (d_ij[0][j] < 0 && d_ij[1][j] < 0 && d_ij[2][j] < 0 && d_ij[3][j] < 0 && d_ij[4][j] < 0 && d_ij[5][j] < 0)
            {
                return true;
            }
        }

        // Check if any aabb edge intersects any of the frustum faces
        for (int i = 0; i < 6; ++i) // frustum face
        {
            for (const glm::ivec2& e : e_jk) // edges
            {
                if (d_ij[i][e.x] * d_ij[i][e.y] < 0) // edge points on opposite sides of the face
                {
                    PGA3D x_jk = (aabb_points[e.x] & aabb_points[e.y]); // edge line

                    const float e0 = (f_mn[i][0] & x_jk)[AO];
                    const float e1 = (f_mn[i][1] & x_jk)[AO];
                    const float e2 = (f_mn[i][2] & x_jk)[AO];
                    const float e3 = (f_mn[i][3] & x_jk)[AO];

                    if (e0 != 0 && e1 != 0 && e2 != 0 && e3 != 0)
                    {
                        const bool sb = std::signbit(e0);
                        if (std::signbit(e1) == sb && std::signbit(e2) == sb && std::signbit(e3) == sb)
                        {
                            return true;
                        }
                    }
                }
            }
        }

        return false;
    };

    std::vector<bool> intersecting(aabbs.size());
    for (int aabb_id = 0; aabb_id < aabbs.size(); ++aabb_id)
    {
        intersecting[aabb_id] = intersect(aabbs[aabb_id]);
    }
    return intersecting;
}

bool intersectFrustum(const glm::mat4x4& view_projection, const AABB& aabb)
{
    std::vector<AABB> aabbs(1, aabb);
    return intersectFrustum(view_projection, aabbs)[0];
}

// Indices for drawing the lines of the box / frustum
static const uint8_t indices[] = {0, 1, 1, 2, 2, 3, 3, 0,
                                  4, 5, 5, 6, 6, 7, 7, 4,
                                  0, 4, 1, 5, 2, 6, 3, 7};

void drawFrustum(const glm::mat4x4& view_projection, const glm::vec3& color)
{
    // http://donw.io/post/frustum-point-extraction
    glm::mat4x4 clip_to_world = glm::inverse(view_projection);

    glm::vec4 frustum[8] =
    {
        {-1, -1, -1, 1}, // 0
        { 1, -1, -1, 1}, // 1
        { 1,  1, -1, 1}, // 2
        {-1,  1, -1, 1}, // 3
        {-1, -1,  1, 1}, // 4
        { 1, -1,  1, 1}, // 5
        { 1,  1,  1, 1}, // 6
        {-1,  1,  1, 1}  // 7
    };

    for (int i = 0; i < 8; i++)
    {
        frustum[i] = clip_to_world * frustum[i];
        frustum[i] /= frustum[i].w;
    }

    glColor3fv(&color.r);
    glEnableClientState(GL_VERTEX_ARRAY);
    glVertexPointer(4, GL_FLOAT, 0, frustum);
    glDrawElements(GL_LINES, 24, GL_UNSIGNED_BYTE, indices);
    glDisableClientState(GL_VERTEX_ARRAY);
}

void drawAABB(const AABB& box, const glm::vec3& color)
{
    const float v[] = {box.min.x, box.min.y, box.min.z, 1,  // 0
                       box.max.x, box.min.y, box.min.z, 1,  // 1
                       box.max.x, box.max.y, box.min.z, 1,  // 2
                       box.min.x, box.max.y, box.min.z, 1,  // 3
                       box.min.x, box.min.y, box.max.z, 1,  // 4
                       box.max.x, box.min.y, box.max.z, 1,  // 5
                       box.max.x, box.max.y, box.max.z, 1,  // 6
                       box.min.x, box.max.y, box.max.z, 1}; // 7

    glColor3fv(&color.r);
    glEnableClientState(GL_VERTEX_ARRAY);
    glVertexPointer(4, GL_FLOAT, 0, v);
    glDrawElements(GL_LINES, 24, GL_UNSIGNED_BYTE, indices);
    glDisableClientState(GL_VERTEX_ARRAY);
}

AABB aabb;

glm::mat4x4 left_camera_pose;
glm::mat4x4 left_projection;

glm::mat4x4 right_camera_pose;
glm::mat4x4 right_projection;

float left_aspect = 1.2;
bool left_is_current = true;

Viewport viewport;

glm::vec2 cursor_position          = glm::vec2(0, 0);
glm::vec2 previous_cursor_position = glm::vec2(0, 0);

void render(GLFWwindow* window)
{
    const glm::mat4x4 left_view            = glm::inverse(left_camera_pose);
    const glm::mat4x4 left_view_projection = left_projection * left_view;

    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    // Draw line separating the viewports
    {
        glMatrixMode(GL_MODELVIEW);
        glLoadIdentity();
        glMatrixMode(GL_PROJECTION);
        glLoadIdentity();

        glViewport(viewport.x, viewport.y, viewport.width, viewport.height);

        const float x = (left_aspect / (left_aspect + 1)) * 2 - 1;

        glColor3f(1, 1, 1);
        glBegin(GL_LINES);
        glVertex2f(x, -1);
        glVertex2f(x, 1);
        glEnd();
    }

    const glm::vec3 aabb_color = intersectFrustum(left_view_projection, aabb) ? glm::vec3(0, 1, 0) : glm::vec3(1, 0, 0);

    // Draw into left viewport
    {
        Viewport left_viewport;
        left_viewport.x      = viewport.x;
        left_viewport.y      = viewport.y;
        left_viewport.width  = (left_aspect / (left_aspect + 1)) * viewport.width;
        left_viewport.height = viewport.height;
        glViewport(left_viewport.x, left_viewport.y, left_viewport.width, left_viewport.height);

        glMatrixMode(GL_MODELVIEW);
        glLoadMatrixf(&left_view[0][0]);
        glMatrixMode(GL_PROJECTION);
        glLoadMatrixf(&left_projection[0][0]);

        // Draw aabb
        drawAABB(aabb, aabb_color);

        // Draw origin
        glBegin(GL_LINES);
        glColor3f(1, 0, 0); glVertex3f(0, 0, 0); glVertex3f(1, 0, 0);
        glColor3f(0, 1, 0); glVertex3f(0, 0, 0); glVertex3f(0, 1, 0);
        glColor3f(0, 0, 1); glVertex3f(0, 0, 0); glVertex3f(0, 0, 1);
        glEnd();
    }

    // Draw into right viewport
    {
        Viewport right_viewport;
        right_viewport.x      = (left_aspect / (left_aspect + 1)) * viewport.width;
        right_viewport.y      = viewport.y;
        right_viewport.width  = (1 - left_aspect / (left_aspect + 1)) * viewport.width;
        right_viewport.height = viewport.height;
        glViewport(right_viewport.x, right_viewport.y, right_viewport.width, right_viewport.height);

        const glm::mat4x4 model_view = glm::inverse(right_camera_pose);

        glMatrixMode(GL_MODELVIEW);
        glLoadMatrixf(&model_view[0][0]);
        glMatrixMode(GL_PROJECTION);
        glLoadMatrixf(&right_projection[0][0]);

        // Draw aabb
        drawAABB(aabb, aabb_color);

        // Draw origin
        glBegin(GL_LINES);
        glColor3f(1, 0, 0); glVertex3f(0, 0, 0); glVertex3f(1, 0, 0);
        glColor3f(0, 1, 0); glVertex3f(0, 0, 0); glVertex3f(0, 1, 0);
        glColor3f(0, 0, 1); glVertex3f(0, 0, 0); glVertex3f(0, 0, 1);
        glEnd();

        // Draw left camera coordiante system
        const glm::vec4 o =     left_camera_pose * glm::vec4(0, 0, 0, 1);
        const glm::vec4 x = o + left_camera_pose * glm::vec4(1, 0, 0, 0);
        const glm::vec4 y = o + left_camera_pose * glm::vec4(0, 1, 0, 0);
        const glm::vec4 z = o + left_camera_pose * glm::vec4(0, 0, 1, 0);

        glBegin(GL_LINES);
        glColor3f(1, 0, 0); glVertex3fv(&o.x); glVertex3fv(&x.x);
        glColor3f(0, 1, 0); glVertex3fv(&o.x); glVertex3fv(&y.x);
        glColor3f(0, 0, 1); glVertex3fv(&o.x); glVertex3fv(&z.x);
        glEnd();

        // Draw frustum from left_camera
        drawFrustum(left_view_projection, glm::vec3(1, 1, 1));
    }

    glfwSwapBuffers(window);
}

void framebufferSizeCallback(GLFWwindow* window, int width, int height)
{
    viewport.x      = 0;
    viewport.y      = 0;
    viewport.width  = width;
    viewport.height = height;

    glViewport(0, 0, width, height);

    render(window);
}

glm::vec2 windowToFramebufferCoordinates(GLFWwindow* window, double x, double y)
{
    int window_width, window_height;
    glfwGetWindowSize(window, &window_width, &window_height);

    int fb_width, fb_height;
    glfwGetFramebufferSize(window, &fb_width, &fb_height);

    return glm::vec2(fb_width * x / window_width, fb_height * (1 - y / window_height));
}

void mouseButtonCallback(GLFWwindow* window, int button, int action, int mods)
{
    double x, y;
    glfwGetCursorPos(window, &x, &y);
    previous_cursor_position = windowToFramebufferCoordinates(window, x, y);

    left_is_current = previous_cursor_position.x < (left_aspect / (left_aspect + 1)) * viewport.width;
}

void cursorPosCallback(GLFWwindow* window, double x, double y)
{
    cursor_position                  = windowToFramebufferCoordinates(window, x, y);
    const glm::vec2 cursor_direction = cursor_position - previous_cursor_position;
    previous_cursor_position         = cursor_position;

    if (viewport.width == 0 || viewport.height == 0)
    {
        return;
    }

    const bool left_pressed = glfwGetMouseButton(window, GLFW_MOUSE_BUTTON_LEFT) == GLFW_PRESS;
    if (left_pressed)
    {
        glm::mat4x4& camera_pose = left_is_current ? left_camera_pose : right_camera_pose;

        const float cursor_distance         = glm::length(cursor_direction);
        const float minimum_cursor_distance = 0.1;
        if (cursor_distance < minimum_cursor_distance)
        {
            return;
        }

        glm::mat3x3 rotation;
        for (int x = 0; x < 3; ++x)
        {
            for (int y = 0; y < 3; ++y)
            {
                rotation[x][y] = camera_pose[x][y];
            }
        }

        const glm::vec3 ex = glm::vec3(1, 0, 0);
        const glm::vec3 ey = glm::vec3(0, 1, 0);

        const float viewport_diagonal = glm::length(glm::vec2(viewport.width, viewport.height));
        const float rotation_angle    = 7 * std::atan(cursor_distance / viewport_diagonal);

        const glm::vec3   rotation_axis = glm::normalize(cursor_direction.x * (rotation * ey) - cursor_direction.y * (rotation * ex));
        const glm::mat4x4 dR            = glm::toMat4(glm::angleAxis(rotation_angle, rotation_axis));

        const glm::vec3   pivot = camera_pose[3];
        const glm::mat4x4 t1    = glm::translate(glm::mat4x4(1), -pivot);
        const glm::mat4x4 t2    = glm::translate(glm::mat4x4(1), pivot);

        camera_pose = glm::inverse(t2 * dR * t1) * camera_pose;

        render(window);
    }
}

void keyCallback(GLFWwindow* window, int key, int scancode, int action, int mods)
{
    if (action == GLFW_PRESS || action == GLFW_REPEAT)
    {
        const bool shift_pressed = (mods & GLFW_MOD_SHIFT);

        bool should_render = true;

        if (key == GLFW_KEY_A || key == GLFW_KEY_D || key == GLFW_KEY_W || key == GLFW_KEY_S || key == GLFW_KEY_Q || key == GLFW_KEY_E)
        {
            const float distance    = shift_pressed ? 5 : 1;
            glm::vec3   translation = glm::vec3(0, 0, 0);
            switch (key)
            {
                case GLFW_KEY_A: translation.x -= distance; break;
                case GLFW_KEY_D: translation.x += distance; break;
                case GLFW_KEY_W: translation.z -= distance; break;
                case GLFW_KEY_S: translation.z += distance; break;
                case GLFW_KEY_Q: translation.y -= distance; break;
                case GLFW_KEY_E: translation.y += distance; break;
            }

            glm::mat4x4& camera_pose = left_is_current ? left_camera_pose : right_camera_pose;
            camera_pose = glm::translate(camera_pose, translation);
        }
        else if (key == GLFW_KEY_SPACE)
        {
            right_camera_pose = left_camera_pose;
        }
        else
        {
            should_render = false;
        }

        if (should_render)
        {
            render(window);
        }
    }
}

int main(int argc, char** argv)
{
    if (!glfwInit())
    {
        std::cerr << "Could not initialize GLFW\n";
        return EXIT_FAILURE;
    }

    GLFWwindow* window;
    {
        const int height = 600;
        const int width  = height * left_aspect + height;

        window = glfwCreateWindow(width, height, "Frustum Culling", nullptr, nullptr);
        if (!window)
        {
            std::cerr << "Could not create a glfw window " << width << " " << height << "\n";
            glfwTerminate();
            return EXIT_FAILURE;
        }
    }

    glfwSetFramebufferSizeCallback(window, framebufferSizeCallback);
    glfwSetMouseButtonCallback(window, mouseButtonCallback);
    glfwSetCursorPosCallback(window, cursorPosCallback);
    glfwSetKeyCallback(window, keyCallback);

    // Initialize rendering
    glfwMakeContextCurrent(window);

    glEnable(GL_DEPTH_TEST);
    glClearColor(0.1, 0.1, 0.1, 1);

    aabb.min = glm::vec3(-3, -5, -4);
    aabb.max = glm::vec3(6.5, 6, 4.5);

    left_camera_pose  = glm::inverse(glm::lookAt(glm::vec3(1, 0.5, 20), 0.5f * (aabb.min + aabb.max), glm::vec3(0, 1, 0)));
    right_camera_pose = glm::inverse(glm::lookAt(glm::vec3(10, 5, 30),  0.5f * (aabb.min + aabb.max), glm::vec3(0, 1, 0)));

    left_projection  = glm::perspective(static_cast<float>(45 * M_PI / 180), left_aspect, 5.0f, 30.0f);
    right_projection = glm::perspective(static_cast<float>(60 * M_PI / 180), 1.0f,        0.05f, 100.0f);

    // Render
    {
        int width, height;
        glfwGetFramebufferSize(window, &width, &height);
        framebufferSizeCallback(window, width, height);
    }
    while (!glfwWindowShouldClose(window))
    {
        glfwWaitEvents();
    }

    glfwTerminate();

    return EXIT_SUCCESS;
}

CMakeLists.txt

cmake_minimum_required(VERSION 3.15)
project(frustum_culling)

set(CMAKE_CXX_STANDARD 14)
set(CMAKE_CXX_STANDARD_REQUIRED ON)

find_package(glm CONFIG REQUIRED)

find_package(glfw3 3.3 CONFIG REQUIRED)

find_package(OpenGL REQUIRED)
add_library(opengl INTERFACE)
target_link_libraries(opengl INTERFACE ${OPENGL_gl_LIBRARY})
target_compile_definitions(opengl INTERFACE -DGL_SILENCE_DEPRECATION)

add_executable(frustum_culling frustum_culling.cpp)
target_link_libraries(frustum_culling glfw glm::glm opengl)
  • Left side shows the view from the frustum
  • Right side shows outside views
  • Clicking into the left side moves the frustum
  • Clicking into the right side moves the outside view
  • Use WASDQE to translate the frustum or the outside view (depending on which side was last clicked)

This is also part of frustum_culling.cpp, but this forum has a character limit :slight_smile:

I modified it slightly to add the constants AO, A0, etc

////////////////////////////////////////////////////////////////////////////////////////////////////
// Code generator from https://bivector.net/tools.html
////////////////////////////////////////////////////////////////////////////////////////////////////

// 3D Projective Geometric Algebra
// Written by a generator written by enki.
#include <array>
#include <cmath>
#include <iostream>
#include <stdio.h>
#include <string>

static const char* basis[] = {"1", "e0", "e1", "e2", "e3", "e01", "e02", "e03", "e12", "e31", "e23", "e021", "e013", "e032", "e123", "e0123"};

class PGA3D
{
public:
    PGA3D() { std::fill(mvec, mvec + sizeof(mvec) / 4, 0.0f); }
    PGA3D(float f, int idx = 0)
    {
        std::fill(mvec, mvec + sizeof(mvec) / 4, 0.0f);
        mvec[idx] = f;
    }
    float&       operator[](size_t idx) { return mvec[idx]; }
    const float& operator[](size_t idx) const { return mvec[idx]; }
    PGA3D        log()
    {
        int n = 0;
        for (int i = 0, j = 0; i < 16; i++)
            if (mvec[i] != 0.0f)
            {
                n++;
                printf("%s%0.7g%s", (j > 0) ? " + " : "", mvec[i], (i == 0) ? "" : basis[i]);
                j++;
            };
        if (n == 0) printf("0");
        printf("\n");
        return *this;
    }
    PGA3D log(const std::string& name)
    {
        printf("%s ", name.c_str());
        return this->log();
    }

    PGA3D Conjugate();
    PGA3D Involute();
    float norm();
    float inorm();
    PGA3D normalized();

private:
    float mvec[16];
};

static const int AO    = 0;
static const int A0    = 1;
static const int A1    = 2;
static const int A2    = 3;
static const int A3    = 4;
static const int A01   = 5;
static const int A02   = 6;
static const int A03   = 7;
static const int A12   = 8;
static const int A31   = 9;
static const int A23   = 10;
static const int A021  = 11;
static const int A013  = 12;
static const int A032  = 13;
static const int A123  = 14;
static const int A0123 = 15;

const PGA3D eO    = PGA3D(1, 0);
const PGA3D e0    = PGA3D(1, 1);
const PGA3D e1    = PGA3D(1, 2);
const PGA3D e2    = PGA3D(1, 3);
const PGA3D e3    = PGA3D(1, 4);
const PGA3D e01   = PGA3D(1, 5);
const PGA3D e02   = PGA3D(1, 6);
const PGA3D e03   = PGA3D(1, 7);
const PGA3D e12   = PGA3D(1, 8);
const PGA3D e31   = PGA3D(1, 9);
const PGA3D e23   = PGA3D(1, 10);
const PGA3D e021  = PGA3D(1, 11);
const PGA3D e013  = PGA3D(1, 12);
const PGA3D e032  = PGA3D(1, 13);
const PGA3D e123  = PGA3D(1, 14);
const PGA3D e0123 = PGA3D(1, 15);

//***********************
// PGA3D.Reverse : res = ~a
// Reverse the order of the basis blades.
//***********************
inline PGA3D operator~(const PGA3D& a)
{
    PGA3D res;
    res[0]  = a[0];
    res[1]  = a[1];
    res[2]  = a[2];
    res[3]  = a[3];
    res[4]  = a[4];
    res[5]  = -a[5];
    res[6]  = -a[6];
    res[7]  = -a[7];
    res[8]  = -a[8];
    res[9]  = -a[9];
    res[10] = -a[10];
    res[11] = -a[11];
    res[12] = -a[12];
    res[13] = -a[13];
    res[14] = -a[14];
    res[15] = a[15];
    return res;
};

//***********************
// PGA3D.Dual : res = !a
// Poincare duality operator.
//***********************
inline PGA3D operator!(const PGA3D& a)
{
    PGA3D res;
    res[0]  = a[15];
    res[1]  = a[14];
    res[2]  = a[13];
    res[3]  = a[12];
    res[4]  = a[11];
    res[5]  = a[10];
    res[6]  = a[9];
    res[7]  = a[8];
    res[8]  = a[7];
    res[9]  = a[6];
    res[10] = a[5];
    res[11] = a[4];
    res[12] = a[3];
    res[13] = a[2];
    res[14] = a[1];
    res[15] = a[0];
    return res;
};

//***********************
// PGA3D.Conjugate : res = a.Conjugate()
// Clifford Conjugation
//***********************
inline PGA3D PGA3D::Conjugate()
{
    PGA3D res;
    res[0]  = this->mvec[0];
    res[1]  = -this->mvec[1];
    res[2]  = -this->mvec[2];
    res[3]  = -this->mvec[3];
    res[4]  = -this->mvec[4];
    res[5]  = -this->mvec[5];
    res[6]  = -this->mvec[6];
    res[7]  = -this->mvec[7];
    res[8]  = -this->mvec[8];
    res[9]  = -this->mvec[9];
    res[10] = -this->mvec[10];
    res[11] = this->mvec[11];
    res[12] = this->mvec[12];
    res[13] = this->mvec[13];
    res[14] = this->mvec[14];
    res[15] = this->mvec[15];
    return res;
};

//***********************
// PGA3D.Involute : res = a.Involute()
// Main involution
//***********************
inline PGA3D PGA3D::Involute()
{
    PGA3D res;
    res[0]  = this->mvec[0];
    res[1]  = -this->mvec[1];
    res[2]  = -this->mvec[2];
    res[3]  = -this->mvec[3];
    res[4]  = -this->mvec[4];
    res[5]  = this->mvec[5];
    res[6]  = this->mvec[6];
    res[7]  = this->mvec[7];
    res[8]  = this->mvec[8];
    res[9]  = this->mvec[9];
    res[10] = this->mvec[10];
    res[11] = -this->mvec[11];
    res[12] = -this->mvec[12];
    res[13] = -this->mvec[13];
    res[14] = -this->mvec[14];
    res[15] = this->mvec[15];
    return res;
};

//***********************
// PGA3D.Mul : res = a * b
// The geometric product.
//***********************
inline PGA3D operator*(const PGA3D& a, const PGA3D& b)
{
    PGA3D res;
    res[0]  = b[0] * a[0] + b[2] * a[2] + b[3] * a[3] + b[4] * a[4] - b[8] * a[8] - b[9] * a[9] - b[10] * a[10] - b[14] * a[14];
    res[1]  = b[1] * a[0] + b[0] * a[1] - b[5] * a[2] - b[6] * a[3] - b[7] * a[4] + b[2] * a[5] + b[3] * a[6] + b[4] * a[7] + b[11] * a[8] + b[12] * a[9] + b[13] * a[10] + b[8] * a[11] + b[9] * a[12] + b[10] * a[13] + b[15] * a[14] - b[14] * a[15];
    res[2]  = b[2] * a[0] + b[0] * a[2] - b[8] * a[3] + b[9] * a[4] + b[3] * a[8] - b[4] * a[9] - b[14] * a[10] - b[10] * a[14];
    res[3]  = b[3] * a[0] + b[8] * a[2] + b[0] * a[3] - b[10] * a[4] - b[2] * a[8] - b[14] * a[9] + b[4] * a[10] - b[9] * a[14];
    res[4]  = b[4] * a[0] - b[9] * a[2] + b[10] * a[3] + b[0] * a[4] - b[14] * a[8] + b[2] * a[9] - b[3] * a[10] - b[8] * a[14];
    res[5]  = b[5] * a[0] + b[2] * a[1] - b[1] * a[2] - b[11] * a[3] + b[12] * a[4] + b[0] * a[5] - b[8] * a[6] + b[9] * a[7] + b[6] * a[8] - b[7] * a[9] - b[15] * a[10] - b[3] * a[11] + b[4] * a[12] + b[14] * a[13] - b[13] * a[14] - b[10] * a[15];
    res[6]  = b[6] * a[0] + b[3] * a[1] + b[11] * a[2] - b[1] * a[3] - b[13] * a[4] + b[8] * a[5] + b[0] * a[6] - b[10] * a[7] - b[5] * a[8] - b[15] * a[9] + b[7] * a[10] + b[2] * a[11] + b[14] * a[12] - b[4] * a[13] - b[12] * a[14] - b[9] * a[15];
    res[7]  = b[7] * a[0] + b[4] * a[1] - b[12] * a[2] + b[13] * a[3] - b[1] * a[4] - b[9] * a[5] + b[10] * a[6] + b[0] * a[7] - b[15] * a[8] + b[5] * a[9] - b[6] * a[10] + b[14] * a[11] - b[2] * a[12] + b[3] * a[13] - b[11] * a[14] - b[8] * a[15];
    res[8]  = b[8] * a[0] + b[3] * a[2] - b[2] * a[3] + b[14] * a[4] + b[0] * a[8] + b[10] * a[9] - b[9] * a[10] + b[4] * a[14];
    res[9]  = b[9] * a[0] - b[4] * a[2] + b[14] * a[3] + b[2] * a[4] - b[10] * a[8] + b[0] * a[9] + b[8] * a[10] + b[3] * a[14];
    res[10] = b[10] * a[0] + b[14] * a[2] + b[4] * a[3] - b[3] * a[4] + b[9] * a[8] - b[8] * a[9] + b[0] * a[10] + b[2] * a[14];
    res[11] = b[11] * a[0] - b[8] * a[1] + b[6] * a[2] - b[5] * a[3] + b[15] * a[4] - b[3] * a[5] + b[2] * a[6] - b[14] * a[7] - b[1] * a[8] + b[13] * a[9] - b[12] * a[10] + b[0] * a[11] + b[10] * a[12] - b[9] * a[13] + b[7] * a[14] - b[4] * a[15];
    res[12] = b[12] * a[0] - b[9] * a[1] - b[7] * a[2] + b[15] * a[3] + b[5] * a[4] + b[4] * a[5] - b[14] * a[6] - b[2] * a[7] - b[13] * a[8] - b[1] * a[9] + b[11] * a[10] - b[10] * a[11] + b[0] * a[12] + b[8] * a[13] + b[6] * a[14] - b[3] * a[15];
    res[13] = b[13] * a[0] - b[10] * a[1] + b[15] * a[2] + b[7] * a[3] - b[6] * a[4] - b[14] * a[5] - b[4] * a[6] + b[3] * a[7] + b[12] * a[8] - b[11] * a[9] - b[1] * a[10] + b[9] * a[11] - b[8] * a[12] + b[0] * a[13] + b[5] * a[14] - b[2] * a[15];
    res[14] = b[14] * a[0] + b[10] * a[2] + b[9] * a[3] + b[8] * a[4] + b[4] * a[8] + b[3] * a[9] + b[2] * a[10] + b[0] * a[14];
    res[15] = b[15] * a[0] + b[14] * a[1] + b[13] * a[2] + b[12] * a[3] + b[11] * a[4] + b[10] * a[5] + b[9] * a[6] + b[8] * a[7] + b[7] * a[8] + b[6] * a[9] + b[5] * a[10] - b[4] * a[11] - b[3] * a[12] - b[2] * a[13] - b[1] * a[14] + b[0] * a[15];
    return res;
};

//***********************
// PGA3D.Wedge : res = a ^ b
// The outer product. (MEET)
//***********************
inline PGA3D operator^(const PGA3D& a, const PGA3D& b)
{
    PGA3D res;
    res[0]  = b[0] * a[0];
    res[1]  = b[1] * a[0] + b[0] * a[1];
    res[2]  = b[2] * a[0] + b[0] * a[2];
    res[3]  = b[3] * a[0] + b[0] * a[3];
    res[4]  = b[4] * a[0] + b[0] * a[4];
    res[5]  = b[5] * a[0] + b[2] * a[1] - b[1] * a[2] + b[0] * a[5];
    res[6]  = b[6] * a[0] + b[3] * a[1] - b[1] * a[3] + b[0] * a[6];
    res[7]  = b[7] * a[0] + b[4] * a[1] - b[1] * a[4] + b[0] * a[7];
    res[8]  = b[8] * a[0] + b[3] * a[2] - b[2] * a[3] + b[0] * a[8];
    res[9]  = b[9] * a[0] - b[4] * a[2] + b[2] * a[4] + b[0] * a[9];
    res[10] = b[10] * a[0] + b[4] * a[3] - b[3] * a[4] + b[0] * a[10];
    res[11] = b[11] * a[0] - b[8] * a[1] + b[6] * a[2] - b[5] * a[3] - b[3] * a[5] + b[2] * a[6] - b[1] * a[8] + b[0] * a[11];
    res[12] = b[12] * a[0] - b[9] * a[1] - b[7] * a[2] + b[5] * a[4] + b[4] * a[5] - b[2] * a[7] - b[1] * a[9] + b[0] * a[12];
    res[13] = b[13] * a[0] - b[10] * a[1] + b[7] * a[3] - b[6] * a[4] - b[4] * a[6] + b[3] * a[7] - b[1] * a[10] + b[0] * a[13];
    res[14] = b[14] * a[0] + b[10] * a[2] + b[9] * a[3] + b[8] * a[4] + b[4] * a[8] + b[3] * a[9] + b[2] * a[10] + b[0] * a[14];
    res[15] = b[15] * a[0] + b[14] * a[1] + b[13] * a[2] + b[12] * a[3] + b[11] * a[4] + b[10] * a[5] + b[9] * a[6] + b[8] * a[7] + b[7] * a[8] + b[6] * a[9] + b[5] * a[10] - b[4] * a[11] - b[3] * a[12] - b[2] * a[13] - b[1] * a[14] + b[0] * a[15];
    return res;
};

//***********************
// PGA3D.Vee : res = a & b
// The regressive product. (JOIN)
//***********************
inline PGA3D operator&(const PGA3D& a, const PGA3D& b)
{
    PGA3D res;
    res[15] = 1 * (a[15] * b[15]);
    res[14] = -1 * (a[14] * -1 * b[15] + a[15] * b[14] * -1);
    res[13] = -1 * (a[13] * -1 * b[15] + a[15] * b[13] * -1);
    res[12] = -1 * (a[12] * -1 * b[15] + a[15] * b[12] * -1);
    res[11] = -1 * (a[11] * -1 * b[15] + a[15] * b[11] * -1);
    res[10] = 1 * (a[10] * b[15] + a[13] * -1 * b[14] * -1 - a[14] * -1 * b[13] * -1 + a[15] * b[10]);
    res[9]  = 1 * (a[9] * b[15] + a[12] * -1 * b[14] * -1 - a[14] * -1 * b[12] * -1 + a[15] * b[9]);
    res[8]  = 1 * (a[8] * b[15] + a[11] * -1 * b[14] * -1 - a[14] * -1 * b[11] * -1 + a[15] * b[8]);
    res[7]  = 1 * (a[7] * b[15] + a[12] * -1 * b[13] * -1 - a[13] * -1 * b[12] * -1 + a[15] * b[7]);
    res[6]  = 1 * (a[6] * b[15] - a[11] * -1 * b[13] * -1 + a[13] * -1 * b[11] * -1 + a[15] * b[6]);
    res[5]  = 1 * (a[5] * b[15] + a[11] * -1 * b[12] * -1 - a[12] * -1 * b[11] * -1 + a[15] * b[5]);
    res[4]  = 1 * (a[4] * b[15] - a[7] * b[14] * -1 + a[9] * b[13] * -1 - a[10] * b[12] * -1 - a[12] * -1 * b[10] + a[13] * -1 * b[9] - a[14] * -1 * b[7] + a[15] * b[4]);
    res[3]  = 1 * (a[3] * b[15] - a[6] * b[14] * -1 - a[8] * b[13] * -1 + a[10] * b[11] * -1 + a[11] * -1 * b[10] - a[13] * -1 * b[8] - a[14] * -1 * b[6] + a[15] * b[3]);
    res[2]  = 1 * (a[2] * b[15] - a[5] * b[14] * -1 + a[8] * b[12] * -1 - a[9] * b[11] * -1 - a[11] * -1 * b[9] + a[12] * -1 * b[8] - a[14] * -1 * b[5] + a[15] * b[2]);
    res[1]  = 1 * (a[1] * b[15] + a[5] * b[13] * -1 + a[6] * b[12] * -1 + a[7] * b[11] * -1 + a[11] * -1 * b[7] + a[12] * -1 * b[6] + a[13] * -1 * b[5] + a[15] * b[1]);
    res[0]  = 1 * (a[0] * b[15] + a[1] * b[14] * -1 + a[2] * b[13] * -1 + a[3] * b[12] * -1 + a[4] * b[11] * -1 + a[5] * b[10] + a[6] * b[9] + a[7] * b[8] + a[8] * b[7] + a[9] * b[6] + a[10] * b[5] - a[11] * -1 * b[4] - a[12] * -1 * b[3] - a[13] * -1 * b[2] - a[14] * -1 * b[1] + a[15] * b[0]);
    return res;
};

//***********************
// PGA3D.Dot : res = a | b
// The inner product.
//***********************
inline PGA3D operator|(const PGA3D& a, const PGA3D& b)
{
    PGA3D res;
    res[0]  = b[0] * a[0] + b[2] * a[2] + b[3] * a[3] + b[4] * a[4] - b[8] * a[8] - b[9] * a[9] - b[10] * a[10] - b[14] * a[14];
    res[1]  = b[1] * a[0] + b[0] * a[1] - b[5] * a[2] - b[6] * a[3] - b[7] * a[4] + b[2] * a[5] + b[3] * a[6] + b[4] * a[7] + b[11] * a[8] + b[12] * a[9] + b[13] * a[10] + b[8] * a[11] + b[9] * a[12] + b[10] * a[13] + b[15] * a[14] - b[14] * a[15];
    res[2]  = b[2] * a[0] + b[0] * a[2] - b[8] * a[3] + b[9] * a[4] + b[3] * a[8] - b[4] * a[9] - b[14] * a[10] - b[10] * a[14];
    res[3]  = b[3] * a[0] + b[8] * a[2] + b[0] * a[3] - b[10] * a[4] - b[2] * a[8] - b[14] * a[9] + b[4] * a[10] - b[9] * a[14];
    res[4]  = b[4] * a[0] - b[9] * a[2] + b[10] * a[3] + b[0] * a[4] - b[14] * a[8] + b[2] * a[9] - b[3] * a[10] - b[8] * a[14];
    res[5]  = b[5] * a[0] - b[11] * a[3] + b[12] * a[4] + b[0] * a[5] - b[15] * a[10] - b[3] * a[11] + b[4] * a[12] - b[10] * a[15];
    res[6]  = b[6] * a[0] + b[11] * a[2] - b[13] * a[4] + b[0] * a[6] - b[15] * a[9] + b[2] * a[11] - b[4] * a[13] - b[9] * a[15];
    res[7]  = b[7] * a[0] - b[12] * a[2] + b[13] * a[3] + b[0] * a[7] - b[15] * a[8] - b[2] * a[12] + b[3] * a[13] - b[8] * a[15];
    res[8]  = b[8] * a[0] + b[14] * a[4] + b[0] * a[8] + b[4] * a[14];
    res[9]  = b[9] * a[0] + b[14] * a[3] + b[0] * a[9] + b[3] * a[14];
    res[10] = b[10] * a[0] + b[14] * a[2] + b[0] * a[10] + b[2] * a[14];
    res[11] = b[11] * a[0] + b[15] * a[4] + b[0] * a[11] - b[4] * a[15];
    res[12] = b[12] * a[0] + b[15] * a[3] + b[0] * a[12] - b[3] * a[15];
    res[13] = b[13] * a[0] + b[15] * a[2] + b[0] * a[13] - b[2] * a[15];
    res[14] = b[14] * a[0] + b[0] * a[14];
    res[15] = b[15] * a[0] + b[0] * a[15];
    return res;
};

//***********************
// PGA3D.Add : res = a + b
// Multivector addition
//***********************
inline PGA3D operator+(const PGA3D& a, const PGA3D& b)
{
    PGA3D res;
    res[0]  = a[0] + b[0];
    res[1]  = a[1] + b[1];
    res[2]  = a[2] + b[2];
    res[3]  = a[3] + b[3];
    res[4]  = a[4] + b[4];
    res[5]  = a[5] + b[5];
    res[6]  = a[6] + b[6];
    res[7]  = a[7] + b[7];
    res[8]  = a[8] + b[8];
    res[9]  = a[9] + b[9];
    res[10] = a[10] + b[10];
    res[11] = a[11] + b[11];
    res[12] = a[12] + b[12];
    res[13] = a[13] + b[13];
    res[14] = a[14] + b[14];
    res[15] = a[15] + b[15];
    return res;
};

//***********************
// PGA3D.Sub : res = a - b
// Multivector subtraction
//***********************
inline PGA3D operator-(const PGA3D& a, const PGA3D& b)
{
    PGA3D res;
    res[0]  = a[0] - b[0];
    res[1]  = a[1] - b[1];
    res[2]  = a[2] - b[2];
    res[3]  = a[3] - b[3];
    res[4]  = a[4] - b[4];
    res[5]  = a[5] - b[5];
    res[6]  = a[6] - b[6];
    res[7]  = a[7] - b[7];
    res[8]  = a[8] - b[8];
    res[9]  = a[9] - b[9];
    res[10] = a[10] - b[10];
    res[11] = a[11] - b[11];
    res[12] = a[12] - b[12];
    res[13] = a[13] - b[13];
    res[14] = a[14] - b[14];
    res[15] = a[15] - b[15];
    return res;
};

//***********************
// PGA3D.smul : res = a * b
// scalar/multivector multiplication
//***********************
inline PGA3D operator*(const float& a, const PGA3D& b)
{
    PGA3D res;
    res[0]  = a * b[0];
    res[1]  = a * b[1];
    res[2]  = a * b[2];
    res[3]  = a * b[3];
    res[4]  = a * b[4];
    res[5]  = a * b[5];
    res[6]  = a * b[6];
    res[7]  = a * b[7];
    res[8]  = a * b[8];
    res[9]  = a * b[9];
    res[10] = a * b[10];
    res[11] = a * b[11];
    res[12] = a * b[12];
    res[13] = a * b[13];
    res[14] = a * b[14];
    res[15] = a * b[15];
    return res;
};

//***********************
// PGA3D.muls : res = a * b
// multivector/scalar multiplication
//***********************
inline PGA3D operator*(const PGA3D& a, const float& b)
{
    PGA3D res;
    res[0]  = a[0] * b;
    res[1]  = a[1] * b;
    res[2]  = a[2] * b;
    res[3]  = a[3] * b;
    res[4]  = a[4] * b;
    res[5]  = a[5] * b;
    res[6]  = a[6] * b;
    res[7]  = a[7] * b;
    res[8]  = a[8] * b;
    res[9]  = a[9] * b;
    res[10] = a[10] * b;
    res[11] = a[11] * b;
    res[12] = a[12] * b;
    res[13] = a[13] * b;
    res[14] = a[14] * b;
    res[15] = a[15] * b;
    return res;
};

//***********************
// PGA3D.sadd : res = a + b
// scalar/multivector addition
//***********************
inline PGA3D operator+(const float& a, const PGA3D& b)
{
    PGA3D res;
    res[0]  = a + b[0];
    res[1]  = b[1];
    res[2]  = b[2];
    res[3]  = b[3];
    res[4]  = b[4];
    res[5]  = b[5];
    res[6]  = b[6];
    res[7]  = b[7];
    res[8]  = b[8];
    res[9]  = b[9];
    res[10] = b[10];
    res[11] = b[11];
    res[12] = b[12];
    res[13] = b[13];
    res[14] = b[14];
    res[15] = b[15];
    return res;
};

//***********************
// PGA3D.adds : res = a + b
// multivector/scalar addition
//***********************
inline PGA3D operator+(const PGA3D& a, const float& b)
{
    PGA3D res;
    res[0]  = a[0] + b;
    res[1]  = a[1];
    res[2]  = a[2];
    res[3]  = a[3];
    res[4]  = a[4];
    res[5]  = a[5];
    res[6]  = a[6];
    res[7]  = a[7];
    res[8]  = a[8];
    res[9]  = a[9];
    res[10] = a[10];
    res[11] = a[11];
    res[12] = a[12];
    res[13] = a[13];
    res[14] = a[14];
    res[15] = a[15];
    return res;
};

//***********************
// PGA3D.ssub : res = a - b
// scalar/multivector subtraction
//***********************
inline PGA3D operator-(const float& a, const PGA3D& b)
{
    PGA3D res;
    res[0]  = a - b[0];
    res[1]  = -b[1];
    res[2]  = -b[2];
    res[3]  = -b[3];
    res[4]  = -b[4];
    res[5]  = -b[5];
    res[6]  = -b[6];
    res[7]  = -b[7];
    res[8]  = -b[8];
    res[9]  = -b[9];
    res[10] = -b[10];
    res[11] = -b[11];
    res[12] = -b[12];
    res[13] = -b[13];
    res[14] = -b[14];
    res[15] = -b[15];
    return res;
};

//***********************
// PGA3D.subs : res = a - b
// multivector/scalar subtraction
//***********************
inline PGA3D operator-(const PGA3D& a, const float& b)
{
    PGA3D res;
    res[0]  = a[0] - b;
    res[1]  = a[1];
    res[2]  = a[2];
    res[3]  = a[3];
    res[4]  = a[4];
    res[5]  = a[5];
    res[6]  = a[6];
    res[7]  = a[7];
    res[8]  = a[8];
    res[9]  = a[9];
    res[10] = a[10];
    res[11] = a[11];
    res[12] = a[12];
    res[13] = a[13];
    res[14] = a[14];
    res[15] = a[15];
    return res;
};

inline float PGA3D::norm() { return sqrt(std::abs(((*this) * Conjugate()).mvec[0])); }
inline float PGA3D::inorm() { return (!(*this)).norm(); }
inline PGA3D PGA3D::normalized() { return (*this) * (1 / norm()); }

// A rotor (Euclidean line) and translator (Ideal line)
static PGA3D rotor(float angle, PGA3D line) { return cos(angle / 2.0f) + sin(angle / 2.0f) * line.normalized(); }
static PGA3D translator(float dist, PGA3D line) { return 1.0f + dist / 2.0f * line; }

// PGA is plane based. Vectors are planes. (think linear functionals)
//static PGA3D e0(1.0f, 1), e1(1.0f, 2), e2(1.0f, 3), e3(1.0f, 4);

// A plane is defined using its homogenous equation ax + by + cz + d = 0
static PGA3D plane(float a, float b, float c, float d) { return a * e1 + b * e2 + c * e3 + d * e0; }

// PGA points are trivectors.
//static PGA3D e123 = e1 ^ e2 ^ e3, e032 = e0 ^ e3 ^ e2, e013 = e0 ^ e1 ^ e3, e021 = e0 ^ e2 ^ e1;

// A point is just a homogeneous point, euclidean coordinates plus the origin
static PGA3D point(float x, float y, float z) { return e123 + x * e032 + y * e013 + z * e021; }

// for our toy problem (generate points on the surface of a torus)
// we start with a function that generates motors.
// circle(t) with t going from 0 to 1.
static PGA3D circle(float t, float radius, PGA3D line)
{
    return rotor(t * 2.0f * M_PI, line) * translator(radius, e1 * e0);
}

// a torus is now the product of two circles.
static PGA3D torus(float s, float t, float r1, PGA3D l1, float r2, PGA3D l2)
{
    return circle(s, r2, l2) * circle(t, r1, l1);
}

// and to sample its points we simply sandwich the origin ..
static PGA3D point_on_torus(float s, float t)
{
    PGA3D to = torus(s, t, 0.25f, e1 * e2, 0.6f, e1 * e3);
    return to * e123 * ~to;
}

////////////////////////////////////////////////////////////////////////////////////////////////////