About 3D PGA Cheat Sheet

meuns · October 27, 2019, 8:12pm

Hi,

I implement a few things of PGA using galgebra and I’m trying to write some tests for each line of the 3D PGA cheat sheet.

I need some insights about norms. Most of the time we are computing line norms and we know how to handle ideal and euclidean lines. But when we takes the norm of something else, the process is very blurry. I can’t use \|M\|^2 = \tilde{M}M and hope for a scalar…

The distance of a point to a line given by \|\hat{P}\vee\hat{l}\|. How should I compute the norm of a plane ?

The distance of a point to a plane is given by \|\hat{P}\wedge\hat{p}\|_{\infty} which is always the pseudo-scalar scaled by the oriented distance. How should I compute the norm of the pseudo-scalar ?

Regards

chakravala · October 28, 2019, 3:42am

The Euclidean norm of a multivector is just simply the standard vector norm applied to the entire set of coefficient values in the multivector.

\displaystyle || M||_\infty^2 = \sum_{g=0}^n ||\langle M\rangle_g||_\infty^2 = \sum_{g=0}^n \sum_{k=1}^{n\choose g} \langle M\rangle_{g,k}^2

Note that \langle M\rangle_{g,k} should be the k-th scalar coefficient without the \Lambda^g basis.

ninepoints · October 28, 2019, 4:16am

@meuns
The course notes has a helpful section regarding the normalization of a motor (pay attention to the section where the closed form of the exponential and logarithm of a bivector and motor are derived specifically, 8.1.3 and 8.1.6.

It mentions normalizing a bivector, but a similar trick (I believe) can be used for motors. It should be self-evident that M\tilde{M} will be of the form u + v\mathbf{I} (M occupies the even subalgebra). Thus, we need to divide by \sqrt{u + v\mathbf{I}} to normalize M properly.

The first trick is to use the fact that \mathbf{I} behaves like a dual number to complete the square:

\begin{aligned}\sqrt{u + v\mathbf{I}} &= \sqrt{u'^2 + 2u'v'\mathbf{I} + v'^2\mathbf{I}^2} \\ &=\sqrt{(u' + v'\mathbf{I})^2}\\ &= u' + v'\mathbf{I} \end{aligned}

where we have added zero (\mathbf{I}^2) in the first step along with the following change of variables:

\begin{aligned} u' &= \sqrt{u} \\ v' &= \frac{v}{2\sqrt{u}} \end{aligned}

Then, we just need to divide by the new quantity. If the new quantity is of the form s + t\mathbf{I}, you can multiply this to u' + v'\mathbf{I}, set equal to one, and solve for s and t in terms of u and v.

\begin{aligned} 1 &= (u' + v'\mathbf{I})(s + t\mathbf{I}) \\ &= u's + (v's + u't)\mathbf{I} \\ \Rightarrow s &= \frac{1}{u'} = \frac{1}{\sqrt{u}} \\ \Rightarrow t &= -\frac{v'}{u'^2} = -\frac{v}{2u\sqrt{u}} \end{aligned}

Multiplying to your original motor M by \frac{1}{\sqrt{u}} - \frac{v}{2u\sqrt{u}}\mathbf{I} should result in a normalized M such that M\tilde{M} = 1.

Disclaimer: I haven’t implemented this yet, but it’s on my list.

edit:

Oh! By the way, an important fact is that if you construct the motor by exponentiating a bivector as outlined in the course notes, the normalization is already accounted for (the Euclidean axis is normalized via the sine and cosine transcendentals, while the polar (ideal) axis is scaled by a dual which extinguishes itself to perform the translation portion of the screw motion. In my case, I will most likely be creating motors in this fashion although I intend to implement the motor normalization eventually (just for numerical stability).

ninepoints · October 28, 2019, 4:45am

@meuns I’ve added motor normalization code here: https://github.com/jeremyong/gal/blob/master/public/gal/pga.hpp#L270

You can see a test demonstrating that it works here:

github.com

jeremyong/gal/blob/master/test/test_pga.cpp#L90


        }


        auto m_norm = compute([](auto m) { return m * ~m; }, m2);
        printf("m_norm: %f, %f\n", m_norm[0], m_norm[1]);
        CHECK_EQ(m_norm[0], doctest::Approx(1));
        CHECK_EQ(m_norm[1], doctest::Approx(0));
    }


    SUBCASE("normalize-motor")
    {
        motor<> m{{293.2, -39.3, 59.3, -1.04, 434.3, 23.0, 72.874}};
        m.normalize();
        auto m_norm = compute([](auto m) { return m * ~m; }, m);
        printf("m_norm: %f, %f\n", m_norm[0], m_norm[1]);
        CHECK_EQ(m_norm[0], doctest::Approx(1));
        CHECK_EQ(m_norm[1], doctest::Approx(0));
    }
}


TEST_SUITE_END();

Thanks for reminding me to get to this

cgunn3 · October 28, 2019, 7:28pm

These are good and important questions, @meuns. I am assuming you are interested in the euclidean case. The non-euclidean case can be handled similarly but there are different special cases.

Good news: for a blade b, b^2 is a scalar (since it’s equal to \pm b \widetilde{b} and this product collapses as a product of squares of 1-vectors, which are by definition scalars).

Since points and planes are both simple (can be written as blades), their norms can be easily calculated in this way. For a euclidean plane p, p^2>0 and \lVert p \rVert = \sqrt{p^2}. Similarly for a euclidean point P, P^2 < 0 and \lVert P \rVert = \sqrt{-P^2}.

A pseudo-scalar is a little different since its “norm” actually is a signed magnitude (while a norm is always non-negative). This is because the vector space of pseudo-scalars is 1-dimensional; once you have defined the unit pseudo-scalar \bf{I} (typically as the product of all the basis 1-vectors) then any pseudo-scalar can be written as \alpha \bf{I} for some non-zero real number \alpha. Then define \lVert \alpha \bf{I} \rVert_\infty = \alpha. Here we use the \infty subscript since this is the so-called “ideal” norm (defined for ideal elements satisfying x^2=0, such as the pseudo-scalar).

There are also ideal norms for ideal points and planes, whose standard norms are 0 and hence cannot be normalized. See Sec. 7.1 of the course notes for details.