## Not famous yet

Calm down, everyone. I didn’t discover a new constant, as tantalizingly suggested (tongue-in-cheek) in my previous post.

Recap: While writing programs to calculate pi, I experimented with a circle-drawing algorithm called Bresenham’s Circle Algorithm. It turns out that if you pretend that it’s a real circle, add up the length of its perimeter, and divide by its diameter, you get a number that’s a constant but not equal to pi. It’s about 3.31371. I ended yesterday’s post by calling this number $\pi_b$ and getting all excited about my upcoming fame. Wow, a new mathematical constant!

Here’s the update. Today, I have an exact value for $\pi_b$ – it’s equal to $8(\sqrt2-1)$. As your calculator will demonstrate, that’s about 3.31371. I was right! A new constant!

Except it’s not new. Or all that interesting.

See, it turns out that if you decompose a Bresenham circle (which is not a circle at all, but a polygon) of a given radius into the line segments that make it up, you can always rearrange it into a shape that’s very close to an octagon, which encloses the true circle of the same radius. And the octagon gets closer and closer to perfect as the radius gets larger. And the perimeter of a perfect octagon wrapped around a circle of radius r is equal to $16r(\sqrt2 - 1)$. Which, when divided by twice the radius, equals “$\pi_b$“.

The same reasoning could produce other kinds of pi, such as “$\pi_{square}$“, which would be equal to 4.

However, none of these other kinds of pi would be of any use in determining the value of the classical pi, which is what I was after. They’re just constants, like, say, 42.

On the bright side, my conjecture that $\sqrt2$ was baked into the Bresenham Circle somehow is definitely true. Note how it’s part of the perimeter formula, and also note all the 45-degree angles inside an octagon. There’s also a nice explanation for the observation that the number of points generated by the algorithm is always equal to $4r\sqrt2$. You’re welcome to look for it if you’re so inclined.

And that concludes this year’s Pi Day festivities!

## A Different Kind of Pi

In my last post I left my readers (Hi Mom!) with a teaser about a possible improvement to my pi-calculating code. Let’s see how it panned out.

My naïve arc-drawing code used two square-root calculations per point – pretty slow. I figured that there must be a quicker version out there, and sure enough, the gold standard in translating curves to pixels is something called Bresenham’s Algorithm. (Edit: the Wikipedia article is awful. Here’s a better one.) Here’s my code, cribbed from this page:

function result = bresarc(r)
% BRESARC(R) - draw a quarter-arc of a circle using Bresenham's algorithm

r = abs(round(r));
x = -r;
y = 0;
err = 2 - 2 * r;

idx = 1;
result = zeros(round(sqrt(2) * r ), 2);

while x <= 0
result(idx, [1 2]) = [-x y];
idx = idx + 1;
t = err;
if t > x
x = x + 1;
err = err + 2 * x + 1;
end
if t <= y
y = y + 1;
err = err + 2 * y + 1;
end
end
end


It’s about four times faster than my arcwalk function, and produces a better approximation to the circle, with fewer points.

Bresenham arc, compared to arc and simplified Bresenham arc

Here’s an oddity. Look at line 10, where I pre-calculate the number of points to return in the result. It turns out (I discovered by experimentation) that this algorithm always generates $radius*\sqrt2$ points for some reason. Something this obvious, in an algorithm that’s very well known, should probably already be common knowledge – but I couldn’t find any evidence of it by Googling for a few minutes. Perhaps my 15 minutes of fame are on the way!

Let’s see how it does in my pi-calculating algorithm. The original form, slightly updated from last time (better memory management):

clc
clear

maxpow = 5;
epsilon = sqrt(2);

tic

for r = 10 .^ (1:maxpow)
pts = arcwalk(r);                       % make the arc
pts = dpsimp(pts, epsilon);             % simplify it
pts = pts(2:end,:) - pts(1:end-1,:);    % vector subtraction
s = sum(hypot(pts(:,1), pts(:,2)));     % add up the lengths
end

toc


The output:

For radius 10, pi is about 3.164823
Elapsed time is 41.948334 seconds.


Swapping out arcwalk for bresarc, we get:

For radius 10, pi is about 3.170979
Elapsed time is 29.530892 seconds.


The killer that remains is the DP polyline simplification routine dpsimp. Hmmm…. the new arc algorithm is closer to the circle. How would it work if I just calculated the arc length of the output of the bresarc function directly, without attempting to simplify it? It will be longer than the true arc (I think; proof is left as an exercise for the reader), but… how much more? Let’s find out.

After commenting out line 11 above, we get:

For radius 10, pi is about 3.297056
Elapsed time is 0.283709 seconds.


Holy crap! 100 times faster! Farewell, favorite recursive algorithm.

Hmmm. Is that value converging? Let’s see. Bumping up maxpower from 5 to 7 and prettying up the output:

For radius       10, pi is about 3.297056274848
Elapsed time is 29.052004 seconds.


It kinda looks like it. So what have we discovered here? A new constant? It would be the ratio of the perimeter length of a Bresenham circle to its diameter. Call it pi-sub-b, or $\pi_b$. For now, let’s assume that $\pi_b\approx3.31371$. We know that this number is transcendental like pi, because pi_b is made up of line segments with two possible lengths: 1 and sqrt(2) (which is also transcendental).

An obvious hack: we can use the (presumed) fact that pi_b is a constant to calculate pi. We figure out the difference between pi and pi_b, then subtract it from the calculated value of pi_b:

EDU>> pib = 3.31371

pib =

3.31371

EDU>> pidiff = pib - pi

pidiff =

0.172117346410207

EDU>> 

Let’s see how it works:

clc
clear

maxpow = 7;
pidiff = 0.172117346410207;

tic

for r = 10 .^ (1:maxpow)
pts = bresarc(r);                       % make the arc
pts = pts(2:end,:) - pts(1:end-1,:);    % vector subtraction
s = sum(hypot(pts(:,1), pts(:,2)));     % add up the lengths
end

toc

For radius       10, pi is about 3.124938928438
Elapsed time is 29.198759 seconds.

Of course, that’s absolutely cheating. I baked the value for pi into the program, as part of the pidiff fudge factor.

But the possibility exists that there’s a more direct route to this result. Pi, the square root of two, and pi_b may all be related in a deep way. Almost certainly the Bresenham algorithm has $\sqrt2$ baked into it, as seen at the beginning of this post.

Consider a Bresenham circle. It has a perimeter length. The length won’t change if we rearrange the line segments it’s made of. It is made of horizontal lines of length 1, vertical lines of length 1, and diagonal lines of length $\sqrt2$. It must be possible to rearrange these into a square with the corners cut off.

Ah, the ancient dream of squaring the circle! My 15 minutes of fame approaches!

## Playing with Pi

My old college buddy Mark VandeWettering tweeted a challenge right before Pi Day: he was going to write code to calculate a value for pi in a novel way.  What were the rest of us slobs going to do?

I had a pretty good idea how he was going to do it, as he’d recently posted something really cool about a way to use the Mandelbrot set to find pi.  (Sure enough that’s what he did.)

I decided to play around with using MATLAB (my current favorite programming environment) to see how many interesting ways I could calculate pi.  I stayed up all night, had fun, learned a lot, and came up with a couple of things worthy of posting.  (At least by the standards of this blog.  Mark doesn’t have anything to worry about.)

To warm up, I wrote a Monte Carlo simulation, with a small wrinkle: I got a load of random numbers from random.org and used those, then compared it to MATLAB’s rand function.  No noticeable difference.

clc
clear

%% Happy Pi Day 2011!
%
% This program uses a Monte Carlo technique to estimate the value of pi.
% Using a list of random numbers between 0 and 1, it uses every possible
% pair of numbers to create [X,Y] points.  These points are then tested to
% see whether they lie within the unit circle.  The ratio of points in the
% circle to total points should be approximately pi/4.

%% ***** Get 10000 random numbers

% Downloaded from random.org - numbers evenly distributed between 0 and 1
% rcount = length(randnums);

randnums = rand([10000,1]);
rcount = length(randnums);

%% ***** Calculate pi

incircle = 0;
trials = 0;

tic

for ii = 1:rcount-1
for jj = ii:rcount
pt = [randnums(ii), randnums(jj)];
len = norm(pt);     % Calculate the distance from the origin
if len < 1          % It's inside the circle
incircle = incircle + 1;
end
if len ~= 1         % Don't count it if it's ON the circle
trials = trials + 1;
end
end
end

fprintf('Happy Pi Day 2011!  For %d trials, Pi is about %f\n', ...
trials, incircle/trials * 4)

toc


Output:

Happy Pi Day 2011!  For 50004999 trials, Pi is about 3.152321
Elapsed time is 73.081508 seconds.


I noticed that this technique converges excruciatingly slowly, so I figured anything I could come up with would have to be faster.

I looked through the Wikipedia articles about estimating pi and was impressed with Ramanujan’s summation formulae, but I didn’t understand where they come from so it seemed like cheating.

Now, I was able to follow the derivation of Machin’s arctan formula, but doing that with a computer seemed… trivial.

EDU>> 4*(4*atan(1/5)-atan(1/239))

ans =

3.14159265358979

EDU>>

Should I try to do something similar with the arcsin function? Naaah, it was getting late and my talents seem to lie mostly with computer algorithms rather than number theory. Let’s try some naïve brute-force code.

And hey, instead of adding up the area inside the circle like in the Monte Carlo simulation, let’s see if we can figure out how long the arc is. That should take a lot fewer points.

Here’s a brute-force arc-drawing routine:

function result = arcwalk(r)
% ARCWALK(R) - generate pixels for quarter-circle arc of radius R

% preallocate array of points
pts = zeros(2*r+1, 2);

x = r;
y = 0;
idx = 1;

%% Walk from [r,0] to [0,r], moving up or left
while x >= 0

pts(idx,[1 2]) = [x y];

up = norm([x, y+1]);
left = norm([x-1, y]);

if abs(r - up) < abs(r - left)
y = y + 1;
else
x = x - 1;
end

idx = idx + 1;
end

result = pts;

end


Calling arcwalk(10) produces these points:

The arcwalk function's output for radius=10

Kewl, now all I have to do is count up the number of pixels it took, assume that’s approximately how long the arc is, and Bob’s your uncle:

clc
clear

maxpow = 5;

tic

for r = 10 .^ [1:maxpow]
pts = arcwalk(r);
s = length(pts);
end

toc

For radius 10, pi is about 4.200000
Elapsed time is 1.026044 seconds.


Wha… oh. Maybe not quite that naïve. Bonus points if you figure out what’s going on. Hint: I had already figured this out in the arcwalk function, right there on line 5.

What if I simplify that arc? Hey, I can use one of my favorite algorithms, the recursive Ramer-Douglas-Peukert algorithm! I used it last summer while (ahem!) interning at NASA. Here’s my implementation:

%% DPSIMP - quick implementation of the Douglas-Peucker algorithm
% PS = DPSIMP(P, E) - simplify polyline P with epsilon E
%
% DPSIMP uses the Ramer-Douglas-Peucker algorithm to simplify a polyline
% http://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm
%
% P is an N-by-2 matrix containing an X,Y point in each of the N rows.
% E is the value of epsilon (the maximum distance from the simplified
%   polyline to the given polyline).
% PS is an N-by-2 matrix containing the X,Y points of the simplified
%   polyline.

% by Doug Weathers, weathers@nmsu.edu

function result = dpsimp(p, epsilon)

%%
% pldist - find distance from point to line
% The point is p0, and the line runs through p1 and p2.  The points are row
% vectors [x y].
function result = pldist(p0, p1, p2)

% First, put the line in the form Ax + By + C = 0
A = p1(2) - p2(2);
B = p2(1) - p1(1);
C = -A * p1(1) - B * p1(2);

% Then, calculate the orthogonal distance with this simple formula:
result = abs(A * p0(1) + B * p0(2) + C) / sqrt(A^2 + B^2);

end

%%
% Find point with maximum distance from the line that runs between the first
% and last points.

dmax = 0;
index = 0;
[rows, cols] = size(p);

% debug
% Draw a line between the endpoints.  Delete it later if it's not
% simplified enough.  This lets you watch the recursion in action.
% h = line([p(1,1) p(end,1)], [p(1,2) p(end,2)], 'LineStyle','-', 'Marker','x');
% pause

for i = 2:rows-1
d = pldist(p(i,:), p(1,:), p(end,:));
if d > dmax
index = i;
dmax = d;
end
end

% Is there a point farther away than epsilon?  If so, we need to simplify
% further.

if dmax >= epsilon

% Yes, there's a point far enough away from the line to be included in
% the simplified version.
% Split the line into two segments and call this routine recursively on
% each segment.
recresult1 = dpsimp(p(1:index,:), epsilon);
recresult2 = dpsimp(p(index:end,:), epsilon);

% Connect the simplified polylines back together!  Only keep one copy
% of the point at the index.
result = [ recresult1(1:end-1,:); recresult2 ];

% debug
% The line we drew earlier needs to be erased
% delete(h);

else

% None of the interior points of the polyline are far enough away from
% the line joining the endpoints to worry about.
% We're done!  We can skip all the intervening points!
result = [ p(1,:); p(end,:) ];
end

end


Here’s how it looks when applied to the arcwalk(10) pixels, with an epsilon of 1:

Quarter-arc, arcwalk pixels, and the DP simplified polyline

Now, how to pick a good value for epsilon? I played with a few values and settled on the square root of 2. It might be worth trying to optimize this some day.

Let’s try to calculate pi with this technique. (Observe how my comments are getting briefer as time passes.)

clc
clear

maxpow = 5;
epsilon = sqrt(2);

tic

for r = 10 .^ (1:maxpow)
pts = arcwalk(r);                       % make the arc
spts = dpsimp(pts, epsilon);            % simplify it
path = spts(2:end,:) - spts(1:end-1,:); % vector subtraction
s = sum(hypot(path(:,1), path(:,2)));   % add up the lengths
end

toc

For radius 10, pi is about 3.164823
Elapsed time is 42.576056 seconds.


Reasonably accurate, but pretty slow. I suspect the problem is with recursion in MATLAB – it’s really inefficient.

Let’s see if we can speed it up and get a better answer by using some calculus. I looked up the formula for the length of the curve of a function, plugged the equation for our quarter-circle into it, simplified it, and ended up with the following code:

clc
clear

tic

pts = 1000000;
x = linspace(0, sqrt(2)/2, pts);
dx = x(2)-x(1);
s = sum(sqrt(1 ./ (1 - x.^2))) .* dx;

fprintf('\nFor %d points, pi is about %f\n', pts, 4*s)

toc

For 1000000 points, pi is about 3.141596
Elapsed time is 0.078690 seconds.


Stupidly fast! A million points, forsooth. I calculated only one-eighth of the circle instead of one quarter, but the big time-saver was the opportunity to vectorize the code. In MATLAB you can perform element-by-element operations on arrays by preceding the operator with a dot. It’s worth doing. No loops – I made an array of all the x-values, then calculated a y-value for each x in one operation.

After this I went to bed and tried to ignore the birds greeting the sun. It wasn’t too hard.

Note: during the writing of this post, I got to wondering about my naïve arc-drawing algorithm. Is there a better one available? Of course there is, and it’s got some serious potential of speeding up my brute-force code. More later – right now I need to see if I can get to sleep before the birds wake up.

## Too many keywords to ignore

I just had to post this.

Way too many keywords to ignore:

- computer
- Apple
- astronomy
- Lego
- Antikythera Mechanism

I think they only missed Star Trek and Ratchet and Clank.

## Pictures and video of runway dedication

I’m supposed to be studying, and I haven’t gotten to this part yet, but I thought I’d post some of this great video from the Spaceport runway dedication.

Official Virgin Galactic footage (contains amazing video from the air, you gotta watch this):

Jeff Foust’s photos:

Photo gallery of the Spaceport America event

My photos (only on Facebook for the moment, sorry):

Spaceport runway dedication

OK, back to studying.

## Day 1 of ISPCS 2010

I blew off a week of school to attend the 2010 International Symposium for Personal and Commercial Spaceflight.  It was well worth it.  Thanks to the good people at the New Mexico Spacegrant Consortium, who not only put on the symposium but also paid for my ticket.  If you’re a student, look into what they can do for you.

