SimpleVec wrap - Embedded.com

SimpleVec wrap

Jack wraps up his discussion of SimpleVec.cpp and demonstrates some useful things you can do with vectors.

In the past few columns, I've embarked on a project to give you a set of C++ files that support mathematical operations on vectors and matrices. This month, I'll wrap up my first offering by collecting the functions into a proper linkable file.

Let's be clear on one thing at the outset: the kinds of vectors we're talking about here are the mathematical entities that physicists and engineers might recognize as useful to represent such practical, real-world entities as velocity, force, momentum, magnetic field, and so forth. Specifically, we're talking about entities that behave like vectors under mathematical operations such as addition, subtraction, and multiplication. These rules, after all, were defined to mimic the way the real, physical properties behave. We are not talking about mere arrays of things, even of numbers, as they are supported by the STL library class, vector . Despite its unfortunately chosen name, that class doesn't support any kind of vector arithmetic.

My first tiny step into the development of a useful library has been a file called SimpleVec.cpp (posted on www.embedded.com/code/2007code.htm). It contains a collection of the most primitive vector math operations with no frills. It's written in very much the style you might find in an old Fortran library–in fact, my old Fortran library, which I used to carry from job to job like a master mechanic carries his toolbox (maybe you've noticed the column's title?). Although written in C++, SimpleVec.cpp would feel at home in a C environment, and that fact is deliberate. I chose to start with functions that could easily be converted to other languages, including C and even assembly language (been there, done that).

For the sake of simplicity and performance, the functions in SimpleVec.cpp include almost no error checking. There is no protection from array bounds violations, zero or negative size parameters, or vectors of unequal lengths. This, again, is very much in the tradition of such routines. Assume the users know what they're doing and write the software so that the functions are not called with ridiculous arguments. If this bothers you, it's easy enough to add asserts, which will at least make sure you know when a function has been misused.

This month I'll wrap up SimpleVec.cpp . Then I'll move on to show you some of the neat things you can do with vectors.

More on options
In my last column, I discussed some of the e-mails I received on this topic, mostly from people who would have done things differently. I explained the disparity by noting that there are many options to consider, which affect the “-ilities” of the functions: things like readability, maintainability, ease of use, reliability, and efficiency. Choices made to improve one of these attributes tend to adversely affect others. The decisions are not clear cut, which is precisely why I tend to make different choices on different occasions, even after all these years.

Case in point: almost all vectors representing real, physical vectors are 3-vectors, meaning that they have three components like x , y , z or length, width, height. Naturally, we'd like to operate on these vectors rather efficiently, because we use them a lot. On the other hand, some useful vectors have different lengths, ranging from the two-dimensional vectors expressing motion in a plane, to the four-vectors made famous by Dr. Einstein and used in GPS calculations, to entities representing hundreds, or even thousands, of individual measurements of some quantity and used in dozens of operations such as statistical analysis or guidance and navigation. We need to be able to support operations on them, also.

Of course, if you have functions that support n -vectors, where n is some arbitrary size, you can always invoke them with n = 3. I've done this many times in the past, but physicists are a lazy lot, and it gets a little tiresome to write the parameter list like (.., .., .., 3) over and over. Aside from the bother (and the slight performance hit), there's always the risk that you might use the wrong number, leading to really obscure and hard-to-find bugs.

In my earliest versions of SimpleVec.cpp , I tried to support both kinds of vectors with two sets of functions: one for 3-vectors and one for n -vectors, where n is some arbitrary number of elements. I pointed out that the functions are, after all, quite small, so we are hardly stressing the memory capacity of a modern PC to provide both sets.

The down side is, now you need two sets of names, as well. And if the data types can be different, you need sets of functions for each type. In my old Fortran days, I can't tell you how many times I maintained functions like addv , addv3 , daddv (for double precision), daddv3 , and so forth. I didn't like it then, and I still don't. Later on, when we're dealing with a C++ class for vectors, we'll be able to use name overloading, but that comes later.

Now, C++ supports a concept that's simple to implement and seems to solve the problems permanently: the default parameter. We can define a size parameter and let it default to three. If we want to do operations on an n -vector, we can do so, simply by passing the size in as a parameter. On the other hand, we can leave it out entirely, and it will default to three. That's the choice I've made in the current version of SimpleVec.cpp .

If, on the other hand, your own choice of “-ilities” happens to stress performance, you should be aware that using a default parameter can make a significant hit on the ability of the compiler to optimize. All the vector operations do things in loops, such as:

for(i = 0; i < n;="" ++i)="" do="" something="" with="" a[i];="">   

If n is an input variable, the compiler has no choice but to use indexed addressing on each access to a[ ] . On the other hand, if the compiler knows the value of n , and knows that it's small, the compiler can unroll the loop into successive statements, like:

    Do something with a[0];    Do something with a[1];    Do something with a[2];   

