Friday, April 6, 2012

Angle Between Intersecting Planes

In a practical sense, there are a few different ways that you can go about determining the angle between two planes.  But before we get to specific methods and tools, let's try to lay the ground work  let's try to understand what we really mean by an angle between two planes.

Definition:  Angle Between Intersecting Planes

Two non-parallel planes intersect at a line.  The angle between intersecting planes is measured perpendicularly to this line of intersection.  In more detail, we can say, for planes P1 and P2 with line of intersection L:
  1. The vertex of the angle between the planes is any given point on the line of intersection (that is, just pick one).
  2. One arm of the angle is coincident with P1 and the other arm is coincident with P2.
  3. Both arms of the angle are perpendicular to L.
Discussion of Definition

We note in passing that there is no need for the two arms of the angle to have any certain length nor is there a need for the lengths to be constrained by any relationship concerning relative lengths  they may differ.  Items 1 and 2 (above) should come as no surprise, but not everyone will immediately understand item 3.  Indeed, if you consider only the two trivial cases  (i) two coincident planes (technically excluded from discussion here) or (ii) perpendicular planes  item 3 appears superfluous.  To see that it is necessary, consider  planes which meet at a 45° angle (as an example).  Imagine putting a 45-45-90 triangle in place so that two of its edges are in contact with a plane each.  Suppose you made one of those planes rest on the top of the triangle and allowed that plane to rotate freely around the line of intersection (L).  Now, keeping the vertex of the triangle (the vertex on the line of intersection) in place and rotating the triangle around that point, you will find that the plane which is resting on the triangle will lift up  the angle between the planes will increase.  By the time you have rotated the triangle so that the vertical edge is against the "resting" plane, the angle between the planes will be 90°.
Here is a video where I demonstrate this scenario:


Perpendicular Angle Between Planes vs. Non-perpendicular Angle Between Planes
(Criteria 1, 2, and 3) vs. (Criteria 1 and 2)

In the above scenario we used the same triangle at different angles and changed the angle between the planes.  Now let's leave the planes be and change the triangle.  Let's draw the orthographic views of planes (represented as rectangles) intersecting at angle i, according to our criteria and draw two triangles in place:  one at 90° to the line of intersection and one at θ to the line of intersection.


I've been lazy in my drawing as the angle for β that I really want is in the plane of the triangle that's on an angle and not the angle of the orthographic projection of that triangle (which is what I have just drawn).  But bear with me.  The angle of intersection is given as i and we take the dimensions L and θ to be given.  From this information we can find formulae for R, h, and β.  Here are some formulae we can take directly off of the drawing:

\[\sin \theta = \frac{R}{L}\] \[\tan \beta = \frac{h}{L}\] \[\tan i = \frac{h}{R}\]

Manipulating these equations gives

\[R = L \sin \theta \] \[h = L \sin \theta \tan i \] \[\tan \beta = \sin \theta \tan i \]

Observe that β  i.

Calculations Without a Scientific Calculator

You can obtain a rise and run for β without a scientific calculator by using a little measurement and geometry.  Suppose you are able to obtain a rise and run for each of angles i and θ.  You might do this by measuring either a shop drawing or a physical situation.  Choose the run arbitrarily and measure the rise which corresponds to it.  The hypotenuse (for θ) can either be measured or calculated using the Pythagorean theorem.  The formulae for the rise and run of β are simply

rise of β = (rise of θ) × (rise of i)

run of β = (hypotenuse of θ) × (run of i)

Saturday, March 31, 2012

Defining a User Coordinate System

In AutoCAD you have User Coordinate Systems (UCS) and the World Coordinate System (WCS). The WCS is the system which your model information is stored in. The UCS is what you draw in. When you first begin a drawing, the current UCS is the same as the WCS.

First, let's familiarize ourselves with how the WCS looks. I've drawn a pyramid in AutoCAD and shown it from several different angles:





The WCS is displayed in all four pictures (and, naturally, is exactly the same in each picture).  One of the techniques for dealing with 3D when you are "doing the math" (perhaps in your own code) and want to make sure the math (or code) you're writing matches the geometry you're working with is known as "the right hand rule."  (We're going to use this for some math coming up in a few paragraphs.) For each picture, look at the WCS and orient your right hand (which should be open face initially) so that your fingers are pointing in the direction of the x-axis and your palm is open toward the direction of the y-axis.  (Don't think about your arm, just your hand. The arm will just go wherever it needs to go, though it may look funny to onlookers.)  Curl your fingers into your palm and stick your thumb out. If you've done it correctly, your thumb should be pointing in the direction of the z-axis.  If you can do this with all four pictures, you're well on your way to being able to understand the rest of the math in this post.

In math notation, we might write the relationship above as something like z = x × y.  The relationships which apply to the WCS (and all UCSs) are

  • x = y × z
  • y = z × x
  • z = x × y

Two of the best ways to change the UCS in AutoCAD are

  • Rotate around a chosen axis:  this is one of the more intuitive ways to adjust the current UCS when you don't have "pickable points" in the plane you want the new xy-plane to be in.
  • Pick three points:  this is a way to put the current UCS in line with a face or in line with somewhere you want to put a face.  Lining up with an existing face of a 3D drawing only requires to you to pick three of the corners of the face.  Depending on the order you pick those points in, you will get different directions for the axes, but the xy-plane of the new UCS will coincide with the three points you pick.
So, how do you define a UCS based on a three points?  Three points define a plane, but that isn't enough to define a UCS.  You need to know where the origin is, which direction the new x-axis is to go, which way the new y-axis is to go, and you also need to pick a direction for the z-axis that is consistent with the way the z-axis relates to the xy-plane in the WCS.  Note that if we just solved the equation of a plane (Ax + By + Cz + D = 0) for the three points we start with, we would not know which way the z-axis goes.  Is it above the plane? below? and why?

To obtain consistent definitions for these axes we use the cross product (which I tried to sneak in earlier:  z = x × y).  Here's the definition of the cross product:

        a × b = (a2 b3 - b2 a3, a3 b1 - b3 a1, a1 b2 - a2 b1),

where a = (a1, a2, a3) and b = (b1, b2, b3). Think about these vectors as direction vectors as opposed to points.  One of the important properties of the cross product is that it produces a vector which is perpendicular to both of the vectors that form the product. That is, a × b is perpendicular to a and to b. To determine which way a × b is pointing, point your fingers (using your right hand) in the direction of a with your palm facing in a direction such that you can curl your fingers toward b.  Your thumb is pointing in the direction of a × b.

Here's how a UCS is defined based on three points. The first point (P0) indicates the position of the new origin. The second point (P1) defines the direction from the first point the x-axis goes in. So, the direction of the new x-axis is x = P1 − P0. The third point (P2) determines not only the plane, but also which side of the x-axis the y-axis will go (within the defined plane). At the same time, we are also defining, albeit indirectly, the direction of the z-axis by means of the right hand rule. We define our UCS accordingly:
  • Origin: P0
  • x-axis: x = P1 − P0
  • z-axis: z = x × (P2 − P0)
  • y-axis: y = z × x
Notice the z-axis is defined before the y-axis as well as the order of factors which produce y. Also observe that the order in which you select the points affects which way your axes are oriented, even though it will not affect what plane is the xy-plane of your new UCS. Given any three points, there are three possible origins - so I need to specify which is the origin. With that selected, I have two options for which will indicate the direction of my x-axis. So, there are six possible UCSs given the same three points in different orders.

Depending on the use of your UCS (if you are programming one in your own software) you may wish to normalize the direction vectors. Normalized vectors (which have length equal to unity) have some useful properties which may save some computation in later calculations. Normalizing is straightforward and simply requires you to divide the value of each component of the vector by the current length of the vector.

Saturday, February 25, 2012

Best Fit Circle: find the center using Excel

Finding the center of a best fit circle depends on minimizing the same function we were concerned with in finding the radius, except that we are going to view it as a function of the center instead of as a function of the radius:

where a and b are the x- and y-coordinates, respectively, of the center and r is given by
What we have here is a function of two variables. It looks like three variables until you realize that r is calculated in terms of the other two. So, we can do a three dimensional plot and see what the scoop is. I used Maxima to do this and obtained a very good view of the surface near the best fit center of the points I have been using in all of my investigations of this problem. Here is the 3D plot of SSE(a,b):

What we are most encouraged to see in this graph is that it looks very smooth and it looks like there is exactly one point that is the lowest point. This lowest point is where the SSE function is minimized and constitutes the best center of the circle. (It might be that there are a few local minima somewhat close together that we could see if we zoomed up really tight to the bottom and we are probably happy with any of these as the "answer". Welcome to numerics.)

These formulae can be used in Excel. Designate two cells for each of the values a and b. You don't know what these are, but start with some guesses for these. You will reference these guesses in your Excel formulae. Put your points in consecutive rows after the pattern (x, y, se, R) where se references the x and y for that line as well as the values for a, b, and r. R will only reference x, y for that line and a, b from above. r above is the average of all the R values in the rows (don't include the 1/n in the R)--you may want to create a cell to contain this average and reference it in your se columns. Use absolute references for a, b, and r (if you have a cell for it) so you can copy and paste the formula easily. Make a sum formula at the bottom of your SSE column and it represents your SSE function as above. You want to use the Excel solver now. The SSE cell is the cell you tell it to minimize and the a and b cells you designate as the cells to be modified. The solver will tweak with the a and b values in an attempt to make SSE as small as possible. (The instructions about r and R might seem circular until you actually implement them. Follow through to the end and you'll see it really isn't circular.) Don't try too hard to follow the instructions--try to do the likely intent (as always).

For a Maxima approach see here.

Monday, February 6, 2012

The Polygon Worksheet

I have put together a worksheet in Excel to demonstrate a technique for calculating the vertices of a regular polygon.  You can download the worksheet here.

Overview

It demonstrates an application of vector rotation (and, by the way, complex number multiplication makes a good way to remember how to do vector rotation).  The motivating principle is to demonstrate how you can reduce calls to trigonometric functions while calculating vertices of a regular polygon.  Note that the spreadsheet does not give such (time) savings because it uses both methods: (1) direct calls to trig functions for each vertex in columns B and C and (2) vector rotation in columns E and F.  Also, the error of the vector rotation results is given in columns H and I.  You will notice that the error becomes larger as the vertex number increases.  The error occurs in the vector rotation method, not the trig function method (which is our control, since we believe these results represent the best approximation we can achieve using the floating point precision available in Excel).

The only values you should need to adjust are the number of sides cell, the radius cell, and the initial angle cell.  Everything else follows from these values.  A few notes are in order:
  1. The polygon that is determined, is inscribed inside of the circle with the given radius, and centered around the origin. 
  2. All angles are specified and calculated in radian measure.
  3. If fewer than 50 vertices are needed, the rows will repeat.
  4. If you wanted the center somewhere other than the origin, just apply a translation to the points.  So, if you want center (a, b), then each point P = (x, y), becomes P' = (x + a, y + b).  If you've made it this far into this post, you are probably more than capable.
Why Try to Reduce Calls to Trig Functions?

As the user of calculation software and scientific calculators, you may be wondering why we would want to reduce our use of trig functions.  The main reason would be for time-critical applications.  This particular spreadsheet is obviously not one of those, but it demonstrates the accuracy of the method.  So, how do I know that trigonometric functions take a long time to perform.  First of all, let me debunk some myths:
  1. Computers and calculators have tables of numbers in them that they use to determine values of "complicated" functions.  So wrong!  Computers and calculators compute/calculate these values!
  2. The calculations that computers and calculators do are just based on tables and interpolation.  This is rarely true.  Sometimes a programmer will use a method like this when he knows that the precision needs of his application will be met this way and will be much faster.  But it is actually more programming effort and would only be done when you didn't need very accurate results but you needed them really, really fast.
  3. Computers and calculators have hardware circuits that make the complex calculations almost as fast as the simple ones.  Not so.  Fact:  Multiplication is done with the fmult instruction.  Addition is done with the fadd instruction.  fmult takes way, way more computer cycles to compute than fadd.  fsincos takes way, way longer than fmult.  It is an unavoidable consequence of the nature of these calculations.