This post contains my personal highlights of the symposium.  For the facts, you probably can’t beat Jeff Foust’s coverage of Day 1 and Day 2.

For me, it started on Tuesday.  I was asked to give a 5-minute talk on my Fachaba project.  I sat next to Rick Homans, the new executive director of Spaceport America.  He doesn’t know it yet, but he’s my future boss.

I got lots of interest and questions after our portion of the program.  Mark Severance, NASA’s ISS National Lab Education Projects Manager, wanted to talk to me about the payload, which to his delight I had in my backpack.  I’m not sure how my payload would be appropriate for ISS, since an inertial tracker would only tell us the position of the ISS, which we already know.  Plus, it wouldn’t work anyway – free fall = no forces = no accelerations = no readings from an accelerometer.  He was just interested in what we had done.  I get the idea that he loves his job.

Richard Mains, from NASA’s CRuSR program, made a beeline for me after the presentation.  His team is working on figuring out how to make the best use of the upcoming family of reusable suborbital vehicles to do research.  The traditional payload is expensive, generally not reusable, and takes two years to deliver.  We did ours in six weeks for \$2500 and got it back intact.  He was very interested.  We swapped business cards.

Wednesday was Day 1 of the ISPCS conference proper, held at the New Mexico Farm and Ranch Heritage Museum.  I missed the first panel – I was tapped to assemble Virgin Galactic’s model of WhiteKnightTwo carrying SpaceShipTwo.  It was a somewhat daunting task – it came in two huge crates with no instructions, had a wingspan of around two meters, and probably cost more than a flight on the real vehicle!  But I figured it all out and got it assembled and on stage (with some help from the Two Aarons).  After that I settled in to enjoy the conference.

On a whim I pulled out my iPod Touch and got online and began posting tweets of items of interest, using the hashtag #ISPCS.  A couple of others were doing the same thing, notably Jeff Foust, Jeff Krukin, @spacecom, and Kevin Russell.  I picked up some followers (hi, folks!) and thank-yous (you’re welcome!).

After the panels we had a wonderful reception at the Hotel Encanto, put on by the Spaceport Sweden people.  They know how to throw a party!

My main takeaway from Day 1: I haven’t been paying enough attention to Masten Space Systems and Sierra Nevada Corporation.

Coverage of Day 2 will have to wait — I have a midterm tomorrow.  See you then!

## Introducing Doug to NASA

This term’s crop of NASA interns (me, Charles, and Elise) will be giving presentations to the NASA employees later today.  This is to introduce us to the staff.  (I don’t know why we’re being introduced after three weeks of gainful employment.  Must be a NASA thing.)

Here’s my presentation. I did it using prezi.com, which I highly recommend.  Instructions:

• Click on the Play arrow at the bottom of the picture.
• Use the “More” menu at the bottom right to go to full screen.
• Once you’re done looking at the first screen, press the right arrow key to go to the next stop along the path.
• Press the left arrow key to navigate backwards along the path.
• If you see something interesting, you can click on it.  Pressing the left/right arrows will put you back on the path.
• You can click and drag the canvas to scroll around.
• You can zoom in and out with the scroll wheel or the up/down arrows.
• Press the Escape key to get your browser window back.

Note: it’s written in Flash.  Be informed.  Be safe.

I’ll let you know how it was received.  <crosses fingers>

## Successful flight

Sorry I’m late on this update, but the rocket with our payload flew on May 4.  It was a perfect flight, and I got the payload back the next day.  In perfect condition.

I had a great time parading around the engineering department of NMSU, showing off the payload to everyone.  Everyone was impressed, including the department head, Dr. Burton.  He wants the data.

The way the payload worked was to begin recording data as soon as it was turned on, and to continue recording until it’s turned off.  Each time it’s turned on, a new data file gets created.  It stores about one megabyte of data per minute, and the flight lasts about 15 minutes.  So I expected to see one new data file on the storage card, and I expected it to be 15 megabytes or larger.

And that’s exactly what I saw when I stuck the memory card into the reader.  One new file, 17.5 megabytes.

