tak0kadaの何でもノート

発声練習、生存確認用。

医学関連は 医学ノート

CGALのサンプルを実行する(その3)

http://doc.cgal.org/latest/Point_set_processing_3/index.html#Point_set_processing_3NormalOrientation を少し変更して実行した。

法線推定

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/pca_estimate_normals.h>
#include <CGAL/mst_orient_normals.h>
#include <CGAL/property_map.h>
#include <CGAL/IO/read_off_points.h>
#include <CGAL/IO/write_xyz_points.h>

#include <cstring> // to use strcmp
#include <string>
#include <utility>
#include <list>
#include <fstream>
#include <iostream>
#include <cstdlib> // to use EXIT_XXX macros

typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;
// Point with normal vector stored in a std::pair.
typedef std::pair<Point, Vector> PointVectorPair;
// Concurrency
#ifdef CGAL_LINKED_WITH_TBB
    typedef CGAL::Parallel_tag Concurrency_tag;
#else
    typedef CGAL::Sequential_tag Concurrency_tag;
#endif

int main(int argc, char* argv[])
{
    std::list<PointVectorPair> points;

    if ( argc != 3 || !std::strcmp(argv[1], argv[2]) )
    {
        std::cerr << "usage: ./a.out input.off output.off" << std::endl
                  << "       input and output must have different names" << std::endl;
        return EXIT_FAILURE;
    }
    const std::string input_filename = argv[1];
    const std::string output_filename = argv[2];

    std::ifstream input(input_filename.c_str());
    if (!input ||
        !CGAL::read_off_points(input
                             , std::back_inserter(points)
                             , CGAL::First_of_pair_property_map<PointVectorPair>()))
    {
        std::cerr << "Error: cannot read " << input_filename << std::endl;
        return EXIT_FAILURE;
    }

    // estimate normals
    const int nb_neighbors = 18;
    CGAL::pca_estimate_normals<Concurrency_tag>(points.begin(), points.end(),
                               CGAL::First_of_pair_property_map<PointVectorPair>(),
                               CGAL::Second_of_pair_property_map<PointVectorPair>(),
                               nb_neighbors);

    // orient normals
    std::list<PointVectorPair>::iterator unoriented_points_begin =
        CGAL::mst_orient_normals(points.begin(), points.end(),
                                 CGAL::First_of_pair_property_map<PointVectorPair>(),
                                 CGAL::Second_of_pair_property_map<PointVectorPair>(),
                                 nb_neighbors);

    std::ofstream output(output_filename.c_str());
    if (!output ||
        !CGAL::write_xyz_points_and_normals(output
                                         , points.begin(), points.end()
                                         , CGAL::First_of_pair_property_map<PointVectorPair>()
                                         , CGAL::Second_of_pair_property_map<PointVectorPair>()))
    {
        std::cerr << "Error: cannot write " << output_filename << std::endl;
        return EXIT_FAILURE;
    }<img src="http://f.st-hatena.com/images/fotolife/k/kuyata/20160422/20160422014516.png" width=350>

    return EXIT_SUCCESS;
}

コンパイルclang++ -std=c++11 -I/home/vagrant/cgal/cgal/usr/local/include/ -I/usr/include/eigen3/ -I/usr/include/ -lCGAL -lboost_thread -lboost_system -lgmp -lmpfr

  • rglで可視化する
camel <- read.table("outputcamel.xyz")
body <- 50 * camel[,1:3] #(適当に拡大する)
normal <- camel[,4:6]

plot3d(rgl)
# 9770頂点の内900頂点に法線をプロットする 
for (i in 1:900) {
    ori <- body[10*i,]
    iro <- body[10*i,] + normal[10*i,]
    x <- rbind(ori,iro)
    segments3d(as.matrix(x), add=TRUE, c='blue')
}

この段階では形は保存されているように見えるが、足の部分の法線は方向がバラバラになっている。

  • 法線だけプロットしたもの

poisson reconstruction

CGALのサンプルを実行する - tak0kadaの何でもノートと同様にした

結果

またしてもうまく行かず、後ろ足が切れてしまった。

分析

元データcamel.offには法線データがないのでkitten.xyzとoni.xyzで分析してみる。

# データ読み込み、ctrlが元データ、testが推定データ
ctrl <- read.table("ctrl.xyz")<img src="http://f.st-hatena.com/images/fotolife/k/kuyata/20160422/20160422014516.png" width=350>
test <- read.table("test.xyz")

# 座標を比較する
# TRUE
any(ctrl[, 1:3] == test[, 1:3])

# 法線のずれを調べる
diff_norm <- apply(test[,4:6] - ctrl[,4:6] ,1 ,function(x){x[1]^2 + x[2]^2 + x[3]^2})

# tagcloudを使ってグラデーションをつけることができる
# colorは差が小さと1→青、大きいと0→赤に対応(分かりにくいが)
library(tagcloud)
color <- (max(diff_norm) - diff_norm) / (max(diff_norm) - min(diff_norm))
plot3d(ctrl[,1:3], col = smoothPalette(color, pal="RdBu"))

比較的形が保存されているように見えるネコでもメッシュが異なる。鬼は口と角でメッシュの違いが大きい。青と赤の点で表示したのは法線の差の絶対値の2乗で、ネコはほぼ均一、鬼はメッシュの違いが表れた口と角でズレが大きい。元のデータで相当数の頂点が含まれているがメッシュ化後では頂点数がかなり減少している。

法線の差をプロットすると、凹凸になっている部分で外側への成分が過小評価されているのが分かる。