Шаг 1: Readimage

Пример читает в изображении 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");

Шаг 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: выполните итерации, чтобы исследовать восстановление

Получившиеся изменения образа с каждой итерацией. Чтобы исследовать эволюцию восстановления изображений, можно сделать развертку на шагах: сделайте набор итераций, смотрите результат, и затем возобновите итерации от того, где они были остановлены. Для этого входное изображение должно быть передано как часть массива ячеек. Например, запустите первый набор итераций путем передачи в {BlurredNoisy} вместо BlurredNoisy как входной параметр изображения.

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

В этом случае выход, luc1_cell, становится массивом ячеек. Ячейка выход состоит из четырех числовых массивов, где первым является BlurredNoisy отобразите, вторым является восстановленное изображение класса double, третий массив является результатом one-last итерации, и четвертый массив является внутренним параметром выполненного с помощью итераций набора. Второй числовой массив выходного массива ячеек, отобразите luc1_cell{2}, идентично выходному массиву изображения Шага 3, luc1, за возможным исключением их класса (ячейка выход всегда дает восстановленное изображение класса double).

Чтобы возобновить итерации, возьмите выход из предыдущего вызова функции, массив ячеек 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");

Следующая часть этого примера исследует weight и subsample введите параметры deconvlucy функция, с помощью симулированного звездообразного изображения (для простоты & скорости).

Шаг 6: создайте демонстрационное изображение

Пример создает черное/белое изображение четырех звезд.

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 изображение ("хорошие" пиксели, расположенные в пунктирных линиях) и нули в ребрах ("плохие" пиксели - те, которые не получают сигнал).

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

Чтобы уменьшать вызов, сопоставленный с границами, примените edgetaper функция с данным 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");

Шаг 8: обеспечьте МАССИВ ВЕСОВ

Веса алгоритма каждое пиксельное значение согласно массиву весов при восстановлении изображения. В нашем примере только используются значения центральных пикселей (где 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 может восстановить субдискретизируемое изображение, учитывая более прекрасный произведенный PSF (более прекрасный subsample \times. Симулировать плохо разрешенное изображение и 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");

Восстановите субдискретизируемое изображение, BinnedImage, с помощью субдискретизируемого PSF, BinnedPSF. Заметьте, что 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- времена более прекрасная сетка). Восстановленное изображение (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");