歪んだ画像を復元することは、画像処理タスクで最も興味深い重要な問題の1つです。理論的および実用的な観点からです。 特殊なケースは、不適切な焦点とぼやけのためにぼやけています-これらの欠陥は、あなたがよく知っているものですが、修正するのは非常に困難です-それらは記事のトピックによって選択されたものです。 残りの歪み（ノイズ、不適切な露出、歪み）により、人類は効果的に対処することを学び、適切なツールがすべての自尊心のある写真エディターにあります。

なぜグリースと焦点ぼけを除去するために実質的に何もないのか（不鮮明なマスクはカウントされない）-これは原則的に不可能でしょうか？ 実際、それは可能です-対応する数学的装置は約70年前に開発され始めましたが、他の多くの画像処理アルゴリズムと同様に、これらすべては最近になって広く使用されています。 ここで、すごい効果のデモンストレーションとして、いくつかの写真：



私は拷問されたレナを使わなかったが、ヴェネツィアの私の写真を見つけた。 右の画像は、左から正直に取得され、48ビット形式のようなトリックを使用することなく（この場合、元の画像が100％復元されます）-左は人為的にぼやけた最も一般的なPNGです。 結果は印象的ですが...実際にはそれほど単純ではありません。 猫の下で理論と実際の結果の詳細なレビュー。
注意、たくさんのPNG画像！

はじめに

遠くから始めましょう。 多くの人々は、ぼかしは不可逆的な操作であり、情報は完全に失われると信じています。 各ピクセルがスポットに変わり、すべてが混ざり合い、大きなぼかし半径で画像全体に均一な色が得られます。 これは完全に真実ではありません-すべての情報はいくつかの法律に従って単に再配布され、いくつかの留保によって明確に復元することができます。 例外は、ぼかし半径の幅を持つ画像のエッジのみです-完全な復元はそこでは不可能です。

1次元の場合の小さな例を使用して、これを「指で」実演します。値を持つピクセルの行があると想像してください。
x ₁ | x ₂ | x ₃ | x ₄ ...-オリジナル画像

歪みの後、各ピクセルの値は左の値と合計されます。 x ' _i = x _i + x _i-1 。 理論的には、2で割る必要がありますが、簡単にするためにこれを省略します。 その結果、ピクセル値を持つぼやけた画像ができます：
x ₁ + x ₀ | x ₂ + x ₁ | x ₃ + x ₂ | x ₄ + x ₃ ...-ぼやけた画像

次に、スキームに従ってチェーンに従って値を順次減算して復元しようとします-1番目は2番目のピクセルから、2番目は3番目から、3番目は4番目から、というようになります。
x ₁ + x ₀ | x ₂ -x ₀ | x ₃ + x ₀ | x ₄ -x ₀ ...-復元された画像

その結果、ぼやけた画像の代わりに、元の画像が得られ、そのピクセルには、符号が交互に変わる未知の定数x ₀が追加されました。 これはすでにはるかに優れています-この定数は視覚的に選択でき、値x ₁にほぼ等しいと仮定でき、隣接ピクセルの値ができる限り「ジャンプ」するなどの基準で自動的に選択できます。 ただし、ノイズを追加するとすぐにすべてが変化します（常に実際の画像にあります）。 説明したスキームでは、各ステップで、総成分へのノイズの寄与が蓄積され、最終的には完全に受け入れられない結果をもたらす可能性がありますが、これまで見てきたように、復元はそのような原始的な方法でも非常に現実的です。

歪みプロセスモデル

それでは、これらの歪みと復元のプロセスについて、より正式で科学的な説明に移りましょう。 フルカラー画像を処理するには、RGBカラーチャンネルごとに必要なすべてのステップを繰り返すだけで十分であるという仮定の下で、グレースケールの白黒画像のみを考慮します。 次の表記法を紹介します。
f（x、y） -歪みのない元の画像
h（x、y）は歪み関数です
n（x、y） -加法性ノイズ
g（x、y）は歪みの結果です。 結果として観察されるもの（ぼやけた画像またはピンぼけした画像）

歪みプロセスのモデルを次のように定式化します。
g（x、y）= h（x、y）* f（x、y）+ n（x、y） （1）

