[Boost-commit] svn:boost r55831 - in sandbox/ggl/other/comparisons: cgal common ggl gpc gtl

Subject: [Boost-commit] svn:boost r55831 - in sandbox/ggl/other/comparisons: cgal common ggl gpc gtl
From: barend.gehrels_at_[hidden]
Date: 2009-08-28 04:54:02

Author: barendgehrels
Date: 2009-08-28 04:54:00 EDT (Fri, 28 Aug 2009)
New Revision: 55831
URL: http://svn.boost.org/trac/boost/changeset/55831

Log:
Added geometry library performance comparison CPP files
Added:
   sandbox/ggl/other/comparisons/cgal/cgal_check.cpp   (contents, props changed)
   sandbox/ggl/other/comparisons/common/common.hpp   (contents, props changed)
   sandbox/ggl/other/comparisons/common/read_shapefile.hpp   (contents, props changed)
   sandbox/ggl/other/comparisons/ggl/ggl_check.cpp   (contents, props changed)
   sandbox/ggl/other/comparisons/gpc/gpc_check.cpp   (contents, props changed)
   sandbox/ggl/other/comparisons/gtl/gtl_check.cpp   (contents, props changed)

Added: sandbox/ggl/other/comparisons/cgal/cgal_check.cpp
==============================================================================
--- (empty file)
+++ sandbox/ggl/other/comparisons/cgal/cgal_check.cpp	2009-08-28 04:54:00 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,395 @@
+
+#include "../common/common.hpp"
+#include "../common/read_shapefile.hpp"
+
+#include <iostream>
+
+#include <CGAL/basic.h>
+
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Cartesian.h>
+
+#include <CGAL/Polygon_2.h>
+
+#include <CGAL/ch_graham_andrew.h>
+#include <CGAL/Boolean_set_operations_2.h>
+#include <CGAL/intersections.h>
+#include <CGAL/centroid.h>
+
+
+#ifdef CGAL_USE_GMP
+   #include <CGAL/Gmpq.h>
+#endif
+
+#include <CGAL/MP_Float.h>
+#include <CGAL/Quotient.h>
+
+
+
+// Define the kernel.
+// NOTE: this is also very important for perfomance
+// Besides that, intersections might crash in some situations for some kernels...
+
+//typedef CGAL::Cartesian<float> K; // -> nearly all polygons are invalid for intersection
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K; // fastest, + only K that compiles for MinGW
+//typedef CGAL::Cartesian<double> K; // slower
+//typedef CGAL::Exact_predicates_exact_constructions_kernel K;
+//typedef CGAL::Cartesian<CGAL::Gmpq> K;
+//typedef CGAL::Cartesian<CGAL::Quotient<CGAL::MP_Float> > K; // slowest
+
+
+typedef K::Point_2 POINT;
+typedef CGAL::Polygon_2<K> POLY;
+typedef CGAL::Bbox_2 BOX;
+
+typedef CGAL::Polygon_with_holes_2<K> HOLEY_POLY; // necessary for intersections
+
+
+template <typename T>
+T convert(SHPObject* psShape)
+{
+    const double* x = psShape->padfX;
+    const double* y = psShape->padfY;
+
+    POLY poly;
+
+    // for intersections:
+    // 1: polygon must NOT be closed otherwise: "Cannot contruct a degenerate segment"
+            /*CGAL-intersections-degenerate:
+
+            terminate called after throwing an instance of 'CGAL::Precondition_exception'
+              what():  CGAL ERROR: precondition violation!
+            Expr: ! is_degen
+            File: ../../../../../CGAL-3.4/include/CGAL/Arr_segment_traits_2.h
+            Line: 120
+            Explanation: Cannot contruct a degenerate segment.
+            Aborted*/
+
+    // 2: polygon must be counter-clockwise otherwise: "The polygon has the wrong orientation"
+    // 3: still get an error, now:
+            /*terminate called after throwing an instance of 'CGAL::Precondition_exception'
+              what():  CGAL ERROR: precondition violation!
+            Expr: comp_f(object, parentP->object) != SMALLER
+            File: ../../../../../CGAL-3.4/include/CGAL/Multiset.h
+            Line: 2128
+            Aborted
+            */
+
+
+    // 4: when going to quotient/etc, we come the most far, getting for one polygon:
+            /*
+            CGAL warning: check violation!
+            Expression : is_strictly_simple
+            File       : ../../../../../CGAL-3.4/include/CGAL/Boolean_set_operations_2/Gps_polygon_validation.h
+            Line       : 174
+            Explanation: The polygon is not strictly simple.
+            Refer to the bug-reporting instructions at http://www.cgal.org/bug_report.html
+            terminate called after throwing an instance of 'CGAL::Precondition_exception'
+              what():  CGAL ERROR: precondition violation!
+            Expr: m_traits->is_valid_2_object()(pgn)
+            File: ../../../../../CGAL-3.4/include/CGAL/General_polygon_set_2.h
+            Line: 200
+            */
+
+    for (int v = psShape->nVertices - 1; v > 0; v--)
+    {
+        // CGAL's intersection does not like degenerate segments. So sort them out here.
+        if (! (v > 0
+            && std::abs(x[v - 1] - x[v]) < 1e-10
+            && std::abs(y[v - 1] - y[v]) < 1e-10))
+        {
+            poly.push_back(POINT(x[v], y[v]));
+        }
+    }
+
+    return poly;
+}
+
+
+
+
+int main(int argc, char** argv)
+{
+    compare::version_info("cgal");
+    std::vector<POLY> polygons;
+    std::vector<POLY> ellipses;
+    std::vector<BOX> boxes;
+    //std::vector<BOX> clip_boxes;
+    std::vector<int> ids;
+
+
+
+    read_shapefile(compare::filename(argc, argv), compare::fieldname(), polygons, ids, convert<POLY>);
+
+
+    // Create envelopes
+    for (std::vector<POLY>::const_iterator it = polygons.begin();
+        it != polygons.end();
+        ++it)
+    {
+        BOX box = CGAL::bbox_2(it->vertices_begin(), it->vertices_end(), K());
+        boxes.push_back(box);
+    }
+
+    // Create the ellipses for intersections lateron
+    if (compare::MEASURE_OVERLAY)
+    {
+        for (std::vector<BOX>::const_iterator bit = boxes.begin();
+            bit != boxes.end();
+            ++bit)
+        {
+            double cx = 0.5 * (bit->xmin() + bit->xmax());
+            double cy = 0.5 * (bit->ymin() + bit->ymax());
+
+            double dx = bit->xmax() - bit->xmin();
+            double dy = bit->ymax() - bit->ymin();
+
+            double a1 = compare::OVERLAY_ELLIPSE_FACTOR1 * 0.5 * dx;
+            double b1 = compare::OVERLAY_ELLIPSE_FACTOR1 * 0.5 * dy;
+            double a2 = compare::OVERLAY_ELLIPSE_FACTOR2 * 0.5 * dx;
+            double b2 = compare::OVERLAY_ELLIPSE_FACTOR2 * 0.5 * dy;
+
+            double angle = 0.0;
+
+            POLY ellipse;
+
+            // Create COUNTER-clockwise ellipse (therefore the minus in angle), WITHOUT last point
+            // (therefore the -1)
+            for (int v = 0; v < compare::OVERLAY_ELLIPSE_COUNT - 1; v++, angle -= compare::delta)
+            {
+                if (v % 2 == 0)
+                {
+                    ellipse.push_back(POINT(cx + a1 * sin(angle), cy + b1 * cos(angle)));
+                }
+                else
+                {
+                    ellipse.push_back(POINT(cx + a2 * sin(angle), cy + b2 * cos(angle)));
+                }
+            }
+
+            ellipses.push_back(ellipse);
+        }
+    }
+
+
+    /***
+    if (compare::MEASURE_CLIP)
+    {
+        for (std::vector<BOX>::const_iterator bit = boxes.begin();
+            bit != boxes.end();
+            ++bit)
+        {
+            double cx = 0.5 * (bit->xmin() + bit->xmax());
+            double cy = 0.5 * (bit->ymin() + bit->ymax());
+
+            double a = compare::OVERLAY_ELLIPSE_FACTOR * 0.5 * (bit->xmax() - bit->xmin());
+            double b = compare::OVERLAY_ELLIPSE_FACTOR * 0.5 * (bit->ymax() - bit->ymin());
+
+            double angle = 45.0 * boost::math::constants::pi<double>() / 180.0;
+
+                BOX box;
+                double angle = 225.0 * ggl::math::d2r;
+                ggl::set<ggl::min_corner, 0>(box, cx + a * sin(angle));
+                ggl::set<ggl::min_corner, 1>(box, cy + b * cos(angle));
+                angle = 45.0 * ggl::math::d2r;
+                ggl::set<ggl::max_corner, 0>(box, cx + a * sin(angle));
+                ggl::set<ggl::max_corner, 1>(box, cy + b * cos(angle));
+
+                .push_back(box);
+
+
+            POLY clip;
+
+            // Create COUNTER-clockwise ellipse (therefore the minus in angle), WITHOUT last point
+            // (therefore the -1)
+            for (int v = 0; v < 4; v++, angle -= compare::delta)
+            {
+                clip.push_back(POINT(cx + a * sin(angle), cy + b * cos(angle)));
+            }
+
+            BOX box = CGAL::bbox_2(clip.vertices_begin(), clip.vertices_end(), K());
+            clip_boxes.push_back(box);
+        }
+    }
+    ***/
+
+
+    compare::wait();
+
+    if (compare::MEASURE_AREA)
+    {
+        double area = 0;
+        boost::timer t;
+        for (int i = 0; i < compare::AREA_COUNT; i++)
+        {
+            for (std::vector<POLY>::const_iterator it = polygons.begin();
+                it != polygons.end();
+                it++)
+            {
+                area += CGAL::to_double(it->area());
+            }
+        }
+        compare::report_area(t, polygons.size(), area);
+    }
+
+    if (compare::MEASURE_CENTROID)
+    {
+        double sum_x = 0, sum_y = 0;
+        boost::timer t;
+        for (int i = 0; i < compare::CENTROID_COUNT; i++)
+        {
+            for (std::vector<POLY>::const_iterator it = polygons.begin();
+                it != polygons.end();
+                ++it)
+            {
+                POINT centroid = CGAL::centroid(it->vertices_begin(), it->vertices_end(), K());
+                sum_x += CGAL::to_double(centroid.x());
+                sum_y += CGAL::to_double(centroid.y());
+            }
+        }
+        compare::report_centroid(t, polygons.size(), sum_x, sum_y);
+    }
+
+
+    if (compare::MEASURE_CONVEX_HULL)
+    {
+        double area = 0.0;
+        boost::timer t;
+        for (std::vector<POLY>::const_iterator it = polygons.begin(); it != polygons.end(); it++)
+        {
+            POLY hull;
+            CGAL::ch_graham_andrew(it->vertices_begin(), it->vertices_end(), std::back_inserter(hull));
+            if (compare::HULL_AREA)
+            {
+                area += CGAL::to_double(hull.area());
+            }
+
+        }
+        compare::report_hull(t, polygons.size(), area);
+    }
+
+
+    /***
+    if (compare::MEASURE_CLIP)
+    {
+        bool first = true;
+        double area1 = 0, area2 = 0;
+
+        boost::timer t;
+
+        std::vector<POLY>::const_iterator cit = clip_boxes.begin();
+        for (std::vector<POLY>::iterator it = polygons.begin();
+            it != polygons.end() && cit != clip_boxes.end();
+            ++it, ++cit)
+        {
+            area1 += CGAL::to_double(it->area());
+
+            std::list<HOLEY_POLY> pv;
+            try
+            {
+                CGAL::intersection(*it, *cit, std::back_inserter(pv));
+                double a = 0;
+                for (std::list<HOLEY_POLY>::const_iterator pit = pv.begin(); pit != pv.end(); ++pit)
+                {
+                    a += CGAL::to_double(pit->outer_boundary().area());
+                }
+                area2 += a;
+                //std::cout << " " << a << std::endl;
+            }
+            catch(std::exception const& e)
+            {
+            }
+            catch(...)
+            {
+            }
+
+        }
+        compare::report_clip(t, polygons.size(), area1, area2);
+    }
+    ***/
+
+
+    #ifndef SKIP_OVERLAY
+    if (compare::MEASURE_OVERLAY)
+    {
+        double area1 = 0, area2 = 0;
+
+        boost::timer t;
+
+        for (int i = 0; i < compare::OVERLAY_COUNT; i++)
+        {
+        int k = 0;
+        std::vector<POLY>::const_iterator eit = ellipses.begin();
+        for (std::vector<POLY>::iterator it = polygons.begin();
+            it != polygons.end() && eit != ellipses.end();
+            ++it, ++eit, ++k)
+        {
+            area1 += CGAL::to_double(it->area());
+
+            std::list<HOLEY_POLY> pv;
+            try
+            {
+                CGAL::intersection(*it, *eit, std::back_inserter(pv));
+                double a = 0;
+                for (std::list<HOLEY_POLY>::const_iterator pit = pv.begin(); pit != pv.end(); ++pit)
+                {
+                    a += CGAL::to_double(pit->outer_boundary().area());
+                }
+                area2 += a;
+            }
+            catch(std::exception const& e)
+            {
+            }
+            catch(...)
+            {
+            }
+
+        }
+        }
+        compare::report_overlay(t, polygons.size(), area1, area2);
+    }
+    #endif
+
+
+    if (compare::MEASURE_WITHIN)
+    {
+        int count = 0;
+        int count_boundary = 0;
+        boost::timer t;
+        for (int e = 0; e < boxes.size(); e++)
+        {
+            const BOX& b = boxes[e];
+            POINT p((b.xmin() + b.xmax()) / 2.0, (b.ymin() + b.ymax()) / 2.0);
+
+            //compare::debug_within(e, count);
+
+            std::vector<BOX>::const_iterator bit = boxes.begin();
+            int k = 0;
+            for (std::vector<POLY>::const_iterator it = polygons.begin();
+                it != polygons.end() && bit != boxes.end();
+                it++, bit++, k++)
+            {
+                if (p.x() > bit->xmin() && p.x() < bit->xmax() && p.y() > bit->ymin() && p.y() < bit->ymax())
+                {
+                    switch(CGAL::bounded_side_2(it->vertices_begin(), it->vertices_end(), p, K()))
+                    {
+                        case CGAL::ON_BOUNDED_SIDE :
+                            count++;
+                            //std::cout << e << " IN " << k << std::endl;
+                            break;
+                        case CGAL::ON_BOUNDARY:
+                            count_boundary++;
+                            break;
+                        case CGAL::ON_UNBOUNDED_SIDE:
+                            break;
+                     }
+                }
+            }
+        }
+        compare::report_within(t, polygons.size(), count, count_boundary);
+    }
+
+    return 0;
+}

