main.c

/*
    References:
    1)http://www.gaalop.de/dhilden_data/CLUScripts/CandG_PrePrint.pdf
    2)http://www2.montgomerycollege.edu/departments/planet/planet/Numerical_Relativity/GA-SIG/Conformal%20Geometry%20Papers/Cambridge/Covarient%20Approach%20to%20Geometry%20Using%20Geometric%20Algebra.pdf
    3)http://peeterjoot.com/archives/geometric-algebra/bivector.pdf
    4)Geometric Algebra for Computer Science: An Object-Oriented Approach to Geometry by Leo Dorst, Daniel Fontijne and Stephen Mann
    5)Multivectors and Clifford Algebra in Electrodynamics By Bernard Jancewicz

    The following code has been written to the best of my understanding of the above concept.
    However, I'm not sure of my math in the function where I compute the inner product of a trivector and a vector.
    Formula from the 5th reference on page 24: (b ^ c ^ d) = b ^ c(d . a) + d ^ b(c . a) + c ^ d(b . a)
    Doing the above for 10 components for 5 times each was painful and prone to errors.
    Sorry for the ugly long names, tried my best to make the code self explanatory.
    Will be more than happy to explain any part which is unclear in the code.
*/

#include <stdio.h>

// vec3. To be read as "x y z"
typedef float vec3[3];

void vec3_extend_vec3(vec3 res, vec3 a, vec3 b);
void vec3_add_vec3(vec3 res, vec3 a, vec3 b);

// Conformal geometric algebra 5d vector notation of 3d vector. To be read as "e1 e2 e3 e0 einf" where e0 is point at origin and einf is point at infinity.
typedef float vec5[5];

// Conformal geometric algebra 10 bivector notation. To be read as "e1e2 e1e3 e1e0 e1einf e2e3 e2e0 e2einf e3e0 e3einf e0einf" where eij = ei ^ ej.
// Conformal geometric algebra 10 trivector notation. To be read as "e1e2e3 e1e2e0 e1e2einf e1e3e0 e1e3einf e1e0einf e2e3e0 e2e3einf e2e0einf e3e0einf"
typedef float vec10[10];

void vec3_to_vec5(vec5 res, vec3 v);
void vec5_outer_vec5(vec10 res, vec5 a, vec5 b);

// bivector ^ vector = trivector (res has 10 components)
void vec10_outer_vec5(vec10 res, vec10 a, vec5 b);

// trivector . vector = bivector (res has 10 components)
void vec10_inner_vec5(vec10 res, vec10 a, vec5 b);

// bivector . bivector = scalar
float vec10_inner_vec10(vec10 a, vec10 b);

void create_sphere(vec5 sphere, vec3 sphere_pos_v3, float rad);
void create_line(vec5 line, vec3 line_pos_v3, vec3 line_dir_v3);

/*
    point_pair(bivector) = line(trivector) . sphere(vector)
    sign(scalar) = point_pair . point_pair
*/
void compute_intersection(vec10 line, vec5 sphere);

int main(void)
{
    vec3 sphere_pos_v3;
    printf("Enter sphere pos in vec3: ");
    scanf("%f %f %f", &sphere_pos_v3[0], &sphere_pos_v3[1], &sphere_pos_v3[2]);

    float rad;
    printf("Enter sphere radius: ");
    scanf("%f", &rad);

    vec5 sphere;
    create_sphere(sphere, sphere_pos_v3, rad);

    vec3 line_pos_v3, line_dir_v3;
    printf("Enter line pos in vec3: ");
    scanf("%f %f %f", &line_pos_v3[0], &line_pos_v3[1], &line_pos_v3[2]);
    printf("Enter line dir in vec3: ");
    scanf("%f %f %f", &line_dir_v3[0], &line_dir_v3[1], &line_dir_v3[2]);

    vec3 line_dir_ext;
    vec3_extend_vec3(line_dir_ext, line_dir_v3, sphere_pos_v3);
    vec3 line_endpos_v3;
    vec3_add_vec3(line_endpos_v3, line_pos_v3, line_dir_ext);

    vec10 line;
    create_line(line, line_pos_v3, line_endpos_v3);

    compute_intersection(line, sphere);

    return 0;
}

void vec3_extend_vec3(vec3 res, vec3 a, vec3 b)
{
    res[0] = a[0] * b[0];
    res[1] = a[1] * b[1];
    res[2] = a[2] * b[2];
}

void vec3_add_vec3(vec3 res, vec3 a, vec3 b)
{
    res[0] = a[0] + b[0];
    res[1] = a[1] + b[1];
    res[2] = a[2] + b[2];
}

void vec3_to_vec5(vec5 res, vec3 v)
{
    res[0] = v[0];
    res[1] = v[1];
    res[2] = v[2];
    res[3] = 1.0f;
    res[4] = (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]) / 2.0f;
}

