点群から代表点を抽出する

# 0. 概要
点群が得られたとする。この点群から代表的な点はどこかを探したい。よくある方法は平均値を取る方法である。しかしこの方法は外れ値に弱い。

また、座標として不連続の場合あり得ない座標（ゴルフコースでいう山林的な場所）を示す場合がある。そこで、全点群の距離を計算し、その距離が最も小さい点群を選ぼうとする。この場合、$O = N^2$となりとてつもない計算量になる。よってこの手法は現実的ではない。

では、どうしたらよいか？それをインターネットで調査したものの全く見つからなかった。というか検索ワードが分からなかった。

その後、5分くらい考えたすえ以下の方法を思いついた。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/95636/763f8996-63cc-9591-a67e-0d454cf425c4.png)

アルゴリズムも簡単で、以下のものを実装する。
1) 全点群をグリッドで分ける
2) グリッド内で点群が多いものを抽出
3) 抽出した点群でまた1)を繰り返す
4) グリッド内に含まれる点群が全て1以下になったら終了条件である。
5) 複数グリッドに1つの点群が含まれる場合は、全点群の平均位置から近い点群を代表点とする

# 1. 実装
## 全探索

最も安定性のある解を出すが、計算量が$O = N^2$になってしまうプログラムである。

``` python
def repPoint_full_search(_x, _y):
  data_num = len(_x)
  distance = np.zeros(data_num)
  for i in range(data_num):
    for j in range(data_num):
      distance[i] += np.abs(_x[i] - _x[j]) + np.abs(_y[i] - _y[j])
  return distance

x = np.random.rand(1000)
y = np.random.rand(1000)
plt.figure(figsize=(12, 10))
plt.scatter(x, y)
dis = repPoint_full_search(x, y)
num = np.argmin(dis)
plt.scatter([x[num]], [y[num]], c='pink')
plt.xlabel("z[0]")
plt.ylabel("z[1]")
plt.savefig("repPoint.png")
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/95636/85df95d0-d758-1284-24d6-5b79b985cad8.png)

1000x1000程度なので20秒程度で終わる。しかし、5000x5000とかになるともう終わらない。
計算量は`25, 000, 000`にもなり、`1ms`で1サイクルの計算が終わるとしても`6.94h`の計算がかかる

## グリッドサーチ法
結果から説明する。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/95636/67d89ace-ee93-b37c-bda6-8770774a5570.png)

赤色がグリッドサーチによる点群の代表点である。ピンク色が全探索の結果の代表点である。
残念ながら同じ点群を指すことは出来ていないが、少ない計算量で近しいところを見つけてこれている。

次に、処理手順を説明する。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/95636/869aac84-b8f7-5cef-0e86-3b31f8e24ca3.png)

全点群データから密集度が最も高いグリッドを見つける。ここだとピンク色である。
因みに、緑色の点が平均位置である。なお、この後はこれを繰り返していく
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/95636/d5426dd8-aa08-941d-53d4-49ad5413662a.png)
2エポック目

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/95636/f8e80d48-17ef-1b67-19b3-b6c376397234.png)
3エポック目
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/95636/32b513fa-bee9-dab5-d942-2f844bbfa0d7.png)
4エポック目

ここで点群が1つに絞られ、1つの結果が返ってきている

4x4のグリッドで探索を行った途中結果を以下に示す。

```
[2. 4. 5. 2.]
[ 7. 12.  9.  1.]
[ 7. 22. 14.  5.]
[0. 4. 6. 0.]
[5]
0.7743798569440619
Minimum grid:
5
--------------------
[1. 3. 6. 1.]
[0. 0. 0. 2.]
[2. 0. 2. 2.]
[1. 2. 0. 0.]
[14]
0.3999928541358709
Minimum grid:
14
--------------------
[0. 0. 0. 1.]
[0. 0. 0. 1.]
[0. 0. 2. 0.]
[1. 1. 0. 0.]
[6]
0.39153543531607427
Minimum grid:
6
--------------------
[0. 0. 0. 1.]
[0. 0. 0. 0.]
[0. 0. 0. 0.]
[1. 0. 0. 0.]
[0, 15]
0.4133443602306934
0.3697265104014552
Minimum grid:
15
--------------------
```

最終的にグリッドの対角のみが残り、平均位置に最も近い方が採用される仕組みとなっている

コードを以下に示す。

``` python
def grid_search(_x, _y, _mx, _my, _grid_num = 4):
  data_num = len(_x)
  distance = np.zeros(data_num)
  ##
  tmp_x = _x
  tmp_y = _y
  center_x_all_data = _mx
  center_y_all_data = _my
  eps = 0.00000001
  # can optimize yet
  min_x = np.min(tmp_x) - eps
  max_x = np.max(tmp_x) + eps
  min_y = np.min(tmp_y) - eps
  max_y = np.max(tmp_y) + eps
  split_num = _grid_num
  step_x = (max_x - min_x) / split_num
  step_y = (max_y - min_y) / split_num
  #
  grid_count = np.zeros(split_num * split_num)
  grid_center_x = np.zeros(split_num * split_num)
  grid_center_y = np.zeros(split_num * split_num)
  
  # phase 1: Deviding grids
  for i in range(data_num):
    x_num = math.floor ( ( (tmp_x[i] - min_x) / step_x) )
    y_num = math.floor( ( (tmp_y[i] - min_y) / step_y) )
    grid_num = int(x_num + ( y_num * split_num )) 
    grid_count[grid_num] += 1
    grid_center_x[grid_num] += tmp_x[i]
    grid_center_y[grid_num] += tmp_y[i]

  for i in range(_grid_num):
    j = (_grid_num - 1) - i
    print(grid_count[j*_grid_num:(j*_grid_num)+_grid_num])

  # phase 2: finding density grid
  maxIndex = [i for i, x in enumerate(grid_count) if x == max(grid_count)]
  print(maxIndex)

  tmp_cost = 99999999
  tmp_idx = -1
  for i in range(len(maxIndex)):
    idx = maxIndex[i]  
    cx = grid_center_x[idx] / grid_count[idx]
    cy = grid_center_y[idx] / grid_count[idx]
    r = np.abs(cx - center_x_all_data) + np.abs(cy - center_y_all_data)
    print(r)
    if r < tmp_cost:
      tmp_cost = r
      tmp_idx = i
  
  print("Minimum grid:")  
  print(maxIndex[tmp_idx])
  print("--------------------")
  gird_idx = maxIndex[tmp_idx]

  # phase : choosing point data
  ext_tmp_x = []
  ext_tmp_y = []
  for i in range(data_num):
    x_num = math.floor ( ( (tmp_x[i] - min_x) / step_x) )
    y_num = math.floor( ( (tmp_y[i] - min_y) / step_y) )
    grid_num = int(x_num + ( y_num * split_num )) 
    if(grid_num == gird_idx):
      ext_tmp_x.append(tmp_x[i])
      ext_tmp_y.append(tmp_y[i])

  return (ext_tmp_x, ext_tmp_y)


