[Boost-commit] svn:boost r55832 - sandbox/ggl/other/comparisons/geos

Subject: [Boost-commit] svn:boost r55832 - sandbox/ggl/other/comparisons/geos
From: barend.gehrels_at_[hidden]
Date: 2009-08-28 04:56:51

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

Log:
Added geometry library performance comparison CPP files
Added:
   sandbox/ggl/other/comparisons/geos/geos_check.cpp   (contents, props changed)

Added: sandbox/ggl/other/comparisons/geos/geos_check.cpp
==============================================================================
--- (empty file)
+++ sandbox/ggl/other/comparisons/geos/geos_check.cpp	2009-08-28 04:56:50 EDT (Fri, 28 Aug 2009)
@@ -0,0 +1,431 @@
+// Copyright Barend Gehrels, 1995-2009, Geodan Holding B.V. Amsterdam, the Netherlands.
+// Copyright Bruno Lalande 2008
+// 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 <iostream>
+#include <iterator>
+#include <vector>
+#include <algorithm>
+#include <vector>
+#include <sstream>
+#include <iomanip>
+#include <geos/geom/PrecisionModel.h>
+#include <geos/geom/GeometryFactory.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/Point.h>
+#include <geos/geom/LinearRing.h>
+#include <geos/geom/LineString.h>
+#include <geos/geom/Polygon.h>
+#include <geos/geom/GeometryCollection.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/CoordinateSequence.h>
+#include <geos/geom/CoordinateArraySequence.h>
+#include <geos/geom/IntersectionMatrix.h>
+
+#include <geos/util/GeometricShapeFactory.h>
+#include <geos/util/GEOSException.h>
+#include <geos/util/IllegalArgumentException.h>
+#include <geos/opLinemerge.h>
+#include <geos/opPolygonize.h>
+
+#include <geos/simplify/DouglasPeuckerSimplifier.h>
+
+#include <geos/io/WKTWriter.h>
+
+
+
+#include "../common/common.hpp"
+#include "../common/read_shapefile.hpp"
+
+using namespace geos;
+using namespace geos::geom;
+
+GeometryFactory *global_factory;
+
+
+template <typename T>
+T convert(SHPObject* psShape)
+{
+    const double* x = psShape->padfX;
+    const double* y = psShape->padfY;
+
+    // We will use a coordinate list to build the linearring
+    CoordinateSequence *cl = new CoordinateArraySequence();
+
+    for (int v = 0; v < psShape->nVertices; v++)
+    {
+        cl->add(Coordinate(x[v], y[v]));
+    }
+    LinearRing *lr = global_factory->createLinearRing(cl);
+    return global_factory->createPolygon(lr, NULL);
+}
+
+
+
+
+
+
+void write_envelope(const Envelope* e)
+{
+    std::cout << e->getMinX() << "," << e->getMinY()
+        << "," << e->getMaxX() << "," << e->getMaxY() << std::endl;
+}
+
+int main(int argc, char** argv)
+{
+    try
+    {
+        compare::version_info("geos");
+
+        PrecisionModel *pm = new PrecisionModel(geos::geom::PrecisionModel::FLOATING);
+        global_factory = new GeometryFactory(pm, -1);
+        delete pm;
+
+        geos::io::WKTWriter writer; // for checking things
+
+
+        std::vector<Geometry*> polygons;
+        std::vector<Geometry*> ellipses;
+        std::vector<int> ids;
+        read_shapefile(compare::filename(argc, argv), compare::fieldname(), polygons, ids, convert<Geometry*>);
+
+
+        // Create envelopes
+        std::vector<const Envelope*> envelopes;
+        std::vector<Geometry*> boxes; // decreased, for intersection
+        {
+            boost::timer t;
+            for (std::vector<Geometry*>::const_iterator it = polygons.begin();
+                it != polygons.end(); ++it)
+            {
+                Geometry* poly = *it;
+                const Envelope* e = poly->getEnvelopeInternal();
+
+                Envelope* b = new Envelope(*e);
+                envelopes.push_back(e);
+
+                b->expandBy(- compare::CLIP_DISTANCE);
+
+                Geometry* g = global_factory->toGeometry(b);
+                boxes.push_back(g);
+            }
+            //std::cout << "ENVELOPE " << polygons.size() << " in " << t.elapsed() << "s" << std::endl;
+        }
+
+        // Create the ellipses for intersections lateron
+        if (compare::MEASURE_OVERLAY)
+        {
+            for (std::vector<const Envelope*>::const_iterator eit = envelopes.begin();
+                    eit != envelopes.end();
+                    ++eit)
+            {
+                double factor = 0.9;
+                double cx = 0.5 * ((*eit)->getMinX() + (*eit)->getMaxX());
+                double cy = 0.5 * ((*eit)->getMinY() + (*eit)->getMaxY());
+
+                double dx = (*eit)->getMaxX() - (*eit)->getMinX();
+                double dy = (*eit)->getMaxY() - (*eit)->getMinY();
+
+                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;
+
+                CoordinateSequence *cl = new CoordinateArraySequence();
+
+                for (int v = 0; v < compare::OVERLAY_ELLIPSE_COUNT - 1; v++, angle += compare::delta)
+                {
+                    if (v % 2 == 0)
+                    {
+                        cl->add(Coordinate(cx + a1 * sin(angle), cy + b1 * cos(angle)));
+                    }
+                    else
+                    {
+                        cl->add(Coordinate(cx + a2 * sin(angle), cy + b2 * cos(angle)));
+                    }
+                }
+                // Close it (when we use the calculation we get a "non closed" exception)
+                cl->add(cl->front());
+
+                LinearRing *lr = global_factory->createLinearRing(cl);
+                ellipses.push_back(global_factory->createPolygon(lr, NULL));
+            }
+        }
+
+
+        compare::wait();
+
+
+        if (compare::MEASURE_AREA)
+        {
+            double area = 0;
+            boost::timer t;
+            for (int i = 0; i < compare::AREA_COUNT; i++)
+            {
+                for (std::vector<Geometry*>::const_iterator it = polygons.begin();
+                    it != polygons.end(); ++it)
+                {
+                    Geometry* poly = *it;
+                    area += poly->getArea();
+                }
+            }
+            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<Geometry*>::const_iterator it = polygons.begin();
+                    it != polygons.end();
+                    ++it)
+                {
+                    Geometry* centroid = (*it)->getCentroid();
+                    Point* c = (Point*) centroid;
+                    sum_x += c->getX();
+                    sum_y += c->getY();
+                    global_factory->destroyGeometry(centroid);
+                }
+            }
+            compare::report_centroid(t, polygons.size(), sum_x, sum_y);
+        }
+
+
+        if (compare::MEASURE_CONVEX_HULL)
+        {
+            double area = 0;
+            boost::timer t;
+            for (std::vector<Geometry*>::const_iterator it = polygons.begin();
+                it != polygons.end(); ++it)
+            {
+                Geometry* hull = (*it)->convexHull();
+                area += hull->getArea();
+                global_factory->destroyGeometry(hull);
+            }
+            compare::report_hull(t, polygons.size(), area);
+        }
+
+
+        if (compare::MEASURE_OVERLAY)
+        {
+            bool first = true;
+            bool first_exception = true;
+            int exception_count = 0;
+            double area1 = 0, area2 = 0;
+            std::cout << std::setprecision(16);
+            std::string exceptions;
+
+            boost::timer t;
+            for (int i = 0; i < compare::OVERLAY_COUNT; i++)
+            {
+                int k = 0;
+                std::vector<Geometry*>::const_iterator eit = ellipses.begin();
+                for (std::vector<Geometry*>::const_iterator it = polygons.begin();
+                        it != polygons.end() && eit != ellipses.end();
+                        ++it, ++eit, ++k)
+                {
+                    //if (k == 768)
+                    {
+                        bool has_exception = true;
+                        Geometry* poly = *it;
+                        Geometry* ellipse = *eit;
+
+                        /***
+                        POLYGON 768 has a problem when factor=1 (assertion failed), I'll report this.
+                        std::cout << k << ": select ST_Area(intersection"
+                            << "(GeometryFromText('"
+                            <<  writer.write(poly)
+                            << "'),GeometryFromText('" << writer.write(ellipse)
+                            << "')));"
+                            << std::endl;
+                        ***/
+
+                        Geometry *g = NULL;
+                        try
+                        {
+                            g = poly->intersection(ellipse);
+
+                            if (compare::OVERLAY_AREA)
+                            {
+                                area1 += poly->getArea();
+                                area2 += g->getArea();
+                                //std::cout << k << ": " << g->getArea() << std::endl;
+                            }
+
+                            has_exception = false;
+                        }
+                        catch(const std::exception& e)
+                        {
+                            exceptions += e.what();
+                            exceptions += "\n";
+                            exception_count++;
+
+                            if (first_exception)
+                            {
+                                std::cout << e.what() << std::endl;
+
+
+                                first_exception = false;
+                            }
+                        }
+                        catch(...)
+                        {
+                            exception_count++;
+                        }
+                        /*if (first && ! has_exception)
+                        {
+                            first = false;
+                            std::cout << "Intersected: " << k << " n0: " << poly->getNumPoints()
+                                << " -> " << g->getNumGeometries() << "("
+                                << g->getNumPoints()
+                                << ")" << std::endl;
+                            first = false;
+                        }*/
+
+                        global_factory->destroyGeometry(g);
+                    }
+                }
+            }
+            compare::report_overlay(t, polygons.size(), area1, area2);
+            if (exception_count > 0)
+            {
+                std::cout << exception_count << " exceptions:" << std::endl;
+                    //<< exceptions << std::endl;
+            }
+        }
+
+
+        if (compare::MEASURE_SIMPLIFY)
+        {
+            int count1 = 0, count2 = 0;
+            double length1 = 0.0, length2 = 0.0;
+            boost::timer t;
+            for (std::vector<Geometry*>::const_iterator it = polygons.begin();
+                it != polygons.end();
+                ++it)
+            {
+                Geometry* poly = *it;
+                count1 += poly->getNumPoints();
+                std::auto_ptr<Geometry> simplified = geos::simplify::DouglasPeuckerSimplifier::simplify(poly, compare::SIMPLIFY_DISTANCE);
+                count2 += simplified->getNumPoints();
+
+                if (compare::SIMPLIFY_LENGTH)
+                {
+                    length1 += ((Polygon *) poly)->getLength();
+                    length2 += ((Polygon *) simplified.get())->getLength();
+                }
+
+            }
+            compare::report_simplify(t, polygons.size(), length1, length2, count1, count2);
+        }
+
+        if (compare::MEASURE_TOUCH)
+        {
+            int count = 0;
+            int count_non_disjoint = 0;
+            boost::timer t;
+            int e1 = 0;
+            for (std::vector<Geometry*>::const_iterator it1 = polygons.begin();
+                it1 != polygons.end(); ++it1, e1++)
+            {
+                int e2 = 0;
+                for (std::vector<Geometry*>::const_iterator it2 = polygons.begin();
+                    it2 != polygons.end(); ++it2, e2++)
+                {
+                    Envelope env1 = *envelopes[e1];
+                    Envelope env2 = *envelopes[e2];
+                    env1.expandBy(2.0);
+                    env2.expandBy(2.0);
+                    if (env1.intersects(env2))
+                    {
+                        count_non_disjoint++;
+                        //std::string r = (*it1)->relate(*it2)->toString();
+                        if (! (*it1)->disjoint(*it2))
+                        {
+                            /*
+                            std::cout << ids[e1] << " - " << ids[e2]
+                                //<< " " << r
+                                << std::endl;
+                            */
+                            count++;
+                        }
+                    }
+                }
+            }
+            compare::report_touch(t, polygons.size(), count, count_non_disjoint);
+        }
+
+
+        if (compare::MEASURE_WITHIN)
+        {
+            int count = 0;
+            boost::timer t;
+            for (int e = 0; e < envelopes.size(); e++)
+            {
+                Coordinate c;
+                envelopes[e]->centre(c);
+
+                Point *p = global_factory->createPoint(c);
+
+                int k = 0;
+                std::vector<const Envelope*>::const_iterator bit = envelopes.begin();
+                for (std::vector<Geometry*>::const_iterator it = polygons.begin();
+                    it != polygons.end() && bit != envelopes.end();
+                    ++it, ++bit, k++)
+                {
+                    /***
+                    if ((e == 85 && k == 85)
+                        || e == 99 && k == 99
+                        )
+                    {
+                        std::cout << e << "," << k << std::endl;
+                        write_envelope(envelopes[e]);
+                        std::cout << p->getX() << "," << p->getY() << std::endl;
+                    }
+                    ***/
+
+                    // First check bbox (like all other tests) (though GEOS does this also internally)
+                    // (perf.difference: ca 1%)
+                    if ((*bit)->contains(c))
+                    {
+                        Geometry* poly = *it;
+                        if (p->within(poly))
+                        {
+                            //std::cout << e << " IN " << k << std::endl;
+                            count++;
+                        }
+                    }
+                }
+                global_factory->destroyGeometry(p);
+            }
+            compare::report_within(t, polygons.size(), count, -1);
+        }
+
+
+        for (int i = 0; i < polygons.size(); i++)
+        {
+            global_factory->destroyGeometry(polygons[i]);
+        }
+        for (int i = 0; i < ellipses.size(); i++)
+        {
+            global_factory->destroyGeometry(ellipses[i]);
+        }
+        for (int i = 0; i < boxes.size(); i++)
+        {
+            global_factory->destroyGeometry(boxes[i]);
+        }
+    }
+    catch(std::exception const& e)
+    {
+        std::cout << e.what() << std::endl;
+    }
+
+    return 0;
+}
+