Boost.Geometry (and more)

SQLite, Ado.Net, SpatiaLite

2011-12-04T14:16:00.010+01:00

SQLite, Ado.Net, SpatiaLite

SQLite is a small but great database. SpatiaLite is its spatial extension. ADO.NET is the method to access databases in .NET.

This blog shows how to combine them, because I did found all information on the web but not a complete and up-to-date story.

SQLite and Ado.Net

An introduction (up and running in 3 minutes) is here. It is good but not up-to-date. In the meantime the ADO.NET driver is transferred from sourceforge/phxsoftware to sqlite.org, which is of course a good sign.

There is only one piece of software which you really need and that is the ADO.NET driver. Furthermore there are some optional (but very useful tools) which I describe below.

ADO.NET driver

The current ADO.NET driver page is here, and more specificly you can download the drivers from this page. Download for example the sqlite-netFx40-setup-bundle-x86-2010-1.0.77.0.exe (10.25 MiB) (at least I did) and install it. The installation folder will contain the driver, and also a manual and a sample database (northwindEF.db).

Setting up the project

Create a new project (or use an existing) in C# Express or Visual Studio (I'm using 2010). I created a Console Project.

Most important is of course to add a reference. Typically you need the System.Data.SQLite.dll from the installation folder (typically "c:\Program Files (x86)\System.Data.SQLite\2010\bin" ). You can either reference it from there or copy it to your project folder and reference it.

Writing the code.

See also the original 3 minutes introduction for some code. Here a different sample is presented. I like complete code and I also include some extra optional pieces which I needed (creating a new database from scratch, and turning on foreign keys).

So a complete program looks like:

// Copyright (c) 2011 Barend Gehrels 
using System; 
using System.Data.Common; 
using System.Data.SQLite; 

namespace SQLiteAndAdoDotNet 
{ 
    class Program 
    { 
        const string mydb = @"c:\temp\mydatabase.db"; 

        static void Main(string[] args) 
        { 
            bool isNew = false;

            if (!System.IO.File.Exists(mydb)) 
            { 
                // Create a new SQLite database if it does not exist
                SQLiteConnection.CreateFile(mydb); 
                isNew = true; 
            } 

            using (SQLiteConnection connection = new SQLiteConnection(@"Data Source=" + mydb)) 
            { 
                connection.Open(); 

                // Recommended action for any SQLite database
                ExecuteStatement("PRAGMA foreign_keys = ON", connection); 

                if (isNew) 
                { 
                    // Create a table for a new database 
                    ExecuteStatement(@" 
                        create table mytable 
                            ( 
                                id integer primary key autoincrement, 
                                myfield text, 
                                mydate integer -- dates can be stored as either text, real or integer 
                            ) 
                        ", connection); 
                } 

                // Insert something 
                ExecuteStatement("insert into mytable(myfield,mydate) values('my row',date('now'))", 
                    connection); 

                // Retrieve something 
                using (SQLiteCommand command = new SQLiteCommand("select *,datetime(mydate) as dt_ok from mytable", 
                            connection)) 
                { 
                    using (DbDataReader reader = command.ExecuteReader())
                    { 
                        while (reader.Read()) 
                        { 
                            System.Console.WriteLine(String.Format( 
                                "Row id={0} field={1} date={2}, dt={3}", 
                                reader["id"], reader["myfield"], 
                                reader["mydate"], reader["dt_ok"])); 
                        } 
                    } 
                } 
            } 
            System.Console.WriteLine("Done"); 
        } 

        private static void ExecuteStatement(string statement, SQLiteConnection connection)
        { 
            using (SQLiteCommand command = new SQLiteCommand(statement, connection))
            { 
                command.ExecuteNonQuery(); 
            } 
        } 
    } 
}

I factored out one method to avoid repetitions. In another blog I'll show you a different approach.

The foreign key pragma is of course not necessary in this short sample. But in a "real" database, with a datamodel and foreign keys, you will need it. Otherwise cascaded deletes will fail unexpectedly!

One caveat

Beware for datetime fields, stored as real or integer! SQLite stores date types in any field, which is reasonable. But if you select such a date(time) field in the normal way you will get only the year back! So alas you have to cast it with the SQLite datetime function (described here). Dates stored as text do not have this caveat.

Useful tools

I like the SQLite Studio which can be downloaded freely from here. You can hardly without such a tool and this one is really cool. If you are coming from SQL Server you might also like a converter, and this one is good, free, and with source.

Using SpatialLite

Using SpatialLite is pretty much the same, from ADO.NET perspective. But of course you have to download some extra libraries. The SpatialLite library of course, but you also need at least two extra libraries. So here we go (see also StackOverflow):

download the SpatiaLite DLL from this page (libspatialite-win-x86-2.3.1.zip)
download proj-win as well from the same page (proj-win-x86-4.6.1.zip)
download iconv as well from the same page (libiconv-win-x86-1.9.2.zip)
download geos as well from the same page (geos-win-x86-3.1.1.zip) (I added geos later, after the comment of rmales, for me it was not necessary because it was found in my path, via PostGIS.)

But them all into your binary folder (Debug and Release). For our test project we also need a spatial database. You can create one from a shapefile (using the spatialite-gui from the same page) or download a test database from SpatialLite, here. I downloaded and extracted world.7z because I like to see world countries (see my other blogs).

Then, before you have to use it, you have to load the extension. That is pretty easy with an SQL statement, select load_extension('libspatialite-1.dll'). If it finds this DLL and the other two, it will run fine (and otherwise it will mention which is missing). Alas I could not get it working in sqlitestudio.

Then below is the program I used, and it runs fine. It is quite similar to the one above but I prefer to list it completely again. In next blog(s) we will do other interesting and maybe surprising things with it.

// Copyright (c) 2011 Barend Gehrels 
using System; 
using System.Data.Common; 
using System.Data.SQLite; 

namespace SQLiteAndAdoDotNet 
{ 
    class Program 
    { 
        const string mydb = @"c:\bdata\blogs\data\world.sqlite"; 

        static void Main(string[] args) 
        { 
            using (SQLiteConnection connection = new SQLiteConnection(@"Data Source=" + mydb)) 
            { 
                connection.Open(); 

                // Recommended action for any SQLite database
                ExecuteStatement("PRAGMA foreign_keys = ON", connection); 

                // Required action for any SpatiaLite database
                ExecuteStatement("select load_extension('libspatialite-1.dll')", connection); 

                // Retrieve some cities in Europe 
                using (SQLiteCommand command = new SQLiteCommand(
                        @" 
                            select wup_aggl,AsText(Geometry) as wkt,
                                X(Geometry) as x,Y(Geometry) as y
                            from Cities 
                            where MBRContains(BuildMBR(-2,40,10,80),Geometry) = 1 
                        ", connection)) 
                { 
                    using (DbDataReader reader = command.ExecuteReader())
                    { 
                        while (reader.Read()) 
                        { 
                            System.Console.WriteLine(String.Format( 
                                "City: {0} location: {1}, {2} wkt: {3}", 
                                reader["wup_aggl"], reader["x"], reader["y"], reader["wkt"])); 
                        } 
                    } 
                } 
            } 
            System.Console.WriteLine("Done"); 
        } 

        private static void ExecuteStatement(string statement, SQLiteConnection connection)
        { 
            using (SQLiteCommand command = new SQLiteCommand(statement, connection))
            { 
                command.ExecuteNonQuery(); 
            } 
        } 

    } 
}

Useful tools

Of course you can use Quantum GIS to perfectly visualize your SpatiaLite database. The SpatiaLite page also contains many useful links.

Summary

Above two complete programs to use SQLite and/or SpatiaLite in your .NET application. This is the C# version, VB.NET will work either.

Linestring/polygon intersection

2011-11-19T18:14:00.000+01:00

Linestring/polygon intersection

After a long while a blog again, a short one this time. I'm currently implementing intersections linestring/polygon for Boost.Geometry. Well, the basics (calculating intersection points) are already there for a long time, but for this combination the intersection points should be followed in another way and the correct pieces should be outputted.

When I do these things I normally check the results with both PostGIS and SQL Server. However, I'm encountering something weird in both of them.

The testcase

My testcase is looking like this:

SQL Server

In SQL Server the intersection is correctly done, but the linestring is reversed! I'm using this query:

with viewy as
(
   select
      geometry::STGeomFromText('POLYGON((1 1,1 3,3 3,3 1,1 1))', 0) as p,
      geometry::STGeomFromText('LINESTRING(2 2,1 2,1 3,2 3)', 0) as q
)
select
   p.STIntersection(q).STLength() as len,
   p.STIntersection(q).STAsText() as wkt
from viewy;

and this is my output:

3 LINESTRING (2 3, 1 3, 1 2, 2 2)

You see, the linestring is reversed! It is not only reversed with this linestring, but with all linestrings. Even if they are, for example, complete inside the polygon. If I reverse the polygon (from clockwise to counter clockwise), the output is still reversed. Mysterious...

PostGIS

Also a surprise in PostGIS.The query is looking there as following:

with viewy as
(
   select
      ST_GeomFromText('POLYGON((1 1,1 3,3 3,3 1,1 1))') as p,
      ST_GeomFromText('LINESTRING(2 2,1 2,1 3,2 3)') as q
)
select 
   ST_Length(ST_Intersection(p, q)) as len,
   ST_AsText(ST_Intersection(p, q)) as wkt 
from viewy;

and this is here the result:

3;"MULTILINESTRING((1 2,1 3),(1 3,2 3),(2 2,1 2))"

I'm getting a multilinestring back! With respect to contents, it is OK. But it is surprising...

OGC

I did not look it up but I don't think this specific behaviour is specified. So yes, all implementations will be correct but...

Boost.Geometry

At the moment of writing, technically more difficult cases as this is one are not completely finished. But I can already tell that the linestring will not be reversed and that the output will consist of only one linestring...

Spherical Side Formula

2011-06-24T23:56:00.006+02:00

Spherical Side Formula

This blog describes an algorithm and a formula to determine the side of one point relative to a segment (two other points), all three points located on a sphere.

The side informarmation is necessary for Boost.Geometry. It is used, among others, in the within algorithm, to check if a point is inside a polygon (point in polygon). With the spherical version of side, it can be checked if a point on a sphere (the earth) is inside a spherical polygon (e.g. a polygon measured in latitude, longitude).

It is also known as orientation (see also here and here)

The Cartesian side is shown below. Point p is on the right side of segment p1-p2:

A visualization of ths spherical side is shown below, also here, point p is right of segment p1-p2.

Note that if the cartesian side formula would be used on lat-long coordinates, this specific point would result on the left side.

Former approach

Boost.Geometry's former approach for spherical side used the useful aviation formulae of Williams, listed here. That did work but I was not completely satisfied about it:

it takes two * (3 cos, 3 sin, 1 atan) for the course, and then another 2 sin, and 1 asin, so in total 17 goniometric functions
it is non defined at the poles
I don't understand it fully, it seemed to me that there should be an easier way

An obvious alternative is translate to 3D cartesian coordinates and why not? After working this out a bit, it appeared that it indeed requires less goniometric functions. These goniometric functions are less performant and less precise, so if they can be avoided, the better. Also the new approach is also relatively easy, from mathematical perspective.

New approach

The new approach also uses sin and cos, but less: for each of the points: 2 sin and 2 cos, in total 12 goniometric functions. Less then the former 17. More important, by design of Boost.Geometry, these might be precalculated such that for a spherical polygon it is not necessary to recalculate this again and again.

Besides that, the algorithm is easier to follow.

The algorithm

Let us calculate a plane through the two points of the segment, and the center of the sphere. This is the plane also contains the Great Circle. Let us then calculate at which side of the plane the point is located, either left of the plane, or right of the plane, or on the plane itself. The side with respect to the plane is also the side of the point with respect to the segment on the sphere.

To calculate a plane, all three points are transformed from spherical coordinates (longitude, latitude) to cartesian coordinates (x, y, z). Described a.o. here with clear images. We can use the unit sphere, it is not necessary to take the radius of the Earth (or another sphere or sphere-like object) into account. Amsteram, for example, then lies somewhere here: (0.614161755 0.042946374 0.788010754).

Converting can conveniently be done using Boost.Geometry's transform algorithm (if p1 is stored in cs::spherical_equatorial):

        typedef model::point<coordinate_type, 3, cs::cartesian> point3d;
        point3d c1, c2, c3;
        transform(p1, c1);
        transform(p2, c2);
        transform(p, c3);

Then we define a plane through the center of the Earth (0,0,0) and the two points. We use the generalized plane equation for that (a x + b y + c z + d = 0), described e.g. here. This is done by taking the cross product of the two vectors (a vector is here: a mathematical vector). These two vectors are normally called v and w. Because we use (0,0,0) the points are the vectors, we don't have to subtract anything to create a vector. So v is equivalent to c1, w is equivalent to c2.

The cross product is the normal vector, called n, normal to the plane, so perpendicular to the great circle of the two points. So we can take the cross product:

vector_type const n = geometry::cross_product(v, w);

We now have our generalized plane equation. To make that explicit, we can use the conventional symbols a, b, c:

        coordinate_type const a = get<0>(n);
        coordinate_type const b = get<1>(n);
        coordinate_type const c = get<2>(n);

For the generalized plane equation, we should calculate d by subsituting one of the points going through the plane. Let us take (0,0,0)! By taking 0,0,0, all terms a x, b y, c z will be zero, so d is zero. Therefore we don't have to calculate d.

double const d = 0;

The generalized plane equation has the property that, if you substitute another point, it is either on the plane (then the result = 0), or at one side of the plane (then the result is positive), or at another side of the plane (then the result is negative). Just like the Cartesian side, described above. And that is exactly what we need:

distance_measure = a * get<0>(c3) + b * get<1>(c3) + c * get<2>(c3) + d;

and we have the result:

int side
           = distance_measure > 0 ? 1 // left
            : distance_measure < 0 ? -1 // right
            : 0;

Spherical Side Formula

We used convenient vector calculation here, making it a mathematical exercise. Of course we could write them all out and get a lot of goniometric functions back.

One big note here: geographers (as I am) always refer to coordinates on a sphere by defining 0 at the equator. But mathematicians define 0 at the upper pole of a sphere. Watch out, it is different. This blog assumes the geographic convention, so assumes a spherical equatorial coordinate system.

Let's present it as mathematic formulae (thanks to mathcast):

The first formula converts to x,y,z. Be sure to enter all angles in radians. Also, be sure that λ is longitude (and comes first in Boost.Geometry) and δ is latitude (there are many conventions for this but I'm using them from here, and not from here and here and here). Latitude is 0 at the equator, 90 at the North Pole (again: mathematicians use 0 at the North Pole, 90 at the aquator; that is known as the spherical (polar) coordinate system, we now refer to the spherical equatorial coordinate system).

The second formula is the cross product and defines the terms of the generalized plane equation. The third formula calculates the distance measure (dm). Because (2) and (3) do not have repetitive clauses, they could be combined into (4) and then, if wished, again, with the conversion, into (5).

The sixth pseudo-formula then defines the result, it is an interpretation of the distance measure dm.

Ellipsoidal Earth

The Earth is not spherical but approximates an ellipsoid, it is flattened at the poles. But the described spherical side formula works well on Earth. Let me try to prove that (note: I'm a geographer, not a mathematician).

A location defined by (λ, δ) essentially means a vector, either on a sphere, or on an ellipsoid. These vectors always have the same direction, so the same on a sphere as on an ellipsoid. Their lengths differ (on an ellipsoid they are often shorter) but the length of a vector does not influence the plane it encompasses. Therefore the plane formed by two locations and the center of the Earth is exactly the same plane as it is using a perfect sphere.

So: two points (a segment), on a sphere or on an ellipsoid, define the same plane.

The same applies for the third point. It is also a vector, having a direction. That vector, either on a sphere or on an ellipsoid, is located at the same side of the plane.

Q.E.D.

Code

Because of the formulae above, I don't give (pseudo)code. It is quite easy. I created a small Excel sheet containing the formula is here. It contains this example (on this useful site): gcmap The formula is also in Boost.Geometry.

Summary

Using some high-school mathematics I presented an algorithm and a formula to calculate at which side a point is with resepect to a segment, on a sphere or on the Earth. Google did not find this for me. So I hope that it will be of some use for some people. At least, it is useful for Boost.Geometry

Debugger Visualizers

2011-05-18T22:40:00.004+02:00

Debugger Visualizers

When debugging in Visual Studio, watching the contents of a variable can be of great help. But as soon as a variable is of a custom type, it might be inconvenient.

For example Boost.Geometry. You can examine the contents of a point or a polygon in the Debugger Watch Window. But the representation of a point type is cumbersome: "{...}". If you click that open, you will see m_values and a meaningless pointer. Still awkward. If you click that m_values open, you see it again, in the Values column the same pointer again. Still no coordinates. If you click that pointer open, yes, you finally see the values! But it takes three steps (clicks).

If you examine a polygon, you have to do these three steps for each point in each of its rings!

You will be tired of that soon, especially if you know it can be helped using the Visual Studio Debugger Visualizers. It seems this method is not documented by Microsoft, so it is a bit of hacking, but descriptions can be found on the web (references below).

This is the recipe how to do it:

go to the folder "c:\Program Files (x86)\Microsoft Visual Studio 10.0\Common7\Packages\Debugger" (usually there, path may change, this path is for VS 2010)
make a backup of the file autoexp.dat because we are going to change it
open autoexp.dat in an editor (it is (on Windows 7) convenient to edit in another path, because saving requires Administrator rights)
go to the bottom, just before [hresult]
enter the cryptic lines shown below this bullet list
save the file (if you did edit it in another path, copy it back, giving Administrator rights)

These are the lines to add:

boost::geometry::model::d2::point_xy<*,*>{
   preview ( #("(x=", $e.m_values[0], " ,y=", $e.m_values[1], ")") )
}

It is a bit abstruse. In preview ( #( ) ) you have to enter member variables of your class (which should be called $e for enigmatic reasons). You can mix them up with strings. All is comma separated. Note that the parentheses should be placed as displayed.

After this exercise, restart Visual Studio and debug again and... you will see a much better appearance of our point in the Watch window:

We don't have to click at all. Our point is just there, x- and y- coordinates specified. Oh yes, 4.56 is now 4.55999 but that is something different, it has to do with floating point precision. We ignore that in this blog.

Definition for Boost.Geometry

At this moment I have point_xy but also model::point defined. I could have more but this is already very useful. Linestrings or polygons automatically use these point presentations, and because std::vectors are visualized by Microsoft (in that same file autoexp.dat), it is all very clear in the Watch Window.

So the section I defined looks like this:

; BOOST.GEOMETRY

boost::geometry::model::point<*,2,*>{
   preview ( #("(", $e.m_values[0], " , ", $e.m_values[1], ")") )
}

boost::geometry::model::point<*,3,*>{
   preview ( #("(", $e.m_values[0], " , ", $e.m_values[1], " , ", $e.m_values[2], ")") )
}

boost::geometry::model::d2::point_xy<*,*>{
   preview ( #("(x=", $e.m_values[0], " ,y=", $e.m_values[1], ")") )
}

I have it running in Visual Studio C++ 2010 Express, and also in Visual Studio C++ 2005 Express. Debugging visualizers is a bit tricky, if you make an error you might get an exception window. But in version 2010 it is better than it was in version 2005.

References

There is more than this. More is explained on these pages:

And in .NET (but this is actually another solution, implementing a custom visualizer in code):

http://geekswithblogs.net/technetbytes/archive/2008/06/11/122792.aspx

SqlGeometry and Boxes

2011-05-06T21:45:00.007+02:00

SqlGeometry and Boxes

I wanted to write a blog about partitioning (calling an algorithm using an on-the-fly quadtree or octree) and how it is implemented in Boost.Geometry, and how it could also be used / implemented for the SqlGeometry type. But somehow that blog will be too large and I blog here just about a box/rectangle in SqlGeometry. It is not too difficult, but I just googled and did not find other blogs saying the same things, so I hope it adds something.

STEnvelope

The OGC/ISO STEnvelope function, supported by spatial databases and many geometry libraries, returns the envelope of a geometry. An envelope is: a bounding box or bounding rectangle, also known as an axis aligned bounding box (aabb), a bbox, a minimum bounding rectangle (mbr) or an extent.

More in detail, STEnvelope returns a Geometry object describing a rectangle. OGC does not know the notion of a box. That is sometimes a bit inconvenient, but it is how it is. Boost.Geometry created a Box Concept, because Boost.Geometry is strongly typed, and envelopes returning polygons with always four or five points would confuse the non-OGC users of our library.

I blogged earlier about the spatial extent of a query in SQL Server, here.

A box is a simple geometry, with (in 2D) a minimum x and y and a maximum x and y. But how to get that from a geometry in a spatial database... I mean, the envelope returns a polygon containing five points (assuming it is closed), is the first point lower left? And is that guaranteed? The OGC specifications say: The minimum bounding box for this Geometry, returned as a Geometry. The polygon is defined by the corner points of the bounding box [(MINX, MINY), (MAXX, MINY), (MAXX, MAXY), (MINX, MAXY), (MINX, MINY)]. So yes, it should be guaranteed.

Let's look at it in more detail.

SQL Server

In SQL Server it indeed seems the case that the first point of the polygon returned by STEnvelope is the lower left point. Then, the third point must be the upper right point. Because, if you consider the polygon clockwise, it is the third point, and if it is counterclockwise, it is also the third point.

We use this query:

select geometry::STGeomFromText('POLYGON((0 0,2 5,4 1,0 0))',0).STEnvelope().STPointN(3).STAsText()

... to get the third point (which should be upper right), and it indeed gives me: POINT(4 5)

PostGIS

In PostGIS, the first point of the polygon returned by ST_Envelope is the lower left point as well. Great.

As you (maybe) know, SQLGeometry uses methods, PostGIS uses functions. So the corresponding nested function calls would be:

select ST_AsText(ST_PointN((ST_Envelope(ST_GeomFromText('POLYGON((0 0,2 5,4 1,0 0))',0))), 3))

But... no result, we get a null value back. Looking at the OGC specs, this is correct: NumPoints is defined for a LineString and not for a Polygon. SqlServer is somewhat more relaxed here (as is Boost.Geometry), returning the total number of points of a polygon (including interior rings, if any).

OK, behaviour is correct, so we try again, now inserting the ST_Boundary function (making a linestring of the boundary of a polygon). Our new function call is:

select ST_AsText(ST_PointN(ST_Boundary((ST_Envelope(ST_GeomFromText('POLYGON((0 0,2 5,4 1,0 0))',0)))), 3))

... and we have our upper right point. Using a 1 for the last value (the index of ST_PointN is one-based) would return the lower left point.

Boost.Geometry

Boost.Geometry, as said, returns a box (conforming a Box Concept) for the envelope (or assigns one). So this is more or less outside the discussion... Having a box you can get the minimum-corner coordinates (might be 2D, might be 3D, or more) and the maximum-corner coordinates.

SqlGeometry

I blogged about the useful SqlGeometry type here. The SqlGeometry type is polymorphic, it can be any OGC geometry. But, therefore, it cannot be a Rectangle (because, remember, a rectangle is not an OGC Geometry). So if it contains a Rectangle, it is a Polygon.

Using C#, it would be useful to get the coordinates from a Polygon returned by STEnvelope and that goes like this. Note that (in the light of defensive programming) we don't even assume that the first point is the lower left point. We add this code as a function (we could add it to a class GeometryExtensions, but for now we don't). It is called BoxToCoordinates, see the code below.

And, for symmetry, and because it is convenient, or often necessary (last reason most important ;-) ) we add a counterpart function as well, which creates a rectangle for you given four coordinate values. We call that one CoordinatesToBox, see code below.

The code

OK, this was some preparation for next blog, not too difficult, and now we can convert an STEnvelope result to coordinate values, and back. The full code is displayed below. Between BoxToCoordinates and CoordinatesToBox, there will be some more interesting calculations in the next blog. Stay tuned.

using System;
using Microsoft.SqlServer.Types;

namespace BlogSqlGeometryBox
{
    class Program
    {
        public static SqlGeometry FromWkt(string text, int srid)
        {
           return SqlGeometry.STGeomFromText(new System.Data.SqlTypes.SqlChars(text), srid);
        }

        public static SqlGeometry CoordinatesToBox(double x1, double y1, double x2, double y2, int srid)
        {
           SqlGeometryBuilder builder = new SqlGeometryBuilder();
           builder.SetSrid(srid);
           builder.BeginGeometry(OpenGisGeometryType.Polygon);
           builder.BeginFigure(x1, y1);
           builder.AddLine(x1, y2);
           builder.AddLine(x2, y2);
           builder.AddLine(x2, y1);
           builder.AddLine(x1, y1); // Yes, we close it neatly
           builder.EndFigure();
           builder.EndGeometry();
           return builder.ConstructedGeometry;
        }

        private static void MakeFirstSmaller<T>(ref T a, ref T b) where T : IComparable
        {
           if (a.CompareTo(b) > 0)
           {
               // Exchange the two values
               T t = a;
               a = b;
               b = t;
           }
        }

        public static void BoxToCoordinates(SqlGeometry box, out double x1, out double y1, out double x2, out double y2)
        {
           SqlGeometry lower_left = box.STPointN(1);
           SqlGeometry upper_right = box.STPointN(3);
           x1 = lower_left.STX.Value;
           y1 = lower_left.STY.Value;
           x2 = upper_right.STX.Value;
           y2 = upper_right.STY.Value;

           // Defensive programming:
           MakeFirstSmaller(ref x1, ref x2);
           MakeFirstSmaller(ref y1, ref y2);
        }

        static void Main(string[] args)
        {
           SqlGeometry triangle = FromWkt("POLYGON((0 0,2 5,4 1,0 0))", 0);
           SqlGeometry box = triangle.STEnvelope();

           double x1, y1, x2, y2;
           BoxToCoordinates(box, out x1, out y1, out x2, out y2);

           // Do something with these coordinates. Let's create a larger box
           x1 -= 1;
           y1 -= 1;
           x2 += 1;
           y2 += 1;

           box = CoordinatesToBox(x1, y1, x2, y2, 0);

           System.Console.WriteLine(box.STAsText().ToSqlString().ToString());
        }
    }
}

Writing:

POLYGON ((-1 -1, -1 6, 5 6, 5 -1, -1 -1))

SqlGeometry types and LINQ

2011-04-29T18:53:00.007+02:00

SqlGeometry types and LINQ

Last year I started to use the Microsoft .NET SqlGeometry types and I like them very much. Thanks to Bert Temme who attended us (at Geodan) on their existance. See also this page with an introduction.

In short, the SqlGeometry is the type used for the SQL Server Spatial geometry objects. But it can also be used outside of SQL Server environments: it is contained in a pure .NET assembly. The download page says: "This component can be installed separately from the server to allow client applications to use these types outside of the server."

So the SqlGeometry type is comparable to other Geometry libraries, such as Boost.Geometry, GEOS and NTS (the JTS for .NET).

SqlGeometry functionality follows the OGC/ISO specifications quite closely. This means there is not much more than OGC functionality. Sometimes there is the wish to do other things than that. For example: select only the polygons from a geometry collection. Or: select only the polygons larger than X square meters from a MultiPolygon. There comes LINQ into place.

Installing the assemblies to use SqlGeometry

To use SQLGeometry, you need to have a reference to Microsoft.SqlServer.Types, which can be found in C:\Program Files\Microsoft SQL Server\100\SDK\Assemblies (on my (64 bit) machine in c:\Program Files (x86)\..., indicating 32 bits). If not there, download it from here (there is a 64 bits version there).
Note that you have to add the reference manually, by browsing on the file system, it does not appear in the reference tab (at least not on my system, using Microsoft Visual C# 2010 Express).

After that, just add:

using Microsoft.SqlServer.Types;

The geometry collection

A geometry collection is one geometry consisting of several sub-geometries, which possibly can have different types (e.g. linestrings and polygons).

A geometry collection is sometimes not convenient. For example: MapServer cannot display them. And there are probably more GIS packages which do not display geometry collections, because most mapping software is based on layers, and most layers are defined either on a layer with points (or multi-points), or a layer with lines (or multi-linestrings), or on a layer with polygons (or multi-polygons), but in most cases not on a layer with polygons and lines, and in even more cases not on a layer with geometry collections.

Geometry collections can come into existence during an intersection. If you overlay, for example, two valid polygons which share a part of their boundaries in opposite directions, that shared boundary will form a linestring, while the rest of the intersection (if any) will form a polygon. The polygon and the linestring form together a geometry collection. The next SQL statement generates a (probably unwanted) geometry collection:

select geometry::STGeomFromText('POLYGON((0 0,0 8,8 8,8 4,4 4,4 0,0 0))', 0) 
    .STIntersection(geometry::STGeomFromText('POLYGON((4 0,4 8,8 8,8 0,4 0))', 0))

The result is: GEOMETRYCOLLECTION (POLYGON ((4 4, 8 4, 8 8, 4 8, 4 4)), LINESTRING (4 4, 4 0))

Because geometry collections are often not displayed, functionality to remove all linestrings from a geometry collection would be very useful, or crucial. SqlGeometry does not offer this functionality, but it is possible to build it. It can be built by creating WKT's of all sub-geometries, and parsing and recombining them, but LINQ offers a more elegant way.

Selecting polygons from a geometry collection

We design our algorithm as:

create a list of geometries from the geometry collection
remove the linestrings and points from this list
build a Polygon or a MultiPolygon from the polygons kept in the list

Creating the list

Creating a list of geometries can be done using the STGeometryN and STNumGeometries functions, which are available. We can iterate through the single geometries and add them to the list, using:

var list = new List<SqlGeometry>(); 
for (int i = 1; i <= geometry.STNumGeometries().Value; i++)
{ 
    list.Add(geometry.STGeometryN(i)); 
}

Note that STGeometryN uses a one-based index, instead of the usual zero-based index.

Because we think this functionality is general, and will often be reused, we make an extension method of it (I really like them!). So we make a static class called GeometryExtensions and we add a static method containing the piece above, and we can just call:

IList<SqlGeometry> list = geom.ToList();

on any variable geom of type SqlGeometry. See source code below for the extension method.

Creating a geometry from a list

Of course, we want to have the reverse as well, making an SqlGeometry from a list of polygons. The key is that an STUnion is required to make a valid result. If two polygons of the list mutually overlap, this overlap will be resolved by STUnion. Even if there is no overlap, the STUnion has to be called because there is no STAdd function available. Alternatively we could make use of the SqlGeometryBuilder type which is available in the SqlGeometry library. However, using STUnion guarantees a valid result. The next loop does this for us:

SqlGeometry result = null; 
foreach (SqlGeometry g in list)
{ 
   result = result == null ? g : result.STUnion(g);
}

Using LINQ

So having programmed step 1 and step 3 as extension methods, and having a convenient list, we can now implement step 2 using LINQ:

var result = (from g in geom.ToList() where g.STGeometryType().ToString().ToLower() == "polygon" select g).ToSqlGeometry();

The magic. This one line statement with LINQ selects the polygons from an on-the-fly list, and put them via another list back into an SqlGeometry. In the listing below it is formatted as three lines.

The complete listing

using System; 
using System.Collections.Generic; 
using System.Linq; 

using Microsoft.SqlServer.Types; 

namespace SqlGeometryBlog1 
{ 
    public static class GeometryExtensions
    { 
        // Converts a multi-geometry to a list of single-geometries 
        public static List<SqlGeometry> ToList(this SqlGeometry geometry)
        { 
          var list = new List<SqlGeometry>(); 
          for (int i = 1; i <= geometry.STNumGeometries().Value; i++)
          { 
              list.Add(geometry.STGeometryN(i)); 
          } 
          return list; 
        } 

        // Converts a list of geometries to a multi-geometry 
        public static SqlGeometry ToSqlGeometry(this IEnumerable<SqlGeometry> list)
        { 
          SqlGeometry result = null; 
          foreach (SqlGeometry g in list)
          { 
              result = result == null ? g : result.STUnion(g);
          } 
          return result; 
        } 

        // For convenience, we add two extra methods 
        public static SqlGeometry FromWkt(string text, int srid)
        { 
          return SqlGeometry.STGeomFromText(new System.Data.SqlTypes.SqlChars(text), srid);
        } 
        public static String ToWkt(this SqlGeometry geometry)
        { 
          return geometry.STAsText().ToSqlString().Value; 
        } 
    } 

    class Program 
    { 
        static void Main(string[] args) 
        { 
          var geom1 = GeometryExtensions.FromWkt("POLYGON((0 0,0 8,8 8,8 4,4 4,4 0,0 0))", 0); 
          var geom2 = GeometryExtensions.FromWkt("POLYGON((4 0,4 8,8 8,8 0,4 0))", 0); 

          var geom = geom1.STIntersection(geom2); 
          System.Console.WriteLine(geom.STIsValid().Value); 
          System.Console.WriteLine(geom.ToWkt()); 

          var list = geom.ToList(); 
          var selection = from g in list where g.STGeometryType().ToString().ToLower() == "polygon" select g; 
          geom = selection.ToSqlGeometry(); 
          System.Console.WriteLine(geom.ToWkt()); 
        } 
    } 
}

The output

True
GEOMETRYCOLLECTION (POLYGON ((4 4, 8 4, 8 8, 4 8, 4 4)), LINESTRING (4 4, 4 0))
POLYGON ((4 4, 8 4, 8 8, 4 8, 4 4))

Extent of a SQL Server Spatial table

2011-04-08T23:10:00.002+02:00

Extent of a SQL Server Spatial table

PostGIS does have a convenient spatial aggregation function: ST_Extent. SQL Server Spatial does not have that functionality. This little blog shows how to achieve its effect.

PostGIS

Given a table (e.g. called world1) with a geometry column (e.g. called geom), a call to the PostGIS function ST_Extent will give a neat box containing all geometries. For example:

select ST_Extent(geom) from world1;

returns a bounding box:

BOX(-180 -89.9,180 83.674733)

SQL Server Spatial

SQL Server Spatial follows OGC closely and therefore does not implement this functionality. Neither it has implemented it as a so-called extended geometry method. I've seen various creative solutions on the web, e.g. here and here and here. Most of these efforts are based on iterating through the rows, which is possible but often inconvenient (because of the required loop) and (also therefore) probably less performant. I will not say that my solution is performing better, but at least it is one SQL statement.

Be prepared on use of with, my current favorite SQL clause. I described it earlier here.

So what happens: we define a sort of on-the-fly view of all the envelopes (select geom.MakeValid().STEnvelope() as envelope from world1). (The call to MakeValid() is only necessary if there might be invalid geometries in the table). Then we ask for the corner points of that viewy thingy. An envelope is a rectangle, always containing five points, where the fifth point is closing point and therefore the same as the first point. So either point 1 and 3, or point 2 and 4 are the opposite corners of a rectangle. So we take point 1 and point 3 (as done in first creative solution I referenced above). We define a second on-the-fly view of this, using the first view as input. To get both diagonally opposite points, we union them (union all is required). Note that I'm referring to an SQL union here, not to a spatial union.

Those on-the-fly views can be (very conveniently) implemented using with. I name them then viewy but that's up to you. So envelope_viewy is the first one and corner_viewy is the second one.

The complete SQL statement is then, for a column called geom in a table called world:

with
  envelope_viewy as 
  (
    select geom.MakeValid().STEnvelope() as envelope from world1
  ),
  corner_viewy as 
  (
    select envelope.STPointN(1) as point from envelope_viewy 
    union all select envelope.STPointN(3) from envelope_viewy
  )
select min(point.STX) as min_x,min(point.STY) as min_y,max(point.STX) as max_x,max(point.STY) as max_y 
from corner_viewy

And voilà, you will get:

min_x    min_y    max_x    max_y
-180    -89.9    180    83.674733

The same! (Of course, I used the same table, world1, used from my shapefile research-still-to-be-finished).

You can also create a stored procedure of this (e.g. called ST_Extent).

PostGIS remark

Writing this blog, I noticed a strange thing of the PostGIS ST_Extent. It only works for 2D, there is explicitly written in the documentation that it does not work for 3D (use a ST_Extent3D for that). But it delivers a 3D box. And indeed it is explicitly typed as a 3D box. I wonder why.

Leaving Geodan

2011-03-31T19:13:00.000+02:00

Leaving Geodan...

As of tomorrow, April 1, I'm not employed with Geodan anymore. I've worked there for nearly 17 years and it has been a fantastic time. However, this year it is time for a change. I will continue with software development, now on my own, as a freelance programmer of Open Source, Geographical, Geometric, .NET and C++ software.

This change is also the reason that there were no new blogs last month on this blogspot. There are still some blog-answers to be given and blog-things to be continued. Will happen. I was quite busy finishing everything, handing over projects, with the release of Boost.Geometry, and with the new start-up.

Of course I will continue with my Open Source activities, of which Boost.Geometry is by far the most important. It is close to release now, in Boost 1.47 it will be included.

And I also take part in a charity/Open Source project, together with Geodan, I will blog about this soon. It is nice to stay related to that great company.

Next Monday I will start in Leiden for a four months contract.

Adios, all my great former Geodan-collegues, thanks for all!

Dissolving the Pentagram

2011-02-19T21:37:00.002+01:00

Dissolving the Pentagram

A week ago Denis asked the GGL-list if the internal segments of a self-intersecting polygon could be dropped. He attached this figure:

I answered that dissolve should do this. But it does not in this case. Dissolve can remove self-intersections, but under certain circumstances and even then it is not always perfect. So this is a reason to remove dissolve to the extensions folder - it should not yet be part of Boost.Geometry as it is released, hopefully in Boost 1.47

Pentagram

A pentagram is a five pointed star which can be drawn with a linestring of five segments. This is a Well-Known text presentation of a pentagram:

POLYGON((5 0,2.5 9,9.5 3.5,0.5 3.5,7.5 9,5 0))

As I often do, I compare the behaviour of algorithms on this WKT with SQL Server, PostGIS, and other packages (in this case Quantum GIS and SVG).

SVG

Boost.Geometry can create SVG images. This is really useful and I often use it in test programs. SVG's can then be depicted in FireFox. (In Google Chrome dropping the second SVG let it crash...).

SVG has thousands of properties, and one of them is the fill-rule. It can be set to winding and to non-zero. This delivers two different presentations:


winding	non-zero

SQL Server

In SQL Server the pentagram is drawn as such:

So it is using the winding rule.

The polygon is not valid, of course, it has five self-intersection points. SQL Server has a very good method to make polygons valid. It is obviously called MakeValid(). It is a SQL Server extension, not an OGC Method. I wrote "very good" and I mean that, we used it a lot in real-world problems and it has always worked very well.

With this specific pentagram MakeValid() creates a multi-polygon containing five triangles. So it directly uses the winding rule to create a multi-geometry. Sounds logical. If we calculate the area we get 17.7316840044235, the summed area of the five triangles. Calculating the area of a non-valid polygon is not possible in SQL Server. The statement is:


select geometry::STGeomFromText('POLYGON((5 0,2.5 9,9.5 3.5,0.5 3.5,7.5 9,5 0))',0)
.MakeValid().STArea()

PostGIS

PostGIS does not have a drawing engine. In PostGIS we can calculate the area of an invalid polygon and it gives us: 33.5. There is (AFAIK) not yet a function to make polygons valid, though it is announced. The good old trick is to use a buffer with buffer-distance 0, so we do this:

select ST_Area(ST_Buffer(ST_GeomFromText(
'POLYGON((5 0,2.5 9,9.5 3.5,0.5 3.5,7.5 9,5 0))',0),0))

and we get an area of 25.6158419936921. If we draw the zero-buffered result we get a filled polygon.

Boost.Geometry

Boost.Geometry can calculate the area of an invalid polygon and it also gives us: 33.5. The simple always-working Surveyor's formula just calculates this for us.

If we use the algorithm dissolve we indeed get a polygon with all internal corners dropped in this case. The area of the dissolved polygon is 25.6158412 and this is the right area. So the Surveyor's formula actually cannot be used with invalid polygons, SQL Server is right to prohibit this.

The SVG files displayed above are both created using Boost.Geometry so no picture here.

Quantum GIS

Quantum GIS has a great extension, QuickWKT, with which WKT geometries can be drawn. It uses the winding rule. It can calculate the area of the polygons (using the Info tool we can get it) and we get 33.5. Right, the Surveyor's formula.

What should happen

What should happen is not that simple. In these four programs, only SQL Server was consistent by prohibiting area calculation and creating a multi-geometry which looks like the image it depicts. All other programs and libraries were consistent, using the Surveyor's formula and creating a closed five-pointed polygon.

The reasoning of Boost.Geometry is that it takes direction into account. So if we draw arrows aside of all linestrings, we get:

And here we see, in the triangle with the smaller arrows, that the direction is not consistent. That triangle should not be generated. Of this whole pentagram several polygons can be drawn with consistent directions, and the largest of them is the whole polygon. So the dissolve function is expected to generate a five-pointed star here. And it does.

There is also a pentagon in the center with consistent directions. Because it has the same orientation as the outer polygon, it should be discarded. And it is. If it would have been orientated counterclockwise, a hole should have been formed.

So there are reasons to say that this pentagram indeed should be converted into a solid five-pointed star. There is no winding rule or non-zero rule involved there, it is just derived by the directions of the segments created by the intersection.

But now the sample

Anyway, using the sample of Denis Pesotsky, Boost.Geometry's results were not that clever. It generates two polygons which are indeed as expected, with consistent directions. The outer polygon does not have a consistent direction and therefore should not be generated, correct. But there is one polygon with consistent direction missing, in the lower right. That one is discarded but should not have been...

Because this, and because we have to think more about the behaviour, the dissolve algorithm has to be moved to the extensions folder. That means it will not be part of the initial release of Boost.Geometry, in the Release branch.

This is the WKT of the input polygon:

POLYGON((55 10, 141 237, 249 23, 21 171, 252 169, 24 89, 266 73, 55 10))

The white pieces are by the winding rule, and in SQL Server it is looking the same (no picture here). SQL Server's MakeValid() creates 9 polygons here.

PostGIS, using the buffer, creates one polygon, looking like this:

And that polygon is indeed the largest polygon which can be formed, following all arrows in clockwise direction. So this is the polygon that was expected by the dissolve algorithm...

Anyway, it is not the result that Denis asked for.

Useful usage

Dissolve can be very useful so I end this blog by just listing some of the testcases where it makes sense:

Inaccurate corners (extra intersection):

Messy corners:

Keyholed holes:

Conclusions

Due to shortcomings in the dissolve algorithm and vagueness in behaviour in general we have to move the dissolve function to the extensions folder. It has to be thought over more thoroughly.

Qt World Mapper with Boost.Geometry

2011-02-13T23:46:00.005+01:00

Qt World Mapper with Boost.Geometry

Last year I created a simple world mapper using WxWidgets. For various reasons I decided to make the Qt equivalent as well.

Both Qt and WxWidgets are well-known cross-platform windows frameworks. Qt is used in Quantum GIS.

The World Mapper I created with WxWidgets did read world countries from a WKT (well-known text) file, and displayed them. Moving the mouse did highlight the country. For Qt it is similar (the highlighting is not yet been done).

The nice thing is that, by setting Boost.Geometry in this context, its force is shown:

Look below to lines marked with // Adapt. The QPointF and QPolygonF are just registered. From now on they act as normal Boost.Geometry geometries. So they can be used in area calculation, distance, transform, etc. We only use transform in this example.

So also look at lines commented with // This is the essention. Here a QPolygonF is declared, a ring (a normal Boost.Geometry part of a polygon, part of a multi-polygon) is transformed using the transformer to a QPolygonF. And the last one is of course displayed without any problem.

Delivering this application:

Complete C++ Code

#include <fstream> 

#include <QtGui> 
#include <QWidget> 
#include <QObject> 
#include <QPainter> 

#include <boost/foreach.hpp> 

#include <boost/geometry/geometry.hpp> 
#include <boost/geometry/geometries/register/point.hpp>
#include <boost/geometry/geometries/register/ring.hpp>

#include <boost/geometry/multi/multi.hpp> 
#include <boost/geometry/extensions/algorithms/selected.hpp>
#include <boost/geometry/extensions/gis/io/wkt/wkt.hpp> 

// Adapt a QPointF such that it can be handled by Boost.Geometry
BOOST_GEOMETRY_REGISTER_POINT_2D_GET_SET(QPointF, double, cs::cartesian, x, y, setX, setY)

// Adapt a QPolygonF as well. 
// A QPolygonF has no holes (interiors) so it is similar to a Boost.Geometry ring
BOOST_GEOMETRY_REGISTER_RING(QPolygonF) 

typedef boost::geometry::model::d2::point_xy<double> point_2d; 
typedef boost::geometry::model::multi_polygon 
    < 
        boost::geometry::model::polygon<point_2d> 
    > country_type; 


class WorldMapper : public QWidget
{ 
public: 
    WorldMapper(std::vector<country_type> const& countries, boost::geometry::model::box<point_2d> const& box)
        : m_countries(countries) 
        , m_box(box) 
    { 
        setPalette(QPalette(QColor(200, 250, 250)));
        setAutoFillBackground(true); 
    } 

protected: 
    void paintEvent(QPaintEvent*) 
    { 
        map_transformer_type transformer(m_box, this->width(), this->height()); 

        QPainter painter(this); 
        painter.setBrush(Qt::green); 

        BOOST_FOREACH(country_type const& country, m_countries)
        { 
          typedef boost::range_value<country_type>::type polygon_type;
          BOOST_FOREACH(polygon_type const& polygon, country)
          { 
           typedef boost::geometry::ring_type<polygon_type>::type ring_type;
           ring_type const& ring = boost::geometry::exterior_ring(polygon);

           // This is the essention: 
           // Directly transform from a multi_polygon (ring-type) to a QPolygonF 
           QPolygonF qring; 
           boost::geometry::transform(ring, qring, transformer); 

           painter.drawPolygon(qring); 
          } 
       } 
    } 

private: 
    typedef boost::geometry::strategy::transform::map_transformer
        < 
          point_2d, QPointF, 
          true, true 
        > map_transformer_type; 

     std::vector<country_type> const& m_countries;
     boost::geometry::model::box<point_2d> const& m_box;
}; 


class MapperWidget : public QWidget
{ 
public: 
     MapperWidget(std::vector<country_type> const& countries, boost::geometry::model::box<point_2d> const& box, QWidget *parent = 0)
         : QWidget(parent) 
     { 
         WorldMapper* mapper = new WorldMapper(countries, box);

         QPushButton *quit = new QPushButton(tr("Quit")); 
         quit->setFont(QFont("Times", 18, QFont::Bold)); 
         connect(quit, SIGNAL(clicked()), qApp, SLOT(quit())); 

         QVBoxLayout *layout = new QVBoxLayout; 
         layout->addWidget(mapper); 
         layout->addWidget(quit); 
         setLayout(layout); 
     } 

}; 


template <typename Geometry, typename Box>
inline void read_wkt(std::string const& filename, std::vector<Geometry>& geometries, Box& box)
{ 
    std::ifstream cpp_file(filename.c_str()); 
    if (cpp_file.is_open()) 
    { 
        while (! cpp_file.eof() ) 
        { 
          std::string line; 
          std::getline(cpp_file, line); 
          if (! line.empty()) 
          { 
              Geometry geometry; 
              boost::geometry::read_wkt(line, geometry); 
              geometries.push_back(geometry); 
              boost::geometry::combine(box, boost::geometry::make_envelope<Box>(geometry)); 
          } 
        } 
    } 
} 

int main(int argc, char *argv[])
{ 
    std::vector<country_type> countries; 
    boost::geometry::model::box<point_2d> box; 
    boost::geometry::assign_inverse(box); 
    read_wkt("../data/world.wkt", countries, box); 

    QApplication app(argc, argv); 
    MapperWidget widget(countries, box); 
    widget.setWindowTitle("Boost.Geometry for Qt - Hello World!"); 
    widget.setGeometry(50, 50, 800, 500);
    widget.show(); 
    return app.exec(); 
}

Conclusions

With less than 160 lines of code, Qt and Boost.Geometry we showed how it can quite conveniently combine together.
The source is shown above and in SVN it can be found here.

Actually this application (and its WxWidgets counterpart) was for me the reason to start the shapefile research, which final conclusions still have to be drawn.

Spikes: research with ESRI

2011-02-13T13:37:00.003+01:00

Spikes: research with ESRI

A short follow-up on my previous blog about spikes and their causes. My collegue Bert recoded the scenario in C# using ESRI. Thanks Bert! The full code is shown below.

The answer is this: with ESRI the spikes are not there. The program produces:

area: 295.44447138901, perimeter: 92.5058167188044

select geometry::STGeomFromText ('POLYGON((180956.436999999
313716.998,180953.978700001 313701.0152,180926.420299999
313732.510699999,180954.359000001
313720.504000001,180956.436999999 313716.998,180956.436999999
313716.998))', 28992)

And this result is right, area and perimeter are without spikes and if you visualize this WKT (e.g. with SQL Server) you get the right picture...

Interesting is that all coordinates seems to be rounded on one millimeter or a tenth of it. That is probably the result of the fuzzy tolerance policy that ESRI uses. Personally I'm normally a bit sceptic about that approach: though it can avoid spikes and slivers, it also brings some fuzzyness into the result. It is sometimes quite hard to determine a tolerance which is sufficient for cleaning and sufficient for not destroying features... Anyway, in this scenario it works well.

It seems interesting to check, if you use rounding in e.g. Boost.Geometry, would the spikes be gone... The answer is (of course), sometimes yes, sometimes no. The rounding is applied on a larger grid, so the effect remains the same. Rounding is not enough. To get rid of spikes, in this scenario, the spike removal functionality is advised. Which uses a gap distance, which can be seen as a sort of fuzzy tolerance...

C# Code for "difference after cut" with ESRI

using System; 
using System.Collections.Generic; 
using System.ComponentModel; 
using System.Data; 
using System.Drawing; 
using System.Text; 
using System.Windows.Forms; 
using ESRI.ArcGIS.Geometry; 
using ESRI.ArcGIS.esriSystem; 
using Microsoft.SqlServer.Types; 
using System.Data.SqlTypes; 

namespace SpikesResearchWithEsri 
{ 
    public partial class Form1 : Form
    { 
        public Form1() 
        { 
          InitializeComponent(); 
        } 

        private void InitializeLicense()
        { 
          AoInitialize aoi = new AoInitializeClass(); 

          //Additional license choices can be included here. 
          esriLicenseProductCode productCode = 
          esriLicenseProductCode.esriLicenseProductCodeArcEditor; 
          if (aoi.IsProductCodeAvailable(productCode) == 
          esriLicenseStatus.esriLicenseAvailable) 
          { 
          aoi.Initialize(productCode); 
          } 
        } 


        private void button1_Click(object sender, EventArgs e)
        { 


          ESRI.ArcGIS.RuntimeManager.Bind(ESRI.ArcGIS.ProductCode.Desktop); 
          InitializeLicense(); 
          IPolygon parcel = this.GetParcel(); 
          IPolyline line = this.GetCutLine(); 

          IPolygon bufferedLine=(IPolygon)((ITopologicalOperator)line).Buffer(0.1); 
          
          IGeometry result=((ITopologicalOperator)parcel).Difference(bufferedLine); 
          IPolygon4 p = (IPolygon4)result; 

          // Big part: 
          IPolygon firstPolygon = this.GetPolygon(p, true);
          // Small part: 
          IPolygon secondPolygon = this.GetPolygon(p, false);


          // use if we do a cutpolygon 
          //IGeometry leftPolygon; 
          //IGeometry rightPolygon; 
          //((ITopologicalOperator)parcel).Cut(line, out leftPolygon, out rightPolygon); 
          //IPolyline cutline=(IPolyline)((ITopologicalOperator)(ICurve)parcel).Intersect(line, esriGeometryDimension.esriGeometryNoDimension); 
          //IPolygon firstPolygon = (IPolygon)leftPolygon; 
          //IPolygon secondPolygon = (IPolygon)rightPolygon; 
          //this.textBox1.Text=("from: " + cutline.FromPoint.X.ToString() + ", " + 
          // cutline.FromPoint.Y.ToString() + ", " + 
          // "to: " + cutline.ToPoint.X.ToString() + ", " + 
          // cutline.ToPoint.Y.ToString()) + Environment.NewLine; 

          this.textBox1.Text += "Small polygon stats: " + System.Environment.NewLine; 
          this.textBox1.Text += GetStats(secondPolygon); 
          this.textBox1.Text += GetWkt(secondPolygon); 

          // now calculate the difference between original and the big part (result should be the small part) 
          IPolygon4 pol2= (IPolygon4)((ITopologicalOperator)parcel).Difference(firstPolygon); 
          this.textBox1.Text += "Result Small polygon stats: " + System.Environment.NewLine; 
          this.textBox1.Text += GetStats(pol2); 
          this.textBox1.Text += GetWkt(pol2); 
        } 

        private string GetStats(IPolygon pol)
        { 
          double area = ((IArea)pol).Area; 
          double perimeter = ((ICurve)pol).Length; 
          return "area: " + area.ToString() + ", perimeter: " + perimeter.ToString() + System.Environment.NewLine; 
        } 


        private string GetWkt(IPolygon pol)
        { 
          string res = GeomHelper.aoGeomToWkt(pol, false, 28992);
          res = res.Substring(0, res.Length - 1);
          string print = "select geometry::STGeomFromText ('" + res + "', 28992)" + Environment.NewLine; 
          return print; 
        } 

        public IPolygon GetPolygon(IPolygon4 input, bool first)
        { 
          IGeometryBag exteriorRings = input.ExteriorRingBag; 
          IEnumGeometry exteriorRingsEnum = exteriorRings as IEnumGeometry;
          exteriorRingsEnum.Reset(); 

          IRing firstRing = exteriorRingsEnum.Next() as IRing;
          IPointCollection pcFirst = (IPointCollection)firstRing; 

          IRing secondRing = exteriorRingsEnum.Next() as IRing;
          IPointCollection pcSecond = (IPointCollection)secondRing; 

          IPolygon p = new PolygonClass(); 

          if(first)((IPointCollection)p).AddPointCollection(pcFirst); 
          if(!first)((IPointCollection)p).AddPointCollection(pcSecond); 
          return p; 
        } 


        private IPolygon GetParcel() 
        { 
          string wkt = "POLYGON ((180956.437 313716.998, 180954.359 313720.504, 180906.303 313741.156, 180896.82 313744.272, 180890.456 313748.17, 180831.49 313797.135, 180809.668 313816.098, 180776.418 313842.336, 180771.483 313846.102, 180746.934 313843.764, 180731.738 313840.386, 180717.32 313833.763, 180710.827 313828.438, 180704.982 313820.515, 180699.267 313811.813, 180698.488 313801.162, 180700.826 313793.889, 180705.539 313791.06, 180706.021 313790.771, 180710.156 313789.208, 180766.157 313768.042, 180825.515 313744.013, 180878.636 313724.271, 180945.593 313695.884, 180952.801 313693.358, 180956.437 313716.998))"; 
          IGeometryFactory factory = new GeometryEnvironmentClass();
          SqlGeometry g = SqlGeometry.Parse(wkt); 
          IGeometry geom2 = new PolygonClass(); 
          int countout; 
          factory.CreateGeometryFromWkbVariant(g.STAsBinary().Value, out geom2, out countout);
          geom2.SpatialReference = this.GetProjectedSpatialReference(28992);
          return (IPolygon) geom2; 
        } 

        private IPolyline GetCutLine() 
        { 
          string line = "LINESTRING(180955 313700,180920 313740)"; 
          IGeometryFactory factory = new GeometryEnvironmentClass();
          SqlGeometry g = SqlGeometry.Parse(line); 
          IGeometry geom2 = new PolylineClass(); 
          geom2.SpatialReference = this.GetProjectedSpatialReference(28992);
          int countout; 
          
          factory.CreateGeometryFromWkbVariant(g.STAsBinary().Value, out geom2, out countout);
          return (IPolyline)geom2; 

        } 


        private ISpatialReference GetProjectedSpatialReference(int pcsType)
        { 
          ISpatialReferenceFactory pSRF = new SpatialReferenceEnvironmentClass();
          IProjectedCoordinateSystem m_ProjectedCoordinateSystem = pSRF.CreateProjectedCoordinateSystem(pcsType);
          ISpatialReference spatialReference = (ISpatialReference)m_ProjectedCoordinateSystem; 
          return spatialReference; 
        } 

    } 
}

Precision, the cause of spikes

2011-01-30T21:30:00.002+01:00

Precision, the cause of spikes

This is a follow-up on my last blog. The spikey effects we have seen are perfectly explainable and understandable.

What we did was:

cut a parcel in two pieces

create a cutline and make a (nearly linear) cutting polygon of it (using buffer)
subtract that cutline from the geometry, using STDifference

keep one piece, using STGeometryN
subtract that piece from the parcel, using STDifference

The result should be the discarded piece, but we got back two free spikes extra.

This short blog will explain this a little more and, as a side note, use just another SQL syntax for the same result.

Side note: with

The query from my previous blog can be rewritten much more elegantly using the with syntax:

with myquery as 
(
  select
    geometry::STGeomFromWKB(0x0103000000010000001A00000023DBF97EE316064146B6F3FDD32513415A643BDFD216064175931804E225134196438B6C52150641C976BE9F34261341F6285C8F06150641022B871641261341F853E3A5D3140641E17A14AE50261341B81E85EBFB120641A4703D8A142713414E6210584D120641AC1C5A64602713414E621058431106414E621058C9271341A01A2FDD1B11064121B07268D8271341F4FDD478571006411904560ECF271341448B6CE7DD0F06418195438BC1271341F6285C8F6A0F06413BDF4F0DA72713410E2DB29D360F06416F1283C091271341E5D022DB070F0641F6285C0F7227134160E5D022DA0E06416F1283404F271341448B6CE7D30E0641F853E3A52427134154E3A59BE60E06411904568E07271341643BDF4F0C0F0641D7A3703DFC2613414A0C022B100F064125068115FB26134191ED7C3F310F0641B6F3FDD4F42613414C378941F11006414A0C022BA0261341EC51B81ECC1206413BDF4F0D40261341022B87167514064125068115F1251341B4C876BE8C160641C74B37897F25134121B07268C6160641508D976E7525134123DBF97EE316064146B6F3FDD3251341,28992)
      as parcel,
    geometry::STGeomFromText('LINESTRING(180955 313700,180920 313740)', 28992).STBuffer(0.1)
      as cutline
)
select parcel.STDifference((select parcel.STDifference(cutline).STGeometryN(2) from myquery))
from myquery

The with clause is useful to avoid repeating pieces (as having the WKB twice), and is essential to do recursive queries. You define an inlined viewy thingy (here "myquery") and use it.

Both SQL Server and PostGreSQL support with queries. See also this link which explains with in depth.

Floating point precision

Floating points are really awful. Everyone knows that if you divide one by six (1/6) and multiply the result with six ((1/6) * 6) you should get one back. However, with floating points this is not always the case. You might get 0.9999999 back. If you try it, you might also get the correct 1.0 back. That is (according to this book) because intermediate results are stored in the floating point register of the CPU. So it is by chance that it is right. In reality it will be wrong.

See also this link and scroll to "Rounding errors can play havoc with math-intense programs". Adding ten 0.1 values will result into 0.99999999999999989. I repeated this and indeed this is the case with a double. With a float the result is 1.0000001192092896. With ttmath it is 1.0, but as soon as you do the same trick with eleven 1/11 additions, you might get 0.9999999999999999999999999999999 back, or 1.0, depending on which precision you specify. Using Boost.Rational you get always 1/1, guaranteed, because it calculates in another way. It uses fractions and keeps fractions all the time. That is great in this use case.

What then caused the spikes

Now that we realize again that floating points are awful, we investigate how imprecision causes the spikes. The damage happens in the first phase. The first STDifference is used to cut the parcel into two pieces.

See the picture below. The grid is the so-called floating point grid. The smallest difference in a 32 bits floating point value (1.175e-38) is depicted. Any point (with float) must be located on one of the gridpoints. (In reality the FP-grid is irregular, but the idea will be clear).

So in this case (in most cases), the intersection point happens to fall between the gridpoints... All the vertices of both input polygons are located on the gridpoints, of course. But the real intersection point does not. It is rounded to one of the four gridpoints. It is just by chance to which one it is rounded. So the rounded intersection point is located on a gridpoint.

Suppose it rounds to the point labeled "3". We then get a resulting polygon as partly depicted here:

The resulting polygon (in blue) is too small. There is a small piece of the original polygon which is not covered by this piece, at the north side, upperleft of point 3. Subtract this resulting polygon from the original and of course, you will get the spike. There is also a similar small piece at the east side but that will not create a visual effect in this use case.

There is no obvious solution. It is the same as adding ten times 0.1. You get too much or too little. These spikes are too little in the first phase, and after the second difference, the two resulting spikes are a bit too much. This explains that sometimes you get zero, sometimes you get one and sometimes you get two spikes. It is by chance. If the rounding would go to point "1" in the picture above, there would be no spike there.

The spikes are not the result of flawed logic, or wrong decisions somewhere. It is just the rounding.

Use ttmath if you want to minimize the chance on it. Or work around it, e.g. by buffering the result with a small value e.g. -0.000001. Boost.Geometry now has some algorithms to find spikes and remove them.

By the way, the nice trick on cutting a parcel into two pieces was introduced in our project by our partners from B3Partners. The project where this study is from extracted is the Portaal Natuur en Landschap

Conclusions

In this use case, spikes are the geometric visualization of floating point imprecision.

SQL Server, PostGIS, STDifference and Precision

2011-01-28T21:45:00.011+01:00

SQL Server, PostGIS, STDifference and Precision

Today I encountered a problem related to the OGC difference function. I will describe this problem in this blog. I created a sort of minimal sample, without database, and present that here.

The problem

The problem I encountered is precision: both SQL Server (where I encountered it) and PostGIS (where I verified it) do not handle the situation described here sufficiently. The results is a polygon with two spikes which should not be there. The cause must be imprecision of floating point types. The message is: be warned!

The situation is as following: there is a parcel. That parcel is cut into two pieces (I will explain below how). One piece is thrown away, the other piece is kept. That kept piece is subtracted from the original parcel. What you would expect is having the original piece that was earlier thrown away. But what is happening is that you get one or two spikes extra...

The parcel and the cutline are shown here:

Cutting into two pieces

We cut (for this experiment) the polygon into two pieces by buffering the cutline, and subtracting (with STDifference) that from the parcel. The cutline is shown in the figure above as the dotted red line. Then we keep one of the two resulting polygons using the STGeometryN function. By trying we get the left piece which we want to keep by index 2.

So this piece of SQL is, for SQL Server:

select geometry::STGeomFromWKB(0x0103000000010000001A00000023DBF97EE316064146B6F3FDD32513415A643BDFD216064175931804E225134196438B6C52150641C976BE9F34261341F6285C8F06150641022B871641261341F853E3A5D3140641E17A14AE50261341B81E85EBFB120641A4703D8A142713414E6210584D120641AC1C5A64602713414E621058431106414E621058C9271341A01A2FDD1B11064121B07268D8271341F4FDD478571006411904560ECF271341448B6CE7DD0F06418195438BC1271341F6285C8F6A0F06413BDF4F0DA72713410E2DB29D360F06416F1283C091271341E5D022DB070F0641F6285C0F7227134160E5D022DA0E06416F1283404F271341448B6CE7D30E0641F853E3A52427134154E3A59BE60E06411904568E07271341643BDF4F0C0F0641D7A3703DFC2613414A0C022B100F064125068115FB26134191ED7C3F310F0641B6F3FDD4F42613414C378941F11006414A0C022BA0261341EC51B81ECC1206413BDF4F0D40261341022B87167514064125068115F1251341B4C876BE8C160641C74B37897F25134121B07268C6160641508D976E7525134123DBF97EE316064146B6F3FDD3251341,28992)
.STDifference(geometry::STGeomFromText('LINESTRING(180955 313700,180920 313740)', 28992).STBuffer(0.1)).STGeometryN(2)

And for PostGIS it is:

select ST_GeometryN(ST_Difference(
ST_SetSRID('0103000000010000001A00000023DBF97EE316064146B6F3FDD32513415A643BDFD216064175931804E225134196438B6C52150641C976BE9F34261341F6285C8F06150641022B871641261341F853E3A5D3140641E17A14AE50261341B81E85EBFB120641A4703D8A142713414E6210584D120641AC1C5A64602713414E621058431106414E621058C9271341A01A2FDD1B11064121B07268D8271341F4FDD478571006411904560ECF271341448B6CE7DD0F06418195438BC1271341F6285C8F6A0F06413BDF4F0DA72713410E2DB29D360F06416F1283C091271341E5D022DB070F0641F6285C0F7227134160E5D022DA0E06416F1283404F271341448B6CE7D30E0641F853E3A52427134154E3A59BE60E06411904568E07271341643BDF4F0C0F0641D7A3703DFC2613414A0C022B100F064125068115FB26134191ED7C3F310F0641B6F3FDD4F42613414C378941F11006414A0C022BA0261341EC51B81ECC1206413BDF4F0D40261341022B87167514064125068115F1251341B4C876BE8C160641C74B37897F25134121B07268C6160641508D976E7525134123DBF97EE316064146B6F3FDD3251341'
 ::geometry, 28992)
, (ST_Buffer(ST_GeomFromText('LINESTRING(180955 313700,180920 313740)',28992),0.1))), 2)

Note that the 28992 is just an SRID of a coordinate system (the Dutch), it probably works (or better stated: does not work) with any coordinate system. The result is this:

Subtracting one of the pieces from the original

After cutting we subtract the piece we kept from the original parcel. This will result in the other piece. In fact we get the area covered by the small buffer as well, but that is hardly visible, and is not important for this story.

This is how it should look like:

And this is the complete query in SQL Server:

select
geometry::STGeomFromWKB(0x0103000000010000001A00000023DBF97EE316064146B6F3FDD32513415A643BDFD216064175931804E225134196438B6C52150641C976BE9F34261341F6285C8F06150641022B871641261341F853E3A5D3140641E17A14AE50261341B81E85EBFB120641A4703D8A142713414E6210584D120641AC1C5A64602713414E621058431106414E621058C9271341A01A2FDD1B11064121B07268D8271341F4FDD478571006411904560ECF271341448B6CE7DD0F06418195438BC1271341F6285C8F6A0F06413BDF4F0DA72713410E2DB29D360F06416F1283C091271341E5D022DB070F0641F6285C0F7227134160E5D022DA0E06416F1283404F271341448B6CE7D30E0641F853E3A52427134154E3A59BE60E06411904568E07271341643BDF4F0C0F0641D7A3703DFC2613414A0C022B100F064125068115FB26134191ED7C3F310F0641B6F3FDD4F42613414C378941F11006414A0C022BA0261341EC51B81ECC1206413BDF4F0D40261341022B87167514064125068115F1251341B4C876BE8C160641C74B37897F25134121B07268C6160641508D976E7525134123DBF97EE316064146B6F3FDD3251341,28992)
.STDifference
(
  (
 select 
 geometry::STGeomFromWKB(0x0103000000010000001A00000023DBF97EE316064146B6F3FDD32513415A643BDFD216064175931804E225134196438B6C52150641C976BE9F34261341F6285C8F06150641022B871641261341F853E3A5D3140641E17A14AE50261341B81E85EBFB120641A4703D8A142713414E6210584D120641AC1C5A64602713414E621058431106414E621058C9271341A01A2FDD1B11064121B07268D8271341F4FDD478571006411904560ECF271341448B6CE7DD0F06418195438BC1271341F6285C8F6A0F06413BDF4F0DA72713410E2DB29D360F06416F1283C091271341E5D022DB070F0641F6285C0F7227134160E5D022DA0E06416F1283404F271341448B6CE7D30E0641F853E3A52427134154E3A59BE60E06411904568E07271341643BDF4F0C0F0641D7A3703DFC2613414A0C022B100F064125068115FB26134191ED7C3F310F0641B6F3FDD4F42613414C378941F11006414A0C022BA0261341EC51B81ECC1206413BDF4F0D40261341022B87167514064125068115F1251341B4C876BE8C160641C74B37897F25134121B07268C6160641508D976E7525134123DBF97EE316064146B6F3FDD3251341,28992)
 .STDifference(geometry::STGeomFromText('LINESTRING(180955 313700,180920 313740)', 28992).STBuffer(0.1)).STGeometryN(2)
 )
)

And in PostGIS:


select ST_Perimeter(ST_Difference(
ST_SetSRID('0103000000010000001A00000023DBF97EE316064146B6F3FDD32513415A643BDFD216064175931804E225134196438B6C52150641C976BE9F34261341F6285C8F06150641022B871641261341F853E3A5D3140641E17A14AE50261341B81E85EBFB120641A4703D8A142713414E6210584D120641AC1C5A64602713414E621058431106414E621058C9271341A01A2FDD1B11064121B07268D8271341F4FDD478571006411904560ECF271341448B6CE7DD0F06418195438BC1271341F6285C8F6A0F06413BDF4F0DA72713410E2DB29D360F06416F1283C091271341E5D022DB070F0641F6285C0F7227134160E5D022DA0E06416F1283404F271341448B6CE7D30E0641F853E3A52427134154E3A59BE60E06411904568E07271341643BDF4F0C0F0641D7A3703DFC2613414A0C022B100F064125068115FB26134191ED7C3F310F0641B6F3FDD4F42613414C378941F11006414A0C022BA0261341EC51B81ECC1206413BDF4F0D40261341022B87167514064125068115F1251341B4C876BE8C160641C74B37897F25134121B07268C6160641508D976E7525134123DBF97EE316064146B6F3FDD3251341'::geometry,28992)
,
(select ST_GeometryN(ST_Difference(
ST_SetSRID('0103000000010000001A00000023DBF97EE316064146B6F3FDD32513415A643BDFD216064175931804E225134196438B6C52150641C976BE9F34261341F6285C8F06150641022B871641261341F853E3A5D3140641E17A14AE50261341B81E85EBFB120641A4703D8A142713414E6210584D120641AC1C5A64602713414E621058431106414E621058C9271341A01A2FDD1B11064121B07268D8271341F4FDD478571006411904560ECF271341448B6CE7DD0F06418195438BC1271341F6285C8F6A0F06413BDF4F0DA72713410E2DB29D360F06416F1283C091271341E5D022DB070F0641F6285C0F7227134160E5D022DA0E06416F1283404F271341448B6CE7D30E0641F853E3A52427134154E3A59BE60E06411904568E07271341643BDF4F0C0F0641D7A3703DFC2613414A0C022B100F064125068115FB26134191ED7C3F310F0641B6F3FDD4F42613414C378941F11006414A0C022BA0261341EC51B81ECC1206413BDF4F0D40261341022B87167514064125068115F1251341B4C876BE8C160641C74B37897F25134121B07268C6160641508D976E7525134123DBF97EE316064146B6F3FDD3251341'::geometry,28992)
,(ST_Buffer(ST_GeomFromText('LINESTRING(180955 313700,180920 313740)', 28992),0.1))),2)
)))

Results

To my surprise I got two spikes in SQL Server. To verify this, I looked into PostGIS. I got two spikes as well. First I thought it was caused by precision, going to WKT and back. But also if you use WKB, or if you use the geometry from the underlying database (not shown here), you will get the same. So it is the imprecision of the handling. When varying the cutline a bit (changing 180955 to 180960), I get one spike in SQL Server and still two in PostGIS.

The spikes are shown here, bordering the blue result:

The perimeter of the result is interesting. The result, without spikes, should be 92.5 meter. But including the spikes it is 151.792871568541 (in SQL Server, using the function .STLength() and 151.792880859545 (in PostGIS, using the function ST_Perimeter(...).

The area is everywhere about 295.4 square meter, with or without spikes, indicating that the spikes have an area about zero (so those are real spikes).

Why would you do such a thing... Creating a parcel, cutting, throwing away one, subtracting and getting the other one...

Well, I was just testing our application. Our application has WebGIS functionality to create parcels. And to cut parcels into two pieces. And to throw away pieces. And to fill-up space which is left from a parcel and a piece. So just this scenario. And to my annoyance, often you get spikes! This project is done with SQL Server, and we get spikes. But this little research shows that, would we have done it with PostGIS, we would have had resulting spikes as well.

It is as a simple mathematical exercise. We cut 5 into two pieces, 3 and 2. We throw away the 2 and keep the 3. We subtract 3 from our 5. And we get the 2 back, of course. In geometry we subtract (difference) twice but that is just because the cut is not an OGC operation. So we trick it.

Note that I didn't get this just once. I got it several times for different geometries. Note also that it can easily be solved using a STBuffer of -0.00001 on the resulting polygon, which will eliminate the spikes and have no visual impact (it has an impact though...).

With Boost.Geometry

Of course I was curious what Boost.Geometry would do with this.

The buffer of Boost.Geometry is not completely implemented, so I create the buffer in SQL Server, saving the WKB.

select geometry::STGeomFromText('LINESTRING(180955 313700,180920 313740)', 28992).STBuffer(0.1).STAsBinary()

and draft some C++ code:

template <typename T> 
void test_difference_parcel_precision() 
{
    namespace bg = boost::geometry; 
    typedef bg::model::d2::point_xy<T> point_type; 
    typedef bg::model::polygon<point_type> polygon_type;
    typedef std::vector<boost::uint8_t> byte_vector;

    polygon_type parcel, buffer; 

    { 
        byte_vector wkb; 
        bg::hex2wkb("0103000000010000001A00000023DBF97EE316064146B6F3FDD32513415A643BDFD216064175931804E225134196438B6C52150641C976BE9F34261341F6285C8F06150641022B871641261341F853E3A5D3140641E17A14AE50261341B81E85EBFB120641A4703D8A142713414E6210584D120641AC1C5A64602713414E621058431106414E621058C9271341A01A2FDD1B11064121B07268D8271341F4FDD478571006411904560ECF271341448B6CE7DD0F06418195438BC1271341F6285C8F6A0F06413BDF4F0DA72713410E2DB29D360F06416F1283C091271341E5D022DB070F0641F6285C0F7227134160E5D022DA0E06416F1283404F271341448B6CE7D30E0641F853E3A52427134154E3A59BE60E06411904568E07271341643BDF4F0C0F0641D7A3703DFC2613414A0C022B100F064125068115FB26134191ED7C3F310F0641B6F3FDD4F42613414C378941F11006414A0C022BA0261341EC51B81ECC1206413BDF4F0D40261341022B87167514064125068115F1251341B4C876BE8C160641C74B37897F25134121B07268C6160641508D976E7525134123DBF97EE316064146B6F3FDD3251341", std::back_inserter(wkb)); 
        bg::read_wkb(wkb.begin(), wkb.end(), parcel); 
    } 
    { 
        byte_vector wkb; 
        bg::hex2wkb(
"01030000000100000083000000000032FCD716064100009E998F25134100000706D81606410040A4998F2513410000DA0FD816064100C0E6998F2513410000A819D81606410080659A8F25134100806A23D816064100C0209B8F25134100801E2DD81606410080189C8F2513410000BE36D816064100404D9D8F25134100004540D816064100C0BE9E8F2513410000AF49D816064100806DA08F2513410000F752D8160641008059A28F2513410000195CD816064100C082A48F25134100800F65D81606410080E9A68F2513410000D66DD816064100408EA98F25134100006876D816064100C070AC8F2513410000C17ED8160641000091AF8F2513410080DC86D816064100C0EFB28F25134100009E8ED816064100C081B68F2513410080EC95D816064100803ABA8F2513410080C79CD8160641000018BE8F25134100002FA3D8160641008017C28F251341008022A9D8160641000037C68F2513410080A1AED8160641000074CA8F2513410000ACB3D81606410040CCCE8F251341000042B8D816064100403DD38F251341000062BCD81606410000C5D78F25134100000DC0D8160641000061DC8F251341000042C3D816064100000FE18F251341000001C6D81606410080CCE58F251341008049C8D8160641004097EA8F25134100001BCAD816064100006DEF8F251341008075CBD816064100804BF48F251341008058CCD8160641004030F98F2513410000C4CCD8160641000019FE8F2513410080B7CCD81606410080030390251341008032CCD81606410000ED0790251341000035CBD81606410000D40C902513410080BEC9D81606410040B511902513410000CFC7D816064100408F1690251341008065C5D816064100005F1B90251341008082C2D81606410080222090251341000025BFD81606410080D7249025134100004DBBD816064100807B29902513410080FAB6D816064100800C2E9025134100002DB2D816064100C08732902513410080E3ACD81606410000EB369025134100801EA7D81606410000343B902513410000DEA0D81606410080603F902513410080209AD816064100406E43902513410080209AC015064100406E43302613410080FC92C015064100004F473026134100008B8BC01506410040F64A302613410000D083C015064100C0634E302613410000D17BC0150641008097513026134100009273C015064100409154302613410000186BC015064100C050573026134100806762C01506410000D6593026134100808559C01506410000215C3026134100007650C01506410000315E3026134100003E47C015064100800660302613410000E23DC01506410000A1613026134100006734C015064100800063302613410080D12AC015064100C024643026134100002621C015064100800D653026134100006917C015064100C0BA653026134100809F0DC015064100402C66302613410000CE03C015064100006266302613410000F9F9BF15064100C05B6630261341000026F0BF1506410040196630261341000058E6BF15064100809A6530261341008095DCBF1506410040DF64302613410080E1D2BF1506410080E76330261341000042C9BF15064100C0B262302613410000BBBFBF1506410040416130261341000051B6BF1506410080925F30261341000009ADBF1506410080A65D302613410000E7A3BF15064100407D5B302613410080F09ABF150641008016593026134100002A92BF15064100C071563026134100009889BF15064100408F533026134100003F81BF15064100006F503026134100802379BF1506410040104D3026134100006271BF15064100407E49302613410080136ABF1506410080C5453026134100803863BF1506410000E841302613410000D15CBF1506410080E83D302613410080DD56BF1506410000C9393026134100805E51BF15064100008C35302613410000544CBF15064100C03331302613410000BE47BF15064100C0C22C3026134100009E43BF15064100003B28302613410000F33FBF15064100009F23302613410000BE3CBF1506410000F11E302613410000FF39BF1506410080331A302613410080B637BF15064100C06815302613410000E535BF150641000093103026134100808A34BF1506410080B40B302613410080A733BF15064100C0CF063026134100003C33BF1506410000E7013026134100804833BF1506410080FCFC2F2613410080CD33BF150641000013F82F2613410000CB34BF15064100002CF32F26134100804136BF15064100C04AEE2F26134100003138BF15064100C070E92F26134100809A3ABF1506410000A1E42F26134100807D3DBF1506410080DDDF2F2613410000DB40BF150641008028DB2F2613410000B344BF150641008084D62F26134100800549BF1506410080F3D12F2613410000D34DBF150641004078CD2F26134100801C53BF150641000015C92F2613410080E158BF1506410000CCC42F2613410000225FBF15064100809FC02F2613410080DF65BF15064100C091BC2F2613410080DF65D716064100C091BC8F2513410080036DD71606410000B1B88F25134100007574D716064100C009B58F2513410000307CD716064100409CB18F25134100002F84D7160641008068AE8F25134100006E8CD716064100C06EAB8F2513410000E894D71606410040AFA88F2513410080989DD716064100002AA68F25134100807AA6D71606410000DFA38F25134100008AAFD71606410000CFA18F2513410000C2B8D71606410080F99F8F25134100001EC2D716064100005F9E8F251341000099CBD71606410080FF9C8F25134100802ED5D71606410040DB9B8F2513410000DADED71606410080F29A8F251341000097E8D71606410040459A8F251341008060F2D716064100C0D3998F251341000032FCD716064100009E998F251341", std::back_inserter(wkb)); 
        bg::read_wkb(wkb.begin(), wkb.end(), buffer); 
    } 
    bg::correct(parcel); 
    bg::correct(buffer); 

    std::vector<polygon_type> pieces; 
    bg::difference(parcel, buffer, pieces); 

    std::vector<polygon_type> filled_out; 
    bg::difference(parcel, pieces.back(), filled_out); 

    std::cout << std::setprecision(16) << bg::wkt(filled_out.front()) << std::endl; 
    std::cout << bg::area(filled_out.front()) << std::endl; 
    std::cout << bg::perimeter(filled_out.front()) << std::endl; 
}

We are using WKB here to not loose precision. I run this with three different numeric types (float, double and ttmath high precision). And that is the essence of Boost.Geometry! Not one type! Many types! All types!

int main(int, char* [])
{ 
    test_difference_parcel_precision<float>(); 
    test_difference_parcel_precision<double>(); 
    test_difference_parcel_precision<ttmath_big>(); 
    return 0; 
}

Using these three types we get:

float: one spike, perimeter 136.28
double: two spikes, perimeter 151.793
ttmath: zero spikes, correct result, perimeter 92.50548050416428971188559805054187252

Conclusions

SQL Server is not perfect, it can cause spikey features
PostGIS is not perfect, it can cause spikey features
The cause is not WKT, nor it is WKB
Boost.Geometry, using double, is not perfect because it causes the same spikey features
Boost.Geometry, using ttmath, gives the correct results
Boost.Geometry, using float, gives one spikey feature
Changing the cutline a bit causes only one spikey feature within SQL Server, still two in PostGIS, and zero in Boost.Geometry (with all numeric types).
Be aware of floating point precision

Free Shapefiles of Countries of the World (4)

2011-01-09T23:59:00.005+01:00

Free Shapefiles of Countries of the World (4)

There was the suggestion to look at gadm (Global Administrative Areas), thanks for this suggestion. This blog handles this new source in detail. Note that some sources are apparently hard to find, i.e. not found immediately by googling to "world country shapefiles" or similar. I'm glad that this blog appears now in the first 10 results, so hope to provide by this an index to world country shapefiles.

Look also at part 1, part 2, and part 3 of this blog.

Large

This newly found shapefile is very interesting. The first observation is that it is rather large, the zipfile is 158MB and takes several minutes of downloading. Unzipped the triplet (shp,shx,dbf) is almost half a gigabyte. When we zoom in, we see how this is caused, it is a shapefile taken from a grid:

The image above is a piece of the island Elba, the magenta line is the kml-ified shapefile displayed in Google Earth (yes, I like Google Earth).

I don't know if I like all those cornery borders... At least, it is not necessary to keep them, why are they not dropped or generalized away?

Anyway, if we zoom out a little again, the borders seem to be reasonably accurate at the first sight, here is the whole island Elba:

Vatican City and Liechtenstein

In last blog I described the small countries in and around Italy. So I shortly redo this for this new datasource. The gadm contains really the most accurate approximatino of Vatican City, compared with Google Earth as still our Single Point of Truth...

It is much better than all other sources found, and it therefore raises an interesting question: what is that little block where it differs, just to the southwest of Saint Peter's Square? Who is right there? And who is wrong? Interesting questions in the neighbourhood of the heart of the Roman Church... But we mean it geographically. Is Google failing here, as in Nicaragua and Costa Rica, nearly causing a war? Or is it the gadm.org website, a particular initiative providing freely available shapefiles for academic and non-commercial use (e.g. blogging)? According to Wikipedia, there indeed is a border dispute near Vatican City, but that is a little more to the east...

And look, we see a similar issue with Liechtenstein, probably a disputed valley...

But for the rest also here this shapefile is relatively accurate. For its size, accuracty is expected, of course...

Monaco

In last blog we saw that the borders of Monaco were correct in only one of the compared shapefiles... Alas, wrong again here:

It seems hard to create a correct shapefile... What is the source of all those errors? Who is really right here? Of course Google! For all certainty I checked this at Wikipedia (because I've never been in Monaco). Yes, the yellow line of Google is following Monaco's portrayal within the Wiki page of Monaco. So alas, gadm is totally wrong here.

The Netherlands

I'm living there, in that little country which is so difficult to be stored in free shapefiles... There is really a large lake, only bordered from the sea by a dike. That lake should be displayed, and it is not. If very small and uninhabited islands are displayed, the lake should really be there. And it is not the only error, in Zeeland the (smaller) water bodies behind the dikes are omitted as well.

The Caspian Sea

This was another requirement (CS1). It is not as bad as it was in some other shapefiles, where the Caspian Sea was divided between surrounding countries. But the Caspian Sea is included as a separate country in gadm... Interesting. Do they issue passports?

Validity

Finally there is validity. All countries should be modelled as multi-polygons and be valid in both SQL Server and PostGIS. We've already seen that this requirement is met in only one caes, until now (11). And also for this case, 13, alas, there are many multi-polygons invalid in SQL Server. In PostGIS it is OK, though the query runs for a long time due to the excessive size of the geometries. The values are post-editted in the previous blog.

Actually it is nice to show a variant of the query, using the with syntax:

with valid as (select ST_IsValid(the_geom) as v from world13)
select v,count(1) from valid group by v;

Conclusion

As always I'm glad with new shapefiles. This one is quite accurate, and in many aspects the best until now. But also, it is way too large, too cornery, and contains errors, as all others do...

So I think a good and freely available shapefile of the countries of the world is still very welcome.

If there are more shapefiles, I'm still interested, but especially I'm interested in a the Google country borders in any format? Are they available?

Free Shapefiles of Countries of the World (3)

2010-12-27T23:28:00.008+01:00

Free Shapefiles of Countries of the World (3)

This is part three of the comparison of free shapefiles of countries of the world. After feedback, new sources, blogs and tweets (thanks everyone) on my previous blogs (part 1 and part 2), I decided to show some things more.

Edited: Note that part 4 handles an additional shapefile, referred as 13 here.

Let's now look at the north of Italy, where several countries are situates: Italy itself, France, Switzerland, Austria, Germany, Slovenia, Croatia, and some very small countries:

San Marino
Liechtenstein
Monaco
Vatican City

And I can't really believe what I saw today. Vatican City has been moved by some shapefile creators...

Vatican City

With respect to Vatican City, all shapefiles that I found are wrong. Below an image created with Google Earth showing boundaries. The yellow one is Google's, so not from a shapefile. It is accurate. I really want a shapefile with Google's country borders! Didn't found it until now, even not in KML-format. All others displayed, 3, 5, 6, 7, 9 are not accurate, but acceptable. I know that Vatican City is a very small country, that on a worldscale map you could not expect too much detail, so it might make sense that they display the country border only roughly. So these are all reckoned sufficient... with respect to the requirements, see more about that below.

But I was really amazed by the other ones, in the zoom out below. Those (1, 2, 8, 12) are all wrong. (Remember that all numbers refer to numbered sources from previous blog, completely listed below). About 1: we already knew that it was far too rough, and it is actually understandable: it exaggarates small countries to show them on a worldscale map. Understandable. But 2, 8 and 12 are completely wrong. Especially from 12 (source: eurostat, europa.eu) you would not expect this. What has happened here...

Some sources (4, 11) do not have Vatican City at all.

Monaco

For Monaco it is a little better, in the sence that countries are not moved by 24 kilometers. But most shapefiles are more or less wrong, 12 being the only exception.

The Yellow line is Google again, and 12 is the cyan one following it. All the rest is wrong.

Liechtenstein and San Marino

The other small countries, San Marino and Liechtenstein, are both displayed more or less acceptable in all but 1, 2 and 8 (all omitted in the pictures - 8 is for Liechtenstein OK). Number 12 (cyan) also here follows Google the most closely.

Based on this pictures, 12 (from Eurostat) is following Google's country borders the most closely, with an amazing exception being Vatican City, where it is totally wrong. The rest of the sources shows not a real winner here.

Elba

Not a country, but also compared, the island Elba is another sample where the vertices of source 12 (Eurostat) follow the borders of the island quite closely. Cyan is the 3M version, blue the 10M, darker (and paler) blue the 20M. In the 60M version Elba is omitted. I agree with this all. The other borders, magenta here below, are all worse. And I omitted 1,2 and 8 here because they are even more wrong.

Requirements

In the last two blogs I made remarks about the quality, but I didn't really listed them as formal requirements. So let's start with that now. They are composed on the fly, and maybe a bit subjective, I know that.

NL1: The five inhabited islands in the Waddensea in the Netherlands should be present
NL2: The inner lake in the Netherlands should be included
NL3: The polders in the Netherlands should be visualized correctly, so absence of non-existing polders and presence of polders existing since 1968 (sometimes you have to make requirements quite explicitly...)
CS1: The Caspian Sea should be there
DB1: Countries should be one multi_polygon per country, not separate polygons
DB2: All (multi)geometries should be valid within SQL Server 2008
DB2: All (multi)geometries should be valid within PostGIS
SM1: Vatican City should not be moved and (probably) not be omitted
SM2: Monaco's borders should follow real borders
SM3: San Marino's borders should follow real borders
SM4: Liechtenstein's borders should follow real borders
SM5: Elba's borders should be followed

In the table below the requirement list is crossed against the shapefiles.

NL1

NL2

NL3

CS1

DB1

DB2

DB3

SM1

SM2

SM3

SM4

SM5

Enlarged

Very poor

26 km off

Very poor

Sufficient

Very poor

Sufficient

Poor

Sufficient

Good

Sufficient

polygons

Omitted

Poor

Sufficient

Poor

Sufficient

Poor

Sufficient

Poor

Sufficient

Poor

Sufficient

Poor

Sufficient

Poor

Sufficient

Poor

26 km off

Poor

Sufficient

Very poor

Sufficient

Poor

Sufficient

Poor

Sufficient

Poor

Omitted

Sufficient

Omitted

omitted

Poor

Sufficient

Good

Omitted

12 (3m)

24 km off

Excellent

12 (10m)

24 km off

Good

12 (20m)

24 km off

Good

Sufficient

12 (60m)

24 km off

Sufficient

Poor

Omitted

Edited: 13

as separate
country

Excellent

Poor

Excellent

All sources compared by image

The images below show the same picture, using all shapefiles, as we did in the previous blog, but now showing these four countries, Italy, and parts of the surrounding countries.

	(1) from aprs file world.shp
	(2) from brothersoft, serving blue marble data file World_countries_shp.shp
	(3) from Geodan file WorldCountries.shp
	(4) from mapping hacks file world_borders.shp
	(5) from thematic mapping file TM_WORLD_BORDERS-0.3.shp
	(6) from openmap here, called cntry02 file cntry02.shp
	(7) from this site ESRI data file cntry08.shp
	(8) from vdstech via this blog file world.shp
	(9) from natural earth file 10m-admin-0-countries.shp
	(9b) from natural earth file 50m_admin_0_countries.shp
	(9c) from natural earth file 110m_admin_0_countries.shp
	(10). Empty
	(11) from huebler file world_adm0.shp
	(12 1:3m) from eurostat file cntr_rg_03m_2006.SHP
	(12 1:10m) from eurostat file cntr_rg_10m_2006.SHP
	(12 1:20m) from eurostat file cntr_rg_20m_2006.SHP
	(12 1:60m) from eurostat file cntr_rg_60m_2006.SHP

Conclusion

Again thanks to commentor Alectora who suggested Eurostat, it is a nice set of four datasets in different scales. The 1:3 million version is excellent in most measures, but it is a pity that Vatican is displaced and the IJsselmeer in the Netherlands is omitted. And that some polygons are not valid...

I still think a good and free shapefile of the world is still welcome.

Doxygen and QuickBook

2010-12-24T22:48:00.006+01:00

Doxygen and QuickBook

Doxygen and QuickBook are both great tools for generating documentation from C++ code. But their cooperation is quite complex... let's say, a challenge.

Quickbook

QuickBook (QBK) is based on BoostBook, which used to be the documentation tool for Boost Libraries. BoostBook is on its turn based on DocBook.

With QuickBook (and also BoostBook) you can make good looking library documentation describing headers, examples, syntax highlighting, etc. Look at the QBK website for more information. For example, I like the callouts.

QuickBook input files (extension .qbk) are plain text files with Wiki style commands. However, when using QuickBook alone it is not maintainable to document hundreds of classes and functions, you really should have a link to the source code. And that is what Doxygen is doing.

Doxygen

Doxygen is a standalone tool, accepting C++ source code (it accepts other languages as well) with JavaDoc style comments to document function parameters, etc. Look at the website for a manual. While Doxygen does a great job, the generated documentation is often considered a bit too technical. In the review report was included:

Using Doxygen alone as a documentation tool has its (known) shortcomings when it comes to generic libraries.

And the individual reviews mentioned things as : "I think doxygen is convenient for developers but not very useful for users, because of the high structure/content ratio." and "I really did not like the doxygen generated content. I found it difficult to find the real documentation among all the doxygen generated fluff." and "I have a general problem with doxygen generated docu with regard to generic libraries. It not easy to separate the important classes and fuctions from the implementation details. GGL helps with this through Grouping by modules and a overview Page, but it could be even better. I think the Modules should them self be further grouped, and the most important should be highlighted.".

This does not mean that Doxygen is wrong or bad, and maybe we did not spend enough time to configure it right. An individual review also mentioned: "Doc: Excellent. Looks very nice. Only an index is missing - but extensive Doxygenation including Doxygen comments are far above what passes for documentation in many packages". We did our best at least.

Anyway, people do like QuickBook generated documentation, so we, with Boost.Geometry, are moving to QuickBook. But we want to keep Doxygen there to link with the source code.

Doxygen -> QuickBook via XSLT

It is certainly possible to let Doxygen and QuickBook work together. Doug Gregor wrote here:

(...) If I haven't scared you yet, I'll put it in a diagram:

C++ sources -> Doxygen XML -> BoostBook XML -\
-> Merged BoostBook XML
Quickbook sources -> BoostBook XML ----------/

Then, for the output:

Merged /--> HTML
BoostBook --> DocBook XML ------> FO --------> PDF
XML \--> Man pages

Look complex? It is. (...)

So it is complex. Especially the XSLT file to convert from Doxygen XML to BoostBook XML proved complex and sometimes incomplete, we (Mateusz) spent time on it and achieved good documentation, but it was really complex. And in the end, you don't get QuickBook, you get BoostBook... Suppose you want something different, for example adding "complexity" to the documentation... It means changing the XSLT files again.

What is in the docs?

I like to see in documentation, per function or per class:

What it does
The header where the function or class is found
Its parameters (for a function)
Its return type (for a function)
Its limitations
An example (this is very important)
An image (very important for a geometry-library, probably not for all libraries) to show what it is doing
Behaviour in various combinations
Its complexity

The cplusplus documentation has for every function or method a little example, I like that. It is great for programmers who program by example. Examples need to be compilable examples, they need to be extracted from real compiled code, otherwise old constructs or other obsolete things or typos will stay there unnoticed. You see that often.

About behaviour: within Boost.Geometry, one function (e.g. distance) can do many things (distance between two points, between point-and-polygon, etc). We also have multiple dimensions (2D, 3D) and coordinate systems. The behaviour of the functions might be different in different circumstances. Therefore, we want to have a possibility to document these differences explicitly.

Besides this, we often have two overloads, one with a so-called strategy, giving the library user the tools to change default calculation methods, or one without, so sticking to the defaults. Besides these differences, these overloads are the same. So documentation is the same. Where does the documentation come from? From the sources. Do we want to have the same documentation duplicated for two functions? No. Or probably not.

Don't repeat yourself, and macro's

Doxygen does have an overload option, but it seems to work only for overloaded member functions. We are overloading free functions, which are often heavily templated. I didn't have success to get Doxygen's overload construct running with our code.

Anyway, Doxygen provides a great utility: an alias. This is a sort of a macro, expanding to whatever you specify. So when we think we repeat ourself (e.g. writing for a template parameter "Coordinate system (e.g. cs::cartesian)", which will be found at many places) we can introduce an alias there, param_macro_coorsystem, and we can use it in the documentation.

The result looks cryptic, but at least you don't repeat yourself and if you want to change that sentence, you do it in one place.

To make the picture complete: QuickBook also provides a great utility: a macro. This is doing about the same: it is a placeholder which is replaced by quickbook by the actual text.

So we can actually chose: write plain English, take an alias, take a macro. Or mix them.

Repeating yourself is usually not so good, but in documentation, it cannot always be avoided. Refactoring all repetitions from comments and document fragments will result in a very cryptic developer-defined language. Actually, it is already looking cryptic now, read on...

Doxygen -> QuickBook via parser

As described above, the XSLT process to go from Doxygen to BoostBook was not completely satisfactory for us, with all our requirements to documentation. It probably can be done, but requires high XSLT skills, and who nowadays still writes XSLT? There are a lot of libraries parsing XML, in C++ (and in other languages), and for a C++ programmer it is quite easy to parse the Doxygen XML and output QuickBook. So that is what is done... We developed a new tool, doxygen_xml2qbk. It goes from the XML's generated by Doxygen to QBK.

The tool might be more useful in general (for more Boost libraries), but currently it does its job for Boost.Geometry.

It is currently (lightly) based on RapidXML (great library, by the way, and super fast).

The markup

So we document functions like this, I take here the area function:

/*!
\brief \brief_calc{area}
\ingroup area
\details \details_calc{area}. \details_default_strategy

So these lines use macros to say that it is calculating the area (brief), and again that it is doing that (detail). The differences in brief and detail are Doxygen differences, sometimes convenient, sometimes not. Then a remark-macro is added that this function takes the default strategy. That is the case with a lot of functions, therefore it is a macro, to avoid repetition. Then it continues:

\tparam Geometry \tparam_geometry
\param geometry \param_geometry
\return \return_calc{area}

So here the parameters are described. There is a template-parameter (Doxygen's tparam) and a parameter (Doxygen's param). Those are similar, geometry! Not literally the same, but the parameter geometry (lower case) is of the type Geometry (Camel Case). This is the case for most of the functions, accepting sometimes one, sometimes two geometries. Therefore macro's (\tparam_geometry) are introduced here... Below we will see that they can be handled in the same line! Though not everyone is convinced of that feature. It might go away or be made flexible.

Avoiding repetitions: it the same for the return value, after Doxygen's \return. It would be possible to write here: "the function area calculates the area of the input geometry". For the function length we will write (about) the same comment (replace area with length). For perimeter the same. Et cetera. So, to avoid repetition, we created the cryptic macro \return_calc, with a parameter a literal string (in this case "area").

We continue:

\qbk{example,area_polygon}
\qbk{example,area_polygon_spherical}

So this is a Doxygen Alias (\qbk) to generate an XML node called <qbk.example>, it generates two of them, one per line. Such an XML node is recognized by our new tool doxygen_xml2qbk. And that tool generates, as expected, an example section. So these two lines generates this QuickBook syntax:

[heading Examples]
[area_polygon]
[area_polygon_spherical]

Where area_polygon can be found in one of the examples, using a QuickBook construct to show a syntax highlighted example.

We continue:

\qbk{behavior,__0dim__:[qbk_ret 0]}
\qbk{behavior,__1dim__:[qbk_ret 0]}
\qbk{behavior,__2dim__:[qbk_ret the area]}
\qbk{behavior,__cart__:[qbk_ret the area] __cs_units__}
\qbk{behavior,__sph__:[qbk_ret the area] __sph1__}
\qbk{behavior,__rev__:[qbk_ret the negative area]}

This all goes to a xml.behavior node, and then to a QuickBook behavior section, a table describing the behaviour for 0 dimensions (point), 1 dimension (linear), etc. All this is macro'd, where __sph__ is also a QuickBook macro, and [qbk_ret] is also a QuickBook macro. They can be nested as well. qbk_ret just says "returns". It has the same length, the only advantage is that you don't repeat yourself in words but now in macros.

OK, using two macro systems through each other, some having a specific meaning for our converter tool, this is quite complex and I realize that. But on the other hand, if everything is written in words here, it is probably unmaintainable.

The results

Let's finally show the results here, as screen-dumps.

Doxygen generates this piece:

To be honest, we must say that examples (using syntax highlighting) within Doxygen are certainly possible as well.

With Doxygen-conversion-QuickBook the following piece is generated:

So we see what our \qbk alias does here... It generates the sections behavior, complexity, examples, and might generate more.

I like it, despite its complexity.

Intersections (2), "recursive polygons"

2010-12-22T22:53:00.006+01:00

Intersections (2), "recursive polygons"

The Boost.Geometry review report said:

Testing: several reviewers mentioned the need for a thorough testing framework allowing the verification of the correctness of the algorithms in a wide range of use cases. Different test strategies need to be employed, such as high volume and random tests, known border case tests, tests using different numeric precision types, etc.

This blog will discuss a part of this, a program called recursive_polygon.cpp, which is a test in the category high volume an random, using different numeric precision types.

Recursively made polygons

In this testprogram, polygons are created in a recursive way. In the first step two polygons are created. They are created at a random place and either a box, or a triangle (in this case: a box with one coordinate omitted). The first version of the program was called recursive_boxes because it then did only boxes. Those two polygons are unioned and the result is either a polygon, or a multi-polygon. In the first step, it is probably a multi-polygon.

After this, two other polygons like this, based on random coordinates, are unioned and result in another multi-polygon.

In the second step, the results of these two earlier unioned multi-polygons are again unioned. This will deliver a third multi-polygon. And so the process goes on, each time creating more complex multi-polygons, and after a a while holes are generated, self-tangencies.

The figure below shows the idea more clearly.

In this figure, a field of 3x3 is used and intersections are done up to level 3.

The test and recursive structure and sequence can be compared to a genealogical Ahnentafel, the unions to marriages, the number 15 to the proband, having two parents, four grandparents, eight great-grandparents, et cetera.

Checking results

While we can run the program and see the results (the program optionally creates an SVG file), we have to have a mechanism to check if the results are correct.

Luckily we are doing geometry, we are doing mathematics, we are doing set theory, so we can use that to check our results.

Checks are done in each step, so during the processing of two polygons. Besides a union (u), an intersection (i) is done. The area is calculated of the original polygons (p and q), of the union, and of the intersection. And it always must be that the area of the union equals to the area of both input polygons, minus the area of the intersection. So: A(u) == A(p) + A(q) - A(i). A simple check, all done by the library itself.

So we can run the check for thousands of times, and for various levels, field sizes, boxes and triangles. It now also runs for counter clockwise polygons and open polygons.

It is not the case that it was without errors the first time... There have been many errors, especially in the self-tangencies, points were two polygons in a multi-polygon touch each other. That all have been solved. So running this test was absolutely valuable.

Seven levels

After seven levels in a field of 10x10 we might get the next results:

At the left the two input polygons (green and blue), at the right the union in red. After these seven levels in the 10x10 field, the unions are that large that they nearly completely fill up the whole area.

Last month I did most tests in this configuration, 10x10 and up to 7 levels. 7 levels results in 255 tests, this is 2⁸-1

Twelve levels

Twelve levels results in 2¹²-1 unions and shows in a 100x100 field like this:

and we can go further and further, but for Boost.Geometry the tests do not give problems anymore. It all runs fine. This test (so 8191 unions and intersections ending in this complexity) runs in two seconds on my current machine (actually, to be precise, it is doing the union twice because the overlay function is reused in other tests as well).

Unit test

It is a sort of a unit test, because the program checks itself. I also use it as a unit test. But note that the program is based on random input, and that a 30-level test resuls in 2³⁰-1 unions and intersections (~1G). So it runs for a while, and including this in the standard Boost test suite will not be appreciated. So this recursive_polygons can be run manually, and there are parameters to specify the number of levels, et cetera.

By the way, there are many real unit tests within Boost.Geometry, there are currently 118 source files using the (great!) Boost.Test library, and some of them use polygons which are created by this testprogram as their input.

Conclusion

This recursive polygon unit test has been very valuable and, besides that, such tests are nice to program. The whole program, links are given above, is rather short, and it is of course Open Source. This test might be useful for other libraries as well.

Intersections (1)

2010-12-19T18:14:00.006+01:00

Intersections (1)

The boolean operations are probably the most interesting parts of the Boost.Geometry library.

Boolean operations are, based on two polygons:

intersection, returning polygons containing areas present in both input polygons
union, returning polygons containing all areas present in either one of the input polygons
difference, returning polygons containing areas present in one of the input polygons, but not in the other
symmetric difference, returning polygons containing areas present in one of either input polygons, but not in both

In other words, boolean operations are spatial implementations of mathematical set theoretic operations. Based on two polygons A and B:

intersection: A and B, also written as A ∩ B
union: A or B, also written as A ∪ B
difference: A and not B, also written as A \ B
symmetric difference: A xor B, also written as A ∆ B

The difference operation (also known as relative complement) is the only operation of these which is not commutative, A \ B results in another geometry than B \ A.

The basics

This blog will not explain the whole algorithm, that would be way too long. Other of my future blogs or a future article will explain more details. The algorithm used in Boost.Geometry is an adapted version of the Weiler-Atherton algorithm. The basics are quite simple, and explained here, and also in adapted versions as from Greiner-Hormann. Note that I have adapted the algorithm myself, such that intersection points are not inserted in the originals but stored separately. Which is essential because we don't want to change the input polygons, and we neither want to copy the input polygons. I have made more adaptations, such that all situations including self-tangencies (in some articles called degenerations) and robustness are handled.

So only the very basics in this blog. There are two input polygons, A and B (here: green and blue), both oriented clockwise. The arrows indicate orientation at the starting points. Intersection points (orange) are calculated. During the calculation of the intersection points, traversal information is gathered: for an intersection (usually: go to the right), would it continue following the green or the blue outline... For example, at point 0, above, it is figured out that for an intersection the green line should be followed, and for a union the blue line. Both until the next intersection point is encountered. That requires sorting which is used heavily within the algorithm.

Turning right for intersections or left for unions is also visualized below, from an older part of the GGL-doc:

How traversal information is figured out, and how intersection points are sorted, will be explained in another blog, maybe much later.

So intersection points and turn information are basically all that is necessary up to this point. The traverse function uses this information, and adds points to the accumulated output polygon. Added points can be either intersection points or original vertices. The output intersection then looks like this (slightly shifted by hand to show it better):

Note that it can become (often much) more complicated, holes, intersections at vertices, self-tangencies, polygons-within-holes, multi-polygons, etc.

Here the process for the intersection of a polygon with a polygon with two holes:

Note that the numbers of the intersection points are shown in the order found, and not in the order of traversal. It is not possible to indicate intersection points with an order of traversal, because that order differs for intersection and union (check this in the picture above). The (manually added) red arrows depict the traversal for intersection, starting with the red "0".

Note also that holes are not a big addition to the algorithm, if they are oriented counterclockwise in clockwise polygons (and vice versa). This is an assumption that is made within Boost.Geometry. If this is the case, traversal can take its routes always in the same way. The union operation would create two holes, which are different from the original holes and stored automatically in the correct orientation.

All boolean algorithms are the same

I mean to say here: intersection, union, and (symmetric) difference are all the same.

We have seen above that for intersection and union the same information is gathered. Only one variable, the traversal-direction, determines what is the output. So easy, we have one algorithm in which we calculate intersection or union.

Now for the difference, it turns out that if we reverse one of the polygons (e.g. green one), such that it is oriented counter-clockwise, the difference is calculated automatically by calculating the intersection of that reversed A and the still normal B. That is not that strange, because it follows from a mathematical rule. The difference of A and B is the same as the intersection of the complement of A with B. And the complement of A (A^c) is just its reverse. Voila, we do another function with the same algorithm.

The complement (from older GGL-doc) of polygon A is the whole world but A.

The symmetric difference, A xor B, is just adding the two differences: it is also defined as A and not B unioned by B and not A.

Concluding this: all operations intersection, union, difference and symmetric difference are calculated with this one algorithm. The algorithm consists of several sub algorithms, relevant for all operations, and I hope to describe them later in more detail.

Reversing

How do we take the reverse in the algorithm? In the past, I cheated and reversed the polygon before doing the algorithm. That was easy, but has two drawbacks: 1) the input must be copyable and 2) it takes time to reverse. So we don't do this anymore. We just iterate through it in reverse way. That is easy enough, using an iterator walking reversely though the vertices. Or, in the case of Boost.Geometry, we use a reversible_range, as explained in previous blog. It is all implemented using specializations and meta-programming.

The difference is shown below. On the left side the result, again shifted slightly to make things more clear. The right side shows the same but in a different way: the intersection of the complement of A (green) with B (blue).

The one algorithm handles all geometries

We have seen above that one algorithm is enough and that we can walk in reverse through the polygons to support differences. What if we have counterclockwise polygons? The answer is the same, we walk through it using our reversible ranges, which are in fact Boost.Range reverse_range adaptors. So we calk backwards through it during getting intersection points, and later when vertices are necessary for the output, and actually everywhere where it is used.

So we can handle clockwise and counterclockwise input polygons, and we can even handle one clockwise, one counterclockwise, and the output can also be clockwise or counterclockwise at wish. That last feature is done by appending points during traversal to the back of the output rings, or inserting them at the front.

The algorithm also handles combinations of rectangles and polygons, because a rectangle can be seen as a polygon. The only difference is that the calculation of intersection points is done more efficient. Furthermore, multi-polygons can be handled. Note that the intersection of two polygons can result in a multi-polygon.

So this seems a litany of what the algorithm can handle, and it does not stop here. Because Boost.Geometry algorithms are independant on their input point types, and also independant on their input polygon types, and other geometry types, the algorithm can handle any geometry type based on any point type. So it can handle 4-byte-floats but also 8-byte-doubles or high-precision types like TTMath. It can handle boost::geometry polygons but also ESRI polygons, provided those polygons are adapted to Boost.Geometry.

Another aspect of this algorithm is that it is coordinate system agnostic. It can handle cartesian polygons, but can also calculate intersections of spherical polygons. This is possible because the algorithm is based on sub-algorithms as calculating the intersection points, and finding the side (left/right) of a point with respect to two other points. If those algorithms are provided (currently they are not...) the intersection algorithm will handle them.

The fastest...

Finally, the implementation of boolean operations in Boost.Geometry is also the most efficient intersection algorithm which can be found. If people don't believe this, or want to test it, or help benchmarking: there are Open Source benchmarks. Other people are more independant as I am, help would be welcome.

Below intersection of two multi-polygons with self-tangencies, intersection is again shifted slightly. This comes from one of the robustness tests.

Final remarks

Using the reversible_range for (symmetric) difference instead of reversing geometries, in the implementaiton, was done yesterday. So all this stuff above was planned but, until yesterday, not implemented like it is today. So another step forwards is made.

The review report (from now more than a year ago) said:

Boolean operations: while the library provides a set of Boolean operations those seem to be not complete and/or robust in terms of clockwise/counterclockwise geometries, closed/open polygons. Robust Boolean operations are a strong requirement, so this should be fixed as reported by at least one reviewer.

This point is now dealt with. There are still a few tweaks to be done but they will disappear coming weeks. I didn't discuss robustness in this blog, and open polygons neither, but these issues are already addressed too.

Don't think you can implement boolean operations easily, though the basic algorithm is quite simple... There are about 100 cases for extraction of turn information from two intersecting segments. Besides that there are also dozens of cases necessary for self-tangencies, in GIS not so ubiquitous, but in games certainly important. I started the algorithm in January 2008 and rewrote it two times after that, one time before the review, and another time after it because there were issues with the implementation, self-tangencies, and robustness. Now it should be ready for production.

The implementation is heavily based on C++, generic programming, templates, specialization and meta-programming. I don't think it is easily translatable into another language (but D). Also in C#, where specializations (at runtime) are possible according to one of my previous blogs, metaprogramming is not possible, it is a compile-time mechanism.

Range adaptors

2010-12-13T22:57:00.002+01:00

Range adaptors

I seem to like doing blogs in series... (Or ranges ;-) .) See my previous blog about ranges here. This is a short bonus.

Since Boost 1.43 released in May 2010, just before BoostCon'10, the Boost.Range library has Range Adaptors. Range adaptors enable users to change range behaviour on the fly, for example traverse it in reverse order.

Reversing through geometries

See this small and working snippet.

using namespace boost;
geometry::model::linestring<P> ls;
geometry::read_wkt("linestring(1 1,2 2,3 3,4 4)", ls);
std::cout << geometry::wkt(ls) << std::endl;
std::cout << geometry::wkt(ls | adaptors::reversed) << std::endl;

This will print out:

LINESTRING(1 1,2 2,3 3,4 4)
LINESTRING(4 4,3 3,2 2,1 1)

So the original one and the reversed one. Works great, and great syntax.

To let this work, I had to adapt the Boost.Range reverse_range to Boost.Geometry, such that Boost.Geometry knows the adapted range is a linestring. In this case it is a linestring, but it might also be a ring (also a range) or a multi-point (also a range), etc. So I adapted it like this:

namespace traits
{

template<typename Geometry>
struct tag<boost::range_detail::reverse_range<Geometry> >
{
typedef typename geometry::tag<Geometry>::type type;
};

}

Just meta-programming, if a reverse_range is specified, it defines its tag as it is for the unreversed geometry. This should work for all these geometries because all of them are very light-weight, essentially only requiring the tag metafunction in namespace traits. (Well, one exception, the linear_ring (ring for short) also might optionally define point order and closure.)

We do the same for the other Boost.Range adaptors, where relevant, and then have filtered, sliced, strided, uniqued, reversed all out of the box. Great! I already liked Boost.Range, and this new feature is cool either.

Why ranges are cool, revisited

We see here again why ranges can be preferred above iterators, and make some other statements:

First we state in another way (than in the previous blog) that ranges are like normal objects:
Ranges are first class citizens, while iterators should be implementation details

And then we make a statement about the cool adaptors we encountered in Boost.Range:
A range allows overloading operators, in contrast with iterators

reverse_view

Before ranges were there in May 2010, I developed a small utility class called reversed_view, where ranges are traversed in either forward or backward direction.

This reversed_view class is still there, but now greatly simplified, the reverse one expressed in terms of the reversed_range adaptor. It is now defined like this:

enum iterate_direction { iterate_forward, iterate_reverse };

template <typename Range, iterate_direction Direction>
struct reversible_view {};

template <typename Range>
struct reversible_view<Range, iterate_forward>
{
typedef identity_view<Range> type;
};

template <typename Range>
struct reversible_view<Range, iterate_reverse>
{
typedef boost::range_detail::reverse_range<Range> type;
};

All meta-functions. Seeming empty classes, but of great value. In the first version it were not meta-functions but normal classes, derived from something or defining something. Therefore they needed an own constructor etc. These are only typedef's, do not define any code, therefore simpler and therefore better. In this case they exist because templated typedefs are not there.

closeable_view

Boost.Geometry also has a closeable view but alas I didn't find this in Boost.Range adaptors, for obvious reasons (the adaptors strip, replace or adapt, but do not add things). So the closeable_view is still there, but now also expressed as meta-functions, the actively-closing one is a real implementation, the already-closed on is another type definition.

Range, return types and references

2010-12-12T13:32:00.000+01:00

Range, return types and references

During BoostCon 2009 the keynote speaker was Andrei Alexandrescu, who wrote the book Modern C++ Design in 2001. Andrei held a fabulous presentation with the name "Iterators Must Go", in which he point out several serious drawbacks of iterators. Iterators have been the basic mechanism to use the Standard Template Library classes, so this speach was revolutionary, and enjoyable. Andrei not only critized iterators, he also gave an alternative: ranges. Ranges are for him the basic piece.

The keynote was the forenote for an article that Andrei wrote, now with the title "On Iteration". The title is changed into a more subtile, but the message is more or less the same. I quote: "STL iterators are marred by lack of safety, difficulty of usage, difficulty of definition" .

Before BoostCon '09 with its famous keynote, I already knew the Boost.Range library quite well. Most of the algorithms in Boost.Geometry are based on ranges. We use Boost.Range everywhere. Ranges are very nice. A std::vector is a range, but a std::pair of iterators is also handled as a range, and a Boost.Array is also handled as a range, and you can create your own ranges. Yes, like iterators, but ranges are more elegant.

Andrei Alexandrescu mentioned Boost.Range during his keynote, of course. He said that Boost.Range was a good start, but should be worked out way further. Until now, that has not been realized as far as I know. But even without changes, Boost.Ranges are far more convenient than iterators.

Iterators

People reading this blog will know iterators and otherwise will have stopped reading here. However, I like to explain the very basics shortly because of the point that Andrei made: "iterators are unnecessary difficult", I've always found that. This is how standard iterators work:

std::vector<int> a;
// fill a
for (std::vector<int>::const_iterator it = a.begin(); it != a.end(); ++it)
// do something with *it

When I first worked with the std::library, somewhere in 199X, I refused to do this and wrote:

for (int i = 0; i < a.size(); i++)
// do something with a[i]

This is half the size, easier to write, easier to read. But it is (much) less efficient and should therefore not be done like that.

Note also that C++ programmers always write ++it and not it++ because it is slightly (and detectably) more efficient.

Boost.Ranges

A Boost.Range is a sort of view on a std::vector, or on any other iterable instance. Boost.Range still defines and works with iterators. It is another level of abstraction above the std:: library. With iterators you iterate like this (still vector a):

for (boost::range_iterator<std::vector<int> const>::type it = boost::begin(a); it != boost::end(a); ++it)
// do something with *it

Hey, this is even much longer! Yes, in this sample this is true and in general it is true, but ranges are normally (at least within Boost.Geometry) used in a template environment, where it becomes:

for (typename boost::range_iterator<Range const>::type it = boost::begin(a); it != boost::end(a); ++it)

Not shorter but more generic: a Range can be anything.

There are shortcuts as well: with BOOST_FOREACH or BOOST_AUTO. They are macro's, so often avoided, but they make things more readable and can save template metaprogramming, see below. But first we will see why ranges are so convenient.

From Iterators to Ranges

Boost.Geometry supports geometry algorithms such as distance. They are defined in a generic way, so

template <typename Geometry1, typename Geometry2>
double distance(Geometry1 const& g1, Geometry2 const& g2)

(The real version does not have the double there but a metafunction). With distance users can calculate the distance between two points, or a point and a line (linestring), or a point and a polygon, etc.

The earlier version of Boost.Geometry (then called GGL), before going to ranges, could calculate the distance to of a part of a linestring by specifying two iterators as well:

template <typename Geometry1, typename Iterator>
double distance(Geometry1 const& g1, Iterator it1, Iterator it2)

This is great functionality, but it is completely inconvenient to implement:

all other functions are generic, having two arguments, why does this one has three arguments?
all other functions are reversable, the distance from a point to a line, and from a line to a point, both result into the same implementation. For iterator support we would need to create an overload having the iterators first
there are more versions (for example: with a strategy), we would need those overloads for those two, leading into another combinatorial explosion (there are more algorithms too...)

When we discussed the GGL previews on the list, Mathias Gaunard suggested to use ranges instead of Boost.Geometries and we took this over and were very glad with this change.

This is quite an argument for ranges instead of iterators in a template-library environment, and this is why I emphasis this at this point.

So let me state: Ranges are compatible with other objects, while iterator-pairs are not.

References

After this lengthy introduction we consider how iterators and ranges are passed.

Iterators are normally passed by value. The main reason for this is explained nicely in StackOverflow: by passing by value you can create local copies, like v.begin().

Ranges are normally passed by reference. Ranges might be non-copyable containers, so passing by reference is required. Besides the non-copyable issue, a std::vector is a range. You definitely don't want to pass a std::vector by value.

Boost.Geometry's polygon

Boost.Geometry contains models of geometries like linestring and ring. They are defined as std::vector or a std::deque, but handled in the code as a range. Any linestring or ring fulfilling the Range Concept can be handled by Boost.Geometry, at least that is the basic idea. A polygon is more more complex because it might contain holes, in the definition of OGC and Boost.Geometry. So a polygon has an exterior ring and zero or more interior rings, defining the holes. The exterior ring fulfills the Ring Concept, the interior ring is a Range of Rings, so a Range of Ranges, a collection or rings.

Boost.Geometry contains a free function exterior_ring which returns the exterior ring by reference. And it also contains a free function interior_rings where the collection of interior rings are returned by reference. Both functions have a const and a non-const version. Of course they don't return the rings by value: polygons can be huge geometries, I've recently processed a polygon of 6 megabytes (in WKT-format). They should never be returned by value.

Other polygons

With Boost.Geometry, geometry classes of other libraries can be adapted, such that Boost.Geometry can run algorithms on them as if it was one of its own geometry classes. There are examples showing this with libgd's points, WxWidgets's points, Qt's polygons, etc. With traits classes these geometry classes can be adapted.

However, until today, polygons were required to return the exterior ring and interior rings by reference. But what if those geometries do not use ranges?

We encountered this adapting Boost.Polygon, our sister library within the Boost family. Boost.Polygon defines a polygon-with-holes-concept, where iterators begin_holes and end_holes are returned.

So after successfully adapting Boost.Polygon's point, rectangle and polygon (without holes, so we call it a ring) to our library, I started with the polygon. It is no problem that Boost.Polygon do not use the range concept internally. Because we can build a sort of proxy, simulating a range. So the adaptation of a Boost.Polygon polygon-with-data consists of defining:

an iterator to iterate through holes
a ring-proxy with an iterator-pair (point-begin, point-end)
an interior-ring-collection-proxy with an iterator-pair (holes-begin, holes-end).

But there was one problem:

Returning local references?

Local variables can (of course) not be returned as references. So you can create ring-proxies, but they cannot be returned... That is to say, they cannot be returned as a reference. Of course, they can be returned by value.

There are several solutions to this problem:

Change the Boost.Geometry polygon concept to work with iterators, so let it return iterators. This would be quite a lot of rework, and sounds like going ten years back. So discarded.
Change Boost.Polygon's polygon concept. But that is of course not wished and not possible. Boost.Polygon is now an incorporated library. Besides that, Boost.Polygon is an example here, many polygon libraries will have polygon types having no ranges internally. For example, shapelib's geometry is just a set of pointers... Boost.Geometry should handle any geometry. So discarded
Let Boost.Geometry's polygon concept return ranges always as values. This is possible, but the return value for its own polygon should then not be the ring but a proxy containing the reference of the ring.
Let Boost.Geometry's polygon concept return ranges sometimes as value, sometimes as references

I've chatted with my companion Bruno Lalande about this. He had the opinion that we definitely not should go back to iterators. He suggested the fourth solution. Let me quote Bruno: in Boost.Geometry, generic functions like exterior_ring() or interior_ring() that only forward the job to the actual adapted function, shouldn't make any assumption on their return type. Currently they seem to be generic in this regard since they get the return type by calling a metafunction, but they're not fully generic because they arbitrarily add a "&" to it. Boost.Range returns ranges by reference and by value.

Andrei Alexandrescu returns ranges in his article, on various places. We might change our own polygon type, returning proxies to ranges, and returning them by value. This would be solution 3. However, as long as this is not necessary, we prefer solution 4, to avoid an unecessary proxy in between, and to enable range-based polygon-types to return references to ranges.

Exterior rings should sometimes be returned as values, and sometimes as references.

So I changed the Boost.Geometry library into range-return-type agnostic behaviour.

So let me summarize: Ranges can be returned as values or as references

More metaprogramming

To implement solution 4 mentioned above, Boost.Geometry's polygon needed two characters extra: two ampersands. So its traits function ring_type now reads:
template
<
    typename Point, bool ClockWise, bool Closed,
    template<typename, typename> class PointList,
    template<typename, typename> class RingList,
    template<typename> class PointAlloc, template<typename> class RingAlloc
>
struct ring_type
<
    model::polygon
        <
            Point, ClockWise, Closed,
            PointList, RingList, PointAlloc, RingAlloc
        >
>
{
    typedef typename model::polygon
        <
            Point, ClockWise, Closed,
            PointList, RingList,
            PointAlloc, RingAlloc
        >::ring_type& type;
};

where the ampersand is added... Meaning, the ring-returning-type is now a reference!

Of course some more changes were necessary, but not more. The free function exterior_ring now uses a meta-function to define its return type. That meta-function uses the traits-function ring_type, so returns sometimes a reference, sometimes a value. For const and non-const, another metafunction is created:
template <typename T>
struct ensure_const_reference
{
    typedef typename mpl::if_
        <
            typename boost::is_reference<T>::type,
            typename boost::add_reference
                <
                    typename boost::add_const
                        <
                            typename boost::remove_reference<T>::type
                        >::type
                >::type,
            T
        >::type type;
};

I like metaprogramming. If the type was a reference, it removes it, adds a const, and adds a reference again. If it was not a reference, it keeps it unchanged. This is type calculation, metaprogramming.

Iterating through const/non const collections of values/references

At every place where Boost.Range iterators were used to iterate through interior rings, a change was necessary. This was possible creating another metafunction, but Bruno suggested to use BOOST_AUTO here. This avoids new metafunctions, and makes all iterations more readable. This works perfectly for both const and non const iterators, so in fact it is all much more simple than before. Loops are now written like this:

for (BOOST_AUTO(it, boost::begin(interior_rings(polygon)));
it != boost::end(interior_rings(polygon));
++it)

This BOOST_AUTO word, which anticipates the C++0x auto keyword, saves defining and declaring metafunctions all the time. It is a macro (yack) but very valuable here, also because if we change things another time, this can stay as it is now.

One thing is noteworthy, there are two calls to interior_rings returning (in the return-by-value case) two copies. Because these copies are suspected to be proxies, containing references to the original (single) object, the iterator-behaviour will be all right, they are both pointing to the same original range.

However, better safe then sorry, better defensive, so we should write here:

BOOST_AUTO(irings, interior_rings(polygon));
for (BOOST_AUTO(it, boost::begin(irings)); it != boost::end(irings); ++it)

which might be also more efficient in some cases.

Conclusions

Ranges are great.
Ranges are compatible with other objects, while iterator-pairs are not
Ranges can be returned as values or as references
Boost.Geometry's polygon concept uses ranges, allowing adaption of other polygon types using implementation of range proxies

Free Shapefiles of Countries of the World (2)

2010-12-08T23:27:00.004+01:00

Free Shapefiles of Countries of the World (2)

Note: read also part 3 and part 4 of this blog.

After input on my recent blog (thanks!) and after that I find two more, another post about this.

To judge the quality of the dataset, I use my own country (the Netherlands) as the primary measure. I've a mental map about my country, and know how it should look like. In many of the found datasets it looks completely destroyed. Therefore, personally I would discard them at once. Other datasets quite precisely depict my country. Maybe it is not only because I live here, but also about the topography: it has a large inner lake, it has some islands, and it has a distinctive shape.

So here the previous datasets, and four more, and some variations, come again...

	(1) from aprs file world.shp from 2002 244 countries Missing a polder existing since 1968
	(2) from brothersoft, serving blue marble data file World_countries_shp.shp from 1996 239 countries All shape is destroyed
	(3) from Geodan file WorldCountries.shp from 2003 253 countries The Netherlands are depicted correctly
	(4) from mapping hacks file world_borders.shp from 2004 251 countries, 3784 polygons Showing a polder which is never realized
	(5) from thematic mapping file TM_WORLD_BORDERS-0.3.shp from 2008 246 countries Showing a polder which is never realized
	(5b) also from thematic mapping but the simplified version file TM_WORLD_BORDERS_SIMPL-0.3.shp also from 2008 also 246 countries Badly simplified
	(6) from openmap here, called cntry02 file cntry02.shp from 2002 251 countries Showing a polder which is never realized, missing three islands
	(7) from this site ESRI data file cntry08.shp from 2008 249 countries Missing three (inhabited) islands in the North. For the rest it looks like a good map on this scale
	(8) from vdstech via this blog file world.shp from 2001 251 countries Badly simplified
	(9) from natural earth file 10m-admin-0-countries.shp from 2009 251 countries Missing the inner lake the IJsselmeer
	(9b) from natural earth file 50m_admin_0_countries.shp from 2010 240 countries Missing the inner lake the IJsselmeer
	(9c) from natural earth file 110m_admin_0_countries.shp from 2010 177 countries Too much simplified
	(10). I wanted to do VMAP0 / VMAP1 here. However, my little research concerns a (one) shapefile of countries of the world. Vmap0 is great but to let users figure out all the detailed files and merge them for four parts of the world is not the purpose I started the research with... so no picture here. Thanks for the tip anyway.
	(11) from huebler file world_adm0.shp from 2002 209 countries Missing the inner lake the IJsselmeer

Based on this series, (3) would probably be the best choice, though a bit heavily detailed. Surprisingly, it is created in the Netherlands... So we cannot judge its quality from this series. Furthermore, (7) would be a reasonable choice, but we saw in the previous blog that (7) is wrong with respect to the Caspian sea. The new dataset (9b) is looking quite good to me for a map on world scale, but alas it misses an important lake. We might see more in another blog, for example showing the countries around the Caspian sea, or the countries in the former republic of Yugoslavia.

Conclusion

I think a good and free shapefile of the world is still welcome.

Free Shapefile of Countries of the World

2010-12-02T23:32:00.001+01:00

Free Shapefile of Countries of the World

Note: read also part 2 and part 3 of this blog.

I'm always surprised that it is hard to find a good shapefile with world countries. I need country vector data for Boost.Geometry sample data. There is some there but I want to have another set. I want to have it as WKT (Well-Known Text) and I can use a shapefile as input.

Asking Google for "free world map shapefile" you get:

aprs (1), with also a modified version
blue marble requires registration, but (probably) the same set can be downloaded from the next entry
brothersoft, serving blue marble data (2)
mapcruzin (leading to no countries but several interesting other shapefiles)

Via my company, Geodan, you can also download a worldmap, here (requiring mailing your contact information) (3)

Via free gis data, here, there is a hit on ESRI, leading to this site, leading to annual subscriptions etc.. I skipped this one.

Via similar terms ("countries") we also find:

mapping hacks, serving three world files (4)
actually originating from this site, which has a more actual version (5)
a file from openmap here, called cntry02 (6)
and finally (via "cntry08") we find ESRI data here (7)

I've downloaded these shapefiles. It gives me (file dates):

is from 2002, modified in 2009 for Antarctica
is from 1996
is from 2003
is from 2004
is from 2008
is from 2002
is from 2008

The Netherlands
Let's first show my country, the Netherlands:

All borders through each other look quite messy... But let's concentrate on the data.
Southern Flevoland, a polder, already existing since 1968, is still not there on maps from 1 (aprsworld, green) and 2 (blue marble, so blue). Of course it is present on the map of my company 3 (geodan, red), because Geodan of course takes care for its own country. Southern Flevoland is also present on all newer maps (5, pink, 6 orange, 7 gray). But on 5 and 6 it has a planned but not realized polder (Markerwaard) included. So best for the Netherlands are 3 and 7.

Of the five inhabited islands along the Waddenzee, 7 depicts only two. So here 3 is the best choice (of course, it is our country, and note that for the rest this is not a commercial talk)

For the whole world, 6 and 7 are quite similar in nearly all aspects

Uruguay
Let's now look at a country not part of Europe or US. I select Uruguay:

The blue vectors are quite rough and shifted to the west and north. Green lines are looking nice for a map on this scale, but deviate a bit from all other maps.

Let's now take a look in Google Earth. For this, we need to create a KML file... hmm, shp2kml cannot be downloaded (blank website saying Missing ID). OK, we use PostGIS then, it can create KML and WKT (we need it below).

shp2pgsql -s4326 world.shp world1 > world1.sql
etc

So I created a database called blog and executed these SQL files.
Five minutes later (PostGIS is great) I do:
psql "-F " -A -n -t -q -Upostgres -dblog "-cselect '<Placemark>',ST_AsKml(the_geom),'</Placemark>' from world2 where name='Uruguay'" -otworld2.kml

Doing this for all tables and reworking the KML's a bit (adding headers etc) gives me:

Considering Google Maps as our Single Point of Truth, we can discard the blue lines from Blue Marble, and might think that the green lines are also wrong. The rest is following the borders more or less but we need to zoom in to judge it better.

This is the new detail on a border:

The Yellow line is Google's line. The red line is following it, probably a bit too much. For a global scale, the other ones (4,5,6,7) are doing well here.

Caspian Sea
Looking at global level, we see immediately one thing, on which we zoom in here: it is surprising that the Caspian sea is not included in map 6 and 7. So the country borders are just across the sea... Not so good.

Not so good. So let's discard 6 and 7. We also see that the red lines (3) are too detailed, also here. Discard 3. Blue (2) and green (1) was already discarded before. So we keep 4 and 5, which are roughly the same. So I decide to keep the more actualy one: 5.

Decision
So we conclude, based on a few samples (and some more but not described here), that the file TM_WORLD_BORDERS-0.3.zip from Thematic Mapping is the best shapefile to use for worldscale countries (even though it contains the Markerwaard).

Country number and validity
We did not consider attributes, and our research was still rough.

An important aspect we want to consider is geometric validity. And, even more important, the number of countries, it can differ over years, but the differences are sometimes still surprising.

Dataset	Invalid (SQL Server)	Invalid (PostGIS)	#Countries	Year
1	1	1	244	2002 / 2009
2	4	3	239	1996
3	34	5	253	2003
4	5	6 polygons	251 (3784 polygons (this dataset contains polygonsp; not multi-polygons)	2004
5	5	4	246	2008
6	3	0	251	2002
7	2	7	249	2008

About the countries themselves, I just glanced through the differences using a query, comparing ISO3 country codes, getting this table:

We see some inconsistent and some missing ISO codes, and some countries not in the one but in the other, and vice versa.

Conclusion

I think a good and free shapefile of the world is still welcome.

Delphi TDateTimePicker with WeekNumber

2010-11-16T13:56:00.001+01:00

This is not about C++ or C# or Geometry but about Delphi.

I used to have a TDateTimePicker in my app showing the week number. This trick is described on several places, e.g. here and (in Dutch) here

However, since Windows 7 (or probably since Windows Vista) this behaves differently, resulting in an ugly calendar where the right side is cut off :-(

The problem was that the actual control (SysMonthCal32) is, since Windows 7, put inside another control (DropDown), which has the same size as the month calendar. So you can enlarge the month calendar, but it only results in things not being displayed. So the parent control (DropDown) should be resized here...

Getting the parent is easy with Windows (GetParent), but it should still work using XP of course. So we check if the parent is really a "DropDown". There might be other ways to do it, but I'm using GetClassName now and that works.

So the adapted solution is here:


procedure TMyForm.BeginDateTimePickerDropDown(Sender: TObject);
const
  MCM_GETMAXTODAYWIDTH  = MCM_FIRST + 21;
var
  Style: LongInt;
  hDTP: THandle;
  r: TRect;
  intTodayWidth: Integer;

  cname: array[0..256] of Char;
begin
  //to get a handle of calendar
  hDTP := DateTime_GetMonthCal((sender as TDateTimePicker).Handle);

  //change a style
  Style := GetWindowLong(hDTP, GWL_STYLE);
  SetWindowLong(hDTP,  GWL_STYLE, Style or MCS_WEEKNUMBERS);

  //now we must change the width for calendar because week numbers shifted all strings
  //1. to get the required rect
  r := Rect(0, 0, 0, 0);
  SendMessage(hDTP, MCM_GETMINREQRECT, 0, Longint(@r));

  //2. to get the maximum width of the "today" string
  intTodayWidth := SendMessage(hDTP, MCM_GETMAXTODAYWIDTH, 0, 0);

  //3. adjust rect width to fit the "today" string
  if intTodayWidth > r.Right then
    r.Right := intTodayWidth;

  // For Win7, the window (class=SysMonthCal32) is automatically inside
  // a parent-window (class=DropDown). Check this. If so, take the parent window.
  // If not so (class=TMainForm), take this window (for XP and lower)
  GetClassName(GetParent(hDtp), cname, sizeof(cname));
  if AnsiSameText(cname, 'DropDown') then
  begin
    hDTP := GetParent(hDTP);

    // To get it perfect (on my machines) is adding this code:
    inc(r.Right, 6);
    inc(r.Bottom, 6);
  end;

  //4. to set new the height and width
  MoveWindow(hDTP, r.Left, r.Top, r.Right-r.Left, r.Bottom-r.Top, true);
end;

Union, Tangencies and Colors

2010-10-26T23:25:00.000+02:00

Union, Tangencies and Colors

I like the picture below very much.

What you see here are, more or less, pink, green and blue. The colors are chosen as such that they always form distinct colors, regardless of how you combine them. There are of course other combinations with these properties, but... this is at least one of them.

When the pink layer is slightly moved to the right, you see better how the image is constructed:

The pink (transparent) layer is the union of the two underlying (multi-)polygons, the green one and the blue one. It has one hole, one kind of indentation and for the rest it is square. The underlying multi-polygons are better seen if the pink one is removed:

So this is clear: there are two multi-polygons here, complete with holes and self-tangencies. These multi-polygons have overlap in some places, no overlap in other places. They nearly fill the area, there are only three "boxes" free. The blue and green colors nicely mix into blue-green. The colors and transparencies are selected for that. SVG definitions are:

Green: style="fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153,204,0);stroke-width:3"

Blue: style="fill-opacity:0.3;fill:rgb(51,51,153);stroke:rgb(51,51,153);stroke-width:3"

I like this combination.

And with the pink style on top of it...

Pink: style="fill-opacity:0.2;stroke-opacity:0.4;fill:rgb(255,0,0);stroke:rgb(255,0,255);stroke-width:8"

...all colors change, however, they are still clearly distinct. Yellow-orange like for green+pink; lila for blue+pink; and brown for the blue+green+pink. So we have seven unique distinct colors now, which happens to be the ideal number of colors for a map.

By the way, this image is created using Boost.Geometry. Boost.Geometry can easily create SVG images. And Boost.Geometry can overlay multi-polygons (already possible for a long time) having all kinds of self-tangencies (new!). A self-tangency in a multi-polygon is a place where two polygon-members touch each other (orange dots in the map below), or where a polygon touches itself (the red dot in the map).

Actually... this could also be a polygon, having five interior rings (holes) which might be mutually tangent and tangent to they exterior ring.

OK, one more image...

Tag Dispatching and Inheritance

2010-10-19T19:24:00.000+02:00

Tag Dispatching and Inheritance

In the previous blog we described Tag Dispatching by Type.

Though very convenient, when many classes (tagged by as many tags) share the same implementation, we get as many dispatch structs with the same implementation (or, better, forwarding to the same struct). This can be avoided by using similarities between classes, or between tags. In Boost.Geometry, for example, some geometries are single while others are multi. That multi-ness can be recognized by the metafunction is_multi which results in true for a multi-polygon and in false for a single polygon. We can include an IsMulti boolean template parameter in the dispatch structure, so that it is selected not only on tag dispatching, but also in is_multi.

This similarity can also be implemented using a tag hierarchy and inheritance.

Let's go back to the fruit example. All citrus fruit species have pulp vesicles. There are many citrus fruits species. And, of course, we don't want to repeat the dispatch structure has_vesicles so many times. So we try to define our tag hierarchy:

struct citrus_tag {};
struct pome_tag {};

struct apple_tag : pome_tag {};
struct pear_tag : pome_tag {};

struct banana_tag {};

struct orange_tag : citrus_tag {};
struct lemon_tag : citrus_tag {};
struct lime_tag : citrus_tag {};

And then we couple the tags, like we did in previous blog:

template <typename T> struct tag {};
template <> struct tag<apple> { typedef apple_tag type; };
template <> struct tag<pear> { typedef pear_tag type; };
template <> struct tag<orange> { typedef orange_tag type; };
template <> struct tag<lime> { typedef lime_tag type; };
template <> struct tag<banana> { typedef banana_tag type; };

The new thing, besides the inheritance in the tags above, is that we now define a tag_cast metafunction:

template
<
    typename Tag, typename BaseTag,
    typename BT2 = void, typename BT3 = void, typename BT4 = void,
    typename BT5 = void, typename BT6 = void, typename BT7 = void
>
struct tag_cast
{
    typedef typename boost::mpl::if_
        <
          typename boost::is_base_of<BaseTag, Tag>::type,
          BaseTag,
          // Try next one in line:
          typename tag_cast<Tag, BT2, BT3, BT4, BT5, BT6, BT7, void>::type
        >::type type;
};

template <typename Tag>
struct tag_cast<Tag, void, void, void, void, void, void, void>
{
    // If not found, take specified tag, so do not cast
    typedef Tag type;
};

We can now check if a tag is a citrus tag, by calling defining a new tag as tag_cast<tag, citrus_tag>::type. This new tag is either a citrus tag, or it is the tag it already was.

And we can dispatch by this tag as well. So it looks like:

namespace dispatch
{
    template <typename T> struct has_vesicles : boost::false_type {};
    template <> struct has_vesicles<citrus_tag> : boost::true_type {};
}

template <typename Fruit>
std::string has_vesicles(Fruit const& fruit)
{
    // (Potentially) go up in hierachy: take tag corresponding to Fruit,
    // downcast to citrus_tag if possible
    typedef typename tag_cast<typename tag<Fruit>::type, citrus_tag>::type tag;

    return std::string("has vesicles: ")
        + (dispatch::has_vesicles<tag>::value ? "true" : "false");
}

This tag_cast, as proposed above, can not only cast to one base-class, but also to a range of base-classes. The definition

 typedef typename tag_cast<typename tag<Fruit>::type, citrus_tag, pome_tag>::type tag;

will result in a citrus_tag if it is a citrus, a pome_tag if it is a pome, and otherwise it results in the input tag. Besides this, it also walks through a tag hierarchy, so if citrus_tag and pome_tag were both derived from, e.g., a rosid_tag, and rosid_tag was specified in this call, it would result in a rosid_tag.

The sample program here shows two base classes, not listed here.

Our main program now displays correctly if the specific fruits have vesicles or not

int main()
{
    using namespace fruit;

    apple a("my apple");
    pear p("my pear");
    orange o("my orange");

    std::cout << has_vesicles(a) << std::endl;
    std::cout << has_vesicles(p) << std::endl;
    std::cout << has_vesicles(o) << std::endl;

    return 0;
}

Back to Boost.Geometry. Until now it is not done, but we might depricate the is_multi meta-function and replace it by this:

struct multi_tag;
struct multi_point_tag : multi_tag {};
struct multi_linestring_tag : multi_tag {};
struct multi_polygon_tag : multi_tag {};

And instead of calling is_multi, we call tag_cast<tag, multi_tag>. We then know if it is a multi and can dispatch on it. The same for linear. This will probably be very convenient, it has to be worked out a little more.

Casting and (multiple) inheritance
The tag_cast as implemented above casts to specified tags. We can not just call the parent-tag. One of the reasons for this is multiple inheritance.
Tag inheritance can be multiple. This is convenient because, of course, there can be different tag classifications for the same set of tags. For example: if a geometry contains segments, or is linear, or is areal. So we might get:

struct multi_tag;
struct multi_point_tag : multi_tag {};
struct multi_linestring_tag : multi_tag, linear_tag {};
struct multi_polygon_tag : multi_tag, areal_tag {};

When calling tag_cast, the results depends on the order of the specified base tags. And so does the dispatching.

Tag Dispatching by Type

2010-10-19T17:32:00.000+02:00

Tag Dispatching by Type

Tag Dispatching (TD) is a Generic Programming (GP) technique for C++. It is well explained in this blog. Google Hit Number one (GH#1) is this site, where the example of the STL function advance is given. Tag Dispatching is not described a lot; GH#7 is from ggl (= Boost.Geometry), as is GH#8..

Nearly all TD descriptions I know of use the example of STL-advance. But in Boost.Geometry TD is done in another way: not by instance, but by type. The difference is explained in this blog.

Suppose you make a generic program doing something with fruit. An apple is normally being eaten differently than a banana or an orange. So if we program a generic function eat, giving it any type of fruit, it should redirect this to an implementation specific to that kind of fruit, so for apple differently than for a banana.

We first implement the tags. Tags are completely empty structures. But (of course) not completely useless.

struct apple_tag {};
struct banana_tag {};
struct orange_tag {};

We then create some sample implementations. Most properties are not used, it is just to show there are (sometimes completely) different classes which can be handled genericly.

struct apple
{
    double radius;

    std::string name;
    apple(std::string const& n) : name(n) {}
};

struct banana
{
    double length;

    std::string name;
    banana(std::string const& n) : name(n) {}
};

All this is quite simple C++ code, though you have to understand templates (generics) and specialization.

OK, now we implement Tag Dispatching by Instance:

namespace dispatch
{
    void eat(apple const& a, apple_tag)
    {
        std::cout << "bite" << std::endl;
    }

    void eat(banana const& b, banana_tag)
    {
        std::cout << "peel" << std::endl;
    }
}

template <typename T>
void eat(T const& fruit)
{
    typename tag<T>::type the_tag;
    dispatch::eat(fruit, the_tag);
}

The function eat at the bottom first declares an instance of the tag. If an apple is entered in this function, it will be an instance of the apple_tag. If a banana is entered, it is a banana_tag. Then it forwards its call to the eat function in the namespace dispatch. There are two versions (overloads) there, one specified with an apple_tag, one with a banana_tag. The compiler selects, based on the tag the right function. Quite easy, and quite powerful, this tag dispatching system.

As said, within Boost.Geometry it is not applied like this. One of the reasons is that the instance of the tag is not necessary at all. If the dispatch::eat function would be implemented within a struct, as a static method, the tag dispatching can be done by type and not by instance. This is still tag dispatching!

So the piece above can be replaced by this Tag Dispatching by Type:

namespace dispatch
{
    template <typename Tag> struct eat {};

    template <> struct eat<apple_tag>
    {
        static void apply(apple const& a)
        {
          std::cout << "bite" << std::endl;
        }
    };

    template <> struct eat<banana_tag>
    {
        static void apply(banana const& b)
        {
          std::cout << "peel" << std::endl;
        }
    };
}

template <typename T>
void eat(T const& fruit)
{
    dispatch::eat<typename tag<T>::type>::apply(fruit);
}

So no instance at all, type only. Boost.Geometry has this structure everywhere, and all static dispatching methods are called apply.

For completeness we list the main program as well here:

int main()
{
    apple a("my apple");
    banana b("my banana");
    eat(a);
    eat(b);

    return 0;
}

It is the same for the instance-version and the type-version.

More advantages of the Tag Dispatching by Type approach are that arguments can be reversed, so you can call distance(point, polygon) and distance(polygon, point) but there is only one implementation with point_tag, polygon_tag necessary. With instances (and function overloads) this is not possible, or at least not that easy. Furthermore, tag dispatching by type can be used to define types, or define constants, as well. Suppose you need to define a constant value dependant on a type. We approach geometry now and want to define if a specific fruit-type is spherical or not. So we define this structure, specialized by tag:

    template <typename Tag> struct spherical {};

    template <> struct spherical<apple_tag>
    {
        static const bool value = true;
    };

    template <> struct spherical<banana_tag>
    {
        static const bool value = false;
    };

We can create a generic structure, which calls this dispatch, just like the free function eat, but now with a struct (or meta-function):

template <typename T>
struct spherical
{
    static const bool value = dispatch::spherical<typename tag<T>::type>::value;
};

Next blog will handle Tag Dispatching and Inheritance.

C# Specializations by Traits IV

2010-07-26T22:00:00.000+02:00

C# Specializations by Traits IV

Geometry

In last three blogs (I'm productive these days in blogs, instead of klocs) I explained how to to specializations in C#. This blog only explains how it begon, give different possibilities to design things, and wraps it up.

Towards specializations

Developing Boost.Geometry, in C++, and doing some more C# this year as well, I wondered how a generic geometry library in C# would look like. Ideally it would be as Boost.Geometry, of course... But Boost.Geometry is based on traits and C++ specializations, partial specializations, tag dispatching. All modern C++ generic programming stuff. Not easily doable in C#.

Calculating the distance between two arbitrary point types would already be something. I tried many methods, and discussed them with my collegue Paul den Dulk, and he helped me, suggesting more methods and some essential pieces. Putting it together, adding the dictionary as the runtime register instead of the compiletime C++ specializations, the point distance function is looking exactly the same as Boost.Geometry now!

Some notes

There are still many differences from this system with C++ specializations, of course.

in C++ it is compile-time and checked at compile-time
therefore in C++ it is very fast
in C# the dictionary is runtime and checked at runtime
because the dictionary lookup and necesary cast it is slower in C#
C++ has partial specializations, don't know (yet) how to do that in C#
C++ has type calculations, the whole template metaprogramming is based on that.
in C# it is only the methods
because in C# it is runtime, it gives some advantages as well: you can register different traits classes when needed
and so, e.g. the ORM can handle the same model in different ways

Timings

As said the dictionary consultation, together with the cast, take some time. Using the geometry case distance, I timed it and it is certainly measurable.

Ten million distance calculations on my computer:

Using interfaces: 0.10 seconds
Using delegates: 0.15 seconds
Using traits: 0.13 seconds
Using a register: 0.42 seconds

So, in fact 4 times slower. It is much in this case, but this is compared to a small calculation (Pythagoras). If applied using GUI, database, disk io, or occasionally (and not 10 million times), the lookup will probably fall away.

Complete program

The program to research, with 10 possible implementations, with interfaces, wrappers, delegates, traits, registers and lambdas, is shown below.

using System;
using System.Diagnostics;

using System.Collections.Generic;


namespace GenericGeometry
{
    interface IPoint
    {
        double x();
        double y();
        void increase();
    }

    delegate double get_x<Point>(Point p);
    delegate double get_y<Point>(Point p);

    public interface ITraitsBase {}

    public interface ITraits<Point> : ITraitsBase
    {
        double get_x(Point p);
        double get_y(Point p);
    }

    class TraitsRegister
    {
        public static IDictionary<Type, ITraitsBase> register;

        public static void Register<PointType>(ITraits<PointType> t)
        {
            if (register == null)
            {
                register = new Dictionary<Type, ITraitsBase>();
            }
            register[typeof(PointType)] = t;
        }
    }

    class LambdaRegister
    {
        public static IDictionary<Type, Tuple<object, object>> register;

        public static void Register<Point>(get_x<Point> gx, get_y<Point> gy)
        {
            if (register == null)
            {
                register = new Dictionary<Type, Tuple<object, object>>();
            }
            register[typeof(Point)] = new Tuple<object, object>(gx, gy);
        }
    }


    // The generic library, with functions as "Distance", working on IPoint or on any point
    class Algorithm
    {
        // (1) Using an interface to access coordinates
        public static double Distance(IPoint p1, IPoint p2)
        {
            double dx = p1.x() - p2.x();
            double dy = p1.y() - p2.y();
            return Math.Sqrt(dx * dx + dy * dy);
        }

        // (7) use delegates to access coordinates of not interfaced point type
        public static double Distance<Point>(Point p1, Point p2, 
            get_x<Point> gx, 
            get_y<Point> gy)
        {
            double dx = gx(p1) - gx(p2);
            double dy = gy(p1) - gy(p2);
            return Math.Sqrt(dx * dx + dy * dy);
        }

        // (8) use traits/accessor to access coordinates of not interfaced point type
        public static double Distance<Point>(Point p1, Point p2, ITraits<Point> traits)
        {
            double dx = traits.get_x(p1) - traits.get_x(p2);
            double dy = traits.get_y(p1) - traits.get_y(p2);
            return Math.Sqrt(dx * dx + dy * dy);
        }

        // (9) use traits/REGISTER accessor to access coordinates of not interfaced point type
        public static double DistanceTR<Point>(Point p1, Point p2)
        {
            ITraits<Point> traits = TraitsRegister.register[typeof(Point)] as ITraits<Point>;
            double dx = traits.get_x(p1) - traits.get_x(p2);
            double dy = traits.get_y(p1) - traits.get_y(p2);
            return Math.Sqrt(dx * dx + dy * dy);
        }

        // (10) use (tupled) lambda/REGISTER accessor to access coordinates of not interfaced point type
        public static double DistanceLR<Point>(Point p1, Point p2)
        {
            var pair = LambdaRegister.register[typeof(Point)] as Tuple<object, object>;
            get_x<Point> gx = pair.Item1 as get_x<Point>;
            get_y<Point> gy = pair.Item2 as get_y<Point>;
            double dx = gx(p1) - gx(p2);
            double dy = gy(p1) - gy(p2);
            return Math.Sqrt(dx * dx + dy * dy);
        }
    }


    // (1) implement IPoint
    class MyPoint1 : IPoint
    {
        private double m_x, m_y;
        public MyPoint1(double x, double y)
        {
            m_x = x;
            m_y = y;
        }
        public double x() { return m_x; }
        public double y() { return m_y; }
    }


    // Non-interfaced class...
    class MyPoint2
    {
        public double x { get; set; }
        public double y { get; set; }
    }

    // (2) Create a wrapper (implementing IPoint)
    class PointWrapper : IPoint
    {
        MyPoint2 m_point;
        public PointWrapper(MyPoint2 p)
        {
            m_point = p;
        }

        public double x() { return m_point.x; }
        public double y() { return m_point.y; }
    }

    // (3) Create a generic wrapper using delegates, implementing IPoint
    class DelPointWrapper<Point> : IPoint
    {
        Point m_point;
        get_x<Point> m_gx;
        get_y<Point> m_gy;

        public DelPointWrapper(Point p, get_x<Point> gx, get_y<Point> gy)
        {
            m_point = p;
            m_gx = gx;
            m_gy = gy;
        }

        public double x() { return m_gx(m_point); }
        public double y() { return m_gy(m_point); }

        public void increase()
        {
            //m_x += 0.000001;
            //m_y += 0.000001;
        }

    }

    // (7) define delegates in a static class
    static class Accessor
    {
        public static get_x<MyPoint2> getX = p => p.x;
        public static get_y<MyPoint2> getY = p => p.y;
    }


    // (8,9) use a sort of Traits class
    class Traits : ITraits<MyPoint2>
    {
        public double get_x(MyPoint2 p)
        {
            return p.x;
        }
        public double get_y(MyPoint2 p)
        {
            return p.y;
        }
    }

    class Program
    {
        const int n = 10000000;

        static void Timed<Point>(String header, Point a, Point b) where Point : IPoint
        {
            double d = 0.0;
            Stopwatch st = new Stopwatch();
            st.Start();
            for (int i = 0; i < n; i++)
            {
                d += Algorithm.Distance(a, b);
            }
            st.Stop();
            System.Console.WriteLine(header + " " + d + " " + st.Elapsed.ToString());
            System.Threading.Thread.Sleep(500);

        }

        // Creating the delegates
        static double point_get_x(MyPoint2 p)
        {
            return p.x;
        }
        static double point_get_y(MyPoint2 p)
        {
            return p.y;
        }

        static void Main(string[] args)
        {

            int option = args.Length > 0 ? Int32.Parse(args[0]) : 10;
            System.Console.Write("{0} ", option);

            // Point without interface
            MyPoint2 p1 = new MyPoint2();
            p1.x = 3;
            p1.y = 5;
            MyPoint2 p2 = new MyPoint2();
            p2.x = 6;
            p2.y = 9;

            switch (option)
            {
                case 1:
                    {
                      // Option 1, point with interface
                      MyPoint1 ip1 = new MyPoint1(3, 5);
                      MyPoint1 ip2 = new MyPoint1(6, 9);
                      Timed("Using interface", ip1, ip2);
                    } break;
                case 2:
                    {
                      // Option 2a:
                      PointWrapper pw1 = new PointWrapper(p1);
                      PointWrapper pw2 = new PointWrapper(p2);
                      Timed("Using wrapper", pw1, pw2);
                    } break;
                case 3:
                    {
                      // Option 2b:
                      DelPointWrapper<MyPoint2> dpw1 = new DelPointWrapper<MyPoint2>(p1, point_get_x, point_get_y);
                      DelPointWrapper<MyPoint2> dpw2 = new DelPointWrapper<MyPoint2>(p2, point_get_x, point_get_y);
                      Timed("Using delegate-wrapper", dpw1, dpw2);
                    } break;
                case 4:
                    {
                      // Option 2c: with inline delegates
                      DelPointWrapper<MyPoint2> idpw1 = new DelPointWrapper<MyPoint2>(p1
                      , delegate(MyPoint2 p) { return p.x; }
                      , delegate(MyPoint2 p) { return p.y; }
                      );
                      DelPointWrapper<MyPoint2> idpw2 = new DelPointWrapper<MyPoint2>(p2
                      , delegate(MyPoint2 p) { return p.x; }
                      , delegate(MyPoint2 p) { return p.y; }
                      );
                      Timed("Using inline delegate-wrapper", idpw1, idpw2);
                    } break;
                case 5:
                    {
                      DelPointWrapper<MyPoint2> ldpw1 = new DelPointWrapper<MyPoint2>(p1, p => p.x, p => p.y);
                      DelPointWrapper<MyPoint2> ldpw2 = new DelPointWrapper<MyPoint2>(p2, p => p.x, p => p.y);
                      Timed("Using lambda-wrapper", ldpw1, ldpw2);
                    } break;
                case 6:
                    {
                      double d = 0.0;
                      Stopwatch st = new Stopwatch();
                      st.Start();
                      for (int i = 0; i < n; i++)
                      {
                      d += Algorithm.Distance(p1, p2, p => p.x, p => p.y);
                      }
                      st.Stop();
                      System.Console.WriteLine("Using delegate-algorithm+lambda " + d + " " + st.Elapsed.ToString());
                    } break;
                case 7:
                    {
                      double d = 0.0;
                      Stopwatch st = new Stopwatch();
                      st.Start();
                      for (int i = 0; i < n; i++)
                      {
                      d += Algorithm.Distance(p1, p2, Accessor.getX, Accessor.getY);
                      }
                      st.Stop();
                      System.Console.WriteLine("Using accessor met static delegates " + d + " " + st.Elapsed.ToString());
                    } break;
                case 8:
                    {
                      Traits traits = new Traits();

                      double d = 0.0;
                      Stopwatch st = new Stopwatch();
                      st.Start();
                      for (int i = 0; i < n; i++)
                      {
                      d += Algorithm.Distance(p1, p2, traits);
                      }
                      st.Stop();
                      System.Console.WriteLine("Using traits " + d + " " + st.Elapsed.ToString());
                    } break;
                case 9:
                    {
                      Traits traits = new Traits();
                      TraitsRegister.Register<MyPoint2>(traits);

                      double d = 0.0;
                      Stopwatch st = new Stopwatch();
                      st.Start();
                      for (int i = 0; i < n; i++)
                      {
                      d += Algorithm.DistanceTR(p1, p2);
                      }
                      st.Stop();
                      System.Console.WriteLine("Using traits/register " + d + " " + st.Elapsed.ToString());
                    } break;
                case 10:
                    {
                      LambdaRegister.Register<MyPoint2>(p => p.x, p => p.y);

                      double d = 0.0;
                      Stopwatch st = new Stopwatch();
                      st.Start();
                      for (int i = 0; i < n; i++)
                      {
                      d += Algorithm.DistanceLR(p1, p2);
                      }
                      st.Stop();
                      System.Console.WriteLine("Using lamba/register " + d + " " + st.Elapsed.ToString());
                    } break;
            }

            System.Console.WriteLine("");
        }
    }
}

And here are the complete timings (for each the smallest of 5 runs):

1 Using interface 50000000 00:00:00.1000114
2 Using wrapper 50000000 00:00:00.1055539
3 Using delegate-wrapper 50000000 00:00:00.1812642
4 Using inline delegate-wrapper 50000000 00:00:00.1843941
5 Using lambda-wrapper 50000000 00:00:00.1810666
6 Using delegate-algorithm+lambda 50000000 00:00:00.1509385
7 Using accessor met static delegates 50000000 00:00:00.1274787
8 Using traits 50000000 00:00:00.1327729
9 Using traits/register 50000000 00:00:00.4192440
10 Using lamba/register 50000000 00:00:00.4194542

Conclusions

Are already done. Specializations in C# are possible. Use them, it makes C# and generics more fun.

C# Specializations by Traits III

2010-07-25T21:28:00.001+02:00

C# Specializations by Traits III

Stormy: Specializations by Traits - Object Relational Model

I like SOCI very much. It is a C++ ORM. Especially the way to go from SQL to model, the object-relational mapping, is very convenient. It is based on specializations.

In the earlier Blogs of this series I have shown that it is possible to do real specializations in C#. This is possible using the combination of traits, a dictionary, and (if necessary) extension methods. This is all explained in previous blogs and I will not explain it again. I will only show how this system is really convenient in Object Relational Mapping (ORM). A hundred-lines prototype is developed. Called Stormy (storm was already occupied).

Use case

In ORM the database (relational) is mapped to objects (models, C# classes). So by nature there are many classes all handled the same way (from SQL to the class instance).

Instead of using the ADO.NET SqlDataReader (delivering just abstract fields and not a model), we want to do things like this:

IEnumerable<Cat> cats = connection.Select<Cat>("select * from cat");

where

Cat is our model, shown below
connection is any connection to a database
the SQL statement is shown, literally. I like that.
so it returns a generic collection of our model

In C# it is not trivial to specialize the Select method per model. But with traits and a register it can be solved.

The model Cat

The model Cat is borrowed from nHibernate (I know nHibernate - it will not map the class below - it wants method to be virtual). Slightly modified (I like this way of specifying attributes):


public class Cat
{
    public Guid Id { get; set; }
    public string Name { get; set; }
    public char Sex { get; set; }
    public float Weight { get; set; }
}

The Object Relational mapping

The ORM is done entirely in code, the way SOCI does it. No XML mappings here. We stay with C#. We show here the implementations Map the Cat from (relational) Database to Object, and vice versa.


public class CatConverter : ITypeConversion<Cat>
{
    public Cat fromSql(SqlDataReader reader)
    {
        Cat cat = new Cat();
        cat.Id = (Guid)reader["CatId"];
        cat.Name = (string)reader["name"];
        cat.Sex = (((string)reader["sex"]) + " ")[0];
        cat.Weight = (float)reader["weight"];
        return cat;
    }

    public void fillParameters(Cat cat, SqlCommand command)
    {
        command.Parameters.Add(new SqlParameter("@name", cat.Name));
        command.Parameters.Add(new SqlParameter("@sex", cat.Sex));
        command.Parameters.Add(new SqlParameter("@weight", cat.Weight));
    }
}

Note that it is mapped using SqlDataReader. It can be made more convenient (eliminating casts) of course, but this is how it is done with SqlDataReader, normally.

That SOCI-like type-converter (I renamed the methods from_base to fromSQL and to_base to fillParameters) is an interface ITypeConversion<T>, and is registered (normally) once per model:

Orm.Register<Cat>(new CatConverter());

Main program

We can implement our scenario, converting from database to Cat and vice versa, using generic programming. Our program:

class Program
{
    static void Main(string[] args)
    {
        Orm.Register<Cat>(new CatConverter());

        Connection connection = new Connection(
            @"
                Data Source=localhost\SQLEXPRESS;
                MultipleActiveResultSets=True; 
                Initial Catalog=stormy;
                Integrated Security=SSPI;
            ");

        var cats = connection.Select<Cat>("select * from cat");
        foreach(Cat cat in cats)
        {
            System.Console.WriteLine("Cat {0} named {1}", cat.Id, cat.Name);
        }

        // Create subselection using linq (not related to SpecializationbyTraits)
        var toms = from cat in cats where cat.Name == "Tom" select cat;
        if (toms.Count() > 0)
        {
            Cat tom = toms.First();
            System.Console.WriteLine("Tom found, ID={0}", tom.Id);
        }
        else
        {
            Cat tom = new Cat();
            tom.Name = "Tom";
            tom.Sex = 'm';
            tom.Weight = 5.2F;
            connection.Execute(CatConverter.InsertSql(), tom);
        }
    }
}

We first register our type conversion specialization. We create a connection to the database (hard coded - alternatively you might get it from e.g. appSettings). We can use any SQL statement to fill any list with any (registered) model. For fun we use LINQ here to do a subselection. We can execute SQL to insert, update or delete things. That InsertSql is specified (for convenience, I prefer it close to its parameters) within the converter as well, it can be anywhere, up to you.

public static string InsertSql()
{
    return "insert into Cat(CatId, Name, Sex, Weight) values(NEWID(), @name, @sex, @weight)";
}

Stormy ORM - wrapping it up

Great isn't it? The Stormy code is about hundred lines now. Of course it is not a full-featured ORM as others are, but in many cases more convenient, more simple to understand, it does not create shadow classes, and it does not force you to learn a new criteria language. If you know SQL, you will like it!

The project including sample can be downloaded here.

We're nearly done with C# Specializations. It is really useful. Can be used everywhere. Makes C# closer to C++.

Next:

geometry...

Stay tuned!