void vec5_outer_vec5(vec10 res, vec5 a, vec5 b)
{
    res[0] = (a[0] * b[1]) - (a[1] * b[0]);
    res[1] = (a[0] * b[2]) - (a[2] * b[0]);
    res[2] = (a[0] * b[3]) - (a[3] * b[0]);
    res[3] = (a[0] * b[4]) - (a[4] * b[0]);

    res[4] = (a[1] * b[2]) - (a[2] * b[1]);
    res[5] = (a[1] * b[3]) - (a[3] * b[1]);
    res[6] = (a[1] * b[4]) - (a[4] * b[1]);

    res[7] = (a[2] * b[3]) - (a[3] * b[2]);
    res[8] = (a[2] * b[4]) - (a[4] * b[2]);

    res[9] = (a[3] * b[4]) - (a[4] * b[3]);
}

void vec10_outer_vec5(vec10 res, vec10 a, vec5 b)
{
    res[0] = (a[0] * b[2]) - (a[1] * b[1]) - (a[4] * b[0]);
    res[1] = (a[0] * b[3]) - (a[2] * b[1]) - (a[5] * b[0]);
    res[2] = (a[0] * b[4]) - (a[3] * b[1]) - (a[6] * b[0]);

    res[3] = (a[1] * b[3]) - (a[2] * b[2]) - (a[7] * b[0]);
    res[4] = (a[1] * b[4]) - (a[3] * b[2]) - (a[8] * b[0]);

    res[5] = (a[2] * b[4]) - (a[3] * b[3]) - (a[9] * b[0]);

    res[6] = (a[4] * b[3]) - (a[5] * b[2]) - (a[7] * b[1]);
    res[7] = (a[4] * b[4]) - (a[6] * b[2]) - (a[8] * b[1]);

    res[8] = (a[5] * b[4]) - (a[6] * b[3]) - (a[9] * b[1]);

    res[9] = (a[7] * b[4]) - (a[8] * b[3]) - (a[9] * b[2]);
}

void vec10_inner_vec5(vec10 res, vec10 a, vec5 b)
{
    res[0] = (a[0] * b[2]) - (a[1] * b[3]) - (a[2] * b[4]);
    res[1] = -(a[0] * b[1]) - (a[3] * b[3]) - (a[4] * b[4]);
    res[2] = -(a[1] * b[1]) - (a[3] * b[2]) - (a[5] * b[4]);
    res[3] = -(a[2] * b[1]) - (a[4] * b[2]) - (a[5] * b[3]);

    res[4] = (a[0] * b[0]) - (a[6] * b[3]) - (a[7] * b[4]);
    res[5] = (a[1] * b[0]) - (a[6] * b[2]) - (a[8] * b[4]);
    res[6] = (a[2] * b[0]) - (a[7] * b[2]) + (a[8] * b[3]);

    res[7] = (a[3] * b[0]) + (a[6] * b[1]) - (a[9] * b[4]);
    res[8] = (a[4] * b[0]) + (a[7] * b[1]) + (a[9] * b[3]);

    res[9] = (a[5] * b[0]) + (a[8] * b[1]) + (a[9] * b[2]);
}

float vec10_inner_vec10(vec10 a, vec10 b)
{
    return (-(a[0] * b[0]) - (a[1] * b[1]) + (a[2] * b[2]) + (a[3] * b[3]) - (a[4] * b[4]) + (a[5] * b[5]) + (a[6] * b[6]) + (a[7] * b[7]) + (a[8] * b[8]) - (a[9] * b[9]));
}

void create_sphere(vec5 sphere, vec3 sphere_pos_v3, float rad)
{
    vec5 sphere_pos_v5;
    vec3_to_vec5(sphere_pos_v5, sphere_pos_v3);

    sphere[0] = sphere_pos_v5[0];
    sphere[1] = sphere_pos_v5[1];
    sphere[2] = sphere_pos_v5[2];
    sphere[3] = sphere_pos_v5[3];
    sphere[4] = sphere_pos_v5[4] - (rad * rad / 2.0f);
}

void create_line(vec10 line, vec3 line_pos_v3, vec3 line_endpos_v3)
{
    vec5 line_pos_v5, line_endpos_v5;
    vec3_to_vec5(line_pos_v5, line_pos_v3);
    vec3_to_vec5(line_endpos_v5, line_endpos_v3);

    vec10 temp;
    vec5_outer_vec5(temp, line_pos_v5, line_endpos_v5);

    vec5 point_at_einf = { 0.0f, 0.0f, 0.0f, 0.0f, 1.0f };
    vec10_outer_vec5(line, temp, point_at_einf);
}

void compute_intersection(vec10 line, vec5 sphere)
{
    vec10 point_pair;
    vec10_inner_vec5(point_pair, line, sphere);
    float sign = vec10_inner_vec10(point_pair, point_pair);
    if (sign > 0)
        printf("True.\n");
    else
        printf("False.\n");
}