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

У меня есть исходное и результирующее изображение. Я знаю, что некоторая матрица свертки была использована на источнике, чтобы получить результат. Можно ли вычислить эту матрицу свертки ? Или, по крайней мере, не точная, но очень похожая.

5 13
algorithm convolution

5 ответов:

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

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

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

Пример:

Вот тестовое изображение Ленны в оттенках серого, уменьшенное до 256×256 пикселей и свернутое с ядром 5×5 в GIMP:

Оригинал ∗ Ядро → Результат

Я буду использовать Python с numpy/scipy для деконволюции: 

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 в реконструкция.

20

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

1

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

1

Я переписал @Ilmari Karonen answer на C / 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
Код не имеет никакого управления памятью, поэтому вы должны очистить свою память !
1

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

Result = Source * Kernel

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

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

1