【Python】点要素の近くの画素をラベリングする【numpy実践例】


# 画素ラベリング(点要素周辺の空間解析)

Pythonでの実装演習
割と初歩的ではあるが便利

## 概要

　基準点1, ..., Nがあって、画素が基準点からの距離に依るときに、どんな空間になっているかを解析する。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/281173/3b782a91-1309-6a53-b92e-19e0b6864f3d.png)

## 使用例

　主に、画像処理や信号処理、シミュレーションなどの分野で使用。「なんかデータが空間的に偏ってるな」と感じるときに使うと、解析がはかどることがある。

- 二等分線図

　二点間の距離$ d_i, d_j $とすると、二等分線図の境界線であれば、明らかに$ d_i = d_j $になる。この性質を使って二等分線図を描くことが出来る。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/281173/38e387cc-cce7-c94c-8f5d-2c56f83714fc.png)

- 無線基地局の分布の解析
- 分子内に存在するエネルギー分布の解析

　上記の現象は、大体の場合は極座標または円柱座標系であり、距離rに依存する量である。信号強度やエネルギー分布を定義できる場合が多いので、出てきたデータを分析するのには都合が良い。
　（※今回は、便利さや処理時間の関係上、2次元でざっくり計算しますが、3次元でも同じです)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/281173/0df36366-3a05-567b-8ac7-c3e51f7128db.png)

## 考え方

### 画素参照

　こういう解析をするときに、ついうっかりやりがちなのが画素間での比較である。例えば、以下のように書きがちだが、Pythonは1画素ずつの参照の処理には時間がかかり、相当な時間がかかることが予想される。

```Python
# 幅, 高さ
width = 500
height = 500

# 参照点情報
random_num = 100
ref_points = np.zeros((random_num, 2))

# 解析
for j in range(0, height):
    for i in range(0, width):
        ...(何か処理を書く)...

```

　解決策として、配列計算はPythonのnumpyの行列演算ライブラリに全てお任せし、並列処理などがなるべくかけられやすいように条件分岐等を入れない形式で表現するのがポイント。

```Python
# 幅, 高さ
width = 500
height = 500

# 参照点情報
random_num = 100
ref_points = np.zeros((random_num, 2))

# 解析
for k in range(0, ref_points):
    ...(何か処理を書く)...

```


### 計算用の行列


* 基準データの反復の度に、minargs_mat(ラベリングデータ)とmin_mat(距離データ)を更新する。ここには、それぞれ最小値のデータが入るようにする。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/281173/d476d621-c8f4-2146-d0c9-428682bfdfd6.png)


* データ距離が同じ場合は、この方式では重なっている領域を配列で管理出来ない。今回は、添え字が最も低い値のデータを参照するか(Priority=True)、どちらなのか分からないとするか(Priority=False)を引数に渡して出力する。

* Priority = Falseになるケースを考慮して、基準データの先頭にダミーデータを入れておき、最後にminargsを-1差し引く。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/281173/d3affb74-69ff-4df2-6904-a9b0d6741290.png)

* coord_mat(座標格納データ)を予め生成する。中身は、x,yの座標が書かれる。

* coord_matを関数によって一つの評価値に変換する(tmp_mat).

* min_matとtmp_matの大小関係を~~比較~~・・・、ではなくて、numpy.signまたはnumpy.signbitによって、1か-1か0で表すようにする。以下、sig_mat, sig_mat_pmの二つを定義する。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/281173/1d032ade-4b08-8b37-3975-bbf977cce258.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/281173/568276ef-bee8-03fb-4092-37254f6a23bf.png)


* 状態遷移図で表すと、次の通り

| sig_mat (a) | Priority | sig_mat_pm(a') | 
|---|-----|---|
|1  |True |1  |
|-1 |True |0  |
|1  |False|1  |
|0  |False|1/2|
|-1 |False|0  |

　前状態のminargs_matをx, k_matをy、新しいminargs_matをuとする。まず、Priority=Trueの場合は、sig_matの変換を数式で行なった後、

```math
u = (1 - a') \times x + a' y 
```

とすれば、変換可能である。Priority=Falseの場合は、上記に加え

```math
v = u - \frac{x+y}{2} \times (1-|a|)
```

として、vを結果とすればよい。値が同じ場合に、a'=1/2なので、その1/2を差し引く演算をvでやる。

### 補足

　離散データをいきなり数式に書いて連続データとして取り扱う、という考え方に違和感を感じる人がいるかもしれない。テストがややこしそう、とかという声もあるだろう。

　そして、行列演算は親切にもこの場合は1、この場合は0、みたいな条件分岐を入れるとかえって重くなるので、離散値で全てを表現し、全て数式で書くやり方をした方が良くなる。こうすることで、高級言語を扱うときにボトルネックとなる「繰り返し動作」を短縮して処理速度を抑制することが出来る。

　・・・まあ、これは「計算無理」という状態から「なんとなく計算できそう」っていう状態に変えるときの話。処理時間が実用上問題になるケースでは、C++やfortranやRustなどと数値計算ライブラリを併せて直接演算するほうが良いだろう。

　その場合も、行列参照のどこがボトルネックになるかを考えることで、処理時間の高速化することは当然必要である。

# プログラム

## 実行結果

左側はラベリングの結果
右側は分布(黒が短く、白が長い)

### N=3

- 二等分線(=空間中のL2ノルムの最小値)

※Priority=True

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/281173/3b852bd4-4f94-9845-5d44-07af1b17bbc4.png)


