[Boost-commit] svn:boost r55833 - in sandbox/ggl/other/comparisons: ggl terralib wykobi

Subject: [Boost-commit] svn:boost r55833 - in sandbox/ggl/other/comparisons: ggl terralib wykobi
From: barend.gehrels_at_[hidden]
Date: 2009-08-28 04:58:00

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

Log:
Added geometry library performance comparison CPP files
Added:
   sandbox/ggl/other/comparisons/ggl/ggl_simplify.cpp   (contents, props changed)
   sandbox/ggl/other/comparisons/terralib/terralib_check.cpp   (contents, props changed)
   sandbox/ggl/other/comparisons/wykobi/wykobi_check.cpp   (contents, props changed)

Added: sandbox/ggl/other/comparisons/ggl/ggl_simplify.cpp
==============================================================================
--- (empty file)
+++ sandbox/ggl/other/comparisons/ggl/ggl_simplify.cpp	2009-08-28 04:57:59 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,88 @@
+
+#include "../common/common.hpp"
+#include "../common/read_shapefile.hpp"
+
+#include <boost/timer.hpp>
+#include <ggl/ggl.hpp>
+#include <ggl/algorithms/perimeter.hpp>
+#include <ggl/geometries/geometries.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 T::point_type POINT;
+        //POINT p;
+        ggl::point_xy<double> p;
+        ggl::set<0>(p, x[v]);
+        ggl::set<1>(p, y[v]);
+        polygon.outer().push_back(p);
+    }
+    return polygon;
+}
+
+
+
+
+
+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> boxes10; // for intersections
+    std::vector<int> ids;
+    read_shapefile(compare::filename(argc, argv), compare::fieldname(), polygons, ids, convert<POLY>);
+
+    double d = 0.01;
+    bool iterate = true;
+    while(iterate)
+    {
+        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, d);
+            count1 += it->outer().size();
+            count2 += p.outer().size();
+            length1 += ggl::perimeter(it->outer());
+            length2 += ggl::perimeter(p.outer());
+        }
+        length1 /= 1000.0;
+        length2 /= 1000.0;
+        double diff = std::abs(length1 - length2);
+        double frac = length2 / length1;
+
+        std::cout << "SIMPLIFY " << polygons.size() << std::endl
+            << " distance: " << d << std::endl
+            << " frac: " << frac << std::endl
+            << " vertices: " << count1 << " -> " << count2 << std::endl
+            << " perimeter: " << length1 << " -> " << length2 << std::endl
+            << " time: "<< t.elapsed() << "s" << std::endl;
+
+        if (frac < 0.99 && d > 0.0001)
+        {
+            d *= 0.95;
+        }
+        else
+        {
+            iterate = false;
+        }
+
+    }
+
+
+
+    return 0;
+}

