最近某細胞の論文の問題で会見が開かれているようで、テレビではその様子が頻繁に放送されています。会見中は大量にフラッシュが焚かれるため「フラッシュの点滅にご注意下さい」という表示が行われますが、これを見て「時間軸でメディアンフィルタを適用すれば軽減できるのでは」と考え検索した自分は
Median Filtering in Constant Timeという論文に遭遇しました。ご覧の通りこれは時間軸のメディアンフィルタとは何の関係も無いのですが、これを少し見た後二分探索を応用することにより本記事のテーマである"順次追加・削除される集合のメジアンを高速に計算する"アルゴリズムの存在に気付きました。
(このアルゴリズムに何か一般的に認知された名称があればご連絡下さい。)
概要
このアルゴリズムは、[0, N - 1]の範囲の整数M個の集合Sからヒストグラムを生成し、集合Sの要素を順に並べた場合の任意の位置の値をO(log₂ N)で求めることができます。さらに、この集合への値の追加/削除もO(log₂ N)で行うことができます。
このアルゴリズムの用途としては、画像のメディアンフィルタが考えられます。メディアンフィルタは画像の各ピクセルを元画像の周辺のピクセルの中央値とするフィルターであり、ノイズ除去フィルタとしては代表的な部類に入ります。メディアンフィルタをヒストグラムの生成によって実装する場合、メディアンを計算する対象とするピクセルのマスクを順次移動して行い、新しくマスクに入ったピクセルをヒストグラムに追加し、逆にマスクから出たピクセルをヒストグラムから取り除くようにして高速化するのが普通ですが、単純なアルゴリズムではこのヒストグラムからメディアンを得るのにO(N)の処理が必要となります。
内容
単純なアルゴリズムでは、求めたい位置の値を検索するために、累計度数を線形探索します。このためO(N)の計算量が必要となりますが、本アルゴリズムでは二分探索を使用します。しかし、累計度数を配列に格納する場合、この配列への値の追加・削除に必要な計算量は平均・最悪O(N)となってしまいます。
そこで、本アルゴリズムでは対象となる整数のヒストグラムだけではなく、これを小さくした整数のヒストグラムを格納します。具体的には元の値$X$のヒストグラムに加え、$\mathrm{floor} \left( X / 2^L \right)$ のヒストグラムを格納します。最初に最も小さいヒストグラム上での位置を求め、その情報を元にそれより1段階大きいヒストグラム上での位置を高速に求めます。これを元の値のヒストグラムまで繰り返せば、求めるべき値を正確に求めることができます。
具体的な動作
例えばN = 8についてこのアルゴリズムを使用することを考えます。例として、以下の様なヒストグラムとなるような集合を考えます。
| 度数分布0 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
|---|---|---|---|---|---|---|---|---|
| 個数 | 35 | 20 | 18 | 40 | 30 | 20 | 10 | 35 |
このアルゴリズムでは、これに加え複数のヒストグラムが生成されます。
| 度数分布1 | 0-1 | 2-3 | 4-5 | 6-7 |
|---|---|---|---|---|
| 個数 | 55 | 58 | 50 | 45 |
| 度数分布2 | 0-3 | 4-7 |
|---|---|---|
| 個数 | 113 | 95 |
| 度数分布3 | 0-7 |
|---|---|
| 個数 | 208 |
ここで下から120番目の値を求めることを考えます。
- 度数分布2を左から順に参照します。0-3の範囲が113個であるため、求める値の範囲が0-3ではなくこの右隣の4-7にあると分かります。
- 度数分布1を4-5から順に参照します。4-5の個数が50であるため、0-5の個数は50 + 113 = 163であり、求める値の範囲は0 - 5と4 - 7の共通部分である4 - 5にあると分かります。
- 度数分布0を4から順に参照します。4の個数が30であるため、0-4の個数は113 + 30 = 143であると分かります。従って求める値の範囲は4 - 5と0 - 4の共通部分である4であると確定します。
この各手順は全て定数時間で行われます。手順の回数はlog₂ Nとなるため、結果的に計算量はO(log₂ N)となります。
ベンチマーク
Emboxenにてstd::sortと速度比較を行うプログラムをご覧頂けます。(若干コードが汚いですが…)