Sunday, September 6, 2015

Plotting a LWPOLYLINE in Maxima

In this post, I want to show how to produce a representation of LWPOLYLINEs in Maxima that is suitable for plotting to the plot2d function. The direct usefulness of such an operation is dubious. Think of this as a worked example to help build understanding of a few mathematical and computing tools.

Recap
The LWPOLYLINE is a 2d representation of polylines that represents each vertex as having x-value, y-value, and bulge value. The bulge value indicates the path of travel from the current point to the next one. Either it will go straight (bulge = 0), veer to the left and back as a CW arc (bulge < 0), or it will veer to the right and back as a CCW arc (bulge > 0).

So, it will
  • bulge left (< 0),
  • go straight (= 0), or
  • bulge right (> 0).
More details on the mathematics of the bulge factor can be found here and here.

The Pieces
First of all, we define our problem and break its solution up into pieces. A polyline is of such a nature that is does not make a good "\(f(x)\)" type function (plotting \(y\) versus \(x\)). However, it lends itself to parametric expression well enough. So, we will have an \(x(s)\) function and a \(y(s)\) function. Ideally, I would like to specify a length along the polyline (measured from the start) and get a point \((x(s), y(s))\). (If you have a surveying background, think \(s = \) station number.)

So, I need to have parametric equations for lines and for arcs that are parameterized for arc length and initial value of the parameter. I then need to stitch these together to make both an \(x\) and a \(y\) piecewise function.

Linear

Parametric equations for lines involve a direction and a starting point. Simply normalize the direction vector to parameterize for length. This looks like
\[(x(s), y(s)) = A + \frac{(B-A)}{||B-A||}s,\]
where \(A\) is the starting point and \(B\) is the final point. We will need to adjust the initial value to account for all of the length of the polyline prior to the initial point on this segment, say, \(s_i\) which changes our equations to
\[(x(s), y(s)) = A + \frac{(B-A)}{||B-A||}(s-s_i).\]

Arc

Parametric equations for arcs use a center, radius, initial angle, direction of rotation, and total angle subtended (some variation of feasible input parameters admitted). Parameterizing for arc length isn't as obvious as it was for lines, so let's start with a basic circle:
\[x(s) = c_x + r \cos{s}\]
\[y(s) = c_y + r \sin{s},\]
where \((c_x, c_y)\) is the center and \(r\) is the radius of the circle. For our reparameterization to work we require that radian measure be used (or else require degrees and do the conversion in the function, but why make work for ourselves and introduce inefficiency?). Recall that the definition of radians is the ratio of arc length to radius (\(\theta = s/r\)). So, our parameter becomes \(s/r\) and so we have

\[x(s) = c_x + r \cos{\left(\frac{s}{r}\right)}\]
\[y(s) = c_y + r \sin{\left(\frac{s}{r}\right)}.\]

To get the right starting point for our function we need the parameter inside the trig functions to be equal to the correct initial angle when the input value of \(s_i\) (length of polyline so up to the start of the arc) is given. We need, then, the initial angle, say \(\theta_i\), plus the angle past initial that the given \(s\) corresponds to. Thus,

\[x(s) = c_x + r \cos{\left(\theta_i + \sigma \frac{s-s_i}{r}\right)}\]
\[y(s) = c_y + r \sin{\left(\theta_i + \sigma \frac{s-s_i}{r}\right)},\]
where \(\sigma\) is 1 or -1 as per the sign of the bulge factor.

Piecewise

Making a piecewise function parameterized for length requires that we keep track of the length of the polyline up to the segment we are about to address. Implementing it in Maxima also requires a bit of acquaintance with the methods available to us for doing so. It is possible to use if-then-else logic to achieve this and obtain something that works fine for plotting purposes. As an example of this, consider a really simple example:

f(x) := if x < 1 then x^2 else x^3;
plot2d(f(x),[x,-1,2]);

Fig 1. If-then-else logic in functions plots fine.
There's another, somewhat handier approach to this problem, namely, the characteristic function. The characteristic function function is a way of expressing conditions as multiplication by 1 or 0 in order to include and exclude the appropriate expression. charfun2(x,a,b) returns 1 if \(x\in [a,b)\), else 0. So, we can modify the above function to be as follows and obtain the same result. 

load(interpol);
f(x) := charfun2(x,-inf,1)*x^2 + charfun2(x,1,inf)*x^3;
plot2d(f(x),[x,-1,2]);

One nice thing about these expressions is that it is possible to create them separately and add them together later and get a relatively normal Maxima expression. (There are some hiccups if you try to differentiate or integrate them.) I can't comment on the performance of the characteristic function approach versus the if-then-else approach, but the characteristic function approach may be more manageable to produce.

The Program
There are a number of details that don't specifically relate to this problem that I will not go into detail about and other details I have given in the two previously mentioned posts on bulge factors. Here is the complete code for producing the parametric equations and plotting them followed by the result for a specific example polyline.



Circle Fitting in Maxima

In previous posts, I've considered the problem of fitting a circle to 2D data points. Fitting problems such as this require:
  1. An equation form to fit to. 
  2. A function to quantify the error, which is to be minimized.
For a circle, the form of equation we want to fit to is
\[r^2=(x-a)^2 + (y-b)^2.\]
Maxima's lsquares package removes the burden of item number 2 for us. But for interest's sake, the error function we want to minimize can be given as
\[\mathrm{r}\left( a,b\right) =\frac{\sum_{i=1}^{n}\sqrt{{\left( {x}_{i}-a\right) }^{2}+{\left( {y}_{i}-b\right) }^{2}}}{n}\]
\[\mathrm{SSE}\left( a,b\right) =\sum_{i=1}^{n}{\left( \mathrm{r}\left( a,b\right) -\sqrt{{\left( {x}_{i}-a\right) }^{2}+{\left( {y}_{i}-b\right) }^{2}}\right) }^{2}.\]
Recalling that \(\sum_{i=1}^{n}(\bar{x}-x_i)^2 = \sum_{i=1}^{n}{x_i}^2 - (\sum_{i=1}^{n}{x_i})/n\), we see that
\[\mathrm{SSE}\left( a,b\right) =\sum_{i=1}^{n}\left({\left( {x}_{i}-a\right) }^{2}+{\left( {y}_{i}-b\right) }^{2}\right)-\frac{{\left(\sum_{i=1}^{n}\sqrt{{\left({x}_{i}-a\right) }^{2}+{\left( {y}_{i}-b\right) }^{2}}\right)}^{2}}{n}.\]

lsquares_estimates()

Given a set of points in matrix form (m), we need only to use the equation of a circle and indicate which are the data point variables and which are the "solve for" variables. In our problem,

lsquares_estimates(m, [x,y], (x-a)^2 + (y-b)^2 = r^2, [a,b,r]);

will suffice. The following function can be used to check the SSE.

sse(a, b, r, pts) := sum((r-sqrt((pts[i][1]-a)^2+(pts[i][2]-b)^2))^2,i,1,length(pts));


Testing for Reasonableness

It is reasonable to ask whether least squares (or other data fitting methods) will reliably produce good results with a particular quality of data for a particular application. This is true for any problem. What if the math works fine, but the likelihood of being able to "math filter" the  likely errors in the data to reach an accurate answer is not so high? In other words, how much inaccuracy in the collected data can you tolerate and still have a reasonable expectation of getting to the right answer? I'm not going to attempt to deal with probabilities (per se) or confidence intervals, but look at a simulation approach that may be able to give you an idea of the reasonableness of the circle fit based on an (intuitively) estimated amount of random error in your data.

Simulation of Circular Arc Data

We will consider the problem of data collected, using a total station, of a circular arc, with an "as-built" scenario particularly in view. As such, auto-generated sample points should simulate the expected behaviors of a rod/prism person stepping out regular sample locations from start to end of an ostensibly circular arc. This can be implemented as random normal perturbations of a circular arc.

One of the first pieces to a realistic perturbed circle is to be able to apply randomization that follows a normal distribution given an expected value and a standard deviation. Maxima has a function that behaves this way. random_normal() can be called with an expected value and standard deviation and optionally with a number of values to generate. Here is a sample use and output of the function, pictured using the histogram function:

load(distrib)$
load(draw)$
histogram(
    random_normal(50.0, 1.0, 1000),
    nclasses = 10,
    xrange = [40,60]);

Fig. 1. Normally distributed data
The standard deviation should be chosen such that roughly two-thirds of the time (≈68%), the sample value will be within one standard deviation of the mean; that is, on the interval \((\bar{x}-\sigma, \bar{x}+\sigma).\) 

Where the radius is concerned, we have the relatively predictable variations of the prism person's actions and the relatively unpredictable variations of the ostensible circular arc, which could turn out to not be an arc at all. In the absence of any sort of meaningful solution to this latter variation, we could possibly just bump up the standard deviation. But I'm not clear on whether this is a meaningful thing to do in assessing the acceptability of the data collection precision. The error should be "small" if the thing we measured conforms well to a circle. In order to check for whether this is "small" you need to produce a standard deviation from the SSE: \[\sigma_{SSE} = \sqrt{\frac{SSE}{n-1}}\] In the scenario we're dealing with here, someone would probably draw the computed arc in a CAD program and compare it to the measured data points and in an "arm-waving" sort of way decide whether it looks good or not.

Where the point sampling distances are concerned, we will assume that the rod person can step out distances (or angle sweeps) within some percentage of that intended. We will also make the pragmatic assumption that the prism person has decided a definite number of samples per given arc and will try to make them evenly spaced. The prism person will try to subdivide the remaining space at each sample point, and two thirds of the time, they will be within some percentage of the distance (or angle sweep) they intended.

Here is the perturbedArc() function with a sample usage:


Fig. 2 - This is a perturbed circular arc. Doesn't look very perturbed does it? If we bumped up the standard deviation of the radius enough, we'd get something more clearly perturbed.
To produce the matrix form of the data I use a simple function I created (see here) called MakeMatrix(). Beyond that, it is a simple as

m: MakeMatrix(pts)$
lsquares_estimates(m, [x,y], (x-a)^2 + (y-b)^2 = r^2, [a,b,r]);
subst(%[1],sse(a,b,abs(r),pts));

Thursday, August 6, 2015

Polylines: Radius-Bulge Turnaround

The dxf representation of LWPOLYLINEs represents arcs using the end points of the arc and a bulge factor which is the tangent of a quarter of the central angle of the arc (see this previous post). We want to interchange between the radius and the bulge factor (both ways). Let's pull up the earlier diagram.

Looking at the diagram, we can produce several useful formulas. Taking the bulge factor as \(b\), observe that
\[b = \tan {\theta \over 4} = {i \over u/2} = {2 i \over u}\tag{1},\]
where we are temporarily ignoring the issue of sign.
The distance \(u\) is simply the distance from \(A\) to \(B\) and these point are always given to us in the context of a polyline whether we are given the bulge factor or the radius. We write this as 
\[u = ||B-A||.\tag{2}\]
Using the Pythagorean theorem we have
\[r^2 = {u^2 \over 4} + a^2 .\tag{3}\]
Using the substitution \(a=r-i\), we obtain
\[r^2= {u^2 \over 4} + (r-i)^2 = {u^2 \over 4} + r^2 - 2 r i + i^2\]
and so
\[2 r i = {u^2 \over 4} + i^2.\tag{4}\]

Radius/bulge direction => Bulge factor

Using equation (3), we can solve for \(a\)
\[a = \sqrt{r^2 - \frac{u^2}{4}}.\]
From this we can obtain \(i = r - a\) which allows us to calculate \(b\) according to equation 1. Additionally, set the sign of \(b\) to be positive for a CCW arc and negative for a CW arc. Alternatively stated, \(b\) is positive if the peak of the arc is a right offset from ray \(\overrightarrow{AB}\) and negative if the peak is a left offset from ray \(\overrightarrow{AB}\). In our diagram above, if point A comes before point B in the order they occur in the polyline, then we have a CW arc (P offset to the left of \(\overrightarrow{AB}\)) and so \(b\) should be negative. A negative bulge value indicates bulging to the left and a positive value indicates bulging to the right.

Bulge factor => Radius/bulge direction

We calculate \(u\) from \(A\) and \(B\) and then substitute equation 1 (solved for \(i\)) into equation 4 to get
\[2 r \frac{b u}{2} = \frac{u^2}{4} + \frac{b^2 u^2}{4}\]
\[4 r b =  u + b^2 u\]
\[r = u \left(\frac{1 + b^2}{4 b}\right).\]
Note that we would actually use the absolute value of \(b\) to obtain a positive radius. We can use the same logic as above for determining CW versus CCW (it's "iff" logic).

Locating the Center Point

Whether we are given the radius or bulge factor, we may want to find the center point. From above, we know we have \(u\), \(i\), \(b\), \(r\), and \(a\) all easily available to us if we need them.

We take our unit vector for the direction of the line segment taken from \(A\) to \(B\) as \(d = (B-A)/u\). Then the normal direction (90° right turn for offsetting) is given by
\[N = (d_y, -d_x).\]
Taking \(\sigma = {\rm signum}(b)\),
\[C = \frac{A+B}{2} - \sigma a N\]
\[P = \frac{A+B}{2} + \sigma i N\].

In my next post, I show a way to implement this in Maxima.

Wednesday, May 6, 2015

>> Imitation Forward Pipe for Common Lisp

Since learning a little F# about a year ago, I took a liking to its forward pipe operator (|>). The semantics of the forward pipe work great with F#'s sequence data type (with so-called, "lazy" computation) and functional style programming. It allows you to express a progression of consecutive steps in the order they will be done, while allowing you to directly feed the results of the last function into the input of the next function, avoiding redundant variable declarations. I suppose you could call these types of consecutive actions "transitive" combinations of operations. Imperative style programming would preserve the apparent sequence of events (meaning, the order you type the operations in, is the order the operations happen in), but it also involves declaring intermediate variables, which perhaps you will only use once.

I'm not an expert on functional style programming and don't intend, right now, to explain or justify its merits, except to say that sometimes it produces cleaner, clearer code. But while I have a current preference for Common Lisp above most other (programming) languages at the present, I really like that forward pipe operator in F#. This morning, I thought it was time to cross that bridge—try to make an imitation forward pipe in Common Lisp.

In F#, if you wanted to turn a into b into c, you might do something like

let a2c x =
    a2b x
    |> b2c

where b2c also takes one parameter. In Lisp, you end up with something more like

(defun a->c (x)
  (b->c (a->b x)))

This is fine for being concise but it doesn't preserve the (apparent) logical ordering of \(a\to b\to c\). We could do that with a change to something more imperative, such as

(defun a->c (x)
  (let ((b (a->b x)))
    (b->c b)))

and this works okay. However, I would like to write

(defun a->c (x)
  (>> 
    (a->b x)
    (b->c)))

Or,

(defun a->e (x)
  (>> :pipe-in
    (a->b x)
    (b->c)
    (c->d :pipe-in 0.0625d0)
    (d->e)))

where we want this last one to expand into

(D->E (C->D (B->C (A->B X)) 0.0625))

I've introduced a keyword so that we can dictate which parameter of the next function the previous result goes into. You can specify any keyword to use as a pipe-in parameter, but make sure you don't use a keyword of one of the functions being called in the sequence—it'll really mess with your mind. Here's the macro definition along with an error condition:



To test that the code works, we can try this at the REPL,

CL-USER> (macroexpand-1 '(>> :pipe-in
 (a->b x)
 (b->c) 
 (c->d :pipe-in 0.0625) 
 (d->e)))
(D->E (C->D (B->C (A->B X)) 0.0625))

This macro does not attempt to handle multiple return values—though such a development might lead to an interesting result. It might be as easy as introducing more keywords.

Monday, April 6, 2015

Dynamic Zoom in LispWorks' CAPI

