1

我发现这篇论文详细介绍了一种快速解决拉普拉斯/泊松问题的方法，尤其是对于图像。作者描述了一种非常简单但功能强大的算法，使他们能够比其他更复杂的算法更快地解决这些问题。这篇论文的关键见解是一个易于计算的离散格林函数的公式。

尽管论文中明确说明了算法的伪代码（第 8 页），但我在计算这个格林函数时遇到了麻烦。我在下面复制了一个示例问题。

论文通过使用卷积关系 A * B = Finv( F(A) ⊙ F(B)) 找到内核，其中⊙是Hadamard乘积。应用这个关系并注意到拉普拉斯核 K_L 和格林函数 V_mono 的卷积等于狄拉克的增量，格林函数 V_mono 可以从其他两个已知量计算出来。

但是我在推导格林函数时一直出错。例如，在下面的代码中，狄拉克 delta 的傅里叶变换中没有零项；但是在拉普拉斯核的傅里叶变换中有一个零项。这使得 Hadamard 将这些量除以找到未定义的格林函数。但是，当将拉普拉斯算子应用于已知图像并使用它在反向过程中找到格林函数时，结果几乎相同。（唯一的功能差异是有问题的 [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 等人命名为“使用绿色函数卷积进行梯度域图像编辑的快速和最优拉普拉斯求解器”。

4

0 回答 0