A quick pass through the file shows that the information is readable and appears to be the right information, in the right order.  Over 300,000 individual readings were taken, and not one of them failed the checksum test (indicating that every record had all its bits recorded correctly with no errors).

Let’s pause for a moment and appreciate the enormity of this achievement.  A lot of hard work by a lot of people, and something me and my friends built got flown into space, did its job, came back with the data, and is ready to fly again.  I am awed and humbled.  Also hooked.  I really want to do this again.

The format of the data file is pretty crude.  With the mad scramble to get the payload delivered on time, I gave up any pretense of calculating anything and just stuffed the raw sensor readings into the file, to be sorted out later.  Well, it’s now “later”, and I’m working on a program to read the data file.

## Fachaba Project presentation

As previously mentioned, I did a presentation about the Fachaba Project at the Southwest Regional Technology Symposium.  Here’s the PowerPoint (and a PDF of same) for your enjoyment.

Note: you’ll probably have to right-click the links below to save them to your hard drive.  Sorry about the inconvenience – my download section is broken at the moment.

PowerPoint file (.pptx)

PDF file (.pdf)

## NASA is not Space, Space is not NASA

(I just watched “Orphans of Apollo“, and it primed me for this rant.)

I keep seeing people make this mistake. Stephen Colbert and Neil DeGrasse Tyson made it, Florida politicians make it, and now a bunch of astronauts have signed an open letter to President Obama that makes it.

NASA does not equal Space.  Space does not equal NASA.

More specifically, NASA is not America’s only way to get people into space.

Back to that letter.  Quoting:

America’s greatness lies in her people: she will always have men and women willing to ride rockets into the heavens. America’s challenge is to match their bravery and acceptance of risk with specific plans and goals worthy of their commitment. NASA must continue at the frontiers of human space exploration in order to develop the technology and set the standards of excellence that will enable commercial space ventures to eventually succeed. Canceling NASA’s human space operations, after 50 years of unparalleled achievement, makes that objective impossible.

See the mistake?  Canceling NASA’s human space operations won’t stop America’s men and women from riding rockets into the heavens.  All they need to do is buy a ticket from Virgin Galactic, or SpaceX, or XCOR.  All of which are American vehicles, designed and built in the US of A.  And then there’s Russia, Chinese, Japan, the ESA, and India, all of whom would likely be glad to fly Americans and American cargo into space.

You know how people lament “Where’s my flying car?”  We should be lamenting “Where’s my PanAm shuttle to orbit?”  or “Where’s my vacation on the moon?”

That was what NASA should have been doing during those 50 years of “unparalleled achievement”.  Developing those technologies they promised us!

What did we get instead?  The Space Shuttle.  Commercially useless at a billion bucks per launch, and a turnaround time measured in months.  Compare that to an airplane, and you can see how NASA has wasted our time and money.  I could weep.

So, slowly, finally, we’re starting to get a commercial space industry that’s focused on commercial viability, not on whatever NASA’s priority is.  (My opinion: they’re a government bureacracy, focused on keeping their jobs, so they literally cannot create a space transportation system that takes fewer than 30,000 people to operate. )

Now it’s time to get NASA out of the business of transporting cargo and humans into space, and refocus them on what they do best.  Big, grand, nearly impossible tasks like expanding humanity into the solar system.  We need to know how to live on other worlds.  We need to know how to get there.  We need an infrastructure to support us, like climbers on the way to the summit of Mt. Everest need a base camp, and supply caches, and maps, and guides.  NASA can build all of these things – if it stops messing around with uneconomical rockets and blowing its budget on things that private industry can do far cheaper and faster.

This is exactly what the Obama administration is advocating.  Good for them for listening to their advisors and taking an unpopular stance.

NASA has failed us by not developing economical transport into orbit.  We need to get NASA out of the way of American engineers and entrepreneurs, and let them get us into orbit.  And focus NASA on the big, not-yet-commercial exploration and development tasks that only NASA can do.

Edit: Searching for something else, I found this PBS Radio transcript of a quick story featuring my personal hero and role model, Jeff Greason.  He also thinks this way.  In fact, I probably got it from him.  He was on the Augustine Commission.

Posted in | 1 Comment