Recall that multiplication is dependant on addition.  Several additions actually.  So it should come as no surprise that by the time you put together a circuit that does the job that several additions and small multiplications can do, you end up with a circuit that takes longer to complete than one that just does addition.

The same is true of trig functions.  You know how you calculate them?  Here are the magic formulas:


The dots (...) indicate that it goes on and on and on and on and on, until the changes are small enough that you aren't changing the digits that mater to you. You may need more terms than are listed to get the accuracy you need. These formulas above are not written in the most computationally efficient manner, but you can rest assured that any circuitry or coding or combination thereof, which computes the values of the above functions will take longer than 4 multiplications and 2 addition/subtractions will (those are the operations involved in doing vector rotation).
Enter Vector Rotation

The key to how this method helps us, is that we can use a few trig function calls at the beginning and reuse the results several times over. We compute the direction vector for the central angle, θ, that each side covers. For a hexagon, this angle is 360°/6 = 60°, or π/3.
We need to determine the value of the central angle using direct calls to trig functions, but we will be able to reuse those values.  We also call trig functions to get our initial position.  From there we use our vector rotation method.  So, E14 and F14 use the values in E11 and F11 (renamed mx and my) and those of E13 and F13 to find their values:

E14 =E13*mx-F13*my
F14 =E13*my+F13*mx

Cells in rows below 14 do the same thing.  They still reference mx and my, which are the cosine and sine of our central angle, θ, respectively, and they reference the cells directly above them.  The concept for each row is:  rotate the point in the row above me around the origin by the central angle and tell me the x and y values for that rotated point.

Proof of the Formulas

Suppose we are given a point (x, y) and denote the distance from the origin as r and the angle it makes with the positive axis as i.  We wish to rotate (x, y) counter-clockwise about the origin by angle θ.  The diagram below illustrates:
Our new point (x', y') can be simplified:

                  
                      
                      
                 
                     
                     
If you look back at the formulas for E14 and F14, you may observe that mx and my correspond with sin θ and cos θ and E13 and F13 correspond with x and y.  That's all there is to it. 

Relationship with Complex Number Multiplication

If we represent the new point (x', y') as x' + y'j (where j is used to indicate the square root of -1), we can say:

This will produce an equivalent result to what we have above.

Saturday, January 14, 2012

Best Fit Circle: find radius given center

Given a center (a, b) of a circle, we want to find the best fit radius, r, to a set of given points Pi = (xi , yi), for i = 1 to n. By addressing the relationship between best radius for a given set of points, we will be able to make r dependent on a and b, rather than being something we guess at independently. (This is ground work for another post.)

In other words, we have a set of points and a candidate center for what we think is a good approximation to a circle (or a circular arc). We want to find the best radius for that center with those points and have a means of quantifying how good of a fit we have. If we can find this, we can decide which of a set of candidate centers is best.

The usual thing to do when looking for a best fit shape is to minimize the sum of squares of the errors. If we knew the radius we might calculate the sum of squares of the errors as


So, for each point, we find out how far it is from the candidate center, find the difference between that distance and the radius to get the error, and then square the error. And then we add all of those squared errors together to get SSE(r). We want to minimize SSE(r), which is a simple problem in differential calculus. When SSE’(r) = 0, SSE(r) is at a minimum, maximum, or perhaps a point of inflection. Then





Observe that the best fit radius is just the average distance from the candidate center to each point. Seems sensible. Also, observe that
and so, by the second derivative test, SSE(r) is at a minimum.

To see it all put together into a practical solution, see Best Fit Circle: find the center using Excel.

Friday, January 6, 2012

Field Measurement of Circular Arcs

There are a number of simple cases for field measurement of a circle. A full circle can be measured across its diameter or around its circumference to determine its area. (To use the circumference to determine the area, first use the formula C = 2πr to determine the radius and then use the radius to determine the area which gives A = C2/(4π).) Half circles and quarter circles are similarly straightforward.

