Open3Dを用いた3D処理入門｜点群処理、メッシュ化

Open3Dの概要とインストール

Open3Dの概要

Open3Dは3Dデータを取り扱うソフトウェアの開発にあたって用いられるオープンソースのライブラリです。OpenCVはC++やPythonから使用することができます。

Open3Dには下記のような特徴があります。

・3D data structures
・3D data processing algorithms
・Scene reconstruction
・Surface alignment
・3D visualization
・Physically based rendering (PBR)
・3D machine learning support with PyTorch and TensorFlow
・GPU acceleration for core 3D operations
・Available in C++ and Python

Open3Dのインストール

Open3Dは下記を実行することでPyPIから入手することができます。

$ pip install open3d

OpenCVの実行例

当節では以下、上記のチュートリアルを元にOpen3Dを用いた点群処理について確認します。基本的にチュートリアルと同じコードを用いていますが、ライブラリの読み込みや重複するコードなどが省略されているのでそのまま動くようにコードを追記しました。

Open3Dを用いた点群処理

以下、上記に基づいてOpen3Dを用いた点群処理のサンプルコードについて確認します。

import numpy as np
import open3d as o3d


ply_point_cloud = o3d.data.PLYPointCloud()
pcd = o3d.io.read_point_cloud(ply_point_cloud.path)

print(pcd)
print(np.asarray(pcd.points))
print(np.asarray(pcd.points).shape)

o3d.visualization.draw_geometries([pcd],
                                  zoom=0.3412,
                                  front=[0.4257, -0.2125, -0.8795],
                                  lookat=[2.6172, 2.0475, 1.532],
                                  up=[-0.0694, -0.9768, 0.2024])

・実行結果

PointCloud with 196133 points.
[[0.65234375 0.84686458 2.37890625]
 [0.65234375 0.83984375 2.38430572]
 [0.66737998 0.83984375 2.37890625]
 ...
 [2.00839925 2.39453125 1.88671875]
 [2.00390625 2.39488506 1.88671875]
 [2.00390625 2.39453125 1.88793314]]
 (196133, 3)

実行したスクリプトは「o3d.data.PLYPointCloudやo3d.io.read_point_cloudを元にpcdに保存したオブジェクトが196,133個の点を持ち、np.asarray(pcd.points)によってNumPy形式で取得できる」と大まかに理解すると良いと思います。

また、o3d.io.read_point_cloudは「ファイルからの点群読み込み」、o3d.visualization.draw_geometriesは「点群の可視化」に用いられるメソッドであることも合わせて抑えておくと良いです。

Open3Dを用いたメッシュ化

以下、上記に基づいてOpen3Dを用いたメッシュ化のサンプルコードについて確認します。

import numpy as np
import open3d as o3d


knot_mesh = o3d.data.KnotMesh()
mesh = o3d.io.read_triangle_mesh(knot_mesh.path)

print(mesh)

print('Vertices:')
print(np.asarray(mesh.vertices))
print(np.asarray(mesh.vertices).shape)

print('Triangles:')
print(np.asarray(mesh.triangles))
print(np.asarray(mesh.triangles).shape)

o3d.visualization.draw_geometries([mesh])

・実行結果

TriangleMesh with 1440 points and 2880 triangles.
Vertices:
[[  4.51268387  28.68865967 -76.55680847]
 [  7.63622284  35.52046967 -69.78063965]
 [  6.21986008  44.22465134 -64.82303619]
 ...
 [-22.12651634  31.28466606 -87.37570953]
 [-13.91188431  25.4865818  -86.25827026]
 [ -5.27768707  23.36245346 -81.43279266]]
(1440, 3)
Triangles:
[[   0   12   13]
 [   0   13    1]
 [   1   13   14]
 ...
 [1438   11 1439]
 [1439   11    0]
 [1439    0 1428]]
(2880, 3)

また、上記のコードで可視化を行う前に下記のようにmesh.compute_vertex_normals()を実行することで法線ベクトル(Surface Normal)の推定を行うことができることも合わせて抑えておくと良いです。

import numpy as np
import open3d as o3d


