【Python】点集合の地図・空間を作る【numpy実践例】

# 点集合データによる地図・空間の構築

※前回の補足的な内容です

前回の記事で、ラベリングは「距離の近さ」を使って、画像のラベリングを行なった。

https://qiita.com/belre/items/7570051d793049ed8734

けれども、例えば無線の基地局がたくさんある場合など、画素があらゆる基準点から影響を与える場合、前回の内容だと厳密性に欠けるので、いろんな基準点を考えた時の地図や空間を作りたい。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/281173/5d2f202c-8b3d-1c12-35af-584528af88f5.png)

イメージとしては、基準点1～3から無線の電波が飛んでいて、それが距離に依存するモデルである場合に、どんな信号空間になるかを検討する、といった使い方をする。

# 前回のプログラムからの変更点

## 参照点に対する評価値を事前に計算する

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/281173/a39fe879-2a14-be15-a5bc-03596acd4bed.png)

dist_matという、基準点に対する評価値のデータを最初に生成しておきます。この評価値は、前回は例えばL1ノルム, L2ノルム、マフィンティンポテンシャルを設定しましたが、今回は各点に指数分布による減衰を表す空間としてi番目の基準点では以下の信号が放たれているものとします。

```math
f_i(x, y) = \exp(- \frac{\sqrt{(x-x_i)^2+(y-y_i)^2}}{\lambda}) - a_0
```

## 情報落ちの防止

　微妙な数値を足し算、引き算する場合には、この問題が必ず付きまといます。例えば、0.001を1000回加算した後、1を加算するものとすると、有効桁数が２桁だったとしますと、1から先に計算した場合

```math
1 + 0.001 + 0.001 + ... + 0.001 = 1
```

という結果になってしまいますが、0.001を1000回加算するときから考えると

```math
0.001 + 0.001 + ... + 0.001 + 1 = 2
```

ときちんと計算されるわけです。

　現代ではあまり意識することは少なそうに見えますが、numpyのfloat64を含む浮動小数点にはIEEEで決められている数値の表現があり、float64やもっと高精度なfloat128, float256でも許容される有効桁数が１０進数で定義されているわけです（これはpythonに限らず、浮動小数点が仕様化されている言語では同様の問題を検討する必要があります）。

　したがって、加算回数が不定で増大する可能性がある場合には、有効桁数を必ず意識して計算をしなければ、本来想定した有効桁よりも信頼できない値が返される可能性が出てきます。

### 多次元配列版のnp.sortを使ってソートをしてから加算する


考え方としては、定義されたdist_matを＋、-に分割して分けた後、その数値をnp.sortを使って並べてあげれば良いです。これによって、dist_mat_p、dist_mat_mの絶対値の小さな値が先に計算されるようになります。そして、最後にdist_mat_pとdist_mat_mを全部足し合わせます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/281173/09b2e2a3-5fe6-10bd-3f49-8b20875ffff8.png)


ソースコードのその部分だけを切り取ると、こんな感じです。

```Python
dist_mat_p = np.sort(np.multiply((dist_mat > 0).astype(np.int64), dist_mat), axis=0)
dist_mat_m = np.sort(np.multiply((dist_mat < 0).astype(np.int64), dist_mat), axis=0)[::-1]
```

　まず、dist_mat>0とすると、PythonではTrue, Falseが出てくるので、これをintに変換して1/0の形式に変換します。その後、dist_mapを行列積で計算すると、正の成分と負の成分をくくりだすことができます。

　最後に、np.sortを使ってソートをしますが、この時axis=0を指定します。添字の２番目、３番目はいわゆる２次元データを管理していて、１番目が基準点の添字を表しているので、１番目の基準点でソートするのが、今回は適切となります。

# プログラム

## 実行結果

### ランダムな点の場合

![スクリーンショット 2021-03-13 23.20.15.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/281173/e3b66f4f-8751-31fa-b414-09599cc561bf.png)

### 四角形の格子形状の場合

![スクリーンショット 2021-03-13 23.21.17.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/281173/690c6311-c63a-bd4c-d533-1997178b1f74.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_all_elements(shape, ref_points, **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_all_elements):
            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

    # onesは前もって用意する
    ones = np.ones((height, width))

    # 参照点に対して評価値を事前に算出
    dist_mat = np.zeros((ref_points.shape[0], coord_mat.shape[0], coord_mat.shape[1]))
    for k in range(1, ref_points.shape[0]):
        dist_mat[k] = eval_func(ref_points[k], coord_mat[:, :])

    # 情報落ちを避けるためソートを行う
    dist_mat_p = np.sort(np.multiply((dist_mat > 0).astype(np.int64), dist_mat), axis=0)
    dist_mat_m = np.sort(np.multiply((dist_mat < 0).astype(np.int64), dist_mat), axis=0)[::-1]

    # データ点に対して計算
    #+
    allpix_mat_p = np.zeros((height, width))
    for k in range(0, ref_points.shape[0]):
        allpix_mat_p[:,:] += dist_mat_p[k, :, :]

    #-
    allpix_mat_m = np.zeros((height, width))
    for k in range(0, ref_points.shape[0]):
        allpix_mat_m[:,:] += dist_mat_m[k, :, :]

    # 総和
    allpix_mat = allpix_mat_p + allpix_mat_m
    
    m_obj.measure_time()

    return allpix_mat

# 幅, 高さ
width = 500
height = 500

# 参照点情報

random_num = 100
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()))

# 指数ポテンシャル
def exponent(r, x):
    dist0 = 50.0
    dist = np.linalg.norm(r - x, axis=2)
    eval = np.exp(-dist/dist0)
    return eval

plt.figure()
dist_mat = fig_all_elements((width, height), ref_points, eval_func=exponent)
plt.imshow(dist_mat, cmap='gray')
for ref in ref_points:
    plt.scatter(ref[0], ref[1], c='red')
plt.show()

# 四角格子
ref_width = 10
ref_height = 10
ref_points_2 = np.zeros((ref_width * ref_height, 2))
cnt = 0
for j in range(0, ref_height):
    for i in range(0, ref_width):
        ref_points_2[cnt] = np.array( ((i + 0.5) * width / ref_width, (j + 0.5) * height / ref_height))
        cnt += 1

plt.figure()
dist_mat = fig_all_elements((width, height), ref_points_2, eval_func=exponent)
plt.imshow(dist_mat, cmap='gray')
for ref in ref_points_2:
    plt.scatter(ref[0], ref[1], c='red')
plt.show()

```