# Archive for January, 2012

## Fast Calculation of Voronoi Polygons in PovRay– Applied

Posted by smathermather on January 30, 2012

In a further exploration of using PovRay to do fast calculation of Voronoi polygons, let’s look to a real stream system as an example.  Here’s where the magic comes out, and where medial axes are found.

First the simple but real example.

Here’s the povray code:

```
#include "transforms.inc"
#version 3.6;

#include "edge_coords1.inc"

#declare Camera_Location = <2258076, 10, 659923>;
#declare Camera_Lookat   = <2258076, 0, 659923>;
#declare Camera_Size = 1800;

camera {
orthographic
location Camera_Location
look_at Camera_Lookat
right Camera_Size*x
up Camera_Size*y
}

background {color <0.5, 0.5, 0.5>}

//union {

#declare Rnd_1 = seed (1153);

#declare LastIndex = dimension_size(edge_coords, 1)-2;
#declare Index = 0;
#while(Index <= LastIndex)

cone {
<0.0, -5, 0.0>, 30,  <0.0, 5, 0.0>, 0.0
hollow
material {
texture {
pigment {
rgb <rand(Rnd_1)*2 + 1, rand(Rnd_1)*5 + 2, rand(Rnd_1)*5 +5>
}
finish {
diffuse 0.6
brilliance 1.0
}
}
interior {
ior 1.3
}
}
translate edge_coords[Index]
}

#declare Index = Index + 1;

#end

```
```
#declare edge_coords = array[16842]
{<2256808.71000000000, 0.00000000000, 660094.46988100000> ,
<2256810.06170000000, 0.00000000000, 660092.99580200000> ,
<2256811.41308000000, 0.00000000000, 660091.52172400000> ,
<2256811.68965000000, 0.00000000000, 660091.21988700000> ,
<2256813.46393000000, 0.00000000000, 660090.29698900000> ,
<2256815.23820000000, 0.00000000000, 660089.37376200000> ,
<2256817.01248000000, 0.00000000000, 660088.45086400000> ,
<2256818.78675000000, 0.00000000000, 660087.52763700000> ,
<2256820.56103000000, 0.00000000000, 660086.60473900000> ,
<2256822.33530000000, 0.00000000000, 660085.68151200000> ,
<2256824.10958000000, 0.00000000000, 660084.75861400000> ,
<2256825.88352000000, 0.00000000000, 660083.83538800000> ,
<2256827.65780000000, 0.00000000000, 660082.91248900000> ,
<2256829.43207000000, 0.00000000000, 660081.98926300000> ,
<2256831.20635000000, 0.00000000000, 660081.06636400000> ,

```

## Fast Calculation of Voronoi Polygons in PovRay

Posted by smathermather on January 20, 2012

Yet further abuse to follow in the application of PovRay for spatial analyses– be forwarned… .

Calculating Voronoi polygons is a useful tool for calculating the initial approximation of a medial axis. We densify the vertices on the polygon for which we want a medial axis, and then calculate Voronoi polygons on said vertices. I’ll confess– I’ve been using ArcGIS for this step. There. I said it. My name is smathermather, and I use ArcGIS. It’s a darn useful tool, and, and, I’m not ashamed to say that I use it all the time.

But, I do like the idea of doing this with non-proprietary, open source tools. I’ve contemplated using the approach that Regina Obe blogs about using PL/R + PostGIS to calculate them, but the installation of PL/R on my system failed my initial Yak Shaving Test.

So, with some google-fu, I stumbled upon articles on using 3D approaches to calculate 2D Voronoi diagrams, and my brain went into high gear. Right cones of equal size, viewed orthogonally generate lovely Voronoi diagrams, the only caveat being you need to know what your maximum distance you want to calculate in your diagram. For my use cases, this isn’t a limitation at all.

And so we abuse PovRay, yet again, to do our dirty work in analyses. I’ll undoubtedly have some follow up posts on this, but in the mean time, some diagrams and really dirty code:

```
global_settings {
ambient_light rgb <0.723364, 0.723368, 0.723364>
}

camera {
orthographic
location < 0.0, 5, 0>
right x * 3
up y * 3
look_at < 0.0, 0.0, 0.0>
}

light_source {
< 0.0, 100, 0>
rgb <0.999996, 1.000000, 0.999996>
parallel
point_at < 0.0, 0.0, 0.0>
}

light_source {
< 0.0, 100, 100>
rgb <0.999996, 1.000000, 0.999996>
parallel
point_at < 0.0, 0.0, 0.0>
}

light_source {
< 100, 100, 0>
rgb <0.999996, 1.000000, 0.999996>
parallel
point_at < 0.0, 0.0, 0.0>
}

cone {
<0.0, -0.5, 0.0>, 0.5,  <0.0, 0.5, 0.0>, 0.0
hollow
material {
texture {
pigment {
rgb <0.089677, 0.000000, 1.000000>
}
finish {
diffuse 0.6
brilliance 1.0
}
}
interior {
ior 1.3
}
}
}

cone {
<0.0, -0.5, 0.0>, 0.5,  <0.0, 0.5, 0.0>, 0.0
hollow
material {
texture {
pigment {
rgb <1, 0.000000, 1.000000>
}
finish {
diffuse 0.6
brilliance 1.0
}
}
interior {
ior 1.3
}
}
translate <-0.5,0,-0.5>
}

cone {
<0.0, -0.5, 0.0>, 0.5,  <0.0, 0.5, 0.0>, 0.0
hollow
material {
texture {
pigment {
rgb <1, 0.000000, 000000>
}
finish {
diffuse 0.6
brilliance 1.0
}
}
interior {
ior 1.3
}
}
translate <-0.5,0,0.5>
}

cone {
<0.0, -0.5, 0.0>, 0.5,  <0.0, 0.5, 0.0>, 0.0
hollow
material {
texture {
pigment {
rgb <1, 1.000000, 000000>
}
finish {
diffuse 0.6
brilliance 1.0
}
}
interior {
ior 1.3
}
}
translate <0.3,0,0.3>
}

cone {
<0.0, -0.5, 0.0>, 0.5,  <0.0, 0.5, 0.0>, 0.0
hollow
material {
texture {
pigment {
rgb <0, 1.000000, 000000>
}
finish {
diffuse 0.6
brilliance 1.0
}
}
interior {
ior 1.3
}
}
translate <0.4,0,-0.2>
}

cone {
<0.0, -0.5, 0.0>, 0.5,  <0.0, 0.5, 0.0>, 0.0
hollow
material {
texture {
pigment {
rgb <1, 1.000000, 1.000000>
}
finish {
diffuse 0.6
brilliance 1.0
}
}
interior {
ior 1.3
}
}
translate <0,0,-.7>
}
```

Posted in Analysis, Optics, Other, POV-Ray | Tagged: , , , | 2 Comments »

## What is the center line of a complex polygon? New approach.

Posted by smathermather on January 17, 2012

Today, I just show pictures of my process for filtering the medial axis of extraneous elements. I’ll let the pics tell the story for now:

Original Geometry, with ST_Segmentize run to densify the vertices

Medial Axis (Calculated outside PostGIS... at this time)

Density of vertex points in Medial Axis-- proportional to how closely it matches the direction of the local geometry

Medial Axis filtered based on density of vertices. Red segments are the removed ones, blue are the retained. It's not perfect, but it's very cheap computationally... .

scale axis transform 3

scale axis transform 2
scale axis transform 1

previous post

original

Posted in Analysis, Database, PostGIS, PostgreSQL, SQL | Tagged: , , , | 3 Comments »

## What is the center line of a complex polygon? (cont. again!)

Posted by smathermather on January 14, 2012

