手ぶれ画像 (640x480)＝
推定画像 × 
(PSF:24x24)
W社にいた2007年当時、福留治隆氏と約1年間シミュレーションを行い、相澤清晴先生の元で月1回発表し指導を受ける形で文書化したこの研究は、 特許申請、学会発表もせず、社内ですぐに利用できない技術であった。現在でも明らかにW社の知的財産であるが、既にそれから6年も経過し、 すでに同種の技術を偶に目にするようになった。すでに技術的に新規でないだろう。それゆえ、W社に感謝を込めて、ここにそのまま提示する。 いつまでも、書かない出さないでは、恐らくこれを誰もよく知らないだろうから、逆に皆に申し訳ない。文章は当時の私によるが筆頭著者を 福留氏にしたのは、彼が実際にLucy法を自力で再発見したからである。自分達の編みだした方法が本当に役にたつことに気が付き、 それが1970年代に発見された方法と同じであると知るのは、素晴しい驚きの経験であった。こういうことがあると仕事はやめられない。
最終ソース deconv149.tgz 2008/02/06 97kB
149のバグフィックス deconv150.tgz 2016/05/31 98kB Xの使用法にバグ、x.c を修正。Vine Linux 6.3で動作するようにフィックスしました。ご迷惑をお掛けして申しわけありません。(2016/5/31) この修正に反応された方に、jpg にするconvertが変わったようで、以前は大丈夫だったyuvの -size 指定が効かないです、すみません。 画像の変化だけはできますので、今はお許しを。jpg にする処理を別にしないとだめか？(2016/6/1)
Linux と X-window と C言語プログラムの開発環境(gcc)で動作するデコンボリューションプログラムです。コマンド convert を使用しますので、 使えることを確認してください。開発環境が異なる場合、Makefile を修正し、make が成功すれば、コマンド deconv ファイル名.拡張子 で動作 します。複数ファイルも処理します。(deconv -quiet *.JPG など) 拡張子がJPG又はjpgなら画像はJPEGとして入出力変換します。デコンボリューション 処理(階層型バイリニアPSF)を行います。YUVに変換した入力画像と、修正中の推定画像と、4隅のPSFをX-window表示しながら変化し、3階層x30回の くり返し後に、PSFの小画像をファイル名_psf.jpgとともに、推定画像を ファイル名_est.jpg として出力します。入力ファイルはそのまま残します。
最終ソースのdeconv.cは、deconv.c22 から輝度方向の非線型を自動的に補正することを入れようとした段階のソースで、そのため十分な性能がでない 場合が多いようです。c22 からコピーし、PSFサイズを deconv.h の #define H, #define V を画像に合わせて変更して make します。画角の小さい 画像では、4隅のPSFを仮定するのは無駄で、単一の PSFを仮定する deconv.c0 を deconv.c にコピーして開始して下さい。 その例、3階層30-30-30回の処理前後の画像を示します。この画像(640x480)の処理(PSF:24x24)は、通常のPC(AMD Athlon(tm) 64 X2 Dual Core Processor 5000+ )で 2分45秒 で結果を得ます。
手ぶれ画像 (640x480)＝
推定画像 × (PSF:24x24)
F(X,Y)= ∫PSF(x,y)Tree(X-x,Y-y)dxdy、∫PSF(k,l)dkdl ＝1
F(I,J)= ΣPSF(k,l)Tree(I-k,J-l)、ΣPSF(k,l)＝1
すき間から見える光源像は逆さではないのに、このX-x, Y-y がなぜマイナスかは、すき間から光源像の上の方を見る地面の場所は、少し下の 場所であるからである。そして面白いことに、我々が見る全ての画像は、元の本当の画像でなく、劣化した画像であると考えることができる。
オーディオのインパルス応答が過去にしか影響を受けないという因果性のため、インパルス応答が時刻が正の右片側にしか伸びないのと違って、 画像ではPSFは、上下左右の 2次元に拡がっている。ここまでは、オーディオのFIRフィルタの単純な拡張と考えられるが、ここで重要な違いが 導入される。PSFの各要素は非負に限定されることである(*)。そのため、高域通過フィルタHPFはないという特徴をもつ。画像のぼけは常に低域 通過フィルタLPFである。またそれからくる重要なこととして、画像劣化のPSFによるLPFフィルタの逆フィルタはHPFの空間周波数特性をもつ必要 があり、画像回復はPSFをもつFIRフィルタではできないということである。
第1の関門として、「劣化した観測画像とPSF画像とを与えて、元の画像を作れるか」という問題がある。答えは、ある方法が存在して、その解が あるかどうかはPSFによるという答えになるが、解が存在することが多い。そして Richardson-Lucy 法[1][2]というある逐次的方法によって解を 求めることが現実に可能であること、それを知ったとしても、第2の関門として、「観測画像だけから元の真の画像を得ることが可能か」という 問いがある。ここでは、真の画像と真のPSFを同時に求めることになる。これを Blind deconvolution という。
Blind deconvolution も Richardson-Lucy法で試すことができる。画像の推定の処理とPSFの推定の処理の繰り返しで多少、第1の関門より難しいが、 PSFなど天体画像でもない限り、通常は調べようがないから、この問題のほうが現実に多く存在し、そして、驚くべきことにこれも解が求まること が多いのである。
(*) PSFの要素が非負の値をもつことが Richardson-Lucy 法において絶対的に必要なことのように見えるのは、その観測画像と劣化画像の画素ごと の比を使うからである。光の古典的な波の理論においても、光が何かによって拡散されるとき、それは正のPSFをもつという制限は、過剰な制限 のように見えるが、現状ではこれ以外、あまり成功した方法はない。しかしながら、波の振幅の2乗が光の強さであるなら、波の振幅で補正しない限り、 その2乗である、強度の補正では限界があるだろう。そして、そういう拡張が可能になれば、オーディオにおいても Richardson-Lucy法による 逆フィルタが可能と考えられる。
まず、それまでの繰り返し法による画像回復の処理を、簡単に記述しよう。劣化フィルタをHとする。Hは、通常のフィルタである。PSFである h(k,l)による2次元FIRフィルタ処理によって推定画像Y(i,j)から劣化画像f(i,j)を作成する。clip(x)は xを[0,255]に制限する関数である。
H： f(i,j)= clip(Σ_kl Y(i-k,j-l) h(k,l))
観測画像x(i,j)と劣化画像f(i,j)との差e(i,j)にh(k,l)を掛けて推定画像Y(i-k,j-l)に足し込む。この処理が、r回目の推定画像から、r+1回の 推定画像を作成するとき、
H^T: Y^(r+1) (i-k,j-l)= Y^r (i-k,j-l) + e(i,j) h(k,l)
e(i,j)= x(i,j) - f(i,j)
この H^T 処理は、画像の周辺領域を例外として、通常のフィルタ処理に変形できる。
H^T： Y^(r+1) (i,j)= Y^r (i,j) + Σ_kl e(i+k,j+l) h(k,l)
e(i,j)= x(i,j) - f(i,j)
このようなとき、エラー画像の正画素について、PSFへの掛け込みを行うPSF修正が有効であった。
h^(r+1)(k,l)= h^r(k,l)*Σ_ij y(i-k,j-l) e(i,j) (e(i,j)＞0)
H^Tの改良として、h(k,l)に比例して推定画像 y(i-k,j-l) に e(i,j)h(k,l) を分配するとき、y(i-k,j-l)/f(i,j)をさらに掛けることが 有効であった。それは、明るい物体画像のリンギングが暗い領域に強く出ることを抑制したいという思いだけだったと、福留氏はいう。 それは、Richardson-Lucy 法として知られるデコンボリューションの方法の再発見であったことを後述する。
H^T: Y^(r+1) (i-k,j-l)= Y^r (i-k,j-l) + e(i,j) h(k,l) y(i-k,j-l)/f(i,j)
Φを観測画像、Ψを推定画像、pをPSF、χ、ξをそれぞれ2次元ベクトルとするとき、次の関係がある。
Φ(χ)＝∫Ψ(ξ)p(χ - ξ)dξ .......(1)
Ψ(ξ)＝∫Φ(χ)g(ξ - χ)dχ .......(2)
(1) は通常のPSFによる画像劣化を表す。(2)はその逆を行う画像回復である。しかし、(2)の画像回復のためのPSF g は未知である。 Φ、pからΨ、gを求めることが通常の逆フィルタ処理であり、Φだけから、Ψ、p、gを求めることが blind-deconvolution である。 (しかし、(2) の式を満たす、非負を要素にもつコンパクトな PSF である g は存在しない。)
画像Φ、Ψを0〜1の間をとる確率、PSFであるpとgを条件付き確率 p(χ - ξ)= P(χ|ξ), g(ξ - χ) = P(ξ|χ) として、Bayes の定理、 P(A｜B)＝ P(B｜A) P(A) / P(B) は、次式を導く。
g(ξ - χ)= p(χ - ξ) Ψ(ξ) / Φ(χ) .......(3)
(3)を(2)に代入時、分母のΦ(χ)は現在のΦ(χ)、分子のΦ(χ)に観測値Φ^0(χ)を使い、Ψ(ξ)への掛け込み漸化式を得る。 そして、PSFの漸化式も、画像とPSFの対称性から得る。
Ψ(ξ)＝Ψ(ξ) ∫Φ^0(χ)/Φ(χ) p(χ - ξ) dχ .......(5)
p(ξ)＝ p(ξ) ∫Φ^0(χ)/Φ(χ) Ψ(χ - ξ) dχ .......(6)
(5)を累加表現にすることができ、式(7) が福留の式に一致する。
Ψ(ξ)＝Ψ(ξ) + ∫(Φ^0(χ) - Φ(χ)) Ψ(ξ)/Φ(χ) p(χ - ξ) dχ .......(7)
PSFの更新の式(6)も、次の累加型漸化式(8)にすることができる。
p(ξ)＝ p(ξ) + ∫(Φ^0(χ) - Φ(χ)) Ψ(χ - ξ)/Φ(χ) p(ξ) dχ .......(8)
PSFの漸化式の別解として、(8)と異なる、次式(9)が有効であった。
p(ξ)＝ p(ξ) ∫ (Φ^0(χ) - Φ(χ)(＞0)) Ψ(χ-ξ) dχ .......(8)
この漸化の初期値は、PSF には平坦、Ψには入力観測画像Φ^0を与えた。推定画像の周辺境界の外に、PSF幅の1/2の領域には観測画像 の画素値を延長するか、領域の外の観測画像がある場合それを用いた。H＾Tの処理がフィルタ型でなく、分配型であるとき処理は使用 した画素に閉じている。
歴史的な文献[1],[2]は、現在、Web公開され一般に無料で入手可能である。(題名でGoogle検索して下さい。)
[1] William Hadley Richardson, "Bayesian-Based Iterative Method of Image Restoration", Jounal of the Optical Society of America Vol. 62, Num. 1. January 1972.
[2] L. B. Lucy, "An iterative technique for the rectification of observed distributions", The Astromical Jounal Vol. 79, Num. 6 June 1974, (American Astronomical Society. Provided by NASA Astrophysics Data System).