Gradient inversion with the mirror-flipped gradient operator¶
import numpy as np
import matplotlib
import matplotlib.pyplot as p
import math
%matplotlib inline
import scipy.ndimage.filters as filters
from scipy import misc
from numpy import pi, sqrt, cos, sin
matplotlib.rcParams['figure.figsize'] = (20.0, 10.0)
We will test our routine on a test image (a Matlab example image). I chose this image intentionally to be non-square, so that I make sure to get the sequence of dimensions right everywhere.
img = np.sum(misc.imread('cell.png'),2).astype(float) # cell is one of the Matlab demo images
# img = np.sum(misc.imread('cameraman.png'),2).astype(float) # another Matlab demo image
# img = np.sum(misc.imread('SheppLogan_Phantom_128.png'),2).astype(float) # cell is one of the Matlab demo images
img.shape, img.dtype
p.imshow(img,cmap=p.get_cmap('gray'))
#define the pixel size
[Nyimg,Nximg] = img.shape
print("Size of img: ",Nximg,"x",Nyimg,"pixels")
Next, lets define a mirror padded version of our input image. It will be this image that we are recovering. We do the mirror padding, in order to ensure periodic boundary conditions. This allows us to use FFTs to invert the gradient.
#define the pixel size
imgPad = np.lib.pad(img, ((0, Nyimg),(0,Nximg)), 'symmetric')
[Ny,Nx] = imgPad.shape
print("Size of imgPad: ",Nx,"x",Ny,"pixels")
p.imshow(imgPad,cmap=p.get_cmap('gray')), p.colorbar()
# define the real space grid
grid2 = np.mgrid[0:Ny:1,0:Nx:1]
x = grid2[1]
y = grid2[0]
dx = x[0][1]-x[0][0]
dy = y[1][0]-y[0][0]
print("Pixel size: ",dx,dy)
We will do the gradient computation in reciprocal space, but in the forward direction this can certainly also be done more quickly in real space. But in order to make sure that the forward and inverse calculation correspond to one another, we need to ensure that we define the gradient filters as follows: $$\frac{\partial f}{\partial x} = \frac{f(x+dx,y)-f(f-dx,y)}{2 dx} = IFT\left[\tilde{f}\left( k_x,k_y\right) \cdot i \sin(k_x) \right]$$ $$\frac{\partial f}{\partial y} = \frac{f(x,y+dy)-f(f,y-dy)}{2 dy} = IFT\left[\tilde{f}\left( k_x,k_y\right) \cdot i \sin(k_y) \right]$$ where $\tilde{f}\left( k_x,k_y\right)$ is the Fourier transform of $f(x,y)$ so that $$f(x,y) = \sum_{x=1}^{N_x} \sum_{y=1}^{Ny} \tilde{f}\left( k_x,k_y\right) \exp\left[ik_x x+ik_y y \right] $$
Using these definitions we define the following differentiating operator $$ D[f] = \left(\frac{\partial}{\partial x} + i\frac{\partial}{\partial y} \right)f(x,y) = IFT\left[\tilde{f}\left( k_x,k_y\right) \cdot \left(i \sin(k_x) -\sin(k_y) \right) \right] = g(x,y)$$ The integrating inverse operator of $D[f]$ is then $$ D^{-1}[g] = \left(\frac{\partial}{\partial x} + i\frac{\partial}{\partial y} \right)^{-1}g(x,y) = \frac{1}{N_x N_y}IFT\left[\frac{\tilde{g}\left( k_x,k_y\right)}{i \sin(k_x) -\sin(k_y)} \right]$$ where again $\tilde{f}\left( k_x,k_y\right)$ is the Fourier transform of $g(x,y)$
# construct reciprocal space
kx1 = np.fft.fftshift(2*pi*np.fft.fftfreq(Nx,dx))
ky1 = np.fft.fftshift(2*pi*np.fft.fftfreq(Ny,dy))
kx = np.tile(kx1.reshape(1,Nx),(Ny,1))
ky = np.tile(ky1.reshape(Ny,1),(1,Nx))
# compute the sin(kx,ky) because this is the reciprocal space filter that corresponds to
# the finite difference derivative operator df/dx = [f(x+dx)-f(x-dx)]/2dx
kxs = np.sin(kx)
kys = np.sin(ky)
# finally we need to define k=isin(kx)-sin(ky) and ik=1/k
k = (1j*kxs-kys)
kmod = sqrt(kx**2+ky**2)
ksmod = kxs**2+kys**2
# ik = np.where(ksmod > 1e-8,(-1j*kxs-kys)/ksmod,1)
ik = np.where(ksmod > 1e-8,1/k,1)
# we also want to apply some slight low-pass filtering in the inversion step,
# in order to not enhance noise
kmax = 0.5*np.amax(kmod)
# ik = np.where(kmod>kmax, ik*0.5*(1+np.cos(pi*(kmod-kmax)/(np.amax(kmod)-kmax))), ik)
# Now, let's define the functions D and D^{-1}:
# Since we have padded these images, we can safely wrap around the edges.
def D(f):
fmean = np.mean(f)
Df = 0.5*(np.roll(f,1,axis=1)-np.roll(f,-1,axis=1))+1j*0.5*(np.roll(f,1,axis=0)-np.roll(f,-1,axis=0))
# Just to convince ourselves that above math is correct, we can, alternatively
# also do this calculation using Fourier transforms (will be a bit slower though):
# D = np.fft.ifft2(np.fft.fft2(f)*k)
return Df, fmean
def iD(g,gmean):
iDg = np.fft.fft2(g)*ik
iDg[0][0] = gmean*Nx*Ny
iDg = np.real(np.fft.ifft2(iDg))
return iDg
The above operator can now also be used to compute the second order derivatives!
# Test the above routine:
[imgD,imgMean] = D(imgPad)
print("Mean of image: ",imgMean)
p.subplot(2,2,1)
p.imshow(np.real(imgD),cmap=p.get_cmap('gray')), p.colorbar(), p.title("df/dx")
p.subplot(2,2,2)
p.imshow(np.imag(imgD),cmap=p.get_cmap('gray')), p.colorbar(), p.title("df/dy")
print("Fraction of nonzero elements in dfdx: ",np.count_nonzero(np.real(imgD))/(Nx*Ny))
print("Fraction of nonzero elements in dfdy: ",np.count_nonzero(np.imag(imgD))/(Nx*Ny))
# Test the inversion !!!
D2 = iD(imgD,imgMean)
p.subplot(2,2,3)
p.imshow(D2,cmap=p.get_cmap('gray')), p.colorbar(), p.title("recovered from gradient")
p.subplot(2,2,4)
p.imshow(imgPad,cmap=p.get_cmap('gray')), p.colorbar(), p.title("original")
# Construct a symmetrized version of the gradient for computing the second derivative
imgDs = np.lib.pad(imgD[0:Nyimg,0:Nximg], ((0, Nyimg),(0,Nximg)), 'symmetric')
# Now apply the operator a second time
[imgD2s,gradMean] = D(imgDs)
print("Mean of image: ",gradMean)
p.subplot(1,2,1)
p.imshow(np.real(imgD2s),cmap=p.get_cmap('gray')), p.colorbar(), p.title("real part")
p.subplot(1,2,2)
p.imshow(np.imag(imgD2s),cmap=p.get_cmap('gray')), p.colorbar(), p.title("imaginary part")
# And this is what happens if we do not symmetrize the gradient beforehand
[imgD2,gradMean] = D(imgD)
print("Mean of image: ",gradMean)
p.subplot(1,2,1)
p.imshow(np.real(imgD2),cmap=p.get_cmap('gray')), p.colorbar(), p.title("real part")
p.subplot(1,2,2)
p.imshow(np.imag(imgD2),cmap=p.get_cmap('gray')), p.colorbar(), p.title("imaginary part")
