PostGISラスタで画像処理、続けてPNGファイルに出力するテンプレート

前回までの続きです。説明は後にして、先にクエリのテンプレートを示します。↓

copy (
    select ST_AsPNG(ST_MapAlgebra
    (
        rast, 1, 'st_range4ma(float[][][], int[][], text[])' :: regprocedure,
        '8BUI', 'FIRST', NULL, 4, 4)
    )
    from adbadge_tall
)
to program '"D:/AppsPortable/GitPortable/2.6.2/App/Git/usr/bin/xxd" -r -p > "出力先PNGファイル"';

昨日も書きましたが、今回の中心部分はPostgres Online Journal掲載の ↓ をなぞってます。そこにあった「PostGISのビルトイン関数を使う例」を、テーブル名だけ変えてcopy (...) to programの括弧内に入れたのが、上のテンプレート。COPYコマンド全体は、昨日の最後に行ったPNGファイルへの出力と同じです。

Writing PostGIS raster Map Algebra Callback Functions in PLV8（Postgres OnLine Journal, 2014.1.27）

上のテンプレート中、ST_MapAlgebraというPostGISの関数でいわゆる「ラスタ演算」の外型を決め、それに渡すコールバック関数でラスタ内のピクセル一つ一つに対し演算を行います。その結果を直ちにPNGファイルに出力することで、クエリ以外の作業がなくて楽。

↓ 画像一つ目が元ラスタ、二つ目が上記クエリ（コールバック関数にST_Range4maを使用）の結果。ST_Range4maは「最大値－最小値」つまりレンジを返す関数で、本来は画像処理よりもデータ分析むけです。処理結果を見てのとおり実用的な画像処理にはなりませんが、とりあえず今回はPostgres Online Journalに合わせてこれで。

ST_MapAlgebraおよびコールバック関数（ST_Range4ma以外にも色々あり）の詳細は下記を参照。ST_MapAlgebraにはPostGIS 2.2でマスク機能が追加されましたが、日本語マニュアルの対象は2.2.0devの開発版だったため説明に入っていません。その他にも細かい違いがあるかも。

PostGISの関数のソースを見て、PL/v8への置き換えを考える

コールバック関数ST_Range4maには古いバージョン（PostGIS 2.0から）と新しい版（2.1で追加）があり、今日使ったのは後者。PL/pgSQLで書かれているので ↓ のようにクエリでソースを確認できます。

\pset format unaligned
\pset tuples_only on -- 関数の定義文だけ見るため、設定変更

select prosrc from pg_proc
where proname = 'st_range4ma'
    and proargnames[1] = 'value'; -- PostGIS 2.1で追加された方

DECLARE
        _value double precision[][][];
        min double precision;
        max double precision;
        x int;
        y int;
        z int;
        ndims int;
BEGIN
        min := 'Infinity'::double precision;
        max := '-Infinity'::double precision;

        ndims := array_ndims(value);
        -- add a third dimension if 2-dimension
        IF ndims = 2 THEN
                _value := _st_convertarray4ma(value);
        ELSEIF ndims != 3 THEN
                RAISE EXCEPTION 'First parameter of function must be a 3-dimension array';
        ELSE
                _value := value;
        END IF;

        -- raster
        FOR z IN array_lower(_value, 1)..array_upper(_value, 1) LOOP
                -- row
                FOR y IN array_lower(_value, 2)..array_upper(_value, 2) LOOP
                        -- column
                        FOR x IN array_lower(_value, 3)..array_upper(_value, 3) LOOP
                                IF _value[z][y][x] IS NULL THEN
                                        IF array_length(userargs, 1) > 0 THEN
                                                _value[z][y][x] = userargs[array_lower(userargs, 1)]::double precision;
                                        ELSE
                                                CONTINUE;
                                        END IF;
                                END IF;

                                IF _value[z][y][x] < min THEN
                                        min := _value[z][y][x];
                                END IF;
                                IF _value[z][y][x] > max THEN
                                        max := _value[z][y][x];
                                END IF;
                        END LOOP;
                END LOOP;
        END LOOP;

        IF max = '-Infinity'::double precision OR min = 'Infinity'::double precision THEN
                RETURN NULL;
        END IF;

        RETURN max - min;
END;

3次元配列をぐるぐるループしていますが、単純化すれば、最後の行「max - min」のとおり対象ピクセルおよび周辺の最大値から最小値を引いて返します。「周辺」の範囲を決めるのがST_MapAlgebraに渡した最後の引数二つで、先のクエリでの「4, 4」は縦・横4ピクセル分の隣接範囲。この値を大きくするほど元ラスタの線的な部分が太り（冒頭の画像のとおり）、また計算量も多くなります。

で、このPL/pgSQLストアドをPL/v8で置き換えたら断然速いよ！というのが、今回ならったPostgres Online Journalの記事の趣旨。確かに上のソースを見ると、多次元配列の最大値・最小値を返すだけならJavaScriptの方が簡単に書けるし、v8スクリプトエンジンの処理速度にも期待できます。

また一般的にラスタ演算のコールバック関数は配列への処理が中心なので、今回うまく行けば「PostGISのラスタ演算はPL/v8で」という方針が立つかも。実際やってみた結果は明日以降。引っ張ってるわけではなく^^;テストが何通りかあって記事が長くなるためです。