Added: sandbox/ggl/other/comparisons/terralib/terralib_check.cpp
==============================================================================
--- (empty file)
+++ sandbox/ggl/other/comparisons/terralib/terralib_check.cpp	2009-08-28 04:57:59 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,272 @@
+#include <TeVersion.h>
+#include <TeDatabase.h>
+#include <TeGeometryAlgorithms.h>
+#include <TeOverlay.h>
+
+#include "TeException.h"
+
+#include "../common/common.hpp"
+#include "../common/read_shapefile.hpp"
+
+
+template <typename T>
+T convert(SHPObject* psShape)
+{
+    const double* x = psShape->padfX;
+    const double* y = psShape->padfY;
+
+    TeLinearRing ring;
+    for (int v = 0; v < psShape->nVertices; v++)
+    {
+        TeCoord2D pt (x[v], y[v]);
+        ring.add (pt);
+    }
+    TePolygon p;
+    p.add(ring);
+    return p;
+}
+
+
+template <typename T>
+T convert_line(SHPObject* psShape)
+{
+    const double* x = psShape->padfX;
+    const double* y = psShape->padfY;
+
+    TeLine2D line;
+    for (int v = 0; v < psShape->nVertices; v++)
+    {
+        TeCoord2D pt (x[v], y[v]);
+        line.add (pt);
+    }
+    return line;
+}
+
+
+int main(int argc, char** argv)
+{
+    compare::version_info("terralib");
+
+    std::vector<TeLine2D> lines;
+    std::vector<TePolygon> polygons;
+    std::vector<TeBox> boxes;
+    std::vector<int> ids;
+
+    std::vector<TePolygonSet> ellipses;
+    std::vector<TePolygonSet> polygonsets;
+
+    read_shapefile(compare::filename(argc, argv), compare::fieldname(), polygons, ids, convert<TePolygon>);
+
+
+    // Read also polygons as (for simplify) (note that thse could be extracted from polygon rings
+    // However, performance is not measured so we take it easy here.
+    read_shapefile(compare::filename(argc, argv), compare::fieldname(), lines, ids, convert_line<TeLine2D>);
+
+    // Create envelops
+    for (std::vector<TePolygon>::const_iterator it = polygons.begin();
+            it != polygons.end(); it++)
+    {
+        TeBox box = it->box();
+        boxes.push_back(box);
+    }
+
+    // Create polygon sets
+    for (std::vector<TePolygon>::const_iterator it = polygons.begin();
+            it != polygons.end(); it++)
+    {
+        TePolygonSet set;
+        set.add(*it);
+        polygonsets.push_back(set);
+    }
+
+    // Create the ellipses for intersections lateron
+    if (compare::MEASURE_OVERLAY)
+    {
+        for (std::vector<TeBox>::const_iterator bit = boxes.begin();
+                bit != boxes.end();
+                ++bit)
+        {
+            double factor = 0.9;
+            double cx = 0.5 * (bit->x1() + bit->x2());
+            double cy = 0.5 * (bit->y1() + bit->y2());
+
+            double dx = bit->x2() - bit->x1();
+            double dy = bit->y2() - bit->y1();
+
+            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;
+
+            TeLinearRing ring;
+
+            for (int v = 0; v < compare::OVERLAY_ELLIPSE_COUNT - 1; v++, angle += compare::delta)
+            {
+                if (v % 2 == 0)
+                {
+                    ring.add(TeCoord2D(cx + a1 * sin(angle), cy + b1 * cos(angle)));
+                }
+                else
+                {
+                    ring.add(TeCoord2D(cx + a2 * sin(angle), cy + b2 * cos(angle)));
+                }
+            }
+            ring.add(ring.first());
+
+            TePolygon poly;
+            TePolygonSet ellipse;
+            poly.add(ring);
+            ellipse.add(poly);
+            ellipses.push_back(ellipse);
+        }
+    }
+
+    compare::wait();
+
+
+    if (compare::MEASURE_AREA)
+    {
+        double area = 0;
+        boost::timer t;
+        for (int i = 0; i < compare::AREA_COUNT; i++)
+        {
+            for (std::vector<TePolygon>::const_iterator it = polygons.begin();
+                    it != polygons.end(); ++it)
+            {
+                area += TeGeometryArea(*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<TePolygon>::const_iterator it = polygons.begin();
+                it != polygons.end();
+                ++it)
+            {
+                TeCoord2D centroid = TeFindCentroid(*it);
+                sum_x += centroid.x();
+                sum_y += centroid.y();
+            }
+        }
+        compare::report_centroid(t, polygons.size(), sum_x, sum_y);
+    }
+
+
+    if (compare::MEASURE_CONVEX_HULL)
+    {
+        double area = 0;
+        boost::timer t;
+        for (std::vector<TePolygon>::const_iterator it = polygons.begin();
+                it != polygons.end(); ++it)
+        {
+            TePolygon hull = TeConvexHull(*it);
+            if (compare::HULL_AREA)
+            {
+                area += TeGeometryArea(hull);
+            }
+        }
+        compare::report_hull(t, polygons.size(), area);
+    }
+
+
+    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<TePolygonSet>::const_iterator eit = ellipses.begin();
+            for (std::vector<TePolygonSet>::const_iterator it = polygonsets.begin();
+                it != polygonsets.end() && eit != ellipses.end();
+                ++it, ++eit, ++k)
+            {
+                if (compare::OVERLAY_AREA)
+                {
+                    area1 += TeGeometryArea(*it);
+                }
+
+                TePolygonSet v;
+                TeOVERLAY::TeIntersection(*it, *eit, v);
+                if (compare::OVERLAY_AREA)
+                {
+                    double a = 0.0;
+                    for (int i = 0; i < v.size(); i++)
+                    {
+                        a += TeGeometryArea(v[i]);
+                    }
+                    area2 += a;
+                }
+            }
+        }
+        compare::report_overlay(t, polygons.size(), area1, area2);
+    }
+
+
+    if (compare::MEASURE_SIMPLIFY)
+    {
+        // Modifies the line! But OK, rest is done on polygon
+        int count1 = 0, count2 = 0;
+        double length1 = 0.0, length2 = 0.0;
+        boost::timer t;
+        for (std::vector<TeLine2D>::iterator it = lines.begin(); it != lines.end(); ++it)
+        {
+            count1 += it->size();
+            if (compare::SIMPLIFY_LENGTH)
+            {
+                length1 += TeLength(*it);
+            }
+
+            TeLineSimplify(*it, compare::SIMPLIFY_DISTANCE, compare::SIMPLIFY_DISTANCE * 100);
+
+            count2 += it->size();
+            if (compare::SIMPLIFY_LENGTH)
+            {
+                length2 += TeLength(*it);
+            }
+        }
+        compare::report_simplify(t, polygons.size(), length1, length2, count1, count2);
+    }
+
+
+    if (compare::MEASURE_WITHIN)
+    {
+        int count = 0;
+        boost::timer t;
+        for (unsigned int e = 0; e < boxes.size(); e++)
+        {
+            TeCoord2D c = boxes[e].center();
+            TePoint p(c);
+
+            std::vector<TeBox>::const_iterator bit = boxes.begin();
+            int k = 0;
+            for (std::vector<TePolygon>::const_iterator it = polygons.begin();
+                it != polygons.end() && bit != boxes.end();
+                ++it, ++bit, ++k)
+            {
+                if (TeWithin(c, *bit) && TeWithin(c, *it))
+                {
+                    //std::cout << e << " IN " << k << std::endl;
+                    count++;
+                }
+            }
+        }
+        compare::report_within(t, polygons.size(), count, -1);
+    }
+
+    return 0;
+}
+

