ルーシー・リチャードソン アルゴリズムを使用したイメージのブレ除去

手順 1: イメージの読み取り

この例では RGB イメージを読み取り、256 X 256 X 3 にトリミングします。関数 deconvlucy は任意の次元の配列を処理できます。

I = imread("board.tif");
I = I(50+(1:256),2+(1:256),:);
figure;
imshow(I);
title("Original Image");
text(size(I,2),size(I,1)+15, ...
    "Image courtesy of courtesy of Alexander V. Panasyuk, Ph.D.", ...
    "FontSize",7,"HorizontalAlignment","right");
text(size(I,2),size(I,1)+25, ...
    "Harvard-Smithsonian Center for Astrophysics", ...
    "FontSize",7,"HorizontalAlignment","right");

手順 2: ブレとノイズのシミュレーション

カメラの移動や焦点が合わないためにブレる可能性がある現実のイメージをシミュレートします。また、ランダム外乱によるノイズがイメージに含まれる場合もあります。この例では、実際のイメージをガウス フィルターで (imfilter を使用して) 畳み込むことによって、ブレをシミュレートします。ガウス フィルターは点像分布関数 PSF で表されます。

PSF = fspecial("gaussian",5,5);
Blurred = imfilter(I,PSF,"symmetric","conv");
figure;
imshow(Blurred);
title("Blurred");

この例では、分散のガウス ノイズ V を、ブレを含んだイメージに (imnoise を使用して) 追加することでノイズをシミュレートします。ノイズ分散 V は、後で、アルゴリズムの減衰パラメーターの定義に使用されます。

V = .002;
BlurredNoisy = imnoise(Blurred,"gaussian",0,V);
figure;
imshow(BlurredNoisy);
title("Blurred & Noisy");

手順 3: ブレとノイズを含むイメージの復元

PSF を指定し、反復回数 5 (既定値は 10) のみを使用して、ブレとノイズを含むイメージを復元します。出力は、入力イメージと同じタイプの配列です。

luc1 = deconvlucy(BlurredNoisy,PSF,5);
figure;
imshow(luc1);
title("Restored Image, NUMIT = 5");

手順 4: 復元を調査するための反復

結果として得られるイメージは、各反復と共に変化します。イメージの復元の変化を調べるには、逆畳み込みを行います。一連の反復を行い、その結果を確認したら、中止した場所から反復を再開します。これを行うには、入力イメージを cell 配列の一部として渡す必要があります。たとえば、入力イメージ パラメーターとして BlurredNoisy の代わりに {BlurredNoisy} を渡すことで、反復の最初のセットを開始します。

luc1_cell = deconvlucy({BlurredNoisy},PSF,5);

その場合、出力 luc1_cell は cell 配列になります。セル出力は 4 つの数値配列から構成されます。その 1 つ目の配列は BlurredNoisy イメージ、2 つ目の配列は double クラスの復元されたイメージ、3 つ目の配列は最後の 1 つ前の反復の結果、4 つ目の配列は反復されたセットの内部パラメーターです。出力 cell 配列の 2 つ目の数値配列であるイメージ luc1_cell{2} は手順 3 の出力配列であるイメージ luc1 と同じですが、クラスに例外が発生する可能性があります (セルの出力は、必ずクラス double の復元されたイメージです)。

反復を再開するには、前の関数呼び出しからの出力である cell 配列 luc1_cell を取得し、関数 deconvlucy に渡します。既定の反復回数 (NUMIT = 10) を使用します。復元されたイメージは、合計 15 回の反復の結果です。

luc2_cell = deconvlucy(luc1_cell,PSF);
luc2 = im2uint8(luc2_cell{2});
figure;
imshow(luc2);
title("Restored Image, NUMIT = 15");

手順 5: 減衰によるノイズ増幅の制御

最新のイメージ luc2 は、15 回の反復の結果です。5 回の反復時点の前の結果に比べてシャープになっていますが、イメージに "染みが付いた" ように見えます。染みは実際のどの構造にも対応しませんが (実際のイメージと比較して)、データ内のノイズをできる限り近似しようとする場合に生じます。

ノイズ増幅を制御するには、DAMPAR パラメーターを指定して減衰オプションを使用します。DAMPAR は入力イメージと同じクラスでなければなりません。アルゴリズムにより、ノイズに比べて差異が小さい領域でのモデルの変化が抑制されます。この例で使用される DAMPAR は、ノイズの 3 つの標準偏差に等しくなります。イメージがよりスムーズであることに注目してください。

DAMPAR = im2uint8(3*sqrt(V));
luc3 = deconvlucy(BlurredNoisy,PSF,15,DAMPAR);
figure;
imshow(luc3);
title("Restored Image with Damping, NUMIT = 15");

