空間フィルターを使った天体写真の画質改善

ようかん@http://blogs.yahoo.co.jp/youkan2000/ 氏から素材画像など必要なデータ一式を提供していただいたので、以前の記事に書いた以下の事を試してみた。

RegiStax で処理する前に空間周波数 65/mm 以上の成分を取り除くような処理を行うともっと少ないスタック数で同じ結果が得られるかもしれない。（中略）プログラムは比較的簡単に作れるから、そのうち試してみたいと思うが、問題は素材の準備だ・・・。

途中いくつか新しい発見があり、何回か作業をやり直しているので時間がかかってしまったが、当初の目論見に近い結果が得られた。

/*-----------  -----------*/

背景と仮説

以前の記事で「ようかんのあやしい実験ブログ」の記事「月面写真のコンポジット処理の実験 (FUJIFILM X-M1 編)」に掲載されていた画像を使い、スタックする画像数の違いによるコントラストの空間周波数分布の変化を調べた。

画像数が 64 を超えて増加すると空間周波数 65/mm 以上のコントラストが次第に低下するが、目で見た画質の変化はほとんど感じられない。この周波数領域に意味のある画像信号が存在しないためだと考えられた。これは次のことからも裏付けられる。

F ナンバー毎の光学系の分解能を MTF カーブで示したもので、当該月面画像の撮影に使われた F 15.4 では空間周波数 65/mm 以上の信号成分の減衰が大きい。そのため画像センサーに結像した画像にその空間周波数領域の成分がほとんど含まれていなかった可能性がある。意味のある画像信号が存在しない空間周波数領域の信号は捨ててしまっても画質に影響は無い。

半分の周波数範囲を捨てても意味のある画像信号に影響が無いなら、ノイズ電力が半減するので S/N 比が 3db 向上する。スタックで同じ効果を得るには画像数を２倍にする必要があるからこの違いは大きい。

空間フィルターの設計

カットオフ特性が鋭いほど効果が高いはずなので、そのようなフィルターの設計が比較的容易な FIR 方式を使った。FIR は周囲の画素の移動平均を計算する方式のフィルターである。

Fig. 4 はオーディオ信号のような１次元の FIR に関するものだ。空間フィルターでは２次元になるが、本質的にやることは同じである。ただし２次元では向きを考える必要がある。例えば「目標の周波数-振幅特性を作る」場合、通過帯域を正方形にすると縦横方向と斜め方向で特性が違ってくる。円形のパターンにすれば等方性の特性になる。

鋭いカットオフ特性を得ようとすると移動平均の重みに負の値が出てくる。計算結果も負の値を含む場合が出てくるが、画像の場合「負の明るさ」は表現できない。レタッチをしている感覚だと「負の明るさ」は値ゼロの黒レベルにしてしまえば良さそうに思えるが、それだと信号波形をクリップすることになり具合が悪い。波形が歪むことで元々無かった周波数の信号成分が発生してしまう。これはノイズで、フィルターの特性を大きく損なう原因になる。

今回は計算結果の最大-最小値が画像データの表現できる範囲に収まるよう、全体の黒レベルとコントラストを調整するようにした。このようにして作ったフィルターにテスト信号を入れて特性を確認した。比較のため、ノイズ除去の目的で使われる場合も有る「ガウスぼかし」の特性も一緒にプロットしている。

FIR フィルターは浮動小数点で計算し、結果を画像データの形式に変換した。「float」と表示したプロットが浮動小数点の内部データで、狙った特性が出ている。これを「tiff16」形式に変換すると高い空間周波数のレベルが若干上昇する。変換する形式が「tiff8」では高い空間周波数のレベル上昇が大きい。どうやらこれは切り捨て誤差ないしは量子化雑音、と言うことらしい。

以前の記事を書いているときはスタック前の画像にフィルターを適用するよう考えていた。しかしこれは処理量が多いので、スタック済の画像に適用するよう変更した。処理対象画像数が数十～数百倍違うのでこの差は大きい。

