CGAL Hello World してみる - monman53のぶろぐ

Computational Geometry Algorithms Library (CGAL) 便利そうなので入門してみる．

これをやっておかないとCGALの哲学がわからなそうな気がしたのでちゃんとメモ．このページを超ざっくり説明できたらと思う．

そのまえにCGALをインストールしよう．公式のダウンロードページの指示に従えば問題なくインストールできた．

ではHelloWorldに戻ろう．

1. 3つの点と1つの線分

points_and_segment.cpp

#include <iostream>
#include <CGAL/Simple_cartesian.h>
typedef CGAL::Simple_cartesian<double> Kernel;
typedef Kernel::Point_2 Point_2;
typedef Kernel::Segment_2 Segment_2;
int main()
{
  Point_2 p(1,1), q(10,10);
  std::cout << "p = " << p << std::endl;
  std::cout << "q = " << q.x() << " " << q.y() << std::endl;
  std::cout << "sqdist(p,q) = " 
            << CGAL::squared_distance(p,q) << std::endl;
  
  Segment_2 s(p,q);
  Point_2 m(5, 9);
  
  std::cout << "m = " << m << std::endl;
  std::cout << "sqdist(Segment_2(p,q), m) = "
            << CGAL::squared_distance(s,m) << std::endl;
  std::cout << "p, q, and m ";
  switch (CGAL::orientation(p,q,m)){
  case CGAL::COLLINEAR: 
    std::cout << "are collinear\n";
    break;
  case CGAL::LEFT_TURN:
    std::cout << "make a left turn\n";
    break;
  case CGAL::RIGHT_TURN: 
    std::cout << "make a right turn\n";
    break;
  }
  std::cout << " midpoint(p,q) = " << CGAL::midpoint(p,q) << std::endl;
  return 0;
}

とりあえずコンパイルしてみる．

$ g++ points_and_segment.cpp -lCGAL
$ ./a.out
p = 1 1
q = 10 10
sqdist(p,q) = 162
m = 5 9
sqdist(Segment_2(p,q), m) = 8
p, q, and m make a left turn
 midpoint(p,q) = 5.5 5.5

このプログラムでは点の定義からその出力，2点の距離(sqrtは時間がかかるためCGALではよく二乗距離を使うようだ)，線分と点の距離，三点の位置関係，中点を求める様子を確認できる．

surprising.cpp

#include <iostream>
#include <CGAL/Simple_cartesian.h>
typedef CGAL::Simple_cartesian<double> Kernel;
typedef Kernel::Point_2 Point_2;
int main()
{
  {
    Point_2 p(0, 0.3), q(1, 0.6), r(2, 0.9);
    std::cout << (CGAL::collinear(p,q,r) ? "collinear\n" : "not collinear\n");   
  }
  {
    Point_2 p(0, 1.0/3.0), q(1, 2.0/3.0), r(2, 1);
    std::cout << (CGAL::collinear(p,q,r) ? "collinear\n" : "not collinear\n");   
  }  
  {
    Point_2 p(0,0), q(1, 1), r(2, 2);
    std::cout << (CGAL::collinear(p,q,r) ? "collinear\n" : "not collinear\n");   
  }  
  return 0;
}

not collinear
not collinear
collinear

collinearは3点が同一直線状に乗っているか判定する関数である．いずれの結果も"collinear"になりそうだが，前2つは丸め誤差によって"not collinear"と判定されている．

exact.cpp

#include <iostream>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <sstream>
typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
typedef Kernel::Point_2 Point_2;
int main()
{
  Point_2 p(0, 0.3), q, r(2, 0.9);
  {
    q  = Point_2(1, 0.6);
    std::cout << (CGAL::collinear(p,q,r) ? "collinear\n" : "not collinear\n");
  }
  
  {
    std::istringstream input("0 0.3   1 0.6   2 0.9");
    input >> p >> q >> r;
    std::cout << (CGAL::collinear(p,q,r) ? "collinear\n" : "not collinear\n");
  }
  
  {
    q = CGAL::midpoint(p,r);
    std::cout << (CGAL::collinear(p,q,r) ? "collinear\n" : "not collinear\n");   
  }
 
  return 0;
}

not collinear
collinear
collinear

今回の"点"は先の"点"と定義される"kernel"が異なるため別物である．これを用いれば正しい判定が行われるが，そのためには数値の扱いに注意しなければならない．もちろん計算速度などの兼ね合いでどのkernelを使うべきか選ぶべきである．たくさんあるCGALのパッケージは，それぞれサポートするkernelが明記されているからその都度確認されたい．

2. 凸包と点の列

点集合の凸法を例に，関数へのデータの入出力のやりかたを見ていく．

array_convex_hull_2.cpp

#include <iostream>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/convex_hull_2.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_2 Point_2;
int main()
{
  Point_2 points[5] = { Point_2(0,0), Point_2(10,0), Point_2(10,10), Point_2(6,5), Point_2(4,1) };
  Point_2 result[5];
  Point_2 *ptr = CGAL::convex_hull_2( points, points+5, result );
  std::cout <<  ptr - result << " points on the convex hull:" << std::endl;
  for(int i = 0; i < ptr - result; i++){
    std::cout << result[i] << std::endl;
  }
  return 0;
}

これをコンパイルするにはGMPへのリンクが必要なようだ．

$ g++11  main.cpp -lCGAL -lgmp
$ ./a.out
3 points on the convex hull:
0 0
10 0
10 10

配列を用いる場合は，convex_hull_2関数に入力のポインタと出力のポインタを渡す．そして最後の凸法の点の後のポインタが返ってくる．

vector_convex_hull_2.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/convex_hull_2.h>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_2 Point_2;
typedef std::vector<Point_2> Points;
int main()
{
  Points points, result;
  points.push_back(Point_2(0,0));
  points.push_back(Point_2(10,0));
  points.push_back(Point_2(10,10));
  points.push_back(Point_2(6,5));
  points.push_back(Point_2(4,1));
  CGAL::convex_hull_2( points.begin(), points.end(), std::back_inserter(result) );
  std::cout << result.size() << " points on the convex hull" << std::endl;
  return 0;
}

3 points on the convex hull

STLのvectorを使った場合はこのようにやる．back_inserterとか使えるようになりたいぞ!

どっちを使っていくべきかはその人の趣味に委ねられているようだ．後者はメモリの確保が不要になっているね．

3. Kernels と Traits について

ここらへんから知りたい内容が始まった．

先の凸法で用いたconvex_hull_2のドキュメントをながめていると，別のバージョンのconvex_hull_2があることに気づく．

template<class InputIterator , class OutputIterator , class Traits >
OutputIterator 
convex_hull_2(InputIterator first,
              InputIterator beyond,
              OutputIterator result,
              const Traits & ch_traits)

Traitsという引数が1つ増えている．まだあまり理解していないが，内部で実行されるアルゴリズムに使う基本的な処理を指定できるっぽい．ここらへん使いこなせると相当便利そう．

There are two obvious questions: What can be used as argument for this template parameter? And why do we have template parameters at all? To answer the first question, any model of the CGAL concept Kernel provides what is required by the concept ConvexHullTraits_2. As for the second question, think about an application where we want to compute the convex hull of 3D points projected into the yz plane. Using the class Projection_traits_yz_3 this is a small modification of the previous example.

4. Concepts と Models について

なんかわからない... Conceptというのは要件のことで，Modelはその要件を満たすもの?

書かれている例を見るとなんとなく分かった気になる．

5 Further Reading

"The C++ Standard Library, A Tutorial and Reference" by Nicolai M. Josuttis from Addison-Wesley や "Generic Programming and the STL" by Matthew H. Austern が勧められていた．チュートリアルも続くようだが，現時点でこのHello Worldしかない．

おわりに

最後の方は超適当になってしまったが，最後の方ほど大事なので後で加筆したい．CGALってどのくらい使われているのだろう...