Deblurring Images Using the Lucy-Richardson Algorithm

本例展示如何使用Lucy-Richardson演算法去模糊影象。 當點擴散函式PSF（模糊運算子）是已知的時候，它可以被有效地使用，但是很少或沒有資訊可用於噪聲。 模糊和噪聲影象通過迭代，加速，衰減的Lucy-Richardson演算法恢復。 您可以使用光學系統的特性作為輸入引數來提高影象恢復的質量。

Step 1: Read Image

該示例讀取RGB影象並將其裁剪為256×256×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');

Step 2: Simulate a Blur and Noise

模擬由於相機運動或對焦不足而可能模糊的真實影象。 由於隨機干擾，影象也可能會產生噪音。 該示例通過將高斯濾波器與真實影象進行卷積來模擬模糊（使用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');

Step 3: Restore the Blurred and Noisy Image

恢復提供PSF的模糊和噪聲影象，並僅使用5次迭代（預設值為10）。 輸出是與輸入影象相同型別的陣列。

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

Step 4: Iterate to Explore the Restoration

每次迭代產生的影象都會改變。 要調查影象恢復的演變過程，您可以按照以下步驟進行反褶積：執行一組迭代，檢視結果，然後從停止的位置恢復迭代。 為此，輸入影象必須作為單元陣列的一部分傳遞。 例如，通過傳遞{BlurredNoisy}而不是BlurredNoisy作為輸入影象引數來啟動第一組迭代。

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

在這種情況下，輸出luc1_cell變成單元陣列。 單元輸出包含四個數值陣列，其中第一個是BlurredNoisy影象，第二個是恢復的類double影象，第三個陣列是第一個倒數第一個迭代的結果，第四個陣列是內部引數 的迭代集合。 輸出單元陣列的第二個陣列陣列影象luc1_cell {2}與步驟3的輸出陣列（影象luc1）完全相同，可能除了它們的類別（單元格輸出始終會給出類別double的恢復映像）。

要恢復迭代，請從前一個函式呼叫的輸出中取出單元陣列luc1_cell，並將其傳遞給去卷積函式。 使用預設的迭代次數（NUMIT = 10）。 恢復的影象是總共15次迭代的結果。

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

Step 5: Control Noise Amplification by Damping

最新的影象luc2是15次迭代的結果。 雖然它比5次迭代的早期結果更清晰，但影象呈現出“斑點”外觀。 這些斑點不符合任何實際結構（將其與真實影象進行比較），而是將資料中的噪音過於緊密地結合起來。

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

本示例的下一部分使用模擬的星形影象（為了簡單和快速）探討了解卷函式的WEIGHT和SUBSMPL輸入引數。

Step 6: Create Sample Image

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');

Step 7: Simulate a Blur

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

現在模擬一個只能觀察部分恆星影象的相機（只能看到模糊）。 建立一個加權函式陣列WEIGHT，其中包含模糊影象中心部分（位於虛線內的“好”畫素）和邊緣零點（“壞”畫素 - 那些不接收訊號的畫素）。

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

為了減少與邊界相關的振鈴，請將邊界限制器功能應用於給定的PSF。

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');

Step 8: Provide the WEIGHT Array

該演算法在恢復影象時根據WEIGHT陣列對每個畫素值進行加權。 在我們的例子中，只使用中心畫素的值（WEIGHT = 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');

Step 9: Provide a finer-sampled PSF

在給定更精細的取樣PSF（更細SUBSMPL時間）的情況下，deconvlucy可以恢復欠取樣影象。 為了模擬較差解析的影象和PSF，該示例在每個維度中將Blurred影象和原始PSF（一箇中的兩個畫素）分組。

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（在SUBSMPL次更精細的網格上定義）。 重建影象（luc6）可以更準確地解析恆星的位置。 請注意它如何在影象右下角的兩顆星星之間分配力量。 這提示存在兩個明亮的物體，而不是一個，就像在以前的修復中那樣。

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