import numpy as np
import scipy
from scipy.fftpack import ifftshift, fftshift, ifft2, fft2
from math import atan, atan2, ceil
# Utility functions for Ptychography experiment
[docs]def circ(X, Y, x0, y0, R):
# R radius of circular aperture
return heaviside(-((X-x0)**2 + (Y-y0)**2 - R**2))
[docs]def heaviside(M):
return (1 + np.sign(M))/2
[docs]def ellipsis(X, Y, x0, y0, a, b):
return heaviside(-(((X-x0)/a)**2 + ((Y-y0)/b)**2 - 1))
[docs]def prop_nf_pa(im, l, delta_z, dx, dy):
"""
Propagates a 2D-wave field within a paraxial approximation
into the optical near field.
The validity of the small-angle approximation is enforced.
If it is not obeyed, the function returns an error.
If sampling of the phase chirp function is not sufficient,
the field is embedded into a wider field of fiew, filling
new pixels, so that a continuous boundary conditions are
obeyed. If the field has been enlarged prior to propagation,
it is cut back to the original size after propagation.
Note that this cropping generally destroys the conservation
of the total fluence in the array.
Original Matlab file author K. Giewekemeyer
Original file from sftp://klaus@login.roentgen.physik.uni-goettingen.de/
home/klaus/Documents/Uni/Daten/Messzeiten/0909_BESSY/analysis_at_home/matlab/main/tools
last modified on 07.07.2010.
last modified and checked 13.10.2010.
Parameters
----------
im : array_like
Incident 2D wave field
l : number
Wave length
delta_z : number
Propagation distance
dx, dy : int
Pixel widths in horizontal and vertical direction
Returns
-------
array_like
Propagated wave field
"""
# TODO not tested
Nx = im.shape[1]
Ny = im.shape[0]
k = 2*np.pi/l
# Enforce paraxial approximation, according to a sufficient
# criterion due to Goodman, Introduction to Fourier Optics,
# 2005, p.68. Note that this criterion can be overly restrictive.
if delta_z == 0:
alpha_th = np.Inf
else:
alpha_th = (4/np.pi * l/delta_z)**0.25
alpha = max([atan(Nx/dx*(2*delta_z)), atan(Ny/dy*(2*delta_z))])
if alpha > alpha_th:
raise Exception('Small-angle approximation is violated')
# Sampling criterion for phase chirp.
px = abs(l*delta_z/(Nx*dx*dx))
py = abs(l*delta_z/(Ny*dy*dy))
f = max([px,py])
# beta = 1 ensures Nyquist sampling. Higher beta's make result better, but
# computationally more expensive.
beta = 1.0
if f > 1:
# This is not yet tested
print ('Probe-FOV is enlarged by factor %3.2f for propagation' % beta*f)
Nx_l = ceil(beta*f)*Nx
dqx = 2*np.pi/(Nx_l*dx)
Ny_l = ceil(beta*f)*Ny
dqy = 2*np.pi/(Ny_l*dy)
# If sampling criterion was obeyed in one direction, result gets even
# better.
# Again enforce paraxial approximation
alpha_th = (4/np.pi*l/delta_z)**0.25
alpha = max([atan(Nx_l/dx*(2*delta_z)), atan(Ny_l/dy*(2*delta_z))])
if alpha > alpha_th:
raise Exception('Small-angle approximation is violated')
rangeNx = np.arange(Nx_l)
rangeNy = np.arange(Ny_l)
Qx, Qy = np.meshgrid((rangeNx-(Nx_l//2))*dqx, (rangeNy-(Ny_l//2))*dqy)
# Paraxial propagator
kappa = -0.5*ifftshift(Qx**2 + Qy**2)/k
# Center in new array
Cx = (Nx_l//2)+1
Cy = (Ny_l//2)+1
# Embedding of image into larger field.
im_l = np.zeros(Ny_l, Nx_l)
# Filling up new pixels with boundary values (continuous boundary
# conditions).
# bottom rows
for i in range(Cy-(Ny//2)-1):
im_l[i,Cx-(Nx//2)-1:Cx+ceil(Nx/2.0)-1] = im[0,:]
# top rows
for i in range(Cy+ceil(Ny/2.0)-1,Ny_l):
im_l[i,Cx-(Nx//2)-1:Cx+ceil(Nx/2.0)-1] = im[-1,:]
# left columns
for i in range(Cx-(Nx//2)-1):
im_l[Cy-(Ny//2)-1:Cy+ceil(Ny/2.0)-1,i] = im[:,0]
# right columns
for i in range(Cx+ceil(Nx/2.0)-1,Nx_l):
im_l[Cy-(Ny//2)-1:Cy+ceil(Ny/2.0)-1,i] = im[:,-1]
# bottom left edge
im_l[0:Cy-(Ny//2)-1,0:Cx-(Nx//2)-1] = im[0,0]
# bottom right edge
im_l[0:Cy-(Ny//2)-1,Cx+ceil(Nx/2.0)-1:] = im[0,-1]
# top left edge
im_l[Cy+ceil(Ny/2.0)-1:Ny_l,0:Cx-(Nx//2)-1] = im[-1,0]
# top right edge
im_l[Cy+ceil(Ny/2.0)-1:Ny_l,Cx+ceil(Nx/2.0)-1:] = im[-1,-1]
# center: actual field
im_l[Cy-(Ny//2)-1:Cy+ceil(Ny/2.0)-1, Cx-(Nx//2)-1:Cx+ceil(Nx/2.0)-1] = im
# propagation
Im_l = fftshift(ifft2(fft2(ifftshift(im_l)) * np.exp(1j*kappa*delta_z)))
# decrease field of view
res = Im_l[Cy-(Ny//2)-1:Cy+ceil(Ny/2.0)-1, Cx-(Nx//2)-1:Cx+ceil(Ny/2.0)-1]
return res
else:
dqx = 2*np.pi/(Nx*dx)
dqy = 2*np.pi/(Ny*dy)
rangeNx = np.arange(Nx)
rangeNy = np.arange(Ny)
Qx, Qy = np.meshgrid((rangeNx-(Nx//2))*dqx, (rangeNy-(Ny//2)*dqy))
kappa = -0.5*ifftshift(Qx**2 + Qy**2)/k
res = fftshift(ifft2(fft2(ifftshift(im)) * np.exp(1j*kappa*delta_z)))
return res
[docs]def im2hsv(im, th):
"""
TODO: what this code does
Original Matlab file author K. Giewekemeyer (based on routine
pie_simulation_code.m by M. Dierolf et al.)
"""
abs_im = np.abs(im)
abs_im -= np.amin(abs_im)
abs_im /= np.amax(abs_im + np.spacing(1.))
imhsv = np.zeros((im.shape[0], im.shape[1], 3))
imhsv[:,:,0] = np.remainder(np.angle(im), 2*np.pi) / (2*np.pi)
imhsv[:,:,1] = 1
tmp = (1.0 + np.sign(abs_im-th)) / 2.0
tmp = np.clip((abs_im - tmp*abs_im + tmp*th) / th, 0, 1)
imhsv[:,:,2] = tmp
return imhsv