No Image– a very simple solution for calculating viewshed

Not much to say here. I had a bit of inspiration, remembering a Christoph Hormann tutorial http://www.imagico.de/pov/icons.html on making transparent icons in povray and realized that what he was describing was the final piece in the povray for viewsheds milleau: no_image option for an object forces it to interact with the scene in every way, but to be seen by the camera. I’ll post the revised code coming up, but it’s pretty simple. Add the following option inside the code to the tree object:
no_image .

http://www.povray.org/documentation/view/3.6.1/328/

Another Blog

Just as I’ve started getting good traffic on this blog, a collegue and I have started another one called Anantasesha, largely on the GIS stack we’re putting together for internal use, but also may include some interesting other bits of code and musings. Check it out:
http://anantasesha.wordpress.com

POV-ray for GIS Visualization

This is a “wouldn’t it be nice” post– i.e. a I-don’t-have-time-to-write-it-and-everything-else-on-my-list post. Wouldn’t it be nice if in addition to returning GeoRSS, WFS, OGC, KML, etc. GeoServer returned a POV-Ray file? It could be returned as part of a simple GET post like this:

http://localhost:8080/geoserver/wms?bbox=2100506,575982,2290532,739026&styles=&Format=application/povray&request=GetMap&version=1.1.1&layers=base:contours&width=800&height=644&srs=EPSG:102722&camera=2195519,657504,1000&look_at=2195029,657004,500

Vector symbolization could be a POV-Ray approximation of the WFS render styles using coordinates and POV primitives. Also, a PHP script could consume this return, and run a command line version of POV to render the file in a webserver implementation. This second part I’ll probably actually write with just a DEM and image overlay, but it would be cool to have the vector symbolization as well as part of another GeoServer protocol.  Of course, I recognize this is not a OGC protocol… .  Alternatively, a PHP script to consume WKT and convert it to a POV-Ray file might be easier than incorporating it in GeoServer.

PostGIS Gripe—Limits to Postgre’s B-tree indexing

Is this a gripe about PostGIS or PostgreSQL?  Hard to say.  At my office, we are migrating to a PostGIS—GeoServer—Client-of-your-choice stack.  As usual, I have to try the hard stuff first.  So, rather than keeping the 2ft contour lines for our 7 county area in shapefiles and doing some indexing on them there, I’m attempting to load them into PostGIS all in a single table.  Make any sense?  Maybe not.  Honestly, I really don’t know.  In doing so though, PostGIS has done pretty well.  I loaded all 9.6 GB of 2ft contours from one county into the database, served it up directly in uDig, and then through GeoServer to OpenLayers and Google Earth.  For viewing the whole county, it was slow, but I’d never let a user view the 2ft contour dataset all in one go.  Zoomed in, however, it did pretty well, in fact, I’d say considering there was no caching turned on for OpenLayers, it did splendidly.

First, how I loaded it:

I used shp2pgsql.  My contours were in a bunch of little chunks derived from a state DEM dataset, so I iterated through:

Create the table:

shp2pgsql -s 102722 N2110610 base.contours > N2110610.sql

psql -U postgres -d CM -f N2110610.sql
Then iterate through the appends (note the -a flag)

shp2pgsql -s 102722 N2110615 base.contours -a > N2110615.sql
psql -U postgres -d CM -f N2110615.sql …

we can script this—I have to confess to having done this in excel with worksheet functions since my Bash shell scripting skills help me none in the Windows world, but a friend told me about FOR…IN…DO in batch scripting, so next time… .

So, loading was easy, and all worked well.  Viewing was fine.  What is my gripe?  B-tree Indexing.  What if I wanted to ensure I had only unique geometric records in my database?  I might initiate the table like this instead:

CREATE TABLE “base”.”contours” (gid serial PRIMARY KEY,
“objectid” int8,
“type” int8,
“elevation” smallint,
“shape_len” numeric);
SELECT AddGeometryColumn(’base’,'contours’,'the_geom’,'102722′,’MULTILINESTRING’,2);
alter table base.contours add constraint unique_geom unique (the_geom);

Adding a unique constraint to the geometry creates an implicit B-tree index, and therein lies the problem– B-tree indices are limited in size to around 2700 bytes.  Our index exceeds that manyfold.  According to the Postgre documentation “Currently, only B-tree indexes can be declared unique.”

This probably is fine for most datasets, but it really limits the use of a uniqueness constraint on geometry columns of large size.  Since this dataset won’t require updates, my work-around (I hope) will be to load the data in the database, perform a select distinct and dump it into another table to clean up any possible duplicates, but it sure would be nice if I could constrain the table in the first place… .

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)