I found this paper detailing a method for fast solution of Laplace/Poisson problems, especially for images. The authors describe a very simple yet powerful algorithm allowing them to solve these problems orders of magnitude faster than other, more complex algorithms. The key insight of the paper is a formulation of a discrete Green's function which is easily computable.
Despite pseudocode for the algorithm being explicitly stated in the paper (pg. 8), I am having trouble computing this Green's function. I have reproduced an example problem below.
The paper finds the kernel by using the convolution relationship A * B = Finv( F(A) ⊙ F(B)), where ⊙ is the Hadamard product. Applying this relationship and noting that the convolution of the Laplacian kernel K_L and the Green's function V_mono equals Dirac's delta, the Green's function V_mono can be computed from the other two known quantities.
But I keep getting an error in deriving the Green's function. For example, in the code below, there are no zero entries in the Fourier transform of Dirac's delta; but there is one zero entry in the Fourier Transform of the Laplacian kernel. This makes Hadamard division of these quantities to find Green's function undefined. However, when applying the Laplacian to a known image, and using that to find the Green's function in the reverse process, the results are nearly equivalent. (The only functional difference is the problematic [0,0] entry has a 10^18 entry in the reverse computed entry, while the analytically computed one has an infinite value there.)
Since this is the crux of the paper, I am unable to apply the reconstruction algorithm. How might I go about fixing my algorithm?
#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)))
P.S. For citation purposes, the paper is named 'Fast and Optimal Laplacian Solver for Gradient-Domain Image Editing using Green Function Convolution' by Beaini et al.
