👟 😘 🖐🏽 ぼやけた画像を回復する別のアルゴリズム 👦🏻 ⛵️ 🏎️

良い一日。画像のデコンボリューション方法については多くのことが言われていますが、追加するものは他にないようです。ただし、以前のアルゴリズムよりも優れた新しいアルゴリズムが常に存在します。少し前まで、低メモリコストで線形収束率を持ち、安定しており、並列化が可能な反復アルゴリズムが説明されていました。そしてしばらくすると、二次収束にまで改善されました。会う：（高速）反復収縮しきい値アルゴリズム。

問題の声明

画像のデコンボリューションアルゴリズムに関する長いシリーズの記事が公開されました。サイクルは、線形代数方程式のシステムとそれらの正則化を解くための古典的な方法を考慮しました。問題のステートメントの観点から説明した方法は、実際には古典的な方法と違いはありませんが、悪魔は詳細にあります。数学の魔法を始めましょう。これは、受信した破損したノイズの多い画像の古典的な表現です。

m a t h b f y = m a t h b f A x + m a t h b f n

$\ mathbf {y} = \ mathbf {Ax} + \ mathbf {n}$

マトリックス

m a t h b f A

$\ mathbf {A}$ は畳み込み行列であり、物理的な観点から、元のデータと破損したデータのポイント間の関係を示します。画像の場合、この行列は畳み込みカーネルを変換することで取得できます。ベクトル

m a t h b f x

$\ mathbf {x}$ そして

m a t h b f y

$\ mathbf {y}$ ベクトル化された形式の元の破損した画像です。つまり、画像内のすべてのピクセル列が次々に連結されます。

この問題の定式化の最初の問題：割り当てられたメモリの量。画像のサイズが640x480の場合、タスクは、サイズが307200x1の要素とサイズが307200x307200要素の行列の画像の2つのベクトルを保存します。畳み込み行列はスパースですが、そのサイズは、スパース形式であっても大きくなり、多くのメモリを占有します。畳み込み行列を含む積は、カーネルを使用した2次元畳み込みに置き換えることができます。これにより、メモリに大きな行列を保存する必要がなくなり、行列を保存する問題が解決します。

m a t h b f A

$\ mathbf {A}$ 。
元の画像にできるだけ近い画像を検索するには、損傷した画像と目的の画像の間の残差の2番目のノルムの古典的な形式で最適化された関数を記述します。

m i n_{h a t m a t h b f x} m i d m i d m a t h b f y - A h a t x m i d m i d_{2}

$min _ {\ hat {\ mathbf {x}}} \ mid \ mid \ mathbf {y-A \ hat {x}} \ mid \ mid_2$

従来の方法とは異なり

この問題に関する声明は最も人気のあるものの1つですが、パフォーマンスを改善するために補足します。使用されるアプローチはマジョライゼーション最大化と呼ばれ、元のタスクを非常によく似た別のタスクに置き換えることで構成されます。新しいタスクには2つのプロパティがあります。

タスクははるかに簡単です
現在の点を除くすべての点で、残差は元の問題より大きい

このアクションは、アルゴリズムの各反復で発生します。とても複雑に聞こえますが、実際に見えるよりもずっと複雑です。結果の反復の最小化関数

k

$k$ 次のように書かれています。

m a t h b f J_{k} （ x ） = m i d m i d m a t h b f y A h a t x m i d m i d_{2} + （ m a t h b f h a t x - h a t x_{k} ）^{T} （ \ア ル フ ァ m a t h b f I - m a t h b f A^{T} A ） （ m a t h b f h a t x - h a t x_{k} ）

$\ mathbf {J_k（x）} = \ mid \ mid \ mathbf {yA \ hat {x}} \ mid \ mid_2 +（\ mathbf {\ hat {x}-\ hat {x_k}}）^ T（\アルファ\ mathbf {I}-\ mathbf {A ^ TA}）（\ mathbf {\ hat {x}-\ hat {x_k}}）$

右項の行列は正定です。これは、

a l p h a > m a x e i g （ m a t h b f A^{T} A ）

