Skip to content

GitLab

Commit e3916c41 authored by p.jhagema's avatar p.jhagema
Browse files

example script and data of a magnesium wire measured at PETRA III, P05 nano.... 

example script and  data of a magnesium wire measured at PETRA III, P05 nano. This is the supplement for the paper 'A phase retrieval framework to directly reconstruct the projected refractive index''
parent 88a13878
Pipeline #265803 passed with stage
in 53 seconds
Show whitespace changes
Inline Side-by-side
examples/MagnesiumWireRefractivePhaseRetrieval/algorithm_options.m 0 → 100755
View file @ e3916c41
%%
p.data_size = [2048 2048];
p.roiSize = [1700 1700]*1;
p.recSize = [2048 2048]*2;
p.positions = [519 521]-133; %mm slider positions, - sliderpos1 + defocus distance
% 4070 = 509 mm distance to goldi, 133 focal length
p.num_distances = 1;
% p.num_images = 1720;
p.z02 = 19.661e9;
p.lambda = 1.2398/11; %nm
p.det_pixel = 6500;
p.xcorr(1:4) = 51;
p.scale_fac = 1;
p = calc_Fresnel_geometry(p);
p.M./p.M(1)
p.xcorr
%%
p.fadeSet = fadeoutImage();
p.fadeSet.method = 'rectangle';
p.fadeSet.ellipseSize = [0.49 0.49];
p.fadeSet.windowShift = [0 0];
p.fadeSet.transitionLength = 20;
p.fadeSet.numSegments = 100;
% p.fadeSet.fadeToVal = 1;
%%
if(do_refractive) %refractive settings
%% magnitude
p.PMset = PFresnelMagnitude.defaultSettings();
p.PMset.projectionOrder = 'averagedRefractive';
% p.PMset.projectionOrder = 'cyclicRefractive';
% PMset.projectionOrder = 'sequential';
% tmp = true(p.data_size);
% tmp = padToSize(tmp, p.recSize, 0);
% p.PMset.applicationMask = repmat(tmp, [1 1 1]);
p.PMset.useGPU = 1;
p.PMset.singlePrecision = true;
p.PMset.representation = 'refractive';
p.PMset.doMErrors = 1;
p.PMset.plotIterates = 0;
p.PMset.plotInterval = 100;
%% refractive
p.refractiveSet = PRefractive.defaultSettings();
p.refractiveSet.ImagNegativeScaleFactor = 0;
p.refractiveSet.deltaOverBeta = [];
p.refractiveSet.minReal = -100;
p.refractiveSet.maxReal = 0;
p.refractiveSet.minImag = -0;
p.refractiveSet.maxImag = 1e-3;
p.refractiveSet.fwhmImag = 2;
p.refractiveSet.fwhmReal = 2;
p.refractiveSet.periodicBoundaries = 0;
p.refractiveSet.roiSize = p.roiSize*1.1;
p.refractiveSet.recSize = p.recSize;
%% support adaption
p.suppAdaptSet = PSupportAdaption.defaultSettings();
p.suppAdaptSet.doRefractive = 1;
p.suppAdaptSet.plotIterates = 1;
p.suppAdaptSet.plotInterval = 10;
p.suppAdaptSet.valueOutsideSupp = 0;
% p.suppAdaptSet.supportThreshold = 0.06;
p.suppAdaptSet.supportThreshold = 0.04;
p.suppAdaptSet.firstAdaption = 10;
p.suppAdaptSet.dilateValue = 5;
p.suppAdaptSet.erodeValue = 5;
% p.suppAdaptSet.
%% algorithm options
p.algSet = projectionAlgorithm.defaultSettings();
p.algSet.doRefractive = true;
p.algSet.doMErrors = 0;
p.algSet.roiSize = p.roiSize;
p.algSet.doConvergenceHistory = 1;
p.algSet.iterations = 2000;
p.algSet.useGPU = 1;
p.algSet.plotIterates = 1;
p.algSet.addMomentum = 0;
p.algSet.plotInterval = 100;
p.algSet.algorithm = 'NAG';%%
p.algSet.abortThreshold = -7.0;
p.algSet.eta = 1;
p.algSet.gamma = 1;
p.algSet.accelerationFrequency = 1;
end
%% amp/pha settings
if(~do_refractive)
%% magnitude
p.PMset = PFresnelMagnitude.defaultSettings();
p.PMset.projectionOrder = 'averaged';
p.PMset.useGPU = 1;
p.PMset.singlePrecision = true;
p.PMset.doMErrors = 1;
p.PMset.plotIterates = 0;
p.PMset.plotInterval = 10;
%% constraints
p.conSet = PConstraints.defaultSettings();
p.conSet.upperAmplitude = 1.0101;
p.conSet.lowerAmplitude = 0.553;
p.conSet.smoothAmplitudesBy = 2;
p.conSet.smoothPhasesBy = 2;
p.conSet.projectOnNegativity = 0;
p.conSet.unwrapPhases = 0;
%% support adaption
p.suppAdaptSet = PSupportAdaption.defaultSettings();
p.suppAdaptSet.plotIterates = 1;
p.suppAdaptSet.plotInterval = 50;
p.suppAdaptSet.valueOutsideSupp = 1;
p.suppAdaptSet.supportThreshold = 0.01;
p.suppAdaptSet.firstAdaption = 50;
p.suppAdaptSet.dilateValue = 5;
p.suppAdaptSet.erodeValue = 5;
%% algorithm options
p.algSet = projectionAlgorithm.defaultSettings();
p.algSet.doMErrors = 0;
p.algSet.roiSize = p.roiSize;
p.algSet.doConvergenceHistory = 1;
p.algSet.iterations = 2000;
p.algSet.useGPU = 1;
p.algSet.plotIterates = 1;
p.algSet.addMomentum = 1;
p.algSet.plotInterval = 100;
p.algSet.algorithm = 'NAG';%%
p.algSet.abortThreshold = -10.0;
p.algSet.eta = 1.1;
p.algSet.gamma = 0.9;
p.algSet.accelerationStartAt = 2;
p.algSet.accelerationFrequency = 1;
end
\ No newline at end of file
examples/MagnesiumWireRefractivePhaseRetrieval/calc_Fresnel_geometry.m 0 → 100644
View file @ e3916c41
% this function calculates the common parameters to calculate the Fresnel
% transformation for a multi distance propagation data set.
% Author: JH 58408 mJD
function p = calc_Fresnel_geometry(p)
if ~isfield(p, 'num_distances')
p.num_distances = numel(p.positions);
end
if ~isfield(p, 'scale_fac')
p.scale_fac = 1;
end
% initialize cell arrays
p.M = zeros(1, p.num_distances);
p.F_calc = zeros(1, p.num_distances);
p.z01 = zeros(1, p.num_distances);
p.z12 = zeros(1, p.num_distances);
p.z_eff = zeros(1, p.num_distances);
p.dxeff = zeros(1, p.num_distances);
if(~isfield(p, 'xcorr'))
p.xcorr = zeros(p.num_distances, 1);
end
% calculate parameters for each distance
for k=1:p.num_distances
p.z01(1, k) = (p.positions(1, k) + p.xcorr(1, k))*1e6;
p.z12(1, k) = p.z02-p.z01(1, k);
p.M(1, k) = p.z02/p.z01(1, k);
p.z_eff(1, k) = p.z12(1, k)/p.M(1, k);
p.dxeff(1, k) = (p.det_pixel* 1/p.scale_fac)/p.M(1, k);
p.F_calc(1, k)= p.dxeff(1, k).^2/(p.z_eff(1, k)*p.lambda);
end
p.F_rescaled = p.F_calc .* (p.M ./ max(p.M)).^2;
end
\ No newline at end of file
examples/MagnesiumWireRefractivePhaseRetrieval/eval_unwrap.m 0 → 100644
View file @ e3916c41
% wrapper function to evaluate phaseunwrapper.
% TO DO: add syntax for single argument function call
function [unwrap] = eval_unwrap(p_psi, a_psi, seed)
% p_psi = gather(angle(psi));
% a_psi = gather(abs(psi));
chunk_size = 2048;
num_steps = ceil(size(p_psi,1)/chunk_size);
num_bins = 2000;
max_num_workers = 16;
ppool = gcp('nocreate');
if(isempty('ppool'))
ppool = parpool(max_num_workers);
pctRunOnAll(maxNumCompThreads(2))
end
shifts = zeros(num_steps*num_steps, 4);
if num_steps > 1
for jj = 1:num_steps
for ii = 1:num_steps
shifts(ii + (jj-1)*num_steps, :) = ...
[ 1+(ii-1)*chunk_size, ...
(ii-1)*chunk_size+chunk_size,...
1+(jj-1)*chunk_size, ...
(jj-1)*chunk_size+chunk_size];
end
end
else
shifts(1, :) = [1 size(a_psi,1) 1 size(a_psi,2)];
end
target_roi = [1 chunk_size 1 chunk_size];
pha_psi = cell(size(shifts,1),1);
abs_psi = cell(size(shifts,1),1);
tmp = zeros(chunk_size);
if num_steps > 1
for ii = 1:size(shifts,1)
tmp = double(transplant_image_region(a_psi, tmp, shifts(ii,:), target_roi));
abs_psi{ii} = ((tmp));
tmp = double(transplant_image_region(p_psi, tmp, shifts(ii,:), target_roi));
pha_psi{ii} = ((tmp));
end
phase_tiles = cell(size(shifts,1),1);
if nargin < 3
seed = [chunk_size/2 chunk_size/2];
end
for ii = 1:size(shifts,1)
% phase_tiles{ii} = robustunwrap_cpp(seed,...
% pha_psi(shifts(ii,1):shifts(ii,2), shifts(ii,3):shifts(ii,4)),...
% abs_psi(shifts(ii,1):shifts(ii,2), shifts(ii,3):shifts(ii,4)), 2000);
phase_tiles{ii} = robustunwrap_cpp(seed,...
pha_psi{ii},...
abs_psi{ii}, num_bins);
end
phase_tiles = reshape(phase_tiles, num_steps, num_steps);
unwrap = cell2mat(phase_tiles);
else
seed = [size(p_psi,1)/2 size(p_psi,2)/2];
unwrap = robustunwrap_cpp(seed,...
double(p_psi),...
double(a_psi), num_bins);
end
end
\ No newline at end of file
examples/MagnesiumWireRefractivePhaseRetrieval/magnesium_wire_refractive_reco.m 0 → 100644
View file @ e3916c41
%% example for iterative refractive phaseretrieval
% This script is a supplement to the paper: "A phase retrieval framework to
% directly reconstruct the projected refractive index". Here we demonstrate the
% reconstruction of the holograms of the magnesium wire which have been obtained
% at P05.
%%
% change your current folder to the folder of this script
addpath(genpath('../../functions'));
addpath('.')
set(0,'DefaultFigureColor','w')
cmap = ((colormap('bone(512)')));
set(0,'DefaultFigureColormap', cmap);
% set defaults for figure
cmd = sprintf(['axis image; colorbar; set(gcf,''PaperPosition'',[0 0 5 5], ''PaperSize'', [5 5]);',...
' set(gca,''xtick'',[],''ytick'',[], ''XColor'', [0 0 0], ''YColor'', [0 0 0]);']);
set(groot,{'DefaultAxesXColor','DefaultAxesYColor','DefaultAxesZColor'},{'k','k','k'})
set(groot,'defaultAxesCreateFcn',@(ax,~)set(ax.Toolbar,'Visible','off'))
set(0,'DefaultImageCreateFcn', cmd)
set(groot,'defaultLineLineWidth',1.5)
set(0,'defaultAxesFontSize',15)
clear cmap cmd fig_cmd;
addToolsToToolbar();
close 1
%% ok lets look at the hologramm we are going to reconstruct
data = load('mg_wire_holo.mat');
holo = data.holo;
figure(1)
imagesc(holo)
colorbarSet = addColorbarLabel;
colorbarSet.Interpreter = 'latex';
addColorbarLabel('$I/I_0$')
% Here we see the flat-field corrected hologram padded to a size of
% 4096x4096 pixels. The rays outside are a 'continuation' of the hologram to
% mitigate gradients in the intensity. This square root of this array is
% used as input for the reconstruction algorithms
p.holo = holo;
% p is a general parameter struct
%% refractive reco
do_refractive = 1;
% algorithm_options is a seperate script which sets the parameters of the
% reconstruction. It takes care of calculating the propagation and algoritm specific
% parameters for the reconstruction. This populates p further.
algorithm_options
ref_reco = reconstruction_refractive(p);
%% plot reconstruction
% the reconstruction is ends with a projection on M, thus the sample
% constraints are not perfectly satisfied.
figure
imagesc(real(ref_reco.x))
addColorbarLabel('$\phi$ (rad)')
title('reconstructed \delta')
figure
imagesc(imag(ref_reco.x))
addColorbarLabel('$\beta$')
title('reconstructed \beta')
%% conventional reco
do_refractive = 0;
algorithm_options
conv_reco = reconstruction_conventional(p);
%% plot reconstruction
% the reconstruction is ends with a projection on M, thus the sample
% constraints are not perfectly satisfied.
figure
imagesc(angle(conv_reco.x))
addColorbarLabel('$\phi$ (rad)')
title('reconstructed phase angle')
% unwrap phases
tmp = conv_reco.x;
% tmp = pha
pha = (gpuArray(eval_unwrap(gather(cropToCenter(angle(tmp), [1700 1700])),...
gather(cropToCenter(abs(tmp), [1700 1700])), [850 850])));
figure()
imagesc(pha)
addColorbarLabel('$\phi$ (rad)')
title('reconstructed phase angle unwrapped')
figure
imagesc(abs(conv_reco.x))
addColorbarLabel('$A$')
title('reconstructed amplitude')
examples/MagnesiumWireRefractivePhaseRetrieval/reconstruction_conventional.m 0 → 100644
View file @ e3916c41
function reco = reconstruction_conventional(p)
% suppguess add the regions above and below the wire to the reconstruction.
% This is again for mitigation of boundary artifacts.
load( 'supp_guess.mat')
%% projectors for conventional reconstruchtion in amplitude and phase
guess = gpuArray(complex(ones(p.recSize.*p.scale_fac, 'single')));
PSuppAdapt = PSupportAdaption(gpuArray(p.supp), p.suppAdaptSet);
PM = PFresnelMagnitude(gpuArray(single(sqrt(p.holo))),...
p.F_calc, p.PMset);
PS = PConstraints(p.conSet) * PSuppAdapt;
%% 1st pass
starttime = datetime('now');
p.algSet.iterations = 1000;
projAlg = projectionAlgorithm(PM, PS, guess, p.algSet);
reco = projAlg.execute(p.algSet.iterations);
seq_conv_hist{1} = reco.stats.convergenceHistory;
%% we refinie the support from the final reconstruchtion
if(1)
p.supp = (angle(reco.x)) < -0.2;
p.supp = padToSize(cropToCenter(p.supp, p.data_size), p.recSize, 0);
se = strel('disk', 70);
p.supp = p.supp | imdilate(supp, se);
p.supp = bwareafilt(p.supp, 1, 'largest');
se = strel('disk', 10);
p.supp = imerode(p.supp, se);
se = strel('disk', 10);
p.supp = imdilate(p.supp, se);
p.supp = imfill(p.supp, 'holes');
p.supp = bwareafilt(p.supp, 1, 'largest');
p.supp = imgaussfilt(single(p.supp), 10/2.35);
figure(1002)
imagesc(p.supp)
end
%%
if(1)
PSupp = PSupport(gpuArray(p.supp), 1);
p.algSet.eta = 1.2;
p.algSet.gamma = 0.9;
% set or unset phase unwrapping, be warned it take time...
p.conSet.unwrapPhases = 1;
if(p.conSet.unwrapPhases)
p.conSet.projectOnNegativity = 1;
end
PS = PConstraints(p.conSet) * PSupp;% * PStat;
p.algSet.iterations = 1000;
projAlg = projectionAlgorithm(PM, PS, guess, p.algSet);
reco = projAlg.execute(p.algSet.iterations);
seq_conv_hist{2} = reco.stats.convergenceHistory;
end
%% done
starttime - datetime('now')
end
\ No newline at end of file
examples/MagnesiumWireRefractivePhaseRetrieval/reconstruction_refractive.m 0 → 100755
View file @ e3916c41
function reco = reconstruction_refractive(p)
% suppguess add the regions above and below the wire to the reconstruction.
% This is again for mitigation of boundary artifacts.
load('supp_guess.mat');
%% define the projectors for the algorithm
% Plain initialisation for the guess with zeros in imaginary and real part.
% This is equivalent to vaccum.
guess = gpuArray(complex(zeros(p.recSize.*p.scale_fac, 'single')));
% projectors in the sample plane: refractive constraints and self refining
% support
PRef = PRefractive(p.refractiveSet);
PSuppAdapt = PSupportAdaption(gpuArray(zeros(p.recSize)), p.suppAdaptSet);
% projectors can be concatuated
PS = PRef * PSuppAdapt;
% magnitude constraint
p.PMset.applicationMask = [];
PM = PFresnelMagnitude(gpuArray(single(sqrt(p.holo))), p.F_calc, p.PMset);
%% 1st pass - find a support
starttime = datetime('now');
p.algSet.iterations = 1000;
projAlg = projectionAlgorithm(PM, PS, guess, p.algSet);
reco = projAlg.execute(p.algSet.iterations);
%% we refinie the support from the final reconstruchtion
seq_conv_hist{1} = reco.stats.convergenceHistory;
if(1)
p.supp = real(reco.x) < -0.2;
p.supp = padToSize(cropToCenter(p.supp, p.data_size), p.recSize, 0);
se = strel('disk', 70);
p.supp = p.supp | imdilate(supp, se);
p.supp = bwareafilt(p.supp, 1, 'largest');
se = strel('disk', 80);
p.supp = imerode(p.supp, se);
se = strel('disk', 80);
p.supp = imdilate(p.supp, se);
p.supp = imfill(p.supp, 'holes');
p.supp = bwareafilt(p.supp, 1, 'largest');
p.supp = imgaussfilt(single(p.supp), 10/2.35);
figure(1002)
imagesc(p.supp)
end
%% 2nd pass - keep support fixed. just iterate
if(1)
% define support Projector with constant support
PSupp = PSupport(gpuArray(p.supp), 0);
% set refractive constraints
p.refractiveSet.minReal = -80;
p.refractiveSet.minImag = -0.5e-3;
p.refractiveSet.maxImag = 500e-3;
% parameters for nesterov
p.algSet.eta = 1.2;
p.algSet.gamma = 0.99;
PRef = PRefractive(p.refractiveSet);
PS = PRef * PSupp;% * PStat;
PM = PFresnelMagnitude(gpuArray(single(sqrt(p.holo))), p.F_calc, p.PMset);
p.algSet.iterations = 10000;
p.algSet.addMomentum = 1;
projAlg = projectionAlgorithm(PM, PS, guess, p.algSet);
reco = projAlg.execute(p.algSet.iterations);
seq_conv_hist{2} = reco.stats.convergenceHistory;
end
%% done
starttime - datetime('now')
end
\ No newline at end of file
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