3次元点群の曲率を計算

この記事では

MATLABの点群処理系のToolboxに、曲率を計算する機能が無かったので簡易的なものを作成してみました。

環境

MATLAB R2024bもしくはR2024
Computer Vision Toolbox

曲率計算手法

まじめにやろうとすると、局所点群に対して球にフィッティングすることになるかと思いますが、簡易的な方法として、固有値をベースに下記の計算で曲率とするケースもあります。

手順

1. 点群データの取得
   点群データを取得し、各点の近傍点を探索します。近傍点の探索には$k$-近傍法や半径検索法が一般的に用いられます。

2. 共分散行列の構築
   各点の近傍点を用いて共分散行列$\mathbf{C}$を構築します。共分散行列は次のように計算されます：

   $$
   \mathbf{C} = \frac{1}{N} \sum_{i=1}^{N} (\mathbf{x}_i - \bar{\mathbf{x}})(\mathbf{x}_i - \bar{\mathbf{x}})^\top
   $$

   ここで、$\mathbf{x}_i$は$i$番目の近傍点の座標、$\bar{\mathbf{x}}$は近傍点の重心、$N$は近傍点の数です。

3. 固有値と固有ベクトルの計算
   共分散行列$\mathbf{C}$の固有値$\lambda_1, \lambda_2, \lambda_3$（$\lambda_1 \geq \lambda_2 \geq \lambda_3$）と対応する固有ベクトルを計算します。

4. 曲率の計算
   曲率は以下のように最小固有値$\lambda_3$を用いて計算されます。

   $$
   \text{曲率} = \frac{\lambda_3}{\lambda_1 + \lambda_2 + \lambda_3}
   $$

MATLABコード

下記のように書くことが可能です。

ptCloud = pcread('teapot.ply');
% 近傍計算に用いる点数
numNeighbors = 10;

% 無効点・重複点を除去
[ptCloud_val, isValid] = ptCloud.removeInvalidPoints;
locs = unique(ptCloud_val.Location,"rows");
ptCloud_val = pointCloud(locs);

% 近傍点のインデックスを取得
ind = ptCloud_val.multiQueryKNNSearchImpl(locs,numNeighbors); % R2024bまではこちら
% ind = findNearestNeighbors(ptCloud_val, locs, numNeighbors); % R2025aからはこちら

% 各点の近傍点座標の抽出
neighPts = reshape(locs(ind(:),:),numNeighbors,[],3);
neighPts = permute(neighPts,[1 3 2]);

% 共分散行列の計算
n_neighPts = neighPts - mean(neighPts,1);
covs = pagemtimes(permute(n_neighPts,[2 1 3]), n_neighPts) ./ (numNeighbors-1);

% 固有値計算
eigVals = pageeig(covs);

% 曲率計算
curvatures = squeeze(eigVals(1,1,:) ./ sum(eigVals,1)); % 一つ目が最小固有値になっている

結果の確認

曲率に応じて色を変えて表示します。

cmap = turbo(256)
pcshow(locs,cmap(round(curvatures*3*255+1),:),'MarkerSize',10);

今回出している曲率の最大値は$\lambda_3=\lambda_2=\lambda_1$となるときで、このとき曲率は1/3となります。カラースケールが0-255なので3*255を掛けています。

それっぽい結果になっています。値が小さいところでもっと変化付けてもいいかもしれませんｎ。。。

まとめ

3次元点群の曲率を直接計算する関数が無かったので、簡易的なものを作成してみました。点群の局所特徴量として固有値をベースにしたものはほかにもたくさんあるようなので、また紹介できればと思います。