ヒストグラム平坦化は、暗部などで細部が見えにくい場合、そのような細部をはっきり見せる目的で使われるアルゴリズムです。基本的な考え方は、通常山があるピクセル分布のヒストグラムをなるべく平坦に直すことで、見えやすくするということです。
ただ、この平坦化を r,g,b チャンネル別々に実行すると、多くの場合色相がずれます。そのため、その色相のずれをなるべくなくそうといくつかの改良版アルゴリズムが開発されています。しかし、先日お示ししたように、この色相のずれという欠点が、じつは褪色フィルム画像の修復に役立つと指摘しておきました。本稿ではこの色相のずれを生む古典的ヒストグラム平坦化アルゴリズムがどのようなものかを、ImageJ API のソースコードを使って解説します。
以前にも、この原理について解説したサイトを紹介しておきましたが、ここでもそのリンクを紹介しておきます。
ただ、具体的にどのようにアルゴリズムを実装しているのでしょうか？要するに、山のあるヒストグラムをなるべく平坦になるように目指して並べ替えるということを行います。
今、8bit の画像を例えに考えてみましょう。8 bit 画像にはそれぞれのチャンネルごとに、 0 ~ 255 の 256通りの明るさのデータを取ります。ブラックポイント (最暗点) が 0, ホワイトポイント (最明点) が 255 の値を取ります。そして、0 ~ 255 までの値に対してヒストグラムを取るということは、例えば 0 の値を取るピクセルは何個、1 の値を取るピクセルは何個、2 の値を取るピクセルは何個、というようなデータが、255 の明るさまで取得されます。
これをどうやってなるべく平坦化するのでしょうか。例えば、今画像全体に、25600 ピクセルあるとします。すると、1つの明るさ当たりの平均ピクセル数は、100ピクセルとなります。とはいえ、実際に各明るさごとに入っているピクセル数はバラバラです。これを 100 ピクセルに均すということはできません。
そこで、ピクセル数を均さない代わりに、各明るさごとのピクセル数に応じて、明るさの階調の間隔を調整します。つまり明るさの値を変えてしまうのです。そのために、各明るさごとに、明るさが 0 をスタートポイントとしたときに累積ピクセル数を求めます。例えば明るさが 128 の場合、明るさが 0 ~ 128 までのそれぞれのピクセル数を足し合わせたものが累積ピクセル数になります。
例えば、仮に明るさ 128 の累積ピクセル数が、12800 だったとします。次に 129 の累積ピクセル数が、12801 だとします。つまり、129の値を取るピクセルが1つしかないわけです。その場合は、129 の値を、128に変更して、128 のピクセルに合併させてしまいます。逆に、128 の累積ピクセル数が、12800 のときに、129 の累積ピクセル数が、14800 だとします。つまり平均 100 ピクセルのところを、129 の値を取るピクセルが、その 20 倍の 2000ピクセルあるということです。この場合は、129 の値をもっと明るくして、元々 128 のピクセルとの間隔を開けます。その場合は、128 の値ももっと下げる必要があるでしょう。つまり累積ピクセル数の差が少ない隣接した明るさのピクセル同士はなるべく諧調の間隔を小さくし (場合によると同じ値にして合併させ)、累積ピクセル数の差が大きいピクセル同士は、ピクセル値を変更して諧調の間隔をあけることで、なるべく全体的にヒストグラムが均されるよう目指します。 これが上のリンク先の図3 の右下のグラフです。要は明度ヒストグラムにおける或る明度区間に属するピクセルに関して、明度 0 からその区間までの累積度数を基に明度の値を割り当てなおす、ということをやっているわけです。
では、具体的なアルゴリズムを見てみたいと思います。ImageJ の Contrast Enhanser API を取り上げます。
API の説明はこちらです。
元々は、ImageJ のプラグインとして書かれています。このコードは、ヒストグラム平坦化だけではないので、ヒストグラム平坦化の部分のみ抜粋してみてみましょう。
基本は equalize というユーザ定義関数ですが、同名の関数が3つあります。一つは ImagePlus (= ImageJ 上でほぼファイルに相当) を引数として取るもの、他は ImageProcessor (= ImageJ 上でファイル中の RGB データの入った実体に相当) を引数として取るものです。このうち一つは、histogram 変数が引数として指定されていない場合、それを引数として取れるようにするためのダミー的関数です。
そして、ImagePlus を引数として取る関数から、ImageProcessor を引数として取る関数を呼び出しています。この後者がアルゴリズムの核心部分が書かれた関数です。
------------
private void equalize(ImageProcessor ip, int histogram) {
ip.resetRoi();
if (ip instanceof ShortProcessor) { // Short
max = 65535;
range = 65535;
} else { //bytes
max = 255;
range = 255;
}
double sum;
sum = getWeightedValue(histogram, 0);
for (int i=1; i<max; i++)
sum += 2 * getWeightedValue(histogram, i);
sum += getWeightedValue(histogram, max);
double scale = range/sum;
int lut = new int[range+1];
lut[0] = 0;
sum = getWeightedValue(histogram, 0);
for (int i=1; i<max; i++) {
double delta = getWeightedValue(histogram, i);
sum += delta;
lut[i] = (int)Math.round(sum*scale);
sum += delta;
}
lut[max] = max;
applyTable(ip, lut);
}
-------------------
この equalize 関数は、ヒストグラム平坦化実行対象の imageProcessor (以下 ip と略)と、ヒストグラムデータの入った配列を引数として取ります。この関数を呼び出す前に、ip に対し、getHistogram() メソッドを適用し、ヒストグラムデータの入った histogram という配列を引数として与えます。なお、ImageJ の場合、getHistogram() メソッドは 32bit 浮動小数点画像では無効なので、8bit もしくは 16bit 画像しかこのアルゴリズムを適用できません。
まず、もし ROI が設定されている (=範囲指定されている) 場合は、リセットし消去します。
次に、画像が 16bit なのか 8bit なのかを判定し、16bit の場合は、明るさの最大値 (max)、およびヒストグラムの区間数 (range) に 65535 を、そうでなければ 255 を与えます。
次にピクセルの累積数を計算するための sum という変数を初期化し、まずヒストグラムの値が 0 (ブラックポイント) に属するピクセル数を入力します。この時、getWeightedValue(histogram, 0) というユーザ定義関数を使っていますが、これは古典的アルゴリズムを使う場合は、そのまま histogram[0] のピクセル数の値を返します。なお色相の変化が起こりにくいアルゴリズムを使う場合は、上の値の平方根を取りますが、この記事の関心は、褪色の修復を可能にする古典的アルゴリズムですので、とりあえず、
getWeightedValue(histogram, i)
と出てきたら、明度 i に属するピクセル数、つまり histogram[i] を取得するものと解釈してください。
その後、ホワイトポイントの1歩手前まで、ピクセル数を2倍して足し合わせていきます。ホワイトポイントに関しては、2倍せず、そのままピクセル数を足します。この理由については後述しますが、基本的には、画像の総ピクセル数 x 2 の値の近似値を求めることになります。但し、2 倍することはヒストグラム平坦化アルゴリズム上必然的に必要なことではなく、いわばこのプログラムの工夫の部分です。
その後、ヒストグラムの区間数を画像の総ピクセル数 x 2 近似値で割りますが、いわば平均値の約2倍の逆数を取ることになります。この値を scale という変数に入れます。
次に元の画像の明るさを変換するための対照表 (Look up table = Lut) を初期化します。この対照表は配列となっており、その要素数は、明度の階調 (256 もしくは 65356) と同じです。lut[i] には明度 i のピクセルをどの明度に変換するか、この後計算した変換後の明度値を入れていきます。そしてまず lut[0] に 0 の値を入れます。これはブラックポイントのデータは明度変換を行わない (0 のまま) ということを意味します。
再び累積値を格納する変数 sum を初期化し、ホワイトポイントの一歩手前まで、for ループで以下の作業を反復します。
まず変数 delta に、明度 i のピクセル数を格納します。
そして、sum にデルタを足し合わせることでそこまでのピクセル数 (x 2) の累積値を計算します。
ここで、lut[i] の値を計算します。この値は、ここまでの累積値に scale を掛けて、整数化した値です。この scale とはいわば各区間のピクセル数平均値 (実はそれをさらに2倍していますが) の逆数です。この意味はどういうことかと考えると (以下平均値を 2 倍していることは一時的に棚上げして説明します) ... 仮に、ヒストグラムの各区間のピクセル数がすべて同じで平均値に一致すると考えましょう。仮にこの平均ピクセル数が 100 だとすると (その場合 8bit 画像なら総ピクセル数は 25600 ピクセルになる)、明度 100 の区間の場合、そこまでの累積ピクセル数は、10000 になるはずです。そこで、ピクセル数の平均値の逆数、つまり 1 / 100 を掛けて、Lut 上の変更する明度の値を求めると、100 になります。ということは、元の明度値が 100 に対するLut の値はそのまま 100、つまり元の 100 の値は変化なし、ということになります。
ところがこの画像がシャドウ域に寄っているとします。そうだとすると、 100 の値のところの累積ピクセル数は、20000 になっているかもしれません。そこに平均値の逆数を掛けると、200 になります。つまり、100 の明度のピクセルを、200 の明度に変更するということになります。これは、シャドウ域の偏りをなるべく直すために、シャドウ域にあるピクセルを明るくするということを意味します。
逆に、画像がハイライト域に寄っている場合は、100 の明度までの累積値が、10000 を下回る数になっているはずです。例えば、5000になっているかもしれません。それに平均値の逆数を掛けると 50 になりますので、100 の明度を 50 に直す、つまり全般的に暗くして平坦化をなるべく実現しようとする、ということです。
このようにして、0 から当該区間までの累積度数に基づいて、当該区間の明度値を Lut を使って振りなおしていきます。
ただ、このコードでは scale に平均値の逆数ではなく平均値の2倍の逆数を取っています。これは lut の値を計算した後さらに、また delta 値を累積値に加算しているためです。これは値を移動する場合でもなるべく区間の中間値に変更させようという意図 (工夫) があるためだと思います。ですので scale は平均値の2倍の逆数を取るのです。また同様の意図から、scale 値を計算するときの sum を計算するときに、最初のブラックポイントと最後のホワイトポイントに関しては、ピクセル数を 2 倍せずに加算しています。
そして最後にこうやって作成した lut を ip に適用して画像を変換します。
これを、R, G, B 3チャンネルに対して実行しますので、各区間に入るピクセルの頻度は異なるものの、ある特定の明度まで (例えば 0 〜 100) の累積頻度は、多少のずれはあるものの概ね揃うことになります。結果として、R, G, B の各ヒストグラムは、全く同じになるわけではありませんが、概ね似た形になります。これは褪色補正の基本原理に合致しますので、褪色補正効果が現れる、ということになります。
なお、ここでヒストグラムの形をある程度揃えたとしても、しっかり R, G, B の違いが残っている、ということが重要で、著しく褪色が進んだ画像の場合は R, G, B の違いがほとんどなくなりほぼモノクロ画像になる場合があります。この場合は色情報が失われたということを意味しますので、着色以外に修復手段がないということになります。
一方、「改良」版のヒストグラム平坦化は、なるべく R, G, B の差を維持したまま (つまり色を変化させないまま)、ヒストグラムの平坦化を行おうとしますので、全く褪色補正効果が現れません。但し、単純に暗い部分をはっきり見せる、という目的のためには、色相が変わらない (R, G, B 差が維持されたままの) 方が好都合です。Photoshop はおそらくこちらの「改良」版のヒストグラム平坦化を実装しているものと思われますので、Photoshop では褪色補正効果がないのです。また、ImageMagick も最近のバージョンでは「改良」版を実装しています。
なお、褪色画像に対して、画像によっては非常に効果的な修復効果が得られますが、その理由の一つには、Lut を使って画像の値を変換しているということがあると思います。それにより他の方法 (例えばトーンカーブを使うなど) で変換するよりも大胆な画像変換が可能になるからです。
また、このアルゴリズムは褪色がある程度進んだ画像では有効ですが、褪色があまり進んでいない画像に適用すると結構色がおかしくなります。褪色が軽度の場合は別の補正手段を使った方が良いでしょう。
なお、以上から明らかなように、このアルゴリズムは諧調併合やトーンジャンプが必然的に起こります。そこで、オリジナルが 16bit 画像であると、元々の諧調間隔が狭い分、ヒストグラム平坦化の結果もおおむね滑らかになりますので、トーンジャンプについても問題にならないレベルにとどまります。従ってこのアルゴリズムを使って、補正を行うつもりなら、なるべくフィルムは 16bit でスキャンして保存することが望ましいことになります。
-----------
なお、褪色補正には色相のシフトを回避するアルゴリズムの「改良」はむしろ不都合で意味がありませんが、余談として、この「改良」について触れたいと思います。ImageJ のこのプログラムの「改良」版アルゴリズムでは、累積ピクセル数の平方根を取ることで色相のシフトを軽減するようになっています。しかし、他に L*a*b* や HSL, LCh に変換して L 値にのみヒストグラム平坦化を適用することで色相のシフトを回避するアルゴリズムを採用しているケースもあるようです。この辺りはプログラムによってどのように実装するかで結果が微妙に異なります。
-----------
[関連サイト]
