Fast Calculation of Voronoi Polygons in PovRay– Applied

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

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>
}

Landscape Position: Conclusion? (part 2)

From earlier post:

“I’ve managed to pilot most of a fast high resolution landscape position workflow with PovRay as my magic tool. The final steps I hope to pipe through PostGIS Raster. In the meantime a screenshot and description: blues are riparian, raw ocre, etc upland categories, grey is mostly flat lake plain and mid slopes, all derived from just a high res DEM input (no hydro lines as input, so it works on areas where you don’t know where the streams may be). There will be more categories in the final, in the meantime, welcome to December.”

I didn’t use PostGIS Raster, but the incredibly useful gdal_calc.py to finish out my landscape position calculations.  I will endeavor to post my code to github soon, but the parts and pieces include using gdal and PovRay.  PovRay helps with the sampling (really!) of nearby neighbors in the raster DEM, and gdal does the averaging and differencing of those to get relative landscape position.  I spent some time yesterday scratching my head over how to show all the new landscape position information on a readable and useable map, and after discussion with collegues, decided to use it to divide the world into two categories– riparian & lowland + headwaters & upland (not reflected yet in the labels).  To find out more about landscape position, follow this link, or better yet this one.  (BTW, the green is park land, blue is riparian/lowland/stream networks, purple is the basin boundary).

Riparian map based on landscape position, calculated using PovRay and GDAL.

Landscape Position: Conclusion?

I’ve managed to pilot most of a fast high resolution landscape position workflow with PovRay as my magic tool. The final steps I hope to pipe through PostGIS Raster. In the meantime a screenshot and description: blues are riparian, raw ocre, etc upland categories, grey is mostly flat lake plain and mid slopes, all derived from just a high res DEM input (no hydro lines as input, so it works on areas where you don’t know where the streams may be). There will be more categories in the final, in the meantime, welcome to December.

Landscape Position and McNab Indices (cont.)

I typed that last one too quickly– too many typos, but my wife says I’m not supposed to revise blogs, but move on… .

So, for clarity, let’s talk a little more about McNab indices.  Field-derived McNab indices are a measure of average angle from the observer to the horizon (mesoscale landform index), or from the observer to another field person a set distance away, e.g. 10m or 18m, or whatever the plot size or settled upon (minor landform index).  Typically these are done in the field at a discreet number of directions, e.g. 4 cardinal directions, or 4 cardinal plus 4 ordinal (NE,NW SE, SW, SE) directions. The landform image in this post and the last is a calculation of mesoscale landform, which is harder to do in a classic GIS than minor landform (I’ll have a follow-up post on minor landform, probably using ArcGIS).

To calculate the values for mesoscale landform computationally, we require calculating the angle to the horizon for a certain number of directions for each point in the landscape.  This has the potential to be computationally intensive and complicated.  If, however, we conceive of the problem differently, as a 3D calculation we can perform in PovRay, arriving at the answer is simplified by the coding already done by the PovRay programmers.

Essentially, McNab mesoscale indices are a field proxy for the steradian of a site.

What does that mean?  Well, essentially, a map of steradians is a proxy for the shading of a white uneven surface on a cloudy day– or how much diffuse light is available to a given site.  Used in combination with site aspect, this is enough information to determine most of of the light conditions of a site, which is why McNab indices in combination with other factors are a good predictor of site productivity, and correlated with different plant communities across the landscape.

How does an uneven white surface shade on a cloudy day? Steeper areas with more of the sky occluded are darker while wide open spaces, like the bottom of river valleys and the tops of ridges and plateaus are brighter.  If you want to witness this effect, look to snow on the ground on a cloudy day (and what a great winter to do it).  (The only difference with snow is subsurface scattering which evens out the light quite a bit.)  You’ll notice the divots in the snow to be darker than the peaks, and the edges of the divots to be darkest of all.

The question then, is how to we compute this within PovRay?  We could use radiance as a global illumination model, but the calculation of inter-reflectivity that is at the core of radiance, while an important real factor in the landscape, would fail to replicate the original mesoscale landform index.  Instead, we set up a series of lights in a dome to illuminate the whole sky sphere, blur them a bit, and call it a day, a technique developed for HDRI rendering by Ben Weston, whose code I borrowed heavily.  The more lights the element of the landscape is exposed to the brighter, and vice versa, essentially making the brightness of an element proportional to the exposure to the sky sphere.  Unlike Ben, I used a simple white sphere to get even lighting.

Rendered at 16-bit resolution, we have a possible range of values from 0-65535.  Let’s assume that a linear transformation of these values will result in values representing the proportion of the sky sphere.  From there, transformation to steradians and then to solid angles in degrees is trivial.  Once it’s solid angles in degrees, it represents the same kind of value as a mesoscale landform index would give us (more later…).

Povray Viewshed with Trees (finally) (cont.)

Ok, your average terrain-only based viewshed (view is from a road to the southeast, viewshed is in cyan):

Note that based on these estimates, we should be able to see most of the valley walls from this little slice of road.  I don’t buy that.  That section of road is pretty closed in with trees.  Let’s add them:

As you may see, just a little cyan in the southeast mostly along the road.  Here’s where you can see the trees occluding our view of the viewshed a bit.  But, it’s a first good stab, even with such boring results… .  I’ll have to get that section of road in here that peeks out over the gorge.  But first, I’ll have to see if I can get tree stem locations from the lidar, so we can trade up for some veracity in our analysis.  Stay tuned.

Povray Viewshed with Trees (finally)

Povray Viewshed with Trees: still some tweaking to do, especially tweaking that represents the real location and size of trees on the landscape, but if we cannot yet have veracity, best to have some verisimilitude.  Also, long term, it’d be nice if the light sources would be effected by the trees, but not the camera.  The trees occlude the ability of the camera to see what the full viewshed is.  Not clear what I mean?  Give me a few days, something nagging the back of my brain says the solution is trivial.  In the mean time, here’s the code:

Now for the include files:

light_posits.inc:

tree_coords5.inc (can you tell I’m not using version control):

Images to come of the output.  It is rendering as we speak.  I didn’t like the Povray 3.7 beta versions I ran at work on a much faster computer, so now I’m waiting for it to finish on my wife’s laptop… .

Povray Landscape with Trees

Just a short teaser post until I remember to bring the code home with me, but I’ve placed trees within the bounds of a canopy layer determined from some county LiDAR dataset.  The tree locations, heights, and rotation are set using pseudo-random numbers, which looks a lot better than a monoculture of the same tree.  Eventually, estimating stem locations and canopy heights to do this would be preferable, but one step at a time… .  We’ll first look for a reasonable result, and follow with an accurate one.  For now I get a pretty realistic canopy over a gorge area for a creek, so below are the first, middle-ish, and last frames in an animation I did of that canopy.

Close in:

rendering partway up…

(Check out the creek running through the gorge)

And the final frame:

POV-Ray for viewsheds (with trees?)

In some of my first posts, I discussed the possibility of using povray for viewshed analyses, since it is a more flexible tool, and can better handle complex analyses, like terrain + vegetation, something which most GIS tools cannot.  In the end though, I just produced a simple terrain based viewshed analysis.  Now, I’m getting ready to go deeper.  Now, I’m ready to put some 28,000 trees on our viewshed terrain.

So, why not just use a digital surface model instead (an elevation model that shows us the top of all vegetation, buildings, ground, etc.)?  Well, trees, especially if there is a single row of them, are not going to necessarily block all views, so a digital elevation (or terrain) model in conjunction with an estimate of tree locations is a far better solution for accurate estimations of viewsheds.

So how do we get there?  I started with a dense LiDAR dataset with 3 returns.  I performed my analysis only on points that were either return 2 or 3, which tells me I very likely am dealing with canopy vegetation.  Then I take a moving average using (yes, I know) ArcGIS’ Spatial Analyst Point Statistics tool to determine canopy heights.

Now, to make things easy, rather than figuring out where the trees really are, I’m going to then generate random points approximately 30 feet apart, and use the elevation from my raster from the Point Statistics tool to determine the heights of my canopy.

To save on parsing time in povray, I’ll generate 28,000 non-unique trees across the landscape and see how things look.

First the code for the placement of the 28,000 random trees (written in BASH):

Now we have some random tree locations.  We make an include file of it (called “tree_coords2.inc”) that looks like this:

Now we need some trees.  I don’t want my machine to use too much memory, so a good tree “include” would be the ticket.  Pov-Tree seems the ticket here, and I’ll even use one of the pre-built trees, the Linden.

I forgot before to include the pov file for rendering:

So here’s the canopy close up:

Which doesn’t look too bad… .

And now far out:

So now we have a start for putting trees on our landscape (coming next…).

Follow-up (#1) for Camera Calibration in POV-Ray (Porro-Koppe Principle in Virtual)

Ok, so projecting an image isn’t hard at all in POV-Ray.  I’m not sure what I was thinking.  In my first post on camera calibration in POV-Ray, I suggested that we could analytically solve for lens and topographic distortion in POV-Ray.  I haven’t gotten far as I don’t have a work project to take advantage of this yet, but I projected a false color infrared image from a Landsat image over Argentina onto the walls of virtual building.  The walls would be replaced with topography, the camera with an orthographic camera, and some lenses would be in-between, but at least it is a start.  It is basically just a modified version of some example code from Boris van Schooten.  Just so long as a photon scene does something similar, the rest of the problem shouldn’t be too hard to solve.  With my schedule as it is right now, about 3 hours of coding is probably about 3 months though… .

Projected Satellite Image

union {
polygon {5, <0,0,0>, <0,1,0>,<1,1,0>,<1,0,0>,<0,0,0>}
disc {<0.5, 1, 0>, <0, 0, 1>, 0.5}
pigment {
image_map {jpeg “landsat-argentina.jpg” interpolate 2 filter all 1.0 }
scale <1,1.5,1>
}
finish{ambient 1.0}
translate x*-0.5
scale <20,20,1>

translate z*13.9
}