Шаг 1: Чтение изображения

Пример читается в изображении RGB и высевает его равным 256 256 3. The 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 image, second - восстановленное изображение класса double, третий массив - результат итерации «один перед последним», а четвертый массив - внутренний параметр итерационного набора. Второй числовой массив выхода массива ячеек, изображение luc1_cell{2}, идентичен выход массиву шага 3, изображение luc1, за возможным исключением их класса (выход камеры всегда дает восстановленное изображение класса double).

Чтобы возобновить итерации, возьмите выход от предыдущего вызова функции, cell-array 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 должен быть того же класса, что и входное изображение. Алгоритм демпфирует изменения в модели в областях, где различия малы по сравнению с шумом. The DAMPAR используемый здесь равен 3 стандартным отклонениям шума. Заметьте, что изображение более плавное.

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

Следующая часть этого примера исследует WEIGHT и SUBSMPL входные параметры функции 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');

Теперь симулируйте камеру, которая может наблюдать только часть изображений звезд (видно только размытие). Создайте массив весовых функций WEIGHT, который состоит из таковых в центральной части изображения 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: Обеспечьте мАССИВ ВЕСОВ

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

Шаг 9: Обеспечьте более мелкую выборку PSF

deconvlucy может восстановить недостаточно дискретизированное изображение, учитывая более мелкую выборку PSF (более мелкую по времени SUBSMPL). Чтобы симулировать плохо разрешенное изображение и 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 (заданной на более мелкой сетке 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');