この例の以下の部分では、(簡略化と高速化のために) シミュレーションを実行した星のイメージを使用して、関数 deconvlucy の入力パラメーター weight および subsample について考察します。

手順 6: サンプル イメージの作成

以下の例では、4 つの星の白黒イメージを作成します。

I = zeros(32);
I(5,5) = 1;
I(10,3) = 1;
I(27,26) = 1;
I(29,25) = 1;
figure;
imshow(1-I,[],"InitialMagnification","fit");
ax = gca;
ax.Visible = "on";
ax.XTickLabel = [];
ax.YTickLabel = [];
ax.XTick = [7 24];
ax.XGrid = "on";
ax.YTick = [5 28];
ax.YGrid = "on";
title("Data");

手順 7: ブレのシミュレーション

以下の例は、ガウス フィルターである PSF を作成し、実際のイメージと畳み込むことで、星のイメージのブレをシミュレートします。

PSF = fspecial("gaussian",15,3);
Blurred = imfilter(I,PSF,"conv","sym");

ここで、星のイメージの部分のみ確認できるカメラをシミュレートします (ブレのみが確認されます)。重み関数配列 WT を作成します。これは、Blurred イメージの中心部分にある 1 (破線内にある "良好な" ピクセル) と、エッジにあるゼロ (信号を受信しない "無効" ピクセル) で構成されます。

WT = zeros(32);
WT(6:27,8:23) = 1;
CutImage = Blurred.*WT;

境界に関連付けられたリンギングを減らすには、PSF が指定された関数 edgetaper を適用します。

CutEdged = edgetaper(CutImage,PSF);
figure;
imshow(1-CutEdged,[],"InitialMagnification","fit");
ax = gca;
ax.Visible = "on";
ax.XTickLabel = [];
ax.YTickLabel = [];
ax.XTick = [7 24];
ax.XGrid = "on";
ax.YTick = [5 28];
ax.YGrid = "on";
title("Observed");

手順 8: 配列 WEIGHT の提供

アルゴリズムは、イメージの復元中に、重み配列に従って各ピクセル値を重み付けします。この例では、"無効" ピクセル値を最適化から除外すると同時に、中心のピクセルの値のみを使用します (ここでは WT = 1)。ただし、アルゴリズムは、カメラの視点のエッジを超えて、信号パワーをこれらの "無効" ピクセルの位置に入れることができます。解決された星の位置の精度に注意してください。

luc4 = deconvlucy(CutEdged,PSF,300,0,WT);
figure;
imshow(1-luc4,[],"InitialMagnification","fit");
ax = gca;
ax.Visible = "on";
ax.XTickLabel = [];
ax.YTickLabel = [];
ax.XTick = [7 24];
ax.XGrid = "on";
ax.YTick = [5 28];
ax.YGrid = "on";
title("Restored");

手順 9: 詳細にサンプリングされた PSF の提供

deconvlucy は、細かくサンプリングされた (subsample 倍細かい) PSF が指定された、アンダーサンプリングされたイメージを復元できます。不完全に解決されたイメージおよび PSF をシミュレートするために、以下の例では、各次元において、Blurred イメージと元の PSF を、2 つのピクセルを 1 つにして分類します。

Binned = squeeze(sum(reshape(Blurred,[2 16 2 16])));
BinnedImage = squeeze(sum(Binned,2));
Binned = squeeze(sum(reshape(PSF(1:14,1:14),[2 7 2 7])));
BinnedPSF = squeeze(sum(Binned,2));
figure;
imshow(1-BinnedImage,[],"InitialMagnification","fit");
ax = gca;
ax.Visible = "on";
ax.XTick = [];
ax.YTick = [];
title("Binned Observed");

アンダーサンプリングされた PSF である BinnedPSF を使用して、アンダーサンプリングされたイメージである BinnedImage を復元します。luc5 イメージは 3 つの星のみを区別することに注目してください。

luc5 = deconvlucy(BinnedImage,BinnedPSF,100);
figure;
imshow(1-luc5,[],"InitialMagnification","fit");
ax = gca;
ax.Visible = "on";
ax.XTick = [];
ax.YTick = [];
title("Poor PSF");

以下の例では、アンダーサンプリングされたイメージ (BinnedImage) を復元します。今回はより細かい PSF (subsample 倍細かいグリッド上で定義された PSF) を使用します。再構成イメージ (luc6) では、星の位置の精度が高くなります。イメージの右下隅にある 2 つの星の間でパワーを分散する方法に注目してください。これは、前の復元と同様に、1 つではなく、2 つの明るいオブジェクトが存在する可能性を示しています。

luc6 = deconvlucy(BinnedImage,PSF,100,[],[],[],2);
figure;
imshow(1-luc6,[],"InitialMagnification","fit");
ax = gca;
ax.Visible = "on";
ax.XTick = [];
ax.YTick = [];
title("Fine PSF");