I didn’t remember that I had 4 previous post on this topic, but I’ve gotten a little obsessed with this problem– how to (with relatively little computational cost) kind the centerline of complex polygons, so we can extract flow lines from streams, label long complex polygons, easily generate the centerline parcel boundary line for river bound parcels, etc..  FYI, here are my 4 previous posts on this topic:

scale axis transform 2
scale axis transform 1

previous post

original

I’ll confess, I became somewhat enamored of the scale axis transform as a method for tackling this problem.  In short, the scale axis transform identifies the medial axis of the polygon using Amenta and Kolluri’s finite union of balls to identify the initial centerline of the polygon, and then tries to fix the deficiencies in the robustness of the medial axis.

The problem, of course, is the degree to which bumps in the shape result in extra little lines hanging off our medial axis.  This does not represent an intuitive medial axis.  This is not what we would (in this case) call the centerline of the stream.

The scale axis transform suggests expanding these balls by a multiplicative factor, and deriving the new axis from that scaled version:

This step can easily be done in PostGIS once you have the medial axis:

```
CREATE TABLE rr_expand_r3_12 AS

SELECT
ST_Buffer(b.the_geom, ST_Distance(b.the_geom, ST_Union(a.the_geom)) * 1.2 )
FROM medial_axis b, polyline a
GROUP BY <a href="http://a.id/" target="_blank">a.id</a>, b.the_geom
;

```

I will play a bit more with this, but the initial results I got around areas with islands (torus like topology) were disappointing.  Also, at least as I have this implemented in PostGIS, this is a computationally intensive process.  I think though that I have a good (enough) alternative, which will be the topic of my next post.

Posted in Analysis, Database, PostGIS, PostgreSQL, SQL | Tagged: , , , | 1 Comment »

## PostGIS Cartographic Effects– Cartoonify Nearly Coincident Lines

Posted by smathermather on January 10, 2012

In my previous post, a long 24-hours ago, I proposed some automatic modification of line for cartographic reasons. I had some flaws in my code. The points were over-rotated by 45 degrees. Can you spot why? Tip: it’s a basic trigonometric mistake. Here’s the corrected code (though there may be a better way):

```CREATE TABLE movedRoad AS
SELECT
ST_Translate(b.the_geom,
0.5 * ST_Distance(b.the_geom, ST_Union(a.the_geom)) * (cos(ST_Azimuth(b.the_geom, ST_ClosestPoint(ST_Union(a.the_geom), b.the_geom))) - 3.14159/4),
0.5 * ST_Distance(b.the_geom, ST_Union(a.the_geom)) * (sin(ST_Azimuth(b.the_geom, ST_ClosestPoint(ST_Union(a.the_geom), b.the_geom))) - 3.14159/4) )
AS the_geom
GROUP BY a.id, b.the_geom
;
```

An alternate approach is to only move those points that are too close, ala:

```CREATE TABLE movedRoadAgain AS
SELECT
ST_Translate(b.the_geom,
(100 - ST_Distance(b.the_geom, ST_Union(a.the_geom))) * (cos(ST_Azimuth(b.the_geom, ST_ClosestPoint(ST_Union(a.the_geom), b.the_geom))) - 3.14159/4) ,
(100 - ST_Distance(b.the_geom, ST_Union(a.the_geom)))* (sin(ST_Azimuth(b.the_geom, ST_ClosestPoint(ST_Union(a.the_geom), b.the_geom))) - 3.14159/4) ) AS the_geom
GROUP BY a.id, b.the_geom
HAVING ST_Distance(b.the_geom, ST_Union(a.the_geom)) < 100

UNION ALL

SELECT
b.the_geom AS the_geom
FROM bridle_points b, apt a
GROUP BY a.id, b.the_geom
HAVING ST_Distance(b.the_geom, ST_Union(a.the_geom)) >= 100
;
```

## PostGIS Cartographic Effects– Cartoonify Nearly Coincident Lines

Posted by smathermather on January 9, 2012

I’m still working on this query, but I thought I’d post what I’ve done so far. My intent is to produce scale-dependent exaggeration of the distances between quasi-parallel lines. The reason for this is so that lines such as street lines which are nearly coincident at a particular viewing scale can be spread from each other, much in the same way great cartography lies a little to display a relatively correct map. The nice thing about digital cartography is that as we zoom in, we can make the lie a smaller and smaller, until it is just barely a fib :).

My thought was to start with the simplest case– move one line away from another (stationary) line. First step is to find the closest point on the stationary line, determine the distance and angle to that point, and then translate the point outward by some scale factor along that angle. Don’t take my code as gospel, though, I think there’s some deeper logical error, but see what you think:

ST_ShortestLine gives us the shortest line from a given point to the road from which we are trying to offset. This is just to help visualize the problem:

```CREATE TABLE shortLine AS

SELECT
ST_ShortestLine(b.the_geom, ST_Union(a.the_geom)) AS the_geom
GROUP BY a.id, b.the_geom
;
```

Now, we try to translate the point outward at the same angle, by half again the distance to the closest point:

```CREATE TABLE movedRoad AS
SELECT
ST_Translate(b.the_geom,
0.5 * ST_Distance(b.the_geom, ST_Union(a.the_geom)) * cos(ST_Azimuth(b.the_geom, ST_ClosestPoint(ST_Union(a.the_geom), b.the_geom))),
0.5 * ST_Distance(b.the_geom, ST_Union(a.the_geom)) * sin(ST_Azimuth(b.the_geom, ST_ClosestPoint(ST_Union(a.the_geom), b.the_geom))) )
AS the_geom
GROUP BY a.id, b.the_geom
;
```

Resulting in something looking like this (red dashed line is the road we are moving away from, the dotted blue are the points making up the line we are moving, and the dotted brown are the newly transformed points, with ST_ShortestLine also visualized in thin brown):

Posted in Analysis, Cartography, Database, PostGIS, PostgreSQL, SQL, Trail Curation, Trails | Tagged: , , , , , , | 1 Comment »

## There’s an App for That, Alternate

Posted by smathermather on January 2, 2012

codinggeekette, or Sarah Dutkiewicz and her husband introduced me to Linux and Open Source many years ago.  Running counter-culture is the genius that is @sadukie, so now she’s a Microsoft MVP.  We’ll forgive her that, since she does crazy stuff like setting up .net environments on Linux and other fun stuff.  Anyway, Sarah has a post on phone apps she likes and has proposed in her post (quasi) simu-blogging on the topic.  I’ll bite.  I’ll focus mine on the limited mapping apps I have on my iPod, since my smartphone solution is an iPod touch and a Motorola flip phone.

App number one: Tiltmeter.

Tiltmeter is simply that– it logs at a specified interval the 2D axis tilt of your iPod or iPhone, and can show that visually or e-mail the log to you.  Oh, where were you Tiltmeter, when I was building camera arrays, putting bandpass filters on them, and launching them on a 12-foot balloon.

It was a low cost, few band, hyperspectral attempt at calculating red-edge effects for productivity estimations.  A very early attempt at something similar to Grassroots Mapping.  Let’s just say, a tiltmeter would have really made the georeferencing easier… .

App number B: Ride the City.

Ride the City (Desktop web interface)

Ride the city is a very nice point-to-point routing solution for a desktop or mobile client.  It’s a discreet app on the mobile side.  Other than not doing looped (recreational) mapping support, it’s a pretty cool application.  Great for commuter bicycling.

App Gamma: PDF Maps

PDF Maps is a tool from Avenza Systems for viewing Geo PDFs.  I’ll confess, since I have an iPod touch, and therefore no GPS, I haven’t run this one through the ringer (and there are other likely choices), but it could be a really useful tool in concert, e.g., with the Geo PDFs from the USGS.

Posted in Other | Tagged: | 1 Comment »