Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

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

この例では、ルーシー・リチャードソン アルゴリズムを使用してイメージのブレを除去する方法を示します。このアルゴリズムは、点像分布関数 (PSF) (ブレの演算子) が既知であり、ノイズについて入手できる情報がわずかであるか、まったくない場合にも効果的に使用できます。ブレとノイズを含むイメージは、減衰を考慮した高速のルーシー・リチャードソン アルゴリズムを反復することで復元します。光学システムの特性を、イメージの復元の品質向上に役立つ入力パラメーターとして使用できます。

手順 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");

Figure contains an axes object. The axes object with title Original Image contains 3 objects of type image, text.

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

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

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

Figure contains an axes object. The axes object with title Blurred contains an object of type image.

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

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

Figure contains an axes object. The axes object with title Blurred & Noisy contains an object of type image.

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

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

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

Figure contains an axes object. The axes object with title Restored Image, NUMIT = 5 contains an object of type image.

手順 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");

Figure contains an axes object. The axes object with title Restored Image, NUMIT = 15 contains an object of type image.

手順 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");

Figure contains an axes object. The axes object with title Restored Image with Damping, NUMIT = 15 contains an object of type image.

この例の以下の部分では、(簡略化と高速化のために) シミュレーションを実行した星のイメージを使用して、関数 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");

Figure contains an axes object. The axes object with title Data contains an object of type image.

手順 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");

Figure contains an axes object. The axes object with title Observed contains an object of type image.

手順 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");

Figure contains an axes object. The axes object with title Restored contains an object of type image.

手順 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");

Figure contains an axes object. The axes object with title Binned Observed contains an object of type image.

アンダーサンプリングされた 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");

Figure contains an axes object. The axes object with title Poor PSF contains an object of type image.

以下の例では、アンダーサンプリングされたイメージ (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");

Figure contains an axes object. The axes object with title Fine PSF contains an object of type image.

参考

| | |

関連するトピック