Monday, 10 September 2012

Halton sequences, at last

A while ago I posted on generating random point sets, and Sean Gillies suggested Halton sequences as a way of generating nice-looking point distributions.  I finally got around to looking into this idea. Indeed he was correct - Halton sequences provide a very pleasing appearance for point distributions (much better than purely random distributions).
The above shows the first 1000 elements of the Halton sequence with bases (2,3).  Larger bases exhibit progressively more coherence, which may or may not be desirable depending on application. Here's H(7,9):
 The code to generate Halton sequences is pretty trivial - see the Wikipedia article for some pseudo-code.  (Although I have to say that the explanation of their derivation is pretty poor.  I've seen this sad trend on some other Wikipedia mathematical articles as well - in particular the one for Barnes surface interpolation.  I suppose the rejoinder would be to get in there and improve it!)

Keying off something else that Sean showed on his blog, I also played around with using Halton sequences to generate quasi-random sets of polygons.  Here's one example, using a slightly different approach to Sean's:

And the same for H(5,7), which looks slightly nicer I think:


Might come in handy someday for generating test data, or possibly as a visual texture.

The complexity of simple

Recently Sean Gillies and others have been looking at the concept of geometric simplicity as implemented in JTS/GEOS, specifically in the case of polygons.  The TL;DR is that JTS has a bug because it reports the following polygon as simple:

POLYGON ((39.14513303 23.75100816, 97.05372209 40.54550088, 105.26527014 48.13302213, 100.91752685 58.43336815, 71.56081448 83.55114801, 60.71189168 86.25316099, 62.00469808 75.1478176, 83.16310007 42.82071673, 92.82305862 37.19175582, 95.99401129 26.47051246, 106.22054482 15.51975192, 39.14513303 23.75100816))






This polygon is invalid, because the shell ring self-intersects.  This means the ring is non-simple - but it's not so clear-cut for the polygon itself.

In fact, reporting the polygon as simple isn't a bug - it was a deliberate design decision.  The original rationale was:

  • the OGC Simple Features Specification in section 2.1.10 defining polygon geometries states that "Polygons are simple geometries".  (This is consistent with the requirement that valid polygons contain no self-intersections).  
  • JTS usually follows the principle that computed results are only well-defined when the input geometries are valid.  This is because it can be expensive to check for validity in order to handle invalid inputs gracefully.  Validity testing is factored out into a separate method (isValid) in order to allow the client to decide when to incur the code of executing it.
For these reasons the initial JTS implementation of isSimple returned true for all polygonal inputs, thus avoiding the cost of checking for self-intersections.

However, perhaps it would be more useful to actually carry out some testing of simplicity.  This raises the question of what exactly are the semantics of isSimple for polygons?  Looking to the SFS for guidance, it has the following somewhat imprecise general definition of simplicity:
the Geometry has no anomalous geometric points, such as self intersection or self tangency. The description of each instantiable geometric class will include the specific conditions that cause an instance of that class to be classified as not simple.
For polygons the spec simply states:
Polygons are simple geometries
(In contrast, for linear geometry the specification has a rigorous mathematical definition of simplicity.)

Apart from the "always simple" semantics, there are at least two other options:
  1. Polygons are simple if their component rings are simple (i.e. have no self-intersections)
  2. Polygons are simple if they have no topological anomalies.  This is equivalent to checking whether the polygon is valid.
Because option #2 simply replicates the semantics of isValid, it seems more useful to adopt option #1 to provide more possibilities for analyzing polygon topology (but I'm definitely open to suggestions about why #2 would be better!)

This change to the semantics of isSimple has been implemented in JTS and will appear in the 1.13 release.  The code will also be extended to handle GeometryCollections, with the semantics that they are simple if all their components are simple.

We're all in deep map now

A few thoughts on the Atlantic article on GooMaps...

  • I love the term "the deep map".  I'm not doing GIS any more - I'm doing deep mapping!
  • Yay!  Geospatial conflation gets a mention in the MSM!
  • The article name-checks Borge's map.  The original story pointed out the paradox of a 1:1 scale map - but this is only paradoxical in 2-space.  In the infinite-dimension virtual world it's possible to have maps that are more detailed than the physical world (in that they can represent relationships across time and other more abstract dimensions)
  • That is some great synergy between StreetView data acquisition and Google's work on self-driving cars
  • The author is convinced that Google has an unassailable lead on building the virtual map of the world.  But without knowing what's going on at Apple (and perhaps Microsoft) I wouldn't be so sure about that.  (And even Google isn't as transparent as this realtime view of OSM edits.)
As always, James Fee has a pithy comment and a great picture. The good news is you get a ping pong table in the office.  The bad news is that 2000 more map edits just came in, so get back to work!