Using the LispWorks' CAPI, I have implemented a simple program to demonstrate the basic elements in the implementation of a 2D dynamic zoom feature. The program is basically uninteresting except for the dynamic zoom implementation. You can start with the code immediately, or skip to the other side for a few comments and refer back.



Every time the user clicks the left button (:button-1), the mouse location is trapped and turned into a "model" location for a line starting position. If the user presses the right button (:button-3) the common terminal point for all of the lines is changed—the same principle of storing a "model" point is followed. Basically it is a spider with as many legs as you want originating from whatever point you select.

When the program commences, the top-left corner is (0, 0) and view scale is 1:1. The rolling of the scroll wheel of the mouse will change both of these parameters by scaling relative to the current location of the mouse pointer. Since the scroll event does not appear to provide a means to get the current mouse pointer location, I am tracking the mouse :motion "gesture" (basically an "event") and storing the pointer location, in screen coordinates, as is moves.

The change-start (I guess I should have said "add-start"—ah, scope creep) and change-center handler functions do basically the same thing. Using the current mouse position (these two handlers do receive mouse position information), determine the model point corresponding to the screen position of the mouse. The top-left corner is an actual ("model") location. So the model location that corresponds to the given screen coordinates is the top-left corner plus the model distance between the mouse pointer and the top left corner of the screen. Thus, \(M' = (M/s + TL)\), where \(M'\) is the model location. We manipulate the same formula to obtain screen coordinates (\(M\)) from the model coordinates when we display to the screen.

The part that is a little trickier is to adjust the top-left corner when dynamic zoom is applied so that you can maintain the mapping between model and screen coordinates. This is an important part of change-scale's job. Our basic formula is
$$TL_i + A_i = TL_f + A_f,$$ where \(TL_i, TL_f\) are the initial and final top-left corners (i.e., before and after zoom) and \(A_i, A_f\) are the initial and final distances between the model position of the mouse and the top-left corner. In other words, zooming does not change the mouse location in either screen or model coordinates. Hence,
$$TL_f = TL_i + A_i - A_f = TL_i + M {{s_f - s_i}\over{s_f s_i}},$$where \(s_i, s_f\) are the initial and final scale factors.

LispWorks Personal Edition for non-commercial use is available here: http://www.lispworks.com/downloads/.

Saturday, December 6, 2014

What is the Sine of 18°?

To find the exact values of sine (and cosine) of a given angle, you need to use identities. This particular value is interesting in that it takes a fair bit of work to arrive at this value analytically, and yet the value is expressible as a fairly simple radical. The first thing to notice is that 90 = 5(18). You probably know the half angle identity, but not a fifth angle identity.  And for good reason—the fifth angle identity has multiple solutions; it's a fifth degree polynomial.

I did this calculation manually the first time through—some years ago. It goes a lot easier with a computer algebra system. Today, I'm going to use Maxima to solve the problem.

An identity which comes up in the development of the complex numbers, gives a convenient way of expressing \(\sin {n\theta}\) in terms of combinations of \(\sin \theta\) and \(\cos \theta\).

De Moivre's formula:
\[r^n (\cos \theta + i \sin \theta)^n = r^n (\cos {n\theta} + i \sin {n\theta}).\]
So, if we take \(r=1\) and \(n=5\), we can do a binomial expansion and easily isolate the \(sin \theta\) and \(\cos \theta\). In a Maxima cell, I type

eq: (cos(%theta) + %i*sin(%theta))^5 = cos(5*%theta) + %i*sin(5*%theta);
realpart(eq);
imagpart(eq);

After executing, I obtain the fifth angle identities (if you want to call them that).
\[\mathrm{cos}\left( 5\,\theta\right) =5\,\mathrm{cos}\left( \theta\right) \,{\mathrm{sin}\left( \theta\right) }^{4}-10\,{\mathrm{cos}\left( \theta\right) }^{3}\,{\mathrm{sin}\left( \theta\right) }^{2}+{\mathrm{cos}\left( \theta\right) }^{5}\]
\[\mathrm{sin}\left( 5\,\theta\right) ={\mathrm{sin}\left( \theta\right) }^{5}-10\,{\mathrm{cos}\left( \theta\right) }^{2}\,{\mathrm{sin}\left( \theta\right) }^{3}+5\,{\mathrm{cos}\left( \theta\right) }^{4}\,\mathrm{sin}\left( \theta\right) \]
To eliminate \(\cos^2 \theta\) from the sine identity we use