$\ alpha> maxeig（\ mathbf {A ^ TA}）$ 。最小化関数は、それぞれが非負の2つの量の合計で構成されます。また、その時点で

m a t h b f h a t x = h a t x_{k}

$\ mathbf {\ hat {x} = \ hat {x} _k}$ 2番目の項はゼロに等しいため、リストの2番目の条件が満たされます。

ソリューションを検索する

解決策の検索は、古典的な方法で開始します。括弧を開き、勾配を書き、ゼロに等しくします。これは上の勾配であることに注意してください

k

$k$ その反復。たくさんの数学に注意してください！

m a t h b f J_{k} （ h a t x ） = （ m a t h b f y A h a t x ）^{T} （ m a t h b f y A h a t x ） + （ m a t h b f h a t x - h a t x_{k} ）^{T} c d o t （ a l p h a m a t h b f I - m a t h b f A^{T} A ） c d o t （ m a t h b f h a t x - h a t x_{k} ）

$\ mathbf {J_k（\ hat {x}）} =（\ mathbf {yA \ hat {x}}）^ T（\ mathbf {yA \ hat {x}}）+（\ mathbf {\ hat {x} -\ hat {x} _k}）^ T \ cdot（\ alpha \ mathbf {I}-\ mathbf {A ^ TA}）\ cdot（\ mathbf {\ hat {x}-\ hat {x} _k}）$

m a t h b f J_{k} （ h a t x ） = m a t h b f y^{T} y + m a t h b f h a t x^{T} A^{T} A h a t x - 2 m a t h b f y^{T} A h a t x - - m a t h b f h a t x^{T} （ a l p h a m a t h b f I - m a t h b f A^{T} A ） h a t x + m a t h b f h a t x_{k}^{T} （ a l p h a m a t h b f I - m a t h b f A^{T} A ） h a t x_{k} - 2 m a t h b f h a t x_{k}^{T} （ a l p h a m a t h b f I - m a t h b f A^{T} A ） h a t x

$\ mathbf {J_k（\ hat {x}）} = \ mathbf {y ^ Ty} + \ mathbf {\ hat {x} ^ TA ^ TA \ hat {x}} -2 \ mathbf {y ^ TA \ hat {x}}-\\-\ mathbf {\ hat {x} ^ T（\ alpha \ mathbf {I}-\ mathbf {A ^ TA}）\ hat {x}} + \ mathbf {\ hat {x} _k ^ T（\ alpha \ mathbf {I}-\ mathbf {A ^ TA}）\ hat {x} _k} -2 \ mathbf {\ hat {x} _k ^ T（\ alpha \ mathbf {I}-\ mathbf {A ^ TA}）\ hat {x}}$

f r a c m a t h b f J_{k} （ h a t x ） d e l t a m a t h b f h a t x = m a t h b f A^{T} A h a t x - 2 m a t h b f A^{T} y + m a t h b f 2 a l p h a m a t h b f I h a t x - m a t h b f A^{T} A h a t x - 2 m a t h b f （ a l p h a m a t h b f I - m a t h b f A^{T} A ） h a t x_{k}

$\ frac {\ mathbf {J_k（\ hat {x}）}} {\ delta \ mathbf {\ hat {x}}} = \ mathbf {A ^ TA \ hat {x}} -2 \ mathbf {A ^ Ty} + \ mathbf {2 \ alpha \ mathbf {I} \ hat {x}-\ mathbf {A ^ TA} \ hat {x}}-2 \ mathbf {（\ alpha \ mathbf {I}-\ mathbf { A ^ TA}）\ hat {x} _k}$

その結果、デリバティブは次の形式になります。

d e l t a J_{k} （ x ） = 2 m a t h b f A^{T} y - 2 （ a l p h a m a t h b f I A^{T} A ） m a t h b f h a t x_{k} + 2 a l p h a m a t h b f h a t x

$\ delta J_k（x）= 2 \ mathbf {A ^ Ty} -2（\ alpha \ mathbf {IA ^ TA}）\ mathbf {\ hat {x} _k} +2 \ alpha \ mathbf {\ hat {x }}$

