特に画像のラプラス/ポアソン問題を高速に解決する方法を詳述しているこの論文を見つけました。著者は、他のより複雑なアルゴリズムよりも桁違いに速くこれらの問題を解決できる、非常に単純でありながら強力なアルゴリズムについて説明しています。この論文の重要な洞察は、簡単に計算できる離散グリーン関数の定式化です。
アルゴリズムの疑似コードが論文に明示的に記載されているにもかかわらず ( 8 ページ)、このグリーン関数の計算に問題があります。以下に問題の例を再現しました。
この論文では、畳み込み関係 A * B = Finv( F(A) ⊙ F(B)) を使用してカーネルを見つけます。ここで、⊙ はアダマール積です。この関係を適用し、ラプラシアン カーネル K_L とグリーンの関数 V_mono の畳み込みがディラックのデルタに等しいことに注意すると、グリーンの関数 V_mono は他の 2 つの既知の量から計算できます。
しかし、グリーン関数を導出する際にエラーが発生し続けます。たとえば、次のコードでは、ディラックのデルタのフーリエ変換にゼロ エントリはありません。しかし、ラプラシアン カーネルのフーリエ変換にはゼロ エントリが 1 つあります。これにより、これらの量のアダマール除算が行われ、未定義のグリーン関数が検出されます。ただし、ラプラシアンを既知の画像に適用し、それを使用して逆のプロセスでグリーン関数を見つけると、結果はほぼ同等になります。(唯一の機能上の違いは、問題のある [0,0] エントリには逆計算されたエントリに 10^18 エントリがあるのに対し、分析的に計算されたエントリには無限の値があることです。)
これがこの論文の核心であるため、再構成アルゴリズムを適用することはできません。アルゴリズムを修正するにはどうすればよいですか?
#Import Test Image
import numpy as np
from scipy import misc
from scipy import fftpack
from scipy import signal
import matplotlib.pyplot as plt
img = misc.face(gray=True)#misc.face()[:,:,0]
pad=4
def pad_with(vector, pad_width, iaxis, kwargs):
pad_value = kwargs.get('padder', 10)
vector[:pad_width[0]] = pad_value
vector[-pad_width[1]:] = pad_value
img=np.pad(img, pad, pad_with, padder=0)#img is padded to allow for proper integration later
#Establish Dirac Kernel
# Dirac Kernel Specified by Paper
kernelPaper=np.zeros(img.shape)
kernelPaper[1,1]=1
# Dirac Kernel Specified by Link
kernelDirac=np.asarray([[0,0,0],[0,1,0],[0,0,0]],dtype="float")
#kernel = np.outer(signal.gaussian(2, 3), signal.gaussian(2, 3))
sz = (img.shape[0] - kernelDirac.shape[0], img.shape[1] - kernelDirac.shape[1]) # total amount of padding
kernelDirac = np.pad(kernelDirac, (((sz[0]+1)//2, sz[0]//2), ((sz[1]+1)//2, sz[1]//2)), 'constant')
kernelDirac = fftpack.ifftshift(kernelDirac)
#Establish Laplace Kernel
# Laplace Kernel Specified by Paper
kernelLaplacePaper=np.zeros(img.shape)
kernelLaplacePaper[np.tile(range(3),(3,1)),np.transpose(np.tile(range(3),(3,1)))]=np.asarray([[0,-1,0],[-1,4,-1],[0,-1,0]],dtype="float")
# Laplace Kernel Specified by Link
kernelLaplace=np.asarray([[0,-1,0],[-1,4,-1],[0,-1,0]],dtype="float")
#kernel = np.outer(signal.gaussian(2, 3), signal.gaussian(2, 3))
sz = (img.shape[0] - kernelLaplace.shape[0], img.shape[1] - kernelLaplace.shape[1]) # total amount of padding
kernelLaplace = np.pad(kernelLaplace, (((sz[0]+1)//2, sz[0]//2), ((sz[1]+1)//2, sz[1]//2)), 'constant')
kernelLaplace = fftpack.ifftshift(kernelLaplace)
#Establish Test Case: Apply Laplacian to Known Image
filtered = np.real(fftpack.ifft2(fftpack.fft2(img) * fftpack.fft2(kernelLaplace)))
##plt.imshow(filtered, vmin=0, vmax=255)
##plt.show()
# ************Desired Result: Green's Function *************
# Green's Function Specified by Paper
green_F_Paper = fftpack.fft2(kernelDirac)/fftpack.fft2(kernelLaplace) #This is actually F(Green's Function)
# Green's Function Specified by Inverse Comparison
green_F_ComputeThroughReversal = fftpack.fft2(img)/fftpack.fft2(filtered)
print("Difference between paper Green's Function kernel and Green's Function computed using original image: ")
print(np.sum(abs(green_F_ComputeThroughReversal-green_F_Paper)))
print("Difference between paper Green's Function kernel and Green's Function computed using original image without [0,0] entry: ")
greenDiff=green_F_ComputeThroughReversal-green_F_Paper
greenDiff[0,0]=0
print(np.sum(abs(greenDiff)))
# Laplace Kernel Specified by Inverse Comparison
kernelLaplace_ComputeThroughReversal = fftpack.ifft2(fftpack.fft2(kernelDirac)/green_F_ComputeThroughReversal)
print("Difference between paper Laplacian kernel and formulation found online: ")
print(np.sum(abs(kernelLaplace_ComputeThroughReversal-kernelLaplace)))
#Desired Result: Solving Laplacian
Ic=(fftpack.ifft2(fftpack.fft2(filtered) * fftpack.fft2(green_F_ComputeThroughReversal))).real #checking to see if reverse operation is preserved
Ic=Ic-Ic[0,0]
#Ic=Ic[pad:-pad,pad:-pad]
print("Difference between solved Laplacian and original image using kernels computed using original image: ")
print(np.sum(abs(Ic-img)))
Ic_Paper=(fftpack.ifft2(fftpack.fft2(filtered) * fftpack.fft2(green_F_Paper))).real #checking to see if reverse operation is preserved
Ic_Paper=Ic_Paper-Ic_Paper[0,0]
#Ic=Ic[pad:-pad,pad:-pad]
print("Difference between reconstructed image and original image using paper formulation: ")
print(np.sum(abs(Ic_Paper-img)))
PS 引用の目的で、この論文は、Beaini et al. による「グリーン関数畳み込みを使用した勾配領域画像編集のための高速かつ最適なラプラシアン ソルバー」と名付けられています。