subst([cos(%theta)=sqrt(1-sin(%theta)^2)],%);
expand(%);

and obtain
\[\mathrm{sin}\left( 5\,\theta\right) =16\,{\mathrm{sin}\left( \theta\right) }^{5}-20\,{\mathrm{sin}\left( \theta\right) }^{3}+5\,\mathrm{sin}\left( \theta\right). \]
Taking \(\theta=18°\), using \(x\) to stand in for \(\sin 18°\), and rearranging, we have
\[0 = 16\,{x}^{5}-20\,{x}^{3}+5\,x - 1.\]
Factoring this one manually was made easier by recognizing \(x = 1\) as a root. After that I tried forming a square free polynomial and found that there was in fact a squared polynomial involved. In general, this helps in a couple ways (potentially):
  1. Reduces the degree of the polynomial to be factored.
  2. Sometimes produces fully reduced factors as a by-product of forming the square free polynomial.
    1. This happens when each repeated factor is of different degree; that is, occurs a different number of times.
(Square free means no factors left that are raised to some power—each factor occurs exactly once.) Doing this manually involves taking the greatest common divisor between a polynomial and it's derivative. But today, let's let Maxima do it.

factor(0=16*x^5-20*x^3+5*x - 1);
\[0=\left( x-1\right) \,{\left( 4\,{x}^{2}+2\,x-1\right) }^{2}.\]
At this point, you really do have to be exceedingly lazy to not continue with the rest of the work manually. It's just a quadratic formula to be applied now (we discard \(x=1\) as an artifact). To use a computer algebra system for the rest is...well...okay, it is Saturday [at original publication], so...

solve(0=4*x^2+2*x-1,[x]);
\[[x=-\frac{\sqrt{5}+1}{4},x=\frac{\sqrt{5}-1}{4}]\]
We reject the negative solution (which is \(\sin {-54°}\) and, interestingly, is the negative of half the Fibonacci number) on the grounds that the range of the sine function on our interval is positive. Hence,
\[\sin 18° = \frac {\sqrt{5} - 1}{4}.\]

Saturday, October 4, 2014

Merging Uniquishly in Common Lisp

I've been working on some code with dxf files and I needed to have some generalized merging. In particular, for my first function, I wanted to merge headers without creating duplicate entries. My reader sorts the header values, so when I go to merge I just need to have a simple merge that looks at each of two lists and collects the lowest value and only advances the list that had an element "consumed" (as in, I collected it). If the values are equal, it should collect one instance, and consume from both lists (i.e., advance—take the cdr of—both lists).

The approach I've taken does not guarantee a list of unique items, only that if the original lists are sorted and contain no duplicates, then the resulting list will be sorted and contain no duplicates.  It is also non-destructive, returning a new list

(Confession: Actually, I did it the hard way the first time. Custom code that can only do one thing. Now I'm going back and amending my ways as an ardent Lisper.)

I went with a loop macro for one reason: the collect key word. This is one of the things that keeps me coming back to the loop macro. I like it overall, but when I'm struggling with the loop macro, I sometimes want to use the do* macro.  But I dislike the way do* lends itself to consing and a final call to reverse.  Hence, I come back whimpering to the loop macro.

The code:



Here's some sample output from an EMACS session:



The custom code solution was 23 lines and not so pleasant looking. The tri-test follows the custom of using a negative value to indicate a < b, positive for a > b, and 0 for a = b.  For numeric data, subtraction is the natural choice for a tri-test.

I tried using some for keywords inside of if statements and my horse (compiler) wouldn't make the jump.  But that's okay, I used a for statement at the beginning to initialize some variables (nexta and nextb) which I could then setf within do clauses as necessary.

The more I use the loop macro, the more I like it. Eventually, maybe I will discover all of the things that cause it to choke, stop doing those things, and then I will like it even more.