Added: sandbox/ggl/other/comparisons/common/common.hpp
==============================================================================
--- (empty file)
+++ sandbox/ggl/other/comparisons/common/common.hpp	2009-08-28 04:54:00 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,335 @@
+// Generic Geometry Library - COMPARISONS
+//
+// Copyright Barend Gehrels, 2008, 2009, Geodan Holding B.V. Amsterdam, the Netherlands.
+// Copyright Bruno Lalande 2008, 2009
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef COMPARE_COMMON_HPP_INCLUDED
+#define COMPARE_COMMON_HPP_INCLUDED
+
+#include <string>
+#include <iostream>
+#include <iomanip>
+#include <sstream>
+
+#include <boost/timer.hpp>
+#include <boost/math/constants/constants.hpp>
+
+//#define OVERLAY_SYMDIFF
+
+
+/*
+To put results into a (e.g. MySQL) database:
+
+drop table if exists measurements;
+create table measurements(id integer auto_increment primary key,
+    lib varchar(16), comp varchar(16), dset varchar(16), alg varchar(16), val double, cnt integer, sec double);
+
+To get the results, do something like:
+select lib,min(cnt) as c,round(min(val),5) as a,round(min(sec),2) as s from measurements
+group by lib
+order by s
+
+
+To get the results in the form which can be pasted into the Wiki-page:
+select alg,concat('||', lib, '||',round(min(val),5),'||',round(min(sec),3),'||') from measurements
+group by alg,lib
+order by alg,lib
+
+
+*/
+
+
+//#define NLP4 // if defined, DUTCH dataset (private not for sale) else PUBLIC US COUNTY dataset
+
+// BEGIN DEFINES for testing ggl-sorting only
+//#define USE_SMOOTH_SORT
+//#define USE_MERGE_SORT
+//#define CONVEX_HULL_USE_INSERTION
+#if defined(USE_MERGE_SORT)
+static const int thread_count = 1;
+#endif
+// END
+
+
+// #define OUTPUT_TEXT // if defined, outputs in text format. If not, outputs SQL insert statements
+
+
+
+namespace compare
+{
+
+    std::string g_name;
+    std::string g_compiler;
+    std::string g_dataset;
+
+    std::string filename(int argc, char** argv)
+    {
+        return argc > 1 ? argv[1] :
+#ifdef NLP4
+         "c:/data/spatial/shape/nl/nlp4_r98"
+#else
+        // Can be downloaded for free from http://www.weather.gov/geodata/catalog/county/html/county100k.htm
+         "c:/data/spatial/shape/c100000jan2902"
+#endif
+    ;
+    }
+
+    std::string fieldname()
+    {
+    // todo: argc,argv
+#ifdef NLP4
+        return "pc4nr";
+#else
+        return "fips";
+#endif
+    }
+
+
+#ifdef NLP4
+    const int AREA_COUNT = 1000;
+    const double AREA_DIVISOR = 1000000.0; // km^2
+
+    const double LENGTH_DIVISOR = 1000.0;
+
+    const double SIMPLIFY_DISTANCE = 91.7; // results in simplified length of 99%, vertex reduction of ca. 50%
+
+    const int CENTROID_COUNT = 1000;
+
+    const double CLIP_DISTANCE = 10.0;
+
+    const int OVERLAY_COUNT = 10;
+    const int CLIP_COUNT = 10;
+
+    const double INTEGER_FACTOR = 100.0;
+
+#else
+    const int AREA_COUNT = 10;
+
+    const double AREA_DIVISOR = 1.0;
+    const double LENGTH_DIVISOR = 1.0;
+
+    const double SIMPLIFY_DISTANCE = 0.0008975; // results in simplified length of 99%, vertex reduction of ca. 10%
+
+    const int CENTROID_COUNT = 10;
+
+
+    const double CLIP_DISTANCE = 0.1;
+
+    const int OVERLAY_COUNT = 1;
+    const int CLIP_COUNT = 1;
+
+    const double INTEGER_FACTOR = 10000.0;
+
+#endif
+
+
+//#define OVERLAY_UNION
+
+    const int OVERLAY_ELLIPSE_COUNT = 101;
+    const double delta = boost::math::constants::pi<double>() * 2.0 / (compare::OVERLAY_ELLIPSE_COUNT - 1);
+    const double OVERLAY_ELLIPSE_FACTOR1 = 1.1; // 1.1
+    const double OVERLAY_ELLIPSE_FACTOR2 = 0.2; // 0.2
+    const double CLIP_FACTOR = 0.9;
+
+    const bool MEASURE_AREA = true;
+    const bool MEASURE_CENTROID = true;
+    const bool MEASURE_CLIP = true;
+    const bool MEASURE_CONVEX_HULL = true;
+    const bool MEASURE_OVERLAY = true;
+    const bool MEASURE_SIMPLIFY = true;
+    const bool MEASURE_TOUCH = false; // currently only for GEOS, not further worked out
+    const bool MEASURE_WITHIN = true;
+
+    // Variables to check the results
+    const bool CLIP_AREA = true;
+    const bool HULL_AREA = true;
+    const bool OVERLAY_AREA = true;
+    const bool SIMPLIFY_LENGTH = true;
+
+
+
+    // Wait a while, after reading shapefiles/generating ellipses, etc
+    // until all buffers/memory allocation/swapping is done.
+    inline void wait(bool first = true)
+    {
+        int const s = first ? 3 : 1;
+        boost::timer t;
+        while (t.elapsed() < s) ;
+    }
+
+
+
+    template <typename T1, typename T2>
+    inline void report_common(double s, int n, const std::string& algorithm, T1 value1, T2 value2 = 0)
+    {
+        #ifdef OUTPUT_TEXT
+        std::cout << algorithm << " " << n
+            << " total: " << value1
+            << " other: " << value2
+            << " time: "<< s << "s" << std::endl;
+        #else
+
+        std::cout << "insert into measurements(lib, comp, dset, alg, val, cnt, sec)"
+            << " values("
+            << "'" << g_name << "'"
+            << ", '" << g_compiler << "'"
+            << ", '" << g_dataset << "'"
+            << ", '" << algorithm << "'"
+            << ", " << value1
+            << ", " << n
+            << ", " << s
+            << ");" << std::endl;
+        #endif
+
+
+        wait(false);
+    }
+
+    inline void report_area(const boost::timer& t, int n, double area)
+    {
+        double s = t.elapsed();
+        area /= double(compare::AREA_COUNT);
+        area /= compare::AREA_DIVISOR;
+        report_common(s, n, "AREA", area, 0);
+    }
+
+    inline void report_centroid(const boost::timer& t, int n, double sum_x, double sum_y)
+    {
+        double s = t.elapsed();
+        sum_x /= double(compare::CENTROID_COUNT);
+        sum_y /= double(compare::CENTROID_COUNT);
+        sum_x /= double(n);
+        sum_y /= double(n);
+        report_common(s, n, "CENTROID", sum_x, sum_y);
+    }
+
+    inline void report_hull(const boost::timer& t, int n, double area)
+    {
+        double s = t.elapsed();
+        report_common(s, n, "HULL", area / compare::AREA_DIVISOR, 0);
+    }
+    inline void report_simplify(const boost::timer& t, int n, double length1, double length2, int count1, int count2)
+    {
+        double s = t.elapsed();
+        length1 /= compare::LENGTH_DIVISOR;
+        length2 /= compare::LENGTH_DIVISOR;
+
+#ifdef OUTPUT_TEXT
+        std::cout
+            << "original : length=" << length1 << ", count=" << count1
+            << " -> "
+            << "simplified : length=" << length2 << ", count=" << count2
+            << std::endl;
+#endif
+
+        report_common(s, n, "SIMPLIFY", length2, length1);
+    }
+
+    inline void report_within(const boost::timer& t, int n, int count, int count_boundary)
+    {
+        double s = t.elapsed();
+            //<< " " << count << " (" << count_boundary << ")"
+        report_common(s, n, "WITHIN", count, count_boundary);
+    }
+
+    inline void report_touch(const boost::timer& t, int n, int count, int count_box)
+    {
+        double s = t.elapsed();
+            //<< " " << count << " (" << count_boundary << ")"
+        report_common(s, n, "TOUCH", count, count_box);
+    }
+
+    inline void report_overlay(const boost::timer& t, int n, double area1, double area2)
+    {
+        double s = t.elapsed();
+
+        area1 /= compare::OVERLAY_COUNT;
+        area2 /= compare::OVERLAY_COUNT;
+        area1 /= compare::AREA_DIVISOR;
+        area2 /= compare::AREA_DIVISOR;
+        report_common(s, OVERLAY_ELLIPSE_COUNT,
+#if defined(OVERLAY_UNION)
+            "UNION"
+#elif defined(OVERLAY_SYMDIFF)
+            "SYM_DIFFERENCE"
+#else
+            "INTERSECTION"
+#endif
+            , area2, area1);
+    }
+
+    inline void report_clip(const boost::timer& t, int n, double area1, double area2)
+    {
+        double s = t.elapsed();
+
+        area1 /= compare::CLIP_COUNT;
+        area2 /= compare::CLIP_COUNT;
+        area1 /= compare::AREA_DIVISOR;
+        area2 /= compare::AREA_DIVISOR;
+        report_common(s, n, "CLIP", area2, area1);
+    }
+
+
+
+    inline void debug_within(int i, int count)
+    {
+        std::cout << ".";
+        if ((i + 1) % 70 == 0)
+        {
+            std::cout << count << std::endl;
+        }
+    }
+
+
+    void version_info(const std::string& name)
+    {
+        g_name = name;
+
+        #ifdef _STLPORT_VERSION
+            std::cout << "STLPORT " << _STLPORT_VERSION << std::endl;
+        #endif
+
+        std::ostringstream out;
+
+        #ifdef _MSC_VER
+            out << "MSC " << _MSC_VER;
+        #endif
+
+        #if defined(__GNUC__)
+            out << "GCC " << __GNUC__ << "." << __GNUC_MINOR__ << "." << __GNUC_PATCHLEVEL__;
+        #endif
+
+        g_compiler = out.str();
+
+        #ifdef NLP4
+            g_dataset = "NLP4";
+        #else
+            g_dataset = "US Counties";
+        #endif
+
+        #ifdef OUTPUT_TEXT
+        std::cout << g_compiler << std::endl;
+        std::cout << g_dataset << std::endl;
+        #endif
+
+
+        #if defined(CONVEX_HULL_USE_INSERTION)
+            std::cout << "insertions" << std::endl;
+        #elif defined(USE_SMOOTH_SORT)
+            std::cout << "smoothsort::sort" << std::endl;
+        #elif defined(USE_MERGE_SORT)
+            std::cout << "MergeSort::sort (" << thread_count << " thread"
+                << (thread_count == 1 ? "" : "s") << ")" << std::endl;
+        #else
+            //std::cout << "std::sort" << std::endl;
+        #endif
+        std::cout << std::setprecision(16) << std::endl;
+    }
+}
+
+
+#endif // COMPARE_COMMON_HPP_INCLUDED
+

