In []: 
from scipy import ndimage
import scipy
In [61]: 
def _airy_func(rr, amplitude=1.0, width=1.0):
    """
    For a simple radially symmetric airy function, returns the value at a given
    (normalized) radius
    """
    return amplitude * (2.0 * scipy.special.j1(rr/width) / (rr/width))**2
In [60]: 
che=ndimage.imread('../figures/CheHigh.jpg',flatten=True)
imshow(che,cmap='gray')
Out[60]:
<matplotlib.image.AxesImage at 0x7f5e6f989890>
In [85]: 
X=arange(-50,50)*1.0+1.e-5
x,y=array(meshgrid(X,X))

rr=hypot(x,y)
psf=_airy_func(rr,width=2.0)

figure(figsize=(10,5))
subplot(121), imshow(psf,cmap='gray')
subplot(122), plot(psf[50,:])
Out[85]:
(<matplotlib.axes.AxesSubplot at 0x7f5e6df4a8d0>,
 [<matplotlib.lines.Line2D at 0x7f5e6decd650>])
In [79]: 
che_bis=fftconvolve(che,psf)

figure(figsize=(10,5))
subplot(121), imshow(che,cmap='gray')
subplot(122), imshow(che_bis,cmap='gray')
Out[79]:
(<matplotlib.axes.AxesSubplot at 0x7f5e6e733110>,
 <matplotlib.image.AxesImage at 0x7f5e6e6a7cd0>)
In [80]: 
figure(figsize=(10,5))
subplot(121), plot(che[2000,:])
subplot(122), plot(che_bis[2000,:])
Out[80]:
(<matplotlib.axes.AxesSubplot at 0x7f5e6e612d90>,
 [<matplotlib.lines.Line2D at 0x7f5e6e5240d0>])
In [82]: 
psf=_airy_func(rr,width=10.0)

figure(figsize=(10,5))
subplot(121)
imshow(psf)
subplot(122)
plot(psf[50,:])
Out[82]:
[<matplotlib.lines.Line2D at 0x7f5e6e2b9bd0>]
In [83]: 
che_bis=fftconvolve(che,psf)

figure(figsize=(10,5))
subplot(121), imshow(che,cmap='gray')
subplot(122), imshow(che_bis,cmap='gray')
Out[83]:
(<matplotlib.axes.AxesSubplot at 0x7f5e6e1ffb50>,
 <matplotlib.image.AxesImage at 0x7f5e6e181750>)
In [84]: 
figure(figsize=(10,5))
subplot(121), plot(che[2000,:])
subplot(122), plot(che_bis[2000,:])
Out[84]:
(<matplotlib.axes.AxesSubplot at 0x7f5e6e0edf50>,
 [<matplotlib.lines.Line2D at 0x7f5e6dff2b10>])