knot_mesh = o3d.data.KnotMesh()
mesh = o3d.io.read_triangle_mesh(knot_mesh.path)

mesh.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh])

・実行結果

外れ値の除去

SOR(Statistical Outlier Removal)

import numpy as np
import open3d as o3d


def display_inlier_outlier(cloud, ind):
    inlier_cloud = cloud.select_by_index(ind)
    outlier_cloud = cloud.select_by_index(ind, invert=True)

    print("Showing outliers (red) and inliers (gray): ")
    outlier_cloud.paint_uniform_color([1, 0, 0])
    inlier_cloud.paint_uniform_color([0.8, 0.8, 0.8])
    o3d.visualization.draw_geometries([inlier_cloud, outlier_cloud],
                                      zoom=0.3412,
                                      front=[0.4257, -0.2125, -0.8795],
                                      lookat=[2.6172, 2.0475, 1.532],
                                      up=[-0.0694, -0.9768, 0.2024])

                                      

print("Load a ply point cloud, print it, and render it")
sample_pcd_data = o3d.data.PCDPointCloud()
pcd = o3d.io.read_point_cloud(sample_pcd_data.path)
print(np.asarray(pcd.points).shape)
o3d.visualization.draw_geometries([pcd],
                                  zoom=0.3412,
                                  front=[0.4257, -0.2125, -0.8795],
                                  lookat=[2.6172, 2.0475, 1.532],
                                  up=[-0.0694, -0.9768, 0.2024])                       

voxel_down_pcd = pcd.voxel_down_sample(voxel_size=0.02)


print("Statistical oulier removal")
cl, ind = voxel_down_pcd.remove_statistical_outlier(nb_neighbors=20,
                                                    std_ratio=2.0)
display_inlier_outlier(voxel_down_pcd, ind)

・実行結果

voxel_down_pcd.remove_statistical_outlierによってSOR(Statistical Outlier Removal)が実行されていることを確認しておくと良いです。nb_neigborsは平均のdistanceを計算するにあたって用いる点の数、std_ratioは$\mu_i$の値を取り除く際に用いる標準偏差です。

SORのアルゴリズムについては下記で詳しくまとめましたので合わせて参照ください。

Radius Outlier Removal

import numpy as np
import open3d as o3d



def display_inlier_outlier(cloud, ind):
    inlier_cloud = cloud.select_by_index(ind)
    outlier_cloud = cloud.select_by_index(ind, invert=True)

    print("Showing outliers (red) and inliers (gray): ")
    outlier_cloud.paint_uniform_color([1, 0, 0])
    inlier_cloud.paint_uniform_color([0.8, 0.8, 0.8])
    o3d.visualization.draw_geometries([inlier_cloud, outlier_cloud],
                                      zoom=0.3412,
                                      front=[0.4257, -0.2125, -0.8795],
                                      lookat=[2.6172, 2.0475, 1.532],
                                      up=[-0.0694, -0.9768, 0.2024])

                                      
print("Load a ply point cloud, print it, and render it")
sample_pcd_data = o3d.data.PCDPointCloud()
pcd = o3d.io.read_point_cloud(sample_pcd_data.path)
print(np.asarray(pcd.points).shape)
o3d.visualization.draw_geometries([pcd],
                                  zoom=0.3412,
                                  front=[0.4257, -0.2125, -0.8795],
                                  lookat=[2.6172, 2.0475, 1.532],
                                  up=[-0.0694, -0.9768, 0.2024])

voxel_down_pcd = pcd.voxel_down_sample(voxel_size=0.02)


print("Radius oulier removal")
cl, ind = voxel_down_pcd.remove_radius_outlier(nb_points=16, radius=0.05)
display_inlier_outlier(voxel_down_pcd, ind)

・実行結果

voxel_down_pcd.remove_radius_outlierによってRadius Outlier Removalが実行されていることを確認しておくと良いです。それぞれの点の周りのradiusで指定した半径の中にnb_pointsで指定した点の数が存在しない場合は外れ値と見なし取り除きます。

SORとRadius Outlier Removalの同時実行

SORとRadius Outlier Removalのそれぞれの挙動を確認するにあたって、下記のプログラムを作成しました。

