Пусть е быть исходное изображение, г быть фильтруют один, а ч фильтр применяется к е, так что:
f * h = g
Передача тха t в частотную область:
F.H = G, so H = G/F
Проблема в том, что инвертирование F ОЧЕНЬ чувствительно к шуму.
Как осуществить это в MATLAB:
close all;
f = imread('cameraman.tif');
[x,y] = size(f);
figure,imshow(f);
h = fspecial('motion', 20, 40); % abitrary filter just for testing the algorithm
F = fft2(f);
H = fft2(h,x,y);
G = F.*H;
g = ifft2(G); % the filtered image
figure, imshow(g/max(g(:)));
% Inverting the original image
epsilon = 10^(-10);
small_values = find(abs(F)<epsilon);
F(small_values) = epsilon;
F_i = ones(x,y)./F;
H_calculated = G.*F_i;
h_calculated = ifft2(H_calculated);
% remove really small values to try to infer the original size of h
r = sum(h_calculated,1)<epsilon;
c = sum(h_calculated,2)<epsilon;
h_real = h_calculated(~r,~c);
% Calculate error
% redo the filtering with the found filter
figure,g_comp = ifft2(fft2(f).*fft2(h_real,x,y));
imshow(g_comp/max(g_comp(:)));
rmse = sqrt(mean(mean((double(g_comp) - double(g)).^2,2),1))
редактировать: Просто объяснить эпсилон часть:
Это может быть, что некоторые значения в F равны нулю, или очень близко к нулю. Если мы попытаемся инвертировать F этими малыми значениями, у нас будут проблемы с бесконечностью. Простым способом решения этого является усечение каждого значения в F, которое меньше сколь угодно малого предела, epsilon на коде.
Математически, что было сделано это:
For all F < epsilon, F = epsilon
Если у вас есть обработанное изображение и исходное изображение, то это очень просто определить коэффициенты фильтра, и оттуда можно классифицировать ядро фильтра. –
да, у меня есть оба, но как я могу определить коэффициенты? с матлабом? – NaN
Вы можете deconvolve, взяв БПФ обоих изображений, разделите, а затем возьмите обратный БПФ результата. См. http://www.mathworks.com/matlabcentral/fileexchange/5465-fast-deconvolution –