2014年3月13日木曜日

boost::filesystem v3 で windows 環境にて日本語ファイルを扱う備忘録

Default encoding under windows を参考にやってみました。 コンパイラが VC8 だったので、最初、手を抜いてソースコードを UTF-8 でファイルに保存して実行してみたら
#include <boost/locale.hpp>
#include <boost/filesystem/path.hpp>
#include <boost/filesystem/fstream.hpp>
#include <boost/format.hpp>

int main()
{
  // Create and install global locale
  std::locale::global(boost::locale::generator().generate(""));
  //std::locale::global(boost::locale::generator().generate("ja_JP.UTF-8"));
  // Make boost.filesystem use it
  boost::filesystem::path::imbue(std::locale());
  // Now Works perfectly fine with UTF-8!
  boost::filesystem::ofstream hello("こんにちわ.txt"); 
  return 0;
}
出力されたファイル名が
こんにちめEtxt
なんじゃこれ??? あれこれ、試しているうちに、どうやら UTF-8 で保存したのがいけないみたいだという事で
#include <boost/locale.hpp>
#include <boost/filesystem/path.hpp>
#include <boost/filesystem/fstream.hpp>
#include <boost/format.hpp>

int main()
{
  // Create and install global locale
  std::locale::global(boost::locale::generator().generate(""));
  //std::locale::global(boost::locale::generator().generate("ja_JP.UTF-8"));
  // Make boost.filesystem use it
  boost::filesystem::path::imbue(std::locale());
  // Now Works perfectly fine with UTF-8!
  //boost::filesystem::ofstream hello("こんにちわ.txt"); 
                                   //  こ          ん          に          ち          わ
  boost::filesystem::ofstream hello("\xE3\x81\x93\xE3\x82\x93\xE3\x81\xAB\xE3\x81\xA1\xE3\x82\x8F.txt"); 
  return 0;
}
と正攻法で書きなおしたら、ちゃんと動作しました。 また、くだらない事に時間を費やしてしまった(´・ω・`)

2014年3月12日水曜日

boost::polygon self intersection problem 顛末

この前のポストの問題ですが、面白い事がわかったので、顛末を書きます。 OpenGLES では、描画が float 精度までで、データは double 精度で保持しておりました。 これを boost::polygon にて、trapezoid にかけて台形化処理するのですが、その時に double から float に static_cast で精度を落とします。そうすると、自己交差の無い simple な図形が精度変換の過程で自己交差を持つ図形に変身するケースがあると判明しました。  じゃあ事前に float 精度に落として 自己交差のある状態の図形は float 精度で自己交差を解消した simple な状態にしてから double に戻せば問題無いだろうと思って、やってみましたが、何故かうまく行かない。詳しく突っ込んで調べていないですが、もしかしたら CPU に依存する話なのかもしれません。  余計な事をすると遅くなるので、なんとかしたかったですが、アプリが落ちてしまっては、元も子もありません。もう少し無駄な計算は省こうという事で、こんな感じで落ち着きました。
    void to_simple( Polygon& polygon ) {
      size_t len = polygon.self_.coords_.size();
      Point ss = polygon.self_.coords_[len-1];
      for( size_t i = 0; i < len-2; ++i ) {
        Point se = polygon.self_.coords_[i];
        Point es = polygon.self_.coords_[i+1];
        for( size_t j = i + 2; j < len; ++j ) {
          Point ee = polygon.self_.coords_[j];
          if( i != j ) {
            if( 
                 (ss.x() == es.x() && ss.y() == es.y() )
              || (ss.x() == ee.x() && ss.y() == ee.y() )
              || (se.x() == es.x() && se.y() == es.y() )
              || (se.x() == ee.x() && se.y() == ee.y() )
            ) {
            } else {
              double xx, yy;
              if( k_4t_close( xx, yy, ss.x(), ss.y(), se.x(), se.y(), es.x(), es.y(), ee.x(), ee.y() ) ) {
                se = Point( xx, yy );
                ee = Point( xx, yy );
                //if( i < j ) { // i < j は自明でしたねorz
                  polygon.self_.coords_.insert( polygon.self_.coords_.begin() + j, ee );
                  polygon.self_.coords_.insert( polygon.self_.coords_.begin() + i, se );
                //} else {
                //  polygon.self_.coords_.insert( polygon.self_.coords_.begin() + i, se );
                //  polygon.self_.coords_.insert( polygon.self_.coords_.begin() + j, ee );
                //}
                len += 2;
              }
            }
          }
          es = ee;
        }
        ss = se;
      }
    }

2014年3月7日金曜日

boost::polygon self intersection problem

ここ2つ前からの投稿の問題は解決しました。この前に geos というライブラリでもテストコードを実行させてみたところ
TopologyException: side location conflict at -752.63088031450525 -148.60182192355231
という例外が発生しました。boost::geometry はロバストに出来ていて、自分が例外キャッチするのを忘れていただけで、ちゃんと
Boost.Geometry Overlay invalid input exception
という例外が送出されていました。  で、どういう問題が発生してたか boost::geometryのコードを読んでみたところ、図形が Simple でない=自己交差を起こしている場合に例外が送出されていました。 なので、自己交差を解消してやれば、問題をクリアできそうな感触を得ました。 残念ながら、自己交差を解消する関数は、見当たらなかったので自作。スタイルが古いですが、古いコンパイラを使っているので auto とか使えなかったんです。 list::insert してるあたりがダサいですが、古い環境なんで、このままです(>_<) 結局、自己交差を解消してやれば、boost::polygon のコードも通りました。
#include <boost/cstdint.hpp>
#include <boost/polygon/polygon.hpp>
#include <vector>
#include <deque>
#include <cmath>
#include <iomanip>

namespace bpl = boost::polygon;
typedef double point_value_type;
typedef bpl::segment_data<point_value_type> Segment;
typedef bpl::polygon_with_holes_data<point_value_type> Polygon;
typedef bpl::polygon_data<point_value_type> Ring;
typedef bpl::polygon_traits<Polygon>::point_type Point;
typedef std::vector<Point> PointSet;
typedef std::vector<Polygon> PolygonSet;
typedef std::list<Segment> SegmentSet;

#define ACCY 1e-6

Polygon make_polygon( const point_value_type* pts, size_t count ) {
  Polygon polygon;
  polygon.self_.coords_.resize( count );
  for( size_t i = 0; i < count; ++i ) {
    polygon.self_.coords_[i].x( *pts++ );
    polygon.self_.coords_[i].y( *pts++ );
  }
  return polygon;
}

void construct( PolygonSet& polygons, Polygon polygon ) {
  using namespace bpl::operators;
  polygons |= polygon; 
}

void construct( SegmentSet& segments, const point_value_type* pts, size_t count ) {
  point_value_type ax = pts[count*2-2];
  point_value_type ay = pts[count*2-1];
  for( size_t i = 0; i < count; ++i ) {
    point_value_type bx = *pts++;
    point_value_type by = *pts++;
    segments.push_back( Segment( Point(ax, ay), Point(bx, by) ) );
    ax = bx;  ay = by;
  }
}

Polygon make_polygon(const SegmentSet& segments ) {
  Polygon polygon;
  size_t count = segments.size();
  polygon.self_.coords_.resize( count );
  size_t i = 0;
  SegmentSet::const_iterator s = segments.begin();
  SegmentSet::const_iterator e = segments.end();
  while( s != e ) {
    polygon.self_.coords_[i++] = bpl::low(*s++);
  }
  return polygon;
}


std::ostream& operator << (std::ostream& os, const Point& p) {
  os << std::setprecision(10);
  os << "{ " << bpl::x(p) << ", " << bpl::y(p) << " },\n";
  return os;
}

std::ostream& operator << (std::ostream& os, const SegmentSet& segments) {
  SegmentSet::const_iterator s = segments.begin();
  SegmentSet::const_iterator e = segments.end();
  while( s != e ) {
    os << bpl::low( *s++ );
  }
  return os;
}

//! 4点交点(線分による完全交差版)を取得する
/*!
 @attention 交点が存在しない場合、x,y は未定義である
 @retval true 交点が存在
 @retval false 交点は存在しない
*/
bool k_4t_close( 
 double& x, //!< [out] 交点X座標値 
 double& y, //!< [out] 交点Y座標値
 double xk, //!< [in] 線分KLの始点X座標値
 double yk, //!< [in] 線分KLの始点Y座標値
 double xl, //!< [in] 線分KLの終点X座標値
 double yl, //!< [in] 線分KLの終点Y座標値
 double xm, //!< [in] 線分MNの始点X座標値
 double ym, //!< [in] 線分MNの始点Y座標値
 double xn, //!< [in] 線分MNの終点X座標値
 double yn   //!< [in] 線分MNの終点Y座標値
) {
 long double xlk = xl - xk;
 long double ylk = yl - yk;
 long double xnm = xn - xm;
 long double ynm = yn - ym;
 long double xmk = xm - xk;
 long double ymk = ym - yk;
 long double det = xnm * ylk - ynm * xlk;
 if( std::fabs( det ) < ACCY ) return false;
 long double s = (xnm * ymk - ynm * xmk ) / det;
 long double t = (xlk * ymk - ylk * xmk ) / det;
 if( s < 0.0 || 1.0 < s || t < 0.0 || 1.0 < t ) return false;
 x = xk + xlk * s;
 y = yk + ylk * s;
 return true;
}

void resolve_segment( SegmentSet& segments ) {
  SegmentSet::iterator is = segments.begin();
  SegmentSet::iterator ie = segments.end();
  while( is != ie ) {
    SegmentSet::iterator cs = segments.begin();
    SegmentSet::iterator ce = segments.end();
    while( cs != ce ) {
      if( cs != is ) {
        if( bpl::intersects( *cs, *is, false ) ) {
          //std::cout << "X" << std::endl;
          double xx, yy;
          Point a = bpl::low( *cs );
          Point b = bpl::high( *cs );
          Point c = bpl::low( *is );
          Point d = bpl::high( *is );
          k_4t_close( xx, yy,
            bpl::x( a ), bpl::y( a ), bpl::x( b ), bpl::y( b ),
            bpl::x( c ), bpl::y( c ), bpl::x( d ), bpl::y( d )
          );
          Point n( xx, yy );
          //std::cout << a << n << b << std::endl;
          *cs = Segment( n, b );
          cs = segments.insert( cs, Segment( a, n ) );
          *is = Segment( n, d );
          is = segments.insert( is, Segment( c, n ) );
        } else {
          //std::cout << "o" << std::endl;
        }
      }
      ++cs;
    }
    ++is;
  }
}


int main() {
  point_value_type pts[] = {
   -162.277344, -684.941406,
   -154.300781, -684.707031,
   -147.291016, -768.253906,
   -148.601562, -752.630859,
   -156.675781, -753.292969,
   -155.425781, -768.546875,
  };
  /*
  point_value_type pts2[] = {
   -155.425781, -768.546875 ,
   -162.277344, -684.941406 ,
   -154.300781, -684.707031 , 
   -148.6018219, -752.6308803 ,
   -147.291016, -768.253906 ,
   -148.601562, -752.630859 ,
   -148.6018219, -752.6308803 ,
   -156.675781, -753.292969 ,
  };
  */
  
  SegmentSet segments;
  construct( segments, pts, 6 );

  std::cout << segments;
  
  resolve_segment( segments );
  
  std::cout << "===== RESOLVE ====\n";
  std::cout << segments;
  
  PolygonSet polygons;
  //construct( polygons, make_polygon( pts, 6 ) );
  construct( polygons, make_polygon( segments ) );
  
  return 0;
}


追記:あれ? insert のあたり、よく考えたら、まずいですかね・・・。

2014年3月6日木曜日

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 というと、複雑な点数を簡略化する(点を間引く)操作を言うので、自己交差のないポリゴンにばらす関数があれば良いのだろうか。

2014年3月5日水曜日

boost::polygon はまった

三角形分割するのに、boost::polygon を使ってたんですが、残念ながら実数版は実装されているもののサポート対象外みたいで、ロバスト性が保証されていないようです。 Android をはじめとする GLES では、Tesselation(三角形分割)がサポートされていないので、Trapezoid (台形)化してポリゴンを塗り潰してたんですけども、残念です。 以下、削ったコード。
#include <boost/cstdint.hpp>
#include <boost/polygon/polygon.hpp>
#include <vector>

using boost::int32_t;
using boost::uint8_t;

namespace bpl = boost::polygon;
typedef float point_value_type;
typedef bpl::polygon_with_holes_data<point_value_type> Polygon;
typedef bpl::polygon_data<point_value_type> Ring;
typedef bpl::polygon_traits<Polygon>::point_type Point;
typedef std::vector<Point> PointSet;
typedef std::vector<Polygon> PolygonSet;

void set_points( Polygon& polygon, const point_value_type* pts, int32_t count ) {
  polygon.self_.coords_.resize( count );
  for( int32_t i = 0; i < count; ++i ) {
    polygon.self_.coords_[i].x( *pts++ );
    polygon.self_.coords_[i].y( *pts++ );
  }
}

void construct( PolygonSet& polygons, const point_value_type* pts, int32_t count ) {
  using namespace bpl::operators;
  Polygon polygon;
  set_points( polygon, pts, count );
  polygons |= polygon;  //// SIGSEGV
}


int main() {
  float pts[] = {
   -30162.277344, -21684.941406,
   -30154.300781, -21684.707031,
   -30147.291016, -21768.253906,
   -30148.601562, -21752.630859,
   -30156.675781, -21753.292969,
   -30155.425781, -21768.546875,
  };
  
  PolygonSet polygons;
  construct( polygons, pts, 6 );

  
  return 0;
}


追記:どんな図形か調べてみた。ま、往々にしてありがちな図形やね。平行な直線が折り返しているとダメなのであろうか?普通、平行な直線は交点の計算ができないので、交点を求めたりはしない。これをやろうとすると泥沼に陥る事が多い。
2014/03/06 追記: 整数の同じ構成にしてテストを行うと、ちゃんとパスする。
#ifndef SEGF
  typedef int32_t point_value_type;
#else
  typedef float point_value_type;
#endif
にして
  point_value_type pts[] = {
#ifndef SEGF
   -22277344,  -4941406,
   -14300781,  -4707031,
    -7291016, -88253906,
    -8601562, -72630859,
   -16675781, -73292969,
   -15425781, -88546875,
#else
   -22.277344,  -4.941406,
   -14.300781,  -4.707031,
    -7.291016, -88.253906,
    -8.601562, -72.630859,
   -16.675781, -73.292969,
   -15.425781, -88.546875,
#endif
  };
ここは、難しいようで、 polygon_arbitrary_format.hpp にもデバッグ用のコードがコメントアウトされる形で埋め込まれていた。 両者で実行して、違いを見比べてみたところ セーフバージョン
processEvent_
loop
current Y -72630881
scanline size 2
( -16675781, -73292969), ( -8601822, -72630881); ( -14300781, -4707031), ( -8601822, -72630881); 
found element in scanline 1
finding elements in tree
first iter y is -7.26309e+007
loop2
loop2
counts_from_scanline size 2
aggregating
loop3
loop3
落ちるバージョン
processEvent_
loop
current Y -72.6309
scanline size 2
( -16.6758, -73.293), ( -8.60156, -72.6309); ( -14.3008, -4.70703), ( -7.29102, -88.2539); 
found element in scanline 1
finding elements in tree
first iter y is -72.6309
loop2
counts_from_scanline size 1
aggregating
loop3
loop3
更に、loop2 の条件判定文のコメントを polygon_arbitrary_format.hpp に追加して追試してみたところ、
        if(iter != scanData_.end())
          std::cout << "first iter y is " << iter->first.evalAtX(x_) << "\n";
        while(iter != scanData_.end() &&
              ((iter->first.pt.x() == x_ && iter->first.pt.y() == currentY) ||
               (iter->first.other_pt.x() == x_ && iter->first.other_pt.y() == currentY))) {
                //iter->first.evalAtX(x_) == (high_precision)currentY) {
          std::cout << "loop2\n";
          elementIters.push_back(iter);
          counts_from_scanline.push_back(std::pair<std::pair<std::pair<Point, Point>, int>, active_tail_arbitrary*>
                                         (std::pair<std::pair<Point, Point>, int>(std::pair<Point, Point>(iter->first.pt,
                                                                                                          iter->first.other_pt),
                                                                                  iter->first.count),
                                          iter->second));
          ++iter;
          if( iter != scanData_.end() ) {
            std::cout << "CHK loop2 " << iter->first.pt.x() << " == " << x_ << " && " << iter->first.pt.y() << " == " << currentY << "\n";
            std::cout << "CHK loop2 " << iter->first.other_pt.x() << " == " << x_ << " && " << iter->first.other_pt.y() << " == " << currentY << "\n";
          }
        }
セーフバージョンと落ちるバージョンでは、そもそもscanline の構成が異なっている事がわかった
processEvent_
loop
current Y -72630881
scanline size 2
( -16675781, -73292969), ( -8601822, -72630881); ( -14300781, -4707031), ( -8601822, -72630881); 
found element in scanline 1
finding elements in tree
first iter y is -7.26309e+007
loop2
CHK loop2 -14300781 == -8601822 && -4707031 == -72630881
CHK loop2 -8601822 == -8601822 && -72630881 == -72630881
loop2
counts_from_scanline size 2
aggregating
processEvent_
loop
current Y -72.6309
scanline size 2
( -16.6758, -73.293), ( -8.60156, -72.6309); ( -14.3008, -4.70703), ( -7.29102, -88.2539); 
found element in scanline 1
finding elements in tree
first iter y is -72.6309
loop2
CHK loop2 -14.3008 == -8.60156 && -4.70703 == -72.6309
CHK loop2 -7.29102 == -8.60156 && -88.2539 == -72.6309
counts_from_scanline size 1
aggregating
loop3
loop3
浮動小数点型の演算は、難しいものだ。やはり、ロバスト性を確保するのは並大抵のことではない感じだ。