import numpy as np
import open3d as o3d


def display_inlier_outlier(cloud, ind1, ind2):
    inlier_cloud = cloud.select_by_index(np.intersect1d(ind1, ind2))
    outlier_cloud = cloud.select_by_index(np.union1d(ind1, ind2), invert=True)
    outlier_cloud1 = cloud.select_by_index(np.setdiff1d(ind1, ind2)) # Radius Outlier Removal
    outlier_cloud2 = cloud.select_by_index(np.setdiff1d(ind2, ind1)) # SOR
    print(inlier_cloud)
    print("RadiusOutlierRemoval_only: ", outlier_cloud1)
    print("SOR_only: ", outlier_cloud2)
    print("Union_OutlierRemoval: ",outlier_cloud)

    print("Showing outliers (red) and inliers (gray): ")
    outlier_cloud1.paint_uniform_color([1, 0, 0])  # Radius Outlier Removal: red
    outlier_cloud2.paint_uniform_color([0, 1, 0])  # SOR_only: green
    outlier_cloud.paint_uniform_color([0, 0, 1])  # Union_OutlierRemoval: blue
    inlier_cloud.paint_uniform_color([0.8, 0.8, 0.8])
    o3d.visualization.draw_geometries([inlier_cloud, outlier_cloud, outlier_cloud1, outlier_cloud2],
                                      zoom=0.3412,
                                      front=[0.4257, -0.2125, -0.8795],
                                      lookat=[2.6172, 2.0475, 1.532],
                                      up=[-0.0694, -0.9768, 0.2024])


ply_point_cloud = o3d.data.PLYPointCloud()
pcd = o3d.io.read_point_cloud(ply_point_cloud.path)

print(pcd)

voxel_down_pcd = pcd.voxel_down_sample(voxel_size=0.02)
print(voxel_down_pcd)


print("Oulier Removal")
cl, ind1 = voxel_down_pcd.remove_statistical_outlier(nb_neighbors=20, std_ratio=2.0)
cl, ind2 = voxel_down_pcd.remove_radius_outlier(nb_points=16, radius=0.05)
print(len(ind1))
print(len(ind2))
display_inlier_outlier(voxel_down_pcd, np.array(ind1), np.array(ind2))

・実行結果

PointCloud with 196133 points.
PointCloud with 27012 points.
Oulier Removal
26136
25696
PointCloud with 25675 points.
RadiusOutlierRemoval_only:  PointCloud with 461 points.
SOR_only:  PointCloud with 21 points.
Union_OutlierRemoval:  PointCloud with 855 points.

SORのみで取り除かれるサンプルを緑、Radius Outlier Removalのみで取り除かれる点を赤、双方の手法で取り除かれるサンプルを青で表示しました。

点群オブジェクトの読み込み、作成、保存

以下、PUNetの10000_poissonのduck.xyz(スペース区切りのテキストファイル)のply化とplyファイルの読み込みについて確認します。

import numpy as np
import open3d as o3d

points = np.loadtxt("data/duck.xyz", delimiter=" ")
print(points.shape)

pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(points)
pcd.colors = o3d.utility.Vector3dVector(np.ones_like(points))
print(pcd)
o3d.io.write_point_cloud("data/duck.ply", pcd)

・実行結果

(10000, 3)
PointCloud with 10000 points.

上記を実行することでduck.xyzに基づいてduck.plyを作成することができます。o3d.geometry.PointCloudクラスに基づいて作成したpcdオブジェクトに対し、pcd.pointsで点の位置、pcd.colorsで点の色を指定している点に着目しておくと良いです。plyファイルの読み込みについては下記のようなプログラムを実行すれば良いです。

import open3d as o3d

pcd = o3d.io.read_point_cloud("data/duck.ply")
print(pcd)

・実行結果

PointCloud with 10000 points.

また、作成したplyファイルはmeshlabなどを用いることで可視化することができます。

$ meshlab data/duck.ply

・実行結果

meshlabのインストール

前項で使用したmeshlabは下記のようなコマンドでインストールすることができます(Ubuntu環境)。

$ apt install meshlab