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
}

19th Century Camera Calibration for Remote Sensing in PovRay (or Porro-Koppe Principle in Virtual)

A complete analytic solution for geometric distortions in remote sensing (ahem, ignoring atmosphere, of course)

Ok, so here is another thought experiment. This time, it will take a few months before I have time to write the code for this, but maybe the thought experiment will inspire someone. A major problem of photogrammetry, remote sensing, and computer vision is correction of lens distortion. Thanks to almost a century of working with this problem, and recent developments in computer vision, your average geek can now calibrate her camera with freely available code. see the Camera Calibration Toolbox for code written in Matlab– there are also links to other free camera calibration projects.

I have a real interest in using cheap off-the-shelf cameras for solving small remote sensing problems, and it would be nice to be able to correct distortions in the images with my favorite computational tool, PovRay.

So, how to do this? Well, let’s not think of this as an empirical problem like most of the corrections of today. The Camera Calibration Toolbox, for example, requires image inputs from the camera to correct for the distortion. What if, instead, we arrive at a completely analytic solution using the known information about the array of lenses, their geometry and index of refraction to correct for distortion? What if we created a virtual array of lenses identical to our camera, and pass our image back through the lenses to cancel the lens distortion effects? This would be the virtual version of an late 19th century technique discussed by Clarke and Fryer, 1998:

For mapping applications the earliest solutions to the problems associated with large radial lens distortions were by direct optical correction whereby the image was re-projected through the camera and lens system which had captured it. This system was termed the Porro-Koppe Principle after the scientists who perfected it in the latter part of the 19th century. In this manner the geometric distortions in the image were canceled.

How do we build such a system? Well, we start with a photon scene, ala Henrik Wann Jensen, and take advantage of the constructive geometry lens set ala Don Barron, and project our image back through a virtual version of our lens set, to be captured on a flat surface with an orthographic camera in PovRay. But here’s an opportunity– If our original scene captured by the camera was not flat, and we knew its three dimensional properties from other information (stereo pair or lidar), we could project the scene back on to it’s 3D geometry, then capture with an orthographic camera, and then we could correct for terrain distortion and camera distortion all in one go. Fun stuff huh– all implemented in free software, a complete analytic solution for geometric distortions in remote sensing.

Just one problem– I’m not quite sure how to project an image in PovRay. Oh, I could do it brute force and create a square transparent pane of color for each pixel (which I may end up doing), but if anyone has a better idea, I’m open to it.

T.A. Clarke and J.G. Fryer, Photogrammetric Record, 16(91): 51-66, April 1998 found at: http://www.vision.caltech.edu/bouguetj/calib_doc/htmls/ref.html)

POV-Ray for viewsheds code

global_settings {
max_trace_level 1
}

#include “shapes.inc”
#include “colors.inc”
#include “textures.inc”

height_field {
png “c:\temp\povray\n2225625.png” // 16-bit integer digital elevation
water_level 0.0

texture {
pigment {
image_map { png “C:\temp\povray\n2225625_drape.png” map_type 0 interpolate 2 once }  // Aerial
rotate x*90
}
}
scale <5000, 65536, 5000>  // Scale to real world size
translate <2224868,0,625008>  // Translate to location in Ohio state plane
}

/// Orthographic Camera
camera {
orthographic
location <2227368.0, 6000, 627508>
look_at <2227368.0, 0, 627508>
right 5000*x //x size of view  (feet)
up 5000*y // y size of view  (feet)
}

// Observer Points

light_source { <2227885, 720, 625464> color Cyan }
light_source { <2228247, 727, 626277> color Yellow }
light_source { <2228091, 729, 626761> color Magenta }

Follow on to POV-Ray for GIS Analyses–metametacode, no code yet, sorry!

Ok, here’s the basic idea– we have an orthographic scene with a height_field object scaled to real world units. The observers in the viewshed become points of light, thus for each observer, we “light” the areas visible from each observer, render, and boom, viewshed created. In addition, if we have three or fewer view points, we can render them in, say cyan, yellow, and magenta, and thus tell which observer points can view a given location. But I digress. I keep leaving the actual code at work so I’ll start out by giving you the

Povray Viewshed Calculation Meta-Metacode:

In a Povray height_field, input elevation images are scaled 1 unit in all 3 dimensions. In our case, the input image is 5000 feet on a side in real world units. So, we multiply x and y dimension by 5000 to scale it. Since the unit values in the Z-direction are scaled to one as well, and they were input as 16-bit, that means that povray effectively divides the data by 65536, so we multiply by 65536 to get back to elevation in feet. Now our x, y, and z dimensions are scaled appropriately to each other.

Since we’re working with geographic data, well put it back in the real world (just in case we want to put other data in our system as well), in this case the Ohio State Plane North NAD83 HARN (feet) projection. So, translate to location in real space (the state plane centroid of the image) and set an orthographic camera above it with an image plane equal to the x and y dimensions. The orthographic camera prevents distortion from perspective, so we retain a flat map. (Side note: I’m contemplating using an orthographic camera and DEM to correct for terrain distortions in uncorrected imagery (and maybe shading as well), but that will take more thought, and may become another post).

Ok. Now to calculate viewshed from a point, place a point light in the scene at the location in question. Render, and boom, we have a viewshed from a point. Want a few points, add a few more point lights. Want to constrain the effect in the x, y, or z direction, add some opaque baffles (I haven’t done this in my code yet), or directional lights.

Finally, just for kicks, we can drape an aerial image on top, and our viewshed lights will only light the parts of interest, something I’ve never seen done in a GIS. But again, (see previous post) the main point is to be able to simulate vegetation, buildings and other comprehensive aspects of the scene, not just a digital terrain (elevation) model. If that’s all we wanted, a GIS would suffice. So, next step, add buildings and trees. Long term, I’d like to do this in Physically Based Renderer (PBRT) or (now that I know it exists…) LuxRender. Heck with physically based rendering, we can do some inverse modeling of physical/chemical characteristics, or WiFi placement optimizations. But again, I digress… .

POV-Ray for GIS Analyses

Well, I tried something that so far has been a real success. I did a viewshed analysis in POV-Ray. Nothing significant in that analysis– it’s just a topography based analysis, just like one might do in GRASS or ArcGIS. But what is interesting, is now I can add buildings (also not hard in a GIS) and trees (which can be quite hard to add). Suddenly the viewshed including trees need not be a solid wall of trees, but something partially penetrable and complex. I’ll do this in part with some new lidar data I have… and update here. In the next couple days, I’ll post the code and results I have so far. The code is very simple and not terribly flexible yet– just a test of concept, but we’ll see where we can go with this.

PostGIS Onesie

Nothing more appropriate for the baby elephant than to be displayed on a onesie.

Well, I have to put it out there: is this an original, the one and only PostGIS onesie? My son was born recently, and as he was sleeping quietly, I surfed to the PostGIS website. One look at that baby elephant icon, and another to our newborn, and I realized I had the perfect inspiration, so time to go make my own onesie, and start a FOSS GIS blog.

FYI, I’m working on a PostGIS/GeoServer/OpenLayers set of services for my day job, but this probably doesn’t count as work.