pages tagged computational-geometryHiking and Codinghttp://chrismgray.github.com//tags/computational-geometry/Hiking and Codingikiwiki2011-12-02T22:05:50ZConvex Hullshttp://chrismgray.github.com//posts/convex-hulls/2011-12-02T22:05:50Z2011-11-30T18:06:31Z
<div id="outline-container-1" class="outline-2">
<h2 id="sec-1">Convex Hulls Three Ways</h2>
<div class="outline-text-2" id="text-1">
<p>
Whenever I watch cooking competition shows, they always have chefs
presenting a foodstuff cooked in multiple different ways. Today,
I'm doing that with algorithms.
</p>
<p>
The algorithm in question today is the <a href="http://en.wikipedia.org/wiki/Convex_hull">convex hull</a> algorithm. In
order of implementation complexity, and descending order of
theoretical running time, there is the Jarvis March, Graham's Scan,
and Chan's Algorithm. All three are implemented in Clojure in my
<a href="https://github.com/chrismgray/convex-hull">github repository</a>.
</p>
</div>
<div id="outline-container-1-1" class="outline-3">
<h3 id="sec-1-1">Jarvis March</h3>
<div class="outline-text-3" id="text-1-1">
<p>
The simplest of the algorithms, the Jarvis March was also one of
the first output-sensitive computational geometry algorithms. In a
nutshell, you first find a point that you know to be on the convex
hull, and then you find the next point by looking at all the rest
of the points and determining which one has a segment that has the
property that all the rest of the points are on one side of it.
You repeatedly find the next point using this procedure until you
get back to the first point. There is basically no way to
implement this algorithm that does not have a running time of
\(O(hn)\), where \(h\) is the number of points on the convex hull
and \(n\) is the number of input points.
</p>
<p>
The one implementation detail of the Jarvis March that is
interesting is that whenever you see the concept of "finding the
next thing" given some previous information, the Clojure
implementation should almost always be lazy. It turns out that
implementing Jarvis March lazily will help in implementing Chan's
Algorithm, so keep that in mind.
</p>
</div>
</div>
<div id="outline-container-1-2" class="outline-3">
<h3 id="sec-1-2">Graham's Scan</h3>
<div class="outline-text-3" id="text-1-2">
<p>
Graham's Scan is one of the algorithms I remember most vividly from
the undergraduate computational geometry course that I took. The
professor, Godfried Toussaint, always referred to it as the "three
coins" algorithm, so I have kept up that tradition in some of my
function names in my implementation.
</p>
<p>
The algorithm first makes a polygon of the input points by sorting
them by angle about the bottom-most point. Then it goes around the
polygon with a stack, pushing and popping points as it goes. If
that sounds familiar, it should – it's the same idea as what I was
talking about when I brought up the idea of <a href="http://chrismgray.github.com/posts/parsing-polygons/">parsing polygons</a> a week
and a half ago.
</p>
<p>
Thus, I used the same polygon-parsing monad in my implementation as
when I computed the visibility polygons last week. It still works
just as well.
</p>
<p>
Since the points must be sorted, Graham's Scan takes \(\Theta(n
\log n)\). Sorting can be reduced to computing convex hulls, so
computing convex hulls has a \(\Omega(n \log n)\) lower bound,
meaning that this algorithm is optimal.
</p>
<p>
But Chan's algorithm is better. Weird huh?
</p>
</div>
</div>
<div id="outline-container-1-3" class="outline-3">
<h3 id="sec-1-3">Chan's Algorithm</h3>
<div class="outline-text-3" id="text-1-3">
<p>
I must confess that I had always been a little intimidated by
Chan's Algorithm. It was invented by Timothy Chan, who has a
well-earned reputation for being a genius, so I thought it would be
really complicated. It's not. There is a decent <a href="http://en.wikipedia.org/wiki/Chan's_algorithm">description</a> of it
on Wikipedia, so I won't go into the details. The gist is that you
combine the two previous algorithms that I discussed. The Jarvis
March needs to be modified so that the points can be input as a
list of smaller convex hulls, and the next point on the convex hull
is found by doing a binary search on the smaller convex hulls. But
that is really the hardest part about the algorithm. I have put
the whole thing below, because I think it's pretty beautiful.
</p>
<pre class="example">(defn chans-algorithm
"Finds the convex hull of a set of points by
the algorithm known as 'Chan's Algorithm'."
[pts]
(let [bottom-pt (apply min-key :y pts)]
(loop [m 3] ; start out with triangles
(let [partitioned-pts (partition m m [] pts)
c-hs (map grahams-scan partitioned-pts)
potential-ch (take m (apply jarvis-march pts c-hs))]
(if (= bottom-pt (last potential-ch)) ; assumes jarvis-march returns bottom-pt last
potential-ch
(recur (min (* m m) (count pts))))))))
</pre>
<p>
The great thing about Chan's Algorithm is that it is also
output-sensitive. But instead of being \(O(nh)\) (which is
\(O(n^2)\) in the worst case), it is \(O(n \log h)\), which is at
least as good as Graham's Scan, but often better. It is also quite
simple to implement, given decent implementations of Jarvis March
and Graham's Scan.
</p>
</div>
</div>
</div>
<div id="outline-container-2" class="outline-2">
<h2 id="sec-2">Conclusion</h2>
<div class="outline-text-2" id="text-2">
<p>
Convex hull algorithms are great. If I was ever to teach a
computational geometry course (admittedly that's looking like a long
shot now), I might start and finish the course with them. The
progression from the ultra-simple Jarvis March to the
more-complicated Chan's Algorithm is really nice, and there are
interesting new things to talk about the whole way. They also show
that computational geometry is not so hard to do in a functional
style. In fact, using laziness is what makes the implementation of
Chan's Algorithm so simple. So this might make a nice talk to give
people who are already into functional programming as well.
</p>
<p>
The next thing I have in mind for this project is to animate the
algorithms. Viewing the output of algorithms is already pretty easy
using Processing, but I would like to be able to see them as they
are operating. It would be great if I could do that without
changing the code too much. I have a couple of ideas, but I'm not
sure if they'll work yet.
</p>
<p>
Also, it is slightly embarrassing to admit, but my blogging system
seems to not support putting images in posts. So I am going to have
to figure out how to work around (or even fix) that before I can
show any results.
</p>
</div>
</div>
Visibility Polygonshttp://chrismgray.github.com//posts/visibility-polygon/2011-12-01T07:45:06Z2011-11-21T19:44:43Z
<div id="outline-container-1" class="outline-2">
<h2 id="sec-1">Visibility Polygon</h2>
<div class="outline-text-2" id="text-1">
<p>
I've put up some preliminary code for finding visibility polygons on
my <a href="https://github.com/chrismgray/visibility-polygon">github repo</a>. I have only tested it with one polygon, but so far
things are looking good. My idea from yesterday of using monads to
"parse" the polygon seems to be paying off. Here is the relevant
code:
</p>
<pre class="example">(defn- add-new-pt [poly]
(fn [pt stack]
(when (or (empty? stack) ; First two points are guaranteed to be visible
(empty? (rest stack))
(visible? pt (first poly) (first stack)))
[pt (rest poly) (cons (first poly) stack)])))
(defn- pop-stack [poly]
(fn [pt stack]
(let [the-seg (seg/new-seg pt (first poly))
top-seg (seg/new-seg (first stack) (second stack))]
(when (pt/left-turn? (second stack) (first stack) (first poly))
(if (seg/intersection-on-seg? the-seg top-seg)
[pt poly (cons (seg/intersection the-seg top-seg) stack)]
[pt poly (rest stack)])))))
(defn- skip-pt [poly]
(fn [pt stack]
(let [the-seg (seg/new-seg pt (first stack))
poly-seg (seg/new-seg (first poly) (second poly))]
(when (not (pt/left-turn? (second stack) (first stack) (first poly)))
(if (seg/intersection-on-seg? the-seg poly-seg)
[pt (cons (seg/intersection the-seg poly-seg) (rest poly)) stack]
[pt (rest poly) stack])))))
(defn- all-conditions [poly]
(with-monad polygon-parser-m
(m-plus (add-new-pt poly) (pop-stack poly) (skip-pt poly))))
</pre>
<p>
This defines three actions to be performed when the polygon and
stack are in certain configurations. The first executes when the
next point on the polygon doesn't obscure any of the stack. In that
case, it removes the point from <code>poly</code> and puts it on <code>stack</code>. The
second executes by popping points off the stack until the next point
on the polygon no longer obscures the stack. The final condition
activates when the next point on the polygon is hidden by the
stack. It skips the points of <code>poly</code> until it finds one that would
be visible. We combine the three conditions using <code>m-plus</code>, which
is defined to try the conditions in order until one returns
something that is not <code>nil</code>.
</p>
<p>
Is this solution any better than a functional programming solution
without monads? I think it is. First, the conditions are easy to
see and very explicitly laid out. Using a <code>cond</code> to accomplish the
same thing is certainly possible, but gets more complicated as the
number of conditions grows. Also, we are explicitly managing state
in this solution. Doing so in an ad-hoc manner would be much more
difficult. In fact, I am not really sure how I would loop until
there are no more points in <code>poly</code> in a solution without monads
(especially since not all the conditions consume a point from
<code>poly</code>). With monads, it is quite simple:
</p>
<pre class="example">(defn- visibility-polygon-helper [pt poly]
((with-monad polygon-parser-m
(m-until empty? all-conditions poly)) pt []))
</pre>
<p>
Doing something like that with a <code>reduce</code> seems like it would veer
off into unreadability fairly quickly.
</p>
<p>
Are there downsides to this solution? I think there is at least
one. That is, the conditions must be functions that return functions.
This makes them a bit more confusing than they really should be. I
needed to do this so that I could have a test in the <code>m-until</code>
function above. On the whole, though, this is a fairly cosmetic
gripe, and can be hidden by using helper functions like the one
above.
</p>
</div>
</div>
<div id="outline-container-2" class="outline-2">
<h2 id="sec-2">On the horizon</h2>
<div class="outline-text-2" id="text-2">
<p>
I would like to test the code quite a bit more. However, doing so
by drawing out polygons by hand and then figuring out their
visibility polygons is quite tedious. Therefore, I need some sort
of GUI to be able to draw polygons and then I will be able to
see whether the computed visibility polygon makes sense.
</p>
<p>
I am thinking that it might be nice to use the <a href="https://github.com/rosado/clj-processing"><code>clj-processing</code> library</a> to do this. Unfortunately, that library is currently only
using Clojure 1.2. Since I am using Clojure 1.3, that might be a
problem. So I might need to do some porting from 1.2 to 1.3.
However, doing so should give me a good idea of how processing
works, and could give a nice tool for more interactive geometry
programs.
</p>
</div>
</div>
Parsing Polygonshttp://chrismgray.github.com//posts/parsing-polygons/2011-12-01T07:44:52Z2011-11-21T00:58:28Z
<div id="outline-container-1" class="outline-2">
<h2 id="sec-1">Functional programming</h2>
<div class="outline-text-2" id="text-1">
<p>
I've gotten a bit obsessed by functional programming over the last
couple of years. I have even seriously thought about writing a book
that looks at computational geometry (which is the subject in which
I was trained) in the light of functional programming. Most of the
standard computational-geometry texts these days approach the
writing of code from an imperative standpoint.
</p>
<p>
I think functional programming is a good way to think about
computational geometry for a couple of reasons. First, most
problems in computational geometry can be expressed in a purely
functional manner. That is, the answers to the problems are usually
the same given the same inputs. Secondly, I have bought into the
idea that parallel and distributed computing are much easier when
state is not being modified. Functional programming forces you to
think this way. While many problems in computational geometry are
inherently sequential, not all of them are. Using a programming
method that allows for the easy addition of more cores seems like a
good practice when a problem is encountered that is easily made
parallel.
</p>
<p>
One of the things I have been thinking about as a final chapter for
the book that I would like to write is the introduction of monads as
a way to "parse" polygons. Monads are commonly used (well, commonly
in the functional programming world) to parse strings. They
generally scan through with the aid of a stack. There are many
algorithms that scan through the points of a polygon with a stack --
the example that I have thought about the most is the visibility
polygon.
</p>
</div>
</div>
<div id="outline-container-2" class="outline-2">
<h2 id="sec-2">Visibility polygons</h2>
<div class="outline-text-2" id="text-2">
<p>
A <i>visibility polygon</i> is the subset of a polygon that can be "seen"
from a point inside the polygon. Here, we regard the segments of
the polygon as opaque walls.
</p>
<p>
The standard algorithm for computing the visibility polygon is not
too difficult to implement in an imperative language. In fact, I
<a href="http://www.win.tue.nl/~cgray/ipelets.html">implemented</a> it already in C++ for the drawing editor <a href="http://tclab.kaist.ac.kr/ipe/">ipe</a>. (Though I
think that code is now obsolete since the release of ipe 7.)
</p>
<p>
I've already tipped my hand at how I would think about implementing
such an algorithm in a functional language. I would use a variant
of a monadic parser to go through the polygon and find the parts
that are visible. (That last sentence doesn't give much of a hint
about how it's done, but the whole algorithm would take too much
space to describe. Suffice it to say that a stack is kept with the
"currently visible" portion of the polygon always on the stack.)
</p>
</div>
</div>
<div id="outline-container-3" class="outline-2">
<h2 id="sec-3">Monads</h2>
<div class="outline-text-2" id="text-3">
<p>
I used monads <a href="http://chrismgray.github.com/posts/time-series-2/">earlier</a> to make a framework for creating and comparing
time-series. However, I only used a monad that was already defined
by someone else. I think I might need to define a monad on my own
this time – a slightly daunting idea. I definitely need to read
the rest of <a href="http://intensivesystems.s3-website-us-east-1.amazonaws.com/tutorials/monads_201.html">Jim Duey's tutorial</a> to make sense of how to do it.
</p>
</div>
</div>