Plane-Based Geometric Algebra: A Practical Guide for GPU Computing

[00:00:00]hello my name is leo durst and there are some exciting developments in the computational representation of euclidean geometry that i'd like to share with you in this talk it's about pga plane based geometric algebra it is an algebra but it's a practical algebra it's a way of organizing your computations in geometry so that primitives and operators are nicely integrated in a way that actually makes them very suitable for implementation on gpus we have universal operators called motors that can move any primitive or motor in exactly the same way exactly the same instruction which enormously saves codes these operators can be generated directly from the geometry so for instance if i want a rotation operation that rotates around the line l with an angle v and i can just do the exponential of that line as an object times phi we have simple expressions to combine geometric constructions so if you have a point and a line and you would like to plane through both then the dot product of the line and the point is that plane all this is coordinate free and can be specified quite generally in less code just because the universal operators require only one line to state what they should do which moreover can be implemented to execute more efficiently than homogeneous coordinates or dual quaternions as a geometric algebra the pga is organized rather unusually to what you're used to normally you would define geometric objects first and then define the operators here it's the other way around so like any geometry the motions that symmetries define what the geometry is we generate the motions as multiple reflections uh for leocalydian geometry that means multiple reflections in planes because two parallel planes give a translation two intersecting planes give a rotation around the intersection line and that means that algebraically we're going to have to let the geometric algebra work on the space of which planes are the vectors it's a linear space with a dog product we call it pga for plane based geometric algebra when we look at that in detail in a few slides we'll see that the fact that we want reflections motivate defining a specific product the geometric product that makes it a very straightforward thing to do and then we find that the geometric product is also capable of intersecting elements we already had the planes as reflectors now we get the planes to intersect to become lines and points so those are the constituents of making objects in the geometry and in this way from motion to objects we see that the motion representation dictates the object representation the objects can then in turn be used to parameterize the motions like we've seen that the rotation around the line is just the exponential of a line and the consequence is that all the elements that we need all the things we need to do euclidean geometry are present in an integrated way the little demo might might illustrate this so we're going to see a tricycle move um the tricycle is determined by the fact that it has a rear axle and a front axle that can be steered over this angle when they intersect they create the point around which the rear axle rotates when it moves and in geometric algebra that'll be just an operator that computes the intersection product then the exponential to make the motion and then the motor the motion operator can be applied to any element of the body of the bicycle the vertices the edges the planes the texture and they all move in exactly the same way and when we see it move you'll see that as the rotation angle changes from straight to left to right there is just no singularity if this intersection is an ideal point at infinity the exponential of an ideal point it just will turn out to be a translation so we're first going to see how we have to represent motions and this is the basic principle behind it we are going to map our present world into a representational space in which the motions are represented as multiple reflections that means that in that space they are orthogonal transformations in order to preserve the dot product and we are helped by a theorem by carton and joe donais that if we have an n-dimensional symmetric bi-linear space then all the orthogonal transformations in that space can be generated from reflections so once we're into orthogonal transformations this is the general principle reflections are represented normally in by this formula here in linear algebra let's see how it comes about you have a vector space you have a vector there's a plane in this space with a normal vector a you'd like to reflect this vector in the plane what do you do well a raw is being played by this projection of x onto a and you have to subtract that twice from x to get to the reflected vector so this is the operation x minus 2 times and this is the projection operator if this is a unit vector just this x dot a times a if it's not a unit vector we have to divide by a norm of a squared now this is a formula that you cannot really use to organize your motions it's just too involved but that is due to the occurrence of the dot product the dot product is symmetrical and the idea by clifford was to replace the dot product by something that you say it's the symmetrical part of a deeper product and a deeper product is called the geometric product it is denoted by just the space and so if we take xa plus ax that's clearly symmetrical divided by 2 then that is what the dot product is okay so if we substitute this equation into that we find that here we have that replacement we now are assuming properties for this new product that is distributivity and so on so that we can expand these terms and we get three terms the first one is x the second one is this ah but a squared divided by the norm of a squared the norm of a squared is a dot a and you see here that this is the new product squared so this equals one and this term can these terms cancel and then we have one term left over well by what i just said you can see that the inverse of a vector a is a divided by a squared because if i multiply this by a i get 1 again and so this is the inverse and we are left with a much more simple formula than this one namely a sandwiching of the thing we want to reflect between the reflector and its inverse with a minus sign otherwise we get this reflection now why is this handy we want to do more than one reflection because we want to make all the orthogonal transformations and if we start with this one and we put it inside another reflection then things get to organize the minus sign cancels here i have ba and here i have ba inverse really because the inverse of a product is the product of the inverses in opposite order you know if i put on my socks and then my shoes then i have to first take my shoes off and then take my socks off to get back to the normal situation okay but we also know from geometry that if we do a double reflection of x in two planes and this is the first reflection in the blue plane the second reflection in the r in the green plane then what i've effectively done is i've ended up at the place that is achieved by rotating this vector in the common line of the planes and if a is the normal vector of the blue plane and b of the green plane and they make an angle phi over 2 then i will end up at an angle phi so i have a rotation operator it is again sandbridging and now without the minus sign because they got cancelled so an operator in its inverse and it's a rotation in 3d around the line through the origin therefore it must be isomorphic to a quaternion and indeed it is but quaternions we're used to are complex things and this is all completely real so that's nice in our real vector algebra there is a geometric product which when applied in a sandwiching operator gives us the same as having done weird things like complex numbers for rotations so that's the principle we are going to do multiple reflections but not just in the origin so now we have to find out if this is the general principle by which geometric algebra can generate reflections using vectors then we need to be in a space in which our vectors what in our real world are planes two planes two parallel planes translation two intersecting planes and rotation and that'll be handy because as we see we will have a really nice representation of orthogonal transformations in linear algebra you wouldn't even think of this because orthogonal transformations are a difficult kind of matrix for us they will be a simple kind of operator now you already know how to represent a plane as a vector a homogeneous coordinate says something like that you take the normal vector of the plane and you slap on something about the distance to the origin and this is made in saturated if you take the dot product of this with the representation homogeneous representation of a point you retrieve the plane equation x dot n equals the distance at the origin we put a plus sign here too because we can choose our units so it's three euclidean dimensions to indicate the normal vector and one extra dimension here with basis vector e that does the embedding so we're in a four dimensional space and we therefore expect that uh aquatic to create object on a that four reflections will generate all our orthogonal transformations to make reflections we need a metric and we need to specify therefore what it is to take dot products in the euclidean part of the space well that's easy e one dot e one is one e one or d two equals 0 etc but we should also say what the extra dimension does and we make it orthogonal so e dot e 1 e 2 and e 3 are all zero it's really an orthogonal part of space um and we have to do one little special thing we don't want to see epsilon squared when we start squaring planes because we don't know how to interpret those so epsilon dot epsilon or e dot e i should say should be zero and that's it so we have a metric space with this metric it's three euclidean dimensions and one dimension that squares to zero now let's see whether it works so we are going to do a double reflection in two parallel planes with normal vector n normal vector n and the second one is at a distance delta over two because we know that as we reflect we double the distance right so just like when we reflected into planes we doubled the angle and if we perform this geometric product then we see we get n squared that's here we get e and delta over 2. we know that let's say that we have normalized these planes to uh to be planes with unit normal vector then n squared equals 1 and then this is precisely the translation vector over which we would like to move and because so this is already our translation motor translation operator we can also write it in exponential form because if we see this as a series we can develop this as one plus this thing plus this thing squared but this thing squared is zero because e squared is zero and i have a minus sign here because i swap those two terms that's the conventional way of representing it same thing with rotation at the origin i take a plane in direction normal vector in direction e1 i take a plane with a tilted orientation over phi over two i perform the product and i get e1 e1 times the cosine e2e1 times the sine i know that e1 e1 equals 1 so this is a scalar and i don't know how to simplify this so i just give it another name i call it i it's sort of related to that plane that i'm rotating in this basis was chosen completely adapted to the plane right it'll be in general situation normally but there's something special about this i uh it was e1 e2 if i multiply this out then i see that i can because the dot product between e2 and e1 equals zero e1 e2 and e2e1 have a minus sign relative to each other and then i can put the brackets back in i get e1 squared e2 squared they're both one so because of that minus sign i have i squared equals minus and that means that i can also make an exponential notation for this thing because as i develop this exponential in a taylor series i got terms with i terms with i squared which gives me a minus sign i to the third which gives me i times the minus sign all precisely right to allow me to write this expression just as you're used to in complex numbers now i've chosen a special basis normally that plane would sit a skew from the regular world coordinates that i want to introduce so e1 and e2 will be linear combinations of e1 world e2 world e3 world that means that i will get mixed terms from in world coordinates e2 e3 e3 e1 and e1 e2 so this really is a vector with components on this basis the basis elements each square to minus one in the way that i've just demonstrated and so what i have here is a unit quaternion but i've made it as the multiplication of two real planes the operators we have should be applied in a sandwiching manner so an object a will be translated or rotated by a motor m in this manner if we now have a composite object that consists of a and b then something really nice happens because of that sandwiching form here we see the composition of the moved a with the move b with your composition operator it's a geometric product and because this and this cancel due to the associativity of the geometric product we can actually move the motion completely outside so we could also have made the composite object a b first and then moved it now and this is really the crucial reason of why we have the motion representation we have we have a two-sided representation that has this property as long as we make compositions out of geometric products or linear combinations of geometric products we will always have that things transform in a covariant manner whether we move them first and then compose them or compose them and then move them it makes no difference that means we don't need a method for every kind of object that we construct lines transform in exactly the same way at the points and in the proper way you know if you would have made the line out of points and move those you would have gotten exactly the same motion as moving the line directly that means that there's no need for a lot of code because all things move the same and that prevents a lot of bugs and in fact this structure is very suitable to put in the compiler and have it generate coordinate code for the things that you want to move you will know that it is correct okay so those are the motions now let's go to the primitives we've used the symmetric part of the geometric product already when we defined the geometric product in terms of the dot product but that leaves the part left over and that's the outer product or the wage product which we denote like this and it's the anti-symmetric part of the geometric product and it turns out that this is precisely the right way of spanning things in the regular geometric algebra around the origin but that in our representation of planes it's actually an intersection operator so the line common to two planes will be the wedge product of the two and it may be an ideal line if the planes are parallel then we get a line at the infinity which will be very useful certainly not exceptional it's just an element of the algebra and a point common to three planes should be the wedge product the intersection of three of these planes so let's make a point we take a plane with normal in e1 direction at the location x1 normally 2 at location x2 normal e3 location x3 that should be the plane at the location the the point at the location x1 x2 x3 when we multiply it out we get something like this and we indeed see something that has component one we see something that is a multiplication of x1 x2 and x3 but we see that the basis on which these are expressed are these tri vectors as we call them composed of these different things for an implementation if this is still homogeneous coordinates but the integration with the rest of the system enforces that we should see it as a trivet on a tri-vector basis so that means that we're getting elements that reside on a basis of scalar of the vectors those are the planes we have made lines through the origin you know when we made the quaternion we saw that we get things like e1 e2 e2 e3 e3 e1 so those are euclidean bi vectors we also get components that involve that extra dimension e which is indicated by a knot here to indicate the pattern we have just seen that we made points on a tri-vector basis and we can go one further we can also have what's called a pseudoscaler which is basically the volume element of the space all of them involved and what we've also done is we've made a rotation operator as a quaternion that was in this space we've made a translator which involves these elements and if we make a general motion we would multiply this with that and we get any element out of this total group and that is as we'll see a dual quaternion so for a three-dimensional space we have a 16-dimensional basis but we're using everything we even recognize everything we even see some things that overlap in usage like line elements are part of making a rotation operator they're all related by the geometric product sometimes specialized to the dot and the wedge for intersection and projection but they're all consistent and if we multiply out all the elements we get multiplication table but what i want to emphasize is that it seems like a lot it seems like we're introducing way too much but since all these elements are already things that you're using in your code this structure is there it's just hidden and you have to reinvent it all the time now we've made points by intersecting three planes we can also make lines by joining two points that is another operator it's called the join and it is defined in a dual way to the wedge operator it's like a boolean way of saying that if you intersect two orthogonal complements and you take the orthogonal complement then you have the connection so this is also an operator and we can show that we can write this out as a linear combination of geometric products so this is also a covariant operator i'm allowed to join the two points and make a line and move that line either by moving the two points or rejoining them or by just moving the line directly if you do it more carefully you find that these homogeneous factors that we have ignored so far are actually very useful you can make measures that clearly have an oriented nature with odd angles in the sense that if the line and the other line meet in this way intersect in this way well they don't intersect but we get something that is proportional to the distance and the angle of their direction vectors and as we move this line to the other side the distance changes if we change the angle of the relative orientation and the distance changes we can make consistent orientation measures so that you don't really need to guess the signs afterwards you can set them up consistently and the ga computes with them for you we also need to measure things if we use the regular norm that we're used to in geometric algebra on these elements that contain e the thing that squares to zero we find that we only get the norm of the euclidean part so we introduce an extra norm which can measure the size of the e part because that also turns out to be relevant in some computations and stephen the koenig really recently found a very nice property that if you take an edge loop you know see a set of connected edges by connecting points consistently from one to the other by the join operator and getting weighted lines then the contour length of this edge loop is the sum of the line lengths not surprisingly and the area is the infinity norm of the summed line objects some line by vectors and this even generalizes to 3d where you have a facet mesh where you connect your points in a consistent oriented manner in facets then the area is the sum of the sizes of these elements of these area of these plane elements and the volume is the infinity norm of the summed elements so this is very nice and we're discovering more simple operators to use in this manner now we have all these elements uh for instance here we have a red plane and a red line that's not so clear to see but that's here all these other elements are completely defined in a coordinate freeway in terms of these basic things l and p this bi vector and vector in the geometric algebra so for instance in the intersection point that's somewhere here is the wedge product of the two now this is a vector that's a bi vectors together tri-vector must be a point it turns out to be intersection point if i want the perpendicular line to display through that intersection point then that is l hp divided by p if i want the projection of the line onto the plane that is this operator p dot l divided by l if i want the projection of the plane onto the line which is the blue line we would normally not call that a projection then it's just l dot p divided by so there's a straightforward structure which can be learned and understood and we've made a cheat sheets of the operations we're beginning to find that work really nicely to fill needs all right that was geometry now let's see how far we can get into physics first we look at the line parameterization that we have with bivectors may recall that we had this in the situation where we looked at the origin and we encountered our first quaternion and we had this as e1 e2 in a general geometric algebra but in pga we would recognize e1 e2 as e1 way g2 as the mean of two planes so this is really the line through the origin so we have a rotation through the origin characterized by its line if we now use covariance that means we go to an arbitrary line that is a distance b from the line l0 then we can apply that operator to that translation to the operator because of covariance we can lift it into the exponential we see that we get the line l and so we have a full parameterization of an arbitrary rotation around an arbitrary line just as the exponential of that line then going to general motors we need an extra translation because it's not just rotations we then get the exponential of a general bi vector as we would call it algebraically geometrically we would more see that like a screw and then we have a dual quaternion but again it is generated by in principle four reflections or two of these operators that we've seen here before and they all are the form of the exponential of the invariants so here the line is invariant the exponential of the line is the rotation operator if a direction is invariant the exponential of the invariant is the translation operator and these are precisely the constructions we used in our demo so we had two lines in 2d they wedged together into a point a point is the invariant and the exponential of the invariant is then the motion operator that we can use on the tricycle to move it around in the orbits no exceptions it's just done so motion is characterized by bi-vectors as they're exponential um and bi vectors are therefore important to describe it now when we then go to the classical physics description of say kinematics then we need the momentum the momentum is a point moving with a certain velocity and if we make the bi vector the line of the point with its velocity weighted with the mass add those over all the points of our object we get the total momentum now when we implement that in pga and we have it moving with a bivector b then we find several components one is how the centroid translates and the other one is the internal rotation a purely euclidean thing this is the euclidean part of that and this is the other part and it's weighted by this i which is now the pseudoscaler of pta so it's e naught one two three a thing that squares to zero looking at dynamics we need to talk about forces forces are attached in a certain direction to a point if we write that out we find again a natural split into something that is a force and a torque then moving to the euler equations and the newton equations you know for translation and for rotation we actually find those as two components of one bi-vector equation which just says that the derivative of momentum is the force and we are now working on integrating those simultaneously in one go by rk-4 methods and you may expect a write-up about that fairly soon so that's the structure of geometric algebra pga in a nutshell we have a four dimensional space with a 16 dimensional basis these basis elements are all recognizable as either geometrical objects or as things that's composed to operators in that are isomorphic to efficient representations we already knew for these things but they're now completely in a real manner integrated with the primitives that they either work on are defined by and stephen the kerning found that you can arrange these things in a somewhat strange fashion where the things that you need for something are always in groups of four closely together in your memory this makes for really efficient gpu implementation of this structure it's been being executed by jeremy young in his program klein and it's he played some more tricks stephen he changed the basis vectors slightly to prevent all sorts of extraneous signs that might not be the same for each group of four and great additional efficiency out of that now if you're already somewhat familiar with geometric algebra you may know that we've been talking about conformal geometric algebra a lot that is a geometric algebra that doesn't have planes as vectors but spheres as vectors and if you do multiple reflections in spheres then you get something that is conformal geometry but planes are also spheres there's spheres through infinity and therefore euclidean geometry is included in this and not too long ago that was the only way we knew how to do it now the restriction to just using planes is a nice subset of cga so it is completely compatible and if you have a problem in which you do not require round elements like the spheres and the tangents and the circles and so on that are in cga then you can use pga for your application if you then suddenly have need of a sphere you just add one extra dimension the extra dimension that is relative to cga from pga and you have those as well and we have some applications in robotics for that this relationship between pga and cga is nicely being written up by chris doran and i give the link here okay maybe you want to know more um i wrote a tutorial pga for cs for computer science with especially computer science audience in mind and some demos we have the cheat sheets and we put presentations from various conferences that are worth seeing again all on this website bivector.net you can play in javascript with pga and many of the other geometric algebras and get a feeling yourself for how short the code is and how surprisingly efficiently it works thank you