Skip to content
GitLab
Commit 9e9d77b7 authored by alexander.dornheim's avatar alexander.dornheim
Browse files

Updated phase 

Correct experiment in Near_field_Siemens_experimental_in is Krueger
parent 7f1b912c
Hide whitespace changes
Inline Side-by-side
proxtoolbox/Problems/Phase/Goettingen_data_processor.py
View file @ 9e9d77b7
from pathlib import Path
def Goettingen_data_processor(config):
# Parameters of the forward problem
......@@ -22,171 +24,172 @@ def Goettingen_data_processor(config):
experiment = config['experiment']
if (experiment == 'Krueger'):
# Sample: NTT-AT, ATN/XREESCO-50HC
# 500 nm thick Ta structure: amplitude transmission 0.9644,
# phase shift 0.4rad at E = 17.5 keV
# parameters see below
# empty waveguide (WG) beam
WG = np.loadtxt('../InputData/Phase/WG_beam.mat')
# hologram
print('Loading data hologram_not-normalized.mat')
I_exp = np.loadtxt('../InputData/Phase/hologram_not-normalized.mat')
I_exp[np.isnan(I_exp)] = 1
# The following is NOT used because it is not clear how the ``normalized" hologram
# is obtained. I suspect it is obtained by simply dividing out the empty beam
# data in the obervation plane. But this is incorrect. Despite the
# efforts of the Hohage team to get around this error by ad hoc regularization,
# we take the following exact approach: back propagate the empty beam, and
# divide this from the unknown object in the object plane. This approach is true
# to the imaging model and does not introduce any approximations. It does, however,
# expose the fact that we do not know the phase of the beam.
# clear I_exp
# load 'data/hologram.mat'
## single image reconstruction
# total number of elements
N = I_exp.size
#number of pixels
Ny = I_exp.shape[0]; Nx = I_exp.shape[1]
config['Ny']=Ny; config['Nx']=Nx
##########################
# Experimental parameters
##########################
# energy in keV
E = 17.5
# wavelength [m]
lambd = 12.398/E*1E-10
k = 2*pi/lambd
# distance source-sample [m]
z1 = 7.48e-3
# distance sample-detector [m]
z2 = 3.09
# effective distance of detector
z_eff = z1*z2/(z1+z2)
# effective magnification factor
M = (z1+z2)/z1
#pixel size in detector plane
dx_det = 55e-6
dy_det = 55e-6
# magnified coordinate system in detector plane
# [X_det,Y_det] = meshgrid(dx_det*[0:1:Nx-1],dy_det*[0:1:Ny-1]);
# demagnifiedpixel size in sample plane
dx = dx_det/M
dy = dy_det/M
# coordinate system in sample plane
# [X,Y] = meshgrid(dx*((1:1:Nx)-floor(Nx/2)-1),dy*((1:1:Ny)-floor(Ny/2)-1));
# magnified coordinate system in detector plane
# [X_det,Y_det] = meshgrid(dx_det*((1:1:Nx)-floor(Nx/2)-1),dy_det*((1:1:Ny)-floor(Ny/2)-1));
# grid conversion in q-space
dqx = 2*pi/(Nx*dx)
dqy = 2*pi/(Ny*dy)
[Qx,Qy] = np.meshgrid(dqx*(range(Nx)-floor(Nx/2)-1),dqy*(range(Ny)-floor(Ny/2)-1))
##################################
# Prepare data from ProxToolbox:
##################################
# Fresnel propagator:
config['use_farfield_formula']=false
config['Nx']=Nx
config['Ny']=Ny
config['fresnel_nr'] = 1*2*pi*Nx # Fresnel number
config['magn']=1 # magnification factor (should this be M above, or
# is it okay since the demagnified pixel size is used?)
kappa = sqrt(k**2-(Qx**2+Qy**2))
config['FT_conv_kernel'] = fftshift(exp(-1j*kappa*z_eff))
config['beam'] = abs(ifft2(fft2(sqrt(WG*7.210116465530644e-04))/config['FT_conv_kernel']))
# the scaling factor of the waveguide array above was
# reverse engineered from the scaling that was apparently
# used in the preparation of the data in the file hologram.mat
# I don't want to use this hologram as the data, preferring instead
# to keep the experimental data as close to the actual observation
# as possible. Also, I am disagreeing with what my colleagues in
# physics think is the correct way to compensate for a structured
# empty beam. According to my reverse engineering, it appears that
# they have divided the EMPTY BEAM MEASUREMENT WG from the object
# IN THE OBJECT PLANE. It is true that you need to do the empty beam
# correction in the object plane, though not with the empty beam, but
# rather the reverse propagated empty beam, as I have done above. The
# reason being that the empty beam is measured in the measurement plane,
# so it does not make sense to divide the object&beam in the object
# plane by the empty beam in the measurement plane - you should divide
# by the empty beam propagated back to the object plane. The difficulty in
# this is that we do not know the phase of the empty beam. The above
# formulation sets the phase of the beam to zero at the object plane.
# The beam is saved as the square root of the empty beam to conform with
# the projeciton operations working on the amplitude instead of the
# magnitude.
# problem data:
config['data'] = I_exp
config['rt_data'] = sqrt(I_exp)
config['data_zeros'] = config['data']==0
config['norm_rt_data'] =sqrt(sum(config['data']))
# use the abs_illumination field to represent the
# support constraint.
config['abs_illumination'] = np.ones(size(I_exp))
config['support_idx'] = config['abs_illumination~=0'].nonzero
config['product_space_dimension'] = 1
config['supp_phase'] = []
# start values
config['u_0']=ifft2(fft2(config['rt_data']*config['magn'])/config['FT_conv_kernel'])/config['beam']
'''
figure(1);
subplot(2,2,1)
imagesc(Nx*dx*linspace(0,1,Ny)*1E6,Nx*dx*linspace(0,1,Nx)*1E6,config['beam);
axis equal;
axis tight;
xlabel('x [\mum]');
ylabel('y [\mum]');
colormap(gray); colorbar; title('empty beam - detector plane')
subplot(2,2,2)
imagesc(Nx*dx*linspace(0,1,Ny)*1E6,Nx*dx*linspace(0,1,Nx)*1E6,config['data);
axis equal;
axis tight;
xlabel('x [\mum]');
ylabel('y [\mum]');
colormap(gray); colorbar; title('uncorrected near field observation')
subplot(2,2,3)
imagesc(Nx*dx*linspace(0,1,Ny)*1E6,Nx*dx*linspace(0,1,Nx)*1E6,(abs(config['u_0)));
axis image; colormap(gray); colorbar;
xlabel('x [\mum]');
ylabel('y [\mum]');
title('initial guess amplitude');
# # pattern
subplot(2,2,4);
imagesc(Nx*dx*linspace(0,1,Ny)*1E6,Nx*dx*linspace(0,1,Nx)*1E6,(angle(config['u_0))');
axis image; colormap(gray); colorbar;
xlabel('x [\mum]');
ylabel('y [\mum]');
title('initial guess phase');
caxis([-0.9 -0.4]);
'''
if (experiment == 'Krueger'):
print('here')
# Sample: NTT-AT, ATN/XREESCO-50HC
# 500 nm thick Ta structure: amplitude transmission 0.9644,
# phase shift 0.4rad at E = 17.5 keV
# parameters see below
# empty waveguide (WG) beam
WG = np.loadtxt('../InputData/Phase/WG_beam.mat')
# hologram
print('Loading data hologram_not-normalized.mat')
I_exp = np.loadtxt('../InputData/Phase/hologram_not-normalized.mat')
I_exp[np.isnan(I_exp)] = 1
# The following is NOT used because it is not clear how the ``normalized" hologram
# is obtained. I suspect it is obtained by simply dividing out the empty beam
# data in the obervation plane. But this is incorrect. Despite the
# efforts of the Hohage team to get around this error by ad hoc regularization,
# we take the following exact approach: back propagate the empty beam, and
# divide this from the unknown object in the object plane. This approach is true
# to the imaging model and does not introduce any approximations. It does, however,
# expose the fact that we do not know the phase of the beam.
# clear I_exp
# load 'data/hologram.mat'
## single image reconstruction
# total number of elements
N = I_exp.size
#number of pixels
Ny = I_exp.shape[0]; Nx = I_exp.shape[1]
config['Ny']=Ny; config['Nx']=Nx
##########################
# Experimental parameters
##########################
# energy in keV
E = 17.5
# wavelength [m]
lambd = 12.398/E*1E-10
k = 2*pi/lambd
# distance source-sample [m]
z1 = 7.48e-3
# distance sample-detector [m]
z2 = 3.09
# effective distance of detector
z_eff = z1*z2/(z1+z2)
# effective magnification factor
M = (z1+z2)/z1
#pixel size in detector plane
dx_det = 55e-6
dy_det = 55e-6
# magnified coordinate system in detector plane
# [X_det,Y_det] = meshgrid(dx_det*[0:1:Nx-1],dy_det*[0:1:Ny-1]);
# demagnifiedpixel size in sample plane
dx = dx_det/M
dy = dy_det/M
# coordinate system in sample plane
# [X,Y] = meshgrid(dx*((1:1:Nx)-floor(Nx/2)-1),dy*((1:1:Ny)-floor(Ny/2)-1));
# magnified coordinate system in detector plane
# [X_det,Y_det] = meshgrid(dx_det*((1:1:Nx)-floor(Nx/2)-1),dy_det*((1:1:Ny)-floor(Ny/2)-1));
# grid conversion in q-space
dqx = 2*pi/(Nx*dx)
dqy = 2*pi/(Ny*dy)
[Qx,Qy] = np.meshgrid(dqx*(range(Nx)-floor(Nx/2)-1),dqy*(range(Ny)-floor(Ny/2)-1))
##################################
# Prepare data from ProxToolbox:
##################################
# Fresnel propagator:
config['use_farfield_formula']=false
config['Nx']=Nx
config['Ny']=Ny
config['fresnel_nr'] = 1*2*pi*Nx # Fresnel number
config['magn']=1 # magnification factor (should this be M above, or
# is it okay since the demagnified pixel size is used?)
kappa = sqrt(k**2-(Qx**2+Qy**2))
config['FT_conv_kernel'] = fftshift(exp(-1j*kappa*z_eff))
config['beam'] = abs(ifft2(fft2(sqrt(WG*7.210116465530644e-04))/config['FT_conv_kernel']))
# the scaling factor of the waveguide array above was
# reverse engineered from the scaling that was apparently
# used in the preparation of the data in the file hologram.mat
# I don't want to use this hologram as the data, preferring instead
# to keep the experimental data as close to the actual observation
# as possible. Also, I am disagreeing with what my colleagues in
# physics think is the correct way to compensate for a structured
# empty beam. According to my reverse engineering, it appears that
# they have divided the EMPTY BEAM MEASUREMENT WG from the object
# IN THE OBJECT PLANE. It is true that you need to do the empty beam
# correction in the object plane, though not with the empty beam, but
# rather the reverse propagated empty beam, as I have done above. The
# reason being that the empty beam is measured in the measurement plane,
# so it does not make sense to divide the object&beam in the object
# plane by the empty beam in the measurement plane - you should divide
# by the empty beam propagated back to the object plane. The difficulty in
# this is that we do not know the phase of the empty beam. The above
# formulation sets the phase of the beam to zero at the object plane.
# The beam is saved as the square root of the empty beam to conform with
# the projeciton operations working on the amplitude instead of the
# magnitude.
# problem data:
config['data'] = I_exp
config['rt_data'] = sqrt(I_exp)
config['data_zeros'] = config['data']==0
config['norm_rt_data'] =sqrt(sum(config['data']))
# use the abs_illumination field to represent the
# support constraint.
config['abs_illumination'] = np.ones(size(I_exp))
config['support_idx'] = config['abs_illumination~=0'].nonzero
config['product_space_dimension'] = 1
config['supp_phase'] = []
# start values
config['u_0']=ifft2(fft2(config['rt_data']*config['magn'])/config['FT_conv_kernel'])/config['beam']
'''
figure(1);
subplot(2,2,1)
imagesc(Nx*dx*linspace(0,1,Ny)*1E6,Nx*dx*linspace(0,1,Nx)*1E6,config['beam);
axis equal;
axis tight;
xlabel('x [\mum]');
ylabel('y [\mum]');
colormap(gray); colorbar; title('empty beam - detector plane')
subplot(2,2,2)
imagesc(Nx*dx*linspace(0,1,Ny)*1E6,Nx*dx*linspace(0,1,Nx)*1E6,config['data);
axis equal;
axis tight;
xlabel('x [\mum]');
ylabel('y [\mum]');
colormap(gray); colorbar; title('uncorrected near field observation')
subplot(2,2,3)
imagesc(Nx*dx*linspace(0,1,Ny)*1E6,Nx*dx*linspace(0,1,Nx)*1E6,(abs(config['u_0)));
axis image; colormap(gray); colorbar;
xlabel('x [\mum]');
ylabel('y [\mum]');
title('initial guess amplitude');
# # pattern
subplot(2,2,4);
imagesc(Nx*dx*linspace(0,1,Ny)*1E6,Nx*dx*linspace(0,1,Nx)*1E6,(angle(config['u_0))');
axis image; colormap(gray); colorbar;
xlabel('x [\mum]');
ylabel('y [\mum]');
title('initial guess phase');
caxis([-0.9 -0.4]);
'''
proxtoolbox/Problems/Phase/Near_field_Siemens_experimental_in.py
View file @ 9e9d77b7
......@@ -37,7 +37,7 @@ new_config = {
## What type of measurements are we working with?
## Options are: 'single diffraction', 'diversity diffraction',
## 'ptychography', and 'complex'
'experiment' : 'Kruger',
'experiment' : 'Krueger',
## Next we move to things that most of our users will know
## better than we will. Some of these may be overwritten in the
......
proxtoolbox/Problems/Phase/phase.py
View file @ 9e9d77b7
......@@ -34,20 +34,24 @@ class Phase(Problem):
module = __import__(self.config['data_filename'])
data_processor = getattr(module, self.config['data_filename'])
data_processor(self.config)
#reshape and rename the data
self.config['data_sq'] = self.config['data'];
self.config['data'] = self.config['rt_data'];
tmp = self.config['data'].shape;
if(tmp[0]==1 or tmp[1]==1):
self.config['data_sq'] = self.config['data_sq'].reshape((self.config['Nx'],self.config['Ny']));
#the projection algorithms work with the square root of the measurement:
self.config['data'] = self.config['data'].reshape((self.config['Nx'],self.config['Ny']));
self.config['normM']=self.config['norm_rt_data']; #previously (in matlab) norm_data
if 'Nz' not in self.config:
self.config['Nz'] = 1;
# reshape and rename the data
if('data_sq' in self.config):
pass # already put data into required format, skip
else:
self.config['data_sq'] = self.config['data']
self.config['data'] = self.config['rt_data']
if('norm_data' in data):
config['norm_data_sq']= config['norm_data']
self.config['norm_data']=self.config['norm_rt_data']
if(tmp[0]==1 or tmp[1]==1):
self.config['data_sq'] = self.config['data_sq'].reshape((self.config['Nx'],self.config['Ny']))
#the prox algorithms work with the square root of the measurement:
self.config['data'] = self.config['data'].reshape((self.config['Nx'],self.config['Ny']))
if 'Nz' not in self.config:
self.config['Nz'] = 1
#If method_config[formulation is does not exist, i.e. not specified in
......
Supports Markdown
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