Having done that, the compiler can then replace references to the indexed array to direct access to the known addresses of the elements of a[ ] .In the many e-mails I've received on the subject of a vector library, I particularly liked the advice of one reader, who I'll paraphrase here: “You shouldn't give a what other people think. You should do what you think is right and let them decide whether to follow suit.”Words to live by, and he's exactly right. After all, it's not my job to write your software for you; it's not my job to give you a file so universal that you can copy it verbatim, secure in the knowledge that it will meet your needs to a T. As I hope I've made clear by now, your own weighting of the various “-ilities” can require you to choose different approaches. I've done my job if I've explained the options to you sufficiently so that you can make wise choices.'Nuff said.New functions
In wrapping up SimpleVec.cpp , I've added a few more functions. One glaring hole in my previous file has been the lack of vector I/O. We could do arithmetic with vectors, but we couldn't input them or show the results. In functions vGet( ) and vPut( ) , I've provided this functionality. The functions are shown in Listing 1, and it's clear that they are very primitive. In the past, I've written much more powerful versions, but these functions are intimately related to the whole user interface. For that reason, you'll probably want to tailor the functions to your application.


Before leaving the I/O routines, I should mention that the function vGet( ) provides no protection whatever for badly formed inputs. Mostly, that behavior comes right from the iostream class itself. Give each input a valid numeric value, and everything works fine. Give it some ridiculous value like 3cf.12q, and strange things will happen, up to and including going into an infinite loop. For production use, a proper input routine should be able to handle any input one throws at it. If your application involves accepting vector inputs from the console, you will definitely want to beef up the error checking considerably.Comparison operators
The other functions I've added allow us to compare two vectors (of the same length, please) for magnitude. This capability requires a little explanation. Since a vector has more than one element, it's certainly possible that one element of one vector is greater than the same element of another vector, while a different element may be less than its counterpart. In such a situation, what could it mean to say:

The only reasonable answer is to compare, not the individual elements of the vectors, but the magnitudes, or absolute values. Notationally, we write:

where this absolute value (norm) is computed from:

The one exception is the comparison for equality and its converse, inequality. Here, we say that two vectors are not equal unless every element is equal, term by term.

There are six comparison operators all told, implemented as three tests and their Boolean complements. Two examples are shown in Listing 2.

The code in the case of vIsLess( ) bears comment. Even though the size parameter defaults to three, as you can see, we must still pass n to nested functions, since they can't otherwise know which size to use.

The final form
Until now, I've shown SimpleVec.cpp as a single file, consisting of a collection of functions plus a test driver. The final step in giving you a proper file is to collect the functions up into a separate file, and create a matching header file. That process is straightforward enough. Space prohibits me from giving you the files in printed form. Instead, you can find them on the www.embedded.com/code/2007code.htm.

Solving problems with vectors
Now that we have a way to use vectors, what shall we use them for ? Well, you wanna build a flight simulator? You can use vectors to compute the relative wind. From that, you can calculate, then sum, the forces of lift, drag, thrust, and gravity. You wanna control an industrial robot? Vectors and matrices let you look at the geometry of the arm segments, and rotations of the joints, to compute and control the state of the end effector.

Designing a flight simulator is a bit much for this column. But there are still some neat things–fundamental things–I can show you. Consider the pair of vectors shown in Figure 1.

Remember the definition of a dot product. It's a scalar, whose value is:


(1)

