Давайте посмотрим на код в вопросе для одного канала, предполагая, что img - это изображение в оттенках серого - здесь все может применяться для каждого канала, поэтому нам не нужно повторять все три раза:
for i in range(25):
img = ndimage.convolve(img, kernel)
img = img / 16 # normalize
Мы доберемся до удаления свертки через минуту.Сначала давайте упростим применяемую операцию.
Упростим обработку
Вышеуказанное идентично (с точностью до цифры):
kernel = kernel / 16 # Normalize
for i in range(25):
img = ndimage.convolve(img, kernel)
Это верно до тех пор, пока img не является целочисленным типом, где происходит отсечение и / или округление.В общем, с * сверткой и C некоторой константой
g = C (f * h) = f * (C h)
Далее мы знаем, что применение свертки 25 раз аналогично применению свертки один раз с составным ядром
g = (((f * h) * h) * h) * h = f * (h * h * h * h)
Как мы можем получить составное ядро?Применение свертки к изображению, которое является всеми нулями и с 1 в среднем пикселе, возвращает ядро снова, так что
delta = np.zeros(kernel.shape)
delta[delta.shape[0]//2, delta.shape[1]//2] = 1
kernel2 = ndimage.convolve(delta, kernel)
kernel2 == kernel # is true everywhere, up to numerical precision
Таким образом, следующий код находит ядро, которое используется для сглаживания изображения ввопрос:
kernel = np.array([[1.0,2.0,1.0],
[2.0,4.0,2.0],
[1.0,2.0,1.0]]) / 16
delta = np.zeros(((kernel.shape[0]-1)*25+1, (kernel.shape[1]-1)*25+1))
delta[delta.shape[0]//2, delta.shape[1]//2] = 1
for i in range(25):
delta = ndimage.convolve(delta, kernel)
kernel = delta
Это ядро будет очень похоже на ядро Гаусса из-за центральной предельной теоремы.
Теперь мы можем получить тот же результат, что и в вопросе с одной сверткой:
output = ndimage.convolve(img, kernel)
Инвертирование свертки
Процесс обратной фильтрации называется деконволюцией.В теории это очень тривиальный процесс, но на практике он очень сложен из-за шума, неточного знания ядра и т. Д.
Мы знаем, что можем вычислить свертку через область Фурье:
output = np.convolve(img, kernel, mode='wrap')
совпадает с
output = np.real(np.fft.ifft2( np.fft.fft2(img) * np.fft.fft2(np.fft.ifftshift(kernel)) ))
(при условии, что kernel соответствует размеру img, нам, как правило, сначала нужно заполнить его нулями).Любые различия между результатами работы в пространственной и частотной областях вызваны тем, как изображение расширяется за его границы при использовании convolve.Метод Фурье предполагает периодическое граничное условие, поэтому я использовал режим 'wrap' для свертки здесь.
Обратная операция - это просто деление в области Фурье:
img = np.real(np.fft.ifft2( np.fft.fft2(output) / np.fft.fft2(np.fft.ifftshift(kernel)) ))
Чтобы это работало, нам нужно знать точные значения kernel, и в процессе не должно быть шума.Для output, вычисленного, как указано выше, это должно дать точный теоретический результат
Однако некоторые ядра могут быть точно равны нулю для некоторых частотных компонентов (то есть np.fft.fft2(np.fft.ifftshift(kernel)) содержит нули).Эти частоты не могут быть восстановлены, и деление на 0 приведет к значениям NaN, которые распространятся по всему изображению в обратном преобразовании, обратное изображение будет всеми NaN.
Для ядра Гаусса нет нулей,так что этого не должно случиться.Однако будет много частот, которые почти равны нулю.Поэтому преобразование Фурье output также будет иметь очень малое значение для этих элементов.Обратный процесс - это деление очень маленького значения на другое очень маленькое значение, что приводит к проблемам с числовой точностью.
И вы можете видеть, как этот процесс, если будет очень мало шума, будет сильноусилить этот шум, так что выходной сигнал почти полностью задается этим шумом.
Деконволюция Винера включает регуляризацию для предотвращения этих проблем с шумом и неточностью чисел.По сути, вы предотвращаете деление на очень маленькие числа, добавляя положительное значение к преобразованию Фурье kernel. Википедия имеет хорошее описание деконволюции Винера.
Демо
Я использую MATLAB с DIPimage 3, чтобы сделать быструю демонстрацию (гораздо меньше усилий для меня, чем стрельбадо Python и выяснить, как все это сделать там).Это код:
% Get data
a = readim('https://i.stack.imgur.com/OfSx2.png');
a = a{1}; % Keep only red channel
% Create kernel
kernel = [1.0,2.0,1.0
2.0,4.0,2.0
1.0,2.0,1.0] / 16;
tmp = deltaim((size(kernel)-1)*25+1);
for ii=1:25
tmp = convolve(tmp,kernel,'add zeros');
end
kernel = tmp;
% Apply convolution
b = convolve(a,kernel,'periodic');
% Find inverse operation
% 1- pad stuff so image and kernel have the same size
% we maintain the periodic boundary condition for image b
b = repmat(b,ceil(imsize(kernel)./imsize(b)));
kernel = extend(kernel,imsize(b));
% 2- apply something similar to Wiener deconvolution
c = ift(ft(b)/(ft(kernel)+1e-6),'real'); % Not exactly Wiener, but there's no noise!
% 3- undo padding
c = cut(c,imsize(a),'top left');
Это вывод, верхняя треть - входное изображение, средняя треть - размытое изображение, нижняя треть - выходное изображение:
Здесь важно отметить, что я использовал периодическое граничное условие для начальной свертки, которое соответствует тому, что происходит в преобразовании Фурье.Другие граничные условия будут вызывать артефакты в обратном преобразовании вблизи краев.Поскольку размер ядра больше, чем у образа, весь образ будет одним большим артефактом, и вы ничего не сможете восстановить.Также обратите внимание, что, чтобы заполнить ядро нулями до размера изображения, мне пришлось повторить изображение, так как ядро больше, чем изображение.Репликация изображения снова соответствует периодическому граничному условию, наложенному преобразованием Фурье.Оба этих трюка можно было бы проигнорировать, если бы входное изображение было намного больше, чем ядро свертки, как и следовало ожидать в нормальной ситуации.
Также обратите внимание, что без регуляризации в деконволюции выводом является весь NaNпотому что мы делим очень маленькие значения на очень маленькие значения.Преобразование Фурье ядра содержит много почти нулей, потому что размытие довольно сильное.
Наконец, обратите внимание, что добавление даже небольшого количества шума к размытому изображению сделает невозможным деконволюциюизображение таким образом, что текст может быть прочитан.Обратное преобразование будет выглядеть очень хорошо, но текстовые штрихи будут искажены настолько, что буквы больше не будут легко распознаваться:
Приведенный выше код использует DIPimage 3, который еще не имеет официального двоичного файла для установки, его необходимо собрать из source .Для запуска кода с использованием DIPimage 2.xa необходимо внести несколько изменений:
Граничное условие должно быть установлено с помощью dip_setboundary, вместо того, чтобы передавать его непосредственно в convolveфункция.Строки 'add zeros' и 'periodic' являются граничными условиями.
Функции ft и ift используют симметричную нормализацию, каждая из которых умножает свой вывод на 1/sqrt(prod(imsize(image))), тогда какв DIPimage 3 нормализация - это более распространенное умножение на 1/prod(imsize(image)) для ift и 1 для ft.Это означает, что преобразование Фурье kernel должно быть умножено на sqrt(prod(imsize(kernel))), чтобы соответствовать результату DIPimage 3:
c = real(ift(ft(b)/((ft(kernel)*sqrt(prod(imsize(kernel))))+1e-6)));