Re: [boost] Geometry and spatial indexes, my opinion

Subject: Re: [boost] Geometry and spatial indexes, my opinion
From: Brandon Kohn (blkohn_at_[hidden])
Date: 2008-10-09 10:25:14

--------------------------------------------------
From: "Patrick Mihelich" <patrick.mihelich_at_[hidden]>
Sent: Thursday, October 09, 2008 6:23 AM
To: <boost_at_[hidden]>
Subject: Re: [boost] Geometry and spatial indexes, my opinion

> On Thu, Oct 9, 2008 at 1:48 AM, Bruno Lalande 
> <bruno.lalande_at_[hidden]>wrote:
>
>> Stop me if I'm wrong but what you're discussing now is how to
>> implement the point class provided by an ideal geometry library,
>> right? If so, I think it's not the most important concern since, as
>> I've already said, the primary goal is not to provide a point class
>> but to define how points are accessed and defined in general. I think
>> that, by talking about tuples, Brandon was rather pointing out that
>> points should be accessible the same way as tuples: get<N>(p). Which
>> is quite a general consensus now, IMO.
>
>
> OK, perhaps I misinterpreted Brandon's comment. I have been getting caught
> up on the March discussions, so apologies if I rehash old debates. If he 
> was
> simply advocating compile-time indexing via get<N>(p) or similar syntax,
> then I completely agree. This should be required by the Point concept (or
> Coordinate concept? I do not quite understand the distinction).
>
> On the other hand, I'm worried to see that the geometry::point
> implementation (at least the version included with SpatialIndexes) 
> supports
> only compile-time indexing. I think we must also have a RuntimeIndexable
> concept, to which a particular point type may or may not conform.
> RuntimeIndexable would require an indexing operator[], so it can only be
> modeled by homogeneous points. Iterators via begin() and end() could be
> useful as well. I have a couple of reasons for wanting this, as usual
> motivated by high dimensional spaces...
>

Actually in the library I've been working on the definition of the point 
type is of lesser importance. All of my algorithms are designed to access 
point coordinates via a traits specialization which details how the point's 
coordinates are accessed. It works like this (please excuse email formatting 
:)

I saw someone mention a fusion vector so I'll use that.

First I have a macro which makes for an easy way to define the basic traits 
of the point:

//! Macro for point type which does not have embedded traits - User Defined 
Points
#define BOOST_DEFINE_USER_POINT_TRAITS( Point, Coordinate, Dimension ) \
template <> \
struct point_traits< Point > \
{ \
   BOOST_STATIC_ASSERT( Dimension > 0 ); \
   typedef Point point_type; \
   typedef dimension_traits<Dimension> dimension_type; \
   typedef coordinate_traits<Coordinate>::coordinate_type coordinate_type;\
};

So we need:

BOOST_DEFINE_USER_POINT_TRAITS(fusion::vector2<double, double>,double,2);
BOOST_DEFINE_USER_POINT_TRAITS(fusion::vector2<int, int>,int,2);

Next you need to specialize the access traits for the coordinate framework 
which the point implicitly describes.

//! Access traits for points in Cartesian coordinates
template <typename Coordinate>
struct cartesian_access_traits<fusion::vector2<Coordinate, Coordinate> >
{
   typedef Point point_type;
   typedef typename point_traits<point_type>::coordinate_type 
coordinate_type;
   typedef typename point_traits<point_type>::dimension_type dimension_type;

   template <unsigned int D>
   static inline coordinate_type get(const point_type& p)
   {
       BOOST_MPL_ASSERT_MSG
       (
          ( dimension_traits< D >::value >= 0 &&
            dimension_traits< D >::value < dimension_type::value )
          , POINT_GET_CALLED_WITH_INDEX_OUT_OF_BOUNDS
          , ( dimension_traits< D > )
       );
       return fusion::at_c< D >::( p );
   }

   //! factory traits.
   static inline point_type construct(const coordinate_type& x,
                                                       const 
coordinate_type& y)
   { return point_type( x, y ); }

};

Then the type may be passed directly into algorithms:

//! Function to find the cross product between two vectors formed by the 
specified three points which share an endpoint at A.
template <typename Point>
point_traits<Point>::coordinate_type cross_product( const Point& A,
                                                                             
       const Point& B,
                                                                             
       const Point& C )
{
   typedef cartesian_access_traits<Point > access_traits;
   boost::function_requires<Point2DConcept<Point> >();
   boost::function_requires
   <
      CartesianCoordinateAccessorConcept<access_traits>
   >();

   typename point_traits< Point >::coordinate_type cross =
   (
      ( access_traits::get<0>( B ) - access_traits::get<0>( A ) )  *
      ( access_traits::get<1>( C ) - access_traits::get<1>( A ) )  -
      ( access_traits::get<0>( C ) - access_traits::get<0>( A ) )  *
      ( access_traits::get<1>( B ) - access_traits::get<1>( A ) )
   );
   return cross;
}

//! Orientation test to check if point C is left, collinear, or right of the 
line formed by A-B.
enum orientation_type
{
        oriented_right     = -1,
        oriented_collinear = 0,
        oriented_left      = 1
};
template <typename PrecisionPolicy, typename Point>
orientation_type get_orientation( const Point& A,
                                                    const Point& B,
                                                    const Point& C,
                                                    const PrecisionPolicy& 
compare )
{
        typedef typename point_traits<Point>::coordinate_type 
coordinate_type;
        coordinate_type cross = cross_product( A, B, C );
        coordinate_type zero = 0;
        if( compare.less_than( cross, zero ) )
        {
                return oriented_right;
        }
        else if( compare.greater_than( cross, zero ) )
        {
                return oriented_left;
        }
        else
        {
                return oriented_collinear;
        }
}

So you could do the following:

typedef fusion::vector2<double,double> point;
point p1( 0., 0. );
point p2( 1., 0. );
point p3( 2., 0. );
orientation_type oType = get_orientation( p1, p2, p3,
      fraction_tolerance_comparison_policy<double>(1e-10) );
assert( oType == oriented_collinear );

Or with ints:

typedef fusion::vector2<int,int> point;
point p1( 0, 0 );
point p2( 1, 1 );
point p3( 2, 2 );
orientation_type oType = get_orientation( p1, p2, p3,
       direct_comparison_policy() );
assert( oType == oriented_collinear );

When I was talking about access like tuple, I meant using compile time 
indexing via the access traits. One of the main goals of my library is to 
facilitate use with legacy geometry code (proprietary in my case) simply by 
specializing access traits. Using access traits results in little constraint 
on the interface of the underlying point type. This frees one to develop 
algorithms for a public library while maintaining the ability to use said 
algorithms with the legacy code (without worrying about translating into a 
boost::point type.) The same idea is used segments and polylines/polygons.

> First, some algorithms require runtime indexing. kd-trees are an 
> interesting
> example. When constructing a kd-tree over a data set, we need some 
> heuristic
> to decide which dimension to split on at each tree node. In 2 or 3
> dimensions, usually you simply cycle through the dimensions in order, for
> which in principle only compile-time indexing is necessary. But this
> heuristic is not so appropriate in high dimensions; then it is usual to
> split each tree node on the dimension with the greatest variance in the
> data, which of course is unknown at compile time.
>
> Second, there are efficiency and code size concerns. Consider the 
> Euclidean
> distance function. For 2D/3D points, it is nice to use compile-time 
> indexing
> and metaprogramming to guarantee tight, loop-less code. For high 
> dimensional
> points, however, beyond some dimension
> (BOOST_GEOMETRY_MAX_UNROLL_DIMENSION?) unrolling the loop is no longer
> desirable, and we would rather just do the runtime iteration. If the 
> points
> are RuntimeIndexable, we can choose an appropriate distance algorithm at
> compile time depending on the points' dimensionality.
>

I presume the reason one would not want loop-less code at high dimensions is 
due to bloat? I do not understand why runtime looping would be faster in 
this case. It would seem that even in the largest cases the bloat would 
still be small compared to memory sizes today. Runtime indexing can 
certainly be accommodated using a tag dispatch on the access traits. I think 
more research needs to be done to quantify the differences between the two 
to see if it is worth the maintenance overhead.

Cheers,

Brandon

> I do think you are on the right track, and I look forward to seeing 
> preview
> 3.
>
> Cheers,
> Patrick
> _______________________________________________
> Unsubscribe & other changes: 
> http://lists.boost.org/mailman/listinfo.cgi/boost
>