(Here I'm using bold-faced type to represent vectors; italics to represent their scalar magnitude.)

Suppose I want to know the angle between a and b . I can find it directly from Equation 1. Solving for θ gives:


(2)

See those two terms in parentheses on the right? The term a /a represents the vector a , divided by its magnitude. That's a vector of unit length, directed along a . The term b /b is similar. Now, remember that I've defined one function that takes the dot product of vectors, and also one that converts them to unit vectors. Writing the software is a no-brainer:

    a = vUnitize(a);    b = vUnitize(b);    theta = acos(Dot(a, b));   

In geometry problems, we often need to know the component of a vector that's parallel to a second vector. That's the same thing as the projection of one vector onto the other. Mathematically, the projection of a onto b is equal to:


(3)

To get this, we base our solution, again, on Equation 1. The formula is almost exactly the same, except that we only unitize b .


(4)

This code is even simpler:

    b = vUnitize(b);    a_sub_b = Dot(a, b);   

The cross product
You think that's cool? Wait'll you see what we can do with the cross product!

Let's begin with a parallelogram defined by two vectors, as in Figure 2.

In school, you learned that the area of a parallelogram is equal to the base times the height. However, a more convenient measure for us is:


(5)

But this is precisely the definition of the vector cross product, or at least its scalar magnitude. Therefore we can write the area as, simply:


(6)

Again, the code to effect this computation is trivial:

    Cross(a, b, c);    S = vNorm(c);   

But wait, there's more. Recall that the cross product is a vector product; its result is another vector, which has both magnitude and direction. Often, the direction is useful to us as well. You might recall that the direction of the result vector is perpendicular to both argument vectors, and has a direction given by the right-hand rule. That is, for the case a × b , you should curl your fingers (right hand, please) in the direction that would rotate a into b . Your thumb now points in the direction of the cross product. One consequence of this definition is that b × a is not the same thing as a × b , but its negative; the cross product operator is not commutative.

In Figure 3, I've tried to show two vectors in three dimensions, and their cross product. Recall that three points in space define a triangle, which in turn defines a plane. As with the parallelogram, we can define the three points, and therefore the triangle and the plane, by giving the two defining vectors.

First, consider the area of the triangle. Since it's half a parallelogram, you will not be surprised to learn that the area of the triangle bounded by a and b is:


(7)

This is precisely half the length of the vector:


(8)

But wait, there's still more. Because n must be perpendicular to both defining vectors, it is also perpendicular to the plane that they define. Is there any value to knowing the orientation of a surface in space? You bet there is. You want to render it in a graphical engine and take out hidden lines? Simply look at the direction of n . If it's pointed away from your line of sight, you can't see the surface. If it's pointed toward you, its orientation relative to your line of sight and to the light sources tells you what sort of shading you should apply to it.

Remember, the vector n is a faithful representation of the area of the triangle, or rather, twice that area. Chances are good that this is the first time anyone ever told you that an area can have direction as well as magnitude. If that surprises you, you may be shocked to learn that area can also be negative. In Figure 3, I've shown the vector n pointing upwards, because that's the direction of a × b , as given by the right hand rule. But swap the positions of a and b , and n should point downwards. Voila! Negative area.

Is there any value to this fact? You bet there is, and in it lies the solution to an eminently practical problem that has kept surveyors thinking for thousands of years. Consider the polygon in Figure 4, which I'm showing in two dimensions. Because all the edges in the polygon lie in the plane of the paper, all “area” vectors lie in or out of the plane. That is, positive or negative.

Imagine that we visit the vertices of the polygon, traversing along its defining edges, in a counterclockwise fashion. Note the little blue arrows in Figure 4. At each vertex, there is a vector like v1 from the “vantage point” at O to the vertex. (We don't absolutely have to define such a vantage point, but the process is much easier to visualize if we do).

As we traverse the polygon from, say, point A to point B , each pair of vectors defines a triangular area. Because the rotation as seen from O is clockwise, that is, left to right, the area is represented by a vector directed downwards. In other words, it's a negative area, colored yellow in the figure. Similarly, as we go further around, from B to C and back to A , we're sweeping in the other direction, sweeping out positive areas, colored blue. Where the two sets of triangles overlap, the difference of the areas (and therefore the sum of the cross products) gives us areas such as the green one on the right. In short, the arithmetic sum of all the cross products gives us the area inside the polygon, which is a nice thing to know. Mathematically, the area is:


(9)

If there's any trick to the computation at all, it's remembering to close the loop by going from the last point back to the first one.

Back in my college days, I worked as a rodman for a surveying company (if you don't know what a rodman is, count your blessings). In the map room area in the back of the office, there was this marvelous machine called a planimeter, which calculated the area of a property of any shape, mechanically. What the surveyor did was to draw a plan of the property of interest, to scale. Then he simply ran the cursor of the planimeter around the boundary of the drawing, much as I've shown it in Figure 4. After this simple process, he read the area off of little dials.

At the time, I thought that the planimeter was something very close to magic. Now that I know the math, I can see what it was doing. Thanks to an ingenious set of wheels and linkages, it was summing, on the fly, a whole bunch of incremental cross products. Just as in Equation 9, the sum of all the cross products gave the area.

Neat. Also, as I said, eminently practical. As a down-to-earth (literally), practical application of vector analysis, it's hard to beat the planimeter.

Coming attractions
We now have a working set of vector math functions. They get the job done, but they ain't purty. Fortunately, by using the features of C++ more effectively, we can do wonderful things to make their use more palatable. In my next column, I'll begin to develop a true vector class, that in the end will make writing math expressions involving vectors as easy as those for scalars.

See you then.

Jack Crenshaw is a senior systems engineer at General Dynamics and the author of Math Toolkit for Real-Time Programming. He holds a PhD in physics from Auburn University. E-mail him at . For more information about Jack

Leave a Reply

This site uses Akismet to reduce spam. Learn how your comment data is processed.