QGIS3.x系向けプラグイン「MagicWand」を作成しました

2020.1.23 Ver1.3について

Ver1.2まではちょっとした修正だったのですが、1.3は処理速度が大幅に改善される修正が含まれており、詳細を記載しておきます。

ポリゴン生成プロセスの修正

従来は、色しきい値内となったピクセル1つに対し矩形が1つ生成されていました。この仕組みでは、画素数の増加につれ指数関数的に矩形の数が増え、処理速度が著しく低下していました。

具体的にどの処理に時間がかかるかといえば、ディゾルブです。矩形の数が多ければ多いほどディゾルブの所要時間が増加していました。ならば処理すべき矩形の数を減らせばよかろう、ということで、以下の画像のような処理へ変更しました。

ご覧の通り、色しきい値内となったピクセルが横軸方向に連続する場合、（イメージ的には）連結してひとつの矩形として扱うようになっています。これにより矩形の数が大幅に削減され、処理速度が大幅に短縮されました。精度を高めてもそれなりの速度で動きます。

はじめに

わたしは仕事柄GISで農地を見る機会が多く、オルソ画像からほ場のポリゴンを作ったりする事もあります。基本的にポリゴンの頂点をマウスのクリックでカチカチして作ります、めんどくさいです。
たとえばグラフィック系のソフトでよくある「自動選択」機能（MagicWandとも言う）みたいに、画像からポリゴンを作れたらなぁと思ってプラグインを作り始めました。QGIS3.x系は日本語情報が少なく、ライブラリに慣れるまで時間がかかりましたが、ざっくり2週間とちょっとで完成しました。本記事ではプラグインの機能や制作過程などを掲載します。

MagicWand

####2019/8/6追記
QGISプラグインリポジトリに追加されました。
Magic Wand/QGIS Python Plugins Repository

MagicWand/GitHub
##機能

• クリックした点の色によるマップキャンバスの分析、ポリゴン生成
• 精度、色しきい値の調節（それぞれ10段階）
• クリックした点の周辺のみ、もしくはキャンバス全体の地物生成か選択可（Single Mode）

開発過程

##開発方針
以下のプロセスで前述のような機能を実装しようと検討。

##問題点
実装自体は着々と進んでいきましたが、いくつか問題が発生しました。

###1.画像の分析方法
####ライブラリについて
QGISプラグインで使えるライブラリはGDAL/OGRやnumpyくらいです。つまり、画像処理の鉄板ライブラリOpenCVが使えません。使えるようにする方法はあるようですが、プラグインの使用者も同様にOpenCVを導入しなければならなくなるため断念。numpy一本で画像処理を行います。

参考：QGISでOpenCVリベンジ - Qiita

####色差について
コンピュータ上の色は（最もポピュラーな形式で）Red, Green, Blue, Alpha（RGBA)でそれぞれ256段階で表現されます。たとえばある色と別の色の差の大きさは、単純にはその2つの色のRGBAのベクトルごとの絶対値の合計で表現出来ます（ユークリッド距離）。たとえば[R,G,B,A]=[63, 127, 191, 255]と[191, 127, 63, 255]の色差は、128+0+128+0=256となります。しかしながら、人間が感じる色の差はそれほど単純には描写出来ないようで、様々な計算方法があります。

参考：色差 - Wikipedia

###2.処理速度
キャンバスを画像化して、1ピクセルごと判定のうえ矩形を生成し、すべての矩形を結合する…という処理が軽いはずがありませんでした。私の環境（画面解像度1920x1080）で生成される画像は概ね1100x500というサイズであり、550000画素となります。当然の結果ですね。

###3.矩形の出来栄え
しきい値内に収まったピクセルで矩形を作成し結合すると、大きな矩形が出来上がる訳ですが、以下のように虫食いが発生してしまいます。しきい値を大きくすれば解決しますがそれだけ精度が下がります。またたまたましきい値内となったピクセルで、小さなノイズ矩形が発生してしまいます。