歪んだ画像を復元するタスクは、元の画像の最適な近似f '（x、y）を見つけることです。 各コンポーネントをより詳細に検討してください。 f（x、y）とg（x、y）では、すべてが明確です。 しかし、関数h（x、y ）については、いくつかの言葉を言う必要があります-どんな感じですか？ 歪みの過程で、元の画像の各ピクセルは、ピンぼけの場合はスポットになり、単純なグリースの場合はセグメントになります。 または、反対に、歪んだ画像の各ピクセルは元の画像の特定の近傍のピクセルから「組み立てられる」ということもできます。 これらはすべて互いに重なり合っており、その結果、画像が歪んでしまいます。 これは、1つのピクセルがスメアまたは収集される法則に従って、歪み関数と呼ばれます。 他の同義語は、PSF（点広がり関数、つまり点の分布関数）、歪曲演算子の核、カーネルなどです。 この関数の次元は、原則として、画像自体の次元よりも小さいです。たとえば、「指で」の例の最初の検査では、関数の次元は2でした。 各ピクセルは2つで構成されていました。

歪み関数

典型的な歪み関数がどのように見えるか見てみましょう。 以下では、Matlabツールを使用します.Matlabツールは、そのような目的のために既に標準になっています。これには、さまざまな画像処理実験に必要なすべてのものが含まれており（だけでなく）、アルゴリズム自体に集中して、すべてのルーチン作業を関数ライブラリにシフトできます。 ただし、このためにはパフォーマンスを支払う必要があります。 そこで、PSFに戻って、その種類の例を示します。


fspecial関数によるガウスぼかしの場合のPSF（「ガウス」、30、8）。


fspecial機能（「モーション」、40、45）を使用した潤滑の場合のPSF。

歪み関数を別の関数（この場合は画像）に適用する操作は、畳み込みと呼ばれます。 元の画像の一部の領域は、歪んだ画像の1ピクセルに縮小されます。 通常の乗算​​と混同しないように、演算子「*」で指定します！ 数学的には、次元M x Nの画像fと次元mxnの歪み関数hに対して、次のように記述されます。

（2）

ここで、a =（m-1）/ 2、b =（n-1）/ 2です。 畳み込みの反対の操作はデコンボリューションと呼ばれ、このような問題の解決は非常に重要です。

ノイズモデル

ノイズの原因となる最後の項、式（1）のn（x、y）を考慮する必要があります。 デジタルセンサーのノイズの原因は大きく異なる可能性がありますが、主な原因は熱変動と暗電流です。 ノイズ値は、ISO値、マトリックスタイプ、ピクセルサイズ、温度、電磁干渉などの多くの要因の影響も受けます。ほとんどの場合、ノイズはガウス（平均と分散の2つのパラメーターによって決定されます）であり、相関もありません画像とともに、ピクセルの座標に依存しません。 最後の3つの仮定は、今後の作業にとって非常に重要です。

畳み込み定理

回復問題の初期定式化に戻りましょう-ノイズを忘れずに、何らかの方法で畳み込みを逆にする必要があります。 式（2）から、 g （x、y）からf（x、y ）を取得することはそれほど容易ではないことがわかります。「正面」と呼ばれるものを決定すると、巨大な方程式系が得られます。 しかし、 フーリエ変換は私たちの助けになります。詳細については説明しません。このトピックについてはすでに多くのことが述べられています。 そのため、空間領域での畳み込み演算は、周波数領域での通常の乗算​​（さらに、行列乗算ではなく、要素ごとの乗算）と同等であると述べるような畳み込み定理があります。 したがって、逆畳み込み演算は、周波数領域での除算と同等です。つまり、次のように記述できます。
（3）

ここで、 H（u、v）、F（u、v）は対応する関数のフーリエ画像です。 したがって、式（1）の歪みプロセスは、周波数領域で次のように書き換えることができます。
（4）

逆フィルタリング

すぐにこの等式をH（u、v）に分割し、元の画像の次の推定値F ^{^} （u、v）を取得するよう求めます。
（5）
これは逆フィルタリングと呼ばれますが、実際にはほとんど機能しません。 なぜそうですか？ この質問に答えるために、式（5）の最後の項を見てみましょう。関数H（u、v）がゼロまたはゼロに近い値を取る場合、この項の寄与が支配的になります。 これはほとんどの場合、実際の例で見られます-これを説明するために、フーリエ変換後のスペクトルがどのように見えるか思い出してみましょう。

元の画像を撮り、


それをグレースケールに変換し、Matlabを使用して、スペクトルを取得します。

% Load image I = imread('image_src.png'); figure(1); imshow(I); title(' '); % Convert image into grayscale I = rgb2gray(I); % Compute Fourier Transform and center it fftRes = fftshift(fft2(I)); % Show result figure(2); imshow(mat2gray(log(1+abs(fftRes)))); title('FFT -   ( )'); figure(3); imshow(mat2gray(angle(fftRes))); title('FFT -  '); 

その結果、振幅スペクトルと位相スペクトルの2つのコンポーネントが得られます。 ちなみに、多くの人はこの段階を忘れています。 振幅スペクトルは、対数目盛で示されていることに注意してください。 その値は大きく異なります-中心は数百万の最大値であり（数百万のオーダー）、中心から遠ざかるにつれて急速にほぼゼロまで減少します。 このため、逆フィルタリングはゼロまたはほぼゼロのノイズ値でのみ機能します。 次のスクリプトを使用して、実際にこれを示しましょう。

 % Load image I = im2double(imread('image_src.png')); figure(1); imshow(I); title(' '); % Blur image Blurred = imfilter(I, PSF,'circular','conv' ); figure(2); imshow(Blurred); title(' '); % Add noise noise_mean = 0; noise_var = 0.0; Blurred = imnoise(Blurred, 'gaussian', noise_mean, noise_var); % Deconvolution figure(3); imshow(deconvwnr(Blurred, PSF, 0)); title(''); 

  noise_var = 0.0000001 noise_var = 0.000005 

非常に小さなノイズの追加でさえ、著しい干渉につながり、それがこの方法の実際の応用を大きく制限することが明らかに見られます。

デコンボリューションの既存のアプローチ

しかし、画像内のノイズの存在を考慮に入れるアプローチがあります-最も有名なものの1つは、Wienerフィルターです。 彼は、画像とノイズをランダムプロセスと見なし、これらの値の標準偏差が最小になるように、歪みのない画像fの推定値f 'を見つけます。 この偏差の最小値は、周波数領域の関数で達成されます。
（6）
この結果は、1942年にWienerによって得られました。 ここでは詳細な結論を述べませんが、興味がある人はここでそれを見ることができます 。 ここでの関数Sは、それぞれノイズと元の画像のエネルギースペクトルを示しています-これらの量はほとんど知られていないため、S _n / S _{fの部分は} 、信号対雑音比としてほぼ記述できる定数Kに置き換えられます。

次の方法は、「リンクによる平滑化最小二乗フィルタリング」、他の名前：「Tikhonovフィルタリング」、「Tikhonov正則化」です。 彼の考えは、対応する最適化問題のさらなる解決策を用いて、問題を行列形式で定式化することです。 このソリューションは次のように書かれています。
（7）
ここで、 yは正則化パラメーター、 P（u、v）はラプラス演算子（3 * 3行列）のフーリエ変換です。

別の興味深いアプローチは、リチャードソン（1972）とルーシー（1974）によって独立して提案されました。 このメソッドは、「ルーシー・リチャードソンメソッド」と呼ばれます。 その際立った特徴は、最初の3つとは異なり、非線形であることです。これにより、最高の結果が得られる可能性があります。 2番目の機能-メソッドは反復的であるため、反復を停止するための基準で問題が発生します。 基本的な考え方は、画像がポアソン分布に従うと想定される最尤法を使用することです。 計算式は、フーリエ変換を使用せずに非常に簡単です-すべてが空間領域で行われます：
（8）
ここで、シンボル「*」は、前述のように、畳み込み演算を示します。 この方法は、天文写真を処理するプログラムで広く使用されています。デフォトコンボリューション（フォトエディタのようにシャープでないマスクの代わりに）を使用することが事実上の標準です。 例はAstra Imageです。ここにデコンボリューションの例があります。 この方法の計算の複雑さは非常に大きく、反復回数に応じて平均的な写真を処理するのに、何時間も、さらには何日もかかることがあります。

