Обратная свертка изображения

Ответ 1

В принципе, да. Просто преобразуйте оба изображения в частотное пространство с помощью FFT и разделите FFT изображения результата на изображение исходного изображения. Затем примените обратный БПФ, чтобы получить приближение ядра свертки.

Чтобы понять, почему это работает, обратите внимание, что свертка в пространственной области соответствует умножению в частотной области, и поэтому деконволюция аналогично соответствует делению в частотной области. В обычной деконволюции можно разделить БПФ свернутого изображения на изображение ядра, чтобы восстановить исходное изображение. Однако, поскольку свертка (например, умножение) является коммутативной операцией, роли ядра и источника могут быть произвольно обмениваться: свертывание источника ядром в точности совпадает с сверткой ядра исходным кодом.

Однако, как и другие ответы, это вряд ли даст точную реконструкцию вашего ядра по тем же причинам, что и обычная деконволюция, как правило, не будет точно воссоздать исходное изображение: округление и отсечение вызовут шум в процессе, и возможно, чтобы свертка полностью уничтожила некоторые частоты (умножив их на нуль), и в этом случае эти частоты не могут быть восстановлены.

Тем не менее, если ваше исходное ядро ​​имело ограниченный размер (поддержка), реконструированное ядро, как правило, должно иметь один острый пик вокруг начала координат, приближаясь к исходному ядру, окруженный низким уровнем шума. Даже если вы не знаете точный размер исходного ядра, не должно быть слишком сложно извлечь этот пик и отказаться от остальной части реконструкции.

Пример:

Здесь тестовое изображение Lenna в оттенках серого, уменьшено до 256 и раз; 256 пикселей и свернуто с ядром 5 и times; 5 в GIMP:

* →

Я использую Python с numpy/scipy для deconvolution:

from scipy import misc
from numpy import fft

orig = misc.imread('lena256.png')
blur = misc.imread('lena256blur.png')
orig_f = fft.rfft2(orig)
blur_f = fft.rfft2(blur)

kernel_f = blur_f / orig_f         # do the deconvolution
kernel = fft.irfft2(kernel_f)      # inverse Fourier transform
kernel = fft.fftshift(kernel)      # shift origin to center of image
kernel /= kernel.max()             # normalize gray levels
misc.imsave('kernel.png', kernel)  # save reconstructed kernel

Результирующее 256-кратное и 256-кратное изображение ядра и масштабирование области 7 и 7 пикселей вокруг его центра показаны ниже:

Сравнивая реконструкцию с исходным ядром, вы можете видеть, что они выглядят довольно похожими; действительно, применение обрезания в пределах от 0,5 до 0,68 к реконструкции приведет к восстановлению исходного ядра. Слабая рябь, окружающая пик в реконструкции, является шумом из-за округления и краевых эффектов.

Я не совсем уверен, что вызывает крестообразный артефакт, который появляется при реконструкции (хотя я уверен, что кто-то, у кого больше опыта с этими вещами, может сказать вам), но, с моей точки зрения, что он имеет какое-то отношение к краям изображения. Когда я свернул исходное изображение в GIMP, я сказал, чтобы он обрабатывал края, расширяя изображение (по существу копируя самые внешние пиксели), тогда как деконволюция FFT предполагает, что края изображения обтекают. Это вполне может ввести ложные корреляции вдоль осей x и y в реконструкции.

Ответ 2

Ну, оценка может быть получена, если известна верхняя граница размера матрицы свертки. Если это N, выберите N * N точек и попытайтесь решить систему линейных уравнений против коэффициентов свертки на основе данных источника и адресата. Учитывая округление цветовых компонентов, система не решит, но с линейным программированием вы сможете минимизировать общее смещение от ожидаемых значений путем небольших изменений этих коэффициентов.

Ответ 3

Вы можете попробовать выполнить deconvolution с исходным изображением в качестве ядра. Но результаты могут быть непредсказуемыми - деконволюция - очень неустойчивый процесс из-за шума, краевых эффектов, ошибок округления и т.д.

