image - 画像の逆畳み込み

Question

ソースと結果の画像があります。結果を得るためにソースでいくつかの畳み込み行列が使用されていることは知っています。この畳み込み行列は計算できますか? または、少なくとも正確なものではありませんが、非常に似ています。

Answer 1

原則として、はい。FFTを使用して両方の画像を周波数空間に変換し、結果画像の FFT をソース画像の FFT で割るだけです。次に、逆 FFT を適用して畳み込みカーネルの近似値を取得します。

これが機能する理由を理解するために、空間ドメインでの畳み込みは周波数ドメインでの乗算に対応し、デコンボリューションは周波数ドメインでの除算に同様に対応することに注意してください。通常のデコンボリューションでは、畳み込み画像の FFT をカーネルの FFT で割り、元の画像を復元します。ただし、畳み込み (乗算と同様) は可換演算であるため、カーネルとソースの役割を任意に交換できます。カーネルによるソースの畳み込みは、ソースによるカーネルの畳み込みとまったく同じです。

ただし、他の回答が指摘しているように、通常のデコンボリューションが一般に元の画像を正確に再構築しないのと同じ理由で、これはカーネルの正確な再構築をもたらす可能性は低いです:丸めとクリッピングはプロセスにノイズを導入し、それは可能です畳み込みが一部の周波数を完全に消去するため (それらにゼロを乗算することにより)、その場合、それらの周波数を再構築することはできません。

とはいえ、元のカーネルのサイズが制限されている (サポートされている) 場合、再構築されたカーネルは通常、原点の周りに 1 つの鋭いピークがあり、元のカーネルに近似し、低レベルのノイズに囲まれているはずです。元のカーネルの正確なサイズがわからない場合でも、このピークを抽出して残りの再構成を破棄するのは難しくありません。

例：

これはグレースケールのLenna テスト画像で、256×256 ピクセルに縮小され、GIMP で 5×5 カーネルで畳み込まれています。

※ →

デコンボリューションには numpy/scipy で Python を使用します。

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 軸に沿って疑似相関が生じる可能性があります。

Answer 2

これはデコンボリューションの古典的な問題です。畳み込み行列と呼ばれるものは、通常「カーネル」と呼ばれます。畳み込み演算は、多くの場合、星印 '*' で示されます (乗算と混同しないように!)。この表記の使用

Result = Source * Kernel

FFT を使用した上記の回答は正しいですが、ノイズが存在する場合は FFT ベースのデコンボリューションを実際に使用することはできません。それを行う正しい方法は、Richardson-Lucy デコンボリューションを使用することです ( https://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolutionを参照) 。

実装は非常に簡単です。この回答では、Matlab の実装例も提供しています。

Answer 3

畳み込み行列のサイズの上限がわかっている場合は、見積もりを作成できます。N の場合、N*N ポイントを選択し、ソースと宛先のデータに基づいて、畳み込み係数に対して線形方程式系を解こうとします。色成分の丸めを考えると、システムは解決しませんが、線形計画法を使用すると、これらの係数を少し変更することで、期待値からの合計オフセットを最小限に抑えることができます。

Answer 4

@Ilmari Karonen の回答を fftw3 を使用して C/C++ に書き直しました。

循環シフト機能

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

コードにはメモリ管理がないため、メモリをクリーンアップする必要があります。