Added: sandbox/ggl/other/comparisons/common/read_shapefile.hpp
==============================================================================
--- (empty file)
+++ sandbox/ggl/other/comparisons/common/read_shapefile.hpp	2009-08-28 04:54:00 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,81 @@
+// Generic Geometry Library - COMPARISONS
+//
+// Copyright Barend Gehrels, 2008, 2009, Geodan Holding B.V. Amsterdam, the Netherlands.
+// Copyright Bruno Lalande 2008, 2009
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef COMPARE_READ_SHAPEFILE_HPP_INCLUDED
+#define COMPARE_READ_SHAPEFILE_HPP_INCLUDED
+
+
+
+// If this does NOT compile:
+// -> download shapelib from http://shapelib.maptools.org/
+#include "shapefil.h"
+
+#include <vector>
+#include <string>
+
+template <typename T, typename F>
+void read_shapefile(std::string const& filename, std::string const & idfieldname,
+                    std::vector<T>& polygons, std::vector<int>& ids,
+                    F functor)
+{
+    try
+    {
+        std::string shape_filename = filename + ".shp";
+        std::string dbase_filename = filename + ".dbf";
+
+        SHPHandle handle = SHPOpen(shape_filename.c_str(), "rb");
+        if (handle <= 0)
+        {
+            throw std::string("File " + shape_filename + " not found");
+        }
+        DBFHandle dbf = DBFOpen(dbase_filename.c_str(), "rb+");
+        if (dbf <= 0)
+        {
+            throw std::string("File " + dbase_filename + " not found");
+        }
+
+        int nShapeType, nEntities;
+        double adfMinBound[4], adfMaxBound[4];
+        SHPGetInfo(handle, &nEntities, &nShapeType, adfMinBound, adfMaxBound );
+
+        int nr_index = DBFGetFieldIndex(dbf, idfieldname.c_str());
+
+        for (int i = 0; i < nEntities; i++)
+        {
+            SHPObject* psShape = SHPReadObject(handle, i );
+
+            // Process only polygons, and from them only single-polygons without holes
+            if (psShape->nSHPType == SHPT_POLYGON && psShape->nParts == 1)
+            {
+                polygons.push_back(functor(psShape));
+                ids.push_back(DBFReadIntegerAttribute(dbf, i, nr_index));
+            }
+            SHPDestroyObject( psShape );
+        }
+        SHPClose(handle);
+        DBFClose(dbf);
+    }
+    catch(std::exception const& e)
+    {
+        throw e;
+    }
+    catch(const std::string& s)
+    {
+        throw s;
+    }
+    catch(...)
+    {
+        throw std::string("Unknown exception ('...')");
+    }
+}
+
+
+
+
+#endif // COMPARE_READ_SHAPEFILE_HPP_INCLUDED
+