##解決方法
###1.画像の分析方法
すべてnumpyで実装しました。キャンバスから出力した画像をnumpyのndarrayに変換し、クリックした点をピクセルごとに比較。差がしきい値未満か判定し、しきい値内に収まっているピクセルはTrue,それ以外をFalseとし、二値化させています。以下は画像処理を担っているImageAnalyzerクラスの全コードです。

image_analyzer.py

import numpy as np

class ImageAnalyzer:
    def __init__(self, image):
        self.image = image

    def to_ndarray(self, resize_multiply):
        scaled_img = self.resize(self.image, resize_multiply).convertToFormat(4)

        width = scaled_img.width()
        height = scaled_img.height()

        ptr = scaled_img.constBits()
        ptr.setsize(scaled_img.byteCount())
        arr_rgba = np.array(ptr).reshape(height, width, 4)
        arr_rgb = np.delete(arr_rgba, 3, 2)
        return arr_rgb
        #arr_rgb structure
        #img = x1y1 x2y1 ... xny1
        #      x1y2 x2y2 ... xny2
        #            ...
        #      x1yn x2yn ... xnyn
        #then ndarray is [[x1y1, x2y1 ... xny1],
        #                 [x1y2, x2y2 ... xny2],
        #                 [x1yn, x2yn ... xnyn]]
        #xnyn = [blue, green, red]

    def resize(self, image, resize_multiply):
        scaled_img = image.scaled(image.width() * resize_multiply, image.height() * resize_multiply, True, False)
        return scaled_img

    def to_binary(self, point, resize_multiply=0.2, threshold=50):
        red, green, blue = self.get_rgb(point)
        img_ndarray = self.to_ndarray(resize_multiply)
        abs_ndarray = abs(img_ndarray - [blue, green, red])
        sum_ndarray = abs_ndarray.sum(axis=2)
        max_ndarray = abs_ndarray.max(axis=2)
        true_index = sum_ndarray + max_ndarray * 0.5 < threshold
        return true_index

    def get_rgb(self, point):
        pixelColor = self.image.pixelColor(point.x(), point.y())
        red_value = pixelColor.red()
        green_value = pixelColor.green()
        blue_value = pixelColor.blue()
        return (red_value, green_value, blue_value)

• to_ndarray()で画像（QImage）をndarrayに変換しています。QImageはRGBAそれぞれの値を持ちますが、Alphaはすべてのピクセルで同値であるとし、値を捨てています（計算量低減のため）。
• to_binary()はndarrayを入力すると、すべての要素を基準色（ここではself.get_rgb(point)）と比較し、しきい値内となる要素のインデックスを返しています。

####色差について
以下のように判定しています。

abs_ndarray = abs(img_ndarray - [blue, green, red])
sum_ndarray = abs_ndarray.sum(axis=2)
max_ndarray = abs_ndarray.max(axis=2)
true_index = sum_ndarray + max_ndarray * 0.5 < threshold

• absは色差の要素ごと絶対値です
• sumは要素ごと絶対値の総和です。
• maxは要素ごと絶対値の最大値です。

true_index = sum_ndarray + max_ndarray * 0.5 < threshold

これが判定式になっています。sumはともかくなぜmaxが入っているのか。たとえば青と水色のsum, maxは以下のように表せます。

青=[0,0,255], 水色=[127,127,255] -> sum=254, max=127

ここで、sumが同じ値になる別の色と比較してみると

青=[0,0,255], ほぼ黒=[0,0,1] -> sum=254, max=254
青=[0,0,255], ほぼ黃=[0,254,255] -> sum=254, max=254

果たして、水色と黄色はどちらが青色に近いでしょうか。人によって様々な意見があると思われますが、本プラグインで主に見る事になるのは航空写真になると思います。もちろんノイズ等があり均一な画像ではないはずです。そのうえで、同系統の色を検出する必要があります。つまり、青を基準色とするなら、黄色よりも水色を検出したい訳です。という事で、判定式にはmaxを入れてあります。実際にはもっと厳密に計算すべきとは思いますが、それなりに思い通り検出してくれるので妥協しています。

---

2.処理速度

numpyの最適化

numpyはforループを使うと遅くなる、らしいです。というのも、numpy自体が持つメソッドを活用する事で、C言語などで演算を行えるから、との事です。