FIR フィルターを構築するのに当初使っていた画像処理ライブラリーが tiff16 形式をサポートしていないので別のツールで tiff8 形式に変換していたのだが、Fig. 5 の結果を見て考え直した。Wavelet のような強調処理はノイズも誇張してしまう場合が多いので、そう言った処理を行う対象画像のノイズはなるべく減らしておいた方がいい。しかし tiff16 形式の読み書きともにちゃんとサポートしている道具は意外と少ないので、実際どの程度影響するのか調べてみた。（tiff8 でもあまり影響が無ければそれで誤魔化したい）

ようかん氏に提供していただいた中にスタック済 tiff16 画像データがあったので、これをそのまま Wavelet 処理した場合と一旦 tiff8 形式に変換してから Wavelet を行った場合を比較したものだ。 tiff8 形式に変換すると高い空間周波数領域のレベルが明らかに上昇している。

仕上がり画像を別々に見ても分からないが、並べて比べると違いが分かる。ただし違うと分かるだけで、どちらが良いかまでは判別できない程度の違いである。しかし数値にすると明らかに差がある。一旦 tiff8 に変換してしまうと高い空間周波数領域の特性の違いがノイズに埋もれて判別できなくなる。画質に対しても悪影響が無いとは思えないので、tiff16 をちゃんとサポートしているモジュールを探して使うようにした。

結果の検証

実際のスタック済画像にフィルターを適用して Wavelet 処理を行った画像を使って検証を行う。フィルターを適用せずに Wavelet 処理を行った画像と比較して変化を確認するのだ。どちらも同じ Wavelet 設定: 「フィルターを適用していない状態でノイズが目立ち始める限界の設定」を使っている。

フィルターは、 X-M1 の画像センサー解像度の半分よりも高い空間周波数の成分を除去するようにカットオフを設定した。数値では狙った通りの効果が得られたが、プロットを見ると意味のある信号を少し削ってしまっているようだ。

問題は画像の仕上がり具合だ。同じ条件で比較するため、フィルターを適用した画像についてノイズが目立ち始めるまで Wavelet を強くしていったのだが、その際のノイズの出方がずいぶん違う。

上は過剰に Wavelet を適用したものだ。ふつう砂嵐的なノイズの出方になるが、フィルターを適用してあると偽像 (artifact) のような感じになる。砂粒の大きさが２倍になっただけなのだが、１つの砂粒を構成するピクセルが増えるため丸みを帯びてくる。そのため「ノイズが目立ち始める」ポイントが把握し辛い。言い換えればノイズが目立ちにくくなるようだ。

フィルターを適用した画像は、適用していない同じスタック画像数のものよりも１段画像数の多いものに近い仕上がりになっているように見える。スタック画像数２倍の効果が得られているようだが、同時に画質が少し劣化してしまったようだ。エッジが鈍くなったり小さなクレーターが大きくなったりしてしまっている。やはりカットオフの設定が低過ぎたようだ。

カットオフをセンサー解像度の 70% にして比較をやり直してみた。

これなら意味のある信号の損失はほとんど無いはずだが、同時にノイズ除去効果も減少しているはずだ。

フィルターを適用したことによる画質劣化は無いと言ってよいだろう。またノイズ除去効果はこちらの方が優れている様に見える。画像がシャープになったのでその様な印象を受けるのだろう。スタック画像数 36 枚のものにフィルターを適用すると、私の眼にはフィルター無しの 64 枚より優れているように見えるし、フィルターあり 16 枚のでもフィルター無し 64 枚に肉薄している（これは少し苦しいか？）。

重要なのは目で見た印象なので、これで良いと言うことだ。

ちなみにフィルターを適用した n=4～36 の画像は、 Wavelet のレイヤー１以外は使う必要が無かった。また n=64 の画像ではハイライト部分が白トビ気味になるのを抑えるためレイヤー２以上に -5 を設定している。

課題

今回は機能の実現を優先したため、パフォーマンスは二の次になってしまった。そのため処理能力に課題を残している。

今回 32bit 版 Python を使ったが、トリミングしていない X-M1 の画像 (16Mピクセル) はメモリー不足で扱えない。 3,200x2,700、つまり約 8M ピクセルが上限で、これは何とか今回の三日月が収まる大きさだ。これまでの比較は VGA サイズ (640x480) の画像を切り出して行った。処理のし易さを優先したデータの持ち方にしているのが原因だと考えられ、これを直せば同じ処理環境でも 24M ピクセル位まで対応できるのではないかと予想している。