Added: sandbox/ggl/other/comparisons/ggl/ggl_check.cpp
==============================================================================
--- (empty file)
+++ sandbox/ggl/other/comparisons/ggl/ggl_check.cpp	2009-08-28 04:54:00 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,491 @@
+// Generic Geometry Library - COMPARISONS
+//
+// Copyright Barend Gehrels, 2008, 2009, Geodan Holding B.V. Amsterdam, the Netherlands.
+// Copyright Bruno Lalande 2008, 2009
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+
+
+//#define GGL_DEBUG_INTERSECTION
+// #define OUTPUT_POSTGIS
+// #define OUTPUT_SQLSERVER
+//#define OUTPUT_ORACLE
+//#define OUTPUT_WKT
+
+
+
+#include "../common/common.hpp"
+#include "../common/read_shapefile.hpp"
+
+#include <ggl/ggl.hpp>
+#include <ggl/geometries/geometries.hpp>
+
+#include <ggl/io/oracle/write_oracle.hpp>
+
+
+
+template <typename T>
+T convert(SHPObject* psShape)
+{
+    const double* x = psShape->padfX;
+    const double* y = psShape->padfY;
+    T polygon;
+    for (int v = 0; v < psShape->nVertices; v++)
+    {
+        typename ggl::point_type<T>::type p;
+        ggl::set<0>(p, x[v]);
+        ggl::set<1>(p, y[v]);
+        polygon.outer().push_back(p);
+    }
+    return polygon;
+}
+
+
+#if defined(OUTPUT_POSTGIS) || defined(OUTPUT_SQLSERVER) || defined(OUTPUT_ORACLE)
+
+template <typename Polygon>
+void output_for_sql(int k, Polygon const& polygon1, Polygon const& polygon2, bool& first, bool last)
+{
+#ifdef OVERLAY_UNION
+    const std::string operation = "union";
+#else
+    const std::string operation = "intersection";
+#endif
+
+
+
+    if (first)
+    {
+
+#if defined(OUTPUT_ORACLE)
+        std::cout << "create table t_comp(k integer, a mdsys.sdo_geometry, b mdsys.sdo_geometry);"
+            << std::endl << std::endl
+            << "declare g1 mdsys.sdo_geometry;" << std::endl
+            << "declare g2 mdsys.sdo_geometry;" << std::endl
+            << "begin" << std::endl;
+
+#else
+        std::cout << "drop table if exists t_comp;" << std::endl;
+        std::cout << "create table t_comp(k integer, a geometry, b geometry);" << std::endl;
+#endif
+
+        first = false;
+    }
+
+    #if defined(OUTPUT_POSTGIS)
+    std::cout << "insert into t_comp values(" << k << ", "
+        << "GeometryFromText('" <<  ggl::wkt(polygon1)
+        << "'),GeometryFromText('" << ggl::wkt(polygon2)
+        << "'));"
+        << std::endl;
+
+    // The query is:
+    // select sum(area(intersection(a,b))) from t_comp
+    #endif
+
+    #if defined(OUTPUT_SQLSERVER)
+    std::cout << "insert into t_comp values(" << k
+        << ", geometry::STGeomFromText('" << ggl::wkt(polygon1)
+        << "', 0), geometry::STGeomFromText('" << ggl::wkt(polygon2)
+        << "', 0));"
+        << std::endl;
+
+
+    // The query is:
+    // select sum(a.STIntersection(b).STArea()) as sum_area from t_comp
+    // where k <> 1911 (for 101 points)
+    // 101: reported 00:00:10 (10 seconds)
+    // 1001: reported 00:01:55, area: 161.239301428967 (GGL: 161.2393008999)
+    #endif
+
+    #if defined(OUTPUT_ORACLE)
+    std::cout
+        << " g1 := " << ggl::sdo_oracle(polygon1, 32662) << ";" << std::endl
+        << " g2 := " << ggl::sdo_oracle(polygon2, 32662) << ";" << std::endl
+        << " execute immediate 'insert into t_comp(k,a,b) values(" << k << ",:gm1,:gm2)' using g1,g2;" << std::endl;
+
+
+    // The query is:
+    // select sum(a.STIntersection(b).STArea()) as sum_area from t_comp
+    // where k <> 1911 (for 101 points)
+    // 101: reported 00:00:10 (10 seconds)
+    // 1001: reported 00:01:55, area: 161.239301428967 (GGL: 161.2393008999)
+    #endif
+
+    if (last)
+    {
+        #if defined(OUTPUT_ORACLE)
+        std::cout
+            << "end;" << std::endl
+            << "/" << std::endl;
+        #endif
+    }
+
+
+
+}
+#endif
+
+
+int main(int argc, char** argv)
+{
+    compare::version_info("ggl");
+
+    typedef ggl::point_xy<double> POINT;
+    typedef ggl::polygon<POINT> POLY;
+    typedef ggl::box<POINT> BOX;
+
+    std::vector<POLY> polygons;
+    std::vector<BOX> boxes;
+    std::vector<BOX> clip_boxes; // for clips
+    std::vector<POLY> ellipses; // for intersections/unions
+    std::vector<int> ids;
+    try
+    {
+        read_shapefile(compare::filename(argc, argv), compare::fieldname(), polygons, ids, convert<POLY>);
+    }
+    catch(std::exception const& e)
+    {
+        std::cout << e.what() << std::endl;
+        return 1;
+    }
+    catch(std::string const* s)
+    {
+        std::cout << s << std::endl;
+        return 1;
+    }
+    catch(...)
+    {
+        std::cout << "Unknown exception" << std::endl;
+        return 1;
+    }
+
+
+    // Create envelopes
+    for (std::vector<POLY>::const_iterator it = polygons.begin();
+        it != polygons.end();
+        ++it)
+    {
+        BOX b;
+        ggl::envelope(*it, b);
+        boxes.push_back(b);
+    }
+
+    // Create the ellipses for intersections lateron
+    if (compare::MEASURE_OVERLAY || compare::MEASURE_CLIP)
+    {
+
+        int k = 0;
+        for (std::vector<BOX>::const_iterator bit = boxes.begin();
+            bit != boxes.end();
+            ++bit, ++k)
+        {
+            double cx = 0.5 * (ggl::get<ggl::max_corner, 0>(*bit) + ggl::get<ggl::min_corner, 0>(*bit));
+            double cy = 0.5 * (ggl::get<ggl::max_corner, 1>(*bit) + ggl::get<ggl::min_corner, 1>(*bit));
+
+            double dx = ggl::get<ggl::max_corner, 0>(*bit) - ggl::get<ggl::min_corner, 0>(*bit);
+            double dy = ggl::get<ggl::max_corner, 1>(*bit) - ggl::get<ggl::min_corner, 1>(*bit);
+
+            if (compare::MEASURE_OVERLAY)
+            {
+                double a1 = compare::OVERLAY_ELLIPSE_FACTOR1 * 0.5 * dx;
+                double b1 = compare::OVERLAY_ELLIPSE_FACTOR1 * 0.5 * dy;
+                double a2 = compare::OVERLAY_ELLIPSE_FACTOR2 * 0.5 * dx;
+                double b2 = compare::OVERLAY_ELLIPSE_FACTOR2 * 0.5 * dy;
+
+                POLY ellipse;
+                ellipse.outer().reserve(compare::OVERLAY_ELLIPSE_COUNT);
+                double angle = 0.0; //45.0 * ggl::math::d2r; //0.0;
+                for (int i = 0; i < compare::OVERLAY_ELLIPSE_COUNT - 1; i++, angle += compare::delta)
+                {
+                    if (i % 2 == 0)
+                    {
+                        ellipse.outer().push_back(POINT(cx + a1 * sin(angle), cy + b1 * cos(angle)));
+                    }
+                    else
+                    {
+                        ellipse.outer().push_back(POINT(cx + a2 * sin(angle), cy + b2 * cos(angle)));
+                    }
+                }
+                ellipse.outer().push_back(ellipse.outer().front());
+                ellipses.push_back(ellipse);
+            }
+
+            if (compare::MEASURE_CLIP)
+            {
+                double a = compare::CLIP_FACTOR * 0.5 * dx;
+                double b = compare::CLIP_FACTOR * 0.5 * dy;
+
+                BOX box;
+                double angle = 225.0 * ggl::math::d2r;
+                ggl::set<ggl::min_corner, 0>(box, cx + a * sin(angle));
+                ggl::set<ggl::min_corner, 1>(box, cy + b * cos(angle));
+                angle = 45.0 * ggl::math::d2r;
+                ggl::set<ggl::max_corner, 0>(box, cx + a * sin(angle));
+                ggl::set<ggl::max_corner, 1>(box, cy + b * cos(angle));
+
+                clip_boxes.push_back(box);
+
+                /*if (k < 5)
+                {
+                    std::cout << ggl::wkt(box) << " " << ggl::area(box) << std::endl;
+                }*/
+
+            }
+
+        }
+    }
+
+    compare::wait();
+
+
+    if (compare::MEASURE_AREA)
+    {
+        double area = 0;
+        boost::timer t;
+        for (int i = 0; i < compare::AREA_COUNT; i++)
+        {
+            for (std::vector<POLY>::const_iterator it = polygons.begin();
+                it != polygons.end();
+                ++it)
+            {
+                area += ggl::area(*it);
+            }
+        }
+        compare::report_area(t, polygons.size(), area);
+    }
+
+
+    if (compare::MEASURE_CENTROID)
+    {
+        double sum_x = 0, sum_y = 0;
+        boost::timer t;
+        for (int i = 0; i < compare::CENTROID_COUNT; i++)
+        {
+            for (std::vector<POLY>::const_iterator it = polygons.begin();
+                it != polygons.end();
+                ++it)
+            {
+                POINT centroid;
+#if defined(CENTROID_WITH_CATCH)
+                try
+#endif
+                {
+                    ggl::centroid(*it, centroid);
+                    sum_x += centroid.x();
+                    sum_y += centroid.y();
+                }
+#if defined(CENTROID_WITH_CATCH)
+                catch(std::exception const& e)
+                {
+                }
+                catch(...)
+                {
+                }
+#endif
+            }
+        }
+        compare::report_centroid(t, polygons.size(), sum_x, sum_y);
+    }
+
+    if (compare::MEASURE_CONVEX_HULL)
+    {
+        double area = 0.0;
+        boost::timer t;
+        for (std::vector<POLY>::const_iterator it = polygons.begin(); it != polygons.end(); ++it)
+        {
+            POLY::ring_type ring;
+            //std::cout << ggl::as_wkt<POLY>(*it) << std::endl;
+            ggl::convex_hull(*it, std::back_inserter(ring));
+            if (compare::HULL_AREA)
+            {
+                area += fabs(ggl::area(ring));
+                /*POLY p;
+                p.outer() = ring;
+                std::cout << ggl::as_wkt<POLY>(p) << " ";
+                std::cout << ggl::area(it->outer()) << " -> " << ggl::area(ring)  << std::endl;
+                */
+            }
+
+        }
+        compare::report_hull(t, polygons.size(), area);
+    }
+
+    if (compare::MEASURE_OVERLAY)
+    {
+        bool first = true;
+        double area1 = 0.0, area2 = 0.0;
+
+        boost::timer t;
+        for (int i = 0; i < compare::OVERLAY_COUNT; i++)
+        {
+            int k = 0;
+            std::vector<POLY>::const_iterator eit = ellipses.begin();
+            for (std::vector<POLY>::const_iterator it = polygons.begin();
+                it != polygons.end() && eit != ellipses.end();
+                ++it, ++eit, ++k)
+            {
+                //if (k == 1590)
+                {
+
+#if defined(OUTPUT_WKT)
+                    if (i == 0)// && k == 264)
+                    {
+                        std::cout
+                            << ggl::wkt(*it) << std::endl
+                            << ggl::wkt(*eit) << std::endl;
+                    }
+#endif
+
+                    if (compare::OVERLAY_AREA)
+                    {
+                        area1 += ggl::area(*it);
+                    }
+
+
+
+                    POLY p;
+                    std::vector<POLY> v;
+                    ggl::intersection<POLY>(*eit, *it, std::back_inserter(v));
+
+                    if (compare::OVERLAY_AREA) // && k != 1911)
+                    {
+                        double a = 0.0;
+                        for (std::vector<POLY>::const_iterator pit = v.begin(); pit != v.end(); ++pit)
+                        {
+                            a += ggl::area(*pit);
+                        }
+                        area2 += a;
+                        //std::cout << "insert into se_ggl values(" << k << ", " << a << ");" << std::endl;
+                    }
+
+#if defined(OUTPUT_WKT)
+                    if (i == 0)// && k == 264)
+                    {
+                        for (std::vector<POLY>::const_iterator pit = v.begin(); pit != v.end(); ++pit)
+                        {
+                            std::cout << ggl::wkt(*pit) << std::endl;
+                        }
+                        //return 0;
+                    }
+#endif
+
+
+#if defined(OUTPUT_POSTGIS) || defined(OUTPUT_SQLSERVER) || defined(OUTPUT_ORACLE)
+                    if (i == 0)
+                    {
+                        output_for_sql(k, *it, *eit, first, it + 1 == polygons.end());
+                    }
+#endif
+                }
+            }
+
+        }
+        compare::report_overlay(t, polygons.size(), area1, area2);
+    }
+
+    if (compare::MEASURE_CLIP)
+    {
+        bool first = true;
+        double area1 = 0.0, area2 = 0.0;
+
+       boost::timer t;
+        for (int i = 0; i < compare::CLIP_COUNT; i++)
+        {
+            int k = 0;
+            std::vector<BOX>::const_iterator bit = clip_boxes.begin();
+            for (std::vector<POLY>::const_iterator it = polygons.begin();
+                it != polygons.end() && bit != clip_boxes.end();
+
+                ++it, ++bit,  ++k)
+            {
+                if (compare::CLIP_AREA)
+                {
+                    area1 += ggl::area(*it);
+                }
+
+                POLY p;
+                std::vector<POLY> v;
+                ggl::intersection<POLY>(*bit, *it, std::back_inserter(v));
+
+                if (compare::CLIP_AREA)
+                {
+                    double a = 0.0;
+                    for (std::vector<POLY>::const_iterator pit = v.begin(); pit != v.end(); ++pit)
+                    {
+                        a += ggl::area(*pit);
+                    }
+                    area2 += a;
+                }
+            }
+        }
+        compare::report_clip(t, polygons.size(), area1, area2);
+    }
+
+
+    if (compare::MEASURE_SIMPLIFY)
+    {
+        int count1 = 0, count2 = 0;
+        double length1 = 0.0, length2 = 0.0;
+        boost::timer t;
+        for (std::vector<POLY>::const_iterator it = polygons.begin(); it != polygons.end(); ++it)
+        {
+            POLY p;
+            ggl::simplify(*it, p, compare::SIMPLIFY_DISTANCE);
+            count1 += it->outer().size();
+            count2 += p.outer().size();
+            if (compare::SIMPLIFY_LENGTH)
+            {
+                length1 += ggl::perimeter(it->outer());
+                length2 += ggl::perimeter(p.outer());
+            }
+
+        }
+        compare::report_simplify(t, polygons.size(), length1, length2, count1, count2);
+    }
+
+
+    if (compare::MEASURE_WITHIN)
+    {
+        int count = 0;
+        boost::timer t;
+        for (int e = 0; e < boxes.size(); e++)
+        {
+            const BOX& b = boxes[e];
+            POINT p((b.min_corner().x() + b.max_corner().x()) / 2.0,
+                (b.min_corner().y() + b.max_corner().y()) / 2.0);
+
+            //compare::debug_within(e, count);
+
+            std::vector<BOX>::const_iterator bit = boxes.begin();
+            int k = 0;
+            for (std::vector<POLY>::const_iterator it = polygons.begin();
+                it != polygons.end() && bit != boxes.end();
+                ++it, ++bit, ++k)
+            {
+                /***
+                if ((e == 85 && k == 85)
+                    || e == 99 && k == 99
+                    )
+                {
+                    std::cout << e << "," << k << std::endl;
+                    std::cout << std::setprecision(16) << ggl::wkt(p) << std::endl;
+                    std::cout << std::setprecision(16) << ggl::wkt(*it) << std::endl;
+                }
+                ***/
+                if (ggl::within(p, *bit) && ggl::within(p, *it))
+                {
+                    //std::cout << e << " IN " << k << std::endl;
+                    count++;
+                }
+            }
+        }
+        compare::report_within(t, polygons.size(), count, -1);
+    }
+
+
+
+    return 0;
+}