Ответ 4

Я переписал @Ilmari Karonen ответ на C/С++, используя fftw3 для кого-то, кто может найти его удобным:

Функция кругового сдвига

template<class ty>
void circshift(ty *out, const ty *in, int xdim, int ydim, int xshift, int yshift)
{
  for (int i =0; i < xdim; i++) 
  {
    int ii = (i + xshift) % xdim;
    for (int j = 0; j < ydim; j++) 
    {
      int jj = (j + yshift) % ydim;
      out[ii * ydim + jj] = in[i * ydim + j];
    }
  }
}

Теперь основной код

int width = 256;
int height = 256;

int index = 0;

MyStringAnsi imageName1 = "C://ka4ag.png";    
MyStringAnsi imageName2 = "C://KyPu2.png";

double * in1 = new double[width * height];
fftw_complex * out1 = new fftw_complex[width * height]; 

double * in2 = new double[width * height];
fftw_complex * out2 = new fftw_complex[width * height]; 

MyUtils::MyImage * im1 = MyUtils::MyImage::Load(imageName1, MyUtils::MyImage::PNG);
MyUtils::MyImage * im2 = MyUtils::MyImage::Load(imageName2, MyUtils::MyImage::PNG);

for (int i = 0; i < width * height; i++)
{
    in1[i] = ((im1->Get(i).r / (255.0 * 0.5)) - 1.0);
    in2[i] = ((im2->Get(i).r / (255.0 * 0.5)) - 1.0);
}


fftw_plan dft_plan1 = fftw_plan_dft_r2c_2d(width, height, in1, out1, FFTW_ESTIMATE);    
fftw_execute(dft_plan1);
fftw_destroy_plan(dft_plan1);

fftw_plan dft_plan2 = fftw_plan_dft_r2c_2d(width, height, in2, out2, FFTW_ESTIMATE);    
fftw_execute(dft_plan2);
fftw_destroy_plan(dft_plan2);

fftw_complex * kernel = new fftw_complex[width * height];   

for (int i = 0; i < width * height; i++)
{
    std::complex<double> c1(out1[i][0], out1[i][1]);
    std::complex<double> c2(out2[i][0], out2[i][1]);

    std::complex<double> div = c2 / c1;

    kernel[i][0] = div.real();
    kernel[i][1] = div.imag();
}

double * kernelOut = new double[width * height];

fftw_plan dft_planOut = fftw_plan_dft_c2r_2d(width, height, kernel, kernelOut, FFTW_ESTIMATE);
fftw_execute(dft_planOut);
fftw_destroy_plan(dft_planOut);

double * kernelShift = new double[width * height];

circshift(kernelShift, kernelOut, width, height, (width/2), (height/2));

double maxKernel = kernelShift[0];
for (int i = 0; i < width * height; i++)
{
    if (maxKernel < kernelShift[i]) maxKernel = kernelShift[i]; 
}

for (int i = 0; i < width * height; i++)
{
    kernelShift[i] /= maxKernel; 
}

uint8 * res = new uint8[width * height];
for (int i = 0; i < width * height; i++)
{                   
   res[i] = static_cast<uint8>((kernelShift[i]+ 1.0) * (255.0 * 0.5));
}

//now in res is similar result as in @Ilmari Karonen, but shifted by +128

В коде нет управления памятью, поэтому вы должны очистить свою память!

Ответ 5

Это классическая проблема деконволюции. То, что вы назвали матрицей свертки, обычно называется "ядром". Операцию свертки часто обозначают звездочкой '*' (не путать с умножением!). Используя это обозначение

Result = Source * Kernel

Ответы выше, используя FFT, являются правильными, но вы не можете использовать деконволюцию на основе FFT при наличии шума. Правильный способ сделать это - использовать деконволюцию Ричардсона-Люси (см. https://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution)

Это довольно просто реализовать. Этот ответ также дает пример реализации Matlab: Будет ли деконволюция Ричардсона-Люси работать для восстановления скрытого ядра?