Geometry Boost List: [ggl] Point/Box in Box Test (Border Cases - 3D)

Subject: [ggl] Point/Box in Box Test (Border Cases - 3D)
From: Frank Glinka (glinkaf)
Date: 2011-07-07 09:47:03

Hi Barend!

Am 03.07.2011 16:36, schrieb Barend Gehrels:
> Hi Frank,
>
> I updated the trunk, and moved the within_code from within_util to
> within.hpp.

Thank you Barend! I used it and it works as described/expected. A thing
I should notice is that the performance of the original 'within'
implementation is negatively impacted by this change.

I attached my basic ad-hoc micro-benchmark if you are interested in
doing your own shot. Results with VS2010 in x86 Release and /O2 (run on
Core 2 Duo 2.4 GHz):

Run with old within: 5.72521s
Run with updated within: 23.2128s
Run with within_code: 25.2096s
Run with intersects: 6.86401s
Run with custom coveredBy: 4.29001s (based on outside (>) check solely)

As the difference is quite significant, it might be a thought to
consider providing specialized implementations of 'within' and
'coveredBy' instead of using the 'outside, onBorder, inside' predicate
as an abstraction in both versions.

> It stays in namespace detail until there is complete agreement on the
> interface (parameter, strategy, function name), so maybe you will have
> to rename something in the future, but the functionality is there. I now
> believe that it is the best to adapt to PostGIS and Oracle here, they
> call it "coveredBy". See
> http://postgis.refractions.net/documentation/manual-1.5/ST_CoveredBy.html,
> defining it as "Returns 1 (TRUE) if no point in Geometry A is outside
> GeometryB", which is what I described. So covered_by would be a
> convenient and (I think) intuitive name.

I like that name and staying compliant to established terms is always
good. I guess you noticed that this specification would allow for a
simpler implementation than currently provided by within_code!? (Which
is able to distinguish between 'inside' and 'on border' while the
specification only requires 'on border' - which probably causes the
previously mentioned performance overhead.)

Best regards,
Frank
-------------- next part --------------
// GeometryTest.cpp : Micro-benchmark for certain boost.geometry functions.
// Author: Frank Glinka, University of Muenster, Germany

#include <vector>
#include <iostream>

#include <boost\thread\xtime.hpp>
#include <boost/geometry/algorithms/within.hpp>
#include <boost/geometry/strategies/strategies.hpp>  // It seems this header must be included, otherwise the register 
                                                     // macros fail to do their job.
#include <boost/geometry/algorithms/intersects.hpp>
#include <boost/geometry/geometries/register/point.hpp>
#include <boost/geometry/geometries/register/box.hpp>

// Create custom vector type.
struct Point3D {
  float x, y, z;

  Point3D() : x(0.0f), y(0.0f), z(0.0f) {}

  Point3D(float a, float b, float c) : x(a), y(b), z(c) {}
};

// Create custom space type.
class Box3D {
public:
  Box3D() : origin(Point3D(0.0f, 0.0f, 0.0f)), size(Point3D(0.0f, 0.0f, 0.0f)) {
  }

  Box3D(Point3D a, Point3D b) : origin(a), size(b) {
  }

  Point3D origin, size;
};

// Register custom types.
BOOST_GEOMETRY_REGISTER_POINT_3D(Point3D, float, boost::geometry::cs::cartesian, x, y, z)
BOOST_GEOMETRY_REGISTER_BOX(Box3D, Point3D, origin, size);

// Helper class for time measurements.
class Timer {
private:
  boost::xtime begin, end;

public:
  Timer() {
    start();
  }

  inline void start() {
    boost::xtime_get(&begin, boost::TIME_UTC);
    end = begin;
  }

  inline double stop() {
    boost::xtime_get(&end, boost::TIME_UTC);
    end.sec -= begin.sec;
    end.nsec -= begin.nsec;
    if(end.nsec < 0) {
      --end.sec;
      end.nsec += 1000000000;
    }

    return duration();
  }

  inline double duration() {
    return ((double)end.sec) + ((double) end.nsec / 1000000000);
  }
};

int main(int argc, char *argv[], char *envp[]) {
  // Experiment settings.
  const size_t batch_size = 1000000;
  const size_t iterations = 1000;
  
  // Create data structures (set of points).
  const Box3D box(Point3D(0,0,0), Point3D(batch_size/2, batch_size/2, batch_size/2));
  const Box3D boxTwo(Point3D(1,1,1), Point3D(2, 2, 2));
  std::vector<Point3D> points;
  points.reserve(batch_size);
  for(size_t i = 0; i < batch_size; ++i) {
    points.push_back(Point3D((float)i, (float)i, (float)i));
  }

  for(int i = 0; i < 10; i++) {
    // Point within box performance test...
    int containsCounter = 0;
    Timer timer = Timer();
    for(size_t j = 0; j < iterations; j++) {
      for(size_t i = 0; i < batch_size; i++) {
        if(boost::geometry::within(points[i], box)) {
          ++containsCounter;
        }
      }
    }
    std::cout << "Point within box testrun took " << timer.stop() << " seconds. Contains: " << containsCounter << std::endl;
  }

  {
    // Point within_code box performance test...
    int containsCounter = 0;
    Timer timer = Timer();
    for(size_t j = 0; j < iterations; j++) {
      for(size_t i = 0; i < batch_size; i++) {
        if(boost::geometry::detail::within_code(points[i], box) != -1) {
          ++containsCounter;
        }
      }
    }
    std::cout << "Point within_code box testrun took " << timer.stop() << " seconds. Contains: " << containsCounter << std::endl;
  }

  {
    // Point intersection box performance test...
    int containsCounter = 0;
    Timer timer = Timer();
    for(size_t j = 0; j < iterations; j++) {
      for(size_t i = 0; i < batch_size; i++) {
        if(boost::geometry::intersects(points[i], box)) {
          ++containsCounter;
        }
      }
    }
    std::cout << "Point intersects box testrun took " << timer.stop() << " seconds. Contains: " << containsCounter << std::endl;

  }

  {
    // Manual point coveredBy box performance test (<)...
    int containsCounter = 0;
    Timer timer = Timer();
    for(size_t j = 0; j < iterations; j++) {
      for(size_t i = 0; i < batch_size; i++) {
        if(points[i].x < box.origin.x || points[i].x > box.size.x || points[i].y < box.origin.y || points[i].y > box.size.y || points[i].z < box.origin.z || points[i].z > box.size.z) {
          continue;
        } else {
          ++containsCounter;
        }
      }
    }
    std::cout << "Manual point coveredBy box testrun (<) took " << timer.stop() << " seconds. Contains: " << containsCounter << std::endl;
  }
   
  //{
  //  // Manual point coveredBy box performance test (>=)...
  //  int containsCounter = 0;
  //  Timer timer = Timer();
  //  for(size_t j = 0; j < iterations; j++) {
  //    for(size_t i = 0; i < batch_size; i++) {
  //      if(points[i].x >= box.origin.x && points[i].x <= box.size.x && points[i].y >= box.origin.y && points[i].y <= box.size.y && points[i].z >= box.origin.z && points[i].z <= box.size.z) {
  //        ++containsCounter;
  //      }
  //    }
  //  }
  //  std::cout << "Manual point coveredBy box performance testrun (>=) took " << timer.stop() << " seconds. Contains: " << containsCounter << std::endl;
  //}

  // Clean up data structures (set of points).
  points.clear();
  std::vector<Point3D>(points).swap(points);

  return 0;
}