Added: sandbox/ggl/other/comparisons/gpc/gpc_check.cpp
==============================================================================
--- (empty file)
+++ sandbox/ggl/other/comparisons/gpc/gpc_check.cpp	2009-08-28 04:54:00 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,289 @@
+
+
+#include "../common/common.hpp"
+#include "../common/read_shapefile.hpp"
+
+
+#include <cstdlib>
+#include <cmath>
+
+extern "C"
+{
+#include "gpc.h"
+}
+
+#include <iomanip>
+
+
+
+template <typename P>
+inline void allocate_polygon(P& polygon, int n)
+{
+    polygon.hole = (int*) malloc(1 * sizeof(int));
+    polygon.contour =  (gpc_vertex_list*) malloc(1 * sizeof(gpc_vertex_list));
+    polygon.contour[0].vertex = (gpc_vertex*) malloc(n * sizeof(gpc_vertex));
+
+    polygon.num_contours = 1;
+    polygon.hole[0] = 0;
+    polygon.contour[0].num_vertices = n;
+}
+
+
+template <typename T>
+T convert(SHPObject* psShape)
+{
+    const double* x = psShape->padfX;
+    const double* y = psShape->padfY;
+
+    int n = psShape->nVertices;
+
+    T polygon;
+
+    allocate_polygon(polygon, n);
+
+    // Add it in CLOCKWISE direction
+    for (int v = 0; v < n; v++, x++, y++)
+    {
+        polygon.contour[0].vertex[v].x = *x;
+        polygon.contour[0].vertex[v].y = *y;
+
+    }
+    return polygon;
+}
+
+
+template <typename P, typename E>
+void envelope_for_gpc(P const& polygon, E& envelope)
+{
+    envelope.first.x = 1e300;
+    envelope.first.y = 1e300;
+    envelope.second.x = -1e300;
+    envelope.second.y = -1e300;
+    int n = polygon.contour[0].num_vertices;
+    const gpc_vertex* v = polygon.contour[0].vertex;
+    for (int i = 0; i < n; i++, v++)
+    {
+        if (v->x < envelope.first.x) envelope.first.x = v->x;
+        if (v->y < envelope.first.y) envelope.first.y = v->y;
+        if (v->x > envelope.second.x) envelope.second.x = v->x;
+        if (v->y > envelope.second.y) envelope.second.y = v->y;
+    }
+}
+
+template <typename P>
+inline double area_for_gpc(P const& polygon, bool close = false)
+{
+    double a = 0.0;
+    const int n = polygon.num_contours;
+    for (int c = 0; c < n; c++)
+    {
+        int n = polygon.contour[c].num_vertices;
+        const gpc_vertex* v = polygon.contour[c].vertex;
+        const gpc_vertex* w = v++;
+        for (int i = 1; i < n; i++, w = v++)
+        {
+            a += w->x * v->y - v->x * w->y;
+        }
+        if (close)
+        {
+            // add closing point
+            v = polygon.contour[c].vertex;
+            a += w->x * v->y - v->x * w->y;
+        }
+    }
+    return std::abs(a / 2.0);
+}
+
+
+
+template <typename B, typename P>
+inline void box_as_polygon(B const& box, P& polygon)
+{
+    int n = 5;
+
+    allocate_polygon(polygon, 5);
+
+    polygon.contour[0].vertex[0].x = box.first.x;
+    polygon.contour[0].vertex[0].y = box.first.y;
+
+    polygon.contour[0].vertex[1].x = box.first.x;
+    polygon.contour[0].vertex[1].y = box.second.y;
+
+    polygon.contour[0].vertex[2].x = box.second.x;
+    polygon.contour[0].vertex[2].y = box.second.y;
+
+    polygon.contour[0].vertex[3].x = box.second.x;
+    polygon.contour[0].vertex[3].y = box.first.y;
+
+    polygon.contour[0].vertex[4].x = box.first.x;
+    polygon.contour[0].vertex[4].y = box.first.y;
+}
+
+
+
+
+
+int main(int argc, char** argv)
+{
+    compare::version_info("gpc");
+
+    typedef gpc_vertex POINT;
+    typedef gpc_polygon POLY;
+    typedef std::pair<POINT, POINT> BOX;
+
+    std::vector<POLY> polygons;
+    std::vector<BOX> boxes;
+    std::vector<BOX> boxes10; // for clips
+    std::vector<POLY> ellipses; // for intersections/unions
+    std::vector<int> ids;
+    try
+    {
+        read_shapefile(compare::filename(argc, argv), compare::fieldname(), polygons, ids, convert<POLY>);
+    }
+    catch(std::exception const& e)
+    {
+        std::cout << e.what() << std::endl;
+        return 1;
+    }
+    catch(std::string const* s)
+    {
+        std::cout << s << std::endl;
+        return 1;
+    }
+    catch(...)
+    {
+        std::cout << "Unknown exception" << std::endl;
+        return 1;
+    }
+
+
+    // Create envelopes
+    for (std::vector<POLY>::const_iterator it = polygons.begin();
+        it != polygons.end();
+        ++it)
+    {
+        BOX b;
+        envelope_for_gpc(*it, b);
+        boxes.push_back(b);
+
+        // Decrease for clipping
+        b.first.x += compare::CLIP_DISTANCE;
+        b.first.y += compare::CLIP_DISTANCE;
+        b.second.x -= compare::CLIP_DISTANCE;
+        b.second.y -= compare::CLIP_DISTANCE;
+
+        boxes10.push_back(b);
+    }
+
+
+    // Create the ellipses for intersections lateron
+    if (compare::MEASURE_OVERLAY)
+    {
+        ellipses.resize(boxes.size());
+
+        int i = 0;
+        for (std::vector<BOX>::const_iterator bit = boxes.begin();
+            bit != boxes.end();
+            ++bit, i++)
+        {
+            double cx = 0.5 * (bit->first.x + bit->second.x);
+            double cy = 0.5 * (bit->first.y + bit->second.y);
+
+            double dx = bit->second.x - bit->first.x;
+            double dy = bit->second.y - bit->first.y;
+
+            double a1 = compare::OVERLAY_ELLIPSE_FACTOR1 * 0.5 * dx;
+            double b1 = compare::OVERLAY_ELLIPSE_FACTOR1 * 0.5 * dy;
+            double a2 = compare::OVERLAY_ELLIPSE_FACTOR2 * 0.5 * dx;
+            double b2 = compare::OVERLAY_ELLIPSE_FACTOR2 * 0.5 * dy;
+
+            double angle = 0.0;
+
+
+            POLY& ellipse = ellipses[i];
+            allocate_polygon(ellipse, compare::OVERLAY_ELLIPSE_COUNT);
+
+            for (int v = 0; v < compare::OVERLAY_ELLIPSE_COUNT - 1; v++, angle += compare::delta)
+            {
+                if (v % 2 == 0)
+                {
+                    ellipse.contour[0].vertex[v].x = cx + a1 * sin(angle);
+                    ellipse.contour[0].vertex[v].y = cy + b1 * cos(angle);
+                }
+                else
+                {
+                    ellipse.contour[0].vertex[v].x = cx + a2 * sin(angle);
+                    ellipse.contour[0].vertex[v].y = cy + b2 * cos(angle);
+                }
+            }
+            // close
+            int last = compare::OVERLAY_ELLIPSE_COUNT - 1;
+            ellipse.contour[0].vertex[last].x = ellipse.contour[0].vertex[0].x;
+            ellipse.contour[0].vertex[last].y = ellipse.contour[0].vertex[0].y;
+        }
+    }
+
+
+    compare::wait();
+
+
+    // OVERLAY
+    if (compare::MEASURE_OVERLAY)
+    {
+        bool first = true;
+        double area1 = 0.0, area2 = 0.0;
+
+        std::cout << std::setprecision(16);
+
+        boost::timer t;
+        for (int i = 0; i < compare::OVERLAY_COUNT; i++)
+        {
+
+            int k = 0;
+            //std::vector<BOX>::const_iterator bit = boxes10.begin();
+            std::vector<POLY>::iterator eit = ellipses.begin();
+            for (std::vector<POLY>::iterator it = polygons.begin();
+                it != polygons.end() && eit != ellipses.end();
+                ++it, ++eit, ++k)
+            {
+                if (compare::OVERLAY_AREA)
+                {
+                    area1 += area_for_gpc(*it);
+                }
+
+                POLY& subject = *it;
+                POLY& clip = *eit;
+                POLY result;
+
+                //box_as_polygon(*bit, clip);
+                gpc_polygon_clip(
+#if defined(OVERLAY_UNION)
+                    GPC_UNION
+#elif defined(OVERLAY_SYMDIFF)
+                    GPC_XOR
+#else
+                    GPC_INT
+#endif
+                    , &subject, &clip, &result);
+
+                if (compare::OVERLAY_AREA)
+                {
+                    area2 += area_for_gpc(result, true);
+                }
+
+                //gpc_free_polygon(&clip);
+                gpc_free_polygon(&result);
+            }
+        }
+        compare::report_overlay(t, polygons.size(), area1, area2);
+    }
+
+
+    for (int i = 0; i < polygons.size(); i++)
+    {
+        gpc_free_polygon(&(polygons[i]));
+        gpc_free_polygon(&(ellipses[i]));
+    }
+
+    return 0;
+}

