boost::geometry でもハマった

同じ図形を geometry でも処理してみたら、見事に落ちた。 ロバストへの道は険しいのであった(´・ω・`)

// Boost.Geometry (aka GGL, Generic Geometry Library)

// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
// Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
// Copyright (c) 2009-2012 Mateusz Loskot, London, UK.

// 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 <algorithm> // for reverse, unique
#include <iostream>
#include <string>

#include <boost/geometry/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/geometries/adapted/c_array.hpp>
#include <boost/geometry/multi/geometries/multi_polygon.hpp>

BOOST_GEOMETRY_REGISTER_C_ARRAY_CS(cs::cartesian)


int main(void)
{
    using namespace boost::geometry;

    typedef model::d2::point_xy<double> point_2d;
    typedef model::polygon<point_2d> polygon_2d;
    typedef model::box<point_2d> box_2d;

    // Define a polygon and fill the outer ring.
    // In most cases you will read it from a file or database
    polygon_2d poly;
    {
#ifndef SEGV
      const double coor[][2] = {
            {2.0, 1.3}, {2.4, 1.7}, {2.8, 1.8}, {3.4, 1.2}, {3.7, 1.6},
            {3.4, 2.0}, {4.1, 3.0}, {5.3, 2.6}, {5.4, 1.2}, {4.9, 0.8}, {2.9, 0.7},
            {2.0, 1.3} // closing point is opening point
            };
#else
        const double coor[][2] = {
          { -162.277344, -684.941406 },
          { -154.300781, -684.707031 },
          { -147.291016, -768.253906 },
          { -148.601562, -752.630859 },
          { -156.675781, -753.292969 },
          { -155.425781, -768.546875 }
            };
#endif
      assign_points(poly, coor);
    }

    // Polygons should be closed, and directed clockwise. If you're not sure if that is the case,
    // call the correct algorithm
    correct(poly);

  
    // Polygons can be streamed as text
    // (or more precisely: as DSV (delimiter separated values))
    std::cout << dsv(poly) << std::endl;
  
    typedef std::vector<polygon_2d> polygon_list;
    polygon_list v;

    try {
      intersection(poly, poly, v);
    } catch(const std::exception& e) {
      std::cerr << e.what() << std::endl;
    }

    std::cout << "Clipped output polygons" << std::endl;
    for (polygon_list::const_iterator it = v.begin(); it != v.end(); ++it)
    {
        std::cout << dsv(*it) << std::endl;
    }
  
    return 0;
}

今回は vc 2008 だったのだ… /D SEGV を付けないケースでは、落ちる事はない。

cl /EHsc /GR /MD test.cpp /D SEGV

2014/03/07 追記: geos で selfIntersection にかけてみたが、やはりダメだった。

TopologyException: side location conflict at -752.63088031450525 -148.60182192355231

追記：例外処理しないといけないですね。キャッチされない例外で落ちるのは、あたりまえですね。修正しました。

Boost.Geometry Overlay invalid input exception

になりました。 追記：boost::geometry has_self_intersections.hpp でスローされる例外のようで、自己交差の無い Simple な Polygon が入力として必要なようである。用語が難しいのだが、一般に、図形の Simplify というと、複雑な点数を簡略化する（点を間引く）操作を言うので、自己交差のないポリゴンにばらす関数があれば良いのだろうか。