コールバック関数ST_Range4maのPL/v8版を作る

やっと本題に来ました。昨日の最後にソースを見たST_Range4ma（PostGISに元からある、PL/pgSQLのストアド）と同じ配列演算をPL/v8で定義すると ↓ こんな感じ。元の関数では細かい条件分岐が入っていましたが（第1引数の次元が想定外の場合など）、今回は省略して配列への処理だけで比べる方針。なので後ほど、同様に簡略化したストアドをPL/pgSQLでも作ります。

中身は、配列の最大値－最小値を出すだけで1行。引数のうち第2・3は今回使わないので、テスト時に省略できるようデフォルト値のNULLを入れました。

create or replace function plv8_range4ma
(
    val float[][][],
    pos int[][] default NULL :: int[][],
    userargs text[] default NULL :: text[]
)
returns float language plv8 immutable as
$$
    return Math.max.apply(null, val) - Math.min.apply(null, val);
$$;

-- Check
select plv8_range4ma(
    '{ { { 1, 2, 3 }, { 4, 5, 6 }, { 4, 5, 6 }, { 7, 8, 9 } } }');
+---------------+
| plv8_range4ma |
+---------------+
|             8 |
+---------------+
(1 row)

関数の動作チェックのクエリは、ラスタと同様の3次元配列を渡してます。次元は「ラスタバンド → 南北 → 東西」で、上の場合は「単バンド、南北4ピクセル、東西3ピクセル」という極小のラスタと同じ。この中で最大値=9、最小値=1から戻り値は8になります。

ところで、昨年9月23日・24日にPL/v8ならではの高速な配列Typed Arrayを試しましたが、今回のような多次元配列には未対応なので使えません。あえて第一引数をTyped Arrayにして多次元配列を渡すと ↓ エラー。というわけで今回はJavaScriptの平凡な配列演算 vs PL/pgSQLとの比較になります。

create or replace function plv8_range4ma_tarray
(
    val plv8_float8array, -- Typed Array
    pos int[][] default NULL :: int[][],
    userargs text[] default NULL :: text[]
)
returns float language plv8 immutable as
$$
    return Math.max.apply(null, val) - Math.min.apply(null, val);
$$;

-- in vain
select plv8_range4ma_tarray(
    '{ { { 1, 2, 3 }, { 4, 5, 6 }, { 4, 5, 6 }, { 7, 8, 9 } } }');
ERROR:  NULL element, or multi-dimension array not allowed in external array type

PL/pgSQLの単純化バージョンも作る

前項の最初に書いたとおり、PL/v8と同じ配列演算だけに単純化したST_Range4maを作ると ↓ こんな感じ。元のは配列の次元やピクセル値を細かくチェックし条件分岐してました。それを省き、あと配列のサイズをループ内で参照していたのを最初の一回だけに変更。これで処理時間がどのくらい変わるのかも（本題ではないですが）後で見ます。

create or replace function simple_range4ma
(
    val float[][][],
    int[][] default NULL :: int[][],
    text[] default NULL :: text[] -- PL/v8と違い、引数名を省略できる
)
returns float language plpgsql immutable as
$$
declare
    nras int := array_length(val, 1);
    nrow int := array_length(val, 2);
    ncol int := array_length(val, 3);
    ras int;
    row int;
    col int;
    min float;
    max float;
    tmp float;
begin
    for ras in 1 .. nras loop
        for row in 1 .. nrow loop
            for col in 1 .. ncol loop
                tmp := val[ras][row][col];
                if min is null or tmp < min then min := tmp; end if;
                if max is null or max < tmp then max := tmp; end if;
            end loop;
        end loop;
    end loop;
    return max - min;
end;    
$$;

SQLでもコールバック関数を作る

今回ならったPostgres Online Journalの記事に載っていたので、同様に作って ↓ 比較に使います。配列を行に展開するunnest関数に多次元配列を渡すと、一気に全次元をキャンセルして一列になる。それを普通の集約関数のmaxとminに渡して差を取るだけ。unnest関数の豪快?な仕様がこんな形で使えるとは。

create or replace function sql_range4ma
(
    val float[][][],
    int[][] default NULL :: int[][],
    text[] default NULL :: text[]
)
returns float language sql immutable as
$$
    select max(v) - min(v) from unnest(val) as v;
$$;

unnestだけの動作は ↓ この通り。配列の定義文をtext[]にキャストして渡してますが、これは「一次元配列」への変換でなく配列一般としての意味しか持ちません。次元数はあくまで定義文で決まり、ここでは3次元配列。これがunnest関数で一気にスカラの列になったという様子です。

select unnest(
    '{ { { 1, 2, 3 }, { 4, 5, 6 }, { 4, 5, 6 }, { 7, 8, 9 } } }' :: text[]);
+--------+
| unnest |
+--------+
| 1      |
| 2      |
| 3      |
| 4      |
| 5      |
| 6      |
| 4      |
| 5      |
| 6      |
| 7      |
| 8      |
| 9      |
+--------+
(12 rows)

テストの全体

以上でテストの準備が整いました。コールバック関数は次の4種類。それぞれを用いた場合の出力ファイルも念のため別々にして確認するので、ファイルに付ける接頭辞も下記のとおり決めました。

PL/pgSQL（PostGISに元からある版）→ 出力ファイル st_range4ma_xx.png
PL/pgSQL（単純化した自作版）→ simple_range4ma_xx.png
PL/v8 → plv8_range4ma_xx.png
SQL → sql_range4ma_xx.png

ラスタ演算に用いる「周辺の範囲」は次の5通り。数値は「南北,東西」それぞれのピクセル数で、これが大きいほど各ピクセルの値を出すのに多くのピクセルを使う、従って全体の計算量が増えます。それと処理時間の関連を見るため5通りにしました。ただし最初の「0, 0」は真っ黒な画像になり（ピクセルの値が全てゼロ）画像処理としては無意味。

0, 0（全ピクセルがゼロになる。対象ピクセル＝周辺ピクセルで、どこでも最大値＝最小値だから）
1, 1
2, 2
3, 3
4, 4

以上、4つの関数×5つの範囲の計20通りごとに ↓ のクエリに入れ、まず各1回実行して出力PNGを確認。クエリは昨日作ったテンプレートと同じで、PostGISラスタの画像処理＆COPYコマンドによるバイナリファイル出力を一度に行うもの。この結果、周辺範囲が同じならどの関数でも同じ画像で、演算自体は正常でした。

copy (
    select ST_AsPNG(ST_MapAlgebra
    (
        rast, 1, 'コールバック関数名(float[][][], int[][], text[])',
        '8BUI', 'FIRST', NULL, 周辺範囲・南北, 周辺範囲・東西)
    )
    from adbadge_tall
)
to program '"D:/AppsPortable/GitPortable/2.6.2/App/Git/usr/bin/xxd" -r -p > "出力先PNGファイルパス"';

参考まで、出力した全PNGファイルをZIPで ↓ 置いておきます。真っ黒い画像は、開くアプリがPNGの透過色を認識する場合は「全ピクセルが透過」扱いになります。

output_range4ma_xx.png.zip（1.13MB, 20files）

20通リを各1回では処理時間を議論しにくいので、複数回行って比較した結果を明日以降書きます。