もう一つは処理速度だ。VGA サイズの画像の処理に 12秒かかっている。16M ピクセルの画像だと 10分位かかるはずだ。移動平均の係数行列の対称性を利用すると簡単に約 1/2 にできる。またアルゴリズムで頑張れば、マルチコアのような力技に頼らなくても更に 1/2～1/4、最終的に１分台に持ち込むことは可能ではないかと考えられる。なお計算量の多い部分は C のモジュールで行っているので Python を使ったことによる影響は２倍程度、フィルター処理を C または C++ で書くと１分を切れるかもしれない。

代替手段について

フィルターを使わなくても同じような効果の得られる機能が標準的なワークフローで使っているツールの中に隠れているかもしれない。もしそう言った機能で代替可能なら敢えてフィルターを使う必要は無い。あるいかそう言った機能があっても今回のフィルターを使うべきなら、そのことを明らかにしておく必要があるだろう。

以前の記事を書いたときは代替手段を探しきれなかったが、今回改めて探すと２つ候補が見つかった。

１つは現像ソフト「RAW FILE CONVERTER EX powered by SILKYPIX」の「ノイズリダクション」設定である。

これまで使った画像は「ノイズ除去」の設定がゼロで現像されているが、これを強くするとどうなるかと言うことである。

RegiStax を使う前の、現像しただけの画像による比較だ。どうやら「ノイズ除去」の実体はぼかしらしい。ノイズが減る効果よりも画質が甘くなる弊害の方が大きい。もっとノイズの多い画像だと効き方が違うのかもしれないが、少なくとも今回の月面画像に使うべき機能ではない。

もう一つ見つかったのは RegiStax の Wavelet で「Wavelet filter」で「Gaussian」を選んだ時に使える「Denoise」設定である。

Wavelet 処理と同時にノイズ除去を行うと、かなり高い効果が得られる場合があるかもしれない。しかし少し試してみた感触では、これもガウスぼかしに近い挙動のような印象を受けた。従って FIR フィルターを併用する余地がありそうだ。

今回の Wavelet filter は Default だ。ようかん氏から設定ファイルを提供していただいたので、それをそのまま使っている。 Gaussian に切り替えると Wavelet の効き方が全然変わってしまうので、一から設定をやり直すことが必要になる。こうなると定量的な比較うには相当な作業量が要るので、これ以上踏み込まないことにした。

FIR フィルターの有効性と限界

今回のように FIR フィルターが効果的なのは、望遠鏡の分解能がカメラのセンサー解像度を下回っている場合の様だ。今回の例では Fig. 2 に示したように望遠鏡の分解能はセンサー解像度の大体半分くらいになる。

仮に２倍テレコンバーターを使わずに F8 だったと仮定すると、この場合の分解能はセンサー解像度とほぼ同じになるので意味のある信号の損失無しに LPF を適用することは不可能だ。カメラのセンサー上に結像する画像の大きさが半分になることを考えれば、何が起こるか容易に想像できるだろう。

別の言い方をすると、今回の月面画像はオーバーサンプリング状態になっていた。オーバーサンプリングされた信号に LPF を適用して S/N 比を稼ぐのはデジタル信号処理でよく使われるテクニックである。この時多くの場合シャープなカットオフ特性を持つ FIR フィルターが使われる。

大気の揺らぎが問題になる場合、望遠鏡の分解能ではなく揺らぎで大きくなった星像の大きさを基準に考えるべきかもしれない。ただしこれは主に大口径の望遠鏡の場合の問題だ。私には縁のない世界の話である。

今回のまとめ

F16 程度の光学系と 16 メガピクセル APS-C センサーの組み合わせで月面写真を撮る場合、FIR 型空間ローパスフィルターを利用するとスタック（コンポジット）画像数を２倍にしたのと同じ効果が得られることが分かった。光学系とセンサーの組み合わせについて有効な条件の範囲を予測できるが、被写体が異なる場合にどのような差異があるか不明である。これを明らかにするには更に実写で試す必要があるだろう。

今回の記事は良質な元画像が多数あったことにより可能になった。文末ながら、画像など多数のデータを提供していただいた「ようかん」氏に謝意を表したい。

ありがとうございます。