true_points = np.where(self.bin_index)
func = lambda x, y, size: self.rect_geo(x, y, size)
np_func = np.frompyfunc(func,3,1)
size_multiply = self.map_canvas.width() / self.bin_index.shape[1]
geos = np_func(true_points[1], true_points[0], size_multiply)

矩形を生成する処理の一部抜粋です。二値化された配列からジオメトリを生成する処理です。ここは当初forループを実装していたのですが、パフォーマンス向上のためnumpyのメソッドを用いてforループを使わないよう実装しました。

参考：Pythonのリストの全要素に任意の関数をapplyする最速の方法 - Qiita

画像の圧縮

numpyの最適化を実施しても劇的な改善はありませんでした。そりゃ、最大で50万もの矩形を生成していれば計算量自体が膨大ですから当然でした。キャンバスを画像として出力して、それをピクセル単位で分析、ひとつのピクセルと同じ位置・サイズで矩形を作成、そして結合というプロセスでしたね。つまり画像サイズと計算量は指数の関係にあります。したがって画像を縮小すれば改善されると考えました。

image_analyzer.py

    def resize(self, image, resize_multiply):
        scaled_img = image.scaled(image.width() * resize_multiply, image.height() * resize_multiply, True, False)
        return scaled_img

ImageAnalyzerクラスは、init()で画像を要求します。その画像をresize()で指定された倍率で縮小しています。

さて画像を圧縮しましたが、矩形を生成するPolygonMakerクラスはその事をまだ知りません。すると以下の画像のようになります。

右下の緑地のポリゴンを作ろうと思ったのに、位置がずれてサイズが小さくなっています。キャンバス画像が勝手に縮小されているのだから当然ですね。以下はPolygonMakerクラスの一部抜粋です。

polygon_maker.py

    #make rectangle geometry by pointXY on Pixels
    def rect_geo(self, x, y, size_multiply):
        point1 = self.map_canvas.getCoordinateTransform().toMapPoint(x * size_multiply, y * size_multiply)
        point2 = self.map_canvas.getCoordinateTransform().toMapPoint((x + 1) * size_multiply, (y + 1) * size_multiply)
        geo = QgsGeometry.fromRect(QgsRectangle(point1.x(), point1.y(), point2.x(), point2.y()))
        return geo

何をやっているかというと、画像の縮小の逆演算です。画像が縮小された分だけ矩形を拡大しています。ちなみに縮小時はImageAnalyzerにおいて、その縮小率を明示的に与えますが、PolygonMakerでは縮小率を、画像とマップキャンバスのサイズ比率から自動計算しています。

結果として、この画像圧縮によってやっと実用的な処理速度となりました（平均50秒くらいかかっていたのが数秒にまで短縮）。以上をまとめると以下の画像のとおりとなります。

3.矩形の出来栄え

• 矩形内部のノイズは、矩形をバッファで拡大して埋めています。Ver1.2からは、プロセッシングツール「穴の削除」により、一定サイズ以下のノイズを削除しています。
• 矩形外部のノイズ処理は、一定値よりも小さい矩形を削除する事で実装。

polygon_maker.py

for feature in single_features:
    if single_mode and not feature.geometry().contains(self.map_canvas.getCoordinateTransform().toMapPoint(point.x(), point.y())):
        continue
    #一定面積よりも小さい矩形は追加しない
    if feature.geometry().area() < minimum_area * noise_multiply:
        continue
    #一定値でバッファ、同値で縮小する事で矩形内部のノイズを削除
    output_geo = feature.geometry().buffer(1 * buffer_dist, 1).buffer(-1 * buffer_dist, 2).simplify(torelance)
    output_feature = QgsFeature()
    output_feature.setGeometry(output_geo)
    output_features.append(output_feature)

そして完成

最終的に処理の流れは以下のとおりとなりました。

色差の判定方法や処理速度等、まだまだ改善の余地はありますが、コード自体もそれなりによくまとまっていて個人的には満足しています。以下のQiita記事で、QGISプラグイン全般についてまとめてありますので、併せてご覧ください。

QGIS3.x系向けプラグイン作成手順や各種処理の実装方法について - Qiita