28 Sep 2008
The code for spatial predicates in PostGIS is largely dependent on the GEOS topology library, because doing predicate calculations in generality is hard, and GEOS already exists. However, moving geometry from PostGIS into GEOS format incurs a cost. And not all predicate algorithms are hard. Point-in-polygon tests, for example, are relatively easy.
So, for performance reasons, a couple years ago, Mark Leslie (at Refractions at the time) implemented a point-in-polygon test directly in the PostGIS library. Whenever ST_Contains(), ST_Intersects(), etc, are called, the geometries are first checked to see if they are points and polygons, and if they are, GEOS is avoided and the calculation is done in PostGIS.
But, why stop at one shortcut?
What if you are testing hundreds (or thousands) of points against one, or a small number, of polygons? Why iterate through all the segments of the polygon for every point, to carry out the test? By indexing the segments of the polygon, you can reduce the computational effort of doing a point-in-polygon test from O(NumberOfPolygonEdges) to O(log(NumberOfPolygonEdges)). However, for the index to be effective, you have to re-use it for each new point, not re-build it for every polygon/point pair. That means it has to be cached between function calls.
A bit over a year ago, Mark implemented a caching version of the point-in-polygon shortcut, with an indexing algorithm from Martin Davis, and that shortcut currently resides in the 1.3 release series. However, it has two drawbacks.
- First, it only works for POLYGON/POINT combinations, and most people working with polygons actually have MULTIPOLYGON/POINT (though their multi-polygons usually only have one member).
- Second, it leaks a lot of memory during processing, so the postgres process size can get very large while the computation is running (PostgreSQL retrieves the memory at the end, so no permanent damage is done).
Last week I took a couple days to delve deeply into what this shortcut was doing, and made the following improvements.
- Removed all the memory leaks, so the process size remains constant throughout the run.
- Added support for MULTIPOLYGON types.
- Improved the caching logic slightly, so that segment indexes are only built when a polygon has been seen two times in a row, otherwise using a standard non-indexed version of the point-in-polygon algorithm.
The point-in-polygon caching shortcut is extremely effective. Using the un-cached code, a spatial join of 8000 points to 80 polygons, where there are an average of 100 points per polygons, takes about 30 seconds on my workstation. With the caching segment indexes, the same join returns in 6 seconds.
The improved code is currently on trunk only, but I will back-port it into 1.3.X next week, so it will be available in the next points release.
Update: These changes have now been back-ported to 1.3 and will be generally available in 1.3.4 when it is released.
25 Sep 2008
I encourage all you lucky bastards folks heading to FOSS4G 2008 in Cape Town next week to upload your photos to Flickr using the tag foss4g2008. We got a great collection of photos from folks last year under the foss4g2007 tag.
24 Sep 2008
It looks like the 2008 edition of FOSS4G (going on next week!) in Cape Town, South Africa is going to crack 500 attendees. This is a great showing for an event which has had to overcome substantial distance from the FOSS4G centers of gravity in North America and Europe. And, supporting the theme of “An Option for Developing Countries”, a large proportion of registrants are South African. So lots of folks will be learning about open source geo-spatial next week, and taking it back to their communities.
As usual, I am suffering non-buyers-remorse as an exciting event that I won’t be attending approaches. However, with the family so wee (1.5 and 4), going away for well over a week just isn’t on right now.
23 Sep 2008
Don’t let the Winter Olympics in Vancouver steal all the thunder for 2010, make your own mark, by organizing FOSS4G, the 2010 edition!
Organizing an international conference is more fun than a barrel of spit and will mark you permanently, leaving you a quivering shell of a man it’s a great learning experience, too!
22 Sep 2008
I wrote previously about how using the Intel “Integrated Performance Primitives” (IPP) and the Intel compiler (ICC) helped to improve performance in a Mapserver raster rendering system by about 40%.
I was feeling pretty proud, until the client pointed out that the outputs, while fast were, um, sort of wrong:
What was going on here?!
Well, the problem looked to be in the JPEG writing, and Mapserver writes its JPEG via libgd by default, so I looked at the JPEG code in libgd for clues. There was little to find. LibJPEG is provides a very simple interface to writing out images, and the code in GD (gd_jpeg.c) was identical to the sample image compression code in the Intel IJG code base (cjpeg.c).
I only had two clues to go on: first, the artefacts would go away, if we created the JPEGs in a “progressive” (interlaced) mode; second, the Intel cjpeg utility could create correct JPEGs in non-interlaced mode.
After taking a brain-break of a couple days, I found the problem while grepping the code and watching TV. Setting up a JPEG compressor requires that you provide the compressor with a buffer-fill function of your own design, so you can provide input from all kinds of sources. The implementation of the buffer-fill routine in the Intel example code had a switch in it, that checked if the mode was progressive or not!
In the case of non-progressive JPEG, the Intel function would allow the compressor to only partially read the buffer if it wanted. Then it would move the un-read bytes to the front of the buffer, and fill in the rest with new bytes. The standard JPEG library expects the compressor to always read the full buffer.
Patching GD to use this partial buffer logic, and JPEG compression worked, I was getting what looked like good ouputs! So, I put the patch in place and told the client and, um, there were still problems:
Gah! The artefacts were smaller, but they were still there. However, having my previous solution in hand, it was easy to find this bug. The artefacts were in a nice blocky patterns, just like the internal TIFF tile structure, which pointed to the JPEG decompressor having the same buffer behavior as the compressor.
So this time I opened up the GDAL library (that Mapserver uses to read the TIFFs) and found the JPEG reading code in the TIFF driver. Sure enough, the Intel decompressor only partially read the buffer, while the standard decompressor was expected to read the whole thing every time.
A small patch to the GDAL TIFF driver, and another patch to the GDAL JPEG driver (for completeness) and now we finally had IPP-accellerated JPEG input and output without artefacts.
Moral of the story: Although the Intel IJG is API-compatible with the standard libjpeg, the slightly different expected behavior of the empty_output_buffer and fill_input_buffer functions means that most existing code will require slight modifications before it can use the Intel libraries – they cannot just be dropped into place as a binary-only update.
