Tellus画像処理の練習メモ：緑地を切り出してGeoJSONを生成してみる

Tellus-APiの方はだいたい理解できたので、次はJupyterNoteでの画像処理に挑戦してみました。
Pythonでの画像処理は初めてなので、[pashango2](https://qiita.com/pashango2)さんのエントリ「[PIL/Pillow チートシート](https://qiita.com/pashango2/items/145d858eff3c505c100a)」を参考にしました。なお、緑地抽出に適した処理手順はわからないので、色々と施行錯誤しながら緑地抽出を試みてみました。

<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">TellusOSを使って、衛星画像の画像処理を練習中です。光学画像から緑地を抽出し、GeoJSONを生成する処理を試してみました。<a href="https://t.co/QZr92WxBAQ">https://t.co/QZr92WxBAQ</a><a href="https://twitter.com/hashtag/Tellus?src=hash&amp;ref_src=twsrc%5Etfw">#Tellus</a> <a href="https://t.co/3xqLna9jTW">pic.twitter.com/3xqLna9jTW</a></p>&mdash; t_mat (@t_mat) <a href="https://twitter.com/t_mat/status/1112747084141068289?ref_src=twsrc%5Etfw">April 1, 2019</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

##画像処理

####1.元画像の読み込み
ASNARO-01「東京北部2016/05/05」から、適度に緑地があるが場所（東京メトロ王子駅付近）をスクリーンショットし、その画像（/tellus_data/に保存されている）をJupyterNoteに読み込みました。

```python
from PIL import Image, ImageFilter ,ImageOps,ImageMath,ImageChops
import numpy as np
import matplotlib.pyplot as plt
im = Image.open('../../tellus_data/20190331-222445.tiff','r')
plt.imshow(np.array(im))
im.save('../data/20190331-222445_0.png')
```
![001.png](https://qiita-image-store.s3.amazonaws.com/0/55735/dc7f776a-56ce-c3d4-2ef3-374faca94578.png)

####2.画素の強調
元画像がくすんでいるので、同じ画像を4回ソフトライトして画素を強調してみました。

```python
def _blend_f(img1, img2, func):
    blend_eval = "convert(func(float(a), float(b)), 'L')"
    bands = [
        ImageMath.eval(
            blend_eval,
            a=a,
            b=b,
            func=func
        )
        for a, b in zip(img1.split(), img2.split())
    ]
    return Image.merge(img1.mode, bands)

def _soft_light(a, b):
    _cl = (a / 255) ** ((255 - b) / 128) * 255
    _ch = (a / 255) ** (128 / b) * 255
    return _cl * (b < 128) + _ch * (b >= 128)
im2=_blend_f(im, im, _soft_light)
plt.imshow(np.array(im2))
im3=_blend_f(im2, im2, _soft_light)
plt.imshow(np.array(im3))
```
![03.png](https://qiita-image-store.s3.amazonaws.com/0/55735/6fc01118-db09-8dbe-0cfc-6fd36f5b3d83.png)

####3.ポスタライズ
ポスタライズで色を単純化しました。

```python
im4=im3.convert("RGB")
im5=ImageOps.posterize(im4, 4)
plt.imshow(np.array(im5))
im5.save('../data/20190331-222445_2.png')
```
![04.png](https://qiita-image-store.s3.amazonaws.com/0/55735/31dab0e7-2168-cffc-2b41-e2ae4d2fff1e.png)

####4.ガンマ補正
ガンマ補正で照度、彩度を調整しました。

```python
def gamma_table(gamma_r, gamma_g, gamma_b, gain_r=1.0, gain_g=1.0, gain_b=1.0):
    r_tbl = [min(255, int((x / 255.) ** (1. / gamma_r) * gain_r * 255.)) for x in range(256)]
    g_tbl = [min(255, int((x / 255.) ** (1. / gamma_g) * gain_g * 255.)) for x in range(256)]
    b_tbl = [min(255, int((x / 255.) ** (1. / gamma_b) * gain_b * 255.)) for x in range(256)]
    return r_tbl + g_tbl + b_tbl

im6=im5.point(lambda x: x * 2.0)
plt.imshow(np.array(im6))
im6.save('../data/20190331-222445_3.png')
```
![05.png](https://qiita-image-store.s3.amazonaws.com/0/55735/79ffae36-9428-e611-c639-ea5e2a58e03f.png)

####5.画素の強調
もう一度ソフトライトを行いました。

```python
im7=_blend_f(im6, im6, _soft_light)
plt.imshow(np.array(im7))
im7.save('../data/20190331-222445_4.png')
```
![07.png](https://qiita-image-store.s3.amazonaws.com/0/55735/e0821476-c681-f175-5baa-ea055c506ea4.png)

####6.色抽出その1
OpenCVをロードし、緑系の色のピクセルを抽出しました。

```python
!pip install opencv-python
import cv2
image_file = '../data/20190331-222445_4.png'
im8 = cv2.imread(image_file)
hsv = cv2.cvtColor(im8, cv2.COLOR_BGR2HSV)
lower = np.array([32, 64, 32])
upper = np.array([180, 255, 120])
img_mask = cv2.inRange(hsv, lower, upper)
im9 = cv2.bitwise_and(im8, im8, mask=img_mask)
plt.imshow(np.array(im9))
cv2.imwrite('../data/20190331-222445_5.png', im9)
```
![08.png](https://qiita-image-store.s3.amazonaws.com/0/55735/b205c92d-12a0-3863-9c03-59bc75deceab.png)

####7.色抽出その2
上記では黄緑系統の色が抽出されなかったので、再度、色抽出を行いました。

```python
lower = np.array([32, 64, 32])
upper = np.array([200, 255, 120])
img_mask = cv2.inRange(im8, lower, upper)
im10 = cv2.bitwise_and(im8, im8, mask=img_mask)
plt.imshow(np.array(im10))
cv2.imwrite('../data/20190331-222445_6.png', im10)
```
![09.png](https://qiita-image-store.s3.amazonaws.com/0/55735/2c673cdc-3e92-dfb4-d8ab-29f847e1fcac.png)

####8.色抽出結果の合成
6、7で抽出した画像を合成しました。

```python
im11 = Image.open('../data/20190331-222445_5.png','r')
im12 = Image.open('../data/20190331-222445_6.png','r')
im13=ImageChops.lighter(im11, im12)
plt.imshow(np.array(im13))
im13.save('../data/20190331-222445_7.png')
```
![10.png](https://qiita-image-store.s3.amazonaws.com/0/55735/71e13955-8c69-7337-782e-fc2178b03db7.png)

####9.ピクセルの再抽出
画像のピクセルを検査し、RGBのうちG値が最も高いピクセルを抽出しました。

```python
size=im13.size
im14 = Image.new('RGB',size)
for x in range(size[0]):
    for y in range(size[1]):
        r,g,b = im13.getpixel((x,y))
        if( g > r and g > b):
            im14.putpixel((x,y),(r,g,b))
        else:
            im14.putpixel((x,y),(0,0,0))
plt.imshow(np.array(im14))
im14.save('../data/20190331-222445_8.png')
```
![11.png](https://qiita-image-store.s3.amazonaws.com/0/55735/b01b99e2-759a-7070-ef0f-5ef918fec8e7.png)

####10.収縮・膨張
収縮・膨張により、ピクセル上のごみを処理し、その後、ガンマ補正を行いました。

```python
im15=im14.filter(ImageFilter.MinFilter()) 
im16=im15.filter(ImageFilter.MaxFilter()) 
plt.imshow(np.array(im16))
im16.save('../data/20190331-222445_9.png')
im17=im16.point(lambda x: x * 2.0)
plt.imshow(np.array(im17))
im17.save('../data/20190331-222445_10.png')
```
![12.png](https://qiita-image-store.s3.amazonaws.com/0/55735/92f9dea9-2a06-27ef-6996-286a97a15dd6.png)

####11.減色処理
色を黒を含め、4色に減色処理しました。

```python
def sub_color(src, K):
    Z = src.reshape((-1,3))
    Z = np.float32(Z)
    criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 10, 1.0)
    ret, label, center = cv2.kmeans(Z, K, None, criteria, 10, cv2.KMEANS_RANDOM_CENTERS)
    center = np.uint8(center)
    res = center[label.flatten()]
    return res.reshape((src.shape))
im19 = cv2.imread('../data/20190331-222445_11.png') 
im20 = sub_color(im19, 4)
plt.imshow(np.array(im20))
cv2.imwrite('../data/20190331-222445_12.png', im20)
```
![13.png](https://qiita-image-store.s3.amazonaws.com/0/55735/a82f493e-9973-b9aa-4966-30701ea1f7bc.png)

####12.色の分離
11で緑3色に減色したので、各色のピクセルを抽出して分離し、2値画像に変換しました。
なお、この処理は使い慣れたJavaで実施しました。

<table>
<tr>
<td><img src="https://qiita-image-store.s3.amazonaws.com/0/55735/8cb68255-e12f-5205-31df-c87e5b464c46.png" width="280px" /></td>
<td><img src="https://qiita-image-store.s3.amazonaws.com/0/55735/e295b3e2-c97a-e6a5-50c3-155bfafa93d2.png" width="280px" /></td>
<td><img src="https://qiita-image-store.s3.amazonaws.com/0/55735/4ef4fb3a-92eb-c764-bd9a-7d60fe1bb5eb.png" width="280px" /></td>
</tr></table>

####13.ベクター化
上記の2値画像を[potrace](https://www.google.com/search?q=potrace&oq=potrace&aqs=chrome..69i57j35i39j0l4.246j0j7&sourceid=chrome&ie=UTF-8)でSVGに変換しました。

<img src="https://qiita-image-store.s3.amazonaws.com/0/55735/efb891bb-a65f-ede3-f69d-6bd288fe1ec2.png" width="380px" />

####14.SVG->GeoJSON生成
[svg2geojson](https://github.com/Phrogz/svg2geojson)により、SVGをGeoJSONに変換しました。
Tellusに生成したGeoJSONを読み込ませた画面を示しますが、そこそこうまく抽出できています。ただ、感覚的ですが、5％ほど、緑地の抽出漏れがある感じです。

<img src="https://qiita-image-store.s3.amazonaws.com/0/55735/60e26eb7-ff2e-5408-4fad-5fdf775aa819.jpeg" width="640px"; />
<img src="https://qiita-image-store.s3.amazonaws.com/0/55735/d3dba431-0bd0-d06e-9a1d-8b14ad8b0297.jpeg" width="640px"; />

##雑感
- 光学画像からでもそこそこの精度で緑地が抽出でき、GeoJSONを生成できることがわかりました。また、林相と画像の色相に関係性があるようなら、林相区分等の自動処理も可能かもと思いました。
- 自動処理で緑地を抽出するとなると、処理手順や抽出の閾値を色々と検討・工夫する必要があると感じました。
- ベクター化やGeoJSON生成には外部ライブラリを使用しましたが、独自に処理を実装するか機能拡張した方がスムーズかもと感じました。