CGALでランダム点をドロネー三角形分割してParaViewで観てみる

期待計算量 $O(n\log n)$ のドロネー三角形分割を実装しようと思ったけど難しそうなので，とりあえずCGALという情報幾何ライブラリに任せちゃうことにした．

CGALのインストール

https://www.cgal.org/download/linux.html 簡単．

ランダム点の生成

このページを参考にした．生成器を作ってそこから取る感じだった．

#include <CGAL/Simple_cartesian.h>
#include <CGAL/point_generators_2.h>
#include <vector>

using namespace CGAL;
using K = Simple_cartesian<double>;
using Point = K::Point_2;
using Creator = Creator_uniform_2<double,Point>;

int main() {
    // Create point set. Prepare a vector for 1000 points.
    std::vector<Point> points;
    points.reserve(1000);
    // Create 1000 points from [-1,1)x[-1,1)
    Random_points_in_square_2<Point, Creator> g(1.0);
    cpp11::copy_n(g, 1000, std::back_inserter(points));
    return 0;
}

ドロネー三角形分割

大変苦戦したが，なんとか書けた．ただ分割するだけではなく頂点にインデックスの情報を持たせたかったのですこし大変．

#include <CGAL/Simple_cartesian.h>
#include <CGAL/point_generators_2.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <vector>

using K         = CGAL::Simple_cartesian<double>;
using Vb        = CGAL::Triangulation_vertex_base_with_info_2<int, K>;
using Tds       = CGAL::Triangulation_data_structure_2<Vb>;
using Delaunay  = CGAL::Delaunay_triangulation_2<K, Tds>;
using Point     = Delaunay::Point;
using Creator   = CGAL::Creator_uniform_2<double,Point>;


int main() {
    // Create point set. Prepare a vector for 1000 points.
    std::vector<Point> points_;
    std::vector<std::pair<Point, int>> points;
    points.reserve(1000);
    // Create 1000 points from [-1,1)x[-1,1)
    CGAL::Random_points_in_square_2<Point, Creator> g(1.0);
    CGAL::cpp11::copy_n(g, 1000, std::back_inserter(points_));
    for(int i=0;i<1000;i++){
        points.push_back({points_[i], i});
    }

    // Delaunay triangulation
    Delaunay t;
    t.insert(points.begin(), points.end());

    // output
    int n_vertices = t.number_of_vertices();
    int n_faces    = t.number_of_faces();

    std::cout << "# vtk DataFile Version 3.0" << std::endl;
    std::cout << "title" << std::endl;
    std::cout << "ASCII" << std::endl;
    std::cout << "DATASET UNSTRUCTURED_GRID" << std::endl;

    std::cout << "POINTS " << n_vertices << " float" << std::endl;
    for(auto point : points_){
        std::cout << point << " 0" << std::endl;
    }

    std::cout << "CELLS " << n_vertices + n_faces << " " << n_vertices*2 + n_faces*4 << std::endl;
    for(int i=0;i<n_vertices;i++){
        std::cout << "1 " << i << std::endl;
    }
    for(auto it = t.finite_faces_begin();it != t.finite_faces_end();it++){
        std::cout << "3 " << it->vertex(0)->info() << " "
                          << it->vertex(1)->info() << " "
                          << it->vertex(2)->info() << std::endl;
    }

    std::cout << "CELL_TYPES " << n_vertices + n_faces << std::endl;
    for(int i=0;i<n_vertices;i++){
        std::cout << "1" << std::endl;
    }
    for(int i=0;i<n_faces;i++){
        std::cout << "5" << std::endl;
    }

    return 0;
}

ParaViewで見てみる

最近ParaViewを使ってるのでvtkフォーマットで結果を出力してみた．そのためにインデックスもたせたのもある． f:id:monman53:20181116051824p:plain できていると信じたい...

終わりに

疲れましたが，次からは使えそうです．

monman53のぶろぐ

いろいろ載せるよ