But suppose you want to measure a circular arc which has an unknown central angle, θ, such that 0 < θ < 180°. There are three things that are normally easy to measure for such circles and they only require a tape measure:

  1. 1. chord distance (also called the run of the arc – denoted u)
  2. 2. rise of the arc (i)
  3. 3. arc length (s) – arc length measurements are easy for existing items if there is a (normally vertical) surface to hold the tape against; not as easy for proposed items unless approximation is acceptable

The following diagram illustrates:

There are number of reasons why you might be interested in the rise, run approach to defining a circular arc. If you are laying out an arc, it may be infeasible to run a tape measure around the center point of the arc – perhaps due to the size or interfering objects. Also, measuring unusual angles requires specialized equipment (such as a transit) which requires set up time. (If you are using a total station this whole discussion is moot, so we are assuming we’re trying to do something without it.) The rise/run approach is also the approach used in some blue prints for defining curved walls. Carpenters and linear measurements are on good working terms, and we like to leverage this where we can.

But whatever the reason you might want to use this approach, you can calculate both ways. First of all, let’s list all of the relevant equations:


 (1)

 (2)
                     
 (3)
         
(4)



Now, let’s suppose we know u and i. We determine the remaining variables from these.
Manipulating equation (4):



(5)

So, we need a and θ (or at least cos θ) to find r, which we do by manipulating equations (1) and (3).
Manipulating equation (3):
(6)
Manipulating equation (1):

(7)
Substitute (6) and (7) into equation (5):






(8)


Now that we have r we can use that number in equation (7) to find θ, etc. I omit all of the rest of the gory details, but here are the final results for solving the unknowns in terms of different combinations of known values (note that θ is understood to be in radian measure and the reader is left to discern what order to calculate the different values in).


I have not addressed how to solve the cases where s and either u, i, or a are given as they are more difficult.  Perhaps a future post will address these cases.

Saturday, December 31, 2011

Field Measurement of Quadrilaterals Using Only a Tape Measure

Trapezoids
A common nonrectangular area that needs to me measured in the field is a trapezoid. Sometimes shapes are approximated by a trapezoid, if a rectangle is not considered a reasonable approximation or not easily visualized. The well-known formula for the area of a trapezoid is

                        clip_image002

where b1 and b2 are the lengths of the two parallel sides and h is the distance between them (measured perpendicularly to the parallel sides, of course). It is noteworthy that the (b1+b2)/2 can be interpreted as the average of the lengths of the parallel sides. This is the length of a line referred to as the median of the trapezoid. In the field, it is often more convenient/faster to measure the median and the height (h, also called the altitude) than to measure both bases and the height. If the location of the median can be “eye-balled” with sufficient precision for the purposes being met, this may increase productivity by decreasing the number of measurements needed. For more information on the trapezoid and the median, see a description here.
General Case Quadrilateral
There are (at least) two approaches to measuring the area quadrilaterals when nothing is known about the internal angles or “parallelness” of nonadjacent lines. The following diagram displays both methods.

We can determine the area of this shape by making 5 measurements and produce an exact picture, or by taking 3 measurements and getting an approximate area. Both are “area by triangles” methods.
Area by Triangles – Heron Method
If we measure the sides AB, BC, CD, DA, and one of the diagonals (either AC or BD), we can use Heron’s theorem to determine the area of the triangles on either side of the chosen diagonal. We can use the law of cosines and law of sines to determine the angles. Alternatively, a basic AutoCAD drawing using temporary construction-line circles would allow you to determine these results without manual calculations. This would give us a full description of the area. This method also generalizes well to an arbitrary number of sides (but requiring a lot of measurements--don't use this method).
Area by Triangles – By Altitude
On the other hand, we could measure the length of a diagonal (say, AC), leave the measuring tape in place (or replace with a string line) and use a second tape measure to measure from the other vertices (B and D) to the diagonal. As long as the accuracy/precision of the result desired is not too high, we can “eye-ball” perpendicular to the diagonal to get the heights of the triangles to an acceptable accuracy. The area calculation is obvious from there. This method involves fewer measurements, minimal extra equipment, and is quicker to calculate, which may make it more efficient for a quadrilateral. It does not generalize to an n-sided figure as easily as the previous method, but could be done with additional string lines used simultaneously.