この状況での次の古典的なステップは、勾配をゼロに等しくし、結果の式から目的のベクトルを表現することです

m a t h b f h a t x

$\ mathbf {\ hat {x}}$ 。記録済み

m a t h b f h a t x

$\ mathbf {\ hat {x}}$ として定義する

m a t h b f h a t x_{n + 1}

$\ mathbf {\ hat {x} _ {n + 1}}$ または、次の反復のソリューションとして。

d e l t a J_{k} （ x ） = 2 m a t h b f A^{T} y - 2 （ a l p h a m a t h b f I A^{T} A ） m a t h b f h a t x_{k} + 2 a l p h a m a t h b f h a t x = 0

$\ delta J_k（x）= 2 \ mathbf {A ^ Ty} -2（\ alpha \ mathbf {IA ^ TA}）\ mathbf {\ hat {x} _k} +2 \ alpha \ mathbf {\ hat {x }} = 0$

m a t h b f A^{T} y + a l p h a m a t h b f I m a t h b f h a t x_{k} - m a t h b f A^{T} A m a t h b f h a t x_{k} = a l p h a m a t h b f h a t x

$\ mathbf {A ^ Ty} + \ alpha \ mathbf {I} \ mathbf {\ hat {x} _k}-\ mathbf {A ^ TA} \ mathbf {\ hat {x} _k} = \ alpha \ mathbf { \ hat {x}}$

m a t h b f h a t x_{k + 1} = m a t h b f h a t x_{k} + f r a c 1 a l p h a c d o t m a t h b f A^{T} （ y A h a t x_{k} ）

$\ mathbf {\ hat {x} _ {k + 1}} = \ mathbf {\ hat {x} _k} + \ frac {1} {\ alpha} \ cdot \ mathbf {A ^ T（yA \ hat { x_k}）}$

結果の式は「Landweber iteration」と呼ばれます。これは、反復間の基本的な式です。これを使用すると、線形速度での反復ごとの残差の減少を保証できます。

基本的なアルゴリズムには、「ソフトしきい値」と呼ばれる追加のステップが含まれており、ソリューションのスパース性を想定しています。このステップは、設定された値より小さいモジュロを求めているベクトルのすべての値をゼロに等しくします。これにより、回復結果に対するノイズの影響が軽減されます。

それはどのように見えますか

このソリューションからわかるように、好みに合わせて調整する2つのパラメーターがあります。

\ラ ム ダ

$\ラムダ$ 近似の精度に責任を負い、収束をしきい値に制限します。

a l p h a

$\ alpha$ 作業の安定性とアルゴリズムの収束速度を担当します。

m a t h b f h a t x_{k + 1} = s o f t （ m a t h b f h a t x_{k} + f r a c 1 a l p h a c d o t m a t h b f A^{T} （ y A h a t x_{k} ） 、 l a m b d a / a l p h a ）

$\ mathbf {\ hat {x} _ {k + 1}} = soft（\ mathbf {\ hat {x} _k} + \ frac {1} {\ alpha} \ cdot \ mathbf {A ^ T（yA \ hat {x_k}）}、\ lambda / \ alpha）$

アルゴリズムの変更

その後、アルゴリズムは、共役勾配法に理論的に類似した追加の用語によって改善されました。各ステップで、前の2回の反復の結果の差を追加して、アルゴリズムの収束を2次に増やします。アルゴリズムの最終手順は次のとおりです。

パラメータを設定する $\ lambda、\ alpha$
セット $t_1 = 1、t_0 = 1、\ mathbf {y_ {temp} = \ hat {x} _0 = A ^ Ty}$
基本的なアルゴリズムを繰り返す $\ mathbf {\ hat {x} _ {k + 1}} = soft（\ mathbf {y_ {temp}} + \ frac {1} {\ alpha} \ cdot \ mathbf {A ^ T（y-Ay_ { temp}）}、\ lambda / \ alpha）$
一時変数を更新する
$t_0 = t_1; t_1 = 0.5 + \ sqrt {0.25 + t_1 ^ 2}$
共役成分を変数に追加します $\ mathbf {y_ {temp} = \ hat {x} _k} + \ frac {t_0-1} {t_1} \ cdot（\ mathbf {\ hat {x} _k- \ hat {x} _ {k-1 }}）$

