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;
+}