検討中の最後のメソッド、またはむしろ、現在積極的に開発および開発されているメソッドのファミリー全体は、ブラインドデコンボリューションです。 以前のすべての方法では、PSFの歪み関数は正確に既知であると想定されていましたが、実際にはそうではなく、通常、PSFは可視歪みの性質によってのみ既知です。 ブラインドデコンボリューションは、これを考慮に入れるための単なる試みです。 原則は非常に簡単です。詳細に入らない場合、PSFの最初の近似が選択され、デコンボリューションはいずれかの方法を使用して実行されます。その後、品質レベルが何らかの基準によって決定され、それに基づいてPSF関数が洗練され、目的の結果が得られるまで反復が繰り返されます

練習する

さて、理論で、すべて-練習に移りましょう、画像にリストされた方法を人工的なぼかしとノイズと比較することから始めましょう。

 % Load image I = im2double(imread('image_src.png')); figure(1); imshow(I); title(' '); % Blur image PSF = fspecial('disk', 15); Blurred = imfilter(I, PSF,'circular','conv' ); % Add noise noise_mean = 0; noise_var = 0.00001; Blurred = imnoise(Blurred, 'gaussian', noise_mean, noise_var); figure(2); imshow(Blurred); title(' '); estimated_nsr = noise_var / var(Blurred(:)); % Restore image figure(3), imshow(deconvwnr(Blurred, PSF, estimated_nsr)), title('Wiener'); figure(4); imshow(deconvreg(Blurred, PSF)); title('Regul'); figure(5); imshow(deconvblind(Blurred, PSF, 100)); title('Blind'); figure(6); imshow(deconvlucy(Blurred, PSF, 100)); title('Lucy'); 

結果：


ウィナーフィルター


チホノフ正則化


ルーシーリチャードソンフィルター


ブラインドデコンボリューション

おわりに

そして、最初の部分の最後で、実際の画像の例を少し触れます。 それ以前は、すべての歪みは人為的なものでした。もちろん、走って勉強するのに適していますが、実際の写真でこれがどのように機能するかを見るのは非常に興味深いです。 Canon 500D DSLRで撮影したこのような写真の一例を次に示します。



次に、簡単なスクリプトを実行します。

 % Load image I = im2double(imread('IMG_REAL.PNG')); figure(1); imshow(I); title(' '); %PSF PSF = fspecial('disk', 8); noise_mean = 0; noise_var = 0.0001; estimated_nsr = noise_var / var(I(:)); I = edgetaper(I, PSF); figure(2); imshow(deconvwnr(I, PSF, estimated_nsr)); title(''); 

そして、次の結果が得られます。



ご覧のように、画像に新しい詳細が表示され、明瞭度ははるかに高くなりましたが、干渉は対照的な境界で「リンギング」の形で現れました。

実際のグリースを使用した例-その実装では、カメラを三脚に取り付け、比較的遅いシャッタースピードを設定し、シャッターを放した瞬間に均一な動きが得られました。


スクリプトはほぼ同じで、PSFタイプのみが「モーション」になります。

 % Load image I = im2double(imread('IMG_REAL_motion_blur.PNG')); figure(1); imshow(I); title(' '); %PSF PSF = fspecial('motion', 14, 0); noise_mean = 0; noise_var = 0.0001; estimated_nsr = noise_var / var(I(:)); I = edgetaper(I, PSF); figure(2); imshow(deconvwnr(I, PSF, estimated_nsr)); title(''); 

結果：


再び、品質が著しく向上しました-窓と車のフレームが区別可能になりました。 アーティファクトは、以前のデフォーカスされた例とはすでに異なっています。

これで面白いと最初の部分を終了します。
第2部では、実際の画像の処理の問題に焦点を当てます-PSFの構築と評価、より複雑で高度なデコンボリューション技術、リンギングなどの欠陥の除去方法、既存のソフトウェアのレビューと比較などを検討します。

PS少し前に、新しいバージョンのPhotoshopでぼやけた写真を修正することに関する記事がハブで公開されました
グリースを除去するための同様の技術（おそらくPhotoshopで使用されるもの）を使用したい場合は、 このリンクを使用してアプリケーションのデモバージョンをダウンロードし、回復例を参照し、操作の原則について読んでください。

文学

ゴンザレスR.、ウッズR.デジタル画像処理
Gonzalez R.、Woods R.、Eddins S.MATLABでのデジタル画像処理

UPD： 続行するリンク

  - 

ウラジミール・ユジコフ