Commit e3b08731 authored by s.gretchko's avatar s.gretchko
Browse files

Adding Ptychography Experiment - Part I: loading data

parent 1697d123
import SetProxPythonPath
from proxtoolbox.experiments.ptychography.ptychographyExperiment import PtychographyExperiment
Pty = PtychographyExperiment(debug = False)
from proxtoolbox.Utilities.cell import Cell, isCell
class DataProcessor:
Abstract class for data processors.
The role of the data processor is to serve as an interface
between a data set and the algorithm that will work on the
data. In essence, it contains all the data that is relevent
for a given experiment. It provides the initial iterate for
the algorithm. Derived classes need to implement the loadData()
method, which is called during the Experiment's initialization
def __init__(self): = None
self.data_sq = None
self.norm_data = None
self.norm_data_sq = None
self.rt_data = None
self.norm_rt_data = None
self.u0 = None
self.truth = None
self.truth_dim = None
self.norm_truth = None
self.product_space_dimension = None
self.sets = None # number of sets
def loadData(self, experiment):
Load or generate the dataset that will be used for
a given experiment. Create the initial iterate.
experiment: an instance of the Experiment class
Experiment object that will use this data processor. The
experiment object contains all the parameters that
describe a given experiment.
raise NotImplementedError("This is just an abstract interface")
def reshapeData(self, Nx, Ny, Nz):
Reshape data based on the given arguments. This method is called
during the initialization of the Experiment object, after the call
to loadData().
Nx, Ny, Nz: int
The new dimensions to be used.
This is a default implementation that is used so far for phase experiments
TODO: too many data fields... See if we can reduce the number of data fields
Could use a variable that tells us which kind of data we are dealing with
# reshape and rename the data
if self.data_sq is None:
# need to put data into required format
self.data_sq = = self.rt_data
if self.norm_data is not None:
self.norm_data_sq = self.norm_data
self.norm_data = self.norm_rt_data
# The following lines of code correspond to the Cone_and_Sphere branch. However,
# the case where the data is a cell is not treated (in Matlab implementation)
# TODO: Need to check is this is OK.
if not isCell(
tmp =
if tmp[0] == 1 or tmp[1] == 1:
self.data_sq = self.data_sq.reshape(Nx, Ny)
# the prox algorithms work with the square root of the measurement: =, Ny)
from proxtoolbox.experiments.dataProcessor import DataProcessor
from proxtoolbox import algorithms
from proxtoolbox import proxoperators
from proxtoolbox.proxoperators.proxoperator import ProxOperator
......@@ -77,11 +76,12 @@ class Experiment(metaclass = ExperimentMetaClass):
Nx = None,
Ny = None,
Nz = None,
noise = 'none',
noise = None,
algorithm = None,
numruns = 1,
MAXIT = 1000,
TOL = 1e-8,
TOL2 = None,
lambda_0 = 0.85, # Initial relaxation value for first iterations
lambda_max = 0.85, # Final relaxation value for later iterations
lambda_switch = 20, # Timescale over which to change relaxation parameter
......@@ -151,10 +151,11 @@ class Experiment(metaclass = ExperimentMetaClass):
self.numruns = numruns
self.TOL = TOL
self.TOL2 = TOL2
self.lambda_0 = lambda_0
self.lambda_max = lambda_max
self.lambda_switch = lambda_switch
self.save_output = save_output
self.json_output = json_output
self.matlab_output = matlab_output
......@@ -188,6 +189,9 @@ class Experiment(metaclass = ExperimentMetaClass):
self.inverse_propagator = None
self.amplitude = None #TODO Quick fix for now. Needed for P_amp prox operator
self.run_number = 1
self.debug = debug
def initialize(self):
......@@ -226,7 +230,8 @@ class Experiment(metaclass = ExperimentMetaClass):
self.propagator = getattr(proxoperators, self.propagator)
self.inverse_propagator = getattr(proxoperators, self.inverse_propagator)
self.TOL2 = 1e-20
if self.TOL2 is None:
self.TOL2 = 1e-20
# TODO in Cone_and_Sphere estimation of the parameters above is done differently (set_fattener)
# New matlab code:
This diff is collapsed.
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
def circ(X, Y, x0, y0, R):
# R radius of circular aperture
return heaviside(-((X-x0)**2 + (Y-y0)**2 - R**2))
def heaviside(M):
return (1 + np.sign(M))/2
def ellipsis(X, Y, x0, y0, a, b):
return heaviside(-(((X-x0)/a)**2 + ((Y-y0)/b)**2 - 1))
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 s
last modified on 07.07.2010.
last modified and checked 13.10.2010.
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
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
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
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
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment