README.md

Previous part: Geometric algebra in 2D and 3D

Disclaimer: there will be NO proofs or rigorous explanations here, just some results and notes on how to use the library

Projective geometric algebra in 2D

Lines

We start with lines. As you can remember, the equation for a line in 2D looks like this:

$$ Ax + By + C = 0 $$

In order to do operations on lines, as elements of the algebra, let's convert it to the form where it will be linear combination of basis vectors:

$$ A e_1 + B e_2 + C e_0 $$

Here we introduce new kind of basis vector $e_0$. It's special, because $e_0^2 = 0$ by definition. Introduction of $e_0$ allows us to make projective geometry out of usual rules of geometric algebra, which we learned in the previous blog post.

For example, $e_1 - e_2 - e_0$ represents line $x - y - 1 = 0$

Let's see what elements this algebra consists of:

       1 
 e1   e2   e0
e20  e01  e12
    e120

As elements of the algebra, we can add lines together:

$$ (A_1 e_1 + B_1 e_2 + C_1 e_0) + (A_2 e_1 + B_2 e_2 + C_2 e_0) = (A_1 + A_2) e_1 + (B_1 + B_2) e_2 + (C_1 + C_2) e_0 $$

Sum of two lines will be some other line, which goes through their intersection point.

For example, $A = e_1 + e_2 + 2 e_0$, $B = 5 e_1 - 5 e_2$:

One thing to note here: $A$ and $\alpha A$ represent the same line! So representation of the line in the algebra is not unique.

In code I define struct that represents line:

struct line2 {
    float32 a, b, c;
};

Basis "vectors" for lines

Let's see what basis lines $e_1$, $e_2$, and $e_0$ represent.

$$ e_1 \mapsto x = 0 $$

This is line Oy.

$$ e_2 \mapsto y = 0 $$

This is line Ox.

$$ e_0 \mapsto 1 = 0 $$

Wait a minute, this is nonsense! But if we pause for a minute, this is very useful representation for our purposes, because it works out exactly as the line at infinity!

You have to imagine line that goes around whole plane at infinite distance, being perpendicular to every other "real" line.

Points

Beside lines we also want to represent points. Turns out, outer product of two lines results in intersection point!

$$ \begin{aligned} L_1 \wedge L_2 & = (A_1 e_1 + B_1 e_2 + C_1 e_0) \wedge (A_2 e_1 + B_2 e_2 + C_2 e_0) = \\ & = (A_1 B_2 - A_2 B_1) e_{12} + (B_1 C_2 - B_2 C_1) e_{20} + (C_1 A_2 - C_2 A_1) e_{01} \end{aligned} $$

This operation is often called meet because it returns a point, where two lines meet.

This is exactly the representation of homogeneous coordinates in 2D:

$$ (x, y, w) = (x / w, y / w, 1) \mapsto xe_{20} + ye_{01} + we_{12} $$

In code points are represented just like you expect:

struct point2 {
    float32 x, y, w;
};

and you can get an intersection of two lines by using outer product:

point2 outer(line2, line2);

But what if we want to intersect two parallel lines? Let's calculate intersection point of $x - y - 1 = 0$ and $x - y + 1 = 0$.

$$ (e_1 - e_2 - e_0) \wedge (e_1 - e_2 + e_0) = -2 e_{20} - 2 e_{01} \mapsto (-2, -2, 0) $$

auto line_1 = make_line2(1, -1, -1); // x - y - 1 = 0
auto line_2 = make_line2(1, -1,  1); // x - y + 1 = 0

auto intersection = outer(line_1, line_2);

Remember how we learned that in homogeneous coordinates (x, y, 1) represent point, while (x, y, 0) represent vectors? Turns out we rediscovered those vectors! In projective geometric algebra vectors are points at infinity:

Duality

Duality is somewhat new concept for us, in geometric algebra each line have a dual point, and each point have a dual line. Usually dual of a multivector is defined by multiplying by something, but I just replace basis vectors with corresponding dual basis:

$$ \begin{aligned} e_1 & \leftrightarrow e_{20} \\ e_2 & \leftrightarrow e_{01} \\ e_0 & \leftrightarrow e_{12} \\ \end{aligned} $$

We see that dual of a line is just point, where normal vector points to, or just normal vector, if line goes through the origin.

Let's see a couple of examples:

The most amazing feat of duality, is that everything can be turn around, but math will work out just the same. For example, in the Projective geometric algebra Wiki basis elements are swapped onto duals comparing to our implementation, but every result still holds.

Connecting two points with a line

Now when we know and understand duality, we can connect two points with a line.

For example, let's draw two points: $(-3, 2)$ and $(3, 1)$. I also draw their dual lines:

The trick goes like this: we see that dual lines cross at some point, let's find this point:

auto point_1 = make_point2(-3, 2);
auto point_2 = make_point2( 3, 1);

auto dual_line_1 = dual(point_1);
auto dual_line_2 = dual(point_2);

auto intersection = outer(dual_line_1, dual_line_2);

Now let's find dual of that intersection point:

auto dual_of_intersection = dual(intersection);

By finding intersection (point that connects two lines) in dual space, we also found line that connects two points in real space. If you ask me, this is kinda magic.

In projective geometric algebra this whole operation is very useful, so they gave it a name regressive product, and write like this:

$$ A \vee B = \mathrm{dual}(\mathrm{dual}(A) \wedge \mathrm{dual}(B)) $$

This operation is often called join because it joins two points by a line.

In code you can use it by calling this function:

line2 join(point2, point2);

Projections

In projective geometric algebra we can define inner and outer product not only between vectors, but between any two objects!

$$ A \wedge B = \sum_{i=0}^n \sum_{j=0}^n \langle AB \rangle_{i + j} $$

$$ A \cdot B = \sum_{i=0}^n \sum_{j=0}^n \langle AB \rangle_{\lvert i - j \rvert} $$

The $\langle A \rangle_i$ is called "i-th blade". This operation strips multivector from all its elements, but elements of i-th dimensions.

For example, $\langle A \rangle_0$ is just scalar part of the multivector, and $\langle A \rangle_2$ is just bivector part (elements with $e_{12}$, $e_{20}$, and $e_{01}$).

Now let's see what the inner product of a line and a point will result to:

$$ (Ae_1 + Be_2 + Ce_0) \cdot (xe_{20} + ye_{01} + we_{12}) = $$

$$ = -Bw e_1 + Aw e_2 + (Bx - Ay) e_0 $$

This is a line, that goes through a point, but perpendicular to the original line!

Now we can use it to make two projections. First, let's observe that if we find intersection of those two lines, we find projection of a point to a line:

auto line = make_line2(-1, 2, 3);
auto point = make_point2(-1, 2);

auto projected_point = outer(inner(line, point), line);

But if we decided to find the perpendicular line again, we could find so called projection of a line onto a point:

auto projected_line = inner(inner(point, line), point);

In projective geometric algebra all kinds of projections can be found by applying this formula:

$$ \mathrm{project} A \mathrm{onto} B = (B \cdot A)B $$

but in code, in order to match types, I changed it to inner(inner(B, A), B) and outer(inner(B, A), B). You also can use functions:

point2 project(point2, line2);
line2 project(line2, point2);

Reflections

In geometric algebra reflections through the line are computed by the formula:

$$ T' = L T L $$

Let's reflect a line $L_1 = A_1 e_1 + B_1 e_2 + C_1 e_0$ through another line $L_2 = A_2 e_1 + B_2 e_2 + C_2 e_0$.

$$ L_2 L_1 L_2 = $$

$$ = (A_2 e_1 + B_2 e_2 + C_2 e_0) (A_1 e_1 + B_1 e_2 + C_1 e_0) (A_2 e_1 + B_2 e_2 + C_2 e_0) = $$

$$ = ((A_2 A_1 + B_2 B_1) + (A_2 B_1 - B_2 A_1) e_{12} + (B_2 C_1 - C_2 B_1) e_{20} + (C_2 A_1 - A_2 C_1) e_{01}) (A_2 e_1 + B_2 e_2 + C_2 e_0) $$

$$ \begin{aligned} = (A_1 A_2 A_2 + B_1 A_2 B_2 + B_1 A_2 B_2 - A_1 B_2 B_2) & e_1 + \\ (A_1 A_2 B_2 + B_1 B_2 B_2 - B_1 A_2 A_2 + A_1 A_2 B_2) & e_2 + \\ (A_1 A_2 C_2 + B_1 B_2 C_2 - C_1 B_2 B_2 + B_1 B_2 C_2 + A_1 A_2 C_2 - C_1 A_2 A_2) & e_0 + \\ (B_1 A_2 C_2 - A_1 B_2 C_2 + C_1 A_2 B_2 - B_1 A_2 C_2 + A_1 B_2 C_2 - C_1 A_2 B_2) & e_{120} \end{aligned} $$

Note that coefficient besides $e_{120}$ is equal to 0.

$$ L_2 L_1 L_2 = S e_1 + T e_2 + R e_0 $$

In code we can't just elliminate one of the coefficient, so element that we got is like a line, but with pseudoscalar. This object is often called flector because it reflects things.

struct flector {
    float32 _e1, _e2, _e0, _e120;
};

But because we know for sure that coefficient besides $e_{120}$ will be equal to 0, we can just discard it in code.

line2 get_line2(flector f) {
    line2 result;
    result._e1 = f._e1;
    result._e2 = f._e2;
    result._e0 = f._e0;
    return result;
}

auto L1 = make_line2(-1, 1, -2);
auto L2 = make_line2(-1, 3, -2);

auto L_prime = get_line2(L2 * L1 * L2);

Rotations

Just like a line with pseudoscalar is a flector, likewise a point with a scalar is a special kind of object, called rotor or motor.

Let's multiply two lines, not with outer product, but with geometric product:

$$ \begin{aligned} (A_1 e_1 + B_1 e_2 & + C_1 e_0) (A_2 e_1 + B_2 e_2 + C_2 e_0) = \\ = (A_1 A_2 & + B_1 B_2 + C_1 C_2) + \\ (A_1 B_2 & - B_1 A_2) e_{12} + \\ (C_1 A_2 & - C_2 A_1) e_{01} + \\ (B_1 C_2 & - B_2 C_1) e_{20} \end{aligned} $$

This object will rotate objects around the point, just like quaternion does, twice the angle between original two lines. Also just like quaternions, you can apply rotation using this formula:

$$ A' = M^\dagger A M $$

Where $M^\dagger$ is conjugated motor.

For example, I made two lines 45 degrees apart each other, multiplied to make a motor, and then rotated a line by it:

auto line_1 = make_line2(1, -1, -1);
auto line_2 = make_line2(0,  1, -1);

auto motor = line_1 * line_2;

auto line_3 = make_line2(1, 0, 2);
auto line_4 = to_line2(conjugated(motor) * line_3 * motor);

You can rotate points just like lines:

auto line_1 = make_line2(1, -1, -1);
auto line_2 = make_line2(0,  1, -1);

auto motor = line_1 * line_2;

auto point_1 = make_point2(-2, 1);
auto point_2 = to_point2(conjugated(motor) * point_1 * motor);

Why some people call motor a motor, and not a rotor, is because in projective geometric algebra you can use it to make translations and rotations at the same time. You can achieve translation by making motor out of two parallel lines. Then it will translate any object by twice the distance between those two lines.

For example:

auto line_1 = make_line2(1,  1, -1);
auto line_2 = make_line2(1,  1,  1);

auto motor = line_1 * line_2;

auto point_1 = make_point2(-1, 2);
auto point_2 = to_point2(conjugated(motor) * point_1 * motor);

And that's it for 2D case. You can find my library which can do everything I have shown here and the examples all working here.

References:

JustSlavic/README.md