次に、アルゴリズムの並列性について説明します。各反復でしきい値を計算するサイクルは、OpenCLとCUDAの両方に含まれる要素ごとの積（アダマールの積）と比較演算で表現できるため、簡単に並列化できます。

CPUとGPUでの計算用の2つのバージョンでMatlabにアルゴリズムを実装しました。将来的には、GPUの方が便利なため、GPUのバージョンのコードのみを使用し始めました。コードを以下に示します。

ここにコード

function [x] = fista_gpu(y,H,Ht,lambda,alpha,Nit) x=Ht(y); y_k=x; t_1=1; T=lambda/alpha; for k = 1:Nit x_old=x; x=soft_gpu((y_k+(1/alpha)*Ht(yH(y_k))), T); t_0=t_1-1; t_1=0.5+sqrt(0.25+t_1^2); y_k=x+(t_0/t_1)*(x-x_old); end end function [y] = soft_gpu(x,T) eq=ge(abs(x),abs(T)); y=eq.*sign(x).*(abs(x)-abs(T)); end

練習に近づき、画像のサイズに応じてアルゴリズムのパフォーマンスをテストしましょう。

私のラップトップ構成

IntelコアI3-4005U
NVIDIA GeForce 820M

左側は、CPUとGPUのピクセル数に応じた、アルゴリズムの100回の繰り返しに対する依存性のグラフです。

結果からわかるように、GPUアルゴリズムはタスクのサイズに実質的に依存せず、非常に大きな画像の場合、GPUメモリによってのみ制限されます。 2000x2000の0.5秒（メモリからのアンロードを除く）はまともな結果になりますか？
右側は、アルゴリズムが繰り返し時間に費やした時間の依存性です。反復回数が増えると、GPUのアルゴリズムはランタイムの増加に反応し始めます。最も可能性が高いのは、コードの記述エラー、またはGPUの周波数を強制的に下げることによるものです。ほとんどの場合、100回の反復で十分です。

次のグラフは、基本アルゴリズムと改善されたアルゴリズムの比較、およびアルゴリズムの収束の依存性を示しています。

a l p h a

$\ alpha$ そして

\ラ ム ダ

$\ラムダ$ パラメータ。

左側のグラフからわかるように、改良されたアルゴリズムは基本的なアルゴリズムよりもはるかに速く収束し、100回の反復後、実際には精度の向上は見られません。 500回目の反復による基本的なアルゴリズムのみが、許容できる結果を示しています。
中央のグラフは、 $\ラムダ$ アルゴリズムの収束のしきい値が変更され、パフォーマンスが制限されます。信号が非常にノイズが多い場合、これは必要な妥協案です。
右側のグラフは、 $\ alpha$ アルゴリズムははるかにゆっくり収束します。前のケースと同様に、畳み込み行列の「品質」と収束率の間には妥協点があります。

そして最後に、FHD画像のアルゴリズムの例、

ソース画像

甘やかされて育った写真

復元された画像

ご覧のとおり、結果は非常に良好です。特に15秒しかかからない場合は、1.5がアルゴリズムの操作で、13.5がGPUメモリから画像をアンロードします。
アルゴリズムに使用されるすべてのコードはGitHubにアップロードされます。

中古文学

I.セレズニック。（2009）概要：スパース信号の復元。
A. BeckおよびM. Teboulle制約付き全変動画像のノイズ除去およびブレ除去の問題に対する高速勾配ベースのアルゴリズムIEEEトランザクションオンイメージプロセッシング、VOL。 18、NO。 2009年11月11日
PL CombettesおよびJ.-C. 信号処理におけるペスケ近位分割法

ぼやけた画像を回復する別のアルゴリズム

問題の声明

従来の方法とは異なり

ソリューションを検索する

アルゴリズムの変更

中古文学

More articles: