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

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

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

Алгоритм взвешивает каждое значение пикселя в соответствии с массивом 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');

Шаг 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');