- 空間中のL1ノルムの最小値

　最小化距離をL2ノルムではなくL1ノルムに変更した場合。
　途中のラベリングの位置が離散化されていることがわかる。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/281173/f47f220e-49e1-29af-7014-495889c5e6fa.png)




### N=100

- 二等分線(=空間中のL2ノルムの最小値)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/281173/33a07e72-cfb3-e653-c931-41bcf20d7052.png)



- 空間中のL1ノルムの最小値

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/281173/8f3b60a8-24cc-8813-f9c6-4989190f7e2e.png)



- マフィンティンモデル

マフィンティンモデルは基準点と対象とする点の距離が近いときのみを最小化する関数で、ポテンシャル・エネルギーの分布として代用されるモデル。先ほどと性質が違うように見えるが、このように距離rに依存する点の分布も解析可能。

距離閾値は固定で20を設定. 

Priorityは、各点のダブりを-1に設定して、再計算することを想定。ダブっているところと、粒子が存在しないところは-1に設定している。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/281173/e7d3adf5-e5f6-cdfd-343f-11676751a382.png)




## ソースコード

```Python


import sys
import numpy as np
import matplotlib.pyplot as plt
import time

# 時間計測クラス
class measure:
    start = time.time()
    isdisplay = True

    def __init__(self):
        self.start = time.time()

    def measure_time(self):
        elapsed_time = time.time() - self.start
        if self.isdisplay == True:
            print ("elapsed_time:{0}".format(elapsed_time) + "[sec]")
        self.start = time.time()

# 位置によって最小の値になるような画像情報を生成する
def fig_map(shape, ref_points, priority=True, **kwargs):
    m_obj=measure()

    width = shape[0]
    height = shape[1]

    eval_func = lambda r, x : np.linalg.norm(r - x, ord=2, axis=2)
    if 'eval_func' in kwargs:
        if type(kwargs['eval_func']) is type(fig_map):
            eval_func = kwargs['eval_func']

    # 距離最小値となる画素番号を抽出
    coord_mat = np.array([[ [i, j] for i in range(0, width)] for j in range(0, height)])
    minargs_mat = np.zeros((height, width), dtype='int64')
    min_mat = np.ones((height, width)) * sys.float_info.max

    # index=0の点を最初に追加する
    # 追加した情報は、priority==Trueの場合-1に、それ以外は>=0に置き換えるようにする.
    ref_points_add0 = np.vstack( (np.array((0, 0)), ref_points))
    
    # onesは前もって用意する
    ones = np.ones((height, width))

    # データ点に対して計算
    for k in range(1, ref_points_add0.shape[0]):
        # ref番号の添え字の配列を生成
        k_mat = ones * k

        # 評価値を現在の座標点と参照点から演算
        tmp_mat = eval_func(ref_points_add0[k], coord_mat[:, :])

        # 符号を計算.
        # priorityがTrueの場合は0も+とする
        sig_mat = None 
        if priority == True:
            sig_mat = (-2 * (np.signbit(min_mat - tmp_mat)  - 1/2)).astype(np.int64)
        else:
            sig_mat = np.sign(min_mat - tmp_mat)

        # ±のみを考慮した値の算出 (※同じ場合は1/2が出力されるので、補正項を足す)
        sig_mat_pm = (sig_mat + 1) / 2

        # ラベリング画像
        if priority == True:
            minargs_mat = ((1 - sig_mat_pm) * minargs_mat + sig_mat_pm * k_mat)
        else:
            minargs_mat_tmp = ((1 - sig_mat_pm) * minargs_mat + sig_mat_pm * k_mat)
            minargs_mat = minargs_mat_tmp - ((minargs_mat + k_mat)/2) * (1-np.abs(sig_mat))

        minargs_mat = minargs_mat.astype(np.int64)

        # 距離画像
        min_mat = (1 - sig_mat_pm) * min_mat + sig_mat_pm * tmp_mat

    m_obj.measure_time()
    
    return min_mat, (minargs_mat-1)

# 幅, 高さ
width = 500
height = 500

# 参照点情報
random_num = 3
ref_points = np.zeros((random_num, 2))
for i in range(0, random_num):
    ref_points[i] = np.array((width * np.random.rand(), height * np.random.rand()))

# 二等分線図
plt.figure()
dist_mat, label_mat = fig_map((width, height), ref_points, True)
plt.subplot(121)
plt.imshow(label_mat, cmap='plasma')
plt.subplot(122)
plt.imshow(dist_mat, cmap='gray')
plt.show()

# L1ノルムで分割
plt.figure()
dist_mat, label_mat = fig_map((width, height), ref_points, True, eval_func=lambda r, x : np.linalg.norm(r - x, ord=1, axis=2))
plt.subplot(121)
plt.imshow(label_mat, cmap='plasma')
plt.subplot(122)
plt.imshow(dist_mat, cmap='gray')
plt.show()

# マフィンティン・ポテンシャル
def muffin(r, x):
    thresh = 20.0
    dist = np.linalg.norm(r - x, axis=2)
    dist = np.sign(dist - thresh * np.ones((x.shape[0], x.shape[1])))
    return dist

plt.figure()
dist_mat, label_mat = fig_map((width, height), ref_points, False, eval_func=muffin)
plt.subplot(121)
plt.imshow(label_mat)
plt.subplot(122)
plt.imshow(dist_mat, cmap='gray')
plt.show()


```