x = np.random.normal(0, 1, 100)
y = np.random.normal(0, 1, 100)
ex = x
ey = y
mx = np.average(x)
my = np.average(y)
carray = ['blue', 'pink', 'red', 'yellow', 'black', 'blue']
for i in range(5):
  plt.clf()
  plt.figure(figsize=(12, 10))
  plt.scatter(ex, ey, c=carray[i])
  plt.scatter([mx], [my], c='green')
  plt.xlabel("z[0]")
  plt.ylabel("z[1]")
  plt.ylim(np.min(ey), np.max(ey) )
  plt.xlim(np.min(ex), np.max(ex) )
  plt.xticks(np.linspace(np.min(ex), np.max(ex), 5))
  plt.yticks(np.linspace(np.min(ey), np.max(ey), 5))
  ex, ey = grid_search(ex, ey, mx, my)
  plt.scatter(ex, ey, c=carray[i+1])
  plt.grid()
  plt.savefig("./gird/repPoint" + str(i) +".png")
  if len(ex) == 1:
    break

# all data
plt.clf()
plt.figure(figsize=(12, 10))
plt.scatter(x, y)
dis = repPoint_full_search(x, y)
num = np.argmin(dis)
plt.scatter([x[num]], [y[num]], c='pink')
plt.scatter([ex], [ey], c='red')
plt.xlabel("z[0]")
plt.ylabel("z[1]")
plt.ylim(np.min(y), np.max(y) )
plt.xlim(np.min(x), np.max(x) )
plt.xticks(np.linspace(np.min(x), np.max(x), 5))
plt.yticks(np.linspace(np.min(y), np.max(y), 5))
plt.grid()
plt.savefig("./gird/repPoint.png")
```

### 100, 000 x 100, 000
なお、この方法だと`100, 000`の点群でも1秒もかからずに処理が終わる。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/95636/94b1a807-869d-f119-9203-f7b74ac679c5.png)
Epoch 1

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/95636/788982a4-05e2-ceae-2870-77cdcbd6a40f.png)
Epoch 2

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/95636/f034b746-d3c1-8008-ffa9-9888721260c5.png)
Epoch 3

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/95636/e6e7506a-4b24-3601-23d4-ee32b830318e.png)
Epoch 4

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/95636/53509f57-c148-9ac1-e244-40b7a3ebc25b.png)
Epoch 5

このように10万個の点群、すなわち愚直にやれば`10, 000, 000, 000`サイクルを要する場合でも
`10万×5回(×3)`程度のサイクルで終わる。約7000倍の効率化である。