Added: sandbox/ggl/other/comparisons/wykobi/wykobi_check.cpp
==============================================================================
--- (empty file)
+++ sandbox/ggl/other/comparisons/wykobi/wykobi_check.cpp	2009-08-28 04:57:59 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,211 @@
+#include <iostream>
+#include <string>
+
+#include "../common/common.hpp"
+#include "../common/read_shapefile.hpp"
+
+#include <boost/timer.hpp>
+
+#include <wykobi.hpp>
+#include <wykobi_algorithm.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++)
+    {
+        typedef typename T::PointType POINT;
+        POINT p;
+        p[0] = x[v];
+        p[1] = y[v];
+        polygon.push_back(p);
+    }
+    return polygon;
+}
+
+
+
+int main(int argc, char** argv)
+{
+    compare::version_info("wykobi");
+
+    typedef wykobi::point2d<double> POINT;
+    typedef wykobi::polygon<double, 2> POLY;
+    typedef wykobi::rectangle<double> BOX;
+
+    std::vector<POLY> polygons;
+    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 b = wykobi::aabb(*it);
+        boxes.push_back(b);
+    }
+
+
+    // Create the smaller boxes for clips lateron
+    if (compare::MEASURE_CLIP)
+    {
+
+        int k = 0;
+        for (std::vector<BOX>::const_iterator bit = boxes.begin();
+            bit != boxes.end();
+            ++bit, ++k)
+        {
+            BOX const& bx = *bit;
+            POINT c;
+            c.x = (bx[0].x + bx[1].x) / 2.0;
+            c.y = (bx[0].y + bx[1].y) / 2.0;
+
+
+            double a = compare::CLIP_FACTOR * 0.5 * (bx[1].x - bx[0].x);
+            double b = compare::CLIP_FACTOR * 0.5 * (bx[1].y - bx[0].y);
+
+            BOX box;
+
+            double angle = 225.0 * wykobi::PIDiv180;
+            box[0].x = c.x + a * sin(angle);
+            box[0].y = c.y + b * cos(angle);
+
+            angle = 45.0 * wykobi::PIDiv180 ;
+            box[1].x = c.x + a * sin(angle);
+            box[1].y = c.y + b * cos(angle);
+
+            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 += wykobi::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 = wykobi::centroid(*it);
+                sum_x += centroid.x;
+                sum_y += centroid.y;
+            }
+        }
+        compare::report_centroid(t, polygons.size(), sum_x, sum_y);
+    }
+
+    if (compare::MEASURE_CLIP)
+    {
+        bool first = true;
+        double area1 = 0.0, area2 = 0.0;
+
+        std::cout << std::setprecision(16);
+
+        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 += wykobi::area(*it);
+                }
+
+
+                POLY out;
+                wykobi::algorithm::sutherland_hodgman_polygon_clipper<POINT>
+                    clipper(*bit, *it, out);
+                if (compare::CLIP_AREA)
+                {
+                    area2 += wykobi::area(out);
+                }
+            }
+        }
+        compare::report_clip(t, polygons.size(), area1, area2);
+    }
+
+
+    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 out;
+            wykobi::algorithm::convex_hull_graham_scan<POINT> graham(it->begin(), it->end(), std::back_inserter(out));
+            if (compare::HULL_AREA)
+            {
+                area += wykobi::area(out);
+            }
+        }
+        compare::report_hull(t, polygons.size(), area);
+    }
+
+    if (compare::MEASURE_WITHIN)
+    {
+        int count = 0;
+        boost::timer t;
+        for (int e = 0; e < boxes.size(); e++)
+        {
+            BOX const& b = boxes[e];
+            POINT p;
+            p.x = (b[0].x + b[1].x) / 2.0;
+            p.y = (b[0].y + b[1].y) / 2.0;
+
+            //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 (wykobi::point_in_rectangle(p, *bit)
+                    && wykobi::point_in_polygon(p, *it))
+                    //&& wykobi::point_in_polygon_winding_number(p, *it))
+                {
+                    count++;
+                }
+            }
+        }
+        compare::report_within(t, polygons.size(), count, -1);
+    }
+
+
+    return 0;
+}