Added: sandbox/ggl/other/comparisons/gtl/gtl_check.cpp
==============================================================================
--- (empty file)
+++ sandbox/ggl/other/comparisons/gtl/gtl_check.cpp	2009-08-28 04:54:00 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,301 @@
+// Generic Geometry Library - COMPARISONS
+//
+// Copyright Barend Gehrels, 2008, 2009, Geodan Holding B.V. Amsterdam, the Netherlands.
+// Copyright Bruno Lalande 2008, 2009
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+
+
+
+#include "../common/common.hpp"
+#include "../common/read_shapefile.hpp"
+
+#include <boost/polygon/polygon.hpp>
+
+
+
+
+// using "double" for polygons with GTL-Booleans
+// - causes a crash in MinGW [[still AUG 25]]
+// - causes a segmentation fault in Linux Ubuntu (gcc 4.3.2)
+// - causes a crash in MSVC2008 (An unhandled win32 exception occured in gtl_check2008.exe [7544].) [[still AUG 25]]
+
+//#define BOOLEANS_DOUBLE
+
+#ifdef BOOLEANS_DOUBLE
+typedef double type_for_booleans;
+const int GTL_INTEGER_FACTOR = 1.0;
+#else
+typedef int type_for_booleans;
+const int GTL_INTEGER_FACTOR = compare::INTEGER_FACTOR;
+#endif
+
+const int GTL_INTEGER_FACTOR_SQR = GTL_INTEGER_FACTOR * GTL_INTEGER_FACTOR;
+
+
+
+// "Factor" is necessary to round coordinates
+template <typename P, typename IP>
+void convert_polygon(P const& double_polygon, IP& integer_polygon)
+{
+    typedef boost::polygon::point_data<type_for_booleans> point_type;
+    std::vector<point_type> points;
+
+    typename boost::polygon::polygon_traits<P>::iterator_type it = double_polygon.begin();
+    int n = double_polygon.size();
+    for(unsigned int i = 0; i < n; i++, ++it)
+    {
+         points.push_back(point_type(
+                GTL_INTEGER_FACTOR * it->x(),
+                GTL_INTEGER_FACTOR * it->y()));
+    }
+
+    boost::polygon::polygon_data<type_for_booleans> polygon;
+    polygon.set(points.begin(), points.end());
+    integer_polygon.insert(polygon);
+}
+
+
+
+template <typename T>
+inline T convert(SHPObject* psShape)
+{
+    const double* x = psShape->padfX;
+    const double* y = psShape->padfY;
+
+    std::vector<boost::polygon::point_data<double> > points;
+
+    for (int v = 0; v < psShape->nVertices; v++)
+    {
+        boost::polygon::point_data<double> pt(x[v], y[v]);
+        points.push_back(pt);
+    }
+
+    boost::polygon::polygon_data<double> polygon;
+    polygon.set(points.begin(), points.end());
+
+    return polygon;
+}
+
+
+
+int main(int argc, char** argv)
+{
+    try
+    {
+        compare::version_info("gtl");
+
+        typedef boost::polygon::polygon_data<double> POLY;
+        typedef boost::polygon::rectangle_data<double> BOX;
+        typedef boost::polygon::point_data<double> POINT;
+
+        std::vector<POLY> polygons;
+        std::vector<POLY> ellipses;
+        std::vector<BOX> boxes;
+        std::vector<int> ids;
+
+        read_shapefile(compare::filename(argc, argv), compare::fieldname(), polygons, ids, convert<POLY>);
+
+
+        // Create envelopes
+        for (std::vector<POLY>::const_iterator it = polygons.begin();
+            it != polygons.end();
+            ++it)
+        {
+            BOX b;
+            boost::polygon::extents(b, *it);
+            boxes.push_back(b);
+        }
+
+
+        // Create the ellipses for intersections lateron
+        if (compare::MEASURE_OVERLAY)
+        {
+            for (std::vector<BOX>::const_iterator bit = boxes.begin();
+                bit != boxes.end();
+                ++bit)
+            {
+                double x1 = bit->get(boost::polygon::HORIZONTAL).low();
+                double y1 = bit->get(boost::polygon::VERTICAL).low();
+                double x2 = bit->get(boost::polygon::HORIZONTAL).high();
+                double y2 = bit->get(boost::polygon::VERTICAL).high();
+
+                double cx = 0.5 * (x1 + x2);
+                double cy = 0.5 * (y1 + y2);
+
+                double a1 = compare::OVERLAY_ELLIPSE_FACTOR1 * 0.5 * (x2 - x1);
+                double b1 = compare::OVERLAY_ELLIPSE_FACTOR1 * 0.5 * (y2 - y1);
+                double a2 = compare::OVERLAY_ELLIPSE_FACTOR2 * 0.5 * (x2 - x1);
+                double b2 = compare::OVERLAY_ELLIPSE_FACTOR2 * 0.5 * (y2 - y1);
+
+                double angle = 0.0;
+
+                std::vector<boost::polygon::point_data<double> > points;
+
+                for (int v = 0; v < compare::OVERLAY_ELLIPSE_COUNT - 1; v++, angle += compare::delta)
+                {
+                    if (v % 2 == 0)
+                    {
+                        points.push_back(boost::polygon::point_data<double>(cx + a1 * sin(angle), cy + b1 * cos(angle)));
+                    }
+                    else
+                    {
+                        points.push_back(boost::polygon::point_data<double>(cx + a2 * sin(angle), cy + b2 * cos(angle)));
+                    }
+                }
+                points.push_back(points.front());
+
+                POLY ellipse;
+                ellipse.set(points.begin(), points.end());
+                ellipses.push_back(ellipse);
+            }
+        }
+
+
+        // Create integer-versions of the polygons
+        // This is done BEFORE to not influence the performance test
+        // by converting inside measurements.
+        typedef boost::polygon::polygon_set_data<type_for_booleans> INT_POLY;
+        std::vector<INT_POLY> integer_polygons;
+        std::vector<INT_POLY> integer_ellipses;
+
+        for (std::vector<POLY>::const_iterator it = polygons.begin();
+            it != polygons.end();
+            ++it)
+        {
+            INT_POLY ip;
+            convert_polygon(*it, ip);
+            integer_polygons.push_back(ip);
+
+            // To check if conversion is OK
+            //std::cout << boost::polygon::area(*it) << " " <<  (boost::polygon::area(ip) / GTL_INTEGER_FACTOR_SQR) << std::endl;
+        }
+
+        for (std::vector<POLY>::const_iterator it = ellipses.begin();
+            it != ellipses.end();
+            ++it)
+        {
+            INT_POLY ip;
+            convert_polygon(*it, ip);
+            integer_ellipses.push_back(ip);
+        }
+
+
+        compare::wait();
+
+        if (compare::MEASURE_AREA)
+        {
+            double area = 0;
+            boost::timer t;
+            for (int i = 0; i < compare::AREA_COUNT; i++)
+            {
+                for (std::vector<POLY>::const_iterator it = polygons.begin();
+                    it != polygons.end();
+                    ++it)
+                {
+                    area += boost::polygon::area(*it);
+                }
+            }
+            compare::report_area(t, polygons.size(), area);
+        }
+
+        if (compare::MEASURE_OVERLAY)
+        {
+            double area1 = 0.0, area2 = 0.0;
+
+            boost::timer t;
+            for (int i = 0; i < compare::OVERLAY_COUNT; i++)
+            {
+                int k = 0;
+                std::vector<INT_POLY>::const_iterator eit = integer_ellipses.begin();
+                for (std::vector<INT_POLY>::const_iterator it = integer_polygons.begin();
+                    it != integer_polygons.end() && eit != integer_ellipses.end();
+                    ++it, ++eit, ++k)
+                {
+                    // * means intersection, + means union, ^ means XOR
+
+#if defined(OVERLAY_UNION)
+                    INT_POLY overlay = (*it) + (*eit);
+#elif defined(OVERLAY_SYMDIFF)
+                    INT_POLY overlay = (*it) ^ (*eit);
+#else
+                    //INT_POLY overlay = (*it) * (*eit); // does not compile anymore,
+                    // so se do it like this:
+                    INT_POLY overlay = boost::polygon::polygon_set_view<INT_POLY, INT_POLY, 1>(*it, *eit);
+#endif
+                    BOX b;
+                    boost::polygon::extents(b, overlay);
+                    if (compare::OVERLAY_AREA)
+                    {
+                        area1 += boost::polygon::area(*it) / GTL_INTEGER_FACTOR_SQR;
+                        //area2 += boost::polygon::area(overlay) / GTL_INTEGER_FACTOR_SQR;
+                        {
+                            typedef std::list<boost::polygon::polygon_with_holes_data<type_for_booleans> > list;
+                            double a = 0.0;
+                            list polys;
+                            overlay.get(polys);
+                            for (list::const_iterator it = polys.begin(); it != polys.end(); ++it)
+                            {
+                                a += boost::polygon::area(*it);
+                            }
+                            area2 += a / GTL_INTEGER_FACTOR_SQR;
+
+
+
+                            //std::cout << "insert into se_gtl values(" << k << ", " << (a / GTL_INTEGER_FACTOR_SQR) << ");" << std::endl;
+
+                            //    CREATING QUERY
+                            // select a.*,b.*,abs(a.area-b.area) as diff from se_ggl a left join se b on a.k=b.k order by diff desc
+                        }
+                    }
+                }
+            }
+            compare::report_overlay(t, polygons.size(), area1, area2);
+        }
+
+
+        if (compare::MEASURE_WITHIN)
+        {
+            int count = 0;
+            boost::timer t;
+            for (int e = 0; e < boxes.size(); e++)
+            {
+                const BOX& b = boxes[e];
+                POINT p;
+                boost::polygon::center(p, b);
+
+                //compare::debug_within(e, count);
+
+                std::vector<BOX>::const_iterator bit = boxes.begin();
+                for (std::vector<POLY>::const_iterator it = polygons.begin();
+                    it != polygons.end() && bit != boxes.end();
+                    ++it, ++bit)
+                {
+                    if (p.x() > bit->get(boost::polygon::HORIZONTAL).low()
+                        && p.x() < bit->get(boost::polygon::HORIZONTAL).high()
+                        && p.y() > bit->get(boost::polygon::VERTICAL).low()
+                        && p.y() < bit->get(boost::polygon::VERTICAL).high())
+                    {
+                        if (boost::polygon::contains(*it, p))
+                        {
+                            count++;
+                        }
+                    }
+                }
+            }
+            compare::report_within(t, polygons.size(), count, -1);
+        }
+    }
+    catch(std::exception const& e)
+    {
+        std::cout << e.what() << std::endl;
+    }
+    catch(...)
+    {
+        std::cout << "Exception encountered" << std::endl;
+    }
+
+    return 0;
+}