【MATLAB】3次元点群の局所特徴量（固有値ベースのいくつか）

## 概要

この記事では、MATLABを使用して3次元点群の局所特徴量を可視化し、その物理的な解釈を実際の地物に対応付けて確認します。対象とする特徴量は、固有値ベースのLinearity、Planarity、Scattering、Eigenentropyです。これらの特徴量は、3次元点群の幾何学的特徴をつかむのに効果的です。

## 使用環境

- **使用言語**: MATLAB
- **バージョン**: R2025a Prerelease
- **必要なツールボックス**: Computer Vision Toolbox

## 固有値ベースの局所特徴量

以下の特徴量を対象にしました。他にも固有値をベースにした特徴量はありますが、役割が似通っていたり解釈性が良くない（というか自分にとって理解しづらい）ものは除きました。

1. **Linearity**
   - **算出式**: $$\text{Linearity} = \frac{e1 - e2}{e1}$$
   - **解釈**: 点群が直線状に分布しているかを示し、値が1に近いほど直線的です。


2. **Planarity**
   - **算出式**: $$ \text{Planarity} = \frac{e2 - e3}{e1} $$
   - **解釈**: 点群が平面上に分布しているかを示し、値が1に近いほど平面的です。


3. **Scattering**
   - **算出式**: $$ \text{Scattering} = \frac{e3}{e1} $$
   - **解釈**: 点群が空間的に散らばっているかを示し、値が1に近いほど散らばっています。


4. **Eigenentropy**
   - **算出式**: $$\text{Eigenentropy} = - (e1 \cdot \log(e1) + e2 \cdot \log(e2) + e3 \cdot \log(e3)) $$
   - **解釈**: 点群の情報エントロピーを示し、値が大きいほど情報が多様です。

これらの特徴量は各点においてそれらの近傍点から計算され、値域は全て0－1です。

## MATLABコード
[先の曲率を計算したときの記事](https://qiita.com/shun-kusano/items/ec211baaec66f1d0aaa8)で、eigValsとなっているのが各点の固有値（3x1xNの配列）です。一行目が最小固有値に対応しています。このeigValsをスタートポイントとして、下記のコードで、上記4つの特徴量の計算ができます。
```matlab
eigVals = squeeze(eigVals); % 3xNに

% extract eigenvalues
e1 = eigVals(3,:);
e2 = eigVals(2,:);
e3 = eigVals(1,:);
    
% Calculate features
% Calculate and store the linearity feature
features.linearity = (e1 - e2) ./ e1; % Linearity

% Calculate and store the planarity feature
features.planarity = (e2 - e3) ./ e1; % Planarity

% Calculate and store the scattering feature
features.scattering = e3 ./ e1; % Scattering

% Calculate and store the eigenentropy feature
sum_eig = e1+e2+e3;
ne1 = e1 ./ sum_eig; % 値域を[0 1]にするために正規化
ne2 = e2 ./ sum_eig;
ne3 = e3 ./ sum_eig;
eigenentropy = - (ne1 .* log(ne1) + ne2 .* log(ne2) + ne3 .* log(ne3)); % Eigenentropy
eigenentropy = max(0, min(1, eigenentropy)); % 念のためクリッピング（不要かも）
features.eigenentropy = eigenentropy;
```

## データセットの説明

この記事で使用するデータセットは、東京都が公表している点群データです。[G空間](https://www.geospatial.jp/ckan/dataset/tokyopc-23ku-2024)から会社のオフィスがある赤阪ガーデンシティ周辺を狙ってダウンロードしました。いったんいくつかの区画をダウンロードし結合したものから、めぼしい範囲で切り取ってます。

```matlab
pcshow(ptCloud)
```
![2025-02-07_14-54-40.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202697/44355377-95d0-507e-b892-9ecfc38cdbe0.gif)

データに含まれる点数は約383万点でした。
読み込んでから特徴量計算が終わるまで第12世代のi7-1270Pで16秒程度かかりました。


## 可視化の結果と解釈

MATLABを用いて、各特徴量を可視化した結果を以下に示します。各特徴量の大きさは下のカラーバーで色付けしています・・・と言えればよかったのですが、図を作った後で正確に0-1に合わせていないことに気づいたので、おおよそそうなっていると思ってみてください。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202697/02f56088-6eaa-264a-2e65-22fbbda4626c.png)

- **Linearity**: 高い値は直線的な構造を示していて、建物のエッジや歩道橋などを表しているようです。
![Linearity.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202697/7ab702cf-de9a-70ec-dc1f-980b6355b119.gif)



- **Planarity**: 高い値は平面的な構造を示していて、道路や建物の屋上を表しています。
![Planarity.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202697/ccbaf126-86b7-aeb7-783c-9743eaf99bfc.gif)



- **Scattering**: 高い値は散らばった構造を示し、基本体にはランダムな構造に反応しているので植生を表しますが、建物の中でも複雑な構造の部分には反応しているのだと思います。実際どのような構造に？というのは良く観察しないとわからないですね。
![Scattering.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202697/d6e0f511-e3c2-d9e5-ffdf-595e9dc86bc9.gif)



- **Eigenentropy**: 乱雑さを表すので基本的にはScatteringと似た傾向を示していると思います。ただLinearityとPlanarityの区別がこっちの方が付きやすいように見えます。面白いですね。

![eigenentropy.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202697/49ab64e6-4597-d086-472f-418522287788.gif)

## まとめ
固有値ベースの3次元点群の局所特徴量について計算、可視化して物理的な解釈と実際の地物がどのように対応するのか見てみました。実際の地物の形状は複雑でセオリー通りに解釈しきれないものも多いとは思いますが、おおよその傾向をつかむには使える特徴量だと思います。それぞれの大きさでもって機械学習にかけて地物の分類なども比較的とっつきやすいかなと思いました。

また、今回は近傍点の数を20にしましたが、特に理由は無いです。適応的にこのあたりのパラメータも変化させることでより安定した結果を得ることはできそうです。

博士の時似たようなこと（レーダ画像から地物を判読するための処理手法の研究）してたな…。